TABLE OF CONTENTS
ABINIT/dfpt_scfcv [ Functions ]
NAME
dfpt_scfcv
FUNCTION
Conducts set of passes or overall iterations of preconditioned conjugate gradient algorithm to converge wavefunctions to optimum and optionally to compute mixed derivatives of energy.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (XG, DRH, MB, XW, MT, SPr) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband*mkmem*nsppol)=pw coefficients of GS wavefunctions at k. cgq(2,mpw1*nspinor*mband*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q. cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k projected with non-local projectors: cprj=<p_i|Cnk> cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q projected with non-local projectors: cprjq=<p_i|Cnk+q> cpus= cpu time limit in seconds doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset eew=2nd derivative of Ewald energy (hartree) efrhar=Contribution from frozen-wavefunction, hartree energy, to the second-derivative of total energy. efrkin=Contribution from frozen-wavefunction, kinetic energy, to the second-derivative of total energy. efrloc=Contribution from frozen-wavefunction, local potential, to the second-derivative of total energy. efrnl=Contribution from frozen-wavefunction, non-local potential, to the second-derivative of total energy. efrx1=Contribution from frozen-wavefunction, xc core correction(1), to the second-derivative of total energy. efrx2=Contribution from frozen-wavefunction, xc core correction(2), to the second-derivative of total energy. eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree) eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree) eii=2nd derivative of pseudopotential core energy (hartree) evdw=DFT-D semi-empirical part of 2nd-order total energy fermie=fermi energy (Hartree) hdr <type(hdr_type)>=the header of wf, den and pot files idir=direction of the current perturbation indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points indsy1(4,nsym1,natom)=indirect indexing array for atom labels ipert=type of the perturbation irrzon1(nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data for RF symmetries istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced planewave coordinates at k kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points. kxc(nfftf,nkxc)=exchange and correlation kernel (see rhotoxc.f) mgfftf=maximum size of 1D FFTs for the "fine" grid (see NOTES in respfn.F90) mkmem =number of k points treated by this node (GS data) mkqmem =number of k+q points which can fit in memory (GS data); 0 if use disk mk1mem =number of k points which can fit in memory (RF data); 0 if use disk mpert=maximum number of ipert mpw=maximum dimensioned size of npw for wfs at k. mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs). nattyp(ntypat)= # atoms of each type. nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point, for each polarization ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>) nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90) ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid (see NOTES in respfn.F90) nkpt=number of k points in the full BZ nkpt_rbz=number of k points in the reduced BZ for this perturbation nkxc=second dimension of the kxc array. mpi_enreg=informations about MPI parallelization my_natom=number of atoms treated by current processor npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point nspden=number of spin-density components nspinor=number of spinorial components of the wavefunctions nsym1=number of symmetry elements in space group consistent with perturbation n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used otherwise, nfftf occkq(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2) at each k+q point of the reduced Brillouin zone. occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2) at each k point of the reduced Brillouin zone. paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS pawang <type(pawang_type)>=paw angular mesh and related data pawang1 <type(pawang_type)>=pawang datastr. containing only symmetries preserving the perturbation pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data pertcase=fuill index of the perturbation phnons1(2,nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic transl. phases, for RF symmetries ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information ph1df(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information for the "fine" grid prtbbb=if 1, band-by-band decomposition (also dim of d2bbb) psps <type(pseudopotential_type)>=variables related to pseudopotentials qphon(3)=reduced coordinates for the phonon wavelength rhog(2,nfftf)=array for Fourier transform of GS electron density rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3. rprimd(3,3)=dimensional primitive translations in real space (bohr) symaf1(nsym1)=anti(ferromagnetic) part of symmetry operations symrc1(3,3,nsym1)=symmetry operations in reciprocal space symrl1(3,3,nsym1)=symmetry operations in real space in terms of primitive translations usecprj= 1 if cprj, cprjq arrays are stored in memory useylmgr = 1 if ylmgr array is allocated useylmgr1= 1 if ylmgr1 array is allocated ddk<wfk_t>=ddk file vpsp1(cplex*nfftf)=first-order derivative of the ionic potential vtrial(nfftf,nspden)=GS potential (Hartree). vxc(nfftf,nspden)=Exchange-Correlation GS potential (Hartree) wtk_rbz(nkpt_rbz)=weight for each k point in the reduced Brillouin zone xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm*useylmgr)= gradients of real spherical harmonics at k ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm*useylmgr1)= gradients of real spherical harmonics at k+q
OUTPUT
blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
cg1_active(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF
wavefunctions at k,q. They are orthogonalized to the active.
d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
second order derivatives
d2lo(2,mpert,3,mpert)=local contributions to the 2DTEs
d2nl(2,mpert,3,mpert)=non-local contributions to the 2DTEs
d2ovl(2,mpert,3,mpert*usepaw)=1st-order change of WF overlap contributions to the 2DTEs
eberry=energy associated with Berry phase
edocc=correction to 2nd-order total energy coming from changes of occupation
eeig0=0th-order eigenenergies part of 2nd-order total energy
ehart01=inhomogeneous 1st-order Hartree part of 2nd-order total energy
for strain perturbation only (zero otherwise, and not used)
ehart1=1st-order Hartree part of 2nd-order total energy
eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues (hartree)
ek0=0th-order kinetic energy part of 2nd-order total energy.
ek1=1st-order kinetic energy part of 2nd-order total energy.
eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy
elpsp1=1st-order local pseudopot. part of 2nd-order total energy.
enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy.
enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy.
eovl1=1st-order change of wave-functions overlap, part of 2nd-order energy
PAW only - Eq(79) and Eq(80) of PRB 78, 035105 (2008)
epaw1=1st-order PAW on-site part of 2nd-order total energy.
etotal=total energy (sum of 7 contributions) (hartree)
exc1=1st-order exchange-correlation part of 2nd-order total energy.
gh1c_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK>
gh0c1_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(0)}|\Psi^{(1)}>
The wavefunction is orthogonal to the active space (for metals). It is not
coherent with cg1.
resid(mband*nkpt_rbz*nsppol)=residuals for each band over all k points
of the reduced Brillouin zone, and spins
residm=maximum value from resid array (except for nbdbuf highest bands)
conv_retcode=return code, 0 if convergence was achieved.
SIDE EFFECTS
cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=updated wavefunctions (ortho. to occ. states);
initialized= if 0 the initialization of the RF run is not yet finished
mpi_enreg=informations about MPI parallelization
rhog1(2,nfftf)=array for Fourier transform of RF electron density
rhor1(cplex*nfftf,nspden)=array for RF electron density in electrons/bohr**3.
=== if psps%usepaw==1
pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
PARENTS
dfpt_looppert
CHILDREN
ab7_mixing_deallocate,ab7_mixing_new,ab7_mixing_use_disk_cache,appdig calcdensph,destroy_efield,dfpt_etot,dfpt_newvtr,dfpt_nselt,dfpt_nstdy dfpt_nstpaw,dfpt_rhofermi,dfpt_rhotov,dfpt_vtorho,dfptff_bec,dfptff_die dfptff_ebp,dfptff_edie,dfptff_initberry,fftdatar_write_from_hdr,fourdp getcut,hdr_update,metric,newfermie1,paw_an_free,paw_an_init paw_an_nullify,paw_an_reset_flags,paw_ij_free,paw_ij_init paw_ij_nullify,paw_ij_reset_flags,pawcprj_alloc,pawcprj_free pawcprj_getdim,pawdenpot,pawdij,pawdijfr,pawmknhat,pawnhatfr pawrhoij_alloc,pawrhoij_free,qmatrix,rf2_getidirs,scprqt,status,symdij timab,wfk_close,wrtout,xmpi_barrier,xmpi_isum,xmpi_sum,xmpi_wait
SOURCE
188 #if defined HAVE_CONFIG_H 189 #include "config.h" 190 #endif 191 192 #include "abi_common.h" 193 194 subroutine dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,& 195 & dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset,& 196 & d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,& 197 & ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,eloc0,elpsp1,& 198 & enl0,enl1,eovl1,epaw1,etotal,evdw,exc1,fermie,gh0c1_set,gh1c_set,hdr,idir,indkpt1,& 199 & indsy1,initialized,ipert,irrzon1,istwfk_rbz,& 200 & kg,kg1,kpt_rbz,kxc,mgfftf,mkmem,mkqmem,mk1mem,& 201 & mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,nattyp,nband_rbz,ncpgr,& 202 & nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,npwarr,npwar1,nspden,& 203 & nsym1,n3xccc,occkq,occ_rbz,& 204 & paw_an,paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,pawrhoij1,pawtab,& 205 & pertcase,phnons1,ph1d,ph1df,& 206 & prtbbb,psps,qphon,resid,residm,rhog,rhog1,& 207 & rhor,rhor1,rprimd,symaf1,symrc1,symrl1,& 208 & usecprj,useylmgr,useylmgr1,ddk_f,vpsp1,vtrial,vxc,& 209 & wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,conv_retcode,& 210 & kramers_deg,& 211 & cg_mq,cg1_mq,cg1_active_mq,docckde_mq,eigen_mq,eigen1_mq,gh0c1_set_mq,gh1c_set_mq,& 212 & kg1_mq,npwar1_mq,occk_mq,resid_mq,residm_mq,rhog1_pq,rhog1_mq,rhor1_pq,rhor1_mq) 213 214 use defs_basis 215 use defs_datatypes 216 use defs_abitypes 217 use m_ab7_mixing 218 use m_efield 219 use m_errors 220 use m_profiling_abi 221 use m_wfk 222 use m_xmpi 223 use m_nctk 224 use m_hdr 225 #ifdef HAVE_NETCDF 226 use netcdf 227 #endif 228 229 use m_cgtools, only : mean_fftr 230 use m_fstrings, only : int2char4, sjoin 231 use m_geometry, only : metric 232 use m_time, only : abi_wtime, sec2str 233 use m_io_tools, only : open_file 234 use m_exit, only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in 235 use m_mpinfo, only : iwrite_fftdatar 236 use m_kg, only : getcut 237 use m_ioarr, only : ioarr, fftdatar_write_from_hdr, fort_denpot_skip 238 use m_pawang, only : pawang_type 239 use m_pawrad, only : pawrad_type 240 use m_pawtab, only : pawtab_type 241 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags 242 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags 243 use m_pawfgrtab,only : pawfgrtab_type 244 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_io, pawrhoij_get_nspden 245 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_getdim 246 use m_pawdij, only : pawdij, pawdijfr, symdij 247 use m_pawfgr, only : pawfgr_type 248 use m_rf2, only : rf2_getidirs 249 250 !This section has been created automatically by the script Abilint (TD). 251 !Do not modify the following lines by hand. 252 #undef ABI_FUNC 253 #define ABI_FUNC 'dfpt_scfcv' 254 use interfaces_14_hidewrite 255 use interfaces_18_timing 256 use interfaces_32_util 257 use interfaces_53_ffts 258 use interfaces_54_abiutil 259 use interfaces_65_paw 260 use interfaces_67_common 261 use interfaces_72_response, except_this_one => dfpt_scfcv 262 !End of the abilint section 263 264 implicit none 265 266 !Arguments ------------------------------------ 267 type(dataset_type),intent(in) :: dtset 268 type(pseudopotential_type),intent(in) :: psps 269 integer,intent(in) :: cplex,dim_eig2rf,idir,ipert,mgfftf,mk1mem,mkmem,mkqmem 270 integer,intent(in) :: mpert,mpw,mpw1,my_natom,n3xccc,ncpgr,nfftf 271 integer,intent(in) :: mpw1_mq !-q duplicate 272 integer,intent(in) :: nkpt,nkpt_rbz,nkxc,nspden 273 integer,intent(in) :: nsym1,pertcase,prtbbb,usecprj,useylmgr,useylmgr1 274 logical,intent(in) :: kramers_deg 275 integer,intent(inout) :: initialized 276 ! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise 277 integer,intent(in) :: atindx(dtset%natom) 278 integer,intent(out) :: blkflg(3,mpert,3,mpert) 279 integer,intent(in) :: indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom) 280 integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 281 integer,intent(in) :: istwfk_rbz(nkpt_rbz) 282 integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem),nattyp(psps%ntypat) 283 integer,intent(in) :: nband_rbz(nkpt_rbz*dtset%nsppol) 284 integer,intent(in) :: npwar1(nkpt_rbz),npwarr(nkpt_rbz) 285 integer,optional,intent(in) :: npwar1_mq(nkpt_rbz) !-q duplicate 286 integer,optional,intent(in) :: kg1_mq(3,mpw1_mq*mk1mem)! 287 integer,intent(in) :: symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1) 288 integer,intent(out) :: conv_retcode 289 real(dp),intent(in) :: cpus,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,eii 290 real(dp),intent(out) :: eberry,edocc,eeig0,ehart01,ehart1,ek0,ek1,eloc0,elpsp1,enl0 291 real(dp),intent(out) :: enl1,eovl1,epaw1,etotal,evdw,exc1,residm 292 real(dp),optional,intent(out) :: residm_mq !-q duplicate 293 real(dp),intent(inout) :: fermie 294 real(dp),intent(in) :: qphon(3) 295 ! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise 296 integer,intent(in) :: ngfftf(18) 297 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*dtset%nsppol) 298 real(dp),intent(inout) :: cg1(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol) 299 real(dp),intent(out) :: cg1_active(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf) 300 real(dp),intent(out) :: gh1c_set(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf) 301 real(dp),intent(out) :: gh0c1_set(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf) 302 real(dp),intent(in) :: cgq(2,mpw1*dtset%nspinor*dtset%mband*mkqmem*dtset%nsppol) 303 real(dp),optional,intent(inout) :: cg1_mq(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol) !start -q duplicates 304 real(dp),optional,intent(out) :: cg1_active_mq(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf)! 305 real(dp),optional,intent(out) :: gh1c_set_mq(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf) ! 306 real(dp),optional,intent(out) :: gh0c1_set_mq(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf) ! 307 real(dp),optional,intent(in) :: cg_mq(2,mpw1_mq*dtset%nspinor*dtset%mband*mkqmem*dtset%nsppol) ! 308 real(dp),optional,intent(in) :: eigen_mq(dtset%mband*nkpt_rbz*dtset%nsppol) ! 309 real(dp),optional,intent(in) :: docckde_mq(dtset%mband*nkpt_rbz*dtset%nsppol) ! 310 real(dp),optional,intent(out) :: eigen1_mq(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol) ! 311 real(dp),optional,intent(in) :: occk_mq(dtset%mband*nkpt_rbz*dtset%nsppol) ! 312 real(dp),optional,intent(out) :: resid_mq(dtset%mband*nkpt_rbz*nspden) !end 313 real(dp),intent(out) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb) 314 real(dp),intent(out) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) 315 real(dp),intent(out) :: d2ovl(2,3,mpert,3,mpert*psps%usepaw) 316 real(dp),intent(in) :: dielt(3,3) 317 real(dp),intent(in) :: doccde_rbz(dtset%mband*nkpt_rbz*dtset%nsppol) 318 real(dp),intent(in) :: docckqde(dtset%mband*nkpt_rbz*dtset%nsppol) 319 real(dp),intent(in) :: eigen0(dtset%mband*nkpt_rbz*dtset%nsppol) 320 real(dp),intent(out) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol) 321 real(dp),intent(in) :: eigenq(dtset%mband*nkpt_rbz*dtset%nsppol) 322 real(dp),intent(in) :: kpt_rbz(3,nkpt_rbz),kxc(nfftf,nkxc) 323 real(dp),intent(in) :: nhat(nfftf,dtset%nspden) 324 real(dp),intent(in) :: occ_rbz(dtset%mband*nkpt_rbz*dtset%nsppol) 325 real(dp),intent(in) :: occkq(dtset%mband*nkpt_rbz*dtset%nsppol) 326 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom),ph1df(2,3*(2*mgfftf+1)*dtset%natom) 327 real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 328 real(dp),intent(out) :: resid(dtset%mband*nkpt_rbz*nspden) 329 real(dp),intent(in) :: rhog(2,nfftf),rhor(nfftf,nspden),rprimd(3,3) 330 real(dp),intent(inout) :: rhog1(2,nfftf),rhor1(cplex*nfftf,nspden),xred(3,dtset%natom) 331 real(dp),optional,intent(inout) :: rhog1_pq(2,nfftf),rhor1_pq(cplex*nfftf,nspden) !+q/-q duplicates 332 real(dp),optional,intent(inout) :: rhog1_mq(2,nfftf),rhor1_mq(cplex*nfftf,nspden) ! 333 real(dp),target,intent(in) :: vtrial(nfftf,nspden) 334 real(dp),intent(in) :: vpsp1(cplex*nfftf),vxc(nfftf,nspden) 335 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xccc3d1(cplex*n3xccc) 336 real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 337 real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm) 338 real(dp),intent(in) :: ylmgr(mpw*mkmem,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 339 real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3+6*((ipert-dtset%natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 340 real(dp),intent(in) :: zeff(3,3,dtset%natom) 341 type(pawcprj_type),intent(in) :: cprj(dtset%natom,dtset%nspinor*dtset%mband*mkmem*dtset%nsppol*usecprj) 342 type(pawcprj_type),intent(in) :: cprjq(dtset%natom,dtset%nspinor*dtset%mband*mkqmem*dtset%nsppol*usecprj) 343 type(datafiles_type),intent(in) :: dtfil 344 type(hdr_type),intent(inout) :: hdr 345 type(pawang_type),intent(in) :: pawang,pawang1 346 type(pawfgr_type),intent(in) :: pawfgr 347 type(paw_an_type),intent(in) :: paw_an(my_natom*psps%usepaw) 348 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw) 349 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 350 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw) 351 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*psps%usepaw) 352 type(pawrhoij_type),intent(inout) :: pawrhoij1(my_natom*psps%usepaw) 353 type(pawtab_type), intent(in) :: pawtab(psps%ntypat*psps%usepaw) 354 type(MPI_type),intent(inout) :: mpi_enreg 355 type(wfk_t),intent(inout) :: ddk_f(4) 356 357 !Local variables------------------------------- 358 !scalars 359 integer,parameter :: level=12,response=1 360 integer :: afford,bantot_rbz,choice,cplex_rhoij,dbl_nnsclo 361 integer :: has_dijfr,iatom,ider,idir_dum,idir_paw1,ierr,iexit,errid,denpot 362 integer :: iprcel,iscf10_mod,iscf_mod,ispden,ispmix 363 integer :: istep,itypat,izero,lmn2_size,me,mgfftdiel,mvdum 364 integer :: nfftdiel,nfftmix,nfftotf,nhat1grdim,npawmix,npwdiel,nspden_rhoij,nstep,nzlmopt 365 integer :: optene,optfr,option,optres,prtfor,quit,quit_sum,qzero 366 integer :: my_quit,quitsum_request,timelimit_exit,varid,ncerr,ncid 367 integer ABI_ASYNC :: quitsum_async 368 integer :: rdwrpaw,spaceComm,sz1,sz2,usexcnhat,Z_kappa 369 integer :: dbl_nnsclo_mq,ifft !-q duplicate for dbl_nnsclo 370 !integer :: pqmq ! pqmq = indicator for potential mixing 371 logical :: need_fermie1,paral_atom,use_nhat_gga 372 real(dp) :: wtime_step,now,prev 373 real(dp) :: born,born_bar,boxcut,deltae,diffor,diel_q,dum,ecut,ecutf,elast 374 real(dp) :: epawdc1_dum,evar,fe1fixed,fermie1,gsqcut,qphon_norm,maxfor,renorm,res2,res3,residm2 375 real(dp) :: ucvol,vxcavg,elmag1 376 real(dp) :: res2_mq,fe1fixed_mq,elast_mq 377 real(dp) :: eberry_mq,edocc_mq,eeig0_mq,ehart01_mq,ehart1_mq,ek0_mq,ek1_mq,eloc0_mq,elpsp1_mq,enl0_mq 378 real(dp) :: enl1_mq,eovl1_mq,epaw1_mq,exc1_mq,fermie1_mq,deltae_mq,elmag1_mq 379 character(len=500) :: msg 380 character(len=fnlen) :: fi1o 381 !character(len=fnlen) :: fi1o_vtk 382 integer :: prtopt 383 type(ab7_mixing_object) :: mix 384 type(efield_type) :: dtefield 385 !arrays 386 integer :: ngfftmix(18) 387 integer,allocatable :: dimcprj(:),pwindall(:,:,:) 388 integer,pointer :: my_atmtab(:) 389 real(dp) :: dielar(7) 390 real(dp) :: favg(3),gmet(3,3),gprimd(3,3),q_cart(3),qphon2(3),qred2cart(3,3) 391 real(dp) :: rmet(3,3),tollist(12),tsec(2) 392 real(dp) :: zeff_red(3),zeff_bar(3,3) 393 real(dp) :: intgden(dtset%nspden,dtset%natom),dentot(dtset%nspden) 394 !real(dp) :: zdmc_red(3),zdmc_bar(3,3),mean_rhor1(1) !dynamic magnetic charges and mean density 395 real(dp),allocatable :: dielinv(:,:,:,:,:) 396 real(dp),allocatable :: fcart(:,:),nhat1(:,:),nhat1gr(:,:,:),nhatfermi(:,:),nvresid1(:,:),nvresid2(:,:) 397 real(dp),allocatable :: qmat(:,:,:,:,:,:),resid2(:),rhog2(:,:),rhor2(:,:),rhorfermi(:,:) 398 real(dp),allocatable :: susmat(:,:,:,:,:),vhartr1(:),vxc1(:,:) 399 real(dp),allocatable :: vhartr1_tmp(:,:) 400 real(dp),allocatable,target :: vtrial1(:,:),vtrial2(:,:) 401 real(dp),allocatable :: vtrial1_pq(:,:),vtrial1_mq(:,:),rhorfermi_mq(:,:) 402 real(dp),allocatable :: nvresid1_mq(:,:),vxc1_mq(:,:),vhartr1_mq(:) 403 real(dp),pointer :: vtrial1_tmp(:,:) 404 type(pawcprj_type),allocatable :: cprj1(:,:) 405 type(paw_an_type),allocatable :: paw_an1(:) 406 type(paw_ij_type),allocatable :: paw_ij1(:) 407 type(pawrhoij_type),allocatable :: pawrhoijfermi(:) 408 409 ! ********************************************************************* 410 411 DBG_ENTER("COLL") 412 413 call timab(120,1,tsec) 414 call timab(154,1,tsec) 415 416 call status(0,dtfil%filstat,iexit,level,'init ') 417 418 ! intel 18 really needs this to be initialized 419 maxfor = zero 420 421 ! enable time limit handler if not done in callers. 422 if (enable_timelimit_in(ABI_FUNC) == ABI_FUNC) then 423 write(std_out,*)"Enabling timelimit check in function: ",trim(ABI_FUNC)," with timelimit: ",trim(sec2str(get_timelimit())) 424 end if 425 426 !Parallelism data 427 spaceComm=mpi_enreg%comm_cell 428 me=mpi_enreg%me_kpt 429 paral_atom=(my_natom/=dtset%natom) 430 my_atmtab=>mpi_enreg%my_atmtab 431 432 _IBM6("XLF in dfpt_scfcv") 433 434 !Save some variables from dataset definition 435 ecut=dtset%ecut 436 ecutf=ecut;if (psps%usepaw==1) ecutf=dtset%pawecutdg 437 iprcel=dtset%iprcel 438 tollist(1)=dtset%tolmxf;tollist(2)=dtset%tolwfr 439 tollist(3)=dtset%toldff;tollist(4)=dtset%toldfe 440 tollist(6)=dtset%tolvrs;tollist(7)=dtset%tolrff 441 nfftotf=product(ngfftf(1:3)) 442 nstep=dtset%nstep 443 iscf_mod=dtset%iscf 444 iscf10_mod=mod(iscf_mod,10) 445 446 qzero=0; if(qphon(1)**2+qphon(2)**2+qphon(3)**2 < tol14) qzero=1 447 448 need_fermie1=((qzero==1.and.dtset%frzfermi==0.and.nstep>0).and.& 449 & (dtset%occopt>=3.and.dtset%occopt<=8).and. & 450 & (ipert<=dtset%natom.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or.ipert==dtset%natom+5)) 451 452 !The value of iscf must be modified if ddk perturbation, see dfpt_looppert.f 453 if (ipert==dtset%natom+1.or.ipert==dtset%natom+10.or.ipert==dtset%natom+11) iscf_mod=-3 454 455 !Compute different geometric tensor, as well as ucvol, from rprimd 456 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 457 458 !Compute large sphere cut-off gsqcut 459 qphon2(:)=zero;if (psps%usepaw==1) qphon2(:)=qphon(:) 460 call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,qphon2,ngfftf) 461 462 !Some variables need to be initialized/nullify at start 463 quit=0 ; dbl_nnsclo=0 ; elast=zero; conv_retcode = -1 464 optres=merge(0,1,abs(iscf_mod)<10) 465 usexcnhat=0 466 !This might be taken away later 467 edocc=zero ; eeig0=zero ; ehart01=zero ; ehart1=zero ; ek0=zero ; ek1=zero 468 eloc0=zero ; elpsp1=zero ; enl0=zero ; enl1=zero ; eovl1=zero; exc1=zero 469 deltae=zero ; fermie1=zero ; epaw1=zero ; eberry=zero ; elmag1=zero 470 elast_mq=zero 471 dbl_nnsclo_mq=0 472 !This might be taken away later 473 edocc_mq=zero ; eeig0_mq=zero ; ehart01_mq=zero ; ehart1_mq=zero ; ek0_mq=zero ; ek1_mq=zero 474 eloc0_mq=zero ; elpsp1_mq=zero ; enl0_mq=zero ; enl1_mq=zero ; eovl1_mq=zero; exc1_mq=zero 475 deltae_mq=zero ; fermie1_mq=zero ; epaw1_mq=zero ; eberry_mq=zero ; elmag1_mq=zero 476 res2_mq=zero 477 478 !Examine tolerance criteria, and eventually print a line to the output 479 !file (with choice=1, the only non-dummy arguments of scprqt are 480 !nstep, tollist and iscf - still, diffor,res2,prtfor,fcart are here initialized to 0) 481 choice=1 ; prtfor=0 ; diffor=zero ; res2=zero 482 ABI_ALLOCATE(fcart,(3,dtset%natom)) 483 484 call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,& 485 & etotal,favg,fcart,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),& 486 & 1,iscf_mod,istep,kpt_rbz,maxfor,& 487 & mvdum,mpi_enreg,nband_rbz,nkpt_rbz,& 488 & nstep,occ_rbz,0,prtfor,0,& 489 & quit,res2,resid,residm,response,& 490 & tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode) 491 492 !Allocations/initializations for PAW only 493 if(psps%usepaw==1) then 494 usexcnhat=maxval(pawtab(:)%usexcnhat) 495 use_nhat_gga=(dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0) 496 ! 1st-order compensation density 497 ABI_ALLOCATE(nhat1,(cplex*nfftf,dtset%nspden)) 498 nhat1=zero 499 ! Projections of 1-st order WF on nl projectors 500 ABI_DATATYPE_ALLOCATE(cprj1,(dtset%natom,dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*usecprj)) 501 if (usecprj==1.and.mk1mem/=0) then 502 !cprj ordered by atom-type 503 ABI_ALLOCATE(dimcprj,(dtset%natom)) 504 call pawcprj_getdim(dimcprj,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 505 call pawcprj_alloc(cprj1,0,dimcprj) 506 ABI_DEALLOCATE(dimcprj) 507 end if 508 ! 1st-order arrays/variables related to the PAW spheres 509 ABI_DATATYPE_ALLOCATE(paw_an1,(my_natom)) 510 ABI_DATATYPE_ALLOCATE(paw_ij1,(my_natom)) 511 call paw_an_nullify(paw_an1) 512 call paw_ij_nullify(paw_ij1) 513 514 has_dijfr=0;if (ipert/=dtset%natom+1.and.ipert/=dtset%natom+10) has_dijfr=1 515 call paw_an_init(paw_an1,dtset%natom,dtset%ntypat,0,dtset%nspden,& 516 & cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1,& 517 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 518 519 call paw_ij_init(paw_ij1,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,0,dtset%natom,& 520 & dtset%ntypat,dtset%typat,pawtab,& 521 & has_dij=1,has_dijhartree=1,has_dijfr=has_dijfr,& 522 & mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom) 523 else 524 ABI_ALLOCATE(nhat1,(0,0)) 525 ABI_DATATYPE_ALLOCATE(cprj1,(0,0)) 526 ABI_DATATYPE_ALLOCATE(paw_an1,(0)) 527 ABI_DATATYPE_ALLOCATE(paw_ij1,(0)) 528 end if ! PAW 529 530 !Various allocations (potentials) 531 ABI_ALLOCATE(vhartr1,(cplex*nfftf)) 532 ABI_ALLOCATE(vtrial1,(cplex*nfftf,nspden)) 533 if(.not.kramers_deg) then 534 ABI_ALLOCATE(vhartr1_mq,(cplex*nfftf)) 535 ABI_ALLOCATE(vtrial1_pq,(cplex*nfftf,nspden)) 536 ABI_ALLOCATE(vtrial1_mq,(cplex*nfftf,nspden)) 537 end if 538 ! TODO: for non collinear case this should always be nspden, in NCPP case as well!!! 539 ABI_ALLOCATE(vxc1,(cplex*nfftf,nspden*(1-usexcnhat))) ! Not always needed 540 vtrial1_tmp => vtrial1 ! this is to avoid errors when vtrial1_tmp is unused 541 542 if (.not.kramers_deg) then 543 ABI_ALLOCATE(vxc1_mq,(cplex*nfftf,nspden*(1-usexcnhat))) 544 end if 545 546 !Several parameters and arrays for the SCF mixing: 547 !These arrays are needed only in the self-consistent case 548 if (iscf_mod>0.or.iscf_mod==-3) then 549 ABI_ALLOCATE(nvresid1,(cplex*nfftf,dtset%nspden)) 550 if (nstep==0) nvresid1=zero 551 if ((dtset%getddb .ne. 0 .or. dtset%irdddb .ne.0) .and. qzero .ne. 1) then 552 ABI_ALLOCATE(nvresid2,(cplex*nfftf,dtset%nspden)) 553 if (nstep==0) nvresid2=zero 554 end if 555 if (.not.kramers_deg) then 556 ABI_ALLOCATE(nvresid1_mq,(cplex*nfftf,dtset%nspden)) 557 if (nstep==0) nvresid1_mq=zero 558 end if 559 else 560 ABI_ALLOCATE(nvresid1,(0,0)) 561 if(.not.kramers_deg) then 562 ABI_ALLOCATE(nvresid1_mq,(0,0)) 563 end if 564 end if 565 if(nstep>0 .and. iscf_mod>0) then 566 dielar(1)=dtset%diecut;dielar(2)=dtset%dielng 567 dielar(3)=dtset%diemac;dielar(4)=dtset%diemix 568 dielar(5)=dtset%diegap;dielar(6)=dtset%dielam 569 dielar(7)=dtset%diemix;if (dtset%iscf>=10) dielar(7)=dtset%diemixmag 570 ! Additional allocation for mixing within PAW 571 npawmix=0 572 if(psps%usepaw==1) then 573 do iatom=1,my_natom 574 itypat=pawrhoij1(iatom)%itypat 575 lmn2_size=pawtab(itypat)%lmn2_size 576 pawrhoij1(iatom)%use_rhoijres=1 577 sz1=pawrhoij1(iatom)%cplex*lmn2_size;sz2=pawrhoij1(iatom)%nspden 578 ABI_ALLOCATE(pawrhoij1(iatom)%rhoijres,(sz1,sz2)) 579 do ispden=1,pawrhoij1(iatom)%nspden 580 pawrhoij1(iatom)%rhoijres(:,ispden)=zero 581 end do 582 ABI_ALLOCATE(pawrhoij1(iatom)%kpawmix,(pawtab(itypat)%lmnmix_sz)) 583 pawrhoij1(iatom)%lmnmix_sz=pawtab(itypat)%lmnmix_sz 584 pawrhoij1(iatom)%kpawmix=pawtab(itypat)%kmix 585 npawmix=npawmix+pawrhoij1(iatom)%nspden*pawtab(itypat)%lmnmix_sz*pawrhoij1(iatom)%cplex 586 end do 587 end if 588 denpot = AB7_MIXING_POTENTIAL 589 if (dtset%iscf > 10) denpot = AB7_MIXING_DENSITY 590 if (psps%usepaw==1.and.dtset%pawmixdg==0) then 591 ispmix=AB7_MIXING_FOURRIER_SPACE;nfftmix=dtset%nfft;ngfftmix(:)=dtset%ngfft(:) 592 else 593 ispmix=AB7_MIXING_REAL_SPACE;nfftmix=nfftf;ngfftmix(:)=ngfftf(:) 594 end if 595 if (iscf10_mod == 5 .or. iscf10_mod == 6) then 596 call ab7_mixing_new(mix, iscf10_mod, denpot, cplex, & 597 & nfftf, dtset%nspden, npawmix, errid, msg, dtset%npulayit) 598 else 599 call ab7_mixing_new(mix, iscf10_mod, denpot, max(cplex, ispmix), & 600 & nfftmix, dtset%nspden, npawmix, errid, msg, dtset%npulayit) 601 end if 602 if (errid /= AB7_NO_ERROR) then 603 MSG_ERROR(msg) 604 end if 605 if (dtset%mffmem == 0) then 606 call ab7_mixing_use_disk_cache(mix, dtfil%fnametmp_fft) 607 end if 608 end if ! iscf, nstep 609 610 !Here, allocate arrays for computation of susceptibility and dielectric matrix or for TDDFT 611 if( (nstep>0 .and. iscf_mod>0) .or. iscf_mod==-1 ) then 612 ! Here, for TDDFT, artificially set iprcel . Also set a variable to reduce the memory needs. 613 afford=1 614 if(iscf_mod==-1) then 615 iprcel=21 616 afford=0 617 end if 618 npwdiel=1 619 mgfftdiel=1 620 nfftdiel=1 621 ! Now, performs allocation 622 ! CAUTION : the dimensions are still those of GS, except for phnonsdiel 623 ABI_ALLOCATE(dielinv,(2,npwdiel*afford,nspden,npwdiel,nspden)) 624 ABI_ALLOCATE(susmat,(2,npwdiel*afford,nspden,npwdiel,nspden)) 625 end if 626 627 !Initialize Berry-phase related stuffs 628 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 629 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 630 ABI_ALLOCATE(pwindall,(max(mpw,mpw1)*mkmem,8,3)) 631 call dfptff_initberry(dtefield,dtset,gmet,kg,kg1,dtset%mband,mkmem,mpi_enreg,& 632 & mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,occ_rbz,pwindall,rprimd) 633 ! calculate inverse of the overlap matrix 634 ABI_ALLOCATE(qmat,(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)) 635 call qmatrix(cg,dtefield,qmat,mpw,mpw1,mkmem,dtset%mband,npwarr,nkpt,dtset%nspinor,dtset%nsppol,pwindall) 636 else 637 ABI_ALLOCATE(pwindall,(0,0,0)) 638 ABI_ALLOCATE(qmat,(0,0,0,0,0,0)) 639 end if 640 641 call timab(154,2,tsec) 642 643 !###################################################################### 644 !PERFORM ELECTRONIC ITERATIONS 645 !###################################################################### 646 647 !Offer option of computing 2nd-order total energy with existing 648 !wavefunctions when nstep<=0, else do nstep iterations 649 !Note that for non-self-consistent calculations, this loop will be exited 650 !after the first call to dfpt_vtorho 651 652 !Pass through the first routines even when nstep==0 653 !write(std_out,*) 'dfpt_scfcv, nstep=', max(1,nstep) 654 655 quitsum_request = xmpi_request_null; timelimit_exit = 0 656 657 do istep=1,max(1,nstep) 658 659 ! Handle time limit condition. 660 if (istep == 1) prev = abi_wtime() 661 if (istep > 1) then 662 now = abi_wtime() 663 wtime_step = now - prev 664 prev = now 665 call wrtout(std_out,sjoin("dfpt_scfcv: previous iteration took ",sec2str(wtime_step))) 666 667 if (have_timelimit_in(ABI_FUNC)) then 668 if (istep > 2) then 669 call xmpi_wait(quitsum_request,ierr) 670 if (quitsum_async > 0) then 671 write(msg,"(3a)")"Approaching time limit ",trim(sec2str(get_timelimit())),". Will exit istep loop in dfpt_scfcv." 672 MSG_COMMENT(msg) 673 call wrtout(ab_out, msg, "COLL") 674 timelimit_exit = 1 675 exit 676 end if 677 end if 678 679 my_quit = 0; if (now - get_start_time() + 2.15 * wtime_step > get_timelimit()) my_quit = 1 680 call xmpi_isum(my_quit,quitsum_async,spacecomm,quitsum_request,ierr) 681 end if 682 end if 683 684 ! ###################################################################### 685 ! The following steps are done once 686 ! ---------------------------------------------------------------------- 687 if (istep==1)then 688 689 ! PAW only: compute frozen part of 1st-order compensation density 690 ! and frozen part of psp strengths Dij 691 ! ---------------------------------------------------------------------- 692 if (psps%usepaw==1) then 693 optfr=0 694 idir_paw1 = idir 695 if (ipert==dtset%natom+11) then 696 call rf2_getidirs(idir,idir_dum,idir_paw1) 697 end if 698 call pawdijfr(cplex,gprimd,idir_paw1,ipert,my_natom,dtset%natom,nfftf,ngfftf,nspden,& 699 & psps%ntypat,optfr,paw_ij1,pawang,pawfgrtab,pawrad,pawtab,qphon,& 700 & rprimd,ucvol,vpsp1,vtrial,vxc,xred,& 701 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 702 703 if ((iscf_mod>=0.or.usexcnhat==0).and.(dtset%pawstgylm/=0)) then 704 ider=0;if ((ipert<=dtset%natom).and.(use_nhat_gga)) ider=1 705 call pawnhatfr(ider,idir_paw1,ipert,my_natom,dtset%natom,nspden,psps%ntypat,& 706 & pawang,pawfgrtab,pawrhoij,pawtab,rprimd,& 707 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 708 end if 709 end if 710 ! PAW only: we sometimes have to compute 1st-order compensation density 711 ! and eventually add it to density from 1st-order WFs 712 ! ---------------------------------------------------------------------- 713 nhat1grdim=0 714 ABI_ALLOCATE(nhat1gr,(0,0,0)) 715 if (psps%usepaw==1.and.ipert/=dtset%natom+1.and.ipert/=dtset%natom+10) then 716 call timab(564,1,tsec) 717 nhat1grdim=0;if (dtset%xclevel==2) nhat1grdim=usexcnhat*dtset%pawnhatxc 718 ider=2*nhat1grdim;izero=0 719 if (nhat1grdim>0) then 720 ABI_DEALLOCATE(nhat1gr) 721 ABI_ALLOCATE(nhat1gr,(cplex*nfftf,dtset%nspden,3*nhat1grdim)) 722 end if 723 call pawmknhat(dum,cplex,ider,idir_paw1,ipert,izero,gprimd,my_natom,dtset%natom,& 724 & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1,& 725 & pawrhoij1,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,& 726 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 727 if (dtfil%ireadwf/=0.and.dtset%get1den==0.and.dtset%ird1den==0.and.initialized==0) then 728 rhor1(:,:)=rhor1(:,:)+nhat1(:,:) 729 call fourdp(cplex,rhog1,rhor1(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 730 end if 731 call timab(564,2,tsec) 732 end if 733 ! Set initial guess for 1st-order potential 734 ! ---------------------------------------------------------------------- 735 call status(istep,dtfil%filstat,iexit,level,'get vtrial1 ') 736 option=1;optene=0;if (iscf_mod==-2) optene=1 737 call dfpt_rhotov(cplex,ehart01,ehart1,elpsp1,exc1,elmag1,gsqcut,idir,ipert,& 738 & dtset%ixc,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,nhat,nhat1,nhat1gr,nhat1grdim,& 739 & nkxc,nspden,n3xccc,optene,option,dtset%paral_kgb,dtset%qptn,& 740 & rhog,rhog1,rhor,rhor1,rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,& 741 & nvresid1,res2,vtrial1,vxc,vxc1,xccc3d1,dtset%ixcrot) 742 743 if(.not.kramers_deg) then 744 vtrial1_pq=vtrial1 !save trial potential at +q 745 !rhor1_mq=rhor1 746 !rhog1_mq=rhog1 747 !get initial guess for vtrial1 at -q 748 do ifft=1,nfftf 749 vtrial1_mq(2*ifft-1,1)=+vtrial1(2*ifft-1,1) 750 vtrial1_mq(2*ifft-1,2)=+vtrial1(2*ifft-1,2) 751 vtrial1_mq(2*ifft ,1)=-vtrial1(2*ifft ,1) 752 vtrial1_mq(2*ifft ,2)=-vtrial1(2*ifft ,2) 753 vtrial1_mq(2*ifft-1,3)= vtrial1(2*ifft ,4) !Re[V^12] 754 vtrial1_mq(2*ifft ,3)= vtrial1(2*ifft-1,4) !Im[V^12],see definition of v(:,4) cplex=2 case 755 vtrial1_mq(2*ifft ,4)= vtrial1(2*ifft-1,3) !Re[V^21]=Re[V^12] 756 vtrial1_mq(2*ifft-1,4)= vtrial1(2*ifft ,3) !Re[V^21]=Re[V^12] 757 end do 758 end if 759 760 ! For Q=0 and metallic occupation, initialize quantities needed to 761 ! compute the first-order Fermi energy 762 ! ---------------------------------------------------------------------- 763 if (need_fermie1) then 764 ABI_ALLOCATE(rhorfermi,(cplex*nfftf,nspden)) 765 if(.not.kramers_deg) then 766 ABI_ALLOCATE(rhorfermi_mq,(cplex*nfftf,nspden)) 767 end if 768 if (psps%usepaw==1.and.usexcnhat==0) then 769 ABI_ALLOCATE(nhatfermi,(cplex*nfftf,nspden)) 770 else 771 ABI_ALLOCATE(nhatfermi,(0,0)) 772 end if 773 ABI_DATATYPE_ALLOCATE(pawrhoijfermi,(my_natom*psps%usepaw)) 774 if (psps%usepaw==1) then 775 cplex_rhoij=max(cplex,dtset%pawcpxocc) 776 nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb) 777 call pawrhoij_alloc(pawrhoijfermi,cplex_rhoij,nspden_rhoij,dtset%nspinor,& 778 & dtset%nsppol,dtset%typat,pawtab=pawtab,mpi_atmtab=mpi_enreg%my_atmtab,& 779 & comm_atom=mpi_enreg%comm_atom) 780 end if 781 782 call dfpt_rhofermi(cg,cgq,cplex,cprj,cprjq,& 783 & doccde_rbz,docckqde,dtfil,dtset,eigenq,eigen0,eigen1,fe1fixed,gmet,gprimd,idir,& 784 & indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,mkmem,mkqmem,mk1mem,mpi_enreg,& 785 & mpw,mpw1,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,ngfftf,nhatfermi,nkpt_rbz,npwarr,npwar1,& 786 & nspden,dtset%nsppol,nsym1,occkq,occ_rbz,& 787 & paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoijfermi,pawtab,& 788 & phnons1,ph1d,dtset%prtvol,psps,rhorfermi,rmet,rprimd,symaf1,symrc1,symrl1,& 789 & ucvol,usecprj,useylmgr1,vtrial,vxc,wtk_rbz,xred,ylm,ylm1,ylmgr1) 790 if (.not.kramers_deg) then 791 call dfpt_rhofermi(cg,cg_mq,cplex,cprj,cprjq,& 792 & doccde_rbz,docckde_mq,dtfil,dtset,eigen_mq,eigen0,eigen1_mq,fe1fixed_mq,gmet,gprimd,idir,& 793 & indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1_mq,kpt_rbz,dtset%mband,mkmem,mkqmem,mk1mem,mpi_enreg,& 794 & mpw,mpw1_mq,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,ngfftf,nhatfermi,nkpt_rbz,npwarr,npwar1_mq,& 795 & nspden,dtset%nsppol,nsym1,occk_mq,occ_rbz,& 796 & paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoijfermi,pawtab,& 797 & phnons1,ph1d,dtset%prtvol,psps,rhorfermi_mq,rmet,rprimd,symaf1,symrc1,symrl1,& 798 & ucvol,usecprj,useylmgr1,vtrial,vxc,wtk_rbz,xred,ylm,ylm1,ylmgr1) 799 end if 800 ! call calcdensph(gmet,mpi_enreg,dtset%natom,nfftf,ngfftf,nspden,& 801 !& dtset%ntypat,ab_out,dtset%ratsph,rhorfermi,rprimd,dtset%typat,ucvol,xred,& 802 !& -1,cplex) !dbg 803 804 end if 805 806 end if ! End the condition of istep==1 807 808 ! ###################################################################### 809 ! The following steps are done at every iteration 810 ! ---------------------------------------------------------------------- 811 812 if (psps%usepaw==1)then 813 ! Computation of "on-site" 2nd-order energy, first-order potentials, first-order densities 814 nzlmopt=0;if (istep==2.and.dtset%pawnzlm>0) nzlmopt=-1 815 if (istep>2) nzlmopt=dtset%pawnzlm 816 call paw_an_reset_flags(paw_an1) ! Force the recomputation of on-site potentials 817 call paw_ij_reset_flags(paw_ij1,self_consistent=.true.) ! Force the recomputation of Dij 818 option=0;if (dtset%iscf>0.and.dtset%iscf<10.and.nstep>0) option=1 819 call status(istep,dtfil%filstat,iexit,level,'call pawdenpot') 820 call pawdenpot(dum,epaw1,epawdc1_dum,ipert,dtset%ixc,my_natom,dtset%natom,& 821 & dtset%nspden,psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an1,paw_an,paw_ij1,pawang,& 822 & dtset%pawprtvol,pawrad,pawrhoij1,dtset%pawspnorb,pawtab,dtset%pawxcdev,& 823 & dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp, & 824 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 825 826 ! First-order Dij computation 827 call status(istep,dtfil%filstat,iexit,level,'call pawdij ') 828 call timab(561,1,tsec) 829 if (has_dijfr>0) then 830 !vpsp1 contribution to Dij already stored in frozen part of Dij 831 ABI_ALLOCATE(vtrial1_tmp,(cplex*nfftf,nspden)) 832 vtrial1_tmp=vtrial1 833 do ispden=1,min(dtset%nspden,2) 834 vtrial1_tmp(:,ispden)=vtrial1_tmp(:,ispden)-vpsp1(:) 835 end do 836 else 837 vtrial1_tmp => vtrial1 838 end if 839 call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,dtset%natom,& 840 & nfftf,nfftotf,dtset%nspden,psps%ntypat,paw_an1,paw_ij1,pawang,& 841 & pawfgrtab,dtset%pawprtvol,pawrad,pawrhoij1,dtset%pawspnorb,pawtab,& 842 & dtset%pawxcdev,qphon,dtset%spnorbscl,ucvol,dtset%charge,vtrial1_tmp,vxc1,xred,& 843 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 844 if (has_dijfr>0) then 845 ABI_DEALLOCATE(vtrial1_tmp) 846 end if 847 848 call status(istep,dtfil%filstat,iexit,level,'call symdij ') 849 call symdij(gprimd,indsy1,ipert,my_natom,dtset%natom,nsym1,psps%ntypat,0,& 850 & paw_ij1,pawang1,dtset%pawprtvol,pawtab,rprimd,symaf1,symrc1, & 851 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,& 852 & qphon=qphon) 853 call timab(561,2,tsec) 854 end if ! end usepaw section 855 856 ! ###################################################################### 857 ! The following steps are done only when nstep>0 858 ! ---------------------------------------------------------------------- 859 call status(istep,dtfil%filstat,iexit,level,'loop istep ') 860 861 if(iscf_mod>0.and.nstep>0)then 862 write(msg, '(a,a,i4)' )ch10,' ITER STEP NUMBER ',istep 863 call wrtout(std_out,msg,'COLL') 864 end if 865 866 ! For Q=0 and metallic occupation, calculate the first-order Fermi energy 867 if (need_fermie1) then 868 call newfermie1(cplex,fermie1,fe1fixed,ipert,istep,dtset%ixc,my_natom,dtset%natom,& 869 & nfftf,nfftotf,nhatfermi,nspden,dtset%ntypat,dtset%occopt,paw_an,paw_an1,paw_ij1,pawang,& 870 & dtset%pawnzlm,pawrad,pawrhoij1,pawrhoijfermi,pawtab,dtset%pawxcdev,& 871 & dtset%prtvol,rhorfermi,ucvol,psps%usepaw,usexcnhat,vtrial1,vxc1,dtset%xclevel,& 872 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 873 if (.not.kramers_deg) then 874 !fermie1_mq is updated as well at "-q" 875 call newfermie1(cplex,fermie1_mq,fe1fixed_mq,ipert,istep,dtset%ixc,my_natom,dtset%natom,& 876 & nfftf,nfftotf,nhatfermi,nspden,dtset%ntypat,dtset%occopt,paw_an,paw_an1,paw_ij1,pawang,& 877 & dtset%pawnzlm,pawrad,pawrhoij1,pawrhoijfermi,pawtab,dtset%pawxcdev,& 878 & dtset%prtvol,rhorfermi_mq,ucvol,psps%usepaw,usexcnhat,vtrial1_mq,vxc1_mq,dtset%xclevel,& 879 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 880 end if 881 end if 882 883 ! No need to continue and call dfpt_vtorho, when nstep==0 884 if(nstep==0) exit 885 886 ! #######################e1magh############################################### 887 ! Compute the 1st-order density rho1 from the 1st-order trial potential 888 ! ---------------------------------------------------------------------- 889 890 call dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,& 891 & dbl_nnsclo,dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,dtset%qptn,edocc,& 892 & eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,enl0,enl1,fermie1,gh0c1_set,gh1c_set,& 893 & gmet,gprimd,idir,indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,& 894 & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,& 895 & nhat1,nkpt_rbz,npwarr,npwar1,res2,nspden,dtset%nsppol,nsym1,dtset%ntypat,nvresid1,& 896 & occkq,occ_rbz,optres,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,& 897 & pawrhoij1,pawtab,phnons1,ph1d,dtset%prtvol,psps,pwindall,qmat,resid,residm,rhog1,& 898 & rhor1,rmet,rprimd,symaf1,symrc1,symrl1,ucvol,usecprj,useylmgr1,ddk_f,& 899 & vtrial,vtrial1,wtk_rbz,xred,ylm,ylm1,ylmgr1) 900 if (.not.kramers_deg) then 901 rhor1_pq=rhor1 !at this stage rhor1_pq contains only one term of the 1st order density at +q 902 rhog1_pq=rhog1 !same for rhog1_pq 903 !get the second term related to 1st order wf at -q 904 call dfpt_vtorho(cg,cg_mq,cg1_mq,cg1_active_mq,cplex,cprj,cprjq,cprj1,& 905 & dbl_nnsclo_mq,dim_eig2rf,doccde_rbz,docckde_mq,dtefield,dtfil,dtset,-dtset%qptn,edocc_mq,& 906 & eeig0_mq,eigen_mq,eigen0,eigen1_mq,ek0_mq,ek1_mq,eloc0_mq,enl0_mq,enl1_mq,fermie1_mq,gh0c1_set_mq,gh1c_set_mq,& 907 & gmet,gprimd,idir,indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1_mq,kpt_rbz,dtset%mband,& 908 & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1_mq,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,& 909 & nhat1,nkpt_rbz,npwarr,npwar1_mq,res2_mq,nspden,dtset%nsppol,nsym1,dtset%ntypat,nvresid1_mq,& 910 & occk_mq,occ_rbz,optres,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,& 911 & pawrhoij1,pawtab,phnons1,ph1d,dtset%prtvol,psps,pwindall,qmat,resid_mq,residm_mq,rhog1_mq,& 912 & rhor1_mq,rmet,rprimd,symaf1,symrc1,symrl1,ucvol,usecprj,useylmgr1,ddk_f,& 913 & vtrial,vtrial1_mq,wtk_rbz,xred,ylm,ylm1,ylmgr1) 914 !reconstruct the +q and -q densities, this might bug if fft parallelization is used, todo... 915 do ifft=1,nfftf 916 rhor1_pq(2*ifft-1,:) = half*(rhor1(2*ifft-1,:)+rhor1_mq(2*ifft-1,:)) 917 rhor1_pq(2*ifft ,:) = half*(rhor1(2*ifft ,:)-rhor1_mq(2*ifft ,:)) 918 rhor1_mq(2*ifft-1,:) = rhor1_pq(2*ifft-1,:) 919 rhor1_mq(2*ifft ,:) =-rhor1_pq(2*ifft ,:) 920 end do 921 rhor1=rhor1_pq 922 call fourdp(cplex,rhog1,rhor1(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 923 call fourdp(cplex,rhog1_mq,rhor1_mq(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 924 end if 925 926 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 927 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 928 929 ! calculate \Omega E \cdot P term 930 if (ipert<=dtset%natom) then 931 ! phonon perturbation 932 call dfptff_ebp(cg,cg1,dtefield,eberry,dtset%mband,mkmem,& 933 & mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat) 934 else if (ipert==dtset%natom+2) then 935 ! electric field perturbation 936 call dfptff_edie(cg,cg1,dtefield,eberry,idir,dtset%mband,mkmem,& 937 & mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat,rprimd) 938 end if 939 end if 940 941 ! SPr: don't remove the following comments for debugging 942 ! call calcdensph(gmet,mpi_enreg,dtset%natom,nfftf,ngfftf,nspden,& 943 !& dtset%ntypat,ab_out,dtset%ratsph,rhor1,rprimd,dtset%typat,ucvol,xred,& 944 !& idir+1,cplex) 945 ! write(*,*) ' n ( 1,2)',intgden(1,1),' ',intgden(1,2) 946 ! write(*,*) ' mx( 1,2)',intgden(2,1),' ',intgden(2,2) 947 ! write(*,*) ' my( 1,2)',intgden(3,1),' ',intgden(3,2) 948 ! write(*,*) ' mz( 1,2)',intgden(4,1),' ',intgden(4,2) 949 ! call dfpt_etot(dtset%berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,& 950 !& efrloc,efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,& 951 !& enl0,enl1,epaw1,etotal,evar,evdw,exc1,elmag1,ipert,dtset%natom,optene) 952 ! write(*,*) 'SPr: ek1=',ek1,' exc1=',exc1,' elmag1=',elmag 953 ! if (ipert==dtset%natom+5) then 954 ! !calculate 1st order magnetic potential contribution to the energy 955 ! call dfpt_e1mag(e1mag,rhor1,rhog1); 956 ! endif 957 958 ! ###################################################################### 959 ! Skip out of step loop if non-SCF (completed) 960 ! ---------------------------------------------------------------------- 961 962 ! Indeed, nstep loops have been done inside dfpt_vtorho 963 if (iscf_mod<=0 .and. iscf_mod/=-3) exit 964 965 ! ###################################################################### 966 ! In case of density mixing , compute the total 2nd-order energy, 967 ! check the exit criterion, then mix the 1st-order density 968 ! ---------------------------------------------------------------------- 969 970 if (iscf_mod>=10) then 971 optene = 1 ! use double counting scheme 972 call dfpt_etot(dtset%berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,& 973 & efrloc,efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,& 974 & enl0,enl1,epaw1,etotal,evar,evdw,exc1,ipert,dtset%natom,optene) 975 call timab(152,1,tsec) 976 choice=2 977 call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,& 978 & etotal,favg,fcart,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),& 979 & 1,iscf_mod,istep,kpt_rbz,maxfor,& 980 & mvdum,mpi_enreg,nband_rbz,nkpt_rbz,& 981 & nstep,occ_rbz,0,prtfor,0,& 982 & quit,res2,resid,residm,response,& 983 & tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode) 984 call timab(152,2,tsec) 985 986 if (istep==nstep) quit=1 987 ! If criteria in scprqt say to quit, then exit the loop over istep 988 quit_sum=quit 989 call xmpi_sum(quit_sum,spaceComm,ierr) 990 991 if (quit_sum>0) exit 992 call status(istep,dtfil%filstat,iexit,level,'call newrho ') 993 ! INSERT HERE CALL TO NEWRHO3 : to be implemented 994 if (psps%usepaw==1) then 995 MSG_BUG("newrho3 not implemented: use potential mixing!") 996 end if 997 initialized=1 998 end if 999 1000 ! ###################################################################### 1001 ! Compute the new 1st-order potential from the 1st-order density 1002 ! ---------------------------------------------------------------------- 1003 1004 if (ipert<dtset%natom+10) then 1005 optene=1 1006 call status(istep,dtfil%filstat,iexit,level,'call dfpt_rhotov ') 1007 call dfpt_rhotov(cplex,ehart01,ehart1,elpsp1,exc1,elmag1,gsqcut,idir,ipert,& 1008 & dtset%ixc,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,& 1009 & nspden,n3xccc,optene,optres,dtset%paral_kgb,dtset%qptn,rhog,rhog1,rhor,rhor1,& 1010 & rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,nvresid1,res2,vtrial1,vxc,vxc1,xccc3d1,dtset%ixcrot) 1011 end if 1012 1013 ! ###################################################################### 1014 ! In case of potential mixing , compute the total 2nd-order energy, 1015 ! check the exit criterion, then mix the 1st-order potential 1016 ! ---------------------------------------------------------------------- 1017 1018 if (iscf_mod<10) then 1019 1020 ! PAW: has to compute here the "on-site" 2nd-order energy 1021 if (psps%usepaw==1) then 1022 nzlmopt=0;if (istep==1.and.dtset%pawnzlm>0) nzlmopt=-1 1023 if (istep>1) nzlmopt=dtset%pawnzlm 1024 call paw_an_reset_flags(paw_an1) ! Force the recomputation of on-site potentials 1025 option=2 1026 call pawdenpot(dum,epaw1,epawdc1_dum,ipert,dtset%ixc,my_natom,dtset%natom,dtset%nspden,& 1027 & psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an1,paw_an,paw_ij1,pawang,dtset%pawprtvol,& 1028 & pawrad,pawrhoij1,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,& 1029 & dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp,& 1030 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1031 end if 1032 1033 optene = 0 ! use direct scheme 1034 call dfpt_etot(dtset%berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,& 1035 & efrloc,efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,& 1036 !& enl0,enl1,epaw1,etotal,evar,evdw,exc1,elmag1,ipert,dtset%natom,optene) 1037 ! !debug: compute the d2E/d-qd+q energy, should be equal to the one from previous line 1038 ! if(.not.kramers_deg) then 1039 ! call dfpt_etot(dtset%berryopt,deltae_mq,eberry_mq,edocc_mq,eeig0_mq,eew,efrhar,efrkin,& 1040 !& efrloc,efrnl,efrx1,efrx2,ehart1_mq,ek0_mq,ek1_mq,eii,elast_mq,eloc0_mq,elpsp1_mq,& 1041 !& enl0_mq,enl1_mq,epaw1_mq,etotal_mq,evar_mq,evdw,exc1_mq,elmag1_mq,ipert,dtset%natom,optene) 1042 ! end if 1043 & enl0,enl1,epaw1,etotal,evar,evdw,exc1,ipert,dtset%natom,optene) 1044 1045 call timab(152,1,tsec) 1046 choice=2 1047 call status(istep,dtfil%filstat,iexit,level,'print info ') 1048 call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,& 1049 & etotal,favg,fcart,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),& 1050 & 1,iscf_mod,istep,kpt_rbz,maxfor,& 1051 & mvdum,mpi_enreg,nband_rbz,nkpt_rbz,& 1052 & nstep,occ_rbz,0,prtfor,0,& 1053 & quit,res2,resid,residm,response,& 1054 & tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode) 1055 ! !debug: print the information about residuals at "-q" 1056 ! if(.not.kramers_deg) then 1057 ! call scprqt(choice,cpus,deltae_mq,diffor,dtset,eigen0,& 1058 !& etotal_mq,favg,fcart,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),& 1059 !& 1,iscf_mod,istep,kpt_rbz,maxfor,& 1060 !& mvdum,mpi_enreg,nband_rbz,nkpt_rbz,& 1061 !& nstep,occ_rbz,0,prtfor,0,& 1062 !& quit,res2_mq,resid_mq,residm_mq,response,& 1063 !& tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode) 1064 ! endif 1065 call timab(152,2,tsec) 1066 1067 ! If criteria in scprqt say to quit, then exit the loop over istep 1068 quit_sum=quit 1069 call xmpi_sum(quit_sum,spaceComm,ierr) 1070 if (quit_sum>0) exit 1071 1072 ! TODO 1073 ! Better error handling is the SCF cycle goes bananas: 1074 ! Write a BIG warning in the output file and save the wavefunctions. 1075 ! so that we can restart. 1076 if(iscf_mod/=-3)then 1077 ! Note that nvresid1 and vtrial1 are called vresid and vtrial inside this routine 1078 call dfpt_newvtr(cplex,dbl_nnsclo,dielar,dtset,etotal,pawfgr%fintocoa,& 1079 & initialized,iscf_mod,ispmix,istep,mix,pawfgr%coatofin,& 1080 & mpi_enreg,my_natom,nfftf,nfftmix,ngfftf,ngfftmix,npawmix,pawrhoij1,& 1081 & qphon,rhor1,rprimd,psps%usepaw,nvresid1,vtrial1) 1082 if (.not.kramers_deg) then 1083 !same problem as with density reconstruction, TODO proper fft parallelization... 1084 do ifft=1,nfftf 1085 vtrial1_mq(2*ifft-1,1)=+vtrial1(2*ifft-1,1) 1086 vtrial1_mq(2*ifft-1,2)=+vtrial1(2*ifft-1,2) 1087 vtrial1_mq(2*ifft ,1)=-vtrial1(2*ifft ,1) 1088 vtrial1_mq(2*ifft ,2)=-vtrial1(2*ifft ,2) 1089 vtrial1_mq(2*ifft-1,3)= vtrial1(2*ifft ,4) !Re[V^12] 1090 vtrial1_mq(2*ifft ,3)= vtrial1(2*ifft-1,4) !Im[V^12],see definition of v(:,4) cplex=2 case 1091 vtrial1_mq(2*ifft ,4)= vtrial1(2*ifft-1,3) !Re[V^21]=Re[V^12] 1092 vtrial1_mq(2*ifft-1,4)= vtrial1(2*ifft ,3) !Re[V^21]=Re[V^12] 1093 end do 1094 end if 1095 initialized=1 1096 end if 1097 end if 1098 1099 ! ###################################################################### 1100 ! END MINIMIZATION ITERATIONS 1101 ! Note that there are different "exit" instructions within the loop 1102 ! ###################################################################### 1103 end do ! istep 1104 1105 ! Avoid pending requests if itime == ntime. 1106 call xmpi_wait(quitsum_request,ierr) 1107 if (timelimit_exit == 1) istep = istep - 1 1108 1109 !SP : Here read the _DDB file and extract the Born effective charge and 1110 ! dielectric constant. 1111 ! The idea is to supress the divergence due to a residual Born effective charge 1112 ! by renormalizing the v_hart1. For this, the difference between the ionic 1113 ! Z_kappa and the Born effective charge divided by the dielectric constant is used. 1114 ! --------------------------------------------------------------------------------- 1115 if ((dtset%getddb .ne. 0 .or. dtset%irdddb .ne.0) .and. qzero .ne. 1) then 1116 ABI_ALLOCATE(rhor2,(cplex*nfftf,nspden)) 1117 ABI_ALLOCATE(resid2,(dtset%mband*nkpt_rbz*nspden)) 1118 ABI_ALLOCATE(rhog2,(2,nfftf)) 1119 ABI_ALLOCATE(vtrial2,(cplex*nfftf,nspden)) 1120 1121 Z_kappa = nint(psps%ziontypat(dtset%typat(ipert))) ! Charge ionic from the psp 1122 qred2cart = two_pi*gprimd 1123 q_cart = MATMUL(qred2cart,qphon) 1124 q_cart = q_cart/SQRT(dot_product(q_cart,q_cart)) 1125 diel_q = dot_product(MATMUL(dielt,q_cart),q_cart) 1126 zeff_bar = SUM(zeff(:,:,:),DIM=3)/dtset%natom 1127 zeff_red = MATMUL(zeff_bar(:,:),rprimd(:,idir))/two_pi 1128 qphon_norm = SQRT(dot_product(qphon,qphon)) 1129 q_cart = MATMUL(qred2cart,qphon) 1130 born_bar = dot_product(q_cart,zeff_red(:)) 1131 zeff_red = MATMUL(zeff(:,:,ipert),rprimd(:,idir))/two_pi 1132 born = dot_product(q_cart,zeff_red(:)) 1133 1134 ! To avoid problem of divergence (0/0) we add a small value to qphon 1135 qphon2 = qphon + tol6 1136 renorm = (1-(qphon2(idir)*Z_kappa-(born-born_bar)/diel_q)/(qphon2(idir)*Z_kappa-born/diel_q)) 1137 1138 vtrial2(:,1) = vtrial1(:,1) -renorm*vhartr1 1139 1140 call dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,& 1141 & dbl_nnsclo,dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,dtset%qptn,edocc,& 1142 & eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,enl0,enl1,fermie1,gh0c1_set,gh1c_set,& 1143 & gmet,gprimd,idir,indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,& 1144 & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,& 1145 & nhat1,nkpt_rbz,npwarr,npwar1,res3,nspden,dtset%nsppol,nsym1,dtset%ntypat,nvresid2,& 1146 & occkq,occ_rbz,optres,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,& 1147 & pawrhoij1,pawtab,phnons1,ph1d,dtset%prtvol,psps,pwindall,qmat,resid2,residm2,rhog2,& 1148 & rhor2,rmet,rprimd,symaf1,symrc1,symrl1,ucvol,usecprj,useylmgr1,ddk_f,& 1149 & vtrial,vtrial2,wtk_rbz,xred,ylm,ylm1,ylmgr1,1) 1150 1151 write(msg,'(a)') ' '//NEW_LINE('A')//'& 1152 & ---------------------------------' 1153 call wrtout(ab_out,msg,'COLL') 1154 write(msg,'(a,a)')' The charge sum rule is activated'//NEW_LINE('A')//'& 1155 & ---------------------------------' 1156 call wrtout(ab_out,msg,'COLL') 1157 write(msg,'(a,i4)') ' Z_ion (psp):',Z_kappa 1158 call wrtout(ab_out,msg,'COLL') 1159 write(msg,'(a,f12.8)') ' Residual Born effective charge: ',born 1160 call wrtout(ab_out,msg,'COLL') 1161 write(msg,'(a,f12.8)') ' Renormalisation: ',renorm 1162 call wrtout(ab_out,msg,'COLL') 1163 if (renorm > 0.01 ) then 1164 write(msg,'(a,a)')' WARNING: The renormalisation seems large (> 0.01).'//NEW_LINE('A')//'& 1165 & You might consider increasing the k-point grid.' 1166 MSG_WARNING(msg) 1167 call wrtout(ab_out,msg,'COLL') 1168 end if 1169 write(msg,'(a)') ' ' 1170 call wrtout(ab_out,msg,'COLL') 1171 1172 ABI_DEALLOCATE(nvresid2) 1173 ABI_DEALLOCATE(rhor2) 1174 ABI_DEALLOCATE(resid2) 1175 ABI_DEALLOCATE(rhog2) 1176 ABI_DEALLOCATE(vtrial2) 1177 end if 1178 1179 if (iscf_mod>0.or.iscf_mod==-3) then 1180 ABI_DEALLOCATE(nvresid1) 1181 if (.not.kramers_deg) then 1182 ABI_DEALLOCATE(nvresid1_mq) 1183 end if 1184 end if 1185 1186 !###################################################################### 1187 !Additional steps after SC iterations 1188 !---------------------------------------------------------------------- 1189 1190 call timab(160,1,tsec) 1191 1192 !Compute Dynamic magnetic charges (dmc) in case of rfphon, 1193 !and magnetic susceptibility in case of rfmagn from first order density 1194 !(results to be comapred to dmc from d2e) 1195 !SPr deb 1196 !if (ipert<=dtset%natom.and.dtset%nspden>=2) then 1197 ! 1198 ! mpi_comm_sphgrid=mpi_enreg%comm_fft 1199 ! call mean_fftr(rhor1(:,1),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid) 1200 ! write(*,*) ' Mean 1st order density: ', mean_rhor1 1201 ! call mean_fftr(rhor1(:,2),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid) 1202 ! if (dtset%nspden==2) then 1203 ! write(*,*) ' 1st order m_z : ', mean_rhor1 1204 ! else !nspden==4 1205 ! write(*,*) ' 1st order m_x : ', mean_rhor1 1206 ! call mean_fftr(rhor1(:,3),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid) 1207 ! write(*,*) ' 1st order m_y : ', mean_rhor1 1208 ! call mean_fftr(rhor1(:,4),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid) 1209 ! write(*,*) ' 1st order m_z : ', mean_rhor1 1210 ! endif 1211 ! 1212 !endif 1213 1214 1215 !Eventually close the dot file, before calling dfpt_nstdy 1216 if ((ipert==dtset%natom+2.and.sum((dtset%qptn(1:3))**2)<=1.0d-7.and.& 1217 & (dtset%berryopt/=4 .and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.& 1218 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17)).or.& 1219 & ipert==dtset%natom+10.or.ipert==dtset%natom+11) then 1220 call wfk_close(ddk_f(1)) 1221 end if 1222 if ((ipert==dtset%natom+10 .and. idir>3) .or. ipert==dtset%natom+11) then 1223 call wfk_close(ddk_f(2)) 1224 end if 1225 if (ipert==dtset%natom+11) then 1226 call wfk_close(ddk_f(3)) 1227 if(idir>3) call wfk_close(ddk_f(4)) 1228 end if 1229 1230 !Deallocate the no more needed arrays 1231 if (iscf_mod>0.and.nstep>0) then 1232 call ab7_mixing_deallocate(mix) 1233 end if 1234 if( (nstep>0 .and. iscf_mod>0) .or. iscf_mod==-1 ) then 1235 ABI_DEALLOCATE(dielinv) 1236 ABI_DEALLOCATE(susmat) 1237 end if 1238 if(allocated(rhorfermi)) then 1239 ABI_DEALLOCATE(rhorfermi) 1240 end if 1241 if(allocated(rhorfermi_mq)) then 1242 ABI_DEALLOCATE(rhorfermi_mq) 1243 end if 1244 if(allocated(nhatfermi)) then 1245 ABI_DEALLOCATE(nhatfermi) 1246 end if 1247 if(allocated(pawrhoijfermi)) then 1248 call pawrhoij_free(pawrhoijfermi) 1249 ABI_DATATYPE_DEALLOCATE(pawrhoijfermi) 1250 end if 1251 if(psps%usepaw==1) then 1252 if (mk1mem/=0.and.usecprj==1) then 1253 call pawcprj_free(cprj1) 1254 end if 1255 do iatom=1,my_natom 1256 if (pawfgrtab(iatom)%nhatfr_allocated>0) then 1257 ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr) 1258 end if 1259 pawfgrtab(iatom)%nhatfr_allocated=0 1260 end do 1261 if (nstep>0.and.iscf_mod>0) then 1262 do iatom=1,my_natom 1263 pawrhoij1(iatom)%lmnmix_sz=0 1264 pawrhoij1(iatom)%use_rhoijres=0 1265 ABI_DEALLOCATE(pawrhoij1(iatom)%kpawmix) 1266 ABI_DEALLOCATE(pawrhoij1(iatom)%rhoijres) 1267 end do 1268 end if 1269 end if ! PAW 1270 ABI_DATATYPE_DEALLOCATE(cprj1) 1271 ABI_DEALLOCATE(nhat1gr) 1272 1273 call timab(160,2,tsec) 1274 call timab(150,1,tsec) 1275 1276 if (psps%usepaw==0.and.dtset%userie/=919.and. & 1277 & (ipert==dtset%natom+3.or.ipert==dtset%natom+4)) then 1278 call status(0,dtfil%filstat,iexit,level,'enter dfpt_nselt ') 1279 call dfpt_nselt(blkflg,cg,cg1,cplex,& 1280 & d2bbb,d2lo,d2nl,ecut,dtset%ecutsm,dtset%effmass_free,& 1281 & gmet,gprimd,gsqcut,idir,& 1282 & ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,dtset%mband,mgfftf,& 1283 & mkmem,mk1mem,mpert,mpi_enreg,psps%mpsang,mpw,mpw1,& 1284 & dtset%natom,nband_rbz,nfftf,ngfftf,& 1285 & nkpt_rbz,nkxc,dtset%nloalg,& 1286 & npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,& 1287 & nsym1,dtset%ntypat,occ_rbz,& 1288 & dtset%paral_kgb,ph1d,dtset%prtbbb,psps,dtset%qptn,rhog,& 1289 & rhor,rhor1,rmet,rprimd,symrc1,dtset%typat,ucvol,& 1290 & wtk_rbz,xred,ylm,ylm1,ylmgr,ylmgr1) 1291 end if 1292 1293 !Use of NSTPAW3 for NCPP (instead of DFPT_NSELT/DFPT_NSTDY) can be forced with userie=919 1294 !MT oct. 2015: this works perfectly on all automatic tests 1295 if(ipert<=dtset%natom+4)then 1296 if (psps%usepaw==1.or.dtset%userie==919) then 1297 call status(0,dtfil%filstat,iexit,level,'enter dfpt_nstpaw ') 1298 call dfpt_nstpaw(blkflg,cg,cgq,cg1,cplex,cprj,cprjq,docckqde,doccde_rbz,dtfil,dtset,d2lo,d2nl,d2ovl,& 1299 & eigenq,eigen0,eigen1,eovl1,gmet,gprimd,gsqcut,idir,indkpt1,indsy1,ipert,irrzon1,istwfk_rbz,& 1300 & kg,kg1,kpt_rbz,kxc,mgfftf,mkmem,mkqmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,ncpgr,& 1301 & nfftf,ngfftf,nhat,nhat1,nkpt_rbz,nkxc,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,& 1302 & nsym1,n3xccc,occkq,occ_rbz,paw_an,paw_an1,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrad,& 1303 & pawrhoij,pawrhoij1,pawtab,phnons1,ph1d,ph1df,psps,rhog,rhor,rhor1,rmet,rprimd,symaf1,symrc1,& 1304 & symrl1,ucvol,usecprj,psps%usepaw,usexcnhat,useylmgr1,vhartr1,vpsp1,vtrial,vtrial1,vxc,wtk_rbz,& 1305 & xccc3d1,xred,ylm,ylm1,ylmgr1) 1306 else 1307 call status(0,dtfil%filstat,iexit,level,'enter dfpt_nstdy ') 1308 if (dtset%nspden==4) then 1309 call dfpt_nstdy(atindx,blkflg,cg,cg1,cplex,dtfil,dtset,d2bbb,d2lo,d2nl,eigen0,eigen1,gmet,& 1310 & gsqcut,idir,indkpt1,indsy1,ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mkmem,mk1mem,mpert,mpi_enreg,& 1311 & mpw,mpw1,nattyp,nband_rbz,nfftf,ngfftf,nkpt,nkpt_rbz,nkxc,npwarr,npwar1,nspden,& 1312 & dtset%nsppol,nsym1,occ_rbz,ph1d,psps,rhor1,rmet,rprimd,symrc1,ucvol,& 1313 & wtk_rbz,xred,ylm,ylm1,rhor=rhor,vxc=vxc) 1314 else 1315 call dfpt_nstdy(atindx,blkflg,cg,cg1,cplex,dtfil,dtset,d2bbb,d2lo,d2nl,eigen0,eigen1,gmet,& 1316 & gsqcut,idir,indkpt1,indsy1,ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mkmem,mk1mem,mpert,mpi_enreg,& 1317 & mpw,mpw1,nattyp,nband_rbz,nfftf,ngfftf,nkpt,nkpt_rbz,nkxc,npwarr,npwar1,nspden,& 1318 & dtset%nsppol,nsym1,occ_rbz,ph1d,psps,rhor1,rmet,rprimd,symrc1,ucvol,& 1319 & wtk_rbz,xred,ylm,ylm1) 1320 end if 1321 end if 1322 end if 1323 1324 call timab(150,2,tsec) 1325 call timab(160,1,tsec) 1326 1327 !calculate Born effective charge and store it in d2lo 1328 if ((dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 1329 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17).and.& 1330 & ipert<=dtset%natom) then 1331 call dfptff_bec(cg,cg1,dtefield,dtset%natom,d2lo,idir,ipert,dtset%mband,mkmem,& 1332 & mpw,mpw1,mpert,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat,rprimd) 1333 blkflg(:,dtset%natom+2,:,1:dtset%natom)=1 1334 end if 1335 1336 !calculate dielectric tensor and store it in d2lo 1337 if ((dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 1338 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17).and.& 1339 & ipert==dtset%natom+2) then 1340 call dfptff_die(cg,cg1,dtefield,d2lo,idir,ipert,dtset%mband,mkmem,& 1341 & mpw,mpw1,mpert,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat,rprimd) 1342 blkflg(:,dtset%natom+2,:,dtset%natom+2)=1 1343 end if 1344 1345 !If SCF convergence was not reached (for nstep>0), 1346 !print a warning to the output file (non-dummy arguments: nstep, 1347 !residm, diffor - infos from tollist have been saved inside ) 1348 !Set also the value of conv_retcode 1349 choice=3 1350 call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,& 1351 & etotal,favg,fcart,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),& 1352 & 1,iscf_mod,istep,kpt_rbz,maxfor,& 1353 & mvdum,mpi_enreg,nband_rbz,nkpt_rbz,& 1354 & nstep,occ_rbz,0,prtfor,0,& 1355 & quit,res2,resid,residm,response,& 1356 & tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode) 1357 1358 !Update the content of the header (evolving variables) 1359 bantot_rbz = sum(nband_rbz(1:nkpt_rbz*dtset%nsppol)) 1360 call hdr_update(hdr,bantot_rbz,etotal,fermie,& 1361 & residm,rprimd,occ_rbz,pawrhoij1,xred,dtset%amu_orig(:,1),& 1362 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab ) 1363 1364 !Optionally provide output of charge density and/or potential in real space, 1365 !as well as analysis of geometrical factors (bond lengths and bond angles). 1366 !Warnings : 1367 !- core charge is excluded from the charge density; 1368 !- the potential is the INPUT vtrial. 1369 1370 if(ipert==dtset%natom+5.or.ipert<=dtset%natom)then 1371 prtopt=1 1372 if(ipert==dtset%natom+5) then 1373 prtopt=idir+1; 1374 call calcdensph(gmet,mpi_enreg,dtset%natom,nfftf,ngfftf,nspden,& 1375 & dtset%ntypat,ab_out,dtset%ratsph,rhor1,rprimd,dtset%typat,ucvol,xred,& 1376 & prtopt,cplex,intgden=intgden,dentot=dentot) 1377 !debug: write out the vtk first-order density components 1378 ! call appdig(pertcase,dtfil%fnameabo_den,fi1o_vtk) 1379 ! call printmagvtk(mpi_enreg,cplex,nspden,nfftf,ngfftf,rhor1,rprimd,adjustl(adjustr(fi1o_vtk)//"_PQ")) 1380 ! call printmagvtk(mpi_enreg,cplex,nspden,nfftf,ngfftf,rhor1,rprimd,adjustl(adjustr(fi1o_vtk)//"_MQ")) 1381 !SPr: add calculation of the contributions to susceptibility from all attomic spheres 1382 end if 1383 end if 1384 1385 if (iwrite_fftdatar(mpi_enreg)) then 1386 if (dtset%prtden>0) then 1387 rdwrpaw=0 1388 call appdig(pertcase,dtfil%fnameabo_den,fi1o) 1389 ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij 1390 call fftdatar_write_from_hdr("first_order_density",fi1o,dtset%iomode,hdr,& 1391 ngfftf,cplex,nfftf,dtset%nspden,rhor1,mpi_enreg) 1392 end if 1393 1394 ! first order potentials are always written because the eph code requires them 1395 ! the files are small (much much smaller that 1WFK, actually we should avoid writing 1WFK) 1396 rdwrpaw=0 1397 call appdig(pertcase,dtfil%fnameabo_pot,fi1o) 1398 ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij 1399 call fftdatar_write_from_hdr("first_order_potential",fi1o,dtset%iomode,hdr,& 1400 ngfftf,cplex,nfftf,dtset%nspden,vtrial1,mpi_enreg) 1401 1402 ! output files for perturbed potential components: vhartr1,vpsp1,vxc 1403 ! NB: only 1 spin for these 1404 if (dtset%prtvha > 0) then 1405 rdwrpaw=0 1406 ABI_ALLOCATE(vhartr1_tmp, (cplex*nfftf, dtset%nspden)) 1407 vhartr1_tmp = zero 1408 vhartr1_tmp(:,1) = vhartr1(:) 1409 call appdig(pertcase,dtfil%fnameabo_vha,fi1o) 1410 ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij 1411 call fftdatar_write_from_hdr("first_order_vhartree",fi1o,dtset%iomode,hdr,& 1412 ngfftf,cplex,nfftf,dtset%nspden,vhartr1_tmp,mpi_enreg) 1413 ABI_DEALLOCATE(vhartr1_tmp) 1414 end if 1415 1416 1417 ! vpsp1 needs to be copied to a temp array - intent(inout) in fftdatar_write_from_hdr though I do not know why 1418 ! if (dtset%prtvpsp > 0) then 1419 ! rdwrpaw=0 1420 ! call appdig(pertcase,dtfil%fnameabo_vpsp,fi1o) 1421 ! ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij 1422 ! call fftdatar_write_from_hdr("first_order_vpsp",fi1o,dtset%iomode,hdr,& 1423 ! ngfftf,cplex,nfftf,1,vpsp1,mpi_enreg) 1424 ! end if 1425 1426 if (dtset%prtvxc > 0) then 1427 rdwrpaw=0 1428 call appdig(pertcase,dtfil%fnameabo_vxc,fi1o) 1429 ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij 1430 call fftdatar_write_from_hdr("first_order_vxc",fi1o,dtset%iomode,hdr,& 1431 ngfftf,cplex,nfftf,dtset%nspden,vxc1,mpi_enreg) 1432 end if 1433 1434 1435 ! Add rhog1(G=0) to file 1436 if (mpi_enreg%me_g0 == 1) then 1437 if (dtset%iomode == IO_MODE_ETSF) then 1438 #ifdef HAVE_NETCDF 1439 NCF_CHECK(nctk_open_modify(ncid, nctk_ncify(fi1o), xmpi_comm_self)) 1440 ncerr = nctk_def_one_array(ncid, nctkarr_t('rhog1_g0', "dp", "two"), varid=varid) 1441 NCF_CHECK(ncerr) 1442 NCF_CHECK(nctk_set_datamode(ncid)) 1443 NCF_CHECK(nf90_put_var(ncid, varid, rhog1(:,1))) 1444 NCF_CHECK(nf90_close(ncid)) 1445 #endif 1446 else 1447 ! Handle Fortran files. 1448 if (open_file(fi1o, msg, newunit=ncid, form='unformatted', status='old', action="readwrite") /= 0) then 1449 MSG_ERROR(msg) 1450 end if 1451 if (fort_denpot_skip(ncid, msg) /= 0) MSG_ERROR(msg) 1452 write(ncid) rhog1(:,1) 1453 close(ncid) 1454 end if 1455 end if 1456 1457 end if ! iwrite_fftdatar(mpi_enreg) 1458 1459 !All procs waiting here... 1460 if(mpi_enreg%paral_kgb==1)then 1461 call timab(61,1,tsec) 1462 call xmpi_barrier(spaceComm) 1463 call timab(61,2,tsec) 1464 end if 1465 1466 !Deallocate arrays 1467 ABI_DEALLOCATE(fcart) 1468 ABI_DEALLOCATE(vtrial1) 1469 if (.not.kramers_deg) then 1470 ABI_DEALLOCATE(vhartr1_mq) 1471 ABI_DEALLOCATE(vxc1_mq) 1472 ABI_DEALLOCATE(vtrial1_pq) 1473 ABI_DEALLOCATE(vtrial1_mq) 1474 end if 1475 ABI_DEALLOCATE(vhartr1) 1476 ABI_DEALLOCATE(vxc1) 1477 ABI_DEALLOCATE(pwindall) 1478 ABI_DEALLOCATE(qmat) 1479 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 1480 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 1481 call destroy_efield(dtefield) 1482 if(allocated(mpi_enreg%kpt_loc2ibz_sp)) then 1483 ABI_DEALLOCATE(mpi_enreg%kpt_loc2ibz_sp) 1484 end if 1485 end if 1486 1487 if(psps%usepaw==1) then 1488 call paw_an_free(paw_an1) 1489 call paw_ij_free(paw_ij1) 1490 end if 1491 ABI_DATATYPE_DEALLOCATE(paw_an1) 1492 ABI_DATATYPE_DEALLOCATE(paw_ij1) 1493 ABI_DEALLOCATE(nhat1) 1494 1495 call status(0,dtfil%filstat,iexit,level,'exit') 1496 1497 call timab(160,2,tsec) 1498 call timab(120,2,tsec) 1499 1500 DBG_EXIT("COLL") 1501 1502 end subroutine dfpt_scfcv