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