Line data Source code
1 : !!****m* ABINIT/m_common
2 : !! NAME
3 : !! m_common
4 : !!
5 : !! FUNCTION
6 : !! This module gathers routines used by higher-level procedures.
7 : !! Mainly printing routines.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, AF, GMR, LBoeri, MT)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_common
24 :
25 : use defs_basis
26 : use m_errors
27 : use m_abicore
28 : use m_exit
29 : use m_fftcore
30 : use m_fock
31 : use m_io_tools
32 : #if defined DEV_YP_VDWXC
33 : use m_xc_vdw
34 : #endif
35 : use netcdf
36 : use m_nctk
37 : use m_crystal
38 : use m_wfk
39 : use m_ebands
40 : use m_hdr
41 : use m_xmpi
42 : use m_dtset
43 : use m_xpapi
44 : use m_yaml
45 : use m_invars2
46 : use m_dtset
47 :
48 : use m_rcpaw, only : rcpaw_type
49 : use m_fstrings, only : indent, endswith, sjoin, itoa
50 : use m_electronpositron, only : electronpositron_type
51 : use m_energies, only : energies_type
52 : use m_pair_list, only : pair_list
53 : use m_geometry, only : mkrdim, metric
54 : use m_kg, only : getcut
55 : use m_parser, only : parsefile, ab_dimensions
56 : use m_invars1, only : invars0, invars1m, indefo
57 : use m_time, only : timab, time_set_papiopt
58 : use defs_abitypes, only : MPI_type
59 : use defs_datatypes, only : pspheader_type
60 : use m_pspheads, only : inpspheads, pspheads_comm
61 : use m_kpts, only : kpts_timrev_from_kptopt
62 : use m_dft_energy, only : entropy
63 :
64 : use m_xg, only : SPACE_CR,SPACE_C
65 : use m_hamiltonian, only : gs_hamiltonian_type
66 : use m_chebfi2, only : chebfi_memInfo
67 : use m_lobpcg2, only : lobpcg_memInfo
68 : use m_invovl, only : invovl_ompgpu_static_mem,invovl_ompgpu_work_mem
69 : use m_gemm_nonlop, only : gemm_nonlop_ompgpu_static_mem,gemm_nonlop_ompgpu_work_mem
70 : use m_gemm_nonlop_projectors, only : gemm_nonlop_split_choice23
71 : use m_getghc, only : getghc_ompgpu_work_mem
72 : use, intrinsic :: iso_c_binding, only : c_size_t
73 :
74 : #if defined(HAVE_GPU)
75 : use m_gpu_toolbox
76 : #endif
77 :
78 : implicit none
79 :
80 : private
81 : !!***
82 :
83 : public :: scprqt
84 : public :: setup1
85 : public :: prteigrs
86 : public :: prtene
87 : public :: get_dtsets_pspheads ! Parse input file, get list of pseudos for files file and build list of datasets
88 : ! pseudopotential headers, maxval of dimensions needed in outvars
89 : public :: ebands_from_file ! Build an ebands_t object from file. Supports Fortran and netcdf files
90 : public :: crystal_from_file ! Build a crystal_t object from netcdf or Fortran file with Header
91 : public :: get_gemm_nonlop_ompgpu_blocksize ! Set OpenMP GPU GEMM nonlop blocksize
92 :
93 : !!***
94 :
95 : contains
96 : !!***
97 :
98 : !!****f* ABINIT/scprqt
99 : !! NAME
100 : !! scprqt
101 : !!
102 : !! FUNCTION
103 : !! Conducts printing inside the scfcv.F90 routine, according to the value of choice.
104 : !! Also checks the convergence with respect to the different criteria.
105 : !! Eventually send a signal to quit the SCF cycle.
106 : !!
107 : !! INPUTS
108 : !! choice= if 1 => called at the initialisation of scfcv.f
109 : !! if 2 => called during the loop in scfcv.f
110 : !! if 3 => called at the end of scfcv.f
111 : !! cpus=cpu time limit in seconds
112 : !! deltae=change in energy between the previous and present SCF cycle
113 : !! diffor=maximum absolute change in component of fcart between present and previous SCF cycle.
114 : !! difmagsph=maximum absolute change in magnetization between present and previous SCF cycle.
115 : !! dtset <type(dataset_type)>=all input variables in this dataset
116 : !! | chkexit= if non-zero, check whether the user wishes to exit
117 : !! | enunit=parameter determining units of output energies
118 : !! | ionmov=governs the movement of atoms (see help file)
119 : !! | kptopt=option for the generation of k points
120 : !! | mband=maximum number of bands
121 : !! | natom=number of atoms in cell.
122 : !! | nnsclo_now=number of non-self-consistent loops for the current vtrial
123 : !! | (often 1 for SCF calculation, =nstep for non-SCF calculations)
124 : !! | nsppol=1 for unpolarized, 2 for spin-polarized
125 : !! | occopt=option for occupancies
126 : !! | prtxml=1 if values have to be stored in an XML file.
127 : !! | prteig=
128 : !! | prtstm=print STM input variable
129 : !! | prtvol= control print volume
130 : !! | usedmatpu=DFT+U: number of SCF steps keeping occ. matrix fixed
131 : !! | usefock=1 if Fock operator is present (hence possibility of a double loop)
132 : !! | usepawu=0 if no DFT+U; /=0 if DFT+U
133 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
134 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
135 : !! etotal=total energy (hartree)
136 : !! favg(3)=average of forces (ha/bohr)
137 : !! fcart(3,natom)=cartesian forces (hartree/bohr)
138 : !! fermie=fermi energy (Hartree) / for electrons thermalized in the conduction bands when occopt==9
139 : !! fermih=fermi energy (Hartree) for holes thermalized in the VB when occopt==9
140 : !! fname_eig=filename for printing of the eigenenergies
141 : !! fock <type(fock_type)>=quantities for the fock operator (optional argument)
142 : !! character(len=fnlen) :: filnam1=character strings giving input file name
143 : !! initGS= 1 if one GS SCF cycle has already be done
144 : !! iscf=( <= 0 =>non-SCF), >0 => SCF)
145 : !! iscf =1 => determination of the largest eigenvalue of the SCF cycle
146 : !! iscf =2 => SCF cycle, simple mixing
147 : !! iscf =3 => SCF cycle, anderson mixing
148 : !! iscf =5 => SCF cycle, CG based on estimations of gradients of the energy
149 : !! iscf =6 => SCF cycle, CG based on true minimization of the energy
150 : !! iscf =-3, although non-SCF, the energy is computed, so print it here.
151 : !! istep=number of the SCF iteration (needed if choice=2)
152 : !! istep_fock_outer=number of outer SCF iteration in the double loop approach
153 : !! istep_mix=number of inner SCF iteration in the double loop approach
154 : !! kpt(3,nkpt)=reduced coordinates of k points.
155 : !! maxfor=maximum absolute value of fcart
156 : !! maxmagsph=maximum absolute value of magnetization among all atoms.
157 : !! moved_atm_inside: if==1, the atoms are allowed to move.
158 : !! mpi_enreg=information about MPI parallelization
159 : !! nband(nkpt*nsppol)=number of bands at each k point, for each polarization
160 : !! nkpt=number of k points
161 : !! nstep=number of steps expected in iterations.
162 : !! occ(mband*nkpt*nsppol)=occupation number for each band at each k point.
163 : !! optres=0 if the residual (res2) is a POTENTIAL residual
164 : !! 1 if the residual (res2) is a DENSITY residual
165 : !! prtfor=1 only if forces have to be printed (0 otherwise)
166 : !! prtxml=1 if XML file has to be output
167 : !! res2=square of the density/potential residual
168 : !! resid(mband*nkpt*nsppol)=residuals for each band over all k points and spins
169 : !! residm=maximum value from resid array (except for nbdbuf highest bands)
170 : !! in Wavelets mode, it is used as the maximum value for the gradient norm.
171 : !! response= if 0, GS case, if 1, RF case.
172 : !! tollist(12)=tolerance list. Presently, the following are defined :
173 : !! tollist(1)=tolmxf ; tollist(2)=tolwfr ; tollist(3)=toldff
174 : !! tollist(4)=toldfe ; tollist(5)=toleig ; tollist(6)=tolvrs
175 : !! tollist(7)=tolrff ; tollist(9)=toldmag
176 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
177 : !! vxcavg=mean of the vxc potential
178 : !! wtk(nkpt)=weight assigned to each k point.
179 : !! xred(3,natom)=reduced dimensionless atomic coordinates
180 : !!
181 : !! OUTPUT
182 : !! quit= 0 if the SCF cycle is not finished; 1 otherwise.
183 : !! conv_retcode=Only if choice==3, != 0 if convergence is not achieved.
184 : !!
185 : !! SOURCE
186 :
187 83892 : subroutine scprqt(choice,cpus,deltae,diffor,maxmagsph,difmagsph,dtset,&
188 83892 : & eigen,etotal,favg,fcart,fermie,fermih,fname_eig,filnam1,initGS,&
189 83892 : & iscf,istep,istep_fock_outer,istep_mix,kpt,maxfor,moved_atm_inside,mpi_enreg,&
190 83892 : & nband,nkpt,nstep,occ,optres,&
191 83892 : & prtfor,prtxml,quit,res2,resid,residm,response,tollist,usepaw,&
192 83892 : & vxcavg,wtk,xred,conv_retcode,&
193 : & electronpositron, fock,rcpaw) ! optional arguments)
194 :
195 : !Arguments ------------------------------------
196 : !scalars
197 : integer,intent(in) :: choice,initGS,iscf,istep,istep_fock_outer,istep_mix
198 : integer,intent(in) :: moved_atm_inside,nkpt,nstep
199 : integer,intent(in) :: optres,prtfor,prtxml,response,usepaw
200 : integer,intent(out) :: quit,conv_retcode
201 : real(dp),intent(in) :: cpus,deltae,diffor,etotal,fermie,fermih,maxfor,res2,residm,maxmagsph,difmagsph
202 : real(dp),intent(in) :: vxcavg
203 : character(len=fnlen),intent(in) :: fname_eig,filnam1
204 : type(electronpositron_type),pointer,optional :: electronpositron
205 : type(fock_type),pointer,optional :: fock
206 : type(MPI_type),intent(in) :: mpi_enreg
207 : type(dataset_type),intent(in) :: dtset
208 : type(rcpaw_type),intent(inout),optional,pointer :: rcpaw
209 : !arrays
210 : integer,intent(in) :: nband(nkpt*dtset%nsppol)
211 : real(dp),intent(in) :: eigen(dtset%mband*nkpt*dtset%nsppol),favg(3)
212 : real(dp),intent(in) :: fcart(3,dtset%natom),kpt(3,nkpt)
213 : real(dp),intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol)
214 : real(dp),intent(in) :: resid(dtset%mband*nkpt*dtset%nsppol),tollist(12)
215 : real(dp),intent(in) :: wtk(nkpt),xred(3,dtset%natom)
216 :
217 : !Local variables-------------------------------
218 : !scalars
219 : integer,parameter :: master=0
220 : !integer,save :: toldfe_ok,toldff_ok,tolrff_ok,ttoldfe,ttoldff,ttolrff,ttolvrs,ttolwfr
221 : integer,save :: toldfe_ok,toldff_ok,tolrff_ok,toldmag_ok,ttoldfe,ttoldff,ttolrff,ttolvrs,ttolwfr,ttoldmag
222 : integer :: iatom,iband,iexit,ikpt,ii,ishift,isppol,my_rank
223 : integer :: nband_index,nband_k,nnsclohf
224 : integer :: openexit,option,tmagnet,usefock
225 : #if defined DEV_YP_VDWXC
226 : integer :: ivdw
227 : #endif
228 : !real(dp),save :: toldfe,toldff,tolrff,tolvrs,tolwfr,vdw_df_threshold
229 : real(dp),save :: toldfe,toldff,tolrff,tolvrs,tolwfr,vdw_df_threshold, toldmag
230 : real(dp) :: diff_e,diff_f,magnet,rhodn,rhoup
231 : logical :: noquit,use_dpfft
232 : character(len=500) :: message, message2, message3
233 : character(len=2) :: format_istep
234 : !character(len=5) :: format_magnet
235 : character(len=8) :: colname
236 : character(len=1) :: firstchar
237 : type(yamldoc_t) :: ydoc
238 : !arrays
239 167784 : real(dp) :: residm_band(dtset%mband,dtset%nsppol), f_tmp(3)
240 :
241 : ! *********************************************************************
242 :
243 : DBG_ENTER("COLL")
244 :
245 83892 : my_rank = mpi_enreg%me_cell
246 :
247 83892 : quit=0; conv_retcode=0
248 83892 : usefock=dtset%usefock
249 83892 : nnsclohf=dtset%nnsclohf
250 83892 : use_dpfft = .False.
251 :
252 83892 : tmagnet=0
253 60648 : if(response==0.and.(iscf>0.or.iscf==-3).and.dtset%nsppol==2.and.dtset%occopt>2) tmagnet=1
254 83892 : if((iscf>0.or.iscf==-3).and.(dtset%nsppol==2.or.dtset%nspden>1)) tmagnet=2
255 :
256 83892 : ishift=0
257 1013037 : residm_band = zero
258 175634 : do isppol=1, dtset%nsppol
259 2072859 : do ikpt=1, nkpt
260 12782935 : do iband=1, nband(ikpt+(isppol-1)*nkpt)
261 10793968 : ishift = ishift+1
262 12691193 : residm_band(iband, isppol) = max (resid(ishift), residm_band(iband, isppol))
263 : end do
264 : end do
265 : end do
266 :
267 93152 : select case (choice)
268 : case (1)
269 : ! choice= if 1 => called at the initialisation of scfcv.f
270 : ! Examine tolerance criteria
271 : ! NB: The tests on tolwfr and the presence of tolerances in the SCF case are
272 : ! also done at the level of the parser in chkinp.
273 9260 : tolwfr=tollist(2)
274 9260 : toldff=tollist(3)
275 9260 : toldfe=tollist(4)
276 9260 : tolvrs=tollist(6)
277 9260 : tolrff=tollist(7)
278 9260 : vdw_df_threshold=tollist(8)
279 9260 : toldmag=tollist(9)
280 9260 : ttolwfr=0 ; ttoldff=0 ; ttoldfe=0 ; ttolvrs=0; ttolrff=0; ttoldmag=0 ;
281 9260 : if(abs(tolwfr)>tiny(zero))ttolwfr=1
282 9260 : if(abs(toldff)>tiny(zero))ttoldff=1
283 9260 : if(abs(tolrff)>tiny(zero))ttolrff=1
284 9260 : if(abs(toldfe)>tiny(zero))ttoldfe=1
285 9260 : if(abs(tolvrs)>tiny(zero))ttolvrs=1
286 9260 : if(abs(toldmag)>tiny(zero))ttoldmag=1
287 : ! If non-scf calculations, tolwfr must be defined
288 9260 : if(ttolwfr /= 1 .and. (iscf<0 .and. iscf/=-3) )then
289 : write(message,'(a,a,a,es14.6,a,a)')&
290 0 : 'when iscf <0 and /= -3, tolwfr must be strictly',ch10,&
291 0 : 'positive, while it is ',tolwfr,ch10,&
292 0 : 'Action: change tolwfr in your input file and resubmit the job.'
293 0 : ABI_ERROR(message)
294 : end if
295 : ! toldff only allowed when prtfor==1
296 : ! FIXME: this test should be done on input, not during calculation
297 9260 : if((ttoldff == 1 .or. ttolrff == 1) .and. prtfor==0 )then
298 0 : ABI_ERROR('toldff only allowed when prtfor=1!')
299 : end if
300 : ! If SCF calculations, one and only one of these can differ from zero
301 : !if( (iscf>0 .or. iscf==-3) .and.(ttolwfr==1.and.ttoldff+ttoldfe+ttolvrs+ttolrff>1) &
302 : if( (iscf>0 .or. iscf==-3) .and.(ttolwfr==1.and.ttoldff+ttoldfe+ttolvrs+ttolrff+ttoldmag>1) &
303 : .and. (ttolwfr==0.and.ttoldff+ttoldfe+ttolvrs+ttolrff+ttoldmag/=1) ) then
304 : !.and. (ttolwfr==0.and.ttoldff+ttoldfe+ttolvrs+ttolrff/=1) ) then
305 : write(message,'(6a,es14.6,a,es14.6,a,es14.6,a,a,es14.6,a,a,a)' )&
306 : & 'For the SCF case, one and only one of the input tolerance criteria ',ch10,&
307 : & 'toldff, tolrff, toldfe, toldmag or tolvrs ','must differ from zero, while they are',ch10,&
308 : & 'toldff=',toldff,', tolrff=',tolrff,', toldfe=',toldfe,ch10,&
309 : !& 'and tolvrs=',tolvrs,' .',ch10,&
310 : & 'toldmag=',toldmag,' and tolvrs=',tolvrs,' .',ch10,&
311 : & 'Action: change your input file and resubmit the job.'
312 : ABI_ERROR(message)
313 : end if
314 :
315 9260 : if (dtset%usewvl == 1) then
316 0 : write(colname, "(A)") "grdnorm "
317 : else
318 9260 : write(colname, "(A)") "residm "
319 : end if
320 9260 : if (nstep>0 .and. (iscf>=0 .or.iscf==-3) .and. dtset%prtstm==0) then
321 8622 : if(tmagnet==2 .and. response==0)then
322 779 : if (prtfor==0) then
323 667 : if (optres==0) then
324 455 : write(message, '(4a)' ) ch10,&
325 910 : ' iter Etot(hartree) deltaE(h) ',colname, ' vres2 difmagsph maxmagsph'
326 : else
327 212 : write(message, '(4a)' ) ch10,&
328 424 : ' iter Etot(hartree) deltaE(h) ',colname, ' nres2 difmagsph maxmagsph'
329 : end if
330 : else
331 112 : if (optres==0) then
332 83 : write(message, '(4a)' ) ch10,&
333 166 : ' iter Etot(hartree) deltaE(h) ',colname, ' vres2 diffor maxfor difmagsph maxmagsph'
334 : else
335 29 : write(message, '(4a)' ) ch10,&
336 58 : ' iter Etot(hartree) deltaE(h) ',colname, ' nres2 diffor maxfor difmagsph maxmagsph'
337 : end if
338 : end if
339 : else
340 7843 : if(response==0)then
341 5625 : if (prtfor==0) then
342 2709 : if (optres==0) then
343 2110 : write(message, '(4a)' ) ch10,&
344 4220 : ' iter Etot(hartree) deltaE(h) ', colname, ' vres2'
345 : else
346 599 : write(message, '(4a)' ) ch10,&
347 1198 : ' iter Etot(hartree) deltaE(h) ', colname, ' nres2'
348 : end if
349 : else
350 2916 : if (optres==0) then
351 2349 : write(message, '(4a)' ) ch10,&
352 4698 : ' iter Etot(hartree) deltaE(h) ',colname,' vres2 diffor maxfor '
353 : else
354 567 : write(message, '(4a)' ) ch10,&
355 1134 : ' iter Etot(hartree) deltaE(h) ',colname,' nres2 diffor maxfor '
356 : end if
357 : end if
358 : else
359 2218 : if(tmagnet==2 .and. iscf>0 )then
360 108 : if (optres==0) then
361 108 : write(message, '(4a)' ) ch10,&
362 216 : ' iter 2DEtotal(Ha) deltaE(Ha) ', colname, ' vres2 difmagsph maxmagsph'
363 : else
364 0 : write(message, '(4a)' ) ch10,&
365 0 : ' iter 2DEtotal(Ha) deltaE(Ha) ', colname, ' nres2 difmagsph maxmagsph'
366 : end if
367 : else
368 2110 : if (optres==0) then
369 2110 : write(message, '(4a)' ) ch10,&
370 4220 : ' iter 2DEtotal(Ha) deltaE(Ha) ', colname, ' vres2 '
371 : else
372 0 : write(message, '(4a)' ) ch10,&
373 0 : ' iter 2DEtotal(Ha) deltaE(Ha) ', colname, ' nres2 '
374 : end if
375 : end if
376 : end if
377 : end if
378 :
379 8622 : ydoc = yamldoc_open('BeginCycle')
380 : ! If wfoptalg=1 or 111 or 112, we should write mdeg_filter
381 : call ydoc%add_ints("iscf, nstep, nline, wfoptalg", &
382 43110 : [dtset%iscf, dtset%nstep, dtset%nline, dtset%wfoptalg], dict_key="solver")
383 : call ydoc%add_reals("tolwfr, toldff, toldfe, tolvrs, tolrff, toldmag", & ! , vdw_df_threshold", &
384 : [tolwfr, toldff, toldfe, tolvrs, tolrff, toldmag], & !, vdw_df_threshold], &
385 : !call ydoc%add_reals("tolwfr, toldff, toldfe, tolvrs, tolrff", & ! , vdw_df_threshold", &
386 : ! [tolwfr, toldff, toldfe, tolvrs, tolrff], & !, vdw_df_threshold], &
387 60354 : real_fmt="(es8.2)", dict_key="tolerances", ignore=zero)
388 :
389 8622 : call ydoc%write_and_free(ab_out, newline=.False.)
390 8622 : call wrtout(ab_out, message)
391 : end if
392 :
393 : case (2)
394 :
395 : ! Examine tolerance criteria
396 65372 : tolwfr=tollist(2)
397 65372 : toldff=tollist(3)
398 65372 : toldfe=tollist(4)
399 65372 : tolvrs=tollist(6)
400 65372 : tolrff=tollist(7)
401 65372 : vdw_df_threshold=tollist(8)
402 65372 : toldmag=tollist(9)
403 65372 : ttolwfr=0 ; ttoldff=0 ; ttoldfe=0 ; ttolvrs=0; ttolrff=0; ttoldmag=0;
404 : !ttolwfr=0 ; ttoldff=0 ; ttoldfe=0 ; ttolvrs=0; ttolrff=0;
405 65372 : if(abs(tolwfr)>tiny(0.0_dp))ttolwfr=1
406 65372 : if(abs(toldff)>tiny(0.0_dp))ttoldff=1
407 65372 : if(abs(tolrff)>tiny(0.0_dp))ttolrff=1
408 65372 : if(abs(toldfe)>tiny(0.0_dp))ttoldfe=1
409 65372 : if(abs(tolvrs)>tiny(0.0_dp))ttolvrs=1
410 65372 : if(abs(toldmag)>tiny(0.0_dp))ttoldmag=1
411 :
412 : ! Conduct printing. If extra output follows, then put a blank line into the output here
413 73468 : if (dtset%prtvol>=10) call wrtout([std_out, ab_out], ' ')
414 :
415 : ! Calculate up and down charge and magnetization
416 65372 : if(tmagnet==1 ) then
417 0 : rhoup = zero
418 0 : rhodn = zero
419 0 : nband_index = 1
420 0 : do isppol=1,dtset%nsppol
421 0 : do ikpt=1,nkpt
422 0 : nband_k=nband(ikpt+(isppol-1)*nkpt)
423 0 : do iband=1,nband_k
424 0 : if(isppol==1) rhoup = rhoup + wtk(ikpt)*occ(nband_index)
425 0 : if(isppol==2) rhodn = rhodn + wtk(ikpt)*occ(nband_index)
426 0 : nband_index = nband_index + 1
427 : end do
428 : end do
429 : end do
430 0 : magnet = abs(rhoup - rhodn)
431 : end if
432 :
433 65372 : if (prtxml == 1) then
434 26 : write(ab_xml_out, "(A)", advance = "NO") ' <scfcvStep'
435 26 : write(message, "(es22.10)") etotal
436 26 : write(ab_xml_out, "(A,A,A)", advance = "NO") ' eTotal="', trim(message) ,'"'
437 26 : write(message, "(es20.8)") deltae
438 26 : write(ab_xml_out, "(A,A,A)", advance = "NO") ' deltaETotal="', trim(message) ,'"'
439 26 : write(message, "(es20.8)") residm
440 26 : write(ab_xml_out, "(A,A,A)", advance = "NO") ' maxResid="', trim(message) ,'"'
441 26 : write(message, "(es20.8)") res2
442 26 : if (optres == 0) then
443 26 : write(ab_xml_out, "(A,A,A)", advance = "NO") ' potResid="', trim(message) ,'"'
444 : else
445 0 : write(ab_xml_out, "(A,A,A)", advance = "NO") ' denResid="', trim(message) ,'"'
446 : end if
447 26 : if (tmagnet== 1) then
448 0 : write(message, "(es20.8)") magnet
449 0 : write(ab_xml_out, "(A,A,A)", advance = "NO") ' magn="', trim(message) ,'"'
450 : end if
451 26 : if (prtfor == 1) then
452 0 : write(message, "(es20.8)") diffor
453 0 : write(ab_xml_out, "(A,A,A)", advance = "NO") ' deltaForces="', trim(message) ,'"'
454 0 : write(message, "(es20.8)") maxfor
455 0 : write(ab_xml_out, "(A,A,A)", advance = "NO") ' maxForces="', trim(message) ,'"'
456 : end if
457 26 : write(ab_xml_out, "(A)") " />"
458 : end if
459 :
460 : ! Print total (free) energy (hartree) and other convergence measures
461 65372 : if(dtset%prtstm==0)then
462 65369 : format_istep='i3'
463 65369 : if(istep>99)format_istep='i5'
464 65369 : if(istep>9999)format_istep='i7'
465 : !if(tmagnet==1)then
466 65369 : if(tmagnet==2 .and. response==0)then
467 : ! if(magnet<10)then
468 : ! format_magnet='f6.3)'
469 : ! else if(magnet<100)then
470 : ! format_magnet='f6.2)'
471 : ! else
472 : ! format_magnet='f6.1)'
473 : ! end if
474 7534 : if (prtfor==0) then
475 : !write(message, '(a,'//format_istep//',1p,g22.14,3es9.2,0p,'//format_magnet ) &
476 : write(message, '(a,'//format_istep//',1p,g22.14,5es10.3)' ) &
477 : !' ETOT',istep,etotal,deltae,residm,res2,magnet
478 6720 : ' ETOT',istep,etotal,deltae,residm,res2,difmagsph,maxmagsph
479 : else
480 : write(message, '(a,'//format_istep//',1p,g22.14,7es10.3)' ) &
481 814 : ' ETOT',istep,etotal,deltae,residm,res2,diffor,maxfor,difmagsph,maxmagsph
482 : end if
483 : else
484 57835 : firstchar=' '
485 57835 : if (response/=0.and.istep==1) firstchar="-"
486 57835 : if (response==0) then
487 39275 : if (prtfor==0) then
488 : write(message, '(2a,'//format_istep//',1p,g22.14,3es10.3)' ) &
489 19597 : firstchar,'ETOT',istep,etotal,deltae,residm,res2
490 : else
491 : write(message, '(2a,'//format_istep//',1p,g22.14,5es10.3)' ) &
492 19678 : firstchar,'ETOT',istep,etotal,deltae,residm,res2,diffor,maxfor
493 : end if
494 : else
495 18560 : if(tmagnet==2 .and. iscf>0 )then
496 : write(message, '(2a,'//format_istep//',1p,g22.14,1x,5es10.3)' ) &
497 1568 : firstchar,'ETOT',istep,etotal,deltae,residm,res2,difmagsph,maxmagsph
498 : else
499 : write(message, '(2a,'//format_istep//',1p,g22.14,1x,3es10.3)' ) &
500 16992 : firstchar,'ETOT',istep,etotal,deltae,residm,res2
501 : end if
502 : end if
503 : end if
504 : !if (etot_yaml_doc%stream%length /= 0) call etot_yaml_doc%add_tabular_line(' '//message(6:))
505 65369 : call wrtout(ab_out,message)
506 :
507 65369 : if(mpi_enreg%paral_pert==1) then
508 337 : call wrtout(std_out, message,'PERS')
509 65032 : elseif(mpi_enreg%paral_pert==0) then
510 65032 : call wrtout(std_out, message,'COLL')
511 : end if
512 :
513 : end if ! dtset%prtstm==0
514 :
515 : ! Print positions/forces every step if dtset%prtvol>=10 and iscf>0 or -3 and GS case
516 65372 : if (dtset%prtvol>=10.and.(iscf>=0.or.iscf==-3).and.response==0.and.dtset%prtstm==0) then
517 1657 : call wrtout(ab_out," ")
518 :
519 : ! Print up and down charge and magnetization
520 1657 : if(tmagnet==1) then
521 : write(message,'(a,f11.6,a,f11.6,a,f10.6)')&
522 0 : ' #electrons spin up=',rhoup,', spin down=',rhodn,', magnetization=',magnet
523 0 : call wrtout([std_out, ab_out], message)
524 : end if
525 :
526 : ! Moreover, print atomic positions if dtset%ionmov==4, and moved_atm_inside==1
527 1657 : if (dtset%ionmov==4 .and. moved_atm_inside==1)then
528 255 : call wrtout([std_out, ab_out], ' reduced coordinates :')
529 255 : do iatom=1,dtset%natom
530 170 : write(message, '(i5,1x,3es21.11)' ) iatom,xred(:,iatom)
531 595 : call wrtout([std_out, ab_out], message)
532 : end do
533 : end if
534 :
535 : ! Slightly change favg for printing reasons
536 1657 : if (prtfor>0) then
537 529 : f_tmp(:)=favg(:)
538 529 : if(abs(favg(1))<1.0d-13)f_tmp(1)=zero
539 529 : if(abs(favg(2))<1.0d-13)f_tmp(2)=zero
540 529 : if(abs(favg(3))<1.0d-13)f_tmp(3)=zero
541 529 : write(message, '(a,3es10.2)' )' cartesian forces (ha/bohr); non-corrected avg=',f_tmp(:)
542 1587 : call wrtout([std_out, ab_out], message)
543 2599 : do iatom=1,dtset%natom
544 8280 : f_tmp(:)=fcart(:,iatom)
545 2070 : if(abs(fcart(1,iatom))<1.0d-13)f_tmp(1)=zero
546 2070 : if(abs(fcart(2,iatom))<1.0d-13)f_tmp(2)=zero
547 2070 : if(abs(fcart(3,iatom))<1.0d-13)f_tmp(3)=zero
548 2070 : write(message, '(i5,1x,3es21.11)' ) iatom,f_tmp(:)
549 6739 : call wrtout([std_out, ab_out], message)
550 : end do
551 : end if
552 :
553 : end if
554 :
555 : ! Print eigenvalues every step if dtset%prtvol>=10 and GS case
556 65372 : if (my_rank == master .and. (dtset%prtvol>=10 .and. response==0 .and. dtset%tfkinfunc==0 .and. dtset%usewvl==0)) then
557 1525 : option=1
558 : call prteigrs(eigen,dtset%enunit,fermie,fermih,fname_eig,ab_out,iscf,kpt,dtset%kptopt,dtset%mband,&
559 1525 : nband,dtset%nbdbuf,nkpt,dtset%nnsclo,dtset%nsppol,occ,dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk)
560 : call prteigrs(eigen,dtset%enunit,fermie,fermih,fname_eig,std_out,iscf,kpt,dtset%kptopt,dtset%mband,&
561 1525 : nband,dtset%nbdbuf,nkpt,dtset%nnsclo,dtset%nsppol,occ,dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk)
562 : end if
563 :
564 65372 : if(response==0)then
565 46812 : write(message, '(a,1p,e15.7,a)' ) ' scprqt: <Vxc>=',vxcavg,' Hartree'
566 46812 : call wrtout(std_out,message)
567 : end if
568 :
569 : ! Check whether exiting was required by the user.
570 65372 : openexit=1 ; if(dtset%chkexit==0) openexit=0
571 65372 : call exit_check(cpus,filnam1,iexit,ab_out,mpi_enreg%comm_cell,openexit)
572 65372 : if (iexit/=0) quit=1
573 :
574 : ! In special cases, do not quit even if convergence is reached
575 : noquit=((istep<nstep).and.(usepaw==1).and.(dtset%usepawu/=0).and.&
576 : (dtset%usedmatpu/=0).and.(istep<=abs(dtset%usedmatpu)).and.&
577 65372 : (dtset%usedmatpu<0.or.initGS==0))
578 :
579 : ! Additional stuff for electron/positron
580 65372 : if (present(electronpositron)) then
581 46812 : if (associated(electronpositron)) then
582 562 : if (electronpositron%istep_scf==1) then
583 79 : toldff_ok=0;tolrff_ok=0;toldfe_ok=0; toldmag_ok=0;
584 : end if
585 : end if
586 : end if
587 :
588 : ! Stopping criteria in the SCF case
589 65372 : if(iscf>1 .or. iscf==-3 .or. iscf == 0) then
590 : ! Here treat the vdw_df_threshold criterion : if the change of energy is less than
591 : ! input vdw_df_threshold, trigger the calculation of vdW interactions
592 : ! write(message,'(1x,a,e10.3,1x,a,e10.3,1x,l1,a)') &
593 : ! & '[vdW-DF][DEBUG] deltae=',deltae,'vdw_df_threshold=',vdw_df_threshold, &
594 : ! & (abs(deltae)<vdw_df_threshold),ch10
595 : ! call wrtout(std_out,message)
596 : #if defined DEV_YP_VDWXC
597 : call xc_vdw_trigger( (abs(deltae)<vdw_df_threshold) )
598 : #endif
599 : ! Here treat the tolwfr criterion: if maximum residual is less than
600 : ! input tolwfr, stop steps (exit loop here)
601 65325 : if (ttolwfr == 1 .and. (ttolvrs+ttoldfe+ttoldff+ttolrff+ttoldmag==0) .and. .not. noquit) then
602 : !if (ttolwfr == 1 .and. (ttolvrs+ttoldfe+ttoldff+ttolrff==0) .and. .not. noquit) then
603 15165 : if (residm < tolwfr) then
604 1633 : if (dtset%usewvl == 0) then
605 1633 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a)' )ch10, &
606 3266 : ' At SCF step',istep,' max residual=',residm,' < tolwfr=',tolwfr,' =>converged.'
607 : else
608 0 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a)' )ch10, &
609 0 : ' At SCF step',istep,' max grdnorm=',residm,' < tolwfr=',tolwfr,' =>converged.'
610 : end if
611 4899 : call wrtout([std_out, ab_out], message)
612 1633 : quit=1
613 : else
614 13532 : use_dpfft = residm < tol7
615 : end if
616 15165 : if(present(rcpaw)) then
617 917 : if (associated(rcpaw).and.residm < (tolwfr)**third) then
618 0 : if(rcpaw%updatepaw(1)==0.and.rcpaw%updatepaw(2)==0) then
619 0 : rcpaw%updatepaw(:)=istep
620 : endif
621 : endif
622 : endif
623 : end if
624 :
625 : ! Here treat the toldff criterion: if maximum change of fcart is less than
626 : ! input toldff twice consecutively, stop steps (exit loop here)
627 65325 : if (ttoldff==1) then
628 16395 : if (istep==1) then
629 2673 : toldff_ok=0
630 13722 : else if (diffor < toldff) then
631 4354 : toldff_ok=toldff_ok+1
632 : ! add warning for forces which are 0 by symmetry. Also added Matteo check below that the wave
633 : ! functions are relatively converged as well
634 4354 : if (diffor < tol12) then
635 1 : write (message,'(3a)') ' toldff criterion is satisfied, but your forces are suspiciously low.', ch10,&
636 2 : ' Check if the forces are 0 by symmetry: in that case you can not use the toldff convergence criterion!'
637 1 : ABI_WARNING(message)
638 1 : if (maxfor < tol16 .and. res2 > tol9) tolrff_ok=0
639 : end if
640 : else
641 9368 : toldff_ok=0
642 9368 : use_dpfft = diffor < tol6
643 : end if
644 16395 : if(present(rcpaw)) then
645 3147 : if (associated(rcpaw).and.diffor < (toldff)**third) then
646 0 : if(rcpaw%updatepaw(1)==0.and.rcpaw%updatepaw(2)==0) then
647 0 : rcpaw%updatepaw(:)=istep
648 : endif
649 : endif
650 : endif
651 :
652 16395 : if(toldff_ok>=2 .and..not.noquit)then
653 2141 : if (ttolwfr==0) then
654 2135 : write(message, '(a,a,i5,a,a,a,es11.3,a,es11.3)' ) ch10, &
655 2135 : ' At SCF step',istep,', forces are converged : ',ch10,&
656 4270 : ' for the second time, max diff in force=',diffor,' < toldff=',toldff
657 6405 : call wrtout([std_out, ab_out], message)
658 2135 : quit=1
659 6 : else if (ttolwfr==1 .and. residm < tolwfr )then
660 2 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a,a,es11.3,a,es11.3)' ) ch10, &
661 2 : ' At SCF step',istep,', max residual=',residm,' < tolwfr=',tolwfr,' AND forces are converged : ',ch10,&
662 4 : ' for the second time, max diff in force=',diffor,' < toldff=',toldff
663 6 : call wrtout([std_out, ab_out], message)
664 2 : quit=1
665 : end if
666 : end if
667 : end if
668 :
669 : ! HEREEEE
670 65325 : if (ttoldmag==1) then
671 30 : if (istep==1) then
672 2 : toldmag_ok=0
673 28 : else if (difmagsph < toldmag) then
674 4 : toldmag_ok=toldmag_ok+1
675 : ! add warning for magnetizations which are 0 . Check below that the wave
676 : ! functions are relatively converged as well
677 4 : if (maxmagsph < tol8) then
678 0 : write (message,'(3a)') ' toldmag criterion is satisfied, but your magnetizations are suspiciously low.', ch10,&
679 0 : ' Check if the magnetizations are 0 : in that case you can not use the toldmag convergence criterion!'
680 0 : ABI_WARNING(message)
681 : ! HEREE ??
682 : end if
683 4 : if (res2 > tol9) toldmag_ok=0
684 : else
685 24 : toldmag_ok=0
686 24 : use_dpfft = difmagsph < tol6
687 : end if
688 :
689 30 : if(toldmag_ok>=2 .and..not.noquit)then
690 2 : if (ttolwfr==0) then
691 2 : write(message, '(a,a,i5,a,a,a,es11.3,a,es11.3)' ) ch10, &
692 2 : ' At SCF step',istep,', magnetizations are converged : ',ch10,&
693 4 : ' for the second time, max diff in magnetization=',difmagsph,' < toldmag=',toldmag
694 6 : call wrtout([std_out, ab_out], message)
695 2 : quit=1
696 0 : else if (ttolwfr==1 .and. residm < tolwfr )then
697 0 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a,a,es11.3,a,es11.3)' ) ch10, &
698 0 : ' At SCF step',istep,', max residual=',residm,' < tolwfr=',tolwfr,' AND magnetizations are converged : ',ch10,&
699 0 : ' for the second time, max diff in magnetization=',difmagsph,' < toldmag=',toldmag
700 0 : call wrtout([std_out, ab_out], message)
701 0 : quit=1
702 : end if
703 : end if
704 : end if
705 :
706 : ! Here treat the tolrff criterion: if maximum change of fcart is less than
707 : ! input tolrff times fcart itself twice consecutively, stop steps (exit loop here)
708 65325 : if (ttolrff==1) then
709 867 : if (istep==1) then
710 132 : tolrff_ok=0
711 : ! 27/7/2009: added test for absolute value of maxfor, otherwise if it is 0 this never exits the scf loop.
712 735 : else if (diffor < tolrff*maxfor .or. (maxfor < tol16 .and. diffor < tol16)) then
713 255 : tolrff_ok=tolrff_ok+1
714 : ! Thu Mar 12 19:01:40 MG: added additional check on res2 to make sure the SCF cycle is close to convergence.
715 : ! Needed for structural relaxations otherwise the stress tensor is wrong and the relax algo makes wrong moves.
716 255 : if (maxfor < tol16 .and. res2 > tol9) tolrff_ok=0
717 : else
718 480 : tolrff_ok=0
719 480 : use_dpfft = diffor < tolrff * maxfor * five
720 : end if
721 867 : if(present(rcpaw)) then
722 134 : if (associated(rcpaw).and.(diffor < ( tolrff*maxfor)**third.or.(maxfor < tol6 .and. diffor < tol6))) then
723 0 : if(rcpaw%updatepaw(1)==0.and.rcpaw%updatepaw(2)==0) then
724 0 : rcpaw%updatepaw(:)=istep
725 : endif
726 : endif
727 : endif
728 867 : if(tolrff_ok>=2 .and. (.not.noquit))then
729 119 : if (ttolwfr==0) then
730 113 : write(message, '(a,a,i5,a,a,a,es11.3,a,es11.3,a)' ) ch10, &
731 113 : ' At SCF step',istep,', forces are sufficiently converged : ',ch10,&
732 113 : ' for the second time, max diff in force=',diffor,&
733 226 : ' is less than < tolrff=',tolrff, ' times max force'
734 339 : call wrtout([std_out, ab_out], message)
735 113 : quit=1
736 6 : else if (ttolwfr==1 .and. residm < tolwfr) then
737 2 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a,a,es11.3,a,es11.3,a)' ) ch10, &
738 2 : ' At SCF step',istep,', max residual=',residm,' < tolwfr=',tolwfr,' AND forces are sufficiently converged : ',ch10,&
739 2 : ' for the second time, max diff in force=',diffor,&
740 4 : ' is less than < tolrff=',tolrff, ' times max force'
741 6 : call wrtout([std_out, ab_out], message)
742 2 : quit=1
743 : end if
744 : end if
745 : end if
746 :
747 : ! Here treat the toldfe criterion: if the change of energy is less than
748 : ! input toldfe twice consecutively, stop steps (exit loop here)
749 65325 : if (ttoldfe==1) then
750 8893 : if (istep==1) then
751 1263 : toldfe_ok=0
752 7630 : else if (abs(deltae)<toldfe) then
753 2173 : toldfe_ok=toldfe_ok+1
754 : else
755 5457 : toldfe_ok=0
756 5457 : use_dpfft = abs(deltae) < tol8
757 : end if
758 8893 : if(present(rcpaw)) then
759 2483 : if (associated(rcpaw).and.abs(deltae) < (toldfe)**third) then
760 7 : if(rcpaw%updatepaw(1)==0.and.rcpaw%updatepaw(2)==0) then
761 6 : rcpaw%updatepaw(:)=istep
762 : endif
763 : endif
764 : endif
765 : ! Fock : tolwfr not taken into account
766 8893 : if(usefock/=0.and.nnsclohf>=2) then
767 197 : if (toldfe_ok==2 .and. (.not.noquit))then
768 : write(message, '(a,i3,a,i3,a,a,a,es11.3,a,es11.3)' ) &
769 11 : ' Outer loop step',istep_fock_outer,' - inner step',istep_mix,' - frozen Fock etot converged : ',ch10,&
770 22 : ' for the second time, diff in etot=',abs(deltae),' < toldfe=',toldfe
771 33 : call wrtout([std_out, ab_out], message)
772 11 : quit=1
773 : end if
774 : ! No Fock : take into account tolwfr if ttolwfr/=0
775 8696 : else if(toldfe_ok>=2 .and. (.not.noquit))then
776 1060 : if (ttolwfr==0) then
777 1054 : write(message, '(a,a,i5,a,a,a,es11.3,a,es11.3)' ) ch10, &
778 1054 : ' At SCF step',istep,', etot is converged : ',ch10,&
779 2108 : ' for the second time, diff in etot=',abs(deltae),' < toldfe=',toldfe
780 3162 : call wrtout([std_out, ab_out], message)
781 1054 : quit=1
782 6 : else if (ttolwfr==1 .and. residm < tolwfr) then
783 2 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a,a,es11.3,a,es11.3)' ) ch10, &
784 2 : ' At SCF step',istep,', max residual=',residm,' < tolwfr=',tolwfr,' AND etot is converged : ',ch10,&
785 4 : ' for the second time, diff in etot=',abs(deltae),' < toldfe=',toldfe
786 6 : call wrtout([std_out, ab_out], message)
787 2 : quit=1
788 : end if
789 : end if
790 8893 : if(usefock==1 .and. nnsclohf>1)then
791 197 : if(istep_mix==1 .and. (.not.noquit))then
792 : ! The change due to the update of the Fock operator is sufficiently small. No need to meet it a second times.
793 60 : if (abs(deltae)<toldfe) then
794 : write(message, '(a,i3,a,i3,a,a,a,es11.3,a,es11.3)' ) &
795 2 : ' Outer loop step',istep_fock_outer,' - inner step',istep_mix,' - etot converged : ',ch10,&
796 4 : ' update of Fock operator yields diff in etot=',abs(deltae),' < toldfe=',toldfe
797 6 : call wrtout([std_out, ab_out], message)
798 2 : fock%fock_common%fock_converged=.true.
799 2 : quit=1
800 : endif
801 : endif
802 : !TODO: separate messages: if HF is imposing a continuation of the loop, then abs(deltae) is actually not > toldfe
803 197 : if(istep_mix==nnsclohf .and. quit==0)then
804 : write(message, '(a,i3,a,i3,a,a,a,es11.3,a,es11.3)' ) &
805 45 : ' Outer loop step',istep_fock_outer,' - inner step',istep_mix,' - frozen Fock etot NOT converged : ',ch10,&
806 90 : ' diff in etot=',abs(deltae),' > toldfe=',toldfe
807 135 : call wrtout([std_out, ab_out], message)
808 : endif
809 : endif
810 :
811 : ! Here treat the vdw_df_threshold criterion for non-SCF vdW-DF
812 : ! calculations: If input vdw_df_threshold is lesss than toldfe
813 : ! then the vdW-DF is triggered once selfconsistency criteria is
814 : ! reached for the first time.
815 : ! write(message,'(1x,a,e10.3,1x,a,e10.3,1x,l1,a)') &
816 : ! & '[vdW-DF][DEBUG] deltae=',deltae,'vdw_df_threshold=',vdw_df_threshold, &
817 : ! & (abs(deltae)<toldfe),ch10
818 : ! call wrtout(std_out,message)
819 : #if defined DEV_YP_VDWXC
820 : ivdw = 0
821 : if ( toldfe > vdw_df_threshold ) then
822 : ivdw = ivdw + 1
823 : end if
824 : call xc_vdw_trigger((toldfe_ok==1 .and. toldfe>vdw_df_threshold))
825 : if ( ivdw == 2) then
826 : quit=1
827 : end if
828 : #endif
829 : end if
830 :
831 : ! Here treat the tolvrs criterion: if density/potential residual (squared)
832 : ! is less than input tolvrs, stop steps (exit loop here)
833 65325 : if (ttolvrs==1 .and. .not. noquit) then
834 23901 : if (ttolwfr==0) then
835 23652 : if (res2 < tolvrs) then
836 1804 : if (optres==0) then
837 1562 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a)' ) ch10,&
838 3124 : ' At SCF step',istep,' vres2 =',res2,' < tolvrs=',tolvrs,' =>converged.'
839 : else
840 242 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a)' ) ch10,&
841 484 : ' At SCF step',istep,' nres2 =',res2,' < tolvrs=',tolvrs,' =>converged.'
842 : end if
843 5412 : call wrtout([std_out, ab_out], message)
844 1804 : quit=1
845 : else
846 21848 : use_dpfft = res2 < tol5
847 : end if
848 249 : else if (ttolwfr==1 .and. residm < tolwfr) then
849 18 : if (res2 < tolvrs) then
850 14 : if (optres==0) then
851 2 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,e10.2,a,e10.2,a)' ) ch10,&
852 2 : ' At SCF step',istep,' max residual=',residm,' < tolwfr=',tolwfr,' AND vres2 =',res2,&
853 4 : & ' < tolvrs=',tolvrs,' =>converged.'
854 : else
855 12 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,e10.2,a,e10.2,a)' ) ch10,&
856 12 : ' At SCF step',istep,' max residual=',residm,' < tolwfr=',tolwfr,' AND nres2 =',res2,&
857 24 : & ' < tolvrs=',tolvrs,' =>converged.'
858 : end if
859 42 : call wrtout([std_out, ab_out], message)
860 14 : quit=1
861 : else
862 4 : use_dpfft = res2 < tol5
863 : end if
864 : end if
865 23901 : if(present(rcpaw)) then
866 7582 : if (associated(rcpaw).and.res2 < (tolvrs)**third) then
867 0 : if(rcpaw%updatepaw(1)==0.and.rcpaw%updatepaw(2)==0) then
868 0 : rcpaw%updatepaw(:)=istep
869 : endif
870 : endif
871 : endif
872 : end if
873 :
874 65325 : if (quit==1.and.noquit) then
875 0 : write(message, '(a,a,a)' ) ch10, &
876 0 : ' SCF cycle will continue as it is in an initialization stage',' (occ. matrix was kept constant)...'
877 0 : call wrtout([std_out, ab_out], message)
878 : end if
879 :
880 : end if
881 :
882 : ! Activate FFT in double-precision.
883 65325 : if (use_dpfft) then
884 18033 : if (fftcore_mixprec == 1) call wrtout(std_out, " Approaching convergence. Activating FFT in double-precision")
885 18033 : ii = fftcore_set_mixprec(0)
886 : end if
887 :
888 : ! Additional stuff for rcpaw
889 65372 : if(present(rcpaw)) then
890 14329 : if(associated(rcpaw)) then
891 21 : if(rcpaw%updatepaw(2)>0.and.rcpaw%frocc) rcpaw%updateocc=rcpaw%updatepaw(2)
892 21 : if(rcpaw%updatetnc==0)rcpaw%updatetnc=rcpaw%updatepaw(2)
893 21 : if(quit==1.and.(rcpaw%istep<rcpaw%updatepaw(2).or.rcpaw%updatepaw(2)==0))quit=0
894 : endif
895 : endif
896 :
897 : case (3)
898 : ! If wavefunction convergence was not reached (for nstep>0) print a warning and return conv_retcode
899 :
900 : conv_retcode = 0
901 9260 : if(nstep>0) then
902 9225 : if (.not. converged()) then
903 1932 : conv_retcode = 1
904 :
905 1932 : if(iscf>=1 .or. iscf==-3 .or. iscf == 0)then
906 1817 : write(message, '(a,a,a,a,i5,a)' ) ch10,&
907 1817 : ' scprqt: WARNING -',ch10,&
908 3634 : ' nstep=',nstep,' was not enough SCF cycles to converge;'
909 :
910 1817 : write(std_out,'(6a,i0,3a)')ch10,&
911 1817 : "--- !ScfConvergenceWarning",ch10,&
912 1817 : "message: |",ch10,&
913 1817 : ' nstep ',nstep,' was not enough SCF cycles to converge.',ch10,&
914 3634 : "..."
915 : !ABI_WARNING_CLASS(message, "ScfConvergenceWarning")
916 : else
917 115 : write(message, '(a,a,a,a,i5,a)' ) ch10,&
918 115 : ' scprqt: WARNING -',ch10,&
919 230 : ' nstep=',nstep,' was not enough non-SCF iterations to converge;'
920 :
921 115 : write(std_out,'(8a)')ch10,&
922 115 : "--- !NscfConvergenceWarning",ch10,&
923 115 : "message: |",ch10,TRIM(indent(message)),ch10,&
924 230 : "..."
925 : !ABI_WARNING_CLASS(message, "NScfConvergenceWarning")
926 : end if
927 5796 : call wrtout([std_out, ab_out], message)
928 :
929 1932 : if (ttolwfr==1 .and. residm > tolwfr) then
930 681 : if (dtset%usewvl == 0) then
931 : write(message, '(a,es11.3,a,es11.3,a)' ) &
932 681 : ' maximum residual=',residm,' exceeds tolwfr=',tolwfr,ch10
933 :
934 : write(message2, '(a,es11.3,2a)' ) &
935 681 : ' maximum residual each band. tolwfr= ',tolwfr,ch10,&
936 1362 : ' iband, isppol, individual band residuals (max over all k-points):'
937 681 : call wrtout(std_out, message2)
938 1524 : do isppol = 1, dtset%nsppol
939 9281 : do iband = 1, dtset%mband
940 7757 : write(message3, '(2i6, es11.3)') iband, isppol, residm_band(iband,isppol)
941 8600 : call wrtout(std_out,message3)
942 : end do
943 : end do
944 :
945 : else
946 : write(message, '(a,es11.3,a,es11.3,a)' ) &
947 0 : ' maximum grdnorm=',residm,' exceeds tolwfr=',tolwfr,ch10
948 : end if
949 :
950 1251 : else if (ttoldff==1) then
951 : write(message, '(a,es11.3,a,es11.3,a)' ) &
952 460 : ' maximum force difference=',diffor,' exceeds toldff=',toldff,ch10
953 :
954 791 : else if (ttoldmag==1) then
955 : write(message, '(a,es11.3,a,es11.3,a)' ) &
956 0 : ' maximum magnetization difference=',difmagsph,' exceeds toldmag=',toldmag,ch10
957 :
958 791 : else if (ttolrff==1) then
959 : write(message, '(a,es11.3,a,es11.3,a)' ) &
960 11 : ' maximum force difference=',diffor,' exceeds tolrff*maxfor=',tolrff*maxfor,ch10
961 :
962 780 : else if (ttoldfe==1) then
963 : write(message, '(a,es11.3,a,es11.3,a)' ) &
964 250 : ' maximum energy difference=',abs(deltae),' exceeds toldfe=',toldfe,ch10
965 :
966 530 : else if(ttolvrs==1)then
967 530 : if (optres==0) then
968 : write(message, '(a,es11.3,a,es11.3,a)' ) &
969 174 : ' potential residual=',res2,' exceeds tolvrs=',tolvrs,ch10
970 : else
971 : write(message, '(a,es11.3,a,es11.3,a)' ) &
972 356 : ' density residual=',res2,' exceeds tolvrs=',tolvrs,ch10
973 : end if
974 : end if
975 5796 : call wrtout([std_out, ab_out], message)
976 :
977 1932 : if (prtxml == 1) then
978 0 : write(ab_xml_out, "(A)", advance = "NO") ' <status cvState="Failed"'
979 : end if
980 :
981 : else
982 : ! Convergence is OK
983 7293 : if (prtxml == 1) then
984 4 : write(ab_xml_out, "(A)", advance = "NO") ' <status cvState="Ok"'
985 : end if
986 : end if ! test for convergence reached or not
987 :
988 9225 : if (prtxml == 1) then
989 4 : if (ttoldfe == 1) then
990 4 : if (ttolwfr==0) then
991 4 : write(ab_xml_out, "(A)") ' stop-criterion="toldfe" />'
992 : else
993 0 : write(ab_xml_out, "(A)") ' stop-criterion="toldfe+tolwfr" />'
994 : end if
995 0 : else if (ttoldff == 1) then
996 0 : if (ttolwfr==0) then
997 0 : write(ab_xml_out, "(A)") ' stop-criterion="toldff" />'
998 : else
999 0 : write(ab_xml_out, "(A)") ' stop-criterion="toldff+tolwfr" />'
1000 : end if
1001 0 : else if (ttoldmag == 1) then
1002 0 : if (ttolwfr==0) then
1003 0 : write(ab_xml_out, "(A)") ' stop-criterion="toldmag" />'
1004 : else
1005 0 : write(ab_xml_out, "(A)") ' stop-criterion="toldmag+tolwfr" />'
1006 : end if
1007 0 : else if (ttolrff == 1) then
1008 0 : if (ttolwfr==0) then
1009 0 : write(ab_xml_out, "(A)") ' stop-criterion="tolrff" />'
1010 : else
1011 0 : write(ab_xml_out, "(A)") ' stop-criterion="tolrff+tolwfr" />'
1012 : end if
1013 0 : else if (ttolvrs == 1) then
1014 0 : if (ttolwfr==0) then
1015 0 : write(ab_xml_out, "(A)") ' stop-criterion="tolvrs" />'
1016 : else
1017 0 : write(ab_xml_out, "(A)") ' stop-criterion="tolvrs+tolwfr" />'
1018 : end if
1019 0 : else if (ttolwfr == 1) then
1020 0 : write(ab_xml_out, "(A)") ' stop-criterion="tolwfr" />'
1021 : else
1022 0 : write(ab_xml_out, "(A)") ' />'
1023 : end if
1024 : end if
1025 :
1026 : ! If enabled, output a YAML document with the ETOT iterations
1027 : !if (etot_yaml_doc%stream%length > 0) call etot_yaml_doc%write_and_free(ab_out)
1028 : end if ! nstep == 0 : no output
1029 :
1030 : case default
1031 0 : write(message, '(a,i0,a)' )' choice = ',choice,' is not an allowed value.'
1032 149264 : ABI_BUG(message)
1033 : end select
1034 :
1035 : ! Additional stuff for the Fock+SCF cycle
1036 83892 : if (present(fock)) then
1037 53730 : if (associated(fock)) then
1038 1498 : fock%fock_common%scf_converged=(quit==1)
1039 : ! At present, the decision that the Fock loop is converged is not taken here
1040 1498 : if (.not.fock%fock_common%fock_converged)quit=0
1041 : end if
1042 : end if
1043 :
1044 : ! Additional stuff for the two-component DFT SCF cycle (electrons+positron)
1045 83892 : if (present(electronpositron)) then
1046 53730 : if (associated(electronpositron)) then
1047 579 : electronpositron%scf_converged=(quit==1)
1048 579 : if (dtset%positron<0) then
1049 530 : diff_e=abs(etotal-electronpositron%etotal_prev)
1050 530 : diff_f=abs(maxfor-electronpositron%maxfor_prev)
1051 : end if
1052 579 : if (choice==1) then
1053 0 : ttoldff=0;ttoldfe=0
1054 0 : if(abs(dtset%postoldff)>tiny(0.0_dp))ttoldff=1
1055 0 : if(abs(dtset%postoldfe)>tiny(0.0_dp))ttoldfe=1
1056 0 : if (dtset%positron<0.and.ttoldff+ttoldfe/=1.and.iscf>0) then
1057 0 : ABI_ERROR('one and only one of toldff or toldfe must differ from zero !')
1058 : end if
1059 : end if
1060 579 : if (choice==2) then
1061 562 : if (dtset%positron<0.and.istep<=nstep) then
1062 519 : if (electronpositron%scf_converged) then
1063 81 : if (electronpositron%istep/=electronpositron%nstep) then
1064 : if ((.not.noquit).and.&
1065 75 : & (diff_e<electronpositron%postoldfe.or.diff_f<electronpositron%postoldff).and.&
1066 : & (mod(electronpositron%calctype,2)==0.or.(dtset%positron>-20.and.dtset%positron/=-2))) then
1067 5 : if (diff_e<electronpositron%postoldfe) then
1068 5 : write(message, '(2a,i5,5a,es11.3,a,es11.3)' ) ch10, &
1069 5 : & ' At SCF step',istep,', the difference between',ch10,&
1070 5 : & ' etotal from electronic calculation and etotal from positronic calculation',ch10,&
1071 10 : & ' is converged : diff(etot_el-etot_pos)=',diff_e,' < postoldfe=',electronpositron%postoldfe
1072 : else
1073 0 : write(message, '(2a,i5,5a,es11.3,a,es11.3)' ) ch10, &
1074 0 : & ' At SCF step',istep,', the difference between',ch10,&
1075 0 : & ' max. force from electronic calculation and max. force from positronic calculation',ch10,&
1076 0 : & ' is converged : diff(maxfor_el-maxfor_pos)=',diff_f,' < postoldff=',electronpositron%postoldff
1077 : end if
1078 15 : call wrtout([std_out, ab_out], message)
1079 : else
1080 70 : quit=0
1081 : end if
1082 : end if
1083 : end if
1084 : end if
1085 : end if
1086 579 : if (choice==3) then
1087 17 : if (dtset%positron<0.and.nstep>0)then
1088 11 : if (diff_e>=electronpositron%postoldfe.and.abs(dtset%postoldfe)>tiny(0.0_dp)) then
1089 5 : write(message, '(4a,i5,5a,es11.3,a,es11.3)' ) ch10,&
1090 5 : & ' scprqt: WARNING -',ch10,&
1091 5 : & ' posnstep=',dtset%posnstep,' was not enough SCF cycles to converge difference between',ch10,&
1092 5 : & ' etotal from electronic calculation and etotal from positronic calculation;',ch10,&
1093 10 : & ' diff=',diff_e,' exceeds postoldfe=',electronpositron%postoldfe
1094 15 : call wrtout([std_out, ab_out], message)
1095 : end if
1096 11 : if (diff_f>=electronpositron%postoldff.and.abs(dtset%postoldff)>tiny(0.0_dp)) then
1097 1 : write(message, '(4a,i5,5a,es11.3,a,es11.3)' ) ch10,&
1098 1 : & ' scprqt: WARNING -',ch10,&
1099 1 : & ' posnstep=',dtset%posnstep,' was not enough SCF cycles to converge difference between',ch10,&
1100 1 : & ' max. force from electronic calculation and max. force from positronic calculation;',ch10,&
1101 2 : & ' diff=',diff_e,' exceeds postoldff=',electronpositron%postoldff
1102 3 : call wrtout([std_out, ab_out], message)
1103 : end if
1104 : end if
1105 : end if
1106 : end if
1107 : end if
1108 :
1109 167784 : call flush_unit(ab_out)
1110 :
1111 : DBG_EXIT("COLL")
1112 :
1113 : contains
1114 :
1115 9225 : logical function converged()
1116 :
1117 : ! LB-02/01/2017:
1118 : ! This code avoids evaluation of undefined variables (which could happen in respfn, apparently)
1119 : logical :: loc_conv
1120 9225 : loc_conv = .true.
1121 9225 : if (ttolwfr==1) then
1122 2839 : if (residm > tolwfr) loc_conv=.false.
1123 : end if
1124 9225 : if (ttoldff==1) then
1125 2673 : if (diffor > toldff) loc_conv=.false.
1126 : end if
1127 9225 : if (ttoldmag==1) then
1128 2 : if (difmagsph > toldmag) loc_conv=.false.
1129 : end if
1130 9225 : if (ttolrff==1) then
1131 132 : if (diffor > tolrff*maxfor .and. maxfor > tol16) loc_conv=.false.
1132 : end if
1133 9225 : if (ttoldfe==1) then
1134 1266 : if (abs(deltae) > toldfe) loc_conv=.false.
1135 : end if
1136 9225 : if (ttolvrs==1) then
1137 2337 : if (res2 > tolvrs) loc_conv=.false.
1138 : end if
1139 9225 : converged = loc_conv
1140 :
1141 9225 : end function converged
1142 :
1143 : end subroutine scprqt
1144 : !!***
1145 :
1146 : !!****f* ABINIT/setup1
1147 : !! NAME
1148 : !! setup1
1149 : !!
1150 : !! FUNCTION
1151 : !! Call near top of main routine to handle setup of various arrays,
1152 : !! filenames, checking of input data, etc.
1153 : !!
1154 : !! INPUTS
1155 : !! acell(3)=length scales (bohr)
1156 : !! ecut_eff=effective energy cutoff (hartree) for planewave basis sphere
1157 : !! ecutc_eff=- PAW only - effective energy cutoff (hartree) for the coarse grid
1158 : !! natom=number of atoms
1159 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1160 : !! ngfftc(18)=contain all needed information about 3D FFT for the coarse grid
1161 : !! nkpt=number of k points
1162 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1163 : !! ntypat=number of types of atoms
1164 : !! response=0 if called by gstate, =1 if called by respfn
1165 : !! rprim(3,3)=dimensionless real space primitive translations
1166 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
1167 : !!
1168 : !! OUTPUT
1169 : !! bantot=total number of bands at all k points
1170 : !! gmet(3,3)=metric for reciprocal space inner products (bohr^-2)
1171 : !! gprimd(3,3)=dimens. primitive translations for reciprocal space (bohr**-1)
1172 : !! gsqcut_eff=Fourier cutoff on G^2 for "large sphere" of radius double
1173 : !! gsqcutc_eff=(PAW) Fourier cutoff on G^2 for "large sphere" of radius double for the coarse FFT grid
1174 : !! that of the basis sphere--appropriate for charge density rho(G),
1175 : !! Hartree potential, and pseudopotentials, corresponding to ecut_eff
1176 : !! rmet(3,3)=real space metric (bohr**2)
1177 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
1178 : !! ucvol=unit cell volume (bohr^3)
1179 : !!
1180 : !! NOTES
1181 : !! SHOULD BE CLEANED !
1182 : !!
1183 : !! SOURCE
1184 :
1185 12872 : subroutine setup1(acell,bantot,dtset,ecut_eff,ecutc_eff,gmet,&
1186 : & gprimd,gsqcut_eff,gsqcutc_eff,ngfft,ngfftc,nkpt,nsppol,&
1187 : & response,rmet,rprim,rprimd,ucvol,usepaw)
1188 :
1189 : !Arguments ------------------------------------
1190 : !scalars
1191 : type(dataset_type),intent(in) :: dtset
1192 : integer,intent(in) :: nkpt,nsppol
1193 : integer,intent(in) :: response,usepaw
1194 : integer,intent(out) :: bantot
1195 : real(dp),intent(in) :: ecut_eff,ecutc_eff
1196 : real(dp),intent(out) :: gsqcut_eff,gsqcutc_eff,ucvol
1197 : !arrays
1198 : integer,intent(in) :: ngfft(18),ngfftc(18)
1199 : real(dp),intent(in) :: acell(3),rprim(3,3)
1200 : real(dp),intent(out) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1201 : real(dp),intent(out) :: rprimd(3,3)
1202 :
1203 : !Local variables-------------------------------
1204 : !scalars
1205 : integer :: ikpt,isppol
1206 : real(dp) :: boxcut,boxcutc
1207 : character(len=500) :: msg
1208 : !arrays
1209 : real(dp) :: k0(3)
1210 :
1211 : ! ************************************************************************
1212 :
1213 : ! Compute bantot
1214 6436 : bantot=0
1215 13481 : do isppol=1,nsppol
1216 144578 : do ikpt=1,nkpt
1217 138142 : bantot=bantot+dtset%nband(ikpt+(isppol-1)*nkpt)
1218 : end do
1219 : end do
1220 :
1221 6436 : if(dtset%nqpt>1.or.dtset%nqpt<0) then
1222 : write(msg,'(a,i0,5a)')&
1223 0 : 'nqpt =',dtset%nqpt,' is not allowed',ch10,'(only 0 or 1 are allowed).',ch10,'Action: correct your input file.'
1224 0 : ABI_ERROR(msg)
1225 : end if
1226 :
1227 : ! Compute dimensional primitive translations rprimd
1228 6436 : call mkrdim(acell,rprim,rprimd)
1229 :
1230 : ! Obtain dimensional translations in reciprocal space gprimd,
1231 : ! metrics and unit cell volume, from rprimd.
1232 : ! Also output rprimd, gprimd and ucvol
1233 6436 : call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol)
1234 :
1235 : ! Get boxcut for given acell, gmet, ngfft, and ecut_eff
1236 : ! (center at 000 for groundstate, center at q for respfn):
1237 : ! boxcut=ratio of basis sphere diameter to fft box side
1238 6436 : k0(:)=0.0_dp
1239 6436 : if(response==1 .and. dtset%nqpt==1)then
1240 2408 : k0(:)=dtset%qptn(:)
1241 1806 : call wrtout([std_out, ab_out], ' setup1 : take into account q-point for computing boxcut.')
1242 : end if
1243 6436 : if (usepaw==1) then
1244 1427 : write(msg,'(2a)') ch10,' Coarse grid specifications (used for wave-functions):'
1245 4281 : call wrtout([std_out, ab_out], msg)
1246 1427 : call getcut(boxcutc,ecutc_eff,gmet,gsqcutc_eff,dtset%iboxcut,ab_out,k0,ngfftc)
1247 1427 : write(msg,'(2a)') ch10,' Fine grid specifications (used for densities):'
1248 4281 : call wrtout([std_out, ab_out], msg)
1249 1427 : call getcut(boxcut,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,ab_out,k0,ngfft)
1250 : else
1251 5009 : call getcut(boxcut,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,ab_out,k0,ngfft)
1252 5009 : gsqcutc_eff=gsqcut_eff
1253 : end if
1254 :
1255 : ! Check that boxcut>=2 if dtset%intxc=1; otherwise dtset%intxc must be set=0
1256 6436 : if (boxcut<2.0_dp.and.dtset%intxc==1) then
1257 : write(msg, '(a,es12.4,a,a,a,a,a)' )&
1258 0 : 'boxcut= ',boxcut,' is < 2.0 => intxc must be 0;',ch10,&
1259 0 : 'Need larger ngfft to use intxc=1.',ch10,&
1260 0 : 'Action: you could increase ngfft, or decrease ecut, or put intxcn=0.'
1261 0 : ABI_ERROR(msg)
1262 : end if
1263 :
1264 6436 : end subroutine setup1
1265 : !!***
1266 :
1267 : !!****f* ABINIT/prteigrs
1268 : !! NAME
1269 : !! prteigrs
1270 : !!
1271 : !! FUNCTION
1272 : !! Print out eigenvalues band by band and k point by k point.
1273 : !! If option=1, do it in a standard way, for self-consistent calculations.
1274 : !! If option=2, print out residuals and eigenvalues, in a format
1275 : !! adapted for nonself-consistent calculations, within the loops.
1276 : !! If option=3, print out eigenvalues, in a format
1277 : !! adapted for nonself-consistent calculations, at the end of the job.
1278 : !! If option=4, print out derivatives of eigenvalues (same format as option==3, except header that is printed)
1279 : !! If option=5, print out Fan contribution to zero-point motion correction to eigenvalues (averaged)
1280 : !! (same format as option==3, except header that is printed)
1281 : !! If option=6, print out DDW contribution to zero-point motion correction to eigenvalues (averaged)
1282 : !! (same format as option==3, except header that is printed)
1283 : !! If option=7, print out Fan+DDW contribution to zero-point motion correction to eigenvalues (averaged)
1284 : !! (same format as option==3, except header that is printed)
1285 : !!
1286 : !! INPUTS
1287 : !! eigen(mband*nkpt*nsppol)=eigenvalues (hartree)
1288 : !! or, if option==4, diagonal of derivative of eigenvalues
1289 : !! or, if option==5...7, zero-point motion correction to eigenvalues (averaged)
1290 : !! enunit=choice parameter: 0=>output in hartree; 1=>output in eV;
1291 : !! 2=> output in both hartree and eV
1292 : !! fermie=fermi energy (Hartree) / for electrons thermalized in the conduction bands when occopt==9
1293 : !! fermih=fermi energy (Hartree) for holes thermalized in the VB when occopt==9
1294 : !! fname_eig=filename of printing the eigenenergies
1295 : !! iout=unit number for formatted output file
1296 : !! iscf=option for self-consistency
1297 : !! kptns(3,nkpt)=k points in reduced coordinates
1298 : !! kptopt=option for the generation of k points
1299 : !! mband=maximum number of bands
1300 : !! nband(nkpt)=number of bands at each k point
1301 : !! nbdbuf= number of buffer bands
1302 : !! nkpt=number of k points
1303 : !! nnsclo_now=number of non-self-consistent loops for the current vtrial
1304 : !! (often 1 for SCF calculation, =nstep for non-SCF calculations)
1305 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1306 : !! occ(maxval(nband(:))*nkpt*nsppol)=occupancies for each band and k point
1307 : !! occopt=option for occupancies
1308 : !! option= (see above)
1309 : !! prteig=control print eigenenergies
1310 : !! prtvol=control print volume and debugging
1311 : !! resid(mband*nkpt*nsppol)=residuals (hartree**2)
1312 : !! tolwfr=tolerance on band residual of wf, hartrees**2 (needed when option=2)
1313 : !! vxcavg=average of vxc potential
1314 : !! wtk(nkpt)=k-point weights
1315 : !!
1316 : !! OUTPUT
1317 : !! (only writing)
1318 : !!
1319 : !! SOURCE
1320 :
1321 12817 : subroutine prteigrs(eigen,enunit,fermie,fermih,fname_eig,iout,iscf,kptns,kptopt,mband,nband,&
1322 12817 : & nbdbuf,nkpt,nnsclo_now,nsppol,occ,occopt,option,prteig,prtvol,resid,tolwfr,vxcavg,wtk)
1323 :
1324 : !Arguments ------------------------------------
1325 : !scalars
1326 : integer,intent(in) :: enunit,iout,iscf,kptopt,mband,nbdbuf,nkpt,nnsclo_now,nsppol
1327 : integer,intent(in) :: occopt,option,prteig,prtvol
1328 : real(dp),intent(in) :: fermie,fermih,tolwfr,vxcavg
1329 : character(len=*),intent(in) :: fname_eig
1330 : !arrays
1331 : integer,intent(in) :: nband(nkpt*nsppol)
1332 : real(dp),intent(in) :: eigen(mband*nkpt*nsppol),kptns(3,nkpt)
1333 : real(dp),intent(in) :: occ(mband*nkpt*nsppol),resid(mband*nkpt*nsppol)
1334 : real(dp),intent(in) :: wtk(nkpt)
1335 :
1336 : !Local variables-------------------------------
1337 : !scalars
1338 : integer,parameter :: nkpt_max=50
1339 : integer :: band_index,iband,ienunit,ii,ikpt,isppol,nband_index,nband_k,nkpt_eff,tmagnet,tmetal,temp_unit
1340 : real(dp) :: convrt,magnet,residk,rhodn,rhoup
1341 : character(len=2) :: ibnd_fmt,ikpt_fmt
1342 : character(len=7) :: strunit1,strunit2
1343 : character(len=39) :: kind_of_output
1344 : character(len=500) :: msg
1345 :
1346 : ! *************************************************************************
1347 :
1348 12817 : if (enunit<0.or.enunit>2) then
1349 0 : ABI_BUG(sjoin('enunit must be 0, 1 or 2. Argument was:', itoa(enunit)))
1350 : end if
1351 :
1352 12817 : if (prteig > 0) then
1353 8821 : call wrtout(iout, sjoin(' prteigrs : about to open file ', fname_eig))
1354 8821 : if (open_file(fname_eig, msg, newunit=temp_unit, status='unknown', form='formatted') /= 0) then
1355 0 : ABI_ERROR(msg)
1356 : end if
1357 8821 : rewind(temp_unit) ! always rewind disk file and print latest eigenvalues
1358 : end if
1359 :
1360 12817 : kind_of_output= ' Eigenvalues '
1361 12817 : if(option==4) kind_of_output=' Expectation of eigenvalue derivatives'
1362 12817 : if(option==5) kind_of_output=' Fan corrections to eigenvalues at T=0'
1363 12817 : if(option==6) kind_of_output=' DDW corrections to eigenvalues at T=0'
1364 12817 : if(option==7) kind_of_output=' Fan+DDW corrs to eigenvalues at T=0'
1365 :
1366 12817 : nkpt_eff=nkpt
1367 :
1368 : !write(msg,'(a,5i5)')' prtvol,iscf,kptopt,nkpt_eff,nkpt_max ',prtvol,iscf,kptopt,nkpt_eff,nkpt_max
1369 : !call wrtout(iout,msg)
1370 :
1371 12817 : if( (prtvol==0.or.prtvol==1) .and. (iscf/=-2 .or. kptopt>0) .and. nkpt_eff>nkpt_max)nkpt_eff=nkpt_max
1372 12817 : if( (prtvol==0.or.prtvol==1) .and. (iscf/=-2 .or. kptopt>0) .and. nkpt_eff>1 .and. iout==ab_out)nkpt_eff=1
1373 :
1374 12817 : if(option==1 .or. (option>=3 .and. option<=7))then
1375 :
1376 37149 : do ienunit=0,1
1377 :
1378 24766 : if (enunit==1 .and. ienunit==0)cycle
1379 24488 : if (enunit==0 .and. ienunit==1)cycle
1380 : ! Print eigenvalues in hartree for enunit=0 or 2
1381 : ! The definition of two different strings is quite ridiculous. Historical reasons ...
1382 :
1383 14935 : if (ienunit==0)then
1384 12105 : convrt=one
1385 12105 : strunit1='hartree'
1386 12105 : strunit2='hartree'
1387 : end if
1388 14935 : if (ienunit==1)then
1389 2830 : convrt=Ha_eV
1390 2830 : strunit1=' eV '
1391 2830 : strunit2='eV '
1392 : end if
1393 :
1394 14935 : band_index=0
1395 :
1396 14935 : if(ienunit==0)then ! XG20140730 I do not know why this is only done when ienunit==0
1397 12105 : tmetal=0
1398 12105 : if(option==1 .and. occopt>=3 .and. occopt<=8)tmetal=1
1399 2972 : tmagnet=0
1400 2972 : if(tmetal==1 .and. nsppol==2)then
1401 290 : tmagnet=1
1402 290 : rhoup = 0._dp
1403 290 : rhodn = 0._dp
1404 290 : nband_index = 1
1405 870 : do isppol=1,nsppol
1406 5310 : do ikpt=1,nkpt
1407 4440 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1408 76024 : do iband=1,nband_k
1409 71004 : if(isppol==1) rhoup = rhoup + wtk(ikpt)*occ(nband_index)
1410 71004 : if(isppol==2) rhodn = rhodn + wtk(ikpt)*occ(nband_index)
1411 75444 : nband_index = nband_index + 1
1412 : end do
1413 : end do
1414 : end do
1415 290 : magnet = abs(rhoup - rhodn)
1416 : end if
1417 : end if
1418 :
1419 14935 : if(iscf>=0 .and. (ienunit==0 .or. option==1))then
1420 13076 : if (occopt == 9) then
1421 : write(msg, '(3a,f10.5,a,f10.5,3a,f10.5)' ) &
1422 2 : ' Fermi energy for thermalized electrons and holes (',trim(strunit2),') =',&
1423 4 : convrt*fermie,', ',convrt*fermih,' Average Vxc (',trim(strunit2),')=',convrt*vxcavg
1424 : else
1425 : write(msg, '(3a,f10.5,3a,f10.5)' ) &
1426 13074 : ' Fermi (or HOMO) energy (',trim(strunit2),') =',convrt*fermie,' Average Vxc (',trim(strunit2),')=',convrt*vxcavg
1427 : end if
1428 13076 : call wrtout(iout,msg)
1429 13076 : if (prteig > 0) call wrtout(temp_unit,msg)
1430 : end if
1431 :
1432 : ! if( (iscf>=0 .or. iscf==-3) .and. ienunit==0)then ! This is the most correct
1433 14935 : if(iscf>=0 .and. ienunit==0)then ! For historical reasons
1434 10808 : if(tmagnet==1)then
1435 : write(msg, '(a,es16.8,a,a,es16.8,a,es16.8)' )&
1436 284 : ' Magnetization (Bohr magneton)=',magnet,ch10,&
1437 568 : ' Total spin up =',rhoup,' Total spin down =',rhodn
1438 284 : call wrtout(iout,msg)
1439 284 : if (prteig > 0) call wrtout(temp_unit,msg)
1440 : end if
1441 : end if
1442 :
1443 : ! Loop over spins (suppress spin data if nsppol not 2)
1444 44837 : do isppol=1,nsppol
1445 :
1446 17519 : ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9"
1447 17519 : if (nsppol==2.and.isppol==1) then
1448 2584 : write(msg, '(4a,'//ikpt_fmt//',2x,a)' )trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points, SPIN UP:'
1449 14935 : else if (nsppol==2.and.isppol==2) then
1450 2584 : write(msg, '(4a,'//ikpt_fmt//',2x,a)' )trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points, SPIN DOWN:'
1451 : else
1452 12351 : write(msg, '(4a,'//ikpt_fmt//',2x,a)' )trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points:'
1453 : end if
1454 17519 : call wrtout(iout,msg)
1455 17519 : if (prteig > 0) call wrtout(temp_unit,msg)
1456 :
1457 17519 : if(ienunit==0)then
1458 13889 : if(option>=4 .and. option<=7)then
1459 556 : msg = ' (in case of degenerate eigenvalues, averaged derivative)'
1460 556 : call wrtout(iout,msg)
1461 556 : if (prteig > 0) call wrtout(temp_unit,msg)
1462 : end if
1463 : end if
1464 :
1465 194964 : do ikpt=1,nkpt
1466 152679 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1467 152679 : ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9"
1468 152679 : ibnd_fmt="i3" ; if(nband_k>=1000)ibnd_fmt="i6" ; if(nband_k>=1000000)ibnd_fmt="i9"
1469 152679 : if(ikpt<=nkpt_eff)then
1470 : write(msg, '(a,'//ikpt_fmt//',a,'//ibnd_fmt//',a,f9.5,a,3f8.4,a)' ) &
1471 266116 : & ' kpt#',ikpt,', nband=',nband_k,', wtk=',wtk(ikpt)+tol10,', kpt=',kptns(1:3,ikpt)+tol10,' (reduced coord)'
1472 66529 : call wrtout(iout,msg)
1473 66529 : if (prteig > 0) call wrtout(temp_unit,msg)
1474 166370 : do ii=0,(nband_k-1)/8
1475 616773 : write(msg, '(8(f10.5,1x))' ) (convrt*eigen(iband+band_index), iband=1+ii*8,min(nband_k,8+ii*8))
1476 99841 : call wrtout(iout,msg)
1477 166370 : if (prteig > 0) call wrtout(temp_unit,msg)
1478 : end do
1479 66529 : if(ienunit==0 .and. option==1 .and. occopt>=3 .and. occopt<=8)then
1480 13848 : write(msg, '(5x,a,'//ikpt_fmt//')' ) ' occupation numbers for kpt#',ikpt
1481 13848 : call wrtout(iout,msg)
1482 41762 : do ii=0,(nband_k-1)/8
1483 203474 : write(msg, '(8(f10.5,1x))' ) (occ(iband+band_index),iband=1+ii*8,min(nband_k,8+ii*8))
1484 41762 : call wrtout(iout,msg)
1485 : end do
1486 : end if
1487 :
1488 : else
1489 86150 : if(ikpt==nkpt_eff+1)then
1490 3524 : write(msg, '(a,a)' )' prteigrs : prtvol=0 or 1, do not print more k-points.',ch10
1491 3524 : call wrtout(iout,msg)
1492 : end if
1493 86150 : if (prteig > 0) then
1494 : write(msg, '(a,'//ikpt_fmt//',a,'//ibnd_fmt//',a,f9.5,a,3f8.4,a)' ) &
1495 210972 : & ' kpt#',ikpt,', nband=',nband_k,', wtk=',wtk(ikpt)+tol10,', kpt=',kptns(1:3,ikpt)+tol10,' (reduced coord)'
1496 52743 : call wrtout(temp_unit,msg)
1497 126778 : do ii=0,(nband_k-1)/8
1498 481557 : write(msg, '(8(f10.5,1x))' ) (convrt*eigen(iband+band_index),iband=1+ii*8,min(nband_k,8+ii*8))
1499 126778 : call wrtout(temp_unit,msg)
1500 : end do
1501 : end if
1502 : end if
1503 170198 : band_index=band_index+nband_k
1504 : end do ! do ikpt=1,nkpt
1505 : end do ! do isppol=1,nsppol
1506 :
1507 : end do ! End loop over Hartree or eV
1508 :
1509 434 : else if(option==2)then
1510 :
1511 : band_index=0
1512 894 : do isppol=1,nsppol
1513 :
1514 460 : if(nsppol==2)then
1515 52 : if(isppol==1)write(msg, '(2a)' ) ch10,' SPIN UP channel '
1516 52 : if(isppol==2)write(msg, '(2a)' ) ch10,' SPIN DOWN channel '
1517 52 : call wrtout(iout,msg)
1518 52 : if(prteig>0) call wrtout(temp_unit,msg)
1519 : end if
1520 :
1521 15169 : do ikpt=1,nkpt
1522 14275 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1523 14275 : ikpt_fmt="i5" ; if(nkpt>=10000)ikpt_fmt="i7" ; if(nkpt>=1000000)ikpt_fmt="i9"
1524 :
1525 14275 : if(ikpt<=nkpt_eff)then
1526 : write(msg, '(1x,a,'//ikpt_fmt//',a,f9.5,2f9.5,a)' ) &
1527 3190 : & 'Non-SCF case, kpt',ikpt,' (',(kptns(ii,ikpt),ii=1,3),'), residuals and eigenvalues='
1528 3190 : call wrtout(iout,msg)
1529 3190 : if (prteig > 0) then
1530 : write(msg, '(1x,a,'//ikpt_fmt//',a,f9.5,2f9.5,a)' ) &
1531 3185 : & 'Non-SCF case, kpt',ikpt,' eig(',(kptns(ii,ikpt),ii=1,3),') '
1532 3185 : call wrtout(temp_unit,msg)
1533 : end if
1534 8606 : do ii=0,(nband_k-1)/8
1535 39588 : write(msg, '(1p,8e10.2)' )(resid(iband+band_index),iband=1+8*ii,min(8+8*ii,nband_k))
1536 8606 : call wrtout(iout,msg)
1537 : end do
1538 9715 : do ii=0,(nband_k-1)/6
1539 40697 : write(msg, '(1p,6e12.4)' )(eigen(iband+band_index),iband=1+6*ii,min(6+6*ii,nband_k))
1540 6525 : call wrtout(iout,msg)
1541 9715 : if (prteig > 0) call wrtout(temp_unit,msg)
1542 : end do
1543 : else
1544 11085 : if(ikpt==nkpt_eff+1)then
1545 226 : write(msg, '(a,a)' )' prteigrs : prtvol=0 or 1, do not print more k-points.',ch10
1546 226 : call wrtout(iout,msg)
1547 : end if
1548 11085 : if (prteig > 0) then
1549 9786 : write(msg, '(1x,a,i5,a,f9.5,2f9.5,a)' )'Non-SCF kpt',ikpt,' eig(',(kptns(ii,ikpt),ii=1,3),') '
1550 9786 : call wrtout(temp_unit,msg)
1551 28082 : do ii=0,(nband_k-1)/6
1552 103144 : write(msg, '(1p,6e12.4)' )(eigen(iband+band_index),iband=1+6*ii,min(6+6*ii,nband_k))
1553 28082 : call wrtout(temp_unit,msg)
1554 : end do
1555 : end if
1556 : end if
1557 :
1558 : ! Don't include the buffer in the output.
1559 14275 : if (nbdbuf>0) then
1560 99907 : residk=maxval(resid(band_index+1:band_index+nband_k-nbdbuf))
1561 : else
1562 31473 : residk=maxval(resid(band_index+1:band_index+nband_k))
1563 : end if
1564 14275 : if (residk>tolwfr) then
1565 : write(msg, '(1x,a,2i5,a,1p,e13.5)' ) &
1566 709 : ' prteigrs : nnsclo,ikpt=',nnsclo_now,ikpt,' max resid (excl. the buffer)=',residk
1567 709 : call wrtout(iout,msg)
1568 : end if
1569 :
1570 14735 : band_index=band_index+nband_k
1571 : end do
1572 : end do
1573 434 : call wrtout(iout," ")
1574 :
1575 : else
1576 0 : ABI_BUG(sjoin('option:', itoa(option),', is not allowed.'))
1577 : end if
1578 :
1579 12817 : if (prteig > 0) close (temp_unit)
1580 :
1581 12817 : end subroutine prteigrs
1582 : !!***
1583 :
1584 : !!****f* ABINIT/prtene
1585 : !!
1586 : !! NAME
1587 : !! prtene
1588 : !!
1589 : !! FUNCTION
1590 : !! Print components of total energy in nice format
1591 : !!
1592 : !! INPUTS
1593 : !! dtset <type(dataset_type)>=all input variables in this dataset
1594 : !! | berryphase
1595 : !! | kptopt
1596 : !! | occopt
1597 : !! | positron=option for electron-positron calculation
1598 : !! | tphysel="physical" electronic temperature with FD occupations
1599 : !! | tsmear=smearing energy or temperature (if metal)
1600 : !! energies <type(energies_type)>=values of parts of total energy
1601 : !! iout=unit number to which output is written
1602 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
1603 : !!
1604 : !! OUTPUT
1605 : !! (only writing)
1606 : !!
1607 : !! SOURCE
1608 :
1609 5277 : subroutine prtene(dtset,energies,iout,usepaw)
1610 :
1611 : !Arguments ------------------------------------
1612 : !scalars
1613 : integer,intent(in) :: iout,usepaw
1614 : type(dataset_type),intent(in) :: dtset
1615 : type(energies_type),intent(inout) :: energies
1616 :
1617 : !Local variables-------------------------------
1618 : !scalars
1619 : integer :: ipositron,optdc
1620 : logical :: directE_avail,testdmft,write_entropy=.false.,write_totalxc=.false.,write_epaw_core=.false.
1621 : real(dp) :: eent,enevalue,etotal,etotaldc,exc_semilocal,el_temp
1622 : ! Do not modify the length of these strings
1623 : character(len=14) :: eneName
1624 : character(len=500) :: info,msg
1625 : type(yamldoc_t) :: edoc,dc_edoc,pawcore_edoc,sdoc,ftxcdoc
1626 : !arrays
1627 : !character(len=10) :: EPName(1:2)=(/"Positronic","Electronic"/)
1628 :
1629 : ! *************************************************************************
1630 :
1631 1376 : directE_avail=(usepaw==0.or.dtset%pawspnorb==0.or.dtset%pawcpxocc==2.or.dtset%kptopt==1.or.dtset%kptopt==2)
1632 :
1633 : !============= Evaluate some parts of the energy ===========
1634 :
1635 5277 : optdc=-1;ipositron=merge(0,2,dtset%positron==0)
1636 5277 : if (abs(energies%e_ewald)<1.e-15_dp.and.abs(energies%e_hartree)<1.e-15_dp) ipositron=1
1637 5277 : call energies%eval_eint(dtset,usepaw,optdc,etotal,etotaldc)
1638 :
1639 5277 : call entropy(dtset,energies)
1640 5277 : eent=energies%e_entropy
1641 :
1642 : ! If DMFT is used and DMFT Entropy is not computed, then do not print
1643 : ! non interacting entropy
1644 : testdmft=(dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(dtset%upawu(:,1))>=tol8.or. &
1645 20547 : & sum(dtset%jpawu(:,1))>tol8).and.dtset%dmft_entropy==0)
1646 : if (dtset%usedmft==1.and.dtset%dmft_triqs_entropy==1.and.dtset%dmft_triqs_compute_integral>0 &
1647 5277 : & .and.(dtset%dmft_solv==6.or.dtset%dmft_solv==7)) testdmft=.false.
1648 5277 : if(testdmft) eent=zero
1649 :
1650 5277 : etotal = etotal + eent
1651 5277 : etotaldc = etotaldc + eent
1652 :
1653 : !============= Printing of Etotal by direct scheme ===========
1654 :
1655 5277 : if (dtset%icoulomb == 1) then
1656 0 : eneName = "Ion-ion energy"
1657 : else
1658 5277 : eneName = "Ewald energy"
1659 : end if
1660 5277 : enevalue = energies%e_ewald
1661 :
1662 :
1663 5277 : if (optdc==0.or.optdc==2) then
1664 :
1665 5277 : if (directE_avail) then
1666 5276 : info = 'Components of total free energy in Hartree'
1667 5276 : if(testdmft) info = 'Components of total energy in Hartree'
1668 5276 : edoc = yamldoc_open('EnergyTerms', info=trim(adjustl(info)), width=20, real_fmt='(es21.14)')
1669 5276 : call edoc%add_real('kinetic', energies%e_kinetic)
1670 5276 : if(abs(energies%e_extfpmd)>tiny(zero)) call edoc%add_real('extfpmd',energies%e_extfpmd)
1671 5276 : if (ipositron/=1) then
1672 5271 : exc_semilocal=energies%e_xc+energies%e_hybcomp_E0-energies%e_hybcomp_v0+energies%e_hybcomp_v
1673 : ! XG20181025 This should NOT be a part of the semilocal XC energy, but treated separately.
1674 : ! At present, there is still a problem with the variational formulation for the Fock term with PAW.
1675 : ! So, for the time being, keep it inside.
1676 5271 : if(usepaw==1)exc_semilocal=exc_semilocal+energies%e_fock
1677 5271 : call edoc%add_real('hartree', energies%e_hartree)
1678 5271 : call edoc%add_real('xc', exc_semilocal)
1679 5271 : call edoc%add_real(eneName, enevalue)
1680 5271 : call edoc%add_real('psp_core', energies%e_corepsp)
1681 : #if defined DEV_YP_VDWXC
1682 : if ( (dtset%vdw_xc > 0) .and. (dtset%vdw_xc < 10) .and. (xc_vdw_status()) ) then
1683 : call edoc%add_real('VdWaals_df', energies%e_xc_vdw)
1684 : end if
1685 : #endif
1686 : end if
1687 5276 : call edoc%add_real('local_psp', energies%e_localpsp)
1688 5276 : if (usepaw==0) then
1689 3901 : if(abs(energies%e_fock0)<tol8)then
1690 3812 : call edoc%add_real('non_local_psp', energies%e_nlpsp_vfock)
1691 : else
1692 89 : call edoc%add_real('non_local_psp+x', energies%e_nlpsp_vfock-energies%e_fock0)
1693 : endif
1694 : else
1695 1375 : if (dtset%use_rcpaw/=0) then
1696 3 : call edoc%add_real('PAW spherical_terms', energies%paw%epaw)
1697 3 : call edoc%add_real('PAW core', energies%paw%epaw_core)
1698 1372 : else if (dtset%paw_add_core==1) then
1699 1 : call edoc%add_real('PAW spherical_terms', energies%paw%epaw-energies%paw%epaw_core)
1700 1 : call edoc%add_real('PAW core', energies%paw%epaw_core)
1701 : else
1702 1371 : call edoc%add_real('PAW spherical_terms', energies%paw%epaw)
1703 : end if
1704 : !!!XG20181025 Does not work (yet)...
1705 : !!!if(abs(energies%e_nlpsp_vfock)>tol8)then
1706 : !!! write(msg, '(a,es21.14)' )' Fock-type term = ',energies%e_nlpsp_vfock
1707 : !!! call wrtout(iout,msg)
1708 : !!! write(msg, '(a,es21.14)' ) ' -frozen Fock en.= ',-energies%e_fock0
1709 : !!! call wrtout(iout,msg)
1710 : !!!endif
1711 : end if
1712 56256 : if (ANY(ABS(dtset%nucdipmom)>tol8)) then
1713 24 : call edoc%add_real('nucl. magn. dipoles',energies%e_nucdip)
1714 : end if
1715 5276 : if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).and.ipositron/=1) then
1716 4 : call edoc%add_real('VdWaals_dft_d', energies%e_vdw_dftd)
1717 : end if
1718 5276 : if (dtset%nzchempot>=1) then
1719 7 : call edoc%add_real('chem_potential', energies%e_chempot)
1720 : end if
1721 5276 : if (dtset%usedmft==1) then
1722 152 : call edoc%add_real('interaction', energies%e_hu)
1723 152 : call edoc%add_real('-double_counting', -energies%e_dc)
1724 : end if
1725 5276 : if(dtset%occopt>=3.and.dtset%occopt<=8.and.ipositron==0) then
1726 1276 : call edoc%add_real('internal', etotal-eent)
1727 1276 : if(.not.testdmft) then
1728 1159 : call edoc%add_real('-kT*entropy', eent)
1729 : end if
1730 3994 : else if (ipositron/=0) then
1731 14 : if (dtset%occopt>=3.and.dtset%occopt<=8) then
1732 6 : call edoc%add_real('-kT*entropy', eent)
1733 : end if
1734 : !write(msg, '(3a,es21.14,a)' ) &
1735 : ! ' >>> ',EPName(ipositron),' E= ',etotal-energies%e0_electronpositron -energies%e_electronpositron,ch10
1736 : !call wrtout(iout,msg)
1737 : !write(msg, '(3a,es21.14,2a,es21.14)' ) &
1738 : ! ' ',EPName(3-ipositron),' ener.= ',energies%e0_electronpositron,ch10,&
1739 : ! ' EP interaction E= ' ,energies%e_electronpositron
1740 : !call wrtout(iout,msg)
1741 14 : if(ipositron == 1) then
1742 5 : call edoc%add_real('positronic', etotal - energies%e0_electronpositron-energies%e_electronpositron)
1743 5 : call edoc%add_real('electronic', energies%e0_electronpositron)
1744 : else
1745 9 : call edoc%add_real('electronic', etotal- energies%e0_electronpositron-energies%e_electronpositron)
1746 9 : call edoc%add_real('positronic', energies%e0_electronpositron)
1747 : end if
1748 14 : call edoc%add_real('electron_positron_interaction', energies%e_electronpositron)
1749 : end if
1750 : if ((dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. &
1751 5276 : dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17) .and.ipositron/=1) then
1752 25 : call edoc%add_real('electric', energies%e_elecfield)
1753 25 : call edoc%add_real('kohn_sham', etotal-energies%e_elecfield)
1754 : end if
1755 5276 : call edoc%add_real('total_energy', etotal)
1756 :
1757 : else
1758 : write(msg, '(9a)' ) &
1759 1 : ' COMMENT: ',ch10,&
1760 1 : ' "Direct" decomposition of total free energy cannot be printed out !!!',ch10,&
1761 1 : ' PAW contribution due to spin-orbit coupling cannot be evaluated',ch10,&
1762 1 : ' without the knowledge of imaginary part of Rhoij atomic occupancies',ch10,&
1763 2 : ' (computed only when pawcpxocc=2).'
1764 1 : call wrtout(iout,msg)
1765 : end if
1766 : end if
1767 : !============= Printing of Etotal by double-counting scheme ===========
1768 :
1769 5277 : if (optdc>=1) then
1770 :
1771 1501 : info = '"Double-counting" decomposition of free energy'
1772 1501 : if(testdmft) info = '"Double-counting" decomposition of internal energy'
1773 : dc_edoc = yamldoc_open('EnergyTermsDC', info=trim(adjustl(info)), &
1774 1501 : width=20, real_fmt="(es21.14)")
1775 1501 : call dc_edoc%add_real('band_energy', energies%e_eigenvalues)
1776 1501 : if(abs(energies%e_extfpmd)>tiny(0.0_dp)) then
1777 7 : call dc_edoc%add_real('extfpmd_dc',energies%edc_extfpmd)
1778 : end if
1779 1501 : if (ipositron/=1) then
1780 : !write(msg, '(2(a,es21.14,a),a,es21.14)' ) &
1781 : ! ' '//eneName//' =',enevalue,ch10,&
1782 : ! ' PspCore energy = ',energies%e_corepsp-energies%e_corepspdc,ch10,&
1783 : ! ' Dble-C XC-energy= ',-energies%e_hartree+energies%e_xc-energies%e_xcdc -energies%e_fock0 + &
1784 : ! energies%e_hybcomp_E0-energies%e_hybcomp_v0
1785 : !call wrtout(iout,msg)
1786 1496 : call dc_edoc%add_real(eneName, enevalue)
1787 1496 : call dc_edoc%add_real('psp_core', energies%e_corepsp-energies%e_corepspdc)
1788 : call dc_edoc%add_real('xc_dc', -energies%e_hartree+energies%e_xc-energies%e_xcdc - energies%e_fock0 + &
1789 1496 : energies%e_hybcomp_E0-energies%e_hybcomp_v0)
1790 : end if
1791 : if ((dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. &
1792 1501 : dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17).and.ipositron/=1) then
1793 9 : call dc_edoc%add_real('electric_field', energies%e_elecfield)
1794 : end if
1795 1501 : if (usepaw==1) then
1796 1374 : if (dtset%use_rcpaw/=0) then
1797 3 : call dc_edoc%add_real('PAW spherical_terms', energies%paw%epaw_dc)
1798 3 : call dc_edoc%add_real('PAW core dc', energies%paw%epaw_core_dc)
1799 1371 : else if (dtset%paw_add_core==1) then
1800 1 : call dc_edoc%add_real('PAW spherical_terms', energies%paw%epaw_dc-energies%paw%epaw_core)
1801 1 : call dc_edoc%add_real('PAW core', energies%paw%epaw_core)
1802 : else
1803 1370 : call dc_edoc%add_real('PAW spherical_terms', energies%paw%epaw_dc)
1804 : end if
1805 : end if
1806 1501 : if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).and.ipositron/=1) then
1807 0 : call dc_edoc%add_real('VdWaals_dft_d', energies%e_vdw_dftd)
1808 : end if
1809 1501 : if (dtset%nzchempot>=1) then
1810 0 : call dc_edoc%add_real('chem_potential', energies%e_chempot)
1811 : end if
1812 1501 : if (dtset%usedmft==1) then
1813 152 : call dc_edoc%add_real('interaction', energies%e_hu)
1814 152 : call dc_edoc%add_real('-double_counting', -energies%e_dc)
1815 : end if
1816 1501 : if(dtset%occopt>=3.and.dtset%occopt<=8.and.ipositron==0) then
1817 689 : if(.not.testdmft) then
1818 : !write(msg, '(a,es21.14,a,a,a,es21.14)' ) &
1819 : ! ' >>>>> Internal E= ',etotaldc-eent,ch10,ch10,&
1820 : ! ' -kT*entropy = ',eent
1821 : !call wrtout(iout,msg)
1822 572 : call dc_edoc%add_real('internal', etotaldc-eent)
1823 572 : call dc_edoc%add_real('-kT*entropy', eent)
1824 : else
1825 117 : call dc_edoc%add_real('internal', etotaldc-eent)
1826 : end if
1827 806 : else if (ipositron/=0) then
1828 12 : if (dtset%occopt>=3 .and. dtset%occopt<=8) then
1829 6 : call dc_edoc%add_real('-kT*entropy', eent)
1830 : end if
1831 : !write(msg, '(a,es21.14,4a,es21.14,a)' ) &
1832 : ! ' - EP dble-ct En.= ',-energies%edc_electronpositron,ch10,&
1833 : ! ' >>> ',EPName(ipositron),' E= ',etotaldc-energies%e0_electronpositron -energies%e_electronpositron,ch10
1834 : !call wrtout(iout,msg)
1835 : !write(msg, '(3a,es21.14,2a,es21.14)' ) &
1836 : ! ' ',EPName(3-ipositron),' ener.= ',energies%e0_electronpositron,ch10,&
1837 : ! ' EP interaction E= ' ,energies%e_electronpositron
1838 : !call wrtout(iout,msg)
1839 12 : call dc_edoc%add_real('electron_positron_dc', -energies%edc_electronpositron)
1840 12 : if(ipositron == 1) then
1841 5 : call dc_edoc%add_real('positronic', etotaldc-energies%e0_electronpositron-energies%e_electronpositron)
1842 5 : call dc_edoc%add_real('electronic', energies%e0_electronpositron)
1843 : else
1844 7 : call dc_edoc%add_real('electronic', etotaldc-energies%e0_electronpositron-energies%e_electronpositron)
1845 7 : call dc_edoc%add_real('positronic', energies%e0_electronpositron)
1846 : end if
1847 12 : call dc_edoc%add_real('electron_positron_interaction', energies%e_electronpositron)
1848 : end if
1849 :
1850 1501 : write(msg, '(a,es21.14)' ) ' >>>> Etotal (DC)= ',etotaldc
1851 : !call wrtout(iout,msg)
1852 1501 : call dc_edoc%add_real('total_energy_dc', etotaldc)
1853 : end if
1854 :
1855 : !======= Additional printing ==========
1856 :
1857 5277 : if (usepaw==0.and.optdc==0) then
1858 3774 : call edoc%add_real('total_energy_eV', etotal*Ha_eV)
1859 3774 : call edoc%add_real('band_energy', energies%e_eigenvalues)
1860 : end if
1861 :
1862 5277 : if ((optdc==0.or.optdc==2).and.(.not.directE_avail)) then
1863 : !write(msg, '(a,a,es18.10)' ) ch10,' Band energy (Ha)= ',energies%e_eigenvalues
1864 : !call wrtout(iout,msg)
1865 1 : call edoc%add_real('band_energy', energies%e_eigenvalues)
1866 : end if
1867 :
1868 5277 : if (usepaw==1) then
1869 1376 : if ((optdc==0.or.optdc==2).and.(directE_avail)) then
1870 1375 : call edoc%add_real('total_energy_eV', etotal*Ha_eV)
1871 : end if
1872 1376 : if (optdc>=1) then
1873 : !if (optdc==1) write(msg, '(a,a,es21.14)' ) ch10,' >Total DC energy in eV = ',etotaldc*Ha_eV
1874 : !if (optdc==2) write(msg, '(a,es21.14)' ) ' >Total DC energy in eV = ',etotaldc*Ha_eV
1875 : !call wrtout(iout,msg)
1876 1374 : call dc_edoc%add_real('total_energy_dc_eV', etotaldc*Ha_eV)
1877 : end if
1878 : end if
1879 :
1880 5277 : if( dtset%icoulomb/=1.and.abs(dtset%cellcharge(1))>tol8) then
1881 : write(msg, '(6a)' ) &
1882 23 : ch10,' Calculation was performed for a charged system with PBC',&
1883 23 : ch10,' You may consider including the monopole correction to the total energy',&
1884 46 : ch10,' The correction is to be divided by the dielectric constant'
1885 23 : call wrtout(iout,msg)
1886 23 : call edoc%add_real('monopole_correction', energies%e_monopole)
1887 23 : call edoc%add_real('monopole_correction_eV', energies%e_monopole*Ha_eV)
1888 : end if
1889 :
1890 : !Print total energy including PAW core contribution
1891 5277 : if (usepaw==1) then
1892 1376 : if (dtset%paw_add_core==0.and.dtset%use_rcpaw==0) then
1893 1372 : if (abs(energies%paw%epaw_core)>tiny(zero).or.abs(energies%paw%epaw_core_dc)>tiny(zero)) then
1894 494 : write_epaw_core=.true.
1895 494 : info = 'Components of total free energy, including PAW core contributions'
1896 : pawcore_edoc = yamldoc_open('EnergyTermsWithPAWCore', info=trim(adjustl(info)), &
1897 494 : width=20, real_fmt='(es21.14)')
1898 494 : call pawcore_edoc%add_real('Total energy', etotal+energies%paw%epaw_core)
1899 494 : if(optdc>=1) call pawcore_edoc%add_real('Total energy DC', etotaldc+energies%paw%epaw_core)
1900 494 : call pawcore_edoc%add_real('Total energy (eV)', (etotal+energies%paw%epaw_core)*Ha_eV)
1901 494 : if(optdc>=1) call pawcore_edoc%add_real('Total energy DC (eV)', (etotaldc+energies%paw%epaw_core)*Ha_eV)
1902 : end if
1903 : end if
1904 : end if
1905 :
1906 : !======== In case other sources of entropies than the non-interacting entropy =========
1907 : !============= of the Kohn-Sham states come into play, print the details ==============
1908 5277 : if(dtset%occopt>=3.and.dtset%occopt<=8) then
1909 1282 : if(abs(energies%entropy)>tiny(zero).and.abs(energies%entropy-energies%entropy_ks)>tiny(zero)) then
1910 11 : write_entropy=.true.
1911 : sdoc = yamldoc_open('EntropyTerms', info='Components of total entropy', &
1912 11 : & width=20, real_fmt="(es21.14)") ! in kB units
1913 11 : call sdoc%add_real('noninteracting',energies%entropy_ks) ! Noninteracting entropy = Entropy of the Kohn-Sham states
1914 11 : if(abs(energies%entropy_xc)>tiny(zero)) call sdoc%add_real('xc',energies%entropy_xc)
1915 11 : if(usepaw==1.and.abs(energies%paw%entropy_paw)>tiny(zero)) call sdoc%add_real('spherical_terms',energies%paw%entropy_paw)
1916 11 : if(abs(energies%entropy_extfpmd)>tiny(zero)) call sdoc%add_real('extfpmd',energies%entropy_extfpmd)
1917 11 : if(abs(energies%entropy_imp)>tiny(zero)) call sdoc%add_real('impurity',energies%entropy_imp)
1918 11 : call sdoc%add_real('total_entropy',energies%entropy) ! Total entropy energy
1919 : end if
1920 : end if
1921 :
1922 : !======== In case finite-temperature exchange-correlation functionals are used =========
1923 : !=================== write the total exchange-correlation components ===================
1924 : !=============================== For testing purposes only =============================
1925 5277 : if(abs(energies%entropy_xc)>tiny(zero)) then
1926 5 : write_totalxc=.false. ! For testing purposes only, set write_totalxc=.true.
1927 : if(write_totalxc) then
1928 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
1929 : ftxcdoc = yamldoc_open('FTXCEnergyTerms', info='Components of total xc energy in Hartree', &
1930 : & width=20, real_fmt="(es21.14)")
1931 : if(usepaw==1) then
1932 : ! For now, only finite-temperature xc functionals contribute to entropy_paw.
1933 : ! We may introduce 'energies%entropy_pawxc' in the future.
1934 : call ftxcdoc%add_real('xc',energies%e_xc)
1935 : call ftxcdoc%add_real('spherical_terms_xc',energies%paw%epaw_xc)
1936 : call ftxcdoc%add_real('internal_xc',energies%e_xc+energies%paw%epaw_xc)
1937 : call ftxcdoc%add_real('-kT*entropy_xc',-el_temp*(energies%entropy_xc+energies%paw%entropy_paw))
1938 : call ftxcdoc%add_real('free_xc',energies%e_xc+energies%paw%epaw_xc-el_temp*(energies%entropy_xc+energies%paw%entropy_paw))
1939 : else
1940 : call ftxcdoc%add_real('internal_xc',energies%e_xc)
1941 : call ftxcdoc%add_real('-kT*entropy_xc',-el_temp*energies%entropy_xc)
1942 : call ftxcdoc%add_real('free_xc',energies%e_xc-el_temp*energies%entropy_xc)
1943 : end if
1944 : end if
1945 : end if
1946 :
1947 : ! Write components of total energies in Yaml format.
1948 5277 : call edoc%write_and_free(iout)
1949 5277 : if(optdc >= 1) call dc_edoc%write_and_free(iout)
1950 5277 : if (write_epaw_core) call pawcore_edoc%write_and_free(iout)
1951 5277 : if(write_entropy) call sdoc%write_and_free(iout)
1952 5277 : if(write_totalxc) call ftxcdoc%write_and_free(iout)
1953 :
1954 5277 : end subroutine prtene
1955 : !!***
1956 :
1957 : !!****f* ABINIT/get_dtsets_pspheads
1958 : !! NAME
1959 : !! get_dtsets_pspheads
1960 : !!
1961 : !! FUNCTION
1962 : !! Parse input file, get list of pseudos for files file and build list of datasets
1963 : !! pseudopotential headers, maxval of dimensions needed in outvars
1964 : !!
1965 : !! INPUTS
1966 : !! input_path: Input filename specifed on the command line. zero lenght if files file syntax is used.
1967 : !! Mainly used to check whether pseudos are defined in the input to avoid entering the files file
1968 : !! branch that prompts for pseudos.
1969 : !! path: Input Filename
1970 : !! comm: MPI communicator
1971 : !!
1972 : !! OUTPUT
1973 : !! lenstr= the length of the resulting string.
1974 : !! ndtset= the number of declared datasets.
1975 : !! string= contains on output the content of the file, ready for parsing.
1976 : !! dtsets(0:ndtset): List of datasets
1977 : !! dmatpuflag=flag controlling the use of an initial density matrix in PAW+U (max. value over datasets)
1978 : !! mx<ab_dimensions>=datatype storing the maximal dimensions.
1979 : !! pspheads(npsp)=<type pspheader_type>=all the important information from the
1980 : !! pseudopotential file headers, as well as the psp file names
1981 : !!
1982 : !! SOURCE
1983 :
1984 1440 : subroutine get_dtsets_pspheads(input_path, path, ndtset, lenstr, string, timopt, dtsets, pspheads, mx, dmatpuflag, comm)
1985 :
1986 : !Arguments ------------------------------------
1987 : !scalars
1988 : integer,intent(out) :: lenstr, ndtset
1989 : type(ab_dimensions),intent(out) :: mx
1990 : character(len=strlen), intent(out) :: string
1991 : character(len=*),intent(in) :: input_path, path
1992 : integer,intent(in) :: comm
1993 : integer,intent(out) :: timopt, dmatpuflag
1994 : !arrays
1995 : type(dataset_type),allocatable,intent(out) :: dtsets(:)
1996 : type(pspheader_type),allocatable,intent(out):: pspheads(:)
1997 :
1998 : !Local variables-------------------------------
1999 : !scalars
2000 : integer :: ipsp,ios, me, ndtset_alloc, nprocs
2001 : integer :: istatr,istatshft, papiopt, npsp, ii, idtset, msym, usepaw
2002 : character(len=fnlen) :: filpsp
2003 : character(len=500) :: msg
2004 : !arrays
2005 1440 : integer,allocatable :: mband_upper_(:)
2006 : real(dp) :: ecut_tmp(3,2,10),tsec(2)
2007 1440 : real(dp),allocatable :: zionpsp(:)
2008 1440 : character(len=fnlen), allocatable :: pspfilnam_(:), pseudo_paths(:)
2009 : !************************************************************************
2010 :
2011 1440 : me = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
2012 :
2013 : ! Read the file, stringify it and return the number of datasets.
2014 1440 : call parsefile(path, lenstr, ndtset, string, comm)
2015 :
2016 1440 : ndtset_alloc = ndtset; if (ndtset == 0) ndtset_alloc=1
2017 195608 : ABI_MALLOC(dtsets, (0:ndtset_alloc))
2018 :
2019 1440 : timopt = 1; if (xmpi_paral==1) timopt = 0
2020 :
2021 : ! Continue to analyze the input string, get upper dimensions, and allocate the remaining arrays.
2022 : call invars0(dtsets, istatr, istatshft, lenstr, msym, mx%natom, mx%nimage, mx%ntypat, &
2023 1440 : ndtset, ndtset_alloc, npsp, pseudo_paths, papiopt, timopt, string, comm)
2024 :
2025 : ! Enable PAPI timers
2026 1440 : call time_set_papiopt(papiopt)
2027 :
2028 9848 : dtsets(:)%timopt = timopt
2029 : dtsets(0)%timopt = 1
2030 1440 : if (xmpi_paral == 1) dtsets(0)%timopt = 0
2031 :
2032 1440 : call timab(timopt,5,tsec)
2033 :
2034 : ! Initialize pspheads, that contains the important information
2035 : ! from the pseudopotential headers, as well as the psp filename
2036 1440 : call timab(102,1,tsec)
2037 1440 : call timab(1021,3,tsec)
2038 :
2039 1440 : usepaw = 0
2040 6270 : ABI_MALLOC(pspheads, (npsp))
2041 1440 : if (npsp > 10) then
2042 0 : ABI_BUG('ecut_tmp is not well defined.')
2043 : end if
2044 131040 : ecut_tmp = -one
2045 :
2046 3390 : pspheads(:)%usewvl = dtsets(1)%usewvl
2047 :
2048 1440 : if (me == 0) then
2049 3486 : ABI_MALLOC(pspfilnam_, (npsp))
2050 :
2051 1162 : if (len_trim(pseudo_paths(1)) == 0) then
2052 : ! Enter Legacy `files file` mode --> Read the name of the psp file from files file.
2053 :
2054 : ! Catch possible mistake done by user (input without pseudos and `abinit t01.in` syntax)
2055 : ! else the code starts to prompt for pseudos and execution gets stuck
2056 2 : if (len_trim(input_path) /= 0) then
2057 0 : ABI_ERROR("`pseudos` variable must be specified in input when the code is invoked with the `abinit t01.in` syntax")
2058 : end if
2059 :
2060 : ! Finish to read the "file" file completely, as npsp is known,
2061 6 : do ipsp=1,npsp
2062 4 : write(std_out,'(/,a)' )' Please give name of formatted atomic psp file (and finish with a newline character)'
2063 4 : read (std_in, '(a)' , iostat=ios ) filpsp
2064 : ! It might be that a file name is missing
2065 4 : if (ios /= 0) then
2066 : write(msg, '(5a)' )&
2067 0 : 'There are not enough names of pseudopotentials provided in the files file.',ch10,&
2068 0 : 'Action: check first the variable ntypat (and/or npsp) in the input file;',ch10,&
2069 0 : 'if they are correct, complete your files file.'
2070 0 : ABI_ERROR(msg)
2071 : end if
2072 4 : pspfilnam_(ipsp) = trim(filpsp)
2073 6 : write(std_out,'(a,i0,2a)' )' For atom type ',ipsp,', psp file is ',trim(filpsp)
2074 : end do ! ipsp
2075 :
2076 : else
2077 : ! Get pseudopotential paths from input file.
2078 3886 : pspfilnam_ = pseudo_paths
2079 2726 : do ipsp=1,npsp
2080 2726 : write(std_out,'(a,i0,2a)' )' For atom type ',ipsp,', psp file is ',trim(pspfilnam_(ipsp))
2081 : end do
2082 : end if
2083 :
2084 : ! Now read the psp headers
2085 1162 : call inpspheads(pspfilnam_, npsp, pspheads, ecut_tmp)
2086 1162 : ABI_FREE(pspfilnam_)
2087 :
2088 2732 : if (minval(abs(pspheads(1:npsp)%pspcod - 7)) == 0) usepaw=1
2089 2732 : if (minval(abs(pspheads(1:npsp)%pspcod - 17)) == 0) usepaw=1
2090 : end if ! me == 0
2091 :
2092 1440 : ABI_FREE(pseudo_paths)
2093 :
2094 : ! Communicate pspheads to all processors
2095 1440 : call pspheads_comm(npsp, pspheads, usepaw)
2096 :
2097 : ! If (all) pspcod are 7 then this is a PAW calculation. Initialize (default) the value of ratsph
2098 9848 : do idtset=0,ndtset_alloc
2099 8408 : dtsets(idtset)%usepaw = usepaw
2100 9848 : if (usepaw == 0) then
2101 15042 : dtsets(idtset)%ratsph(:)=two
2102 : else
2103 : ! Note that the following coding assumes that npsp=ntypat for PAW, which is true as of now (XG20101024).
2104 : ! dtsets(idtset)%ratsph(1:npsp)=token%pspheads(1:npsp)%pawheader%rpaw
2105 3903 : do ipsp=1,npsp
2106 3903 : dtsets(idtset)%ratsph(ipsp) = pspheads(ipsp)%pawheader%rpaw
2107 : end do
2108 : endif
2109 : end do
2110 :
2111 : ! Take care of other dimensions, and part of the content of dtsets that is or might be needed early.
2112 4320 : ABI_MALLOC(zionpsp, (npsp))
2113 3390 : do ii=1,npsp
2114 3390 : zionpsp(ii) = pspheads(ii)%zionpsp
2115 : end do
2116 :
2117 4320 : ABI_MALLOC(mband_upper_, (0:ndtset_alloc))
2118 :
2119 : ! Get MAX dimension over datasets
2120 : call invars1m(dmatpuflag, dtsets, ab_out, lenstr, mband_upper_, mx,&
2121 1440 : msym, ndtset, ndtset_alloc, string, npsp, zionpsp, comm)
2122 :
2123 1427 : ABI_FREE(zionpsp)
2124 1427 : call timab(1021,2,tsec)
2125 1427 : call timab(1022,3,tsec)
2126 :
2127 : ! Provide defaults for the variables that have not yet been initialized.
2128 1427 : call indefo(dtsets, ndtset_alloc, nprocs)
2129 :
2130 1427 : call timab(1022,2,tsec)
2131 1427 : call timab(1023,3,tsec)
2132 :
2133 : ! Perform some global initialization, depending on the value of
2134 : ! pseudopotentials, parallelism variables, or macro input variables
2135 1427 : call macroin(dtsets, ecut_tmp, lenstr, ndtset_alloc, string)
2136 :
2137 : ! If all the pseudopotentials have the same pspxc, override the default value for dtsets 1 to ndtset
2138 3364 : if (minval(abs((pspheads(1:npsp)%pspxc - pspheads(1)%pspxc)))==0) then
2139 8382 : dtsets(1:ndtset_alloc)%ixc = pspheads(1)%pspxc
2140 : end if
2141 :
2142 : ! Call the main input routine.
2143 1427 : call invars2m(dtsets,ab_out,lenstr,mband_upper_,msym,ndtset,ndtset_alloc,npsp,pspheads,string, comm)
2144 :
2145 1427 : call macroin2(dtsets, ndtset_alloc)
2146 :
2147 1427 : mx%mband = dtsets(1)%mband
2148 8382 : do ii=1,ndtset_alloc
2149 8382 : mx%mband = max(dtsets(ii)%mband, mx%mband)
2150 : end do
2151 :
2152 1427 : call timab(1023,2,tsec)
2153 1427 : call timab(102,2,tsec)
2154 :
2155 1427 : ABI_FREE(mband_upper_)
2156 :
2157 2880 : end subroutine get_dtsets_pspheads
2158 : !!***
2159 :
2160 : !!****f* ABINIT/ebands_from_file
2161 : !! NAME
2162 : !! ebands_from_file
2163 : !!
2164 : !! FUNCTION
2165 : !! Build and ebands_t object from file. Supports Fortran and netcdf files
2166 : !! provided they have a Abinit header and obviously GS eigenvalues
2167 : !!
2168 : !! INPUTS
2169 : !! path: File name.
2170 : !! comm: MPI communicator.
2171 : !!
2172 : !! OUTPUT
2173 : !!
2174 : !! SOURCE
2175 :
2176 :
2177 0 : type(ebands_t) function ebands_from_file(path, comm) result(new)
2178 :
2179 : !Arguments ------------------------------------
2180 : !scalars
2181 : character(len=*),intent(in) :: path
2182 : integer,intent(in) :: comm
2183 :
2184 : !Local variables-------------------------------
2185 : !scalars
2186 : integer :: ncid, fform
2187 0 : type(hdr_type) :: hdr
2188 : !arrays
2189 0 : real(dp),pointer :: gs_eigen(:,:,:)
2190 :
2191 : ! *************************************************************************
2192 :
2193 : ! NOTE: Assume file with header. Must use wfk_read_eigenvalues to handle Fortran WFK
2194 0 : if (endswith(path, "_WFK") .or. endswith(path, "_WFK.nc")) then
2195 0 : call wfk_read_eigenvalues(path, gs_eigen, hdr, comm)
2196 0 : call new%from_hdr(hdr, maxval(hdr%nband), gs_eigen)
2197 :
2198 0 : else if (endswith(path, ".nc")) then
2199 0 : NCF_CHECK(nctk_open_read(ncid, path, comm))
2200 0 : call hdr%ncread(ncid, fform)
2201 0 : ABI_CHECK(fform /= 0, "fform == 0")
2202 0 : ABI_MALLOC(gs_eigen, (hdr%mband, hdr%nkpt, hdr%nsppol))
2203 0 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "eigenvalues"), gs_eigen))
2204 0 : call new%from_hdr(hdr, maxval(hdr%nband), gs_eigen)
2205 0 : NCF_CHECK(nf90_close(ncid))
2206 : else
2207 0 : ABI_ERROR(sjoin("Don't know how to construct crystal structure from: ", path, ch10, "Supported extensions: _WFK or .nc"))
2208 : end if
2209 :
2210 0 : ABI_FREE(gs_eigen)
2211 0 : call hdr%free()
2212 :
2213 0 : end function ebands_from_file
2214 : !!***
2215 :
2216 : !!****f* ABINIT/crystal_from_file
2217 : !! NAME
2218 : !! crystal_from_file
2219 : !!
2220 : !! FUNCTION
2221 : !! Build crystal_t object from netcdf file
2222 : !!
2223 : !! INPUTS
2224 : !!
2225 : !! OUTPUT
2226 : !!
2227 : !! SOURCE
2228 :
2229 0 : type(crystal_t) function crystal_from_file(path, comm) result(new)
2230 :
2231 : !Arguments ------------------------------------
2232 : !scalars
2233 : character(len=*),intent(in) :: path
2234 : integer,intent(in) :: comm
2235 :
2236 : !Local variables-------------------------------
2237 : !scalars
2238 : integer :: fform, ncid
2239 0 : type(hdr_type) :: hdr
2240 : ! *************************************************************************
2241 :
2242 0 : if (endswith(path, ".nc")) then
2243 0 : NCF_CHECK(nctk_open_read(ncid, path, comm))
2244 0 : call new%ncread(ncid)
2245 0 : NCF_CHECK(nf90_close(ncid))
2246 : else
2247 : ! Assume file with Abinit header
2248 : ! TODO: Should add routine to read crystal from structure without hdr
2249 0 : call hdr%from_fname(path, fform, comm)
2250 0 : ABI_CHECK(fform /= 0, "fform == 0")
2251 0 : new = hdr%get_crystal()
2252 0 : call hdr%free()
2253 : end if
2254 :
2255 0 : end function crystal_from_file
2256 : !!***
2257 :
2258 :
2259 : !!****f* ABINIT/get_gemm_nonlop_ompgpu_blocksize
2260 : !! NAME
2261 : !! get_gemm_nonlop_ompgpu_blocksize
2262 : !!
2263 : !! FUNCTION
2264 : !! Check GPU memory constraints for many ABINIT routines (lobpcg,chebfi,getghc)
2265 : !! and compute a size of GEMM nonlop block that is small enough to fit in GPU memory.
2266 : !! Print its estimation of memory consumption in output.
2267 : !!
2268 : !! INPUTS
2269 : !! gs_hamk : hamiltonian structure (used by getghc)
2270 : !! ndat : size of batching, usually matching bandpp or nblock_lobpcg
2271 : !! npw : number of planewaves
2272 : !! nband : number of bands
2273 : !! nspinor : number of spin-orbits
2274 : !! paral_kgb : 1 if KGB parallelism is enabled, 0 otherwise
2275 : !! optfor : 1 if forces computation is enabled, 0 otherwise
2276 : !! optstr : 1 if stress computation is enabled, 0 otherwise
2277 : !! wfoptalg : Which diago algorithm if used:
2278 : !! -1: none (GEMM nonlop isn't used for diago)
2279 : !! 111: CHEBFI2
2280 : !! 112: SLICE
2281 : !! 114: LOBPCG2
2282 : !! other: Only account for getghc
2283 : !! gpu_option : If GPU is enabled (expected to be ABI_GPU_OPENMP for now)
2284 : !! blocksize : if higher than 0, only print memory estimation and exit
2285 : !!
2286 : !! OUTPUT
2287 : !! nfft_blocks : Number of blocks fourwf computation will be divided into
2288 : !! blocksize : Size of MPI tasks blocks to be used in GEMM nonlop
2289 : !! nblocks : Number of MPI blocks to be used in GEMM nonlop
2290 : !!
2291 : !! SOURCE
2292 0 : subroutine get_gemm_nonlop_ompgpu_blocksize(ikpt,gs_hamk,ndat,nband,nspinor,nspden,paral_kgb,&
2293 : & npband,optfor,optstr,wfoptalg,gpu_option,use_distrib,&
2294 : & blocksize,nblocks,nfft_blocks,warn_on_fail,disable_output)
2295 :
2296 : integer,intent(in) :: ikpt,ndat,nband,nspinor,nspden,paral_kgb,npband,optfor,optstr,wfoptalg,gpu_option
2297 : logical,intent(in) :: use_distrib
2298 : logical,intent(in),optional :: warn_on_fail,disable_output
2299 : type(gs_hamiltonian_type),intent(in) :: gs_hamk
2300 : integer,intent(inout) :: blocksize,nfft_blocks
2301 : integer,intent(out) :: nblocks
2302 :
2303 : integer(kind=c_size_t) :: nonlop_smem,invovl_smem,getghc_wmem,invovl_wmem,nonlop_wmem,gs_ham_smem,updrho_wmem,prep_nonlop_wmem
2304 : integer(kind=c_size_t) :: sum_mem,sum_bandpp_mem,sum_other_mem,free_mem,localMem,fourwf_smem,fourwf_wmem,fourwf_mem,hegvd_mem
2305 : integer :: icplx,space,i,ndat_try,rank,nprocs,ndgxdt,blockdim,max_slices,npw,npw_fft,signs,nprojs,itypat
2306 : integer, target :: t_fft(3)
2307 : logical :: print_and_exit,l_warn_on_fail,l_disable_output,fixed_blocksize,fixed_fourwf_blocks
2308 : integer(kind=c_size_t) :: chebfiMem(2),lobpcgMem(2)
2309 : character(len=500) :: message
2310 :
2311 : ! *********************************************************************
2312 :
2313 0 : free_mem=256*1e9 ! Dummy value
2314 0 : l_warn_on_fail=.false.;if(present(warn_on_fail)) l_warn_on_fail=warn_on_fail
2315 0 : l_disable_output=.false.;if(present(disable_output)) l_disable_output=disable_output
2316 : #ifdef HAVE_GPU
2317 : if(gpu_option /= ABI_GPU_DISABLED) then
2318 : call gpu_get_max_mem(free_mem)
2319 : ! NOTE: Cutting 10% out to be safe
2320 : ! I computed this value using the memory allocation summary that
2321 : ! NVHPC provides after a GPU out-of-memory crash, on a run NVIDIA A100 with 80GB.
2322 : ! 8.5% is the amount of memory that wasn't reported allocated nor freed.
2323 : ! I add extra 1.5% to account for minor buffers eventually allocated.
2324 : ! This amount may either be hidden allocations for CUDA/cuBLAS/cuFFT/cuSOLVER,
2325 : ! or memory allocated outside ABINIT.
2326 : free_mem = 0.9 * free_mem
2327 : end if
2328 : #else
2329 : ABI_UNUSED(gpu_option)
2330 : #endif
2331 :
2332 0 : if(gpu_option /= ABI_GPU_OPENMP) then
2333 : ! No distribution is attempted outside of OpenMP GPU. User is already warned in chkinp
2334 0 : blocksize=1; nblocks=0; nfft_blocks=1
2335 0 : return
2336 : end if
2337 :
2338 0 : rank = xmpi_comm_rank(xmpi_world); nprocs = xmpi_comm_size(xmpi_world)
2339 0 : if ( gs_hamk%istwf_k == 2 ) then ! Real only
2340 0 : space = SPACE_CR
2341 0 : icplx = 2
2342 : else ! complex
2343 0 : space = SPACE_C
2344 0 : icplx = 1
2345 : end if
2346 :
2347 0 : npw=gs_hamk%npw_k
2348 0 : npw_fft=gs_hamk%npw_fft_k
2349 0 : ndat_try=ndat
2350 0 : blockdim=npband*ndat
2351 0 : ndgxdt=0
2352 0 : if(optfor>0) ndgxdt=ndgxdt+3
2353 0 : if(optstr>0) ndgxdt=ndgxdt+6
2354 0 : signs=2
2355 : !wfoptalg==-1 means we're in forstr
2356 0 : if(wfoptalg==-1) signs=1
2357 0 : t_fft(1) = gs_hamk%ngfft(3);
2358 0 : t_fft(2) = gs_hamk%ngfft(2);
2359 0 : t_fft(3) = gs_hamk%ngfft(1);
2360 0 : nprojs=0
2361 0 : do itypat=1,gs_hamk%ntypat
2362 0 : nprojs = nprojs + count(gs_hamk%indlmn(3,:,itypat)>0)*gs_hamk%nattyp(itypat)
2363 : end do
2364 :
2365 0 : nonlop_smem=0; invovl_smem=0; getghc_wmem=0; invovl_wmem=0; nonlop_wmem=0; gs_ham_smem=0
2366 0 : updrho_wmem=0; prep_nonlop_wmem=0; sum_mem=0; sum_bandpp_mem=0; sum_other_mem=0;
2367 0 : localMem=0; fourwf_smem=0; fourwf_wmem=0; fourwf_mem=0; hegvd_mem=0
2368 0 : chebfiMem(:)=0; lobpcgMem(:)=0
2369 :
2370 : !HEGVD work memory estimate.
2371 : ! Since *_bufferSize routines from (cu/hip)SOLVER require buffer
2372 : ! to be provided, I measeured the work size given by those routines
2373 : ! on many big cases and guess an approximate value.
2374 : ! hipSolver is eager than cuSolver, hence the extra multiplier
2375 :
2376 0 : if(wfoptalg==111 .or. wfoptalg==11) then
2377 0 : hegvd_mem = int(dp, c_size_t) * nband * nband * 3
2378 0 : else if(wfoptalg==114 .or. wfoptalg==14) then
2379 0 : hegvd_mem = int(dp, c_size_t) * (ndat*3) * (ndat*3) * 3
2380 : end if
2381 : #ifdef HAVE_GPU_CUDA
2382 : if(space == SPACE_C) hegvd_mem = hegvd_mem * 2
2383 : #endif
2384 : #ifdef HAVE_GPU_HIP
2385 : if(space == SPACE_C) hegvd_mem = hegvd_mem * 3
2386 : ! ROCm 7 memory usage was measured to be more or less
2387 : ! on par with CUDA but ROCm 6 was indeed ~9 times higher.
2388 : ! For now, we can't drop ROCm 6 so we keep this workaround.
2389 : if(gpu_get_lib_version_major() < 7) then
2390 : hegvd_mem = hegvd_mem * 9
2391 : end if
2392 : #endif
2393 :
2394 : if(wfoptalg>=0) then
2395 : #ifdef HAVE_GPU
2396 : call gpu_fft_get_estimate_work_size(3, c_loc(t_fft), FFT_Z2Z, ndat, fourwf_smem);
2397 : #endif
2398 : end if
2399 :
2400 0 : nonlop_smem = gemm_nonlop_ompgpu_static_mem(npw_fft, gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%ntypat, max(1,blocksize), ndgxdt, use_distrib)
2401 0 : getghc_wmem = getghc_ompgpu_work_mem(gs_hamk, ndat, max(nfft_blocks,1))
2402 : fourwf_wmem = int(2, c_size_t) * dp * gs_hamk%n4 * gs_hamk%n5 * gs_hamk%n6 &
2403 0 : & * (ndat/max(nfft_blocks,1) + modulo(ndat,max(nfft_blocks,1)))
2404 0 : fourwf_mem = fourwf_wmem+fourwf_smem
2405 :
2406 : nonlop_wmem = gemm_nonlop_ompgpu_work_mem(gs_hamk%istwf_k, ndat, ndgxdt, npw_fft,&
2407 0 : & gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%ntypat, gs_hamk%lmnmax, signs, wfoptalg)
2408 0 : gs_ham_smem = int(2,c_size_t)*npw_fft*size(gs_hamk%ffnl_k,dim=3)*size(gs_hamk%ffnl_k,dim=4) + int(3,c_size_t)*npw_fft
2409 0 : if(associated(gs_hamk%ph3d_k)) gs_ham_smem = gs_ham_smem + int(2,c_size_t) * npw_fft * gs_hamk%matblk
2410 0 : gs_ham_smem = gs_ham_smem*dp
2411 :
2412 : ! Mapped arrays used in mkrho or vtowfk
2413 0 : updrho_wmem = int(2,c_size_t)*npw_fft*ndat*nspinor ! cwavef
2414 0 : if(nspden==4) updrho_wmem = updrho_wmem + int(2,c_size_t)*npw_fft*ndat*2 ! cwavef_x + cwavef_y
2415 0 : updrho_wmem = updrho_wmem*dp
2416 :
2417 : ! Mapped arrays used in prep_nonlop (only paral_kgb==1)
2418 0 : prep_nonlop_wmem = 0
2419 0 : if(paral_kgb==1) prep_nonlop_wmem = int(2,c_size_t)*npw_fft*ndat*nspinor*3*dp ! cwavef_alltoall{1,2} + cwavef
2420 :
2421 0 : if(wfoptalg==111 .or. wfoptalg==112) then
2422 0 : chebfiMem = chebfi_memInfo(nband,icplx*npw*nspinor,space,paral_kgb,icplx*npw*nspinor,blockdim)
2423 0 : invovl_smem = invovl_ompgpu_static_mem(gs_hamk)
2424 0 : invovl_wmem = invovl_ompgpu_work_mem(gs_hamk, ndat)
2425 : end if
2426 0 : if(wfoptalg==114) then
2427 0 : lobpcgMem = lobpcg_memInfo(nband,icplx*npw*nspinor,space,paral_kgb,blockdim)
2428 : end if
2429 0 : localMem = (int(2,c_size_t)*npw*nspinor*nband+3*nband)*kind(1.d0) ! cg, eig, occ, resid in chebfiwf/lobpcgwf
2430 :
2431 : ! Check if arrays outside of GEMM nonlop projectors and ompgpu_fourwf fit in GPU memory
2432 0 : sum_other_mem = gs_ham_smem
2433 :
2434 0 : if(wfoptalg>=0) then
2435 0 : sum_other_mem = sum_other_mem+updrho_wmem+prep_nonlop_wmem
2436 : else
2437 0 : sum_other_mem = sum_other_mem+nonlop_wmem+prep_nonlop_wmem
2438 : end if
2439 :
2440 0 : if(wfoptalg==111) then
2441 0 : sum_other_mem = sum_other_mem + invovl_wmem+invovl_smem+chebfiMem(1)+chebfiMem(2)+localMem
2442 : end if
2443 :
2444 0 : if(wfoptalg==114) then
2445 0 : sum_other_mem = sum_other_mem + lobpcgMem(1)+lobpcgMem(2)+localMem
2446 : end if
2447 :
2448 0 : print_and_exit=.false.
2449 0 : fixed_blocksize=.false.
2450 0 : fixed_fourwf_blocks=.false.
2451 0 : nblocks=0
2452 0 : if(blocksize > 0 .and. nfft_blocks > 0) then
2453 0 : nblocks=max(1,nprocs/blocksize)
2454 0 : print_and_exit=.true.
2455 : else
2456 0 : if(blocksize > 0) then
2457 0 : nblocks=max(1,nprocs/blocksize)
2458 0 : fixed_blocksize=.true.
2459 : else
2460 0 : blocksize=1
2461 0 : if(.not. l_disable_output) write(std_out,*) "Setting GEMM nonlop block number...", new_line('A')
2462 : end if
2463 :
2464 0 : if(nfft_blocks > 0) then
2465 : fixed_fourwf_blocks=.true.
2466 : else
2467 0 : nfft_blocks=1
2468 0 : if(.not. l_disable_output) write(std_out,*) "Setting FFT blocks number...", new_line('A')
2469 : end if
2470 : end if
2471 :
2472 0 : max_slices=max(100,nprocs*2); if(sum_other_mem > free_mem) max_slices=1
2473 : ! How we try to optimize GPU memory consumption:
2474 : ! We work on two variables :
2475 : ! - blocksize : for slicing GEMM nonlop projectors arrays
2476 : ! - nfft_blocks : for slicing ompgpu_fourwf work buffer
2477 : ! At each iteration, we check which routine from GEMM nonlop or ompgpu_fourwf
2478 : ! have the higher memory requirement and increase related variable.
2479 : !
2480 : ! We arbitraly give at least 100 tries to improve GPU memory consumption before eventually
2481 : ! ruling out that use case is too big to run with available GPU memory and abort.
2482 : !
2483 : ! User may hard set slicing for both fourwf and GEMM nonlop, in which case the code will
2484 : ! warn the user about possible GPU memory overpassing instead of aborting.
2485 : !
2486 : ! However, if arrays from other parts of the code already have higher memory requirements,
2487 : ! we fail anyway and advise the user to increase nblock_lobpcg or run on more nodes.
2488 0 : do i=1,max_slices
2489 :
2490 : ! First iteration or user provided parameters to split fourwf and GEMM nonlop
2491 : ! Just measure
2492 0 : if(i==1 .or. print_and_exit) then
2493 0 : if(wfoptalg>=0) then
2494 : #ifdef HAVE_GPU
2495 : call gpu_fft_get_estimate_work_size(3, c_loc(t_fft), FFT_Z2Z, ndat/nfft_blocks, fourwf_smem);
2496 : #endif
2497 0 : getghc_wmem = getghc_ompgpu_work_mem(gs_hamk, ndat, nfft_blocks)
2498 : fourwf_wmem = int(2, c_size_t) * dp * gs_hamk%n4 * gs_hamk%n5 * gs_hamk%n6 &
2499 0 : & * (ndat/nfft_blocks + modulo(ndat,nfft_blocks))
2500 0 : fourwf_mem = fourwf_wmem + fourwf_smem
2501 : end if
2502 : nonlop_smem = gemm_nonlop_ompgpu_static_mem(npw_fft,gs_hamk%indlmn,gs_hamk%nattyp,&
2503 0 : & gs_hamk%ntypat,blocksize,ndgxdt,use_distrib)
2504 : else
2505 : ! Raise fourwf slicing if :
2506 : ! - GEMM nonlop block has been set by user
2507 : ! or
2508 : ! - fourwf memory requirements are higher
2509 : ! - fourwf slicing wasn't set by user
2510 : ! - fourwf is still sliceable
2511 : !
2512 : ! Raise GEMM nonlop blocks otherwise
2513 0 : if(fixed_blocksize .or. &
2514 : & (wfoptalg >= 0 &
2515 : & .and. nonlop_smem < fourwf_mem &
2516 : & .and. fourwf_mem >= getghc_wmem &
2517 : & .and. .not. fixed_fourwf_blocks &
2518 : & .and. nfft_blocks < ndat &
2519 : & .and. ndat_try > 1)) then
2520 : ! Fourwf work memory requirement is higher, split here
2521 0 : if(nfft_blocks == ndat) cycle ! Can't split more than ndat
2522 :
2523 : if(i>1 .and. .not. print_and_exit) then
2524 0 : do while(ndat_try <= (ndat/nfft_blocks + modulo(ndat,nfft_blocks)))
2525 0 : nfft_blocks=nfft_blocks+1
2526 : end do
2527 : ndat_try = (ndat/nfft_blocks + modulo(ndat,nfft_blocks))
2528 : end if
2529 : #ifdef HAVE_GPU
2530 : call gpu_fft_get_estimate_work_size(3, c_loc(t_fft), FFT_Z2Z, ndat/nfft_blocks, fourwf_smem);
2531 : #endif
2532 0 : getghc_wmem = getghc_ompgpu_work_mem(gs_hamk, ndat, nfft_blocks)
2533 : fourwf_wmem = int(2, c_size_t) * dp * gs_hamk%n4 * gs_hamk%n5 * gs_hamk%n6 &
2534 0 : & * (ndat/nfft_blocks + modulo(ndat,nfft_blocks))
2535 0 : fourwf_mem = fourwf_wmem + fourwf_smem
2536 : else
2537 : ! Gemm nonlop static memory requirement is higher, split here
2538 0 : if(i>1 .and. .not. print_and_exit) blocksize = blocksize + 1
2539 0 : if(modulo(nprocs,blocksize)/=0 .and. use_distrib) cycle
2540 0 : if(nprocs < blocksize .and. use_distrib) cycle
2541 : !FIXME : Skipping uneven blocksize <=5 if using MPI distrib, as the amount of GPU per node is even usually
2542 : !For example, with 3 nodes of 4 GPU, we don't want to have a blocksize of 3 as
2543 : !it would generate 4 comms-block, with 2 inter-node comms.
2544 : !While using a blocksize of 4 would generate 3 comms, one for each node, leading to less MPI comms
2545 0 : if(i>1 .and. modulo(blocksize,2)/=0 .and. use_distrib .and. .not. print_and_exit) cycle
2546 0 : if(i>1) nblocks=nprocs/blocksize
2547 :
2548 : nonlop_smem = gemm_nonlop_ompgpu_static_mem(npw_fft,gs_hamk%indlmn,gs_hamk%nattyp,&
2549 0 : & gs_hamk%ntypat,blocksize,ndgxdt,use_distrib)
2550 : end if
2551 : end if
2552 :
2553 : ! Bandpp~ndat sized buffer memory requirements are higher, split there
2554 0 : sum_mem = nonlop_smem + gs_ham_smem
2555 0 : sum_bandpp_mem = getghc_wmem
2556 :
2557 0 : if(wfoptalg>=0) then
2558 0 : sum_mem = sum_mem+getghc_wmem+updrho_wmem+prep_nonlop_wmem+fourwf_smem+hegvd_mem
2559 : else
2560 0 : sum_mem = sum_mem+nonlop_wmem+prep_nonlop_wmem
2561 : end if
2562 :
2563 0 : if(wfoptalg==111 .or. wfoptalg==112) then
2564 0 : sum_mem = sum_mem + invovl_smem+invovl_wmem+chebfiMem(1)+chebfiMem(2)+localMem
2565 0 : sum_bandpp_mem = sum_bandpp_mem + invovl_wmem
2566 : end if
2567 :
2568 0 : if(wfoptalg==114) then
2569 0 : sum_mem = sum_mem + lobpcgMem(1)+lobpcgMem(2)+localMem
2570 : end if
2571 :
2572 0 : if(sum_mem < free_mem .or. print_and_exit) exit
2573 :
2574 : end do
2575 :
2576 : ! Corner case : not enough GPU memory in forstrnps for forces and stress computation.
2577 : ! By default, forces and stress are computed in one gemm_nonlop call using choice==23
2578 : ! This translates to have various arrays sized by ndgxdt == 9 (6 for stress, 3 for forces)
2579 : ! To try circumventing the lack of GPU memory in that case, we may compute stress and forces
2580 : ! separately so arrays will be sized after ndgxdt=6 at most instead.
2581 0 : if(sum_mem > free_mem .and. optfor > 0 .and. optstr > 0 .and. wfoptalg < 0) then
2582 0 : ndgxdt = 6 ! number of derivatives for stress
2583 : nonlop_wmem = gemm_nonlop_ompgpu_work_mem(gs_hamk%istwf_k, ndat, ndgxdt, npw_fft,&
2584 0 : & gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%ntypat, gs_hamk%lmnmax, signs, wfoptalg)
2585 0 : blocksize=1
2586 : ! Same loop as above, simplified to forstrnps use case
2587 0 : do i=1,nprocs
2588 : ! Gemm nonlop static memory requirement is higher, split here
2589 0 : if(i>1 .and. .not. print_and_exit) blocksize = blocksize + 1
2590 0 : if(modulo(nprocs,blocksize)/=0 .and. use_distrib) cycle
2591 : !FIXME : Skipping uneven blocksize <=5 if using MPI distrib, as the amount of GPU per node is even usually
2592 : !For example, with 3 nodes of 4 GPU, we don't want to have a blocksize of 3 as
2593 : !it would generate 4 comms-block, with 2 inter-node comms.
2594 : !While using a blocksize of 4 would generate 3 comms, one for each node, leading to less MPI comms
2595 0 : if(i>1 .and. modulo(blocksize,2)/=0 .and. use_distrib .and. .not. print_and_exit) cycle
2596 0 : if(i>1) nblocks=nprocs/blocksize
2597 :
2598 : nonlop_smem = gemm_nonlop_ompgpu_static_mem(npw_fft,gs_hamk%indlmn,gs_hamk%nattyp,gs_hamk%ntypat,&
2599 0 : & blocksize,ndgxdt, use_distrib)
2600 0 : sum_mem = nonlop_smem + gs_ham_smem + nonlop_wmem + prep_nonlop_wmem
2601 :
2602 0 : if(sum_mem < free_mem) then
2603 0 : gemm_nonlop_split_choice23 = .true.
2604 0 : exit
2605 : end if
2606 : end do
2607 : end if
2608 :
2609 : ! Quickfix : sometimes, we may run out of GPU memory when computing forces/stresses because of fragmentation.
2610 : ! We try to reduce the risk by forcing even more blocking:
2611 0 : if((wfoptalg < 0 .and. (optfor > 0 .or. optstr > 0)) .and. sum_mem > 0.95*free_mem) then
2612 0 : if(blocksize > 5) then
2613 0 : if(.not. gemm_nonlop_split_choice23) then
2614 0 : gemm_nonlop_split_choice23 = .true.
2615 0 : ndgxdt=6
2616 : else
2617 0 : blocksize=blocksize*1.5
2618 0 : blocksize=min(nprojs,blocksize)
2619 : end if
2620 : else
2621 0 : blocksize=blocksize*1.5
2622 : end if
2623 :
2624 : nonlop_smem = gemm_nonlop_ompgpu_static_mem(npw_fft,gs_hamk%indlmn,gs_hamk%nattyp,gs_hamk%ntypat,&
2625 0 : & blocksize,ndgxdt,use_distrib)
2626 0 : sum_mem = nonlop_smem + gs_ham_smem + nonlop_wmem + prep_nonlop_wmem
2627 : end if
2628 :
2629 0 : if(.not. l_disable_output .or. sum_other_mem > free_mem .or. sum_mem > free_mem) then
2630 0 : write(std_out,'(A,I3,A)') "GPU memory consumption estimate per MPI task for K-point ",ikpt,":"
2631 0 : if(blocksize>1) then
2632 0 : if(use_distrib) then
2633 0 : write(std_out,'(A,I3,A,I3,A)') "MPI distribution of GEMM nonlop projectors using ",&
2634 0 : & nblocks, " blocks of ", blocksize, " MPI tasks."
2635 : else
2636 0 : write(std_out,'(A,I3,A)') "Local slicing of GEMM nonlop projectors using ",&
2637 0 : & blocksize, " blocks."
2638 : end if
2639 : end if
2640 0 : if(nfft_blocks>1 .and. wfoptalg>=0) then
2641 0 : write(std_out,'(A,I3,A)') "Local slicing of FFT work array using ",&
2642 0 : & nfft_blocks, " blocks."
2643 : end if
2644 0 : write(std_out,'(A,F10.3,1x,A)') " Considered available memory : ", real(free_mem)/(1024*1024), "MiB"
2645 0 : write(std_out,'(A)')
2646 0 : write(std_out,'(A)') "| Buffers governed by blocking/slicing |"
2647 0 : write(std_out,'(A)') "|:--------------------------|---------:|-------------:|---------------:|"
2648 0 : write(std_out,'(A,I4,A,F10.2,1x,A)') "| gemm_nonlop_projectors | ", blocksize, " blk | npw,*natom* | ", real(nonlop_smem,dp)/(1024*1024), "MiB |"
2649 0 : if(wfoptalg>=0) then
2650 0 : write(std_out,'(A,I4,A,F10.2,1x,A)') "| fourwf (fofr work array) | ", nfft_blocks, " blk | npw,*bandpp* | ", real(fourwf_wmem,dp)/(1024*1024), "MiB |"
2651 0 : write(std_out,'(A,F10.2,1x,A)') "| xFFT~internal buffers | NA | NA | ", real(fourwf_smem,dp)/(1024*1024), "MiB |"
2652 : end if
2653 :
2654 0 : if(sum_other_mem > free_mem) then
2655 0 : write(std_out,'(A)')
2656 0 : write(std_out,'(A)') "/!\ No slicing attempted as other arrays are too big to fit"
2657 0 : write(std_out,'(A)')
2658 : end if
2659 :
2660 0 : write(std_out,'(A)')
2661 0 : write(std_out,'(A)') "| Static buffers, computed once and permanently on card |"
2662 0 : write(std_out,'(A)') "|:-------------------------|--------------:|----------------:|"
2663 : ! CHEBFI2 or SLICE
2664 0 : if(wfoptalg==111 .or. wfoptalg==112) then
2665 0 : write(std_out,'(A,F10.2,1x,A)') "| invovl (mkinvovl) | natom | ", real(invovl_smem,dp)/(1024*1024), "MiB |"
2666 0 : write(std_out,'(A,F10.2,1x,A)') "| chebfi2 | npw | ", real(chebfiMem(1))/(1024*1024), "MiB |"
2667 : end if
2668 :
2669 : ! LOBPCG2
2670 0 : if(wfoptalg==114) then
2671 0 : write(std_out,'(A,F10.2,1x,A)') "| lobpcg2 | npw | ", real(lobpcgMem(1))/(1024*1024), "MiB |"
2672 : end if
2673 :
2674 0 : write(std_out,'(A,F10.2,1x,A)') "| hamiltonian arrays | npw | ", real(gs_ham_smem)/(1024*1024), "MiB |"
2675 :
2676 0 : write(std_out,'(A)')
2677 0 : write(std_out,'(A)') "| Work buffers (mostly sized after bandpp or nblock_lobpcg) |"
2678 0 : write(std_out,'(A)') "|:-------------------------|--------------:|----------------:|"
2679 : ! getghc (any diago algorithm)
2680 0 : if(wfoptalg>=0) then
2681 0 : if(getghc_wmem /= fourwf_wmem) then
2682 0 : write(std_out,'(A,F10.2,1x,A)') "| gemm_nonlop | bandpp | ", real(getghc_wmem,dp)/(1024*1024), "MiB |"
2683 : end if
2684 0 : write(std_out,'(A,F10.2,1x,A)') "| mkrho~vtowfk_extra | npw,bandpp | ", real(updrho_wmem,dp)/(1024*1024), "MiB |"
2685 0 : write(std_out,'(A,F10.2,1x,A)') "| hegvd | bandpp | ", real(hegvd_mem,dp)/(1024*1024), "MiB |"
2686 : else
2687 0 : write(std_out,'(A,F10.2,1x,A)') "| gemm_nonlop | natom,bandpp | ", real(nonlop_wmem,dp)/(1024*1024), "MiB |"
2688 : end if
2689 0 : if(paral_kgb==1) then
2690 0 : write(std_out,'(A,F10.2,1x,A)') "| prep_nonlop | npw,bandpp | ", real(prep_nonlop_wmem,dp)/(1024*1024), "MiB |"
2691 : end if
2692 :
2693 : ! CHEBFI2 or SLICE
2694 0 : if(wfoptalg==111 .or. wfoptalg==112) then
2695 0 : write(std_out,'(A,F10.2,1x,A)') "| invovl | natom,bandpp | ", real(invovl_wmem,dp)/(1024*1024), "MiB |"
2696 0 : write(std_out,'(A,F10.2,1x,A)') "| chebfi2 (RR buffers) | nband | ", real(chebfiMem(2))/(1024*1024), "MiB |"
2697 0 : write(std_out,'(A,F10.2,1x,A)') "| chebfiwf (cg,resid,eig) | npw,nband | ", real(localMem)/(1024*1024), "MiB |"
2698 : end if
2699 :
2700 : ! LOBPCG2
2701 0 : if(wfoptalg==114) then
2702 0 : write(std_out,'(A,F10.2,1x,A)') "| lobpcg2 (RR buffers) | bandpp | ", real(lobpcgMem(2))/(1024*1024), "MiB |"
2703 0 : write(std_out,'(A,F10.2,1x,A)') "| lobpcgwf (cg,resid,eig) | npw,nband | ", real(localMem)/(1024*1024), "MiB |"
2704 : end if
2705 :
2706 0 : write(std_out,'(A)')
2707 0 : if(sum_other_mem > free_mem) then
2708 0 : write(std_out,'(A,F10.2,1x,A)') "Sum : ", real(sum_other_mem)/(1024*1024), "MiB"
2709 : else
2710 0 : write(std_out,'(A,F10.2,1x,A)') "Sum : ", real(sum_mem)/(1024*1024), "MiB"
2711 : end if
2712 0 : write(std_out,'(A)')
2713 0 : flush(std_out)
2714 :
2715 : end if
2716 :
2717 0 : if(rank==0) then
2718 0 : if(sum_other_mem > free_mem) then
2719 : write(message,'(3a)') &
2720 0 : & ' Your case is too big to fit in GPU memory regardless of possible array optimizations in fourwf and GEMM nonlop.',ch10,&
2721 0 : & ' Action : run on more nodes and/or increase nblock_lobpcg if using LOBPCG.'
2722 0 : ABI_ERROR(message)
2723 : end if
2724 0 : if(sum_mem > free_mem) then
2725 0 : if(l_warn_on_fail) then
2726 0 : ABI_WARNING("It seems the test case you're trying to run is too big to run with given GPU resources !")
2727 : else
2728 0 : ABI_ERROR("It seems the test case you're trying to run is too big to run with given GPU resources !")
2729 : end if
2730 : end if
2731 : end if
2732 : !call xmpi_barrier(xmpi_world)
2733 :
2734 : end subroutine get_gemm_nonlop_ompgpu_blocksize
2735 : !!***
2736 :
2737 : end module m_common
2738 : !!***
|
