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 6918 : subroutine afterscfloop(atindx,atindx1,cg,computed_forces,cprj,cpus,&
281 6918 : & deltae,diffor,difmag,dtefield,dtfil,dtset,eigen,electronpositron,elfr,&
282 13836 : & energies,etotal,extfpmd,favg,fcart,fock,forold,grchempottn,grcondft,&
283 13836 : & gred,gresid,grewtn,grhf,grhor,grvdw,&
284 20754 : & grxc,gsqcut,hdr,indsym,intgres,irrzon,istep,istep_fock_outer,istep_mix,&
285 13836 : & kg,kxc,lrhor,maxfor,maxmag,mcg,mcprj,mgfftf,&
286 13836 : & moved_atm_inside,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfft,ngfftf,ngrvdw,nhat,&
287 13836 : & nkxc,npwarr,nvresid,occ,optres,paw_an,paw_ij,pawang,pawfgr,&
288 6918 : & pawfgrtab,pawrad,pawrhoij,pawtab,pel,pel_cg,ph1d,ph1df,phnons,pion,prtfor,prtxml,&
289 13836 : & psps,pwind,pwind_alloc,pwnsfac,res2,resid,residm,results_gs,&
290 13836 : & rhog,rhor,rprimd,stress_needed,strscondft,strsxc,strten,symrec,synlgr,taug,&
291 13836 : & taur,tollist,usecprj,usevxctau,vhartr,vpsp,vtrial,vxc,vxctau,vxcavg,wvl,&
292 6918 : & 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 359736 : type(crystal_t) :: crystal
374 6918 : type(ebands_t) :: ebands_k
375 6918 : 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 6918 : real(dp),allocatable :: cg1_3(:,:,:),doccde(:),mpibuf(:,:),qphon(:),rhonow(:,:,:),sqnormgrhor(:,:)
385 6918 : 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 6918 : call timab(250,1,tsec)
395 6918 : call timab(251,1,tsec)
396 :
397 : !Compute different geometric tensor, as well as ucvol, from rprimd
398 6918 : 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 13806 : & symrel=dtset%symrel,tnons=dtset%tnons,symafm=dtset%symafm)
405 :
406 20754 : ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
407 497468 : 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 6918 : & hdr%kptrlatt, hdr%nshiftk, hdr%shiftk)
413 6918 : ABI_FREE(doccde)
414 :
415 : !MPI FFT communicator
416 6918 : spaceComm_fft=mpi_enreg%comm_fft
417 :
418 : !Recompute structure factor phases if atomic positions have changed
419 6918 : 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 6918 : 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 6918 : call timab(251,2,tsec)
540 6918 : call timab(252,1,tsec)
541 :
542 : !----------------------------------------------------------------------
543 : !Polarization Calculation, but not orbital magnetism
544 : !----------------------------------------------------------------------
545 :
546 6918 : 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 6918 : call timab(252,2,tsec)
554 6918 : call timab(253,1,tsec)
555 :
556 : !----------------------------------------------------------------------
557 : !Orbital magnetization calculation using PEAD DDK wavefunctions
558 : !----------------------------------------------------------------------
559 :
560 6918 : 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 6918 : 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 6918 : call timab(253,2,tsec)
674 6918 : call timab(254,1,tsec)
675 :
676 : !We use routine mkrho with option=1 to compute kinetic energy density taur (and taug)
677 6918 : 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 6918 : 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 6918 : call timab(254,2,tsec)
726 6918 : 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 6918 : 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 6918 : call timab(255,2,tsec)
925 6918 : call timab(256,1,tsec)
926 :
927 6918 : optfor=0
928 :
929 6918 : if (computed_forces==0.and.dtset%optforces>0.and.dtset%iscf>=0) then
930 3270 : if (dtset%nstep>0.or.dtfil%ireadwf==1) optfor=1
931 : end if
932 :
933 6918 : 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 6333 : if (psps%usepaw==1) then
937 4061 : test_nfgd =any(pawfgrtab(:)%nfgd==0)
938 2797 : test_rfgd =any(pawfgrtab(:)%rfgd_allocated==0)
939 2796 : test_gylmgr=any(pawfgrtab(:)%gylmgr_allocated==0)
940 : if (test_nfgd.or.&
941 : & (test_gylmgr.and.dtset%pawstgylm==1).or.&
942 1417 : & (test_rfgd.and.stress_needed==1.and.dtset%pawstgylm==1).or.&
943 : (test_rfgd.and.dtset%pawstgylm==0)) then
944 664 : optcut=0;optgr0=0;optgr1=dtset%pawstgylm;optgr2=0
945 664 : optrad=1-dtset%pawstgylm;if (stress_needed==1) optrad=1
946 664 : 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 664 : & 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 6333 : & 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 6918 : if (optfor==1) computed_forces=1
980 6918 : if (optfor==1) diffor = MAGIC_UNDEF
981 10614 : if (stress_needed==0) strten = MAGIC_UNDEF
982 13530 : if (computed_forces==0) fcart = MAGIC_UNDEF
983 6918 : if (dtset%prtstm/=0) strten(:)=zero
984 :
985 6918 : call timab(256,2,tsec)
986 6918 : 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 6918 : 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 6918 : & electronpositron=electronpositron, fock=fock)
1000 :
1001 : !output POSCAR and FORCES files, VASP style, for PHON code and friends.
1002 6918 : 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 6918 : 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 6918 : 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 6918 : 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 154 : & 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 6918 : bantot=hdr%bantot
1036 6918 : 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 6901 : 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 6918 : if (mpi_enreg%nproc_fft>1) then
1066 440 : bufsz=2+2*dtset%natom;if (moved_atm_inside==1) bufsz=bufsz+dtset%natom
1067 1320 : ABI_MALLOC(mpibuf,(3,bufsz))
1068 4280 : mpibuf(:,1:dtset%natom)=gred(:,1:dtset%natom)
1069 4280 : mpibuf(:,dtset%natom+1:2*dtset%natom)=fcart(:,1:dtset%natom)
1070 440 : if (moved_atm_inside==1) mpibuf(:,2*dtset%natom+1:3*dtset%natom)=xred(:,1:dtset%natom)
1071 3960 : mpibuf(1:3,bufsz-1:bufsz)=reshape(strten(1:6),(/3,2/))
1072 440 : call xmpi_sum(mpibuf,mpi_enreg%comm_fft,ierr)
1073 4280 : gred(:,1:dtset%natom)=mpibuf(:,1:dtset%natom)/mpi_enreg%nproc_fft
1074 4280 : fcart(:,1:dtset%natom)=mpibuf(:,dtset%natom+1:2*dtset%natom)/mpi_enreg%nproc_fft
1075 440 : if (moved_atm_inside==1) xred(:,1:dtset%natom)=mpibuf(:,2*dtset%natom+1:3*dtset%natom)/mpi_enreg%nproc_fft
1076 3080 : strten(1:6)=reshape(mpibuf(1:3,bufsz-1:bufsz),(/6/))/mpi_enreg%nproc_fft
1077 880 : ABI_FREE(mpibuf)
1078 : end if
1079 :
1080 : !results_gs%energies = energies
1081 6918 : call energies%copy(results_gs%energies)
1082 6918 : results_gs%etotal =etotal
1083 6918 : results_gs%deltae =deltae
1084 6918 : results_gs%diffor =diffor
1085 6918 : results_gs%residm =residm
1086 6918 : results_gs%res2 =res2
1087 88262 : results_gs%fcart(:,:) =fcart(:,:)
1088 88262 : results_gs%gred(:,:) =gred(:,:)
1089 88262 : results_gs%grchempottn(:,:)=grchempottn(:,:)
1090 88262 : results_gs%gresid(:,:)=gresid(:,:)
1091 88262 : results_gs%grewtn(:,:)=grewtn(:,:)
1092 88262 : results_gs%grxc(:,:) =grxc(:,:)
1093 6918 : results_gs%berryopt =dtset%berryopt
1094 27672 : results_gs%pel(1:3) =pel(1:3)
1095 27672 : results_gs%pion(1:3) =pion(1:3)
1096 48426 : results_gs%strten(1:6)=strten(1:6)
1097 88262 : results_gs%synlgr(:,:)=synlgr(:,:)
1098 6918 : results_gs%vxcavg =vxcavg
1099 6994 : if (ngrvdw>0) results_gs%grvdw(1:3,1:ngrvdw)=grvdw(1:3,1:ngrvdw)
1100 6918 : 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 49863 : results_gs%intgres(:,:)=zero
1110 88262 : results_gs%grcondft(:,:)=zero
1111 16395 : 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 6918 : if (dtset%nstep == 0 .and. dtset%occopt>=3.and.dtset%occopt<=8) then
1117 0 : 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 6918 : 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 6918 : call crystal%free()
1130 6918 : call ebands_k%free()
1131 :
1132 6918 : call timab(257,2,tsec)
1133 6918 : 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 6918 : end subroutine afterscfloop
1142 : !!***
1143 :
1144 : end module m_afterscfloop
1145 : !!***
|
