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