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-2024 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_mem*mkmem*nsppol)=planewave coefficients of wavefunctions at k cgq(2,mpw1*nspinor*mband_mem*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q. cg1(2,mpw1*nspinor*mband_mem*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_mem*mkmem*nsppol*usecprj)= wave functions at k projected with non-local projectors cprjq(natom,nspinor*mband_mem*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 mband_mem=number of bands per processor 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 vectornd(with_vectornd*nfftf,3)=nuclear dipole moment vector potential 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 with_vectornd = 1 if vectornd allocated 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) [[cite:Audouze2008]]
NOTES
We perform here the computation of
delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
see PRB 78, 035105 (2008), Eq. (42) [[cite:Audouze2008]]
SOURCE
217 subroutine dfpt_nstpaw(blkflg,cg,cgq,cg1,cplex,cprj,cprjq,docckqde,doccde_rbz,dtfil,dtset,d2lo,d2nl,d2ovl,& 218 & eigenq,eigen0,eigen1,eovl1,gmet,gprimd,gsqcut,idir,indkpt1,indsy1,ipert,irrzon1,istwfk_rbz,& 219 & kg,kg1,kpt_rbz,kxc,mgfftf,mkmem,mkqmem,mk1mem,& 220 & mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,mband_mem_rbz,ncpgr,nfftf,ngfftf,nhat,nhat1,& 221 & nkpt_rbz,nkxc,npwarr,npwar1,nspden,nspinor,nsppol,nsym1,n3xccc,occkq,occ_rbz,& 222 & paw_an,paw_an1,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,& 223 & pawrhoij1,pawtab,phnons1,ph1d,ph1df,psps,rhog,rhor,rhor1,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,& 224 & ucvol,usecprj,usepaw,usexcnhat,useylmgr1,vectornd,vhartr1,vpsp1,vtrial,vtrial1,vxc,& 225 & with_vectornd,wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr1) 226 227 !Arguments ------------------------------- 228 !scalars 229 integer,intent(in) :: cplex,idir,ipert,mgfftf,mkmem,mkqmem,mk1mem,mpert,mpw,mpw1 230 integer,intent(in) :: ncpgr,nfftf,nkpt_rbz,nkxc,nspden,nspinor,nsppol,nsym1 231 integer,intent(in) :: n3xccc,usecprj,usepaw,usexcnhat,useylmgr1 232 integer,intent(in) :: mband_mem_rbz,with_vectornd 233 real(dp),intent(in) :: gsqcut,ucvol 234 real(dp),intent(out) :: eovl1 235 type(datafiles_type),intent(in) :: dtfil 236 type(dataset_type),intent(in) :: dtset 237 type(MPI_type),intent(in) :: mpi_enreg 238 type(pawang_type),intent(in) :: pawang,pawang1 239 type(pawfgr_type),intent(in) :: pawfgr 240 type(pseudopotential_type),intent(in) :: psps 241 !arrays 242 integer,intent(in) :: nattyp(dtset%ntypat),nband_rbz(nkpt_rbz*nsppol) 243 integer,intent(in) :: indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom) 244 integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4)) 245 integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 246 integer,intent(in) :: ngfftf(18),npwarr(nkpt_rbz),npwar1(nkpt_rbz) 247 integer,intent(in) :: symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1) 248 integer,intent(inout) :: blkflg(3,mpert,3,mpert) 249 real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem_rbz*mkmem*nsppol) 250 real(dp),intent(in) :: cgq(2,mpw1*nspinor*mband_mem_rbz*mkqmem*nsppol) 251 real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband_mem_rbz*mk1mem*nsppol) 252 real(dp),intent(in) :: docckqde(dtset%mband*nkpt_rbz*nsppol) 253 real(dp),intent(in) :: doccde_rbz(dtset%mband*nkpt_rbz*nsppol) 254 real(dp),intent(in) :: eigenq(dtset%mband*nkpt_rbz*nsppol),eigen0(dtset%mband*nkpt_rbz*nsppol) 255 real(dp),intent(in) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*nsppol) 256 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt_rbz(3,nkpt_rbz) 257 real(dp),intent(in) :: kxc(nfftf,nkxc),nhat(nfftf,nspden),nhat1(cplex*nfftf,nspden*usepaw) 258 real(dp),intent(in) :: occkq(dtset%mband*nkpt_rbz*dtset%nsppol) 259 real(dp),intent(in) :: occ_rbz(dtset%mband*nkpt_rbz*nsppol) 260 real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4)) 261 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 262 real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom),rhog(2,nfftf) 263 real(dp),intent(in) :: rhor(cplex*nfftf,nspden),rhor1(cplex*nfftf,nspden),rmet(3,3),rprimd(3,3) 264 real(dp),intent(in) :: tnons1(3,nsym1),vhartr1(cplex*nfftf),vtrial1(cplex*nfftf,nspden),vxc(nfftf,nspden) 265 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,dtset%natom) 266 real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 267 real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm) 268 real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 269 real(dp),target,intent(in) :: vpsp1(cplex*nfftf),vtrial(nfftf,nspden),xccc3d1(cplex*n3xccc) 270 real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert) 271 real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert),d2ovl(2,3,mpert,3,mpert*usepaw) 272 real(dp),intent(inout) :: vectornd(with_vectornd*nfftf,3) 273 type(pawcprj_type),intent(in) :: cprj(dtset%natom,nspinor*mband_mem_rbz*mkmem*nsppol*usecprj) 274 type(pawcprj_type),intent(in) :: cprjq(dtset%natom,nspinor*mband_mem_rbz*mkqmem*nsppol*usecprj) 275 type(paw_an_type),intent(in) :: paw_an(:) 276 type(paw_an_type),intent(inout) :: paw_an1(:) 277 type(paw_ij_type),intent(in) :: paw_ij(:) 278 type(paw_ij_type),intent(inout) :: paw_ij1(:) 279 type(pawfgrtab_type),intent(inout) :: pawfgrtab(:) 280 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*usepaw) 281 type(pawrhoij_type),intent(in) :: pawrhoij(:) 282 type(pawrhoij_type),intent(in) :: pawrhoij1(:) 283 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw) 284 285 !Local variables------------------------------- 286 !scalars 287 integer,parameter :: tim_fourwf=18,tim_getgh1c=2,tim_projbd=3,formeig1=1 288 integer :: bd2tot_index,bdtot_index,berryopt,bufsz,choice,cpopt,cplex_rhoij,ddkcase 289 integer :: dimffnl,dimffnl1,dimffnl1_idir1,dimylmgr1,g0term 290 integer :: ia,iatom,iband,ibg,ibgq,ibg1,icg,icgq,icg1,ider,idir0,idir1,idir_cprj 291 integer :: ierr,ii,ikg,ikg1,ikpt,ikpt_me,ilmn,iorder_cprj,ipert1 292 integer :: ispden,isppol,istwf_k,istr,istr1,itypat,jband,jj,kdir1,kpert1,master,mcgq,mcprjq 293 integer :: mdir1,me,mpert1,my_natom,my_comm_atom,my_nsppol,nband_k,nband_kocc,need_ylmgr1 294 integer :: nfftot,nkpg,nkpg1,nkpt_me,npw_,npw_k,npw1_k,nspden_rhoij 295 integer :: nvh1,nvxc1,nzlmopt_ipert,nzlmopt_ipert1,optlocal,optnl 296 integer :: option,opt_gvnlx1,qphase_rhoij,sij_opt,spaceworld,usevnl,wfcorr,ik_ddk 297 integer :: nband_me, iband_me, jband_me, iband_ 298 integer :: do_scprod, do_bcast 299 integer :: startband, endband 300 real(dp) :: arg,doti,dotr,dot1i,dot1r,dot2i,dot2r,dot3i,dot3r,elfd_fact,invocc,lambda,wtk_k 301 logical :: force_recompute,has_dcwf,has_dcwf2,has_drho,has_ddk_file,has_vectornd 302 logical :: is_metal,is_metal_or_qne0,need_ddk_file,need_pawij10 303 logical :: need_wfk,need_wf1,nmxc,paral_atom,qne0,t_exist 304 character(len=500) :: msg 305 character(len=fnlen) :: fiwfddk(3) 306 type(gs_hamiltonian_type) :: gs_hamkq 307 type(rf_hamiltonian_type) :: rf_hamkq 308 type(MPI_type) :: mpi_enreg_seq 309 !arrays 310 integer :: ddkfil(3),my_spintab(2),nband_tmp(1),npwar1_tmp(1) 311 integer,allocatable :: bands_treated_now(:),band_procs(:) 312 integer,allocatable :: jpert1(:),jdir1(:),kg1_k(:,:),kg_k(:,:) 313 integer,pointer :: my_atmtab(:) 314 real(dp) :: dum1(1,1),dum2(1,1),dum3(1,1),epawnst(2),kpoint(3),kpq(3) 315 real(dp) :: sumelfd(2),symfact(3),tsec(2),ylmgr_dum(1,3,1) 316 real(dp),allocatable :: buffer(:),ch1c(:,:,:,:),cs1c(:,:,:,:) 317 real(dp),allocatable :: ch1c_tmp(:,:) 318 real(dp),allocatable :: cs1c_tmp(:,:) 319 real(dp),allocatable :: cwave0(:,:),cwavef(:,:),dcwavef(:,:) 320 real(dp),allocatable :: cg_ddk(:,:,:) 321 real(dp),allocatable :: doccde_k(:),doccde_kq(:) 322 real(dp),allocatable :: dnhat1(:,:),drhoaug1(:,:,:,:) 323 real(dp),allocatable :: drhor1(:,:),drho1wfg(:,:),drho1wfr(:,:,:) 324 real(dp),allocatable :: d2nl_elfd(:,:),dkinpw(:) 325 real(dp),allocatable :: d2nl_k(:,:),d2ovl_drho(:,:,:,:,:),d2ovl_k(:,:) 326 real(dp),allocatable :: eig_k(:),eig_kq(:),eig1_k(:) 327 real(dp),allocatable,target :: e1kbfr_spin(:,:,:,:,:,:),ffnlk(:,:,:,:),ffnl1(:,:,:,:) 328 real(dp),allocatable :: gh1(:,:),gs1(:,:),gvnlx1(:,:),kinpw1(:),kpg_k(:,:),kpg1_k(:,:) 329 real(dp),allocatable :: gvnlx1_tmp(:,:) 330 real(dp),allocatable :: occ_k(:),occ_kq(:),ph3d(:,:,:),ph3d1(:,:,:),rhotmp(:,:),rocceig(:,:) 331 real(dp),allocatable :: vectornd_pac(:,:,:,:,:),vectornd_pac_idir(:,:,:,:) 332 real(dp),allocatable :: ylm_k(:,:),ylm1_k(:,:),ylmgr1_k(:,:,:),vtmp1(:,:),vxc10(:,:) 333 real(dp),allocatable,target :: work(:,:,:),e1kb_work(:,:,:,:) 334 real(dp),pointer :: e1kbfr(:,:,:,:,:),e1kb_ptr(:,:,:,:) 335 real(dp),pointer :: ffnl1_idir1(:,:,:,:),vhartr01(:),vpsp1_idir1(:),xccc3d1_idir1(:) 336 type(pawcprj_type),allocatable :: dcwaveprj(:,:) 337 type(pawcprj_type),allocatable,target :: cwaveprj0(:,:) 338 type(pawcprj_type),pointer :: cwaveprj0_idir1(:,:) 339 type(paw_ij_type),allocatable :: paw_ij10(:,:) 340 type(pawrhoij_type),target,allocatable :: pawdrhoij1(:,:) 341 type(pawrhoij_type),pointer :: pawdrhoij1_unsym(:,:) 342 type(wfk_t) :: ddks(3) 343 344 ! ********************************************************************* 345 346 DBG_ENTER("COLL") 347 348 !Keep track of total time spent in dfpt_nstpaw 349 call timab(566,1,tsec) 350 351 !Not valid for PrintBandByBand 352 if (dtset%prtbbb/=0) then 353 ABI_BUG('not yet valid for prtbbb/=0!') 354 end if 355 356 !NCPP restrictions 357 if (usepaw==0) then 358 ! cprj cannot be used 359 if (usecprj/=0) then 360 ABI_BUG('NCPP: usecprj should be 0!') 361 end if 362 ! d2ovl cannot be used 363 if (size(d2ovl)/=0) then 364 ABI_BUG('NCPP: d2ovl should not be allocated!') 365 end if 366 end if 367 368 !PAW restrictions 369 if (usepaw==1) then 370 ! Test on FFT grid sizes 371 if (pawfgr%nfft/=nfftf) then 372 ABI_BUG('PAW: wrong values for nfft, nfftf!') 373 end if 374 ! Test gradients of cprj 375 if (ipert<=dtset%natom.and.ncpgr/=3) then 376 ABI_BUG('PAW: wrong value of ncpgr for ipert<=natom!') 377 end if 378 if (ipert==dtset%natom+1.and.ncpgr/=1.and.dtset%orbmag==0) then 379 ABI_BUG('PAW: wrong value of ncpgr for ipert=natom+1!') 380 end if 381 if (ipert==dtset%natom+2.and.ncpgr/=3) then 382 ABI_BUG('PAW: wrong value of ncpgr for ipert=natom+2!') 383 end if 384 if ((ipert==dtset%natom+3.or.ipert==dtset%natom+4).and.ncpgr/=1) then 385 ABI_BUG('PAW: wrong value of ncpgr for ipert=natom+3 or 4!') 386 end if 387 ! Test on availability of DijHartree and XC on-site potentials 388 if (mpi_enreg%my_natom>0.and.ipert/=dtset%natom+1) then 389 if (paw_ij1(1)%has_dijhartree==0.or.paw_an1(1)%has_vxc==0) then 390 msg='PAW: paw_ij1%dijhartree and paw=_an1%vxc1 should be allocated !' 391 end if 392 end if 393 end if 394 395 !Set up parallelism 396 master=0;me=mpi_enreg%me_kpt 397 spaceworld=mpi_enreg%comm_cell 398 paral_atom=(mpi_enreg%my_natom/=dtset%natom) 399 my_comm_atom=mpi_enreg%comm_atom 400 my_natom=mpi_enreg%my_natom 401 my_atmtab=>mpi_enreg%my_atmtab 402 my_spintab=mpi_enreg%my_isppoltab 403 my_nsppol=count(my_spintab==1) 404 405 !Fake MPI data to be used in sequential calls to parallel routines 406 call initmpi_seq(mpi_enreg_seq) 407 mpi_enreg_seq%my_natom=dtset%natom 408 409 !Compute effective number of k-points 410 nkpt_me=nkpt_rbz 411 if(xmpi_paral==1)then 412 nkpt_me=0 413 do isppol=1,nsppol 414 do ikpt=1,nkpt_rbz 415 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 416 if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))) nkpt_me=nkpt_me+1 417 end do 418 end do 419 end if 420 421 !Sizes for WF at k+q 422 mcgq=mpw1*nspinor*mband_mem_rbz*mkqmem*nsppol 423 mcprjq=nspinor*mband_mem_rbz*mkqmem*nsppol*usecprj 424 425 ABI_MALLOC(bands_treated_now, (maxval(nband_rbz))) 426 427 !Check ddk files (needed to compute electric field perturbations) 428 ddkfil(:)=0 429 do idir1=1,3 430 ddkcase=idir1+dtset%natom*3 431 call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk(idir1)) 432 t_exist = file_exists(fiwfddk(idir1)) 433 434 if (.not. t_exist) then 435 ! Try netcdf file. 436 t_exist = file_exists(nctk_ncify(fiwfddk(idir1))) 437 if (t_exist) then 438 fiwfddk(idir1) = nctk_ncify(fiwfddk(idir1)) 439 write(msg,"(3a)")"- File: ",trim(fiwfddk(idir1))," does not exist but found netcdf file with similar name." 440 call wrtout(std_out,msg,'COLL') 441 end if 442 end if 443 444 if (t_exist) then 445 ddkfil(idir1)=20+idir1 ! Note the use of unit numbers 21, 22 and 23 446 end if 447 end do 448 has_ddk_file=(any(ddkfil(:)>0)) 449 450 !Define the set of perturbations (j1)=(ipert1,idir1) 451 !The first perturbation must be (j1)=(j2)=(ipert,idir) 452 !because we need to compute <g|H^(j2)-Eps.S^(j2)|u0> first. 453 if (ipert/=dtset%natom+1) then 454 mpert1=0 455 ABI_MALLOC(jpert1,(mpert)) 456 jpert1 = 0 457 if (ipert/=dtset%natom+2.or.has_ddk_file) then 458 mpert1=mpert1+1;jpert1(mpert1)=ipert 459 end if 460 do ipert1=1,mpert 461 if (ipert1/=ipert.and.& 462 & (ipert1<=dtset%natom.or.& 463 & (ipert1==dtset%natom+2.and.has_ddk_file).or.& 464 & ((ipert>dtset%natom.and.ipert/=dtset%natom+5).and.(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)).or. & 465 & ((ipert1==dtset%natom+2).and.has_ddk_file))) then 466 mpert1=mpert1+1;jpert1(mpert1)=ipert1 467 end if 468 end do 469 else 470 mpert1=1 471 ABI_MALLOC(jpert1,(mpert1)) 472 jpert1(1)=dtset%natom+1 473 end if 474 mdir1=3 475 ABI_MALLOC(jdir1,(mdir1)) 476 jdir1(1:3)= (/ (idir1,idir1=1,3) /) 477 jdir1(1)=idir;jdir1(idir)=1 478 479 !Index of strain perturbation, if any 480 istr=idir;if (ipert==dtset%natom+4) istr=idir+3 481 482 !Open ddk WF file(s) in sequential mode 483 if (has_ddk_file) then 484 do kdir1=1,mdir1 485 idir1=jdir1(kdir1) 486 if (ddkfil(idir1)/=0) then 487 write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk(idir1)) 488 call wrtout([std_out, ab_out],msg) 489 call wfk_open_read(ddks(idir1),fiwfddk(idir1),formeig1,dtset%iomode,ddkfil(idir1), xmpi_comm_self) 490 end if 491 end do 492 493 ABI_MALLOC(cg_ddk,(2,mpw1*nspinor*mband_mem_rbz,3)) 494 cg_ddk = zero ! not all may be initialized below if only certain ddk directions are provided 495 end if 496 497 !Zero only portion of matrix to be computed here 498 d2nl(:,:,1:dtset%natom+4,idir,ipert)=zero 499 if (usepaw==1) d2ovl(:,:,1:dtset%natom+4,idir,ipert)=zero 500 501 !Update list of computed matrix elements 502 do kpert1=1,mpert1 503 ipert1=jpert1(kpert1) 504 do kdir1=1,mdir1 505 idir1=jdir1(kdir1) 506 if ((ipert1<=dtset%natom).or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4.or.& 507 & (ipert1==dtset%natom+1.and.((ddkfil(idir1)/=0).or.(dtset%rfdir(idir1)/=0.and.idir1<=idir))).or.& 508 & (ipert1==dtset%natom+2.and.ddkfil(idir1)/=0).or.& 509 & ((ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and.ddkfil(idir1)/=0)) then 510 blkflg(idir1,ipert1,idir,ipert)=1 511 end if 512 end do 513 end do 514 515 !Initialize most of the (1st-order) Hamiltonian 516 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 517 !2) Perform the setup needed for the non-local factors: 518 call init_hamiltonian(gs_hamkq,psps,pawtab,nspinor,nsppol,nspden,dtset%natom,& 519 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,& 520 & paw_ij=paw_ij,mpi_atmtab=my_atmtab,comm_atom=my_comm_atom,mpi_spintab=mpi_enreg%my_isppoltab,& 521 & usecprj=usecprj,nucdipmom=dtset%nucdipmom,gpu_option=dtset%gpu_option) 522 has_vectornd = (with_vectornd .EQ. 1) 523 if(has_vectornd) then 524 ABI_MALLOC(vectornd_pac,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),gs_hamkq%nvloc,3)) 525 ABI_MALLOC(vectornd_pac_idir,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),gs_hamkq%nvloc)) 526 end if 527 528 !Variables common to all perturbations 529 arg=maxval(occ_rbz)-minval(occ_rbz) 530 qne0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2>=tol14) 531 is_metal=((dtset%occopt>=3.and.dtset%occopt<=8).or.(abs(arg)>tol8)) 532 is_metal_or_qne0=((is_metal).or.(qne0)) 533 nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 534 ABI_MALLOC(ch1c,(2,mband_mem_rbz,dtset%mband,nkpt_me)) 535 ABI_MALLOC(ch1c_tmp,(2,mband_mem_rbz)) 536 ch1c(:,:,:,:)=zero 537 nzlmopt_ipert=0;nzlmopt_ipert1=0 538 if (usepaw==1) then 539 ABI_MALLOC(d2ovl_drho,(2,3,mpert,3,mpert)) 540 d2ovl_drho=zero 541 if (dtset%pawnzlm/=0) then 542 nzlmopt_ipert=1;if (dtset%nstep<2) nzlmopt_ipert=-1 543 nzlmopt_ipert1=-1 544 end if 545 if (is_metal_or_qne0) then 546 ABI_MALLOC(cs1c,(2,dtset%mband,mband_mem_rbz,nkpt_me)) 547 ABI_MALLOC(cs1c_tmp,(2,dtset%mband)) 548 cs1c(:,:,:,:)=zero 549 end if 550 end if 551 552 !Force the recomputation of on-site potentials and DijHartree 553 if (usepaw==1) then 554 call paw_an_reset_flags(paw_an1) 555 call paw_ij_reset_flags(paw_ij1,dijhartree=.true.) 556 end if 557 558 559 !LOOP OVER PERTURBATION TYPES (j1) 560 do kpert1=1,mpert1 561 ipert1=jpert1(kpert1) 562 563 ! Flag for use of DDK file 564 need_ddk_file=(has_ddk_file.and.(ipert1==dtset%natom+1.or.ipert1==dtset%natom+2)) 565 566 ! Factor to be applied for electric Field (Eff. charges and piezo. tensor are "minus" d2E) 567 elfd_fact=one 568 if ((ipert <=dtset%natom.or.ipert ==dtset%natom+3.or.ipert ==dtset%natom+4).and. & 569 & (ipert1==dtset%natom+2)) elfd_fact=-one 570 if ((ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and. & 571 & (ipert ==dtset%natom+2)) elfd_fact=-one 572 573 ! We want to compute delta_u^(j1))=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>] 574 ! see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42) 575 has_dcwf=.false.;has_dcwf2=.false.;has_drho=.false. 576 if (usepaw==1) then 577 has_dcwf =(ipert1/=dtset%natom+2) 578 has_dcwf2=(ipert /=dtset%natom+2) 579 has_drho =(has_dcwf.and.ipert1/=dtset%natom+1) 580 end if 581 582 ! Select which WF are needed 583 need_wfk=(.true.) 584 need_wf1=(.true.) 585 586 ! Initialize data for NL 1st-order (j1) hamiltonian 587 call init_rf_hamiltonian(cplex,gs_hamkq,ipert1,rf_hamkq,mpi_spintab=[0,0]) 588 589 ! The following contributions are needed only for non-DDK perturbation: 590 ! - Frozen part of 1st-order Dij 591 ! - Contribution from local potential to dynamical matrix (due to Vxc^(j1)(tild_nc)+VH^(j1)(tild_nZc)) 592 if (ipert/=dtset%natom+1.and.ipert1/=dtset%natom+1) then 593 594 ! Allocations 595 force_recompute=(usepaw==0) ! This is dangerous... 596 nfftot=ngfftf(1)*ngfftf(2)*ngfftf(3) 597 nvxc1=0;if (ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) nvxc1=nspden 598 nvh1=0;if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) nvh1=1 599 ABI_MALLOC(vxc10,(cplex*nfftf,nvxc1)) 600 ABI_MALLOC(vhartr01,(cplex*nfftf*nvh1)) 601 need_pawij10=(usepaw==1) 602 if (need_pawij10) then 603 ABI_MALLOC(paw_ij10,(my_natom,mdir1)) 604 ABI_MALLOC(e1kbfr_spin,(rf_hamkq%dime1kb1,rf_hamkq%dime1kb2,nspinor**2,cplex,mdir1,my_nsppol)) 605 else 606 ABI_MALLOC(paw_ij10,(0,0)) 607 end if 608 609 ! LOOP OVER PERTURBATION DIRECTIONS 610 do kdir1=1,mdir1 611 idir1=jdir1(kdir1) 612 istr1=idir1;if (ipert1==dtset%natom+4) istr1=idir1+3 613 614 ! Get first-order local potential and first-order pseudo core density 615 if (ipert==ipert1.and.idir==idir1.and.(.not.force_recompute)) then 616 vpsp1_idir1 => vpsp1 617 xccc3d1_idir1 => xccc3d1 618 else 619 ABI_MALLOC(vpsp1_idir1,(cplex*nfftf)) 620 ABI_MALLOC(xccc3d1_idir1,(cplex*n3xccc)) 621 if (usepaw==1) then 622 call dfpt_atm2fft(gs_hamkq%atindx,cplex,gmet,gprimd,gsqcut,istr1,ipert1,& 623 & mgfftf,psps%mqgrid_vl,dtset%natom,1,nfftf,ngfftf,dtset%ntypat,& 624 & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 625 & atmrhor1=xccc3d1_idir1,atmvlocr1=vpsp1_idir1,optv_in=1,optn_in=n3xccc/nfftf,optn2_in=1,& 626 & vspl=psps%vlspl) 627 else 628 if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 629 630 !To compute Absolute Deformation Potentials toghether with FxE tensor 631 !the reference has to be the same as in the FxE routines 632 g0term=0; if (dtset%rfstrs_ref==1) g0term=1 633 634 call vlocalstr(gmet,gprimd,gsqcut,istr1,mgfftf,mpi_enreg,psps%mqgrid_vl,dtset%natom,& 635 & nattyp,nfftf,ngfftf,dtset%ntypat,ph1df,psps%qgrid_vl,ucvol,& 636 & psps%vlspl,vpsp1_idir1,g0term=g0term) 637 else 638 call dfpt_vlocal(gs_hamkq%atindx,cplex,gmet,gsqcut,idir1,ipert1,mpi_enreg,psps%mqgrid_vl,& 639 & dtset%natom,nattyp,nfftf,ngfftf,dtset%ntypat,ngfftf(1),ngfftf(2),ngfftf(3),& 640 & ph1df,psps%qgrid_vl,dtset%qptn,ucvol,psps%vlspl,vpsp1_idir1,xred) 641 end if 642 if(psps%n1xccc/=0)then 643 call dfpt_mkcore(cplex,idir1,ipert1,dtset%natom,dtset%ntypat,ngfftf(1),psps%n1xccc,& 644 & ngfftf(2),ngfftf(3),dtset%qptn,rprimd,dtset%typat,ucvol,& 645 & psps%xcccrc,psps%xccc1d,xccc3d1_idir1,xred) 646 end if 647 end if 648 end if 649 650 ! Compute 1st-order non-local factors (Dij^(j1)_fr) 651 if (need_pawij10) then 652 call paw_ij_nullify(paw_ij10(:,idir1)) 653 call paw_ij_init(paw_ij10(:,idir1),cplex,nspinor,dtset%nsppol,dtset%nspden,& 654 & 0,dtset%natom,dtset%ntypat,dtset%typat,pawtab,has_dijfr=1,& 655 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 656 if (ipert/=ipert1.or.idir/=idir1.or.force_recompute) then 657 option=0 658 call pawdijfr(gprimd,idir1,ipert1,my_natom,dtset%natom,nfftf,ngfftf,& 659 & nspden,nsppol,dtset%ntypat,option,paw_ij10(:,idir1),pawang,pawfgrtab,pawrad,& 660 & pawtab,cplex,dtset%qptn,rprimd,ucvol,vpsp1_idir1,vtrial,vxc,xred,& 661 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 662 else 663 do iatom=1,my_natom 664 paw_ij10(iatom,idir1)%has_dijfr=paw_ij1(iatom)%has_dijfr 665 if (paw_ij1(iatom)%has_dijfr==2) paw_ij10(iatom,idir1)%dijfr=paw_ij1(iatom)%dijfr 666 end do 667 end if 668 end if 669 670 ! Get first-order exchange-correlation potential (core-correction contribution only) 671 if (nvxc1>0) then 672 if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 673 option=0 674 call dfpt_mkvxcstr(cplex,idir1,ipert1,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,& 675 & nhat,nhat1,nkxc,nmxc,nspden,n3xccc,option,dtset%qptn,rhor,rhor1,rprimd,& 676 & usepaw,usexcnhat,vxc10,xccc3d1_idir1) 677 else 678 ! Non-collinear magnetism (should the second nkxc be nkxc_cur ?) 679 if (nspden==4) then 680 option=0 681 call dfpt_mkvxc_noncoll(cplex,dtset%ixc,kxc,mpi_enreg,nfftf,ngfftf,dum1,0,dum2,0,dum3,0,nkxc,& 682 & nmxc,nspden,n3xccc,1,option,dtset%qptn,dum1,dum1,rprimd,0,vxc,& 683 & vxc10,xccc3d1_idir1) 684 else 685 call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfftf,ngfftf,dum2,0,dum3,0,nkxc,& 686 & nmxc,nspden,n3xccc,0,dtset%qptn,dum1,rprimd,0,vxc10,xccc3d1_idir1) 687 end if 688 end if 689 end if 690 691 ! Get first-order Hartree potential (metric tensor contribution only) 692 if (nvh1>0) then 693 call hartrestr(gsqcut,idir1,ipert1,mpi_enreg,dtset%natom,& 694 & nfftf,ngfftf,rhog,rprimd,vhartr01) 695 end if 696 697 ! Get Hartree + xc + local contributions to dynamical matrix or elastic tensor 698 if (ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 699 if (usepaw==0) then 700 ! vxc1 is integrated with the total 1st-order density (rhor1 ) 701 ! vpsp1 is integrated with the 1st-order pseudo density (rhor1) 702 call dotprod_vn(cplex,rhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol) 703 call dotprod_vn(cplex,rhor1,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol) 704 else if (usexcnhat/=0) then 705 ! vxc1 is integrated with the total 1st-order density (rhor1 including nhat1) 706 ! vpsp1 is integrated with the 1st-order pseudo density (rhor1 without nhat1) 707 ABI_MALLOC(rhotmp,(cplex*nfftf,1)) 708 rhotmp(:,1)=rhor1(:,1)-nhat1(:,1) 709 call dotprod_vn(cplex,rhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol) 710 call dotprod_vn(cplex,rhotmp,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol) 711 ABI_FREE(rhotmp) 712 else 713 ! vxc1 is integrated with the 1st-order pseudo density (rhor1 without nhat1) 714 ! vpsp1 is integrated with the 1st-order pseudo density (rhor1 without nhat1) 715 ABI_MALLOC(rhotmp,(cplex*nfftf,nspden)) 716 rhotmp(:,:)=rhor1(:,:)-nhat1(:,:) 717 call dotprod_vn(cplex,rhotmp,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol) 718 call dotprod_vn(cplex,rhotmp,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol) 719 ABI_FREE(rhotmp) 720 end if 721 if (nvh1>0) then 722 call dotprod_vn(cplex,rhor1,dot3r,dot3i,nfftf,nfftot,1,2,vhartr01,ucvol) 723 else 724 dot3r=zero ; dot3i=zero 725 end if 726 ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2 727 dotr=dot1r+dot2r+dot3r;doti=dot1i+dot2i+dot3i 728 ! In case ipert = natom+2, these lines compute the local part 729 ! of the Born effective charges from phonon and electric 730 ! field type perturbations, see eq. 43 of X. Gonze and C. Lee, PRB 55, 10355 (1997) [[cite:Gonze1997a]] 731 ! The minus sign is due to the fact that the effective charges 732 ! are minus the second derivatives of the energy 733 if (ipert/=dtset%natom+1) then 734 d2lo(1,idir1,ipert1,idir,ipert)=elfd_fact*dotr 735 d2lo(2,idir1,ipert1,idir,ipert)=elfd_fact*doti 736 end if 737 end if ! ipert1<=natom 738 739 if (ipert/=ipert1.or.idir/=idir1.or.force_recompute) then 740 ABI_FREE(vpsp1_idir1) 741 ABI_FREE(xccc3d1_idir1) 742 end if 743 744 ! End loop on directions 745 end do 746 747 ! Free memory 748 ABI_FREE(vxc10) 749 ABI_FREE(vhartr01) 750 751 else ! ddk perturbation 752 d2lo(1:2,1:mdir1,ipert1,idir,ipert)=zero 753 need_pawij10=.false. 754 ABI_MALLOC(paw_ij10,(0,0)) 755 end if 756 757 ! Prepare RF PAW files for reading and writing if mkmem, mkqmem or mk1mem==0 758 iorder_cprj=0 759 760 ! Allocate arrays used to accumulate density change due to overlap 761 if (has_drho) then 762 ABI_MALLOC(drhoaug1,(cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),mdir1)) 763 ABI_MALLOC(drho1wfr,(cplex*dtset%nfft,dtset%nspden,mdir1)) 764 ABI_MALLOC(pawdrhoij1,(my_natom,mdir1)) 765 if (paral_atom) then 766 ABI_MALLOC(pawdrhoij1_unsym,(dtset%natom,mdir1)) 767 else 768 pawdrhoij1_unsym => pawdrhoij1 769 end if 770 do kdir1=1,mdir1 771 idir1=jdir1(kdir1) 772 drho1wfr(:,:,idir1)=zero 773 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 774 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc) 775 call pawrhoij_alloc(pawdrhoij1(:,idir1),cplex_rhoij,nspden_rhoij,dtset%nspinor,& 776 & dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,use_rhoijp=1,use_rhoij_=0,& 777 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=my_atmtab) 778 if (paral_atom) then 779 call pawrhoij_alloc(pawdrhoij1_unsym(:,idir1),cplex_rhoij,nspden_rhoij,dtset%nspinor,& 780 & dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,& 781 & use_rhoijp=0,use_rhoij_=1) 782 else 783 call pawrhoij_init_unpacked(pawdrhoij1_unsym(:,idir1)) 784 end if 785 end do 786 end if 787 788 ! Initialize shifts for global arrays 789 bdtot_index=0 790 bd2tot_index=0 791 ibg=0;icg=0 792 ibg1=0;icg1=0 793 ibgq=0;icgq=0 794 795 796 ! Has to get 1st-order non-local factors before the loop over spins 797 ! because this needs a communication over comm_atom (=comm_spinkpt) 798 if (need_pawij10) then 799 if (my_nsppol<nsppol) then 800 ABI_MALLOC(e1kb_work,(rf_hamkq%dime1kb1,rf_hamkq%dime1kb2,nspinor**2,cplex)) 801 end if 802 ii=0 803 do isppol=1,nsppol 804 if (my_spintab(isppol)==1) ii=ii+1 805 if (my_spintab(isppol)/=1) e1kb_ptr => e1kb_work 806 do kdir1=1,mdir1 807 idir1=jdir1(kdir1) 808 if (my_spintab(isppol)==1) e1kb_ptr => e1kbfr_spin(:,:,:,:,idir1,ii) 809 call pawdij2e1kb(paw_ij10(:,idir1),isppol,my_comm_atom,e1kbfr=e1kb_ptr,mpi_atmtab=my_atmtab) 810 end do 811 end do 812 if (my_nsppol<nsppol) then 813 ABI_FREE(e1kb_work) 814 end if 815 end if 816 817 ! LOOP OVER SPINS 818 do isppol=1,nsppol 819 820 ikpt_me=0 821 822 ! Rewind (k+G) data if needed 823 ikg=0;ikg1=0 824 825 ! Continue to initialize the GS/RF Hamiltonian 826 call gs_hamkq%load_spin(isppol,with_nonlocal=.true.) 827 if (need_pawij10) then 828 ii=min(isppol,size(e1kbfr_spin,6)) 829 if (ii>0) e1kbfr => e1kbfr_spin(:,:,:,:,:,ii) 830 end if 831 832 ! if vectornd is present, set it up for addition to gs_hamkq and rf_hamkq. 833 ! Note that it must be done for the three Cartesian directions. Also, the following 834 ! code assumes explicitly and implicitly that nvloc = 1. This should eventually be generalized. 835 if(has_vectornd) then 836 call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, & 837 & dtset%nspden, gs_hamkq%nvloc, 3, pawfgr, mpi_enreg, vectornd,vectornd_pac) 838 call gs_hamkq%load_spin(isppol, vectornd=vectornd_pac) 839 vectornd_pac_idir(:,:,:,:)=vectornd_pac(:,:,:,:,idir) 840 call rf_hamkq%load_spin(isppol, vectornd=vectornd_pac_idir) 841 end if 842 843 ! Initialize accumulation of density 844 if (has_drho) drhoaug1(:,:,:,:)=zero 845 846 ! LOOP OVER K-POINTS 847 do ikpt=1,nkpt_rbz 848 849 ! Load dimensions for this k-point 850 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 851 istwf_k=istwfk_rbz(ikpt) 852 npw_k=npwarr(ikpt) 853 npw1_k=npwar1(ikpt) 854 nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me) 855 856 ! Skip loop if this k-point is not to be treated by this proc 857 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 858 bdtot_index=bdtot_index+nband_k 859 bd2tot_index=bd2tot_index+2*nband_k**2 860 cycle ! Skip the rest of the k-point loop 861 end if 862 863 ! Allocate/initialize local arrays and scalars for this k-point 864 ABI_MALLOC(d2nl_k,(2,3)) 865 ABI_MALLOC(d2ovl_k,(2,3)) 866 ABI_MALLOC(eig_k,(nband_k)) 867 ABI_MALLOC(eig_kq,(nband_k)) 868 ABI_MALLOC(eig1_k,(2*nband_k**2)) 869 ABI_MALLOC(occ_k,(nband_k)) 870 d2nl_k(:,:)=zero 871 d2ovl_k(:,:)=zero 872 eig_k (:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 873 eig_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index) 874 eig1_k(:)=eigen1(1+bd2tot_index:2*nband_k**2+bd2tot_index) 875 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 876 nband_kocc=count(abs(occ_k(:))>tol8) 877 kpoint(:)=kpt_rbz(:,ikpt) 878 kpq(:)=kpoint(:);if (ipert1<dtset%natom+3) kpq(:)=kpq(:)+dtset%qptn(1:3) 879 wtk_k=wtk_rbz(ikpt) 880 need_ylmgr1=0;dimylmgr1=0 881 nkpg=0;nkpg1=0 882 ikpt_me=ikpt_me+1 883 if (is_metal) then 884 ! For each pair of active bands (m,n), generates the ratios 885 ! rocceig(m,n)=(occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)) 886 ABI_MALLOC(doccde_k,(nband_k)) 887 ABI_MALLOC(doccde_kq,(nband_k)) 888 ABI_MALLOC(occ_kq,(nband_k)) 889 ABI_MALLOC(rocceig,(nband_k,nband_k)) 890 doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index) 891 doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index) 892 occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index) 893 call occeig(doccde_k,doccde_kq,eig_k,eig_kq,nband_k,dtset%occopt,occ_k,occ_kq,rocceig) 894 end if 895 896 ! Take care of the npw and kg records in WF and DDK files 897 if (need_ddk_file) then 898 do kdir1=1,mdir1 899 idir1=jdir1(kdir1) 900 if (ddkfil(idir1)/=0)then 901 !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt) 902 ik_ddk = indkpt1(ikpt) 903 npw_ = ddks(idir1)%hdr%npwarr(ik_ddk) 904 if (npw_/=npw_k) then 905 write(msg, '(a,i0,a,i0,a,i0,a,a,i0,a,a,i0)')& 906 'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',idir,ch10,& 907 'the number of plane waves in the ddk file is equal to', npw_,ch10,& 908 'while it should be ',npw_k 909 ABI_ERROR(msg) 910 end if 911 912 ! NB: this will fail if the bands are not contiguous. 913 startband = nband_k 914 endband = 1 915 do iband=1,nband_k 916 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) == mpi_enreg%me_kpt) then 917 if (iband < startband) startband = iband 918 if (iband > endband) endband = iband 919 end if 920 end do 921 ! NB: eig_k is band distributed in call to read_band_block, though array has full size, 922 ! only certain columns for my iband are filled, then used below 923 call ddks(idir1)%read_band_block((/startband,endband/),ik_ddk,isppol,xmpio_collective, & 924 & cg_k=cg_ddk(:,:,idir1)) 925 end if ! ddk file is already present 926 end do 927 end if 928 929 ! Allocate arrays used for NL form factors 930 ABI_MALLOC(kg_k,(3,npw_k)) 931 ABI_MALLOC(kg1_k,(3,npw1_k)) 932 ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 933 ABI_MALLOC(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)) 934 if (psps%useylm==1.and.(need_ddk_file.or.ipert1==dtset%natom+1)) need_ylmgr1=1 935 if (psps%useylm==1.and.(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)) need_ylmgr1=1 936 dimylmgr1=max(useylmgr1,need_ylmgr1) 937 ABI_MALLOC(ylmgr1_k,(npw1_k,3,psps%mpsang*psps%mpsang*psps%useylm*dimylmgr1)) 938 939 ! Get plane-wave vectors and related data at k 940 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 941 if (psps%useylm==1) then 942 do jj=1,psps%mpsang*psps%mpsang 943 ylm_k(1:npw_k,jj)=ylm(1+ikg:npw_k+ikg,jj) 944 end do 945 end if 946 947 ! Get plane-wave vectors and related data at k+q 948 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 949 if (psps%useylm==1) then 950 do jj=1,psps%mpsang*psps%mpsang 951 ylm1_k(1:npw1_k,jj)=ylm1(1+ikg1:npw1_k+ikg1,jj) 952 end do 953 if (need_ylmgr1==1.and.useylmgr1/=0) then 954 do jj=1,psps%mpsang*psps%mpsang 955 do ia=1,3 956 ylmgr1_k(1:npw1_k,ia,jj)=ylmgr1(1+ikg1:npw1_k+ikg1,ia,jj) 957 end do 958 end do 959 end if 960 end if 961 962 ! If Ylm gradients at k+q are needed and not in memory, compute them 963 if (need_ylmgr1==1.and.useylmgr1==0) then 964 option=-1;npwar1_tmp(1)=npw1_k;nband_tmp(1)=nband_k 965 !Subtlety: initylmg is called in sequential mode 966 call initylmg(gprimd,kg1_k,kpq,1,mpi_enreg_seq,psps%mpsang,& 967 & npw1_k,nband_tmp,1,npwar1_tmp,nsppol,option,rprimd,& 968 & ylm1_k,ylmgr1_k) 969 end if 970 971 ! Compute (k+G) vectors 972 nkpg=0;if(ipert1<=dtset%natom) nkpg=3*dtset%nloalg(3) 973 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 974 if (nkpg>0) then 975 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 976 end if 977 978 ! Compute (k+q+G) vectors 979 nkpg1=0;if(ipert1<=dtset%natom.or.need_ylmgr1==1) nkpg1=3*dtset%nloalg(3) 980 ABI_MALLOC(kpg1_k,(npw1_k,nkpg1)) 981 if (nkpg1>0) then 982 call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k) 983 end if 984 985 ! Allocate kinetic contributions 986 ABI_MALLOC(dkinpw,(npw_k)) 987 ABI_MALLOC(kinpw1,(npw1_k)) 988 dkinpw=zero 989 ! Compute (1/2) (2 Pi)**2 (k+q+G)**2: 990 ! call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k) 991 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0) 992 993 ! Compute nonlocal form factors ffnl at (k+G), for all atoms 994 ider=0;idir0=0 995 dimffnl=0;if (ipert1<=dtset%natom) dimffnl=1 996 ABI_MALLOC(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 997 if (ipert1<=dtset%natom) then 998 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,& 999 & ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,& 1000 & psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,usepaw,& 1001 & psps%useylm,ylm_k,ylmgr_dum) 1002 end if 1003 1004 ! Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms 1005 ider=1;if (ipert1<=dtset%natom) ider=0 1006 if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)then 1007 dimffnl1=1;if (ider>=1) dimffnl1=2+5*psps%useylm 1008 idir0=0;if (ider>0.and.psps%useylm==1) idir0=-7 1009 else 1010 dimffnl1=1;if (ider>=1) dimffnl1=2+2*psps%useylm 1011 idir0=0;if (ider>0.and.psps%useylm==1) idir0=4 1012 end if 1013 ABI_MALLOC(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat)) 1014 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,& 1015 & ider,idir0,psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,& 1016 & psps%mqgrid_ff,nkpg1,npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,usepaw,& 1017 & psps%useylm,ylm1_k,ylmgr1_k) 1018 1019 ! Extract non-local form factors for H^(j1) 1020 if (ipert1<=dtset%natom) then 1021 ffnl1_idir1 => ffnl1(:,:,:,:) 1022 dimffnl1_idir1=dimffnl1 1023 else 1024 dimffnl1_idir1=1+ider 1025 ABI_MALLOC(ffnl1_idir1,(npw1_k,dimffnl1_idir1,psps%lmnmax,psps%ntypat)) 1026 ii=1;if (psps%useylm==0) ii=1+ider 1027 do itypat=1,psps%ntypat 1028 do ilmn=1,psps%lmnmax 1029 ffnl1_idir1(1:npw1_k,1:ii,ilmn,itypat)=ffnl1(1:npw1_k,1:ii,ilmn,itypat) 1030 end do 1031 end do 1032 end if 1033 1034 ! Load k-dependent part in the Hamiltonian datastructure 1035 ABI_MALLOC(ph3d,(2,npw_k,gs_hamkq%matblk)) 1036 call gs_hamkq%load_k(kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,& 1037 & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.) 1038 if (size(ffnlk)>0) then 1039 call gs_hamkq%load_k(ffnl_k=ffnlk) 1040 else 1041 call gs_hamkq%load_k(ffnl_k=ffnl1) 1042 end if 1043 1044 ! Load k+q-dependent part in the Hamiltonian datastructure 1045 ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead) 1046 call gs_hamkq%load_kprime(kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,& 1047 & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1_idir1,& 1048 & compute_gbound=.true.) 1049 if (qne0) then 1050 ABI_MALLOC(ph3d1,(2,npw1_k,gs_hamkq%matblk)) 1051 call gs_hamkq%load_kprime(ph3d_kp=ph3d1,compute_ph3d=.true.) 1052 end if 1053 1054 ! Load k-dependent part in the 1st-order Hamiltonian datastructure 1055 call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dkinpw) 1056 1057 ! Allocate memory space for one band 1058 if (need_wfk) then 1059 ABI_MALLOC(cwave0,(2,npw_k*nspinor)) 1060 end if 1061 if (need_wf1) then 1062 ABI_MALLOC(cwavef,(2,npw1_k*nspinor)) 1063 end if 1064 ABI_MALLOC(gh1,(2,npw1_k*nspinor)) 1065 nullify(cwaveprj0_idir1) 1066 if (usecprj==1) then 1067 ABI_MALLOC(cwaveprj0,(dtset%natom,nspinor)) 1068 call pawcprj_alloc(cwaveprj0,ncpgr,gs_hamkq%dimcprj) 1069 end if 1070 if (has_dcwf) then 1071 ABI_MALLOC(gs1,(2,npw1_k*nspinor)) 1072 else 1073 ABI_MALLOC(gs1,(0,0)) 1074 end if 1075 1076 ABI_MALLOC(band_procs, (nband_k)) 1077 call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,nband_k,& 1078 & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band) 1079 1080 ! LOOP OVER BANDS 1081 iband_me = 0 1082 do iband=1,nband_k 1083 1084 ! Skip band if not to be treated by this proc 1085 if (xmpi_paral==1) then 1086 if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle 1087 end if 1088 iband_me = iband_me + 1 1089 1090 ch1c_tmp = zero 1091 1092 bands_treated_now = 0 1093 bands_treated_now(iband) = 1 1094 call xmpi_sum(bands_treated_now,mpi_enreg%comm_band,ierr) 1095 1096 ! Extract GS wavefunctions 1097 if (need_wfk) then 1098 cwave0(:,:)=cg(:,1+(iband_me-1)*npw_k*nspinor+icg:iband_me*npw_k*nspinor+icg) 1099 if (usecprj==1) then 1100 call pawcprj_get(gs_hamkq%atindx1,cwaveprj0,cprj,dtset%natom,iband_me,ibg,ikpt,iorder_cprj,& 1101 & isppol,mband_mem_rbz,mkmem,dtset%natom,1,nband_me,nspinor,nsppol,dtfil%unpaw) 1102 ! in distributed cprj memory, no need for these? cg and cprj have same distribution 1103 !& mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 1104 end if 1105 end if 1106 1107 ! Extract 1st-order wavefunctions 1108 if (need_wf1) then 1109 cwavef(:,:)=cg1(:,1+(iband_me-1)*npw1_k*nspinor+icg1:iband_me*npw1_k*nspinor+icg1) 1110 end if 1111 1112 ! LOOP OVER PERTURBATION DIRECTIONS 1113 do kdir1=1,mdir1 1114 idir1=jdir1(kdir1) 1115 istr1=idir1;if(ipert1==dtset%natom+4) istr1=idir1+3 1116 1117 ! Not able to compute if ipert1=(Elect. field) and no ddk WF file 1118 if (ipert1==dtset%natom+2.and.ddkfil(idir1)==0) cycle 1119 1120 ! Extract 1st-order NL form factors derivatives for this idir1 1121 if (dimffnl1_idir1>=2.and.psps%useylm==1.and.ipert1>dtset%natom) then 1122 do itypat=1,psps%ntypat 1123 do ilmn=1,psps%lmnmax 1124 ffnl1_idir1(1:npw1_k,2,ilmn,itypat)=ffnl1(1:npw1_k,1+istr1,ilmn,itypat) 1125 end do 1126 end do 1127 end if 1128 1129 ! Extract ground state projected WF and derivatives in idir1 direction 1130 if (usecprj==1) then 1131 cpopt=1 1132 1133 ! === Atomic displ. perturbation 1134 if (ipert<=dtset%natom) then 1135 ABI_MALLOC(cwaveprj0_idir1,(dtset%natom,nspinor)) 1136 call pawcprj_alloc(cwaveprj0_idir1,1,gs_hamkq%dimcprj) 1137 if (ipert1<=dtset%natom) then 1138 call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=idir1) 1139 else 1140 if (ipert1==dtset%natom+2) then 1141 idir_cprj=idir1;choice=5 1142 end if 1143 if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 1144 idir_cprj=istr1;choice=3 1145 end if 1146 call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1) 1147 call getcprj(choice,cpopt,cwave0,cwaveprj0_idir1,& 1148 & gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,& 1149 & gs_hamkq%kg_kp,gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,& 1150 & gs_hamkq%mgfft,mpi_enreg,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,& 1151 & gs_hamkq%nloalg,gs_hamkq%npw_kp,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,& 1152 & gs_hamkq%ph1d,gs_hamkq%ph3d_kp,gs_hamkq%ucvol,gs_hamkq%useylm) 1153 end if 1154 1155 ! === Wave-vector perturbation 1156 else if (ipert==dtset%natom+1) then 1157 cwaveprj0_idir1 => cwaveprj0 1158 1159 ! == Electric field perturbation 1160 else if (ipert==dtset%natom+2) then 1161 ABI_MALLOC(cwaveprj0_idir1,(dtset%natom,nspinor)) 1162 call pawcprj_alloc(cwaveprj0_idir1,1,gs_hamkq%dimcprj) 1163 if (ipert1==dtset%natom+2) then 1164 call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=idir1) 1165 else 1166 if (ipert1<=dtset%natom) then 1167 idir_cprj=idir1;choice=2 1168 end if 1169 if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 1170 idir_cprj=istr1;choice=3 1171 end if 1172 call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1) 1173 call getcprj(choice,cpopt,cwave0,cwaveprj0_idir1,& 1174 & gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,& 1175 & gs_hamkq%kg_kp,gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,& 1176 & gs_hamkq%mgfft,mpi_enreg,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,gs_hamkq%nloalg,& 1177 & gs_hamkq%npw_kp,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,gs_hamkq%ph1d,& 1178 & gs_hamkq%ph3d_kp,gs_hamkq%ucvol,gs_hamkq%useylm) 1179 end if 1180 1181 ! === Strain perturbation 1182 else if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then 1183 if ((ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and.(istr==istr1)) then 1184 cwaveprj0_idir1 => cwaveprj0 1185 else 1186 ABI_MALLOC(cwaveprj0_idir1,(dtset%natom,nspinor)) 1187 call pawcprj_alloc(cwaveprj0_idir1,ncpgr,gs_hamkq%dimcprj) 1188 if (ipert1<=dtset%natom) then 1189 idir_cprj=idir1;choice=2 1190 end if 1191 if (ipert1==dtset%natom+2) then 1192 idir_cprj=idir1;choice=5 1193 end if 1194 if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 1195 idir_cprj=istr1;choice=3 1196 end if 1197 call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1) 1198 call getcprj(choice,cpopt,cwave0,cwaveprj0_idir1,& 1199 & gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,gs_hamkq%kg_kp,& 1200 & gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,gs_hamkq%mgfft,mpi_enreg,& 1201 & gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,gs_hamkq%nloalg,gs_hamkq%npw_kp,& 1202 & gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,gs_hamkq%ph1d,gs_hamkq%ph3d_kp,& 1203 & gs_hamkq%ucvol,gs_hamkq%useylm) 1204 end if 1205 end if ! ipert 1206 1207 else ! usecprj=0: cwaveprj0_idir1 is not used 1208 cwaveprj0_idir1 => cwaveprj0 1209 end if 1210 1211 ! Eventually compute 1st-order kinetic operator 1212 if (ipert1==dtset%natom+1) then 1213 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir1,0) 1214 else if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then 1215 call kpgstr(dkinpw,dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,gprimd,istr1,kg_k,kpoint,npw_k) 1216 end if 1217 1218 ! Finalize initialization of 1st-order NL hamiltonian 1219 if (need_pawij10) rf_hamkq%e1kbfr => e1kbfr(:,:,:,:,idir1) 1220 1221 ! Read DDK wave function (if ipert1=electric field) 1222 if (ipert1==dtset%natom+2) then 1223 usevnl=1 1224 ABI_MALLOC(gvnlx1,(2,npw1_k*nspinor*usevnl)) 1225 if (need_ddk_file) then 1226 if (ddkfil(idir1)/=0) then 1227 !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt) 1228 !ik_ddk = indkpt1(ikpt) 1229 !call ddks(idir1)%read_bks(iband, ik_ddk, isppol, xmpio_single, cg_bks=gvnlx1) 1230 gvnlx1 = cg_ddk(:,1+(iband_me-1)*npw1_k*nspinor:iband_me*npw1_k*nspinor,idir1) 1231 else 1232 gvnlx1=zero 1233 end if 1234 if (ipert1==dtset%natom+2) then 1235 do ii=1,npw1_k*nspinor ! Multiply ddk by +i (to be consistent with getgh1c) 1236 arg=gvnlx1(1,ii) 1237 gvnlx1(1,ii)=-gvnlx1(2,ii) 1238 gvnlx1(2,ii)=arg 1239 end do 1240 end if 1241 else 1242 gvnlx1=zero 1243 end if 1244 else 1245 usevnl=0 1246 ABI_MALLOC(gvnlx1,(2,npw1_k*nspinor*usevnl)) 1247 end if 1248 1249 ! Get |H^(j2)-Eps_k_i.S^(j2)|u0_k_i> (VHxc-dependent part not taken into account) and S^(j2)|u0> 1250 lambda=eig_k(iband);berryopt=0;optlocal=0 1251 optnl=0;if (ipert1/=dtset%natom+1.or.idir==idir1) optnl=1 1252 opt_gvnlx1=0;if (ipert1==dtset%natom+2) opt_gvnlx1=2 1253 sij_opt=-1;if (has_dcwf) sij_opt=1 1254 if (usepaw==0) sij_opt=0 1255 call getgh1c(berryopt,cwave0,cwaveprj0_idir1,gh1,dum1,gs1,gs_hamkq,gvnlx1,idir1,ipert1,& 1256 & lambda,mpi_enreg,optlocal,optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 1257 if (sij_opt==1.and.optnl==1) gh1=gh1-lambda*gs1 1258 ABI_FREE(gvnlx1) 1259 1260 ! If needed, compute here <delta_u^(j1)_k_i|H^(j2)-Eps_k_i.S^(j2)|u0_k_i> 1261 ! with delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>] 1262 ! (see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42)) 1263 ! This can be rewritten as: 1264 ! -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> 1265 ! The sum over j can be computed with a single call to projbd routine 1266 ! At first call (when j1=j2), ch1c=<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i> is stored 1267 ! For the next calls, it is re-used. 1268 if (has_dcwf.or.(ipert==ipert1.and.idir==idir1.and.usepaw==1)) then 1269 ! note: gvnlx1 used as temporary space 1270 ABI_MALLOC(gvnlx1,(2,npw1_k*nspinor)) 1271 ABI_MALLOC(gvnlx1_tmp,(2,npw1_k*nspinor)) 1272 if (ipert==ipert1.and.idir==idir1) then 1273 option=0;gvnlx1=gh1 1274 else 1275 option=1;gvnlx1=zero 1276 end if 1277 do iband_=1, nband_k 1278 if (bands_treated_now(iband_) == 0) cycle 1279 1280 ! distribute gvnlx1 to my subcomm 1281 gvnlx1_tmp = zero 1282 if (iband_ == iband) then 1283 gvnlx1_tmp = gvnlx1 1284 end if 1285 ! TODO CHECK IF IT IS BAND_PROCS(IBAND_) 1286 !call xmpi_bcast(gvnlx1_tmp, band_procs(iband_), mpi_enreg%comm_band, ierr) 1287 call xmpi_sum(gvnlx1_tmp, mpi_enreg%comm_band, ierr) 1288 if (option == 1) then 1289 ! in case I need to reuse the ch1c (option 1) then load them here 1290 ch1c_tmp(:,1:nband_me) = ch1c(:,1:nband_me,iband_,ikpt_me) 1291 end if 1292 1293 1294 ! Compute -Sum_{j}[<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>.|u0_k+q_j> 1295 call projbd(cgq,gvnlx1_tmp,-1,icgq,0,istwf_k,mcgq,0,nband_me,npw1_k,nspinor,& 1296 & dum1,ch1c_tmp,option,tim_projbd,0,mpi_enreg%me_g0,mpi_enreg%comm_fft) 1297 1298 !sum over all jband by combining the projbd 1299 call xmpi_sum(gvnlx1_tmp,mpi_enreg%comm_band,ierr) 1300 ! keep my own gvnlx 1301 if (iband_ == iband) then 1302 ! if bands are parallelized, I have only projected against bands on my cpu 1303 ! Pc|work> = |work> - Sum_l <psi_{k+q, l}|work> |psi_{k+q, l}> 1304 ! = Sum_nproc_band (|work> - Sum_{my l} <psi_{k+q, l}|work> |psi_{k+q, l}>) - (nproc_band-1) |work> 1305 !TODO: make this a blas call? zaxpy 1306 gvnlx1 = gvnlx1_tmp - (mpi_enreg%nproc_band-1)*gvnlx1 1307 end if 1308 1309 if (option == 0) then 1310 ! save ch1c for all of the iband_ on each proc, for later use. First band index only for my nband_me which matches cgq 1311 ch1c(:,1:nband_me,iband_,ikpt_me) = ch1c_tmp(:,1:nband_me) 1312 end if 1313 end do ! iband_ 1314 ABI_FREE(gvnlx1_tmp) 1315 1316 if (has_dcwf) then 1317 if (ipert==ipert1.and.idir==idir1) gvnlx1=gvnlx1-gh1 1318 dotr=zero;doti=zero 1319 if (abs(occ_k(iband))>tol8) then 1320 ! 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> 1321 call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,gs1,gvnlx1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1322 ! Add contribution to DDB 1323 ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2 1324 ! (-1) factor already present 1325 d2ovl_k(1,idir1)=d2ovl_k(1,idir1)+wtk_k*occ_k(iband)*elfd_fact*dotr 1326 d2ovl_k(2,idir1)=d2ovl_k(2,idir1)+wtk_k*occ_k(iband)*elfd_fact*doti 1327 end if 1328 end if 1329 ABI_FREE(gvnlx1) 1330 end if 1331 1332 ! If needed, compute here <delta_u^(j1)_k_i|H-Eps_k_i.S|u^(j2)_k_i> 1333 ! This is equal to <delta_u^(j1)_k_i|H-Eps_k_i.S|delta_u^(j2)_k_i> (I) 1334 ! +<delta_u^(j1)_k_i|H-Eps_k_i.S|u^paral^(j2)_k_i> (II) 1335 ! (u^paral^(j2)_k_i is the part of u^(j2)_k_i parallel to active space : metals) 1336 ! (I) can be rewritten as: 1337 ! 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) } 1338 ! (II) can be rewritten as: 1339 ! Sum_j{1/2.(occ_kq_j-occ_k_i).Eps1_k,q_ij.<u0_k_i|S^(j1)|u0_k+q_j> } 1340 ! 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> 1341 ! At first call (when j1=j2), cs1c=<u0_k_i|S^(j1)|u0_k+q_j> is stored 1342 ! For the next calls, it is re-used. 1343 if (has_dcwf.and.is_metal_or_qne0) then 1344 ABI_MALLOC(gvnlx1,(2,npw1_k*nspinor)) 1345 ! dotX is local to my proc, and should accumulate sum over all jband, for my iband_me 1346 dotr=zero;doti=zero 1347 ! flag to broadcast the j dependent vector in the band pool, to get full sum over j 1348 do_bcast = 0 1349 ! flag to do scalar product: only needed if we are saving cs1c or if the band is occupied 1350 do_scprod = 0 1351 if ((ipert==ipert1.and.idir==idir1)) then 1352 do_bcast = 1 1353 do_scprod = 1 1354 end if 1355 invocc=zero 1356 if (abs(occ_k(iband))>tol8) then 1357 invocc=two/occ_k(iband) 1358 do_scprod = 1 1359 end if 1360 ! does anyone else need the cgq(j) below? 1361 do iband_ = 1, nband_k 1362 if (bands_treated_now(iband_) > 0 .and. abs(occ_k(iband_))>tol8) then 1363 do_bcast = 1 1364 end if 1365 end do 1366 1367 1368 jband_me = 0 1369 do jband=1,nband_k 1370 if (do_bcast > 0) then 1371 gvnlx1 = zero 1372 if (mpi_enreg%proc_distrb(ikpt,jband,isppol)==me) then 1373 jband_me = jband_me + 1 1374 ! gvnlx1 depends on j only, I have it, and everyone needs it 1375 gvnlx1(:,1:npw1_k*nspinor)=cgq(:,1+npw1_k*nspinor*(jband_me-1)+icgq:npw1_k*nspinor*jband_me+icgq) 1376 end if 1377 ! xmpi bcast the current jband to other procs in band pool 1378 !call xmpi_bcast(gvnlx1, band_procs(jband), mpi_enreg%comm_band, ierr) 1379 call xmpi_sum(gvnlx1, mpi_enreg%comm_band, ierr) 1380 ! Computation of cs1c=<u0_k_i|S^(j1)|u0_k+q_j> 1381 ! do _I_ need to calculate the dot1X? 1382 if (do_scprod > 0) then 1383 call dotprod_g(dot1r,dot1i,istwf_k,npw1_k*nspinor,2,gs1,gvnlx1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1384 if (ipert==ipert1.and.idir==idir1.and.has_dcwf2) then 1385 cs1c(1,jband,iband_me,ikpt_me)=dot1r 1386 cs1c(2,jband,iband_me,ikpt_me)=dot1i 1387 end if 1388 end if 1389 end if ! ipert==ipert1.and.idir==idir1 or some iband for some proc in pool is filled 1390 1391 if (abs(occ_k(iband))>tol8) then 1392 ! Computation of term (I) 1393 if (has_dcwf2) then 1394 arg=eig_kq(jband)-eig_k(iband) 1395 dot2r=cs1c(1,jband,iband_me,ikpt_me) 1396 dot2i=cs1c(2,jband,iband_me,ikpt_me) 1397 dotr=dotr+(dot1r*dot2r+dot1i*dot2i)*arg 1398 doti=doti+(dot1i*dot2r-dot1r*dot2i)*arg 1399 end if 1400 ! Computation of term (II) TODO: the next two ifs could be combined 1401 if (is_metal) then 1402 if (abs(rocceig(jband,iband))>tol8) then 1403 ii=2*jband-1+(iband-1)*2*nband_k 1404 arg=invocc*rocceig(jband,iband)*(eig_k(iband)-eig_kq(jband)) 1405 dot2r=eig1_k(ii) 1406 dot2i=eig1_k(ii+1) 1407 dotr=dotr+arg*(dot1r*dot2r-dot1i*dot2i) 1408 doti=doti+arg*(dot1r*dot2i+dot2r*dot1i) 1409 end if 1410 end if 1411 end if ! occ bands 1412 end do ! jband 1413 1414 dotr=quarter*dotr 1415 doti=quarter*doti 1416 1417 ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) 1418 ! Note2: do not sum over bands here - comm_band is a sub communicator of spacecomm, and a full sum is done later 1419 d2ovl_k(1,idir1)=d2ovl_k(1,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*dotr 1420 d2ovl_k(2,idir1)=d2ovl_k(2,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*doti 1421 ABI_FREE(gvnlx1) 1422 end if 1423 1424 ! Build the matrix element <u0_k_i|H^(j1)-Eps_k_i.S^(j1)|u^(j2)_k,q_i> 1425 ! and add contribution to DDB 1426 if (ipert1/=dtset%natom+1) then 1427 if (abs(occ_k(iband))>tol8) then 1428 call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,gh1,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1429 ! Case ipert1=natom+2 (electric field): 1430 ! gh1 contains H^(j1)|u0_k_i> (VHxc constant) which corresponds 1431 ! to i.d/dk in Eq. (38) of Gonze, PRB 55, 10355 (1997) [[cite:Gonze1997a]]. 1432 ! * if ipert==natom+2, we apply directly Eq. (38) 1433 ! * if ipert/=natom+2, Born effective charges are minus D2E 1434 d2nl_k(1,idir1)=d2nl_k(1,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*dotr 1435 d2nl_k(2,idir1)=d2nl_k(2,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*doti 1436 end if 1437 1438 ! Or compute localisation tensor (ddk) 1439 ! See M. Veithen thesis Eq(2.5) 1440 ! MT jan-2010: this is probably not correctly implemented for PAW !!! 1441 ! missing terms due to S^(1) and S^(2) 1442 else 1443 ! note: gh1 used as temporary space (to store idir ddk WF) 1444 if (idir==idir1) then 1445 gh1=cwavef 1446 else 1447 if (need_ddk_file.and.ddkfil(idir1)/=0) then 1448 !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt) 1449 !ik_ddk = indkpt1(ikpt) 1450 !call ddks(idir1)%read_bks(iband, ik_ddk, isppol, xmpio_single, cg_bks=gh1) 1451 gh1 = cg_ddk(:,1+(iband_me-1)*npw1_k*nspinor:iband_me*npw1_k*nspinor,idir1) 1452 else 1453 gh1=zero 1454 end if 1455 end if 1456 if (abs(occ_k(iband))>tol8) then 1457 call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,gh1,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1458 call dotprod_g(dot1r,dot1i,istwf_k,npw1_k*nspinor,2,cwave0,gh1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1459 call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*nspinor,2,cwavef,cwave0,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1460 dotr=dotr-(dot1r*dot2r-dot1i*dot2i) 1461 doti=doti-(dot1r*dot2i+dot1i*dot2r) 1462 d2nl_k(1,idir1)=d2nl_k(1,idir1)+wtk_k*occ_k(iband)*dotr/(nband_kocc*two) 1463 d2nl_k(2,idir1)=d2nl_k(2,idir1)+wtk_k*occ_k(iband)*doti/(nband_kocc*two) 1464 end if 1465 end if 1466 1467 ! Accumulate here 1st-order density change due to overlap operator changes (if any) 1468 if (has_drho) then 1469 ! Compute here delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>] 1470 ! (see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42)) 1471 ABI_MALLOC(dcwavef,(2,npw1_k*nspinor)) 1472 ABI_MALLOC(dcwaveprj,(dtset%natom,nspinor)) 1473 call pawcprj_alloc(dcwaveprj,0,gs_hamkq%dimcprj) 1474 ! NB: have to call getdc with all band processors to distribute cgq cprjq correctly 1475 call getdc1(iband,band_procs,bands_treated_now,cgq,cprjq,dcwavef,dcwaveprj,& 1476 & ibgq,icgq,istwf_k,mcgq,& 1477 & mcprjq,mpi_enreg,dtset%natom,nband_k,nband_me,npw1_k,nspinor,1,gs1) 1478 1479 if (abs(occ_k(iband))>tol8) then 1480 ! Accumulate 1st-order density due to delta_u^(j1) 1481 option=1;wfcorr=0 1482 call dfpt_accrho(cplex,cwave0,dcwavef,dcwavef,cwaveprj0_idir1,dcwaveprj,& 1483 & lambda,gs_hamkq,iband,idir1,ipert1,isppol,dtset%kptopt,& 1484 & mpi_enreg,dtset%natom,nband_k,1,npw_k,npw1_k,nspinor,occ_k,option,& 1485 & pawdrhoij1_unsym(:,idir1),drhoaug1(:,:,:,idir1),tim_fourwf,wfcorr,wtk_k) 1486 end if 1487 1488 call pawcprj_free(dcwaveprj) 1489 ABI_FREE(dcwaveprj) 1490 ABI_FREE(dcwavef) 1491 end if ! has_drho 1492 1493 if((usecprj==1).and..not.(associated(cwaveprj0_idir1,cwaveprj0)))then 1494 call pawcprj_free(cwaveprj0_idir1) 1495 ABI_FREE(cwaveprj0_idir1) 1496 end if 1497 1498 ! End of loops 1499 end do ! idir1 1500 end do ! iband 1501 1502 ABI_FREE(band_procs) 1503 1504 ! Accumulate contribution of this k-point 1505 d2nl (:,:,ipert1,idir,ipert)=d2nl (:,:,ipert1,idir,ipert)+d2nl_k (:,:) 1506 if (usepaw==1) d2ovl(:,:,ipert1,idir,ipert)=d2ovl(:,:,ipert1,idir,ipert)+d2ovl_k(:,:) 1507 1508 1509 ! Deallocations of arrays used for this k-point 1510 ABI_FREE(gh1) 1511 ABI_FREE(gs1) 1512 if (need_wfk) then 1513 ABI_FREE(cwave0) 1514 end if 1515 if (need_wf1) then 1516 ABI_FREE(cwavef) 1517 end if 1518 ABI_FREE(kg_k) 1519 ABI_FREE(kg1_k) 1520 ABI_FREE(ylm_k) 1521 ABI_FREE(ylm1_k) 1522 ABI_FREE(ylmgr1_k) 1523 ABI_FREE(kpg_k) 1524 ABI_FREE(kpg1_k) 1525 ABI_FREE(d2nl_k) 1526 ABI_FREE(d2ovl_k) 1527 ABI_FREE(eig_k) 1528 ABI_FREE(eig_kq) 1529 ABI_FREE(eig1_k) 1530 ABI_FREE(occ_k) 1531 if (is_metal) then 1532 ABI_FREE(doccde_k) 1533 ABI_FREE(doccde_kq) 1534 ABI_FREE(occ_kq) 1535 ABI_FREE(rocceig) 1536 end if 1537 ABI_FREE(dkinpw) 1538 ABI_FREE(kinpw1) 1539 ABI_FREE(ph3d) 1540 if (allocated(ph3d1)) then 1541 ABI_FREE(ph3d1) 1542 end if 1543 ABI_FREE(ffnlk) 1544 ABI_FREE(ffnl1) 1545 if (ipert1>dtset%natom) then 1546 ABI_FREE(ffnl1_idir1) 1547 end if 1548 nullify(ffnl1_idir1) 1549 if (usecprj==1) then 1550 call pawcprj_free(cwaveprj0) 1551 ABI_FREE(cwaveprj0) 1552 end if 1553 nullify(cwaveprj0_idir1) 1554 ! Shift arrays 1555 bdtot_index=bdtot_index+nband_k 1556 bd2tot_index=bd2tot_index+2*nband_k**2 1557 if (mkmem/=0) then 1558 ibg=ibg+nspinor*nband_me 1559 icg=icg+npw_k*nspinor*nband_me 1560 ikg=ikg+npw_k 1561 end if 1562 if (mkqmem/=0) then 1563 ibgq=ibgq+nspinor*nband_me 1564 icgq=icgq+npw1_k*nspinor*nband_me 1565 end if 1566 if (mk1mem/=0) then 1567 ibg1=ibg1+nspinor*nband_me 1568 icg1=icg1+npw1_k*nspinor*nband_me 1569 ikg1=ikg1+npw1_k 1570 end if 1571 1572 1573 end do ! End loop over K-POINTS 1574 !---------------------------------------------------------------- 1575 1576 ! Transfer 1st-order density change due to overlap; also take into account the spin. 1577 if(has_drho) then 1578 do kdir1=1,mdir1 1579 idir1=jdir1(kdir1) 1580 call fftpac(isppol,mpi_enreg,nspden,cplex*dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),& 1581 & cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),& 1582 & dtset%ngfft,drho1wfr(:,:,idir1),drhoaug1(:,:,:,idir1),1) 1583 end do 1584 end if 1585 1586 end do ! End loop over SPINS 1587 !---------------------------------------------------------------- 1588 1589 ! Free memory used for this type of perturbation 1590 call rf_hamkq%free() 1591 if (has_vectornd) then 1592 ABI_FREE(vectornd_pac_idir) 1593 end if 1594 if (has_drho) then 1595 ABI_FREE(drhoaug1) 1596 end if 1597 if (need_pawij10) then 1598 do kdir1=1,mdir1 1599 idir1=jdir1(kdir1) 1600 call paw_ij_free(paw_ij10(:,idir1)) 1601 end do 1602 ABI_FREE(e1kbfr_spin) 1603 end if 1604 ABI_FREE(paw_ij10) 1605 1606 ! In case of parallelism, sum 1st-order density and occupation matrix over processors 1607 if (has_drho.and.xmpi_paral==1) then 1608 1609 ! Accumulate 1st-order density 1610 call timab(48,1,tsec) 1611 bufsz=cplex*dtset%nfft*nspden*mdir1 1612 ABI_MALLOC(buffer,(bufsz)) 1613 buffer(1:bufsz)=reshape(drho1wfr,(/bufsz/)) 1614 call xmpi_sum(buffer,bufsz,spaceworld,ierr) 1615 drho1wfr(:,:,:)=reshape(buffer(1:bufsz),(/cplex*dtset%nfft,nspden,mdir1/)) 1616 ABI_FREE(buffer) 1617 call timab(48,2,tsec) 1618 1619 ! Accumulate 1st-order PAW occupancies 1620 if (usepaw==1) then 1621 call pawrhoij_mpisum_unpacked(pawdrhoij1_unsym,spaceworld) 1622 end if 1623 1624 end if 1625 1626 ! Compute second part of overlap contribution (due to VHxc^(j2)(tild_n+hat_n)) 1627 if (has_drho) then 1628 1629 ABI_MALLOC(drhor1,(cplex*nfftf,nspden)) 1630 ABI_MALLOC(dnhat1,(cplex*nfftf,nspden)) 1631 1632 ! LOOP OVER PERTURBATION DIRECTIONS 1633 do kdir1=1,mdir1 1634 idir1=jdir1(kdir1) 1635 1636 ! Build and symmetrize 1st-order density change due to change of overlap 1637 ABI_MALLOC(drho1wfg,(2,dtset%nfft)) 1638 call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,& 1639 & nspden,nsppol,nsym1,phnons1,drho1wfg,drho1wfr(:,:,idir1),& 1640 & rprimd,symaf1,symrl1,tnons1) 1641 if (dtset%pawstgylm/=0) then 1642 option=0 1643 call pawnhatfr(option,idir1,ipert1,my_natom,dtset%natom,nspden,dtset%ntypat,& 1644 & pawang,pawfgrtab,pawrhoij,pawtab,rprimd,& 1645 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 1646 end if 1647 call pawmkrho(1,arg,cplex,gprimd,idir1,indsy1,ipert1,& 1648 & mpi_enreg,my_natom,dtset%natom,nspden,nsym1,dtset%ntypat,dtset%paral_kgb,pawang,& 1649 & pawfgr,pawfgrtab,-10001,pawdrhoij1(:,idir1),pawdrhoij1_unsym(:,idir1),pawtab,& 1650 & dtset%qptn,drho1wfg,drho1wfr(:,:,idir1),drhor1,rprimd,symaf1,symrc1,dtset%typat,& 1651 & ucvol,dtset%usewvl,xred,pawang_sym=pawang1,pawnhat=dnhat1,pawrhoij0=pawrhoij) 1652 ABI_FREE(drho1wfg) 1653 1654 ! Compute plane-wave contribution to overlap contribution 1655 ! This is subtle as it is a mix of Eq(79) and Eq(80) of PRB 78, 035105 (2008) [[cite:Audouze2008]] 1656 ! Details: 1657 ! The VH(tild_nZc)^(1) term of Eq(79) is: 1658 ! <VH(tild_nZc)^(j2)|delta_tild_rho^(j1)> = <vpsp1|drhor1-dnhat1> 1659 ! The first term of Eq(80) is: 1660 ! <VHxc^(j2)|delta_tild_rho^(j1)+delta_hat_rho^(j1)> = <vtrial1-vpsp1|drhor1> 1661 ! The addition of these two terms gives: 1662 ! <vtrial1|drhor1>-<vpsp1|dnhat1> 1663 ! And this is more subtle when usexcnhat=0 1664 call dotprod_vn(cplex,drhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vtrial1,ucvol) 1665 if (usexcnhat/=0) then 1666 call dotprod_vn(cplex,dnhat1,dot2r,dot2i,nfftf,nfftot,1 ,2,vpsp1,ucvol) 1667 else 1668 ABI_MALLOC(vtmp1,(cplex*nfftf,nspden)) 1669 do ispden=1,nspden 1670 vtmp1(:,ispden)=vtrial1(:,ispden)-vhartr1(:) 1671 end do 1672 call dotprod_vn(cplex,dnhat1,dot2r,dot2i,nfftf,nfftot,nspden,2,vtmp1,ucvol) 1673 ABI_FREE(vtmp1) 1674 end if 1675 dotr=dot1r-dot2r;doti=dot1i-dot2i 1676 1677 ! Compute on-site contributions to overlap contribution 1678 ! (two last terms of Eq(80) of PRB 78, 035105 (2008)) [[cite:Audouze2008]] 1679 ! (note: Dij^(j2) and Vxc^(j2) are computed for ipert at first call) 1680 call pawdfptenergy(epawnst,ipert,ipert1,dtset%ixc,my_natom,dtset%natom,dtset%ntypat,& 1681 & nzlmopt_ipert,nzlmopt_ipert1,paw_an,paw_an1,paw_ij1,pawang,dtset%pawprtvol,& 1682 & pawrad,pawrhoij1,pawdrhoij1(:,idir1),pawtab,dtset%pawxcdev,dtset%xclevel,& 1683 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 1684 1685 ! Accumulate in 2nd-order matrix: 1686 ! Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2 1687 ! has to take the complex conjugate because we want here Int[VHxc^(j1)^*.delta_rho^(j2)] 1688 dotr=dotr+epawnst(1);doti=-(doti+epawnst(2)) 1689 d2ovl_drho(1,idir1,ipert1,idir,ipert)=elfd_fact*dotr 1690 d2ovl_drho(2,idir1,ipert1,idir,ipert)=elfd_fact*doti 1691 1692 end do ! End loop over perturbation directions 1693 1694 ! Free no more needed memory 1695 ABI_FREE(drhor1) 1696 ABI_FREE(dnhat1) 1697 ABI_FREE(drho1wfr) 1698 do iatom=1,my_natom 1699 if (pawfgrtab(iatom)%nhatfr_allocated>0) then 1700 ABI_FREE(pawfgrtab(iatom)%nhatfr) 1701 end if 1702 pawfgrtab(iatom)%nhatfr_allocated=0 1703 end do 1704 if (paral_atom) then 1705 do kdir1=1,mdir1 1706 idir1=jdir1(kdir1) 1707 call pawrhoij_free(pawdrhoij1_unsym(:,idir1)) 1708 end do 1709 ABI_FREE(pawdrhoij1_unsym) 1710 end if 1711 do kdir1=1,mdir1 1712 idir1=jdir1(kdir1) 1713 call pawrhoij_free(pawdrhoij1(:,idir1)) 1714 end do 1715 ABI_FREE(pawdrhoij1) 1716 end if ! has_drho 1717 1718 ! End loop over perturbations (j1) 1719 end do 1720 1721 !Final deallocations 1722 ABI_FREE(ch1c) 1723 ABI_FREE(ch1c_tmp) 1724 if (usepaw==1.and.is_metal_or_qne0) then 1725 ABI_FREE(cs1c) 1726 ABI_FREE(cs1c_tmp) 1727 end if 1728 call gs_hamkq%free() 1729 if (has_vectornd) then 1730 ABI_FREE(vectornd_pac) 1731 end if 1732 1733 !In case of parallelism, sum over processors 1734 if (xmpi_paral==1)then 1735 call timab(161,1,tsec) 1736 call xmpi_barrier(spaceworld) 1737 call timab(161,2,tsec) 1738 ABI_MALLOC(buffer,(6*mpert*(1+usepaw))) 1739 buffer(1:6*mpert)=reshape(d2nl(:,:,:,idir,ipert),(/6*mpert/)) 1740 if (usepaw==1) buffer(6*mpert+1:6*mpert+6*mpert)=reshape(d2ovl(:,:,:,idir,ipert),(/6*mpert/)) 1741 call timab(48,1,tsec) 1742 call xmpi_sum(buffer,6*mpert*(1+usepaw),spaceworld,ierr) 1743 call timab(48,2,tsec) 1744 d2nl (:,:,:,idir,ipert)=reshape(buffer(1:6*mpert),(/2,3,mpert/)) 1745 if (usepaw==1) d2ovl(:,:,:,idir,ipert)=reshape(buffer(6*mpert+1:6*mpert+6*mpert),(/2,3,mpert/)) 1746 ABI_FREE(buffer) 1747 end if 1748 1749 !Build complete d2ovl matrix 1750 if (usepaw==1) d2ovl(:,:,:,idir,ipert)=d2ovl(:,:,:,idir,ipert)+d2ovl_drho(:,:,:,idir,ipert) 1751 1752 if (usepaw==1) then 1753 ABI_FREE(d2ovl_drho) 1754 end if 1755 1756 !Close the ddk WF files 1757 if (has_ddk_file) then 1758 do kdir1=1,mdir1 1759 idir1=jdir1(kdir1) 1760 if (ddkfil(idir1)/=0) call ddks(idir1)%close() 1761 end do 1762 ABI_FREE(cg_ddk) 1763 end if 1764 ABI_FREE(jpert1) 1765 ABI_FREE(jdir1) 1766 1767 !Symmetrize the phonons contributions, as was needed for the forces in a GS calculation 1768 ABI_MALLOC(work,(2,3,dtset%natom)) 1769 do ipert1=1,dtset%natom 1770 do idir1=1,3 1771 work(:,idir1,ipert1)=d2nl(:,idir1,ipert1,idir,ipert) 1772 end do 1773 end do 1774 call dfpt_sygra(dtset%natom,d2nl(:,:,:,idir,ipert),work,indsy1,ipert,nsym1,dtset%qptn,symrc1) 1775 if (usepaw==1) then 1776 do ipert1=1,dtset%natom 1777 do idir1=1,3 1778 work(:,idir1,ipert1)=d2ovl(:,idir1,ipert1,idir,ipert) 1779 end do 1780 end do 1781 call dfpt_sygra(dtset%natom,d2ovl(:,:,:,idir,ipert),work,indsy1,ipert,nsym1,dtset%qptn,symrc1) 1782 end if 1783 ABI_FREE(work) 1784 1785 !In the case of the strain perturbation time-reversal symmetry will always 1786 !be true so imaginary part of d2nl will be must be set to zero here since 1787 !the symmetry-reduced kpt set will leave a non-zero imaginary part. 1788 if(ipert==dtset%natom+3.or.ipert==dtset%natom+4) then 1789 d2nl(2,:,:,idir,ipert)=zero 1790 if (usepaw==1) d2ovl(2,:,:,idir,ipert)=zero 1791 else 1792 d2nl(2,:,dtset%natom+3:dtset%natom+4,idir,ipert)=zero 1793 if (usepaw==1) d2ovl(2,:,dtset%natom+3:dtset%natom+4,idir,ipert)=zero 1794 end if 1795 1796 1797 !Symmetrize the strain perturbation contributions, as was needed for the stresses in a GS calculation 1798 !if (ipert==dtset%natom+3.or.ipert==dtset%natom+4)then 1799 if (nsym1>1) then 1800 ABI_MALLOC(work,(6,1,1)) 1801 ii=0 1802 do ipert1=dtset%natom+3,dtset%natom+4 1803 do idir1=1,3 1804 ii=ii+1 1805 work(ii,1,1)=d2nl(1,idir1,ipert1,idir,ipert) 1806 end do 1807 end do 1808 call stresssym(gprimd,nsym1,work(:,1,1),symrc1) 1809 ii=0 1810 do ipert1=dtset%natom+3,dtset%natom+4 1811 do idir1=1,3 1812 ii=ii+1 1813 d2nl(1,idir1,ipert1,idir,ipert)=work(ii,1,1) 1814 end do 1815 end do 1816 if (usepaw==1) then 1817 ii=0 1818 do ipert1=dtset%natom+3,dtset%natom+4 1819 do idir1=1,3 1820 ii=ii+1 1821 work(ii,1,1)=d2ovl(1,idir1,ipert1,idir,ipert) 1822 end do 1823 end do 1824 call stresssym(gprimd,nsym1,work(:,1,1),symrc1) 1825 ii=0 1826 do ipert1=dtset%natom+3,dtset%natom+4 1827 do idir1=1,3 1828 ii=ii+1 1829 d2ovl(1,idir1,ipert1,idir,ipert)=work(ii,1,1) 1830 end do 1831 end do 1832 ABI_FREE(work) 1833 end if 1834 end if 1835 !end if 1836 1837 !Must also symmetrize the electric field perturbation response ! 1838 !Note: d2ovl is not symetrized because it is zero for electric field perturbation 1839 if (has_ddk_file) then 1840 ABI_MALLOC(d2nl_elfd,(2,3)) 1841 ! There should not be any imaginary part, but stay general (for debugging) 1842 d2nl_elfd (:,:)=d2nl(:,:,dtset%natom+2,idir,ipert) 1843 do ii=1,3 1844 sumelfd(:)=zero 1845 do ia=1,nsym1 1846 do jj=1,3 1847 if(symrl1(ii,jj,ia)/=0)then 1848 if(ddkfil(jj)==0)then 1849 blkflg(ii,dtset%natom+2,idir,ipert)=0 1850 end if 1851 end if 1852 end do 1853 symfact(1)=dble(symrl1(ii,1,ia)) 1854 symfact(2)=dble(symrl1(ii,2,ia)) 1855 symfact(3)=dble(symrl1(ii,3,ia)) 1856 sumelfd(:)=sumelfd(:)+symfact(1)*d2nl_elfd(:,1) & 1857 & +symfact(2)*d2nl_elfd(:,2)+symfact(3)*d2nl_elfd(:,3) 1858 end do 1859 d2nl(:,ii,dtset%natom+2,idir,ipert)=sumelfd(:)/dble(nsym1) 1860 end do 1861 ABI_FREE(d2nl_elfd) 1862 end if 1863 1864 !Overlap: store the diagonal part of the matrix in the 1865 ! 2nd-order energy non-stationnary expression 1866 eovl1=zero;if (usepaw==1) eovl1=d2ovl(1,idir,ipert,idir,ipert) 1867 1868 ABI_FREE(bands_treated_now) 1869 1870 call destroy_mpi_enreg(mpi_enreg_seq) 1871 call timab(566,2,tsec) 1872 1873 DBG_EXIT("COLL") 1874 1875 end subroutine dfpt_nstpaw
ABINIT/dfpt_nstwf [ Functions ]
NAME
dfpt_nstwf
FUNCTION
This routine computes the non-local contribution to the 2DTE matrix elements, in the non-stationary formulation Only for norm-conserving pseudopotentials (no PAW)
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (XG,AR,MB,MVer,MT, MVeithen) 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_mem*mkmem*nsppol)=planewave coefficients of wavefunctions at k
cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
ddkfil(3)=unit numbers for the three possible ddk files for ipert1
equal to 0 if no dot file is available for this direction
dtset <type(dataset_type)>=all input variables for this dataset
eig_k(mband*nsppol)=GS eigenvalues at k (hartree)
eig1_k(2*nsppol*mband**2)=matrix of first-order eigenvalues (hartree)
gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
icg=shift to be applied on the location of data in the array cg
icg1=shift to be applied on the location of data in the array cg1
idir=direction of the current perturbation
ikpt=number of the k-point
ipert=type of the perturbation
isppol=1 for unpolarized, 2 for spin-polarized
istwf_k=parameter that describes the storage of wfs
kg_k(3,npw_k)=reduced planewave coordinates.
kg1_k(3,npw1_k)=reduced planewave coordinates at k+q, with RF k points
kpt(3)=reduced coordinates of k point
kpq(3)=reduced coordinates of k+q point
mkmem =number of k points treated by this node
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).
nband_k=number of bands at this k point for that spin polarization
npw_k=number of plane waves at this k point
npw1_k=number of plane waves at this k+q point
nsppol=1 for unpolarized, 2 for spin-polarized
occ_k(nband_k)=occupation number for each band (usually 2) for each k.
psps <type(pseudopotential_type)>=variables related to pseudopotentials
rmet(3,3)=real space metric (bohr**2)
ddks(3)<wfk_t>=struct info for for the three possible DDK files for ipert1
wtk_k=weight assigned to the k point.
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
OUTPUT
d2bbb_k(2,3,mband,mband*prtbbb)=band by band decomposition of the second order derivatives, for the present k point, and perturbation idir, ipert d2nl_k(2,3,mpert)=non-local contributions to non-stationary 2DTE, for the present k point, and perturbation idir, ipert
TODO
XG 20141103 The localization tensor cannot be defined in the metallic case. It should not be computed.
SOURCE
1943 subroutine dfpt_nstwf(cg,cg1,ddkfil,dtset,d2bbb_k,d2nl_k,eig_k,eig1_k,gs_hamkq,& 1944 & icg,icg1,idir,ikpt,ipert,isppol,istwf_k,kg_k,kg1_k,kpt,kpq,& 1945 & mband_mem_rbz,mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,nband_k,npw_k,npw1_k,nsppol,& 1946 & occ_k,psps,rmet,ddks,wtk_k,ylm,ylm1) 1947 1948 !Arguments ------------------------------------ 1949 !scalars 1950 integer,intent(in) :: icg,icg1,idir,ikpt,ipert,isppol,istwf_k 1951 integer,intent(in) :: mkmem,mk1mem,mpert,mpw,mpw1,nsppol 1952 integer,intent(in) :: mband_mem_rbz 1953 integer,intent(inout) :: nband_k,npw1_k,npw_k 1954 real(dp),intent(in) :: wtk_k 1955 type(MPI_type),intent(in) :: mpi_enreg 1956 type(dataset_type),intent(in) :: dtset 1957 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 1958 type(pseudopotential_type),intent(in) :: psps 1959 !arrays 1960 integer,intent(in) :: ddkfil(3),kg1_k(3,npw1_k) 1961 integer,intent(in) :: kg_k(3,npw_k) 1962 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*mband_mem_rbz*mkmem*nsppol) 1963 real(dp),intent(in) :: cg1(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem*nsppol) 1964 real(dp),intent(in) :: eig_k(dtset%mband*nsppol),kpt(3),kpq(3),occ_k(nband_k),rmet(3,3) 1965 real(dp),intent(in) :: ylm(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 1966 real(dp),intent(in) :: ylm1(npw1_k,psps%mpsang*psps%mpsang*psps%useylm) 1967 real(dp),intent(inout) :: eig1_k(2*nsppol*dtset%mband**2) 1968 real(dp),intent(out) :: d2bbb_k(2,3,dtset%mband,dtset%mband*dtset%prtbbb) 1969 real(dp),intent(inout) :: d2nl_k(2,3,mpert) 1970 type(wfk_t),intent(inout) :: ddks(3) 1971 1972 !Local variables------------------------------- 1973 !scalars 1974 integer :: berryopt,dimffnl,dimffnl1,dimph3d 1975 integer :: iband,ider,idir1,ipert1,ipw,jband,nband_kocc,nkpg,nkpg1 1976 integer :: ierr, iband_me, jband_me 1977 integer :: npw_disk,nsp,optlocal,optnl,opt_gvnlx1,sij_opt,tim_getgh1c,usevnl 1978 integer :: nddk_needed, startband, endband 1979 logical :: ddk 1980 real(dp) :: aa,dot1i,dot1r,dot2i,dot2r,dot_ndiagi,dot_ndiagr,doti,dotr,lambda 1981 character(len=500) :: msg 1982 type(rf_hamiltonian_type) :: rf_hamkq 1983 !arrays 1984 integer :: ik_ddks(3) 1985 integer :: band_procs(nband_k) 1986 logical :: distrb_cycle(nband_k) 1987 real(dp) :: dum_grad_berry(1,1),dum_gvnlx1(1,1),dum_gs1(1,1),dum_ylmgr(1,3,1),tsec(2) 1988 real(dp),allocatable :: cg_k(:,:),cwave0(:,:),cwavef(:,:),cwavef_da(:,:) 1989 real(dp),allocatable :: cwaveddk(:,:,:) 1990 real(dp),allocatable :: cg_ddk(:,:,:) !2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem*nsppol,3) == 1991 1992 real(dp),allocatable :: cwavef_db(:,:),dkinpw(:),eig2_k(:),ffnl1(:,:,:,:),ffnlk(:,:,:,:) 1993 real(dp),allocatable :: eig2_ddk(:,:) 1994 real(dp),allocatable :: gvnlx1(:,:),kinpw1(:),kpg1_k(:,:),kpg_k(:,:),ph3d(:,:,:) 1995 type(pawcprj_type),allocatable :: dum_cwaveprj(:,:) 1996 1997 ! ********************************************************************* 1998 1999 DBG_ENTER("COLL") 2000 2001 !Not valid for PAW 2002 if (psps%usepaw==1) then 2003 msg=' This routine cannot be used for PAW (use pawnst3 instead) !' 2004 ABI_BUG(msg) 2005 end if 2006 2007 !Keep track of total time spent in dfpt_nstwf 2008 call timab(112,1,tsec) 2009 tim_getgh1c=2 2010 2011 !Miscelaneous inits 2012 ABI_MALLOC(dum_cwaveprj,(0,0)) 2013 ddk=(ipert==dtset%natom+1.or.ipert==dtset%natom+10.or.ipert==dtset%natom+11) 2014 2015 ! filter for bands on this cpu for cg cg1 etc. 2016 distrb_cycle = (mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol) /= mpi_enreg%me_kpt) 2017 call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,nband_k,& 2018 & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band) 2019 2020 !Additional allocations 2021 if (.not.ddk) then 2022 ABI_MALLOC(dkinpw,(npw_k)) 2023 ABI_MALLOC(kinpw1,(npw1_k)) 2024 kinpw1=zero;dkinpw=zero 2025 else 2026 ABI_MALLOC(dkinpw,(0)) 2027 ABI_MALLOC(kinpw1,(0)) 2028 end if 2029 ABI_MALLOC(gvnlx1,(2,npw1_k*dtset%nspinor)) 2030 ABI_MALLOC(eig2_k,(2*nsppol*dtset%mband**2)) 2031 ABI_MALLOC(cwave0,(2,npw_k*dtset%nspinor)) 2032 ABI_MALLOC(cwavef,(2,npw1_k*dtset%nspinor)) 2033 2034 nddk_needed = 0 2035 do idir1=1,3 2036 if (ddkfil(idir1)/=0) nddk_needed = nddk_needed+1 2037 end do 2038 if (nddk_needed > 0) then 2039 ABI_MALLOC(cwaveddk,(2,npw1_k*dtset%nspinor,3)) 2040 !TODO: for the moment avoid indirect indexing of the ddk directions in case not all are present. Here all are allocated and read in 2041 ABI_MALLOC(cg_ddk,(2,mpw1*dtset%nspinor*mband_mem_rbz,3)) 2042 cg_ddk = zero ! not all may be initialized below if only certain ddk directions are provided 2043 2044 ABI_MALLOC(eig2_ddk,(2*dtset%mband**2,3)) 2045 eig2_ddk = zero 2046 end if 2047 2048 !Compute (k+G) vectors 2049 nkpg=0;if (.not.ddk) nkpg=3*gs_hamkq%nloalg(3) 2050 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 2051 if (nkpg>0) then 2052 call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k) 2053 end if 2054 2055 !Compute (k+q+G) vectors 2056 nkpg1=0;if (.not.ddk) nkpg1=3*gs_hamkq%nloalg(3) 2057 ABI_MALLOC(kpg1_k,(npw1_k,nkpg1)) 2058 if (nkpg1>0) then 2059 call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k) 2060 end if 2061 2062 !Compute nonlocal form factors ffnl at (k+G) 2063 dimffnl=0;if (.not.ddk) dimffnl=1 2064 ABI_MALLOC(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 2065 if (.not.ddk) then 2066 ider=0 2067 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamkq%gmet,& 2068 & gs_hamkq%gprimd,ider,ider,psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,& 2069 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,& 2070 & psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,dum_ylmgr) 2071 end if 2072 2073 !Compute nonlocal form factors ffnl1 at (k+q+G) 2074 dimffnl1=0;if (.not.ddk) dimffnl1=1 2075 ABI_MALLOC(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat)) 2076 if (.not.ddk) then 2077 ider=0 2078 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,& 2079 & gs_hamkq%gprimd,ider,ider,psps%indlmn,kg1_k,kpg1_k,kpq,& 2080 & psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,npw1_k,psps%ntypat,& 2081 & psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1,dum_ylmgr) 2082 end if 2083 2084 !Load k-dependent part in the Hamiltonian datastructure 2085 call gs_hamkq%load_k(kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,& 2086 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnlk) 2087 2088 !Load k+q-dependent part in the Hamiltonian datastructure 2089 ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead) 2090 dimph3d=0;if (.not.ddk) dimph3d=gs_hamkq%matblk 2091 ABI_MALLOC(ph3d,(2,npw1_k,dimph3d)) 2092 call gs_hamkq%load_kprime(kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,& 2093 & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,& 2094 & ph3d_kp=ph3d,compute_ph3d=(.not.ddk)) 2095 2096 !Load k-dependent part in the 1st-order Hamiltonian datastructure 2097 call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dkinpw) 2098 2099 !Take care of the npw and kg records 2100 !NOTE : one should be able to modify the rwwf routine to take care 2101 !of the band parallelism, which is not the case yet ... 2102 ik_ddks = 0 2103 do idir1=1,3 2104 if (ddkfil(idir1)/=0)then 2105 ! Read npw record 2106 nsp=dtset%nspinor 2107 ik_ddks(idir1) = ddks(idir1)%findk(kpt) 2108 ABI_CHECK(ik_ddks(idir1) /= -1, "Cannot find kpt") 2109 npw_disk = ddks(idir1)%hdr%npwarr(ik_ddks(idir1)) 2110 if (npw_k /= npw_disk) then 2111 write(unit=msg,fmt='(a,i3,a,i5,a,i3,a,a,i5,a,a,i5)')& 2112 & 'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',idir,ch10,& 2113 & 'the number of plane waves in the ddk file is equal to', npw_disk,ch10,& 2114 & 'while it should be ',npw_k 2115 ABI_BUG(msg) 2116 end if 2117 2118 ! NB: this will fail if the bands are not contiguous. 2119 startband = nband_k 2120 endband = 1 2121 do iband=1,nband_k 2122 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) == mpi_enreg%me_kpt) then 2123 if (iband < startband) startband = iband 2124 if (iband > endband) endband = iband 2125 end if 2126 end do 2127 ! NB: eig_k is band distributed in call to read_band_block, though array has full size, 2128 ! only certain columns for my iband are filled, then used below 2129 call ddks(idir1)%read_band_block((/startband,endband/),ik_ddks(idir1),isppol,xmpio_collective, & 2130 & cg_k=cg_ddk(:,:,idir1), eig_k=eig2_ddk(:,idir1)) 2131 end if ! ddk file is already present 2132 end do ! idir1 2133 2134 if (ipert==dtset%natom+1) then 2135 nband_kocc = 0 2136 do iband = 1,nband_k 2137 if (abs(occ_k(iband)) > tol8) nband_kocc = nband_kocc + 1 2138 nband_kocc = max (nband_kocc, 1) 2139 end do 2140 end if 2141 2142 if(dtset%prtbbb==1)then 2143 ABI_MALLOC(cwavef_da,(2,npw1_k*dtset%nspinor)) 2144 ABI_MALLOC(cwavef_db,(2,npw1_k*dtset%nspinor)) 2145 ABI_MALLOC(cg_k,(2,npw_k*dtset%nspinor*mband_mem_rbz)) 2146 if ((ipert == dtset%natom + 1).or.(ipert <= dtset%natom).or. & 2147 & (ipert == dtset%natom + 2).or.(ipert == dtset%natom + 5)) then 2148 cg_k(:,:) = cg(:,1+icg:icg+mband_mem_rbz*npw_k*dtset%nspinor) 2149 end if 2150 d2bbb_k(:,:,:,:) = zero 2151 end if 2152 2153 !Loop over ALL bands 2154 iband_me = 0 2155 do iband=1,nband_k 2156 2157 ! if band is mine, retrieve it and then broadcast it 2158 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) == mpi_enreg%me_kpt) then 2159 iband_me = iband_me + 1 2160 2161 ! Read ground-state wavefunction for iband 2162 if (dtset%prtbbb==0 .or. ipert==dtset%natom+2) then 2163 cwave0(:,:)=cg(:,1+(iband_me-1)*npw_k*dtset%nspinor+icg:iband_me*npw_k*dtset%nspinor+icg) 2164 else ! prtbbb==1 and ipert<=natom , already in cg_k 2165 cwave0(:,:)=cg_k(:,1+(iband_me-1)*npw_k*dtset%nspinor:iband_me*npw_k*dtset%nspinor) 2166 end if 2167 2168 ! Get first-order wavefunctions for iband 2169 cwavef(:,:)=cg1(:,1+(iband_me-1)*npw1_k*dtset%nspinor+icg1:iband_me*npw1_k*dtset%nspinor+icg1) 2170 ! Get ddk wavefunctions for iband 2171 if(nddk_needed > 0) then 2172 cwaveddk(:,:,:)=cg_ddk(:,1+(iband_me-1)*npw1_k*dtset%nspinor:iband_me*npw1_k*dtset%nspinor,:) 2173 end if 2174 end if 2175 call xmpi_bcast(cwave0, band_procs(iband), mpi_enreg%comm_band, ierr) 2176 call xmpi_bcast(cwavef, band_procs(iband), mpi_enreg%comm_band, ierr) 2177 if(nddk_needed > 0) then 2178 call xmpi_bcast(cwaveddk, band_procs(iband), mpi_enreg%comm_band, ierr) 2179 end if 2180 2181 ! In case non ddk perturbation 2182 if (ipert /= dtset%natom + 1) then 2183 2184 do ipert1=1,mpert 2185 2186 if( ipert1<=dtset%natom .or. ipert1==dtset%natom+2 )then 2187 2188 ! Initialize data for NL 1st-order hamiltonian 2189 call init_rf_hamiltonian(1,gs_hamkq,ipert1,rf_hamkq) 2190 2191 if (((ipert <= dtset%natom).or.(ipert == dtset%natom + 2)) & 2192 & .and.(ipert1 == dtset%natom+2).and. dtset%prtbbb==1) then 2193 call gaugetransfo(cg_k,cwavef,cwavef_db,mpi_enreg%comm_band,distrb_cycle,eig_k,eig1_k,iband,nband_k, & 2194 & dtset%mband,mband_mem_rbz,npw_k,npw1_k,dtset%nspinor,nsppol,mpi_enreg%nproc_band,occ_k) 2195 cwavef(:,:) = cwavef_db(:,:) 2196 end if 2197 2198 ! Define the direction along which to move the atom : 2199 ! the polarisation (ipert1,idir1) is refered as j1. 2200 do idir1=1,3 2201 if (ipert1<=dtset%natom.or.(ipert1==dtset%natom+2.and.ddkfil(idir1)/=0)) then 2202 2203 ! Get |Vnon-locj^(1)|u0> : 2204 ! First-order non-local, applied to zero-order wavefunction 2205 ! This routine gives MINUS the non-local contribution 2206 2207 ! ==== Atomic displ. perturbation 2208 if( ipert1<=dtset%natom )then 2209 lambda=eig_k((isppol-1)*nband_k+iband) 2210 berryopt=1;optlocal=0;optnl=1;usevnl=0;opt_gvnlx1=0;sij_opt=0 2211 call getgh1c(berryopt,cwave0,dum_cwaveprj,gvnlx1,dum_grad_berry,& 2212 & dum_gs1,gs_hamkq,dum_gvnlx1,idir1,ipert1,lambda,mpi_enreg,optlocal,& 2213 & optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 2214 2215 ! ==== Electric field perturbation 2216 else if( ipert1==dtset%natom+2 )then 2217 ! TODO: Several tests fail here ifdef HAVE_MPI_IO_DEFAULT 2218 ! The problem is somehow related to the use of MPI-IO file views!. 2219 !TODO MJV: Check if it works with HAVE_MPI_IO_DEFAULT now. 2220 2221 gvnlx1 = cwaveddk(:,:,idir1) 2222 eig2_k(1+(iband-1)*2*nband_k:iband*2*nband_k) = eig2_ddk(1+(iband-1)*2*nband_k:iband*2*nband_k,idir1) 2223 2224 !write(777,*)"eig2_k, gvnlx1 for band: ",iband,", ikpt: ",ikpt 2225 !do ii=1,2*nband_k 2226 ! write(777,*)eig2_k(ii+(iband-1)) 2227 !end do 2228 !write(777,*)gvnlx1 2229 2230 ! In case of band-by-band, 2231 ! construct the first-order wavefunctions in the diagonal gauge 2232 if (((ipert <= dtset%natom).or.(ipert == dtset%natom + 2)).and.(dtset%prtbbb==1)) then 2233 call gaugetransfo(cg_k,gvnlx1,cwavef_da,mpi_enreg%comm_band,distrb_cycle,eig_k,eig2_k,iband,nband_k, & 2234 & dtset%mband,mband_mem_rbz,npw_k,npw1_k,dtset%nspinor,nsppol,mpi_enreg%nproc_band,occ_k) 2235 gvnlx1(:,:) = cwavef_da(:,:) 2236 end if 2237 ! Multiplication by -i 2238 do ipw=1,npw1_k*dtset%nspinor 2239 aa=gvnlx1(1,ipw) 2240 gvnlx1(1,ipw)=gvnlx1(2,ipw) 2241 gvnlx1(2,ipw)=-aa 2242 end do 2243 end if 2244 2245 ! at this stage if iband is not mine I can cycle 2246 if (mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) then 2247 cycle 2248 end if 2249 2250 ! MVeithen 021212 : 2251 ! 1) Case ipert1 = natom + 2 and ipert = natom + 2: 2252 ! the second derivative of the energy with respect to an electric 2253 ! field is computed from Eq. (38) of X. Gonze, PRB 55 ,10355 (1997) [[cite:Gonze1997a]]. 2254 ! The evaluation of this formula needs the operator $i \frac{d}{dk}. 2255 ! 2) Case ipert1 = natom + 2 and ipert < natom: 2256 ! the computation of the Born effective charge tensor uses 2257 ! the operator $-i \frac{d}{dk}. 2258 if (ipert==dtset%natom+2) gvnlx1(:,:) = -gvnlx1(:,:) 2259 2260 ! <G|Vnl1|Cnk> is contained in gvnlx1 2261 ! construct the matrix element (<uj2|vj1|u0>)complex conjug and add it to the 2nd-order matrix 2262 call dotprod_g(dotr,doti,istwf_k,npw1_k*dtset%nspinor,2,cwavef,gvnlx1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2263 d2nl_k(1,idir1,ipert1)=d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*two*dotr 2264 d2nl_k(2,idir1,ipert1)=d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*two*doti 2265 2266 ! Band by band decomposition of the Born effective charges 2267 ! calculated from a phonon perturbation 2268 if(dtset%prtbbb==1) then ! .and. mpi_enreg%proc_distrb(ikpt,iband,isppol) == mpi_enreg%me_kpt)then 2269 d2bbb_k(1,idir1,iband,iband) = wtk_k*occ_k(iband)*two*dotr 2270 d2bbb_k(2,idir1,iband,iband) = -one*wtk_k*occ_k(iband)*two*doti 2271 end if 2272 2273 end if 2274 end do ! idir 2275 2276 call rf_hamkq%free() 2277 end if ! ipert1<=dtset%natom .or. ipert1==dtset%natom+2 2278 end do ! ipert1 2279 end if ! ipert /= natom +1 2280 2281 ! Compute the localization tensor 2282 2283 if (ipert==dtset%natom+1) then 2284 2285 ipert1=dtset%natom+1 2286 if(dtset%prtbbb==1)then 2287 call gaugetransfo(cg_k,cwavef,cwavef_db,mpi_enreg%comm_band,distrb_cycle,eig_k,eig1_k,iband,nband_k, & 2288 & dtset%mband,mband_mem_rbz,npw_k,npw1_k,dtset%nspinor,nsppol,mpi_enreg%nproc_band,occ_k) 2289 cwavef(:,:) = cwavef_db(:,:) 2290 end if 2291 2292 do idir1 = 1,3 2293 eig2_k(:) = zero 2294 gvnlx1(:,:) = zero 2295 if (idir == idir1) then 2296 gvnlx1(:,:) = cwavef(:,:) 2297 eig2_k(:) = eig1_k(:) 2298 else 2299 if (ddkfil(idir1) /= 0) then 2300 gvnlx1 = cwaveddk(:,:,idir1) 2301 eig2_k(1+(iband-1)*2*nband_k:iband*2*nband_k) = eig2_ddk(1+(iband-1)*2*nband_k:iband*2*nband_k,idir1) 2302 2303 !write(778,*)"eig2_k, gvnlx1 for band: ",iband,", ikpt: ",ikpt 2304 !do ii=1,2*nband_k 2305 ! write(778,*)eig2_k(ii+(iband-1)) 2306 !end do 2307 !write(778,*)gvnlx1 2308 2309 if(dtset%prtbbb==1)then 2310 call gaugetransfo(cg_k,gvnlx1,cwavef_da,mpi_enreg%comm_band,distrb_cycle,eig_k,eig2_k,iband,nband_k, & 2311 & dtset%mband,mband_mem_rbz,npw_k,npw1_k,dtset%nspinor,nsppol,mpi_enreg%nproc_band,occ_k) 2312 2313 gvnlx1(:,:) = cwavef_da(:,:) 2314 end if 2315 2316 end if !ddkfil(idir1) 2317 end if !idir == idir1 2318 2319 ! <G|du/dqa> is contained in gvnlx1 and <G|du/dqb> in cwavef 2320 ! construct the matrix elements <du/dqa|du/dqb> -> dot 2321 ! <u|du/dqa> -> dot1 2322 ! <du/dqb|u> -> dot2 2323 ! and add them to the 2nd-order matrix 2324 2325 call dotprod_g(dotr,doti,istwf_k,npw1_k*dtset%nspinor,2,gvnlx1,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2326 d2nl_k(1,idir1,ipert1)=d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*dotr/(nband_kocc*two) 2327 d2nl_k(2,idir1,ipert1)=d2nl_k(2,idir1,ipert1)+wtk_k*occ_k(iband)*doti/(nband_kocc*two) 2328 2329 2330 ! XG 020216 : Marek, could you check the next forty lines 2331 ! In the parallel gauge, dot1 and dot2 vanishes 2332 if(dtset%prtbbb==1)then 2333 if (mpi_enreg%proc_distrb(ikpt,iband,isppol) == mpi_enreg%me_kpt) then 2334 d2bbb_k(1,idir1,iband,iband)=d2bbb_k(1,idir1,iband,iband)+dotr 2335 d2bbb_k(2,idir1,iband,iband)=d2bbb_k(2,idir1,iband,iband)+doti 2336 end if 2337 dot_ndiagr=zero ; dot_ndiagi=zero 2338 jband_me = 0 2339 do jband = 1,nband_k !compute dot1 and dot2 2340 if (mpi_enreg%proc_distrb(ikpt,jband,isppol) /= mpi_enreg%me_kpt) then 2341 cycle 2342 end if 2343 jband_me = jband_me + 1 2344 2345 if (abs(occ_k(jband)) > tol8) then 2346 dot1r=zero ; dot1i=zero 2347 dot2r=zero ; dot2i=zero 2348 cwave0(:,:)=cg_k(:,1+(jband_me-1)*npw_k*dtset%nspinor:jband_me*npw_k*dtset%nspinor) 2349 2350 call dotprod_g(dot1r,dot1i,istwf_k,npw1_k*dtset%nspinor,2,cwave0,gvnlx1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2351 call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*dtset%nspinor,2,cwavef,cwave0,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2352 2353 dot_ndiagr = dot_ndiagr + dot1r*dot2r - dot1i*dot2i 2354 dot_ndiagi = dot_ndiagi + dot1r*dot2i + dot1i*dot2r 2355 ! this should fill all of the iband but only the local cpu jband indices 2356 d2bbb_k(1,idir1,iband,jband) = d2bbb_k(1,idir1,iband,jband) - & 2357 & (dot1r*dot2r - dot1i*dot2i) 2358 d2bbb_k(2,idir1,iband,jband) = d2bbb_k(2,idir1,iband,jband) - & 2359 & (dot1r*dot2i + dot1i*dot2r) 2360 end if ! occ_k 2361 end do !jband 2362 2363 d2bbb_k(:,idir1,iband,:)=d2bbb_k(:,idir1,iband,:)*wtk_k*occ_k(iband)*half 2364 d2nl_k(1,idir1,ipert1)= & 2365 & d2nl_k(1,idir1,ipert1)-wtk_k*occ_k(iband)*dot_ndiagr/(nband_kocc*two) 2366 d2nl_k(2,idir1,ipert1)=& 2367 & d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*dot_ndiagi/(nband_kocc*two) 2368 end if ! prtbbb==1 2369 2370 end do ! idir1 2371 end if ! Compute localization tensor, ipert=natom+1 2372 2373 end do ! End loop over iband 2374 2375 ! if(dtset%prtbbb==1)then 2376 ! ! complete over jband index 2377 ! call xmpi_sum(d2bbb_k, mpi_enreg%comm_band, ierr) 2378 ! end if 2379 2380 !Final deallocations 2381 ABI_FREE(cwave0) 2382 ABI_FREE(cwavef) 2383 ABI_FREE(eig2_k) 2384 ABI_FREE(gvnlx1) 2385 ABI_FREE(ffnlk) 2386 ABI_FREE(ffnl1) 2387 ABI_FREE(dkinpw) 2388 ABI_FREE(kinpw1) 2389 ABI_FREE(kpg_k) 2390 ABI_FREE(kpg1_k) 2391 ABI_FREE(ph3d) 2392 ABI_FREE(dum_cwaveprj) 2393 if(dtset%prtbbb==1) then 2394 ABI_FREE(cg_k) 2395 ABI_FREE(cwavef_da) 2396 ABI_FREE(cwavef_db) 2397 end if 2398 if (nddk_needed > 0) then 2399 ABI_FREE(cwaveddk) 2400 ABI_FREE(cg_ddk) 2401 ABI_FREE(eig2_ddk) 2402 end if 2403 2404 call timab(112,2,tsec) 2405 2406 DBG_EXIT("COLL") 2407 2408 end subroutine dfpt_nstwf
ABINIT/gaugetransfo [ Functions ]
NAME
gaugetransfo
FUNCTION
This routine allows the passage from the parallel-transport gauge to the diagonal gauge for the first-order wavefunctions
INPUTS
cg_k(2,mpw*nspinor*mband_mem*nsppol)=planewave coefficients of wavefunctions
for a particular k point.
cwavef(2,npw1_k*nspinor)=first order wavefunction for a particular k point
in the parallel gauge
comm=mpi communicator for bands
distrb_cycle=array of logical flags to skip certain bands in parallelization scheme
eig_k(mband*nsppol)=GS eigenvalues at k (hartree)
eig1_k(2*nsppol*mband**2)=matrix of first-order eigenvalues (hartree)
iband=band index of the 1WF for which the transformation has to be applied
mband=maximum number of bands
mband_mem_rbz=maximum number of bands on this cpu
nband_k=number of bands for this k point
npw_k=maximum dimensioned size of npw or wfs at k
npw1_k=number of plane waves at this k+q point
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
occ_k(nband_k)=occupation number for each band (usually 2) for each k
OUTPUT
cwavef_d(2,npw1_k*nspinor)=first order wavefunction for a particular k point
in the diagonal gauge
SOURCE
2444 subroutine gaugetransfo(cg_k,cwavef,cwavef_d,comm,distrb_cycle,eig_k,eig1_k,iband,nband_k, & 2445 & mband,mband_mem_rbz,npw_k,npw1_k,nspinor,nsppol,nproc_band,occ_k) 2446 2447 !Arguments ------------------------------------ 2448 !scalars 2449 integer,intent(in) :: iband,mband,mband_mem_rbz,nband_k,npw1_k,npw_k,nspinor,nsppol 2450 integer,intent(in) :: comm, nproc_band 2451 !arrays 2452 logical, intent(in) :: distrb_cycle(nband_k) 2453 real(dp),intent(in) :: cg_k(2,npw_k*nspinor*mband_mem_rbz),cwavef(2,npw1_k*nspinor) 2454 real(dp),intent(in) :: eig1_k(2*nsppol*mband**2),eig_k(mband*nsppol) 2455 real(dp),intent(in) :: occ_k(nband_k) 2456 real(dp),intent(out) :: cwavef_d(2,npw1_k*nspinor) 2457 2458 !Local variables------------------------------- 2459 !tolerance for non degenerated levels 2460 !scalars 2461 integer :: ierr, jband,jband_me 2462 real(dp),parameter :: etol=1.0d-3 2463 !arrays 2464 real(dp) :: cwave0(2,npw1_k*nspinor),eig1(2) 2465 2466 ! ********************************************************************* 2467 2468 cwavef_d(:,:) = cwavef(:,:) 2469 2470 jband_me = 0 2471 do jband = 1,nband_k !loop over bands 2472 if (distrb_cycle(jband)) cycle 2473 jband_me = jband_me + 1 2474 2475 if ((abs(eig_k(iband)-eig_k(jband)) > etol).and.(abs(occ_k(jband)) > tol8 )) then 2476 2477 cwave0(:,:) = cg_k(:,1+(jband_me-1)*npw_k*nspinor:jband_me*npw_k*nspinor) 2478 2479 eig1(1) = eig1_k(2*jband-1+(iband-1)*2*nband_k) 2480 eig1(2) = eig1_k(2*jband +(iband-1)*2*nband_k) 2481 2482 cwavef_d(1,:)=cwavef_d(1,:) & 2483 & - (eig1(1)*cwave0(1,:)-eig1(2)*cwave0(2,:))/(eig_k(jband)-eig_k(iband)) 2484 cwavef_d(2,:)=cwavef_d(2,:) & 2485 & - (eig1(1)*cwave0(2,:)+eig1(2)*cwave0(1,:))/(eig_k(jband)-eig_k(iband)) 2486 2487 end if 2488 2489 end do !loop over bands 2490 call xmpi_sum(cwavef_d, comm, ierr) 2491 ! here we have summed the cwavef N times (N-1 too many), but the correction is completed over bands 2492 cwavef_d = cwavef_d - dble(nproc_band-1)*cwavef 2493 2494 end subroutine gaugetransfo
ABINIT/m_dfpt_nstwf [ Modules ]
NAME
m_dfpt_nstwf
FUNCTION
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group () This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_dfpt_nstwf 22 23 use defs_basis 24 use m_xmpi 25 use m_errors 26 use m_abicore 27 use m_wfk 28 use m_hamiltonian 29 use m_cgtools 30 use m_nctk 31 use m_dtset 32 use m_dtfil 33 34 use defs_datatypes, only : pseudopotential_type 35 use defs_abitypes, only : MPI_type 36 use m_time, only : timab 37 use m_io_tools, only : file_exists 38 use m_fourier_interpol, only : transgrid 39 use m_geometry, only : stresssym 40 use m_dynmat, only : dfpt_sygra 41 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq, proc_distrb_cycle, proc_distrb_band, proc_distrb_nband 42 use m_hdr, only : hdr_skip 43 use m_occ, only : occeig 44 use m_pawang, only : pawang_type 45 use m_pawrad, only : pawrad_type 46 use m_pawtab, only : pawtab_type 47 use m_paw_an, only : paw_an_type, paw_an_reset_flags 48 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags 49 use m_pawfgrtab,only : pawfgrtab_type 50 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, pawrhoij_nullify, & 51 pawrhoij_init_unpacked, pawrhoij_mpisum_unpacked, pawrhoij_inquire_dim 52 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get, pawcprj_copy 53 use m_pawdij, only : pawdijfr 54 use m_pawfgr, only : pawfgr_type 55 use m_paw_mkrho,only : pawmkrho 56 use m_paw_nhat, only : pawnhatfr 57 use m_paw_dfpt, only : pawdfptenergy 58 use m_kg, only : mkkin, kpgstr, mkkpg 59 use m_fft, only : fftpac 60 use m_spacepar, only : hartrestr, symrhg 61 use m_initylmg, only : initylmg 62 use m_mkffnl, only : mkffnl 63 use m_getgh1c, only : getgh1c, getdc1 64 use m_dfpt_mkrho, only : dfpt_accrho 65 use m_atm2fft, only : dfpt_atm2fft 66 use m_mkcore, only : dfpt_mkcore 67 use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll 68 use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr 69 use m_mklocl, only : dfpt_vlocal, vlocalstr 70 use m_cgprj, only : getcprj 71 72 implicit none 73 74 private