Line data Source code
1 : !!****m* ABINIT/m_gstate
2 : !! NAME
3 : !! m_gstate
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, JYR, MKV, MT, FJ, MB, DJA)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : ! nvtx related macro definition
22 : #include "nvtx_macros.h"
23 :
24 : module m_gstate
25 :
26 : use defs_basis
27 : use defs_rectypes
28 : use m_errors
29 : use m_xmpi
30 : use m_abicore
31 : use libxc_functionals
32 : use m_exit
33 : use m_crystal
34 : use m_scf_history
35 : use m_abimover
36 : use m_wffile
37 : use m_rec
38 : use m_efield
39 : use m_ddb
40 : use m_bandfft_kpt
41 : use m_gemm_nonlop_projectors
42 : use m_xg_nonlop
43 : use m_wfk
44 : use m_nctk
45 : use m_hdr
46 : use m_ebands
47 : use m_dtfil
48 : use m_extfpmd
49 : use m_rcpaw
50 : use m_alloc_hamilt_gpu
51 :
52 : use defs_datatypes, only : pseudopotential_type
53 : use defs_abitypes, only : MPI_type
54 : use m_time, only : timab
55 : use m_matrix, only : matr3inv
56 : use m_io_tools, only : open_file
57 : use m_occ, only : newocc, getnel
58 : use m_ddb_hdr, only : ddb_hdr_type
59 : use m_fstrings, only : strcat, sjoin
60 : use m_geometry, only : fixsym, mkradim, metric
61 : use m_kpts, only : tetra_from_kptrlatt
62 : use m_kg, only : kpgio, getph
63 : use m_fft, only : fourdp
64 : use m_pawang, only : pawang_type
65 : use m_pawrad, only : pawrad_type
66 : use m_pawtab, only : pawtab_type, pawtab_print
67 : use m_pawcprj, only : pawcprj_type,pawcprj_free,pawcprj_alloc, pawcprj_getdim
68 : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
69 : use m_abi2big, only : wvl_occ_abi2big, wvl_setngfft, wvl_setBoxGeometry
70 : use m_energies, only : energies_type
71 : use m_args_gs, only : args_gs_type
72 : use m_results_gs, only : results_gs_type
73 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_free
74 : use m_paw_dmft, only : init_sc_dmft,destroy_sc_dmft,print_sc_dmft,paw_dmft_type,readocc_dmft
75 : use m_paw_sphharm, only : setsym_ylm
76 : use m_paw_occupancies, only : initrhoij
77 : use m_paw_init, only : pawinit,paw_gencond
78 : use m_paw_correlations, only : pawpuxinit
79 : use m_paw_uj, only : pawuj_ini,pawuj_free,pawuj_det, macro_uj_type
80 : use m_data4entropyDMFT, only : data4entropyDMFT_t, data4entropyDMFT_init, data4entropyDMFT_destroy
81 : use m_electronpositron, only : electronpositron_type,init_electronpositron,destroy_electronpositron, &
82 : electronpositron_calctype
83 : use m_scfcv, only : scfcv_t, scfcv_init, scfcv_destroy, scfcv_run
84 : use m_jellium, only : jellium
85 : use m_iowf, only : outwf, outresid
86 : use m_outqmc, only : outqmc
87 : use m_ioarr, only : ioarr,read_rhor
88 : use m_inwffil, only : inwffil
89 : use m_spacepar, only : setsym
90 : use m_mkrho, only : mkrho, initro, prtrhomxmn
91 : use m_initylmg, only : initylmg
92 : use m_pspini, only : pspini
93 : use m_mover, only : mover
94 : use m_mpinfo, only : proc_distrb_cycle
95 : use m_common, only : setup1, prteigrs, prtene
96 : use m_fourier_interpol, only : transgrid
97 : use m_psolver, only : psolver_kernel
98 : use m_paw2wvl, only : paw2wvl, wvl_paw_free
99 : use m_berryphase_new, only : init_e_field_vars,prtefield
100 : use m_wvl_wfs, only : wvl_wfs_set, wvl_wfs_free, wvl_wfs_lr_copy
101 : use m_wvl_rho, only : wvl_initro, wvl_mkrho
102 : use m_wvl_descr_psp, only : wvl_descr_psp_set, wvl_descr_free, wvl_descr_atoms_set, wvl_descr_atoms_set_sym
103 : use m_wvl_denspot, only : wvl_denspot_set, wvl_denspot_free
104 : use m_wvl_projectors, only : wvl_projectors_set, wvl_projectors_free
105 : use m_cgprj, only : ctocprj
106 : use m_nonlop_ylm, only : nonlop_ylm_init_counters,nonlop_ylm_output_counters
107 : use m_fft, only : fft_init_counters,fft_output_counters
108 : use m_pstat, only : pstat_proc
109 :
110 : #if defined(HAVE_GPU_MARKERS)
111 : use m_nvtx_data
112 : #endif
113 :
114 : #if defined HAVE_YAKL
115 : use gator_mod
116 : #endif
117 :
118 : use defs_wvltypes, only : wvl_data,coulomb_operator,wvl_wf_type
119 : #if defined HAVE_BIGDFT
120 : use BigDFT_API, only : wvl_timing => timing,xc_init,xc_end,XC_MIXED,XC_ABINIT,&
121 : local_potential_dimensions,nullify_gaussian_basis, &
122 : copy_coulomb_operator,deallocate_coulomb_operator
123 : #else
124 : use defs_wvltypes, only : coulomb_operator
125 : #endif
126 :
127 : #if defined HAVE_LOTF
128 : use defs_param_lotf, only : lotfparam_init
129 : #endif
130 :
131 : implicit none
132 :
133 : private
134 : !!***
135 :
136 : public :: gstate
137 : !!***
138 :
139 : contains
140 : !!***
141 :
142 : !!****f* m_gstate/gstate
143 : !! NAME
144 : !! gstate
145 : !!
146 : !! FUNCTION
147 : !! Primary routine for conducting DFT calculations by CG minimization.
148 : !!
149 : !! INPUTS
150 : !! args_gs<type(args_gs_type)>=various input arguments for the GS calculation
151 : !! Possibly different from dtset
152 : !! codvsn=code version
153 : !! cpui=initial CPU time
154 : !! itimimage_gstate=counter for calling do loop
155 : !!
156 : !! OUTPUT
157 : !! npwtot(nkpt) = total number of plane waves at each k point
158 : !! results_gs <type(results_gs_type)>=results (energy and its components,
159 : !! forces and its components, the stress tensor) of a ground-state computation
160 : !!
161 : !! SIDE EFFECTS
162 : !! acell(3)=unit cell length scales (bohr)
163 : !! dtfil <type(datafiles_type)>=variables related to files
164 : !! dtset <type(dataset_type)>=all input variables in this dataset
165 : !! | mband =maximum number of bands (IN)
166 : !! | mgfft =maximum single fft dimension (IN)
167 : !! | mkmem =number of k points treated by this processor (IN)
168 : !! | mpw =maximum number of planewaves in basis sphere (large number) (IN)
169 : !! | natom =number of atoms in unit cell (IN)
170 : !! | nfft =(effective) number of FFT grid points (for this processor) (IN)
171 : !! | nkpt =number of k points (IN)
172 : !! | nspden=number of spin-density components (IN)
173 : !! | nsppol=number of channels for spin-polarization (1 or 2) (IN)
174 : !! | nsym =number of symmetry elements in space group
175 : !! iexit= exit flag
176 : !! initialized= 0 for the first GS calculation (not initialized), else 1
177 : !! mpi_enreg=MPI-parallelisation information (some already initialized,
178 : !! some others to be initialized here)
179 : !! occ(mband*nkpt*nsppol) = occupation number for each band and k
180 : !! pawang <type(pawang_type)>=paw angular mesh and related data
181 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
182 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
183 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
184 : !! Before entering the first time in gstate, a significant part of
185 : !! psps has been initialized :
186 : !! the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,
187 : !! ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp
188 : !! All the remaining components of psps are to be initialized in the call
189 : !! to pspini .
190 : !! The next time the code enters gstate, psps might be identical to the
191 : !! one of the previous dtset, in which case, no reinitialisation is scheduled
192 : !! in pspini.f .
193 : !! rprim(3,3)=dimensionless real space primitive translations
194 : !! scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles
195 : !! vel(3,natom)=value of velocity
196 : !! vel_cell(3,3)=value of cell parameters velocity
197 : !! wvl <type(wvl_data)>=all wavelets data
198 : !! xred(3,natom) = reduced atomic coordinates
199 : !!
200 : !! NOTES
201 : !! USE OF FFT GRIDS:
202 : !! =================
203 : !! In case of PAW:
204 : !! ---------------
205 : !! Two FFT grids are used:
206 : !! - A "coarse" FFT grid (defined by ecut)
207 : !! for the application of the Hamiltonian on the plane waves basis.
208 : !! It is defined by nfft, ngfft, mgfft, ...
209 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ...
210 : !! are expressed on this grid.
211 : !! - A "fine" FFT grid (defined) by ecutdg)
212 : !! for the computation of the density inside PAW spheres.
213 : !! It is defined by nfftf, ngfftf, mgfftf, ...
214 : !! Total density, potentials, ...
215 : !! are expressed on this grid.
216 : !! In case of norm-conserving:
217 : !! ---------------------------
218 : !! - Only the usual FFT grid (defined by ecut) is used.
219 : !! It is defined by nfft, ngfft, mgfft, ...
220 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf)
221 : !! are set equal to (nfft,ngfft,mgfft) in that case.
222 : !! In case of wavelets:
223 : !! --------------------
224 : !! - Only the usual FFT grid (defined by wvl_crmult) is used.
225 : !! It is defined by nfft, ngfft, mgfft, ... This is strictly not
226 : !! an FFT grid since its dimensions are not suited for FFTs. They are
227 : !! defined by wvl_setngfft().
228 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf)
229 : !! are set equal to (nfft,ngfft,mgfft) in that case.
230 : !!
231 : !! TODO
232 : !! Not yet possible to use restartxf in parallel when localrdwf==0
233 : !!
234 : !! SOURCE
235 :
236 5629 : subroutine gstate(args_gs,acell,codvsn,cpui,dtfil,dtset,iexit,initialized,&
237 5629 : & itimimage_gstate,mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,&
238 5629 : & psps,results_gs,rprim,scf_history,vel,vel_cell,wvl,xred)
239 :
240 : !Arguments ------------------------------------
241 : !scalars
242 : integer,intent(inout) :: iexit,initialized
243 : integer,intent(in) :: itimimage_gstate
244 : real(dp),intent(in) :: cpui
245 : character(len=8),intent(in) :: codvsn
246 : type(MPI_type),intent(inout) :: mpi_enreg
247 : type(args_gs_type),intent(in) :: args_gs
248 : type(datafiles_type),intent(inout) :: dtfil
249 : type(dataset_type),intent(inout) :: dtset
250 : type(pawang_type),intent(inout) :: pawang
251 : type(pseudopotential_type),intent(inout) :: psps
252 : type(results_gs_type),intent(inout) :: results_gs
253 : type(scf_history_type),target,intent(inout) :: scf_history
254 : type(wvl_data),intent(inout) :: wvl
255 : !arrays
256 : integer,intent(out) :: npwtot(dtset%nkpt)
257 : real(dp),intent(inout) :: acell(3),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
258 : real(dp),intent(inout) :: rprim(3,3),vel(3,dtset%natom),vel_cell(3,3),xred(3,dtset%natom)
259 : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
260 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
261 :
262 : !Local variables-------------------------------
263 : !Define file format for different type of files. Presently,
264 : !only one file format is supported for each type of files, but this might change soon ...
265 : !2 for wavefunction file, new format (version 2.0 and after) (fform) NOT USED
266 : !52 for density rho(r) (fformr)
267 : !102 for potential V(r) file. (fformv) NOT USED
268 : !scalars
269 : logical :: compute_cprj
270 : integer,parameter :: formeig=0, level=101, response=0 ,cplex1=1, master=0, itime0=0
271 : integer :: ndtpawuj=0 ! Cannot use parameter because scfargs points to this! Have to get rid of pointers to scalars!
272 : #if defined HAVE_BIGDFT
273 : integer :: icoulomb
274 : #endif
275 : integer :: accessfil,ask_accurate,bantot,choice,comm_psp,fform
276 : integer :: gnt_option,gscase,iatom,idir,ierr,ii,indx,jj,kk,ios,iorder_cprj,itypat
277 : integer :: ixfh,mband_cprj,mcg,mcprj,me,mgfftf,mpert,mu,my_natom,my_nspinor
278 : integer :: nband_k,nbandtot,nblok,ncprj,ncpgr,nfftf,nfftot,npwmin
279 : integer :: openexit,option,optorth,psp_gencond,conv_retcode
280 : integer :: pwind_alloc,rdwrpaw,comm,tim_mkrho,use_sc_dmft
281 : integer :: cnt,spin,band,ikpt,usecg,usecprj,ylm_option
282 : real(dp) :: cpus,ecore,ecut_eff,ecutdg_eff,etot,fermie,fermih
283 : real(dp) :: gsqcut_eff,gsqcut_shp,gsqcutc_eff,hyb_range_fock,residm,ucvol
284 : logical :: read_wf_or_den,has_to_init,call_pawinit,write_wfk,inv_sij
285 : logical :: is_dfpt=.false.,wvlbigdft=.false.
286 : character(len=500) :: msg
287 : character(len=fnlen) :: dscrpt,filnam,wfkfull_path
288 : real(dp) :: fatvshift
289 292708 : type(crystal_t) :: cryst
290 11258 : type(ebands_t) :: bstruct, ebands, ebands_bz
291 5629 : type(efield_type) :: dtefield
292 : type(electronpositron_type),pointer :: electronpositron
293 5629 : type(hdr_type) :: hdr, hdr_den, hdr_bz
294 : type(extfpmd_type),pointer :: extfpmd => null()
295 : type(rcpaw_type), pointer :: rcpaw => null()
296 61919 : type(macro_uj_type) :: dtpawuj(1)
297 5629 : type(paw_dmft_type) :: paw_dmft
298 : type(pawfgr_type) :: pawfgr
299 5629 : type(recursion_type) ::rec_set
300 : type(wffile_type) :: wff1,wffnew,wffnow
301 5629 : type(ab_xfh_type) :: ab_xfh
302 5629 : type(ddb_type) :: ddb
303 292708 : type(ddb_hdr_type) :: ddb_hdr
304 5629 : type(scfcv_t) :: scfcv_args
305 : !arrays
306 : integer :: itimes(2),ngfft(18),ngfftf(18)
307 16887 : integer,allocatable :: atindx(:),atindx1(:),indsym(:,:,:),dimcprj_srt(:)
308 5629 : integer,allocatable :: irrzon(:,:,:),kg(:,:),nattyp(:),symrec(:,:,:)
309 5629 : integer,allocatable,target :: npwarr(:)
310 5629 : integer,pointer :: npwarr_(:),pwind(:,:,:)
311 : real(dp) :: efield_band(3),gmet(3,3),gmet_for_kg(3,3),gprimd(3,3),gprimd_for_kg(3,3)
312 : real(dp) :: rmet(3,3),rprimd(3,3),rprimd_for_kg(3,3),tsec(2)
313 5629 : real(dp),allocatable :: doccde(:)
314 5629 : real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phnons(:,:,:),resid(:),rhowfg(:,:)
315 5629 : real(dp),allocatable :: rhowfr(:,:),spinat_dum(:,:),start(:,:),work(:)
316 5629 : real(dp),allocatable :: ylm(:,:),ylmgr(:,:,:)
317 : real(dp),contiguous, pointer :: cg(:,:) => null()
318 5629 : real(dp),pointer :: eigen(:),pwnsfac(:,:),rhog(:,:),rhor(:,:)
319 5629 : real(dp),pointer :: taug(:,:),taur(:,:),xred_old(:,:)
320 5629 : type(pawrhoij_type),pointer :: pawrhoij(:)
321 : type(coulomb_operator) :: kernel_dummy
322 5629 : type(pawcprj_type),allocatable :: cprj(:,:)
323 5629 : type(xg_nonlop_t) :: xg_nonlop
324 : ! ***********************************************************************
325 :
326 : DBG_ENTER("COLL")
327 :
328 5629 : call timab(1232,1,tsec)
329 5629 : call timab(1211,3,tsec)
330 :
331 : !###########################################################
332 : !### 01. Initializations XML, MPI, WVL, etc
333 :
334 : !Init MPI data
335 5629 : comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm)
336 :
337 : !Set up MPI information from the dataset
338 5629 : my_natom=mpi_enreg%my_natom
339 :
340 : !Set up information when wavelets are in use
341 5629 : if (dtset%usewvl == 1) then
342 :
343 : ! If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
344 0 : wvlbigdft=(dtset%wvl_bigdft_comp==1)
345 :
346 : ! Default value, to be set-up elsewhere.
347 0 : wvl%descr%h(:) = dtset%wvl_hgrid
348 :
349 : #if defined HAVE_BIGDFT
350 : wvl%descr%paw%usepaw=psps%usepaw
351 : wvl%descr%paw%natom=dtset%natom
352 : #endif
353 :
354 : ! We set the atom-related internal wvl variables.
355 : call wvl_descr_atoms_set(acell, dtset%icoulomb, dtset%natom, &
356 0 : & dtset%ntypat, dtset%typat, wvl%descr)
357 : if(dtset%usepaw==0) then
358 : ! nullify PAW proj_G in NC case:
359 : #if defined HAVE_BIGDFT
360 : ABI_MALLOC(wvl%projectors%G,(dtset%ntypat))
361 : do itypat=1,dtset%ntypat
362 : call nullify_gaussian_basis(wvl%projectors%G(itypat))
363 : end do
364 : #endif
365 : end if
366 :
367 0 : wvl%descr%exctxpar = "OP2P"
368 : end if
369 :
370 5629 : if (me == master .and. dtset%prtxml == 1) then
371 : ! gstate() will handle a dataset, so we output the dataSet markup.
372 1 : write(ab_xml_out, "(A)") ' <dataSet>'
373 : end if
374 :
375 : !Define FFT grid(s) sizes (be careful !)
376 : !See NOTES in the comments at the beginning of this file.
377 5629 : call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
378 :
379 : !Structured debugging if prtvol==-level
380 5629 : if(dtset%prtvol==-level)then
381 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' gstate : enter , debug mode '
382 0 : call wrtout(std_out, msg)
383 : end if
384 :
385 : !###########################################################
386 : !### 02. Calls setup1, kpgio, initylmg
387 :
388 5629 : ecore=zero
389 22516 : results_gs%pel(1:3) =zero
390 58573 : results_gs%grchempottn(:,:)=zero
391 58573 : results_gs%grewtn(:,:)=zero
392 : !MT Feb 2012: I dont know why but grvdw has to be allocated
393 : !when using BigDFT to ensure success on inca_gcc44_sdebug
394 5629 : if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).or.dtset%usewvl==1) then
395 5 : results_gs%ngrvdw=dtset%natom
396 5 : ABI_SFREE(results_gs%grvdw)
397 15 : ABI_MALLOC(results_gs%grvdw,(3,dtset%natom))
398 5653 : results_gs%grvdw(:,:)=zero
399 : end if
400 5629 : call results_gs%energies%init()
401 :
402 : !Set up for iterations
403 : call setup1(acell,bantot,dtset,&
404 : ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,&
405 : ngfftf,ngfft,dtset%nkpt,dtset%nsppol,&
406 5629 : response,rmet,rprim,rprimd,ucvol,psps%usepaw)
407 :
408 : !In some cases (e.g. getcell/=0), the plane wave vectors have
409 : ! to be generated from the original simulation cell
410 5629 : rprimd_for_kg=rprimd
411 5797 : if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=args_gs%rprimd_orig
412 5629 : if (dtset%optcell/=0.and.dtset%imgmov/=0) rprimd_for_kg=args_gs%rprimd_orig
413 5629 : call matr3inv(rprimd_for_kg,gprimd_for_kg)
414 225160 : gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg)
415 :
416 : !Set up the basis sphere of planewaves
417 16887 : ABI_MALLOC(npwarr,(dtset%nkpt))
418 5629 : if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2) then
419 16869 : ABI_MALLOC(kg,(3,dtset%mpw*dtset%mkmem))
420 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet_for_kg,dtset%istwfk,kg, &
421 : & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,&
422 5623 : & dtset%mpw,npwarr,npwtot,dtset%nsppol)
423 : call bandfft_kpt_init1(bandfft_kpt,dtset%istwfk,kg,dtset%mgfft,dtset%mkmem,mpi_enreg,&
424 5623 : & dtset%mpw,dtset%nband,dtset%nkpt,npwarr,dtset%nsppol,gpu_option=dtset%gpu_option)
425 : else
426 6 : ABI_MALLOC(kg,(0,0))
427 12 : npwarr(:) = 0
428 12 : npwtot(:) = 0
429 : end if
430 :
431 : ! Handling GEMM nonlop use
432 : ! Not enabled by default for CPU and CUDA implementations
433 : ! Enabled if using OpenMP GPU offload (only implementation)
434 5629 : gemm_nonlop_use_gemm = .false.
435 :
436 5629 : gemm_nonlop_is_distributed = .false.
437 5629 : if(dtset%gpu_nl_distrib == 1) gemm_nonlop_is_distributed = .true.
438 5629 : if(dtset%gpu_nl_splitsize > 0) gemm_nonlop_block_size = dtset%gpu_nl_splitsize
439 :
440 5629 : if(dtset%gpu_option == ABI_GPU_OPENMP .or. dtset%use_gemm_nonlop == 1) then
441 29 : gemm_nonlop_use_gemm = .true.
442 29 : call init_gemm_nonlop(dtset%gpu_option)
443 : end if
444 :
445 : ! Handle GPU FFT slicing
446 5629 : hamilt_gpu_nfft_blocks = dtset%gpu_nfft_blocks
447 :
448 : !Set up the Ylm for each k point
449 5629 : if ( dtset%tfkinfunc /= 2) then
450 22492 : ABI_MALLOC(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm))
451 22492 : ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm))
452 5623 : if (psps%useylm==1) then
453 1338 : ylm_option=0
454 1338 : if (dtset%prtstm==0.and.dtset%iscf>0.and.dtset%positron/=1) ylm_option=1 ! compute gradients of YLM
455 1338 : if (dtset%berryopt==4 .and. dtset%optstress /= 0 .and. psps%usepaw==1) ylm_option = 1 ! compute gradients of YLM
456 : call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,&
457 : & psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,&
458 1338 : & npwarr,dtset%nsppol,ylm_option,rprimd,ylm,ylmgr)
459 : end if
460 : else
461 6 : ABI_MALLOC(ylm,(0,0))
462 6 : ABI_MALLOC(ylmgr,(0,0,0))
463 : end if
464 :
465 : !SCF history management (allocate it at first call)
466 5629 : if (initialized==0) then
467 : ! This call has to be done before any use of SCF history
468 4231 : usecg=0
469 4231 : if(dtset%extrapwf>0 .or. dtset%imgwfstor==1)usecg=1
470 4231 : call scf_history_init(dtset,mpi_enreg,usecg,scf_history)
471 : end if
472 5629 : has_to_init=(initialized==0.or.scf_history%history_size<0)
473 :
474 5629 : call timab(1211,2,tsec)
475 5629 : call timab(1212,3,tsec)
476 :
477 : !###########################################################
478 : !### 03. Calls pspini
479 :
480 : !Open and read pseudopotential files
481 5629 : comm_psp=mpi_enreg%comm_cell;if (dtset%usewvl==1) comm_psp=mpi_enreg%comm_wvl
482 5653 : if (dtset%nimage>1) psps%mixalch(:,:)=args_gs%mixalch(:,:) ! mixalch can evolve for some image algos
483 5629 : call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,psps,rprimd,comm_mpi=comm_psp)
484 5629 : call timab(1212,2,tsec)
485 5629 : call timab(1211,3,tsec)
486 :
487 : !In case of isolated computations, ecore must set to zero
488 : !because its contribution is counted in the ewald energy as the ion-ion interaction.
489 5629 : if (dtset%icoulomb == 1) ecore = zero
490 :
491 : !WVL - Now that psp data are available, we compute rprimd, acell... from the atomic positions.
492 5629 : if (dtset%usewvl == 1) then
493 0 : call wvl_descr_psp_set(trim(dtfil%filnam_ds(3))//"_OCCUP",dtset%nsppol,psps,dtset%spinat,wvl%descr)
494 : call wvl_setBoxGeometry(dtset%prtvol, psps%gth_params%radii_cf, rprimd, xred, &
495 0 : & wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult)
496 0 : call mkradim(acell,rprim,rprimd)
497 0 : rprimd_for_kg=rprimd
498 : call wvl_denspot_set(wvl%den, psps%gth_params, dtset%ixc, dtset%natom, dtset%nsppol, rprimd, &
499 0 : & wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult, mpi_enreg%comm_wvl, xred)
500 : ! TODO: to be moved in a routine.
501 : #if defined HAVE_BIGDFT
502 : if (wvl%descr%atoms%astruct%geocode == "F") then
503 : icoulomb = 1
504 : else if (wvl%descr%atoms%astruct%geocode == "S") then
505 : icoulomb = 2
506 : else
507 : icoulomb = 0
508 : end if
509 : ! calculation of the Poisson kernel anticipated to reduce memory peak for small systems
510 : call psolver_kernel( wvl%den%denspot%dpbox%hgrids, 1, icoulomb, mpi_enreg%me_wvl, wvl%den%denspot%pkernel , &
511 : & mpi_enreg%comm_wvl, wvl%den%denspot%dpbox%ndims, mpi_enreg%nproc_wvl, dtset%nscforder)
512 : nullify(wvl%den%denspot%pkernelseq%kernel)
513 : !call copy_coulomb_operator(wvl%den%denspot%pkernel,wvl%den%denspot%pkernelseq, "gstate")
514 : ! Associate the denspot distribution into mpi_enreg.
515 : mpi_enreg%nscatterarr => wvl%den%denspot%dpbox%nscatterarr
516 : mpi_enreg%ngatherarr => wvl%den%denspot%dpbox%ngatherarr
517 : mpi_enreg%ngfft3_ionic = wvl%den%denspot%dpbox%n3pi
518 : call wvl_setngfft(mpi_enreg%me_wvl, dtset%mgfft, dtset%nfft, &
519 : & dtset%ngfft, mpi_enreg%nproc_wvl, wvl%den%denspot%dpbox%ndims(1), &
520 : & wvl%den%denspot%dpbox%ndims(2),wvl%den%denspot%dpbox%ndims(3),&
521 : & wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl, 1))
522 : #endif
523 0 : nfftf = dtset%nfft
524 0 : mgfftf = dtset%mgfft
525 0 : ngfftf(:) = dtset%ngfft(:)
526 : ! Recalculate gprimd
527 0 : call matr3inv(rprimd,gprimd)
528 : ! PAW section
529 0 : if(psps%usepaw==1) then
530 : ! Reinitialize Pawfgr with new values of ngfft
531 0 : call pawfgr_destroy(pawfgr)
532 0 : call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
533 : ! fill wvl objects from paw objects
534 : ! wvl%descr%paw%usepaw=dtset%usepaw
535 0 : call paw2wvl(pawtab,wvl%projectors,wvl%descr)
536 : end if
537 : else if (dtset%icoulomb /= 0) then
538 : #if defined HAVE_BIGDFT
539 : if (dtset%ixc < 0) then
540 : call xc_init(wvl%den%denspot%xc, dtset%ixc, XC_MIXED, dtset%nsppol)
541 : else
542 : call xc_init(wvl%den%denspot%xc, dtset%ixc, XC_ABINIT, dtset%nsppol)
543 : end if
544 : #endif
545 : end if
546 :
547 : !Initialize band structure datatype
548 5629 : if (dtset%paral_kgb/=0) then ! We decide to store total npw in bstruct,
549 1809 : ABI_MALLOC(npwarr_,(dtset%nkpt))
550 7449 : npwarr_(:)=npwarr(:)
551 603 : call xmpi_sum(npwarr_,mpi_enreg%comm_bandfft,ierr)
552 : else
553 5026 : npwarr_ => npwarr
554 : end if
555 :
556 5629 : call bstruct%from_dtset(dtset, npwarr_)
557 :
558 5629 : if (dtset%paral_kgb/=0) then
559 603 : ABI_FREE(npwarr_)
560 : end if
561 5629 : nullify(npwarr_)
562 :
563 : !Initialize PAW atomic occupancies
564 5629 : if (scf_history%history_size>=0) then
565 1454 : pawrhoij => scf_history%pawrhoij_last
566 : else
567 14291 : ABI_MALLOC(pawrhoij,(my_natom*psps%usepaw))
568 : end if
569 5629 : if (psps%usepaw==1.and.has_to_init) then
570 : call initrhoij(dtset%pawcpxocc,dtset%lexexch,&
571 : & dtset%lpawu,my_natom,dtset%natom,dtset%nspden,dtset%nspinor,dtset%nsppol,&
572 : & dtset%ntypat,pawrhoij,dtset%pawspnorb,pawtab,cplex1,dtset%spinat,dtset%typat,&
573 1070 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
574 : end if
575 :
576 : !Initialize header
577 5629 : gscase=0
578 : call hdr%init(bstruct,codvsn,dtset,pawtab,gscase,psps,wvl%descr,&
579 5629 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
580 :
581 : !Clean band structure datatype (should use it more in the future !)
582 5629 : call bstruct%free()
583 :
584 : !Update header, with evolving variables, when available
585 : !Here, rprimd, xred and occ are available
586 5629 : etot=hdr%etot ; fermie=hdr%fermie ; fermih=hdr%fermih ; residm=hdr%residm
587 : call hdr%update(bantot,etot,fermie,fermih,&
588 : residm,rprimd,occ,pawrhoij,xred,args_gs%amu,&
589 5629 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
590 :
591 : ! PW basis set: test if the problem is ill-defined.
592 5629 : if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2) then
593 75328 : npwmin=minval(hdr%npwarr(:))
594 5623 : if (dtset%mband > npwmin) then
595 : ! No way we can solve the problem. Abort now!
596 : write(msg,"(2(a,i0),4a)")&
597 0 : "Number of bands nband= ",dtset%mband," > number of planewaves npw= ",npwmin,ch10,&
598 0 : "The number of eigenvectors cannot be greater that the size of the Hamiltonian!",ch10,&
599 0 : "Action: decrease nband or, alternatively, increase ecut"
600 0 : if (dtset%ionmov/=23) then
601 0 : ABI_ERROR(msg)
602 : else
603 0 : ABI_WARNING(msg)
604 : end if
605 :
606 5623 : else if (dtset%mband >= 0.9 * npwmin) then
607 : ! Warn the user
608 : write(msg,"(a,i0,a,f6.1,4a)")&
609 3 : & "Number of bands nband= ",dtset%mband," >= 0.9 * maximum number of planewaves= ",0.9*npwmin,ch10,&
610 3 : & "The problem is ill-defined and the GS algorithm will show numerical instabilities!",ch10,&
611 6 : & "Assume experienced user. Execution will continue."
612 3 : ABI_WARNING(msg)
613 : end if
614 : end if
615 :
616 : !###########################################################
617 : !### 04. Symmetry operations when nsym>1
618 :
619 : !Do symmetry stuff only for nsym>1
620 5629 : if (dtset%usewvl == 0) then
621 5629 : nfftot=ngfft(1)*ngfft(2)*ngfft(3)
622 : else
623 : #if defined HAVE_BIGDFT
624 : nfftot=product(wvl%den%denspot%dpbox%ndims)
625 : #else
626 0 : BIGDFT_NOTENABLED_ERROR()
627 : #endif
628 : end if
629 22516 : ABI_MALLOC(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
630 22516 : ABI_MALLOC(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
631 22516 : ABI_MALLOC(indsym,(4,dtset%nsym,dtset%natom))
632 16887 : ABI_MALLOC(symrec,(3,3,dtset%nsym))
633 133484430 : irrzon(:,:,:)=0
634 200204024 : phnons(:,:,:)=zero
635 1168255 : indsym(:,:,:)=0
636 1373671 : symrec(:,:,:)=0
637 :
638 5629 : if (dtset%nsym>1) then
639 : call setsym(indsym,irrzon,dtset%iscf,dtset%natom,&
640 : & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
641 4486 : & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred)
642 :
643 : ! Make sure dtset%iatfix does not break symmetry
644 4486 : call fixsym(dtset%iatfix,indsym,dtset%natom,dtset%nsym)
645 : else
646 : ! The symrec array is used by initberry even in case nsym = 1
647 14859 : symrec(:,:,1) = 0
648 1143 : symrec(1,1,1) = 1 ; symrec(2,2,1) = 1 ; symrec(3,3,1) = 1
649 : end if
650 5629 : if (dtset%usewvl == 1) then
651 : call wvl_descr_atoms_set_sym(wvl%descr, dtset%efield, irrzon, dtset%nsppol, &
652 0 : & dtset%nsym, phnons, dtset%symafm, dtset%symrel, dtset%tnons, dtset%tolsym)
653 : #if defined HAVE_BIGDFT
654 : wvl%den%symObj = wvl%descr%atoms%astruct%sym%symObj
655 : #endif
656 : end if
657 :
658 : !###########################################################
659 : !### 05. Calls inwffil
660 : ABI_NVTX_START_RANGE(NVTX_INIT_INWFFIL)
661 :
662 : ! if paral_kgb == 0, it may happen that some processors are idle (no entry in proc_distrb)
663 : ! but mkmem == nkpt and this can cause integer overflow in mcg or allocation error.
664 : ! Here we count the number of states treated by the proc. if cnt == 0, mcg is then set to 0.
665 5629 : cnt = 0
666 5629 : nbandtot = 0
667 11799 : do spin=1,dtset%nsppol
668 83582 : do ikpt=1,dtset%nkpt
669 71783 : nband_k = dtset%nband(ikpt + (spin-1) * dtset%nkpt)
670 516066 : do band=1,nband_k
671 516066 : if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, band, band, spin, mpi_enreg%me_kpt)) cnt = cnt + 1
672 : end do
673 77953 : nbandtot = nbandtot + nband_k
674 : end do
675 : end do
676 :
677 5629 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
678 5629 : mcg=dtset%mpw*my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
679 5629 : if (cnt == 0) then
680 42 : mcg = 0
681 42 : write(msg,"(2(a,i0))")"rank: ",mpi_enreg%me, "does not have wavefunctions to treat. Setting mcg to: ",mcg
682 42 : ABI_WARNING(msg)
683 : end if
684 :
685 5629 : if (dtset%usewvl == 0 .and. dtset%mpw > 0 .and. cnt /= 0)then
686 5587 : if (my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol > floor(real(HUGE(0))/real(dtset%mpw) )) then
687 : write (msg,'(9a)')&
688 0 : "Default integer is not wide enough to store the size of the wavefunction array (mcg).",ch10,&
689 0 : "This usually happens when paral_kgb == 0 and there are not enough procs to distribute kpts and spins",ch10,&
690 0 : "Action: if paral_kgb == 0, use nprocs = nkpt * nsppol to reduce the memory per node.",ch10,&
691 0 : "If this does not solve the problem, use paral_kgb 1 with nprocs > nkpt * nsppol and use npfft/npband/npspinor",ch10,&
692 0 : "to decrease the memory requirements. Consider also OpenMP threads."
693 : ii = 0
694 0 : ABI_ERROR_NOSTOP(msg,ii)
695 : write (msg,'(5(a,i0), 2a)')&
696 0 : "my_nspinor: ",my_nspinor, ", mpw: ",dtset%mpw, ", mband: ",dtset%mband,&
697 0 : ", mkmem: ",dtset%mkmem, ", nsppol: ",dtset%nsppol,ch10,&
698 0 : 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS...) compiled in int64 mode'
699 0 : ABI_ERROR(msg)
700 : end if
701 : end if
702 :
703 5629 : if (dtset%imgwfstor==1) then
704 53 : cg => scf_history%cg(:,:,1)
705 53 : eigen => scf_history%eigen(:,1)
706 : else
707 5576 : if(dtset%gpu_option == ABI_GPU_KOKKOS) then
708 : #if defined HAVE_GPU && defined HAVE_YAKL
709 : ABI_MALLOC_MANAGED(cg, (/2,mcg/))
710 : #endif
711 : else
712 16728 : ABI_MALLOC_OR_DIE(cg,(2,mcg), ierr)
713 : end if
714 16728 : ABI_MALLOC(eigen,(dtset%mband*dtset%nkpt*dtset%nsppol))
715 : end if
716 :
717 16887 : ABI_MALLOC(resid,(dtset%mband*dtset%nkpt*dtset%nsppol))
718 894347 : eigen(:)=zero ; resid(:)=zero
719 : !mpi_enreg%paralbd=0 ; ask_accurate=0
720 5629 : ask_accurate=0
721 :
722 : !WVL - Branching, allocating wavefunctions as wavelets.
723 5629 : if (dtset%usewvl == 1) then
724 0 : call wvl_wfs_lr_copy(wvl%wfs, wvl%descr)
725 : ! Create access arrays for wavefunctions and allocate wvl%wfs%psi (other arrays are left unallocated).
726 : call wvl_wfs_set(dtset%strprecon,dtset%spinmagntarget, dtset%kpt, mpi_enreg%me_wvl,&
727 : & dtset%natom, sum(dtset%nband), &
728 : & dtset%nkpt, mpi_enreg%nproc_wvl, dtset%nspinor, dtset%nsppol, dtset%nwfshist, occ, &
729 : & psps, rprimd, wvl%wfs, dtset%wtk, wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult, &
730 0 : & xred)
731 : ! We transfer wavelets information to the hdr structure.
732 : #if defined HAVE_BIGDFT
733 : call local_potential_dimensions(mpi_enreg%me_wvl,wvl%wfs%ks%lzd,wvl%wfs%ks%orbs,wvl%den%denspot%xc,&
734 : & wvl%den%denspot%dpbox%ngatherarr(0,1))
735 : hdr%nwvlarr(1) = wvl%wfs%ks%lzd%Glr%wfd%nvctr_c
736 : hdr%nwvlarr(2) = 7 * wvl%wfs%ks%lzd%Glr%wfd%nvctr_f
737 : #endif
738 : ! Create access arrays for projectors and allocate them.
739 : ! Compute projectors from each atom.
740 : call wvl_projectors_set(mpi_enreg%me_wvl, dtset%natom, wvl%projectors, psps, rprimd, &
741 0 : & wvl%wfs, wvl%descr, dtset%wvl_frmult, xred)
742 : end if
743 :
744 5629 : read_wf_or_den=(dtset%iscf<=0.or.dtfil%ireadwf/=0.or.(dtfil%ireadden/=0.and.dtset%positron<=0))
745 5629 : read_wf_or_den=(read_wf_or_den.and.has_to_init)
746 :
747 : !RECURSION - initialization
748 5629 : if(has_to_init .and. dtset%userec==1) then
749 18 : call InitRec(dtset,mpi_enreg,rec_set,rmet,maxval(psps%indlmn(3,:,:)))
750 : end if
751 :
752 : !LOTF - initialization
753 : #if defined HAVE_LOTF
754 : if(has_to_init .and. dtset%ionmov==23) then
755 : call lotfparam_init(dtset%natom,dtset%lotf_version,1,&
756 : & dtset%lotf_nitex,dtset%lotf_nneigx,&
757 : & dtset%lotf_classic,1,1)
758 : end if
759 : #endif
760 :
761 : !Initialize wavefunctions.
762 5629 : if(dtset%imgwfstor==1 .and. initialized==1)then
763 1012074 : cg(:,:)=scf_history%cg(:,:,1)
764 714 : eigen(:)=scf_history%eigen(:,1)
765 5587 : else if(dtset%tfkinfunc /=2) then
766 : !if(dtset%tfkinfunc /=2) then
767 5581 : wff1%unwff=dtfil%unwff1
768 5581 : optorth=1 !if (psps%usepaw==1) optorth=0
769 5581 : if(psps%usepaw==1 .and. dtfil%ireadwf==1)optorth=0
770 72553 : hdr%rprimd=rprimd_for_kg ! We need the rprimd that was used to generate de G vectors
771 : ABI_NVTX_START_RANGE(NVTX_INIT_INWFFIL2)
772 : call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen,&
773 : & dtset%exchn2n3d,formeig,hdr,dtfil%ireadwf,dtset%istwfk,kg,&
774 : & dtset%kptns,dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,&
775 : & dtset%mpw,dtset%nband,ngfft,dtset%nkpt,npwarr,&
776 : & dtset%nsppol,dtset%nsym,occ,optorth,dtset%symafm,&
777 : & dtset%symrel,dtset%tnons,dtfil%unkg,wff1,wffnow,dtfil%unwff1,&
778 5581 : & dtfil%fnamewffk,wvl)
779 : ABI_NVTX_END_RANGE()
780 72553 : hdr%rprimd=rprimd
781 : end if
782 :
783 5629 : if (psps%usepaw==1.and.dtfil%ireadwf==1)then
784 164 : call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
785 : end if
786 :
787 : ! Now that wavefunctions are initialized, calls of nonlocal operations are possible, so we start the counting (if enabled)
788 5629 : if (dtset%useylm==1.and.dtset%nonlop_ylm_count/=0.and.dtset%paral_kgb==0) then
789 0 : call nonlop_ylm_init_counters()
790 : end if
791 : ! Same for fft counters
792 5629 : if (dtset%fft_count/=0.and.dtset%paral_kgb==0) then
793 0 : call fft_init_counters()
794 : end if
795 : ABI_NVTX_END_RANGE()
796 :
797 : !###########################################################
798 : !### 06. Operations related to restartxf (Old version)
799 :
800 : !TODO: Remove ab_xfh
801 : !Initialize xf history (should be put in inwffil)
802 5629 : ab_xfh%nxfh=0
803 5629 : if(dtset%restartxf>=1 .and. dtfil%ireadwf==1)then
804 :
805 : ! Should exchange the data about history in parallel localrdwf==0
806 0 : if(xmpi_paral==1 .and. dtset%localrdwf==0)then
807 : write(msg, '(a,a,a)' )&
808 0 : & 'It is not yet possible to use non-zero restartxf,',ch10,&
809 0 : & 'in parallel, when localrdwf=0. Sorry for this ...'
810 0 : ABI_BUG(msg)
811 : end if
812 :
813 0 : ABI_MALLOC(ab_xfh%xfhist,(3,dtset%natom+4,2,0))
814 0 : call outxfhist(ab_xfh,dtset%natom,2,wff1,ios)
815 0 : ABI_FREE(ab_xfh%xfhist)
816 :
817 0 : if(ios>0)then
818 : write(msg,'(a,a,a)')&
819 0 : & 'An error occurred reading the input wavefunction file,',ch10,&
820 0 : & 'with restartxf=1.'
821 0 : ABI_ERROR(msg)
822 0 : else if(ios==0)then
823 0 : write(msg, '(a,a,i4,a)' )ch10,&
824 0 : & ' gstate : reading',ab_xfh%nxfh,' (x,f) history pairs from input wf file.'
825 0 : call wrtout([std_out, ab_out], msg)
826 : end if
827 : ! WARNING : should check that restartxf is not negative
828 : ! WARNING : should check that restartxf /= only when dtfil%ireadwf is activated
829 : end if
830 :
831 : !Allocate the xf history array : takes into account the existing
832 : !pairs, minus those that will be discarded, then those that will
833 : !be computed, governed by dtset%ntime, and some additional pairs
834 : !(needed when it will be possible to use xfhist for move.f)
835 5629 : ab_xfh%mxfh=(ab_xfh%nxfh-dtset%restartxf+1)+dtset%ntime+5
836 28145 : ABI_MALLOC(ab_xfh%xfhist,(3,dtset%natom+4,2,ab_xfh%mxfh))
837 2549693 : ab_xfh%xfhist(:,:,:,:) = zero
838 : !WARNING : should check that the number of atoms in the wf file and natom are the same
839 :
840 : !Initialize the xf history array
841 5629 : if(ab_xfh%nxfh>=dtset%restartxf .and. ab_xfh%nxfh>0)then
842 : ! Eventually skip some of the previous history
843 0 : if(dtset%restartxf>=2)then
844 0 : do ixfh=1,dtset%restartxf-1
845 0 : call WffReadSkipRec(ios,1,wff1)
846 : end do
847 : end if
848 :
849 : ! Read and store the relevant history
850 0 : ab_xfh%nxfhr=ab_xfh%nxfh-dtset%restartxf+1
851 0 : call outxfhist(ab_xfh,dtset%natom,3,wff1,ios)
852 : end if
853 :
854 : !Close wff1, if it was ever opened (in inwffil)
855 5629 : if (dtfil%ireadwf==1) then
856 1214 : call WffClose(wff1,ierr)
857 : end if
858 :
859 : !###########################################################
860 : !### 07. Calls setup2
861 :
862 : !Further setup
863 16887 : ABI_MALLOC(start,(3,dtset%natom))
864 5629 : call setup2(dtset,npwtot,start,wvl%wfs,xred)
865 :
866 : !Allocation of previous atomic positions
867 5629 : if (scf_history%history_size>=0) then
868 1454 : xred_old => scf_history%xred_last
869 : else
870 12525 : ABI_MALLOC(xred_old,(3,dtset%natom))
871 : end if
872 48053 : if (has_to_init) xred_old=xred
873 :
874 : !Initialize (eventually) extfpmd object
875 5629 : if(dtset%useextfpmd>=1) then
876 7 : if(extfpmd_chkinp(dtset)) then
877 7 : ABI_MALLOC(extfpmd,)
878 : call extfpmd%init(dtset%mband,hdr%extfpmd_eshift,dtset%extfpmd_nbcut,dtset%extfpmd_nbdbuf,&
879 : & nfftf,dtset%nspden,dtset%nsppol,dtset%nkpt,dtset%occopt,rprimd,dtset%tphysel,&
880 7 : & dtset%tsmear,dtset%useextfpmd,mpi_enreg,dtset%extfpmd_nband,dtset%extfpmd_pawsph)
881 : end if
882 : end if
883 :
884 : !Timing for initialisation period
885 5629 : call timab(1211,2,tsec)
886 5629 : call timab(1213,3,tsec)
887 :
888 :
889 : !###########################################################
890 : !### 08. Compute new occupation numbers
891 :
892 : !Compute new occupation numbers, in case wavefunctions and eigenenergies
893 : !were read from disk, occupation scheme is metallic (this excludes iscf=-1),
894 : !and occupation numbers are required by iscf
895 : if( dtfil%ireadwf==1 .and. &
896 : & (dtset%occopt>=3.and.dtset%occopt<=9) .and. &
897 5629 : & (dtset%iscf>0 .or. dtset%iscf==-3) .and. dtset%positron/=1 ) then
898 :
899 879 : ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
900 : ! Warning : ideally, results_gs%entropy_ks should not be set up here XG 20011007
901 : ! Do not take into account the possible STM bias
902 : call newocc(doccde,eigen,results_gs%energies%entropy_ks,&
903 : & results_gs%energies%e_fermie,results_gs%energies%e_fermih,dtset%ivalence,&
904 : & dtset%spinmagntarget,dtset%mband,dtset%nband,&
905 : & dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,&
906 : & dtset%occopt,dtset%prtvol,dtset%tphysel,dtset%tsmear,dtset%wtk,&
907 293 : & extfpmd=extfpmd,rcpaw=rcpaw)
908 : if (dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(args_gs%upawu(:))>=tol8.or. &
909 965 : & sum(args_gs%jpawu(:))>tol8).and.dtset%dmft_entropy==0) results_gs%energies%entropy_ks=zero
910 :
911 293 : if(associated(extfpmd)) then
912 : ! Get nelect to build density
913 0 : extfpmd%nelect=zero
914 : call extfpmd%compute_nelect(results_gs%energies%e_fermie,dtset%nband,extfpmd%nelect,&
915 0 : & dtset%nkpt,dtset%nspinor,dtset%nsppol,dtset%wtk)
916 : end if
917 293 : ABI_FREE(doccde)
918 :
919 : ! Transfer occupations to bigdft object:
920 293 : if(dtset%usewvl==1 .and. .not. wvlbigdft) then
921 0 : call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,1,wvl%wfs)
922 : ! call wvl_energies_abi2big(results_gs%energies,wvl%wfs,2)
923 : end if
924 :
925 : else
926 : ! Warning : ideally, results_gs%entropy should not be set up here XG 20011007
927 5336 : results_gs%energies%entropy=zero
928 : end if
929 :
930 : !###########################################################
931 : !### 09. Generate an index table of atoms
932 :
933 : !Definition of atindx array
934 : !Generate an index table of atoms, in order for them to be used type after type.
935 16887 : ABI_MALLOC(atindx,(dtset%natom))
936 11258 : ABI_MALLOC(atindx1,(dtset%natom))
937 16887 : ABI_MALLOC(nattyp,(psps%ntypat))
938 5629 : indx=1
939 13620 : do itypat=1,psps%ntypat
940 7991 : nattyp(itypat)=0
941 34079 : do iatom=1,dtset%natom
942 28450 : if(dtset%typat(iatom)==itypat)then
943 13236 : atindx(iatom)=indx
944 13236 : atindx1(indx)=iatom
945 13236 : indx=indx+1
946 13236 : nattyp(itypat)=nattyp(itypat)+1
947 : end if
948 : end do
949 : end do
950 :
951 : !Compute structure factor phases for current atomic pos:
952 5629 : if ((.not.read_wf_or_den).or.(scf_history%history_size>0.and.has_to_init)) then
953 12621 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*dtset%natom))
954 4207 : call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
955 : end if
956 :
957 : !Here allocation of GPU for vtorho calculations
958 : #if defined HAVE_GPU
959 : if (dtset%gpu_option/=ABI_GPU_DISABLED) then
960 : call alloc_hamilt_gpu(atindx1,dtset,gprimd,mpi_enreg,nattyp,npwarr,2,psps,dtset%gpu_option)
961 : end if
962 : #endif
963 :
964 : !###########################################################
965 : !### 10. PAW related operations
966 :
967 : !Initialize paw_dmft, even if neither dmft nor paw are used
968 : !write(std_out,*) "dtset%usedmft",dtset%usedmft
969 5629 : use_sc_dmft=dtset%usedmft
970 : ! if(dtset%paral_kgb>0) use_sc_dmft=0
971 : !Should be done inside init_sc_dmft
972 5629 : if ( dtset%usedmft /= 0 .and. dtset%dmft_entropy > 0 .and. dtset%usedmft /= 10) then
973 : call data4entropyDMFT_init(paw_dmft%forentropyDMFT,&
974 : dtset%natom,&
975 : dtset%typat,&
976 : dtset%lpawu,&
977 : dtset%dmft_t2g==1, &
978 : args_gs%upawu,&
979 3 : args_gs%jpawu)
980 : end if
981 : !write(std_out,*) "paw_dmft%use_dmft",paw_dmft%use_dmft
982 :
983 : !PAW: 1- Initialize values for several arrays unchanged during iterations
984 : !2- Initialize data for DFT+U
985 : !3- Eventually open temporary storage file
986 5629 : if(psps%usepaw==1) then
987 : ! 1-
988 1285 : gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2
989 :
990 : ! Test if we have to call pawinit
991 : ! Some gen-cond have to be added...
992 1285 : call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
993 :
994 1285 : if (psp_gencond==1.or.call_pawinit) then
995 534 : call timab(553,1,tsec)
996 534 : gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
997 534 : hyb_range_fock=zero;if (dtset%ixc<0) call libxc_functionals_get_hybridparams(hyb_range=hyb_range_fock)
998 : call pawinit(dtset%effmass_free,gnt_option,gsqcut_shp,hyb_range_fock,dtset%pawlcutd,dtset%pawlmix,&
999 : & psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,&
1000 534 : & pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero)
1001 :
1002 : ! Update internal values
1003 534 : call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
1004 534 : call timab(553,2,tsec)
1005 : #if defined HAVE_BIGDFT
1006 : ! In the PAW+WVL case, copy sij:
1007 : if(dtset%usewvl==1) then
1008 : do itypat=1,dtset%ntypat
1009 : wvl%descr%paw%sij(1:pawtab(itypat)%lmn2_size,itypat)=pawtab(itypat)%sij(:)
1010 : end do
1011 : end if
1012 : #endif
1013 : end if
1014 3094 : psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore)
1015 1285 : call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot)
1016 : ! 2-Initialize and compute data for DFT+U, EXX, or DFT+DMFT
1017 : call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,&
1018 : & is_dfpt,args_gs%jpawu,dtset%lexexch,dtset%lpawu,dtset%nspinor,dtset%ntypat,dtset%optdcmagpawu,pawang,dtset%pawprtvol,&
1019 : & pawrad,pawtab,args_gs%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu,ucrpa=dtset%ucrpa,dmft_orbital=dtset%dmft_orbital(:),&
1020 : & dmft_dc=dtset%dmft_dc,dmft_orbital_filepath=dtset%dmft_orbital_filepath,dmft_yukawa_param=dtset%dmft_yukawa_param,&
1021 1285 : & dmft_yukawa_lambda=dtset%dmft_yukawa_lambda,dmft_yukawa_epsilon=dtset%dmft_yukawa_epsilon)
1022 :
1023 : ! DEBUG:
1024 : !if (me == master) call pawtab_print(Pawtab)
1025 : end if
1026 :
1027 : call init_sc_dmft(dtset,psps%mpsang,paw_dmft,gprimd(:,:),kg(:,:),mpi_enreg,npwarr(:),occ(:),pawang, &
1028 5629 : & pawrad(:),pawtab(:),rprimd(:,:),ucvol,dtfil%unpaw,use_sc_dmft,xred(:,:),ylm(:,:))
1029 5629 : if (paw_dmft%use_dmft == 1) then
1030 43 : if (paw_dmft%myproc == 0) then
1031 22 : call readocc_dmft(paw_dmft,dtfil%filnam_ds(3),dtfil%filnam_ds(4))
1032 : end if
1033 43 : if (paw_dmft%dmft_read_occnd /= 0) then
1034 5 : call xmpi_bcast(paw_dmft%occnd(:,:,:,:,:),0,paw_dmft%spacecomm,ierr)
1035 : end if
1036 43 : call print_sc_dmft(paw_dmft,dtset%pawprtvol)
1037 : end if
1038 :
1039 : ! Initialize (eventually) rcpaw object
1040 5629 : if (dtset%use_rcpaw==1) then
1041 3 : ABI_WARNING("Untested Mode RCPAW")
1042 3 : ABI_MALLOC(rcpaw,)
1043 3 : call rcpaw_init(rcpaw,dtset,psps%filpsp,pawrad,pawtab,psps%ntypat,1,.false.,my_natom,mpi_enreg%comm_atom,mpi_enreg%my_atmtab)
1044 : end if
1045 :
1046 : !###########################################################
1047 : !### 11. Initialize (eventually) electron-positron data and
1048 : !### electric and magnetic field data
1049 :
1050 : !Initialize (eventually) electron-positron data
1051 5629 : nullify (electronpositron)
1052 5629 : if (dtset%positron/=0) then
1053 14 : call init_electronpositron(dtfil%ireadwf,dtset,electronpositron,mpi_enreg,nfftf,pawrhoij,pawtab)
1054 : end if
1055 :
1056 : !###########################################################
1057 : ! Initialisation of cprj
1058 :
1059 : ! xg_nonlop available only for cprj_in_memory=1 and (LOBPCG or Chebfi or Slicing)
1060 : ! cprj_in_memory=2 is used for Congugate Gradient
1061 5629 : if (dtset%cprj_in_memory==1) then
1062 162 : if (dtset%useylm/=1) then
1063 0 : ABI_ERROR('xg_nonlop cannot be used with useylm/=1')
1064 : end if
1065 : call xg_nonlop_init(xg_nonlop,psps%indlmn,my_natom,nattyp,dtset%mkmem,dtset%ntypat,&
1066 : dtset%nspinor,ucvol,dtset%usepaw,dtset%xg_nonlop_option,&
1067 : mpi_enreg%me_band,mpi_enreg%comm_band,mpi_enreg%comm_atom,&
1068 162 : mpi_atmtab=mpi_enreg%my_atmtab)
1069 162 : if (xg_nonlop%paw) then
1070 138 : inv_sij=dtset%wfoptalg==111.or.dtset%wfoptalg==112
1071 138 : call xg_nonlop_make_Sij(xg_nonlop,pawtab,inv_sij=inv_sij)
1072 : else
1073 24 : call xg_nonlop_make_ekb(xg_nonlop,psps%ekb)
1074 : end if
1075 : end if
1076 :
1077 5629 : usecprj=0; mcprj=0;mband_cprj=0
1078 5629 : compute_cprj=.false.
1079 5629 : if (dtset%cprj_in_memory==2) then
1080 : compute_cprj=.true.
1081 : usecprj=1
1082 : else
1083 5629 : if (dtset%usepaw==1) then
1084 1285 : if (associated(electronpositron)) then
1085 13 : if (dtset%positron/=0.and.electronpositron%dimcprj>0) usecprj=1
1086 : end if
1087 1285 : if (dtset%prtnabla>0) usecprj=1
1088 1285 : if (dtset%extrapwf>0) usecprj=1
1089 1285 : if (dtset%pawfatbnd>0)usecprj=1
1090 1285 : if (dtset%prtdos==3) usecprj=1
1091 1285 : if (dtset%usewvl==1) usecprj=1
1092 1285 : if (dtset%nstep==0) usecprj=0
1093 1285 : if (dtset%usefock==1) usecprj=1
1094 : end if
1095 : end if
1096 1270 : if (usecprj==0) then
1097 5570 : ABI_MALLOC(cprj,(0,0))
1098 : end if
1099 : if (usecprj==1) then
1100 59 : mband_cprj=dtset%mband;if (dtset%paral_kgb/=0) mband_cprj=mband_cprj/mpi_enreg%nproc_band
1101 59 : mcprj=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
1102 : !Was allocated above for valgrind sake so should always be true (safety)
1103 59 : if (allocated(cprj)) then
1104 0 : call pawcprj_free(cprj)
1105 0 : ABI_FREE(cprj)
1106 : end if
1107 16096 : ABI_MALLOC(cprj,(dtset%natom,mcprj))
1108 59 : ncpgr=0
1109 59 : if (dtset%usefock==1) then ! Note that compute_cprj = false if usefock/=0
1110 15 : if (dtset%optforces == 1) then
1111 5 : ncpgr = 3
1112 : end if
1113 : ! if (dtset%optstress /= 0) then
1114 : ! ncpgr = 6 ; ctocprj_choice = 3
1115 : ! end if
1116 : end if
1117 177 : ABI_MALLOC(dimcprj_srt,(dtset%natom))
1118 59 : call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
1119 59 : call pawcprj_alloc(cprj,ncpgr,dimcprj_srt)
1120 59 : if (compute_cprj) then
1121 0 : choice = 1 ! no derivative...
1122 0 : idir = 0 ! ...so no direction
1123 0 : iatom = 0 ! all atoms
1124 0 : iorder_cprj = 0 ! ordered by atom types
1125 0 : ncprj = dtset%natom
1126 : ! Compute structure factor phases and large sphere cut-off (gsqcut):
1127 0 : ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
1128 0 : call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
1129 0 : call wrtout(std_out,' Computing cprj from initial wavefunctions (gstate)')
1130 : call ctocprj(atindx,cg,choice,cprj,gmet,gprimd,iatom,idir,&
1131 : & iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,&
1132 : & dtset%mpw,dtset%natom,nattyp,dtset%nband,ncprj,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,&
1133 : & dtset%nsppol,dtset%nsppol,psps%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,dtfil%unpaw,&
1134 0 : & xred,ylm,ylmgr)
1135 0 : call wrtout(std_out,' cprj is computed')
1136 0 : ABI_FREE(ph1d)
1137 : end if
1138 : end if
1139 :
1140 : !Timing for initialisation period
1141 5629 : call timab(1213,2,tsec)
1142 5629 : call timab(1214,3,tsec)
1143 :
1144 : !###########################################################
1145 : !### 12. Operations dependent of iscf value
1146 :
1147 : !Get starting charge density : rhor as well as rhog
1148 : !Also initialize the kinetic energy density
1149 5629 : if (scf_history%history_size>=0) then
1150 1454 : rhor => scf_history%rhor_last
1151 1454 : taur => scf_history%taur_last
1152 : else
1153 16700 : ABI_MALLOC(rhor,(nfftf,dtset%nspden))
1154 16700 : ABI_MALLOC(taur,(nfftf,dtset%nspden*dtset%usekden))
1155 : end if
1156 16887 : ABI_MALLOC(rhog,(2,nfftf))
1157 16887 : ABI_MALLOC(taug,(2,nfftf*dtset%usekden))
1158 :
1159 5629 : if (has_to_init) then
1160 :
1161 : ! === Self-consistent case
1162 4492 : if (dtset%iscf>0 .or. (dtset%iscf==0 .and. dtset%usewvl==1 )) then
1163 :
1164 : ! >>> Initialize charge density
1165 :
1166 3981 : if (dtfil%ireadden/=0.and.dtset%positron<=0) then
1167 : ! Choice 1: read charge density from file
1168 56 : rdwrpaw=psps%usepaw ; if(dtfil%ireadwf/=0) rdwrpaw=0
1169 56 : if (dtset%usewvl==0) then
1170 : call read_rhor(dtfil%fildensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw, &
1171 56 : mpi_enreg,rhor,hdr_den,pawrhoij,comm,check_hdr=hdr,allow_interp=.True.)
1172 56 : results_gs%etotal = hdr_den%etot
1173 56 : call hdr_den%free()
1174 : else
1175 0 : fform=52 ; accessfil=0
1176 0 : if (dtset%iomode == IO_MODE_MPI ) accessfil=4
1177 0 : if (dtset%iomode == IO_MODE_ETSF) accessfil=3
1178 : call ioarr(accessfil,rhor,dtset,results_gs%etotal,fform,dtfil%fildensin,hdr,&
1179 0 : mpi_enreg,ngfftf,cplex1,nfftf,pawrhoij,1,rdwrpaw,wvl%den)
1180 : end if
1181 56 : if (rdwrpaw/=0) then
1182 : call hdr%update(bantot,etot,fermie,fermih,residm,&
1183 : rprimd,occ,pawrhoij,xred,args_gs%amu,&
1184 17 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1185 : end if
1186 : ! Compute up+down rho(G) by fft
1187 56 : call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0)
1188 :
1189 3925 : else if (dtfil%ireadwf/=0) then
1190 : ! Choice 2: obtain charge density from wfs that were read previously
1191 : ! Warning: in PAW, rho does not include the compensation density (added later)
1192 855 : tim_mkrho=1
1193 855 : if (psps%usepaw==1) then
1194 420 : ABI_MALLOC(rhowfg,(2,dtset%nfft))
1195 560 : ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden))
1196 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1197 140 : & rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,extfpmd=extfpmd)
1198 140 : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor)
1199 140 : ABI_FREE(rhowfg)
1200 140 : ABI_FREE(rhowfr)
1201 : else
1202 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1203 715 : & rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,extfpmd=extfpmd)
1204 : end if
1205 :
1206 3070 : else if (dtfil%ireadwf==0.and.dtset%positron/=1) then
1207 : ! Choice 3: crude, but realistic initialisation of the charge density
1208 : ! There is not point to compute it from random wavefunctions
1209 3065 : if (dtset%usewvl == 0) then
1210 : call initro(atindx,dtset%densty,gmet,gsqcut_eff,psps%usepaw,&
1211 : & mgfftf,mpi_enreg,psps%mqgrid_vl,dtset%natom,nattyp,nfftf,&
1212 : & ngfftf,dtset%nspden,psps%ntypat,psps,pawtab,ph1df,&
1213 : & psps%qgrid_vl,rhog,rhor,dtset%spinat,ucvol,psps%usepaw,&
1214 3065 : & dtset%ziontypat,dtset%znucl)
1215 : else
1216 0 : if (dtset%usepaw==0) then
1217 : !Wavelet density corresponds exactly to the wavefunctions,
1218 : !since wavefunctions are taken from diagonalisation of LCAO.
1219 0 : call wvl_mkrho(dtset,irrzon,mpi_enreg,phnons,rhor,wvl%wfs,wvl%den)
1220 : else
1221 : #if defined HAVE_BIGDFT
1222 : call wvl_initro(atindx1,wvl%descr%atoms%astruct%geocode,wvl%descr%h,&
1223 : & mpi_enreg%me_wvl,dtset%natom,nattyp,nfftf,dtset%nspden,psps%ntypat,&
1224 : & wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,&
1225 : & wvl%descr%Glr%d%n3,pawrad,pawtab,psps%gth_params%psppar,rhor,rprimd,&
1226 : & dtset%spinat,wvl%den,dtset%xc_denpos,xred,dtset%ziontypat)
1227 : call wvl_mkrho(dtset,irrzon,mpi_enreg,phnons,rhor,wvl%wfs,wvl%den)
1228 : #endif
1229 : end if ! usepaw
1230 : end if ! usewvl
1231 : ! Update initialized density taking into account jellium slab
1232 3065 : if (dtset%jellslab/=0) then
1233 6 : option=2
1234 18 : ABI_MALLOC(work,(nfftf))
1235 : call jellium(gmet,gsqcut_eff,mpi_enreg,nfftf,ngfftf,dtset%nspden,option,&
1236 6 : & dtset%slabwsrad,rhog,rhor,rprimd,work,dtset%slabzbeg,dtset%slabzend)
1237 6 : ABI_FREE(work)
1238 : end if ! of usejell
1239 :
1240 : end if ! choice for charge density initialization
1241 :
1242 : ! >>> Initialize kinetic energy density
1243 3981 : if (dtset%usekden==1) then
1244 :
1245 73 : if (dtfil%ireadkden/=0.and.dtset%positron<=0) then
1246 : ! Choice 1: read kinetic energy density from file
1247 2 : rdwrpaw=0
1248 : call read_rhor(dtfil%filkdensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw, &
1249 2 : mpi_enreg,taur,hdr_den,pawrhoij,comm,check_hdr=hdr,allow_interp=.True.)
1250 2 : call hdr_den%free()
1251 : ! Compute up+down tau(G) by fft
1252 2 : call fourdp(1,taug,taur(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0)
1253 :
1254 71 : else if (dtfil%ireadwf/=0) then
1255 : ! Choice 2: obtain kinetic energy density from wfs that were read previously
1256 13 : tim_mkrho=1
1257 13 : if (psps%usepaw==1) then
1258 0 : ABI_MALLOC(rhowfg,(2,dtset%nfft))
1259 0 : ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden))
1260 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1261 0 : & rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
1262 0 : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,taug,rhowfr,taur)
1263 0 : ABI_FREE(rhowfg)
1264 0 : ABI_FREE(rhowfr)
1265 : else
1266 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1267 13 : & taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
1268 : end if
1269 :
1270 58 : else if(dtfil%ireadwf==0.and.dtset%positron/=1)then
1271 : ! Choice 3: kinetic energy density initialized to zero (?)
1272 12007059 : taur=zero ; taug=zero
1273 :
1274 : end if ! choice for kinetic energy density initialization
1275 : end if ! usekden
1276 :
1277 : ! === Non self-consistent case
1278 511 : else if ((dtset%iscf==-1.or.dtset%iscf==-2.or.dtset%iscf==-3).and.dtset%positron<=0) then
1279 :
1280 : ! Read density from a disk file (this is mandatory for non-self-consistent calculations)
1281 : ! Note : results_gs%etotal is read here,
1282 : ! and might serve in the tddft routine, but it is contrary to the intended use of results_gs ...
1283 : ! Warning : should check the use of results_gs%e_fermie
1284 : ! Warning : should check the use of results_gs%residm
1285 : ! One might make them separate variables.
1286 :
1287 : ! Read charge density and get Fermi level from hdr_den
1288 511 : rdwrpaw=psps%usepaw
1289 : call read_rhor(dtfil%fildensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw,&
1290 511 : & mpi_enreg,rhor,hdr_den,pawrhoij,comm,check_hdr=hdr)
1291 511 : results_gs%etotal = hdr_den%etot;
1292 511 : results_gs%energies%e_fermie = hdr_den%fermie
1293 511 : results_gs%energies%e_fermih = hdr_den%fermih
1294 : ! Compute up+down rho(G) by fft
1295 511 : call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0)
1296 511 : call hdr_den%free()
1297 :
1298 : ! Read kinetic energy density
1299 511 : if(dtset%usekden==1)then
1300 1 : rdwrpaw=0
1301 : call read_rhor(dtfil%filkdensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw, &
1302 1 : mpi_enreg,taur,hdr_den,pawrhoij,comm,check_hdr=hdr)
1303 1 : call hdr_den%free()
1304 : ! Compute up+down tau(G) by fft
1305 1 : call fourdp(1,taug,taur(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0)
1306 : end if
1307 :
1308 : end if ! self-consistent/non self-consistent
1309 : end if ! has_to_init
1310 :
1311 : !Timing for initialisation period
1312 5629 : call timab(1214,2,tsec)
1313 5629 : call timab(1215,3,tsec)
1314 :
1315 : !###########################################################
1316 : !### 13. If needed, initialize SCF history variables
1317 :
1318 : !If needed, initialize atomic density in SCF history
1319 5629 : if (scf_history%history_size>0.and.has_to_init) then
1320 : ! If rhor is an atomic density, just store it in history
1321 109 : if (.not.read_wf_or_den) then
1322 9343117 : scf_history%atmrho_last(:)=rhor(:,1)
1323 : else
1324 : ! If rhor is not an atomic density, has to compute rho_at(r)
1325 0 : ABI_MALLOC(rhowfg,(2,nfftf))
1326 0 : ABI_MALLOC(rhowfr,(nfftf,1))
1327 0 : ABI_MALLOC(spinat_dum,(3,dtset%natom))
1328 0 : spinat_dum=zero
1329 : call initro(atindx,dtset%densty,gmet,gsqcut_eff,psps%usepaw,mgfftf,mpi_enreg,&
1330 : & psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,1,psps%ntypat,psps,pawtab,&
1331 : & ph1df,psps%qgrid_vl,rhowfg,rhowfr,spinat_dum,ucvol,&
1332 0 : & psps%usepaw,dtset%ziontypat,dtset%znucl)
1333 0 : scf_history%atmrho_last(:)=rhowfr(:,1)
1334 0 : ABI_FREE(rhowfg)
1335 0 : ABI_FREE(rhowfr)
1336 0 : ABI_FREE(spinat_dum)
1337 : end if
1338 : end if
1339 :
1340 5629 : if ((.not.read_wf_or_den).or.(scf_history%history_size>0.and.has_to_init)) then
1341 4207 : ABI_FREE(ph1df)
1342 : end if
1343 :
1344 : !!Electric field: initialization stage
1345 : !!further initialization and updates happen in scfcv.F90
1346 : call init_e_field_vars(dtefield,dtset,gmet,gprimd,kg,&
1347 : & mpi_enreg,npwarr,occ,pawang,pawrad,pawtab,psps,&
1348 5629 : & pwind,pwind_alloc,pwnsfac,rprimd,symrec,xred)
1349 :
1350 5629 : fatvshift=one
1351 :
1352 : !Check whether exiting was required by the user. If found then do not start minimization steps
1353 : !At this first call to chkexi, initialize cpus, if it
1354 : !is non-zero (which would mean that no action has to be taken)
1355 : !Should do this in driver ...
1356 5629 : cpus=dtset%cpus
1357 5629 : if(abs(cpus)>1.0d-5)cpus=cpus+cpui
1358 5629 : openexit=1 ; if(dtset%chkexit==0) openexit=0
1359 5629 : call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit)
1360 :
1361 : !If immediate exit, and wavefunctions were not read, must zero eigenvalues
1362 5629 : if (iexit/=0) eigen(:)=zero
1363 :
1364 : #if defined HAVE_BIGDFT
1365 : if (dtset%usewvl == 1 .and. dtset%timopt==10) then
1366 : call wvl_timing(xmpi_world,'== INITS','PR')
1367 : end if
1368 : #endif
1369 :
1370 5629 : call timab(1215,2,tsec)
1371 5629 : call pstat_proc%print(_PSTAT_ARGS_)
1372 :
1373 5629 : conv_retcode = 0
1374 :
1375 5629 : if (iexit==0) then
1376 :
1377 : ! ###########################################################
1378 : ! ### 14. Move atoms and acell according to ionmov value
1379 :
1380 5629 : call timab(1225,3,tsec)
1381 :
1382 : call scfcv_init(scfcv_args,atindx,atindx1,cg,cprj,cpus,&
1383 : & args_gs%dmatpawu,dtefield,dtfil,dtpawuj,dtset,ecore,eigen,hdr,extfpmd,rcpaw,&
1384 : & indsym,initialized,irrzon,kg,mcg,mcprj,mpi_enreg,my_natom,nattyp,ndtpawuj,&
1385 : & nfftf,npwarr,occ,pawang,pawfgr,pawrad,pawrhoij,&
1386 : & pawtab,phnons,psps,pwind,pwind_alloc,pwnsfac,rec_set,&
1387 : & resid,results_gs,scf_history,fatvshift,&
1388 5629 : & symrec,taug,taur,wvl,ylm,ylmgr,paw_dmft,wffnew,wffnow,xg_nonlop)
1389 :
1390 5629 : call dtfil%init_time(0)
1391 :
1392 455949 : write(msg,'(a,80a)')ch10,('=',mu=1,80)
1393 16887 : call wrtout([std_out, ab_out], msg)
1394 :
1395 5629 : if (dtset%ionmov==0 .or. dtset%imgmov==6) then
1396 :
1397 : ! Should merge this call with the call for dtset%ionmov==4 and 5
1398 5280 : if (dtset%macro_uj==0) then
1399 5270 : itimes(1)=itime0 ; itimes(2)=itimimage_gstate
1400 5270 : call scfcv_run(scfcv_args,electronpositron,itimes,rhog,rhor,rprimd,xred,xred_old,conv_retcode)
1401 : else
1402 : ! Conduct determination of U
1403 10 : call pawuj_drive(scfcv_args,dtset,electronpositron,rhog,rhor,rprimd,xred,xred_old)
1404 : end if
1405 :
1406 : ! ========================================
1407 : ! New structure for geometry optimization
1408 : ! ========================================
1409 349 : else if (dtset%ionmov>50.or.dtset%ionmov<=28) then
1410 :
1411 : ! TODO: return conv_retcode
1412 : call mover(scfcv_args,ab_xfh,acell,args_gs%amu,dtfil,&
1413 349 : & electronpositron,rhog,rhor,rprimd,vel,vel_cell,xred,xred_old,itimimage_gstate=itimimage_gstate)
1414 :
1415 : ! Compute rprim from rprimd and acell
1416 1396 : do kk=1,3
1417 4537 : do jj=1,3
1418 4188 : rprim(jj,kk)=rprimd(jj,kk)/acell(kk)
1419 : end do
1420 : end do
1421 :
1422 : ! =========================================
1423 : ! New structure for geometry optimization
1424 : ! =========================================
1425 :
1426 : else ! Not an allowed option
1427 : write(msg, '(a,i0,2a)' )&
1428 0 : 'Disallowed value for ionmov=',dtset%ionmov,ch10,&
1429 0 : 'Allowed values are: 1,2,3,4,5,6,7,8,9,10,11,12,13,14,20,21,22,23,24,28 and 30'
1430 0 : ABI_BUG(msg)
1431 : end if
1432 :
1433 5629 : call scfcv_destroy(scfcv_args)
1434 :
1435 5629 : call timab(1225,2,tsec)
1436 :
1437 : ! ###########################################################
1438 : ! ### 15. Final operations and output for gstate
1439 :
1440 : end if ! End of the check of hasty exit
1441 :
1442 5629 : call timab(1226,3,tsec)
1443 :
1444 455949 : write(msg, '(80a,a,a,a,a)' ) ('=',mu=1,80),ch10,ch10,&
1445 461578 : & ' ----iterations are completed or convergence reached----',ch10
1446 16887 : call wrtout([std_out, ab_out], msg)
1447 :
1448 : !Mark this GS computation as done
1449 5629 : initialized=1
1450 :
1451 : !Update the header, before using it
1452 : call hdr%update(bantot,results_gs%etotal,results_gs%energies%e_fermie,results_gs%energies%e_fermih,&
1453 : results_gs%residm,rprimd,occ,pawrhoij,xred,args_gs%amu,&
1454 5629 : comm_atom=mpi_enreg%comm_atom,extfpmd_eshift=results_gs%extfpmd_eshift,mpi_atmtab=mpi_enreg%my_atmtab)
1455 :
1456 16887 : ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
1457 449988 : doccde=zero
1458 :
1459 : call ebands%init(bantot, dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
1460 : doccde,eigen,hdr%istwfk,hdr%kptns,hdr%nband,&
1461 : hdr%nkpt,hdr%npwarr,hdr%nsppol,hdr%nspinor,hdr%tphysel,hdr%tsmear,hdr%occopt,hdr%occ,hdr%wtk,&
1462 : hdr%cellcharge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, &
1463 5629 : hdr%kptrlatt, hdr%nshiftk, hdr%shiftk)
1464 :
1465 5629 : ebands%fermie = results_gs%energies%e_fermie
1466 5629 : ebands%fermih = results_gs%energies%e_fermih
1467 5629 : ABI_FREE(doccde)
1468 :
1469 : ! Compute and print the gaps.
1470 5629 : call ebands%report_gap(header="Gap info",unit=std_out,mode_paral="COLL",gaps=results_gs%gaps)
1471 :
1472 5629 : call timab(1226,2,tsec)
1473 5629 : call timab(1227,3,tsec)
1474 :
1475 5629 : if(dtset%nqpt==0)filnam=dtfil%fnameabo_wfk
1476 5629 : if(dtset%nqpt==1)filnam=dtfil%fnameabo_wfq
1477 :
1478 : ! Write wavefunctions file only if convergence was not achieved.
1479 : !write(std_out,*)"conv_retcode", conv_retcode
1480 5629 : write_wfk = .True.
1481 5629 : if (dtset%prtwf==-1 .and. conv_retcode == 0) then
1482 0 : write_wfk = .False.
1483 0 : msg = "GS calculation converged with prtwf=-1 --> Skipping WFK file output"
1484 0 : call wrtout(ab_out, msg)
1485 0 : ABI_COMMENT(msg)
1486 : end if
1487 :
1488 : !To print out the WFs, need the rprimd that was used to generate the G vectors
1489 73177 : hdr%rprimd=rprimd_for_kg
1490 :
1491 5629 : if (write_wfk) then
1492 5629 : call outresid(dtset,dtset%kptns,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,resid)
1493 : call outwf(cg,dtset,psps,eigen,filnam,hdr,kg,dtset%kptns,&
1494 : dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%natom,&
1495 : dtset%nband,dtset%nkpt,npwarr,dtset%nsppol,&
1496 5629 : occ,response,dtfil%unwff2,wvl%wfs,wvl%descr)
1497 :
1498 : ! Generate WFK with k-mesh from WFK containing list of k-points inside pockets.
1499 5629 : if (dtset%getkerange_filepath /= ABI_NOFILE) then
1500 2 : call wfk_klist2mesh(dtfil%fnameabo_wfk, dtset%getkerange_filepath, dtset, comm)
1501 : end if
1502 :
1503 : !SPr: add input variable managing the .vtk file OUTPUT (Please don't remove the next commented line)
1504 : !call printmagvtk(mpi_enreg,cplex1,dtset%nspden,nfftf,ngfftf,rhor,rprimd,'DEN')
1505 : end if
1506 :
1507 5629 : if (dtset%prtwf==2) call outqmc(cg,dtset,eigen,gprimd,hdr,kg,mcg,mpi_enreg,npwarr,occ,psps,results_gs)
1508 :
1509 : !Restore the original rprimd in hdr
1510 73177 : hdr%rprimd=rprimd
1511 :
1512 : ! Generate WFK in full BZ (needed by LOBSTER)
1513 5629 : if (me == master .and. dtset%prtwf == 1 .and. dtset%prtwf_full == 1 .and. dtset%nqpt == 0) then
1514 1 : wfkfull_path = strcat(dtfil%filnam_ds(4), "_FULL_WFK")
1515 1 : if (dtset%iomode == IO_MODE_ETSF) wfkfull_path = nctk_ncify(wfkfull_path)
1516 1 : call wfk_to_bz(filnam, dtset, psps, pawtab, wfkfull_path, hdr_bz, ebands_bz)
1517 1 : call hdr_bz%free(); call ebands_bz%free(); call cryst%free()
1518 : end if
1519 :
1520 5629 : call timab(1227,2,tsec)
1521 5629 : call timab(1228,3,tsec)
1522 :
1523 : call clnup1(acell,dtset,eigen,results_gs%energies%e_fermie,results_gs%energies%e_fermih,&
1524 : & dtfil%fnameabo_dos,dtfil%fnameabo_eig,results_gs%gred,&
1525 : & mpi_enreg,nfftf,ngfftf,occ,dtset%optforces,&
1526 5629 : & resid,rhor,rprimd,results_gs%vxcavg,xred,rcpaw)
1527 :
1528 5629 : if ( (dtset%iscf>=0 .or. dtset%iscf==-3) .and. dtset%prtstm==0) then
1529 5142 : call prtene(dtset,results_gs%energies,ab_out,psps%usepaw)
1530 : end if
1531 :
1532 5629 : call timab(1228,2,tsec)
1533 5629 : call timab(1229,3,tsec)
1534 :
1535 : !write final electric field components HONG
1536 :
1537 : if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt ==7 .or. &
1538 5629 : & dtset%berryopt == 14 .or. dtset%berryopt == 16 .or. dtset%berryopt ==17 ) then ! output final electric field data !!HONG
1539 : if (dtset%berryopt == 4) then
1540 25 : write(msg,'(a,a)') ch10, 'Constant unreduced E calculation - final values:'
1541 : else if (dtset%berryopt == 6 ) then
1542 0 : write(msg,'(a,a)') ch10, 'Constant unreduced D calculation - final values:'
1543 : else if (dtset%berryopt == 14) then
1544 0 : write(msg,'(a,a)') ch10, 'Constant reduced ebar calculation - final values:'
1545 : else if (dtset%berryopt == 16 ) then
1546 0 : write(msg,'(a,a)') ch10, 'Constant reduced d calculation - final values:'
1547 : else if (dtset%berryopt == 17) then
1548 0 : write(msg,'(a,a)') ch10, 'Constant reduced ebar and d calculation - final values:'
1549 : end if
1550 :
1551 75 : call wrtout([std_out, ab_out], msg)
1552 25 : call prtefield(dtset,dtefield,ab_out,rprimd)
1553 25 : call prtefield(dtset,dtefield,std_out,rprimd)
1554 :
1555 : ! To check if the final electric field is below the critical field
1556 100 : do kk = 1, 3
1557 100 : efield_band(kk) = abs(dtset%red_efieldbar(kk))*dtefield%nkstr(kk)
1558 : end do
1559 : ! eg = maxval(eg_dir)
1560 : ! eg_ev = eg*Ha_eV
1561 25 : write(msg,'(a,a,a,a,a,a,a,a,f7.2,a,a)')ch10,&
1562 25 : & ' Please check: COMMENT - ',ch10,&
1563 25 : & ' As a rough estimate,',ch10,&
1564 25 : & ' to be below the critical field, the bandgap of your system',ch10,&
1565 150 : & ' should be larger than ',maxval(efield_band)*Ha_eV,' eV.',ch10
1566 75 : call wrtout([std_out, ab_out], msg)
1567 :
1568 25 : write(msg,'(a)') '--------------------------------------------------------------------------------'
1569 75 : call wrtout([std_out, ab_out], msg)
1570 : end if
1571 :
1572 5629 : call timab(1229,2,tsec)
1573 5629 : call timab(1230,3,tsec)
1574 :
1575 : !In the // case, only master writes the energy and the gradients to the DDB
1576 5629 : if (me==0.and.dtset%nimage==1.and.((dtset%iscf > 0).or.&
1577 : & (dtset%berryopt == -1).or.(dtset%berryopt) == -3)) then
1578 :
1579 : ! DDB dimensions
1580 2803 : if (dtset%iscf > 0) then
1581 2796 : nblok = 2 ! 1st blok = gradients, 2nd blok = energy
1582 : else
1583 7 : nblok = 1 ! 1st blok = gradients
1584 : end if
1585 2803 : mpert = dtset%natom + 6
1586 :
1587 : ! Create header and ddb objects
1588 2803 : dscrpt=' Note : temporary (transfer) database '
1589 : call ddb_hdr%init(dtset,psps,pawtab,dscrpt,nblok,&
1590 2803 : & xred=xred,occ=occ,ngfft=ngfft)
1591 :
1592 2803 : call ddb%init(dtset, nblok, mpert, with_d1E=.true.)
1593 :
1594 : ! Set the electronic polarization
1595 2803 : if ((abs(dtset%berryopt) == 1).or.(abs(dtset%berryopt) == 3)) then
1596 45 : call ddb%set_pel(results_gs%pel, dtset%rfdir, 1)
1597 : end if
1598 :
1599 2803 : if (dtset%iscf > 0) then
1600 :
1601 : ! Set the gradients (forces) in reduced coordinates
1602 2796 : call ddb%set_gred(results_gs%gred, 1)
1603 :
1604 : ! Set the stress tensor
1605 2796 : call ddb%set_strten(results_gs%strten, 1)
1606 :
1607 : ! Set the total energy
1608 2796 : call ddb%set_etotal(results_gs%etotal, 2)
1609 : end if
1610 :
1611 2803 : if (dtset%prtddb==1) then
1612 : ! Write the DDB
1613 2760 : call ddb%write(ddb_hdr, dtfil%fnameabo_ddb, with_psps=0)
1614 : end if
1615 :
1616 : ! Free memory
1617 2803 : call ddb_hdr%free()
1618 2803 : call ddb%free()
1619 :
1620 : end if
1621 :
1622 5629 : call timab(1230,2,tsec)
1623 5629 : call timab(1231,3,tsec)
1624 :
1625 :
1626 5629 : if (dtset%nstep>0 .and. dtset%prtstm==0 .and. dtset%positron/=1) then
1627 : call clnup2(psps%n1xccc,results_gs%gred,results_gs%grchempottn,results_gs%gresid,&
1628 : & results_gs%grewtn,results_gs%grvdw,results_gs%grxc,dtset%iscf,dtset%natom,&
1629 : & results_gs%ngrvdw,dtset%optforces,dtset%optstress,dtset%prtvol,start,&
1630 5586 : & results_gs%strten,results_gs%synlgr,xred)
1631 : end if
1632 :
1633 : ! Write nonlop_ylm_counters (if enabled) in outputs
1634 5629 : if (dtset%useylm==1.and.dtset%nonlop_ylm_count/=0.and.dtset%paral_kgb==0) then
1635 0 : call nonlop_ylm_output_counters(dtset%natom,nbandtot,dtset%ntypat,dtset%typat,mpi_enreg)
1636 : end if
1637 : ! Write fft_counters (if enabled) in output
1638 5629 : if (dtset%fft_count/=0.and.dtset%paral_kgb==0) then
1639 0 : call fft_output_counters(nbandtot,mpi_enreg)
1640 : end if
1641 :
1642 5629 : if(dtset%imgwfstor==1)then
1643 2554229 : scf_history%cg(:,:,1)=cg(:,:)
1644 1749 : scf_history%eigen(:,1)=eigen(:)
1645 : endif
1646 :
1647 : !Deallocate arrays
1648 5629 : ABI_FREE(atindx)
1649 5629 : ABI_FREE(atindx1)
1650 5629 : ABI_FREE(indsym)
1651 5629 : ABI_FREE(npwarr)
1652 5629 : ABI_FREE(nattyp)
1653 5629 : ABI_FREE(resid)
1654 5629 : ABI_FREE(rhog)
1655 5629 : ABI_FREE(start)
1656 5629 : ABI_FREE(symrec)
1657 5629 : ABI_FREE(taug)
1658 5629 : ABI_FREE(ab_xfh%xfhist)
1659 5629 : call pawfgr_destroy(pawfgr)
1660 5629 : if (dtset%cprj_in_memory==1) then
1661 : !if (xg_nonlop%paw) then
1662 : ! call xg_nonlop_destroy_Sij(xg_nonlop)
1663 : !else
1664 : ! call xg_nonlop_destroy_ekb(xg_nonlop)
1665 : !end if
1666 162 : call xg_nonlop_destroy(xg_nonlop)
1667 : end if
1668 :
1669 5629 : if(dtset%imgwfstor==0)then
1670 5576 : if(dtset%gpu_option == ABI_GPU_KOKKOS) then
1671 : #if defined HAVE_GPU && defined HAVE_YAKL
1672 : ABI_FREE_MANAGED(cg)
1673 : #endif
1674 : else
1675 5576 : ABI_FREE(cg)
1676 : end if
1677 5576 : ABI_FREE(eigen)
1678 : else
1679 53 : nullify(cg,eigen)
1680 : endif
1681 :
1682 5629 : if (dtset%usewvl == 0) then
1683 : ! In wavelet case, irrzon and phnons are deallocated by wavelet object.
1684 5629 : ABI_FREE(irrzon)
1685 5629 : ABI_FREE(phnons)
1686 : end if
1687 :
1688 5629 : ABI_FREE(ylm)
1689 5629 : ABI_FREE(ylmgr)
1690 :
1691 5629 : if (scf_history%history_size<0) then
1692 4175 : if (psps%usepaw==1) then
1693 956 : call pawrhoij_free(pawrhoij)
1694 : end if
1695 4175 : ABI_FREE(rhor)
1696 4175 : ABI_FREE(taur)
1697 4175 : ABI_FREE(pawrhoij)
1698 4175 : ABI_FREE(xred_old)
1699 : else
1700 1454 : nullify(rhor,taur,pawrhoij,xred_old)
1701 : end if
1702 :
1703 : !PAW+DMFT
1704 5629 : call destroy_sc_dmft(paw_dmft)
1705 : ! This call should be done inside destroy_sc_dmft
1706 5629 : if ( dtset%usedmft /= 0 .and. dtset%dmft_entropy > 0) then
1707 3 : call data4entropyDMFT_destroy(paw_dmft%forentropyDMFT)
1708 : end if
1709 :
1710 : !Destroy extfpmd datastructure
1711 5629 : if(associated(extfpmd)) then
1712 7 : call extfpmd%destroy()
1713 7 : ABI_FREE(extfpmd)
1714 : end if
1715 :
1716 : !Destroy rcpaw datastructure
1717 5629 : if(associated(rcpaw)) then
1718 3 : call rcpaw_destroy(rcpaw)
1719 3 : ABI_FREE(rcpaw)
1720 : endif
1721 :
1722 : !Destroy electronpositron datastructure
1723 5629 : if (dtset%positron/=0) then
1724 14 : call destroy_electronpositron(electronpositron)
1725 : end if
1726 :
1727 : !Deallocating the basis set.
1728 5629 : if (dtset%usewvl == 1) then
1729 0 : call wvl_projectors_free(wvl%projectors)
1730 0 : call wvl_wfs_free(wvl%wfs)
1731 0 : call wvl_descr_free(wvl%descr)
1732 0 : call wvl_denspot_free(wvl%den)
1733 0 : if(dtset%usepaw == 1) then
1734 0 : call wvl_paw_free(wvl%descr)
1735 : end if
1736 : end if
1737 :
1738 5629 : ABI_FREE(kg)
1739 :
1740 5629 : if (dtset%icoulomb /= 0) then
1741 : call psolver_kernel((/ 0._dp, 0._dp, 0._dp /), 0, dtset%icoulomb, 0, kernel_dummy, &
1742 0 : & 0, dtset%ngfft, 1, dtset%nscforder)
1743 : end if
1744 :
1745 5629 : if (associated(pwind)) then
1746 5629 : ABI_FREE(pwind)
1747 : end if
1748 5629 : if (associated(pwnsfac)) then
1749 5629 : ABI_FREE(pwnsfac)
1750 : end if
1751 5629 : if ((dtset%berryopt<0).or.&
1752 : & (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
1753 : & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17)) then
1754 : if (xmpi_paral == 1) then
1755 82 : ABI_FREE(mpi_enreg%kptdstrb)
1756 : if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
1757 82 : & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then
1758 25 : ABI_FREE(mpi_enreg%kpt_loc2ibz_sp)
1759 : end if
1760 : end if
1761 82 : ABI_SFREE(mpi_enreg%kpt_loc2ibz_sp)
1762 82 : ABI_SFREE(mpi_enreg%kpt_loc2fbz_sp)
1763 82 : ABI_SFREE(mpi_enreg%mkmem)
1764 : end if
1765 : ! deallocate cprj
1766 5629 : if(usecprj==1) then
1767 59 : ABI_FREE(dimcprj_srt)
1768 59 : call pawcprj_free(cprj)
1769 : end if
1770 16589 : ABI_FREE(cprj)
1771 :
1772 : ! deallocate efield
1773 5629 : call dtefield%free()
1774 :
1775 5629 : if (dtset%userec == 1) call CleanRec(rec_set)
1776 :
1777 5629 : call hdr%free()
1778 5629 : call ebands%free()
1779 :
1780 5629 : if (me == master .and. dtset%prtxml == 1) then
1781 : ! The dataset given in argument has been treated, then we output its variables.
1782 : ! call outvarsXML()
1783 : ! gstate() will handle a dataset, so we output the dataSet markup.
1784 1 : write(ab_xml_out, "(A)") ' </dataSet>'
1785 : end if
1786 :
1787 5629 : if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2 .and. dtset%optdriver /= RUNL_GWLS) then
1788 : ! Plane-wave case
1789 5615 : call bandfft_kpt_destroy_array(bandfft_kpt,mpi_enreg)
1790 : end if
1791 :
1792 : !Clean gemm_nonlop work spaces
1793 5629 : if(gemm_nonlop_use_gemm) then
1794 29 : call destroy_gemm_nonlop(dtset%gpu_option)
1795 29 : gemm_nonlop_use_gemm = .false.
1796 : end if
1797 :
1798 : !Clean GPU work spaces
1799 : #if defined HAVE_GPU
1800 : if (dtset%gpu_option/=ABI_GPU_DISABLED) then
1801 : call dealloc_hamilt_gpu(2,dtset%gpu_option)
1802 : end if
1803 : #endif
1804 :
1805 : #if defined HAVE_BIGDFT
1806 : if (dtset%usewvl == 1 .and. dtset%timopt==10) then
1807 : call wvl_timing(xmpi_world,'== WFN OPT','PR')
1808 : end if
1809 : #endif
1810 :
1811 5629 : call timab(1231,2,tsec)
1812 5629 : call timab(1232,2,tsec)
1813 :
1814 : DBG_EXIT("COLL")
1815 :
1816 28145 : end subroutine gstate
1817 : !!***
1818 :
1819 : !!****f* m_gstate/setup2
1820 : !!
1821 : !! NAME
1822 : !! setup2
1823 : !!
1824 : !! FUNCTION
1825 : !! Call within main routine for setup of various arrays.
1826 : !!
1827 : !! INPUTS
1828 : !! dtset <type(dataset_type)>=all input variables in this dataset
1829 : !! | ecut=kinetic energy cutoff for planewave basis (hartree)
1830 : !! | natom=number of atoms in unit cell
1831 : !! | nkpt=number of k points
1832 : !! | wtk(nkpt)=integration weight associated with each k point
1833 : !! | iscf=parameter controlling scf or non-scf choice
1834 : !! npwtot(nkpt)=number of planewaves in basis and boundary at each k point
1835 : !! xred(3,natom)=starting reduced atomic coordinates
1836 : !!
1837 : !! OUTPUT
1838 : !! start(3,natom)=copy of starting xred
1839 : !!
1840 : !! SOURCE
1841 :
1842 5629 : subroutine setup2(dtset,npwtot,start,wfs,xred)
1843 :
1844 : !Arguments ------------------------------------
1845 : !scalars
1846 : type(dataset_type),intent(in) :: dtset
1847 : type(wvl_wf_type),intent(in) :: wfs
1848 : !arrays
1849 : integer,intent(in) :: npwtot(dtset%nkpt)
1850 : real(dp),intent(in) :: xred(3,dtset%natom)
1851 : real(dp),intent(out) :: start(3,dtset%natom)
1852 :
1853 : !Local variables-------------------------------
1854 : !scalars
1855 : integer :: ikpt,npw
1856 : real(dp) :: arith,geom,wtknrm
1857 : character(len=500) :: msg
1858 : ! *************************************************************************
1859 :
1860 5629 : if (dtset%iscf>=0) then
1861 :
1862 : ! Copy coordinates into array start
1863 53898 : start(:,:)=xred(:,:)
1864 :
1865 5118 : if (dtset%usewvl == 0) then
1866 : ! Get average number of planewaves per k point:
1867 : ! both arithmetic and GEOMETRIC averages are desired--
1868 : ! need geometric average to use method of Francis and Payne,
1869 : ! J. Phys.: Condens. Matter 2, 4395-4404 (1990) [[cite:Francis1990]].
1870 : ! Also note: force k point wts to sum to 1 for this averaging.
1871 : ! (wtk is not forced to add to 1 in a case with occopt=2)
1872 5118 : arith=zero
1873 5118 : geom=one
1874 5118 : wtknrm=zero
1875 57650 : do ikpt=1,dtset%nkpt
1876 52532 : npw=npwtot(ikpt)
1877 52532 : wtknrm=wtknrm+dtset%wtk(ikpt)
1878 52532 : arith=arith+npw*dtset%wtk(ikpt)
1879 57650 : geom=geom*npw**dtset%wtk(ikpt)
1880 : end do
1881 :
1882 : ! Enforce normalization of weights to 1
1883 5118 : arith=arith/wtknrm
1884 5118 : geom=geom**(1.0_dp/wtknrm)
1885 :
1886 : end if
1887 :
1888 : ! Ensure portability of output thanks to tol8
1889 5118 : if (dtset%usewvl == 0) then
1890 5118 : write(msg, '(a,2f12.3)' ) '_setup2: Arith. and geom. avg. npw (full set) are',arith+tol8,geom
1891 : else
1892 : #if defined HAVE_BIGDFT
1893 : write(msg, '(a,2I8)' ) ' setup2: nwvl coarse and fine are', &
1894 : & wfs%ks%lzd%Glr%wfd%nvctr_c, wfs%ks%lzd%Glr%wfd%nvctr_f
1895 : #endif
1896 : end if
1897 15354 : call wrtout([std_out, ab_out], msg)
1898 : end if
1899 :
1900 : #if !defined HAVE_BIGDFT
1901 : if (.false.) write(std_out,*) wfs%ks
1902 : #endif
1903 :
1904 5629 : end subroutine setup2
1905 : !!***
1906 :
1907 : !!****f* m_gstate/clnup1
1908 : !! NAME
1909 : !! clnup1
1910 : !!
1911 : !! FUNCTION
1912 : !! Perform "cleanup" at end of execution of gstate routine.
1913 : !!
1914 : !! INPUTS
1915 : !! acell(3)=length scales of primitive translations (bohr)
1916 : !! dosdeltae=DOS delta of Energy
1917 : !! dtset <type(dataset_type)>=all input variables in this dataset
1918 : !! eigen(mband*nkpt*nsppol)=eigenvalues (hartree) for all bands
1919 : !! at each k point
1920 : !! enunit=choice for units of output eigenvalues: 0=>hartree,
1921 : !! 1=> eV, 2=> hartree and eV
1922 : !! fermie=fermi energy (Hartree)
1923 : !! fermih=fermi energy for holes (Hartree) (for occopt 9)
1924 : !! fnameabo_dos=filename of output DOS file
1925 : !! fnameabo_eig=filename of output EIG file
1926 : !! gred(3,natom)=d(E)/d(xred) (hartree)
1927 : !! iatfix(3,natom)=0 if not fixed along specified direction,
1928 : !! 1 if fixed
1929 : !! iscf=parameter controlling scf or non-scf choice
1930 : !! kptopt=option for the generation of k points
1931 : !! kptns(3,nkpt)=k points in terms of recip primitive translations
1932 : !! mband=maximum number of bands
1933 : !! mpi_enreg=information about MPI parallelization
1934 : !! natom=number of atoms in unit cell
1935 : !! nband(nkpt*nsppol)=number of bands
1936 : !! nfft=(effective) number of FFT grid points (for this processor)
1937 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1938 : !! nkpt=number of k points
1939 : !! nspden=number of spin-density components
1940 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1941 : !! nstep=desired number of electron iteration steps
1942 : !! occ(maxval(nband(:))*nkpt*nsppol)=occupancies for each band and k point
1943 : !! occopt=option for occupancies
1944 : !! prtdos= if == 1, will print the density of states
1945 : !! prtfor= if >0, will print the forces
1946 : !! prtstm= input variable prtstm
1947 : !! prtvol=control print volume and debugging
1948 : !! resid(mband*nkpt*nsppol)=squared residuals for each band and k point where
1949 : !! resid(n,k)=|<C(n,k)|(H-e(n,k))|C(n,k)>|^2
1950 : !! rhor(nfft,nspden)=electron density (electrons/bohr^3)
1951 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1952 : !! tphysel="physical" electronic temperature with FD occupations
1953 : !! tsmear=smearing energy or temperature (if metal)
1954 : !! vxcavg=average of vxc potential
1955 : !! wtk(nkpt)=real(dp) array of k-point weights
1956 : !! xred(3,natom)=reduced atomic coordinates
1957 : !!
1958 : !! OUTPUT
1959 : !! (only print and write to disk)
1960 : !!
1961 : !! SOURCE
1962 5629 : subroutine clnup1(acell,dtset,eigen,fermie,fermih, fnameabo_dos,fnameabo_eig,gred,&
1963 5629 : mpi_enreg,nfft,ngfft,occ,prtfor, resid,rhor,rprimd,vxcavg,xred,rcpaw)
1964 :
1965 : !Arguments ------------------------------------
1966 : !scalars
1967 : integer,intent(in) :: nfft, prtfor
1968 : real(dp),intent(in) :: fermie, fermih, vxcavg
1969 : character(len=*),intent(in) :: fnameabo_dos,fnameabo_eig
1970 : type(dataset_type),intent(in) :: dtset
1971 : type(MPI_type),intent(in) :: mpi_enreg
1972 : type(rcpaw_type),pointer,intent(inout) :: rcpaw
1973 : !arrays
1974 : integer,intent(in) :: ngfft(18)
1975 : real(dp),intent(in) :: acell(3)
1976 : real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
1977 : real(dp),intent(in) :: gred(3,dtset%natom)
1978 : real(dp),intent(in) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
1979 : real(dp),intent(in) :: rhor(nfft,dtset%nspden)
1980 : real(dp),intent(in) :: rprimd(3,3)
1981 : real(dp),intent(in) :: xred(3,dtset%natom)
1982 : real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
1983 :
1984 : !Local variables-------------------------------
1985 : !scalars
1986 : integer,parameter :: master=0
1987 : integer :: comm,iatom,ii,iscf_dum,iwfrc,me,nnonsc,option,unitdos
1988 : real(dp) :: entropy,grmax,grsum,maxocc,nelect,tolwf,ucvol
1989 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1990 : character(len=500) :: msg
1991 : character(len=fnlen) filename
1992 : !arrays
1993 5629 : real(dp),allocatable :: doccde(:)
1994 :
1995 : ! ****************************************************************
1996 :
1997 5629 : comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm)
1998 :
1999 5629 : if(dtset%prtstm==0)then ! Write reduced coordinates xred
2000 5626 : write(msg, '(a,i5,a)' )' reduced coordinates (array xred) for',dtset%natom,' atoms'
2001 5626 : call wrtout(ab_out, msg)
2002 18847 : do iatom=1,dtset%natom
2003 13221 : write(msg, '(1x,3f20.12)' ) xred(:,iatom)
2004 18847 : call wrtout(ab_out, msg)
2005 : end do
2006 : end if
2007 :
2008 : !Write reduced gradients if iscf > 0 and dtset%nstep>0 and
2009 5629 : if (dtset%iscf>=0.and.dtset%nstep>0.and.dtset%prtstm==0) then
2010 :
2011 : ! Compute absolute maximum and root mean square value of gradients
2012 5095 : grmax=0.0_dp
2013 5095 : grsum=0.0_dp
2014 17215 : do iatom=1,dtset%natom
2015 53575 : do ii=1,3
2016 : ! To be activated in v5.5
2017 : ! grmax=max(grmax,abs(gred(ii,iatom)))
2018 36360 : grmax=max(grmax,gred(ii,iatom))
2019 48480 : grsum=grsum+gred(ii,iatom)**2
2020 : end do
2021 : end do
2022 5095 : grsum=sqrt(grsum/dble(3*dtset%natom))
2023 :
2024 5095 : write(msg, '(1x,a,1p,e12.4,a,e12.4,a)' )'rms dE/dt=',grsum,'; max dE/dt=',grmax,'; dE/dt below (all hartree)'
2025 5095 : call wrtout(ab_out, msg)
2026 17215 : do iatom=1,dtset%natom
2027 12120 : write(msg, '(i5,1x,3f20.12)' ) iatom,gred(1:3,iatom)
2028 17215 : call wrtout(ab_out, msg)
2029 : end do
2030 :
2031 : end if
2032 :
2033 5629 : if(dtset%prtstm==0)then
2034 :
2035 : ! Compute and write out dimensional cartesian coords and forces:
2036 5626 : call wrtout(ab_out,' ')
2037 :
2038 : ! (only write forces if iscf > 0 and dtset%nstep>0)
2039 5626 : if (dtset%iscf<0.or.dtset%nstep<=0.or.prtfor==0) then
2040 634 : iwfrc=0
2041 : else
2042 4992 : iwfrc=1
2043 : end if
2044 :
2045 5626 : call prtxf(gred,dtset%iatfix,ab_out,iwfrc,dtset%natom,rprimd,xred)
2046 :
2047 : ! Write length scales
2048 5626 : write(msg, '(1x,a,3f16.12,a)' )'length scales=',acell,' bohr'
2049 5626 : call wrtout(ab_out, msg)
2050 22504 : write(msg, '(14x,a,3f16.12,a)' )'=',Bohr_Ang*acell(1:3),' angstroms'
2051 5626 : call wrtout(ab_out, msg)
2052 :
2053 : end if
2054 :
2055 5629 : option=1; nnonsc=0; tolwf=0.0_dp
2056 :
2057 5629 : if(dtset%iscf<0 .and. dtset%iscf/=-3)option=3
2058 5629 : iscf_dum=dtset%iscf
2059 5629 : if(dtset%nstep==0)iscf_dum=-1
2060 :
2061 5629 : if(dtset%tfkinfunc==0)then
2062 5619 : if (me == master) then
2063 : call prteigrs(eigen,dtset%enunit,fermie,fermih,fnameabo_eig,ab_out,&
2064 : & iscf_dum,dtset%kptns,dtset%kptopt,dtset%mband,&
2065 : & dtset%nband,dtset%nbdbuf,dtset%nkpt,nnonsc,dtset%nsppol,occ,&
2066 : & dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwf,&
2067 4371 : & vxcavg,dtset%wtk)
2068 : call prteigrs(eigen,dtset%enunit,fermie,fermih,fnameabo_eig,std_out,&
2069 : & iscf_dum,dtset%kptns,dtset%kptopt,dtset%mband,&
2070 : & dtset%nband,dtset%nbdbuf,dtset%nkpt,nnonsc,dtset%nsppol,occ,&
2071 : & dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwf,&
2072 4371 : & vxcavg,dtset%wtk)
2073 : end if
2074 :
2075 5619 : if (dtset%prteig==1 .and. me == master) then
2076 2694 : filename=trim(fnameabo_eig)//'.nc'
2077 2694 : call write_eig(eigen,fermie,filename,dtset%kptns,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol)
2078 : end if
2079 : end if
2080 :
2081 : !Compute and print location of maximal and minimal density
2082 5629 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2083 5629 : call prtrhomxmn(std_out,mpi_enreg,nfft,ngfft,dtset%nspden,2,rhor,ucvol=ucvol)
2084 5629 : if( dtset%prtvol>1)then
2085 347 : call prtrhomxmn(ab_out,mpi_enreg,nfft,ngfft,dtset%nspden,2,rhor,ucvol=ucvol)
2086 : end if
2087 :
2088 : !If needed, print DOS (unitdos is closed in getnel, occ is not changed if option == 2
2089 5629 : if ((dtset%prtdos==1.or.dtset%prtdos==4) .and. me == master) then
2090 9 : if (open_file(fnameabo_dos,msg, newunit=unitdos, status='unknown', action="write", form='formatted') /= 0) then
2091 0 : ABI_ERROR(msg)
2092 : end if
2093 9 : rewind(unitdos)
2094 9 : maxocc=two/(dtset%nspinor*dtset%nsppol) ! Will not work in the fixed moment case
2095 9 : option=2
2096 27 : ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
2097 : call getnel(doccde,dtset%dosdeltae,eigen,entropy,fermie,fermih,&
2098 : & maxocc,dtset%mband,dtset%nband,nelect,dtset%nkpt,&
2099 : & dtset%nsppol,occ,dtset%occopt,option,dtset%tphysel,&
2100 9 : & dtset%tsmear,unitdos,dtset%wtk,1,dtset%nband(1),rcpaw=rcpaw)!CP: added 1, nband(1) to fit new definition of getnel; parameters only used if
2101 9 : ABI_FREE(doccde)
2102 : end if
2103 :
2104 5629 : end subroutine clnup1
2105 : !!***
2106 :
2107 : !!****f* m_gstate/prtxf
2108 : !! NAME
2109 : !! prtxf
2110 : !!
2111 : !! FUNCTION
2112 : !! Compute and print out dimensional cartesian coordinates and forces.
2113 : !! Note: for x=cartesian coordinates, t=reduced coordinates (xred),
2114 : !! =>
2115 : !! $ x= R t $
2116 : !! =>
2117 : !! $ x(1)=rprimd(1,1) t(1)+rprimd(2,1) t(2)+rprimd(3,1) t(3)$
2118 : !! etc. Also $ t = (R^{-1}) x$ .
2119 : !! To convert gradients, $d(E)/dx(n) = [d(E)/dt(m)] [dt(m)/dx(n)]$
2120 : !! and $ dt(m)/dx(n) = (R^{-1})_{mn} = G_{nm}$ because G is the
2121 : !! inverse transpose of R. Finally then
2122 : !! $d(E)/dx(n) = G_{nm} [d(E)/dt(m)]$.
2123 : !! The vector $d(E)/dt(m)$ for each atom is input in gred
2124 : !! (grad. wrt xred).
2125 : !!
2126 : !! INPUTS
2127 : !! gred(3,natom)=gradients of Etot (hartree) wrt xred(3,natom)
2128 : !! iatfix(3,natom)=1 for each fixed atom along specified
2129 : !! direction, else 0
2130 : !! iout=unit number for output file
2131 : !! iwfrc=controls force output: 0=> no forces output,
2132 : !! 1=>forces out in eV/A and Ha/bohr,
2133 : !! 2=>forces out in Ha/bohr
2134 : !! natom=number of atoms in unit cell
2135 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
2136 : !! xred(3,natom)=relative coordinates of atoms (in terms of prim. transl.)
2137 : !!
2138 : !! OUTPUT
2139 : !! (data written to unit iout)
2140 : !!
2141 : !! SOURCE
2142 :
2143 5626 : subroutine prtxf(gred,iatfix,iout,iwfrc,natom,rprimd,xred)
2144 :
2145 : !Arguments ------------------------------------
2146 : !scalars
2147 : integer,intent(in) :: iout,iwfrc,natom
2148 : !arrays
2149 : integer,intent(in) :: iatfix(3,natom)
2150 : real(dp),intent(in) :: gred(3,natom),rprimd(3,3),xred(3,natom)
2151 :
2152 : !Local variables-------------------------------
2153 : !scalars
2154 : integer :: iatom,mu,unfixd
2155 : real(dp) :: convt,fmax,frms
2156 : character(len=15) :: format_line21
2157 : character(len=15) :: format_line25
2158 : character(len=15) :: format_line
2159 : character(len=500) :: msg
2160 : !arrays
2161 : real(dp) :: favg(3),favg_out(3),ff(3),gprimd(3,3),xx(3)
2162 :
2163 : ! ****************************************************************
2164 :
2165 5626 : format_line21='(i5,1x,3f21.14)'
2166 5626 : format_line25='(i5,1x,3f25.14)'
2167 :
2168 : !Write cartesian coordinates in angstroms
2169 5626 : call wrtout(iout,' cartesian coordinates (angstrom) at end:')
2170 18847 : do iatom=1,natom
2171 13221 : format_line=format_line21
2172 52884 : do mu=1,3
2173 : xx(mu)=(rprimd(mu,1)*xred(1,iatom)+&
2174 : & rprimd(mu,2)*xred(2,iatom)+&
2175 39663 : & rprimd(mu,3)*xred(3,iatom))*Bohr_Ang
2176 52884 : if(xx(mu)>99999 .or. xx(mu)<-9999)format_line=format_line25
2177 : end do
2178 13221 : write(msg,format_line) iatom,xx
2179 18847 : call wrtout(iout, msg)
2180 : end do
2181 :
2182 : !Optionally write cartesian forces in eV/Angstrom (also provide same in hartree/bohr)
2183 5626 : if (iwfrc/=0) then
2184 : ! First, provide results in hartree/bohr
2185 4992 : write(msg, '(a,a)' ) ch10,' cartesian forces (hartree/bohr) at end:'
2186 4992 : call wrtout(iout, msg)
2187 4992 : frms=zero
2188 4992 : fmax=zero
2189 4992 : favg(1)=zero
2190 4992 : favg(2)=zero
2191 4992 : favg(3)=zero
2192 : ! To get cartesian forces from input gradients with respect to
2193 : ! dimensionless coordinates xred, multiply by G and negate
2194 : ! (see notes at top of this subroutine)
2195 4992 : call matr3inv(rprimd,gprimd)
2196 : ! First compute (spurious) average force favg
2197 16562 : do iatom=1,natom
2198 51272 : do mu=1,3
2199 : ff(mu)=-(gprimd(mu,1)*gred(1,iatom)+&
2200 : & gprimd(mu,2)*gred(2,iatom)+&
2201 34710 : & gprimd(mu,3)*gred(3,iatom))
2202 46280 : favg(mu)=favg(mu)+ff(mu)
2203 : end do
2204 : end do
2205 4992 : favg(1) = favg(1)/dble(natom)
2206 4992 : favg(2) = favg(2)/dble(natom)
2207 4992 : favg(3) = favg(3)/dble(natom)
2208 :
2209 : ! Subtract off average force in what follows
2210 : ! (avg is also subtracted off in carfor, called by loopcv,
2211 : ! called by grad)
2212 4992 : unfixd=0
2213 16562 : do iatom=1,natom
2214 11570 : format_line=format_line21
2215 46280 : do mu=1,3
2216 : ff(mu)=-(gprimd(mu,1)*gred(1,iatom)+&
2217 : & gprimd(mu,2)*gred(2,iatom)+&
2218 34710 : & gprimd(mu,3)*gred(3,iatom))-favg(mu)
2219 34710 : if(ff(mu)>99999 .or. ff(mu)<-9999)format_line=format_line25
2220 : ! For rms and max force, include only unfixed components
2221 46280 : if (iatfix(mu,iatom) /= 1) then
2222 30701 : unfixd=unfixd+1
2223 30701 : frms=frms+ff(mu)**2
2224 30701 : fmax=max(fmax,abs(ff(mu)))
2225 : end if
2226 : end do
2227 11570 : write(msg, format_line) iatom,ff
2228 16562 : call wrtout(iout, msg)
2229 : end do
2230 4992 : if ( unfixd /= 0 ) frms = sqrt(frms/dble(unfixd))
2231 :
2232 : ! The average force is obtained from the cancellation of numbers
2233 : ! of typical size unity, so an absolute value lower
2234 : ! than tol14 is meaningless for the output file.
2235 4992 : favg_out(:)=favg(:)
2236 4992 : if(abs(favg_out(1))<tol14)favg_out(1)=zero
2237 4992 : if(abs(favg_out(2))<tol14)favg_out(2)=zero
2238 4992 : if(abs(favg_out(3))<tol14)favg_out(3)=zero
2239 :
2240 4992 : write(msg, '(a,1p,2e14.7,1x,3e11.3,a)' )' frms,max,avg=',frms,fmax,favg_out(1:3),' h/b'
2241 4992 : call wrtout(iout, msg)
2242 :
2243 4992 : if (iwfrc==1) then
2244 :
2245 4992 : write(msg, '(a,a)' )ch10,' cartesian forces (eV/Angstrom) at end:'
2246 4992 : call wrtout(iout, msg)
2247 4992 : convt=Ha_eV/Bohr_Ang
2248 :
2249 : ! Note: subtract off average force
2250 16562 : do iatom=1,natom
2251 11570 : format_line=format_line21
2252 46280 : do mu=1,3
2253 : ff(mu)=(-(gprimd(mu,1)*gred(1,iatom)+&
2254 : & gprimd(mu,2)*gred(2,iatom)+&
2255 34710 : & gprimd(mu,3)*gred(3,iatom))-favg(mu))*convt
2256 46280 : if(ff(mu)>99999 .or. ff(mu)<-9999)format_line=format_line25
2257 : end do
2258 11570 : write(msg, format_line) iatom,ff
2259 16562 : call wrtout(iout, msg)
2260 : end do
2261 24960 : write(msg, '(a,1p,2e14.7,1x,3e11.3,a)' )' frms,max,avg=',convt*frms,convt*fmax,convt*favg_out(1:3),' e/A'
2262 4992 : call wrtout(iout, msg)
2263 :
2264 : end if
2265 : end if
2266 :
2267 5626 : end subroutine prtxf
2268 : !!***
2269 :
2270 : !!****f* m_gstate/clnup2
2271 : !! NAME
2272 : !! clnup2
2273 : !!
2274 : !! FUNCTION
2275 : !! Perform more "cleanup" after completion of iterations.
2276 : !! This subroutine prints out more breakdown of force
2277 : !! information, shifts of atomic positions, and stresses.
2278 : !!
2279 : !! INPUTS
2280 : !! gred(3,natom)=d(E_total)/d(xred) derivatives (hartree)
2281 : !! grchempottn(3,natom)=d(E_chempot)/d(xred) derivatives (hartree)
2282 : !! grewtn(3,natom)=d(E_Ewald)/d(xred) derivatives (hartree)
2283 : !! grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D2 dispersion (hartree)
2284 : !! grxc(3,natom)=d(Exc)/d(xred) derivatives (0 without core charges)
2285 : !! iscf=parameter controlling scf or non-scf iterations
2286 : !! natom=number of atoms in unit cell
2287 : !! ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
2288 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
2289 : !! prtfor= >0 if forces have to be printed (0 otherwise)
2290 : !! prtstr= >0 if stresses have to be printed (0 otherwise)
2291 : !! prtvol=control print volume and debugging output
2292 : !! start(3,natom)=starting coordinates in terms of real space
2293 : !! primitive translations
2294 : !! strten(6)=components of the stress tensor (hartree/bohr^3)
2295 : !! synlgr(3,natom)=d(E_nlpsp)/d(xred) derivatives (hartree)
2296 : !! xred(3,natom)=final coordinates in terms of primitive translations
2297 : !!
2298 : !! OUTPUT
2299 : !! (only print)
2300 : !!
2301 : !! SOURCE
2302 :
2303 5586 : subroutine clnup2(n1xccc,gred,grchempottn,gresid,grewtn,grvdw,grxc,iscf,natom,ngrvdw,&
2304 5586 : & prtfor,prtstr,prtvol,start,strten,synlgr,xred)
2305 :
2306 : !Arguments ------------------------------------
2307 : !scalars
2308 : integer,intent(in) :: iscf,n1xccc,natom,ngrvdw,prtfor,prtstr,prtvol
2309 : !arrays
2310 : real(dp),intent(in) :: gred(3,natom),grchempottn(3,natom),gresid(3,natom)
2311 : real(dp),intent(in) :: grewtn(3,natom),grvdw(3,ngrvdw)
2312 : real(dp),intent(in) :: grxc(3,natom),start(3,natom),strten(6),synlgr(3,natom)
2313 : real(dp),intent(in) :: xred(3,natom)
2314 :
2315 : !Local variables-------------------------------
2316 : character(len=*), parameter :: format01020 ="(i5,1x,3f20.12)"
2317 : !scalars
2318 : integer :: iatom,mu
2319 : real(dp) :: devsqr,grchempot2
2320 : character(len=500) :: msg
2321 : integer :: units(2)
2322 : ! *************************************************************************
2323 :
2324 : !write(std_out,*)' clnup2 : enter '
2325 :
2326 : !Only print additional info for scf calculations
2327 5586 : if (iscf>=0) then
2328 :
2329 5090 : if(prtvol >= 10 .and. prtfor > 0) then
2330 188 : write(msg, '(a,10x,a)' ) ch10, '===> extra information on forces <==='
2331 188 : call wrtout(ab_out,msg)
2332 :
2333 188 : call wrtout(ab_out, ' ewald contribution to reduced grads')
2334 993 : do iatom=1,natom
2335 805 : write(msg,format01020) iatom,(grewtn(mu,iatom),mu=1,3)
2336 993 : call wrtout(ab_out,msg)
2337 : end do
2338 :
2339 3408 : grchempot2=sum(grchempottn(:,:)**2)
2340 188 : if(grchempot2>tol16)then
2341 6 : call wrtout(ab_out, ' chemical potential contribution to reduced grads')
2342 18 : do iatom=1,natom
2343 12 : write(msg,format01020) iatom,(grchempottn(mu,iatom),mu=1,3)
2344 18 : call wrtout(ab_out,msg)
2345 : end do
2346 : end if
2347 :
2348 188 : call wrtout(ab_out,' nonlocal contribution to red. grads')
2349 993 : do iatom=1,natom
2350 805 : write(msg,format01020) iatom,(synlgr(mu,iatom),mu=1,3)
2351 993 : call wrtout(ab_out,msg)
2352 : end do
2353 :
2354 188 : call wrtout(ab_out, ' local psp contribution to red. grads')
2355 188 : if (n1xccc /= 0) then
2356 641 : do iatom=1,natom
2357 567 : write(msg,format01020) iatom,gred(:,iatom) - &
2358 2835 : (grewtn(:,iatom)+grchempottn(:,iatom)+synlgr(:,iatom)+grxc(:,iatom)+gresid(:,iatom))
2359 641 : call wrtout(ab_out,msg)
2360 : end do
2361 : else
2362 352 : do iatom=1,natom
2363 238 : write(msg,format01020) iatom,gred(:,iatom) - &
2364 1190 : (grewtn(:,iatom)+grchempottn(:,iatom)+synlgr(:,iatom)+gresid(:,iatom))
2365 352 : call wrtout(ab_out,msg)
2366 : end do
2367 : end if
2368 :
2369 188 : if (n1xccc /= 0) then
2370 74 : call wrtout(ab_out,' core charge xc contribution to reduced grads')
2371 641 : do iatom=1,natom
2372 567 : write(msg,format01020) iatom,(grxc(mu,iatom),mu=1,3)
2373 641 : call wrtout(ab_out,msg)
2374 : end do
2375 : end if
2376 :
2377 188 : if (ngrvdw == natom) then
2378 0 : call wrtout(ab_out,' Van der Waals DFT-D contribution to reduced grads')
2379 0 : do iatom=1,natom
2380 0 : write(msg,format01020) iatom,(grvdw(mu,iatom),mu=1,3)
2381 0 : call wrtout(ab_out,msg)
2382 : end do
2383 : end if
2384 :
2385 188 : call wrtout(ab_out,' residual contribution to red. grads')
2386 993 : do iatom=1,natom
2387 805 : write(msg,format01020) iatom,(gresid(mu,iatom),mu=1,3)
2388 993 : call wrtout(ab_out,msg)
2389 : end do
2390 :
2391 : end if
2392 :
2393 : ! Compute mean squared deviation from starting coords
2394 5090 : devsqr=zero
2395 17183 : do iatom=1,natom
2396 53462 : do mu=1,3
2397 48372 : devsqr=devsqr+(xred(mu,iatom)-start(mu,iatom))**2
2398 : end do
2399 : end do
2400 :
2401 : ! When shift is nonnegligible then print values
2402 5090 : if (devsqr>1.d-14) then
2403 286 : write(msg, '(a,1p,e12.4,3x,a)' )' rms coord change=',sqrt(devsqr/dble(3*natom)),'atom, delta coord (reduced):'
2404 286 : call wrtout(ab_out,msg)
2405 1198 : do iatom=1,natom
2406 3648 : write(msg, '(1x,i5,2x,3f20.12)' ) iatom, (xred(mu,iatom)-start(mu,iatom),mu=1,3)
2407 1198 : call wrtout(ab_out,msg)
2408 : end do
2409 : end if
2410 :
2411 5090 : if (prtstr > 0) then
2412 : ! Write out stress results
2413 14997 : units = [std_out, ab_out]
2414 4999 : write(msg, '(a,a)' ) ch10,' Cartesian components of stress tensor (hartree/bohr^3)'
2415 4999 : call wrtout(units, msg)
2416 :
2417 4999 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4)
2418 4999 : call wrtout(units, msg)
2419 4999 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5)
2420 4999 : call wrtout(units, msg)
2421 4999 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) ' sigma(3 3)=',strten(3),' sigma(2 1)=',strten(6)
2422 4999 : call wrtout(units, msg)
2423 :
2424 : ! Also output the pressure (minus one third the trace of the stress tensor).
2425 4999 : write(msg, '(a,a,es12.4,a)' ) ch10,&
2426 4999 : '-Cartesian components of stress tensor (GPa) [Pressure=',&
2427 9998 : -(strten(1)+strten(2)+strten(3))*HaBohr3_GPa/3.0_dp,' GPa]'
2428 4999 : call wrtout(units, msg)
2429 :
2430 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2431 4999 : '- sigma(1 1)=',strten(1)*HaBohr3_GPa,&
2432 9998 : ' sigma(3 2)=',strten(4)*HaBohr3_GPa
2433 4999 : call wrtout(units, msg)
2434 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2435 4999 : '- sigma(2 2)=',strten(2)*HaBohr3_GPa,&
2436 9998 : ' sigma(3 1)=',strten(5)*HaBohr3_GPa
2437 4999 : call wrtout(units, msg)
2438 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2439 4999 : '- sigma(3 3)=',strten(3)*HaBohr3_GPa,&
2440 9998 : ' sigma(2 1)=',strten(6)*HaBohr3_GPa
2441 4999 : call wrtout(units, msg)
2442 : end if
2443 :
2444 : end if ! iscf > 0
2445 :
2446 : !write(std_out,*)' clnup2 : exit '
2447 :
2448 5586 : end subroutine clnup2
2449 : !!***
2450 :
2451 : !!****f* m_gstate/pawuj_drive
2452 : !! NAME
2453 : !! pawuj_drive
2454 : !!
2455 : !! FUNCTION
2456 : !! Drive for automatic determination of U
2457 : !! Relevant only in PAW+U context
2458 : !!
2459 : !! INPUTS
2460 : !! atindx(natom)=index table for atoms (see gstate.f)
2461 : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
2462 : !! cpus= cpu time limit in seconds
2463 : !! dtfil <type(datafiles_type)>=variables related to files
2464 : !! dtset <type(dataset_type)>=all input variables for this dataset
2465 : !! | mband=maximum number of bands
2466 : !! | mgfft=maximum size of 1D FFTs for the "coarse" grid (see NOTES below)
2467 : !! | mkmem =number of k points treated by this node.
2468 : !! | mpw=maximum dimensioned size of npw.
2469 : !! | natom=number of atoms in cell.
2470 : !! | nfft=(effective) number of FFT grid points (for this processor)
2471 : !! | for the "coarse" grid (see NOTES below)
2472 : !! | nkpt=number of k points
2473 : !! | nspden=number of spin-density components
2474 : !! | nsppol=1 for unpolarized, 2 for spin-polarized
2475 : !! | nsym=number of symmetry elements in space group
2476 : !! ecore=core psp energy (part of total energy) (hartree)
2477 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
2478 : !! mpi_enreg=information about MPI parallelization
2479 : !! nattyp(ntypat)= # atoms of each type.
2480 : !! npwarr(nkpt)=number of planewaves in basis at this k point
2481 : !! nspinor=number of spinorial components of the wavefunctions
2482 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
2483 : !! | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
2484 : !! pwind(pwind_alloc,2,3) = array used to compute
2485 : !! the overlap matrix smat between k-points (see initberry.f)
2486 : !! pwind_alloc = first dimension of pwind
2487 : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
2488 : !! (see initberry.f)
2489 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
2490 : !! ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
2491 : !!
2492 : !! OUTPUT
2493 : !! resid(mband*nkpt*nsppol)=residuals for each band over all k points and spins
2494 : !!
2495 : !! SIDE EFFECTS
2496 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)=updated wavefunctions.
2497 : !! dtefield <type(efield_type)> = variables related to Berry phase calculations (see initberry.f)
2498 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
2499 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
2500 : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
2501 : !! initialized= if 0 the initialization of the gstate run is not yet finished
2502 : !! irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
2503 : !! nfftf=(effective) number of FFT grid points (for this processor)
2504 : !! for the "fine" grid (see NOTES below)
2505 : !! occ(mband*nkpt*nsppol)=occupation number for each band (often 2) at each k point
2506 : !! pawrhoij(natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
2507 : !! phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
2508 : !! results_gs <type(results_gs_type)>=results (energy and its components,
2509 : !! forces and its components, the stress tensor) of a ground-state computation
2510 : !! (should be made a pure output quantity)
2511 : !! rhog(2,nfftf)=array for Fourier transform of electron density
2512 : !! rhor(nfftf,nspden)=array for electron density in el./bohr**3
2513 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2514 : !! scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles
2515 : !! symrec(3,3,nsym)=symmetry operations in reciprocal space
2516 : !! taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density
2517 : !! taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
2518 : !! wffnew,wffnow=struct info for wf disk files.
2519 : !! wvl <type(wvl_data)>=all wavelets data.
2520 : !! xred(3,natom)=reduced dimensionless atomic coordinates
2521 : !! xred_old(3,natom)= at input, previous reduced dimensionless atomic coordinates
2522 : !! at output, current xred is transferred to xred_old
2523 : !!
2524 : !! SOURCE
2525 :
2526 10 : subroutine pawuj_drive(scfcv_args, dtset,electronpositron,rhog,rhor,rprimd, xred,xred_old)
2527 :
2528 : !Arguments ------------------------------------
2529 : !scalars
2530 : type(scfcv_t), intent(inout) :: scfcv_args
2531 : type(dataset_type),intent(inout) :: dtset
2532 : type(electronpositron_type),pointer :: electronpositron
2533 : !type(wffile_type),intent(inout) :: wffnew,wffnow
2534 : !arrays
2535 : real(dp), intent(inout) :: rprimd(3,3)
2536 : real(dp), pointer :: rhog(:,:),rhor(:,:)
2537 : real(dp), intent(inout) :: xred(3,dtset%natom),xred_old(3,dtset%natom)
2538 :
2539 : !Local variables -------------------------
2540 : !scalars
2541 : integer,parameter :: itime0 = 0
2542 : integer,target :: ndtpawuj=4
2543 : integer :: iuj,conv_retcode
2544 : integer :: itimes(2)
2545 : real(dp) :: ures
2546 : !character(len=500) :: msg
2547 : !arrays
2548 : !real(dp),allocatable :: cgstart(:,:)
2549 10 : type(macro_uj_type),allocatable,target :: dtpawuj(:)
2550 : ! *********************************************************************
2551 :
2552 :
2553 : DBG_ENTER("COLL")
2554 :
2555 10 : if (dtset%macro_uj==0) then
2556 0 : ABI_BUG('Macro_uj must be set !')
2557 : end if
2558 :
2559 80 : ABI_MALLOC(dtpawuj,(0:ndtpawuj))
2560 : !ABI_MALLOC(cgstart,(2,scfcv_args%mcg))
2561 :
2562 10 : call pawuj_ini(dtpawuj,ndtpawuj)
2563 :
2564 : !cgstart=scfcv_args%cg
2565 50 : do iuj=1,ndtpawuj
2566 : ! allocate(dtpawuj(iuj)%rprimd(3,3)) ! this has already been done in pawuj_ini
2567 40 : dtpawuj(iuj)%macro_uj=dtset%macro_uj
2568 40 : dtpawuj(iuj)%pawprtvol=dtset%pawprtvol
2569 40 : dtpawuj(iuj)%diemix=dtset%diemix
2570 40 : dtpawuj(iuj)%diemixmag=dtset%diemixmag
2571 40 : dtpawuj(iuj)%pawujat=dtset%pawujat
2572 40 : dtpawuj(iuj)%nspden=dtset%nspden
2573 560 : dtpawuj(iuj)%rprimd=dtset%rprimd_orig(1:3,1:3,1)
2574 50 : dtpawuj(iuj)%dmatpuopt=dtset%dmatpuopt
2575 : end do
2576 :
2577 10 : iuj=1 !LMac Flag to collect occupancies for unperturbed calculation
2578 10 : dtpawuj(iuj)%iuj=iuj
2579 :
2580 10 : scfcv_args%ndtpawuj=>ndtpawuj
2581 10 : scfcv_args%dtpawuj=>dtpawuj
2582 :
2583 10 : itimes(1)=itime0 ; itimes(2)=1
2584 10 : call scfcv_run(scfcv_args,electronpositron,itimes,rhog,rhor,rprimd,xred,xred_old,conv_retcode)
2585 :
2586 : !Calculate Hubbard U (or J)
2587 10 : call pawuj_det(dtpawuj, ndtpawuj, dtset, scfcv_args%dtfil, ures, scfcv_args%mpi_enreg%comm_cell)
2588 10 : dtset%upawu(dtset%typat(dtset%pawujat),1)=ures/Ha_eV
2589 :
2590 : !Deallocations
2591 60 : do iuj=0,ndtpawuj
2592 60 : call pawuj_free(dtpawuj(iuj))
2593 : end do
2594 :
2595 60 : ABI_FREE(dtpawuj)
2596 : !ABI_FREE(cgstart)
2597 :
2598 : DBG_EXIT("COLL")
2599 :
2600 10 : end subroutine pawuj_drive
2601 : !!***
2602 :
2603 : !!****f* ABINIT/outxfhist
2604 : !! NAME
2605 : !! outxfhist
2606 : !!
2607 : !! FUNCTION
2608 : !! read/write xfhist
2609 : !!
2610 : !! COPYRIGHT
2611 : !! Copyright (C) 2003-2026 ABINIT group (MB)
2612 : !! This file is distributed under the terms of the
2613 : !! GNU General Public License, see ~abinit/COPYING
2614 : !! or http://www.gnu.org/copyleft/gpl.txt .
2615 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
2616 : !!
2617 : !! INPUTS
2618 : !! option =
2619 : !! 1: write
2620 : !! 2: read only nxfh
2621 : !! 3: read xfhist
2622 : !! response =
2623 : !! 0: GS wavefunctions
2624 : !! 1: RF wavefunctions
2625 : !! natom = number of atoms in unit cell
2626 : !! mxfh = last dimension of the xfhist array
2627 : !!
2628 : !! OUTPUT
2629 : !! ios = error code returned by read operations
2630 : !!
2631 : !! SIDE EFFECTS
2632 : !! nxfh = actual number of (x,f) history pairs, see xfhist array
2633 : !! wff2 = structured info for wavefunctions
2634 : !! xfhist(3,natom+4,2,ab_xfh%mxfh) = (x,f) history array, also including
2635 : !! rprim and stress
2636 : !!
2637 : !! SOURCE
2638 :
2639 0 : subroutine outxfhist(ab_xfh,natom,option,wff2,ios)
2640 :
2641 : use defs_basis
2642 : use m_abicore
2643 : use m_abimover
2644 : use m_xmpi
2645 : use m_wffile
2646 : use m_errors
2647 : use netcdf
2648 :
2649 : !Arguments ------------------------------------
2650 : integer ,intent(in) :: natom,option
2651 : integer ,intent(out) :: ios
2652 : type(wffile_type),intent(inout) :: wff2
2653 : type(ab_xfh_type),intent(inout) :: ab_xfh
2654 :
2655 : !Local variables-------------------------------
2656 : integer :: ierr,ixfh,ncid_hdr,spaceComm,xfdim2
2657 0 : real(dp),allocatable :: xfhist_tmp(:)
2658 : character(len=500) :: msg
2659 : !no_abirules
2660 : integer :: ncerr
2661 : integer :: nxfh_id, mxfh_id, xfdim2_id, dim2inout_id, dimr3_id,xfhist_id
2662 : integer :: nxfh_tmp,mxfh_tmp,xfdim2_tmp,dim2inout_tmp
2663 : ! *************************************************************************
2664 :
2665 0 : ncid_hdr = wff2%unwff
2666 0 : xfdim2 = natom+4
2667 :
2668 0 : ios = 0
2669 :
2670 : !### (Option=1) Write out content of all iterations
2671 : !#####################################################################
2672 0 : if ( option == 1 ) then
2673 :
2674 : ! Write the (x,f) history
2675 0 : if (wff2%iomode == IO_MODE_FORTRAN) then
2676 0 : write(unit=wff2%unwff)ab_xfh%nxfh
2677 0 : do ixfh=1,ab_xfh%nxfh
2678 0 : write(unit=wff2%unwff)ab_xfh%xfhist(:,:,:,ixfh)
2679 : end do
2680 :
2681 0 : else if (wff2%iomode == IO_MODE_FORTRAN_MASTER) then
2682 : ! FIXME: should copy the xfhist to other processors, and check that we are on the master to read in this case
2683 : ! if node is master
2684 0 : write(msg, "(A,A,A,A)") ch10, " outxfhist: ERROR -", ch10, &
2685 0 : & 'iomode == -1 (localrdwf ) has not been coded yet for xfhist rereading.'
2686 0 : ABI_ERROR(msg)
2687 :
2688 0 : write(unit=wff2%unwff)ab_xfh%nxfh
2689 0 : do ixfh=1,ab_xfh%nxfh
2690 0 : write(unit=wff2%unwff)ab_xfh%xfhist(:,:,:,ixfh)
2691 : end do
2692 :
2693 : ! insert mpi broadcast here
2694 :
2695 0 : else if(wff2%iomode==IO_MODE_MPI)then
2696 0 : ABI_MALLOC(xfhist_tmp,(3*(natom+4)*2))
2697 0 : spaceComm=xmpi_comm_self
2698 0 : call xderiveWRecInit(wff2,ierr)
2699 0 : call xderiveWrite(wff2,ab_xfh%nxfh,ierr)
2700 0 : call xderiveWRecEnd(wff2,ierr)
2701 0 : do ixfh=1,ab_xfh%nxfh
2702 0 : xfhist_tmp(:)=reshape(ab_xfh%xfhist(:,:,:,ixfh),(/3*(natom+4)*2/))
2703 0 : call xderiveWRecInit(wff2,ierr)
2704 0 : call xderiveWrite(wff2,xfhist_tmp,3*(natom+4)*2,spaceComm,ierr)
2705 0 : call xderiveWRecEnd(wff2,ierr)
2706 : end do
2707 0 : ABI_FREE(xfhist_tmp)
2708 :
2709 0 : else if (wff2%iomode == IO_MODE_NETCDF) then
2710 : ! check if nxfh and xfhist are defined
2711 0 : ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id)
2712 :
2713 0 : if (ncerr /= NF90_NOERR) then
2714 : ! need to define everything
2715 0 : ncerr = nf90_redef (ncid=ncid_hdr)
2716 0 : NCF_CHECK_MSG(ncerr," outxfhist : going to define mode ")
2717 :
2718 0 : ncerr = nf90_def_dim(ncid=ncid_hdr,name="dim2inout",len=2,dimid=dim2inout_id)
2719 0 : NCF_CHECK_MSG(ncerr," outxfhist : define dim2inout")
2720 0 : ncerr = nf90_def_dim(ncid=ncid_hdr,name="mxfh",len=ab_xfh%mxfh,dimid=mxfh_id)
2721 0 : NCF_CHECK_MSG(ncerr," outxfhist : define mxfh")
2722 0 : ncerr = nf90_def_dim(ncid=ncid_hdr,name="nxfh",len=ab_xfh%nxfh,dimid=nxfh_id)
2723 0 : NCF_CHECK_MSG(ncerr," outxfhist : define nxfh")
2724 0 : ncerr = nf90_def_dim(ncid=ncid_hdr,name="xfdim2",len=xfdim2,dimid=xfdim2_id)
2725 0 : NCF_CHECK_MSG(ncerr," outxfhist : define xfdim2")
2726 :
2727 0 : ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="dimr3",dimid=dimr3_id)
2728 0 : NCF_CHECK_MSG(ncerr," outxfhist : inquire dimr3")
2729 :
2730 : ! ab_xfh%xfhist(3,natom+4,2,ab_xfh%mxfh)
2731 : ncerr = nf90_def_var(ncid=ncid_hdr,name="xfhist",xtype=NF90_DOUBLE,&
2732 0 : & dimids=(/dimr3_id,xfdim2_id,dim2inout_id,mxfh_id/),varid=xfhist_id)
2733 0 : NCF_CHECK_MSG(ncerr," outxfhist : define xfhist")
2734 :
2735 : ! End define mode and go to data mode
2736 0 : ncerr = nf90_enddef(ncid=ncid_hdr)
2737 0 : NCF_CHECK_MSG(ncerr," outxfhist : enddef call ")
2738 : else
2739 : ! check that the dimensions are correct
2740 0 : ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id)
2741 0 : NCF_CHECK_MSG(ncerr," outxfhist : inquire nxfh")
2742 : ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=nxfh_id,&
2743 0 : & len=nxfh_tmp)
2744 0 : NCF_CHECK_MSG(ncerr," outxfhist : get nxfh")
2745 0 : ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="xfdim2",dimid=xfdim2_id)
2746 0 : NCF_CHECK_MSG(ncerr," outxfhist : inquire xfdim2")
2747 : ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=xfdim2_id,&
2748 0 : & len=xfdim2_tmp)
2749 0 : NCF_CHECK_MSG(ncerr," outxfhist : get xfdim2")
2750 0 : ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="mxfh",dimid=mxfh_id)
2751 0 : NCF_CHECK_MSG(ncerr," outxfhist : inquire mxfh")
2752 : ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=mxfh_id,&
2753 0 : & len=mxfh_tmp)
2754 0 : NCF_CHECK_MSG(ncerr," outxfhist : get mxfh")
2755 0 : ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="dim2inout",dimid=dim2inout_id)
2756 0 : NCF_CHECK_MSG(ncerr," outxfhist : inquire dim2inout")
2757 : ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=dim2inout_id,&
2758 0 : & len=dim2inout_tmp)
2759 0 : NCF_CHECK_MSG(ncerr," outxfhist : get dim2inout")
2760 :
2761 0 : ncerr = nf90_inq_varid(ncid=ncid_hdr,name="xfhist",varid=xfhist_id)
2762 0 : NCF_CHECK_MSG(ncerr," outxfhist : inquire xfhist")
2763 :
2764 0 : if (mxfh_tmp /= ab_xfh%mxfh .or. dim2inout_tmp /= 2 .or. xfdim2_tmp /= xfdim2) then
2765 0 : write (msg,"(A)") 'outxfhist : ERROR xfhist has bad dimensions in NetCDF file. Can not re-write it.'
2766 0 : ABI_ERROR(msg)
2767 : end if
2768 :
2769 : end if
2770 :
2771 : ! Now fill the data
2772 0 : ncerr = nf90_put_var(ncid=ncid_hdr,varid=xfhist_id,values=ab_xfh%xfhist)
2773 0 : NCF_CHECK_MSG(ncerr," outxfhist : fill xfhist")
2774 :
2775 : ! end NETCDF definition ifdef
2776 : end if ! end iomode if
2777 :
2778 : ! ### (Option=2) Read in number of iterations
2779 : ! #####################################################################
2780 0 : else if ( option == 2 ) then
2781 :
2782 0 : if (wff2%iomode == IO_MODE_FORTRAN) then
2783 0 : read(unit=wff2%unwff,iostat=ios)ab_xfh%nxfh
2784 :
2785 0 : else if (wff2%iomode == IO_MODE_FORTRAN_MASTER) then
2786 : ! FIXME: should copy the xfhist to other processors, and check that we are on the master to read in this case
2787 : ! if node is master
2788 0 : write(msg, "(A,A,A,A)") ch10, " outxfhist: ERROR -", ch10, &
2789 0 : & 'iomode == -1 (localrdwf ) has not been coded yet for xfhist rereading.'
2790 0 : ABI_ERROR(msg)
2791 :
2792 0 : read(unit=wff2%unwff,iostat=ios)ab_xfh%nxfh
2793 :
2794 0 : else if (wff2%iomode == IO_MODE_MPI) then
2795 0 : call xderiveRRecInit(wff2,ierr)
2796 0 : call xderiveRead(wff2,ab_xfh%nxfh,ierr)
2797 0 : call xderiveRRecEnd(wff2,ierr)
2798 :
2799 0 : else if (wff2%iomode == IO_MODE_NETCDF) then
2800 0 : ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id)
2801 0 : NCF_CHECK_MSG(ncerr," outxfhist : inquire nxfh")
2802 : ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=nxfh_id,&
2803 0 : & len=ab_xfh%nxfh)
2804 0 : NCF_CHECK_MSG(ncerr," outxfhist : get nxfh")
2805 : end if
2806 :
2807 : ! ### (Option=3) Read in iteration content
2808 : ! #####################################################################
2809 0 : else if ( option == 3 ) then
2810 0 : if (wff2%iomode == IO_MODE_FORTRAN) then
2811 0 : do ixfh=1,ab_xfh%nxfhr
2812 0 : read(unit=wff2%unwff,iostat=ios)ab_xfh%xfhist(:,:,:,ixfh)
2813 : end do
2814 0 : else if (wff2%iomode == IO_MODE_FORTRAN_MASTER) then
2815 : ! FIXME: should copy the xfhist to other processors, and check that we are on the master to read in this case
2816 : ! if node is master
2817 0 : write(msg, "(A,A,A,A)") ch10, " outxfhist: ERROR -", ch10, &
2818 0 : & 'iomode == -1 (localrdwf ) has not been coded yet for xfhist rereading.'
2819 0 : ABI_ERROR(msg)
2820 :
2821 0 : do ixfh=1,ab_xfh%nxfhr
2822 0 : read(unit=wff2%unwff,iostat=ios)ab_xfh%xfhist(:,:,:,ixfh)
2823 : end do
2824 :
2825 0 : else if (wff2%iomode == IO_MODE_MPI) then
2826 0 : ABI_MALLOC(xfhist_tmp,(3*(natom+4)*2))
2827 0 : spaceComm=xmpi_comm_self
2828 0 : do ixfh=1,ab_xfh%nxfhr
2829 0 : call xderiveRRecInit(wff2,ierr)
2830 0 : call xderiveRead(wff2,xfhist_tmp,3*(natom+4)*2,spaceComm,ierr)
2831 0 : call xderiveRRecEnd(wff2,ierr)
2832 0 : xfhist_tmp(:)=xfhist_tmp(:)
2833 : end do
2834 0 : ABI_FREE(xfhist_tmp)
2835 : end if
2836 :
2837 : ! FIXME: should this be inside the if not mpi as above for options 1 and 2?
2838 : ! it is placed here because the netcdf read is a single operation
2839 0 : if (wff2%iomode == IO_MODE_NETCDF) then
2840 0 : ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id)
2841 0 : NCF_CHECK_MSG(ncerr," outxfhist : inquire nxfh")
2842 : ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=nxfh_id,&
2843 0 : & len=ab_xfh%nxfhr)
2844 0 : NCF_CHECK_MSG(ncerr," outxfhist : get nxfh")
2845 :
2846 0 : ncerr = nf90_inq_varid(ncid=ncid_hdr,varid=xfhist_id,name="xfhist")
2847 0 : NCF_CHECK_MSG(ncerr," outxfhist : inquire xfhist")
2848 : ncerr = nf90_get_var(ncid=ncid_hdr,varid=xfhist_id,values=ab_xfh%xfhist,&
2849 0 : & start=(/1,1,1,1/),count=(/3,natom+4,2,ab_xfh%nxfhr/))
2850 0 : NCF_CHECK_MSG(ncerr," outxfhist : read xfhist")
2851 : end if
2852 :
2853 : else
2854 : ! write(std_out,*)' outxfhist : option ', option , ' not available '
2855 0 : write(msg, "(A,A,A,A,I3,A)") ch10, "outxfhist: ERROR -", ch10, &
2856 0 : & "option ", option, " not available."
2857 0 : ABI_ERROR(msg)
2858 : end if
2859 :
2860 0 : end subroutine outxfhist
2861 : !!***
2862 :
2863 : end module m_gstate
2864 : !!***
|
