Branch data Line data Source code
1 : : !!****m* ABINIT/m_afterscfloop
2 : : !! NAME
3 : : !! m_afterscfloop
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2008-2026 ABINIT group (XG)
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 : : module m_afterscfloop
23 : :
24 : : use defs_basis
25 : : use defs_wvltypes
26 : : use m_energies
27 : : use m_errors
28 : : use m_abicore
29 : : use m_ebands
30 : : use m_efield
31 : : use m_abi_mixing
32 : : use m_hdr
33 : : use m_dtset
34 : : use m_dtfil
35 : : use m_extfpmd
36 : :
37 : : use defs_datatypes, only : pseudopotential_type
38 : : use defs_abitypes, only : mpi_type
39 : : use m_time, only : timab
40 : : use m_xmpi, only : xmpi_sum, xmpi_comm_rank,xmpi_comm_size
41 : : use m_berryphase_new, only : berryphase_new
42 : : use m_geometry, only : xred2xcart, metric
43 : : use m_crystal, only : crystal_t,prtposcar
44 : : use m_results_gs , only : results_gs_type
45 : : use m_electronpositron, only : electronpositron_type, electronpositron_calctype, exchange_electronpositron
46 : : use m_paw_dmft, only : paw_dmft_type
47 : : use m_pawang, only : pawang_type
48 : : use m_pawrad, only : pawrad_type
49 : : use m_pawtab, only : pawtab_type
50 : : use m_pawrhoij, only : pawrhoij_type
51 : : use m_paw_an, only : paw_an_type
52 : : use m_paw_ij, only : paw_ij_type
53 : : use m_pawfgrtab, only : pawfgrtab_type
54 : : use m_pawcprj, only : pawcprj_type,pawcprj_getdim
55 : : use m_pawfgr, only : pawfgr_type
56 : : use m_paw_mkrho, only : pawmkrho
57 : : use m_paw_nhat, only : nhatgrid,wvl_nhatgrid
58 : : use m_paw_occupancies, only : pawmkrhoij
59 : : use m_paw_correlations, only : setnoccmmp
60 : : use m_fock, only : fock_type
61 : : use m_kg, only : getph
62 : : use m_spin_current, only : spin_current
63 : : use m_mkrho, only : mkrho, prtrhomxmn
64 : : use m_elpolariz, only : elpolariz
65 : : use m_orbmag, only : orbmag
66 : : use m_nonlop_test, only : nonlop_test
67 : : use m_common, only : scprqt
68 : : use m_xctk, only : xcden
69 : : use m_forstr, only : forstr
70 : : use m_wvl_rho, only : wvl_mkrho
71 : : use m_wvl_psi, only : wvl_psitohpsi, wvl_tail_corrections
72 : : use m_fourier_interpol, only : transgrid
73 : : use m_xg_nonlop, only : xg_nonlop_t
74 : :
75 : : #ifdef HAVE_BIGDFT
76 : : use m_abi2big
77 : : use BigDFT_API, only : last_orthon, &
78 : : & kswfn_free_scf_data, denspot_free_history,&
79 : : & write_energies, total_energies, XC_potential,&
80 : : & eigensystem_info, applyprojectorsonthefly
81 : : #endif
82 : :
83 : : implicit none
84 : :
85 : : private
86 : : !!***
87 : :
88 : : public :: afterscfloop
89 : : !!***
90 : :
91 : : contains
92 : : !!***
93 : :
94 : : !!****f* ABINIT/afterscfloop
95 : : !! NAME
96 : : !! afterscfloop
97 : : !!
98 : : !! FUNCTION
99 : : !! Perform all calculations needed after the SCF loop, independent of the
100 : : !! call to scfcv (with or without atomic displacements), and exclusive
101 : : !! of print or write purposes, or deallocations.
102 : : !!
103 : : !! INPUTS
104 : : !! atindx(natom)=index table for atoms (see gstate.f)
105 : : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
106 : : !! cg(2,mcg)=wavefunctions (may be read from disk instead of input)
107 : : !! cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn>
108 : : !! cpus= cpu time limit in seconds
109 : : !! deltae=change in energy between the previous and present SCF cycle
110 : : !! dtfil <type(datafiles_type)>=variables related to files
111 : : !! dtset <type(dataset_type)>=all input variables in this dataset
112 : : !! | mband=maximum number of bands
113 : : !! | mgfft=maximum size of 1D FFTs (see NOTES at beginning of scfcv)
114 : : !! | mkmem=maximum number of k points in core memory
115 : : !! | mpw = maximum number of plane waves
116 : : !! | natom=number of atoms in cell
117 : : !! | nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv)
118 : : !! | nkpt=number of k points in Brillouin zone
119 : : !! | nspden=number of spin-density components
120 : : !! | nsppol=1 for unpolarized, 2 for spin-polarized
121 : : !! | nsym=number of symmetries in space group
122 : : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
123 : : !! extfpmd <type(extfpmd_type)>=extended first-principles molecular dynamics type
124 : : !! fock <type(fock_type)>= quantities to calculate Fock exact exchange
125 : : !! grchempottn(3,natom)=d(E_chemical_potential)/d(xred) (hartree)
126 : : !! grcondft(3,natom)=d(E_constrained_DFT)/d(xred) (hartree)
127 : : !! grewtn(3,natom)=d(Ewald)/d(xred) (hartree)
128 : : !! grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D2 dispersion (hartree)
129 : : !! gsqcut=cutoff value on G**2 for (large) sphere inside FFT box.
130 : : !! gsqcut=(boxcut**2)*dtset%ecut/(2._dp*(Pi**2)
131 : : !! hdr <type(hdr_type)>=the header of wf, den and pot files
132 : : !! indsym(4,nsym,natom)=index showing transformation of atom labels
133 : : !! under symmetry operations (computed in symatm)
134 : : !! intgres(nspden,natom)=integrated residuals from constrained DFT. They are also Lagrange parameters, or gradients with respect to constraints.
135 : : !! irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
136 : : !! istep=number of the SCF iteration
137 : : !! istep_fock_outer=number of outer SCF iteration in the double loop approach
138 : : !! istep_mix=number of inner SCF iteration in the double loop approach
139 : : !! kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere
140 : : !! kxc(nfftf,nkxc)=XC kernel
141 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
142 : : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
143 : : !! mgfftf= - PAW only - maximum size of 1D FFTs for the "fine" grid (see NOTES at beginning of scfcv)
144 : : !! ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
145 : : !! moved_atm_inside: if==1, the atoms are allowed to move.
146 : : !! mpi_enreg=information about MPI parallelization
147 : : !! my_natom=number of atoms treated by current processor
148 : : !! n3xccc=dimension of the xccc3d array (0 or nfft).
149 : : !! nattyp(dtset%ntypat)=number of atoms of each type
150 : : !! nfftf= - PAW only - number of FFT grid points for the "fine" grid (see NOTES at beginning of scfcv)
151 : : !! ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
152 : : !! ngfftf(18)= - PAW only - contain all needed information about 3D FFT for the "fine" grid
153 : : !! nhat(nfftf,nspden*psps%usepaw)= -PAW only- compensation density
154 : : !! nkxc=dimension of kxc
155 : : !! npwarr(nkpt)=number of planewaves in basis and on boundary for each k
156 : : !! nvresid(nfftf,nspden)=array for the residual of the density/potential
157 : : !! occ(mband*nkpt*nsppol)=occupancies of bands at various k points
158 : : !! optres=0: the potential residual has been computed in scfcv
159 : : !! 1: the density residual has been computed in scfcv
160 : : !! paw_an(my_natom*usepaw) <type(paw_an_type)>=paw arrays given on angular mesh
161 : : !! paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
162 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
163 : : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
164 : : !! pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
165 : : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
166 : : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
167 : : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
168 : : !! pel(3)=reduced coordinates of the electronic polarization (a. u.)
169 : : !! pel_cg(3) = reduced coordinates of the electronic polarization (a. u.)
170 : : !! computed in the SCF loop
171 : : !! ph1df(2,3*(2*mgfftf+1)*natom)= - PAW only - 1-dim structure factor phases for the "fine" grid
172 : : !! (see NOTES at beginning of scfcv)
173 : : !! phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
174 : : !! pion(3)=reduced coordinates of the ionic polarization (a. u.)
175 : : !! prtfor=1 only if forces have to be printed (0 otherwise)
176 : : !! prtxml=1 if XML file has to be output
177 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
178 : : !! | mpsang= 1+maximum angular momentum
179 : : !! pwind(pwind_alloc,2,3) = array used to compute
180 : : !! the overlap matrix smat between k-points (see initberry.f)
181 : : !! pwind_alloc = first dimension of pwind
182 : : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
183 : : !! res2=density/potential residual (squared)
184 : : !! resid(mband*nkpt*nsppol)=residuals for each band over all k points and spins
185 : : !! residm=maximum value from resid array (except for nbdbuf highest bands)
186 : : !! rhog(2,nfftf)=Fourier transform of total electron density (including compensation density in PAW)
187 : : !! rhor(nfftf,nspden)=total electron density (including compensation density in PAW)
188 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
189 : : !! stress_needed=1 if stresses are needed, 0 otherwise
190 : : !! strscondft(6)=cDFT correction to stress
191 : : !! strsxc(6)=xc correction to stress
192 : : !! symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
193 : : !! tollist(12)=list of tolerances
194 : : !! usecprj=1 if cprj datastructure has been allocated
195 : : !! usevxctau=1 if kinetic energy density contribution has to be included (mGGA)
196 : : !! vhartr(nfftf)=Hartree potential
197 : : !! vpsp(nfftf)=array for holding local psp
198 : : !! vxc(nfftf,nspden)=exchange-correlation potential (hartree) in real space
199 : : !! vxctau(nfft,nspden,4*usevxctau)=(only for meta-GGA) derivative of XC energy density
200 : : !! wrt kinetic energy density (depsxcdtau)
201 : : !! vxcavg=vxc average
202 : : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
203 : : !! xcctau3d(n3xccc*usekden)=(only for meta-GGA): 3D core electron kinetic energy density for XC core correction
204 : : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
205 : : !! ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
206 : : !!
207 : : !! OUTPUT
208 : : !! conv_retcode= Non-zero if convergence is not achieved.
209 : : !! elfr(nfftf,nspden)=electron localization function
210 : : !! grhor(nfft,nspden,3)= gradient of electron density in electrons/bohr**4, real space
211 : : !! lrhor(nfft,nspden)= Laplacian of electron density in electrons/bohr**5, real space
212 : : !! results_gs <type(results_gs_type)>=results (energy and its components,
213 : : !! forces and its components, the stress tensor) of a ground-state computation
214 : : !! (should be made a pure output quantity)
215 : : !! taug(2,nfftf)=Fourier transform of total kinetic energy density
216 : : !! taur(nfftf,nspden)=total kinetic energy density in real space
217 : : !! ==== if forces are required ====
218 : : !! diffor=maximal absolute value of changes in the components of
219 : : !! force between the input and the output.
220 : : !! favg(3)=mean of the forces before correction for translational symmetry
221 : : !! fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
222 : : !! at input, previous value of forces,
223 : : !! at output, new value.
224 : : !! Note : unlike gred, this array has been corrected by enforcing
225 : : !! the translational symmetry, namely that the sum of force
226 : : !! on all atoms is zero.
227 : : !! gred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
228 : : !! gresid(3,natom)=forces due to the residual of the potential
229 : : !! grhf(3,natom)=Hellman-Feynman derivatives of the total energy
230 : : !! grxc(9+3*natom)=d(Exc)/d(xred) if core charges are used
231 : : !! maxfor=maximal absolute value of the output array force.
232 : : !! synlgr(3,natom)=symmetrized gradients of energy due to nonlocal contributions
233 : : !! ==== if stress tensor is required ====
234 : : !! strten(6)=components of the stress tensor (hartree/bohr^3) for the
235 : : !! 6 unique components of this symmetric 3x3 tensor:
236 : : !! Given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1).
237 : : !!
238 : : !! SIDE EFFECTS
239 : : !! computed_forces=1 if forces have been computed, 0 otherwise
240 : : !! dtefield <type(efield_type)> = variables related to Berry phase
241 : : !! and electric field calculations (see initberry.f).
242 : : !! In case dtset%berryopt = 4/6/7/14/16/17, the overlap matrices computed
243 : : !! in this routine are stored in dtefield%smat in order
244 : : !! to be used in the electric field calculation.
245 : : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
246 : : !! energies <type(energies_type)>=all part of total energy.
247 : : !! | entropy(IN)=entropy due to the occupation number smearing (if metal)
248 : : !! | e_localpsp(IN)=local psp energy (hartree)
249 : : !! | e_eigenvalues(IN)=Sum of the eigenvalues - Band energy (Hartree)
250 : : !! | e_chempot(IN)=energy from spatially varying chemical potential (hartree)
251 : : !! | e_ewald(IN)=Ewald energy (hartree)
252 : : !! | e_vdw_dftd(IN)=VdW DFT-D energy
253 : : !! | e_hartree(IN)=Hartree part of total energy (hartree units)
254 : : !! | e_corepsp(IN)=psp core-core energy
255 : : !! | e_hybcomp_E0=energy compensation term for hybrid exchange-correlation energy (hartree) at fixed density
256 : : !! | e_hybcomp_v0=potential compensation term for hybrid exchange-correlation energy (hartree) at fixed density
257 : : !! | e_hybcomp_v=potential compensation term for hybrid exchange-correlation energy (hartree) at self-consistent den
258 : : !! | e_kinetic(IN)=kinetic energy part of total energy.
259 : : !! | e_nlpsp_vfock(IN)=nonlocal psp + potential Fock ACE part of total energy.
260 : : !! | e_xc(IN)=exchange-correlation energy (hartree)
261 : : !! | e_xcdc(IN)=exchange-correlation double-counting energy (hartree)
262 : : !! | paw%e_paw(IN)=PAW spherical part energy
263 : : !! | paw%e_pawdc(IN)=PAW spherical part double-counting energy
264 : : !! | e_elecfield(OUT)=the term of the energy functional that depends explicitely !!HONG
265 : : !! | on the electric field:
266 : : !! | enefield = -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j for fixed E/ebar
267 : : !! | = Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j for fixed D/d
268 : : !! etotal=total energy, might be correct by improved polarization computation
269 : : !! forold(3,natom)=old forces
270 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
271 : : !! ===== if dtset%densfor_pred==3 .and. moved_atm_inside==1 =====
272 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases (coarse grid)
273 : : !! ph1df(2,3*(2*mgfftf+1)*natom)=1-dim structure factor phases (fine PAW grid)
274 : : !! wvl <type(wvl_data)>=all wavelets data.
275 : : !!
276 : : !! NOTES
277 : : !!
278 : : !! SOURCE
279 : :
280 : 6953 : subroutine afterscfloop(atindx,atindx1,cg,computed_forces,cprj,cpus,&
281 : 6953 : & deltae,diffor,difmag,dtefield,dtfil,dtset,eigen,electronpositron,elfr,&
282 : 13906 : & energies,etotal,extfpmd,favg,fcart,fock,forold,grchempottn,grcondft,&
283 : 13906 : & gred,gresid,grewtn,grhf,grhor,grvdw,&
284 : 20859 : & grxc,gsqcut,hdr,indsym,intgres,irrzon,istep,istep_fock_outer,istep_mix,&
285 : 13906 : & kg,kxc,lrhor,maxfor,maxmag,mcg,mcprj,mgfftf,&
286 : 13906 : & moved_atm_inside,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfft,ngfftf,ngrvdw,nhat,&
287 : 13906 : & nkxc,npwarr,nvresid,occ,optres,paw_an,paw_ij,pawang,pawfgr,&
288 : 6953 : & pawfgrtab,pawrad,pawrhoij,pawtab,pel,pel_cg,ph1d,ph1df,phnons,pion,prtfor,prtxml,&
289 : 13906 : & psps,pwind,pwind_alloc,pwnsfac,res2,resid,residm,results_gs,&
290 : 13906 : & rhog,rhor,rprimd,stress_needed,strscondft,strsxc,strten,symrec,synlgr,taug,&
291 : 13906 : & taur,tollist,usecprj,usevxctau,vhartr,vpsp,vtrial,vxc,vxctau,vxcavg,wvl,&
292 [ - + - + ]: 6953 : & xccc3d,xcctau3d,xred,ylm,ylmgr,qvpotzero,conv_retcode,xg_nonlop)
293 : :
294 : : !Arguments ------------------------------------
295 : : !scalars
296 : : integer,intent(in) :: istep,istep_fock_outer,istep_mix
297 : : integer,intent(in) :: mcg,mcprj,mgfftf,moved_atm_inside,my_natom,n3xccc,nfftf,ngrvdw,nkxc
298 : : integer,intent(in) :: optres,prtfor,prtxml,pwind_alloc,stress_needed,usecprj,usevxctau
299 : : integer,intent(inout) :: computed_forces
300 : : real(dp),intent(in) :: cpus,deltae,gsqcut,res2,residm
301 : : real(dp),intent(in) :: qvpotzero
302 : : real(dp),intent(inout) :: diffor,etotal,maxfor,vxcavg
303 : : type(MPI_type),intent(inout) :: mpi_enreg
304 : : type(datafiles_type),intent(in) :: dtfil
305 : : type(dataset_type),intent(inout) :: dtset
306 : : type(efield_type),intent(inout) :: dtefield
307 : : type(electronpositron_type),pointer :: electronpositron
308 : : type(energies_type),intent(inout) :: energies
309 : : type(extfpmd_type),pointer,intent(inout) :: extfpmd
310 : : type(hdr_type),intent(inout) :: hdr
311 : : type(pawang_type),intent(in) :: pawang
312 : : type(pawfgr_type),intent(in) :: pawfgr
313 : : type(pseudopotential_type),intent(in) :: psps
314 : : type(results_gs_type),intent(inout) :: results_gs
315 : : type(wvl_data),intent(inout) :: wvl
316 : : type(fock_type),pointer, intent(inout) :: fock
317 : : type(xg_nonlop_t), intent(inout) :: xg_nonlop
318 : : !arrays
319 : : integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom)
320 : : integer,intent(in) :: indsym(4,dtset%nsym,dtset%natom)
321 : : integer,intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
322 : : integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(dtset%ntypat)
323 : : integer,intent(in) :: ngfft(18),ngfftf(18),npwarr(dtset%nkpt)
324 : : integer,intent(in) :: pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym)
325 : : integer,intent(out) :: conv_retcode
326 : : real(dp),intent(in) :: grchempottn(3,dtset%natom),grewtn(3,dtset%natom),grvdw(3,ngrvdw)
327 : : real(dp),intent(in) :: grcondft(:,:) ! (3,natom) if constrainedDFT otherwise (3,0)
328 : : real(dp),intent(in) :: intgres(:,:) ! (nspden,natom) if constrainedDFT otherwise (nspden,0)
329 : : real(dp),intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
330 : : real(dp),intent(in) :: pwnsfac(2,pwind_alloc)
331 : : real(dp),intent(in) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
332 : : real(dp),intent(in) :: strscondft(6)
333 : : real(dp),intent(in) :: tollist(12),vpsp(nfftf)
334 : : real(dp),intent(inout) :: vtrial(nfftf,dtset%nspden)
335 : : real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
336 : : real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
337 : : real(dp),intent(in) :: maxmag, difmag
338 : : real(dp),intent(inout) :: cg(2,mcg)
339 : : real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
340 : : real(dp),intent(inout) :: forold(3,dtset%natom)
341 : : real(dp),intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw)
342 : : real(dp),intent(inout) :: nvresid(nfftf,dtset%nspden),pel(3)
343 : : real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),pel_cg(3)
344 : : real(dp),intent(inout) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)
345 : : real(dp),intent(inout) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom),pion(3)
346 : : real(dp),intent(inout) :: rprimd(3,3)
347 : : real(dp),intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden),strsxc(6)
348 : : real(dp),intent(inout) :: vhartr(nfftf),vxc(nfftf,dtset%nspden),vxctau(nfftf,dtset%nspden,4*usevxctau)
349 : : real(dp),intent(inout) :: xccc3d(n3xccc),xcctau3d(n3xccc*dtset%usekden),xred(3,dtset%natom)
350 : : real(dp),intent(inout) :: favg(3),fcart(3,dtset%natom),gred(3,dtset%natom)
351 : : real(dp),intent(inout) :: gresid(3,dtset%natom),grhf(3,dtset%natom)
352 : : real(dp),intent(inout) :: grxc(3,dtset%natom),kxc(nfftf,nkxc),strten(6)
353 : : real(dp),intent(inout) :: synlgr(3,dtset%natom)
354 : : real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:),taug(:,:),taur(:,:)
355 : : type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj)
356 : : type(paw_an_type),intent(inout) :: paw_an(my_natom*psps%usepaw)
357 : : type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw)
358 : : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
359 : : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
360 : : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw)
361 : : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
362 : :
363 : : !Local variables-------------------------------
364 : : !scalars
365 : : integer,parameter :: response=0
366 : : integer :: bantot,bufsz,calc_pol_ddk,choice,cplex,ierr,ifft,igrad,ishift,ispden
367 : : integer :: mcg1_3,nfftotf,ngrad,optcut,optfor,optgr0,optgr1,optgr2,optrad,quit,shft
368 : : integer :: spaceComm_fft,tim_mkrho
369 : : logical :: save_cg1_3,test_gylmgr,test_nfgd,test_rfgd
370 : : logical :: remove_inv=.false.,wvlbigdft=.false.
371 : : real(dp) :: c_fermi,dtaur,dtaurzero,ucvol
372 : : character(len=500) :: message
373 [ + + + + : 361556 : type(crystal_t) :: crystal
+ + + + +
+ + + + +
+ + + + ]
374 : 6953 : type(ebands_t) :: ebands_k
375 : 6953 : type(paw_dmft_type) :: paw_dmft
376 : : #if defined HAVE_BIGDFT
377 : : integer :: ia,ii,mband_cprj
378 : : logical :: do_last_ortho,compute_wvl_tail=.false.
379 : : real(dp) :: dum,eexctx,eh,ekin,eloc,enl,eproj,esicdc,evxc,exc,ucvol_local
380 : : #endif
381 : : !arrays
382 : : real(dp) :: gmet(3,3),gprimd(3,3),pelev(3),ptot(3),red_ptot(3),rmet(3,3),tsec(2)
383 : : real(dp) :: dmatdum(0,0,0,0)
384 : 6953 : real(dp),allocatable :: cg1_3(:,:,:),doccde(:),mpibuf(:,:),qphon(:),rhonow(:,:,:),sqnormgrhor(:,:)
385 : 6953 : real(dp),allocatable :: tauwfg(:,:),tauwfr(:,:),vtrial_local(:,:)
386 : : #if defined HAVE_BIGDFT
387 : : integer,allocatable :: dimcprj_srt(:)
388 : : real(dp),allocatable :: hpsi_tmp(:),xcart(:,:)
389 : : #endif
390 : : ! *************************************************************************
391 : :
392 : : DBG_ENTER("COLL")
393 : :
394 : 6953 : call timab(250,1,tsec)
395 : 6953 : call timab(251,1,tsec)
396 : :
397 : : !Compute different geometric tensor, as well as ucvol, from rprimd
398 : 6953 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
399 : : nfftotf=product(ngfftf(1:3))
400 : :
401 : : call crystal%init(dtset%amu_orig(:,1),dtset%spgroup,dtset%natom,dtset%npsp,&
402 : : & psps%ntypat,dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,&
403 : : & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,psps%title,&
404 [ + + + + ]: 13876 : & symrel=dtset%symrel,tnons=dtset%tnons,symafm=dtset%symafm)
405 : :
406 [ - + + - : 20859 : ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
+ - + - -
+ - + ]
407 [ + + ]: 500792 : doccde=zero
408 : : call ebands_k%init(hdr%bantot,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
409 : : & doccde,eigen,hdr%istwfk,hdr%kptns,hdr%nband,&
410 : : & hdr%nkpt,hdr%npwarr,hdr%nsppol,hdr%nspinor,hdr%tphysel,hdr%tsmear,hdr%occopt,hdr%occ,hdr%wtk,&
411 : : & hdr%cellcharge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, &
412 : 6953 : & hdr%kptrlatt, hdr%nshiftk, hdr%shiftk)
413 [ - + ]: 6953 : ABI_FREE(doccde)
414 : :
415 : : !MPI FFT communicator
416 : 6953 : spaceComm_fft=mpi_enreg%comm_fft
417 : :
418 : : !Recompute structure factor phases if atomic positions have changed
419 [ + + ]: 6953 : if (moved_atm_inside==1) then
420 [ - + - - ]: 12 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
421 : 0 : call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
422 : : else
423 [ + + + + ]: 15006 : ph1d(:,:)=ph1df(:,:)
424 : : end if
425 : : end if
426 : :
427 : : !----------------------------------------------------------------------
428 : : !Wavelet case: transform psi to KS orbitals
429 : : !----------------------------------------------------------------------
430 [ - + ]: 6953 : if (dtset%usewvl == 1) then
431 : :
432 : : ! wvlbigdft indicates that the BigDFT workflow will be followed
433 : 0 : wvlbigdft=(dtset%wvl_bigdft_comp==1)
434 : :
435 : : #if defined HAVE_BIGDFT
436 : : ! Transform to KS orbitals
437 : :
438 : : ! Need xcart
439 : : ABI_MALLOC(xcart,(3, dtset%natom))
440 : : call xred2xcart(dtset%natom, rprimd, xcart, xred)
441 : : ucvol_local=product(wvl%den%denspot%dpbox%hgrids)*real(product(wvl%den%denspot%dpbox%ndims),dp)
442 : :
443 : : ! do_last_ortho in case of direct minimization, since
444 : : ! We never diagonalized the hamiltonian and the eigenvalues are unknown.
445 : : if ( wvlbigdft) do_last_ortho=(dtset%iscf==0)
446 : : if (.not.wvlbigdft) do_last_ortho=(.false.)
447 : : if (do_last_ortho) then
448 : : call total_energies(wvl%e%energs, istep, mpi_enreg%me_wvl)
449 : : call write_energies(istep,0,wvl%e%energs,zero,zero,"FINAL")
450 : : if(.not.wvlbigdft) then
451 : : ! If ISCF>10, we exit scfcv at a place where bigdft objects
452 : : ! do not contain the KS potential. Hence, we copy vtrial to wvl%den
453 : : if(dtset%iscf>=10) call wvl_vtrial_abi2big(1,vtrial,wvl%den)
454 : : ! hpsi is lost in hpsitopsi, so we recalculate it (needed for last_orthon).
455 : : call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
456 : : & istep,1,-1,mpi_enreg%me_wvl,dtset%natom,&
457 : : & nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
458 : : & dum,.false.,evxc,wvl,wvlbigdft,xcart,strsxc)
459 : : if (dtset%iscf==0) then
460 : : energies%e_kinetic=ekin ; energies%e_hartree=eh
461 : : energies%e_xc=exc ; energies%e_localpsp=eloc
462 : : energies%e_nlpsp_vfock=enl ; energies%e_exactX=eexctx
463 : : energies%e_sicdc=esicdc ; energies%e_xcdc=evxc
464 : : energies%e_eigenvalues = energies%e_kinetic + energies%e_localpsp &
465 : : & + energies%e_xcdc + two*energies%e_hartree +energies%e_nlpsp_vfock
466 : : end if
467 : : end if
468 : : call last_orthon(mpi_enreg%me_wvl,mpi_enreg%nproc_wvl,istep,wvl%wfs%ks,wvl%e%energs%evsum,.true.)
469 : : if (mpi_enreg%nproc_wvl == 1) nullify(wvl%wfs%ks%psit)
470 : : call eigensystem_info(mpi_enreg%me_wvl,mpi_enreg%nproc_wvl,0.d0,&
471 : : wvl%wfs%ks%Lzd%Glr%wfd%nvctr_c+7*wvl%wfs%ks%Lzd%Glr%wfd%nvctr_f,&
472 : : wvl%wfs%ks%orbs,wvl%wfs%ks%psi)
473 : : ! Copy eigenvalues from BigDFT object to "eigen"
474 : : call wvl_eigen_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,eigen,2,wvl%wfs)
475 : : ! Copy occupations from BigDFT objects to ABINIT
476 : : call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,2,wvl%wfs)
477 : : end if
478 : :
479 : : ! Tail corrections, pending for wvlbigdft==.false.
480 : : ! TODO put it at the end of gstate.
481 : : ! WVL - maybe compute the tail corrections to energy
482 : : compute_wvl_tail=(dtset%tl_radius>tol12.and.wvlbigdft)
483 : : if (compute_wvl_tail) then
484 : : ! Use the tails to improve energy precision.
485 : : call wvl_tail_corrections(dtset, energies, etotal, mpi_enreg, psps, wvl, xcart)
486 : : end if
487 : :
488 : : ! Clean KSwfn parts only needed in the SCF loop.
489 : : call kswfn_free_scf_data(wvl%wfs%ks, (mpi_enreg%nproc_wvl > 1))
490 : : ! Clean denspot parts only needed in the SCF loop.
491 : : call denspot_free_history(wvl%den%denspot)
492 : :
493 : : ! If WF have been modified, change the density according to the KS projection.
494 : : if ( do_last_ortho ) then
495 : :
496 : : ! Density from new orthogonalized WFs
497 : : call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl%wfs, wvl%den)
498 : :
499 : : ! PAW: has to update cprj, rhoij and compensation charge density
500 : : if (psps%usepaw==1) then
501 : : ! 1-Compute cprj
502 : : ABI_MALLOC(hpsi_tmp,(size(wvl%wfs%ks%hpsi)))
503 : : call applyprojectorsonthefly(mpi_enreg%me_wvl,wvl%wfs%ks%orbs,wvl%descr%atoms,wvl%descr%Glr,&
504 : : & xcart,wvl%descr%h(1),wvl%descr%h(2),wvl%descr%h(3),wvl%wfs%ks%lzd%Glr%wfd,&
505 : : & wvl%projectors%nlpsp,wvl%wfs%ks%psi,hpsi_tmp,eproj,&
506 : : & proj_G=wvl%projectors%G,paw=wvl%descr%paw)
507 : : ABI_FREE(hpsi_tmp)
508 : : do ii=1,mcprj
509 : : do ia=1,dtset%natom
510 : : !Note that cprj should be allocated (i.e. usepcrj=1 imposed in scfcv)
511 : : cprj(ia,ii)%cp(:,:)= wvl%descr%paw%cprj(ia,ii)%cp(:,:)
512 : : end do
513 : : end do
514 : : ! 2-Compute rhoij
515 : : ABI_MALLOC(dimcprj_srt,(dtset%natom))
516 : : call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
517 : : mband_cprj=mcprj/(dtset%nspinor*dtset%mkmem*dtset%nsppol)
518 : : paw_dmft%use_sc_dmft=0 ; paw_dmft%use_dmft=0 ! dmft not used here
519 : : call pawmkrhoij(atindx,atindx1,cprj,dimcprj_srt,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,&
520 : : & mcprj,dtset%mkmem,mpi_enreg,dtset%natom,dtset%nband,dtset%nkpt,dtset%nspden,dtset%nspinor,&
521 : : & dtset%nsppol,occ,mpi_enreg%paral_kgb,paw_dmft,pawrhoij,dtfil%unpaw,dtset%usewvl,dtset%wtk)
522 : : ABI_FREE(dimcprj_srt)
523 : : ! 3-Symetrize rhoij, compute nhat and add it to rhor
524 : : call pawmkrho(1,dum,1,gprimd,0,indsym,0,mpi_enreg,my_natom,dtset%natom,dtset%nspden,dtset%nsym,&
525 : : & dtset%ntypat,mpi_enreg%paral_kgb,pawang,pawfgr,pawfgrtab,dtset%pawprtvol,pawrhoij,pawrhoij,&
526 : : & pawtab,(/zero,zero,zero/),rhog,rhor,rhor,rprimd,dtset%symafm,symrec,dtset%typat,ucvol_local,&
527 : : & dtset%usewvl,xred,pawnhat=nhat)
528 : : call wvl_rho_abi2big(1,rhor,wvl%den)
529 : : end if
530 : : end if
531 : :
532 : : ABI_FREE(xcart)
533 : :
534 : : #else
535 : 0 : BIGDFT_NOTENABLED_ERROR()
536 : : #endif
537 : : end if
538 : :
539 : 6953 : call timab(251,2,tsec)
540 : 6953 : call timab(252,1,tsec)
541 : :
542 : : !----------------------------------------------------------------------
543 : : !Polarization Calculation, but not orbital magnetism
544 : : !----------------------------------------------------------------------
545 : :
546 [ + + + + ]: 6953 : if(dtset%berryopt/=0 .AND. dtset%orbmag == 0)then
547 : : call elpolariz(atindx1,cg,cprj,dtefield,dtfil,dtset,etotal,energies%e_elecfield,gprimd,hdr,&
548 : : & kg,dtset%mband,mcg,mcprj,dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,nattyp,dtset%nkpt,&
549 : : & npwarr,dtset%nsppol,psps%ntypat,pawrhoij,pawtab,pel,pel_cg,pelev,pion,&
550 : 85 : & psps,pwind,pwind_alloc,pwnsfac,rprimd,ucvol,usecprj,xred)
551 : : end if
552 : :
553 : 6953 : call timab(252,2,tsec)
554 : 6953 : call timab(253,1,tsec)
555 : :
556 : : !----------------------------------------------------------------------
557 : : !Orbital magnetization calculation using PEAD DDK wavefunctions
558 : : !----------------------------------------------------------------------
559 : :
560 [ + + + + ]: 6953 : if (dtset%berryopt == -2 .AND. dtset%orbmag /= 0) then
561 : 1 : save_cg1_3 = .TRUE.
562 : 1 : mcg1_3 = mcg
563 [ - + + - : 4 : ABI_MALLOC(cg1_3,(2,mcg1_3,3))
+ - + - +
- - + -
+ ]
564 : :
565 : 1 : calc_pol_ddk = 2
566 : : call berryphase_new(atindx1,cg,cg1_3,cprj,dtefield,dtfil,dtset,psps,&
567 : : & gprimd,hdr,psps%indlmn,kg,psps%lmnmax,dtset%mband,mcg,mcg1_3,mcprj,&
568 : : & dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,npwarr,dtset%nsppol,psps%ntypat,&
569 : : & dtset%nkpt,calc_pol_ddk,pawrhoij,pawtab,pel,pelev,pion,ptot,red_ptot,pwind,& !!REC
570 : : & pwind_alloc,pwnsfac,rprimd,save_cg1_3,dtset%typat,ucvol,ab_out,&
571 : 1 : & usecprj,psps%usepaw,xred,psps%ziontypat)
572 : :
573 : : if ( .NOT. ALLOCATED(vtrial_local)) then
574 [ - + + - : 4 : ABI_MALLOC(vtrial_local,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
575 : : end if
576 [ + - - + : 13827 : vtrial_local = vtrial
- - + + +
+ ]
577 : : call orbmag(cg,cg1_3,cprj,crystal,dtfil,dtset,ebands_k,gsqcut,hdr,kg,mcg,mcg1_3,&
578 : : & mcprj,dtset%mkmem,mpi_enreg,dtset%mpw,nfftf,ngfftf,paw_ij,pawfgr,&
579 : 1 : & pawrad,pawtab,psps,usevxctau,vtrial_local,vxctau,ylm,ylmgr)
580 : :
581 [ - + ]: 1 : ABI_FREE(vtrial_local)
582 [ - + ]: 1 : ABI_FREE(cg1_3)
583 : :
584 : : end if
585 : :
586 : : !----------------------------------------------------------------------
587 : : !Gradient and Laplacian of the Density Calculation
588 : : !----------------------------------------------------------------------
589 : :
590 : : !We use routine xcden which get gradient of rhor (grhor), and eventually laplacian of rhor (lrhor).
591 [ + + - + ]: 6953 : if(dtset%prtgden/=0 .or. dtset%prtlden/=0)then
592 : :
593 : : ! Compute gradient of the electron density
594 : 3 : ngrad=2
595 : 3 : cplex=1
596 : 3 : ishift=0
597 [ - + + - : 15 : ABI_MALLOC(rhonow,(nfftf,dtset%nspden,ngrad*ngrad))
- + + - +
- + - + -
- + - + ]
598 [ + + ]: 3 : if(dtset%prtlden/=0)then
599 : 2 : nullify(lrhor)
600 [ - + + - : 8 : ABI_MALLOC(lrhor,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
601 : : end if
602 : 3 : write(message,'(a,a)') ch10, " Compute gradient of the electron density"
603 : 3 : call wrtout(ab_out,message)
604 [ + + ]: 3 : if(dtset%prtlden/=0) then
605 : 2 : write(message,'(a)') " and also Compute Laplacian of the electron density"
606 : 2 : call wrtout(ab_out,message)
607 : : end if
608 : 3 : write(message,'(a)') "--------------------------------------------------------------------------------"
609 : 3 : call wrtout(ab_out,message)
610 : :
611 [ - + ]: 3 : ABI_MALLOC(qphon,(3))
612 [ + + ]: 12 : qphon(:)=zero
613 [ + + ]: 3 : if(dtset%prtlden/=0) then
614 [ - + ]: 2 : call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,qphon,rhor,rhonow,lrhonow=lrhor)
615 : : else
616 : 1 : call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,qphon,rhor,rhonow)
617 : : end if
618 [ - + ]: 3 : ABI_FREE(qphon)
619 : :
620 : : ! Copy gradient which has been output in rhonow to grhor (and free rhonow)
621 : 3 : nullify(grhor)
622 [ - + + - : 15 : ABI_MALLOC(grhor,(nfftf,dtset%nspden,3))
- + + - +
- + - + -
- + - + ]
623 [ + + ]: 7 : do ispden=1,dtset%nspden
624 [ + + ]: 473039 : do ifft=1,nfftf
625 [ + + ]: 1892132 : grhor(ifft,ispden,1:3) = rhonow(ifft,ispden,2:4)
626 : : end do
627 : : end do
628 [ - + ]: 3 : ABI_FREE(rhonow)
629 : :
630 [ + - ]: 3 : if(dtset%prtgden/=0) then
631 : : ! Print result for grhor
632 : 3 : write(message,'(a,a)') ch10, " Result for gradient of the electron density for each direction (1,2,3):"
633 : 3 : call wrtout(ab_out,message)
634 : 3 : write(message,'(a,a,a,a)') ch10," 1rst direction:",ch10,&
635 : 6 : & "--------------------------------------------------------------------------------"
636 : 3 : call wrtout(ab_out,message)
637 [ - + ]: 3 : call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,1),optrhor=2,ucvol=ucvol)
638 : 3 : write(message,'(a)') "--------------------------------------------------------------------------------"
639 : 3 : call wrtout(ab_out,message)
640 : 3 : write(message,'(a,a,a,a)') ch10," 2nd direction:",ch10,&
641 : 6 : & "--------------------------------------------------------------------------------"
642 : 3 : call wrtout(ab_out,message)
643 [ - + ]: 3 : call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,2),optrhor=2,ucvol=ucvol)
644 : 3 : write(message,'(a)') "--------------------------------------------------------------------------------"
645 : 3 : call wrtout(ab_out,message)
646 : 3 : write(message,'(a,a,a,a)') ch10," 3rd direction:",ch10,&
647 : 6 : & "--------------------------------------------------------------------------------"
648 : 3 : call wrtout(ab_out,message)
649 [ - + ]: 3 : call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,3),optrhor=2,ucvol=ucvol)
650 : 3 : write(message,'(a)') "--------------------------------------------------------------------------------"
651 : 3 : call wrtout(ab_out,message)
652 : : end if
653 : :
654 [ + + ]: 3 : if(dtset%prtlden/=0) then
655 : : ! Print result for lrhor
656 : 2 : write(message,'(a,a)') ch10, " Result for Laplacian of the electron density :"
657 : 2 : call wrtout(ab_out,message)
658 : 2 : write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------"
659 : 2 : call wrtout(ab_out,message)
660 [ - + ]: 2 : call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,lrhor,optrhor=3,ucvol=ucvol)
661 : 2 : write(message,'(a)') "--------------------------------------------------------------------------------"
662 : 2 : call wrtout(ab_out,message)
663 : : end if
664 : :
665 : 3 : write(message,'(a)') "--------------------------------------------------------------------------------"
666 : 3 : call wrtout(ab_out,message)
667 : : end if
668 : :
669 : : !----------------------------------------------------------------------
670 : : !Kinetic Energy Density Calculation
671 : : !----------------------------------------------------------------------
672 : :
673 : 6953 : call timab(253,2,tsec)
674 : 6953 : call timab(254,1,tsec)
675 : :
676 : : !We use routine mkrho with option=1 to compute kinetic energy density taur (and taug)
677 [ + + + + ]: 6953 : if(dtset%usekden==0 .and. dtset%prtelf/=0)then
678 : : ! tauX are reused in outscfcv for output
679 : : ! should be deallocated there
680 : 4 : nullify(taug,taur)
681 [ - + + - : 12 : ABI_MALLOC(taug,(2,nfftf))
+ - + - -
+ - + ]
682 [ - + + - : 16 : ABI_MALLOC(taur,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
683 : 4 : tim_mkrho=5
684 : : if(dtset%prtelf/=0) then
685 : 4 : write(message,'(a,a)') ch10, " Compute ELF"
686 : 4 : call wrtout(ab_out,message)
687 : 4 : write(message,'(a)') "--------------------------------------------------------------------------------"
688 : 4 : call wrtout(ab_out,message)
689 : : end if
690 : 4 : write(message,'(a,a)') ch10, " Compute kinetic energy density"
691 : 4 : call wrtout(ab_out,message)
692 : 4 : paw_dmft%use_sc_dmft=0 ! dmft not used here
693 : 4 : paw_dmft%use_dmft=0 ! dmft not used here
694 [ + - ]: 4 : if (psps%usepaw==0) then
695 : : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,&
696 : : & npwarr,occ,paw_dmft,phnons,taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,&
697 [ - + - + ]: 4 : & option=1)
698 : : else
699 [ # # # # : 0 : ABI_MALLOC(tauwfg,(2,dtset%nfft))
# # # # #
# # # ]
700 [ # # # # : 0 : ABI_MALLOC(tauwfr,(dtset%nfft,dtset%nspden))
# # # # #
# # # # #
# # ]
701 : : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,&
702 : : & npwarr,occ,paw_dmft,phnons,tauwfg,tauwfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,&
703 : 0 : & option=1)
704 [ # # # # ]: 0 : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,tauwfg,taug,tauwfr,taur)
705 [ # # ]: 0 : ABI_FREE(tauwfg)
706 [ # # ]: 0 : ABI_FREE(tauwfr)
707 : : end if
708 [ - + ]: 4 : ABI_FREE(taug)
709 : : end if
710 : : !Print result
711 [ + + ]: 6953 : if(dtset%prtkden/=0) then
712 : 45 : write(message,'(a,a)') ch10, "Result for kinetic energy density :"
713 : 45 : call wrtout(ab_out,message)
714 : 45 : write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------"
715 : 45 : call wrtout(ab_out,message)
716 [ - + ]: 45 : call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,taur,optrhor=1,ucvol=ucvol)
717 : 45 : write(message,'(a)') "--------------------------------------------------------------------------------"
718 : 45 : call wrtout(ab_out,message)
719 : : end if
720 : :
721 : : !----------------------------------------------------------------------
722 : : !Electron Localization Function (ELF) Calculation
723 : : !----------------------------------------------------------------------
724 : :
725 : 6953 : call timab(254,2,tsec)
726 : 6953 : call timab(255,1,tsec)
727 : :
728 : : !We use routine xcden to compute gradient of electron density (grhor),
729 : : !NOTE: If GGA is used, gradient of electron density is already computed
730 : : !and it is stored in exchange correlation kernel kxc(:,5:7) (nspden=1) or kxc(:,14:19) (nspden=2).
731 : : !In order to save memory and do not have the same quantity twice
732 : : !in memory we should use kxc.
733 : : !But unfortunately only spin up ans spin down component are stored in kxc.
734 : : !So far we thus use grhor which contains all component (nspden=4)
735 : : !just like rhonow in xcden instead of kxc.
736 : :
737 [ + + ]: 6953 : if((dtset%prtelf/=0))then
738 [ + - ]: 4 : if(dtset%nspden<=2) then
739 : :
740 : 4 : ngrad=2
741 : 4 : cplex=1
742 [ - + ]: 4 : if((cplex*dtset%nfft)/=nfftf)then
743 : 0 : write(message, '(a,a,a,a)' ) ch10,&
744 : 0 : & ' afterscfloop: ERROR -', ch10, &
745 : 0 : & ' The density is complex, ELF analysis cannot be performed.'
746 : 0 : call wrtout(std_out,message)
747 : : ! ABI_ERROR(message)
748 : : end if
749 : :
750 [ + - + - ]: 4 : if((dtset%prtgden==0) .and. (dtset%prtlden==0)) then
751 : : ! Compute gradient of the electron density
752 : 4 : ishift=0
753 [ - + + - : 20 : ABI_MALLOC(rhonow,(nfftf,dtset%nspden,ngrad*ngrad))
- + + - +
- + - + -
- + - + ]
754 : 4 : write(message,'(a,a)') ch10, " Compute gradient of the electron density"
755 : 4 : call wrtout(ab_out,message)
756 [ - + ]: 4 : ABI_MALLOC(qphon,(3))
757 [ + + ]: 16 : qphon(:)=zero
758 : 4 : call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,qphon,rhor,rhonow)
759 [ - + ]: 4 : ABI_FREE(qphon)
760 : : ! Copy gradient which has been output in rhonow to grhor (and free rhonow)
761 [ - + + - : 20 : ABI_MALLOC(grhor,(nfftf,dtset%nspden,3))
- + + - +
- + - + -
- + - + ]
762 [ + + ]: 10 : do ispden=1,dtset%nspden
763 [ + + ]: 135650 : do ifft=1,nfftf
764 [ + + ]: 542566 : grhor(ifft,ispden,1:3) = rhonow(ifft,ispden,2:4)
765 : : end do
766 : : end do
767 [ - + ]: 4 : ABI_FREE(rhonow)
768 : : end if
769 : : ! Compute square norm of gradient of the electron density (|grhor|**2)
770 [ + + ]: 4 : if(dtset%nspden==1)then
771 [ - + + - : 8 : ABI_MALLOC(sqnormgrhor,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
772 [ + + ]: 27642 : do ifft=1,nfftf
773 : 27642 : sqnormgrhor(ifft,1) = zero
774 : : end do
775 [ + - ]: 2 : elseif(dtset%nspden==2)then
776 [ - + + - : 8 : ABI_MALLOC(sqnormgrhor,(nfftf,dtset%nspden+1))
- + + - +
- + - - +
- + ]
777 : : ! because we not only want (total and up quantities, but also down)
778 : : ! Indeed after having token the square norm we can not recover the
779 : : ! down quantity by substracting total and up quantities (as we do for usual densities)
780 [ + + ]: 8 : do ispden=1,dtset%nspden+1
781 [ + + ]: 162008 : do ifft=1,nfftf
782 : 162006 : sqnormgrhor(ifft,ispden) = zero
783 : : end do
784 : : end do
785 : : end if
786 : :
787 [ + + ]: 16 : do igrad=1,3
788 [ + + ]: 30 : do ispden=1,dtset%nspden !total (and up)
789 [ + + ]: 406950 : do ifft=1,nfftf
790 : 406938 : sqnormgrhor(ifft,ispden) = sqnormgrhor(ifft,ispden) + grhor(ifft,ispden,igrad)**2
791 : : end do
792 : : end do
793 [ + + ]: 16 : if(dtset%nspden==2)then
794 [ + + ]: 162006 : do ifft=1,nfftf !down
795 : 162006 : sqnormgrhor(ifft,3) = sqnormgrhor(ifft,3) + (grhor(ifft,1,igrad)-grhor(ifft,2,igrad))**2
796 : : end do
797 : : end if
798 : : end do
799 : :
800 : : ! Compute electron localization function (ELF) (here it is elfr)
801 : :
802 : 4 : nullify(elfr)
803 [ + + ]: 4 : if(dtset%nspden==1)then
804 [ - + + - : 8 : ABI_MALLOC(elfr,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + ]
805 [ + - ]: 2 : elseif(dtset%nspden==2)then
806 [ - + + - : 8 : ABI_MALLOC(elfr,(nfftf,dtset%nspden+1))
- + + - +
- + - - +
- + ]
807 : : ! 1rst is total elf, 2nd is spin-up elf, and 3rd is spin-down elf. (elf_tot /= elf_up + elf_down)
808 : : end if
809 : 4 : c_fermi = 3.0d0/10.0d0*((3.0d0*pi**2)**(2.0d0/3.0d0))
810 : :
811 : : ! First compute total elf
812 : 4 : ispden=1
813 [ + + ]: 81644 : do ifft=1,nfftf
814 : 81640 : dtaurzero = c_fermi*rhor(ifft,ispden)**(5.0d0/3.0d0)
815 : 81640 : dtaur = taur(ifft,ispden)
816 : 81640 : dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden)/rhor(ifft,ispden))
817 : : ! Ensure that dtaur is always positive or zero, as it should be.
818 : : if(dtaur<0.0d0)dtaur=0.0d0
819 : : ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur
820 [ - + ]: 81644 : if(dtaurzero<(1.0d-20*dtaur)) then
821 : 0 : elfr(ifft,ispden) = 0.0d0
822 : : else
823 : 81640 : elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2)
824 : : ! For atomic shell studies we could also add the condition that when (dtaur/dtaurzero)
825 : : ! is very close to zero we actually set it exactly to zero (or --> elfr = 1.0d0
826 : : ! which is the upper limit of elfr.)
827 : : end if
828 : : end do
829 : :
830 : : ! If spin-dependent densities are avalaible, compute spin-dependent elf
831 : : ! and offer the possibility to compute total elf in an alternative approach
832 : : ! (see doc/theory/ELF)
833 [ + + ]: 4 : if(dtset%nspden==2)then
834 : :
835 : : ! alternative approach to the total elf
836 [ + + ]: 2 : if(dtset%prtelf==2)then
837 : 27001 : ispden=1
838 [ + + ]: 27001 : do ifft=1,nfftf
839 : : dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*( rhor(ifft,ispden+1)**(5.0d0/3.0d0) + &
840 : 27000 : & (rhor(ifft,ispden) - rhor(ifft,ispden+1))**(5.0d0/3.0d0) )
841 : 27000 : dtaur = taur(ifft,ispden)
842 : : dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+1)/rhor(ifft,ispden+1)) &
843 : 27000 : & - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+2)/(rhor(ifft,ispden)-rhor(ifft,ispden+1)))
844 : : if(dtaur<0.0d0)dtaur=0.0d0
845 : : ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur
846 [ - + ]: 27001 : if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then
847 : 0 : elfr(ifft,ispden) = 0.0d0
848 : : else
849 : 27000 : elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2)
850 : : end if
851 : : end do
852 : : end if
853 : :
854 : : ! elf_up
855 : 54002 : ispden=2
856 [ + + ]: 54002 : do ifft=1,nfftf
857 : 54000 : dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*rhor(ifft,ispden)**(5.0d0/3.0d0)
858 : 54000 : dtaur = taur(ifft,ispden)
859 : 54000 : dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden)/rhor(ifft,ispden))
860 : : if(dtaur<0.0d0)dtaur=0.0d0
861 : : ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur
862 [ - + ]: 54002 : if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then
863 : 0 : elfr(ifft,ispden) = 0.0d0
864 : : else
865 : 54000 : elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2)
866 : : end if
867 : : end do
868 : :
869 : : ! elf_down
870 : 54002 : ispden=1
871 [ + + ]: 54002 : do ifft=1,nfftf
872 : 54000 : dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*(rhor(ifft,ispden)-rhor(ifft,ispden+1))**(5.0d0/3.0d0)
873 : 54000 : dtaur = taur(ifft,ispden)-taur(ifft,ispden+1)
874 : 54000 : dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+2)/(rhor(ifft,ispden)-rhor(ifft,ispden+1)))
875 : : if(dtaur<0.0d0)dtaur=0.0d0
876 : : ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur
877 [ - + ]: 54002 : if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then
878 : 0 : elfr(ifft,ispden+2) = 0.0d0
879 : : else
880 : 54000 : elfr(ifft,ispden+2) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2)
881 : : end if
882 : : end do
883 : :
884 : : end if !endif dtset%nspden==2
885 : :
886 : : ! Print result for elfr
887 [ - + ]: 4 : call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,elfr,optrhor=4,ucvol=ucvol)
888 : :
889 [ - + ]: 4 : ABI_FREE(grhor)
890 [ - + ]: 4 : ABI_FREE(sqnormgrhor)
891 : :
892 : : else
893 : 0 : message ='ELF is not yet implemented for non collinear spin cases.'
894 : 0 : ABI_WARNING(message)
895 : :
896 [ # # # # : 0 : ABI_MALLOC(elfr,(nfftf,dtset%nspden))
# # # # #
# # # # #
# # ]
897 [ # # ]: 0 : do ispden=1,dtset%nspden
898 [ # # ]: 0 : do ifft=1,nfftf
899 : 0 : elfr(ifft,ispden) = -2.0d0
900 : : end do
901 : : end do
902 : : ! even if elf is not computed we want to finish the abinit run.
903 : : ! To ensure that users won't use the _ELF file which will be produced
904 : : ! we set elf to -2.0 (a meaningless value)
905 : :
906 : : end if ! endif dtset%nspden<=2
907 : :
908 : 4 : write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------"
909 : 4 : call wrtout(ab_out,message)
910 : 4 : write(message,'(a)') " End of ELF section"
911 : 4 : call wrtout(ab_out,message)
912 : :
913 [ + - ]: 4 : if (dtset%usekden==0) then
914 [ - + ]: 4 : ABI_FREE(taur)
915 : : end if
916 : :
917 : : end if !endif prtelf/=0
918 : :
919 : : !######################################################################
920 : : !Compute forces (if they were not computed during the elec. iterations)
921 : : !and stresses (if requested by user)
922 : : !----------------------------------------------------------------------
923 : :
924 : 6953 : call timab(255,2,tsec)
925 : 6953 : call timab(256,1,tsec)
926 : :
927 : 6953 : optfor=0
928 : :
929 [ + + + + : 6953 : if (computed_forces==0.and.dtset%optforces>0.and.dtset%iscf>=0) then
+ + ]
930 [ + + + + ]: 3282 : if (dtset%nstep>0.or.dtfil%ireadwf==1) optfor=1
931 : : end if
932 : :
933 [ + + + + ]: 6953 : if (optfor>0.or.stress_needed>0) then
934 : :
935 : : ! PAW: eventually, compute g_l(r).Y_lm(r) gradients (if not already done)
936 [ + + ]: 6362 : if (psps%usepaw==1) then
937 [ + + + + ]: 4126 : test_nfgd =any(pawfgrtab(:)%nfgd==0)
938 [ + + + + ]: 2852 : test_rfgd =any(pawfgrtab(:)%rfgd_allocated==0)
939 [ + + + + ]: 2851 : test_gylmgr=any(pawfgrtab(:)%gylmgr_allocated==0)
940 : : if (test_nfgd.or.&
941 : : & (test_gylmgr.and.dtset%pawstgylm==1).or.&
942 [ + + + + : 1438 : & (test_rfgd.and.stress_needed==1.and.dtset%pawstgylm==1).or.&
+ + - + -
- - - -
- ]
943 : : (test_rfgd.and.dtset%pawstgylm==0)) then
944 : 668 : optcut=0;optgr0=0;optgr1=dtset%pawstgylm;optgr2=0
945 [ + + ]: 668 : optrad=1-dtset%pawstgylm;if (stress_needed==1) optrad=1
946 [ + - ]: 668 : if (dtset%usewvl==0) then
947 : : call nhatgrid(atindx1,gmet,my_natom,dtset%natom,nattyp,ngfftf,dtset%ntypat,&
948 : : & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,&
949 : : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
950 : 668 : & comm_fft=spaceComm_fft,distribfft=mpi_enreg%distribfft)
951 : : else
952 : : shft=0
953 : : #if defined HAVE_BIGDFT
954 : : shft=wvl%descr%Glr%d%n1i*wvl%descr%Glr%d%n2i*wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,4)
955 : : call wvl_nhatgrid(atindx1,wvl%descr%atoms%astruct%geocode,&
956 : : & wvl%descr%h,wvl%den%denspot%dpbox%i3s,dtset%natom,dtset%natom,&
957 : : & nattyp,psps%ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,&
958 : : & wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,&
959 : : & wvl%den%denspot%dpbox%n3pi,optcut,optgr0,optgr1,optgr2,optrad,&
960 : : & pawfgrtab,pawtab,psps%gth_params%psppar,rprimd,shft,xred)
961 : : #endif
962 : : end if
963 : : end if
964 : : end if
965 : :
966 : : call forstr(atindx1,cg,cprj,diffor,dtefield,dtset,&
967 : : & eigen,electronpositron,energies,favg,fcart,fock,&
968 : : & forold,gred,grchempottn,grcondft,gresid,grewtn,&
969 : : & grhf,grvdw,grxc,gsqcut,extfpmd,indsym,&
970 : : & kg,kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,&
971 : : & n3xccc,nattyp,nfftf,ngfftf,ngrvdw,nhat,nkxc,&
972 : : & npwarr,dtset%ntypat,nvresid,occ,optfor,optres,&
973 : : & paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1df,&
974 : : & psps,rhog,rhor,rprimd,stress_needed,strscondft,strsxc,strten,symrec,synlgr,&
975 [ - + - - : 6362 : & ucvol,usecprj,usevxctau,vhartr,vpsp,vxc,vxctau,wvl,xccc3d,xcctau3d,xred,ylm,ylmgr,qvpotzero,xg_nonlop)
- - - - -
+ ]
976 : : end if
977 : :
978 : : ! Init values with MAGIC_UNDEF if not computed.
979 [ + + ]: 6953 : if (optfor==1) computed_forces=1
980 [ + + ]: 6953 : if (optfor==1) diffor = MAGIC_UNDEF
981 [ + + + + ]: 10697 : if (stress_needed==0) strten = MAGIC_UNDEF
982 [ + + + + : 13597 : if (computed_forces==0) fcart = MAGIC_UNDEF
+ + ]
983 [ + + ]: 6953 : if (dtset%prtstm/=0) strten(:)=zero
984 : :
985 : 6953 : call timab(256,2,tsec)
986 : 6953 : call timab(257,1,tsec)
987 : :
988 : : !If SCF convergence was not reached (for dtset%nstep>0),
989 : : !print a warning to the output file (non-dummy arguments: dtset%nstep,
990 : : !residm, diffor - infos from tollist have been saved inside )
991 : 6953 : choice=3
992 : : call scprqt(choice,cpus,deltae,diffor,maxmag,difmag,dtset,&
993 : : & eigen,etotal,favg,fcart,energies%e_fermie,energies%e_fermih,&
994 : : & dtfil%fnameabo_app_eig,dtfil%filnam_ds(1),&
995 : : & 1,dtset%iscf,istep,istep_fock_outer,istep_mix,dtset%kptns,maxfor,&
996 : : & moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,&
997 : : & dtset%nstep,occ,optres,prtfor,prtxml,quit,&
998 : : & res2,resid,residm,response,tollist,psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode,&
999 : 6953 : & electronpositron=electronpositron, fock=fock)
1000 : :
1001 : : !output POSCAR and FORCES files, VASP style, for PHON code and friends.
1002 [ + + ]: 6953 : if (dtset%prtposcar == 1) then
1003 : 1 : call prtposcar(fcart, dtfil%filnam_ds(4), dtset%natom, dtset%ntypat, rprimd, dtset%typat, ucvol, xred, dtset%znucl)
1004 : : end if ! prtposcar
1005 : :
1006 : : if(allocated(qphon)) then
1007 : : ABI_FREE(qphon)
1008 : : end if
1009 : :
1010 : : !get current operator on wavefunctions
1011 [ + + ]: 6953 : if (dtset%prtspcur == 1) then
1012 : 1 : call spin_current(cg,dtfil,dtset,gprimd,hdr,kg,mcg,mpi_enreg,psps)
1013 : : end if
1014 : :
1015 : : !Electron-positron stuff: if last calculation was a positron minimization,
1016 : : !exchange electron and positron data in order to
1017 : : !get electronic quantities in global variables
1018 [ + + ]: 6953 : if (dtset%positron/=0) then
1019 : 17 : electronpositron%scf_converged=.false.
1020 [ + + + + ]: 17 : if (dtset%positron<0.and.electronpositron_calctype(electronpositron)==1) then
1021 : : call exchange_electronpositron(cg,cprj,dtset,eigen,electronpositron,energies,gred,mcg,mcprj,&
1022 : 6 : & mpi_enreg,my_natom,nfftf,ngfftf,nhat,npwarr,occ,paw_an,pawrhoij,rhog,rhor,strten,usecprj,vhartr)
1023 : : end if
1024 : : end if
1025 : :
1026 : : !If PAW+U and density mixing, has to update nocc_mmp
1027 [ + + + + : 6953 : if (psps%usepaw==1.and.dtset%usepawu/=0.and.(dtset%iscf>0.or.dtset%iscf==-3)) then
+ + ]
1028 : : call setnoccmmp(1,0,dmatdum,0,0,indsym,my_natom,dtset%natom,dtset%natpawu,&
1029 : : & dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,dtset%pawprtvol,&
1030 : : & pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,&
1031 : 155 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1032 : : end if
1033 : :
1034 : : !Update the content of the header (evolving variables)
1035 : 6953 : bantot=hdr%bantot
1036 [ + + ]: 6953 : if (dtset%positron==0) then
1037 : : call hdr%update(bantot,etotal,energies%e_fermie,energies%e_fermih,residm,rprimd,occ,&
1038 : : pawrhoij,xred,dtset%amu_orig(:,1),&
1039 : 6936 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1040 : : else
1041 : : call hdr%update(bantot,electronpositron%e0,energies%e_fermie,energies%e_fermih,residm,rprimd,occ,&
1042 : : pawrhoij,xred,dtset%amu_orig(:,1),&
1043 : 17 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1044 : : end if
1045 : :
1046 : : #ifdef HAVE_LOTF
1047 : : if(dtset%ionmov==23 .and. mpi_enreg%nproc_band>1) then
1048 : : bufsz=2+2*dtset%natom;if (moved_atm_inside==1) bufsz=bufsz+dtset%natom
1049 : : ABI_MALLOC(mpibuf,(3,bufsz))
1050 : : mpibuf(:,1:dtset%natom)=gred(:,1:dtset%natom)
1051 : : mpibuf(:,dtset%natom+1:2*dtset%natom)=fcart(:,1:dtset%natom)
1052 : : if (moved_atm_inside==1) mpibuf(:,2*dtset%natom+1:3*dtset%natom)=xred(:,1:dtset%natom)
1053 : : mpibuf(1:3,bufsz-1:bufsz)=reshape(strten(1:6),(/3,2/))
1054 : : call xmpi_sum(mpibuf,mpi_enreg%comm_band,ierr)
1055 : : gred(:,1:dtset%natom)=mpibuf(:,1:dtset%natom)/mpi_enreg%nproc_band
1056 : : fcart(:,1:dtset%natom)=mpibuf(:,dtset%natom+1:2*dtset%natom)/mpi_enreg%nproc_band
1057 : : if (moved_atm_inside==1) xred(:,1:dtset%natom)=mpibuf(:,2*dtset%natom+1:3*dtset%natom)/mpi_enreg%nproc_band
1058 : : strten(1:6)=reshape(mpibuf(1:3,bufsz-1:bufsz),(/6/))/mpi_enreg%nproc_band
1059 : : ABI_FREE(mpibuf)
1060 : : end if
1061 : : #endif
1062 : :
1063 : : !In case of FFT parallelisation, has to synchronize positions and forces
1064 : : !to avoid numerical noise
1065 [ + + ]: 6953 : if (mpi_enreg%nproc_fft>1) then
1066 [ - + ]: 444 : bufsz=2+2*dtset%natom;if (moved_atm_inside==1) bufsz=bufsz+dtset%natom
1067 [ - + + - : 1332 : ABI_MALLOC(mpibuf,(3,bufsz))
+ - + - -
+ - + -
+ ]
1068 [ + + + + ]: 4300 : mpibuf(:,1:dtset%natom)=gred(:,1:dtset%natom)
1069 [ + + + + ]: 4300 : mpibuf(:,dtset%natom+1:2*dtset%natom)=fcart(:,1:dtset%natom)
1070 [ - + - - : 444 : if (moved_atm_inside==1) mpibuf(:,2*dtset%natom+1:3*dtset%natom)=xred(:,1:dtset%natom)
- - ]
1071 [ + + + + ]: 3996 : mpibuf(1:3,bufsz-1:bufsz)=reshape(strten(1:6),(/3,2/))
1072 : 444 : call xmpi_sum(mpibuf,mpi_enreg%comm_fft,ierr)
1073 [ + + + + ]: 4300 : gred(:,1:dtset%natom)=mpibuf(:,1:dtset%natom)/mpi_enreg%nproc_fft
1074 [ + + + + ]: 4300 : fcart(:,1:dtset%natom)=mpibuf(:,dtset%natom+1:2*dtset%natom)/mpi_enreg%nproc_fft
1075 [ - + - - : 444 : if (moved_atm_inside==1) xred(:,1:dtset%natom)=mpibuf(:,2*dtset%natom+1:3*dtset%natom)/mpi_enreg%nproc_fft
- - ]
1076 [ + + ]: 3108 : strten(1:6)=reshape(mpibuf(1:3,bufsz-1:bufsz),(/6/))/mpi_enreg%nproc_fft
1077 [ - + ]: 888 : ABI_FREE(mpibuf)
1078 : : end if
1079 : :
1080 : : !results_gs%energies = energies
1081 : 6953 : call energies%copy(results_gs%energies)
1082 : 6953 : results_gs%etotal =etotal
1083 : 6953 : results_gs%deltae =deltae
1084 : 6953 : results_gs%diffor =diffor
1085 : 6953 : results_gs%residm =residm
1086 : 6953 : results_gs%res2 =res2
1087 [ + + + + ]: 88545 : results_gs%fcart(:,:) =fcart(:,:)
1088 [ + + + + ]: 88545 : results_gs%gred(:,:) =gred(:,:)
1089 [ + + + + ]: 88545 : results_gs%grchempottn(:,:)=grchempottn(:,:)
1090 [ + + + + ]: 88545 : results_gs%gresid(:,:)=gresid(:,:)
1091 [ + + + + ]: 88545 : results_gs%grewtn(:,:)=grewtn(:,:)
1092 [ + + + + ]: 88545 : results_gs%grxc(:,:) =grxc(:,:)
1093 : 6953 : results_gs%berryopt =dtset%berryopt
1094 [ + + ]: 27812 : results_gs%pel(1:3) =pel(1:3)
1095 [ + + ]: 27812 : results_gs%pion(1:3) =pion(1:3)
1096 [ + + ]: 48671 : results_gs%strten(1:6)=strten(1:6)
1097 [ + + + + ]: 88545 : results_gs%synlgr(:,:)=synlgr(:,:)
1098 : 6953 : results_gs%vxcavg =vxcavg
1099 [ + + + + : 7029 : if (ngrvdw>0) results_gs%grvdw(1:3,1:ngrvdw)=grvdw(1:3,1:ngrvdw)
+ + ]
1100 [ + + ]: 6953 : if (associated(extfpmd)) then
1101 [ - + ]: 7 : if(dtset%extfpmd_prterr==1) then
1102 : : call extfpmd_err(extfpmd,eigen,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,&
1103 [ # # ]: 0 : & dtset%wtk,trim(dtfil%filnam_ds(4))//'_EXTFPMD_ERR')
1104 : : endif
1105 : 7 : results_gs%nelect_extfpmd=extfpmd%nelect
1106 : 7 : results_gs%extfpmd_eshift=extfpmd%eshift
1107 : : end if
1108 : :
1109 [ + + + + ]: 50044 : results_gs%intgres(:,:)=zero
1110 [ + + + + ]: 88545 : results_gs%grcondft(:,:)=zero
1111 [ + + + + : 16486 : if(any(dtset%constraint_kind(:)/=0))then
+ + ]
1112 [ + + + + ]: 606 : results_gs%intgres(1:dtset%nspden,:) =intgres(1:dtset%nspden,:)
1113 [ + + + + ]: 558 : results_gs%grcondft(:,:) =grcondft(:,:)
1114 : : endif
1115 : :
1116 [ + + + + : 6953 : if (dtset%nstep == 0 .and. dtset%occopt>=3.and.dtset%occopt<=8) then
+ - ]
1117 : 2 : results_gs%etotal = results_gs%etotal - dtset%tsmear * results_gs%entropy
1118 : : end if
1119 : :
1120 : : !This call is only for testing purpose:
1121 : : !test of the nonlop routine (DFPT vs Finite Differences)
1122 [ - + ]: 6953 : if (dtset%useria==112233) then
1123 : : call nonlop_test(cg,eigen,dtset%istwfk,kg,dtset%kptns,dtset%mband,mcg,dtset%mgfft,dtset%mkmem,&
1124 : : & mpi_enreg,dtset%mpw,my_natom,dtset%natom,dtset%nband,dtset%nfft,dtset%ngfft,dtset%nkpt,&
1125 : : & dtset%nloalg,npwarr,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,paw_ij,&
1126 : 0 : & pawtab,ph1d,psps,rprimd,dtset%typat,xred)
1127 : : end if
1128 : :
1129 : 6953 : call crystal%free()
1130 : 6953 : call ebands_k%free()
1131 : :
1132 : 6953 : call timab(257,2,tsec)
1133 : 6953 : call timab(250,2,tsec)
1134 : :
1135 : : DBG_EXIT("COLL")
1136 : :
1137 : : #if !defined HAVE_BIGDFT
1138 : : if (.false.) write(std_out,*) vtrial(1,1)
1139 : : #endif
1140 : :
1141 [ - + - + : 6953 : end subroutine afterscfloop
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - + -
+ - + - +
- + - - -
- - - - -
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
1142 : : !!***
1143 : :
1144 : : end module m_afterscfloop
1145 : : !!***
|
