TABLE OF CONTENTS
ABINIT/m_dfptlw_wf [ Modules ]
NAME
m_dfptlw_wf
FUNCTION
FIXME: add description.
COPYRIGHT
Copyright (C) 2022-2024 ABINIT group (MR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
PARENTS
CHILDREN
SOURCE
23 #if defined HAVE_CONFIG_H 24 #include "config.h" 25 #endif 26 27 #include "abi_common.h" 28 29 module m_dfptlw_wf 30 31 use defs_basis 32 use defs_abitypes 33 use defs_datatypes 34 use m_dtset 35 use m_errors 36 use m_profiling_abi 37 use m_hamiltonian 38 use m_cgtools 39 use m_pawcprj 40 use m_pawfgr 41 use m_wfk 42 use m_xmpi 43 use m_getgh1c 44 use m_mklocl 45 46 use m_fstrings, only : itoa, sjoin 47 use m_io_tools, only : file_exists 48 use m_time, only : cwtime 49 use m_kg, only : mkkpg 50 51 implicit none 52 53 public :: dfpt_1wf 54 55 private 56 57 ! ************************************************************************* 58 59 contains
ABINIT/m_dfptlw_wf/dfpt_1wf [ Functions ]
NAME
dfpt_1wf
FUNCTION
Compute the spin, band and kpt resolved contributions to the spatial-dispersion third-order energy derivatives that depend on first-order response functions.
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k cplex: if 1, several magnitudes are REAL, if 2, COMPLEX ddk_f = wf files d2_dkdk_f = wf files dimffnl= third dimension of ffnl_k dtset <type(dataset_type)>=all input variables for this dataset eig1_k(2*nband_k**2)=1st-order eigenvalues at k for i1pert,i1dir eig2_k(2*nband_k**2)=1st-order eigenvalues at k for i2pert,i2dir ffnl_k(dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives for this k point gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k cg1 = first derivative of cg with respect the perturbation i1pert cg2 = first derivative of cg with respect the perturbation i2pert icg=shift to be applied on the location of data in the array cg i1dir,i2dir,i3dir=directions of the corresponding perturbations i1pert,i2pert,i3pert = type of perturbation that has to be computed ikpt=number of the k-point 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. kpt(3)=reduced coordinates of k point mkmem =number of k points treated by this node mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw or wfs at k natom= number of atoms in the unit cell natpert=number of atomic displacement perturbations nband_k=number of bands at this k point for that spin polarization nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other npw_k=number of plane waves at this k point nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nylmgr=second dimension of ylmgr_k occ_k(nband_k)=occupation number for each band (usually 2) for each k. pawfgr <type(pawfgr_type)>=fine grid parameters and related data psps <type(pseudopotential_type)>=variables related to pseudopotentials rmet(3,3)=real space metric (bohr**2) rprimd(3,3) = dimensional primitive translations (bohr) samepert= .true. if i1pert=i2pert and i1dir=i2dir useylmgr= if 1 use the derivative of spherical harmonics vpsp1_i1pertdq(cplex*nfft,nspden,n1dq)= local potential of first-order gradient Hamiltonian for i1pert vpsp1_i2pertdq(cplex*nfft,nspden,n2dq)= local potential of first-order gradient Hamiltonian for i2pert wtk_k=weight assigned to the k point. ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics for the k point ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics for the k point
OUTPUT
d3etot_t(1-5)_k= stationary 1wf contributions to the third-order energy derivatives for kpt
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
134 subroutine dfpt_1wf(cg,cg1,cg2,cplex,ddk_f,d2_dkdk_f,d2_dkdk_f2,& 135 & d3etot_t1_k,d3etot_t2_k,d3etot_t3_k,& 136 & d3etot_t4_k,d3etot_t5_k,dimffnl,dtset,eig1_k,eig2_k,ffnl_k,gs_hamkq,icg,& 137 & i1dir,i2dir,i3dir,i1pert,i2pert,ikpt,isppol,istwf_k,& 138 & kg_k,kpt,mkmem,mpi_enreg,mpw,natom,nband_k,& 139 & n1dq,n2dq,nfft,ngfft,npw_k,nspden,nsppol,nylmgr,occ_k,& 140 & pawfgr,psps,rmet,rprimd,samepert,useylmgr,& 141 & vpsp1_i1pertdq,vpsp1_i2pertdq,& 142 & wtk_k,ylm_k,ylmgr_k) 143 144 use defs_basis 145 146 implicit none 147 148 !Arguments ------------------------------------ 149 !scalars 150 integer,intent(in) :: cplex,dimffnl,i1dir,i1pert,i2dir,i2pert,i3dir 151 integer,intent(in) :: icg,ikpt,isppol,istwf_k 152 integer,intent(in) :: mkmem,mpw,natom,nband_k,n1dq,n2dq,nfft 153 integer,intent(in) :: npw_k,nspden,nsppol,nylmgr 154 integer,intent(in) :: useylmgr 155 real(dp),intent(in) :: wtk_k 156 logical,intent(in) :: samepert 157 type(dataset_type),intent(in) :: dtset 158 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 159 type(MPI_type),intent(in) :: mpi_enreg 160 type(pseudopotential_type),intent(in) :: psps 161 type(wfk_t),intent(inout) :: ddk_f,d2_dkdk_f, d2_dkdk_f2 162 type(pawfgr_type),intent(in) :: pawfgr 163 164 !arrays 165 integer,intent(in) :: kg_k(3,npw_k),ngfft(18) 166 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 167 real(dp),intent(in) :: cg1(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 168 real(dp),intent(in) :: cg2(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 169 real(dp),intent(out) :: d3etot_t1_k(2) 170 real(dp),intent(out) :: d3etot_t2_k(2) 171 real(dp),intent(out) :: d3etot_t3_k(2) 172 real(dp),intent(out) :: d3etot_t4_k(2,n2dq) 173 real(dp),intent(out) :: d3etot_t5_k(2,n1dq) 174 real(dp),intent(in) :: eig1_k(2*nband_k**2),eig2_k(2*nband_k**2) 175 real(dp),intent(in) :: ffnl_k(npw_k,dimffnl,psps%lmnmax,psps%ntypat) 176 real(dp),intent(in) :: kpt(3),occ_k(nband_k) 177 real(dp),intent(in) :: rmet(3,3),rprimd(3,3) 178 real(dp),intent(in) :: vpsp1_i1pertdq(2*nfft,nspden,n1dq) 179 real(dp),intent(in) :: vpsp1_i2pertdq(2*nfft,nspden,n2dq) 180 real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 181 real(dp),intent(in) :: ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 182 183 !Local variables------------------------------- 184 !scalars 185 integer :: berryopt,dimffnlk,dimffnl1,iband,idir,idq,ii,jband,nkpg,nkpg1,nylmgrtmp 186 integer :: offset_cgi,opt_gvnl1,optlocal,optnl,reuse_ffnlk,reuse_ffnl1,sij_opt 187 integer :: size_wf,tim_getgh1c,usepaw,usevnl,useylmgr1 188 real(dp) :: cprodi,cprodr,doti,dotr,dum_lambda,fac,tmpim,tmpre 189 real(dp) :: cpu,wall,gflops 190 logical :: with_nonlocal_i1pert,with_nonlocal_i2pert 191 type(rf_hamiltonian_type) :: rf_hamkq 192 193 !arrays 194 real(dp),allocatable :: cg1_aux(:,:),cg1_ddk(:,:,:),cwave0i(:,:),cwave0j(:,:) 195 real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:) 196 real(dp),allocatable :: dkinpw(:),gv1c(:,:) 197 real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:) 198 real(dp),allocatable :: gvloc1dqc(:,:),gvnl1dqc(:,:) 199 real(dp) :: cj_h1_ci(2),dum_grad_berry(1,1),dum_gs1(1,1),dum_gvnl1(1,1) 200 real(dp),allocatable :: kinpw1(:),kpg_k(:,:),kpg1_k(:,:) 201 real(dp),allocatable :: part_ylmgr_k(:,:,:),ph3d(:,:,:),ph3d1(:,:,:) 202 real(dp),allocatable :: dum_vlocal(:,:,:,:),vlocal1(:,:,:,:),dum_vpsp(:) 203 real(dp),allocatable :: vpsp1(:) 204 type(pawcprj_type),allocatable :: dum_cwaveprj(:,:) 205 206 ! ************************************************************************* 207 208 DBG_ENTER("COLL") 209 210 !Additional definitions 211 tim_getgh1c=0 212 useylmgr1=useylmgr 213 usepaw=dtset%usepaw 214 size_wf= dtset%nspinor*npw_k 215 with_nonlocal_i1pert=.true. ; if (i1pert==natom+2) with_nonlocal_i1pert=.false. 216 with_nonlocal_i2pert=.true. ; if (i2pert==natom+2) with_nonlocal_i2pert=.false. 217 reuse_ffnlk=1 ; if (dtset%ffnl_lw==1) reuse_ffnlk=0 218 reuse_ffnl1=1 ; if (dtset%ffnl_lw==1) reuse_ffnl1=0 219 220 !Additional allocations 221 ABI_MALLOC(cwave0i,(2,size_wf)) 222 ABI_MALLOC(cwave0j,(2,size_wf)) 223 ABI_MALLOC(cwavef1,(2,size_wf)) 224 ABI_MALLOC(cwavef2,(2,size_wf)) 225 ABI_MALLOC(cg1_aux,(2,size_wf)) 226 ABI_MALLOC(gv1c,(2,size_wf)) 227 ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 228 ABI_MALLOC(dum_vpsp,(nfft)) 229 ABI_MALLOC(dum_vlocal,(ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 230 ABI_MALLOC(vpsp1,(cplex*nfft)) 231 ABI_MALLOC(dum_cwaveprj,(0,0)) 232 ABI_MALLOC(part_ylmgr_k,(npw_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 233 part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3,:) 234 235 !------------------------------------T1------------------------------------------------ 236 !q1-gradient of gs Hamiltonian: 237 ! < u_{i,k}^{\lambda1}} | \partial_{gamma} H^{(0)} | u_{i,k}^{\lambda2} > 238 !-------------------------------------------------------------------------------------- 239 240 call cwtime(cpu, wall, gflops, "start") 241 242 !Specific definitions 243 d3etot_t1_k=zero 244 vlocal1=zero 245 dum_lambda=zero 246 berryopt=0;optlocal=0;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 247 248 !Initialize rf Hamiltonian (the k-dependent part is prepared in getgh1c_setup) 249 call init_rf_hamiltonian(cplex,gs_hamkq,natom+1,rf_hamkq,& 250 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 251 & mpi_spintab=mpi_enreg%my_isppoltab) 252 253 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 254 255 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 256 if (dtset%ffnl_lw==0) then 257 ABI_MALLOC(ffnlk,(npw_k,0,psps%lmnmax,psps%ntypat)) 258 ABI_MALLOC(ffnl1,(npw_k,2,psps%lmnmax,psps%ntypat)) 259 ffnl1(:,1,:,:)=ffnl_k(:,1,:,:) 260 ffnl1(:,2,:,:)=ffnl_k(:,1+i3dir,:,:) 261 end if 262 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 263 kpt,kpt,i3dir,natom+1,natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 264 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,part_ylmgr_k,& ! In 265 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,& ! Out 266 reuse_ffnlk=reuse_ffnlk, reuse_ffnl1=reuse_ffnl1) ! Optional 267 268 !LOOP OVER BANDS 269 do iband=1,nband_k 270 271 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 272 273 !Select bks wf1 274 offset_cgi = (iband-1)*size_wf+icg 275 cwavef1(:,:)= cg1(:,1+offset_cgi:size_wf+offset_cgi) 276 cwavef2(:,:)= cg2(:,1+offset_cgi:size_wf+offset_cgi) 277 278 !Compute < g |\partial_{gamma} H^{(0)} | u_{i,k}^{\lambda2} > 279 call getgh1c(berryopt,cwavef2,dum_cwaveprj,gv1c,dum_grad_berry,& 280 & dum_gs1,gs_hamkq,dum_gvnl1,i3dir,natom+1,dum_lambda,mpi_enreg,optlocal,& 281 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 282 283 !Apply the dot product with the ket wf (take into account occupation here) 284 ! < u_{i,k}^{\lambda1}} | \partial_{gamma} H^{(0)} | u_{i,k}^{lambda2}} > 285 call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef1,gv1c, & 286 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 287 288 d3etot_t1_k(1)=d3etot_t1_k(1)+occ_k(iband)*dotr 289 d3etot_t1_k(2)=d3etot_t1_k(2)+occ_k(iband)*doti 290 291 end do !iband 292 293 !Clean rf_hamiltonian 294 call rf_hamkq%free() 295 296 !Deallocations 297 ABI_FREE(kpg_k) 298 ABI_FREE(kpg1_k) 299 ABI_FREE(dkinpw) 300 ABI_FREE(kinpw1) 301 ABI_FREE(ffnlk) 302 ABI_FREE(ffnl1) 303 ABI_FREE(ph3d) 304 305 call cwtime(cpu, wall, gflops, "stop") 306 307 !------------------------------------T2------------------------------------------------ 308 !q-gradient of CB projector x rf Hamiltonian lambda 2: 309 ! < u_{i,k}^{\lambda1}} | \partial_{gamma} Q_k H^{\lambda2} | u_{i,k}^{(0)} > 310 !-------------------------------------------------------------------------------------- 311 312 !Create array for ddk 1wf from file 313 ABI_MALLOC(cg1_ddk,(2,size_wf,nband_k)) 314 do iband=1,nband_k 315 call ddk_f%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=cg1_aux) 316 cg1_ddk(:,:,iband)=cg1_aux(:,:) 317 end do 318 319 call cwtime(cpu, wall, gflops, "start") 320 321 !For \lambda1=\lambda2 T2 is inferred from the cc of T3 322 if (.not.samepert) then 323 324 !Specific definitions 325 d3etot_t2_k=zero 326 327 !LOOP OVER BANDS 328 do iband=1,nband_k 329 330 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 331 332 !Select bks wfs 333 offset_cgi = (iband-1)*size_wf+icg 334 cwavef1(:,:)= cg1(:,1+offset_cgi:size_wf+offset_cgi) 335 336 !LOOP OVER BANDS 337 do jband=1,nband_k 338 339 !Select ddk wf1 340 cg1_aux(:,:)=cg1_ddk(:,:,jband) 341 342 !Load < u_{j,k}^{(0) | H^{\lambda2}+V^{\lambda2}} | u_{i,k}^{(0)} > 343 ii=2*jband-1+(iband-1)*2*nband_k 344 cj_h1_ci(1)=eig2_k(ii) 345 cj_h1_ci(2)=eig2_k(ii+1) 346 347 !Calculate: < u_{i,k}^{lambda1}} | u_{j,k}^{k_{\gamma}} > 348 call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef1,cg1_aux, & 349 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 350 351 !Calculate the contribution to T2 352 cprodr=dotr*cj_h1_ci(1)-doti*cj_h1_ci(2) 353 cprodi=dotr*cj_h1_ci(2)+doti*cj_h1_ci(1) 354 d3etot_t2_k(1)=d3etot_t2_k(1)-cprodr*occ_k(iband) 355 d3etot_t2_k(2)=d3etot_t2_k(2)-cprodi*occ_k(iband) 356 357 end do !jband 358 359 end do !iband 360 361 end if !samepert 362 363 call cwtime(cpu, wall, gflops, "stop") 364 365 !------------------------------------T3------------------------------------------------ 366 !rf Hamiltonian lambda 1 x q-gradient of CB projector 367 ! < u_{i,k}^{(0) | (H^{\lambda1})^{\dagger} \partial_{gamma} Q_k | u_{i,k}^{\lambda2}} > 368 !-------------------------------------------------------------------------------------- 369 370 call cwtime(cpu, wall, gflops, "start") 371 !Specific definitions 372 d3etot_t3_k=zero 373 374 !LOOP OVER BANDS 375 do iband=1,nband_k 376 377 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 378 379 !Select bks wfs 380 offset_cgi = (iband-1)*size_wf+icg 381 cwavef2(:,:)= cg2(:,1+offset_cgi:size_wf+offset_cgi) 382 383 !LOOP OVER BANDS 384 do jband=1,nband_k 385 386 !Select ddk wf1 387 cg1_aux(:,:)=cg1_ddk(:,:,jband) 388 389 !Load (< u_{j,k}^{(0) | H^{\lambda1}+V^{\lambda1}} | u_{i,k}^{(0)} >)^* 390 ii=2*jband-1+(iband-1)*2*nband_k 391 cj_h1_ci(1)=eig1_k(ii) 392 cj_h1_ci(2)=-eig1_k(ii+1) 393 394 !Calculate: < u_{j,k}^{k_{\gamma}} | u_{i,k}^{lambda2}} > 395 call dotprod_g(dotr,doti,istwf_k,size_wf,2,cg1_aux,cwavef2, & 396 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 397 398 !Calculate the contribution to T3 399 cprodr=dotr*cj_h1_ci(1)-doti*cj_h1_ci(2) 400 cprodi=dotr*cj_h1_ci(2)+doti*cj_h1_ci(1) 401 d3etot_t3_k(1)=d3etot_t3_k(1)-cprodr*occ_k(iband) 402 d3etot_t3_k(2)=d3etot_t3_k(2)-cprodi*occ_k(iband) 403 404 end do !jband 405 406 end do !iband 407 408 ABI_FREE(cg1_ddk) 409 ABI_FREE(cg1_aux) 410 ABI_FREE(vpsp1) 411 ABI_FREE(vlocal1) 412 413 if (samepert) then 414 d3etot_t2_k(1)=d3etot_t3_k(1) 415 d3etot_t2_k(2)=-d3etot_t3_k(2) 416 end if 417 418 call cwtime(cpu, wall, gflops, "stop") 419 !------------------------------------T4------------------------------------------------ 420 !q-gradient of rf Hamiltonian lambda 2 421 ! < u_{i,k}^{\lambda1} | H^{\lambda2}_{gamma} | u_{i,k}^{(0)} > 422 !-------------------------------------------------------------------------------------- 423 424 !For \lambda1=\lambda2 T4 is inferred from the cc of T5 425 if (.not.samepert) then 426 427 !Specific definitions and allocations 428 d3etot_t4_k=zero 429 optlocal=1;optnl=1 430 dimffnlk=0 431 if (i2pert/=natom+2) then 432 ABI_MALLOC(vlocal1,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 433 ABI_MALLOC(vpsp1,(2*nfft)) 434 ABI_MALLOC(gvloc1dqc,(2,size_wf)) 435 ABI_MALLOC(gvnl1dqc,(2,size_wf)) 436 end if 437 if (i2pert<=natom) fac=-one 438 if (i2pert==natom+2) fac=one 439 if (i2pert==natom+3.or.i2pert==natom+4) fac=-half 440 if (i2pert<=natom) then 441 nylmgrtmp=3 442 dimffnlk=1 443 dimffnl1=2 444 else if (i2pert==natom+3.or.i2pert==natom+4) then 445 nylmgrtmp=nylmgr 446 dimffnl1=10 447 ABI_FREE(part_ylmgr_k) 448 ABI_MALLOC(part_ylmgr_k,(npw_k,nylmgrtmp,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 449 part_ylmgr_k(:,:,:)=ylmgr_k(:,:,:) 450 end if 451 452 !Do loop to compute both extradiagonal shear-strain components 453 do idq=1,n2dq 454 455 call cwtime(cpu, wall, gflops, "start") 456 457 if (i2pert/=natom+2) then 458 idir=i2dir; if (i2pert==natom+4) idir=idq*3+i2dir 459 !Initialize rf Hamiltonian (the k-dependent part is prepared in getgh1c_setup) 460 call init_rf_hamiltonian(2,gs_hamkq,i2pert,rf_hamkq,& 461 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 462 & mpi_spintab=mpi_enreg%my_isppoltab) 463 464 !Set up local potentials with proper dimensioning 465 !and load the spin-dependent part of the Hamiltonians 466 vpsp1=vpsp1_i2pertdq(:,isppol,idq) 467 call rf_transgrid_and_pack(isppol,nspden,usepaw,2,nfft,nfft,ngfft,& 468 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1) 469 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,& 470 & with_nonlocal=with_nonlocal_i2pert) 471 472 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 473 if (dtset%ffnl_lw==0) then 474 ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,psps%ntypat)) 475 if (dimffnlk==1) ffnlk(:,1,:,:)=ffnl_k(:,1,:,:) 476 ABI_MALLOC(ffnl1,(npw_k,dimffnl1,psps%lmnmax,psps%ntypat)) 477 if (dimffnl1==2) then 478 ffnl1(:,1,:,:)=ffnl_k(:,1,:,:) 479 ffnl1(:,2,:,:)=ffnl_k(:,1+i3dir,:,:) 480 else 481 ffnl1(:,1:dimffnl1,:,:)=ffnl_k(:,1:dimffnl1,:,:) 482 end if 483 end if 484 485 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,i2pert,i3dir, & 486 & dtset%natom,rmet,rprimd,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgrtmp,useylmgr1,kg_k, & 487 & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1, & 488 & reuse_ffnlk=reuse_ffnlk,reuse_ffnl1=reuse_ffnl1) 489 490 end if 491 492 !LOOP OVER BANDS 493 do iband=1,nband_k 494 495 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 496 497 !Select bks wfs 498 offset_cgi = (iband-1)*size_wf+icg 499 cwavef1(:,:)= cg1(:,1+offset_cgi:size_wf+offset_cgi) 500 501 !Perturbation-specific part 502 if (i2pert==natom+2) then 503 if (samepert) then 504 ! Read from d2_dkdk_f 505 call d2_dkdk_f%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=gv1c) 506 else 507 ! Read from d2_dkdk_f2 508 call d2_dkdk_f2%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=gv1c) 509 end if 510 else 511 cwave0i(:,:)= cg(:,1+offset_cgi:size_wf+offset_cgi) 512 513 !Compute < g |H^{\lambda2}}_{\gamma} | u_{i,k}^{(0)} > 514 call getgh1dqc(cwave0i,dum_cwaveprj,gv1c,gvloc1dqc,gvnl1dqc,gs_hamkq, & 515 & idir,i2pert,mpi_enreg,optlocal,optnl,i3dir,rf_hamkq) 516 end if 517 518 !Calculate: < u_{j,k}^{\lambda1} | |H^{\lambda2}}_{\gamma} | u_{i,k}^{(0)} > 519 call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef1,gv1c, & 520 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 521 522 !Calculate the contribution to T4 523 d3etot_t4_k(1,idq)=d3etot_t4_k(1,idq)+dotr*occ_k(iband) 524 d3etot_t4_k(2,idq)=d3etot_t4_k(2,idq)+doti*occ_k(iband) 525 526 end do !iband 527 528 if (i2pert/=natom+2) then 529 530 !Clean the rf_hamiltonian 531 call rf_hamkq%free() 532 533 !Deallocations 534 ABI_FREE(kpg_k) 535 ABI_FREE(kpg1_k) 536 ABI_FREE(dkinpw) 537 ABI_FREE(kinpw1) 538 ABI_FREE(ffnlk) 539 ABI_FREE(ffnl1) 540 ABI_FREE(ph3d) 541 542 end if 543 544 call cwtime(cpu, wall, gflops, "stop") 545 546 !Apply the perturbation-dependent prefactors on T4 547 tmpre=d3etot_t4_k(1,idq); tmpim=d3etot_t4_k(2,idq) 548 if (i2pert<=natom.or.i2pert==natom+2) then 549 d3etot_t4_k(1,idq)=-tmpim 550 d3etot_t4_k(2,idq)=tmpre 551 end if 552 d3etot_t4_k(:,idq)=d3etot_t4_k(:,idq)*fac 553 554 end do !idq 555 556 if (i2pert/=natom+2) then 557 ABI_FREE(gvloc1dqc) 558 ABI_FREE(gvnl1dqc) 559 ABI_FREE(vlocal1) 560 ABI_FREE(vpsp1) 561 end if 562 563 end if !samepert 564 565 !------------------------------------T5------------------------------------------------ 566 !q-gradient of rf Hamiltonian lambda 1 567 ! < u_{i,k}^{(0)} | (H^{\lambda1}_{gamma})^{\dagger} | u_{i,k}^{\lambda2} > 568 !-------------------------------------------------------------------------------------- 569 570 !Specific definitions and allocations 571 d3etot_t5_k=zero 572 optlocal=1;optnl=1 573 dimffnlk=0 574 if (i1pert/=natom+2) then 575 ABI_MALLOC(vlocal1,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 576 ABI_MALLOC(vpsp1,(2*nfft)) 577 ABI_MALLOC(gvloc1dqc,(2,size_wf)) 578 ABI_MALLOC(gvnl1dqc,(2,size_wf)) 579 end if 580 if (i1pert<=natom) fac=-one 581 if (i1pert==natom+2) fac=one 582 if (i1pert==natom+3.or.i1pert==natom+4) fac=-half 583 if (i1pert<=natom) then 584 nylmgrtmp=3 585 dimffnlk=1 586 dimffnl1=2 587 ABI_FREE(part_ylmgr_k) 588 ABI_MALLOC(part_ylmgr_k,(npw_k,nylmgrtmp,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 589 part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3,:) 590 else if (i1pert==natom+3.or.i1pert==natom+4) then 591 nylmgrtmp=nylmgr 592 dimffnl1=10 593 ABI_FREE(part_ylmgr_k) 594 ABI_MALLOC(part_ylmgr_k,(npw_k,nylmgrtmp,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 595 part_ylmgr_k(:,:,:)=ylmgr_k(:,:,:) 596 end if 597 598 !Do loop to compute both extradiagonal shear-strain components 599 do idq=1,n1dq 600 601 call cwtime(cpu, wall, gflops, "start") 602 if (i1pert/=natom+2) then 603 idir=i1dir; if (i1pert==natom+4) idir=idq*3+i1dir 604 !Initialize rf Hamiltonian (the k-dependent part is prepared in getgh1c_setup) 605 call init_rf_hamiltonian(2,gs_hamkq,i1pert,rf_hamkq,& 606 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 607 & mpi_spintab=mpi_enreg%my_isppoltab) 608 609 !Set up local potentials with proper dimensioning 610 !and load the spin-dependent part of the Hamiltonians 611 vpsp1=vpsp1_i1pertdq(:,isppol,idq) 612 call rf_transgrid_and_pack(isppol,nspden,usepaw,2,nfft,nfft,ngfft,& 613 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1) 614 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,& 615 & with_nonlocal=with_nonlocal_i1pert) 616 617 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 618 if (dtset%ffnl_lw==0) then 619 ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,psps%ntypat)) 620 if (dimffnlk==1) ffnlk(:,1,:,:)=ffnl_k(:,1,:,:) 621 ABI_MALLOC(ffnl1,(npw_k,dimffnl1,psps%lmnmax,psps%ntypat)) 622 if (dimffnl1==2) then 623 ffnl1(:,1,:,:)=ffnl_k(:,1,:,:) 624 ffnl1(:,2,:,:)=ffnl_k(:,1+i3dir,:,:) 625 else 626 ffnl1(:,1:dimffnl1,:,:)=ffnl_k(:,1:dimffnl1,:,:) 627 end if 628 end if 629 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,i1pert,i3dir, & 630 & dtset%natom,rmet,rprimd,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgrtmp,useylmgr1,kg_k, & 631 & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1, & 632 & reuse_ffnlk=reuse_ffnlk,reuse_ffnl1=reuse_ffnl1) 633 634 end if 635 636 !LOOP OVER BANDS 637 do iband=1,nband_k 638 639 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 640 641 !Select bks wfs 642 offset_cgi = (iband-1)*size_wf+icg 643 cwavef2(:,:)= cg2(:,1+offset_cgi:size_wf+offset_cgi) 644 645 !Perturbation-specific part 646 if (i1pert==natom+2) then 647 call d2_dkdk_f%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=gv1c) 648 else 649 cwave0i(:,:)= cg(:,1+offset_cgi:size_wf+offset_cgi) 650 651 !Compute < g |H^{\lambda1}}_{\gamma} | u_{i,k}^{(0)} > 652 call getgh1dqc(cwave0i,dum_cwaveprj,gv1c,gvloc1dqc,gvnl1dqc,gs_hamkq, & 653 & idir,i1pert,mpi_enreg,optlocal,optnl,i3dir,rf_hamkq) 654 end if 655 656 !Calculate: < u_{j,k}^{\lambda2} | |H^{\lambda1}}_{\gamma} | u_{i,k}^{(0)} > 657 call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef2,gv1c, & 658 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 659 660 !Calculate the contribution to T5: 661 d3etot_t5_k(1,idq)=d3etot_t5_k(1,idq)+dotr*occ_k(iband) 662 d3etot_t5_k(2,idq)=d3etot_t5_k(2,idq)+doti*occ_k(iband) 663 664 end do !iband 665 666 if (i1pert/=natom+2) then 667 668 !Clean the rf_hamiltonian 669 call rf_hamkq%free() 670 671 !Deallocations 672 ABI_FREE(kpg_k) 673 ABI_FREE(kpg1_k) 674 ABI_FREE(dkinpw) 675 ABI_FREE(kinpw1) 676 ABI_FREE(ffnlk) 677 ABI_FREE(ffnl1) 678 ABI_FREE(ph3d) 679 680 end if 681 682 call cwtime(cpu, wall, gflops, "stop") 683 684 !Apply the perturbation-dependent prefactors on T5 685 tmpre=d3etot_t5_k(1,idq); tmpim=d3etot_t5_k(2,idq) 686 if (i1pert<=natom.or.i1pert==natom+2) then 687 d3etot_t5_k(1,idq)=-tmpim 688 d3etot_t5_k(2,idq)=tmpre 689 end if 690 d3etot_t5_k(:,idq)=d3etot_t5_k(:,idq)*fac 691 692 !Apply now the conjugate complex: 693 !(< u_{j,k}^{\lambda2} | |H^{\lambda1}}_{\gamma} | u_{i,k}^{(0)} >)* 694 tmpim=d3etot_t5_k(2,idq) 695 d3etot_t5_k(2,idq)=-tmpim 696 697 end do !idq 698 699 700 if (i1pert/=natom+2) then 701 ABI_FREE(gvloc1dqc) 702 ABI_FREE(gvnl1dqc) 703 ABI_FREE(vlocal1) 704 ABI_FREE(vpsp1) 705 end if 706 707 if (samepert) then 708 d3etot_t4_k(1,:)=d3etot_t5_k(1,:) 709 d3etot_t4_k(2,:)=-d3etot_t5_k(2,:) 710 end if 711 712 !Scale d3etot_k contributions by the kpt weight 713 d3etot_t1_k(:)=d3etot_t1_k(:)*wtk_k 714 d3etot_t2_k(:)=d3etot_t2_k(:)*wtk_k 715 d3etot_t3_k(:)=d3etot_t3_k(:)*wtk_k 716 d3etot_t4_k(:,:)=d3etot_t4_k(:,:)*wtk_k 717 d3etot_t5_k(:,:)=d3etot_t5_k(:,:)*wtk_k 718 719 !Deallocations 720 ABI_FREE(cwave0i) 721 ABI_FREE(cwave0j) 722 ABI_FREE(cwavef1) 723 ABI_FREE(cwavef2) 724 ABI_FREE(gv1c) 725 ABI_FREE(dum_vpsp) 726 ABI_FREE(dum_vlocal) 727 ABI_FREE(dum_cwaveprj) 728 ABI_FREE(part_ylmgr_k) 729 730 731 DBG_EXIT("COLL") 732 733 end subroutine dfpt_1wf