Branch data 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 : 84096 : subroutine scprqt(choice,cpus,deltae,diffor,maxmagsph,difmagsph,dtset,&
188 : 84096 : & eigen,etotal,favg,fcart,fermie,fermih,fname_eig,filnam1,initGS,&
189 : 84096 : & iscf,istep,istep_fock_outer,istep_mix,kpt,maxfor,moved_atm_inside,mpi_enreg,&
190 : 84096 : & nband,nkpt,nstep,occ,optres,&
191 : 84096 : & prtfor,prtxml,quit,res2,resid,residm,response,tollist,usepaw,&
192 : 84096 : & 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 : 168192 : real(dp) :: residm_band(dtset%mband,dtset%nsppol), f_tmp(3)
240 : :
241 : : ! *********************************************************************
242 : :
243 : : DBG_ENTER("COLL")
244 : :
245 : 84096 : my_rank = mpi_enreg%me_cell
246 : :
247 : 84096 : quit=0; conv_retcode=0
248 : 84096 : usefock=dtset%usefock
249 : 84096 : nnsclohf=dtset%nnsclohf
250 : 84096 : use_dpfft = .False.
251 : :
252 : 84096 : tmagnet=0
253 [ + + + + : 60756 : if(response==0.and.(iscf>0.or.iscf==-3).and.dtset%nsppol==2.and.dtset%occopt>2) tmagnet=1
+ + ]
254 [ + + + + : 84096 : if((iscf>0.or.iscf==-3).and.(dtset%nsppol==2.or.dtset%nspden>1)) tmagnet=2
+ + ]
255 : :
256 : 84096 : ishift=0
257 [ + + + + ]: 1014261 : residm_band = zero
258 [ + + ]: 176042 : do isppol=1, dtset%nsppol
259 [ + + ]: 2073471 : do ikpt=1, nkpt
260 [ + + ]: 12784159 : do iband=1, nband(ikpt+(isppol-1)*nkpt)
261 : 10794784 : ishift = ishift+1
262 : 12692213 : residm_band(iband, isppol) = max (resid(ishift), residm_band(iband, isppol))
263 : : end do
264 : : end do
265 : : end do
266 : :
267 : 93372 : 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 : 9276 : tolwfr=tollist(2)
274 : 9276 : toldff=tollist(3)
275 : 9276 : toldfe=tollist(4)
276 : 9276 : tolvrs=tollist(6)
277 : 9276 : tolrff=tollist(7)
278 : 9276 : vdw_df_threshold=tollist(8)
279 : 9276 : toldmag=tollist(9)
280 : 9276 : ttolwfr=0 ; ttoldff=0 ; ttoldfe=0 ; ttolvrs=0; ttolrff=0; ttoldmag=0 ;
281 [ + + ]: 9276 : if(abs(tolwfr)>tiny(zero))ttolwfr=1
282 [ + + ]: 9276 : if(abs(toldff)>tiny(zero))ttoldff=1
283 [ + + ]: 9276 : if(abs(tolrff)>tiny(zero))ttolrff=1
284 [ + + ]: 9276 : if(abs(toldfe)>tiny(zero))ttoldfe=1
285 [ + + ]: 9276 : if(abs(tolvrs)>tiny(zero))ttolvrs=1
286 [ + + ]: 9276 : if(abs(toldmag)>tiny(zero))ttoldmag=1
287 : : ! If non-scf calculations, tolwfr must be defined
288 [ + + - + ]: 9276 : 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 [ + + + + : 9276 : 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 [ - + ]: 9276 : if (dtset%usewvl == 1) then
316 : 0 : write(colname, "(A)") "grdnorm "
317 : : else
318 : 9276 : write(colname, "(A)") "residm "
319 : : end if
320 [ + + + + : 9276 : if (nstep>0 .and. (iscf>=0 .or.iscf==-3) .and. dtset%prtstm==0) then
+ + ]
321 [ + + + + ]: 8638 : 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 [ + + ]: 7859 : if(response==0)then
341 [ + + ]: 5629 : if (prtfor==0) then
342 [ + + ]: 2713 : if (optres==0) then
343 : 2110 : write(message, '(4a)' ) ch10,&
344 : 4220 : ' iter Etot(hartree) deltaE(h) ', colname, ' vres2'
345 : : else
346 : 603 : write(message, '(4a)' ) ch10,&
347 : 1206 : ' 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 [ + + + + ]: 2230 : 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 [ + - ]: 2122 : if (optres==0) then
369 : 2122 : write(message, '(4a)' ) ch10,&
370 : 4244 : ' 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 : 8638 : 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 [ + + ]: 43190 : [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 [ + + ]: 60466 : real_fmt="(es8.2)", dict_key="tolerances", ignore=zero)
388 : :
389 : 8638 : call ydoc%write_and_free(ab_out, newline=.False.)
390 : 8638 : call wrtout(ab_out, message)
391 : : end if
392 : :
393 : : case (2)
394 : :
395 : : ! Examine tolerance criteria
396 : 65544 : tolwfr=tollist(2)
397 : 65544 : toldff=tollist(3)
398 : 65544 : toldfe=tollist(4)
399 : 65544 : tolvrs=tollist(6)
400 : 65544 : tolrff=tollist(7)
401 : 65544 : vdw_df_threshold=tollist(8)
402 : 65544 : toldmag=tollist(9)
403 : 65544 : ttolwfr=0 ; ttoldff=0 ; ttoldfe=0 ; ttolvrs=0; ttolrff=0; ttoldmag=0;
404 : : !ttolwfr=0 ; ttoldff=0 ; ttoldfe=0 ; ttolvrs=0; ttolrff=0;
405 [ + + ]: 65544 : if(abs(tolwfr)>tiny(0.0_dp))ttolwfr=1
406 [ + + ]: 65544 : if(abs(toldff)>tiny(0.0_dp))ttoldff=1
407 [ + + ]: 65544 : if(abs(tolrff)>tiny(0.0_dp))ttolrff=1
408 [ + + ]: 65544 : if(abs(toldfe)>tiny(0.0_dp))ttoldfe=1
409 [ + + ]: 65544 : if(abs(tolvrs)>tiny(0.0_dp))ttolvrs=1
410 [ + + ]: 65544 : 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 [ + + + + ]: 73640 : if (dtset%prtvol>=10) call wrtout([std_out, ab_out], ' ')
414 : :
415 : : ! Calculate up and down charge and magnetization
416 [ - + ]: 65544 : 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 [ + + ]: 65544 : 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 [ + + ]: 65544 : if(dtset%prtstm==0)then
462 : 65541 : format_istep='i3'
463 [ + + ]: 65541 : if(istep>99)format_istep='i5'
464 [ - + ]: 65541 : if(istep>9999)format_istep='i7'
465 : : !if(tmagnet==1)then
466 [ + + + + ]: 65541 : 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 : 58007 : firstchar=' '
485 [ + + + + ]: 58007 : if (response/=0.and.istep==1) firstchar="-"
486 [ + + ]: 58007 : if (response==0) then
487 [ + + ]: 39375 : if (prtfor==0) then
488 : : write(message, '(2a,'//format_istep//',1p,g22.14,3es10.3)' ) &
489 : 19697 : 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 [ + + + + ]: 18632 : 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 : 17064 : 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 : 65541 : call wrtout(ab_out,message)
506 : :
507 [ + + ]: 65541 : if(mpi_enreg%paral_pert==1) then
508 : 337 : call wrtout(std_out, message,'PERS')
509 [ + - ]: 65204 : elseif(mpi_enreg%paral_pert==0) then
510 : 65204 : 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 [ + + + - : 65544 : 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 [ + + + + : 65544 : 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 [ + + ]: 65544 : if(response==0)then
565 : 46912 : write(message, '(a,1p,e15.7,a)' ) ' scprqt: <Vxc>=',vxcavg,' Hartree'
566 : 46912 : call wrtout(std_out,message)
567 : : end if
568 : :
569 : : ! Check whether exiting was required by the user.
570 [ + + ]: 65544 : openexit=1 ; if(dtset%chkexit==0) openexit=0
571 : 65544 : call exit_check(cpus,filnam1,iexit,ab_out,mpi_enreg%comm_cell,openexit)
572 [ - + ]: 65544 : 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 [ + + + + : 65544 : (dtset%usedmatpu<0.or.initGS==0))
+ + + + +
+ + - +
- ]
578 : :
579 : : ! Additional stuff for electron/positron
580 [ + + ]: 65544 : if (present(electronpositron)) then
581 [ + + ]: 46912 : 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 [ + + - + ]: 65544 : 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 [ + + + + : 65497 : 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 [ + + ]: 15237 : if (residm < tolwfr) then
604 [ + - ]: 1645 : if (dtset%usewvl == 0) then
605 : 1645 : write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a)' )ch10, &
606 : 3290 : ' 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 [ + + ]: 4935 : call wrtout([std_out, ab_out], message)
612 : 1645 : quit=1
613 : : else
614 : 13592 : use_dpfft = residm < tol7
615 : : end if
616 [ + + ]: 15237 : 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 [ + + ]: 65497 : 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 [ + + ]: 65497 : 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 [ + + ]: 65497 : 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 [ + + ]: 65497 : 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 [ + + + + ]: 65497 : if (ttolvrs==1 .and. .not. noquit) then
834 [ + + ]: 24001 : if (ttolwfr==0) then
835 [ + + ]: 23752 : 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 : 21948 : 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 [ + + ]: 24001 : if(present(rcpaw)) then
866 [ - + - - ]: 7682 : 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 [ + + - + ]: 65497 : 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 [ + + ]: 65497 : if (use_dpfft) then
884 [ - + ]: 18161 : if (fftcore_mixprec == 1) call wrtout(std_out, " Approaching convergence. Activating FFT in double-precision")
885 : 18161 : ii = fftcore_set_mixprec(0)
886 : : end if
887 : :
888 : : ! Additional stuff for rcpaw
889 [ + + ]: 65544 : if(present(rcpaw)) then
890 [ + + ]: 14429 : 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 [ + + ]: 9276 : if(nstep>0) then
902 [ + + ]: 9241 : if (.not. converged()) then
903 : 1936 : conv_retcode = 1
904 : :
905 [ + + - + ]: 1936 : if(iscf>=1 .or. iscf==-3 .or. iscf == 0)then
906 : 1821 : write(message, '(a,a,a,a,i5,a)' ) ch10,&
907 : 1821 : ' scprqt: WARNING -',ch10,&
908 : 3642 : ' nstep=',nstep,' was not enough SCF cycles to converge;'
909 : :
910 : 1821 : write(std_out,'(6a,i0,3a)')ch10,&
911 : 1821 : "--- !ScfConvergenceWarning",ch10,&
912 : 1821 : "message: |",ch10,&
913 : 1821 : ' nstep ',nstep,' was not enough SCF cycles to converge.',ch10,&
914 : 3642 : "..."
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 [ + + ]: 5808 : call wrtout([std_out, ab_out], message)
928 : :
929 [ + + + + ]: 1936 : 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 [ + + ]: 1255 : 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 [ - + ]: 795 : 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 [ + + ]: 795 : 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 [ + + ]: 784 : 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 [ + - ]: 534 : else if(ttolvrs==1)then
967 [ + + ]: 534 : 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 : 360 : ' density residual=',res2,' exceeds tolvrs=',tolvrs,ch10
973 : : end if
974 : : end if
975 [ + + ]: 5808 : call wrtout([std_out, ab_out], message)
976 : :
977 [ - + ]: 1936 : 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 [ + + ]: 7305 : 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 [ + + ]: 9241 : 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 [ + + + - ]: 149640 : ABI_BUG(message)
1033 : : end select
1034 : :
1035 : : ! Additional stuff for the Fock+SCF cycle
1036 [ + + ]: 84096 : if (present(fock)) then
1037 [ + + ]: 53834 : 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 [ + + ]: 84096 : if (present(electronpositron)) then
1046 [ + + ]: 53834 : 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 [ + + ]: 168192 : call flush_unit(ab_out)
1110 : :
1111 : : DBG_EXIT("COLL")
1112 : :
1113 : : contains
1114 : :
1115 : 9241 : 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 : 9241 : loc_conv = .true.
1121 [ + + ]: 9241 : if (ttolwfr==1) then
1122 [ + + ]: 2851 : if (residm > tolwfr) loc_conv=.false.
1123 : : end if
1124 [ + + ]: 9241 : if (ttoldff==1) then
1125 [ + + ]: 2673 : if (diffor > toldff) loc_conv=.false.
1126 : : end if
1127 [ + + ]: 9241 : if (ttoldmag==1) then
1128 [ + - ]: 2 : if (difmagsph > toldmag) loc_conv=.false.
1129 : : end if
1130 [ + + ]: 9241 : if (ttolrff==1) then
1131 [ + + - + ]: 132 : if (diffor > tolrff*maxfor .and. maxfor > tol16) loc_conv=.false.
1132 : : end if
1133 [ + + ]: 9241 : if (ttoldfe==1) then
1134 [ + + ]: 1266 : if (abs(deltae) > toldfe) loc_conv=.false.
1135 : : end if
1136 [ + + ]: 9241 : if (ttolvrs==1) then
1137 [ + + ]: 2341 : if (res2 > tolvrs) loc_conv=.false.
1138 : : end if
1139 : 9241 : converged = loc_conv
1140 : :
1141 : 9241 : 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 : 12896 : 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 : 6448 : bantot=0
1215 [ + + ]: 13505 : do isppol=1,nsppol
1216 [ + + ]: 144722 : do ikpt=1,nkpt
1217 : 138274 : bantot=bantot+dtset%nband(ikpt+(isppol-1)*nkpt)
1218 : : end do
1219 : : end do
1220 : :
1221 [ - + ]: 6448 : 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 : 6448 : 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 : 6448 : 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 : 6448 : k0(:)=0.0_dp
1239 [ + + + + ]: 6448 : 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 [ + + ]: 6448 : if (usepaw==1) then
1244 : 1439 : write(msg,'(2a)') ch10,' Coarse grid specifications (used for wave-functions):'
1245 [ + + ]: 4317 : call wrtout([std_out, ab_out], msg)
1246 : 1439 : call getcut(boxcutc,ecutc_eff,gmet,gsqcutc_eff,dtset%iboxcut,ab_out,k0,ngfftc)
1247 : 1439 : write(msg,'(2a)') ch10,' Fine grid specifications (used for densities):'
1248 [ + + ]: 4317 : call wrtout([std_out, ab_out], msg)
1249 : 1439 : 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 [ + + - + ]: 6448 : 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 : 6448 : 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 : 12822 : subroutine prteigrs(eigen,enunit,fermie,fermih,fname_eig,iout,iscf,kptns,kptopt,mband,nband,&
1322 : 12822 : & 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 [ - + ]: 12822 : 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 [ + + ]: 12822 : if (prteig > 0) then
1353 : 8826 : call wrtout(iout, sjoin(' prteigrs : about to open file ', fname_eig))
1354 [ - + ]: 8826 : if (open_file(fname_eig, msg, newunit=temp_unit, status='unknown', form='formatted') /= 0) then
1355 : 0 : ABI_ERROR(msg)
1356 : : end if
1357 : 8826 : rewind(temp_unit) ! always rewind disk file and print latest eigenvalues
1358 : : end if
1359 : :
1360 : 12822 : kind_of_output= ' Eigenvalues '
1361 [ + + ]: 12822 : if(option==4) kind_of_output=' Expectation of eigenvalue derivatives'
1362 [ + + ]: 12822 : if(option==5) kind_of_output=' Fan corrections to eigenvalues at T=0'
1363 [ + + ]: 12822 : if(option==6) kind_of_output=' DDW corrections to eigenvalues at T=0'
1364 [ + + ]: 12822 : if(option==7) kind_of_output=' Fan+DDW corrs to eigenvalues at T=0'
1365 : :
1366 : 12822 : 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 [ + + + + : 12822 : if( (prtvol==0.or.prtvol==1) .and. (iscf/=-2 .or. kptopt>0) .and. nkpt_eff>nkpt_max)nkpt_eff=nkpt_max
+ + ]
1372 [ + + + + : 12822 : 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 [ + + ]: 12822 : if(option==1 .or. (option>=3 .and. option<=7))then
1375 : :
1376 [ + + ]: 37164 : do ienunit=0,1
1377 : :
1378 [ + + + + ]: 24776 : if (enunit==1 .and. ienunit==0)cycle
1379 [ + + + + ]: 24498 : 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 [ + + ]: 14940 : if (ienunit==0)then
1384 : 12110 : convrt=one
1385 : 12110 : strunit1='hartree'
1386 : 12110 : strunit2='hartree'
1387 : : end if
1388 [ + + ]: 14940 : if (ienunit==1)then
1389 : 2830 : convrt=Ha_eV
1390 : 2830 : strunit1=' eV '
1391 : 2830 : strunit2='eV '
1392 : : end if
1393 : :
1394 : 14940 : band_index=0
1395 : :
1396 [ + + ]: 14940 : if(ienunit==0)then ! XG20140730 I do not know why this is only done when ienunit==0
1397 : 12110 : tmetal=0
1398 [ + + + + : 12110 : 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 [ + + + + : 14940 : if(iscf>=0 .and. (ienunit==0 .or. option==1))then
+ - ]
1420 [ + + ]: 13078 : 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 : 13076 : ' Fermi (or HOMO) energy (',trim(strunit2),') =',convrt*fermie,' Average Vxc (',trim(strunit2),')=',convrt*vxcavg
1427 : : end if
1428 : 13078 : call wrtout(iout,msg)
1429 [ + + ]: 13078 : 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 [ + + + + ]: 14940 : if(iscf>=0 .and. ienunit==0)then ! For historical reasons
1434 [ + + ]: 10810 : 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 [ + + ]: 44852 : do isppol=1,nsppol
1445 : :
1446 [ - + - + ]: 17524 : ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9"
1447 [ + + + + ]: 17524 : 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 [ + + + - ]: 14940 : 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 : 12356 : write(msg, '(4a,'//ikpt_fmt//',2x,a)' )trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points:'
1453 : : end if
1454 : 17524 : call wrtout(iout,msg)
1455 [ + + ]: 17524 : if (prteig > 0) call wrtout(temp_unit,msg)
1456 : :
1457 [ + + ]: 17524 : if(ienunit==0)then
1458 [ + + ]: 13894 : if(option>=4 .and. option<=7)then
1459 : 559 : msg = ' (in case of degenerate eigenvalues, averaged derivative)'
1460 : 559 : call wrtout(iout,msg)
1461 [ + + ]: 559 : if (prteig > 0) call wrtout(temp_unit,msg)
1462 : : end if
1463 : : end if
1464 : :
1465 [ + + ]: 194984 : do ikpt=1,nkpt
1466 : 152684 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1467 [ - + - + ]: 152684 : ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9"
1468 [ - + - + ]: 152684 : ibnd_fmt="i3" ; if(nband_k>=1000)ibnd_fmt="i6" ; if(nband_k>=1000000)ibnd_fmt="i9"
1469 [ + + ]: 152684 : if(ikpt<=nkpt_eff)then
1470 : : write(msg, '(a,'//ikpt_fmt//',a,'//ibnd_fmt//',a,f9.5,a,3f8.4,a)' ) &
1471 [ + + ]: 266136 : & ' kpt#',ikpt,', nband=',nband_k,', wtk=',wtk(ikpt)+tol10,', kpt=',kptns(1:3,ikpt)+tol10,' (reduced coord)'
1472 : 66534 : call wrtout(iout,msg)
1473 [ + + ]: 66534 : if (prteig > 0) call wrtout(temp_unit,msg)
1474 [ + + ]: 166380 : do ii=0,(nband_k-1)/8
1475 [ + - + + ]: 616798 : write(msg, '(8(f10.5,1x))' ) (convrt*eigen(iband+band_index), iband=1+ii*8,min(nband_k,8+ii*8))
1476 : 99846 : call wrtout(iout,msg)
1477 [ + + ]: 166380 : if (prteig > 0) call wrtout(temp_unit,msg)
1478 : : end do
1479 [ + + + + : 66534 : 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 : 170208 : 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 [ + + ]: 12822 : if (prteig > 0) close (temp_unit)
1580 : :
1581 : 12822 : 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 : 5281 : 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 [ + + + + : 1380 : 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 [ + + ]: 5281 : optdc=-1;ipositron=merge(0,2,dtset%positron==0)
1636 [ + + + + ]: 5281 : if (abs(energies%e_ewald)<1.e-15_dp.and.abs(energies%e_hartree)<1.e-15_dp) ipositron=1
1637 : 5281 : call energies%eval_eint(dtset,usepaw,optdc,etotal,etotaldc)
1638 : :
1639 : 5281 : call entropy(dtset,energies)
1640 : 5281 : 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 [ + + + + : 20559 : & 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 [ + + - + : 5281 : & .and.(dtset%dmft_solv==6.or.dtset%dmft_solv==7)) testdmft=.false.
- - - - ]
1648 [ + + ]: 5281 : if(testdmft) eent=zero
1649 : :
1650 : 5281 : etotal = etotal + eent
1651 : 5281 : etotaldc = etotaldc + eent
1652 : :
1653 : : !============= Printing of Etotal by direct scheme ===========
1654 : :
1655 [ - + ]: 5281 : if (dtset%icoulomb == 1) then
1656 : 0 : eneName = "Ion-ion energy"
1657 : : else
1658 : 5281 : eneName = "Ewald energy"
1659 : : end if
1660 : 5281 : enevalue = energies%e_ewald
1661 : :
1662 : :
1663 [ + - ]: 5281 : if (optdc==0.or.optdc==2) then
1664 : :
1665 [ + + ]: 5281 : if (directE_avail) then
1666 : 5280 : info = 'Components of total free energy in Hartree'
1667 [ + + ]: 5280 : if(testdmft) info = 'Components of total energy in Hartree'
1668 [ + - ]: 5280 : edoc = yamldoc_open('EnergyTerms', info=trim(adjustl(info)), width=20, real_fmt='(es21.14)')
1669 : 5280 : call edoc%add_real('kinetic', energies%e_kinetic)
1670 [ + + ]: 5280 : if(abs(energies%e_extfpmd)>tiny(zero)) call edoc%add_real('extfpmd',energies%e_extfpmd)
1671 [ + + ]: 5280 : if (ipositron/=1) then
1672 : 5275 : 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 [ + + ]: 5275 : if(usepaw==1)exc_semilocal=exc_semilocal+energies%e_fock
1677 : 5275 : call edoc%add_real('hartree', energies%e_hartree)
1678 : 5275 : call edoc%add_real('xc', exc_semilocal)
1679 : 5275 : call edoc%add_real(eneName, enevalue)
1680 : 5275 : 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 : 5280 : call edoc%add_real('local_psp', energies%e_localpsp)
1688 [ + + ]: 5280 : 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 [ + + ]: 1379 : 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 [ + + ]: 1376 : 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 : 1375 : 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 : : ! Suppress output of e_nucdip in the npband case, due to problematic direct
1713 : : ! call to getghc_nucdip in m_vtowfk in this case. SCF cycle, which depends
1714 : : ! on getghc itself, is ok (parallel bands are distributed properly there
1715 : : ! before getghc_nucdip is called
1716 [ + + + + : 56260 : if (ANY(ABS(dtset%nucdipmom)>tol8) .AND. (dtset%npband .EQ. 1) ) then
+ + + + +
+ ]
1717 : 24 : call edoc%add_real('nucl. magn. dipoles',energies%e_nucdip)
1718 : : end if
1719 [ + + + - ]: 5280 : if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).and.ipositron/=1) then
1720 : 4 : call edoc%add_real('VdWaals_dft_d', energies%e_vdw_dftd)
1721 : : end if
1722 [ + + ]: 5280 : if (dtset%nzchempot>=1) then
1723 : 7 : call edoc%add_real('chem_potential', energies%e_chempot)
1724 : : end if
1725 [ + + ]: 5280 : if (dtset%usedmft==1) then
1726 : 152 : call edoc%add_real('interaction', energies%e_hu)
1727 : 152 : call edoc%add_real('-double_counting', -energies%e_dc)
1728 : : end if
1729 [ + + + + ]: 5280 : if(dtset%occopt>=3.and.dtset%occopt<=8.and.ipositron==0) then
1730 : 1276 : call edoc%add_real('internal', etotal-eent)
1731 [ + + ]: 1276 : if(.not.testdmft) then
1732 : 1159 : call edoc%add_real('-kT*entropy', eent)
1733 : : end if
1734 [ + + ]: 3998 : else if (ipositron/=0) then
1735 [ + + ]: 14 : if (dtset%occopt>=3.and.dtset%occopt<=8) then
1736 : 6 : call edoc%add_real('-kT*entropy', eent)
1737 : : end if
1738 : : !write(msg, '(3a,es21.14,a)' ) &
1739 : : ! ' >>> ',EPName(ipositron),' E= ',etotal-energies%e0_electronpositron -energies%e_electronpositron,ch10
1740 : : !call wrtout(iout,msg)
1741 : : !write(msg, '(3a,es21.14,2a,es21.14)' ) &
1742 : : ! ' ',EPName(3-ipositron),' ener.= ',energies%e0_electronpositron,ch10,&
1743 : : ! ' EP interaction E= ' ,energies%e_electronpositron
1744 : : !call wrtout(iout,msg)
1745 [ + + ]: 14 : if(ipositron == 1) then
1746 : 5 : call edoc%add_real('positronic', etotal - energies%e0_electronpositron-energies%e_electronpositron)
1747 : 5 : call edoc%add_real('electronic', energies%e0_electronpositron)
1748 : : else
1749 : 9 : call edoc%add_real('electronic', etotal- energies%e0_electronpositron-energies%e_electronpositron)
1750 : 9 : call edoc%add_real('positronic', energies%e0_electronpositron)
1751 : : end if
1752 : 14 : call edoc%add_real('electron_positron_interaction', energies%e_electronpositron)
1753 : : end if
1754 : : if ((dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. &
1755 [ + + + - ]: 5280 : dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17) .and.ipositron/=1) then
1756 : 25 : call edoc%add_real('electric', energies%e_elecfield)
1757 : 25 : call edoc%add_real('kohn_sham', etotal-energies%e_elecfield)
1758 : : end if
1759 : 5280 : call edoc%add_real('total_energy', etotal)
1760 : :
1761 : : else
1762 : : write(msg, '(9a)' ) &
1763 : 1 : ' COMMENT: ',ch10,&
1764 : 1 : ' "Direct" decomposition of total free energy cannot be printed out !!!',ch10,&
1765 : 1 : ' PAW contribution due to spin-orbit coupling cannot be evaluated',ch10,&
1766 : 1 : ' without the knowledge of imaginary part of Rhoij atomic occupancies',ch10,&
1767 : 2 : ' (computed only when pawcpxocc=2).'
1768 : 1 : call wrtout(iout,msg)
1769 : : end if
1770 : : end if
1771 : : !============= Printing of Etotal by double-counting scheme ===========
1772 : :
1773 [ + + ]: 5281 : if (optdc>=1) then
1774 : :
1775 : 1505 : info = '"Double-counting" decomposition of free energy'
1776 [ + + ]: 1505 : if(testdmft) info = '"Double-counting" decomposition of internal energy'
1777 : : dc_edoc = yamldoc_open('EnergyTermsDC', info=trim(adjustl(info)), &
1778 [ + - ]: 1505 : width=20, real_fmt="(es21.14)")
1779 : 1505 : call dc_edoc%add_real('band_energy', energies%e_eigenvalues)
1780 [ + + ]: 1505 : if(abs(energies%e_extfpmd)>tiny(0.0_dp)) then
1781 : 7 : call dc_edoc%add_real('extfpmd_dc',energies%edc_extfpmd)
1782 : : end if
1783 [ + + ]: 1505 : if (ipositron/=1) then
1784 : : !write(msg, '(2(a,es21.14,a),a,es21.14)' ) &
1785 : : ! ' '//eneName//' =',enevalue,ch10,&
1786 : : ! ' PspCore energy = ',energies%e_corepsp-energies%e_corepspdc,ch10,&
1787 : : ! ' Dble-C XC-energy= ',-energies%e_hartree+energies%e_xc-energies%e_xcdc -energies%e_fock0 + &
1788 : : ! energies%e_hybcomp_E0-energies%e_hybcomp_v0
1789 : : !call wrtout(iout,msg)
1790 : 1500 : call dc_edoc%add_real(eneName, enevalue)
1791 : 1500 : call dc_edoc%add_real('psp_core', energies%e_corepsp-energies%e_corepspdc)
1792 : : call dc_edoc%add_real('xc_dc', -energies%e_hartree+energies%e_xc-energies%e_xcdc - energies%e_fock0 + &
1793 : 1500 : energies%e_hybcomp_E0-energies%e_hybcomp_v0)
1794 : : end if
1795 : : if ((dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. &
1796 [ + + + - ]: 1505 : dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17).and.ipositron/=1) then
1797 : 9 : call dc_edoc%add_real('electric_field', energies%e_elecfield)
1798 : : end if
1799 [ + + ]: 1505 : if (usepaw==1) then
1800 [ + + ]: 1378 : if (dtset%use_rcpaw/=0) then
1801 : 3 : call dc_edoc%add_real('PAW spherical_terms', energies%paw%epaw_dc)
1802 : 3 : call dc_edoc%add_real('PAW core dc', energies%paw%epaw_core_dc)
1803 [ + + ]: 1375 : else if (dtset%paw_add_core==1) then
1804 : 1 : call dc_edoc%add_real('PAW spherical_terms', energies%paw%epaw_dc-energies%paw%epaw_core)
1805 : 1 : call dc_edoc%add_real('PAW core', energies%paw%epaw_core)
1806 : : else
1807 : 1374 : call dc_edoc%add_real('PAW spherical_terms', energies%paw%epaw_dc)
1808 : : end if
1809 : : end if
1810 [ - + - - ]: 1505 : if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).and.ipositron/=1) then
1811 : 0 : call dc_edoc%add_real('VdWaals_dft_d', energies%e_vdw_dftd)
1812 : : end if
1813 [ - + ]: 1505 : if (dtset%nzchempot>=1) then
1814 : 0 : call dc_edoc%add_real('chem_potential', energies%e_chempot)
1815 : : end if
1816 [ + + ]: 1505 : if (dtset%usedmft==1) then
1817 : 152 : call dc_edoc%add_real('interaction', energies%e_hu)
1818 : 152 : call dc_edoc%add_real('-double_counting', -energies%e_dc)
1819 : : end if
1820 [ + + + + ]: 1505 : if(dtset%occopt>=3.and.dtset%occopt<=8.and.ipositron==0) then
1821 [ + + ]: 689 : if(.not.testdmft) then
1822 : : !write(msg, '(a,es21.14,a,a,a,es21.14)' ) &
1823 : : ! ' >>>>> Internal E= ',etotaldc-eent,ch10,ch10,&
1824 : : ! ' -kT*entropy = ',eent
1825 : : !call wrtout(iout,msg)
1826 : 572 : call dc_edoc%add_real('internal', etotaldc-eent)
1827 : 572 : call dc_edoc%add_real('-kT*entropy', eent)
1828 : : else
1829 : 117 : call dc_edoc%add_real('internal', etotaldc-eent)
1830 : : end if
1831 [ + + ]: 810 : else if (ipositron/=0) then
1832 [ + + ]: 12 : if (dtset%occopt>=3 .and. dtset%occopt<=8) then
1833 : 6 : call dc_edoc%add_real('-kT*entropy', eent)
1834 : : end if
1835 : : !write(msg, '(a,es21.14,4a,es21.14,a)' ) &
1836 : : ! ' - EP dble-ct En.= ',-energies%edc_electronpositron,ch10,&
1837 : : ! ' >>> ',EPName(ipositron),' E= ',etotaldc-energies%e0_electronpositron -energies%e_electronpositron,ch10
1838 : : !call wrtout(iout,msg)
1839 : : !write(msg, '(3a,es21.14,2a,es21.14)' ) &
1840 : : ! ' ',EPName(3-ipositron),' ener.= ',energies%e0_electronpositron,ch10,&
1841 : : ! ' EP interaction E= ' ,energies%e_electronpositron
1842 : : !call wrtout(iout,msg)
1843 : 12 : call dc_edoc%add_real('electron_positron_dc', -energies%edc_electronpositron)
1844 [ + + ]: 12 : if(ipositron == 1) then
1845 : 5 : call dc_edoc%add_real('positronic', etotaldc-energies%e0_electronpositron-energies%e_electronpositron)
1846 : 5 : call dc_edoc%add_real('electronic', energies%e0_electronpositron)
1847 : : else
1848 : 7 : call dc_edoc%add_real('electronic', etotaldc-energies%e0_electronpositron-energies%e_electronpositron)
1849 : 7 : call dc_edoc%add_real('positronic', energies%e0_electronpositron)
1850 : : end if
1851 : 12 : call dc_edoc%add_real('electron_positron_interaction', energies%e_electronpositron)
1852 : : end if
1853 : :
1854 : 1505 : write(msg, '(a,es21.14)' ) ' >>>> Etotal (DC)= ',etotaldc
1855 : : !call wrtout(iout,msg)
1856 : 1505 : call dc_edoc%add_real('total_energy_dc', etotaldc)
1857 : : end if
1858 : :
1859 : : !======= Additional printing ==========
1860 : :
1861 [ + + + + ]: 5281 : if (usepaw==0.and.optdc==0) then
1862 : 3774 : call edoc%add_real('total_energy_eV', etotal*Ha_eV)
1863 : 3774 : call edoc%add_real('band_energy', energies%e_eigenvalues)
1864 : : end if
1865 : :
1866 [ + - + + ]: 5281 : if ((optdc==0.or.optdc==2).and.(.not.directE_avail)) then
1867 : : !write(msg, '(a,a,es18.10)' ) ch10,' Band energy (Ha)= ',energies%e_eigenvalues
1868 : : !call wrtout(iout,msg)
1869 : 1 : call edoc%add_real('band_energy', energies%e_eigenvalues)
1870 : : end if
1871 : :
1872 [ + + ]: 5281 : if (usepaw==1) then
1873 [ + - + + ]: 1380 : if ((optdc==0.or.optdc==2).and.(directE_avail)) then
1874 : 1379 : call edoc%add_real('total_energy_eV', etotal*Ha_eV)
1875 : : end if
1876 [ + + ]: 1380 : if (optdc>=1) then
1877 : : !if (optdc==1) write(msg, '(a,a,es21.14)' ) ch10,' >Total DC energy in eV = ',etotaldc*Ha_eV
1878 : : !if (optdc==2) write(msg, '(a,es21.14)' ) ' >Total DC energy in eV = ',etotaldc*Ha_eV
1879 : : !call wrtout(iout,msg)
1880 : 1378 : call dc_edoc%add_real('total_energy_dc_eV', etotaldc*Ha_eV)
1881 : : end if
1882 : : end if
1883 : :
1884 [ + - + + ]: 5281 : if( dtset%icoulomb/=1.and.abs(dtset%cellcharge(1))>tol8) then
1885 : : write(msg, '(6a)' ) &
1886 : 23 : ch10,' Calculation was performed for a charged system with PBC',&
1887 : 23 : ch10,' You may consider including the monopole correction to the total energy',&
1888 : 46 : ch10,' The correction is to be divided by the dielectric constant'
1889 : 23 : call wrtout(iout,msg)
1890 : 23 : call edoc%add_real('monopole_correction', energies%e_monopole)
1891 : 23 : call edoc%add_real('monopole_correction_eV', energies%e_monopole*Ha_eV)
1892 : : end if
1893 : :
1894 : : !Print total energy including PAW core contribution
1895 [ + + ]: 5281 : if (usepaw==1) then
1896 [ + + + + ]: 1380 : if (dtset%paw_add_core==0.and.dtset%use_rcpaw==0) then
1897 [ + + - + ]: 1376 : if (abs(energies%paw%epaw_core)>tiny(zero).or.abs(energies%paw%epaw_core_dc)>tiny(zero)) then
1898 : 498 : write_epaw_core=.true.
1899 : 498 : info = 'Components of total free energy, including PAW core contributions'
1900 : : pawcore_edoc = yamldoc_open('EnergyTermsWithPAWCore', info=trim(adjustl(info)), &
1901 [ + - ]: 498 : width=20, real_fmt='(es21.14)')
1902 : 498 : call pawcore_edoc%add_real('Total energy', etotal+energies%paw%epaw_core)
1903 [ + + ]: 498 : if(optdc>=1) call pawcore_edoc%add_real('Total energy DC', etotaldc+energies%paw%epaw_core)
1904 : 498 : call pawcore_edoc%add_real('Total energy (eV)', (etotal+energies%paw%epaw_core)*Ha_eV)
1905 [ + + ]: 498 : if(optdc>=1) call pawcore_edoc%add_real('Total energy DC (eV)', (etotaldc+energies%paw%epaw_core)*Ha_eV)
1906 : : end if
1907 : : end if
1908 : : end if
1909 : :
1910 : : !======== In case other sources of entropies than the non-interacting entropy =========
1911 : : !============= of the Kohn-Sham states come into play, print the details ==============
1912 [ + + ]: 5281 : if(dtset%occopt>=3.and.dtset%occopt<=8) then
1913 [ + + + + ]: 1282 : if(abs(energies%entropy)>tiny(zero).and.abs(energies%entropy-energies%entropy_ks)>tiny(zero)) then
1914 : 11 : write_entropy=.true.
1915 : : sdoc = yamldoc_open('EntropyTerms', info='Components of total entropy', &
1916 : 11 : & width=20, real_fmt="(es21.14)") ! in kB units
1917 : 11 : call sdoc%add_real('noninteracting',energies%entropy_ks) ! Noninteracting entropy = Entropy of the Kohn-Sham states
1918 [ + + ]: 11 : if(abs(energies%entropy_xc)>tiny(zero)) call sdoc%add_real('xc',energies%entropy_xc)
1919 [ + + + + ]: 11 : if(usepaw==1.and.abs(energies%paw%entropy_paw)>tiny(zero)) call sdoc%add_real('spherical_terms',energies%paw%entropy_paw)
1920 [ + + ]: 11 : if(abs(energies%entropy_extfpmd)>tiny(zero)) call sdoc%add_real('extfpmd',energies%entropy_extfpmd)
1921 [ - + ]: 11 : if(abs(energies%entropy_imp)>tiny(zero)) call sdoc%add_real('impurity',energies%entropy_imp)
1922 : 11 : call sdoc%add_real('total_entropy',energies%entropy) ! Total entropy energy
1923 : : end if
1924 : : end if
1925 : :
1926 : : !======== In case finite-temperature exchange-correlation functionals are used =========
1927 : : !=================== write the total exchange-correlation components ===================
1928 : : !=============================== For testing purposes only =============================
1929 [ + + ]: 5281 : if(abs(energies%entropy_xc)>tiny(zero)) then
1930 : 5 : write_totalxc=.false. ! For testing purposes only, set write_totalxc=.true.
1931 : : if(write_totalxc) then
1932 : : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
1933 : : ftxcdoc = yamldoc_open('FTXCEnergyTerms', info='Components of total xc energy in Hartree', &
1934 : : & width=20, real_fmt="(es21.14)")
1935 : : if(usepaw==1) then
1936 : : ! For now, only finite-temperature xc functionals contribute to entropy_paw.
1937 : : ! We may introduce 'energies%entropy_pawxc' in the future.
1938 : : call ftxcdoc%add_real('xc',energies%e_xc)
1939 : : call ftxcdoc%add_real('spherical_terms_xc',energies%paw%epaw_xc)
1940 : : call ftxcdoc%add_real('internal_xc',energies%e_xc+energies%paw%epaw_xc)
1941 : : call ftxcdoc%add_real('-kT*entropy_xc',-el_temp*(energies%entropy_xc+energies%paw%entropy_paw))
1942 : : call ftxcdoc%add_real('free_xc',energies%e_xc+energies%paw%epaw_xc-el_temp*(energies%entropy_xc+energies%paw%entropy_paw))
1943 : : else
1944 : : call ftxcdoc%add_real('internal_xc',energies%e_xc)
1945 : : call ftxcdoc%add_real('-kT*entropy_xc',-el_temp*energies%entropy_xc)
1946 : : call ftxcdoc%add_real('free_xc',energies%e_xc-el_temp*energies%entropy_xc)
1947 : : end if
1948 : : end if
1949 : : end if
1950 : :
1951 : : ! Write components of total energies in Yaml format.
1952 : 5281 : call edoc%write_and_free(iout)
1953 [ + + ]: 5281 : if(optdc >= 1) call dc_edoc%write_and_free(iout)
1954 [ + + ]: 5281 : if (write_epaw_core) call pawcore_edoc%write_and_free(iout)
1955 [ + + ]: 5281 : if(write_entropy) call sdoc%write_and_free(iout)
1956 [ - + ]: 5281 : if(write_totalxc) call ftxcdoc%write_and_free(iout)
1957 : :
1958 [ + + ]: 5281 : end subroutine prtene
1959 : : !!***
1960 : :
1961 : : !!****f* ABINIT/get_dtsets_pspheads
1962 : : !! NAME
1963 : : !! get_dtsets_pspheads
1964 : : !!
1965 : : !! FUNCTION
1966 : : !! Parse input file, get list of pseudos for files file and build list of datasets
1967 : : !! pseudopotential headers, maxval of dimensions needed in outvars
1968 : : !!
1969 : : !! INPUTS
1970 : : !! input_path: Input filename specifed on the command line. zero lenght if files file syntax is used.
1971 : : !! Mainly used to check whether pseudos are defined in the input to avoid entering the files file
1972 : : !! branch that prompts for pseudos.
1973 : : !! path: Input Filename
1974 : : !! comm: MPI communicator
1975 : : !!
1976 : : !! OUTPUT
1977 : : !! lenstr= the length of the resulting string.
1978 : : !! ndtset= the number of declared datasets.
1979 : : !! string= contains on output the content of the file, ready for parsing.
1980 : : !! dtsets(0:ndtset): List of datasets
1981 : : !! dmatpuflag=flag controlling the use of an initial density matrix in PAW+U (max. value over datasets)
1982 : : !! mx<ab_dimensions>=datatype storing the maximal dimensions.
1983 : : !! pspheads(npsp)=<type pspheader_type>=all the important information from the
1984 : : !! pseudopotential file headers, as well as the psp file names
1985 : : !!
1986 : : !! SOURCE
1987 : :
1988 : 1446 : subroutine get_dtsets_pspheads(input_path, path, ndtset, lenstr, string, timopt, dtsets, pspheads, mx, dmatpuflag, comm)
1989 : :
1990 : : !Arguments ------------------------------------
1991 : : !scalars
1992 : : integer,intent(out) :: lenstr, ndtset
1993 : : type(ab_dimensions),intent(out) :: mx
1994 : : character(len=strlen), intent(out) :: string
1995 : : character(len=*),intent(in) :: input_path, path
1996 : : integer,intent(in) :: comm
1997 : : integer,intent(out) :: timopt, dmatpuflag
1998 : : !arrays
1999 : : type(dataset_type),allocatable,intent(out) :: dtsets(:)
2000 : : type(pspheader_type),allocatable,intent(out):: pspheads(:)
2001 : :
2002 : : !Local variables-------------------------------
2003 : : !scalars
2004 : : integer :: ipsp,ios, me, ndtset_alloc, nprocs
2005 : : integer :: istatr,istatshft, papiopt, npsp, ii, idtset, msym, usepaw
2006 : : character(len=fnlen) :: filpsp
2007 : : character(len=500) :: msg
2008 : : !arrays
2009 : 1446 : integer,allocatable :: mband_upper_(:)
2010 : : real(dp) :: ecut_tmp(3,2,10),tsec(2)
2011 : 1446 : real(dp),allocatable :: zionpsp(:)
2012 : 1446 : character(len=fnlen), allocatable :: pspfilnam_(:), pseudo_paths(:)
2013 : : !************************************************************************
2014 : :
2015 : 1446 : me = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
2016 : :
2017 : : ! Read the file, stringify it and return the number of datasets.
2018 : 1446 : call parsefile(path, lenstr, ndtset, string, comm)
2019 : :
2020 [ + + ]: 1446 : ndtset_alloc = ndtset; if (ndtset == 0) ndtset_alloc=1
2021 [ - + + - : 196407 : ABI_MALLOC(dtsets, (0:ndtset_alloc))
+ - + - -
+ - + - +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + - ]
2022 : :
2023 : 1446 : timopt = 1; if (xmpi_paral==1) timopt = 0
2024 : :
2025 : : ! Continue to analyze the input string, get upper dimensions, and allocate the remaining arrays.
2026 : : call invars0(dtsets, istatr, istatshft, lenstr, msym, mx%natom, mx%nimage, mx%ntypat, &
2027 : 1446 : ndtset, ndtset_alloc, npsp, pseudo_paths, papiopt, timopt, string, comm)
2028 : :
2029 : : ! Enable PAPI timers
2030 : 1446 : call time_set_papiopt(papiopt)
2031 : :
2032 [ + + ]: 9873 : dtsets(:)%timopt = timopt
2033 : : dtsets(0)%timopt = 1
2034 : 1446 : if (xmpi_paral == 1) dtsets(0)%timopt = 0
2035 : :
2036 : 1446 : call timab(timopt,5,tsec)
2037 : :
2038 : : ! Initialize pspheads, that contains the important information
2039 : : ! from the pseudopotential headers, as well as the psp filename
2040 : 1446 : call timab(102,1,tsec)
2041 : 1446 : call timab(1021,3,tsec)
2042 : :
2043 : 1446 : usepaw = 0
2044 [ - + + - : 6294 : ABI_MALLOC(pspheads, (npsp))
+ - + - -
+ - + - +
+ + + - ]
2045 [ - + ]: 1446 : if (npsp > 10) then
2046 : 0 : ABI_BUG('ecut_tmp is not well defined.')
2047 : : end if
2048 [ + + + + : 131586 : ecut_tmp = -one
+ + ]
2049 : :
2050 [ + + ]: 3402 : pspheads(:)%usewvl = dtsets(1)%usewvl
2051 : :
2052 [ + + ]: 1446 : if (me == 0) then
2053 [ - + + - : 3495 : ABI_MALLOC(pspfilnam_, (npsp))
+ - + - -
+ - + ]
2054 : :
2055 [ + + ]: 1165 : if (len_trim(pseudo_paths(1)) == 0) then
2056 : : ! Enter Legacy `files file` mode --> Read the name of the psp file from files file.
2057 : :
2058 : : ! Catch possible mistake done by user (input without pseudos and `abinit t01.in` syntax)
2059 : : ! else the code starts to prompt for pseudos and execution gets stuck
2060 [ - + ]: 2 : if (len_trim(input_path) /= 0) then
2061 : 0 : ABI_ERROR("`pseudos` variable must be specified in input when the code is invoked with the `abinit t01.in` syntax")
2062 : : end if
2063 : :
2064 : : ! Finish to read the "file" file completely, as npsp is known,
2065 [ + + ]: 6 : do ipsp=1,npsp
2066 : 4 : write(std_out,'(/,a)' )' Please give name of formatted atomic psp file (and finish with a newline character)'
2067 : 4 : read (std_in, '(a)' , iostat=ios ) filpsp
2068 : : ! It might be that a file name is missing
2069 [ - + ]: 4 : if (ios /= 0) then
2070 : : write(msg, '(5a)' )&
2071 : 0 : 'There are not enough names of pseudopotentials provided in the files file.',ch10,&
2072 : 0 : 'Action: check first the variable ntypat (and/or npsp) in the input file;',ch10,&
2073 : 0 : 'if they are correct, complete your files file.'
2074 : 0 : ABI_ERROR(msg)
2075 : : end if
2076 : 4 : pspfilnam_(ipsp) = trim(filpsp)
2077 : 6 : write(std_out,'(a,i0,2a)' )' For atom type ',ipsp,', psp file is ',trim(filpsp)
2078 : : end do ! ipsp
2079 : :
2080 : : else
2081 : : ! Get pseudopotential paths from input file.
2082 [ + - - + : 3895 : pspfilnam_ = pseudo_paths
- - - - +
+ ]
2083 [ + + ]: 2732 : do ipsp=1,npsp
2084 [ + - ]: 2732 : write(std_out,'(a,i0,2a)' )' For atom type ',ipsp,', psp file is ',trim(pspfilnam_(ipsp))
2085 : : end do
2086 : : end if
2087 : :
2088 : : ! Now read the psp headers
2089 : 1165 : call inpspheads(pspfilnam_, npsp, pspheads, ecut_tmp)
2090 [ - + ]: 1165 : ABI_FREE(pspfilnam_)
2091 : :
2092 [ + + + + ]: 2738 : if (minval(abs(pspheads(1:npsp)%pspcod - 7)) == 0) usepaw=1
2093 [ + + + + ]: 2738 : if (minval(abs(pspheads(1:npsp)%pspcod - 17)) == 0) usepaw=1
2094 : : end if ! me == 0
2095 : :
2096 [ - + ]: 1446 : ABI_FREE(pseudo_paths)
2097 : :
2098 : : ! Communicate pspheads to all processors
2099 : 1446 : call pspheads_comm(npsp, pspheads, usepaw)
2100 : :
2101 : : ! If (all) pspcod are 7 then this is a PAW calculation. Initialize (default) the value of ratsph
2102 [ + + ]: 9873 : do idtset=0,ndtset_alloc
2103 : 8427 : dtsets(idtset)%usepaw = usepaw
2104 [ + + ]: 9873 : if (usepaw == 0) then
2105 [ + + ]: 15046 : dtsets(idtset)%ratsph(:)=two
2106 : : else
2107 : : ! Note that the following coding assumes that npsp=ntypat for PAW, which is true as of now (XG20101024).
2108 : : ! dtsets(idtset)%ratsph(1:npsp)=token%pspheads(1:npsp)%pawheader%rpaw
2109 [ + + ]: 3937 : do ipsp=1,npsp
2110 : 3937 : dtsets(idtset)%ratsph(ipsp) = pspheads(ipsp)%pawheader%rpaw
2111 : : end do
2112 : : endif
2113 : : end do
2114 : :
2115 : : ! Take care of other dimensions, and part of the content of dtsets that is or might be needed early.
2116 [ - + + - : 4338 : ABI_MALLOC(zionpsp, (npsp))
+ - + - -
+ - + ]
2117 [ + + ]: 3402 : do ii=1,npsp
2118 : 3402 : zionpsp(ii) = pspheads(ii)%zionpsp
2119 : : end do
2120 : :
2121 [ - + + - : 4338 : ABI_MALLOC(mband_upper_, (0:ndtset_alloc))
+ - + - -
+ - + ]
2122 : :
2123 : : ! Get MAX dimension over datasets
2124 : : call invars1m(dmatpuflag, dtsets, ab_out, lenstr, mband_upper_, mx,&
2125 : 1446 : msym, ndtset, ndtset_alloc, string, npsp, zionpsp, comm)
2126 : :
2127 [ - + ]: 1432 : ABI_FREE(zionpsp)
2128 : 1432 : call timab(1021,2,tsec)
2129 : 1432 : call timab(1022,3,tsec)
2130 : :
2131 : : ! Provide defaults for the variables that have not yet been initialized.
2132 : 1432 : call indefo(dtsets, ndtset_alloc, nprocs)
2133 : :
2134 : 1432 : call timab(1022,2,tsec)
2135 : 1432 : call timab(1023,3,tsec)
2136 : :
2137 : : ! Perform some global initialization, depending on the value of
2138 : : ! pseudopotentials, parallelism variables, or macro input variables
2139 : 1432 : call macroin(dtsets, ecut_tmp, lenstr, ndtset_alloc, string)
2140 : :
2141 : : ! If all the pseudopotentials have the same pspxc, override the default value for dtsets 1 to ndtset
2142 [ + + + - ]: 3374 : if (minval(abs((pspheads(1:npsp)%pspxc - pspheads(1)%pspxc)))==0) then
2143 [ + + ]: 8399 : dtsets(1:ndtset_alloc)%ixc = pspheads(1)%pspxc
2144 : : end if
2145 : :
2146 : : ! Call the main input routine.
2147 : 1432 : call invars2m(dtsets,ab_out,lenstr,mband_upper_,msym,ndtset,ndtset_alloc,npsp,pspheads,string, comm)
2148 : :
2149 : 1432 : call macroin2(dtsets, ndtset_alloc)
2150 : :
2151 : 1432 : mx%mband = dtsets(1)%mband
2152 [ + + ]: 8399 : do ii=1,ndtset_alloc
2153 : 8399 : mx%mband = max(dtsets(ii)%mband, mx%mband)
2154 : : end do
2155 : :
2156 : 1432 : call timab(1023,2,tsec)
2157 : 1432 : call timab(102,2,tsec)
2158 : :
2159 [ - + ]: 1432 : ABI_FREE(mband_upper_)
2160 : :
2161 : 2892 : end subroutine get_dtsets_pspheads
2162 : : !!***
2163 : :
2164 : : !!****f* ABINIT/ebands_from_file
2165 : : !! NAME
2166 : : !! ebands_from_file
2167 : : !!
2168 : : !! FUNCTION
2169 : : !! Build and ebands_t object from file. Supports Fortran and netcdf files
2170 : : !! provided they have a Abinit header and obviously GS eigenvalues
2171 : : !!
2172 : : !! INPUTS
2173 : : !! path: File name.
2174 : : !! comm: MPI communicator.
2175 : : !!
2176 : : !! OUTPUT
2177 : : !!
2178 : : !! SOURCE
2179 : :
2180 : :
2181 : 0 : type(ebands_t) function ebands_from_file(path, comm) result(new)
2182 : :
2183 : : !Arguments ------------------------------------
2184 : : !scalars
2185 : : character(len=*),intent(in) :: path
2186 : : integer,intent(in) :: comm
2187 : :
2188 : : !Local variables-------------------------------
2189 : : !scalars
2190 : : integer :: ncid, fform
2191 : 0 : type(hdr_type) :: hdr
2192 : : !arrays
2193 : 0 : real(dp),pointer :: gs_eigen(:,:,:)
2194 : :
2195 : : ! *************************************************************************
2196 : :
2197 : : ! NOTE: Assume file with header. Must use wfk_read_eigenvalues to handle Fortran WFK
2198 [ # # ]: 0 : if (endswith(path, "_WFK") .or. endswith(path, "_WFK.nc")) then
2199 : 0 : call wfk_read_eigenvalues(path, gs_eigen, hdr, comm)
2200 [ # # # # ]: 0 : call new%from_hdr(hdr, maxval(hdr%nband), gs_eigen)
2201 : :
2202 [ # # ]: 0 : else if (endswith(path, ".nc")) then
2203 [ # # ]: 0 : NCF_CHECK(nctk_open_read(ncid, path, comm))
2204 : 0 : call hdr%ncread(ncid, fform)
2205 [ # # ]: 0 : ABI_CHECK(fform /= 0, "fform == 0")
2206 [ # # # # : 0 : ABI_MALLOC(gs_eigen, (hdr%mband, hdr%nkpt, hdr%nsppol))
# # # # #
# # # # #
# # # # #
# ]
2207 [ # # ]: 0 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "eigenvalues"), gs_eigen))
2208 [ # # # # ]: 0 : call new%from_hdr(hdr, maxval(hdr%nband), gs_eigen)
2209 [ # # ]: 0 : NCF_CHECK(nf90_close(ncid))
2210 : : else
2211 : 0 : ABI_ERROR(sjoin("Don't know how to construct crystal structure from: ", path, ch10, "Supported extensions: _WFK or .nc"))
2212 : : end if
2213 : :
2214 [ # # ]: 0 : ABI_FREE(gs_eigen)
2215 : 0 : call hdr%free()
2216 : :
2217 [ # # # # : 0 : end function ebands_from_file
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
2218 : : !!***
2219 : :
2220 : : !!****f* ABINIT/crystal_from_file
2221 : : !! NAME
2222 : : !! crystal_from_file
2223 : : !!
2224 : : !! FUNCTION
2225 : : !! Build crystal_t object from netcdf file
2226 : : !!
2227 : : !! INPUTS
2228 : : !!
2229 : : !! OUTPUT
2230 : : !!
2231 : : !! SOURCE
2232 : :
2233 : 0 : type(crystal_t) function crystal_from_file(path, comm) result(new)
2234 : :
2235 : : !Arguments ------------------------------------
2236 : : !scalars
2237 : : character(len=*),intent(in) :: path
2238 : : integer,intent(in) :: comm
2239 : :
2240 : : !Local variables-------------------------------
2241 : : !scalars
2242 : : integer :: fform, ncid
2243 : 0 : type(hdr_type) :: hdr
2244 : : ! *************************************************************************
2245 : :
2246 : 0 : if (endswith(path, ".nc")) then
2247 [ # # ]: 0 : NCF_CHECK(nctk_open_read(ncid, path, comm))
2248 : 0 : call new%ncread(ncid)
2249 [ # # ]: 0 : NCF_CHECK(nf90_close(ncid))
2250 : : else
2251 : : ! Assume file with Abinit header
2252 : : ! TODO: Should add routine to read crystal from structure without hdr
2253 : 0 : call hdr%from_fname(path, fform, comm)
2254 [ # # ]: 0 : ABI_CHECK(fform /= 0, "fform == 0")
2255 : 0 : new = hdr%get_crystal()
2256 : 0 : call hdr%free()
2257 : : end if
2258 : :
2259 [ # # # # : 0 : end function crystal_from_file
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
2260 : : !!***
2261 : :
2262 : :
2263 : : !!****f* ABINIT/get_gemm_nonlop_ompgpu_blocksize
2264 : : !! NAME
2265 : : !! get_gemm_nonlop_ompgpu_blocksize
2266 : : !!
2267 : : !! FUNCTION
2268 : : !! Check GPU memory constraints for many ABINIT routines (lobpcg,chebfi,getghc)
2269 : : !! and compute a size of GEMM nonlop block that is small enough to fit in GPU memory.
2270 : : !! Print its estimation of memory consumption in output.
2271 : : !!
2272 : : !! INPUTS
2273 : : !! gs_hamk : hamiltonian structure (used by getghc)
2274 : : !! ndat : size of batching, usually matching bandpp or nblock_lobpcg
2275 : : !! npw : number of planewaves
2276 : : !! nband : number of bands
2277 : : !! nspinor : number of spin-orbits
2278 : : !! paral_kgb : 1 if KGB parallelism is enabled, 0 otherwise
2279 : : !! optfor : 1 if forces computation is enabled, 0 otherwise
2280 : : !! optstr : 1 if stress computation is enabled, 0 otherwise
2281 : : !! wfoptalg : Which diago algorithm if used:
2282 : : !! -1: none (GEMM nonlop isn't used for diago)
2283 : : !! 111: CHEBFI2
2284 : : !! 112: SLICE
2285 : : !! 114: LOBPCG2
2286 : : !! other: Only account for getghc
2287 : : !! gpu_option : If GPU is enabled (expected to be ABI_GPU_OPENMP for now)
2288 : : !! blocksize : if higher than 0, only print memory estimation and exit
2289 : : !!
2290 : : !! OUTPUT
2291 : : !! nfft_blocks : Number of blocks fourwf computation will be divided into
2292 : : !! blocksize : Size of MPI tasks blocks to be used in GEMM nonlop
2293 : : !! nblocks : Number of MPI blocks to be used in GEMM nonlop
2294 : : !!
2295 : : !! SOURCE
2296 : 0 : subroutine get_gemm_nonlop_ompgpu_blocksize(ikpt,gs_hamk,ndat,nband,nspinor,nspden,paral_kgb,&
2297 : : & npband,optfor,optstr,wfoptalg,gpu_option,use_distrib,&
2298 : : & blocksize,nblocks,nfft_blocks,warn_on_fail,disable_output)
2299 : :
2300 : : integer,intent(in) :: ikpt,ndat,nband,nspinor,nspden,paral_kgb,npband,optfor,optstr,wfoptalg,gpu_option
2301 : : logical,intent(in) :: use_distrib
2302 : : logical,intent(in),optional :: warn_on_fail,disable_output
2303 : : type(gs_hamiltonian_type),intent(in) :: gs_hamk
2304 : : integer,intent(inout) :: blocksize,nfft_blocks
2305 : : integer,intent(out) :: nblocks
2306 : :
2307 : : integer(kind=c_size_t) :: nonlop_smem,invovl_smem,getghc_wmem,invovl_wmem,nonlop_wmem,gs_ham_smem,updrho_wmem,prep_nonlop_wmem
2308 : : integer(kind=c_size_t) :: sum_mem,sum_bandpp_mem,sum_other_mem,free_mem,localMem,fourwf_smem,fourwf_wmem,fourwf_mem,hegvd_mem
2309 : : integer :: icplx,space,i,ndat_try,rank,nprocs,ndgxdt,blockdim,max_slices,npw,npw_fft,signs,nprojs,itypat
2310 : : integer, target :: t_fft(3)
2311 : : logical :: print_and_exit,l_warn_on_fail,l_disable_output,fixed_blocksize,fixed_fourwf_blocks
2312 : : integer(kind=c_size_t) :: chebfiMem(2),lobpcgMem(2)
2313 : : character(len=500) :: message
2314 : :
2315 : : ! *********************************************************************
2316 : :
2317 : 0 : free_mem=256*1e9 ! Dummy value
2318 [ # # ]: 0 : l_warn_on_fail=.false.;if(present(warn_on_fail)) l_warn_on_fail=warn_on_fail
2319 [ # # ]: 0 : l_disable_output=.false.;if(present(disable_output)) l_disable_output=disable_output
2320 : : #ifdef HAVE_GPU
2321 : : if(gpu_option /= ABI_GPU_DISABLED) then
2322 : : call gpu_get_max_mem(free_mem)
2323 : : ! NOTE: Cutting 10% out to be safe
2324 : : ! I computed this value using the memory allocation summary that
2325 : : ! NVHPC provides after a GPU out-of-memory crash, on a run NVIDIA A100 with 80GB.
2326 : : ! 8.5% is the amount of memory that wasn't reported allocated nor freed.
2327 : : ! I add extra 1.5% to account for minor buffers eventually allocated.
2328 : : ! This amount may either be hidden allocations for CUDA/cuBLAS/cuFFT/cuSOLVER,
2329 : : ! or memory allocated outside ABINIT.
2330 : : free_mem = 0.9 * free_mem
2331 : : end if
2332 : : #else
2333 : : ABI_UNUSED(gpu_option)
2334 : : #endif
2335 : :
2336 [ # # ]: 0 : if(gpu_option /= ABI_GPU_OPENMP) then
2337 : : ! No distribution is attempted outside of OpenMP GPU. User is already warned in chkinp
2338 : 0 : blocksize=1; nblocks=0; nfft_blocks=1
2339 : 0 : return
2340 : : end if
2341 : :
2342 : 0 : rank = xmpi_comm_rank(xmpi_world); nprocs = xmpi_comm_size(xmpi_world)
2343 [ # # ]: 0 : if ( gs_hamk%istwf_k == 2 ) then ! Real only
2344 : 0 : space = SPACE_CR
2345 : 0 : icplx = 2
2346 : : else ! complex
2347 : 0 : space = SPACE_C
2348 : 0 : icplx = 1
2349 : : end if
2350 : :
2351 : 0 : npw=gs_hamk%npw_k
2352 : 0 : npw_fft=gs_hamk%npw_fft_k
2353 : 0 : ndat_try=ndat
2354 : 0 : blockdim=npband*ndat
2355 : 0 : ndgxdt=0
2356 [ # # ]: 0 : if(optfor>0) ndgxdt=ndgxdt+3
2357 [ # # ]: 0 : if(optstr>0) ndgxdt=ndgxdt+6
2358 : 0 : signs=2
2359 : : !wfoptalg==-1 means we're in forstr
2360 [ # # ]: 0 : if(wfoptalg==-1) signs=1
2361 : 0 : t_fft(1) = gs_hamk%ngfft(3);
2362 : 0 : t_fft(2) = gs_hamk%ngfft(2);
2363 : 0 : t_fft(3) = gs_hamk%ngfft(1);
2364 : 0 : nprojs=0
2365 [ # # ]: 0 : do itypat=1,gs_hamk%ntypat
2366 [ # # # # ]: 0 : nprojs = nprojs + count(gs_hamk%indlmn(3,:,itypat)>0)*gs_hamk%nattyp(itypat)
2367 : : end do
2368 : :
2369 : 0 : nonlop_smem=0; invovl_smem=0; getghc_wmem=0; invovl_wmem=0; nonlop_wmem=0; gs_ham_smem=0
2370 : 0 : updrho_wmem=0; prep_nonlop_wmem=0; sum_mem=0; sum_bandpp_mem=0; sum_other_mem=0;
2371 : 0 : localMem=0; fourwf_smem=0; fourwf_wmem=0; fourwf_mem=0; hegvd_mem=0
2372 : 0 : chebfiMem(:)=0; lobpcgMem(:)=0
2373 : :
2374 : : !HEGVD work memory estimate.
2375 : : ! Since *_bufferSize routines from (cu/hip)SOLVER require buffer
2376 : : ! to be provided, I measeured the work size given by those routines
2377 : : ! on many big cases and guess an approximate value.
2378 : : ! hipSolver is eager than cuSolver, hence the extra multiplier
2379 : :
2380 [ # # ]: 0 : if(wfoptalg==111 .or. wfoptalg==11) then
2381 : 0 : hegvd_mem = int(dp, c_size_t) * nband * nband * 3
2382 [ # # ]: 0 : else if(wfoptalg==114 .or. wfoptalg==14) then
2383 : 0 : hegvd_mem = int(dp, c_size_t) * (ndat*3) * (ndat*3) * 3
2384 : : end if
2385 : : #ifdef HAVE_GPU_CUDA
2386 : : if(space == SPACE_C) hegvd_mem = hegvd_mem * 2
2387 : : #endif
2388 : : #ifdef HAVE_GPU_HIP
2389 : : if(space == SPACE_C) hegvd_mem = hegvd_mem * 3
2390 : : ! ROCm 7 memory usage was measured to be more or less
2391 : : ! on par with CUDA but ROCm 6 was indeed ~9 times higher.
2392 : : ! For now, we can't drop ROCm 6 so we keep this workaround.
2393 : : if(gpu_get_lib_version_major() < 7) then
2394 : : hegvd_mem = hegvd_mem * 9
2395 : : end if
2396 : : #endif
2397 : :
2398 : : if(wfoptalg>=0) then
2399 : : #ifdef HAVE_GPU
2400 : : call gpu_fft_get_estimate_work_size(3, c_loc(t_fft), FFT_Z2Z, ndat, fourwf_smem);
2401 : : #endif
2402 : : end if
2403 : :
2404 : 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)
2405 : 0 : getghc_wmem = getghc_ompgpu_work_mem(gs_hamk, ndat, max(nfft_blocks,1))
2406 : : fourwf_wmem = int(2, c_size_t) * dp * gs_hamk%n4 * gs_hamk%n5 * gs_hamk%n6 &
2407 : 0 : & * (ndat/max(nfft_blocks,1) + modulo(ndat,max(nfft_blocks,1)))
2408 : 0 : fourwf_mem = fourwf_wmem+fourwf_smem
2409 : :
2410 : : nonlop_wmem = gemm_nonlop_ompgpu_work_mem(gs_hamk%istwf_k, ndat, ndgxdt, npw_fft,&
2411 : 0 : & gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%ntypat, gs_hamk%lmnmax, signs, wfoptalg)
2412 : 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
2413 [ # # ]: 0 : if(associated(gs_hamk%ph3d_k)) gs_ham_smem = gs_ham_smem + int(2,c_size_t) * npw_fft * gs_hamk%matblk
2414 : 0 : gs_ham_smem = gs_ham_smem*dp
2415 : :
2416 : : ! Mapped arrays used in mkrho or vtowfk
2417 : 0 : updrho_wmem = int(2,c_size_t)*npw_fft*ndat*nspinor ! cwavef
2418 [ # # ]: 0 : if(nspden==4) updrho_wmem = updrho_wmem + int(2,c_size_t)*npw_fft*ndat*2 ! cwavef_x + cwavef_y
2419 : 0 : updrho_wmem = updrho_wmem*dp
2420 : :
2421 : : ! Mapped arrays used in prep_nonlop (only paral_kgb==1)
2422 : 0 : prep_nonlop_wmem = 0
2423 [ # # ]: 0 : if(paral_kgb==1) prep_nonlop_wmem = int(2,c_size_t)*npw_fft*ndat*nspinor*3*dp ! cwavef_alltoall{1,2} + cwavef
2424 : :
2425 [ # # ]: 0 : if(wfoptalg==111 .or. wfoptalg==112) then
2426 : 0 : chebfiMem = chebfi_memInfo(nband,icplx*npw*nspinor,space,paral_kgb,icplx*npw*nspinor,blockdim)
2427 : 0 : invovl_smem = invovl_ompgpu_static_mem(gs_hamk)
2428 : 0 : invovl_wmem = invovl_ompgpu_work_mem(gs_hamk, ndat)
2429 : : end if
2430 [ # # ]: 0 : if(wfoptalg==114) then
2431 : 0 : lobpcgMem = lobpcg_memInfo(nband,icplx*npw*nspinor,space,paral_kgb,blockdim)
2432 : : end if
2433 : 0 : localMem = (int(2,c_size_t)*npw*nspinor*nband+3*nband)*kind(1.d0) ! cg, eig, occ, resid in chebfiwf/lobpcgwf
2434 : :
2435 : : ! Check if arrays outside of GEMM nonlop projectors and ompgpu_fourwf fit in GPU memory
2436 : 0 : sum_other_mem = gs_ham_smem
2437 : :
2438 [ # # ]: 0 : if(wfoptalg>=0) then
2439 : 0 : sum_other_mem = sum_other_mem+updrho_wmem+prep_nonlop_wmem
2440 : : else
2441 : 0 : sum_other_mem = sum_other_mem+nonlop_wmem+prep_nonlop_wmem
2442 : : end if
2443 : :
2444 [ # # ]: 0 : if(wfoptalg==111) then
2445 : 0 : sum_other_mem = sum_other_mem + invovl_wmem+invovl_smem+chebfiMem(1)+chebfiMem(2)+localMem
2446 : : end if
2447 : :
2448 [ # # ]: 0 : if(wfoptalg==114) then
2449 : 0 : sum_other_mem = sum_other_mem + lobpcgMem(1)+lobpcgMem(2)+localMem
2450 : : end if
2451 : :
2452 : 0 : print_and_exit=.false.
2453 : 0 : fixed_blocksize=.false.
2454 : 0 : fixed_fourwf_blocks=.false.
2455 : 0 : nblocks=0
2456 [ # # # # ]: 0 : if(blocksize > 0 .and. nfft_blocks > 0) then
2457 : 0 : nblocks=max(1,nprocs/blocksize)
2458 : 0 : print_and_exit=.true.
2459 : : else
2460 [ # # ]: 0 : if(blocksize > 0) then
2461 : 0 : nblocks=max(1,nprocs/blocksize)
2462 : 0 : fixed_blocksize=.true.
2463 : : else
2464 : 0 : blocksize=1
2465 [ # # ]: 0 : if(.not. l_disable_output) write(std_out,*) "Setting GEMM nonlop block number...", new_line('A')
2466 : : end if
2467 : :
2468 [ # # ]: 0 : if(nfft_blocks > 0) then
2469 : : fixed_fourwf_blocks=.true.
2470 : : else
2471 : 0 : nfft_blocks=1
2472 [ # # ]: 0 : if(.not. l_disable_output) write(std_out,*) "Setting FFT blocks number...", new_line('A')
2473 : : end if
2474 : : end if
2475 : :
2476 [ # # ]: 0 : max_slices=max(100,nprocs*2); if(sum_other_mem > free_mem) max_slices=1
2477 : : ! How we try to optimize GPU memory consumption:
2478 : : ! We work on two variables :
2479 : : ! - blocksize : for slicing GEMM nonlop projectors arrays
2480 : : ! - nfft_blocks : for slicing ompgpu_fourwf work buffer
2481 : : ! At each iteration, we check which routine from GEMM nonlop or ompgpu_fourwf
2482 : : ! have the higher memory requirement and increase related variable.
2483 : : !
2484 : : ! We arbitraly give at least 100 tries to improve GPU memory consumption before eventually
2485 : : ! ruling out that use case is too big to run with available GPU memory and abort.
2486 : : !
2487 : : ! User may hard set slicing for both fourwf and GEMM nonlop, in which case the code will
2488 : : ! warn the user about possible GPU memory overpassing instead of aborting.
2489 : : !
2490 : : ! However, if arrays from other parts of the code already have higher memory requirements,
2491 : : ! we fail anyway and advise the user to increase nblock_lobpcg or run on more nodes.
2492 [ # # ]: 0 : do i=1,max_slices
2493 : :
2494 : : ! First iteration or user provided parameters to split fourwf and GEMM nonlop
2495 : : ! Just measure
2496 [ # # ]: 0 : if(i==1 .or. print_and_exit) then
2497 [ # # ]: 0 : if(wfoptalg>=0) then
2498 : : #ifdef HAVE_GPU
2499 : : call gpu_fft_get_estimate_work_size(3, c_loc(t_fft), FFT_Z2Z, ndat/nfft_blocks, fourwf_smem);
2500 : : #endif
2501 : 0 : getghc_wmem = getghc_ompgpu_work_mem(gs_hamk, ndat, nfft_blocks)
2502 : : fourwf_wmem = int(2, c_size_t) * dp * gs_hamk%n4 * gs_hamk%n5 * gs_hamk%n6 &
2503 : 0 : & * (ndat/nfft_blocks + modulo(ndat,nfft_blocks))
2504 : 0 : fourwf_mem = fourwf_wmem + fourwf_smem
2505 : : end if
2506 : : nonlop_smem = gemm_nonlop_ompgpu_static_mem(npw_fft,gs_hamk%indlmn,gs_hamk%nattyp,&
2507 : 0 : & gs_hamk%ntypat,blocksize,ndgxdt,use_distrib)
2508 : : else
2509 : : ! Raise fourwf slicing if :
2510 : : ! - GEMM nonlop block has been set by user
2511 : : ! or
2512 : : ! - fourwf memory requirements are higher
2513 : : ! - fourwf slicing wasn't set by user
2514 : : ! - fourwf is still sliceable
2515 : : !
2516 : : ! Raise GEMM nonlop blocks otherwise
2517 [ # # # # : 0 : if(fixed_blocksize .or. &
# # # # #
# # # ]
2518 : : & (wfoptalg >= 0 &
2519 : : & .and. nonlop_smem < fourwf_mem &
2520 : : & .and. fourwf_mem >= getghc_wmem &
2521 : : & .and. .not. fixed_fourwf_blocks &
2522 : : & .and. nfft_blocks < ndat &
2523 : : & .and. ndat_try > 1)) then
2524 : : ! Fourwf work memory requirement is higher, split here
2525 [ # # ]: 0 : if(nfft_blocks == ndat) cycle ! Can't split more than ndat
2526 : :
2527 : : if(i>1 .and. .not. print_and_exit) then
2528 [ # # ]: 0 : do while(ndat_try <= (ndat/nfft_blocks + modulo(ndat,nfft_blocks)))
2529 : 0 : nfft_blocks=nfft_blocks+1
2530 : : end do
2531 : : ndat_try = (ndat/nfft_blocks + modulo(ndat,nfft_blocks))
2532 : : end if
2533 : : #ifdef HAVE_GPU
2534 : : call gpu_fft_get_estimate_work_size(3, c_loc(t_fft), FFT_Z2Z, ndat/nfft_blocks, fourwf_smem);
2535 : : #endif
2536 : 0 : getghc_wmem = getghc_ompgpu_work_mem(gs_hamk, ndat, nfft_blocks)
2537 : : fourwf_wmem = int(2, c_size_t) * dp * gs_hamk%n4 * gs_hamk%n5 * gs_hamk%n6 &
2538 : 0 : & * (ndat/nfft_blocks + modulo(ndat,nfft_blocks))
2539 : 0 : fourwf_mem = fourwf_wmem + fourwf_smem
2540 : : else
2541 : : ! Gemm nonlop static memory requirement is higher, split here
2542 : 0 : if(i>1 .and. .not. print_and_exit) blocksize = blocksize + 1
2543 [ # # # # ]: 0 : if(modulo(nprocs,blocksize)/=0 .and. use_distrib) cycle
2544 [ # # # # ]: 0 : if(nprocs < blocksize .and. use_distrib) cycle
2545 : : !FIXME : Skipping uneven blocksize <=5 if using MPI distrib, as the amount of GPU per node is even usually
2546 : : !For example, with 3 nodes of 4 GPU, we don't want to have a blocksize of 3 as
2547 : : !it would generate 4 comms-block, with 2 inter-node comms.
2548 : : !While using a blocksize of 4 would generate 3 comms, one for each node, leading to less MPI comms
2549 [ # # # # ]: 0 : if(i>1 .and. modulo(blocksize,2)/=0 .and. use_distrib .and. .not. print_and_exit) cycle
2550 : 0 : if(i>1) nblocks=nprocs/blocksize
2551 : :
2552 : : nonlop_smem = gemm_nonlop_ompgpu_static_mem(npw_fft,gs_hamk%indlmn,gs_hamk%nattyp,&
2553 : 0 : & gs_hamk%ntypat,blocksize,ndgxdt,use_distrib)
2554 : : end if
2555 : : end if
2556 : :
2557 : : ! Bandpp~ndat sized buffer memory requirements are higher, split there
2558 : 0 : sum_mem = nonlop_smem + gs_ham_smem
2559 : 0 : sum_bandpp_mem = getghc_wmem
2560 : :
2561 [ # # ]: 0 : if(wfoptalg>=0) then
2562 : 0 : sum_mem = sum_mem+getghc_wmem+updrho_wmem+prep_nonlop_wmem+fourwf_smem+hegvd_mem
2563 : : else
2564 : 0 : sum_mem = sum_mem+nonlop_wmem+prep_nonlop_wmem
2565 : : end if
2566 : :
2567 [ # # ]: 0 : if(wfoptalg==111 .or. wfoptalg==112) then
2568 : 0 : sum_mem = sum_mem + invovl_smem+invovl_wmem+chebfiMem(1)+chebfiMem(2)+localMem
2569 : 0 : sum_bandpp_mem = sum_bandpp_mem + invovl_wmem
2570 : : end if
2571 : :
2572 [ # # ]: 0 : if(wfoptalg==114) then
2573 : 0 : sum_mem = sum_mem + lobpcgMem(1)+lobpcgMem(2)+localMem
2574 : : end if
2575 : :
2576 [ # # ]: 0 : if(sum_mem < free_mem .or. print_and_exit) exit
2577 : :
2578 : : end do
2579 : :
2580 : : ! Corner case : not enough GPU memory in forstrnps for forces and stress computation.
2581 : : ! By default, forces and stress are computed in one gemm_nonlop call using choice==23
2582 : : ! This translates to have various arrays sized by ndgxdt == 9 (6 for stress, 3 for forces)
2583 : : ! To try circumventing the lack of GPU memory in that case, we may compute stress and forces
2584 : : ! separately so arrays will be sized after ndgxdt=6 at most instead.
2585 [ # # # # : 0 : if(sum_mem > free_mem .and. optfor > 0 .and. optstr > 0 .and. wfoptalg < 0) then
# # # # ]
2586 : 0 : ndgxdt = 6 ! number of derivatives for stress
2587 : : nonlop_wmem = gemm_nonlop_ompgpu_work_mem(gs_hamk%istwf_k, ndat, ndgxdt, npw_fft,&
2588 : 0 : & gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%ntypat, gs_hamk%lmnmax, signs, wfoptalg)
2589 : 0 : blocksize=1
2590 : : ! Same loop as above, simplified to forstrnps use case
2591 [ # # ]: 0 : do i=1,nprocs
2592 : : ! Gemm nonlop static memory requirement is higher, split here
2593 [ # # ]: 0 : if(i>1 .and. .not. print_and_exit) blocksize = blocksize + 1
2594 [ # # # # ]: 0 : if(modulo(nprocs,blocksize)/=0 .and. use_distrib) cycle
2595 : : !FIXME : Skipping uneven blocksize <=5 if using MPI distrib, as the amount of GPU per node is even usually
2596 : : !For example, with 3 nodes of 4 GPU, we don't want to have a blocksize of 3 as
2597 : : !it would generate 4 comms-block, with 2 inter-node comms.
2598 : : !While using a blocksize of 4 would generate 3 comms, one for each node, leading to less MPI comms
2599 [ # # # # : 0 : if(i>1 .and. modulo(blocksize,2)/=0 .and. use_distrib .and. .not. print_and_exit) cycle
# # # # ]
2600 : 0 : if(i>1) nblocks=nprocs/blocksize
2601 : :
2602 : : nonlop_smem = gemm_nonlop_ompgpu_static_mem(npw_fft,gs_hamk%indlmn,gs_hamk%nattyp,gs_hamk%ntypat,&
2603 : 0 : & blocksize,ndgxdt, use_distrib)
2604 : 0 : sum_mem = nonlop_smem + gs_ham_smem + nonlop_wmem + prep_nonlop_wmem
2605 : :
2606 [ # # ]: 0 : if(sum_mem < free_mem) then
2607 : 0 : gemm_nonlop_split_choice23 = .true.
2608 : 0 : exit
2609 : : end if
2610 : : end do
2611 : : end if
2612 : :
2613 : : ! Quickfix : sometimes, we may run out of GPU memory when computing forces/stresses because of fragmentation.
2614 : : ! We try to reduce the risk by forcing even more blocking:
2615 [ # # # # : 0 : if((wfoptalg < 0 .and. (optfor > 0 .or. optstr > 0)) .and. sum_mem > 0.95*free_mem) then
# # # # ]
2616 [ # # ]: 0 : if(blocksize > 5) then
2617 [ # # ]: 0 : if(.not. gemm_nonlop_split_choice23) then
2618 : 0 : gemm_nonlop_split_choice23 = .true.
2619 : 0 : ndgxdt=6
2620 : : else
2621 : 0 : blocksize=blocksize*1.5
2622 : 0 : blocksize=min(nprojs,blocksize)
2623 : : end if
2624 : : else
2625 : 0 : blocksize=blocksize*1.5
2626 : : end if
2627 : :
2628 : : nonlop_smem = gemm_nonlop_ompgpu_static_mem(npw_fft,gs_hamk%indlmn,gs_hamk%nattyp,gs_hamk%ntypat,&
2629 : 0 : & blocksize,ndgxdt,use_distrib)
2630 : 0 : sum_mem = nonlop_smem + gs_ham_smem + nonlop_wmem + prep_nonlop_wmem
2631 : : end if
2632 : :
2633 [ # # # # ]: 0 : if(.not. l_disable_output .or. sum_other_mem > free_mem .or. sum_mem > free_mem) then
2634 : 0 : write(std_out,'(A,I3,A)') "GPU memory consumption estimate per MPI task for K-point ",ikpt,":"
2635 [ # # ]: 0 : if(blocksize>1) then
2636 [ # # ]: 0 : if(use_distrib) then
2637 : 0 : write(std_out,'(A,I3,A,I3,A)') "MPI distribution of GEMM nonlop projectors using ",&
2638 : 0 : & nblocks, " blocks of ", blocksize, " MPI tasks."
2639 : : else
2640 : 0 : write(std_out,'(A,I3,A)') "Local slicing of GEMM nonlop projectors using ",&
2641 : 0 : & blocksize, " blocks."
2642 : : end if
2643 : : end if
2644 [ # # # # ]: 0 : if(nfft_blocks>1 .and. wfoptalg>=0) then
2645 : 0 : write(std_out,'(A,I3,A)') "Local slicing of FFT work array using ",&
2646 : 0 : & nfft_blocks, " blocks."
2647 : : end if
2648 : 0 : write(std_out,'(A,F10.3,1x,A)') " Considered available memory : ", real(free_mem)/(1024*1024), "MiB"
2649 : 0 : write(std_out,'(A)')
2650 : 0 : write(std_out,'(A)') "| Buffers governed by blocking/slicing |"
2651 : 0 : write(std_out,'(A)') "|:--------------------------|---------:|-------------:|---------------:|"
2652 : 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 |"
2653 [ # # ]: 0 : if(wfoptalg>=0) then
2654 : 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 |"
2655 : 0 : write(std_out,'(A,F10.2,1x,A)') "| xFFT~internal buffers | NA | NA | ", real(fourwf_smem,dp)/(1024*1024), "MiB |"
2656 : : end if
2657 : :
2658 [ # # ]: 0 : if(sum_other_mem > free_mem) then
2659 : 0 : write(std_out,'(A)')
2660 : 0 : write(std_out,'(A)') "/!\ No slicing attempted as other arrays are too big to fit"
2661 : 0 : write(std_out,'(A)')
2662 : : end if
2663 : :
2664 : 0 : write(std_out,'(A)')
2665 : 0 : write(std_out,'(A)') "| Static buffers, computed once and permanently on card |"
2666 : 0 : write(std_out,'(A)') "|:-------------------------|--------------:|----------------:|"
2667 : : ! CHEBFI2 or SLICE
2668 [ # # ]: 0 : if(wfoptalg==111 .or. wfoptalg==112) then
2669 : 0 : write(std_out,'(A,F10.2,1x,A)') "| invovl (mkinvovl) | natom | ", real(invovl_smem,dp)/(1024*1024), "MiB |"
2670 : 0 : write(std_out,'(A,F10.2,1x,A)') "| chebfi2 | npw | ", real(chebfiMem(1))/(1024*1024), "MiB |"
2671 : : end if
2672 : :
2673 : : ! LOBPCG2
2674 [ # # ]: 0 : if(wfoptalg==114) then
2675 : 0 : write(std_out,'(A,F10.2,1x,A)') "| lobpcg2 | npw | ", real(lobpcgMem(1))/(1024*1024), "MiB |"
2676 : : end if
2677 : :
2678 : 0 : write(std_out,'(A,F10.2,1x,A)') "| hamiltonian arrays | npw | ", real(gs_ham_smem)/(1024*1024), "MiB |"
2679 : :
2680 : 0 : write(std_out,'(A)')
2681 : 0 : write(std_out,'(A)') "| Work buffers (mostly sized after bandpp or nblock_lobpcg) |"
2682 : 0 : write(std_out,'(A)') "|:-------------------------|--------------:|----------------:|"
2683 : : ! getghc (any diago algorithm)
2684 [ # # ]: 0 : if(wfoptalg>=0) then
2685 [ # # ]: 0 : if(getghc_wmem /= fourwf_wmem) then
2686 : 0 : write(std_out,'(A,F10.2,1x,A)') "| gemm_nonlop | bandpp | ", real(getghc_wmem,dp)/(1024*1024), "MiB |"
2687 : : end if
2688 : 0 : write(std_out,'(A,F10.2,1x,A)') "| mkrho~vtowfk_extra | npw,bandpp | ", real(updrho_wmem,dp)/(1024*1024), "MiB |"
2689 : 0 : write(std_out,'(A,F10.2,1x,A)') "| hegvd | bandpp | ", real(hegvd_mem,dp)/(1024*1024), "MiB |"
2690 : : else
2691 : 0 : write(std_out,'(A,F10.2,1x,A)') "| gemm_nonlop | natom,bandpp | ", real(nonlop_wmem,dp)/(1024*1024), "MiB |"
2692 : : end if
2693 [ # # ]: 0 : if(paral_kgb==1) then
2694 : 0 : write(std_out,'(A,F10.2,1x,A)') "| prep_nonlop | npw,bandpp | ", real(prep_nonlop_wmem,dp)/(1024*1024), "MiB |"
2695 : : end if
2696 : :
2697 : : ! CHEBFI2 or SLICE
2698 [ # # ]: 0 : if(wfoptalg==111 .or. wfoptalg==112) then
2699 : 0 : write(std_out,'(A,F10.2,1x,A)') "| invovl | natom,bandpp | ", real(invovl_wmem,dp)/(1024*1024), "MiB |"
2700 : 0 : write(std_out,'(A,F10.2,1x,A)') "| chebfi2 (RR buffers) | nband | ", real(chebfiMem(2))/(1024*1024), "MiB |"
2701 : 0 : write(std_out,'(A,F10.2,1x,A)') "| chebfiwf (cg,resid,eig) | npw,nband | ", real(localMem)/(1024*1024), "MiB |"
2702 : : end if
2703 : :
2704 : : ! LOBPCG2
2705 [ # # ]: 0 : if(wfoptalg==114) then
2706 : 0 : write(std_out,'(A,F10.2,1x,A)') "| lobpcg2 (RR buffers) | bandpp | ", real(lobpcgMem(2))/(1024*1024), "MiB |"
2707 : 0 : write(std_out,'(A,F10.2,1x,A)') "| lobpcgwf (cg,resid,eig) | npw,nband | ", real(localMem)/(1024*1024), "MiB |"
2708 : : end if
2709 : :
2710 : 0 : write(std_out,'(A)')
2711 [ # # ]: 0 : if(sum_other_mem > free_mem) then
2712 : 0 : write(std_out,'(A,F10.2,1x,A)') "Sum : ", real(sum_other_mem)/(1024*1024), "MiB"
2713 : : else
2714 : 0 : write(std_out,'(A,F10.2,1x,A)') "Sum : ", real(sum_mem)/(1024*1024), "MiB"
2715 : : end if
2716 : 0 : write(std_out,'(A)')
2717 : 0 : flush(std_out)
2718 : :
2719 : : end if
2720 : :
2721 [ # # ]: 0 : if(rank==0) then
2722 [ # # ]: 0 : if(sum_other_mem > free_mem) then
2723 : : write(message,'(3a)') &
2724 : 0 : & ' Your case is too big to fit in GPU memory regardless of possible array optimizations in fourwf and GEMM nonlop.',ch10,&
2725 : 0 : & ' Action : run on more nodes and/or increase nblock_lobpcg if using LOBPCG.'
2726 : 0 : ABI_ERROR(message)
2727 : : end if
2728 [ # # ]: 0 : if(sum_mem > free_mem) then
2729 [ # # ]: 0 : if(l_warn_on_fail) then
2730 : 0 : ABI_WARNING("It seems the test case you're trying to run is too big to run with given GPU resources !")
2731 : : else
2732 : 0 : ABI_ERROR("It seems the test case you're trying to run is too big to run with given GPU resources !")
2733 : : end if
2734 : : end if
2735 : : end if
2736 : : !call xmpi_barrier(xmpi_world)
2737 : :
2738 : : end subroutine get_gemm_nonlop_ompgpu_blocksize
2739 : : !!***
2740 : :
2741 : : end module m_common
2742 : : !!***
|
