Branch data Line data Source code
1 : : !!****m* ABINIT/m_scfcv_core
2 : : !! NAME
3 : : !! m_scfcv_core
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (XG, GMR, AR, MKV, MT, FJ, MB)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : ! nvtx related macro definition
23 : : #include "nvtx_macros.h"
24 : :
25 : : module m_scfcv_core
26 : :
27 : : use defs_basis
28 : : use defs_wvltypes
29 : : use defs_rectypes
30 : : use m_xmpi
31 : : use m_abicore
32 : : use m_wffile
33 : : use m_rec
34 : : use m_abi_mixing
35 : : use m_precon
36 : : use m_errors
37 : : use m_efield
38 : : use mod_prc_memory
39 : : use m_nctk
40 : : use m_hdr
41 : : use m_xcdata
42 : : use m_cgtools
43 : : use m_dtfil
44 : : use m_extfpmd
45 : : use m_invovl
46 : : use m_xg_nonlop
47 : : use m_rcpaw
48 : :
49 : : use m_nonlop, only : nonlop_counter
50 : : use defs_datatypes, only : pseudopotential_type
51 : : use defs_abitypes, only : MPI_type
52 : : use m_berryphase_new, only : update_e_field_vars
53 : : use m_dens, only : constrained_dft_t, constrained_dft_ini, constrained_dft_free, calcdenmagsph, calmaxdifmag
54 : : use m_time, only : timab
55 : : use m_fstrings, only : int2char4, sjoin, itoa
56 : : use m_symtk, only : symmetrize_xred
57 : : use m_geometry, only : metric
58 : : use m_fftcore, only : getng, sphereboundary
59 : : use m_time, only : abi_wtime, sec2str
60 : : use m_exit, only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in
61 : : use m_mpinfo, only : destroy_mpi_enreg, iwrite_fftdatar, initmpi_seq, proc_distrb_cycle
62 : : use m_ioarr, only : fftdatar_write_from_hdr
63 : : use m_results_gs , only : results_gs_type
64 : : use m_scf_history, only : scf_history_type, scf_history_init, scf_history_free
65 : : use m_energies, only : energies_type
66 : : use m_electronpositron, only : electronpositron_type, electronpositron_calctype
67 : : use m_pawang, only : pawang_type
68 : : use m_pawrad, only : pawrad_type
69 : : use m_pawtab, only : pawtab_type,pawtab_get_lsize
70 : : use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags
71 : : use m_pawxc, only : pawxc_get_usekden
72 : : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
73 : : use m_pawrhoij, only : pawrhoij_type
74 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, &
75 : : & pawcprj_free, pawcprj_axpby, pawcprj_put, pawcprj_getdim, pawcprj_reorder
76 : : use m_pawdij, only : pawdij, symdij
77 : : use m_pawfgr, only : pawfgr_type
78 : : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags
79 : : use m_paw_dmft, only : paw_dmft_type
80 : : use m_paw_nhat, only : nhatgrid,wvl_nhatgrid,pawmknhat
81 : : use m_paw_tools, only : chkpawovlp
82 : : use m_paw_denpot, only : pawdenpot,paw_relax_core
83 : : use m_paw_atomorb, only : orb_relaxed_core
84 : : use m_paw_occupancies, only : pawmkrhoij
85 : : use m_paw_correlations, only : setnoccmmp,setrhoijpbe0
86 : : use m_paw_mkrho, only : pawmkrho
87 : : use m_paw_uj, only : pawuj_red, macro_uj_type
88 : : use m_paw_dfpt, only : pawgrnl
89 : : use m_fock, only : fock_type, fock_init, fock_destroy, fock_ACE_destroy, fock_common_destroy, &
90 : : fock_BZ_destroy, fock_update_exc, fock_updatecwaveocc
91 : : use m_gemm_nonlop_projectors, only : reset_gemm_nonlop, gemm_nonlop_use_gemm
92 : : use m_gwls_hamiltonian, only : build_vxc
93 : : #if defined HAVE_BIGDFT
94 : : use BigDFT_API, only : cprj_clean,cprj_paw_alloc
95 : : #endif
96 : : use m_outxml, only : out_resultsgs_XML, out_geometry_XML
97 : : use m_kg, only : getcut, getmpw, kpgio, getph
98 : : use m_fft, only : fourdp
99 : : use m_vtorhorec, only : first_rec, vtorhorec
100 : : use m_vtorhotf, only : vtorhotf
101 : : use m_outscfcv, only : outscfcv
102 : : use m_afterscfloop, only : afterscfloop
103 : : use m_extraprho, only : extraprho
104 : : use m_spacepar, only : make_vectornd,setsym
105 : : use m_newrho, only : newrho
106 : : use m_newvtr, only : newvtr
107 : : use m_vtorho, only : vtorho
108 : : use m_setvtr, only : setvtr
109 : : use m_mkrho, only : mkrho
110 : : use m_rhotov, only : rhotov
111 : : use m_forces, only : fresid, forces
112 : : use m_dft_energy, only : energy, entropy
113 : : use m_initylmg, only : initylmg
114 : : use m_rhotoxc, only : rhotoxc
115 : : use m_drivexc, only : check_kxc, xc_need_kden
116 : : use m_odamix, only : odamix
117 : : use m_common, only : scprqt, prtene
118 : : use m_fourier_interpol, only : transgrid
119 : : use m_fock_getghc, only : fock2ACE
120 : : use m_forstr, only : nres2vres
121 : : use m_positron, only : setup_positron
122 : : use m_cgprj, only : ctocprj
123 : : use m_psolver, only : psolver_rhohxc
124 : : use m_paw2wvl, only : paw2wvl_ij, wvl_cprjreorder
125 : : use m_pstat, only : pstat_proc
126 : : use m_pspini, only : pspcor
127 : : use m_ewald, only : ewald
128 : : use m_atm2fft, only : atm2fft
129 : : use m_paw_correlations, only : loc_orbmom_cal
130 : :
131 : : #if defined(HAVE_GPU_MARKERS)
132 : : use m_nvtx_data
133 : : #endif
134 : :
135 : : implicit none
136 : :
137 : : private
138 : : !!***
139 : :
140 : : public :: scfcv_core
141 : : !!***
142 : :
143 : : contains
144 : : !!***
145 : :
146 : : !!****f* ABINIT/scfcv_core
147 : : !! NAME
148 : : !! scfcv_core
149 : : !!
150 : : !! FUNCTION
151 : : !! Self-consistent-field convergence.
152 : : !! Conducts set of passes or overall iterations of preconditioned
153 : : !! conjugate gradient algorithm to converge wavefunctions to
154 : : !! ground state and optionally to compute forces and energy.
155 : : !! This routine is called to compute forces for given atomic
156 : : !! positions or else to do non-SCF band structures.
157 : : !!
158 : : !! INPUTS
159 : : !! atindx(natom)=index table for atoms (see gstate.f)
160 : : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
161 : : !! cpus= cpu time limit in seconds
162 : : !! dmatpawu= fixed occupation matrix of correlated orbitals (DFT+U or DMFT only)
163 : : !! dtfil <type(datafiles_type)>=variables related to files
164 : : !! dtset <type(dataset_type)>=all input variables for this dataset
165 : : !! | mband=maximum number of bands
166 : : !! | mgfft=maximum size of 1D FFTs for the "coarse" grid (see NOTES below)
167 : : !! | mkmem =number of k points treated by this node.
168 : : !! | mpw=maximum dimensioned size of npw.
169 : : !! | natom=number of atoms in cell.
170 : : !! | nfft=(effective) number of FFT grid points (for this processor)
171 : : !! | for the "coarse" grid (see NOTES below)
172 : : !! | nkpt=number of k points
173 : : !! | nspden=number of spin-density components
174 : : !! | nsppol=1 for unpolarized, 2 for spin-polarized
175 : : !! | nsym=number of symmetry elements in space group
176 : : !! ecore=core psp energy (part of total energy) (hartree)
177 : : !! fatvshift=factor to multiply dtset%atvshift
178 : : !! itimes(2)=itime array, contain itime=itimes(1) and itimimage_gstate=itimes(2) from outer loops
179 : : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
180 : : !! mcg=size of wave-functions array (cg) =mpw*my_nspinor*mband*mkmem*nsppol
181 : : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
182 : : !! mpi_enreg=information about MPI parallelization
183 : : !! my_natom=number of atoms treated by current processor
184 : : !! nattyp(ntypat)= # atoms of each type.
185 : : !! ndtpawuj=size of dtpawuj
186 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
187 : : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
188 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
189 : : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
190 : : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and
191 : : !! related data
192 : : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
193 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
194 : : !! | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
195 : : !! pwind(pwind_alloc,2,3) = array used to compute
196 : : !! the overlap matrix smat between k-points (see initberry.f)
197 : : !! pwind_alloc = first dimension of pwind
198 : : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
199 : : !! (see initberry.f)
200 : : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for
201 : : !! each G and k point
202 : : !! ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real
203 : : !! spherical harmonics
204 : : !!
205 : : !! OUTPUT
206 : : !! resid(mband*nkpt*nsppol)=residuals for each band over all k points and spins
207 : : !!
208 : : !! SIDE EFFECTS
209 : : !! cg(2,mcg)=updated wavefunctions; if mkmem>=nkpt, these are kept in a disk file.
210 : : !! cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn>
211 : : !! dtefield <type(efield_type)> = variables related to Berry phase
212 : : !! dtpawuj(ndtpawuj)= data used for the automatic determination of U
213 : : !! (relevant only for PAW+U) calculations (see initberry.f)
214 : : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
215 : : !! electronpositron <type(electronpositron_type)>=quantities for
216 : : !! the electron-positron annihilation
217 : : !! extfpmd <type(extfpmd_type)>=extended first-principles molecular dynamics type
218 : : !! hdr <type(hdr_type)>=the header of wf, den and pot files
219 : : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
220 : : !! initialized= if 0 the initialization of the gstate run is not yet finished
221 : : !! irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
222 : : !! nfftf=(effective) number of FFT grid points (for this processor)
223 : : !! for the "fine" grid (see NOTES below)
224 : : !! occ(mband*nkpt*nsppol)=occupation number for each band (often 2) at each k point
225 : : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
226 : : !! phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
227 : : !! results_gs <type(results_gs_type)>=results (energy and its components,
228 : : !! forces and its components, the stress tensor) of a ground-state
229 : : !! computation (should be made a pure output quantity)
230 : : !! rhog(2,nfftf)=array for Fourier transform of electron density
231 : : !! rhor(nfftf,nspden)=array for electron density in el./bohr**3
232 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
233 : : !! scf_history <type(scf_history_type)>=arrays obtained from previous
234 : : !! SCF cycles
235 : : !! symrec(3,3,nsym)=symmetry operations in reciprocal space
236 : : !! taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic
237 : : !! energy density
238 : : !! taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
239 : : !! wffnew=struct info for wf disk files
240 : : !! wvl <type(wvl_data)>=all wavelets data
241 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
242 : : !! xred_old(3,natom)= at input, previous reduced dimensionless atomic
243 : : !! coordinates at output, current xred is transferred to xred_old
244 : : !! conv_retcode=return code, 0 if convergence was achieved
245 : : !!
246 : : !! NOTES
247 : : !! It is worth to explain THE USE OF FFT GRIDS:
248 : : !! ============================================
249 : : !! In case of PAW:
250 : : !! ---------------
251 : : !! Two FFT grids are used:
252 : : !! - A "coarse" FFT grid (defined by ecut)
253 : : !! for the application of the Hamiltonian on the plane waves basis.
254 : : !! It is defined by nfft, ngfft, mgfft, ...
255 : : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ...
256 : : !! are expressed on this grid.
257 : : !! - A "fine" FFT grid (defined) by ecutdg)
258 : : !! for the computation of the density inside PAW spheres.
259 : : !! It is defined by nfftf, ngfftf, mgfftf, ...
260 : : !! Total density, potentials, ...
261 : : !! are expressed on this grid.
262 : : !! In case of norm-conserving:
263 : : !! ---------------------------
264 : : !! - Only the usual FFT grid (defined by ecut) is used.
265 : : !! It is defined by nfft, ngfft, mgfft, ...
266 : : !! For compatibility reasons, (nfftf,ngfftf,mgfftf)
267 : : !! are set equal to (nfft,ngfft,mgfft) in that case.
268 : : !!
269 : : !! SOURCE
270 : :
271 : 13878 : subroutine scfcv_core(atindx,atindx1,cg,cprj,cpus,dmatpawu,dtefield,dtfil,dtpawuj,&
272 : 13878 : & dtset,ecore,eigen,electronpositron,fatvshift,hdr,extfpmd,rcpaw,indsym,&
273 : 13878 : & initialized,irrzon,itimes,kg,mcg,mcprj,mpi_enreg,my_natom,nattyp,ndtpawuj,nfftf,npwarr,occ,&
274 : 13878 : & paw_dmft,pawang,pawfgr,pawrad,pawrhoij,pawtab,phnons,psps,pwind,&
275 : 13878 : & pwind_alloc,pwnsfac,rec_set,resid,results_gs,rhog,rhor,rprimd,&
276 [ - + ]: 6939 : & scf_history,symrec,taug,taur,wffnew,wvl,xg_nonlop,xred,xred_old,ylm,ylmgr,conv_retcode)
277 : :
278 : : !Arguments ------------------------------------
279 : : !scalars
280 : : integer,intent(in) :: mcg,my_natom,ndtpawuj,pwind_alloc
281 : : integer,intent(inout) :: initialized,nfftf,mcprj
282 : : integer,intent(out) :: conv_retcode
283 : : real(dp),intent(in) :: cpus,fatvshift
284 : : real(dp),intent(inout) :: ecore
285 : : type(MPI_type),intent(inout) :: mpi_enreg
286 : : type(datafiles_type),intent(in) :: dtfil
287 : : type(dataset_type),intent(inout) :: dtset
288 : : type(efield_type),intent(inout) :: dtefield
289 : : type(electronpositron_type),pointer:: electronpositron
290 : : type(hdr_type),intent(inout) :: hdr
291 : : type(extfpmd_type),pointer,intent(inout) :: extfpmd
292 : : type(rcpaw_type), pointer,intent(inout) :: rcpaw
293 : : type(pawang_type),intent(in) :: pawang
294 : : type(pawfgr_type),intent(inout) :: pawfgr
295 : : type(pseudopotential_type),intent(inout) :: psps
296 : : type(recursion_type),intent(inout) :: rec_set
297 : : type(results_gs_type),intent(inout) :: results_gs
298 : : type(scf_history_type),intent(inout) :: scf_history
299 : : type(wffile_type),intent(inout) :: wffnew
300 : : type(wvl_data),intent(inout) :: wvl
301 : : !arrays
302 : : integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom)
303 : : integer,intent(in) :: indsym(4,dtset%nsym,dtset%natom),itimes(2)
304 : : !no_abirules
305 : : integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
306 : : !(nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise)
307 : : integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem)
308 : : integer, intent(in) :: nattyp(psps%ntypat),npwarr(dtset%nkpt),pwind(pwind_alloc,2,3)
309 : : integer, intent(in) :: symrec(3,3,dtset%nsym)
310 : : real(dp), intent(inout) :: cg(2,mcg),dmatpawu(:,:,:,:)
311 : : real(dp), intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
312 : : real(dp), intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
313 : : real(dp), intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
314 : : !(nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise)
315 : : real(dp), intent(in) :: pwnsfac(2,pwind_alloc)
316 : : real(dp), intent(inout) :: rprimd(3,3)
317 : : real(dp), pointer :: rhog(:,:),rhor(:,:)
318 : : real(dp), pointer :: taug(:,:),taur(:,:)
319 : : real(dp), intent(inout) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
320 : : real(dp), intent(inout) :: xred(3,dtset%natom)
321 : : real(dp), intent(inout) :: xred_old(3,dtset%natom)
322 : : real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
323 : : real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
324 : : type(macro_uj_type),intent(inout) :: dtpawuj(0:ndtpawuj)
325 : : type(pawrhoij_type), intent(inout) :: pawrhoij(my_natom*psps%usepaw)
326 : : type(pawrad_type), intent(in) :: pawrad(psps%ntypat*psps%usepaw)
327 : : type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
328 : : type(paw_dmft_type), intent(inout) :: paw_dmft
329 : : type(pawcprj_type),pointer, intent(inout) :: cprj(:,:)
330 : : type(xg_nonlop_t),intent(inout) :: xg_nonlop
331 : : !Local variables -------------------------
332 : : !scalars
333 : : integer,parameter :: level=110,response=0,cplex1=1
334 : : integer :: afford,bantot,choice
335 : : integer :: computed_forces,cplex,cplex_hf,ctocprj_choice,dbl_nnsclo,dielop,dielstrt,dimdmat
336 : : integer :: forces_needed,errid,has_dijhat,has_dijnd,has_dijU,has_vhartree,has_dijfock
337 : : integer :: history_size
338 : : integer :: iatom,ider,idir,ierr,ii,ikpt,impose_dmat,denpot
339 : : integer :: initialized0,iorder_cprj,ipert,ipositron,isave_den,isave_kden,iscf10,ispden
340 : : integer :: ispmix,istep,istep_fock_outer,istep_mix,istep_updatedfock,itypat,izero,lmax_diel,lpawumax,mband_cprj
341 : : #if defined HAVE_BIGDFT
342 : : integer :: mcprj_wvl
343 : : #endif
344 : : integer :: me,me_wvl,mgfftdiel,mgfftf,moved_atm_inside,moved_rhor,my_nspinor,n1xccc
345 : : integer :: n3xccc,ncpgr,nfftdiel,nfftmix,nfftmix_per_nfft,nfftotf,ngrcondft,ngrvdw,nhatgrdim,nk3xc,nkxc
346 : : integer :: npawmix,npwdiel,nremit,nstep,nzlmopt,optcut,optcut_hf,optene,optgr0,optgr0_hf
347 : : integer :: optgr1,optgr2,optgr1_hf,optgr2_hf,option,optrad,optrad_hf,optres,optxc,prtfor,prtxml,quit
348 : : integer :: quit_sum,rdwrpaw,shift,spaceComm,spaceComm_fft,spaceComm_wvl,spaceComm_grid
349 : : integer :: spare_mem,optn
350 : : integer :: stress_needed,sz1,sz2,tim_mkrho,unit_out
351 : : integer :: usecprj,usevxctau,usevxctau_paw,usexcnhat,use_hybcomp
352 : : integer :: my_quit,quitsum_request,timelimit_exit,usecg,wfmixalg,with_vectornd
353 : : integer ABI_ASYNC :: quitsum_async
354 : : real(dp) :: boxcut,compch_fft,compch_sph,deltae,diecut,diffor,ecut
355 : : real(dp) :: ecutf,ecutsus,edum,sdum,elast,etotal,evxc,fermie,fermih,gsqcut,hyb_mixing,hyb_mixing_sr
356 : : real(dp) :: maxfor,res2,residm,ucvol,ucvol_local,val_max
357 : : real(dp) :: val_min,vxcavg,vxcavg_dum,el_temp
358 : : real(dp) :: zion,wtime_step,now,prev,esum,enonlocalpsp !MRM
359 : : real(dp) :: invovl_precision_orig,invovl_dyn_ratio
360 : : ! Margin (in units of the active tolerance) by which invovl_precision is
361 : : ! tightened to its final value ahead of actual SCF convergence.
362 : : real(dp),parameter :: invovl_dyn_prec_safety=0.1_dp
363 : : character(len=10) :: tag
364 : : character(len=500) :: MY_NAME = "scfcv_core"
365 : : character(len=1500) :: msg
366 : : !character(len=500) :: dilatmx_errmsg
367 : : character(len=fnlen) :: fildata
368 : 6939 : type(MPI_type) :: mpi_enreg_diel
369 : : type(xcdata_type) :: xcdata
370 : : type(energies_type), target :: energies
371 : : type(abi_mixing_object) :: mix,mix_mgga
372 : 6939 : type(precon_object) :: precon
373 : : logical,parameter :: VERBOSE=.FALSE.
374 : : logical :: dummy_nhatgr
375 : : logical :: finite_efield_flag=.false.
376 : : logical :: non_magnetic_xc=.false.
377 : : logical :: recompute_cprj=.false.,reset_mixing=.false.
378 : : logical,save :: tfw_activated=.false.
379 : : logical :: wvlbigdft=.false.
380 : : !type(energies_type),pointer :: energies_wvl ! TO BE ACTIVATED LATER
381 : : !arrays
382 : : integer :: ngfft(18),ngfftdiel(18),ngfftf(18),ngfftmix(18),npwarr_diel(1)
383 : : integer :: npwtot_diel(1)
384 : : integer, save :: scfcv_jdtset = 0 ! To simulate iapp behavior
385 : : integer, save :: scfcv_itime = 1 ! To simulate iapp behavior
386 : 6939 : integer,allocatable :: dimcprj(:),dimcprj_srt(:)
387 : 6939 : integer,allocatable :: gbound_diel(:,:),irrzondiel(:,:,:),kg_diel(:,:)
388 : 6939 : integer,allocatable :: l_size_atm(:)
389 : 6939 : integer,allocatable :: indsym_dum(:,:,:),symrec_dum(:,:,:), rmm_diis_status(:,:,:)
390 : 13878 : logical :: eijkl_is_sym(psps%ntypat)
391 : 6939 : logical,pointer :: lmselect_ep(:,:)
392 : : real(dp) :: dielar(7),dphase(3),favg(3),gmet(3,3),gprimd(3,3)
393 : : real(dp) :: kpt_diel(3),pel(3),pel_cg(3),pelev(3),pion(3),ptot(3),qpt(3),red_ptot(3) !!REC
394 : : real(dp) :: rhodum(1),rmet(3,3),strscondft(6),strsxc(6),strten(6),tollist(12)
395 : 13878 : real(dp) :: tsec(2),vnew_mean(dtset%nspden),vres_mean(dtset%nspden)
396 : : real(dp) :: efield_old_cart(3), ptot_cart(3)
397 : : real(dp) :: red_efield2(3),red_efield2_old(3)
398 : : real(dp) :: vpotzero(2)
399 : : real(dp) :: maxmag , difmag
400 : : real(dp) :: dmatdum(0,0,0,0)
401 : 13878 : real(dp) :: orb_mom_atom(10,3,dtset%natom)
402 : : ! red_efield1(3),red_efield2(3) is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009) [[cite:Stengel2009]]
403 : : ! red_efield1(3) for fixed ebar calculation, red_efield2(3) for fixed reduced d calculation, in mixed BC
404 : : ! red_efieldbar_lc(3) is local reduced electric field, defined by Eq.(28) of Nat. Phys. suppl. (2009) [[cite:Stengel2009]]
405 : : ! pbar(3) and dbar(3) are reduced polarization and displacement field,
406 : : ! defined by Eq.(27) and (29) Nat. Phys. suppl. (2009) [[cite:Stengel2009]]
407 : : real(dp) :: dummy01(0),dummy02(0),dummy03(0),dummy04(0),dummy05(0),dummy06(0),dummy07(0),dummy08(0),dummy09(0),dummy010(0)
408 : : real(dp) :: dummy61(6),dummy62(6)
409 : : real(dp),parameter :: k0(3)=(/zero,zero,zero/)
410 : 6939 : real(dp),allocatable :: dielinv(:,:,:,:,:),dtn_pc(:,:)
411 : 6939 : real(dp),allocatable :: fcart(:,:),forold(:,:),gred(:,:),gresid(:,:)
412 : 6939 : real(dp),allocatable :: grchempottn(:,:),grcondft(:,:),grewtn(:,:)
413 : 6939 : real(dp),allocatable :: grhf(:,:),grnl(:),grvdw(:,:),grxc(:,:)
414 : 6939 : real(dp),allocatable :: intgres(:,:),kxc(:,:),nhat(:,:),nhatgr(:,:,:),nvresid(:,:),nvtauresid(:,:)
415 : 6939 : real(dp),allocatable :: ph1d(:,:),ph1ddiel(:,:),ph1df(:,:)
416 : 6939 : real(dp),allocatable :: phnonsdiel(:,:,:),rhowfg(:,:),rhowfr(:,:),shiftvector(:)
417 : 6939 : real(dp),allocatable :: susmat(:,:,:,:,:),synlgr(:,:)
418 : 6939 : real(dp),allocatable :: vectornd(:,:,:),vhartr(:),vpsp(:),vtrial(:,:)
419 : 6939 : real(dp),allocatable :: vxc(:,:),vxc_hybcomp(:,:),vxctau(:,:,:),workr(:,:),xccc3d(:),xcctau3d(:),ylmdiel(:,:)
420 : 6939 : real(dp),allocatable :: intgden(:,:),intgden0(:,:)
421 : 6939 : real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:)
422 : 6939 : type(scf_history_type) :: scf_history_wf
423 : 6939 : type(constrained_dft_t) :: constrained_dft
424 : 6939 : type(paw_an_type),allocatable :: paw_an(:)
425 : 6939 : type(paw_ij_type),allocatable :: paw_ij(:)
426 : 6939 : type(pawfgrtab_type),allocatable,save :: pawfgrtab(:)
427 : 6939 : type(pawrhoij_type),pointer :: pawrhoij_ep(:)
428 : : type(fock_type),pointer :: fock
429 : 6939 : type(pawcprj_type),allocatable, target :: cprj_local(:,:)
430 : :
431 : : ! *********************************************************************
432 : :
433 : : !write(std_out,'(a,5i4)')' scfcv_core, enter : itimes(1:2)=',itimes(1:2)
434 : : DBG_ENTER("COLL")
435 : :
436 : 6939 : call timab(1440,1,tsec)
437 : 6939 : call timab(1441,3,tsec)
438 : :
439 : : ! enable time limit handler if not done in callers.
440 [ + + ]: 6939 : if (enable_timelimit_in(MY_NAME) == MY_NAME) then
441 [ + - ]: 6943 : write(std_out,*)"Enabling timelimit check in function: ",trim(MY_NAME)," with timelimit: ",trim(sec2str(get_timelimit()))
442 : : end if
443 : :
444 : : ! Initialise non_magnetic_xc for rhohxc
445 : 6939 : non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
446 : :
447 : : !######################################################################
448 : : !Initializations - Memory allocations
449 : : !----------------------------------------------------------------------
450 : 6939 : lmax_diel = 0
451 : :
452 : : !MPI communicators
453 [ - + ]: 6939 : if (xmpi_paral==1.and.mpi_enreg%paral_hf==1) then
454 : 0 : spaceComm=mpi_enreg%comm_kpt
455 : : else
456 : 6939 : spaceComm=mpi_enreg%comm_cell
457 : : end if
458 : 6939 : me=xmpi_comm_rank(spaceComm)
459 : 6939 : spaceComm_fft=mpi_enreg%comm_fft
460 : 6939 : spaceComm_wvl=mpi_enreg%comm_wvl
461 : 6939 : me_wvl=mpi_enreg%me_wvl
462 : 6939 : spaceComm_grid=mpi_enreg%comm_fft
463 [ - + ]: 6939 : if(dtset%usewvl==1) spaceComm_grid=mpi_enreg%comm_wvl
464 : 6939 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
465 : :
466 : : !Save some variables from dataset definition
467 : 6939 : nstep=dtset%nstep
468 : 6939 : ecut=dtset%ecut
469 [ + + ]: 6939 : ecutf=ecut; if (psps%usepaw==1) ecutf=dtset%pawecutdg
470 [ + + + + ]: 6939 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) ecutf=dtset%pawecutdg
471 : 6939 : iscf10=mod(dtset%iscf,10)
472 : 6939 : tollist(1)=dtset%tolmxf;tollist(2)=dtset%tolwfr
473 : 6939 : tollist(3)=dtset%toldff;tollist(4)=dtset%toldfe
474 : 6939 : tollist(6)=dtset%tolvrs;tollist(7)=dtset%tolrff
475 : 6939 : tollist(8)=dtset%vdw_df_threshold
476 : 6939 : tollist(9)=dtset%toldmag
477 : 6939 : dielstrt=0
478 : : finite_efield_flag=(dtset%berryopt == 4 .or. &
479 : : & dtset%berryopt == 6 .or. &
480 : : & dtset%berryopt == 7 .or. &
481 : : & dtset%berryopt == 14 .or. &
482 : : & dtset%berryopt == 16 .or. &
483 [ + + ]: 6939 : & dtset%berryopt == 17)
484 : :
485 : : !Get electronic temperature from dtset
486 [ + + + - : 13836 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
- + ]
487 : :
488 : : !Get FFT grid(s) sizes (be careful !)
489 : : !See NOTES in the comments at the beginning of this file.
490 [ + + ]: 131841 : ngfft(:)=dtset%ngfft(:)
491 [ + + ]: 6939 : if (psps%usepaw==1) then
492 [ + + ]: 28899 : mgfftf=pawfgr%mgfft;ngfftf(:)=pawfgr%ngfft(:)
493 : : else
494 : 5418 : mgfftf=dtset%mgfft;ngfftf(:)=ngfft(:)
495 : : end if
496 : :
497 : : !Calculate zion: the total positive charge acting on the valence electrons
498 : 6939 : zion=zero
499 [ + + ]: 27313 : do iatom=1,dtset%natom
500 : 27313 : zion=zion+psps%ziontypat(dtset%typat(iatom))
501 : : end do
502 : :
503 : : !Compute different geometric tensor, as well as ucvol, from rprimd
504 : 6939 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
505 : :
506 : : !Fock: be sure that the pointer is initialized to Null.
507 : 6939 : nullify(fock)
508 : :
509 : : !If Chebishev Filtering or Slicing algo is used, init invovl routine structure
510 : : if((dtset%wfoptalg == 1 .or. dtset%wfoptalg == 111 .or. dtset%wfoptalg == 112) .and. psps%usepaw == 1 &
511 [ + + + + : 6939 : & .and. dtset%cprj_in_memory==0) then
+ + + + ]
512 : 96 : call init_invovl(dtset%nkpt)
513 : : end if
514 : :
515 : : !Special care in case of WVL
516 : : !wvlbigdft indicates that the BigDFT workflow will be followed
517 [ - + - - ]: 6939 : wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
518 : : !if (wvlbigdft) then ! TO BE ACTIVATED LATER
519 : : ! ABI_MALLOC(energies_wvl,)
520 : : !end if
521 : 6939 : ucvol_local = ucvol
522 : : #if defined HAVE_BIGDFT
523 : : if (dtset%usewvl == 1) then
524 : : ! We need to tune the volume when wavelets are used because, not
525 : : ! all FFT points are used.
526 : : ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3)
527 : : ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(product(wvl%den%denspot%dpbox%ndims), dp)
528 : : end if
529 : : #endif
530 : :
531 : : !Some variables need to be initialized/nullified at start
532 : 6939 : nullify(grhor,lrhor,elfr)
533 : 6939 : quit=0 ; dbl_nnsclo=0 ; conv_retcode=0
534 : 6939 : dielop=0 ; strsxc=zero
535 : 6939 : deltae=zero ; elast=zero ;
536 : 6939 : vpotzero(:)=zero
537 : : ! JWZ April 12 2018: Intel 18 compiler seems to require maxfor initialized,
538 : : ! else it dies in scprqt in some scenarios
539 : 6939 : maxfor=zero
540 : : ! Save the user-requested (tightest) invovl_precision: restored at the end
541 : : ! of the SCF loop, and used as the floor of the dynamic-precision schedule.
542 : 6939 : invovl_precision_orig=dtset%invovl_precision
543 : : !
544 : 6939 : results_gs%residm=zero;results_gs%res2=zero
545 : 6939 : results_gs%deltae=zero;results_gs%diffor=zero
546 : 6939 : call energies%init()
547 [ + + + + ]: 6939 : if (dtset%positron/=0.and.initialized/=0) then
548 : 3 : energies%e0_electronpositron =results_gs%energies%e0_electronpositron
549 : 3 : energies%e_electronpositron =results_gs%energies%e_electronpositron
550 : 3 : energies%edc_electronpositron=results_gs%energies%edc_electronpositron
551 : : maxfor=zero
552 : : end if
553 : :
554 : : ! Initialize fermi level.
555 [ + + + + ]: 6939 : if (dtset%nstep==0 .or. dtset%iscf < 0) then
556 : : !if ((dtset%nstep==0 .or. dtset%iscf < 0) .and. dtset%plowan_compute==0) then
557 : 531 : energies%e_fermie = results_gs%energies%e_fermie
558 : 531 : results_gs%fermie = results_gs%energies%e_fermie
559 : 531 : energies%e_fermih = results_gs%energies%e_fermih
560 : 531 : results_gs%fermih = results_gs%energies%e_fermih
561 : : ! End CP addition
562 : : end if
563 : :
564 : 13877 : select case (dtset%usepotzero)
565 : : case (0,1)
566 : 6938 : energies%e_corepsp = ecore / ucvol
567 : 6938 : energies%e_corepspdc = zero
568 : : case (2)
569 : : ! No need to include the PspCore energy since it is already included in the
570 : : ! local pseudopotential (vpsp)
571 : 1 : energies%e_corepsp = zero
572 [ + + - ]: 6939 : energies%e_corepspdc = zero
573 : : end select
574 [ - + ]: 6939 : if(wvlbigdft) energies%e_corepsp = zero
575 [ + + ]: 6939 : if(dtset%icutcoul.ne.3) energies%e_corepsp = zero
576 : :
577 : 6939 : fermie=energies%e_fermie
578 : 6939 : fermih=energies%e_fermih
579 : 6939 : isave_den=0; isave_kden=0 !initial index of density protection file
580 [ + + ]: 6939 : optres=merge(0,1,dtset%iscf<10)
581 : 6939 : usexcnhat=0!;mcprj=0
582 : 6939 : initialized0=initialized
583 [ - + ]: 6939 : if (dtset%tfkinfunc==12) tfw_activated=.true.
584 : 6939 : ipert=0;idir=0;cplex=1
585 : 6939 : istep_mix=1
586 : 6939 : istep_fock_outer=1
587 : 6939 : ipositron=electronpositron_calctype(electronpositron)
588 [ + + ]: 6939 : unit_out=0;if (dtset%prtvol >= 10) unit_out=ab_out
589 [ + + ]: 27756 : nfftotf=product(ngfftf(1:3))
590 : :
591 : 6939 : usecprj=0
592 [ + + ]: 6939 : if (mcprj>0) then
593 : 78 : usecprj=1
594 : : end if
595 : :
596 : : !Stresses and forces flags
597 : 6939 : forces_needed=0;prtfor=0
598 : : if ((dtset%optforces==1.or.dtset%ionmov==4.or.dtset%ionmov==5.or.&
599 [ + + + - : 6939 : & abs(tollist(3))>tiny(0._dp)).or.abs(tollist(7))>tiny(0._dp)) then
+ - + - -
+ ]
600 [ + + + - ]: 3052 : if (dtset%iscf>0.and.nstep>0) forces_needed=1
601 [ + + ]: 3052 : if (nstep==0) forces_needed=2
602 : 3052 : prtfor=1
603 [ + + + + ]: 3887 : else if (dtset%iscf>0.and.dtset%optforces==2) then
604 : 3270 : forces_needed=2
605 : : end if
606 : :
607 : 6939 : stress_needed=0
608 [ + + + + : 6939 : if (dtset%optstress>0.and.dtset%iscf>0.and.dtset%prtstm==0.and. (nstep>0.or.dtfil%ireadwf==1)) stress_needed=1
+ + + + +
+ ]
609 : : if (dtset%optstress>0.and.dtset%iscf>0.and.psps%usepaw==1 &
610 [ + + + + : 6939 : & .and.finite_efield_flag.and.(nstep>0.or.dtfil%ireadwf==1)) stress_needed=1
+ + + + -
+ - - ]
611 : :
612 : : !This is only needed for the tddft routine, and does not
613 : : !correspond to the intended use of results_gs (should be only
614 : : !for output of scfcv_core
615 : 6939 : etotal = results_gs%etotal
616 : :
617 : : !Entering a scfcv_core loop, printing data to XML file if required.
618 [ + + + + ]: 6939 : prtxml=0;if (me==0.and.dtset%prtxml==1) prtxml=1
619 [ + + ]: 6939 : if (prtxml == 1) then
620 : : ! scfcv_core() will handle a scf loop, so we output the scfcv markup.
621 : 4 : write(ab_xml_out, "(A)") ' <scfcvLoop>'
622 : 4 : write(ab_xml_out, "(A)") ' <initialConditions>'
623 : : ! We output the geometry of the dataset given in argument.
624 : : ! xred and rprimd are given independently since dtset only
625 : : ! stores original and final values.
626 : 4 : call out_geometry_XML(dtset, 4, dtset%natom, rprimd, xred)
627 : 4 : write(ab_xml_out, "(A)") ' </initialConditions>'
628 : : end if
629 : :
630 : : !Examine tolerance criteria, and eventually print a line to the output
631 : : !file (with choice=1, the only non-dummy arguments of scprqt are
632 : : !nstep, tollist and iscf - still, diffor and res2 are here initialized to 0)
633 : 6939 : choice=1 ; diffor=zero ; res2=zero
634 [ - + + - : 20817 : ABI_MALLOC(fcart,(3,dtset%natom))
+ - + - -
+ - + ]
635 [ - + + - : 13878 : ABI_MALLOC(gred,(3,dtset%natom))
- + - + ]
636 [ - + + - : 27756 : ABI_MALLOC(intgden,(dtset%nspden,dtset%natom))
- + + - +
- + - - +
- + ]
637 [ - + + - : 20817 : ABI_MALLOC(intgden0,(dtset%nspden,dtset%natom))
- + + - -
+ - + ]
638 [ + + + + ]: 88435 : gred(:,:)=zero
639 [ + + + + ]: 88435 : fcart(:,:)=results_gs%fcart(:,:) ! This is a side effect ...
640 [ + + + + ]: 49960 : intgden(:,:)=zero
641 [ + + + + : 699655 : orb_mom_atom=zero
+ + ]
642 : : !results_gs should not be used as input of scfcv_core
643 : : !HERE IS PRINTED THE FIRST LINE OF SCFCV
644 : :
645 : : call scprqt(choice,cpus,deltae,diffor,maxmag,difmag,dtset,&
646 : : & eigen,etotal,favg,fcart,energies%e_fermie,energies%e_fermih,dtfil%fnameabo_app_eig,&
647 : : & dtfil%filnam_ds(1),initialized0,dtset%iscf,istep,istep_fock_outer,istep_mix,dtset%kptns,&
648 : : & maxfor,moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,nstep,&
649 : : & occ,optres,prtfor,prtxml,quit,res2,resid,residm,response,tollist,&
650 : 6939 : & psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode)
651 : :
652 : : !Various allocations (potentials, gradients, ...)
653 [ - + + - : 20817 : ABI_MALLOC(forold,(3,dtset%natom))
+ - + - -
+ - + ]
654 [ - + + - : 13878 : ABI_MALLOC(grchempottn,(3,dtset%natom))
- + - + ]
655 [ - + + - : 13878 : ABI_MALLOC(grcondft,(3,dtset%natom))
- + - + ]
656 [ - + + - : 13878 : ABI_MALLOC(gresid,(3,dtset%natom))
- + - + ]
657 [ - + + - : 13878 : ABI_MALLOC(grewtn,(3,dtset%natom))
- + - + ]
658 [ + - - + : 20817 : ABI_MALLOC(grnl,(3*dtset%natom))
+ - + - +
- - + -
+ ]
659 [ - + + - : 13878 : ABI_MALLOC(grxc,(3,dtset%natom))
- + - + ]
660 [ - + + - : 13878 : ABI_MALLOC(synlgr,(3,dtset%natom))
- + - + ]
661 [ + - - + : 20817 : ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
+ - + - +
- - + - +
- + ]
662 [ + - - + : 20817 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*dtset%natom))
+ - + - +
- - + -
+ ]
663 [ - + + - : 20817 : ABI_MALLOC(vhartr,(nfftf))
+ - + - -
+ - + ]
664 [ - + + - : 27756 : ABI_MALLOC(vtrial,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
665 [ - + + - : 13878 : ABI_MALLOC(vpsp,(nfftf))
- + - + ]
666 [ - + + - : 20817 : ABI_MALLOC(vxc,(nfftf,dtset%nspden))
- + + - -
+ - + -
+ ]
667 [ + + ]: 6939 : usevxctau=merge(1,0,xc_need_kden(dtset%ixc))
668 [ - + + - : 34746 : ABI_MALLOC(vxctau,(nfftf,dtset%nspden,4*usevxctau))
- + + - +
+ + + + -
+ - + + -
+ - + -
+ ]
669 : :
670 : 6939 : wfmixalg=dtset%fockoptmix/100
671 : 6939 : use_hybcomp=0
672 [ + + ]: 6939 : if(mod(dtset%fockoptmix,100)==11)use_hybcomp=1
673 [ - + + - : 27756 : ABI_MALLOC(vxc_hybcomp,(nfftf,dtset%nspden*use_hybcomp))
+ + + - +
- + + - +
- + ]
674 : :
675 [ + + ]: 6939 : ngrvdw=0;if (dtset%vdw_xc>=5.and.dtset%vdw_xc<=7) ngrvdw=dtset%natom
676 [ + + + - : 20817 : ABI_MALLOC(grvdw,(3,ngrvdw))
+ - + + -
+ - + ]
677 : :
678 : 6939 : ngrcondft=0
679 [ + + + + : 16454 : if(any(dtset%constraint_kind(:)/=0)) ngrcondft=dtset%natom
+ + ]
680 [ - + + - : 27756 : ABI_MALLOC(intgres,(dtset%nspden,ngrcondft))
+ + + - +
- + + - +
- + ]
681 [ + + ]: 6939 : if(ngrcondft/=0)then
682 [ + + + + ]: 606 : intgres(:,:)=zero
683 : : endif
684 : :
685 [ + + + + ]: 88435 : grchempottn(:,:)=zero
686 [ + + + + ]: 88435 : grcondft(:,:)=zero
687 [ + + + + : 169931 : forold(:,:)=zero ; gresid(:,:)=zero ; pel(:)=zero
+ + + + ]
688 : 6939 : strscondft(:)=zero
689 [ + + + + : 369808547 : vtrial(:,:)=zero; vxc(:,:)=zero
+ + + + ]
690 [ + + ]: 6939 : n1xccc=0;if (psps%n1xccc/=0) n1xccc=psps%n1xccc
691 [ + + ]: 6939 : n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf
692 [ + + + - : 20817 : ABI_MALLOC(xccc3d,(n3xccc))
+ - + + -
+ - + ]
693 [ + + + - : 20817 : ABI_MALLOC(xcctau3d,(nfftf*dtset%usekden))
+ - + + -
+ - + ]
694 : :
695 : : !Allocations/initializations for PAW only
696 : 6939 : lpawumax=-1
697 [ + + ]: 6939 : if(psps%usepaw==1) then
698 : : ! Variables/arrays related to the fine FFT grid
699 [ - + + - : 6084 : ABI_MALLOC(nhat,(nfftf,dtset%nspden*psps%usepaw))
- + + - +
- + - - +
- + ]
700 [ + + + + : 212988 : if (nstep==0) nhat=zero
+ + ]
701 [ + + + - : 7656 : ABI_MALLOC(pawfgrtab,(my_natom))
+ - + + -
+ - + - +
+ + ]
702 [ + + ]: 1521 : if (my_natom>0) then
703 : : call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,&
704 : 1373 : & mpi_atmtab=mpi_enreg%my_atmtab)
705 : : call pawfgrtab_init(pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat,&
706 : 1373 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
707 [ - + ]: 1373 : ABI_FREE(l_size_atm)
708 : : end if
709 : 1521 : compch_fft=-1.d5
710 [ + + ]: 3654 : usexcnhat=maxval(pawtab(:)%usexcnhat)
711 [ + + - + : 1521 : if (usexcnhat==0.and.dtset%ionmov==4.and.dtset%iscf<10) then
- - ]
712 : 0 : ABI_ERROR('You cannot simultaneously use ionmov=4 and such a PAW psp file !')
713 : : end if
714 : :
715 : : ! Variables/arrays related to the PAW spheres
716 [ + + + - : 7656 : ABI_MALLOC(paw_ij,(my_natom))
+ - + + -
+ - + + +
+ - ]
717 [ + + + - : 7656 : ABI_MALLOC(paw_an,(my_natom))
+ - + + -
+ - + + +
+ - ]
718 : 1521 : call paw_an_nullify(paw_an)
719 : 1521 : call paw_ij_nullify(paw_ij)
720 [ - + ]: 1521 : has_dijhat=0;if (dtset%iscf==22) has_dijhat=1
721 [ + + - + ]: 1521 : has_vhartree=0; if (dtset%prtvha > 0 .or. dtset%prtvclmb > 0) has_vhartree=1
722 [ + + ]: 1521 : has_dijfock=0; if (dtset%usefock == 1) has_dijfock=1
723 [ + + + + : 17817 : has_dijnd=0;if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1
+ + + + ]
724 [ + + ]: 1521 : has_dijU=merge(0,1,dtset%usepawu>0) !Be careful on this!
725 : 1521 : usevxctau_paw=pawxc_get_usekden(dtset%ixc)
726 : : call paw_an_init(paw_an,dtset%natom,dtset%ntypat,0,0,dtset%nspden,&
727 : : & cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1,&
728 : : & has_vxctau=usevxctau_paw,has_vxc_ex=1,has_vhartree=has_vhartree,&
729 : 1521 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
730 : : call paw_ij_init(paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,&
731 : : & dtset%pawspnorb,dtset%natom,dtset%ntypat,dtset%typat,pawtab,&
732 : : & has_dij=1,has_dijfock=has_dijfock,has_dijhartree=1,has_dijnd=has_dijnd,has_dijso=1,&
733 : : & has_dijhat=has_dijhat,&
734 : : & has_dijU=has_dijU,has_pawu_occ=1,has_exexch_pot=1,nucdipmom=dtset%nucdipmom,&
735 : 1521 : & zora=dtset%zora,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
736 [ - + ]: 1521 : if(dtset%usewvl==1) then
737 : 0 : call paw2wvl_ij(1,paw_ij,wvl%descr)
738 : : end if
739 : 1521 : compch_sph=-1.d5
740 [ - + + - : 4563 : ABI_MALLOC(dimcprj,(dtset%natom))
+ - + - -
+ - + -
+ ]
741 [ - + + - : 3042 : ABI_MALLOC(dimcprj_srt,(dtset%natom))
- + - + ]
742 : 1521 : call pawcprj_getdim(dimcprj ,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'R')
743 : 1521 : call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
744 [ + + ]: 3654 : do itypat=1,dtset%ntypat
745 [ + + ]: 3654 : if (pawtab(itypat)%usepawu/=0) lpawumax=max(pawtab(itypat)%lpawu,lpawumax)
746 : : end do
747 [ + + + - ]: 1521 : if (dtset%usedmatpu/=0.and.lpawumax>0) then
748 [ + - - + ]: 15 : if (2*lpawumax+1/=size(dmatpawu,1).or.2*lpawumax+1/=size(dmatpawu,2)) then
749 : 0 : ABI_BUG('Incorrect size for dmatpawu!')
750 : : end if
751 : : end if
752 : :
753 [ + + + - ]: 1521 : if (dtset%cprj_in_memory==1.and.psps%usepaw==1) then
754 : 142 : call xg_nonlop_init_cplex_alldij(xg_nonlop,paw_ij)
755 : : end if
756 : :
757 [ + + ]: 3654 : eijkl_is_sym=.true.
758 : :
759 : : ! Allocation of projected WF (optional)
760 [ + + ]: 1521 : if (usecprj==1) then
761 : 78 : iorder_cprj=0
762 [ + + ]: 78 : if (dtset%usefock==1) then
763 : 26 : ctocprj_choice = 1
764 [ + + ]: 26 : if (dtset%optforces == 1) then
765 : 16 : ctocprj_choice = 2; ! ncpgr = 3
766 : : end if
767 : : ! if (dtset%optstress /= 0) then
768 : : ! ncpgr = 6 ; ctocprj_choice = 3
769 : : ! end if
770 : : end if
771 : :
772 : : #if defined HAVE_BIGDFT
773 : : if (dtset%usewvl==1) then
774 : : mband_cprj=dtset%mband;if (dtset%paral_kgb/=0) mband_cprj=mband_cprj/mpi_enreg%nproc_band
775 : : mcprj_wvl=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
776 : : ABI_MALLOC(wvl%descr%paw%cprj,(dtset%natom,mcprj_wvl))
777 : : call cprj_paw_alloc(wvl%descr%paw%cprj,0,dimcprj_srt)
778 : : end if
779 : : #endif
780 : : end if
781 : :
782 : : ! Other variables for PAW
783 [ + + ]: 1521 : nullify(pawrhoij_ep);if(associated(electronpositron))pawrhoij_ep=>electronpositron%pawrhoij_ep
784 [ + + ]: 1521 : nullify(lmselect_ep);if(associated(electronpositron))lmselect_ep=>electronpositron%lmselect_ep
785 : : else
786 [ - + - + ]: 5418 : ABI_MALLOC(dimcprj,(0))
787 [ - + ]: 5418 : ABI_MALLOC(dimcprj_srt,(0))
788 [ - + ]: 5418 : ABI_MALLOC(nhat,(0,0))
789 [ - + ]: 5418 : ABI_MALLOC(paw_ij,(0))
790 [ - + ]: 5418 : ABI_MALLOC(paw_an,(0))
791 [ - + - + ]: 5418 : ABI_MALLOC(pawfgrtab,(0))
792 : : end if ! PAW
793 : :
794 : : !Several parameters and arrays for the SCF mixing:
795 : : !These arrays are needed only in the self-consistent case
796 [ + + ]: 6939 : if (dtset%iscf>=0) then
797 : 6428 : dielar(1)=dtset%diecut;dielar(2)=dtset%dielng
798 : 6428 : dielar(3)=dtset%diemac;dielar(4)=dtset%diemix
799 : 6428 : dielar(5)=dtset%diegap;dielar(6)=dtset%dielam
800 [ + + ]: 6428 : dielar(7)=dtset%diemix;if (dtset%iscf>=10) dielar(7)=dtset%diemixmag
801 [ - + + - : 25712 : ABI_MALLOC(nvresid,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
802 [ + + ]: 6428 : if (dtset%iscf>=10) then
803 [ - + + - : 5716 : ABI_MALLOC(nvtauresid,(nfftf,dtset%nspden*dtset%usekden))
+ + + - +
- + + - +
- + - + ]
804 : : else
805 [ - + + - : 19996 : ABI_MALLOC(nvtauresid,(nfftf,dtset%nspden*usevxctau))
+ + + - +
- + + - +
- + - + ]
806 : : end if
807 [ + + ]: 6428 : if (nstep==0) then
808 [ + + + + ]: 279698 : nvresid=zero
809 [ - + - - ]: 20 : nvtauresid=zero
810 : : end if
811 [ - + + - : 19284 : ABI_MALLOC(dtn_pc,(3,dtset%natom))
+ - + - -
+ - + ]
812 : : ! The next arrays are needed if iscf==5 and ionmov==4,
813 : : ! but for the time being, they are always allocated
814 [ - + + - : 12856 : ABI_MALLOC(grhf,(3,dtset%natom))
- + - + ]
815 : : ! Additional allocation for mixing within PAW
816 : 6428 : npawmix=0
817 [ + + ]: 6428 : if(psps%usepaw==1) then
818 [ + + ]: 4546 : do iatom=1,my_natom
819 : 3048 : itypat=pawrhoij(iatom)%itypat
820 : 3048 : pawrhoij(iatom)%use_rhoijres=1
821 : 3048 : sz1=pawrhoij(iatom)%cplex_rhoij*pawtab(itypat)%lmn2_size
822 : 3048 : sz2=pawrhoij(iatom)%nspden
823 [ - + + - : 12192 : ABI_MALLOC(pawrhoij(iatom)%rhoijres,(sz1,sz2))
- + + - +
- + - - +
- + - + ]
824 [ + + ]: 7082 : do ispden=1,pawrhoij(iatom)%nspden
825 [ + + ]: 369211 : pawrhoij(iatom)%rhoijres(:,ispden)=zero
826 : : end do
827 [ - + + - : 9144 : ABI_MALLOC(pawrhoij(iatom)%kpawmix,(pawtab(itypat)%lmnmix_sz))
+ - + - -
+ - + -
+ ]
828 : 3048 : pawrhoij(iatom)%lmnmix_sz=pawtab(itypat)%lmnmix_sz
829 [ + - - + : 188085 : pawrhoij(iatom)%kpawmix=pawtab(itypat)%kmix
- - - - +
+ ]
830 : : npawmix=npawmix+pawrhoij(iatom)%nspden*pawtab(itypat)%lmnmix_sz &
831 : 4546 : & *pawrhoij(iatom)%cplex_rhoij*pawrhoij(iatom)%qphase
832 : : end do
833 : : end if
834 [ + - ]: 6428 : if (dtset%iscf > 0) then
835 : 6428 : denpot = ABI_MIXING_POTENTIAL
836 [ + + ]: 6428 : if (dtset%iscf > 10) denpot = ABI_MIXING_DENSITY
837 [ + + + + : 6428 : if (psps%usepaw==1.and.dtset%pawmixdg==0 .and. dtset%usewvl==0) then
+ - ]
838 : 1190 : ispmix=ABI_MIXING_FOURRIER_SPACE;nfftmix=dtset%nfft;ngfftmix(:)=ngfft(:)
839 : : else
840 : 5238 : ispmix=ABI_MIXING_REAL_SPACE;nfftmix=nfftf;ngfftmix(:)=ngfftf(:)
841 : : end if
842 : : !TRangel: added to avoid segfaults with Wavelets
843 [ + - ]: 6428 : nfftmix_per_nfft=0;if(nfftf>0) nfftmix_per_nfft=(1-nfftmix/nfftf)
844 : : call abi_mixing_new(mix, iscf10, denpot, ispmix, nfftmix, dtset%nspden, npawmix, errid, msg, &
845 : 6428 : & dtset%npulayit, pulayhist_storage=dtset%pulayhiststore)
846 [ - + ]: 6428 : if (errid /= AB7_NO_ERROR) then
847 : 0 : ABI_ERROR(msg)
848 : : end if
849 [ + + ]: 6428 : if (dtset%usekden/=0) then
850 [ - + ]: 73 : if (dtset%useria==12345) then ! This is temporary
851 : : call abi_mixing_new(mix_mgga, iscf10, denpot, ispmix, nfftmix, dtset%nspden, 0, errid, msg, &
852 : 0 : & dtset%npulayit, pulayhist_storage=dtset%pulayhiststore)
853 : : else
854 : : call abi_mixing_new(mix_mgga, 0, denpot, ispmix, nfftmix, dtset%nspden, 0, errid, msg, &
855 : 73 : & dtset%npulayit, pulayhist_storage=dtset%pulayhiststore)
856 : : end if
857 [ - + ]: 73 : if (errid /= AB7_NO_ERROR) then
858 : 0 : ABI_ERROR(msg)
859 : : end if
860 : : end if
861 [ - + ]: 6428 : if (dtset%mffmem == 0) then
862 : 0 : call abi_mixing_use_disk_cache(mix, dtfil%fnametmp_fft)
863 [ # # # # ]: 0 : if (dtset%usekden/=0.and.denpot==ABI_MIXING_DENSITY) &
864 : 0 : & call abi_mixing_use_disk_cache(mix, dtfil%fnametmp_fft_mgga)
865 : : end if
866 : : ! else if (dtset%iscf==0.and.dtset%usewvl==1) then
867 : : ! ispmix=ABI_MIXING_REAL_SPACE;nfftmix=nfftf;ngfftmix(:)=ngfftf(:)
868 : : end if
869 : : else
870 [ - + ]: 511 : ABI_MALLOC(nvresid,(0,0))
871 [ - + - + ]: 511 : ABI_MALLOC(nvtauresid,(0,0))
872 [ - + ]: 511 : ABI_MALLOC(dtn_pc,(0,0))
873 [ - + ]: 511 : ABI_MALLOC(grhf,(0,0))
874 : : end if ! iscf>0
875 : :
876 : : ! Initializing precon-object for chi0 based preconditioning
877 : : call precon%init(dtset, atindx, atindx1, cg, cprj, dimcprj, dtfil, eigen, energies%e_fermie, gmet, gprimd, indsym, &
878 : : & irrzon, kg, mcprj, nattyp, nfftmix, ngfftmix, npwarr, occ, pawang, pawfgr, pawfgrtab, &
879 : 6939 : & pawtab, ph1d, phnons, psps, rhor, rmet, rprimd, symrec, ucvol, usecprj, vxc, xred, ylm)
880 : :
881 : : ! Here initialize the datastructure constrained_dft, for constrained DFT calculations
882 : : ! as well as penalty function constrained magnetization
883 [ + + + + : 16454 : if(any(dtset%constraint_kind(:)/=0).or.dtset%magconon/=0)then
+ + + + ]
884 : : call constrained_dft_ini(dtset%chrgat,constrained_dft,dtset%constraint_kind,dtset%magconon,dtset%magcon_lambda,&
885 : : & mpi_enreg,dtset%natom,nfftf,ngfftf,dtset%nspden,dtset%ntypat,&
886 : 59 : & dtset%ratsm,dtset%ratsph,rprimd,dtset%spinat,dtset%typat,xred,dtset%ziontypat,dtset%znucl,dtset%qgbt,dtset%use_gbt,dtset%spinaxis)
887 : : endif
888 : :
889 : : !Here, allocate arrays for computation of susceptibility and dielectric matrix or for TDDFT
890 : :
891 [ + + + + : 6939 : if( (nstep>0 .and. dtset%iscf>=0) .or. dtset%iscf==-1 ) then !MF
+ + ]
892 : :
893 : : ! Here, for TDDFT, artificially set iprcel . Also set a variable to reduce the memory needs.
894 : 6430 : afford=1
895 [ + + ]: 6430 : if(dtset%iscf==-1) then
896 : : ! dtset%iprcel=21
897 : 22 : afford=0
898 : : end if
899 : :
900 : : ! First compute dimensions
901 [ + + + + ]: 6430 : if(dtset%iprcel>=21 .or. dtset%iscf==-1)then
902 : : ! With dielop=1, the matrices will be computed when istep=dielstrt
903 : : ! With dielop=2, the matrices will be computed when istep=dielstrt and 1
904 : 80 : dielop=1
905 [ + + ]: 80 : if(dtset%iprcel>=41)dielop=2
906 [ + + ]: 80 : if((dtset%iprcel >= 71).and.(dtset%iprcel<=79)) dielop=0 !RSkerker preconditioner do not need the susceptibility matrix
907 : : ! Immediate computation of dielectric matrix
908 : 80 : dielstrt=1
909 : : ! Or delayed computation
910 [ + + ]: 80 : if(modulo(dtset%iprcel,100)>21 .and. modulo(dtset%iprcel,100)<=29)dielstrt=modulo(dtset%iprcel,100)-20
911 [ - + ]: 80 : if(modulo(dtset%iprcel,100)>31 .and. modulo(dtset%iprcel,100)<=39)dielstrt=modulo(dtset%iprcel,100)-30
912 [ + + ]: 80 : if(modulo(dtset%iprcel,100)>41 .and. modulo(dtset%iprcel,100)<=49)dielstrt=modulo(dtset%iprcel,100)-40
913 [ + + ]: 80 : if(modulo(dtset%iprcel,100)>51 .and. modulo(dtset%iprcel,100)<=59)dielstrt=modulo(dtset%iprcel,100)-50
914 [ + + ]: 80 : if(modulo(dtset%iprcel,100)>61 .and. modulo(dtset%iprcel,100)<=69)dielstrt=modulo(dtset%iprcel,100)-60
915 : : ! Get diecut, and the fft grid to be used for the susceptibility computation
916 : 80 : diecut=abs(dtset%diecut)
917 [ - + ]: 80 : if( dtset%diecut<0.0_dp )then
918 : 0 : ecutsus=ecut
919 : : else
920 : 80 : ecutsus= ( sqrt(ecut) *0.5_dp + sqrt(diecut) *0.25_dp )**2
921 : : end if
922 : : ! Impose sequential calculation
923 [ + + + + ]: 1040 : ngfftdiel(1:3)=0 ; ngfftdiel(7)=100 ; ngfftdiel(9)=0; ngfftdiel(8)=dtset%ngfft(8);ngfftdiel(10:18)=0
924 [ + + ]: 80 : if(dtset%iscf==-1)ngfftdiel(7)=102
925 : :
926 : : ! The dielectric stuff is performed in sequential mode; set mpi_enreg_diel accordingly
927 : 80 : call initmpi_seq(mpi_enreg_diel)
928 : : call getng(dtset%boxcutmin,dtset%chksymtnons,ecutsus,gmet,k0,mpi_enreg_diel%me_fft,mgfftdiel,nfftdiel,ngfftdiel,&
929 : : & mpi_enreg_diel%nproc_fft,dtset%nsym,mpi_enreg_diel%paral_kgb,dtset%symrel,dtset%tnons,&
930 : 80 : & gpu_option=dtset%gpu_option)
931 : : ! Update the fft distribution
932 : 80 : call mpi_enreg_diel%distribfft%init_seq('c',ngfftdiel(2),ngfftdiel(3),'all')
933 : :
934 : : ! Compute the size of the dielectric matrix
935 : 80 : kpt_diel(1:3)=(/ 0.0_dp, 0.0_dp, 0.0_dp /)
936 : 80 : call getmpw(diecut,dtset%exchn2n3d,gmet,(/1/),kpt_diel,mpi_enreg_diel,npwdiel,1)
937 : 80 : lmax_diel=0
938 [ - + ]: 160 : if (psps%usepaw==1) then
939 [ # # ]: 0 : do ii=1,dtset%ntypat
940 : 0 : lmax_diel=max(lmax_diel,pawtab(ii)%lcut_size)
941 : : end do
942 : : end if
943 : : else
944 : 6350 : npwdiel=1
945 : 6350 : mgfftdiel=1
946 : 6350 : nfftdiel=1
947 : 6350 : lmax_diel=0
948 : 6350 : afford=0
949 : : end if
950 : :
951 : : ! Now, performs allocation
952 [ + + + - : 38580 : ABI_MALLOC(dielinv,(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden))
- + + - -
+ + - - +
+ - + - +
+ - + -
+ ]
953 [ + + + - : 32150 : ABI_MALLOC(susmat,(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden))
- + + - -
+ + - - +
+ - - + -
+ ]
954 [ - + + - : 19290 : ABI_MALLOC(kg_diel,(3,npwdiel))
+ - + - -
+ - + ]
955 [ + - - + : 25720 : ABI_MALLOC(gbound_diel,(2*mgfftdiel+8,2))
+ - + - +
- + - - +
- + ]
956 [ - + + - : 25720 : ABI_MALLOC(irrzondiel,(nfftdiel**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
- + + - +
- + - - +
- + ]
957 [ - + + - : 25720 : ABI_MALLOC(phnonsdiel,(2,nfftdiel**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
- + + - +
- + - - +
- + ]
958 [ + + + + : 19290 : ABI_MALLOC(ph1ddiel,(2,3*(2*mgfftdiel+1)*dtset%natom*psps%usepaw))
+ - + - +
+ - + - +
- + ]
959 [ - + + - : 25720 : ABI_MALLOC(ylmdiel,(npwdiel,lmax_diel**2))
+ - - - +
- - + - +
- + ]
960 : : ! Then, compute the values of different arrays
961 [ + + ]: 6430 : if(dielop>=1)then
962 : : ! Note : npwarr_diel is dummy, npwtot_diel is dummy
963 : : ! This kpgio call for going from the suscep FFT grid to the diel sphere
964 : 78 : npwarr_diel(1)=npwdiel
965 : :
966 : : call kpgio(diecut,dtset%exchn2n3d,gmet,(/1/),kg_diel,&
967 : : & kpt_diel,1,(/1/),1,'COLL',mpi_enreg_diel,npwdiel,&
968 : 78 : & npwarr_diel,npwtot_diel,dtset%nsppol)
969 : 78 : call sphereboundary(gbound_diel,1,kg_diel,mgfftdiel,npwdiel)
970 : :
971 [ + + + + ]: 78 : if (dtset%nsym>1 .and. dtset%iscf>=0 ) then
972 : : ! Should replace this initialization of irrzondiel and phnonsdiel through setsym by a direct call to irrzg
973 [ - + + - : 168 : ABI_MALLOC(indsym_dum,(4,dtset%nsym,dtset%natom))
- + + - +
- + - - +
- + ]
974 [ - + + - : 126 : ABI_MALLOC(symrec_dum,(3,3,dtset%nsym))
+ - - + -
+ ]
975 : : call setsym(indsym_dum,irrzondiel,dtset%iscf,dtset%natom,&
976 : : & nfftdiel,ngfftdiel,dtset%nspden,dtset%nsppol,dtset%nsym,phnonsdiel,&
977 : 42 : & dtset%symafm,symrec_dum,dtset%symrel,dtset%tnons,dtset%typat,xred)
978 [ - + ]: 42 : ABI_FREE(indsym_dum)
979 [ - + ]: 42 : ABI_FREE(symrec_dum)
980 : : end if
981 [ - + ]: 78 : if (psps%usepaw==1) then
982 : : call getph(atindx,dtset%natom,ngfftdiel(1),ngfftdiel(2),&
983 : 0 : & ngfftdiel(3),ph1ddiel,xred)
984 : : call initylmg(gprimd,kg_diel,kpt_diel,1,mpi_enreg_diel,&
985 : : & lmax_diel,npwdiel,dtset%nband,1,npwarr_diel,dtset%nsppol,0,&
986 : 0 : & rprimd,ylmdiel,rhodum)
987 : : end if
988 : : end if
989 : :
990 [ + + + + ]: 6430 : if(dtset%iprcel>=21 .or. dtset%iscf==-1)then
991 : 80 : call destroy_mpi_enreg(mpi_enreg_diel)
992 : : end if
993 : :
994 : : else
995 : 509 : npwdiel=1
996 : 509 : mgfftdiel=1
997 : 509 : nfftdiel=1
998 : 509 : afford = 0
999 [ - + ]: 509 : ABI_MALLOC(susmat,(0,0,0,0,0))
1000 [ - + ]: 509 : ABI_MALLOC(kg_diel,(0,0))
1001 [ - + ]: 509 : ABI_MALLOC(gbound_diel,(0,0))
1002 [ - + ]: 509 : ABI_MALLOC(irrzondiel,(0,0,0))
1003 [ - + ]: 509 : ABI_MALLOC(phnonsdiel,(0,0,0))
1004 [ - + - + ]: 509 : ABI_MALLOC(ph1ddiel,(0,0))
1005 [ - + ]: 509 : ABI_MALLOC(ylmdiel,(0,0))
1006 : : end if
1007 : :
1008 : 6939 : nkxc=0
1009 : : !TDDFT - For a first coding
1010 [ + + + + ]: 6939 : if (dtset%iscf==-1 .and. dtset%nspden==1) nkxc=2
1011 [ + + + + ]: 6939 : if (dtset%iscf==-1 .and. dtset%nspden==2) nkxc=3
1012 : : !Eventually need kxc-LDA when susceptibility matrix has to be computed
1013 [ + + + + : 6939 : if (dtset%iscf>0.and.modulo(dtset%iprcel,100)>=61.and.(dtset%iprcel<71.or.dtset%iprcel>79)) nkxc=2*min(dtset%nspden,2)-1
+ + ]
1014 : : !Eventually need kxc-LDA for residual forces (when density mixing is selected)
1015 : : if (dtset%iscf>=10.and.dtset%usewvl==0.and.forces_needed>0 .and. &
1016 [ + + + - : 6939 : & abs(dtset%densfor_pred)>=1.and.abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5) then
+ + + - +
+ ]
1017 [ + + + - ]: 1327 : if (dtset%xclevel==1.or.dtset%densfor_pred>=0) nkxc=2*min(dtset%nspden,2)-1
1018 [ + + + + : 1327 : if (dtset%xclevel==2.and.dtset%nspden==1.and.dtset%densfor_pred<0) nkxc=7 ! This is not full kxc for mGGA
- + ]
1019 [ + + + + : 1327 : if (dtset%xclevel==2.and.dtset%nspden==2.and.dtset%densfor_pred<0) nkxc=19 ! This is not full kxc for mGGA
- + ]
1020 : : end if
1021 : : !Eventually need Kxc to precondition the SCF.
1022 [ + + ]: 6939 : if (precon%use_kxc) then
1023 : 1 : nkxc = precon%nkxc
1024 : : end if
1025 [ + + ]: 6939 : if (nkxc>0) then
1026 : 1351 : call check_kxc(dtset%ixc,dtset%optdriver)
1027 : : end if
1028 [ - + + - : 27756 : ABI_MALLOC(kxc,(nfftf,nkxc))
+ + + - +
- + + - +
- + - + ]
1029 : 6939 : call precon%init_kxc(kxc)
1030 : :
1031 : : !This flag will be set to 1 just before an eventual change of atomic
1032 : : !positions inside the iteration, and set to zero when the consequences
1033 : : !of this change are taken into account.
1034 : 6939 : moved_atm_inside=0
1035 : : !This flag will be set to 1 if the forces are computed inside the iteration.
1036 : 6939 : computed_forces=0
1037 : :
1038 [ - + ]: 6939 : if(dtset%wfoptalg==2)then
1039 [ # # # # : 0 : ABI_MALLOC(shiftvector,((dtset%mband+2)*dtset%nkpt))
# # # # #
# # # ]
1040 : 6939 : val_min=-1.0_dp
1041 : 6939 : val_max=zero
1042 : : else
1043 [ - + ]: 6939 : ABI_MALLOC(shiftvector,(1))
1044 : : end if
1045 : :
1046 : : !!PAW+DMFT: allocate structured datatype paw_dmft if dtset%usedmft=1
1047 : : !call init_sc_dmft(dtset%dmftbandi,dtset%dmftbandf,dtset%mband,dtset%nkpt,&
1048 : : !& dtset%nsppol,dtset%usedmft,paw_dmft,dtset%usedmft)
1049 : : !call print_sc_dmft(paw_dmft)
1050 : :
1051 : : !!Electric field initializations: initialize pel_cg(:) and p_ion(:)
1052 : : call update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,&
1053 : : & efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,&
1054 : : & dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,nattyp,ngfft,dtset%nkpt,npwarr,&
1055 : : & dtset%ntypat,pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,&
1056 : : & pwind,pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,&
1057 : 6939 : & 0,quit,istep,ucvol,unit_out,psps%usepaw,xred,ylm,ylmgr)
1058 : :
1059 [ - + ]: 6939 : if (dtset%iscf==22) energies%h0=zero
1060 : :
1061 : 6939 : call timab(1441,2,tsec)
1062 : 6939 : call pstat_proc%print(_PSTAT_ARGS_)
1063 : :
1064 : : !##################################################################
1065 : : !PERFORM ELECTRONIC ITERATIONS
1066 : : !##################################################################
1067 : :
1068 : : !Offer option of computing total energy with existing
1069 : : !wavefunctions when nstep<=0, else do nstep iterations
1070 : : !Note that for non-self-consistent calculations, this loop will be exited
1071 : : !after the first call to vtorho
1072 : : !Pass through the first routines even when nstep==0
1073 : :
1074 : 6939 : quitsum_request = xmpi_request_null; timelimit_exit = 0
1075 : 6939 : istep_updatedfock=0
1076 : :
1077 [ - + + - : 272678 : ABI_ICALLOC(rmm_diis_status, (2, dtset%nkpt, dtset%nsppol))
- + + - +
- + - - +
- + + + +
+ + + ]
1078 : :
1079 : : ! start SCF loop
1080 : : ABI_NVTX_START_RANGE(NVTX_SCF)
1081 [ + + ]: 48465 : do istep=1,max(1,nstep)
1082 : :
1083 [ + + ]: 47688 : if(associated(rcpaw)) then
1084 : 21 : rcpaw%istep=istep
1085 [ + + ]: 42 : eijkl_is_sym=rcpaw%eijkl_is_sym
1086 : : endif
1087 : :
1088 : :
1089 : : ! Handle time limit condition.
1090 [ + + ]: 47688 : if (istep == 1) prev = abi_wtime()
1091 [ + + ]: 47688 : if (istep > 1) then
1092 : 40749 : now = abi_wtime()
1093 : 40749 : wtime_step = now - prev
1094 : 40749 : prev = now
1095 : : call wrtout(std_out, sjoin("{SCF_istep:", itoa(istep-1), ", Vnl|psi>:", itoa(nonlop_counter), &
1096 : 40749 : ", wall_time: '", sec2str(wtime_step), "'} <<< TIME"))
1097 : 40749 : nonlop_counter = 0
1098 : :
1099 [ + + ]: 40749 : if (have_timelimit_in(MY_NAME)) then
1100 [ + + ]: 16 : if (istep > 2) then
1101 : 12 : call xmpi_wait(quitsum_request,ierr)
1102 [ - + ]: 12 : if (quitsum_async > 0) then
1103 [ # # ]: 0 : write(msg,"(3a)")"Approaching time limit ",trim(sec2str(get_timelimit())),". Will exit istep loop in scfcv_core."
1104 : 0 : ABI_COMMENT(msg)
1105 : 0 : call wrtout(ab_out, msg)
1106 : 0 : timelimit_exit = 1
1107 : 0 : exit
1108 : : end if
1109 : : end if
1110 : :
1111 [ - + ]: 16 : my_quit = 0; if (now - get_start_time() + 2.15 * wtime_step > get_timelimit()) my_quit = 1
1112 : 16 : call xmpi_isum(my_quit,quitsum_async,spacecomm,quitsum_request,ierr)
1113 : : end if
1114 : : end if
1115 : :
1116 : 47688 : call timab(1442,1,tsec)
1117 [ + + + + ]: 47688 : if (moved_atm_inside==1 .or. istep==1) then
1118 : : ! ##############################################################
1119 : : ! The following steps are done once for a given set of atomic
1120 : : ! coordinates or for the nstep=1 case
1121 : : ! --------------------------------------------------------------
1122 : :
1123 : : ! Eventually symmetrize atomic coordinates over space group elements:
1124 : 7045 : call symmetrize_xred(dtset%natom,dtset%nsym,dtset%symrel,dtset%tnons,xred,indsym=indsym)
1125 : :
1126 [ + - ]: 7045 : if (dtset%usewvl == 0) then
1127 : : ! Get cut-off for g-vectors
1128 [ + + ]: 7045 : if (psps%usepaw==1) call wrtout(std_out,' FFT (fine) grid used in SCF cycle:')
1129 : 7045 : call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf)
1130 : :
1131 : : ! Compute structure factor phases and large sphere cut-off (gsqcut):
1132 : 7045 : call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
1133 : :
1134 [ + + + + ]: 7045 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
1135 : 1452 : call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
1136 : : else
1137 [ + + + + ]: 7174228 : ph1df(:,:)=ph1d(:,:)
1138 : : end if
1139 : : end if
1140 : :
1141 : : ! If using GEMM nonlop, reset precomputed projectors
1142 [ + + ]: 7045 : if(gemm_nonlop_use_gemm) call reset_gemm_nonlop()
1143 : :
1144 : : ! Initialization of atomic data for PAW
1145 [ + + ]: 7045 : if (psps%usepaw==1) then
1146 : :
1147 : : ! Check for non-overlapping spheres, allow one remittance in case of optimization or MD or image algorithms
1148 : 1521 : nremit=0
1149 [ + + + + ]: 1521 : if(dtset%ionmov>0 .and. itimes(1)==1)nremit=1
1150 [ + + + + ]: 1521 : if(dtset%imgmov>0 .and. itimes(2)==1)nremit=mpi_enreg%my_nimage
1151 : 1521 : call chkpawovlp(dtset%natom,psps%ntypat,dtset%pawovlp,pawtab,rmet,dtset%typat,xred,nremit=nremit)
1152 : :
1153 : : ! Identify parts of the rectangular grid where the density has to be calculated
1154 : 1521 : optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm
1155 : : if ((forces_needed==1).or. &
1156 [ + + + + : 1521 : & (dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0).or. &
+ + + + +
+ - + ]
1157 : : & (dtset%positron/=0.and.forces_needed==2)) then
1158 [ + + - + ]: 721 : optgr1=dtset%pawstgylm;if (stress_needed==1) optrad=1; if (dtset%pawprtwf==1) optrad=1
1159 : : end if
1160 : :
1161 [ + - ]: 1521 : if(dtset%usewvl==0) then
1162 : : call nhatgrid(atindx1,gmet,my_natom,dtset%natom,&
1163 : : & nattyp,ngfftf,psps%ntypat,optcut,optgr0,optgr1,optgr2,optrad,&
1164 : : & pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,&
1165 : : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1166 : 1521 : & comm_fft=spaceComm_fft,distribfft=mpi_enreg%distribfft)
1167 : : else
1168 : : shift=0
1169 : : #if defined HAVE_BIGDFT
1170 : : shift=wvl%descr%Glr%d%n1i*wvl%descr%Glr%d%n2i*wvl%den%denspot%dpbox%nscatterarr(me_wvl,4)
1171 : : call wvl_nhatgrid(atindx1,wvl%descr%atoms%astruct%geocode,&
1172 : : & wvl%descr%h,wvl%den%denspot%dpbox%i3s,dtset%natom,dtset%natom,&
1173 : : & nattyp,psps%ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,&
1174 : : & wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,&
1175 : : & wvl%den%denspot%dpbox%n3pi,optcut,optgr0,optgr1,optgr2,optrad,&
1176 : : & pawfgrtab,pawtab,psps%gth_params%psppar,rprimd,shift,xred)
1177 : : #endif
1178 : : end if
1179 : : end if
1180 : :
1181 : : ! If we are inside SCF cycle or inside dynamics over ions,
1182 : : ! we have to translate the density of previous iteration
1183 : 7045 : moved_rhor=0
1184 : :
1185 [ + - + - : 2791 : if (initialized/=0.and.dtset%usewvl == 0.and.ipositron/=1.and. &
+ + + + ]
1186 [ + + ]: 7045 : & (abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)) then
1187 : 2318 : moved_rhor=1
1188 [ + + ]: 2318 : if (abs(dtset%densfor_pred)==2) then
1189 : 1976 : option=2
1190 [ - + + - : 7904 : ABI_MALLOC(workr,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
1191 : : call fresid(dtset,gresid,mpi_enreg,nfftf,ngfftf,&
1192 : : & psps%ntypat,option,pawtab,rhor,rprimd,&
1193 [ - + ]: 1976 : & ucvol,workr,xred,xred_old,psps%znuclpsp)
1194 [ + + + + ]: 37954426 : rhor=workr
1195 [ - + ]: 1976 : ABI_FREE(workr)
1196 [ + - ]: 342 : else if (abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6) then
1197 : 342 : scf_history%icall=scf_history%icall+1
1198 : : call extraprho(atindx,atindx1,cg,cprj,dtset,gmet,gprimd,gsqcut,&
1199 : : & scf_history%icall,kg,mcg,mcprj,mgfftf,mpi_enreg,psps%mqgrid_vl,&
1200 : : & my_natom,nattyp,nfftf,ngfftf,npwarr,psps%ntypat,pawrhoij,pawtab,&
1201 : : & ph1df,psps,psps%qgrid_vl,rhor,rprimd,scf_history,ucvol,&
1202 [ - + ]: 342 : & psps%usepaw,xred,xred_old,ylm,psps%ziontypat,psps%znuclpsp,extfpmd=extfpmd)
1203 : : end if
1204 [ - + - + ]: 2318 : call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0)
1205 : : end if
1206 [ + + + - : 7045 : if (initialized/=0.and.dtset%usewvl == 0.and.ipositron/=1) then
+ - ]
1207 : : ! In some cases cprj are kept in memory, so we have to update them before the call of vtorho
1208 [ - + ]: 2791 : if (dtset%cprj_in_memory==2) then
1209 : 0 : iatom=0
1210 : 0 : idir=0
1211 : 0 : iorder_cprj=0
1212 : 0 : call wrtout(std_out,' Computing cprj from wavefunctions (scfcv_core)')
1213 : : ABI_NVTX_START_RANGE(NVTX_CTOCPRJ)
1214 : : call ctocprj(atindx,cg,1,cprj,gmet,gprimd,iatom,idir,&
1215 : : & iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,&
1216 : : & dtset%mpw,dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft, dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,&
1217 : : & dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,dtfil%unpaw,&
1218 [ # # ]: 0 : & xred,ylm,ylmgr)
1219 : : ABI_NVTX_END_RANGE()
1220 : 0 : call wrtout(std_out,' cprj is computed')
1221 : : end if
1222 : : end if
1223 : :
1224 : : ! if any nuclear dipoles are nonzero, compute the vector potential in real space (depends on
1225 : : ! atomic position so should be done for nstep = 1 and for updated ion positions
1226 [ + + + + : 90033 : if ( any(abs(dtset%nucdipmom(:,:))>tol8) ) then
+ + + + ]
1227 : 28 : with_vectornd = 1
1228 : : else
1229 : 7017 : with_vectornd = 0
1230 : : end if
1231 [ + + ]: 7045 : if(allocated(vectornd)) then
1232 : 106 : ABI_FREE(vectornd)
1233 : : end if
1234 [ + + + - : 35225 : ABI_MALLOC(vectornd,(with_vectornd*nfftf,dtset%nspden,3))
- + + - +
- + - + +
- + - + -
+ ]
1235 [ + + + + : 17773270 : vectornd=zero
+ + ]
1236 [ + + ]: 7045 : if(with_vectornd .EQ. 1) then
1237 : : call make_vectornd(1,gsqcut,psps%usepaw,mpi_enreg,dtset%natom,nfftf,ngfftf,&
1238 : 28 : & dtset%nspden,dtset%nucdipmom,rprimd,vectornd,xred)
1239 : : endif
1240 : :
1241 : : end if ! moved_atm_inside==1 .or. istep==1
1242 : :
1243 : 47688 : call timab(1442,2,tsec)
1244 : :
1245 : : !Initialize/Update data in the case of an Exact-exchange (Hartree-Fock) or hybrid XC calculation
1246 : 47688 : hyb_mixing=zero;hyb_mixing_sr=zero
1247 [ + + ]: 47688 : if (dtset%usefock==1) then
1248 : 1375 : call timab(1443,1,tsec)
1249 [ + + ]: 1375 : if (istep==1) then
1250 : : ! Initialize data_type fock for the calculation
1251 [ + + ]: 123 : cplex_hf=cplex; if (psps%usepaw==1) cplex_hf=dtset%pawcpxocc
1252 : 123 : call fock_init(atindx,cplex_hf,dtset,fock,gsqcut,kg,mpi_enreg,nattyp,npwarr,pawang,pawfgr,pawtab,rprimd)
1253 : :
1254 [ + + ]: 123 : if (fock%fock_common%usepaw==1) then
1255 : 26 : optcut_hf = 0 ! use rpaw to construct local_pawfgrtab
1256 : 26 : optgr0_hf = 0; optgr1_hf = 0; optgr2_hf = 0 ! dont need gY terms locally
1257 : 26 : optrad_hf = 1 ! do store r-R
1258 : : call nhatgrid(atindx1,gmet,dtset%natom,dtset%natom,nattyp,ngfftf,psps%ntypat,&
1259 : : & optcut_hf,optgr0_hf,optgr1_hf,optgr2_hf,optrad_hf,fock%fock_common%pawfgrtab,pawtab,&
1260 : 26 : & rprimd,dtset%typat,ucvol,xred,typord=1)
1261 : 26 : iatom=-1;idir=0
1262 : : call ctocprj(atindx,cg,ctocprj_choice,cprj,gmet,gprimd,iatom,idir,&
1263 : : & iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,&
1264 : : & dtset%mpw,dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft, dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,&
1265 : : & dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,dtfil%unpaw,&
1266 [ - + ]: 26 : & xred,ylm,ylmgr)
1267 : : end if
1268 [ + + ]: 123 : if(wfmixalg/=0)then
1269 : 4 : spare_mem=0
1270 : : if(spare_mem==1)history_size=wfmixalg ! Not yet coded
1271 : 4 : if(spare_mem==0)history_size=2*(wfmixalg-1)+1
1272 : : ! Specific case of simple mixing : always history_size=1
1273 [ + + ]: 4 : if(wfmixalg==2)history_size=1
1274 : 4 : scf_history_wf%history_size=history_size
1275 : 4 : usecg=2
1276 : 4 : call scf_history_init(dtset,mpi_enreg,usecg,scf_history_wf)
1277 : : end if
1278 : : end if
1279 : :
1280 : : !Fock energy
1281 : 1375 : energies%e_exactX=zero
1282 [ + + + + : 3135 : if (fock%fock_common%optfor) fock%fock_common%forces=zero
+ + ]
1283 : :
1284 : 1375 : call timab(1443,2,tsec)
1285 : :
1286 [ + + + + : 1375 : if (istep==1 .or. istep_updatedfock==fock%fock_common%nnsclo_hf .or. &
+ - + + ]
1287 : : & (fock%fock_common%nnsclo_hf>1 .and. fock%fock_common%scf_converged) ) then
1288 : :
1289 : 1238 : istep_updatedfock=1
1290 : 1238 : fock%fock_common%scf_converged=.false.
1291 : :
1292 : 1238 : call timab(1444,1,tsec)
1293 : :
1294 : : !Possibly mix the wavefunctions from different steps before computing the Fock operator
1295 [ + + + + : 1238 : if(wfmixalg/=0 .and. .not. (wfmixalg==2 .and. abs(scf_history_wf%alpha-one)<tol8) )then
+ - ]
1296 : : call wf_mixing(atindx1,cg,cprj,dtset,istep_fock_outer,mcg,mcprj,mpi_enreg,&
1297 [ - + ]: 43 : & nattyp,npwarr,pawtab,scf_history_wf)
1298 : 43 : istep_fock_outer=istep_fock_outer+1
1299 : :
1300 : 43 : call timab(1444,2,tsec)
1301 : :
1302 : : !DEBUG
1303 : : if(.false.)then
1304 : : !Update the density, from the newly mixed cg and cprj.
1305 : : !Be careful: in PAW, rho does not include the compensation density (added later) !
1306 : : tim_mkrho=6
1307 : : if (psps%usepaw==1) then
1308 : : ABI_MALLOC(rhowfg,(2,dtset%nfft))
1309 : : ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden))
1310 : : ! 1-Compute density from WFs
1311 : : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,&
1312 : : & rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
1313 : : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor)
1314 : : ! 2-Compute rhoij
1315 : : call pawmkrhoij(atindx,atindx1,cprj,dimcprj,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,&
1316 : : & mcprj,dtset%mkmem,mpi_enreg,dtset%natom,dtset%nband,dtset%nkpt,dtset%nspden,dtset%nspinor,&
1317 : : & dtset%nsppol,occ,dtset%paral_kgb,paw_dmft,pawrhoij,dtfil%unpaw,dtset%usewvl,dtset%wtk)
1318 : : ! 3-Symetrize rhoij, compute nhat and add it to rhor
1319 : : ! Note pawrhoij_unsym and pawrhoij are the same, which means that pawrhoij
1320 : : ! cannot be distributed over different atomic sites.
1321 : : cplex=1;ipert=0;idir=0;qpt(:)=zero
1322 : : call pawmkrho(1,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
1323 : : & my_natom,dtset%natom,dtset%nspden,dtset%nsym,dtset%ntypat,&
1324 : : & dtset%paral_kgb,pawang,pawfgr,pawfgrtab,dtset%pawprtvol,pawrhoij,pawrhoij,&
1325 : : & pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,symrec,dtset%typat,ucvol,&
1326 : : & dtset%usewvl,xred,pawnhat=nhat,rhog=rhog,extfpmd=extfpmd)
1327 : : ! 2-Take care of kinetic energy density
1328 : : if(dtset%usekden==1)then
1329 : : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,&
1330 : : & rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
1331 : : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,taug,rhowfr,taur)
1332 : : end if
1333 : : ABI_FREE(rhowfg)
1334 : : ABI_FREE(rhowfr)
1335 : : else
1336 : : write(std_out,*)' scfcv_core : recompute the density after the wf mixing '
1337 : : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
1338 : : & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,extfpmd=extfpmd)
1339 : : if(dtset%usekden==1)then
1340 : : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,taug,taur,&
1341 : : & rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
1342 : : end if
1343 : : end if
1344 : : end if ! false
1345 : : !ENDDEBUG
1346 : :
1347 : : end if
1348 : :
1349 : 1238 : call timab(1445,1,tsec)
1350 : :
1351 : : ! Update data relative to the occupied states in fock
1352 [ - + ]: 1238 : call fock_updatecwaveocc(cg,cprj,dtset,fock,indsym,mcg,mcprj,mpi_enreg,nattyp,npwarr,occ,ucvol)
1353 : :
1354 : 1238 : call timab(1445,2,tsec)
1355 : 1238 : call timab(1446,3,tsec)
1356 : :
1357 : : ! Possibly (re)compute the ACE operator
1358 [ + - ]: 1238 : if(fock%fock_common%use_ACE/=0) then
1359 : : call fock2ACE(cg,cprj,fock,dtset%istwfk,kg,dtset%kptns,dtset%mband,mcg,mcprj,dtset%mgfft,&
1360 : : & dtset%mkmem,mpi_enreg,psps%mpsang,&
1361 : : & dtset%mpw,my_natom,dtset%natom,dtset%nband,dtset%nfft,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspden,&
1362 : : & dtset%nspinor,dtset%nsppol,dtset%ntypat,occ,dtset%optforces,paw_ij,pawtab,ph1d,psps,rprimd,&
1363 [ - + ]: 1238 : & dtset%typat,usecprj,dtset%gpu_option,dtset%wtk,xred,ylm)
1364 : 1238 : energies%e_fock0=fock%fock_common%e_fock0
1365 : : end if
1366 : :
1367 : 1238 : call timab(1446,2,tsec)
1368 : :
1369 : : !Should place a test on whether there should be the final exit of the istep loop.
1370 : : !This test should use focktoldfe.
1371 : : !This should update the info in fock%fock_common%fock_converged.
1372 : : !For the time being, fock%fock_common%fock_converged=.false., so the loop end with the maximal value of nstep always,
1373 : : !except when nnsclo_hf==1 (so the Fock operator is always updated), in which case, the usual exit tests (toldfe, tolvrs, etc)
1374 : : !work fine.
1375 : : !if(fock%fock_common%nnsclo_hf==1 .and. fock%fock_common%use_ACE==0)then
1376 [ + + ]: 1238 : if(fock%fock_common%nnsclo_hf==1) fock%fock_common%fock_converged=.TRUE.
1377 : :
1378 : : !Depending on fockoptmix, possibly restart the mixing procedure for the potential
1379 [ + + ]: 1238 : if(mod(dtset%fockoptmix,10)==1) istep_mix=1
1380 : : else
1381 : 137 : istep_updatedfock=istep_updatedfock+1
1382 : : end if
1383 : :
1384 : : !Used locally
1385 : 1375 : hyb_mixing=fock%fock_common%hyb_mixing ; hyb_mixing_sr=fock%fock_common%hyb_mixing_sr
1386 : : end if ! usefock
1387 : :
1388 : 47688 : call timab(1447,1,tsec)
1389 : :
1390 : : ! Initialize/update data in the electron-positron case
1391 [ + + + + : 47688 : if (dtset%positron<0.or.(dtset%positron>0.and.istep==1)) then
+ + ]
1392 : : call setup_positron(atindx,atindx1,cg,cprj,dtefield,dtfil,dtset,ecore,eigen,&
1393 : : & etotal,electronpositron,energies,fock,forces_needed,gred,gmet,gprimd,&
1394 : : & grchempottn,grcondft,grewtn,grvdw,gsqcut,hdr,extfpmd,initialized0,indsym,istep,istep_mix,kg,&
1395 : : & kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfftf,ngrvdw,nhat,&
1396 : : & nkxc,npwarr,nvresid,occ,optres,paw_ij,pawang,pawfgr,pawfgrtab,&
1397 : : & pawrad,pawrhoij,pawtab,ph1df,ph1d,psps,rhog,rhor,rmet,rprimd,&
1398 : : & stress_needed,strscondft,strsxc,symrec,ucvol,usecprj,usevxctau,vhartr,vpsp,vxc,vxctau,&
1399 [ - + - + : 525 : & xccc3d,xcctau3d,xred,ylm,ylmgr,xg_nonlop)
- + ]
1400 : 525 : ipositron=electronpositron_calctype(electronpositron)
1401 : : end if
1402 : :
1403 : 47688 : call timab(1447,2,tsec)
1404 : 47688 : call timab(1448,3,tsec)
1405 : :
1406 : : if ((moved_atm_inside==1 .or. istep==1).or.&
1407 [ + + + + : 47688 : & (dtset%positron<0.and.istep_mix==1).or.&
+ + + + +
+ + + ]
1408 : : & (mod(dtset%fockoptmix,100)==11 .and. istep_updatedfock==1)) then
1409 : : ! PAW only: we sometimes have to compute compensation density
1410 : : ! and eventually add it to density from WFs
1411 : 7131 : nhatgrdim=0
1412 : 7131 : dummy_nhatgr = .False.
1413 : : if (psps%usepaw==1.and.(dtset%positron>=0.or.ipositron/=1) &
1414 [ + + + + : 7131 : & .and.((usexcnhat==0) &
+ + + + +
+ + + + +
+ + + + +
+ + + ]
1415 : : & .or.(dtset%xclevel==2.and.(dtfil%ireadwf/=0.or.dtfil%ireadden/=0.or.initialized/=0)) &
1416 : : & .or.(dtfil%ireadwf/=0.and.dtfil%ireadden==0.and.initialized==0))) then
1417 : 963 : call timab(558,1,tsec)
1418 [ + + ]: 963 : nhatgrdim=0;if (dtset%xclevel==2) nhatgrdim=usexcnhat*dtset%pawnhatxc
1419 : 963 : ider=2*nhatgrdim;izero=0
1420 [ + + ]: 963 : if (nhatgrdim>0) then
1421 [ - + + - : 140 : ABI_MALLOC(nhatgr,(cplex*nfftf,dtset%nspden,3*nhatgrdim))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1422 : : else
1423 [ - + - + ]: 935 : ABI_MALLOC(nhatgr,(0,0,0))
1424 : 935 : dummy_nhatgr = .True.
1425 : : end if
1426 : : call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,&
1427 : : & nfftf,ngfftf,nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,&
1428 : : & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,k0,rprimd,ucvol_local,dtset%usewvl,xred,&
1429 : : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1430 : : & comm_fft=spaceComm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
1431 : : & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl,&
1432 : 963 : & gpu_thread_limit=dtset%gpu_thread_limit)
1433 [ + + + + : 963 : if (dtfil%ireadwf/=0.and.dtfil%ireadden==0.and.initialized==0) then
+ + ]
1434 [ + + + + ]: 7131206 : rhor(:,:)=rhor(:,:)+nhat(:,:)
1435 [ + - ]: 140 : if(dtset%usewvl==0) then
1436 [ - + - + ]: 140 : call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0)
1437 : : end if
1438 : : end if
1439 : 963 : call timab(558,2,tsec)
1440 : : end if
1441 : :
1442 : : ! The following steps have been gathered in the setvtr routine:
1443 : : ! - get Ewald energy and Ewald forces
1444 : : ! - compute local ionic pseudopotential vpsp
1445 : : ! - possibly compute 3D core electron density xccc3d
1446 : : ! - possibly compute 3D core kinetic energy density
1447 : : ! - possibly compute vxc and vhartr
1448 : : ! - set up vtrial
1449 : :
1450 : 7131 : optene = 4 * optres
1451 [ + + ]: 7131 : if(dtset%iscf==-3) optene=4
1452 [ - + ]: 7131 : if (wvlbigdft) optene = 1 ! VH needed for the WF mixing
1453 : :
1454 [ + + ]: 7131 : if (.not.allocated(nhatgr)) then
1455 [ - + + - : 30840 : ABI_MALLOC(nhatgr,(nfftf,dtset%nspden,3*nhatgrdim))
- + + - -
+ + - - -
+ - - + -
+ - + -
+ ]
1456 : 6168 : dummy_nhatgr = .True.
1457 : : end if
1458 : :
1459 : : ! Compute trial potential
1460 : : ABI_NVTX_START_RANGE(NVTX_SCFCV_SETVTR)
1461 : : call setvtr(atindx1,dtset,energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcut,&
1462 : : & istep,kxc,mgfftf,moved_atm_inside,moved_rhor,mpi_enreg,&
1463 : : & nattyp,nfftf,ngfftf,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,psps%ntypat,&
1464 : : & n1xccc,n3xccc,optene,pawang,pawrad,pawrhoij,pawtab,ph1df,psps,rhog,rhor,&
1465 : : & rmet,rprimd,strsxc,ucvol,usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,&
1466 : : & xccc3d,xred,electronpositron=electronpositron,&
1467 [ - + - + : 7131 : & taur=taur,vxc_hybcomp=vxc_hybcomp,vxctau=vxctau,add_tfw=tfw_activated,xcctau3d=xcctau3d)
- + ]
1468 : : ABI_NVTX_END_RANGE()
1469 : :
1470 : : ! set the zero of the potentials here
1471 [ + + + + ]: 34131 : if(dtset%usepotzero==2) vpsp(:) = vpsp(:) + ecore / ( zion * ucvol )
1472 : :
1473 [ + + ]: 7131 : if(dtset%optdriver==RUNL_GWLS) call build_vxc(vxc,nfftf,dtset%nspden)
1474 : :
1475 [ + + - + : 7131 : if ((nhatgrdim>0.and.nstep>0).or.dummy_nhatgr) then
+ - ]
1476 [ - + ]: 7131 : ABI_FREE(nhatgr)
1477 : : end if
1478 : :
1479 : : ! Recursion Initialisation
1480 [ + + + - ]: 7131 : if(dtset%userec==1 .and. istep==1) then
1481 : 11 : rec_set%quitrec = 0
1482 : : ! --At any step calculate the metric
1483 : 11 : call Init_MetricRec(rec_set%inf,rec_set%nl%nlpsp,rmet,ucvol,rprimd,xred,dtset%ngfft(1:3),dtset%natom,rec_set%debug)
1484 : 11 : call rec_set%mpi%distribfft%free()
1485 : 11 : call rec_set%mpi%distribfft%init('c',rec_set%mpi%nproc_fft,rec_set%ngfftrec(2),rec_set%ngfftrec(3))
1486 : 11 : call rec_set%mpi%distribfft%init('f',rec_set%mpi%nproc_fft,dtset%ngfft(2),dtset%ngfft(3))
1487 [ + + ]: 11 : if(initialized==0) call first_rec(dtset,psps,rec_set)
1488 : : end if
1489 : :
1490 : : ! End the condition of atomic position change or istep==1
1491 : : end if
1492 : :
1493 : 47688 : call timab(1448,2,tsec)
1494 : 47688 : call timab(1449,1,tsec)
1495 : :
1496 : : ! ######################################################################
1497 : : ! The following steps are done at every iteration
1498 : : ! ----------------------------------------------------------------------
1499 : : ! PAW: Compute energies and potentials in the augmentation regions (spheres)
1500 : : ! Compute pseudopotential strengths (Dij quantities)
1501 [ + + ]: 47688 : if (psps%usepaw==1)then
1502 : :
1503 : : ! Local exact exch.: impose occ. matrix if required
1504 [ + + ]: 15474 : if (dtset%useexexch/=0) then
1505 : : call setrhoijpbe0(dtset,initialized0,istep,istep_mix,&
1506 : : & spaceComm,my_natom,dtset%natom,dtset%ntypat,pawrhoij,pawtab,dtset%typat,&
1507 : 21 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1508 : : end if
1509 : :
1510 : : ! Computation of on-site densities/potentials/energies
1511 [ + + + + ]: 15474 : nzlmopt=0;if (istep_mix==2.and.dtset%pawnzlm>0) nzlmopt=-1
1512 [ + + ]: 15474 : if (istep_mix>2) nzlmopt=dtset%pawnzlm
1513 : 15474 : call paw_an_reset_flags(paw_an) ! Force the recomputation of on-site potentials
1514 : 15474 : call paw_ij_reset_flags(paw_ij,self_consistent=.true.) ! Force the recomputation of Dij
1515 [ + + + - ]: 15474 : option=0;if (dtset%iscf>0.and.dtset%iscf<10.and.nstep>0) option=1
1516 : : ABI_NVTX_START_RANGE(NVTX_SCFCV_PAWDENPOT)
1517 : : call pawdenpot(compch_sph,el_temp,&
1518 : : & gprimd,ipert,dtset%ixc,my_natom,dtset%natom,dtset%nspden,psps%ntypat,dtset%nucdipmom,nzlmopt,&
1519 : : & option,paw_an,paw_an,energies%paw,paw_ij,pawang,dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,&
1520 : : & pawtab,dtset%pawxcdev,dtset%spnorbscl,dtset%xclevel,&
1521 : : & dtset%xc_denpos,dtset%xc_taupos,xred,ucvol,psps%znuclpsp,dtset%spinaxis,&
1522 : : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1523 : : & hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,&
1524 : 15474 : & electronpositron=electronpositron,vpotzero=vpotzero,rcpaw=rcpaw,extfpmd=extfpmd)
1525 : : ABI_NVTX_END_RANGE()
1526 : :
1527 : : ! Correct the average potential with the calculated constant vpotzero
1528 : : ! Correct the total energies accordingly
1529 : : ! vpotzero(1) = -beta/ucvol
1530 : : ! vpotzero(2) = -1/ucvol sum_ij rho_ij gamma_ij
1531 : : write(msg,'(a,f14.6,2x,f14.6)') &
1532 [ + + + + ]: 92844 : & ' average electrostatic smooth potential [Ha] , [eV]',SUM(vpotzero(:)),SUM(vpotzero(:))*Ha_eV
1533 : 15474 : call wrtout(std_out, msg)
1534 [ + + + + : 821061146 : vtrial(:,:)=vtrial(:,:)+SUM(vpotzero(:))
+ + ]
1535 [ + + ]: 15474 : if(option/=1)then
1536 : : ! Fix the direct total energy (non-zero only for charged systems)
1537 [ + + ]: 43845 : energies%paw%epaw=energies%paw%epaw-SUM(vpotzero(:))*dtset%cellcharge(1)
1538 : : ! Fix the double counting total energy accordingly (for both charged AND
1539 : : ! neutral systems)
1540 [ + + ]: 43845 : energies%paw%epaw_dc=energies%paw%epaw_dc-SUM(vpotzero(:))*zion+vpotzero(2)*dtset%cellcharge(1)
1541 : : end if
1542 : :
1543 : : ! PAW+U: impose density matrix if required
1544 : : ! not available if usepawu<0 (PAW+U without occupation matrix)
1545 [ + + - + ]: 15474 : if (dtset%usepawu>0.and.(ipositron/=1)) then
1546 : 1876 : impose_dmat=0
1547 [ + + + - : 1876 : if ((istep<=abs(dtset%usedmatpu)).and.(dtset%usedmatpu<0.or.initialized0==0)) impose_dmat=1
+ - ]
1548 [ + + + + ]: 1876 : if (impose_dmat==1.or.dtset%dmatudiag/=0) then
1549 [ + + ]: 106 : dimdmat=0;if (impose_dmat==1) dimdmat=2*lpawumax+1
1550 : : call setnoccmmp(0,dimdmat,&
1551 : : & dmatpawu(1:dimdmat,1:dimdmat,1:dtset%nsppol*dtset%nspinor,1:dtset%natpawu*impose_dmat),&
1552 : : & dtset%dmatudiag,impose_dmat,indsym,my_natom,dtset%natom,dtset%natpawu,&
1553 : : & dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,dtset%pawprtvol,&
1554 : : & pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,&
1555 [ + + + + : 490 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
+ + + + +
+ + + +
+ ]
1556 : : ! Reinitialize mixing if PAW+U and occupation matrix now allowed to change
1557 : : ! For experimental purpose...
1558 [ + - - - : 106 : if ((dtset%userib==1234).and.(istep==abs(dtset%usedmatpu)).and. &
- - - - ]
1559 : 0 : & (dtset%usedmatpu<0.or.initialized0==0)) reset_mixing=.true.
1560 : : end if
1561 : : end if
1562 : :
1563 : : ! Write out unperturbed occupancies to dtpawuj-dataset LMac
1564 [ + + + + : 15474 : if (dtset%usepawu/=0.and.dtset%macro_uj>0.and.istep==1.and.ipositron/=1) then
+ + + - ]
1565 : : call pawuj_red(istep, 0, dtfil, dtset,dtpawuj,fatvshift,my_natom,dtset%natom,dtset%ntypat,&
1566 : 10 : paw_ij,pawrad,pawtab,ndtpawuj,spaceComm, comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1567 : : end if
1568 : : ABI_NVTX_START_RANGE(NVTX_SCFCV_DIJ)
1569 : :
1570 : : ! Dij computation
1571 : 15474 : call timab(561,1,tsec)
1572 : :
1573 : : call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,dtset%natom,nfftf,nfftotf,&
1574 : : & dtset%nspden,psps%ntypat,paw_an,paw_ij,pawang,pawfgrtab,dtset%pawprtvol,&
1575 : : & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,k0,dtset%spnorbscl,&
1576 : : & ucvol_local,dtset%cellcharge(1),vtrial,vxc,xred,dtset%znucl,&
1577 : : & natvshift=dtset%natvshift,atvshift=dtset%atvshift,fatvshift=fatvshift,&
1578 : : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1579 : : & mpi_comm_grid=spaceComm_grid,&
1580 : : & hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,&
1581 : : & electronpositron_calctype=ipositron,&
1582 : : & electronpositron_pawrhoij=pawrhoij_ep,&
1583 : : & electronpositron_lmselect=lmselect_ep,&
1584 : 15474 : & nucdipmom=dtset%nucdipmom,eijkl_is_sym=eijkl_is_sym,spinaxis=dtset%spinaxis)
1585 : :
1586 : : ! Symetrize Dij
1587 : : call symdij(gprimd,indsym,ipert,my_natom,dtset%natom,dtset%nsym,&
1588 : : & psps%ntypat,0,paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,&
1589 : 15474 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1590 [ - + ]: 15474 : if (has_dijhat==1) then
1591 : : call symdij(gprimd,indsym,ipert,my_natom,dtset%natom,dtset%nsym,&
1592 : : & psps%ntypat,1,paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,&
1593 : 0 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1594 : : end if
1595 [ - + ]: 15474 : if(dtset%usewvl==1) then
1596 : 0 : call paw2wvl_ij(3,paw_ij,wvl%descr)
1597 : : end if
1598 : :
1599 : 15474 : call timab(561,2,tsec)
1600 : : ABI_NVTX_END_RANGE()
1601 : : end if
1602 : :
1603 : : ! Now that the perturbation has been applied, we harvest occupancies for the perturbed case: LMac
1604 [ + + + + : 47688 : if (dtset%usepawu/=0.and.dtset%macro_uj>0.and.istep>1.and.ipositron/=1) then
+ + + - ]
1605 : : call pawuj_red(istep, 1, dtfil, dtset,dtpawuj,fatvshift,my_natom,dtset%natom,dtset%ntypat,&
1606 : 136 : paw_ij,pawrad,pawtab,ndtpawuj,spaceComm,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1607 : : end if
1608 : :
1609 : 47688 : call timab(1449,2,tsec)
1610 : :
1611 : : ! No need to continue and call vtorho, when nstep==0
1612 [ + + ]: 47688 : if(nstep==0)exit
1613 : :
1614 : : ! ######################################################################
1615 : : ! The following steps are done only when nstep>0
1616 : : ! ----------------------------------------------------------------------
1617 : 47653 : call timab(1450,1,tsec)
1618 : :
1619 [ + + ]: 47653 : if(dtset%iscf>=0)then
1620 : 47157 : write(msg, '(a,a,i4)' )ch10,' ITER STEP NUMBER ',istep
1621 : 47157 : call wrtout(std_out,msg)
1622 : : end if
1623 : :
1624 : : ! The next flag says whether the xred have to be changed in the current iteration
1625 : 47653 : moved_atm_inside=0
1626 : : ! /< Hack to remove iapp from scfcv_core
1627 : : ! for ionmov 4|5 ncycle=1
1628 : : ! Hence iapp = itime
1629 [ + + ]: 47653 : if ( dtset%jdtset /= scfcv_jdtset ) then
1630 : : ! new dtset -> reinitialize
1631 : 3961 : scfcv_jdtset = dtset%jdtset
1632 : 3961 : scfcv_itime = 0
1633 : : end if
1634 [ + + ]: 47653 : if ( istep==1 ) scfcv_itime = scfcv_itime + 1
1635 [ + + + + : 47653 : if(dtset%ionmov==4 .and. mod(scfcv_itime,2)/=1 .and. dtset%iscf>=0 ) moved_atm_inside=1
+ - ]
1636 [ + + + + : 47653 : if(dtset%ionmov==5 .and. scfcv_itime/=1 .and. istep==1 .and. dtset%iscf>=0) moved_atm_inside=1
+ + + - ]
1637 : : ! /< Hack to remove iapp from scfcv_core
1638 : :
1639 : : ! Thomas-Fermi scheme might use a different toldfe criterion
1640 [ + + ]: 47653 : if (dtset%tfkinfunc>0.and.dtset%tfkinfunc/=2) then
1641 [ + + ]: 277 : tollist(4)=dtset%toldfe;if (.not.tfw_activated) tollist(4)=dtset%tfw_toldfe
1642 : : end if
1643 : :
1644 : : ! The next flag says whether the forces have to be computed in the current iteration
1645 : 47653 : computed_forces=0
1646 [ + + - + : 47653 : if ((dtset%optforces==1 .and. dtset%usewvl == 0).or.(moved_atm_inside==1)) computed_forces=1
- + ]
1647 [ + + ]: 47653 : if (abs(tollist(3))>tiny(0._dp)) computed_forces=1
1648 [ + + ]: 47653 : if (dtset%iscf<0) computed_forces=0
1649 [ + + + + ]: 47653 : if ((istep==1).and.(dtset%optforces/=1)) then
1650 [ - + ]: 3856 : if (moved_atm_inside==1) then
1651 : : write(msg,'(5a)')&
1652 : 0 : & 'Although the computation of forces during electronic iterations',ch10,&
1653 : 0 : & 'was not required by user, it is done (required by the',ch10,&
1654 : 0 : & 'choice of ionmov input parameter).'
1655 : 0 : ABI_WARNING(msg)
1656 : : end if
1657 [ - + ]: 3856 : if (abs(tollist(3))+abs(tollist(7))>tiny(0._dp)) then
1658 : : write(msg,'(5a)')&
1659 : 0 : & 'Although the computation of forces during electronic iterations',ch10,&
1660 : 0 : & 'was not required by user, it is done (required by the',ch10,&
1661 : 0 : & '"toldff" or "tolrff" tolerance criteria).'
1662 : 0 : ABI_WARNING(msg)
1663 : : end if
1664 : : end if
1665 [ + + + + : 47653 : if ((istep==1).and.(dtset%optforces==1).and. dtset%usewvl == 1) then
- + ]
1666 : : write(msg,'(5a)')&
1667 : 0 : & 'Although the computation of forces during electronic iterations',ch10,&
1668 : 0 : & 'was required by user, it has been disable since the tolerence',ch10,&
1669 : 0 : & 'is not on forces (force computation is expensive in wavelets).'
1670 : 0 : ABI_WARNING(msg)
1671 : : end if
1672 : :
1673 : 47653 : call timab(1450,2,tsec)
1674 : :
1675 : : ! ######################################################################
1676 : : ! Compute the density rho from the trial potential
1677 : : ! ----------------------------------------------------------------------
1678 : 47653 : call timab(1451,1,tsec)
1679 : : ! Compute the density from the trial potential
1680 [ + + ]: 47653 : if (dtset%tfkinfunc==0) then
1681 : : if(VERBOSE) call wrtout(std_out,'*. Compute the density from the trial potential (vtorho)')
1682 : :
1683 : : ! Dynamic invovl_precision (ChebFI+PAW only): start loose, tighten back
1684 : : ! to the user-requested value as the active SCF criterion (whichever of
1685 : : ! tolwfr/toldfe/tolvrs was specified) approaches convergence.
1686 : : if (dtset%invovl_dyn_prec/=0 .and. psps%usepaw==1 &
1687 [ + + + - : 47338 : & .and. (dtset%wfoptalg==1.or.dtset%wfoptalg==111.or.dtset%wfoptalg==112)) then
- + - - ]
1688 [ + + ]: 80 : if (istep==1) then
1689 : 6 : dtset%invovl_precision=dtset%invovl_prec_ceil
1690 : : else
1691 [ + + ]: 74 : if (abs(dtset%tolwfr)>tiny(zero)) then
1692 : 18 : invovl_dyn_ratio=residm/dtset%tolwfr
1693 [ + + ]: 56 : else if (abs(dtset%toldfe)>tiny(zero)) then
1694 : 28 : invovl_dyn_ratio=abs(deltae)/dtset%toldfe
1695 [ + - ]: 28 : else if (abs(dtset%tolvrs)>tiny(zero)) then
1696 : 28 : invovl_dyn_ratio=res2/dtset%tolvrs
1697 : : else
1698 : 0 : invovl_dyn_ratio=residm/invovl_precision_orig
1699 : : end if
1700 : : dtset%invovl_precision=min(dtset%invovl_prec_ceil, &
1701 : 74 : & max(invovl_precision_orig,invovl_precision_orig*invovl_dyn_ratio*invovl_dyn_prec_safety))
1702 : : end if
1703 : : end if
1704 : :
1705 : : ABI_NVTX_START_RANGE(NVTX_VTORHO)
1706 : :
1707 : : call vtorho(afford,atindx,atindx1,cg,compch_fft,cprj,cpus,dbl_nnsclo,&
1708 : : & dielop,dielstrt,dmatpawu,dphase,dtefield,dtfil,dtset,&
1709 : : & eigen,electronpositron,energies,etotal,gbound_diel,&
1710 : : & gmet,gprimd,grnl,gsqcut,hdr,extfpmd,indsym,irrzon,irrzondiel,&
1711 : : & istep,istep_mix,itimes,kg,kg_diel,kxc,lmax_diel,mcg,mcprj,mgfftdiel,mpi_enreg,&
1712 : : & my_natom,dtset%natom,nattyp,nfftf,nfftdiel,ngfftdiel,nhat,nkxc,&
1713 : : & npwarr,npwdiel,res2,psps%ntypat,nvresid,occ,&
1714 : : & computed_forces,optres,paw_dmft,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,&
1715 : : & pawrhoij,pawtab,phnons,phnonsdiel,ph1d,ph1ddiel,psps,fock,&
1716 : : & pwind,pwind_alloc,pwnsfac,results_gs,resid,residm,rhog,rhor,rmet,rprimd,&
1717 : : & susmat,symrec,taug,taur,nvtauresid,ucvol_local,usecprj,usevxctau,wffnew,with_vectornd,&
1718 [ - + - + : 47338 : & vectornd,vtrial,vxctau,wvl,xg_nonlop,xred,ylm,ylmgr,ylmdiel, rmm_diis_status,rcpaw)
- + - + ]
1719 : :
1720 : : ABI_NVTX_END_RANGE()
1721 : :
1722 [ + + - + ]: 315 : else if (dtset%tfkinfunc==1.or.dtset%tfkinfunc==11.or.dtset%tfkinfunc==12) then
1723 : : ! CP: occopt 9 not available with Thomas Fermi functionals
1724 : 277 : ABI_WARNING('THOMAS FERMI')
1725 : : call vtorhotf(dtset,energies%e_kinetic,energies%e_nlpsp_vfock,&
1726 : : & energies%entropy_ks,energies%e_fermie,gprimd,grnl,irrzon,mpi_enreg,&
1727 : : & dtset%natom,nfftf,dtset%nspden,dtset%nsppol,dtset%nsym,phnons,&
1728 [ - + - + ]: 277 : & rhog,rhor,rprimd,ucvol,vtrial)
1729 : :
1730 : 277 : residm=zero
1731 : 277 : energies%e_eigenvalues=zero
1732 : : end if
1733 : :
1734 : : ! Recursion method
1735 [ + + ]: 47653 : if(dtset%userec==1)then
1736 : : call vtorhorec(dtset,&
1737 : : & energies%e_kinetic,energies%e_nlpsp_vfock,energies%entropy_ks,energies%e_eigenvalues,&
1738 : : & energies%e_fermie,grnl,initialized,irrzon,nfftf,phnons,&
1739 : 38 : & rhog,rhor,vtrial,rec_set,istep-nstep,rprimd,gprimd)
1740 : 38 : residm=zero
1741 : : end if
1742 : :
1743 : : ! Update Fermi level in energies
1744 : 47653 : results_gs%fermie = energies%e_fermie
1745 : :
1746 [ - + ]: 47653 : if(dtset%wfoptalg==2)then
1747 [ # # ]: 0 : do ikpt=1,dtset%nkpt
1748 : 0 : shiftvector(1+(ikpt-1)*(dtset%mband+2))=val_min
1749 : : shiftvector(2+(ikpt-1)*(dtset%mband+2):ikpt*(dtset%mband+2)-1)=&
1750 [ # # ]: 0 : & eigen((ikpt-1)*dtset%mband+1:ikpt*dtset%mband)
1751 : 0 : shiftvector(ikpt*(dtset%mband+2))=val_max
1752 : : end do
1753 : : end if
1754 : :
1755 : 47653 : call timab(1451,2,tsec)
1756 : :
1757 : : ! ######################################################################
1758 : : ! Skip out of step loop if non-SCF (completed)
1759 : : ! ----------------------------------------------------------------------
1760 : :
1761 : : ! Indeed, nstep loops have been done inside vtorho
1762 [ + + ]: 47653 : if (dtset%iscf<0) exit
1763 : :
1764 : : ! ######################################################################
1765 : : ! In case of density mixing or wavelet handling, compute the total energy
1766 : : ! ----------------------------------------------------------------------
1767 : :
1768 : 47157 : call timab(1452,1,tsec)
1769 [ + + - + ]: 47157 : if (dtset%iscf>=10 .or. wvlbigdft) then
1770 : 14674 : optene = 1 ! use double counting scheme (default)
1771 [ - + - - ]: 14674 : if (wvlbigdft.and.dtset%iscf==0) optene = 0 ! use direct scheme
1772 [ - + ]: 14674 : if (dtset%iscf==22) optene = -1
1773 : :
1774 : : ! Add the Fock contribution to E_xc and E_xcdc if required
1775 [ + + ]: 14674 : if (dtset%usefock==1) then
1776 : 188 : energies%e_fockdc=two*energies%e_fock
1777 : : end if
1778 : :
1779 : : ! if the mixing is the ODA mixing, compute energy and new density here
1780 [ - + ]: 14674 : if (dtset%iscf==22) then
1781 : : call odamix(deltae,dtset,&
1782 : : & elast,energies,etotal,gprimd,gsqcut,kxc,mpi_enreg,&
1783 : : & my_natom,nfftf,ngfftf,nhat,nkxc,psps%ntypat,nvresid,n3xccc,optres,&
1784 : : & paw_ij,paw_an,pawang,pawfgrtab,pawrad,pawrhoij,pawtab,&
1785 : : & red_ptot,psps,rhog,rhor,rprimd,strsxc,ucvol,psps%usepaw,&
1786 : : & usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,xccc3d,xred,&
1787 [ # # # # : 0 : & taur=taur,vxctau=vxctau,add_tfw=tfw_activated)
# # ]
1788 : : end if
1789 : : ! If the density mixing is required, compute the total energy here
1790 : : ! TODO: add nvtauresid if needed (for forces?)
1791 : : call etotfor(atindx1,deltae,diffor,dtefield,dtset,&
1792 : : & elast,electronpositron,energies,&
1793 : : & etotal,favg,fcart,fock,forold,gred,gmet,grchempottn,grcondft,gresid,grewtn,grhf,grnl,grvdw,&
1794 : : & grxc,gsqcut,extfpmd,indsym,kxc,maxfor,mgfftf,mpi_enreg,my_natom,&
1795 : : & nattyp,nfftf,ngfftf,ngrvdw,nhat,nkxc,psps%ntypat,nvresid,n1xccc,n3xccc,&
1796 : : & optene,computed_forces,optres,pawang,pawfgrtab,pawrad,pawrhoij,pawtab,&
1797 : : & ph1df,red_ptot,psps,rhog,rhor,rmet,rprimd,symrec,synlgr,ucvol,&
1798 [ - + - + ]: 14674 : & psps%usepaw,usevxctau,vhartr,vpsp,vxc,vxctau,wvl%descr,wvl%den,xccc3d,xred,rcpaw)
1799 : : !if (wvlbigdft) energies%copy(energies_wvl) ! TO BE ACTIVATED LATER
1800 : : end if
1801 : 47157 : call timab(1452,2,tsec)
1802 : :
1803 : : ! ######################################################################
1804 : : ! In case of density mixing, check the exit criterion
1805 : : ! ----------------------------------------------------------------------
1806 [ + + - + : 47157 : if (dtset%iscf>=10.or.(wvlbigdft.and.dtset%iscf>0)) then
- - ]
1807 : : ! Check exit criteria
1808 : 14674 : call timab(1453,1,tsec)
1809 : 14674 : choice=2
1810 [ + + ]: 14674 : if(paw_dmft%use_dmft==1) then
1811 : 113 : call prtene(dtset,energies,std_out,psps%usepaw)
1812 : : end if
1813 [ + - + + : 14674 : if(response==0.and.(dtset%iscf>0.or.dtset%iscf==-3).and.(dtset%nsppol==2.or.dtset%nspden>1)) then
+ + ]
1814 [ - + + + : 27404 : intgden0=intgden
+ + ]
1815 : : call calcdenmagsph(mpi_enreg,dtset%natom,nfftf,ngfftf,dtset%nspden,&
1816 [ - + ]: 2886 : dtset%ntypat,dtset%ratsm,dtset%ratsph,rhor,rprimd,dtset%typat,xred,1,cplex1,dtset%qgbt,dtset%use_gbt,intgden=intgden)
1817 : : !Compute maximal magnet and maximal difference of magnet
1818 : 2886 : call calmaxdifmag(cplex1,intgden,intgden0,dtset%natom,dtset%nspden,maxmag,difmag)
1819 [ + + - + : 2886 : if (dtset%prt_lorbmag==1 .and. (dtset%nspinor==2) .and. (dtset%nspden==4 ) .and. (dtset%usepawu .ne. 0)) then
- + ]
1820 : : call loc_orbmom_cal(1,0,dmatdum,0,0,indsym,my_natom,dtset%natom,dtset%natpawu,&
1821 : : & dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,pawrad,dtset%pawprtvol,&
1822 : : & pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,dtset%znucl,&
1823 : 8 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,orb_mom_atom=orb_mom_atom,maxmag=maxmag,difmag=difmag)
1824 : : endif
1825 : : endif
1826 : : call scprqt(choice,cpus,deltae,diffor,maxmag,difmag,dtset,&
1827 : : & eigen,etotal,favg,fcart,energies%e_fermie,energies%e_fermih,dtfil%fnameabo_app_eig,&
1828 : : & dtfil%filnam_ds(1),initialized0,dtset%iscf,istep,istep_fock_outer,istep_mix,dtset%kptns,&
1829 : : & maxfor,moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,nstep,&
1830 : : & occ,optres,prtfor,prtxml,quit,res2,resid,residm,response,tollist,&
1831 : : & psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode,&
1832 : 14674 : & electronpositron=electronpositron,fock=fock,rcpaw=rcpaw)
1833 : 14674 : call timab(1453,2,tsec)
1834 : :
1835 : : ! Check if we need to exit the loop
1836 : 14674 : call timab(1454,1,tsec)
1837 [ - + - - : 14674 : if (dtset%tfkinfunc>10.and.(.not.tfw_activated).and.quit==1) then
- - ]
1838 : 0 : quit=0;tfw_activated=.true.;reset_mixing=.true.
1839 : : end if
1840 [ - + - - ]: 14674 : if(dtset%userec==1.and.rec_set%quitrec==2)quit=1
1841 [ + + ]: 14674 : if (istep==nstep) quit=1
1842 : 14674 : quit_sum=quit
1843 : 14674 : call xmpi_sum(quit_sum,spaceComm,ierr)
1844 [ + + ]: 14674 : if (quit_sum>0) quit=1
1845 : 14674 : call timab(1454,2,tsec)
1846 : :
1847 : : ! If criteria in scprqt say to quit, then exit the loop over istep.
1848 [ + + ]: 14674 : if (quit==1) exit
1849 : : end if
1850 : :
1851 : :
1852 : : ! ######################################################################
1853 : : ! Mix the total density (if required)
1854 : : ! ----------------------------------------------------------------------
1855 : 45729 : call timab(1455,1,tsec)
1856 : :
1857 [ + + + - ]: 45729 : if (dtset%iscf>=10 .and.dtset%iscf/=22.and. .not. wvlbigdft ) then
1858 : :
1859 : : ! If LDA dielectric matrix is used for preconditionning, has to update here Kxc
1860 : : if (nkxc>0.and.modulo(dtset%iprcel,100)>=61.and.(dtset%iprcel<71.or.dtset%iprcel>79) &
1861 [ + + - + : 13246 : & .and.((istep==1.or.istep==dielstrt).or.(dtset%iprcel>=100))) then
- - - - -
- - - ]
1862 : 0 : optxc=10
1863 : 0 : call xcdata_init(xcdata,dtset=dtset)
1864 : : ! to be adjusted for the call to rhotoxc
1865 : 0 : nk3xc=1
1866 [ # # # # ]: 0 : if(dtset%icoulomb==0 .and. dtset%usewvl==0) then
1867 [ # # # # ]: 0 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
1868 : : call rhotoxc(edum,sdum,kxc,mpi_enreg,nfftf,&
1869 : : & ngfftf,nhat,psps%usepaw,nhatgr,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
1870 : : & optxc,rhor,rprimd,0,vxc,vxcavg_dum,xccc3d,xcdata,&
1871 [ # # ]: 0 : & add_tfw=tfw_activated,taur=taur,vhartr=vhartr,vxctau=vxctau,xcctau3d=xcctau3d)
1872 [ # # ]: 0 : else if(.not. wvlbigdft) then
1873 : : ! WVL case:
1874 : : call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, &
1875 : : & dtset%icoulomb, dtset%ixc, &
1876 : : & mpi_enreg, nfftf, ngfftf,&
1877 : : & nhat,psps%usepaw,&
1878 : : & dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, &
1879 : : & usexcnhat,psps%usepaw,dtset%usewvl,vhartr, vxc, vxcavg,&
1880 : : & wvl%descr,wvl%den,&
1881 [ # # ]: 0 : & wvl%e,xccc3d,dtset%xclevel,dtset%xc_denpos)
1882 : : end if
1883 : : end if
1884 : :
1885 [ + + + - : 13246 : if(associated(extfpmd).and..not.moved_rhor==1.and.dtset%extfpmd_pawsph==1) then
+ + ]
1886 [ + + ]: 6 : if(istep==1) then
1887 : 1 : dielar(3)=one
1888 : 1 : dielar(4)=one
1889 : 1 : dielar(7)=one
1890 [ + + ]: 5 : else if(istep==2) then
1891 : 1 : istep_mix=1
1892 : 1 : dielar(3)=dtset%diemac
1893 : 1 : dielar(4)=dtset%diemix
1894 : 1 : dielar(7)=dtset%diemix
1895 : : endif
1896 : : endif
1897 : :
1898 [ + + ]: 13246 : if(associated(rcpaw)) then
1899 [ + - + + : 18 : if(istep>=rcpaw%updatepaw(1).and.istep<=rcpaw%updatepaw(2)+1.and.rcpaw%updatepaw(2)>0) then
+ + ]
1900 : 4 : istep_mix=1
1901 : 4 : dielar(3)=one
1902 : 4 : dielar(4)=one
1903 : 4 : dielar(7)=one
1904 : : else
1905 : 14 : dielar(3)=dtset%diemac
1906 : 14 : dielar(4)=dtset%diemix
1907 : 14 : dielar(7)=dtset%diemix
1908 : : endif
1909 : : endif
1910 : :
1911 : : ABI_NVTX_START_RANGE(NVTX_SCFCV_NEWRHO)
1912 : : call newrho(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,dtn_pc,&
1913 : : & dtset,etotal,fcart,pawfgr%fintocoa,&
1914 : : & gmet,grhf,gsqcut,initialized,ispmix,istep_mix,kg_diel,kxc,&
1915 : : & mgfftf,mix,pawfgr%coatofin,moved_atm_inside,mpi_enreg,my_natom,nattyp,nfftf,&
1916 : : & nfftmix,nfftmix_per_nfft,ngfftf,ngfftmix,nkxc,npawmix,npwdiel,nvresid,psps%ntypat,&
1917 : : & n1xccc,pawrhoij,pawtab,ph1df,precon,psps,rhog,rhor,&
1918 : : & rprimd,susmat,psps%usepaw,vtrial,wvl%descr,wvl%den,xred,rcpaw,extfpmd,&
1919 [ + - - - : 13246 : & mix_mgga=mix_mgga,taug=taug,taur=taur,tauresid=nvtauresid)
- - + - -
- - - - +
- + - + -
+ - + -
+ ]
1920 : : ABI_NVTX_END_RANGE()
1921 : :
1922 : : end if ! iscf>=10
1923 : :
1924 : 45729 : call timab(1455,2,tsec)
1925 : :
1926 : : ! ######################################################################
1927 : : ! Additional computation in case of an electric field or electric displacement field
1928 : : ! ----------------------------------------------------------------------
1929 : :
1930 : 45729 : call timab(1456,1,tsec)
1931 : :
1932 : : call update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,&
1933 : : & efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,&
1934 : : & dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,nattyp,ngfft,dtset%nkpt,npwarr,&
1935 : : & dtset%ntypat,pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,&
1936 : : & pwind,pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,&
1937 : 45729 : & 1,quit,istep,ucvol,unit_out,psps%usepaw,xred,ylm,ylmgr)
1938 : :
1939 : 45729 : call timab(1456,2,tsec)
1940 : :
1941 : : ! ######################################################################
1942 : : ! Compute the new potential from the trial density
1943 : : ! ----------------------------------------------------------------------
1944 : :
1945 : 45729 : call timab(1457,1,tsec)
1946 : : if(VERBOSE) call wrtout(std_out,'*. Compute the new potential from the trial density')
1947 : :
1948 : : ! Set XC computation flag
1949 : 45729 : optxc=1
1950 [ + + ]: 45729 : if (nkxc>0) then
1951 : : ! MJV 2017 May 25: you should not be able to get here with iscf < 0
1952 [ - + ]: 12381 : if (dtset%iscf<0) optxc=2
1953 : : if (modulo(dtset%iprcel,100)>=61.and.(dtset%iprcel<71.or.dtset%iprcel>79).and. &
1954 [ + + + - : 12381 : & dtset%iscf<10.and. &
+ - + - +
+ + + ]
1955 : 4 : & (dtset%iprcel>=100.or.istep==1.or.istep==dielstrt)) optxc=2
1956 [ + + ]: 12381 : if (precon%use_kxc) optxc=2 ! Kxc needed for (chi0-based) preconditioning.
1957 [ + + + - : 12381 : if (dtset%iscf>=10.and.dtset%densfor_pred/=0.and.abs(dtset%densfor_pred)/=5) optxc=2
+ - ]
1958 [ + + + + : 12381 : if (optxc==2.and.dtset%xclevel==2.and.nkxc==2*min(dtset%nspden,2)-1) optxc=12
+ + ]
1959 : : end if
1960 : :
1961 [ + - ]: 45729 : if (dtset%iscf/=22) then
1962 : : ! PAW: eventually recompute compensation density (and gradients)
1963 : 45729 : nhatgrdim=0
1964 [ - + ]: 45729 : if ( allocated(nhatgr) ) then
1965 : 0 : ABI_FREE(nhatgr)
1966 : : end if
1967 [ + + ]: 45729 : if (psps%usepaw==1) then
1968 [ + + + + : 14030 : ider=-1;if (dtset%iscf>=10.and.((dtset%xclevel==2.and.dtset%pawnhatxc>0).or.usexcnhat==0)) ider=0
+ + + + ]
1969 [ + + + + : 14030 : if (dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0) ider=ider+2
+ + ]
1970 [ + + ]: 14030 : if (ipositron==1) ider=-1
1971 [ + + ]: 14030 : if (ider>0) then
1972 : 741 : nhatgrdim=1
1973 [ - + + - : 3705 : ABI_MALLOC(nhatgr,(nfftf,dtset%nspden,3))
- + + - +
- + - + -
- + - + -
+ ]
1974 : : else
1975 [ - + - + ]: 13289 : ABI_MALLOC(nhatgr,(0,0,0))
1976 : : end if
1977 [ + + ]: 14030 : if (ider>=0) then
1978 : : ABI_NVTX_START_RANGE(NVTX_SCFCV_PAWKNHAT)
1979 : 9090 : call timab(558,1,tsec)
1980 : 9090 : izero=0
1981 : :
1982 : : call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,nfftf,ngfftf,&
1983 : : & nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,nhatgr,nhat,&
1984 : : & pawrhoij,pawrhoij,pawtab,k0,rprimd,ucvol_local,dtset%usewvl,xred,&
1985 : : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1986 : : & comm_fft=spaceComm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
1987 : 9090 : & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
1988 : :
1989 : 9090 : call timab(558,2,tsec)
1990 : : ABI_NVTX_END_RANGE()
1991 : : end if
1992 : : else
1993 [ - + - + ]: 31699 : ABI_MALLOC(nhatgr,(0,0,0))
1994 : : end if
1995 : :
1996 [ + + ]: 45729 : if(associated(rcpaw)) then
1997 [ + + + - ]: 18 : nzlmopt=0;if (istep_mix==2.and.dtset%pawnzlm>0) nzlmopt=-1
1998 [ + + ]: 18 : if (istep_mix>2) nzlmopt=dtset%pawnzlm
1999 [ - + - - ]: 18 : option=0;if (dtset%iscf>0.and.dtset%iscf<10.and.nstep>0) option=1
2000 : : call paw_relax_core(pawtab,pawrad,pawang,pawrhoij,dtset%ntypat,rcpaw,psps,dtset,&
2001 : : & cplex,nzlmopt,option,ucvol_local,paw_an,my_natom,&
2002 : 18 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,extfpmd=extfpmd)
2003 [ + - + + : 18 : if(istep>=rcpaw%updatepaw(1).and.istep<=rcpaw%updatepaw(2).and.dtset%cprj_in_memory==1)then
- + ]
2004 : 0 : call xg_nonlop_destroy_Sij(xg_nonlop)
2005 : 0 : call xg_nonlop_make_Sij(xg_nonlop,pawtab,inv_sij=dtset%wfoptalg==111)
2006 : : endif
2007 [ + - - + : 18 : if(.not.rcpaw%all_atoms_relaxed.and.any(rcpaw%atm(:)%zcore_orig>0)) then
+ - + - ]
2008 : 18 : optn=n3xccc/nfftf
2009 [ - + - - ]: 18 : if(rcpaw%istep>rcpaw%updatetnc.and.rcpaw%updatetnc>0) optn=0
2010 : : call atm2fft(atindx1,xccc3d,vpsp,dummy01,dummy02,dummy03,dummy04,&
2011 : : & gmet,gprimd,dummy05,dummy06,gsqcut,mgfftf,psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,psps%ntypat,&
2012 : : & 1,0,0,0,optn,1,0,1,psps,pawtab,ph1df,psps%qgrid_vl,dtset%qprtrb,&
2013 : : & dtset%rcut,dummy07,rprimd,dummy61,dummy62,ucvol,psps%usepaw,dummy08,dummy09,dummy010,dtset%vprtrb,psps%vlspl,&
2014 : : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
2015 : 18 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
2016 : 18 : call pspcor(ecore,psps%epsatm,dtset%natom,dtset%ntypat,dtset%typat,psps%ziontypat)
2017 [ + - ]: 18 : if(dtset%usepotzero<2) then
2018 : 18 : energies%e_corepsp = ecore / ucvol
2019 : : endif
2020 : : call ewald(energies%e_ewald,gmet,grewtn,gsqcut,dtset%icutcoul,dtset%natom,ngfftf,dtset%nkpt,dtset%ntypat,&
2021 : 18 : & dtset%rcut,rmet,rprimd,dtset%typat,ucvol,dtset%vcutgeo,xred,psps%ziontypat)
2022 : : endif
2023 : : endif
2024 : :
2025 : : ! Compute new potential from the trial density
2026 [ + + ]: 45729 : optene=2*optres;if(psps%usepaw==1) optene=2
2027 : : ABI_NVTX_START_RANGE(NVTX_SCFCV_RHOTOV)
2028 : : call rhotov(constrained_dft,dtset,energies,gprimd,grcondft,gsqcut,intgres,istep,&
2029 : : & kxc,mpi_enreg,nfftf,ngfftf,nhat,nhatgr,nhatgrdim,nkxc,nvresid,n3xccc,&
2030 : : & optene,optres,optxc,pawang,pawrad,pawrhoij,pawtab,&
2031 : : & rhog,rhor,rprimd,strscondft,strsxc,taur,ucvol_local,psps%usepaw,usexcnhat,usevxctau,&
2032 : : & vhartr,vnew_mean,vpsp,vres_mean,res2,vtrial,vxcavg,vxc,vxctau,nvtauresid,wvl,xccc3d,xcctau3d,xred,&
2033 [ - + - + : 45729 : & electronpositron=electronpositron,vxc_hybcomp=vxc_hybcomp,add_tfw=tfw_activated)
- + ]
2034 : : ABI_NVTX_END_RANGE()
2035 : :
2036 [ + + ]: 45729 : if(associated(rcpaw)) then
2037 [ + - ]: 18 : if(.not.rcpaw%all_atoms_relaxed) then
2038 : : call rcpaw_core_eig(pawtab,pawrad,dtset%ntypat,rcpaw,dtset,&
2039 : : & nfftf,vhartr+vpsp,cplex,ucvol_local,&
2040 : : & gmet,rprimd,xred,ngfftf,my_natom,&
2041 : : & distribfft=mpi_enreg%distribfft,comm_fft=spaceComm_fft,&
2042 [ + - + + ]: 1316370 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
2043 : : endif
2044 : : endif
2045 : :
2046 : : end if
2047 : :
2048 : 45729 : call timab(1457,2,tsec)
2049 : 45729 : call timab(1452,1,tsec)
2050 : :
2051 : : ! This is inside the loop, its not equivalent to the line 1821
2052 [ + + + + : 47153 : if(moved_atm_inside==1) xred_old(:,:)=xred(:,:)
+ + ]
2053 : :
2054 [ + + ]: 45729 : if (dtset%iscf<10) then
2055 : :
2056 : : if(VERBOSE) call wrtout(std_out,'Check exit criteria in case of potential mixing')
2057 : :
2058 : : ! If the potential mixing is required, compute the total energy here
2059 : : ! PAW: has to compute here spherical terms
2060 [ + + ]: 32483 : if (psps%usepaw==1) then
2061 [ + + + - ]: 859 : nzlmopt=0;if (istep_mix==1.and.dtset%pawnzlm>0) nzlmopt=-1
2062 [ + + ]: 859 : if (istep_mix>1) nzlmopt=dtset%pawnzlm
2063 : 859 : call paw_an_reset_flags(paw_an) ! Force the recomputation of on-site potentials
2064 : 859 : option=2
2065 : : ABI_NVTX_START_RANGE(NVTX_SCFCV_PAWDENPOT)
2066 : : call pawdenpot(compch_sph,el_temp,gprimd,&
2067 : : & ipert,dtset%ixc,my_natom,dtset%natom,dtset%nspden,&
2068 : : & psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,&
2069 : : & energies%paw,paw_ij,pawang,dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,&
2070 : : & pawtab,dtset%pawxcdev,dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,&
2071 : : & dtset%xc_taupos,xred,ucvol,psps%znuclpsp,dtset%spinaxis,&
2072 : : & hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,comm_atom=mpi_enreg%comm_atom,&
2073 : 859 : & mpi_atmtab=mpi_enreg%my_atmtab,electronpositron=electronpositron)
2074 : : ABI_NVTX_END_RANGE()
2075 : : end if
2076 : :
2077 : : ! Add the Fock contribution to E_xc and E_xcdc if required
2078 [ + + ]: 32483 : if (dtset%usefock==1) energies%e_fockdc=two*energies%e_fock
2079 : :
2080 : : ABI_NVTX_START_RANGE(NVTX_SCFCV_ETOTFOR)
2081 [ + - ]: 32483 : if (.not.wvlbigdft) then
2082 : : ! TODO: add nvtauresid if needed (for forces?)
2083 : : call etotfor(atindx1,deltae,diffor,dtefield,dtset,&
2084 : : & elast,electronpositron,energies,&
2085 : : & etotal,favg,fcart,fock,forold,gred,gmet,grchempottn,grcondft,gresid,grewtn,grhf,grnl,grvdw,&
2086 : : & grxc,gsqcut,extfpmd,indsym,kxc,maxfor,mgfftf,mpi_enreg,my_natom,&
2087 : : & nattyp,nfftf,ngfftf,ngrvdw,nhat,nkxc,dtset%ntypat,nvresid,n1xccc, &
2088 : : & n3xccc,0,computed_forces,optres,pawang,pawfgrtab,pawrad,pawrhoij,&
2089 : : & pawtab,ph1df,red_ptot,psps,rhog,rhor,rmet,rprimd,symrec,synlgr,ucvol,&
2090 [ - + - + ]: 32483 : & psps%usepaw,usevxctau,vhartr,vpsp,vxc,vxctau,wvl%descr,wvl%den,xccc3d,xred,rcpaw)
2091 : : !& xcctau3d=xcctau3d)
2092 : : end if
2093 : : ABI_NVTX_END_RANGE()
2094 : :
2095 : : end if
2096 : 45729 : call timab(1452,2,tsec)
2097 : :
2098 : : ! ######################################################################
2099 : : ! Check exit criteria in case of potential mixing or direct minimization
2100 : : ! ----------------------------------------------------------------------
2101 [ + + - + : 45729 : if ((dtset%iscf<10.and.(.not.wvlbigdft)) .or. dtset%iscf == 0) then
- + ]
2102 : : ! Check exit criteria
2103 : 32483 : call timab(1453,1,tsec)
2104 : 32483 : choice=2
2105 : : ABI_NVTX_START_RANGE(NVTX_SCFCV_SCPRQT)
2106 [ + - + + : 32483 : if(response==0.and.(dtset%iscf>0.or.dtset%iscf==-3).and.(dtset%nsppol==2.or.dtset%nspden>1)) then
+ + ]
2107 [ - + + + : 44246 : intgden0=intgden
+ + ]
2108 : : call calcdenmagsph(mpi_enreg,dtset%natom,nfftf,ngfftf,dtset%nspden,&
2109 [ - + ]: 4648 : dtset%ntypat,dtset%ratsm,dtset%ratsph,rhor,rprimd,dtset%typat,xred,1,cplex1,dtset%qgbt,dtset%use_gbt,intgden=intgden)
2110 : : !Compute maximal magnet and maximal difference of magnet
2111 : 4648 : call calmaxdifmag(cplex1,intgden,intgden0,dtset%natom,dtset%nspden,maxmag,difmag)
2112 [ + - - - : 4648 : if (dtset%prt_lorbmag==1 .and. (dtset%nspinor==2) .and. (dtset%nspden==4 ) .and. (dtset%usepawu .ne. 0)) then
- - ]
2113 : : call loc_orbmom_cal(1,0,dmatdum,0,0,indsym,my_natom,dtset%natom,dtset%natpawu,&
2114 : : & dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,pawrad,dtset%pawprtvol,&
2115 : : & pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,dtset%znucl,&
2116 : 0 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,orb_mom_atom=orb_mom_atom,maxmag=maxmag,difmag=difmag)
2117 : : endif
2118 : : endif
2119 : :
2120 : : call scprqt(choice,cpus,deltae,diffor,maxmag,difmag,dtset,&
2121 : : & eigen,etotal,favg,fcart,energies%e_fermie,energies%e_fermih,dtfil%fnameabo_app_eig,&
2122 : : & dtfil%filnam_ds(1),initialized0,dtset%iscf,istep,istep_fock_outer,istep_mix,dtset%kptns,&
2123 : : & maxfor,moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,nstep,&
2124 : : & occ,optres,prtfor,prtxml,quit,res2,resid,residm,response,tollist,&
2125 : : & psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode,&
2126 : 32483 : & electronpositron=electronpositron,fock=fock)
2127 : : ABI_NVTX_END_RANGE()
2128 : 32483 : call timab(1453,2,tsec)
2129 : :
2130 : : ! Check if we need to exit the loop
2131 : 32483 : call timab(1454,1,tsec)
2132 [ + + + + : 32483 : if (dtset%tfkinfunc>10.and.(.not.tfw_activated).and.quit==1) then
+ + ]
2133 : 1 : quit=0;tfw_activated=.true.;reset_mixing=.true.
2134 : : end if
2135 [ + + + + ]: 32483 : if (istep==nstep.and.psps%usepaw==1) quit=1
2136 [ + + + + ]: 32483 : if(dtset%userec==1 .and. rec_set%quitrec==2) quit=1
2137 : 32483 : quit_sum=quit
2138 : 32483 : call xmpi_sum(quit_sum,spaceComm,ierr)
2139 [ + + ]: 32483 : if (quit_sum > 0) quit=1
2140 : :
2141 : : ! If criteria in scprqt say to quit, then exit the loop over istep.
2142 [ + + ]: 32483 : if (quit==1) then
2143 [ + + ]: 8876 : do ispden=1,dtset%nspden
2144 [ + + ]: 65320293 : vtrial(:,ispden)=vtrial(:,ispden)+nvresid(:,ispden)+vres_mean(ispden)
2145 : : end do
2146 : : !if (usevxctau==1) then
2147 : : ! do ispden=1,dtset%nspden
2148 : : ! vxctau(:,ispden,1)=vxctau(:,ispden,1)+nvtauresid(:,ispden)
2149 : : ! end do
2150 : : !end if
2151 : 4203 : call timab(1454,2,tsec) ! Due to the exit instruction, two timab calls are needed
2152 : 4203 : exit ! exit the loop over istep
2153 : : end if
2154 : 28280 : call timab(1454,2,tsec) ! Due to the exit instruction, two timab calls are needed
2155 : : end if
2156 : :
2157 : : ! ######################################################################
2158 : : ! Mix the potential (if required) - Check exit criteria
2159 : : ! ----------------------------------------------------------------------
2160 : :
2161 : 41526 : call timab(1458,1,tsec)
2162 [ + + + - ]: 41526 : if (dtset%iscf<10 .and. dtset%iscf>0 .and. .not. wvlbigdft) then
2163 : :
2164 : : if(VERBOSE) call wrtout(std_out,'*. Mix the potential (if required) - Check exit criteria')
2165 : :
2166 : : ! Precondition the residual and forces, then determine the new vtrial
2167 : : ! (Warning: the (H)xc potential may have been subtracted from vtrial)
2168 : :
2169 : : call newvtr(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,&
2170 : : & dtn_pc,dtset,etotal,fcart,pawfgr%fintocoa,&
2171 : : & gmet,grhf,gsqcut,initialized,ispmix,&
2172 : : & istep_mix,kg_diel,kxc,mgfftf,mix,pawfgr%coatofin,&
2173 : : & moved_atm_inside,mpi_enreg,my_natom,nattyp,nfftf,nfftmix,&
2174 : : & ngfftf,ngfftmix,nkxc,npawmix,npwdiel,&
2175 : : & nstep,psps%ntypat,n1xccc,&
2176 : : & pawrhoij,ph1df,precon,psps,rhor,rprimd,susmat,psps%usepaw,&
2177 : : & vhartr,vnew_mean,vpsp,nvresid,vres_mean,vtrial,vxc,xred,&
2178 : : & nfftf,pawtab,rhog,wvl,&
2179 [ + - - - : 28280 : & mix_mgga=mix_mgga,vtau=vxctau,vtauresid=nvtauresid)
- - + - -
- - - - +
- + - + -
+ ]
2180 : :
2181 : : end if ! iscf<10
2182 : :
2183 : : ! ######################################################################
2184 : : ! END MINIMIZATION ITERATIONS
2185 : : ! ######################################################################
2186 : :
2187 : : if(VERBOSE) call wrtout(std_out,'*. END MINIMIZATION ITERATIONS')
2188 : :
2189 : : ! The initialisation of the gstate run should be done when this point is reached
2190 : 41526 : initialized=1
2191 : :
2192 : : ! This is to save the density for restart.
2193 [ + + ]: 41526 : if (iwrite_fftdatar(mpi_enreg)) then
2194 : :
2195 [ + + - + ]: 33205 : if(dtset%prtden<0.or.dtset%prtkden<0) then
2196 : : ! Update the content of the header (evolving variables)
2197 : : ! Don't use parallelism over atoms because only me=0 accesses here
2198 : 36 : bantot=hdr%bantot
2199 [ + - ]: 36 : if (dtset%positron==0) then
2200 : : call hdr%update(bantot,etotal,energies%e_fermie,energies%e_fermih,residm,&
2201 : 36 : & rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1))
2202 : : else
2203 : : call hdr%update(bantot,electronpositron%e0,energies%e_fermie,energies%e_fermih,residm,&
2204 : 0 : & rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1))
2205 : : end if
2206 : : end if
2207 : :
2208 [ + + ]: 33205 : if (dtset%prtden<0) then
2209 [ + - ]: 36 : if (mod(istep-1,abs(dtset%prtden))==0) then
2210 : 36 : isave_den=isave_den+1
2211 : 36 : rdwrpaw=0
2212 : 36 : call int2char4(mod(isave_den,2),tag)
2213 [ - + ]: 36 : ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!')
2214 [ + - ]: 36 : fildata=trim(dtfil%fnametmp_app_den)//'_'//trim(tag)
2215 [ + - ]: 36 : if (dtset%iomode == IO_MODE_ETSF) fildata = nctk_ncify(fildata)
2216 : : call fftdatar_write_from_hdr("density",fildata,dtset%iomode,hdr,ngfftf,cplex1,nfftf,&
2217 [ - + ]: 36 : & dtset%nspden,rhor,mpi_enreg,eigen=eigen)
2218 : : end if
2219 : : end if
2220 : :
2221 [ - + ]: 33205 : if (dtset%prtkden<0) then
2222 [ # # ]: 0 : if (mod(istep-1,abs(dtset%prtkden))==0) then
2223 : 0 : isave_kden=isave_kden+1
2224 : 0 : rdwrpaw=0
2225 : 0 : call int2char4(mod(isave_kden,2),tag)
2226 [ # # ]: 0 : ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!')
2227 [ # # ]: 0 : fildata=trim(dtfil%fnametmp_app_kden)//'_'//trim(tag)
2228 [ # # ]: 0 : if (dtset%iomode == IO_MODE_ETSF) fildata = nctk_ncify(fildata)
2229 : : ! output the Laplacian of density
2230 : : call fftdatar_write_from_hdr("kinedr",fildata,dtset%iomode,hdr,ngfftf,cplex1,nfftf,&
2231 [ # # ]: 0 : & dtset%nspden,taur,mpi_enreg,eigen=eigen)
2232 : : end if
2233 : : end if
2234 : :
2235 : : end if
2236 : :
2237 [ - + ]: 41526 : ABI_FREE(nhatgr)
2238 : :
2239 : 41526 : istep_mix=istep_mix+1
2240 [ + + ]: 41526 : if (reset_mixing) then
2241 : 1 : istep_mix=1;reset_mixing=.false.
2242 : : end if
2243 [ + + ]: 41526 : if (ipositron/=0) electronpositron%istep_scf=electronpositron%istep_scf+1
2244 : :
2245 : 44262 : call timab(1458,2,tsec)
2246 : : end do ! istep
2247 : : ABI_NVTX_END_RANGE()
2248 : :
2249 : : ! Restore the user-requested invovl_precision (dynamic schedule is local to the SCF loop).
2250 : 6939 : dtset%invovl_precision=invovl_precision_orig
2251 : :
2252 [ - + ]: 6939 : ABI_FREE(rmm_diis_status)
2253 [ + + ]: 6939 : ABI_SFREE(nhatgr)
2254 : :
2255 : : if((dtset%wfoptalg == 1 .or. dtset%wfoptalg == 111 .or. dtset%wfoptalg == 112) .and. psps%usepaw == 1 &
2256 [ + + + + : 6939 : & .and. dtset%cprj_in_memory==0) then
+ + + + ]
2257 : 96 : call destroy_invovl(dtset%nkpt,dtset%gpu_option)
2258 : : end if
2259 : :
2260 : : ! Avoid pending requests if itime == ntime.
2261 : 6939 : call xmpi_wait(quitsum_request,ierr)
2262 [ - + ]: 6939 : if (timelimit_exit == 1) istep = istep - 1
2263 : :
2264 : 6939 : call timab(1459,1,tsec)
2265 : :
2266 [ + + ]: 6939 : if (dtset%iscf > 0) then
2267 : 6428 : call abi_mixing_deallocate(mix)
2268 [ + + ]: 6428 : if (dtset%usekden/=0) call abi_mixing_deallocate(mix_mgga)
2269 : : end if
2270 : :
2271 [ + + ]: 6939 : if (dtset%usefock==1)then
2272 [ + + ]: 123 : if(wfmixalg/=0) call scf_history_free(scf_history_wf)
2273 : : end if
2274 : :
2275 [ + + + + ]: 6939 : if (quit==1.and.nstep==1) initialized=1
2276 : :
2277 : : !######################################################################
2278 : : !Case nstep==0: compute energy based on incoming wf
2279 : : !----------------------------------------------------------------------
2280 : :
2281 [ + + ]: 6939 : if(nstep==0) then
2282 : 35 : optene=2*psps%usepaw+optres
2283 : 35 : energies%entropy=results_gs%energies%entropy !MT20070219: entropy is not recomputed in routine energy
2284 [ + - ]: 35 : if (.not.allocated(nhatgr) ) then
2285 [ - + - + ]: 35 : ABI_MALLOC(nhatgr,(0,0,0))
2286 : : end if
2287 : :
2288 : : call energy(cg,compch_fft,constrained_dft,dtset,electronpositron,&
2289 : : & energies,eigen,etotal,gsqcut,extfpmd,indsym,irrzon,kg,mcg,mpi_enreg,my_natom,&
2290 : : & nfftf,ngfftf,nhat,nhatgr,nhatgrdim,npwarr,n3xccc,&
2291 : : & occ,optene,paw_dmft,paw_ij,pawang,pawfgr,pawfgrtab,pawrhoij,pawtab,&
2292 : : & phnons,ph1d,psps,resid,rhog,rhor,rprimd,strsxc,symrec,taug,taur,usexcnhat,&
2293 : : & vhartr,vtrial,vpsp,vxc,wvl%wfs,wvl%descr,wvl%den,wvl%e,xccc3d,xred,ylm,&
2294 [ - + - + : 35 : & add_tfw=tfw_activated,vxctau=vxctau,xcctau3d=xcctau3d,rcpaw=rcpaw)
- + - + ]
2295 : :
2296 [ - + ]: 35 : if (nhatgrdim>0) then
2297 [ # # ]: 0 : ABI_FREE(nhatgr)
2298 : : end if
2299 : :
2300 : : end if ! nstep==0
2301 : :
2302 : : !######################################################################
2303 : : !Additional steps after SC iterations, including force, stress, polarization calculation
2304 : : !----------------------------------------------------------------------
2305 : :
2306 [ + + ]: 6939 : if (dtset%userec==1) then
2307 : 11 : call prtene(dtset,energies,ab_out,psps%usepaw)
2308 : 11 : call prtene(dtset,energies,std_out,psps%usepaw)
2309 : : end if
2310 : :
2311 : : !if (wvlbigdft) call energies_wvl%copy(energies) ! TO BE ACTIVATED LATER
2312 : :
2313 : : !PAW: if cprj=<p_lmn|Cnk> are in memory,
2314 : : !need to reorder them (from atom-sorted to unsorted)
2315 [ + + + + ]: 6939 : if (psps%usepaw==1.and.usecprj==1) then
2316 : 78 : iorder_cprj=1
2317 : 78 : call pawcprj_reorder(cprj,atindx1)
2318 [ + + ]: 78 : if (dtset%positron/=0) then
2319 [ + - ]: 3 : if (electronpositron%dimcprj>0) then
2320 : 3 : call pawcprj_reorder(electronpositron%cprj_ep,atindx1)
2321 : : end if
2322 : : end if
2323 [ - + ]: 78 : if (dtset%usewvl==1) then
2324 : 0 : call wvl_cprjreorder(wvl%descr,atindx1)
2325 : : end if
2326 : : end if
2327 : :
2328 : : !PAW: if cprj=<p_lmn|Cnk> are not in memory,need to compute them in some cases
2329 : : recompute_cprj = psps%usepaw ==1 .and. usecprj==0 .and. &
2330 : : & (dtset%prtwant ==2 .or. &
2331 : : & dtset%prtwant ==3 .or. &
2332 : : & dtset%prtnabla > 0 .or. &
2333 : : & dtset%prtdos ==3 .or. &
2334 : : & dtset%prtdos ==4 .or. &
2335 : : & dtset%berryopt /=0 .or. &
2336 : : & dtset%kssform ==3 .or. &
2337 : : & dtset%pawfatbnd> 0 .or. &
2338 : : & dtset%pawprtwf > 0 .or. &
2339 : : & dtset%plowan_compute > 0 .or. &
2340 [ + + + + : 6939 : & dtset%userid .EQ. 1 )
+ + - + -
+ - + + +
+ + - + +
+ - + -
+ ]
2341 : :
2342 [ + + + + : 88623 : if( (psps%usepaw.EQ.1) .AND. ANY(ABS(dtset%nucdipmom)>tol8) ) recompute_cprj=.TRUE.
+ + + + +
+ ]
2343 [ + + + + ]: 6939 : if(dtset%berryopt == -2 .AND. dtset%orbmag /= 0) recompute_cprj=.TRUE.
2344 : :
2345 [ + + ]: 6939 : if (recompute_cprj) then
2346 : 66 : usecprj=1
2347 : 66 : mband_cprj=dtset%mband/mpi_enreg%nproc_band
2348 : 66 : mcprj=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
2349 [ - + + - : 42316 : ABI_MALLOC(cprj_local,(dtset%natom,mcprj))
- + + - +
- + - - +
- + - + +
+ + + ]
2350 : 66 : ncpgr = 0 ; ctocprj_choice = 1
2351 [ + + ]: 66 : if (finite_efield_flag) then
2352 [ + - - + ]: 9 : if (forces_needed /= 0 .and. stress_needed == 0) then
2353 : 0 : ncpgr = 3 ; ctocprj_choice = 2
2354 [ + - + - ]: 9 : else if (forces_needed /= 0 .and. stress_needed /= 0) then
2355 : 9 : ncpgr = 9 ; ctocprj_choice = 23
2356 [ # # # # ]: 0 : else if (forces_needed == 0 .and. stress_needed /= 0) then
2357 : 0 : ncpgr = 6 ; ctocprj_choice = 3
2358 : : end if
2359 : : end if
2360 [ + + + - ]: 66 : if (dtset%berryopt == -2 .AND. dtset%orbmag /= 0) then
2361 : 1 : ncpgr=3; ctocprj_choice=5
2362 : : end if
2363 : 66 : call pawcprj_alloc(cprj_local,ncpgr,dimcprj)
2364 : 66 : cprj=> cprj_local
2365 : 66 : iatom=0 ; iorder_cprj=1 ! cprj are not ordered
2366 : : call ctocprj(atindx,cg,ctocprj_choice,cprj_local,gmet,gprimd,&
2367 : : & iatom,idir,iorder_cprj,dtset%istwfk,kg,dtset%kptns,&
2368 : : & mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,&
2369 : : & dtset%mpw,dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,&
2370 : : & dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,dtset%nsppol,&
2371 : : & dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,&
2372 : 66 : & dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr)
2373 : : end if
2374 : :
2375 : : ! MRM print final Hartree energy components
2376 : 6939 : enonlocalpsp=energies%e_nlpsp_vfock-2.0d0*energies%e_fock0
2377 : : esum=energies%e_kinetic+energies%e_ewald+energies%e_corepsp+energies%e_hartree+energies%e_xc&
2378 : : &+energies%e_localpsp+enonlocalpsp+energies%e_fock0&
2379 : : &+energies%e_hybcomp_E0+energies%e_hybcomp_v0+energies%e_hybcomp_v+energies%e_vdw_dftd&
2380 : 6939 : &+energies%e_elecfield+energies%e_magfield+energies%paw%epaw
2381 : :
2382 [ + + ]: 6939 : if (me == 0) then
2383 : 5462 : write(std_out,'(a1)')' '
2384 : 5462 : write(std_out,'(a98)')'-------------------------------------------------------------------------------------------------'
2385 : 5462 : write(std_out,'(a,2(es16.6,a))')' Ekinetic = : ',energies%e_kinetic ,' Ha ,',energies%e_kinetic*Ha_eV ,' eV'
2386 : 5462 : write(std_out,'(a,2(es16.6,a))')' Evext_l = : ',energies%e_localpsp ,' Ha ,',energies%e_localpsp*Ha_eV ,' eV'
2387 [ + + ]: 5462 : if(enonlocalpsp>1.0d-6) then
2388 : 2211 : write(std_out,'(a,2(es16.6,a))')' Evext_nl = : ',enonlocalpsp ,' Ha ,',enonlocalpsp*Ha_eV ,' eV'
2389 : : endif
2390 : 5462 : write(std_out,'(a,2(es16.6,a))')' Epsp_core = : ',energies%e_corepsp ,' Ha ,',energies%e_corepsp*Ha_eV ,' eV'
2391 : 5462 : write(std_out,'(a,2(es16.6,a))')' Ehartree = : ',energies%e_hartree ,' Ha ,',energies%e_hartree*Ha_eV ,' eV'
2392 [ + + ]: 5462 : if(dtset%usefock==1) then
2393 : 107 : write(std_out,'(a,2(es16.6,a))')' Efock = : ',energies%e_fock0 ,' Ha ,',energies%e_fock0*Ha_eV ,' eV'
2394 : : endif
2395 : 5462 : write(std_out,'(a,2(es16.6,a))')' Exc_ks = : ',energies%e_xc ,' Ha ,',energies%e_xc*Ha_eV ,' eV'
2396 [ + + ]: 5462 : if(abs(energies%e_vdw_dftd)>1.0d-6) then
2397 : 18 : write(std_out,'(a,2(es16.6,a))')' EvdW-D = : ',energies%e_vdw_dftd ,' Ha ,',energies%e_vdw_dftd*Ha_eV ,' eV'
2398 : : endif
2399 [ + + ]: 5462 : if(abs(energies%e_elecfield)>1.0d-6) then
2400 : 17 : write(std_out,'(a,2(es16.6,a))')' Eefield = : ',energies%e_elecfield,' Ha ,',energies%e_elecfield*Ha_eV ,' eV'
2401 : : endif
2402 [ - + ]: 5462 : if(abs(energies%e_magfield)>1.0d-6) then
2403 : 0 : write(std_out,'(a,2(es16.6,a))')' Emfield = : ',energies%e_magfield ,' Ha ,',energies%e_magfield*Ha_eV ,' eV'
2404 : : endif
2405 [ + + ]: 5462 : if(abs(energies%paw%epaw)>1.0d-6) then
2406 : 1049 : write(std_out,'(a,2(es16.6,a))')' Epaw = : ',energies%paw%epaw ,' Ha ,',energies%paw%epaw*Ha_eV ,' eV'
2407 : : endif
2408 : 5462 : write(std_out,'(a,2(es16.6,a))')' Enn = : ',energies%e_ewald ,' Ha ,',energies%e_ewald*Ha_eV ,' eV'
2409 : 5462 : write(std_out,'(a98)')'-------------------------------------------------------------------------------------------------'
2410 : 5462 : write(std_out,'(a,2(es16.6,a))')' Etot = : ',esum ,' Ha ,',esum*Ha_eV ,' eV'
2411 : 5462 : write(std_out,'(a98)')'-------------------------------------------------------------------------------------------------'
2412 : : end if ! end MRM printing energy components
2413 : :
2414 : 6939 : call timab(1459,2,tsec)
2415 : 6939 : call timab(1460,1,tsec)
2416 : :
2417 : :
2418 : : !SHOULD CLEAN THE ARGS OF THIS ROUTINE
2419 : : call afterscfloop(atindx,atindx1,cg,computed_forces,cprj,cpus,&
2420 : : & deltae,diffor,difmag,dtefield,dtfil,dtset,eigen,electronpositron,elfr,&
2421 : : & energies,etotal,extfpmd,favg,fcart,fock,forold,grchempottn,grcondft,&
2422 : : & gred,gresid,grewtn,grhf,grhor,grvdw,&
2423 : : & grxc,gsqcut,hdr,indsym,intgres,irrzon,istep,istep_fock_outer,istep_mix,&
2424 : : & kg,kxc,lrhor,maxfor,maxmag,mcg,mcprj,mgfftf,&
2425 : : & moved_atm_inside,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfft,ngfftf,ngrvdw,nhat,&
2426 : : & nkxc,npwarr,nvresid,occ,optres,paw_an,paw_ij,pawang,pawfgr,&
2427 : : & pawfgrtab,pawrad,pawrhoij,pawtab,pel,pel_cg,ph1d,ph1df,phnons,pion,prtfor,&
2428 : : & prtxml,psps,pwind,pwind_alloc,pwnsfac,res2,resid,residm,results_gs,&
2429 : : & rhog,rhor,rprimd,stress_needed,strscondft,strsxc,strten,symrec,synlgr,taug,&
2430 : : & taur,tollist,usecprj,usevxctau,vhartr,vpsp,vtrial,vxc,vxctau,vxcavg,wvl,&
2431 [ + + - + : 20817 : & xccc3d,xcctau3d,xred,ylm,ylmgr,dtset%cellcharge(1)*SUM(vpotzero(:)),conv_retcode,xg_nonlop)
- + - + ]
2432 : :
2433 : : !Before leaving the present routine, save the current value of xred.
2434 [ + + + + ]: 88435 : xred_old(:,:)=xred(:,:)
2435 : :
2436 : 6939 : call timab(1460,2,tsec)
2437 : :
2438 : : !######################################################################
2439 : : !All calculations in scfcv_core are finished. Printing section
2440 : : !----------------------------------------------------------------------
2441 : :
2442 : 6939 : call timab(1461,1,tsec)
2443 : :
2444 : : call outscfcv(atindx1,cg,compch_fft,compch_sph,cprj,dimcprj,dmatpawu,dtfil,&
2445 : : & dtset,ecut,eigen,electronpositron,elfr,etotal,&
2446 : : & gmet,gprimd,grhor,hdr,intgres,kg,lrhor,dtset%mband,mcg,mcprj,dtset%mgfft,&
2447 : : & dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,my_natom,dtset%natom,nattyp,&
2448 : : & nfftf,ngfftf,nhat,dtset%nkpt,npwarr,dtset%nspden,&
2449 : : & dtset%nsppol,dtset%nsym,psps%ntypat,n3xccc,occ,paw_dmft,pawang,pawfgr,pawfgrtab,&
2450 : : & pawrad,pawrhoij,pawtab,paw_an,paw_ij,dtset%prtvol,psps,results_gs,&
2451 [ - + - + ]: 6939 : & rhor,rprimd,taur,ucvol,usecprj,vhartr,vpsp,vtrial,vxc,wvl%den,xccc3d,xred,rcpaw)
2452 : :
2453 : 6939 : call timab(1461,2,tsec)
2454 : 6939 : call timab(1462,1,tsec)
2455 : :
2456 : : !Transfer eigenvalues and occupation computed by BigDFT in afterscfloop to eigen.
2457 : : #if defined HAVE_BIGDFT
2458 : : if (dtset%usewvl == 1) then
2459 : : if (dtset%nsppol == 1) then
2460 : : eigen = wvl%wfs%ks%orbs%eval
2461 : : occ = wvl%wfs%ks%orbs%occup
2462 : : else
2463 : : eigen(1:wvl%wfs%ks%orbs%norbu) = wvl%wfs%ks%orbs%eval(1:wvl%wfs%ks%orbs%norbu)
2464 : : eigen(dtset%mband + 1:dtset%mband + wvl%wfs%ks%orbs%norbd) = &
2465 : : & wvl%wfs%ks%orbs%eval(wvl%wfs%ks%orbs%norbu + 1:wvl%wfs%ks%orbs%norb)
2466 : : occ(1:wvl%wfs%ks%orbs%norbu) = wvl%wfs%ks%orbs%occup(1:wvl%wfs%ks%orbs%norbu)
2467 : : occ(dtset%mband + 1:dtset%mband + wvl%wfs%ks%orbs%norbd) = &
2468 : : & wvl%wfs%ks%orbs%occup(wvl%wfs%ks%orbs%norbu + 1:wvl%wfs%ks%orbs%norb)
2469 : : end if
2470 : : end if
2471 : : #endif
2472 : : !need to reorder cprj (from unsorted to atom-sorted)
2473 [ + + + + ]: 6939 : if (psps%usepaw==1.and.usecprj==1) then
2474 : 144 : iorder_cprj=0
2475 : 144 : call pawcprj_reorder(cprj,atindx)
2476 [ + + ]: 144 : if (dtset%positron/=0) then
2477 [ + - ]: 3 : if (electronpositron%dimcprj>0) then
2478 : 3 : call pawcprj_reorder(electronpositron%cprj_ep,atindx)
2479 : : end if
2480 : : end if
2481 : : end if
2482 : : !######################################################################
2483 : : !Deallocate memory and save results
2484 : : !----------------------------------------------------------------------
2485 : :
2486 : 6939 : call prc_mem_free()
2487 : :
2488 [ - + ]: 6939 : ABI_FREE(fcart)
2489 [ - + ]: 6939 : ABI_FREE(gred)
2490 [ - + ]: 6939 : ABI_FREE(forold)
2491 [ - + ]: 6939 : ABI_FREE(grchempottn)
2492 [ - + ]: 6939 : ABI_FREE(grcondft)
2493 [ - + ]: 6939 : ABI_FREE(gresid)
2494 [ - + ]: 6939 : ABI_FREE(grewtn)
2495 [ - + ]: 6939 : ABI_FREE(grnl)
2496 [ - + ]: 6939 : ABI_FREE(grvdw)
2497 [ - + ]: 6939 : ABI_FREE(grxc)
2498 [ - + ]: 6939 : ABI_FREE(intgres)
2499 [ - + ]: 6939 : ABI_FREE(synlgr)
2500 [ - + ]: 6939 : ABI_FREE(ph1d)
2501 [ - + ]: 6939 : ABI_FREE(ph1df)
2502 [ - + ]: 6939 : ABI_FREE(vhartr)
2503 [ - + ]: 6939 : ABI_FREE(vtrial)
2504 [ - + ]: 6939 : ABI_FREE(vpsp)
2505 [ - + ]: 6939 : ABI_FREE(vxc)
2506 [ - + ]: 6939 : ABI_FREE(vxc_hybcomp)
2507 [ - + ]: 6939 : ABI_FREE(vxctau)
2508 [ - + ]: 6939 : ABI_FREE(xccc3d)
2509 [ - + ]: 6939 : ABI_FREE(kxc)
2510 [ - + ]: 6939 : ABI_FREE(shiftvector)
2511 [ - + ]: 6939 : ABI_FREE(dtn_pc)
2512 [ - + ]: 6939 : ABI_FREE(grhf)
2513 [ - + ]: 6939 : ABI_FREE(nvresid)
2514 [ - + ]: 6939 : ABI_FREE(nvtauresid)
2515 [ - + ]: 6939 : ABI_FREE(intgden)
2516 [ - + ]: 6939 : ABI_FREE(intgden0)
2517 : :
2518 : : !Deallocate precon-object
2519 : : !call precon%save(ngfft, 1) !DEBUG
2520 : 6939 : call precon%free()
2521 : :
2522 [ + - ]: 6939 : if(allocated(vectornd)) then
2523 : 6939 : ABI_FREE(vectornd)
2524 : : end if
2525 : :
2526 [ + + ]: 6939 : if(associated(rcpaw)) then
2527 : 3 : call rcpaw_reinit(rcpaw)
2528 : : endif
2529 : :
2530 [ + + + + : 6939 : if((nstep>0.and.dtset%iscf>0).or.dtset%iscf==-1) then
+ + ]
2531 [ - + ]: 6430 : ABI_FREE(dielinv)
2532 : : end if
2533 [ - + ]: 6939 : ABI_FREE(gbound_diel)
2534 [ - + ]: 6939 : ABI_FREE(irrzondiel)
2535 [ - + ]: 6939 : ABI_FREE(kg_diel)
2536 [ - + ]: 6939 : ABI_FREE(phnonsdiel)
2537 [ - + ]: 6939 : ABI_FREE(susmat)
2538 [ - + ]: 6939 : ABI_FREE(ph1ddiel)
2539 [ - + ]: 6939 : ABI_FREE(ylmdiel)
2540 : :
2541 [ + + ]: 6939 : if (psps%usepaw==1) then
2542 [ + + ]: 1521 : if (dtset%iscf>0) then
2543 [ + + ]: 4546 : do iatom=1,my_natom
2544 : 3048 : pawrhoij(iatom)%lmnmix_sz=0
2545 : 3048 : pawrhoij(iatom)%use_rhoijres=0
2546 [ - + ]: 3048 : ABI_FREE(pawrhoij(iatom)%kpawmix)
2547 [ - + ]: 4546 : ABI_FREE(pawrhoij(iatom)%rhoijres)
2548 : : end do
2549 : : end if
2550 : : ! if (recompute_cprj.or.usecprj==1) then
2551 [ + + ]: 1521 : if (recompute_cprj) then
2552 : 66 : usecprj=0;mcprj=0
2553 : 66 : call pawcprj_free(cprj)
2554 [ + - + + : 25800 : ABI_FREE(cprj_local)
- + - + -
+ ]
2555 : : end if
2556 : 1521 : call paw_an_free(paw_an)
2557 : 1521 : call paw_ij_free(paw_ij)
2558 : 1521 : call pawfgrtab_free(pawfgrtab)
2559 [ - + ]: 1521 : if(dtset%usewvl==1) then
2560 : : #if defined HAVE_BIGDFT
2561 : : call cprj_clean(wvl%descr%paw%cprj)
2562 : : ABI_FREE(wvl%descr%paw%cprj)
2563 : : #endif
2564 : 0 : call paw2wvl_ij(2,paw_ij,wvl%descr)
2565 : : end if
2566 : : end if
2567 [ + - + + : 10032 : ABI_FREE(pawfgrtab)
- + - + -
+ - + - +
- + - + -
+ - + ]
2568 [ + - + + : 10032 : ABI_FREE(paw_an)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
2569 [ + - + + : 10032 : ABI_FREE(paw_ij)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
2570 [ - + ]: 6939 : ABI_FREE(nhat)
2571 [ - + ]: 6939 : ABI_FREE(xcctau3d)
2572 [ - + ]: 6939 : ABI_FREE(dimcprj_srt)
2573 [ - + ]: 6939 : ABI_FREE(dimcprj)
2574 : :
2575 : :
2576 : : ! Deallocate exact exchange data at the end of the calculation
2577 [ + + ]: 6939 : if (dtset%usefock==1) then
2578 [ + - ]: 123 : if (fock%fock_common%use_ACE/=0) call fock_ACE_destroy(fock%fockACE)
2579 : 123 : call fock_common_destroy(fock%fock_common)
2580 : 123 : call fock_BZ_destroy(fock%fock_BZ)
2581 : 123 : call fock_destroy(fock)
2582 : 123 : nullify(fock)
2583 : : end if
2584 : :
2585 [ + + ]: 6939 : if (prtxml == 1) then
2586 : : ! We output the final result given in results_gs
2587 : 4 : write(ab_xml_out, "(A)") ' <finalConditions>'
2588 : 4 : call out_resultsgs_XML(dtset, 4, results_gs, psps%usepaw)
2589 : 4 : write(ab_xml_out, "(A)") ' </finalConditions>'
2590 : 4 : write(ab_xml_out, "(A)") ' </scfcvLoop>'
2591 : : end if
2592 : :
2593 : : !Free the datastructure constrained_dft
2594 : 6939 : call constrained_dft_free(constrained_dft)
2595 : :
2596 : 6939 : call timab(1462,2,tsec)
2597 : 6939 : call timab(1440,2,tsec)
2598 : :
2599 : : DBG_EXIT("COLL")
2600 : :
2601 [ - + - + : 20817 : end subroutine scfcv_core
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
- - - - -
- - - - -
- - - - +
- - - - -
- - - - -
- - - - -
+ - - - -
- - - + -
+ - + - +
- + - + -
+ - + - +
+ + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
- - - - -
- + - + -
+ - + - +
- + - + -
+ ]
2602 : : !!***
2603 : :
2604 : : !!****f* ABINIT/etotfor
2605 : : !! NAME
2606 : : !! etotfor
2607 : : !!
2608 : : !! FUNCTION
2609 : : !! This routine is called to compute the total energy and various parts of it.
2610 : : !! The routine computes -if requested- the forces.
2611 : : !!
2612 : : !! INPUTS
2613 : : !! atindx1(natom)=index table for atoms, inverse of atindx
2614 : : !! dtefield <type(efield_type)> = variables related to Berry phase
2615 : : !! dtset <type(dataset_type)>=all input variables in this dataset
2616 : : !! | berryopt = 4: electric field is on -> add the contribution of the
2617 : : !! | - \Omega E.P term to the total energy
2618 : : !! | /= 4: electric field is off
2619 : : !! | bfield = cartesian coordinates of magnetic field in atomic units
2620 : : !! | efield = cartesian coordinates of the electric field in atomic units
2621 : : !! | iatfix(3,natom)=1 for frozen atom along some direction, 0 for unfrozen
2622 : : !! | ionmov=governs the movement of atoms (see help file)
2623 : : !! | densfor_pred=governs the mixed electronic-atomic part of the preconditioner
2624 : : !! | natom=number of atoms in cell.
2625 : : !! | nconeq=number of atomic constraint equations
2626 : : !! | nspden=number of spin-density components
2627 : : !! | nsym=number of symmetry elements in space group
2628 : : !! | occopt=option for occupancies
2629 : : !! | prtvol=integer controlling volume of printed output
2630 : : !! | tsmear=smearing energy or temperature (if metal)
2631 : : !! | typat(natom)=type integer for each atom in cell
2632 : : !! | wtatcon(3,natom,nconeq)=weights for atomic constraints
2633 : : !! extfpmd <type(extfpmd_type)>=extended first-principles molecular dynamics type
2634 : : !! gmet(3,3)=metric tensor for G vecs (in bohr**-2)
2635 : : !! fock <type(fock_type)>= quantities to calculate Fock exact exchange
2636 : : !! grchempottn(3,natom)=grads of spatially-varying chemical potential energy (hartree)
2637 : : !! grcondft(3,natom)=grads of constrained DFT energy (hartree)
2638 : : !! grewtn(3,natom)=grads of Ewald energy (hartree)
2639 : : !! grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree)
2640 : : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
2641 : : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
2642 : : !! kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0
2643 : : !! mgfft=maximum size of 1D FFTs
2644 : : !! mpi_enreg=information about MPI parallelization
2645 : : !! my_natom=number of atoms treated by current processor
2646 : : !! nattyp(ntypat)=number of atoms of each type
2647 : : !! nfft=(effective) number of FFT grid points (for this processor)
2648 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
2649 : : !! ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
2650 : : !! nhat(nfft,nspden*usepaw)= -PAW only- compensation density
2651 : : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
2652 : : !! ntypat=number of types of atoms in unit cell.
2653 : : !! nvresid(nfft,nspden)=potential or density residual
2654 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
2655 : : !! n3xccc=dimension of the xccc3d array (0 or nfft).
2656 : : !! optene=option for the computation of total energy
2657 : : !! (-1=no computation; 0=direct scheme; 1=double-counting scheme)
2658 : : !! optforces=option for the computation of forces
2659 : : !! optres=0 if residual array (nvresid) contains the potential residual
2660 : : !! =1 if residual array (nvresid) contains the density residual
2661 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
2662 : : !! pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
2663 : : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
2664 : : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
2665 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) information.
2666 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
2667 : : !! rhog(2,nfft)=array for Fourier transform of electron density
2668 : : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3
2669 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2670 : : !! symrec(3,3,nsym)=symmetry operations in reciprocal space
2671 : : !! ucvol=unit cell volume
2672 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
2673 : : !! usevxctau=1 if if XC functional depends on kinetic energy density
2674 : : !! vhartr(nfft)=array for holding Hartree potential
2675 : : !! vpsp(nfft)=array for holding local psp
2676 : : !! vxc(nfft,nspden)=array for holding XC potential
2677 : : !! vxctau(nfftf,dtset%nspden,4*usevxctau)]=derivative of XC energy density wrt
2678 : : !! kinetic energy density (metaGGA cases)
2679 : : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
2680 : : !! xcctau3d(n3xccc)=3D core electron kinetic energy density for XC core correction, Ha / bohr^-3 (TODO: check unit)
2681 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
2682 : : !!
2683 : : !! OUTPUT
2684 : : !! deltae=change in total energy between the previous and present SCF cycle
2685 : : !! etotal=total energy (hartree)
2686 : : !! ===== if optforces==1
2687 : : !! diffor=maximum absolute change in component of forces between present and previous SCF cycle.
2688 : : !! favg(3)=mean of fcart before correction for translational symmetry
2689 : : !! fcart(3,natom)=cartesian forces from gred (hartree/bohr)
2690 : : !! gred(3,natom)=symmetrized form of grtn (grads of Etot) (hartree)
2691 : : !! gresid(3,natom)=forces due to the residual of the density/potential
2692 : : !! grhf(3,natom)=Hellman-Feynman derivatives of the total energy
2693 : : !! grxc(3,natom)=d(Exc)/d(xred) derivatives (0 without core charges)
2694 : : !! maxfor=maximum absolute value of force
2695 : : !! synlgr(3,natom)=symmetrized form of grads of Enl (hartree)
2696 : : !!
2697 : : !! SIDE EFFECTS
2698 : : !! Input/Output:
2699 : : !! elast=previous value of the energy,
2700 : : !! needed to compute deltae, then updated.
2701 : : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
2702 : : !! energies <type(energies_type)>=all part of total energy.
2703 : : !! | entropy(IN)=entropy due to the occupation number smearing (if metal)
2704 : : !! | e_localpsp(IN)=local psp energy (hartree)
2705 : : !! | e_eigenvalues(IN)=Sum of the eigenvalues - Band energy (Hartree)
2706 : : !! | e_chempot(IN)=energy from spatially varying chemical potential (hartree)
2707 : : !! | e_ewald(IN)=Ewald energy (hartree)
2708 : : !! | e_vdw_dftd(IN)=VdW DFT-D energy
2709 : : !! | e_hartree(IN)=Hartree part of total energy (hartree units)
2710 : : !! | e_corepsp(IN)=psp core-core energy
2711 : : !! | e_hybcomp_E0(IN)=energy compensation energy for the hybrid functionals at frozen density
2712 : : !! | e_hybcomp_v0(IN)=potential compensation energy for the hybrid functionals at frozen density
2713 : : !! | e_hybcomp_v (IN)=potential compensation energy for the hybrid functionals at self-consistent density
2714 : : !! | e_kinetic(IN)=kinetic energy part of total energy.
2715 : : !! | e_nlpsp_vfock(IN)=nonlocal psp + potential Fock ACE part of total energy.
2716 : : !! | e_nucdip(IN)=energy due to array of nuclear magnetic dipoles
2717 : : !! | e_xc(IN)=exchange-correlation energy (hartree)
2718 : : !! | e_xcdc(IN)=exchange-correlation double-counting energy (hartree)
2719 : : !! | e_elecfield(OUT)=the term of the energy functional that depends explicitly
2720 : : !! | on the electric field: enefield = -ucvol*E*P
2721 : : !! | e_magfield(OUT)=the term of the energy functional that depends explicitly
2722 : : !! | on the magnetic field: e_magfield = -ucvol*E*P
2723 : : !! | e_entropy(OUT)=entropy energy due to the occupation number smearing (if metal)
2724 : : !! | this value is %entropy * dtset%tsmear (hartree).
2725 : : !! | paw%epaw(IN)=PAW spherical part energy
2726 : : !! | paw%epaw_dc(IN)=PAW spherical part double-counting energy
2727 : : !! | paw%epaw_core(IN)=PAW spherical part energy from core electrons
2728 : : !! | paw%epaw_core_dc(IN)=PAW spherical part double-counting energy from core electrons
2729 : : !! ===== if optforces==1
2730 : : !! forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)
2731 : : !! grnl(3*natom)=gradients of Etot due to nonlocal contributions
2732 : : !! Input for norm-conserving psps, output for PAW
2733 : : !! ===== if psps%usepaw==1
2734 : : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
2735 : : !! (gradients of rhoij for each atom with respect to atomic positions are computed here)
2736 : : !!
2737 : : !! NOTES
2738 : : !! In case of PAW calculations:
2739 : : !! All computations are done on the fine FFT grid.
2740 : : !! All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
2741 : : !! All arrays (densities/potentials...) are computed on this fine FFT grid.
2742 : : !! ! Developers have to be careful when introducing others arrays:
2743 : : !! they have to be stored on the fine FFT grid.
2744 : : !! In case of norm-conserving calculations the FFT grid is the usual FFT grid.
2745 : : !!
2746 : : !! SOURCE
2747 : :
2748 : 47157 : subroutine etotfor(atindx1,deltae,diffor,dtefield,dtset,&
2749 : : & elast,electronpositron,energies,&
2750 : 47157 : & etotal,favg,fcart,fock,forold,gred,gmet,grchempottn,grcondft,gresid,grewtn,grhf,grnl,grvdw,&
2751 : 47157 : & grxc,gsqcut,extfpmd,indsym,kxc,maxfor,mgfft,mpi_enreg,my_natom,nattyp,&
2752 : 47157 : & nfft,ngfft,ngrvdw,nhat,nkxc,ntypat,nvresid,n1xccc,n3xccc,optene,optforces,optres,&
2753 : 47157 : & pawang,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,red_ptot,psps,rhog,rhor,rmet,rprimd,&
2754 : 47157 : & symrec,synlgr,ucvol,usepaw,usevxctau,vhartr,vpsp,vxc,vxctau,wvl,wvl_den,xccc3d,xred,rcpaw)
2755 : : !& xcctau3d)
2756 : :
2757 : : !Arguments ------------------------------------
2758 : : !scalars
2759 : : integer,intent(in) :: my_natom,mgfft,n1xccc,n3xccc,nfft,ngrvdw,nkxc,ntypat,optene,optforces
2760 : : integer,intent(in) :: optres,usepaw,usevxctau
2761 : : real(dp),intent(in) :: gsqcut
2762 : : real(dp),intent(inout) :: elast,ucvol
2763 : : real(dp),intent(out) :: deltae,diffor,etotal,maxfor
2764 : : type(MPI_type),intent(in) :: mpi_enreg
2765 : : type(efield_type),intent(in) :: dtefield
2766 : : type(dataset_type),intent(in) :: dtset
2767 : : type(electronpositron_type),pointer :: electronpositron
2768 : : type(energies_type),intent(inout) :: energies
2769 : : type(extfpmd_type),pointer,intent(inout) :: extfpmd
2770 : : type(pawang_type),intent(in) :: pawang
2771 : : type(pseudopotential_type),intent(in) :: psps
2772 : : type(wvl_internal_type), intent(in) :: wvl
2773 : : type(wvl_denspot_type), intent(inout) :: wvl_den
2774 : : type(fock_type),pointer, intent(inout) :: fock
2775 : : type(rcpaw_type), pointer, intent(in) :: rcpaw
2776 : : !arrays
2777 : : integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom)
2778 : : integer,intent(in) :: nattyp(ntypat),ngfft(18),symrec(3,3,dtset%nsym)
2779 : : real(dp),intent(in) :: gmet(3,3),grchempottn(3,dtset%natom),grcondft(3,dtset%natom)
2780 : : real(dp),intent(in) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc)
2781 : : real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),red_ptot(3)
2782 : : real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,dtset%nspden),rmet(3,3)
2783 : : real(dp),intent(in) :: vhartr(nfft),vpsp(nfft),vxc(nfft,dtset%nspden)
2784 : : real(dp),intent(in) :: vxctau(nfft,dtset%nspden,4*usevxctau)
2785 : : real(dp),intent(in) :: xccc3d(n3xccc)
2786 : : !real(dp),intent(in),optional :: xcctau3d(n3xccc)
2787 : : real(dp),intent(inout) :: forold(3,dtset%natom),grnl(3*dtset%natom)
2788 : : real(dp),intent(inout) :: nhat(nfft,dtset%nspden*psps%usepaw)
2789 : : real(dp),intent(inout),target :: nvresid(nfft,dtset%nspden)
2790 : : real(dp),intent(inout) :: xred(3,dtset%natom)
2791 : : real(dp),intent(out) :: favg(3),gred(3,dtset%natom)
2792 : : real(dp),intent(inout) :: fcart(3,dtset%natom)
2793 : : real(dp),intent(inout) :: rprimd(3,3)
2794 : : real(dp),intent(out) :: gresid(3,dtset%natom),grhf(3,dtset%natom)
2795 : : real(dp),intent(inout) :: grxc(3,dtset%natom)
2796 : : real(dp),intent(out) :: synlgr(3,dtset%natom)
2797 : : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
2798 : : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
2799 : : type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw)
2800 : : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
2801 : :
2802 : : !Local variables-------------------------------
2803 : : !scalars
2804 : : integer :: comm_grid,dimnhat,ifft,ipositron,ispden,itypat,optgr,optgr2,option
2805 : : integer :: optnc,optstr,optstr2,iir,jjr,kkr
2806 : : logical :: apply_residual
2807 : : real(dp) :: eenth,ucvol_
2808 : : !arrays
2809 : : real(dp),parameter :: k0(3)=(/zero,zero,zero/)
2810 : : real(dp) :: tsec(2),A(3,3),A1(3,3),A_new(3,3),efield_new(3)
2811 : : real(dp) :: dummy(0),nhat_dum(0,0)
2812 : 47157 : real(dp),allocatable :: vlocal(:,:)
2813 : 47157 : real(dp), ABI_CONTIGUOUS pointer :: resid(:,:)
2814 : :
2815 : : ! *********************************************************************
2816 : :
2817 : 47157 : call timab(80,1,tsec)
2818 : :
2819 : 47157 : ipositron=electronpositron_calctype(electronpositron)
2820 : :
2821 [ + - ]: 47157 : if (optene>-1) then
2822 : :
2823 : 47157 : call entropy(dtset,energies)
2824 : :
2825 : : ! Turn it into an electric enthalpy, refer to Eq.(33) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]],
2826 : : ! the missing volume is added here
2827 : 47157 : energies%e_elecfield=zero
2828 [ + + + - ]: 47157 : if ((dtset%berryopt==4.or.dtset%berryopt==14).and.ipositron/=1) then
2829 [ + + ]: 700 : energies%e_elecfield=-dot_product(dtset%red_efieldbar,red_ptot) !!ebar_i p_i
2830 : 175 : eenth=zero
2831 [ + + ]: 700 : do iir=1,3
2832 [ + + ]: 2275 : do jjr=1,3
2833 : 2100 : eenth=eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !!g^{-1})_ij ebar_i ebar_j
2834 : : end do
2835 : : end do
2836 : 175 : energies%e_elecfield=energies%e_elecfield-eenth*ucvol/(8._dp*pi)
2837 : : end if
2838 : :
2839 : : ! Turn it into an internal energy, refer to Eq.(36) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]],
2840 : : ! but a little different: U=E_ks + (vol/8*pi) * g^{-1})_ij ebar_i ebar_j
2841 [ - + - - ]: 47157 : if ((dtset%berryopt==6.or.dtset%berryopt==16).and.ipositron/=1) then
2842 : : energies%e_elecfield=zero
2843 : 0 : eenth=zero
2844 [ # # ]: 0 : do iir=1,3
2845 [ # # ]: 0 : do jjr=1,3
2846 : 0 : eenth=eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! g^{-1})_ij ebar_i ebar_j
2847 : : end do
2848 : : end do
2849 : 0 : energies%e_elecfield=energies%e_elecfield+eenth*ucvol/(8._dp*pi)
2850 : : end if
2851 : :
2852 : : ! Calculate internal energy and electric enthalpy for mixed BC case.
2853 [ - + - - ]: 47157 : if (dtset%berryopt==17.and.ipositron/=1) then
2854 : 0 : energies%e_elecfield=zero
2855 [ # # # # ]: 0 : A(:,:)=(four_pi/ucvol)*rmet(:,:)
2856 : 0 : A1(:,:)=A(:,:) ; A_new(:,:)=A(:,:)
2857 [ # # ]: 0 : efield_new(:)=dtset%red_efield(:)
2858 : : eenth=zero
2859 [ # # ]: 0 : do kkr=1,3
2860 [ # # ]: 0 : if (dtset%jfielddir(kkr)==1) then ! fixed ebar direction
2861 : : ! step 1 add -ebar*p
2862 : 0 : eenth=eenth-dtset%red_efieldbar(kkr)*red_ptot(kkr)
2863 : : ! step 2 chang to e_new (change e to ebar)
2864 : 0 : efield_new(kkr)=dtset%red_efieldbar(kkr)
2865 : : ! step 3 chang matrix A to A1
2866 [ # # ]: 0 : do iir=1,3
2867 [ # # ]: 0 : do jjr=1,3
2868 [ # # ]: 0 : if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr)
2869 [ # # ]: 0 : if ((iir==kkr .and. jjr/=kkr) .or. (iir/=kkr .and. jjr==kkr)) &
2870 : 0 : & A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr)
2871 [ # # ]: 0 : if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr)
2872 : : end do
2873 : : end do
2874 : 0 : A(:,:)=A1(:,:) ; A_new(:,:)=A1(:,:)
2875 : : end if
2876 : : end do ! end for kkr
2877 [ # # ]: 0 : do iir=1,3
2878 [ # # ]: 0 : do jjr=1,3
2879 : 0 : eenth= eenth+half*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr)
2880 : : end do
2881 : : end do
2882 : 0 : energies%e_elecfield=energies%e_elecfield+eenth
2883 : : end if ! berryopt==17
2884 : :
2885 : : ! Turn it into a magnetic enthalpy, by adding orbital electronic contribution
2886 : 47157 : energies%e_magfield = zero
2887 : : ! if (dtset%berryopt == 5 .and. ipositron/=1) then
2888 : : ! emag = dot_product(mag_cart,dtset%bfield)
2889 : : ! energies%e_magfield = emag
2890 : : ! end if
2891 : :
2892 : : ! Compute total (free)- energy by direct scheme
2893 [ + + ]: 47157 : if (optene==0) then
2894 : : etotal = energies%e_kinetic + energies%e_hartree + energies%e_xc &
2895 : : & + energies%e_localpsp + energies%e_corepsp &
2896 : : & + energies%e_entropy + energies%e_elecfield &
2897 : : & + energies%e_magfield + energies%e_nucdip &
2898 : : & + energies%e_hybcomp_E0 - energies%e_hybcomp_v0 + energies%e_hybcomp_v &
2899 : : & + energies%e_constrained_dft + energies%e_ewald &
2900 : 32483 : & + energies%e_chempot + energies%e_vdw_dftd
2901 : : ! +two*energies%e_fock-energies%e_fock0 ! The Fock energy is already included in the non-local one
2902 : : ! +energies%e_nlpsp_vfock - energies%e_fock0
2903 : :
2904 : : ! See similar section in m_energies.F90
2905 : : ! XG 20181025 This gives a variational energy in case of NCPP with all bands occupied - not yet for metals.
2906 [ + + ]: 32483 : if (usepaw==0) etotal = etotal + energies%e_nlpsp_vfock - energies%e_fock0
2907 : : ! XG 20181025 I was expecting the following to give also a variational energy in case of PAW, but this is not true.
2908 : : ! if (usepaw==1) etotal = etotal + energies%paw%epaw + energies%e_nlpsp_vfock - energies%e_fock0
2909 : : ! XG 20181025 So, the following is giving a non-variational expression ...
2910 [ + + ]: 32483 : if (usepaw==1) etotal = etotal + energies%paw%epaw + energies%e_fock
2911 [ - + ]: 32483 : if (dtset%usedmft/=0) etotal = etotal + energies%e_hu - energies%e_dc
2912 : : end if
2913 : :
2914 : : ! Compute total (free) energy by double-counting scheme
2915 [ + + ]: 47157 : if (optene==1) then
2916 : : etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc &
2917 : : & - energies%e_xcdc + energies%e_corepsp - energies%e_corepspdc- energies%e_fock0 &
2918 : : & + energies%e_entropy + energies%e_elecfield + energies%e_magfield &
2919 : 14674 : & + energies%e_hybcomp_E0 - energies%e_hybcomp_v0 + energies%e_constrained_dft
2920 : 14674 : etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
2921 [ + + ]: 14674 : if (usepaw/=0) etotal = etotal + energies%paw%epaw_dc
2922 [ + + ]: 14674 : if (dtset%usedmft/=0) etotal = etotal + energies%e_hu - energies%e_dc
2923 : : end if
2924 : :
2925 : : ! Additional stuff for electron-positron
2926 [ + + ]: 47157 : if (dtset%positron/=0) then
2927 [ + + ]: 562 : if (ipositron==0) then
2928 : 69 : energies%e_electronpositron =zero
2929 : 69 : energies%edc_electronpositron=zero
2930 : : else
2931 : 493 : energies%e_electronpositron =electronpositron%e_hartree+electronpositron%e_xc
2932 : 493 : energies%edc_electronpositron=electronpositron%e_hartree+electronpositron%e_xcdc
2933 [ + + ]: 493 : if (usepaw==1) then
2934 : 488 : energies%e_electronpositron =energies%e_electronpositron +electronpositron%e_paw
2935 : 488 : energies%edc_electronpositron=energies%edc_electronpositron+electronpositron%e_pawdc
2936 : : end if
2937 : : end if
2938 [ - + ]: 562 : if (optene==0) electronpositron%e0=etotal
2939 [ + - ]: 562 : if (optene==1) electronpositron%e0=etotal-energies%edc_electronpositron
2940 : 562 : etotal=electronpositron%e0+energies%e0_electronpositron+energies%e_electronpositron
2941 : : end if
2942 : :
2943 : : ! Add the extfpmd energy contribution to the internal energy
2944 [ + + ]: 47157 : if (associated(extfpmd)) then
2945 : 90 : energies%e_extfpmd=extfpmd%e_kinetic
2946 : 90 : energies%edc_extfpmd=extfpmd%edc_kinetic
2947 [ - + ]: 90 : if (optene==0) etotal=etotal+energies%e_extfpmd
2948 [ + - ]: 90 : if (optene==1) etotal=etotal+energies%edc_extfpmd
2949 : : end if
2950 : :
2951 : : ! Add the PAW core energy contribution to the internal energy
2952 [ + + ]: 47157 : if(associated(rcpaw)) then
2953 : 21 : energies%paw%epaw_core=rcpaw%ehnzc+rcpaw%ekinc
2954 : 21 : energies%paw%epaw_core_dc=rcpaw%eeigc-rcpaw%edcc+rcpaw%ehnzc
2955 [ + - ]: 21 : if (ipositron/=1) then
2956 [ + + ]: 42 : do itypat=1,dtset%ntypat
2957 : 21 : energies%paw%epaw_core=energies%paw%epaw_core+pawtab(itypat)%exccore*rcpaw%atm(itypat)%mult
2958 : 42 : energies%paw%epaw_core_dc=energies%paw%epaw_core_dc+pawtab(itypat)%exccore*rcpaw%atm(itypat)%mult
2959 : : enddo
2960 : : endif
2961 [ - + ]: 21 : if(optene==0) etotal=etotal+energies%paw%epaw_core
2962 [ + - ]: 21 : if(optene==1) etotal=etotal+energies%paw%epaw_core_dc
2963 : : end if
2964 : :
2965 : : ! Compute energy residual
2966 : 47157 : deltae=etotal-elast
2967 : 47157 : elast=etotal
2968 : : end if !optene/=-1
2969 : :
2970 : 47157 : call timab(80,2,tsec)
2971 : :
2972 : : !------Compute forces-----------------------------------------------------
2973 : :
2974 [ + + ]: 47157 : if (optforces==1) then
2975 : :
2976 : : ! PAW: add gradients due to Dij derivatives to non-local term
2977 [ + + ]: 20737 : if (usepaw==1) then
2978 [ - + + - : 24776 : ABI_MALLOC(vlocal,(nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
2979 [ + + ]: 12665 : do ispden=1,min(dtset%nspden,2)
2980 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vhartr,vlocal,vpsp,vxc)
2981 [ + + ]: 332323958 : do ifft=1,nfft
2982 : 332317764 : vlocal(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)
2983 : : end do
2984 : : end do
2985 : :
2986 [ + + ]: 6194 : if(dtset%nspden==4)then
2987 [ + + ]: 120 : do ispden=3,4
2988 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vlocal,vxc)
2989 [ + + ]: 655480 : do ifft=1,nfft
2990 : 655440 : vlocal(ifft,ispden)=vxc(ifft,ispden)
2991 : : end do
2992 : : end do
2993 : : end if
2994 : 6194 : ucvol_=ucvol
2995 : : #if defined HAVE_BIGDFT
2996 : : if (dtset%usewvl==1) ucvol_=product(wvl_den%denspot%dpbox%hgrids)*real(product(wvl_den%denspot%dpbox%ndims),dp)
2997 : : #endif
2998 : 6194 : dimnhat=0;optgr=1;optgr2=0;optstr=0;optstr2=0
2999 : 6194 : comm_grid=mpi_enreg%comm_fft;if(dtset%usewvl==1) comm_grid=mpi_enreg%comm_wvl
3000 : : call pawgrnl(atindx1,dimnhat,dummy,1,dummy,grnl,gsqcut,mgfft,my_natom, &
3001 : : & dtset%natom, nattyp,nfft,ngfft,nhat_dum,dummy,dtset%nspden,dtset%nsym,ntypat,optgr,optgr2,optstr,optstr2,&
3002 : : & pawang,pawfgrtab,pawrhoij,pawtab,ph1d,psps,k0,rprimd,symrec,dtset%typat,ucvol_,vlocal,vxc,xred, &
3003 : : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,mpi_comm_grid=mpi_enreg%comm_fft,&
3004 : 6194 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,paral_kgb=mpi_enreg%paral_kgb)
3005 [ - + ]: 6194 : ABI_FREE(vlocal)
3006 : : end if
3007 : :
3008 : : apply_residual=(optres==1 .and. dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1 .and. &
3009 [ + + + - : 20737 : & abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5)
+ - + - ]
3010 : :
3011 : : ! If residual is a density residual (and forces from residual asked),
3012 : : ! has to convert it into a potential residual before calling forces routine
3013 : : if (apply_residual) then
3014 [ - + + - : 24056 : ABI_MALLOC(resid,(nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
3015 [ - + ]: 6014 : option=0; if (dtset%densfor_pred<0) option=1
3016 [ + + - + ]: 6014 : optnc=1;if (dtset%nspden==4.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) optnc=2
3017 : : call nres2vres(dtset,gsqcut,usepaw,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,&
3018 : : & nkxc,nvresid,n3xccc,optnc,option,pawang,pawfgrtab,pawrhoij,pawtab,&
3019 : 6014 : & rhor,rprimd,usepaw,resid,xccc3d,xred,vxc)
3020 : : else
3021 : 14723 : resid => nvresid
3022 : : end if
3023 : : call forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,forold,gred,grchempottn,grcondft,gresid,grewtn,&
3024 : : & grhf,grnl,grvdw,grxc,gsqcut,indsym,maxfor,mgfft,mpi_enreg,&
3025 : : & n1xccc,n3xccc,nattyp,nfft,ngfft,ngrvdw,ntypat,pawrad,pawtab,&
3026 : : & ph1d,psps,rhog,rhor,rprimd,symrec,synlgr,dtset%usefock,usevxctau,resid,vxc,vxctau,wvl,wvl_den,xred,&
3027 : 20737 : & electronpositron=electronpositron)
3028 [ + + ]: 20737 : if (apply_residual) then
3029 [ - + ]: 6014 : ABI_FREE(resid)
3030 : : end if
3031 : :
3032 : : ! Returned gred are full symmetrized gradients of Etotal
3033 : : ! wrt reduced coordinates xred, d(Etotal)/d(xred)
3034 : : ! Forces are contained in array fcart
3035 : :
3036 : : else ! if optforces==0
3037 [ + + + + ]: 346052 : fcart=zero
3038 [ + + + + ]: 346052 : gred=zero
3039 : 26420 : favg=zero
3040 : 26420 : diffor=zero
3041 [ + + + + ]: 346052 : gresid=zero
3042 [ + + + + ]: 346052 : grhf=zero
3043 : 26420 : maxfor=zero
3044 [ + + + + ]: 346052 : synlgr=zero
3045 : : end if
3046 : :
3047 : 47157 : call timab(80,2,tsec)
3048 : :
3049 : 47157 : end subroutine etotfor
3050 : : !!***
3051 : :
3052 : : !!****f* ABINIT/wf_mixing
3053 : : !! NAME
3054 : : !! wf_mixing
3055 : : !!
3056 : : !! FUNCTION
3057 : : !! Mixing of wavefunctions in the outer loop of a double loop SCF approach.
3058 : : !! Different algorithms are implemented, depending on the value of wfmixalg.
3059 : : !!
3060 : : !! INPUTS
3061 : : !! atindx1(dtset%natom)=index table for atoms, inverse of atindx
3062 : : !! dtset <type(dataset_type)>=all input variables in this dataset
3063 : : !! istep=number of call the routine (usually the outer loop in the SCF double loop)
3064 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
3065 : : !! mcprj=size of cprj array
3066 : : !! mpi_enreg=information about MPI parallelization
3067 : : !! nattyp(dtset%ntypat)=number of atoms of each type in cell.
3068 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
3069 : : !! pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
3070 : : !!
3071 : : !! SIDE EFFECTS
3072 : : !! cg(2,mcg)= plane wave wavefunction coefficient
3073 : : !! Value from previous SCF cycle is input and stored in some form
3074 : : !! Extrapolated value is output
3075 : : !! cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
3076 : : !! Value from previous SCF cycle is input and stored in some form
3077 : : !! Extrapolated value is output
3078 : : !! scf_history_wf <type(scf_history_type)>=arrays obtained from previous SCF cycles
3079 : : !!
3080 : : !! SOURCE
3081 : :
3082 : 43 : subroutine wf_mixing(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,&
3083 : 43 : & nattyp,npwarr,pawtab,scf_history_wf)
3084 : :
3085 : : use m_cgcprj, only : dotprod_set_cgcprj, dotprodm_sumdiag_cgcprj, lincom_cgcprj, cgcprj_cholesky
3086 : :
3087 : : !Arguments ------------------------------------
3088 : : !scalars
3089 : : integer,intent(in) :: istep,mcg,mcprj
3090 : : type(MPI_type),intent(in) :: mpi_enreg
3091 : : type(dataset_type),intent(in) :: dtset
3092 : : type(scf_history_type),intent(inout) :: scf_history_wf
3093 : : !arrays
3094 : : integer,intent(in) :: atindx1(dtset%natom),nattyp(dtset%ntypat)
3095 : : integer,intent(in) :: npwarr(dtset%nkpt)
3096 : : real(dp), intent(inout) :: cg(2,mcg)
3097 : : type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj)
3098 : : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
3099 : :
3100 : : !Local variables-------------------------------
3101 : : !scalars
3102 : : integer :: hermitian
3103 : : integer :: ibdmix,ibdsp,ibg,ibg_hist,icg,icg_hist
3104 : : integer :: ierr,ikpt,indh,ind_biorthog,ind_biorthog_eff,ind_newwf,ind_residual,inplace
3105 : : integer :: iorder,iset2,isppol,istep_cycle,istep_new,istwf_k,kk,me_distrb,my_nspinor
3106 : : integer :: nband_k,nbdmix,npw_k,nset1,nset2,ntypat
3107 : : integer :: shift_set1,shift_set2,spaceComm_band,spare_mem,usepaw,wfmixalg
3108 : : real(dp) :: alpha,beta
3109 : : complex(dp) :: sum_coeffs
3110 : : !arrays
3111 : 43 : integer,allocatable :: ipiv(:),dimcprj(:)
3112 : : real(dp) :: tsec(2)
3113 : 43 : real(dp),allocatable :: al(:,:),mmn(:,:,:)
3114 : 43 : real(dp),allocatable :: dotprod_res(:,:,:),dotprod_res_k(:,:,:),res_mn(:,:,:),smn(:,:,:)
3115 : 43 : complex(dp),allocatable :: coeffs(:)
3116 : 43 : type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kh(:,:)
3117 : : ! *************************************************************************
3118 : :
3119 : : !DEBUG
3120 : : !write(std_out,*)
3121 : : !write(std_out,*)' wf_mixing : enter, istep= ',istep
3122 : : !call flush(std_out)
3123 : : !write(std_out,*)' istep,scf_history_wf%alpha=',istep,scf_history_wf%alpha
3124 : : !write(std_out,*)' cg(1,1)=',cg(1,1)
3125 : : !write(std_out,*)' scf_history_wf%cg(1,1,1:5)=',scf_history_wf%cg(1,1,1:5)
3126 : : !ABI_MALLOC(cg_ref,(2,mcg))
3127 : : !cg_ref(:,:)=cg(:,:)
3128 : : !ABI_MALLOC(cprj_ref,(dtset%natom,mcprj))
3129 : : !cprj_ref(:,:)=cprj(:,:)
3130 : : ! write(std_out,*)' scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)=',&
3131 : : !& scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)
3132 : : ! call flush(std_out)
3133 : : !ENDDEBUG
3134 : :
3135 [ + - ]: 43 : if (istep==0) return
3136 : :
3137 : 43 : ntypat=dtset%ntypat
3138 : 43 : usepaw=dtset%usepaw
3139 : 43 : wfmixalg=scf_history_wf%wfmixalg
3140 : 43 : nbdmix=dtset%nbandhf
3141 : 43 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
3142 : 43 : me_distrb=mpi_enreg%me_kpt
3143 : 43 : spaceComm_band=xmpi_comm_self
3144 : :
3145 : : spare_mem=0
3146 : 43 : if(scf_history_wf%history_size==wfmixalg-1)spare_mem=1
3147 : :
3148 : : !scf_history_wf%alpha contains dtset%wfmix
3149 : 43 : alpha=scf_history_wf%alpha
3150 : 43 : beta=one-scf_history_wf%alpha
3151 : 43 : icg=0
3152 : 43 : icg_hist=0
3153 : 43 : ibg=0
3154 : 43 : ibg_hist=0
3155 : :
3156 : : !Useful array
3157 [ - + + - : 129 : ABI_MALLOC(dimcprj,(dtset%natom))
+ - + - -
+ - + ]
3158 [ + + ]: 43 : if (usepaw==1) then
3159 : 16 : iorder=0 ! There is no change of ordering in the mixing of wavefunctions
3160 : 16 : call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O')
3161 : : end if
3162 : :
3163 [ + + ]: 43 : if(istep==1)then
3164 [ + + ]: 28 : do indh=1,scf_history_wf%history_size
3165 : 28 : call pawcprj_alloc(scf_history_wf%cprj(:,:,indh),0,dimcprj)
3166 : : end do
3167 : : end if
3168 : :
3169 [ - + + - : 628 : ABI_MALLOC(cprj_k,(dtset%natom,my_nspinor*nbdmix))
- + + - +
- + - - +
- + + + +
+ ]
3170 [ - + + - : 585 : ABI_MALLOC(cprj_kh,(dtset%natom,my_nspinor*nbdmix))
- + + - -
+ - + + +
+ + ]
3171 [ + + ]: 43 : if(usepaw==1) then
3172 : 16 : call pawcprj_alloc(cprj_k,0,dimcprj)
3173 : 16 : call pawcprj_alloc(cprj_kh,0,dimcprj)
3174 : : end if
3175 [ - + + - : 172 : ABI_MALLOC(smn,(2,nbdmix,nbdmix))
- + + - +
- + - - +
- + ]
3176 [ - + + - : 129 : ABI_MALLOC(mmn,(2,nbdmix,nbdmix))
- + + - -
+ - + ]
3177 : :
3178 [ + + ]: 43 : if(wfmixalg>2)then
3179 : 16 : nset1=1
3180 : 16 : nset2=min(istep-1,wfmixalg-1)
3181 [ + + + - : 48 : ABI_MALLOC(dotprod_res_k,(2,1,nset2))
+ - + + -
+ - + ]
3182 [ + + + - : 32 : ABI_MALLOC(dotprod_res,(2,1,nset2))
- + - + ]
3183 [ - + + - : 64 : ABI_MALLOC(res_mn,(2,wfmixalg-1,wfmixalg-1))
- + + - +
- - + -
+ ]
3184 [ + + + + : 216 : dotprod_res=zero
+ + ]
3185 [ + + ]: 16 : if(istep==1)then
3186 [ + + + + : 750 : scf_history_wf%dotprod_sumdiag_cgcprj_ij=zero
+ + ]
3187 : : end if
3188 : : end if
3189 : :
3190 : : !Explanation for the index for the wavefunction stored in scf_history_wf
3191 : : !The reference is the cg+cprj output after the wf optimization at istep 1.
3192 : : !It comes as input to the present routine as cgcprj input at step 2, and is usually found at indh=1.
3193 : :
3194 : : !In the simple mixing case (wfmixalg==2), the reference is never stored, because it is used "on-the-fly" to biothogonalize the
3195 : : !previous input (that was stored in indh=1), then generate the next input, which is stored again in indh=1
3196 : :
3197 : : !When the storage is not spared:
3198 : : !- the values of indh from 2 to wfmixalg store the (computed here) biorthogonalized input cgcprj, then the residual
3199 : : !- the values of indh from wfmixalg+1 to 2*wfmixalg-1 store the biorthogonalized output cgcprj (coming as argument)
3200 : :
3201 : : !First step
3202 [ + + + + : 43 : if (istep==1 .or. (wfmixalg==2 .and. abs(scf_history_wf%alpha-one)<tol8) ) then
- + ]
3203 : :
3204 : 4 : indh=2 ! This input wavefunction is NOT the reference
3205 [ + + ]: 4 : if(wfmixalg==2)indh=1 ! But this does not matter in the simple mixing case that has history_size=1
3206 : :
3207 : : ! Simply store the wavefunctions and cprj. However, nband_k might be different from nbandhf...
3208 : : ! LOOP OVER SPINS
3209 [ + + ]: 8 : do isppol=1,dtset%nsppol
3210 : :
3211 : : ! BIG FAT k POINT LOOP
3212 [ + + ]: 21 : do ikpt=1,dtset%nkpt
3213 : :
3214 : : ! Select k point to be treated by this proc
3215 : 13 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
3216 [ - + ]: 13 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
3217 : :
3218 : 13 : npw_k=npwarr(ikpt)
3219 : :
3220 [ + + + + ]: 32137 : scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
3221 [ + + ]: 13 : if(usepaw==1) then
3222 : : ! scf_history_wf%cprj(:,ibg_hist+1:ibg_hist+my_nspinor*nbdmix,1)=cprj(:,ibg+1:ibg+my_nspinor*nbdmix)
3223 : : call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,iorder,isppol,dtset%mband,&
3224 : : & dtset%mkmem,dtset%natom,nbdmix,nband_k,my_nspinor,dtset%nsppol,0,&
3225 : 4 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
3226 : : call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,&
3227 : : & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
3228 : 4 : & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
3229 : : end if
3230 : :
3231 : : ! Update the counters
3232 : 13 : ibg=ibg+my_nspinor*nband_k
3233 : 13 : ibg_hist=ibg_hist+my_nspinor*nbdmix
3234 : 13 : icg=icg+my_nspinor*nband_k*npw_k
3235 : 17 : icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
3236 : :
3237 : : end do
3238 : : end do
3239 : :
3240 : : else
3241 : : ! From istep==2
3242 : :
3243 : : ! First part of the computation : biorthogonalization, and computation of the residual (possibly, prediction of the next input in the case of simple mixing)
3244 : : ! Index for the wavefunctions stored in scf_history_wf whose scalar products with the argument cgcprj will have to be computed.
3245 : 14 : indh=1 ! This input wavefunction is the reference
3246 [ + + ]: 14 : if(wfmixalg/=2 .and. istep==2)indh=2 ! except for istep=2 in the rmm-diis
3247 : :
3248 [ + - ]: 14 : if(wfmixalg>2)then
3249 : : ! istep inside the cycle defined by wfmixalg, and next index. Then, indices of the wavefunction sets.
3250 : 14 : istep_cycle=mod((istep-2),wfmixalg-1)
3251 : 14 : istep_new=mod((istep-1),wfmixalg-1)
3252 : 14 : ind_biorthog=1+wfmixalg+istep_cycle
3253 : 14 : ind_residual=2+istep_cycle
3254 : 14 : ind_newwf=2+istep_new
3255 : 14 : shift_set1=ind_residual-1
3256 : 14 : shift_set2=1
3257 : : end if
3258 : :
3259 : : ! LOOP OVER SPINS
3260 [ + + ]: 78 : do isppol=1,dtset%nsppol
3261 : :
3262 : : ! BIG FAT k POINT LOOP
3263 [ + + ]: 219 : do ikpt=1,dtset%nkpt
3264 : :
3265 : : ! Select k point to be treated by this proc
3266 : 141 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
3267 [ - + ]: 141 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
3268 : :
3269 : 141 : istwf_k=dtset%istwfk(ikpt)
3270 : 141 : npw_k=npwarr(ikpt)
3271 : :
3272 : : ! Biorthogonalization
3273 : :
3274 [ + + ]: 141 : if(usepaw==1) then
3275 : : call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,iorder,isppol,dtset%mband,&
3276 : : & dtset%mkmem,dtset%natom,nbdmix,nband_k,my_nspinor,dtset%nsppol,0,&
3277 : 60 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
3278 : : call pawcprj_get(atindx1,cprj_kh,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,&
3279 : : & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,my_nspinor,dtset%nsppol,0,&
3280 : 60 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
3281 : : end if !end usepaw=1
3282 : :
3283 : 141 : hermitian=0
3284 [ + + + + ]: 141 : if(wfmixalg==2 .or. istep==2)then
3285 : : call dotprod_set_cgcprj(atindx1,cg,scf_history_wf%cg(:,:,indh),cprj_k,cprj_kh,dimcprj,hermitian,&
3286 : : & 0,0,icg,icg_hist,ikpt,isppol,istwf_k,nbdmix,mcg,mcg,mcprj,mcprj,dtset%mkmem,&
3287 : 105 : & mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn,usepaw)
3288 : : else
3289 : : call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),cg,cprj_kh,cprj_k,dimcprj,hermitian,&
3290 : : & 0,0,icg_hist,icg,ikpt,isppol,istwf_k,nbdmix,mcg,mcg,mcprj,mcprj,dtset%mkmem,&
3291 : 36 : & mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn,usepaw)
3292 : : end if
3293 : :
3294 : : ! Invert S matrix, that is NOT hermitian.
3295 : : ! Calculate M=S^-1
3296 [ + + + + : 7233 : mmn=zero
+ + ]
3297 [ + + ]: 690 : do kk=1,nbdmix
3298 : 690 : mmn(1,kk,kk)=one
3299 : : end do
3300 : :
3301 [ - + + - : 423 : ABI_MALLOC(ipiv,(nbdmix))
+ - + - -
+ - + ]
3302 : : ! The smn is destroyed by the following inverse call
3303 : 141 : call zgesv(nbdmix,nbdmix,smn,nbdmix,ipiv,mmn,nbdmix,ierr)
3304 [ - + ]: 141 : ABI_CHECK(ierr == 0, sjoin('zgesv general inversion routine returned ierr:', itoa(ierr)))
3305 [ - + ]: 141 : ABI_FREE(ipiv)
3306 : :
3307 : : ! The M matrix is used to compute the biorthogonalized set of wavefunctions, and to store it at the proper place
3308 [ + + + + ]: 141 : if(wfmixalg==2 .or. istep==2)then
3309 : 105 : inplace=1
3310 : : call lincom_cgcprj(mmn,scf_history_wf%cg(:,:,indh),cprj_kh,dimcprj,&
3311 : 105 : & icg_hist,inplace,mcg,my_nspinor*nbdmix,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw)
3312 : : else
3313 : 36 : inplace=0
3314 : : call lincom_cgcprj(mmn,cg,cprj_k,dimcprj,&
3315 : : & icg,inplace,mcg,my_nspinor*nbdmix,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw,&
3316 : 36 : & cgout=scf_history_wf%cg(:,:,ind_biorthog),cprjout=scf_history_wf%cprj(:,:,ind_biorthog),icgout=icg_hist)
3317 : : end if
3318 : :
3319 : : ! The biorthogonalised set of wavefunctions is now stored at the proper place
3320 : :
3321 : : ! Finalize this first part of the computation, depending on the algorithm and the step.
3322 : :
3323 [ + + ]: 141 : if(wfmixalg==2)then
3324 : :
3325 : : ! Wavefunction extrapolation, simple mixing case
3326 : : ! alpha contains dtset%wfmix, beta contains one-alpha
3327 : : cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=&
3328 : : & alpha*cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)&
3329 [ + + + + ]: 170560 : & +beta*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)
3330 [ + + ]: 100 : if(usepaw==1) then
3331 [ + + ]: 300 : do ibdmix=1,nbdmix
3332 : 300 : call pawcprj_axpby(beta,alpha,cprj_kh(:,ibdmix:ibdmix),cprj_k(:,ibdmix:ibdmix))
3333 : : end do ! end loop on ibdmix
3334 : : end if
3335 : :
3336 : : ! Back to usual orthonormalization
3337 : : call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,my_nspinor*nband_k,dtset%mkmem,&
3338 : 100 : & mpi_enreg,dtset%natom,nattyp,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw)
3339 : :
3340 : : ! Store the newly extrapolated wavefunctions, orthonormalized, in scf_history_wf
3341 [ + + + + ]: 170560 : scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
3342 [ + + ]: 100 : if(usepaw==1) then
3343 [ + + ]: 300 : do ibdmix=1,nbdmix
3344 : : call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,&
3345 : : & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
3346 : 300 : & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
3347 : : end do ! end loop on ibdmix
3348 : : end if
3349 : :
3350 : : else ! wfmixalg/=2
3351 : : ! RMM-DIIS
3352 : :
3353 [ + + ]: 41 : if (istep==2)then
3354 : : ! Store the argument wf as the reference for all future steps, in scf_history_wf with index 1.
3355 [ + + + + ]: 17837 : scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,1)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
3356 [ - + ]: 5 : if(usepaw==1) then
3357 [ # # ]: 0 : do ibdmix=1,nbdmix
3358 : : call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,1),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,&
3359 : : & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
3360 : 0 : & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
3361 : : end do ! end loop on ibdmix
3362 : : end if
3363 : : end if
3364 : :
3365 : 41 : ind_biorthog_eff=ind_biorthog
3366 [ + + ]: 41 : if(istep==2)ind_biorthog_eff=1 ! The argument wf has not been stored in ind_biorthog
3367 : : ! Compute the residual of the wavefunctions for this istep,
3368 : : ! that replaces the previously stored set of (biorthogonalized) input wavefunctions
3369 : : scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)=&
3370 : : & scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog_eff)&
3371 [ + + + + ]: 124337 : & -scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)
3372 [ - + ]: 41 : if(usepaw==1) then
3373 [ # # ]: 0 : do ibdmix=1,nbdmix
3374 : : call pawcprj_axpby(one,-one,scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog_eff),&
3375 : 0 : & scf_history_wf%cprj(:,ibdmix:ibdmix,ind_residual))
3376 : : end do ! end loop on ibdmix
3377 : : end if
3378 : :
3379 : : ! Compute the new scalar products to fill the res_mn matrix
3380 : : call dotprodm_sumdiag_cgcprj(atindx1,scf_history_wf%cg,scf_history_wf%cprj,dimcprj,&
3381 : : & ibg_hist,icg_hist,ikpt,isppol,istwf_k,nbdmix,mcg,mcprj,dtset%mkmem,&
3382 : : & mpi_enreg,scf_history_wf%history_size,dtset%natom,nattyp,nbdmix,npw_k,nset1,nset2,my_nspinor,dtset%nsppol,ntypat,&
3383 : 41 : & shift_set1,shift_set2,pawtab,dotprod_res_k,usepaw)
3384 : :
3385 [ - + + + : 702 : dotprod_res=dotprod_res+dotprod_res_k
+ + + + ]
3386 : :
3387 : : ! scf_history_wf for index ind_biorthog will contain the extrapolated wavefunctions (and no more the output of the SCF loop).
3388 : : scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog)=&
3389 : : & scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog_eff)+&
3390 [ + + + + ]: 124337 : & (alpha-one)*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)
3391 [ - + ]: 41 : if(usepaw==1) then
3392 [ # # ]: 0 : do ibdmix=1,nbdmix
3393 [ # # ]: 0 : if(ind_biorthog/=ind_biorthog_eff)then
3394 [ # # # # : 0 : scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog)=scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog_eff)
# # # # #
# # # # #
# # ]
3395 : : end if
3396 : : call pawcprj_axpby((alpha-one),one,scf_history_wf%cprj(:,ibdmix:ibdmix,ind_residual),&
3397 : 0 : & scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog))
3398 : : end do ! end loop on ibdmix
3399 : : end if
3400 : :
3401 : : end if
3402 : :
3403 : 141 : ibg=ibg+my_nspinor*nband_k
3404 : 141 : ibg_hist=ibg_hist+my_nspinor*nbdmix
3405 : 141 : icg=icg+my_nspinor*nband_k*npw_k
3406 : 180 : icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
3407 : :
3408 : : end do ! End big k point loop
3409 : : end do ! End loop over spins
3410 : :
3411 : : end if ! istep>=2
3412 : :
3413 [ + + + + ]: 43 : if(wfmixalg>2 .and. istep>1)then
3414 : :
3415 : : !DEBUG
3416 : : ! write(std_out,*)' '
3417 : : ! write(std_out,*)' Entering the residual minimisation part '
3418 : : ! write(std_out,*)' '
3419 : : ! call flush(std_out)
3420 : : !ENDDEBUG
3421 : :
3422 : 14 : call timab(48,1,tsec)
3423 : 14 : call xmpi_sum(dotprod_res,mpi_enreg%comm_kpt,ierr)
3424 : 14 : call timab(48,2,tsec)
3425 : :
3426 [ + + + + ]: 164 : scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,1+shift_set1,1+shift_set2:nset2+shift_set2)=dotprod_res(:,1,1:nset2)
3427 [ + + ]: 64 : scf_history_wf%dotprod_sumdiag_cgcprj_ij(1,1+shift_set2:nset2+shift_set2,1+shift_set1)=dotprod_res(1,1,1:nset2)
3428 [ + + ]: 64 : scf_history_wf%dotprod_sumdiag_cgcprj_ij(2,1+shift_set2:nset2+shift_set2,1+shift_set1)=-dotprod_res(2,1,1:nset2)
3429 : :
3430 : : end if ! wfmixalg>2 and istep>1
3431 : :
3432 [ + + ]: 16 : if(wfmixalg>2 .and. istep>2)then
3433 : :
3434 : : ! Extract the relevant matrix R_mn
3435 : : res_mn(:,1:nset2,1:nset2)=&
3436 [ + + + + : 684 : & scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,1+shift_set2:nset2+shift_set2,1+shift_set2:nset2+shift_set2)
+ + ]
3437 : :
3438 : : !DEBUG
3439 : : ! write(std_out,*)' The matrix res_mn(:,1:nset2,1:nset2) is :'
3440 : : ! write(std_out,*)res_mn(:,1:nset2,1:nset2)
3441 : : ! call flush(std_out)
3442 : : !ENDDEBUG
3443 : :
3444 : : ! Solve R_mn \alpha_n = 1_m
3445 [ - + + - : 36 : ABI_MALLOC(ipiv,(nset2))
+ - + - -
+ - + ]
3446 [ - + + - : 36 : ABI_MALLOC(coeffs,(nset2))
+ - + - -
+ - + ]
3447 [ + + ]: 60 : coeffs(:)=cone
3448 : : ! The res_mn is destroyed by the following inverse call
3449 : 12 : call zgesv(nset2,1,res_mn,wfmixalg-1,ipiv,coeffs,nset2,ierr)
3450 [ - + ]: 12 : ABI_CHECK(ierr == 0, sjoin('zgesv general inversion routine returned ierr:', itoa(ierr)))
3451 [ - + ]: 12 : ABI_FREE(ipiv)
3452 : : ! The coefficients must sum to one
3453 [ + + ]: 60 : sum_coeffs=sum(coeffs)
3454 [ + + ]: 60 : coeffs=coeffs/sum_coeffs
3455 : :
3456 : : !DEBUG
3457 : : ! write(std_out,*)' The coefficients that minimize the residual have been found'
3458 : : ! write(std_out,*)' coeffs =',coeffs
3459 : : ! call flush(std_out)
3460 : : !ENDDEBUG
3461 : : end if ! wfmixalg>2 and istep>2
3462 : :
3463 [ + + ]: 16 : if(wfmixalg>2 .and. istep>1)then
3464 : :
3465 : : ! Find the new "input" wavefunction, bi-orthogonalized, and store it replacing the adequate "old" input wavefunction.
3466 : :
3467 : 14 : icg=0
3468 : 14 : icg_hist=0
3469 : 14 : ibg=0
3470 : 14 : ibg_hist=0
3471 [ - + + - : 42 : ABI_MALLOC(al,(2,nset2))
+ - + - -
+ - + ]
3472 [ + + ]: 14 : if(istep>2)then
3473 [ + + ]: 60 : do iset2=1,nset2
3474 : 60 : al(1,iset2)=real(coeffs(iset2)) ; al(2,iset2)=aimag(coeffs(iset2))
3475 : : end do
3476 : : else
3477 : 2 : al(1,1)=one ; al(2,1)=zero
3478 : : end if
3479 : :
3480 : : !DEBUG
3481 : : ! write(std_out,*)' Overload the coefficients, in order to simulate a simple mixing with wfmix '
3482 : : ! write(std_out,*)' Set al(1,ind_biorthog-3)=one, for ind_biorthog=',ind_biorthog
3483 : : ! write(std_out,*)' This will feed scf_history for set ind_biorthog-3+wfmixalg=',ind_biorthog-3+wfmixalg
3484 : : ! al(:,:)=zero
3485 : : ! al(1,ind_biorthog-3)=one
3486 : : ! call flush(std_out)
3487 : : !ENDDEBUG
3488 : :
3489 : : ! LOOP OVER SPINS
3490 [ + + ]: 28 : do isppol=1,dtset%nsppol
3491 : :
3492 : : ! BIG FAT k POINT LOOP
3493 [ + + ]: 69 : do ikpt=1,dtset%nkpt
3494 : :
3495 : : ! Select k point to be treated by this proc
3496 : 41 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
3497 [ - + ]: 41 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
3498 : :
3499 : 41 : istwf_k=dtset%istwfk(ikpt)
3500 : 41 : npw_k=npwarr(ikpt)
3501 : :
3502 [ + + ]: 41 : if(istep>2)then
3503 : : ! Make the appropriate linear combination (from the extrapolated wfs)
3504 [ + + + + ]: 106500 : cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=zero
3505 [ + + ]: 186 : do iset2=1,nset2
3506 : : cg(1,icg+1:icg+my_nspinor*npw_k*nband_k)=cg(1,icg+1:icg+my_nspinor*npw_k*nband_k)&
3507 : : & +al(1,iset2)*scf_history_wf%cg(1,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)&
3508 [ + + ]: 141926 : & -al(2,iset2)*scf_history_wf%cg(2,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)
3509 : : cg(2,icg+1:icg+my_nspinor*npw_k*nband_k)=cg(2,icg+1:icg+my_nspinor*npw_k*nband_k)&
3510 : : & +al(1,iset2)*scf_history_wf%cg(2,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)&
3511 [ + + ]: 141962 : & +al(2,iset2)*scf_history_wf%cg(1,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)
3512 : : end do
3513 : : else ! One needs a simple copy from the extrapolated wavefunctions
3514 [ + + + + ]: 17837 : cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,1+wfmixalg)
3515 : : end if
3516 : : ! Note the storage in cprj_k. By the way, a simple copy might also be used in case istep=2.
3517 [ - + ]: 41 : if(usepaw==1) then
3518 [ # # ]: 0 : do ibdsp=1,my_nspinor*nbdmix
3519 : 0 : call pawcprj_lincom(al,scf_history_wf%cprj(:,ibdsp,1+wfmixalg:nset2+wfmixalg),cprj_k(:,ibdsp:ibdsp),nset2)
3520 : : end do
3521 : : end if
3522 : :
3523 : : ! Store the newly extrapolated wavefunctions for this k point, still bi-orthonormalized, in scf_history_wf
3524 [ + + + + ]: 124337 : scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_newwf)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)
3525 [ - + ]: 41 : if(usepaw==1) then
3526 : : call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,ind_newwf),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,&
3527 : : & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
3528 : 0 : & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
3529 : : end if
3530 : :
3531 : : ! Back to usual orthonormalization for the cg and cprj_k
3532 : : call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,my_nspinor*nband_k,dtset%mkmem,&
3533 : 41 : & mpi_enreg,dtset%natom,nattyp,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw)
3534 : :
3535 : : ! Need to transfer cprj_k to cprj
3536 [ - + ]: 41 : if(usepaw==1) then
3537 : : call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,iorder,isppol,&
3538 : : & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,&
3539 : 0 : & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb)
3540 : : end if
3541 : :
3542 : 41 : ibg=ibg+my_nspinor*nband_k
3543 : 41 : ibg_hist=ibg_hist+my_nspinor*nbdmix
3544 : 41 : icg=icg+my_nspinor*nband_k*npw_k
3545 : 55 : icg_hist=icg_hist+my_nspinor*nbdmix*npw_k
3546 : :
3547 : : end do ! End big k point loop
3548 : : end do ! End loop over spins
3549 : :
3550 [ + + ]: 14 : if(istep>2)then
3551 [ - + ]: 12 : ABI_FREE(coeffs)
3552 : : end if
3553 [ - + ]: 14 : ABI_FREE(al)
3554 : :
3555 : : end if ! wfmixalg>2 and istep>1
3556 : :
3557 : : !DEBUG
3558 : : ! write(std_out,*)' wf_mixing : exit '
3559 : : ! write(std_out,*)' scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)=',&
3560 : : !& scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)
3561 : : ! write(std_out,*)' cg(1:2,1:2)=',cg(1:2,1:2)
3562 : : ! write(std_out,*)' scf_history_wf%cg(1:2,1:2,1)=',scf_history_wf%cg(1:2,1:2,1)
3563 : : ! ABI_FREE(cg_ref)
3564 : : ! ABI_FREE(cprj_ref)
3565 : : !ENDDEBUG
3566 : :
3567 [ + + ]: 43 : if(usepaw==1) then
3568 : 16 : call pawcprj_free(cprj_k)
3569 : 16 : call pawcprj_free(cprj_kh)
3570 : : end if
3571 [ + - + + : 345 : ABI_FREE(cprj_k)
- + - + -
+ ]
3572 [ + - + + : 345 : ABI_FREE(cprj_kh)
- + - + -
+ ]
3573 [ - + ]: 43 : ABI_FREE(dimcprj)
3574 [ - + ]: 43 : ABI_FREE(mmn)
3575 [ - + ]: 43 : ABI_FREE(smn)
3576 [ + + ]: 43 : if(wfmixalg>2)then
3577 [ - + ]: 16 : ABI_FREE(dotprod_res_k)
3578 [ - + ]: 16 : ABI_FREE(dotprod_res)
3579 [ - + ]: 16 : ABI_FREE(res_mn)
3580 : : end if
3581 : :
3582 [ - + - + : 43 : end subroutine wf_mixing
- + - + -
+ - - - -
- - - + -
- - - - -
- + ]
3583 : : !!***
3584 : :
3585 : : end module m_scfcv_core
3586 : : !!***
|
