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 13836 : subroutine scfcv_core(atindx,atindx1,cg,cprj,cpus,dmatpawu,dtefield,dtfil,dtpawuj,&
272 13836 : & dtset,ecore,eigen,electronpositron,fatvshift,hdr,extfpmd,rcpaw,indsym,&
273 13836 : & initialized,irrzon,itimes,kg,mcg,mcprj,mpi_enreg,my_natom,nattyp,ndtpawuj,nfftf,npwarr,occ,&
274 13836 : & paw_dmft,pawang,pawfgr,pawrad,pawrhoij,pawtab,phnons,psps,pwind,&
275 13836 : & pwind_alloc,pwnsfac,rec_set,resid,results_gs,rhog,rhor,rprimd,&
276 6918 : & 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 6918 : 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 6918 : 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 6918 : integer,allocatable :: dimcprj(:),dimcprj_srt(:)
383 6918 : integer,allocatable :: gbound_diel(:,:),irrzondiel(:,:,:),kg_diel(:,:)
384 6918 : integer,allocatable :: l_size_atm(:)
385 6918 : integer,allocatable :: indsym_dum(:,:,:),symrec_dum(:,:,:), rmm_diis_status(:,:,:)
386 13836 : logical :: eijkl_is_sym(psps%ntypat)
387 6918 : 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 13836 : 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 13836 : 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 6918 : real(dp),allocatable :: dielinv(:,:,:,:,:),dtn_pc(:,:)
407 6918 : real(dp),allocatable :: fcart(:,:),forold(:,:),gred(:,:),gresid(:,:)
408 6918 : real(dp),allocatable :: grchempottn(:,:),grcondft(:,:),grewtn(:,:)
409 6918 : real(dp),allocatable :: grhf(:,:),grnl(:),grvdw(:,:),grxc(:,:)
410 6918 : real(dp),allocatable :: intgres(:,:),kxc(:,:),nhat(:,:),nhatgr(:,:,:),nvresid(:,:),nvtauresid(:,:)
411 6918 : real(dp),allocatable :: ph1d(:,:),ph1ddiel(:,:),ph1df(:,:)
412 6918 : real(dp),allocatable :: phnonsdiel(:,:,:),rhowfg(:,:),rhowfr(:,:),shiftvector(:)
413 6918 : real(dp),allocatable :: susmat(:,:,:,:,:),synlgr(:,:)
414 6918 : real(dp),allocatable :: vectornd(:,:,:),vhartr(:),vpsp(:),vtrial(:,:)
415 6918 : real(dp),allocatable :: vxc(:,:),vxc_hybcomp(:,:),vxctau(:,:,:),workr(:,:),xccc3d(:),xcctau3d(:),ylmdiel(:,:)
416 6918 : real(dp),allocatable :: intgden(:,:),intgden0(:,:)
417 6918 : real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:)
418 6918 : type(scf_history_type) :: scf_history_wf
419 6918 : type(constrained_dft_t) :: constrained_dft
420 6918 : type(paw_an_type),allocatable :: paw_an(:)
421 6918 : type(paw_ij_type),allocatable :: paw_ij(:)
422 6918 : type(pawfgrtab_type),allocatable,save :: pawfgrtab(:)
423 6918 : type(pawrhoij_type),pointer :: pawrhoij_ep(:)
424 : type(fock_type),pointer :: fock
425 6918 : 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 6918 : call timab(1440,1,tsec)
433 6918 : call timab(1441,3,tsec)
434 :
435 : ! enable time limit handler if not done in callers.
436 6918 : if (enable_timelimit_in(MY_NAME) == MY_NAME) then
437 6922 : 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 6918 : non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
442 :
443 : !######################################################################
444 : !Initializations - Memory allocations
445 : !----------------------------------------------------------------------
446 6918 : lmax_diel = 0
447 :
448 : !MPI communicators
449 6918 : if (xmpi_paral==1.and.mpi_enreg%paral_hf==1) then
450 0 : spaceComm=mpi_enreg%comm_kpt
451 : else
452 6918 : spaceComm=mpi_enreg%comm_cell
453 : end if
454 6918 : me=xmpi_comm_rank(spaceComm)
455 6918 : spaceComm_fft=mpi_enreg%comm_fft
456 6918 : spaceComm_wvl=mpi_enreg%comm_wvl
457 6918 : me_wvl=mpi_enreg%me_wvl
458 6918 : spaceComm_grid=mpi_enreg%comm_fft
459 6918 : if(dtset%usewvl==1) spaceComm_grid=mpi_enreg%comm_wvl
460 6918 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
461 :
462 : !Save some variables from dataset definition
463 6918 : nstep=dtset%nstep
464 6918 : ecut=dtset%ecut
465 6918 : ecutf=ecut; if (psps%usepaw==1) ecutf=dtset%pawecutdg
466 6918 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) ecutf=dtset%pawecutdg
467 6918 : iscf10=mod(dtset%iscf,10)
468 6918 : tollist(1)=dtset%tolmxf;tollist(2)=dtset%tolwfr
469 6918 : tollist(3)=dtset%toldff;tollist(4)=dtset%toldfe
470 6918 : tollist(6)=dtset%tolvrs;tollist(7)=dtset%tolrff
471 6918 : tollist(8)=dtset%vdw_df_threshold
472 6918 : tollist(9)=dtset%toldmag
473 6918 : 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 6918 : & dtset%berryopt == 17)
480 :
481 : !Get electronic temperature from dtset
482 13794 : 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 131442 : ngfft(:)=dtset%ngfft(:)
487 6918 : if (psps%usepaw==1) then
488 28500 : 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 6918 : zion=zero
495 27254 : do iatom=1,dtset%natom
496 27254 : zion=zion+psps%ziontypat(dtset%typat(iatom))
497 : end do
498 :
499 : !Compute different geometric tensor, as well as ucvol, from rprimd
500 6918 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
501 :
502 : !Fock: be sure that the pointer is initialized to Null.
503 6918 : 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 6918 : & .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 6918 : 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 6918 : 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 6918 : nullify(grhor,lrhor,elfr)
529 6918 : quit=0 ; dbl_nnsclo=0 ; conv_retcode=0
530 6918 : dielop=0 ; strsxc=zero
531 6918 : deltae=zero ; elast=zero ;
532 6918 : 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 6918 : maxfor=zero
536 : !
537 6918 : results_gs%residm=zero;results_gs%res2=zero
538 6918 : results_gs%deltae=zero;results_gs%diffor=zero
539 6918 : call energies%init()
540 6918 : 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 6918 : 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 13835 : select case (dtset%usepotzero)
558 : case (0,1)
559 6917 : energies%e_corepsp = ecore / ucvol
560 6917 : 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 6918 : energies%e_corepspdc = zero
566 : end select
567 6918 : if(wvlbigdft) energies%e_corepsp = zero
568 6918 : if(dtset%icutcoul.ne.3) energies%e_corepsp = zero
569 :
570 6918 : fermie=energies%e_fermie
571 6918 : fermih=energies%e_fermih
572 6918 : isave_den=0; isave_kden=0 !initial index of density protection file
573 6918 : optres=merge(0,1,dtset%iscf<10)
574 6918 : usexcnhat=0!;mcprj=0
575 6918 : initialized0=initialized
576 6918 : if (dtset%tfkinfunc==12) tfw_activated=.true.
577 6918 : ipert=0;idir=0;cplex=1
578 6918 : istep_mix=1
579 6918 : istep_fock_outer=1
580 6918 : ipositron=electronpositron_calctype(electronpositron)
581 6918 : unit_out=0;if (dtset%prtvol >= 10) unit_out=ab_out
582 27672 : nfftotf=product(ngfftf(1:3))
583 :
584 6918 : usecprj=0
585 6918 : if (mcprj>0) then
586 78 : usecprj=1
587 : end if
588 :
589 : !Stresses and forces flags
590 6918 : forces_needed=0;prtfor=0
591 : if ((dtset%optforces==1.or.dtset%ionmov==4.or.dtset%ionmov==5.or.&
592 6918 : & 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 3883 : else if (dtset%iscf>0.and.dtset%optforces==2) then
597 3270 : forces_needed=2
598 : end if
599 :
600 6918 : stress_needed=0
601 6918 : 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 6918 : & .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 6918 : etotal = results_gs%etotal
609 :
610 : !Entering a scfcv_core loop, printing data to XML file if required.
611 6918 : prtxml=0;if (me==0.and.dtset%prtxml==1) prtxml=1
612 6918 : 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 6918 : choice=1 ; diffor=zero ; res2=zero
627 20754 : ABI_MALLOC(fcart,(3,dtset%natom))
628 13836 : ABI_MALLOC(gred,(3,dtset%natom))
629 27672 : ABI_MALLOC(intgden,(dtset%nspden,dtset%natom))
630 20754 : ABI_MALLOC(intgden0,(dtset%nspden,dtset%natom))
631 88262 : gred(:,:)=zero
632 88262 : fcart(:,:)=results_gs%fcart(:,:) ! This is a side effect ...
633 49863 : intgden(:,:)=zero
634 698342 : 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 6918 : & psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode)
644 :
645 : !Various allocations (potentials, gradients, ...)
646 20754 : ABI_MALLOC(forold,(3,dtset%natom))
647 13836 : ABI_MALLOC(grchempottn,(3,dtset%natom))
648 13836 : ABI_MALLOC(grcondft,(3,dtset%natom))
649 13836 : ABI_MALLOC(gresid,(3,dtset%natom))
650 13836 : ABI_MALLOC(grewtn,(3,dtset%natom))
651 20754 : ABI_MALLOC(grnl,(3*dtset%natom))
652 13836 : ABI_MALLOC(grxc,(3,dtset%natom))
653 13836 : ABI_MALLOC(synlgr,(3,dtset%natom))
654 20754 : ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
655 20754 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*dtset%natom))
656 20754 : ABI_MALLOC(vhartr,(nfftf))
657 27672 : ABI_MALLOC(vtrial,(nfftf,dtset%nspden))
658 13836 : ABI_MALLOC(vpsp,(nfftf))
659 20754 : ABI_MALLOC(vxc,(nfftf,dtset%nspden))
660 6918 : usevxctau=merge(1,0,xc_need_kden(dtset%ixc))
661 34641 : ABI_MALLOC(vxctau,(nfftf,dtset%nspden,4*usevxctau))
662 :
663 6918 : wfmixalg=dtset%fockoptmix/100
664 6918 : use_hybcomp=0
665 6918 : if(mod(dtset%fockoptmix,100)==11)use_hybcomp=1
666 27672 : ABI_MALLOC(vxc_hybcomp,(nfftf,dtset%nspden*use_hybcomp))
667 :
668 6918 : ngrvdw=0;if (dtset%vdw_xc>=5.and.dtset%vdw_xc<=7) ngrvdw=dtset%natom
669 20754 : ABI_MALLOC(grvdw,(3,ngrvdw))
670 :
671 6918 : ngrcondft=0
672 16395 : if(any(dtset%constraint_kind(:)/=0)) ngrcondft=dtset%natom
673 27672 : ABI_MALLOC(intgres,(dtset%nspden,ngrcondft))
674 6918 : if(ngrcondft/=0)then
675 606 : intgres(:,:)=zero
676 : endif
677 :
678 88262 : grchempottn(:,:)=zero
679 88262 : grcondft(:,:)=zero
680 169606 : forold(:,:)=zero ; gresid(:,:)=zero ; pel(:)=zero
681 6918 : strscondft(:)=zero
682 365338468 : vtrial(:,:)=zero; vxc(:,:)=zero
683 6918 : n1xccc=0;if (psps%n1xccc/=0) n1xccc=psps%n1xccc
684 6918 : n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf
685 20754 : ABI_MALLOC(xccc3d,(n3xccc))
686 20754 : ABI_MALLOC(xcctau3d,(nfftf*dtset%usekden))
687 :
688 : !Allocations/initializations for PAW only
689 6918 : lpawumax=-1
690 6918 : if(psps%usepaw==1) then
691 : ! Variables/arrays related to the fine FFT grid
692 6000 : ABI_MALLOC(nhat,(nfftf,dtset%nspden*psps%usepaw))
693 212967 : if (nstep==0) nhat=zero
694 7555 : ABI_MALLOC(pawfgrtab,(my_natom))
695 1500 : if (my_natom>0) then
696 : call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,&
697 1352 : & mpi_atmtab=mpi_enreg%my_atmtab)
698 : call pawfgrtab_init(pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat,&
699 1352 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
700 1352 : ABI_FREE(l_size_atm)
701 : end if
702 1500 : compch_fft=-1.d5
703 3595 : usexcnhat=maxval(pawtab(:)%usexcnhat)
704 1500 : 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 7555 : ABI_MALLOC(paw_ij,(my_natom))
710 7555 : ABI_MALLOC(paw_an,(my_natom))
711 1500 : call paw_an_nullify(paw_an)
712 1500 : call paw_ij_nullify(paw_ij)
713 1500 : has_dijhat=0;if (dtset%iscf==22) has_dijhat=1
714 1500 : has_vhartree=0; if (dtset%prtvha > 0 .or. dtset%prtvclmb > 0) has_vhartree=1
715 1500 : has_dijfock=0; if (dtset%usefock == 1) has_dijfock=1
716 17660 : has_dijnd=0;if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1
717 1500 : has_dijU=merge(0,1,dtset%usepawu>0) !Be careful on this!
718 1500 : 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 1500 : & 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 1500 : & zora=dtset%zora,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
729 1500 : if(dtset%usewvl==1) then
730 0 : call paw2wvl_ij(1,paw_ij,wvl%descr)
731 : end if
732 1500 : compch_sph=-1.d5
733 4500 : ABI_MALLOC(dimcprj,(dtset%natom))
734 3000 : ABI_MALLOC(dimcprj_srt,(dtset%natom))
735 1500 : call pawcprj_getdim(dimcprj ,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'R')
736 1500 : call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
737 3595 : do itypat=1,dtset%ntypat
738 3595 : if (pawtab(itypat)%usepawu/=0) lpawumax=max(pawtab(itypat)%lpawu,lpawumax)
739 : end do
740 1500 : 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 1500 : 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 3595 : eijkl_is_sym=.true.
751 :
752 : ! Allocation of projected WF (optional)
753 1500 : 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 1500 : nullify(pawrhoij_ep);if(associated(electronpositron))pawrhoij_ep=>electronpositron%pawrhoij_ep
777 1500 : 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 6918 : if (dtset%iscf>=0) then
790 6407 : dielar(1)=dtset%diecut;dielar(2)=dtset%dielng
791 6407 : dielar(3)=dtset%diemac;dielar(4)=dtset%diemix
792 6407 : dielar(5)=dtset%diegap;dielar(6)=dtset%dielam
793 6407 : dielar(7)=dtset%diemix;if (dtset%iscf>=10) dielar(7)=dtset%diemixmag
794 25628 : ABI_MALLOC(nvresid,(nfftf,dtset%nspden))
795 6407 : if (dtset%iscf>=10) then
796 5632 : ABI_MALLOC(nvtauresid,(nfftf,dtset%nspden*dtset%usekden))
797 : else
798 19996 : ABI_MALLOC(nvtauresid,(nfftf,dtset%nspden*usevxctau))
799 : end if
800 6407 : if (nstep==0) then
801 279698 : nvresid=zero
802 20 : nvtauresid=zero
803 : end if
804 19221 : 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 12814 : ABI_MALLOC(grhf,(3,dtset%natom))
808 : ! Additional allocation for mixing within PAW
809 6407 : npawmix=0
810 6407 : if(psps%usepaw==1) then
811 4487 : do iatom=1,my_natom
812 3010 : itypat=pawrhoij(iatom)%itypat
813 3010 : pawrhoij(iatom)%use_rhoijres=1
814 3010 : sz1=pawrhoij(iatom)%cplex_rhoij*pawtab(itypat)%lmn2_size
815 3010 : sz2=pawrhoij(iatom)%nspden
816 12040 : ABI_MALLOC(pawrhoij(iatom)%rhoijres,(sz1,sz2))
817 7000 : do ispden=1,pawrhoij(iatom)%nspden
818 367113 : pawrhoij(iatom)%rhoijres(:,ispden)=zero
819 : end do
820 9030 : ABI_MALLOC(pawrhoij(iatom)%kpawmix,(pawtab(itypat)%lmnmix_sz))
821 3010 : pawrhoij(iatom)%lmnmix_sz=pawtab(itypat)%lmnmix_sz
822 186641 : pawrhoij(iatom)%kpawmix=pawtab(itypat)%kmix
823 : npawmix=npawmix+pawrhoij(iatom)%nspden*pawtab(itypat)%lmnmix_sz &
824 4487 : & *pawrhoij(iatom)%cplex_rhoij*pawrhoij(iatom)%qphase
825 : end do
826 : end if
827 6407 : if (dtset%iscf > 0) then
828 6407 : denpot = ABI_MIXING_POTENTIAL
829 6407 : if (dtset%iscf > 10) denpot = ABI_MIXING_DENSITY
830 6407 : 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 5234 : ispmix=ABI_MIXING_REAL_SPACE;nfftmix=nfftf;ngfftmix(:)=ngfftf(:)
834 : end if
835 : !TRangel: added to avoid segfaults with Wavelets
836 6407 : 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 6407 : & dtset%npulayit, pulayhist_storage=dtset%pulayhiststore)
839 6407 : if (errid /= AB7_NO_ERROR) then
840 0 : ABI_ERROR(msg)
841 : end if
842 6407 : 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 6407 : 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 6918 : & 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 16395 : 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 6918 : 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 6409 : afford=1
888 6409 : if(dtset%iscf==-1) then
889 : ! dtset%iprcel=21
890 22 : afford=0
891 : end if
892 :
893 : ! First compute dimensions
894 6409 : 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 6329 : npwdiel=1
938 6329 : mgfftdiel=1
939 6329 : nfftdiel=1
940 6329 : lmax_diel=0
941 6329 : afford=0
942 : end if
943 :
944 : ! Now, performs allocation
945 38454 : ABI_MALLOC(dielinv,(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden))
946 32045 : ABI_MALLOC(susmat,(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden))
947 19227 : ABI_MALLOC(kg_diel,(3,npwdiel))
948 25636 : ABI_MALLOC(gbound_diel,(2*mgfftdiel+8,2))
949 25636 : ABI_MALLOC(irrzondiel,(nfftdiel**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
950 25636 : ABI_MALLOC(phnonsdiel,(2,nfftdiel**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
951 19227 : ABI_MALLOC(ph1ddiel,(2,3*(2*mgfftdiel+1)*dtset%natom*psps%usepaw))
952 25636 : ABI_MALLOC(ylmdiel,(npwdiel,lmax_diel**2))
953 : ! Then, compute the values of different arrays
954 6409 : 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 6409 : 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 6918 : nkxc=0
1002 : !TDDFT - For a first coding
1003 6918 : if (dtset%iscf==-1 .and. dtset%nspden==1) nkxc=2
1004 6918 : if (dtset%iscf==-1 .and. dtset%nspden==2) nkxc=3
1005 : !Eventually need kxc-LDA when susceptibility matrix has to be computed
1006 6918 : 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 6918 : & 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 6918 : if (precon%use_kxc) then
1016 1 : nkxc = precon%nkxc
1017 : end if
1018 6918 : if (nkxc>0) then
1019 1334 : call check_kxc(dtset%ixc,dtset%optdriver)
1020 : end if
1021 27672 : ABI_MALLOC(kxc,(nfftf,nkxc))
1022 6918 : 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 6918 : moved_atm_inside=0
1028 : !This flag will be set to 1 if the forces are computed inside the iteration.
1029 6918 : computed_forces=0
1030 :
1031 6918 : if(dtset%wfoptalg==2)then
1032 0 : ABI_MALLOC(shiftvector,((dtset%mband+2)*dtset%nkpt))
1033 6918 : val_min=-1.0_dp
1034 6918 : val_max=zero
1035 : else
1036 6918 : 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 6918 : & 0,quit,istep,ucvol,unit_out,psps%usepaw,xred,ylm,ylmgr)
1051 :
1052 6918 : if (dtset%iscf==22) energies%h0=zero
1053 :
1054 6918 : call timab(1441,2,tsec)
1055 6918 : 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 6918 : quitsum_request = xmpi_request_null; timelimit_exit = 0
1068 6918 : istep_updatedfock=0
1069 :
1070 272153 : ABI_ICALLOC(rmm_diis_status, (2, dtset%nkpt, dtset%nsppol))
1071 :
1072 : ! start SCF loop
1073 : ABI_NVTX_START_RANGE(NVTX_SCF)
1074 48120 : do istep=1,max(1,nstep)
1075 :
1076 47343 : 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 47343 : if (istep == 1) prev = abi_wtime()
1084 47343 : if (istep > 1) then
1085 40425 : now = abi_wtime()
1086 40425 : wtime_step = now - prev
1087 40425 : prev = now
1088 : call wrtout(std_out, sjoin("{SCF_istep:", itoa(istep-1), ", Vnl|psi>:", itoa(nonlop_counter), &
1089 40425 : ", wall_time: '", sec2str(wtime_step), "'} <<< TIME"))
1090 40425 : nonlop_counter = 0
1091 :
1092 40425 : 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 47343 : call timab(1442,1,tsec)
1110 47343 : 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 7024 : call symmetrize_xred(dtset%natom,dtset%nsym,dtset%symrel,dtset%tnons,xred,indsym=indsym)
1118 :
1119 7024 : if (dtset%usewvl == 0) then
1120 : ! Get cut-off for g-vectors
1121 7024 : if (psps%usepaw==1) call wrtout(std_out,' FFT (fine) grid used in SCF cycle:')
1122 7024 : 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 7024 : call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
1126 :
1127 7024 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
1128 1431 : 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 7024 : if(gemm_nonlop_use_gemm) call reset_gemm_nonlop()
1136 :
1137 : ! Initialization of atomic data for PAW
1138 7024 : 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 1500 : nremit=0
1142 1500 : if(dtset%ionmov>0 .and. itimes(1)==1)nremit=1
1143 1500 : if(dtset%imgmov>0 .and. itimes(2)==1)nremit=mpi_enreg%my_nimage
1144 1500 : 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 1500 : optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm
1148 : if ((forces_needed==1).or. &
1149 1500 : & (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 1500 : 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 1500 : & 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 7024 : moved_rhor=0
1177 :
1178 2791 : if (initialized/=0.and.dtset%usewvl == 0.and.ipositron/=1.and. &
1179 7024 : & (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 7024 : 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 89876 : if ( any(abs(dtset%nucdipmom(:,:))>tol8) ) then
1220 24 : with_vectornd = 1
1221 : else
1222 7000 : with_vectornd = 0
1223 : end if
1224 7024 : if(allocated(vectornd)) then
1225 106 : ABI_FREE(vectornd)
1226 : end if
1227 35120 : ABI_MALLOC(vectornd,(with_vectornd*nfftf,dtset%nspden,3))
1228 11773123 : vectornd=zero
1229 7024 : if(with_vectornd .EQ. 1) then
1230 : call make_vectornd(1,gsqcut,psps%usepaw,mpi_enreg,dtset%natom,nfftf,ngfftf,&
1231 24 : & dtset%nspden,dtset%nucdipmom,rprimd,vectornd,xred)
1232 : endif
1233 :
1234 : end if ! moved_atm_inside==1 .or. istep==1
1235 :
1236 47343 : call timab(1442,2,tsec)
1237 :
1238 : !Initialize/Update data in the case of an Exact-exchange (Hartree-Fock) or hybrid XC calculation
1239 47343 : hyb_mixing=zero;hyb_mixing_sr=zero
1240 47343 : 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 47343 : call timab(1447,1,tsec)
1382 :
1383 : ! Initialize/update data in the electron-positron case
1384 47343 : 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 47343 : call timab(1447,2,tsec)
1397 47343 : call timab(1448,3,tsec)
1398 :
1399 : if ((moved_atm_inside==1 .or. istep==1).or.&
1400 47343 : & (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 7110 : nhatgrdim=0
1405 7110 : dummy_nhatgr = .False.
1406 : if (psps%usepaw==1.and.(dtset%positron>=0.or.ipositron/=1) &
1407 7110 : & .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 942 : call timab(558,1,tsec)
1411 942 : nhatgrdim=0;if (dtset%xclevel==2) nhatgrdim=usexcnhat*dtset%pawnhatxc
1412 942 : ider=2*nhatgrdim;izero=0
1413 942 : if (nhatgrdim>0) then
1414 140 : ABI_MALLOC(nhatgr,(cplex*nfftf,dtset%nspden,3*nhatgrdim))
1415 : else
1416 914 : ABI_MALLOC(nhatgr,(0,0,0))
1417 914 : 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 942 : & gpu_thread_limit=dtset%gpu_thread_limit)
1426 942 : 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 942 : 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 7110 : optene = 4 * optres
1444 7110 : if(dtset%iscf==-3) optene=4
1445 7110 : if (wvlbigdft) optene = 1 ! VH needed for the WF mixing
1446 :
1447 7110 : 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 7110 : & 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 34110 : if(dtset%usepotzero==2) vpsp(:) = vpsp(:) + ecore / ( zion * ucvol )
1465 :
1466 7110 : if(dtset%optdriver==RUNL_GWLS) call build_vxc(vxc,nfftf,dtset%nspden)
1467 :
1468 7110 : if ((nhatgrdim>0.and.nstep>0).or.dummy_nhatgr) then
1469 7110 : ABI_FREE(nhatgr)
1470 : end if
1471 :
1472 : ! Recursion Initialisation
1473 7110 : 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 47343 : call timab(1448,2,tsec)
1487 47343 : 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 47343 : if (psps%usepaw==1)then
1495 :
1496 : ! Local exact exch.: impose occ. matrix if required
1497 15129 : 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 15129 : nzlmopt=0;if (istep_mix==2.and.dtset%pawnzlm>0) nzlmopt=-1
1505 15129 : if (istep_mix>2) nzlmopt=dtset%pawnzlm
1506 15129 : call paw_an_reset_flags(paw_an) ! Force the recomputation of on-site potentials
1507 15129 : call paw_ij_reset_flags(paw_ij,self_consistent=.true.) ! Force the recomputation of Dij
1508 15129 : 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 15129 : & 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 90774 : & ' average electrostatic smooth potential [Ha] , [eV]',SUM(vpotzero(:)),SUM(vpotzero(:))*Ha_eV
1526 15129 : call wrtout(std_out, msg)
1527 767672886 : vtrial(:,:)=vtrial(:,:)+SUM(vpotzero(:))
1528 15129 : if(option/=1)then
1529 : ! Fix the direct total energy (non-zero only for charged systems)
1530 42810 : 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 42810 : 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 15129 : 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 15129 : 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 15129 : 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 15129 : & 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 15129 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1583 15129 : 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 15129 : if(dtset%usewvl==1) then
1589 0 : call paw2wvl_ij(3,paw_ij,wvl%descr)
1590 : end if
1591 :
1592 15129 : 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 47343 : 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 47343 : call timab(1449,2,tsec)
1603 :
1604 : ! No need to continue and call vtorho, when nstep==0
1605 47343 : if(nstep==0)exit
1606 :
1607 : ! ######################################################################
1608 : ! The following steps are done only when nstep>0
1609 : ! ----------------------------------------------------------------------
1610 47308 : call timab(1450,1,tsec)
1611 :
1612 47308 : if(dtset%iscf>=0)then
1613 46812 : write(msg, '(a,a,i4)' )ch10,' ITER STEP NUMBER ',istep
1614 46812 : 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 47308 : moved_atm_inside=0
1619 : ! /< Hack to remove iapp from scfcv_core
1620 : ! for ionmov 4|5 ncycle=1
1621 : ! Hence iapp = itime
1622 47308 : if ( dtset%jdtset /= scfcv_jdtset ) then
1623 : ! new dtset -> reinitialize
1624 3940 : scfcv_jdtset = dtset%jdtset
1625 3940 : scfcv_itime = 0
1626 : end if
1627 47308 : if ( istep==1 ) scfcv_itime = scfcv_itime + 1
1628 47308 : if(dtset%ionmov==4 .and. mod(scfcv_itime,2)/=1 .and. dtset%iscf>=0 ) moved_atm_inside=1
1629 47308 : 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 47308 : 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 47308 : computed_forces=0
1639 47308 : if ((dtset%optforces==1 .and. dtset%usewvl == 0).or.(moved_atm_inside==1)) computed_forces=1
1640 47308 : if (abs(tollist(3))>tiny(0._dp)) computed_forces=1
1641 47308 : if (dtset%iscf<0) computed_forces=0
1642 47308 : if ((istep==1).and.(dtset%optforces/=1)) then
1643 3852 : 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 3852 : 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 47308 : 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 47308 : call timab(1450,2,tsec)
1667 :
1668 : ! ######################################################################
1669 : ! Compute the density rho from the trial potential
1670 : ! ----------------------------------------------------------------------
1671 47308 : call timab(1451,1,tsec)
1672 : ! Compute the density from the trial potential
1673 47308 : 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 46993 : & 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 47308 : 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 47308 : results_gs%fermie = energies%e_fermie
1716 :
1717 47308 : 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 47308 : 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 47308 : if (dtset%iscf<0) exit
1734 :
1735 : ! ######################################################################
1736 : ! In case of density mixing or wavelet handling, compute the total energy
1737 : ! ----------------------------------------------------------------------
1738 :
1739 46812 : call timab(1452,1,tsec)
1740 46812 : if (dtset%iscf>=10 .or. wvlbigdft) then
1741 14329 : optene = 1 ! use double counting scheme (default)
1742 14329 : if (wvlbigdft.and.dtset%iscf==0) optene = 0 ! use direct scheme
1743 14329 : if (dtset%iscf==22) optene = -1
1744 :
1745 : ! Add the Fock contribution to E_xc and E_xcdc if required
1746 14329 : 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 14329 : 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 14329 : & 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 46812 : call timab(1452,2,tsec)
1773 :
1774 : ! ######################################################################
1775 : ! In case of density mixing, check the exit criterion
1776 : ! ----------------------------------------------------------------------
1777 46812 : if (dtset%iscf>=10.or.(wvlbigdft.and.dtset%iscf>0)) then
1778 : ! Check exit criteria
1779 14329 : call timab(1453,1,tsec)
1780 14329 : choice=2
1781 14329 : if(paw_dmft%use_dmft==1) then
1782 113 : call prtene(dtset,energies,std_out,psps%usepaw)
1783 : end if
1784 14329 : 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 14329 : & electronpositron=electronpositron,fock=fock,rcpaw=rcpaw)
1804 14329 : call timab(1453,2,tsec)
1805 :
1806 : ! Check if we need to exit the loop
1807 14329 : call timab(1454,1,tsec)
1808 14329 : 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 14329 : if(dtset%userec==1.and.rec_set%quitrec==2)quit=1
1812 14329 : if (istep==nstep) quit=1
1813 14329 : quit_sum=quit
1814 14329 : call xmpi_sum(quit_sum,spaceComm,ierr)
1815 14329 : if (quit_sum>0) quit=1
1816 14329 : call timab(1454,2,tsec)
1817 :
1818 : ! If criteria in scprqt say to quit, then exit the loop over istep.
1819 14329 : if (quit==1) exit
1820 : end if
1821 :
1822 :
1823 : ! ######################################################################
1824 : ! Mix the total density (if required)
1825 : ! ----------------------------------------------------------------------
1826 45405 : call timab(1455,1,tsec)
1827 :
1828 45405 : 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 12922 : & .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 12922 : 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 12922 : 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 12922 : & mix_mgga=mix_mgga,taug=taug,taur=taur,tauresid=nvtauresid)
1891 : ABI_NVTX_END_RANGE()
1892 :
1893 : end if ! iscf>=10
1894 :
1895 45405 : call timab(1455,2,tsec)
1896 :
1897 : ! ######################################################################
1898 : ! Additional computation in case of an electric field or electric displacement field
1899 : ! ----------------------------------------------------------------------
1900 :
1901 45405 : 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 45405 : & 1,quit,istep,ucvol,unit_out,psps%usepaw,xred,ylm,ylmgr)
1909 :
1910 45405 : call timab(1456,2,tsec)
1911 :
1912 : ! ######################################################################
1913 : ! Compute the new potential from the trial density
1914 : ! ----------------------------------------------------------------------
1915 :
1916 45405 : 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 45405 : optxc=1
1921 45405 : 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 45405 : if (dtset%iscf/=22) then
1933 : ! PAW: eventually recompute compensation density (and gradients)
1934 45405 : nhatgrdim=0
1935 45405 : if ( allocated(nhatgr) ) then
1936 0 : ABI_FREE(nhatgr)
1937 : end if
1938 45405 : if (psps%usepaw==1) then
1939 13706 : ider=-1;if (dtset%iscf>=10.and.((dtset%xclevel==2.and.dtset%pawnhatxc>0).or.usexcnhat==0)) ider=0
1940 13706 : if (dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0) ider=ider+2
1941 13706 : if (ipositron==1) ider=-1
1942 13706 : if (ider>0) then
1943 741 : nhatgrdim=1
1944 3705 : ABI_MALLOC(nhatgr,(nfftf,dtset%nspden,3))
1945 : else
1946 12965 : ABI_MALLOC(nhatgr,(0,0,0))
1947 : end if
1948 13706 : if (ider>=0) then
1949 : ABI_NVTX_START_RANGE(NVTX_SCFCV_PAWKNHAT)
1950 8766 : call timab(558,1,tsec)
1951 8766 : 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 8766 : & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
1959 :
1960 8766 : 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 45405 : 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 45405 : 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 45405 : & electronpositron=electronpositron,vxc_hybcomp=vxc_hybcomp,add_tfw=tfw_activated)
2005 : ABI_NVTX_END_RANGE()
2006 :
2007 45405 : 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 45405 : call timab(1457,2,tsec)
2020 45405 : call timab(1452,1,tsec)
2021 :
2022 : ! This is inside the loop, its not equivalent to the line 1821
2023 46829 : if(moved_atm_inside==1) xred_old(:,:)=xred(:,:)
2024 :
2025 45405 : 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 45405 : call timab(1452,2,tsec)
2068 :
2069 : ! ######################################################################
2070 : ! Check exit criteria in case of potential mixing or direct minimization
2071 : ! ----------------------------------------------------------------------
2072 45405 : 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 41202 : call timab(1458,1,tsec)
2133 41202 : 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 41202 : initialized=1
2162 :
2163 : ! This is to save the density for restart.
2164 41202 : if (iwrite_fftdatar(mpi_enreg)) then
2165 :
2166 32929 : 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 32929 : 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 32929 : 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 41202 : ABI_FREE(nhatgr)
2209 :
2210 41202 : istep_mix=istep_mix+1
2211 41202 : if (reset_mixing) then
2212 1 : istep_mix=1;reset_mixing=.false.
2213 : end if
2214 41202 : if (ipositron/=0) electronpositron%istep_scf=electronpositron%istep_scf+1
2215 :
2216 43917 : call timab(1458,2,tsec)
2217 : end do ! istep
2218 : ABI_NVTX_END_RANGE()
2219 :
2220 6918 : ABI_FREE(rmm_diis_status)
2221 6918 : ABI_SFREE(nhatgr)
2222 :
2223 : if((dtset%wfoptalg == 1 .or. dtset%wfoptalg == 111 .or. dtset%wfoptalg == 112) .and. psps%usepaw == 1 &
2224 6918 : & .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 6918 : call xmpi_wait(quitsum_request,ierr)
2230 6918 : if (timelimit_exit == 1) istep = istep - 1
2231 :
2232 6918 : call timab(1459,1,tsec)
2233 :
2234 6918 : if (dtset%iscf > 0) then
2235 6407 : call abi_mixing_deallocate(mix)
2236 6407 : if (dtset%usekden/=0) call abi_mixing_deallocate(mix_mgga)
2237 : end if
2238 :
2239 6918 : if (dtset%usefock==1)then
2240 123 : if(wfmixalg/=0) call scf_history_free(scf_history_wf)
2241 : end if
2242 :
2243 6918 : if (quit==1.and.nstep==1) initialized=1
2244 :
2245 : !######################################################################
2246 : !Case nstep==0: compute energy based on incoming wf
2247 : !----------------------------------------------------------------------
2248 :
2249 6918 : 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 6918 : 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 6918 : 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 6918 : & dtset%userid .EQ. 1 )
2309 :
2310 88466 : if( (psps%usepaw.EQ.1) .AND. ANY(ABS(dtset%nucdipmom)>tol8) ) recompute_cprj=.TRUE.
2311 6918 : if(dtset%berryopt == -2 .AND. dtset%orbmag /= 0) recompute_cprj=.TRUE.
2312 :
2313 6918 : if (recompute_cprj) then
2314 62 : usecprj=1
2315 62 : mband_cprj=dtset%mband/mpi_enreg%nproc_band
2316 62 : mcprj=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
2317 42284 : ABI_MALLOC(cprj_local,(dtset%natom,mcprj))
2318 62 : ncpgr = 0 ; ctocprj_choice = 1
2319 62 : 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 62 : if (dtset%berryopt == -2 .AND. dtset%orbmag /= 0) then
2329 1 : ncpgr=3; ctocprj_choice=5
2330 : end if
2331 62 : call pawcprj_alloc(cprj_local,ncpgr,dimcprj)
2332 62 : cprj=> cprj_local
2333 62 : 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 62 : & dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr)
2341 : end if
2342 :
2343 : ! MRM print final Hartree energy components
2344 6918 : 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 6918 : &+energies%e_elecfield+energies%e_magfield+energies%paw%epaw
2349 :
2350 6918 : if (me == 0) then
2351 5444 : write(std_out,'(a1)')' '
2352 5444 : write(std_out,'(a98)')'-------------------------------------------------------------------------------------------------'
2353 5444 : write(std_out,'(a,2(es16.6,a))')' Ekinetic = : ',energies%e_kinetic ,' Ha ,',energies%e_kinetic*Ha_eV ,' eV'
2354 5444 : write(std_out,'(a,2(es16.6,a))')' Evext_l = : ',energies%e_localpsp ,' Ha ,',energies%e_localpsp*Ha_eV ,' eV'
2355 5444 : 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 5444 : write(std_out,'(a,2(es16.6,a))')' Epsp_core = : ',energies%e_corepsp ,' Ha ,',energies%e_corepsp*Ha_eV ,' eV'
2359 5444 : write(std_out,'(a,2(es16.6,a))')' Ehartree = : ',energies%e_hartree ,' Ha ,',energies%e_hartree*Ha_eV ,' eV'
2360 5444 : 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 5444 : write(std_out,'(a,2(es16.6,a))')' Exc_ks = : ',energies%e_xc ,' Ha ,',energies%e_xc*Ha_eV ,' eV'
2364 5444 : 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 5444 : 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 5444 : 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 5444 : if(abs(energies%paw%epaw)>1.0d-6) then
2374 1031 : write(std_out,'(a,2(es16.6,a))')' Epaw = : ',energies%paw%epaw ,' Ha ,',energies%paw%epaw*Ha_eV ,' eV'
2375 : endif
2376 5444 : write(std_out,'(a,2(es16.6,a))')' Enn = : ',energies%e_ewald ,' Ha ,',energies%e_ewald*Ha_eV ,' eV'
2377 5444 : write(std_out,'(a98)')'-------------------------------------------------------------------------------------------------'
2378 5444 : write(std_out,'(a,2(es16.6,a))')' Etot = : ',esum ,' Ha ,',esum*Ha_eV ,' eV'
2379 5444 : write(std_out,'(a98)')'-------------------------------------------------------------------------------------------------'
2380 : end if ! end MRM printing energy components
2381 :
2382 6918 : call timab(1459,2,tsec)
2383 6918 : 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 20754 : & 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 88262 : xred_old(:,:)=xred(:,:)
2403 :
2404 6918 : call timab(1460,2,tsec)
2405 :
2406 : !######################################################################
2407 : !All calculations in scfcv_core are finished. Printing section
2408 : !----------------------------------------------------------------------
2409 :
2410 6918 : 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 6918 : & rhor,rprimd,taur,ucvol,usecprj,vhartr,vpsp,vtrial,vxc,wvl%den,xccc3d,xred,rcpaw)
2420 :
2421 6918 : call timab(1461,2,tsec)
2422 6918 : 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 6918 : if (psps%usepaw==1.and.usecprj==1) then
2442 140 : iorder_cprj=0
2443 140 : call pawcprj_reorder(cprj,atindx)
2444 140 : 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 6918 : call prc_mem_free()
2455 :
2456 6918 : ABI_FREE(fcart)
2457 6918 : ABI_FREE(gred)
2458 6918 : ABI_FREE(forold)
2459 6918 : ABI_FREE(grchempottn)
2460 6918 : ABI_FREE(grcondft)
2461 6918 : ABI_FREE(gresid)
2462 6918 : ABI_FREE(grewtn)
2463 6918 : ABI_FREE(grnl)
2464 6918 : ABI_FREE(grvdw)
2465 6918 : ABI_FREE(grxc)
2466 6918 : ABI_FREE(intgres)
2467 6918 : ABI_FREE(synlgr)
2468 6918 : ABI_FREE(ph1d)
2469 6918 : ABI_FREE(ph1df)
2470 6918 : ABI_FREE(vhartr)
2471 6918 : ABI_FREE(vtrial)
2472 6918 : ABI_FREE(vpsp)
2473 6918 : ABI_FREE(vxc)
2474 6918 : ABI_FREE(vxc_hybcomp)
2475 6918 : ABI_FREE(vxctau)
2476 6918 : ABI_FREE(xccc3d)
2477 6918 : ABI_FREE(kxc)
2478 6918 : ABI_FREE(shiftvector)
2479 6918 : ABI_FREE(dtn_pc)
2480 6918 : ABI_FREE(grhf)
2481 6918 : ABI_FREE(nvresid)
2482 6918 : ABI_FREE(nvtauresid)
2483 6918 : ABI_FREE(intgden)
2484 6918 : ABI_FREE(intgden0)
2485 :
2486 : !Deallocate precon-object
2487 : !call precon%save(ngfft, 1) !DEBUG
2488 6918 : call precon%free()
2489 :
2490 6918 : if(allocated(vectornd)) then
2491 6918 : ABI_FREE(vectornd)
2492 : end if
2493 :
2494 6918 : if(associated(rcpaw)) then
2495 3 : call rcpaw_reinit(rcpaw)
2496 : endif
2497 :
2498 6918 : if((nstep>0.and.dtset%iscf>0).or.dtset%iscf==-1) then
2499 6409 : ABI_FREE(dielinv)
2500 : end if
2501 6918 : ABI_FREE(gbound_diel)
2502 6918 : ABI_FREE(irrzondiel)
2503 6918 : ABI_FREE(kg_diel)
2504 6918 : ABI_FREE(phnonsdiel)
2505 6918 : ABI_FREE(susmat)
2506 6918 : ABI_FREE(ph1ddiel)
2507 6918 : ABI_FREE(ylmdiel)
2508 :
2509 6918 : if (psps%usepaw==1) then
2510 1500 : if (dtset%iscf>0) then
2511 4487 : do iatom=1,my_natom
2512 3010 : pawrhoij(iatom)%lmnmix_sz=0
2513 3010 : pawrhoij(iatom)%use_rhoijres=0
2514 3010 : ABI_FREE(pawrhoij(iatom)%kpawmix)
2515 4487 : ABI_FREE(pawrhoij(iatom)%rhoijres)
2516 : end do
2517 : end if
2518 : ! if (recompute_cprj.or.usecprj==1) then
2519 1500 : if (recompute_cprj) then
2520 62 : usecprj=0;mcprj=0
2521 62 : call pawcprj_free(cprj)
2522 25788 : ABI_FREE(cprj_local)
2523 : end if
2524 1500 : call paw_an_free(paw_an)
2525 1500 : call paw_ij_free(paw_ij)
2526 1500 : call pawfgrtab_free(pawfgrtab)
2527 1500 : 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 9973 : ABI_FREE(pawfgrtab)
2536 9973 : ABI_FREE(paw_an)
2537 9973 : ABI_FREE(paw_ij)
2538 6918 : ABI_FREE(nhat)
2539 6918 : ABI_FREE(xcctau3d)
2540 6918 : ABI_FREE(dimcprj_srt)
2541 6918 : ABI_FREE(dimcprj)
2542 :
2543 :
2544 : ! Deallocate exact exchange data at the end of the calculation
2545 6918 : 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 6918 : 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 6918 : call constrained_dft_free(constrained_dft)
2563 :
2564 6918 : call timab(1462,2,tsec)
2565 6918 : call timab(1440,2,tsec)
2566 :
2567 : DBG_EXIT("COLL")
2568 :
2569 20754 : 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 46812 : subroutine etotfor(atindx1,deltae,diffor,dtefield,dtset,&
2717 : & elast,electronpositron,energies,&
2718 46812 : & etotal,favg,fcart,fock,forold,gred,gmet,grchempottn,grcondft,gresid,grewtn,grhf,grnl,grvdw,&
2719 46812 : & grxc,gsqcut,extfpmd,indsym,kxc,maxfor,mgfft,mpi_enreg,my_natom,nattyp,&
2720 46812 : & nfft,ngfft,ngrvdw,nhat,nkxc,ntypat,nvresid,n1xccc,n3xccc,optene,optforces,optres,&
2721 46812 : & pawang,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,red_ptot,psps,rhog,rhor,rmet,rprimd,&
2722 46812 : & 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 46812 : real(dp),allocatable :: vlocal(:,:)
2781 46812 : real(dp), ABI_CONTIGUOUS pointer :: resid(:,:)
2782 :
2783 : ! *********************************************************************
2784 :
2785 46812 : call timab(80,1,tsec)
2786 :
2787 46812 : ipositron=electronpositron_calctype(electronpositron)
2788 :
2789 46812 : if (optene>-1) then
2790 :
2791 46812 : 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 46812 : energies%e_elecfield=zero
2796 46812 : 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 46812 : 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 46812 : 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 46812 : 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 46812 : 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 46812 : 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 14329 : & + energies%e_hybcomp_E0 - energies%e_hybcomp_v0 + energies%e_constrained_dft
2888 14329 : etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
2889 14329 : if (usepaw/=0) etotal = etotal + energies%paw%epaw_dc
2890 14329 : if (dtset%usedmft/=0) etotal = etotal + energies%e_hu - energies%e_dc
2891 : end if
2892 :
2893 : ! Additional stuff for electron-positron
2894 46812 : 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 46812 : 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 46812 : 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 46812 : deltae=etotal-elast
2935 46812 : elast=etotal
2936 : end if !optene/=-1
2937 :
2938 46812 : call timab(80,2,tsec)
2939 :
2940 : !------Compute forces-----------------------------------------------------
2941 :
2942 46812 : 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 345552 : fcart=zero
3006 345552 : gred=zero
3007 26320 : favg=zero
3008 26320 : diffor=zero
3009 345552 : gresid=zero
3010 345552 : grhf=zero
3011 26320 : maxfor=zero
3012 345552 : synlgr=zero
3013 : end if
3014 :
3015 46812 : call timab(80,2,tsec)
3016 :
3017 46812 : 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 : !!***
|
