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