Branch data 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 : 5650 : subroutine gstate(args_gs,acell,codvsn,cpui,dtfil,dtset,iexit,initialized,&
237 : 5650 : & itimimage_gstate,mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,&
238 : 5650 : & 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 [ + + + + : 355950 : type(crystal_t) :: cryst
+ + + + +
+ + + + +
+ + + + +
+ ]
290 : 11300 : type(ebands_t) :: bstruct, ebands, ebands_bz
291 : 5650 : type(efield_type) :: dtefield
292 : : type(electronpositron_type),pointer :: electronpositron
293 : 5650 : type(hdr_type) :: hdr, hdr_den, hdr_bz
294 : : type(extfpmd_type),pointer :: extfpmd => null()
295 : : type(rcpaw_type), pointer :: rcpaw => null()
296 [ + + + + : 62150 : type(macro_uj_type) :: dtpawuj(1)
+ - + - +
- + - + -
+ - + - ]
297 : 5650 : type(paw_dmft_type) :: paw_dmft
298 : : type(pawfgr_type) :: pawfgr
299 : 5650 : type(recursion_type) ::rec_set
300 : : type(wffile_type) :: wff1,wffnew,wffnow
301 : 5650 : type(ab_xfh_type) :: ab_xfh
302 : 5650 : type(ddb_type) :: ddb
303 [ + + + + : 355950 : type(ddb_hdr_type) :: ddb_hdr
+ + + + +
+ + + + +
+ + + + +
+ ]
304 : 5650 : type(scfcv_t) :: scfcv_args
305 : : !arrays
306 : : integer :: itimes(2),ngfft(18),ngfftf(18)
307 : 16950 : integer,allocatable :: atindx(:),atindx1(:),indsym(:,:,:),dimcprj_srt(:)
308 : 5650 : integer,allocatable :: irrzon(:,:,:),kg(:,:),nattyp(:),symrec(:,:,:)
309 : 5650 : integer,allocatable,target :: npwarr(:)
310 : 5650 : 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 : 5650 : real(dp),allocatable :: doccde(:)
314 : 5650 : real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phnons(:,:,:),resid(:),rhowfg(:,:)
315 : 5650 : real(dp),allocatable :: rhowfr(:,:),spinat_dum(:,:),start(:,:),work(:)
316 : 5650 : real(dp),allocatable :: ylm(:,:),ylmgr(:,:,:)
317 : : real(dp),contiguous, pointer :: cg(:,:) => null()
318 : 5650 : real(dp),pointer :: eigen(:),pwnsfac(:,:),rhog(:,:),rhor(:,:)
319 : 5650 : real(dp),pointer :: taug(:,:),taur(:,:),xred_old(:,:)
320 : 5650 : type(pawrhoij_type),pointer :: pawrhoij(:)
321 : : type(coulomb_operator) :: kernel_dummy
322 : 5650 : type(pawcprj_type),allocatable :: cprj(:,:)
323 : 5650 : type(xg_nonlop_t) :: xg_nonlop
324 : : ! ***********************************************************************
325 : :
326 : : DBG_ENTER("COLL")
327 : :
328 : 5650 : call timab(1232,1,tsec)
329 : 5650 : call timab(1211,3,tsec)
330 : :
331 : : !###########################################################
332 : : !### 01. Initializations XML, MPI, WVL, etc
333 : :
334 : : !Init MPI data
335 : 5650 : comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm)
336 : :
337 : : !Set up MPI information from the dataset
338 : 5650 : my_natom=mpi_enreg%my_natom
339 : :
340 : : !Set up information when wavelets are in use
341 [ - + ]: 5650 : 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 [ + + + + ]: 5650 : 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 : 5650 : call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
378 : :
379 : : !Structured debugging if prtvol==-level
380 [ - + ]: 5650 : 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 : 5650 : ecore=zero
389 [ + + ]: 22600 : results_gs%pel(1:3) =zero
390 [ + + + + ]: 58746 : results_gs%grchempottn(:,:)=zero
391 [ + + + + ]: 58746 : 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 [ + + - + ]: 5650 : 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 [ + + + + ]: 5674 : results_gs%grvdw(:,:)=zero
399 : : end if
400 : 5650 : 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 : 5650 : 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 : 5650 : rprimd_for_kg=rprimd
411 [ + + + - : 5818 : if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=args_gs%rprimd_orig
+ + + + ]
412 [ + + - + : 5650 : if (dtset%optcell/=0.and.dtset%imgmov/=0) rprimd_for_kg=args_gs%rprimd_orig
- - - - ]
413 : 5650 : call matr3inv(rprimd_for_kg,gprimd_for_kg)
414 [ + + + + : 226000 : gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg)
+ + ]
415 : :
416 : : !Set up the basis sphere of planewaves
417 [ - + + - : 16950 : ABI_MALLOC(npwarr,(dtset%nkpt))
+ - + - -
+ - + ]
418 [ + - + + ]: 5650 : if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2) then
419 [ - + + - : 16932 : 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 : 5644 : & dtset%mpw,npwarr,npwtot,dtset%nsppol)
423 : : call bandfft_kpt_init1(bandfft_kpt,dtset%istwfk,kg,dtset%mgfft,dtset%mkmem,mpi_enreg,&
424 : 5644 : & 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 : 5650 : gemm_nonlop_use_gemm = .false.
435 : :
436 : 5650 : gemm_nonlop_is_distributed = .false.
437 [ - + ]: 5650 : if(dtset%gpu_nl_distrib == 1) gemm_nonlop_is_distributed = .true.
438 [ - + ]: 5650 : if(dtset%gpu_nl_splitsize > 0) gemm_nonlop_block_size = dtset%gpu_nl_splitsize
439 : :
440 [ + - + + ]: 5650 : 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 : 5650 : hamilt_gpu_nfft_blocks = dtset%gpu_nfft_blocks
447 : :
448 : : !Set up the Ylm for each k point
449 [ + + ]: 5650 : if ( dtset%tfkinfunc /= 2) then
450 [ - + + - : 22576 : ABI_MALLOC(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm))
+ + + - +
- + + - +
- + ]
451 [ - + + - : 22576 : ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm))
+ + + - +
- + + - +
- + ]
452 [ + + ]: 5644 : if (psps%useylm==1) then
453 : 1359 : ylm_option=0
454 [ + - + + : 1359 : if (dtset%prtstm==0.and.dtset%iscf>0.and.dtset%positron/=1) ylm_option=1 ! compute gradients of YLM
+ + ]
455 [ + + + - : 1359 : 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 : 1359 : & 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 [ + + ]: 5650 : if (initialized==0) then
467 : : ! This call has to be done before any use of SCF history
468 : 4252 : usecg=0
469 [ + + + + ]: 4252 : if(dtset%extrapwf>0 .or. dtset%imgwfstor==1)usecg=1
470 : 4252 : call scf_history_init(dtset,mpi_enreg,usecg,scf_history)
471 : : end if
472 [ + + + + ]: 5650 : has_to_init=(initialized==0.or.scf_history%history_size<0)
473 : :
474 : 5650 : call timab(1211,2,tsec)
475 : 5650 : call timab(1212,3,tsec)
476 : :
477 : : !###########################################################
478 : : !### 03. Calls pspini
479 : :
480 : : !Open and read pseudopotential files
481 [ - + ]: 5650 : comm_psp=mpi_enreg%comm_cell;if (dtset%usewvl==1) comm_psp=mpi_enreg%comm_wvl
482 [ + + + + : 5674 : if (dtset%nimage>1) psps%mixalch(:,:)=args_gs%mixalch(:,:) ! mixalch can evolve for some image algos
+ + ]
483 : 5650 : call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,psps,rprimd,comm_mpi=comm_psp)
484 : 5650 : call timab(1212,2,tsec)
485 : 5650 : 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 [ - + ]: 5650 : if (dtset%icoulomb == 1) ecore = zero
490 : :
491 : : !WVL - Now that psp data are available, we compute rprimd, acell... from the atomic positions.
492 [ - + ]: 5650 : 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 [ + + ]: 5650 : if (dtset%paral_kgb/=0) then ! We decide to store total npw in bstruct,
549 [ - + + - : 1821 : ABI_MALLOC(npwarr_,(dtset%nkpt))
+ - + - -
+ - + ]
550 [ + - + + : 7461 : npwarr_(:)=npwarr(:)
+ + ]
551 : 607 : call xmpi_sum(npwarr_,mpi_enreg%comm_bandfft,ierr)
552 : : else
553 : 5043 : npwarr_ => npwarr
554 : : end if
555 : :
556 [ - + ]: 5650 : call bstruct%from_dtset(dtset, npwarr_)
557 : :
558 [ + + ]: 5650 : if (dtset%paral_kgb/=0) then
559 [ - + ]: 607 : ABI_FREE(npwarr_)
560 : : end if
561 : 5650 : nullify(npwarr_)
562 : :
563 : : !Initialize PAW atomic occupancies
564 [ + + ]: 5650 : if (scf_history%history_size>=0) then
565 : 1454 : pawrhoij => scf_history%pawrhoij_last
566 : : else
567 [ + + + - : 14392 : ABI_MALLOC(pawrhoij,(my_natom*psps%usepaw))
+ - + + -
+ - + + +
+ - ]
568 : : end if
569 [ + + + + ]: 5650 : 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 [ - + ]: 1091 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
574 : : end if
575 : :
576 : : !Initialize header
577 : 5650 : gscase=0
578 : : call hdr%init(bstruct,codvsn,dtset,pawtab,gscase,psps,wvl%descr,&
579 : 5650 : 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 : 5650 : call bstruct%free()
583 : :
584 : : !Update header, with evolving variables, when available
585 : : !Here, rprimd, xred and occ are available
586 : 5650 : 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 [ + - - - : 5650 : 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 [ + - + + ]: 5650 : if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2) then
593 [ + + ]: 75489 : npwmin=minval(hdr%npwarr(:))
594 [ - + ]: 5644 : 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 [ + + ]: 5644 : 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 [ + - ]: 5650 : if (dtset%usewvl == 0) then
621 : 5650 : 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 [ - + + - : 22600 : ABI_MALLOC(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
- + + - +
- + - - +
- + - + ]
630 [ - + + - : 22600 : ABI_MALLOC(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
- + + - +
- + - - +
- + - + ]
631 [ - + + - : 22600 : ABI_MALLOC(indsym,(4,dtset%nsym,dtset%natom))
- + + - +
- + - - +
- + ]
632 [ - + + - : 16950 : ABI_MALLOC(symrec,(3,3,dtset%nsym))
+ - + - -
+ - + ]
633 [ + + + + : 133484556 : irrzon(:,:,:)=0
+ + ]
634 [ + + + + : 200204129 : phnons(:,:,:)=zero
+ + ]
635 [ + + + + : 1168504 : indsym(:,:,:)=0
+ + ]
636 [ + + + + : 1373965 : symrec(:,:,:)=0
+ + ]
637 : :
638 [ + + ]: 5650 : 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 [ + + + + ]: 15132 : symrec(:,:,1) = 0
648 : 1164 : symrec(1,1,1) = 1 ; symrec(2,2,1) = 1 ; symrec(3,3,1) = 1
649 : : end if
650 [ - + ]: 5650 : 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 : 5650 : cnt = 0
666 : 5650 : nbandtot = 0
667 [ + + ]: 11841 : do spin=1,dtset%nsppol
668 [ + + ]: 83764 : do ikpt=1,dtset%nkpt
669 : 71923 : nband_k = dtset%nband(ikpt + (spin-1) * dtset%nkpt)
670 [ + + ]: 517374 : do band=1,nband_k
671 [ + + ]: 517374 : if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, band, band, spin, mpi_enreg%me_kpt)) cnt = cnt + 1
672 : : end do
673 : 78114 : nbandtot = nbandtot + nband_k
674 : : end do
675 : : end do
676 : :
677 : 5650 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
678 : 5650 : mcg=dtset%mpw*my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
679 [ + + ]: 5650 : 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 [ + - + - : 5650 : if (dtset%usewvl == 0 .and. dtset%mpw > 0 .and. cnt /= 0)then
+ + ]
686 [ - + - + ]: 5608 : 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 [ + + ]: 5650 : if (dtset%imgwfstor==1) then
704 : 53 : cg => scf_history%cg(:,:,1)
705 : 53 : eigen => scf_history%eigen(:,1)
706 : : else
707 [ + - ]: 5597 : 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 [ + + + - : 16791 : ABI_MALLOC_OR_DIE(cg,(2,mcg), ierr)
+ - + + +
- + - -
+ ]
713 : : end if
714 [ - + + - : 16791 : ABI_MALLOC(eigen,(dtset%mband*dtset%nkpt*dtset%nsppol))
+ - + - -
+ - + ]
715 : : end if
716 : :
717 [ - + + - : 16950 : ABI_MALLOC(resid,(dtset%mband*dtset%nkpt*dtset%nsppol))
+ - + - -
+ - + ]
718 [ + + + + ]: 896704 : eigen(:)=zero ; resid(:)=zero
719 : : !mpi_enreg%paralbd=0 ; ask_accurate=0
720 : 5650 : ask_accurate=0
721 : :
722 : : !WVL - Branching, allocating wavefunctions as wavelets.
723 [ - + ]: 5650 : 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 [ + + + + : 5650 : read_wf_or_den=(dtset%iscf<=0.or.dtfil%ireadwf/=0.or.(dtfil%ireadden/=0.and.dtset%positron<=0))
+ + + + ]
745 : 5650 : read_wf_or_den=(read_wf_or_den.and.has_to_init)
746 : :
747 : : !RECURSION - initialization
748 [ + + + + ]: 5650 : 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 [ + + + + ]: 5650 : if(dtset%imgwfstor==1 .and. initialized==1)then
763 [ + + + + ]: 1012074 : cg(:,:)=scf_history%cg(:,:,1)
764 [ + + ]: 714 : eigen(:)=scf_history%eigen(:,1)
765 [ + + ]: 5608 : else if(dtset%tfkinfunc /=2) then
766 : : !if(dtset%tfkinfunc /=2) then
767 : 5602 : wff1%unwff=dtfil%unwff1
768 : 5602 : optorth=1 !if (psps%usepaw==1) optorth=0
769 [ + + + + ]: 5602 : if(psps%usepaw==1 .and. dtfil%ireadwf==1)optorth=0
770 [ + + + + ]: 72826 : 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 [ - + ]: 5602 : & dtfil%fnamewffk,wvl)
779 : : ABI_NVTX_END_RANGE()
780 [ + + + + ]: 72826 : hdr%rprimd=rprimd
781 : : end if
782 : :
783 [ + + + + ]: 5650 : 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 [ + + - + : 5650 : 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 [ - + - - ]: 5650 : 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 : 5650 : ab_xfh%nxfh=0
803 [ - + - - ]: 5650 : 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 : 5650 : ab_xfh%mxfh=(ab_xfh%nxfh-dtset%restartxf+1)+dtset%ntime+5
836 [ + - - + : 28250 : ABI_MALLOC(ab_xfh%xfhist,(3,dtset%natom+4,2,ab_xfh%mxfh))
+ - + - -
+ + - + -
+ - - + -
+ - + ]
837 [ + + + + : 2556987 : 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 [ + - - + ]: 5650 : 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 [ + + ]: 5650 : 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 [ - + + - : 16950 : ABI_MALLOC(start,(3,dtset%natom))
+ - + - -
+ - + ]
864 : 5650 : call setup2(dtset,npwtot,start,wvl%wfs,xred)
865 : :
866 : : !Allocation of previous atomic positions
867 [ + + ]: 5650 : if (scf_history%history_size>=0) then
868 : 1454 : xred_old => scf_history%xred_last
869 : : else
870 [ - + + - : 12588 : ABI_MALLOC(xred_old,(3,dtset%natom))
+ - + - -
+ - + ]
871 : : end if
872 [ + + + + : 48226 : if (has_to_init) xred_old=xred
+ + ]
873 : :
874 : : !Initialize (eventually) extfpmd object
875 [ + + ]: 5650 : 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 : 5650 : call timab(1211,2,tsec)
886 : 5650 : 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 [ + + + + : 5650 : & (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 : 5357 : 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 [ - + + - : 16950 : ABI_MALLOC(atindx,(dtset%natom))
+ - + - -
+ - + ]
936 [ - + + - : 11300 : ABI_MALLOC(atindx1,(dtset%natom))
- + - + ]
937 [ - + + - : 16950 : ABI_MALLOC(nattyp,(psps%ntypat))
+ - + - -
+ - + -
+ ]
938 : 5650 : indx=1
939 [ + + ]: 13679 : do itypat=1,psps%ntypat
940 : 8029 : nattyp(itypat)=0
941 [ + + ]: 34210 : do iatom=1,dtset%natom
942 [ + + ]: 28560 : if(dtset%typat(iatom)==itypat)then
943 : 13274 : atindx(iatom)=indx
944 : 13274 : atindx1(indx)=iatom
945 : 13274 : indx=indx+1
946 : 13274 : 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 [ + + - + ]: 5650 : if ((.not.read_wf_or_den).or.(scf_history%history_size>0.and.has_to_init)) then
953 [ + - - + : 12684 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*dtset%natom))
+ - + - +
- - + - +
- + ]
954 : 4228 : 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 : 5650 : use_sc_dmft=dtset%usedmft
970 : : ! if(dtset%paral_kgb>0) use_sc_dmft=0
971 : : !Should be done inside init_sc_dmft
972 [ + + + + : 5650 : 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 [ + + ]: 5650 : if(psps%usepaw==1) then
987 : : ! 1-
988 [ + + + + : 1306 : 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 : 1306 : call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
993 : :
994 [ + + + + ]: 1306 : if (psp_gencond==1.or.call_pawinit) then
995 : 541 : call timab(553,1,tsec)
996 : 541 : gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
997 [ + + ]: 541 : 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 : 541 : & pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero)
1001 : :
1002 : : ! Update internal values
1003 : 541 : call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
1004 : 541 : 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 [ + + ]: 3153 : psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore)
1015 : 1306 : 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 [ + - - - : 1306 : & 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 : 5650 : & pawrad(:),pawtab(:),rprimd(:,:),ucvol,dtfil%unpaw,use_sc_dmft,xred(:,:),ylm(:,:))
1029 [ + + ]: 5650 : 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 [ + + ]: 5650 : 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 : 5650 : nullify (electronpositron)
1052 [ + + ]: 5650 : 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 [ + + ]: 5650 : 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 : 5650 : usecprj=0; mcprj=0;mband_cprj=0
1078 : 5650 : compute_cprj=.false.
1079 [ + - ]: 5650 : if (dtset%cprj_in_memory==2) then
1080 : : compute_cprj=.true.
1081 : : usecprj=1
1082 : : else
1083 [ + + ]: 5650 : if (dtset%usepaw==1) then
1084 [ + + ]: 1306 : if (associated(electronpositron)) then
1085 [ - + + + ]: 13 : if (dtset%positron/=0.and.electronpositron%dimcprj>0) usecprj=1
1086 : : end if
1087 [ + + ]: 1306 : if (dtset%prtnabla>0) usecprj=1
1088 [ + + ]: 1306 : if (dtset%extrapwf>0) usecprj=1
1089 [ + + ]: 1306 : if (dtset%pawfatbnd>0)usecprj=1
1090 [ + + ]: 1306 : if (dtset%prtdos==3) usecprj=1
1091 [ - + ]: 1306 : if (dtset%usewvl==1) usecprj=1
1092 [ + + ]: 1306 : if (dtset%nstep==0) usecprj=0
1093 [ + + ]: 1306 : if (dtset%usefock==1) usecprj=1
1094 : : end if
1095 : : end if
1096 [ + + ]: 1291 : if (usecprj==0) then
1097 [ - + - + ]: 5591 : 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 : 5650 : call timab(1213,2,tsec)
1142 : 5650 : 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 [ + + ]: 5650 : if (scf_history%history_size>=0) then
1150 : 1454 : rhor => scf_history%rhor_last
1151 : 1454 : taur => scf_history%taur_last
1152 : : else
1153 [ - + + - : 16784 : ABI_MALLOC(rhor,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
1154 [ - + + - : 16784 : ABI_MALLOC(taur,(nfftf,dtset%nspden*dtset%usekden))
+ + + - +
- + + - +
- + ]
1155 : : end if
1156 [ - + + - : 16950 : ABI_MALLOC(rhog,(2,nfftf))
+ - + - -
+ - + ]
1157 [ + + + - : 16950 : ABI_MALLOC(taug,(2,nfftf*dtset%usekden))
+ - + + -
+ - + ]
1158 : :
1159 [ + + ]: 5650 : if (has_to_init) then
1160 : :
1161 : : ! === Self-consistent case
1162 [ + + - + : 4513 : if (dtset%iscf>0 .or. (dtset%iscf==0 .and. dtset%usewvl==1 )) then
- - ]
1163 : :
1164 : : ! >>> Initialize charge density
1165 : :
1166 [ + + + + ]: 4002 : 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 [ + + ]: 3946 : 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 [ + + ]: 3091 : 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 [ + - ]: 3086 : 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 [ - + - + ]: 3086 : & 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 [ + + ]: 3086 : 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 [ + + ]: 4002 : 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 : 5650 : call timab(1214,2,tsec)
1313 : 5650 : 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 [ + + + + ]: 5650 : 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 [ + + - + ]: 5650 : if ((.not.read_wf_or_den).or.(scf_history%history_size>0.and.has_to_init)) then
1341 [ - + ]: 4228 : 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 : 5650 : & pwind,pwind_alloc,pwnsfac,rprimd,symrec,xred)
1349 : :
1350 : 5650 : 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 : 5650 : cpus=dtset%cpus
1357 [ - + ]: 5650 : if(abs(cpus)>1.0d-5)cpus=cpus+cpui
1358 [ + + ]: 5650 : openexit=1 ; if(dtset%chkexit==0) openexit=0
1359 : 5650 : 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 [ - + - - ]: 5650 : 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 : 5650 : call timab(1215,2,tsec)
1371 : 5650 : call pstat_proc%print(_PSTAT_ARGS_)
1372 : :
1373 : 5650 : conv_retcode = 0
1374 : :
1375 [ + - ]: 5650 : if (iexit==0) then
1376 : :
1377 : : ! ###########################################################
1378 : : ! ### 14. Move atoms and acell according to ionmov value
1379 : :
1380 : 5650 : 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 [ - + - + : 5650 : & symrec,taug,taur,wvl,ylm,ylmgr,paw_dmft,wffnew,wffnow,xg_nonlop)
- + - + ]
1389 : :
1390 : 5650 : call dtfil%init_time(0)
1391 : :
1392 [ + - + + ]: 457650 : write(msg,'(a,80a)')ch10,('=',mu=1,80)
1393 [ + + ]: 16950 : call wrtout([std_out, ab_out], msg)
1394 : :
1395 [ + + + + ]: 5650 : 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 [ + + ]: 5301 : if (dtset%macro_uj==0) then
1399 : 5291 : itimes(1)=itime0 ; itimes(2)=itimimage_gstate
1400 [ - + ]: 5291 : 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 : 5650 : call scfcv_destroy(scfcv_args)
1434 : :
1435 : 5650 : 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 : 5650 : call timab(1226,3,tsec)
1443 : :
1444 [ + + ]: 457650 : write(msg, '(80a,a,a,a,a)' ) ('=',mu=1,80),ch10,ch10,&
1445 [ + - ]: 463300 : & ' ----iterations are completed or convergence reached----',ch10
1446 [ + + ]: 16950 : call wrtout([std_out, ab_out], msg)
1447 : :
1448 : : !Mark this GS computation as done
1449 : 5650 : 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 [ + - - - : 5650 : comm_atom=mpi_enreg%comm_atom,extfpmd_eshift=results_gs%extfpmd_eshift,mpi_atmtab=mpi_enreg%my_atmtab)
- - - + ]
1455 : :
1456 [ - + + - : 16950 : ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
+ - + - -
+ - + ]
1457 [ + + ]: 451177 : 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 [ - + ]: 5650 : hdr%kptrlatt, hdr%nshiftk, hdr%shiftk)
1464 : :
1465 : 5650 : ebands%fermie = results_gs%energies%e_fermie
1466 : 5650 : ebands%fermih = results_gs%energies%e_fermih
1467 [ - + ]: 5650 : ABI_FREE(doccde)
1468 : :
1469 : : ! Compute and print the gaps.
1470 : 5650 : call ebands%report_gap(header="Gap info",unit=std_out,mode_paral="COLL",gaps=results_gs%gaps)
1471 : :
1472 : 5650 : call timab(1226,2,tsec)
1473 : 5650 : call timab(1227,3,tsec)
1474 : :
1475 [ + + ]: 5650 : if(dtset%nqpt==0)filnam=dtfil%fnameabo_wfk
1476 [ + + ]: 5650 : 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 : 5650 : write_wfk = .True.
1481 [ - + - - ]: 5650 : 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 [ + + + + ]: 73450 : hdr%rprimd=rprimd_for_kg
1490 : :
1491 [ + - ]: 5650 : if (write_wfk) then
1492 : 5650 : 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 [ - + ]: 5650 : 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 [ + + ]: 5650 : 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 [ - + - - ]: 5650 : 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 [ + + + + ]: 73450 : hdr%rprimd=rprimd
1511 : :
1512 : : ! Generate WFK in full BZ (needed by LOBSTER)
1513 [ + + + + : 5650 : 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 : 5650 : call timab(1227,2,tsec)
1521 : 5650 : 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 [ - + - + ]: 5650 : & resid,rhor,rprimd,results_gs%vxcavg,xred,rcpaw)
1527 : :
1528 [ + + + + ]: 5650 : if ( (dtset%iscf>=0 .or. dtset%iscf==-3) .and. dtset%prtstm==0) then
1529 : 5163 : call prtene(dtset,results_gs%energies,ab_out,psps%usepaw)
1530 : : end if
1531 : :
1532 : 5650 : call timab(1228,2,tsec)
1533 : 5650 : 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 [ + + ]: 5650 : & 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 : 5650 : call timab(1229,2,tsec)
1573 : 5650 : call timab(1230,3,tsec)
1574 : :
1575 : : !In the // case, only master writes the energy and the gradients to the DDB
1576 [ + + + + : 5650 : 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 [ + + ]: 2821 : if (dtset%iscf > 0) then
1581 : 2814 : nblok = 2 ! 1st blok = gradients, 2nd blok = energy
1582 : : else
1583 : 7 : nblok = 1 ! 1st blok = gradients
1584 : : end if
1585 : 2821 : mpert = dtset%natom + 6
1586 : :
1587 : : ! Create header and ddb objects
1588 : 2821 : dscrpt=' Note : temporary (transfer) database '
1589 : : call ddb_hdr%init(dtset,psps,pawtab,dscrpt,nblok,&
1590 : 2821 : & xred=xred,occ=occ,ngfft=ngfft)
1591 : :
1592 : 2821 : call ddb%init(dtset, nblok, mpert, with_d1E=.true.)
1593 : :
1594 : : ! Set the electronic polarization
1595 [ + + ]: 2821 : 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 [ + + ]: 2821 : if (dtset%iscf > 0) then
1600 : :
1601 : : ! Set the gradients (forces) in reduced coordinates
1602 : 2814 : call ddb%set_gred(results_gs%gred, 1)
1603 : :
1604 : : ! Set the stress tensor
1605 : 2814 : call ddb%set_strten(results_gs%strten, 1)
1606 : :
1607 : : ! Set the total energy
1608 : 2814 : call ddb%set_etotal(results_gs%etotal, 2)
1609 : : end if
1610 : :
1611 [ + + ]: 2821 : if (dtset%prtddb==1) then
1612 : : ! Write the DDB
1613 : 2778 : call ddb%write(ddb_hdr, dtfil%fnameabo_ddb, with_psps=0)
1614 : : end if
1615 : :
1616 : : ! Free memory
1617 : 2821 : call ddb_hdr%free()
1618 : 2821 : call ddb%free()
1619 : :
1620 : : end if
1621 : :
1622 : 5650 : call timab(1230,2,tsec)
1623 : 5650 : call timab(1231,3,tsec)
1624 : :
1625 : :
1626 [ + + + + : 5650 : 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 : 5607 : & results_gs%strten,results_gs%synlgr,xred)
1631 : : end if
1632 : :
1633 : : ! Write nonlop_ylm_counters (if enabled) in outputs
1634 [ + + - + : 5650 : 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 [ - + - - ]: 5650 : 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 [ + + ]: 5650 : 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 [ - + ]: 5650 : ABI_FREE(atindx)
1649 [ - + ]: 5650 : ABI_FREE(atindx1)
1650 [ - + ]: 5650 : ABI_FREE(indsym)
1651 [ - + ]: 5650 : ABI_FREE(npwarr)
1652 [ - + ]: 5650 : ABI_FREE(nattyp)
1653 [ - + ]: 5650 : ABI_FREE(resid)
1654 [ - + ]: 5650 : ABI_FREE(rhog)
1655 [ - + ]: 5650 : ABI_FREE(start)
1656 [ - + ]: 5650 : ABI_FREE(symrec)
1657 [ - + ]: 5650 : ABI_FREE(taug)
1658 [ - + ]: 5650 : ABI_FREE(ab_xfh%xfhist)
1659 : 5650 : call pawfgr_destroy(pawfgr)
1660 [ + + ]: 5650 : 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 [ + + ]: 5650 : if(dtset%imgwfstor==0)then
1670 [ + - ]: 5597 : 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 [ - + ]: 5597 : ABI_FREE(cg)
1676 : : end if
1677 [ - + ]: 5597 : ABI_FREE(eigen)
1678 : : else
1679 : 53 : nullify(cg,eigen)
1680 : : endif
1681 : :
1682 [ + - ]: 5650 : if (dtset%usewvl == 0) then
1683 : : ! In wavelet case, irrzon and phnons are deallocated by wavelet object.
1684 [ - + ]: 5650 : ABI_FREE(irrzon)
1685 [ - + ]: 5650 : ABI_FREE(phnons)
1686 : : end if
1687 : :
1688 [ - + ]: 5650 : ABI_FREE(ylm)
1689 [ - + ]: 5650 : ABI_FREE(ylmgr)
1690 : :
1691 [ + + ]: 5650 : if (scf_history%history_size<0) then
1692 [ + + ]: 4196 : if (psps%usepaw==1) then
1693 : 977 : call pawrhoij_free(pawrhoij)
1694 : : end if
1695 [ - + ]: 4196 : ABI_FREE(rhor)
1696 [ - + ]: 4196 : ABI_FREE(taur)
1697 [ - + ]: 4196 : ABI_FREE(pawrhoij)
1698 [ - + ]: 4196 : ABI_FREE(xred_old)
1699 : : else
1700 : 1454 : nullify(rhor,taur,pawrhoij,xred_old)
1701 : : end if
1702 : :
1703 : : !PAW+DMFT
1704 : 5650 : call destroy_sc_dmft(paw_dmft)
1705 : : ! This call should be done inside destroy_sc_dmft
1706 [ + + + + ]: 5650 : 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 [ + + ]: 5650 : if(associated(extfpmd)) then
1712 : 7 : call extfpmd%destroy()
1713 [ - + - + : 7 : ABI_FREE(extfpmd)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ ]
1714 : : end if
1715 : :
1716 : : !Destroy rcpaw datastructure
1717 [ + + ]: 5650 : if(associated(rcpaw)) then
1718 : 3 : call rcpaw_destroy(rcpaw)
1719 [ - + - + : 3 : ABI_FREE(rcpaw)
- + - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- + - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- + - - -
- - - -
- ]
1720 : : endif
1721 : :
1722 : : !Destroy electronpositron datastructure
1723 [ + + ]: 5650 : if (dtset%positron/=0) then
1724 : 14 : call destroy_electronpositron(electronpositron)
1725 : : end if
1726 : :
1727 : : !Deallocating the basis set.
1728 [ - + ]: 5650 : 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 [ - + ]: 5650 : ABI_FREE(kg)
1739 : :
1740 [ - + ]: 5650 : 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 [ + - ]: 5650 : if (associated(pwind)) then
1746 : 5650 : ABI_FREE(pwind)
1747 : : end if
1748 [ + - ]: 5650 : if (associated(pwnsfac)) then
1749 : 5650 : ABI_FREE(pwnsfac)
1750 : : end if
1751 [ + + + + ]: 5650 : 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 [ + + ]: 5650 : if(usecprj==1) then
1767 [ - + ]: 59 : ABI_FREE(dimcprj_srt)
1768 : 59 : call pawcprj_free(cprj)
1769 : : end if
1770 [ + - + + : 16610 : ABI_FREE(cprj)
- + - + -
+ ]
1771 : :
1772 : : ! deallocate efield
1773 : 5650 : call dtefield%free()
1774 : :
1775 [ + + ]: 5650 : if (dtset%userec == 1) call CleanRec(rec_set)
1776 : :
1777 : 5650 : call hdr%free()
1778 : 5650 : call ebands%free()
1779 : :
1780 [ + + + + ]: 5650 : 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 [ + - + + : 5650 : if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2 .and. dtset%optdriver /= RUNL_GWLS) then
+ + ]
1788 : : ! Plane-wave case
1789 : 5636 : call bandfft_kpt_destroy_array(bandfft_kpt,mpi_enreg)
1790 : : end if
1791 : :
1792 : : !Clean gemm_nonlop work spaces
1793 [ + + ]: 5650 : 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 : 5650 : call timab(1231,2,tsec)
1812 : 5650 : call timab(1232,2,tsec)
1813 : :
1814 : : DBG_EXIT("COLL")
1815 : :
1816 [ - + - + : 28250 : 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 : 5650 : 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 [ + + ]: 5650 : if (dtset%iscf>=0) then
1861 : :
1862 : : ! Copy coordinates into array start
1863 [ + + + + ]: 54071 : start(:,:)=xred(:,:)
1864 : :
1865 [ + - ]: 5139 : 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 : 5139 : arith=zero
1873 : 5139 : geom=one
1874 : 5139 : wtknrm=zero
1875 [ + + ]: 57811 : do ikpt=1,dtset%nkpt
1876 : 52672 : npw=npwtot(ikpt)
1877 : 52672 : wtknrm=wtknrm+dtset%wtk(ikpt)
1878 : 52672 : arith=arith+npw*dtset%wtk(ikpt)
1879 : 57811 : geom=geom*npw**dtset%wtk(ikpt)
1880 : : end do
1881 : :
1882 : : ! Enforce normalization of weights to 1
1883 : 5139 : arith=arith/wtknrm
1884 : 5139 : geom=geom**(1.0_dp/wtknrm)
1885 : :
1886 : : end if
1887 : :
1888 : : ! Ensure portability of output thanks to tol8
1889 [ + - ]: 5139 : if (dtset%usewvl == 0) then
1890 : 5139 : 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 [ + + ]: 15417 : 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 : 5650 : 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 : 5650 : subroutine clnup1(acell,dtset,eigen,fermie,fermih, fnameabo_dos,fnameabo_eig,gred,&
1963 : 5650 : 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 : 5650 : real(dp),allocatable :: doccde(:)
1994 : :
1995 : : ! ****************************************************************
1996 : :
1997 : 5650 : comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm)
1998 : :
1999 [ + + ]: 5650 : if(dtset%prtstm==0)then ! Write reduced coordinates xred
2000 : 5647 : write(msg, '(a,i5,a)' )' reduced coordinates (array xred) for',dtset%natom,' atoms'
2001 : 5647 : call wrtout(ab_out, msg)
2002 [ + + ]: 18906 : do iatom=1,dtset%natom
2003 : 13259 : write(msg, '(1x,3f20.12)' ) xred(:,iatom)
2004 : 18906 : 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 [ + + + + : 5650 : 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 : 5116 : grmax=0.0_dp
2013 : 5116 : grsum=0.0_dp
2014 [ + + ]: 17274 : do iatom=1,dtset%natom
2015 [ + + ]: 53748 : do ii=1,3
2016 : : ! To be activated in v5.5
2017 : : ! grmax=max(grmax,abs(gred(ii,iatom)))
2018 : 36474 : grmax=max(grmax,gred(ii,iatom))
2019 : 48632 : grsum=grsum+gred(ii,iatom)**2
2020 : : end do
2021 : : end do
2022 : 5116 : grsum=sqrt(grsum/dble(3*dtset%natom))
2023 : :
2024 : 5116 : 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 : 5116 : call wrtout(ab_out, msg)
2026 [ + + ]: 17274 : do iatom=1,dtset%natom
2027 : 12158 : write(msg, '(i5,1x,3f20.12)' ) iatom,gred(1:3,iatom)
2028 : 17274 : call wrtout(ab_out, msg)
2029 : : end do
2030 : :
2031 : : end if
2032 : :
2033 [ + + ]: 5650 : if(dtset%prtstm==0)then
2034 : :
2035 : : ! Compute and write out dimensional cartesian coords and forces:
2036 : 5647 : call wrtout(ab_out,' ')
2037 : :
2038 : : ! (only write forces if iscf > 0 and dtset%nstep>0)
2039 [ + + + + : 5647 : if (dtset%iscf<0.or.dtset%nstep<=0.or.prtfor==0) then
+ + ]
2040 : 638 : iwfrc=0
2041 : : else
2042 : 5009 : iwfrc=1
2043 : : end if
2044 : :
2045 : 5647 : call prtxf(gred,dtset%iatfix,ab_out,iwfrc,dtset%natom,rprimd,xred)
2046 : :
2047 : : ! Write length scales
2048 : 5647 : write(msg, '(1x,a,3f16.12,a)' )'length scales=',acell,' bohr'
2049 : 5647 : call wrtout(ab_out, msg)
2050 [ + + ]: 22588 : write(msg, '(14x,a,3f16.12,a)' )'=',Bohr_Ang*acell(1:3),' angstroms'
2051 : 5647 : call wrtout(ab_out, msg)
2052 : :
2053 : : end if
2054 : :
2055 : 5650 : option=1; nnonsc=0; tolwf=0.0_dp
2056 : :
2057 [ + + ]: 5650 : if(dtset%iscf<0 .and. dtset%iscf/=-3)option=3
2058 : 5650 : iscf_dum=dtset%iscf
2059 [ + + ]: 5650 : if(dtset%nstep==0)iscf_dum=-1
2060 : :
2061 [ + + ]: 5650 : if(dtset%tfkinfunc==0)then
2062 [ + + ]: 5640 : 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 : 4389 : & 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 : 4389 : & vxcavg,dtset%wtk)
2073 : : end if
2074 : :
2075 [ + + + + ]: 5640 : if (dtset%prteig==1 .and. me == master) then
2076 [ + - ]: 2695 : filename=trim(fnameabo_eig)//'.nc'
2077 : 2695 : 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 : 5650 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2083 : 5650 : call prtrhomxmn(std_out,mpi_enreg,nfft,ngfft,dtset%nspden,2,rhor,ucvol=ucvol)
2084 [ + + ]: 5650 : 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 [ + + + - ]: 5650 : 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 : 5650 : 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 : 5647 : 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 : 5647 : format_line21='(i5,1x,3f21.14)'
2166 : 5647 : format_line25='(i5,1x,3f25.14)'
2167 : :
2168 : : !Write cartesian coordinates in angstroms
2169 : 5647 : call wrtout(iout,' cartesian coordinates (angstrom) at end:')
2170 [ + + ]: 18906 : do iatom=1,natom
2171 : 13259 : format_line=format_line21
2172 [ + + ]: 53036 : do mu=1,3
2173 : : xx(mu)=(rprimd(mu,1)*xred(1,iatom)+&
2174 : : & rprimd(mu,2)*xred(2,iatom)+&
2175 : 39777 : & rprimd(mu,3)*xred(3,iatom))*Bohr_Ang
2176 [ + - - + ]: 53036 : if(xx(mu)>99999 .or. xx(mu)<-9999)format_line=format_line25
2177 : : end do
2178 : 13259 : write(msg,format_line) iatom,xx
2179 : 18906 : call wrtout(iout, msg)
2180 : : end do
2181 : :
2182 : : !Optionally write cartesian forces in eV/Angstrom (also provide same in hartree/bohr)
2183 [ + + ]: 5647 : if (iwfrc/=0) then
2184 : : ! First, provide results in hartree/bohr
2185 : 5009 : write(msg, '(a,a)' ) ch10,' cartesian forces (hartree/bohr) at end:'
2186 : 5009 : call wrtout(iout, msg)
2187 : 5009 : frms=zero
2188 : 5009 : fmax=zero
2189 : 5009 : favg(1)=zero
2190 : 5009 : favg(2)=zero
2191 : 5009 : 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 : 5009 : call matr3inv(rprimd,gprimd)
2196 : : ! First compute (spurious) average force favg
2197 [ + + ]: 16613 : do iatom=1,natom
2198 [ + + ]: 51425 : do mu=1,3
2199 : : ff(mu)=-(gprimd(mu,1)*gred(1,iatom)+&
2200 : : & gprimd(mu,2)*gred(2,iatom)+&
2201 : 34812 : & gprimd(mu,3)*gred(3,iatom))
2202 : 46416 : favg(mu)=favg(mu)+ff(mu)
2203 : : end do
2204 : : end do
2205 : 5009 : favg(1) = favg(1)/dble(natom)
2206 : 5009 : favg(2) = favg(2)/dble(natom)
2207 : 5009 : 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 : 5009 : unfixd=0
2213 [ + + ]: 16613 : do iatom=1,natom
2214 : 11604 : format_line=format_line21
2215 [ + + ]: 46416 : do mu=1,3
2216 : : ff(mu)=-(gprimd(mu,1)*gred(1,iatom)+&
2217 : : & gprimd(mu,2)*gred(2,iatom)+&
2218 : 34812 : & gprimd(mu,3)*gred(3,iatom))-favg(mu)
2219 [ + - - + ]: 34812 : if(ff(mu)>99999 .or. ff(mu)<-9999)format_line=format_line25
2220 : : ! For rms and max force, include only unfixed components
2221 [ + + ]: 46416 : if (iatfix(mu,iatom) /= 1) then
2222 : 30803 : unfixd=unfixd+1
2223 : 30803 : frms=frms+ff(mu)**2
2224 : 30803 : fmax=max(fmax,abs(ff(mu)))
2225 : : end if
2226 : : end do
2227 : 11604 : write(msg, format_line) iatom,ff
2228 : 16613 : call wrtout(iout, msg)
2229 : : end do
2230 [ + + ]: 5009 : 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 : 5009 : favg_out(:)=favg(:)
2236 [ + + ]: 5009 : if(abs(favg_out(1))<tol14)favg_out(1)=zero
2237 [ + + ]: 5009 : if(abs(favg_out(2))<tol14)favg_out(2)=zero
2238 [ + + ]: 5009 : if(abs(favg_out(3))<tol14)favg_out(3)=zero
2239 : :
2240 : 5009 : write(msg, '(a,1p,2e14.7,1x,3e11.3,a)' )' frms,max,avg=',frms,fmax,favg_out(1:3),' h/b'
2241 : 5009 : call wrtout(iout, msg)
2242 : :
2243 [ + - ]: 5009 : if (iwfrc==1) then
2244 : :
2245 : 5009 : write(msg, '(a,a)' )ch10,' cartesian forces (eV/Angstrom) at end:'
2246 : 5009 : call wrtout(iout, msg)
2247 : 5009 : convt=Ha_eV/Bohr_Ang
2248 : :
2249 : : ! Note: subtract off average force
2250 [ + + ]: 16613 : do iatom=1,natom
2251 : 11604 : format_line=format_line21
2252 [ + + ]: 46416 : do mu=1,3
2253 : : ff(mu)=(-(gprimd(mu,1)*gred(1,iatom)+&
2254 : : & gprimd(mu,2)*gred(2,iatom)+&
2255 : 34812 : & gprimd(mu,3)*gred(3,iatom))-favg(mu))*convt
2256 [ + - + + ]: 46416 : if(ff(mu)>99999 .or. ff(mu)<-9999)format_line=format_line25
2257 : : end do
2258 : 11604 : write(msg, format_line) iatom,ff
2259 : 16613 : call wrtout(iout, msg)
2260 : : end do
2261 [ + + ]: 25045 : 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 : 5009 : call wrtout(iout, msg)
2263 : :
2264 : : end if
2265 : : end if
2266 : :
2267 : 5647 : 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 : 5607 : subroutine clnup2(n1xccc,gred,grchempottn,gresid,grewtn,grvdw,grxc,iscf,natom,ngrvdw,&
2304 : 5607 : & 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 [ + + ]: 5607 : if (iscf>=0) then
2328 : :
2329 [ + + + - ]: 5111 : 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 : 5111 : devsqr=zero
2395 [ + + ]: 17242 : do iatom=1,natom
2396 [ + + ]: 53635 : do mu=1,3
2397 : 48524 : 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 [ + + ]: 5111 : 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 [ + + ]: 5111 : if (prtstr > 0) then
2412 : : ! Write out stress results
2413 [ + + ]: 15048 : units = [std_out, ab_out]
2414 : 5016 : write(msg, '(a,a)' ) ch10,' Cartesian components of stress tensor (hartree/bohr^3)'
2415 : 5016 : call wrtout(units, msg)
2416 : :
2417 : 5016 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4)
2418 : 5016 : call wrtout(units, msg)
2419 : 5016 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5)
2420 : 5016 : call wrtout(units, msg)
2421 : 5016 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) ' sigma(3 3)=',strten(3),' sigma(2 1)=',strten(6)
2422 : 5016 : call wrtout(units, msg)
2423 : :
2424 : : ! Also output the pressure (minus one third the trace of the stress tensor).
2425 : 5016 : write(msg, '(a,a,es12.4,a)' ) ch10,&
2426 : 5016 : '-Cartesian components of stress tensor (GPa) [Pressure=',&
2427 : 10032 : -(strten(1)+strten(2)+strten(3))*HaBohr3_GPa/3.0_dp,' GPa]'
2428 : 5016 : call wrtout(units, msg)
2429 : :
2430 : : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2431 : 5016 : '- sigma(1 1)=',strten(1)*HaBohr3_GPa,&
2432 : 10032 : ' sigma(3 2)=',strten(4)*HaBohr3_GPa
2433 : 5016 : call wrtout(units, msg)
2434 : : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2435 : 5016 : '- sigma(2 2)=',strten(2)*HaBohr3_GPa,&
2436 : 10032 : ' sigma(3 1)=',strten(5)*HaBohr3_GPa
2437 : 5016 : call wrtout(units, msg)
2438 : : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2439 : 5016 : '- sigma(3 3)=',strten(3)*HaBohr3_GPa,&
2440 : 10032 : ' sigma(2 1)=',strten(6)*HaBohr3_GPa
2441 : 5016 : call wrtout(units, msg)
2442 : : end if
2443 : :
2444 : : end if ! iscf > 0
2445 : :
2446 : : !write(std_out,*)' clnup2 : exit '
2447 : :
2448 : 5607 : 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 : : !!***
|
