TABLE OF CONTENTS
ABINIT/dfpt_nstpaw [ Functions ]
NAME
dfpt_nstpaw
FUNCTION
Initially designed for PAW approach, but works also for NCPP. This routine compute the non-stationary expression for the second derivative of the total energy, for a whole row of mixed derivatives (including diagonal terms contributing to non-stationnary 2nd-order total energy). Compared with NC-pseudopotentials, PAW contributions include: - changes of the overlap between 0-order wave-functions, - on-site contributions.
COPYRIGHT
Copyright (C) 2010-2018 ABINIT group (MT, AM) 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
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k cgq(2,mpw1*nspinor*mband*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q. cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF 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 cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q projected with non-local projectors docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree) eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree) eigen1(2*mband*mband*nkpt_rbz*nsppol)=1st-order eigenvalues at k,q (hartree) gmet(3,3)=reciprocal space metric tensor in bohr**-2. gprimd(3,3)=dimensional reciprocal space primitive translations gsqcut=cutoff on (k+G)^2 (bohr^-2) idir=direction of the 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. kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points in the reduced BZ kxc(nfftf,nkxc)=exchange and correlation kernel mgfftf=maximum size of 1D FFTs for the "fine" grid (see NOTES in respfn.F90) mkmem =number of k points treated by this node. mkqmem =number of k+q points treated by this node (GS data). mk1mem =number of k points treated by this node (RF data) mpert =maximum number of ipert mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw or 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 spin 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 ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid nhat(nfft,nspden*nhatdim)= -PAW only- compensation density nhat1(cplex*nfftf,nspden*usepaw)=1st-order compensation charge density (PAW) nkpt_rbz=number of k points in the reduced BZ for this perturbation nkxc=second dimension of the kxc array 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 nsppol=1 for unpolarized, 2 for spin-polarized nsym1=number of symmetry elements in space group consistent with i perturbation n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used otherwise, cplex*nfftf occkq(mband*nkpt_rbz*nsppol)=occupation number for each band at each k+q point of the reduced BZ occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k in the reduced BZ paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS paw_an1(natom) <type(paw_an_type)>=1st-order paw arrays given on angular mesh for the perturbation (j1) paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS paw_ij1(natom) <type(paw_ij_type)>=1st-order paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawang1 <type(pawang_type)>=pawang datastructure containing only the 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 pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data 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 psps <type(pseudopotential_type)>=variables related to pseudopotentials rhor(nfft,nspden)=array for GS electron density in electrons/bohr**3. rhor1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3. rmet(3,3)=real space metric (bohr**2) 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 ucvol=unit cell volume in bohr**3. usecprj= 1 if cprj, cprjq arrays are stored in memory usepaw= 0 for non paw calculation; =1 for paw calculation usexcnhat= -PAW only- flag controling use of compensation density in Vxc useylmgr1= 1 if ylmgr1 array is allocated vhartr1(cplex*nfft)=1-order Hartree potential vpsp1(cplex*nfftf)=first-order derivative of the ionic potential vtrial(nfftf,nspden)=GS potential (Hartree). vtrial1(cplex*nfftf,nspden)= RF 1st-order potential (Hartree). vxc(nfftf,nspden)=XC GS potential wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ 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 ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics at k+q
OUTPUT
blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed) d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs d2ovl(2,3,mpert,3,mpert*usepaw)=overlap contributions to the 2DTEs (PAW only) 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)
NOTES
delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
PARENTS
dfpt_scfcv
CHILDREN
appdig,destroy_hamiltonian,destroy_mpi_enreg,destroy_rf_hamiltonian dfpt_accrho,dfpt_atm2fft,dfpt_mkcore,dfpt_mkvxc,dfpt_mkvxc_noncoll dfpt_mkvxcstr,dfpt_sygra,dfpt_vlocal,dotprod_g,dotprod_vn,fftpac getcprj,getdc1,getgh1c,hartrestr,init_hamiltonian,init_rf_hamiltonian initmpi_seq,initylmg,kpgstr,load_k_hamiltonian,load_k_rf_hamiltonian load_kprime_hamiltonian,load_spin_hamiltonian,mkffnl,mkkin,mkkpg,occeig paw_an_reset_flags,paw_ij_free,paw_ij_init,paw_ij_nullify paw_ij_reset_flags,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawcprj_get pawdfptenergy,pawdij2e1kb,pawdijfr,pawmkrho,pawnhatfr,pawrhoij_alloc pawrhoij_free,pawrhoij_init_unpacked,pawrhoij_mpisum_unpacked,projbd stresssym,symrhg,timab,vlocalstr,wfk_close,wfk_open_read,wfk_read_bks wrtout,xmpi_barrier,xmpi_sum
SOURCE
148 #if defined HAVE_CONFIG_H 149 #include "config.h" 150 #endif 151 152 #include "abi_common.h" 153 154 subroutine dfpt_nstpaw(blkflg,cg,cgq,cg1,cplex,cprj,cprjq,docckqde,doccde_rbz,dtfil,dtset,d2lo,d2nl,d2ovl,& 155 & eigenq,eigen0,eigen1,eovl1,gmet,gprimd,gsqcut,idir,indkpt1,indsy1,ipert,irrzon1,istwfk_rbz,& 156 & kg,kg1,kpt_rbz,kxc,mgfftf,mkmem,mkqmem,mk1mem,& 157 & mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nhat1,& 158 & nkpt_rbz,nkxc,npwarr,npwar1,nspden,nspinor,nsppol,nsym1,n3xccc,occkq,occ_rbz,& 159 & paw_an,paw_an1,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,& 160 & pawrhoij1,pawtab,phnons1,ph1d,ph1df,psps,rhog,rhor,rhor1,rmet,rprimd,symaf1,symrc1,symrl1,& 161 & ucvol,usecprj,usepaw,usexcnhat,useylmgr1,vhartr1,vpsp1,vtrial,vtrial1,vxc,& 162 & wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr1) 163 164 use defs_basis 165 use defs_datatypes 166 use defs_abitypes 167 use m_xmpi 168 use m_errors 169 use m_profiling_abi 170 use m_wfk 171 use m_hamiltonian 172 use m_cgtools 173 use m_nctk 174 175 use m_io_tools, only : file_exists 176 use m_mpinfo, only : destroy_mpi_enreg 177 use m_hdr, only : hdr_skip 178 use m_occ, only : occeig 179 use m_pawang, only : pawang_type 180 use m_pawrad, only : pawrad_type 181 use m_pawtab, only : pawtab_type 182 use m_paw_an, only : paw_an_type, paw_an_reset_flags 183 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags 184 use m_pawfgrtab,only : pawfgrtab_type 185 use m_pawrhoij, only : pawrhoij_type,pawrhoij_alloc,pawrhoij_free,pawrhoij_copy,pawrhoij_nullify, & 186 & pawrhoij_init_unpacked,pawrhoij_mpisum_unpacked,pawrhoij_get_nspden 187 use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_free,pawcprj_get,pawcprj_copy 188 use m_pawdij, only : pawdijfr 189 use m_pawfgr, only : pawfgr_type 190 use m_kg, only : mkkin, kpgstr 191 use m_cgtools, only : dotprod_vn 192 193 !This section has been created automatically by the script Abilint (TD). 194 !Do not modify the following lines by hand. 195 #undef ABI_FUNC 196 #define ABI_FUNC 'dfpt_nstpaw' 197 use interfaces_14_hidewrite 198 use interfaces_18_timing 199 use interfaces_32_util 200 use interfaces_41_geometry 201 use interfaces_51_manage_mpi 202 use interfaces_53_ffts 203 use interfaces_56_recipspace 204 use interfaces_56_xc 205 use interfaces_64_psp 206 use interfaces_65_paw 207 use interfaces_66_nonlocal 208 use interfaces_66_wfs 209 use interfaces_67_common 210 use interfaces_72_response, except_this_one => dfpt_nstpaw 211 !End of the abilint section 212 213 implicit none 214 215 !Arguments ------------------------------- 216 !scalars 217 integer,intent(in) :: cplex,idir,ipert,mgfftf,mkmem,mkqmem,mk1mem,mpert,mpw,mpw1 218 integer,intent(in) :: ncpgr,nfftf,nkpt_rbz,nkxc,nspden,nspinor,nsppol,nsym1 219 integer,intent(in) :: n3xccc,usecprj,usepaw,usexcnhat,useylmgr1 220 real(dp),intent(in) :: gsqcut,ucvol 221 real(dp),intent(out) :: eovl1 222 type(datafiles_type),intent(in) :: dtfil 223 type(dataset_type),intent(in) :: dtset 224 type(MPI_type),intent(in) :: mpi_enreg 225 type(pawang_type),intent(in) :: pawang,pawang1 226 type(pawfgr_type),intent(in) :: pawfgr 227 type(pseudopotential_type),intent(in) :: psps 228 !arrays 229 integer,intent(in) :: nattyp(dtset%ntypat),nband_rbz(nkpt_rbz*nsppol) 230 integer,intent(in) :: indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom) 231 integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4)) 232 integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 233 integer,intent(in) :: ngfftf(18),npwarr(nkpt_rbz),npwar1(nkpt_rbz) 234 integer,intent(in) :: symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1) 235 integer,intent(inout) :: blkflg(3,mpert,3,mpert) 236 real(dp),intent(in) :: cg(2,mpw*nspinor*dtset%mband*mkmem*nsppol) 237 real(dp),intent(in) :: cgq(2,mpw1*nspinor*dtset%mband*mkqmem*nsppol) 238 real(dp),intent(in) :: cg1(2,mpw1*nspinor*dtset%mband*mk1mem*nsppol) 239 real(dp),intent(in) :: docckqde(dtset%mband*nkpt_rbz*nsppol) 240 real(dp),intent(in) :: doccde_rbz(dtset%mband*nkpt_rbz*nsppol) 241 real(dp),intent(in) :: eigenq(dtset%mband*nkpt_rbz*nsppol),eigen0(dtset%mband*nkpt_rbz*nsppol) 242 real(dp),intent(in) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*nsppol) 243 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt_rbz(3,nkpt_rbz) 244 real(dp),intent(in) :: kxc(nfftf,nkxc),nhat(nfftf,nspden),nhat1(cplex*nfftf,nspden*usepaw) 245 real(dp),intent(in) :: occkq(dtset%mband*nkpt_rbz*dtset%nsppol) 246 real(dp),intent(in) :: occ_rbz(dtset%mband*nkpt_rbz*nsppol) 247 real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4)) 248 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 249 real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom),rhog(2,nfftf) 250 real(dp),intent(in) :: rhor(cplex*nfftf,nspden),rhor1(cplex*nfftf,nspden),rmet(3,3),rprimd(3,3) 251 real(dp),intent(in) :: vhartr1(cplex*nfftf),vtrial1(cplex*nfftf,nspden),vxc(nfftf,nspden) 252 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,dtset%natom) 253 real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 254 real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm) 255 real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 256 real(dp),target,intent(in) :: vpsp1(cplex*nfftf),vtrial(nfftf,nspden),xccc3d1(cplex*n3xccc) 257 real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert) 258 real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert),d2ovl(2,3,mpert,3,mpert*usepaw) 259 type(pawcprj_type),intent(in) :: cprj(dtset%natom,nspinor*dtset%mband*mkmem*nsppol*usecprj) 260 type(pawcprj_type),intent(in) :: cprjq(dtset%natom,nspinor*dtset%mband*mkqmem*nsppol*usecprj) 261 type(paw_an_type),intent(in) :: paw_an(:) 262 type(paw_an_type),intent(inout) :: paw_an1(:) 263 type(paw_ij_type),intent(in) :: paw_ij(:) 264 type(paw_ij_type),intent(inout) :: paw_ij1(:) 265 type(pawfgrtab_type),intent(inout) :: pawfgrtab(:) 266 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*usepaw) 267 type(pawrhoij_type),intent(in) :: pawrhoij(:) 268 type(pawrhoij_type),intent(in) :: pawrhoij1(:) 269 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw) 270 271 !Local variables------------------------------- 272 !scalars 273 integer,parameter :: tim_fourwf=18,tim_getgh1c=2,tim_projbd=3,formeig1=1 274 integer :: bd2tot_index,bdtot_index,berryopt,bufsz,choice,cpopt,counter,cplex_rhoij,ddkcase 275 integer :: dimffnl,dimffnl1,dimffnl1_idir1,dimylmgr1 276 integer :: ia,iatom,iband,ibg,ibgq,ibg1,icg,icgq,icg1,ider,idir0,idir1,idir_cprj 277 integer :: ierr,ii,ikg,ikg1,ikpt,ikpt_me,ilmn,iorder_cprj,ipert1 278 integer :: ispden,isppol,istwf_k,istr,istr1,itypat,jband,jj,kdir1,kpert1,master,mcgq,mcprjq 279 integer :: mdir1,me,mpert1,my_natom,my_comm_atom,my_nsppol,nband_k,nband_kocc,need_ylmgr1 280 integer :: nfftot,nkpg,nkpg1,nkpt_me,npw_,npw_k,npw1_k,nspden_rhoij 281 integer :: nvh1,nvxc1,nzlmopt_ipert,nzlmopt_ipert1,optlocal,optnl 282 integer :: option,opt_gvnl1,sij_opt,spaceworld,usevnl,wfcorr,ik_ddk 283 real(dp) :: arg,doti,dotr,dot1i,dot1r,dot2i,dot2r,dot3i,dot3r,elfd_fact,invocc,lambda,wtk_k 284 logical :: force_recompute,has_dcwf,has_dcwf2,has_drho,has_ddk_file 285 logical :: is_metal,is_metal_or_qne0,need_ddk_file,need_pawij10 286 logical :: need_wfk,need_wf1,paral_atom,qne0,t_exist 287 character(len=500) :: msg 288 character(len=fnlen) :: fiwfddk(3) 289 type(gs_hamiltonian_type) :: gs_hamkq 290 type(rf_hamiltonian_type) :: rf_hamkq 291 type(MPI_type) :: mpi_enreg_seq 292 !arrays 293 integer :: ddkfil(3),my_spintab(2),nband_tmp(1),npwar1_tmp(1) 294 integer,allocatable :: jpert1(:),jdir1(:),kg1_k(:,:),kg_k(:,:) 295 integer,pointer :: my_atmtab(:) 296 real(dp) :: dum1(1,1),dum2(1,1),dum3(1,1),epawnst(2),kpoint(3),kpq(3) 297 real(dp) :: sumelfd(2),symfact(3),tsec(2),ylmgr_dum(1,3,1) 298 real(dp),allocatable :: buffer(:),ch1c(:,:,:,:),cs1c(:,:,:,:) 299 real(dp),allocatable :: cwave0(:,:),cwavef(:,:),dcwavef(:,:) 300 real(dp),allocatable :: doccde_k(:),doccde_kq(:) 301 real(dp),allocatable :: dnhat1(:,:),drhoaug1(:,:,:,:) 302 real(dp),allocatable :: drhor1(:,:),drho1wfg(:,:),drho1wfr(:,:,:) 303 real(dp),allocatable :: d2nl_elfd(:,:),dkinpw(:) 304 real(dp),allocatable :: d2nl_k(:,:),d2ovl_drho(:,:,:,:,:),d2ovl_k(:,:) 305 real(dp),allocatable :: eig_k(:),eig_kq(:),eig1_k(:) 306 real(dp),allocatable,target :: e1kbfr_spin(:,:,:,:,:),ffnlk(:,:,:,:),ffnl1(:,:,:,:) 307 real(dp),allocatable :: gh1(:,:),gs1(:,:),gvnl1(:,:),kinpw1(:),kpg_k(:,:),kpg1_k(:,:) 308 real(dp),allocatable :: occ_k(:),occ_kq(:),ph3d(:,:,:),ph3d1(:,:,:),rhotmp(:,:),rocceig(:,:) 309 real(dp),allocatable :: ylm_k(:,:),ylm1_k(:,:),ylmgr1_k(:,:,:),vtmp1(:,:),vxc10(:,:) 310 real(dp),allocatable,target :: work(:,:,:) 311 real(dp),pointer :: e1kbfr(:,:,:,:),e1kb_ptr(:,:,:) 312 real(dp),pointer :: ffnl1_idir1(:,:,:,:),vhartr01(:),vpsp1_idir1(:),xccc3d1_idir1(:) 313 type(pawcprj_type),allocatable :: dcwaveprj(:,:) 314 type(pawcprj_type),allocatable,target :: cwaveprj0(:,:) 315 type(pawcprj_type),pointer :: cwaveprj0_idir1(:,:) 316 type(paw_ij_type),allocatable :: paw_ij10(:,:) 317 type(pawrhoij_type),target,allocatable :: pawdrhoij1(:,:) 318 type(pawrhoij_type),pointer :: pawdrhoij1_unsym(:,:) 319 type(wfk_t) :: ddks(3) 320 321 ! ********************************************************************* 322 323 DBG_ENTER("COLL") 324 325 !Keep track of total time spent in dfpt_nstpaw 326 call timab(566,1,tsec) 327 328 !Not valid for PrintBandByBand 329 if (dtset%prtbbb/=0) then 330 MSG_BUG('not yet valid for prtbbb/=0!') 331 end if 332 333 !NCPP restrictions 334 if (usepaw==0) then 335 ! cprj cannot be used 336 if (usecprj/=0) then 337 MSG_BUG('NCPP: usecprj should be 0!') 338 end if 339 ! d2ovl cannot be used 340 if (size(d2ovl)/=0) then 341 MSG_BUG('NCPP: d2ovl should not be allocated!') 342 end if 343 end if 344 345 !PAW restrictions 346 if (usepaw==1) then 347 ! Test on FFT grid sizes 348 if (pawfgr%nfft/=nfftf) then 349 MSG_BUG('PAW: wrong values for nfft, nfftf!') 350 end if 351 ! Test gradients of cprj 352 if (ipert<=dtset%natom.and.ncpgr/=3) then 353 MSG_BUG('PAW: wrong value of ncpgr for ipert<=natom!') 354 end if 355 if (ipert==dtset%natom+1.and.ncpgr/=1) then 356 MSG_BUG('PAW: wrong value of ncpgr for ipert=natom+1!') 357 end if 358 if (ipert==dtset%natom+2.and.ncpgr/=3) then 359 MSG_BUG('PAW: wrong value of ncpgr for ipert=natom+2!') 360 end if 361 if ((ipert==dtset%natom+3.or.ipert==dtset%natom+4).and.ncpgr/=1) then 362 MSG_BUG('PAW: wrong value of ncpgr for ipert=natom+3 or 4!') 363 end if 364 ! Test on availability of DijHartree and XC on-site potentials 365 if (mpi_enreg%my_natom>0.and.ipert/=dtset%natom+1) then 366 if (paw_ij1(1)%has_dijhartree==0.or.paw_an1(1)%has_vxc==0) then 367 msg='PAW: paw_ij1%dijhartree and paw=_an1%vxc1 should be allocated !' 368 end if 369 end if 370 end if 371 372 !Set up parallelism 373 master=0;me=mpi_enreg%me_kpt 374 spaceworld=mpi_enreg%comm_cell 375 paral_atom=(mpi_enreg%my_natom/=dtset%natom) 376 my_comm_atom=mpi_enreg%comm_atom 377 my_natom=mpi_enreg%my_natom 378 my_atmtab=>mpi_enreg%my_atmtab 379 my_spintab=mpi_enreg%my_isppoltab 380 my_nsppol=count(my_spintab==1) 381 382 !Fake MPI data to be used in sequential calls to parallel routines 383 call initmpi_seq(mpi_enreg_seq) 384 mpi_enreg_seq%my_natom=dtset%natom 385 386 !Compute effective number of k-points 387 nkpt_me=nkpt_rbz 388 if(xmpi_paral==1)then 389 nkpt_me=0 390 do isppol=1,nsppol 391 do ikpt=1,nkpt_rbz 392 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 393 if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))) nkpt_me=nkpt_me+1 394 end do 395 end do 396 end if 397 398 !Sizes for WF at k+q 399 mcgq=mpw1*nspinor*dtset%mband*mkqmem*nsppol 400 mcprjq=nspinor*dtset%mband*mkqmem*nsppol*usecprj 401 402 !Check ddk files (needed to compute electric field perturbations) 403 ddkfil(:)=0 404 do idir1=1,3 405 ddkcase=idir1+dtset%natom*3 406 call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk(idir1)) 407 t_exist = file_exists(fiwfddk(idir1)) 408 409 if (.not. t_exist) then 410 ! Try netcdf file. 411 t_exist = file_exists(nctk_ncify(fiwfddk(idir1))) 412 if (t_exist) then 413 fiwfddk(idir1) = nctk_ncify(fiwfddk(idir1)) 414 write(msg,"(3a)")"- File: ",trim(fiwfddk(idir1))," does not exist but found netcdf file with similar name." 415 call wrtout(std_out,msg,'COLL') 416 end if 417 end if 418 419 if (t_exist) then 420 ddkfil(idir1)=20+idir1 ! Note the use of unit numbers 21, 22 and 23 421 end if 422 end do 423 has_ddk_file=(any(ddkfil(:)>0)) 424 425 !Define the set of perturbations (j1)=(ipert1,idir1) 426 !The first perturbation must be (j1)=(j2)=(ipert,idir) 427 !because we need to compute <g|H^(j2)-Eps.S^(j2)|u0> first. 428 if (ipert/=dtset%natom+1) then 429 mpert1=0 430 ABI_ALLOCATE(jpert1,(mpert)) 431 jpert1 = 0 432 if (ipert/=dtset%natom+2.or.has_ddk_file) then 433 mpert1=mpert1+1;jpert1(mpert1)=ipert 434 end if 435 do ipert1=1,mpert 436 if (ipert1/=ipert.and.& 437 & (ipert1<=dtset%natom.or.& 438 & (ipert1==dtset%natom+2.and.has_ddk_file).or.& 439 & ((ipert>dtset%natom.and.ipert/=dtset%natom+5).and.(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)).or. & 440 & ((ipert1==dtset%natom+2).and.has_ddk_file))) then 441 mpert1=mpert1+1;jpert1(mpert1)=ipert1 442 end if 443 end do 444 else 445 mpert1=1 446 ABI_ALLOCATE(jpert1,(mpert1)) 447 jpert1(1)=dtset%natom+1 448 end if 449 mdir1=3 450 ABI_ALLOCATE(jdir1,(mdir1)) 451 jdir1(1:3)= (/ (idir1,idir1=1,3) /) 452 jdir1(1)=idir;jdir1(idir)=1 453 454 !Index of strain perturbation, if any 455 istr=idir;if (ipert==dtset%natom+4) istr=idir+3 456 457 !Open ddk WF file(s) 458 if (has_ddk_file) then 459 do kdir1=1,mdir1 460 idir1=jdir1(kdir1) 461 if (ddkfil(idir1)/=0) then 462 write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk(idir1)) 463 call wrtout(std_out,msg,'COLL') 464 call wrtout(ab_out,msg,'COLL') 465 call wfk_open_read(ddks(idir1),fiwfddk(idir1),formeig1,dtset%iomode,ddkfil(idir1),spaceworld) 466 end if 467 end do 468 end if 469 470 !Zero only portion of matrix to be computed here 471 d2nl(:,:,1:dtset%natom+4,idir,ipert)=zero 472 if (usepaw==1) d2ovl(:,:,1:dtset%natom+4,idir,ipert)=zero 473 474 !Update list of computed matrix elements 475 do kpert1=1,mpert1 476 ipert1=jpert1(kpert1) 477 do kdir1=1,mdir1 478 idir1=jdir1(kdir1) 479 if ((ipert1<=dtset%natom).or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4.or.& 480 & (ipert1==dtset%natom+1.and.((ddkfil(idir1)/=0).or.(dtset%rfdir(idir1)/=0.and.idir1<=idir))).or.& 481 & (ipert1==dtset%natom+2.and.ddkfil(idir1)/=0).or.& 482 & ((ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and.ddkfil(idir1)/=0)) then 483 blkflg(idir1,ipert1,idir,ipert)=1 484 end if 485 end do 486 end do 487 488 !Initialize most of the (1st-order) Hamiltonian 489 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 490 !2) Perform the setup needed for the non-local factors: 491 call init_hamiltonian(gs_hamkq,psps,pawtab,nspinor,nsppol,nspden,dtset%natom,& 492 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,& 493 & paw_ij=paw_ij,mpi_atmtab=my_atmtab,comm_atom=my_comm_atom,mpi_spintab=mpi_enreg%my_isppoltab,& 494 & usecprj=usecprj,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda) 495 496 !Variables common to all perturbations 497 arg=maxval(occ_rbz)-minval(occ_rbz) 498 qne0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2>=tol14) 499 is_metal=((dtset%occopt>=3.and.dtset%occopt<=8).or.(abs(arg)>tol8)) 500 is_metal_or_qne0=((is_metal).or.(qne0)) 501 ABI_ALLOCATE(ch1c,(2,dtset%mband,dtset%mband,nkpt_me)) 502 ch1c(:,:,:,:)=zero 503 nzlmopt_ipert=0;nzlmopt_ipert1=0 504 if (usepaw==1) then 505 ABI_ALLOCATE(d2ovl_drho,(2,3,mpert,3,mpert)) 506 d2ovl_drho=zero 507 if (dtset%pawnzlm/=0) then 508 nzlmopt_ipert=1;if (dtset%nstep<2) nzlmopt_ipert=-1 509 nzlmopt_ipert1=-1 510 end if 511 if (is_metal_or_qne0) then 512 ABI_ALLOCATE(cs1c,(2,dtset%mband,dtset%mband,nkpt_me)) 513 cs1c(:,:,:,:)=zero 514 end if 515 end if 516 517 !Force the recomputation of on-site potentials and DijHartree 518 if (usepaw==1) then 519 call paw_an_reset_flags(paw_an1) 520 call paw_ij_reset_flags(paw_ij1,dijhartree=.true.) 521 end if 522 523 !LOOP OVER PERTURBATION TYPES (j1) 524 do kpert1=1,mpert1 525 ipert1=jpert1(kpert1) 526 527 ! Flag for use of DDK file 528 need_ddk_file=(has_ddk_file.and.(ipert1==dtset%natom+1.or.ipert1==dtset%natom+2)) 529 530 ! Factor to be applied for electric Field (Eff. charges and piezo. tensor are "minus" d2E) 531 elfd_fact=one 532 if ((ipert <=dtset%natom.or.ipert ==dtset%natom+3.or.ipert ==dtset%natom+4).and. & 533 & (ipert1==dtset%natom+2)) elfd_fact=-one 534 if ((ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and. & 535 & (ipert ==dtset%natom+2)) elfd_fact=-one 536 537 ! We want to compute delta_u^(j1))=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>] 538 ! see PRB 78, 035105 (2008), Eq. (42) 539 has_dcwf=.false.;has_dcwf2=.false.;has_drho=.false. 540 if (usepaw==1) then 541 has_dcwf =(ipert1/=dtset%natom+2) 542 has_dcwf2=(ipert /=dtset%natom+2) 543 has_drho =(has_dcwf.and.ipert1/=dtset%natom+1) 544 end if 545 546 ! Select which WF are needed 547 need_wfk=(.true.) 548 need_wf1=(.true.) 549 550 ! Initialize data for NL 1st-order (j1) hamiltonian 551 call init_rf_hamiltonian(cplex,gs_hamkq,ipert1,rf_hamkq,mpi_spintab=[0,0]) 552 553 ! The following contributions are needed only for non-DDK perturbation: 554 ! - Frozen part of 1st-order Dij 555 ! - Contribution from local potential to dynamical matrix (due to Vxc^(j1)(tild_nc)+VH^(j1)(tild_nZc)) 556 if (ipert/=dtset%natom+1.and.ipert1/=dtset%natom+1) then 557 558 ! Allocations 559 force_recompute=(usepaw==0) ! This is dangerous... 560 nfftot=ngfftf(1)*ngfftf(2)*ngfftf(3) 561 nvxc1=0;if (ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) nvxc1=nspden 562 nvh1=0;if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) nvh1=1 563 ABI_ALLOCATE(vxc10,(cplex*nfftf,nvxc1)) 564 ABI_ALLOCATE(vhartr01,(cplex*nfftf*nvh1)) 565 need_pawij10=(usepaw==1) 566 if (need_pawij10) then 567 ABI_DATATYPE_ALLOCATE(paw_ij10,(my_natom,mdir1)) 568 ABI_ALLOCATE(e1kbfr_spin,(rf_hamkq%dime1kb1,rf_hamkq%dime1kb2,nspinor**2,mdir1,my_nsppol)) 569 else 570 ABI_DATATYPE_ALLOCATE(paw_ij10,(0,0)) 571 end if 572 573 ! LOOP OVER PERTURBATION DIRECTIONS 574 do kdir1=1,mdir1 575 idir1=jdir1(kdir1) 576 istr1=idir1;if (ipert1==dtset%natom+4) istr1=idir1+3 577 578 ! Get first-order local potential and first-order pseudo core density 579 if (ipert==ipert1.and.idir==idir1.and.(.not.force_recompute)) then 580 vpsp1_idir1 => vpsp1 581 xccc3d1_idir1 => xccc3d1 582 else 583 ABI_ALLOCATE(vpsp1_idir1,(cplex*nfftf)) 584 ABI_ALLOCATE(xccc3d1_idir1,(cplex*n3xccc)) 585 if (usepaw==1) then 586 call dfpt_atm2fft(gs_hamkq%atindx,cplex,gmet,gprimd,gsqcut,istr1,ipert1,& 587 & mgfftf,psps%mqgrid_vl,dtset%natom,1,nfftf,ngfftf,dtset%ntypat,& 588 & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 589 & atmrhor1=xccc3d1_idir1,atmvlocr1=vpsp1_idir1,optv_in=1,optn_in=n3xccc/nfftf,optn2_in=1,& 590 & vspl=psps%vlspl) 591 else 592 if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 593 call vlocalstr(gmet,gprimd,gsqcut,istr1,mgfftf,mpi_enreg,psps%mqgrid_vl,dtset%natom,& 594 & nattyp,nfftf,ngfftf,dtset%ntypat,dtset%paral_kgb,ph1df,psps%qgrid_vl,ucvol,& 595 & psps%vlspl,vpsp1_idir1) 596 else 597 call dfpt_vlocal(gs_hamkq%atindx,cplex,gmet,gsqcut,idir1,ipert1,mpi_enreg,psps%mqgrid_vl,& 598 & dtset%natom,nattyp,nfftf,ngfftf,dtset%ntypat,ngfftf(1),ngfftf(2),ngfftf(3),dtset%paral_kgb,& 599 & ph1df,psps%qgrid_vl,dtset%qptn,ucvol,psps%vlspl,vpsp1_idir1,xred) 600 end if 601 if(psps%n1xccc/=0)then 602 call dfpt_mkcore(cplex,idir1,ipert1,dtset%natom,dtset%ntypat,ngfftf(1),psps%n1xccc,& 603 & ngfftf(2),ngfftf(3),dtset%qptn,rprimd,dtset%typat,ucvol,& 604 & psps%xcccrc,psps%xccc1d,xccc3d1_idir1,xred) 605 end if 606 end if 607 end if 608 609 ! Compute 1st-order non-local factors (Dij^(j1)_fr) 610 if (need_pawij10) then 611 call paw_ij_nullify(paw_ij10(:,idir1)) 612 call paw_ij_init(paw_ij10(:,idir1),cplex,nspinor,dtset%nsppol,dtset%nspden,& 613 & 0,dtset%natom,dtset%ntypat,dtset%typat,pawtab,has_dijfr=1,& 614 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 615 if (ipert/=ipert1.or.idir/=idir1.or.force_recompute) then 616 option=0 617 call pawdijfr(cplex,gprimd,idir1,ipert1,my_natom,dtset%natom,nfftf,ngfftf,& 618 & nspden,dtset%ntypat,option,paw_ij10(:,idir1),pawang,pawfgrtab,pawrad,& 619 & pawtab,dtset%qptn,rprimd,ucvol,vpsp1_idir1,vtrial,vxc,xred,& 620 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 621 else 622 do iatom=1,my_natom 623 paw_ij10(iatom,idir1)%has_dijfr=paw_ij1(iatom)%has_dijfr 624 if (paw_ij1(iatom)%has_dijfr==2) paw_ij10(iatom,idir1)%dijfr=paw_ij1(iatom)%dijfr 625 end do 626 end if 627 end if 628 629 ! Get first-order exchange-correlation potential (core-correction contribution only) 630 if (nvxc1>0) then 631 if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 632 option=0 633 call dfpt_mkvxcstr(cplex,idir1,ipert1,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,& 634 & nhat,nhat1,nkxc,nspden,n3xccc,option,dtset%paral_kgb,dtset%qptn,rhor,rhor1,rprimd,& 635 & usepaw,usexcnhat,vxc10,xccc3d1_idir1) 636 else 637 ! Non-collinear magnetism (should the second nkxc be nkxc_cur ?) 638 if (nspden==4) then 639 option=0 640 call dfpt_mkvxc_noncoll(cplex,dtset%ixc,kxc,mpi_enreg,nfftf,ngfftf,dum1,0,dum2,0,dum3,0,nkxc,& 641 & nspden,n3xccc,1,option,dtset%paral_kgb,dtset%qptn,dum1,dum1,rprimd,0,vxc,& 642 & vxc10,xccc3d1_idir1) 643 else 644 call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfftf,ngfftf,dum2,0,dum3,0,nkxc,& 645 & nspden,n3xccc,0,dtset%paral_kgb,dtset%qptn,dum1,rprimd,0,vxc10,xccc3d1_idir1) 646 end if 647 end if 648 end if 649 650 ! Get first-order Hartree potential (metric tensor contribution only) 651 if (nvh1>0) then 652 call hartrestr(gsqcut,idir1,ipert1,mpi_enreg,dtset%natom,& 653 & nfftf,ngfftf,dtset%paral_kgb,rhog,rprimd,vhartr01) 654 end if 655 656 ! Get Hartree + xc + local contributions to dynamical matrix or elastic tensor 657 if (ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 658 if (usepaw==0) then 659 ! vxc1 is integrated with the total 1st-order density (rhor1 ) 660 ! vpsp1 is integrated with the 1st-order pseudo density (rhor1) 661 call dotprod_vn(cplex,rhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol) 662 call dotprod_vn(cplex,rhor1,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol) 663 else if (usexcnhat/=0) then 664 ! vxc1 is integrated with the total 1st-order density (rhor1 including nhat1) 665 ! vpsp1 is integrated with the 1st-order pseudo density (rhor1 without nhat1) 666 ABI_ALLOCATE(rhotmp,(cplex*nfftf,1)) 667 rhotmp(:,1)=rhor1(:,1)-nhat1(:,1) 668 call dotprod_vn(cplex,rhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol) 669 call dotprod_vn(cplex,rhotmp,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol) 670 ABI_DEALLOCATE(rhotmp) 671 else 672 ! vxc1 is integrated with the 1st-order pseudo density (rhor1 without nhat1) 673 ! vpsp1 is integrated with the 1st-order pseudo density (rhor1 without nhat1) 674 ABI_ALLOCATE(rhotmp,(cplex*nfftf,nspden)) 675 rhotmp(:,:)=rhor1(:,:)-nhat1(:,:) 676 call dotprod_vn(cplex,rhotmp,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol) 677 call dotprod_vn(cplex,rhotmp,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol) 678 ABI_DEALLOCATE(rhotmp) 679 end if 680 if (nvh1>0) then 681 call dotprod_vn(cplex,rhor1,dot3r,dot3i,nfftf,nfftot,1,2,vhartr01,ucvol) 682 else 683 dot3r=zero ; dot3i=zero 684 end if 685 ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2 686 dotr=dot1r+dot2r+dot3r;doti=dot1i+dot2i+dot3i 687 ! In case ipert = natom+2, these lines compute the local part 688 ! of the Born effective charges from phonon and electric 689 ! field type perturbations, see eq. 43 of X. Gonze and C. Lee, PRB 55, 10355 (1997) 690 ! The minus sign is due to the fact that the effective charges 691 ! are minus the second derivatives of the energy 692 if (ipert/=dtset%natom+1) then 693 d2lo(1,idir1,ipert1,idir,ipert)=elfd_fact*dotr 694 d2lo(2,idir1,ipert1,idir,ipert)=elfd_fact*doti 695 end if 696 end if ! ipert1<=natom 697 698 if (ipert/=ipert1.or.idir/=idir1.or.force_recompute) then 699 ABI_DEALLOCATE(vpsp1_idir1) 700 ABI_DEALLOCATE(xccc3d1_idir1) 701 end if 702 703 ! End loop on directions 704 end do 705 706 ! Free memory 707 ABI_DEALLOCATE(vxc10) 708 ABI_DEALLOCATE(vhartr01) 709 710 else ! ddk perturbation 711 d2lo(1:2,1:mdir1,ipert1,idir,ipert)=zero 712 need_pawij10=.false. 713 ABI_DATATYPE_ALLOCATE(paw_ij10,(0,0)) 714 end if 715 716 ! Prepare RF PAW files for reading and writing if mkmem, mkqmem or mk1mem==0 717 iorder_cprj=0 718 719 ! Allocate arrays used to accumulate density change due to overlap 720 if (has_drho) then 721 ABI_ALLOCATE(drhoaug1,(cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),mdir1)) 722 ABI_ALLOCATE(drho1wfr,(cplex*dtset%nfft,dtset%nspden,mdir1)) 723 ABI_DATATYPE_ALLOCATE(pawdrhoij1,(my_natom,mdir1)) 724 if (paral_atom) then 725 ABI_DATATYPE_ALLOCATE(pawdrhoij1_unsym,(dtset%natom,mdir1)) 726 else 727 pawdrhoij1_unsym => pawdrhoij1 728 end if 729 do kdir1=1,mdir1 730 idir1=jdir1(kdir1) 731 drho1wfr(:,:,idir1)=zero 732 cplex_rhoij=max(cplex,dtset%pawcpxocc) 733 nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb) 734 call pawrhoij_alloc(pawdrhoij1(:,idir1),cplex_rhoij,nspden_rhoij,dtset%nspinor,& 735 & dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=1,use_rhoij_=0,& 736 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=my_atmtab) 737 if (paral_atom) then 738 call pawrhoij_alloc(pawdrhoij1_unsym(:,idir1),cplex_rhoij,nspden_rhoij,dtset%nspinor,& 739 & dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1) 740 else 741 call pawrhoij_init_unpacked(pawdrhoij1_unsym(:,idir1)) 742 end if 743 end do 744 end if 745 746 ! Initialize shifts for global arrays 747 bdtot_index=0 748 bd2tot_index=0 749 ibg=0;icg=0 750 ibg1=0;icg1=0 751 ibgq=0;icgq=0 752 753 ! Has to 1st-order non-local factors before the loop over spins 754 ! because this needs a communication over comm_atom (=comm_spinkpt) 755 if (need_pawij10) then 756 if (my_nsppol<nsppol) then 757 ABI_ALLOCATE(work,(rf_hamkq%dime1kb1,rf_hamkq%dime1kb2,nspinor**2)) 758 end if 759 ii=0 760 do isppol=1,nsppol 761 if (my_spintab(isppol)==1) ii=ii+1 762 if (my_spintab(isppol)/=1) e1kb_ptr => work 763 do kdir1=1,mdir1 764 idir1=jdir1(kdir1) 765 if (my_spintab(isppol)==1) e1kb_ptr => e1kbfr_spin(:,:,:,idir1,ii) 766 call pawdij2e1kb(paw_ij10(:,idir1),isppol,my_comm_atom,e1kbfr=e1kb_ptr,mpi_atmtab=my_atmtab) 767 end do 768 end do 769 if (my_nsppol<nsppol) then 770 ABI_DEALLOCATE(work) 771 end if 772 end if 773 774 ! LOOP OVER SPINS 775 do isppol=1,nsppol 776 777 ! either this has to be inside the sppol loop or the size of ch1c and cs1c has to be nkpt_me*nsppol 778 ikpt_me=0 779 780 ! Rewind (k+G) data if needed 781 ikg=0;ikg1=0 782 783 ! Continue to initialize the GS/RF Hamiltonian 784 call load_spin_hamiltonian(gs_hamkq,isppol,with_nonlocal=.true.) 785 if (need_pawij10) then 786 ii=min(isppol,size(e1kbfr_spin,5)) 787 if (ii>0) e1kbfr => e1kbfr_spin(:,:,:,:,ii) 788 end if 789 790 ! Initialize accumulation of density 791 if (has_drho) drhoaug1(:,:,:,:)=zero 792 793 ! LOOP OVER K-POINTS 794 do ikpt=1,nkpt_rbz 795 796 ! Load dimensions for this k-point 797 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 798 istwf_k=istwfk_rbz(ikpt) 799 npw_k=npwarr(ikpt) 800 npw1_k=npwar1(ikpt) 801 802 ! Skip loop if this k-point is not to be treated by this proc 803 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 804 bdtot_index=bdtot_index+nband_k 805 bd2tot_index=bd2tot_index+2*nband_k**2 806 cycle ! Skip the rest of the k-point loop 807 end if 808 809 ! Allocate/initialize local arrays and scalars for this k-point 810 ABI_ALLOCATE(d2nl_k,(2,3)) 811 ABI_ALLOCATE(d2ovl_k,(2,3)) 812 ABI_ALLOCATE(eig_k,(nband_k)) 813 ABI_ALLOCATE(eig_kq,(nband_k)) 814 ABI_ALLOCATE(eig1_k,(2*nband_k**2)) 815 ABI_ALLOCATE(occ_k,(nband_k)) 816 d2nl_k(:,:)=zero;d2ovl_k(:,:)=zero 817 eig_k (:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 818 eig_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index) 819 eig1_k(:)=eigen1(1+bd2tot_index:2*nband_k**2+bd2tot_index) 820 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 821 nband_kocc=count(abs(occ_k(:))>tol8) 822 kpoint(:)=kpt_rbz(:,ikpt) 823 kpq(:)=kpoint(:);if (ipert1<dtset%natom+3) kpq(:)=kpq(:)+dtset%qptn(1:3) 824 wtk_k=wtk_rbz(ikpt) 825 need_ylmgr1=0;dimylmgr1=0 826 nkpg=0;nkpg1=0 827 ikpt_me=ikpt_me+1 828 if (is_metal) then 829 ! For each pair of active bands (m,n), generates the ratios 830 ! rocceig(m,n)=(occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)) 831 ABI_ALLOCATE(doccde_k,(nband_k)) 832 ABI_ALLOCATE(doccde_kq,(nband_k)) 833 ABI_ALLOCATE(occ_kq,(nband_k)) 834 ABI_ALLOCATE(rocceig,(nband_k,nband_k)) 835 doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index) 836 doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index) 837 occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index) 838 call occeig(doccde_k,doccde_kq,eig_k,eig_kq,nband_k,dtset%occopt,occ_k,occ_kq,rocceig) 839 end if 840 841 ! Take care of the npw and kg records in WF and DDK files 842 if (need_ddk_file) then 843 do kdir1=1,mdir1 844 idir1=jdir1(kdir1) 845 if (ddkfil(idir1)/=0)then 846 !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt) 847 ik_ddk = indkpt1(ikpt) 848 npw_ = ddks(idir1)%hdr%npwarr(ik_ddk) 849 if (npw_/=npw_k) then 850 write(unit=msg,fmt='(a,i3,a,i5,a,i3,a,a,i5,a,a,i5)')& 851 & 'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',idir,ch10,& 852 & 'the number of plane waves in the ddk file is equal to', npw_,ch10,& 853 & 'while it should be ',npw_k 854 MSG_ERROR(msg) 855 end if 856 857 end if 858 end do 859 end if 860 861 ! Allocate arrays used for NL form factors 862 ABI_ALLOCATE(kg_k,(3,npw_k)) 863 ABI_ALLOCATE(kg1_k,(3,npw1_k)) 864 ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 865 ABI_ALLOCATE(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)) 866 if (psps%useylm==1.and.(need_ddk_file.or.ipert1==dtset%natom+1)) need_ylmgr1=1 867 if (psps%useylm==1.and.(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)) need_ylmgr1=1 868 dimylmgr1=max(useylmgr1,need_ylmgr1) 869 ABI_ALLOCATE(ylmgr1_k,(npw1_k,3,psps%mpsang*psps%mpsang*psps%useylm*dimylmgr1)) 870 871 ! Get plane-wave vectors and related data at k 872 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 873 if (psps%useylm==1) then 874 do jj=1,psps%mpsang*psps%mpsang 875 ylm_k(1:npw_k,jj)=ylm(1+ikg:npw_k+ikg,jj) 876 end do 877 end if 878 879 ! Get plane-wave vectors and related data at k+q 880 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 881 if (psps%useylm==1) then 882 do jj=1,psps%mpsang*psps%mpsang 883 ylm1_k(1:npw1_k,jj)=ylm1(1+ikg1:npw1_k+ikg1,jj) 884 end do 885 if (need_ylmgr1==1.and.useylmgr1/=0) then 886 do jj=1,psps%mpsang*psps%mpsang 887 do ia=1,3 888 ylmgr1_k(1:npw1_k,ia,jj)=ylmgr1(1+ikg1:npw1_k+ikg1,ia,jj) 889 end do 890 end do 891 end if 892 end if 893 894 ! If Ylm gradients at k+q are needed and not in memory, compute them 895 if (need_ylmgr1==1.and.useylmgr1==0) then 896 option=-1;npwar1_tmp(1)=npw1_k;nband_tmp(1)=nband_k 897 !Subtlety: initylmg is called in sequential mode 898 call initylmg(gprimd,kg1_k,kpq,1,mpi_enreg_seq,psps%mpsang,& 899 & npw1_k,nband_tmp,1,npwar1_tmp,nsppol,option,rprimd,& 900 & ylm1_k,ylmgr1_k) 901 end if 902 903 ! Compute (k+G) vectors 904 nkpg=0;if(ipert1<=dtset%natom) nkpg=3*dtset%nloalg(3) 905 ABI_ALLOCATE(kpg_k,(npw_k,nkpg)) 906 if (nkpg>0) then 907 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 908 end if 909 910 ! Compute (k+q+G) vectors 911 nkpg1=0;if(ipert1<=dtset%natom.or.need_ylmgr1==1) nkpg1=3*dtset%nloalg(3) 912 ABI_ALLOCATE(kpg1_k,(npw1_k,nkpg1)) 913 if (nkpg1>0) then 914 call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k) 915 end if 916 917 ! Allocate kinetic contributions 918 ABI_ALLOCATE(dkinpw,(npw_k)) 919 ABI_ALLOCATE(kinpw1,(npw1_k)) 920 dkinpw=zero 921 ! Compute (1/2) (2 Pi)**2 (k+q+G)**2: 922 ! call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k) 923 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0) 924 925 ! Compute nonlocal form factors ffnl at (k+G), for all atoms 926 ider=0;idir0=0 927 dimffnl=0;if (ipert1<=dtset%natom) dimffnl=1 928 ABI_ALLOCATE(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 929 if (ipert1<=dtset%natom) then 930 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,& 931 & ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,& 932 & psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,usepaw,& 933 & psps%useylm,ylm_k,ylmgr_dum) 934 end if 935 936 ! Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms 937 ider=1;if (ipert1<=dtset%natom) ider=0 938 if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)then 939 dimffnl1=1;if (ider>=1) dimffnl1=2+5*psps%useylm 940 idir0=0;if (ider>0.and.psps%useylm==1) idir0=-7 941 else 942 dimffnl1=1;if (ider>=1) dimffnl1=2+2*psps%useylm 943 idir0=0;if (ider>0.and.psps%useylm==1) idir0=4 944 end if 945 ABI_ALLOCATE(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat)) 946 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,& 947 & ider,idir0,psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,& 948 & psps%mqgrid_ff,nkpg1,npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,usepaw,& 949 & psps%useylm,ylm1_k,ylmgr1_k) 950 951 ! Extract non-local form factors for H^(j1) 952 if (ipert1<=dtset%natom) then 953 ffnl1_idir1 => ffnl1(:,:,:,:) 954 dimffnl1_idir1=dimffnl1 955 else 956 dimffnl1_idir1=1+ider 957 ABI_ALLOCATE(ffnl1_idir1,(npw1_k,dimffnl1_idir1,psps%lmnmax,psps%ntypat)) 958 ii=1;if (psps%useylm==0) ii=1+ider 959 do itypat=1,psps%ntypat 960 do ilmn=1,psps%lmnmax 961 ffnl1_idir1(1:npw1_k,1:ii,ilmn,itypat)=ffnl1(1:npw1_k,1:ii,ilmn,itypat) 962 end do 963 end do 964 end if 965 966 ! Load k-dependent part in the Hamiltonian datastructure 967 ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamkq%matblk)) 968 call load_k_hamiltonian(gs_hamkq,kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,& 969 & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.) 970 if (size(ffnlk)>0) then 971 call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnlk) 972 else 973 call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnl1) 974 end if 975 976 ! Load k+q-dependent part in the Hamiltonian datastructure 977 ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead) 978 call load_kprime_hamiltonian(gs_hamkq,kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,& 979 & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1_idir1,& 980 & compute_gbound=.true.) 981 if (qne0) then 982 ABI_ALLOCATE(ph3d1,(2,npw1_k,gs_hamkq%matblk)) 983 call load_kprime_hamiltonian(gs_hamkq,ph3d_kp=ph3d1,compute_ph3d=.true.) 984 end if 985 986 ! Load k-dependent part in the 1st-order Hamiltonian datastructure 987 call load_k_rf_hamiltonian(rf_hamkq,npw_k=npw_k,dkinpw_k=dkinpw) 988 989 ! Allocate memory space for one band 990 if (need_wfk) then 991 ABI_ALLOCATE(cwave0,(2,npw_k*nspinor)) 992 end if 993 if (need_wf1) then 994 ABI_ALLOCATE(cwavef,(2,npw1_k*nspinor)) 995 end if 996 ABI_ALLOCATE(gh1,(2,npw1_k*nspinor)) 997 nullify(cwaveprj0_idir1) 998 if (usecprj==1) then 999 ABI_DATATYPE_ALLOCATE(cwaveprj0,(dtset%natom,nspinor)) 1000 call pawcprj_alloc(cwaveprj0,ncpgr,gs_hamkq%dimcprj) 1001 end if 1002 if (has_dcwf) then 1003 ABI_ALLOCATE(gs1,(2,npw1_k*nspinor)) 1004 else 1005 ABI_ALLOCATE(gs1,(0,0)) 1006 end if 1007 1008 ! LOOP OVER BANDS 1009 do iband=1,nband_k 1010 1011 ! Skip band if not to be treated by this proc 1012 if (xmpi_paral==1) then 1013 if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle 1014 end if 1015 1016 ! Extract GS wavefunctions 1017 if (need_wfk) then 1018 cwave0(:,:)=cg(:,1+(iband-1)*npw_k*nspinor+icg:iband*npw_k*nspinor+icg) 1019 if (usecprj==1) then 1020 call pawcprj_get(gs_hamkq%atindx1,cwaveprj0,cprj,dtset%natom,iband,ibg,ikpt,iorder_cprj,& 1021 & isppol,dtset%mband,mkmem,dtset%natom,1,nband_k,nspinor,nsppol,dtfil%unpaw,& 1022 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 1023 end if 1024 end if 1025 1026 ! Extract 1st-order wavefunctions 1027 if (need_wf1) then 1028 cwavef(:,:)=cg1(:,1+(iband-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1) 1029 end if 1030 1031 ! LOOP OVER PERTURBATION DIRECTIONS 1032 do kdir1=1,mdir1 1033 idir1=jdir1(kdir1) 1034 istr1=idir1;if(ipert1==dtset%natom+4) istr1=idir1+3 1035 1036 ! Not able to compute if ipert1=(Elect. field) and no ddk WF file 1037 if (ipert1==dtset%natom+2.and.ddkfil(idir1)==0) cycle 1038 1039 ! Extract 1st-order NL form factors derivatives for this idir1 1040 if (dimffnl1_idir1>=2.and.psps%useylm==1.and.ipert1>dtset%natom) then 1041 do itypat=1,psps%ntypat 1042 do ilmn=1,psps%lmnmax 1043 ffnl1_idir1(1:npw1_k,2,ilmn,itypat)=ffnl1(1:npw1_k,1+istr1,ilmn,itypat) 1044 end do 1045 end do 1046 end if 1047 1048 ! Extract ground state projected WF and derivatives in idir1 direction 1049 if (usecprj==1) then 1050 cpopt=1 1051 1052 ! === Atomic displ. perturbation 1053 if (ipert<=dtset%natom) then 1054 ABI_DATATYPE_ALLOCATE(cwaveprj0_idir1,(dtset%natom,nspinor)) 1055 call pawcprj_alloc(cwaveprj0_idir1,1,gs_hamkq%dimcprj) 1056 if (ipert1<=dtset%natom) then 1057 call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=idir1) 1058 else 1059 if (ipert1==dtset%natom+2) then 1060 idir_cprj=idir1;choice=5 1061 end if 1062 if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 1063 idir_cprj=istr1;choice=3 1064 end if 1065 call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1) 1066 call getcprj(choice,cpopt,cwave0,cwaveprj0_idir1,& 1067 & gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,& 1068 & gs_hamkq%kg_kp,gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,& 1069 & gs_hamkq%mgfft,mpi_enreg,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,& 1070 & gs_hamkq%nloalg,gs_hamkq%npw_kp,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,& 1071 & gs_hamkq%ph1d,gs_hamkq%ph3d_kp,gs_hamkq%ucvol,gs_hamkq%useylm) 1072 end if 1073 1074 ! === Wave-vector perturbation 1075 else if (ipert==dtset%natom+1) then 1076 cwaveprj0_idir1 => cwaveprj0 1077 1078 ! == Electric field perturbation 1079 else if (ipert==dtset%natom+2) then 1080 ABI_DATATYPE_ALLOCATE(cwaveprj0_idir1,(dtset%natom,nspinor)) 1081 call pawcprj_alloc(cwaveprj0_idir1,1,gs_hamkq%dimcprj) 1082 if (ipert1==dtset%natom+2) then 1083 call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=idir1) 1084 else 1085 if (ipert1<=dtset%natom) then 1086 idir_cprj=idir1;choice=2 1087 end if 1088 if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 1089 idir_cprj=istr1;choice=3 1090 end if 1091 call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1) 1092 call getcprj(choice,cpopt,cwave0,cwaveprj0_idir1,& 1093 & gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,& 1094 & gs_hamkq%kg_kp,gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,& 1095 & gs_hamkq%mgfft,mpi_enreg,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,gs_hamkq%nloalg,& 1096 & gs_hamkq%npw_kp,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,gs_hamkq%ph1d,& 1097 & gs_hamkq%ph3d_kp,gs_hamkq%ucvol,gs_hamkq%useylm) 1098 end if 1099 1100 ! === Strain perturbation 1101 else if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then 1102 if ((ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and.(istr==istr1)) then 1103 cwaveprj0_idir1 => cwaveprj0 1104 else 1105 ABI_DATATYPE_ALLOCATE(cwaveprj0_idir1,(dtset%natom,nspinor)) 1106 call pawcprj_alloc(cwaveprj0_idir1,ncpgr,gs_hamkq%dimcprj) 1107 if (ipert1<=dtset%natom) then 1108 idir_cprj=idir1;choice=2 1109 end if 1110 if (ipert1==dtset%natom+2) then 1111 idir_cprj=idir1;choice=5 1112 end if 1113 if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 1114 idir_cprj=istr1;choice=3 1115 end if 1116 call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1) 1117 call getcprj(choice,cpopt,cwave0,cwaveprj0_idir1,& 1118 & gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,gs_hamkq%kg_kp,& 1119 & gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,gs_hamkq%mgfft,mpi_enreg,& 1120 & gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,gs_hamkq%nloalg,gs_hamkq%npw_kp,& 1121 & gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,gs_hamkq%ph1d,gs_hamkq%ph3d_kp,& 1122 & gs_hamkq%ucvol,gs_hamkq%useylm) 1123 end if 1124 end if ! ipert 1125 1126 else ! usecprj=0: cwaveprj0_idir1 is not used 1127 cwaveprj0_idir1 => cwaveprj0 1128 end if 1129 1130 ! Eventually compute 1st-order kinetic operator 1131 if (ipert1==dtset%natom+1) then 1132 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir1,0) 1133 else if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 1134 call kpgstr(dkinpw,dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,gprimd,istr1,kg_k,kpoint,npw_k) 1135 end if 1136 1137 ! Finalize initialization of 1st-order NL hamiltonian 1138 if (need_pawij10) rf_hamkq%e1kbfr => e1kbfr(:,:,:,idir1) 1139 1140 ! Read DDK wave function (if ipert1=electric field) 1141 if (ipert1==dtset%natom+2) then 1142 usevnl=1 1143 ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor*usevnl)) 1144 if (need_ddk_file) then 1145 if (ddkfil(idir1)/=0) then 1146 !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt) 1147 ik_ddk = indkpt1(ikpt) 1148 call wfk_read_bks(ddks(idir1), iband, ik_ddk, isppol, xmpio_single, cg_bks=gvnl1) 1149 else 1150 gvnl1=zero 1151 end if 1152 if (ipert1==dtset%natom+2) then 1153 do ii=1,npw1_k*nspinor ! Multiply ddk by +i (to be consistent with getgh1c) 1154 arg=gvnl1(1,ii) 1155 gvnl1(1,ii)=-gvnl1(2,ii) 1156 gvnl1(2,ii)=arg 1157 end do 1158 end if 1159 else 1160 gvnl1=zero 1161 end if 1162 else 1163 usevnl=0 1164 ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor*usevnl)) 1165 end if 1166 1167 ! Get |H^(j2)-Eps_k_i.S^(j2)|u0_k_i> (VHxc-dependent part not taken into account) and S^(j2)|u0> 1168 lambda=eig_k(iband);berryopt=0;optlocal=0 1169 optnl=0;if (ipert1/=dtset%natom+1.or.idir==idir1) optnl=1 1170 opt_gvnl1=0;if (ipert1==dtset%natom+2) opt_gvnl1=2 1171 sij_opt=-1;if (has_dcwf) sij_opt=1 1172 if (usepaw==0) sij_opt=0 1173 call getgh1c(berryopt,cwave0,cwaveprj0_idir1,gh1,dum1,gs1,gs_hamkq,gvnl1,idir1,ipert1,& 1174 & lambda,mpi_enreg,optlocal,optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 1175 if (sij_opt==1.and.optnl==1) gh1=gh1-lambda*gs1 1176 ABI_DEALLOCATE(gvnl1) 1177 1178 ! If needed, compute here <delta_u^(j1)_k_i|H^(j2)-Eps_k_i.S^(j2)|u0_k_i> 1179 ! with delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>] 1180 ! (see PRB 78, 035105 (2008), Eq. (42)) 1181 ! This can be rewritten as: 1182 ! -1/2.<u0_k_i|S^(j1)| Sum_{j}[<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>.|u0_k+q_j> 1183 ! The sum over j can be computed with a single call to projbd routine 1184 ! At first call (when j1=j2), ch1c=<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i> is stored 1185 ! For the next calls, it is re-used. 1186 if (has_dcwf.or.(ipert==ipert1.and.idir==idir1.and.usepaw==1)) then 1187 ! note: gvnl1 used as temporary space 1188 ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor)) 1189 if (ipert==ipert1.and.idir==idir1) then 1190 option=0;gvnl1=gh1 1191 else 1192 option=1;gvnl1=zero 1193 end if 1194 ! Compute -Sum_{j}[<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>.|u0_k+q_j> 1195 call projbd(cgq,gvnl1,-1,icgq,0,istwf_k,mcgq,0,nband_k,npw1_k,nspinor,& 1196 & dum1,ch1c(:,1:nband_k,iband,ikpt_me),option,tim_projbd,0,mpi_enreg%me_g0,mpi_enreg%comm_fft) 1197 if (has_dcwf) then 1198 if (ipert==ipert1.and.idir==idir1) gvnl1=gvnl1-gh1 1199 dotr=zero;doti=zero 1200 if (abs(occ_k(iband))>tol8) then 1201 ! Compute: -<u0_k_i|S^(j1)| Sum_{j}[<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>.|u0_k+q_j> 1202 call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,gs1,gvnl1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1203 ! Add contribution to DDB 1204 ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2 1205 ! (-1) factor already present 1206 d2ovl_k(1,idir1)=d2ovl_k(1,idir1)+wtk_k*occ_k(iband)*elfd_fact*dotr 1207 d2ovl_k(2,idir1)=d2ovl_k(2,idir1)+wtk_k*occ_k(iband)*elfd_fact*doti 1208 end if 1209 end if 1210 ABI_DEALLOCATE(gvnl1) 1211 end if 1212 1213 ! If needed, compute here <delta_u^(j1)_k_i|H-Eps_k_i.S|u^(j2)_k_i> 1214 ! This is equal to <delta_u^(j1)_k_i|H-Eps_k_i.S|delta_u^(j2)_k_i> (I) 1215 ! +<delta_u^(j1)_k_i|H-Eps_k_i.S|u^paral^(j2)_k_i> (II) 1216 ! (u^paral^(j2)_k_i is the part of u^(j2)_k_i parallel to active space : metals) 1217 ! (I) can be rewritten as: 1218 ! Sum_j{ 1/4.<u0_k_i|S^(j1)|u0_k+q_j>.<u0_k+q_j|S^(j2)|u0_k_i>.(Eps_k+q_j-Eps_k_i) } 1219 ! (II) can be rewritten as: 1220 ! Sum_j{1/2.(occ_kq_j-occ_k_i).Eps1_k,q_ij.<u0_k_i|S^(j1)|u0_k+q_j> } 1221 ! where Eps1_k,q_ij=<u0_k+q_j|H^(j2)-1/2(Eps_k+q_j-Eps_k_i)S^(j2)|u0_k_i> 1222 ! At first call (when j1=j2), cs1c=<u0_k_i|S^(j1)|u0_k+q_j> is stored 1223 ! For the next calls, it is re-used. 1224 if (has_dcwf.and.is_metal_or_qne0) then 1225 ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor)) 1226 dotr=zero;doti=zero 1227 invocc=zero ; if (abs(occ_k(iband))>tol8) invocc=two/occ_k(iband) 1228 do jband=1,nband_k 1229 ! Computation of cs1c=<u0_k_i|S^(j1)|u0_k+q_j> 1230 if ((ipert==ipert1.and.idir==idir1).or.(abs(occ_k(iband))>tol8)) then 1231 gvnl1(:,1:npw1_k*nspinor)=cgq(:,1+npw1_k*nspinor*(jband-1)+icgq:npw1_k*nspinor*jband+icgq) 1232 call dotprod_g(dot1r,dot1i,istwf_k,npw1_k*nspinor,2,gs1,gvnl1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1233 if (ipert==ipert1.and.idir==idir1.and.has_dcwf2) then 1234 cs1c(1,jband,iband,ikpt_me)=dot1r 1235 cs1c(2,jband,iband,ikpt_me)=dot1i 1236 end if 1237 end if 1238 if (abs(occ_k(iband))>tol8) then 1239 ! Computation of term (I) 1240 if (has_dcwf2) then 1241 arg=eig_kq(jband)-eig_k(iband) 1242 dot2r=cs1c(1,jband,iband,ikpt_me) 1243 dot2i=cs1c(2,jband,iband,ikpt_me) 1244 dotr=dotr+(dot1r*dot2r+dot1i*dot2i)*arg 1245 doti=doti+(dot1i*dot2r-dot1r*dot2i)*arg 1246 end if 1247 ! Computation of term (II) 1248 if (is_metal) then 1249 if (abs(rocceig(jband,iband))>tol8) then 1250 ii=2*jband-1+(iband-1)*2*nband_k 1251 arg=invocc*rocceig(jband,iband)*(eig_k(iband)-eig_kq(jband)) 1252 dot2r=eig1_k(ii);dot2i=eig1_k(ii+1) 1253 dotr=dotr+arg*(dot1r*dot2r-dot1i*dot2i) 1254 doti=doti+arg*(dot1r*dot2i+dot2r*dot1i) 1255 end if 1256 end if 1257 end if 1258 end do 1259 dotr=quarter*dotr;doti=quarter*doti 1260 ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) 1261 d2ovl_k(1,idir1)=d2ovl_k(1,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*dotr 1262 d2ovl_k(2,idir1)=d2ovl_k(2,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*doti 1263 ABI_DEALLOCATE(gvnl1) 1264 end if 1265 1266 ! Build the matrix element <u0_k_i|H^(j1)-Eps_k_i.S^(j1)|u^(j2)_k,q_i> 1267 ! and add contribution to DDB 1268 if (ipert1/=dtset%natom+1) then 1269 if (abs(occ_k(iband))>tol8) then 1270 call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,gh1,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1271 ! Case ipert1=natom+2 (electric field): 1272 ! gh1 contains H^(j1)|u0_k_i> (VHxc constant) which corresponds 1273 ! to i.d/dk in Eq. (38) of Gonze, PRB 55, 10355 (1997). 1274 ! * if ipert==natom+2, we apply directly Eq. (38) 1275 ! * if ipert/=natom+2, Born effective charges are minus D2E 1276 d2nl_k(1,idir1)=d2nl_k(1,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*dotr 1277 d2nl_k(2,idir1)=d2nl_k(2,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*doti 1278 end if 1279 1280 ! Or compute localisation tensor (ddk) 1281 ! See M. Veithen thesis Eq(2.5) 1282 ! MT jan-2010: this is probably not correctly implemented for PAW !!! 1283 ! missing terms due to S^(1) and S^(2) 1284 else 1285 ! note: gh1 used as temporary space (to store idir ddk WF) 1286 if (idir==idir1) then 1287 gh1=cwavef 1288 else 1289 if (need_ddk_file.and.ddkfil(idir1)/=0) then 1290 !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt) 1291 ik_ddk = indkpt1(ikpt) 1292 call wfk_read_bks(ddks(idir1), iband, ik_ddk, isppol, xmpio_single, cg_bks=gh1) 1293 else 1294 gh1=zero 1295 end if 1296 end if 1297 if (abs(occ_k(iband))>tol8) then 1298 call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,gh1,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1299 call dotprod_g(dot1r,dot1i,istwf_k,npw1_k*nspinor,2,cwave0,gh1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1300 call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*nspinor,2,cwavef,cwave0,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1301 dotr=dotr-(dot1r*dot2r-dot1i*dot2i) 1302 doti=doti-(dot1r*dot2i+dot1i*dot2r) 1303 d2nl_k(1,idir1)=d2nl_k(1,idir1)+wtk_k*occ_k(iband)*dotr/(nband_kocc*two) 1304 d2nl_k(2,idir1)=d2nl_k(2,idir1)+wtk_k*occ_k(iband)*doti/(nband_kocc*two) 1305 end if 1306 end if 1307 1308 ! Accumulate here 1st-order density change due to overlap operator changes (if any) 1309 if (has_drho) then 1310 if (abs(occ_k(iband))>tol8) then 1311 ! Compute here delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>] 1312 ! (see PRB 78, 035105 (2008), Eq. (42)) 1313 ABI_ALLOCATE(dcwavef,(2,npw1_k*nspinor)) 1314 ABI_DATATYPE_ALLOCATE(dcwaveprj,(dtset%natom,nspinor)) 1315 call pawcprj_alloc(dcwaveprj,0,gs_hamkq%dimcprj) 1316 call getdc1(cgq,cprjq,dcwavef,dcwaveprj,ibgq,icgq,istwf_k,mcgq,& 1317 & mcprjq,mpi_enreg,dtset%natom,nband_k,npw1_k,nspinor,1,gs1) 1318 1319 ! Accumulate 1st-order density due to delta_u^(j1) 1320 counter=500*iband;option=1;wfcorr=0 1321 call dfpt_accrho(counter,cplex,cwave0,dcwavef,dcwavef,cwaveprj0_idir1,dcwaveprj,& 1322 & lambda,dtfil%filstat,gs_hamkq,iband,idir1,ipert1,isppol,dtset%kptopt,& 1323 & mpi_enreg,dtset%natom,nband_k,1,npw_k,npw1_k,nspinor,occ_k,option,& 1324 & pawdrhoij1_unsym(:,idir1),dtset%prtvol,drhoaug1(:,:,:,idir1),tim_fourwf,wfcorr,wtk_k) 1325 call pawcprj_free(dcwaveprj) 1326 ABI_DEALLOCATE(dcwavef) 1327 ABI_DATATYPE_DEALLOCATE(dcwaveprj) 1328 end if 1329 end if 1330 1331 if((usecprj==1).and..not.(associated(cwaveprj0_idir1,cwaveprj0)))then 1332 call pawcprj_free(cwaveprj0_idir1) 1333 ABI_DATATYPE_DEALLOCATE(cwaveprj0_idir1) 1334 end if 1335 1336 ! End of loops 1337 end do ! idir1 1338 end do ! iband 1339 1340 ! Accumulate contribution of this k-point 1341 d2nl (:,:,ipert1,idir,ipert)=d2nl (:,:,ipert1,idir,ipert)+d2nl_k (:,:) 1342 if (usepaw==1) d2ovl(:,:,ipert1,idir,ipert)=d2ovl(:,:,ipert1,idir,ipert)+d2ovl_k(:,:) 1343 1344 ! Deallocations of arrays used for this k-point 1345 ABI_DEALLOCATE(gh1) 1346 ABI_DEALLOCATE(gs1) 1347 if (need_wfk) then 1348 ABI_DEALLOCATE(cwave0) 1349 end if 1350 if (need_wf1) then 1351 ABI_DEALLOCATE(cwavef) 1352 end if 1353 ABI_DEALLOCATE(kg_k) 1354 ABI_DEALLOCATE(kg1_k) 1355 ABI_DEALLOCATE(ylm_k) 1356 ABI_DEALLOCATE(ylm1_k) 1357 ABI_DEALLOCATE(ylmgr1_k) 1358 ABI_DEALLOCATE(kpg_k) 1359 ABI_DEALLOCATE(kpg1_k) 1360 ABI_DEALLOCATE(d2nl_k) 1361 ABI_DEALLOCATE(d2ovl_k) 1362 ABI_DEALLOCATE(eig_k) 1363 ABI_DEALLOCATE(eig_kq) 1364 ABI_DEALLOCATE(eig1_k) 1365 ABI_DEALLOCATE(occ_k) 1366 if (is_metal) then 1367 ABI_DEALLOCATE(doccde_k) 1368 ABI_DEALLOCATE(doccde_kq) 1369 ABI_DEALLOCATE(occ_kq) 1370 ABI_DEALLOCATE(rocceig) 1371 end if 1372 ABI_DEALLOCATE(dkinpw) 1373 ABI_DEALLOCATE(kinpw1) 1374 ABI_DEALLOCATE(ph3d) 1375 if (allocated(ph3d1)) then 1376 ABI_DEALLOCATE(ph3d1) 1377 end if 1378 ABI_DEALLOCATE(ffnlk) 1379 ABI_DEALLOCATE(ffnl1) 1380 if (ipert1>dtset%natom) then 1381 ABI_DEALLOCATE(ffnl1_idir1) 1382 end if 1383 nullify(ffnl1_idir1) 1384 if (usecprj==1) then 1385 call pawcprj_free(cwaveprj0) 1386 ABI_DATATYPE_DEALLOCATE(cwaveprj0) 1387 end if 1388 nullify(cwaveprj0_idir1) 1389 ! Shift arrays 1390 bdtot_index=bdtot_index+nband_k 1391 bd2tot_index=bd2tot_index+2*nband_k**2 1392 if (mkmem/=0) then 1393 ibg=ibg+nspinor*nband_k 1394 icg=icg+npw_k*nspinor*nband_k 1395 ikg=ikg+npw_k 1396 end if 1397 if (mkqmem/=0) then 1398 ibgq=ibgq+nspinor*nband_k 1399 icgq=icgq+npw1_k*nspinor*nband_k 1400 end if 1401 if (mk1mem/=0) then 1402 ibg1=ibg1+nspinor*nband_k 1403 icg1=icg1+npw1_k*nspinor*nband_k 1404 ikg1=ikg1+npw1_k 1405 end if 1406 1407 end do ! End loop over K-POINTS 1408 1409 ! Transfer 1st-order density change due to overlap; also take into account the spin. 1410 if(has_drho) then 1411 do kdir1=1,mdir1 1412 idir1=jdir1(kdir1) 1413 call fftpac(isppol,mpi_enreg,nspden,cplex*dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),& 1414 & cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),& 1415 & dtset%ngfft,drho1wfr(:,:,idir1),drhoaug1(:,:,:,idir1),1) 1416 end do 1417 end if 1418 1419 end do ! End loop over SPINS 1420 1421 ! Free memory used for this type of perturbation 1422 call destroy_rf_hamiltonian(rf_hamkq) 1423 if (has_drho) then 1424 ABI_DEALLOCATE(drhoaug1) 1425 end if 1426 if (need_pawij10) then 1427 do kdir1=1,mdir1 1428 idir1=jdir1(kdir1) 1429 call paw_ij_free(paw_ij10(:,idir1)) 1430 end do 1431 ABI_DEALLOCATE(e1kbfr_spin) 1432 end if 1433 ABI_DATATYPE_DEALLOCATE(paw_ij10) 1434 1435 ! In case of parallelism, sum 1st-order density and occupation matrix over processors 1436 if (has_drho.and.xmpi_paral==1) then 1437 1438 ! Accumulate 1st-order density 1439 call timab(48,1,tsec) 1440 bufsz=cplex*dtset%nfft*nspden*mdir1 1441 ABI_ALLOCATE(buffer,(bufsz)) 1442 buffer(1:bufsz)=reshape(drho1wfr,(/bufsz/)) 1443 call xmpi_sum(buffer,bufsz,spaceworld,ierr) 1444 drho1wfr(:,:,:)=reshape(buffer(1:bufsz),(/cplex*dtset%nfft,nspden,mdir1/)) 1445 ABI_DEALLOCATE(buffer) 1446 call timab(48,2,tsec) 1447 1448 ! Accumulate 1st-order PAW occupancies 1449 if (usepaw==1) then 1450 call pawrhoij_mpisum_unpacked(pawdrhoij1_unsym,spaceworld) 1451 end if 1452 1453 end if 1454 1455 ! Compute second part of overlap contribution (due to VHxc^(j2)(tild_n+hat_n)) 1456 if (has_drho) then 1457 1458 ABI_ALLOCATE(drhor1,(cplex*nfftf,nspden)) 1459 ABI_ALLOCATE(dnhat1,(cplex*nfftf,nspden)) 1460 1461 ! LOOP OVER PERTURBATION DIRECTIONS 1462 do kdir1=1,mdir1 1463 idir1=jdir1(kdir1) 1464 1465 ! Build and symmetrize 1st-order density change due to change of overlap 1466 ABI_ALLOCATE(drho1wfg,(2,dtset%nfft)) 1467 call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,& 1468 & nspden,nsppol,nsym1,dtset%paral_kgb,phnons1,drho1wfg,drho1wfr(:,:,idir1),& 1469 & rprimd,symaf1,symrl1) 1470 if (dtset%pawstgylm/=0) then 1471 option=0 1472 call pawnhatfr(option,idir1,ipert1,my_natom,dtset%natom,nspden,dtset%ntypat,& 1473 & pawang,pawfgrtab,pawrhoij,pawtab,rprimd,& 1474 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 1475 end if 1476 call pawmkrho(arg,cplex,gprimd,idir1,indsy1,ipert1,& 1477 & mpi_enreg,my_natom,dtset%natom,nspden,nsym1,dtset%ntypat,dtset%paral_kgb,pawang,& 1478 & pawfgr,pawfgrtab,-10001,pawdrhoij1(:,idir1),pawdrhoij1_unsym(:,idir1),pawtab,& 1479 & dtset%qptn,drho1wfg,drho1wfr(:,:,idir1),drhor1,rprimd,symaf1,symrc1,dtset%typat,& 1480 & ucvol,dtset%usewvl,xred,pawang_sym=pawang1,pawnhat=dnhat1,pawrhoij0=pawrhoij) 1481 ABI_DEALLOCATE(drho1wfg) 1482 1483 ! Compute plane-wave contribution to overlap contribution 1484 ! This is subtle as it is a mix of Eq(79) and Eq(80) of PRB 78, 035105 (2008) 1485 ! Details: 1486 ! The VH(tild_nZc)^(1) term of Eq(79) is: 1487 ! <VH(tild_nZc)^(j2)|delta_tild_rho^(j1)> = <vpsp1|drhor1-dnhat1> 1488 ! The first term of Eq(80) is: 1489 ! <VHxc^(j2)|delta_tild_rho^(j1)+delta_hat_rho^(j1)> = <vtrial1-vpsp1|drhor1> 1490 ! The addition of these two terms gives: 1491 ! <vtrial1|drhor1>-<vpsp1|dnhat1> 1492 ! And this is more subtle when usexcnhat=0 1493 call dotprod_vn(cplex,drhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vtrial1,ucvol) 1494 if (usexcnhat/=0) then 1495 call dotprod_vn(cplex,dnhat1,dot2r,dot2i,nfftf,nfftot,1 ,2,vpsp1,ucvol) 1496 else 1497 ABI_ALLOCATE(vtmp1,(cplex*nfftf,nspden)) 1498 do ispden=1,nspden 1499 vtmp1(:,ispden)=vtrial1(:,ispden)-vhartr1(:) 1500 end do 1501 call dotprod_vn(cplex,dnhat1,dot2r,dot2i,nfftf,nfftot,nspden,2,vtmp1,ucvol) 1502 ABI_DEALLOCATE(vtmp1) 1503 end if 1504 dotr=dot1r-dot2r;doti=dot1i-dot2i 1505 1506 ! Compute on-site contributions to overlap contribution 1507 ! (two last terms of Eq(80) of PRB 78, 035105 (2008)) 1508 ! (note: Dij^(j2) and Vxc^(j2) are computed for ipert at first call) 1509 call pawdfptenergy(epawnst,ipert,ipert1,dtset%ixc,my_natom,dtset%natom,dtset%ntypat,& 1510 & nzlmopt_ipert,nzlmopt_ipert1,paw_an,paw_an1,paw_ij1,pawang,dtset%pawprtvol,& 1511 & pawrad,pawrhoij1,pawdrhoij1(:,idir1),pawtab,dtset%pawxcdev,dtset%xclevel,& 1512 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 1513 1514 ! Accumulate in 2nd-order matrix: 1515 ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2 1516 ! has to take the complex conjugate because we want here Int[VHxc^(j1)^*.delta_rho^(j2)] 1517 dotr=dotr+epawnst(1);doti=-(doti+epawnst(2)) 1518 d2ovl_drho(1,idir1,ipert1,idir,ipert)=elfd_fact*dotr 1519 d2ovl_drho(2,idir1,ipert1,idir,ipert)=elfd_fact*doti 1520 1521 end do ! End loop over perturbation directions 1522 1523 ! Free no more needed memory 1524 ABI_DEALLOCATE(drhor1) 1525 ABI_DEALLOCATE(dnhat1) 1526 ABI_DEALLOCATE(drho1wfr) 1527 do iatom=1,my_natom 1528 if (pawfgrtab(iatom)%nhatfr_allocated>0) then 1529 ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr) 1530 end if 1531 pawfgrtab(iatom)%nhatfr_allocated=0 1532 end do 1533 if (paral_atom) then 1534 do kdir1=1,mdir1 1535 idir1=jdir1(kdir1) 1536 call pawrhoij_free(pawdrhoij1_unsym(:,idir1)) 1537 end do 1538 ABI_DATATYPE_DEALLOCATE(pawdrhoij1_unsym) 1539 end if 1540 do kdir1=1,mdir1 1541 idir1=jdir1(kdir1) 1542 call pawrhoij_free(pawdrhoij1(:,idir1)) 1543 end do 1544 ABI_DATATYPE_DEALLOCATE(pawdrhoij1) 1545 end if ! has_drho 1546 1547 ! End loop over perturbations (j1) 1548 end do 1549 1550 !Final deallocations 1551 ABI_DEALLOCATE(ch1c) 1552 if (usepaw==1.and.is_metal_or_qne0) then 1553 ABI_DEALLOCATE(cs1c) 1554 end if 1555 call destroy_hamiltonian(gs_hamkq) 1556 1557 !In case of parallelism, sum over processors 1558 if (xmpi_paral==1)then 1559 call timab(161,1,tsec) 1560 call xmpi_barrier(spaceworld) 1561 call timab(161,2,tsec) 1562 ABI_ALLOCATE(buffer,(6*mpert*(1+usepaw))) 1563 buffer(1:6*mpert)=reshape(d2nl(:,:,:,idir,ipert),(/6*mpert/)) 1564 if (usepaw==1) buffer(6*mpert+1:6*mpert+6*mpert)=reshape(d2ovl(:,:,:,idir,ipert),(/6*mpert/)) 1565 call timab(48,1,tsec) 1566 call xmpi_sum(buffer,6*mpert*(1+usepaw),spaceworld,ierr) 1567 call timab(48,2,tsec) 1568 d2nl (:,:,:,idir,ipert)=reshape(buffer(1:6*mpert),(/2,3,mpert/)) 1569 if (usepaw==1) d2ovl(:,:,:,idir,ipert)=reshape(buffer(6*mpert+1:6*mpert+6*mpert),(/2,3,mpert/)) 1570 ABI_DEALLOCATE(buffer) 1571 end if 1572 1573 !Build complete d2ovl matrix 1574 if (usepaw==1) d2ovl(:,:,:,idir,ipert)=d2ovl(:,:,:,idir,ipert)+d2ovl_drho(:,:,:,idir,ipert) 1575 1576 if (usepaw==1) then 1577 ABI_DEALLOCATE(d2ovl_drho) 1578 end if 1579 1580 !Close the ddk WF files 1581 if (has_ddk_file) then 1582 do kdir1=1,mdir1 1583 idir1=jdir1(kdir1) 1584 if (ddkfil(idir1)/=0)then 1585 call wfk_close(ddks(idir1)) 1586 end if 1587 end do 1588 end if 1589 ABI_DEALLOCATE(jpert1) 1590 ABI_DEALLOCATE(jdir1) 1591 1592 !Symmetrize the phonons contributions, as was needed for the forces in a GS calculation 1593 ABI_ALLOCATE(work,(2,3,dtset%natom)) 1594 do ipert1=1,dtset%natom 1595 do idir1=1,3 1596 work(:,idir1,ipert1)=d2nl(:,idir1,ipert1,idir,ipert) 1597 end do 1598 end do 1599 call dfpt_sygra(dtset%natom,d2nl(:,:,:,idir,ipert),work,indsy1,ipert,nsym1,dtset%qptn,symrc1) 1600 if (usepaw==1) then 1601 do ipert1=1,dtset%natom 1602 do idir1=1,3 1603 work(:,idir1,ipert1)=d2ovl(:,idir1,ipert1,idir,ipert) 1604 end do 1605 end do 1606 call dfpt_sygra(dtset%natom,d2ovl(:,:,:,idir,ipert),work,indsy1,ipert,nsym1,dtset%qptn,symrc1) 1607 end if 1608 ABI_DEALLOCATE(work) 1609 1610 !In the case of the strain perturbation time-reversal symmetry will always 1611 !be true so imaginary part of d2nl will be must be set to zero here since 1612 !the symmetry-reduced kpt set will leave a non-zero imaginary part. 1613 if(ipert==dtset%natom+3.or.ipert==dtset%natom+4) then 1614 d2nl(2,:,:,idir,ipert)=zero 1615 if (usepaw==1) d2ovl(2,:,:,idir,ipert)=zero 1616 else 1617 d2nl(2,:,dtset%natom+3:dtset%natom+4,idir,ipert)=zero 1618 if (usepaw==1) d2ovl(2,:,dtset%natom+3:dtset%natom+4,idir,ipert)=zero 1619 end if 1620 1621 1622 !Symmetrize the strain perturbation contributions, as was needed for the stresses in a GS calculation 1623 !if (ipert==dtset%natom+3.or.ipert==dtset%natom+4)then 1624 if (nsym1>1) then 1625 ABI_ALLOCATE(work,(6,1,1)) 1626 ii=0 1627 do ipert1=dtset%natom+3,dtset%natom+4 1628 do idir1=1,3 1629 ii=ii+1 1630 work(ii,1,1)=d2nl(1,idir1,ipert1,idir,ipert) 1631 end do 1632 end do 1633 call stresssym(gprimd,nsym1,work(:,1,1),symrc1) 1634 ii=0 1635 do ipert1=dtset%natom+3,dtset%natom+4 1636 do idir1=1,3 1637 ii=ii+1 1638 d2nl(1,idir1,ipert1,idir,ipert)=work(ii,1,1) 1639 end do 1640 end do 1641 if (usepaw==1) then 1642 ii=0 1643 do ipert1=dtset%natom+3,dtset%natom+4 1644 do idir1=1,3 1645 ii=ii+1 1646 work(ii,1,1)=d2ovl(1,idir1,ipert1,idir,ipert) 1647 end do 1648 end do 1649 call stresssym(gprimd,nsym1,work(:,1,1),symrc1) 1650 ii=0 1651 do ipert1=dtset%natom+3,dtset%natom+4 1652 do idir1=1,3 1653 ii=ii+1 1654 d2ovl(1,idir1,ipert1,idir,ipert)=work(ii,1,1) 1655 end do 1656 end do 1657 ABI_DEALLOCATE(work) 1658 end if 1659 end if 1660 !end if 1661 1662 !Must also symmetrize the electric field perturbation response ! 1663 !Note: d2ovl is not symetrized because it is zero for electric field perturbation 1664 if (has_ddk_file) then 1665 ABI_ALLOCATE(d2nl_elfd,(2,3)) 1666 ! There should not be any imaginary part, but stay general (for debugging) 1667 d2nl_elfd (:,:)=d2nl(:,:,dtset%natom+2,idir,ipert) 1668 do ii=1,3 1669 sumelfd(:)=zero 1670 do ia=1,nsym1 1671 do jj=1,3 1672 if(symrl1(ii,jj,ia)/=0)then 1673 if(ddkfil(jj)==0)then 1674 blkflg(ii,dtset%natom+2,idir,ipert)=0 1675 end if 1676 end if 1677 end do 1678 symfact(1)=dble(symrl1(ii,1,ia)) 1679 symfact(2)=dble(symrl1(ii,2,ia)) 1680 symfact(3)=dble(symrl1(ii,3,ia)) 1681 sumelfd(:)=sumelfd(:)+symfact(1)*d2nl_elfd(:,1) & 1682 & +symfact(2)*d2nl_elfd(:,2)+symfact(3)*d2nl_elfd(:,3) 1683 end do 1684 d2nl(:,ii,dtset%natom+2,idir,ipert)=sumelfd(:)/dble(nsym1) 1685 end do 1686 ABI_DEALLOCATE(d2nl_elfd) 1687 end if 1688 1689 !Overlap: store the diagonal part of the matrix in the 1690 ! 2nd-order energy non-stationnary expression 1691 eovl1=zero;if (usepaw==1) eovl1=d2ovl(1,idir,ipert,idir,ipert) 1692 1693 call destroy_mpi_enreg(mpi_enreg_seq) 1694 call timab(566,2,tsec) 1695 1696 DBG_EXIT("COLL") 1697 1698 end subroutine dfpt_nstpaw