TABLE OF CONTENTS
ABINIT/dfpt_eltfrloc [ Functions ]
NAME
dfpt_eltfrloc
FUNCTION
Compute the frozen-wavefunction local pseudopotential contribution to the elastic tensor and the internal strain (derivative wrt one cartesian strain component and one reduced-coordinate atomic displacement).
COPYRIGHT
Copyright (C) 2000-2018 ABINIT group (DRH, XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
atindx(natom)=index table for atoms (see gstate.f) gmet(3,3)=reciprocal space metric (bohr^-2) gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1) gsqcut=cutoff on G^2 based on ecut mgfft=maximum size of 1D FFTs mpi_enreg=informations about MPI parallelization mqgrid=dimensioned number of q grid points for local psp spline natom=number of atoms in unit cell nattyp(ntypat)=number of atoms of each type nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ntypat=number of types of atoms ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information qgrid(mqgrid)=q point array for local psp spline fits rhog(2,nfft)=electron density in G space vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
OUTPUT
eltfrloc(6+3*natom,6)=non-symmetrized local pseudopotenial contribution to the elastic tensor and internal strain.
PARENTS
respfn
CHILDREN
ptabs_fourdp,timab,xmpi_sum
SOURCE
49 #if defined HAVE_CONFIG_H 50 #include "config.h" 51 #endif 52 53 #include "abi_common.h" 54 55 56 subroutine dfpt_eltfrloc(atindx,eltfrloc,gmet,gprimd,gsqcut,mgfft,& 57 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,ntypat,ph1d,qgrid,rhog,vlspl) 58 59 60 use defs_basis 61 use defs_abitypes 62 use m_errors 63 use m_profiling_abi 64 use m_xmpi 65 66 use m_mpinfo, only : ptabs_fourdp 67 68 !This section has been created automatically by the script Abilint (TD). 69 !Do not modify the following lines by hand. 70 #undef ABI_FUNC 71 #define ABI_FUNC 'dfpt_eltfrloc' 72 use interfaces_18_timing 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: mgfft,mqgrid,natom,nfft,ntypat 80 real(dp),intent(in) :: gsqcut 81 type(MPI_type),intent(in) :: mpi_enreg 82 !arrays 83 integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18) 84 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom) 85 real(dp),intent(in) :: qgrid(mqgrid),rhog(2,nfft),vlspl(mqgrid,2,ntypat) 86 real(dp),intent(out) :: eltfrloc(6+3*natom,6) 87 88 !Local variables------------------------------- 89 !scalars 90 integer,parameter :: im=2,re=1 91 integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ielt,ieltx,ierr,ig1,ig2,ig3,ii 92 integer :: is1,is2,itypat,jj,ka,kb,kd,kg,me_fft,n1,n2,n3,nproc_fft 93 real(dp),parameter :: tolfix=1.0000001_dp 94 real(dp) :: aa,bb,cc,cutoff,d2g 95 real(dp) :: dd,dg1,dg2,diff,dq 96 real(dp) :: dq2div6,dqdiv6,dqm1,ee,ff,gmag,gsquar 97 real(dp) :: sfi,sfr,term,term1 98 !real(dp) :: ph1_elt,ph2_elt,ph3_elt,phi_elt,phr_elt 99 real(dp) :: term2,term3,term4,term5,vion1,vion2,vion3 100 !arrays 101 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 102 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 103 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 104 real(dp) :: dgm(3,3,6),tsec(2) 105 real(dp),allocatable :: d2gm(:,:,:,:),elt_work(:,:) 106 107 ! ************************************************************************* 108 109 !Define G^2 based on G space metric gmet. 110 ! gsq_elt(i1,i2,i3)=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 111 !& dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 112 !& dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 113 114 !Define dG^2/ds based on G space metric derivative 115 ! dgsqds_elt(i1,i2,i3,is)=dble(i1*i1)*dgm(1,1,is)+dble(i2*i2)*dgm(2,2,is)+& 116 !& dble(i3*i3)*dgm(3,3,is)+& 117 !& dble(i1*i2)*(dgm(1,2,is)+dgm(2,1,is))+& 118 !& dble(i1*i3)*(dgm(1,3,is)+dgm(3,1,is))+& 119 !& dble(i2*i3)*(dgm(2,3,is)+dgm(3,2,is)) 120 121 !Define 2dG^2/ds1ds2 based on G space metric derivative 122 ! d2gsqds_elt(i1,i2,i3,is1,is2)=dble(i1*i1)*d2gm(1,1,is1,is2)+& 123 !& dble(i2*i2)*d2gm(2,2,is1,is2)+dble(i3*i3)*d2gm(3,3,is1,is2)+& 124 !& dble(i1*i2)*(d2gm(1,2,is1,is2)+d2gm(2,1,is1,is2))+& 125 !& dble(i1*i3)*(d2gm(1,3,is1,is2)+d2gm(3,1,is1,is2))+& 126 !& dble(i2*i3)*(d2gm(2,3,is1,is2)+d2gm(3,2,is1,is2)) 127 128 !Real and imaginary parts of phase--statment functions: 129 ! phr_elt(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 130 ! phi_elt(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 131 ! ph1_elt(nri,i1,ia)=ph1d(nri,i1+1+n1+(ia-1)*(2*n1+1)) 132 ! ph2_elt(nri,i2,ia)=ph1d(nri,i2+1+n2+(ia-1)*(2*n2+1)+& 133 !& natom*(2*n1+1)) 134 ! ph3_elt(nri,i3,ia)=ph1d(nri,i3+1+n3+(ia-1)*(2*n3+1)+& 135 !& natom*(2*n1+1+2*n2+1)) 136 ! phre_elt(i1,i2,i3,ia)=phr_elt(ph1_elt(re,i1,ia),ph1_elt(im,i1,ia),ph2_elt(re,i2,ia),& 137 !& ph2_elt(im,i2,ia),ph3_elt(re,i3,ia),ph3_elt(im,i3,ia)) 138 ! phimag_elt(i1,i2,i3,ia)=phi_elt(ph1_elt(re,i1,ia),ph1_elt(im,i1,ia),ph2_elt(re,i2,ia),& 139 !& ph2_elt(im,i2,ia),ph3_elt(re,i3,ia),ph3_elt(im,i3,ia)) 140 141 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 142 me_fft=ngfft(11) ; nproc_fft=ngfft(10) 143 144 !Get the distrib associated with this fft_grid 145 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 146 147 !----- 148 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components 149 !and store for use in inner loop below. 150 ABI_ALLOCATE(d2gm,(3,3,6,6)) 151 152 !Loop over 2nd strain index 153 do is2=1,6 154 kg=idx(2*is2-1);kd=idx(2*is2) 155 do jj = 1,3 156 dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj)) 157 end do 158 ! Loop over 1st strain index, upper triangle only 159 do is1=1,is2 160 ka=idx(2*is1-1);kb=idx(2*is1) 161 d2gm(:,:,is1,is2)=0._dp 162 do jj = 1,3 163 if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 164 & +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj) 165 if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 166 & +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj) 167 if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 168 & +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj) 169 if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 170 & +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj) 171 end do 172 end do !is1 173 end do !is2 174 175 !Zero out array to permit accumulation over atom types below: 176 eltfrloc(:,:)=0.0_dp 177 178 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 179 dqm1=1.0_dp/dq 180 dqdiv6=dq/6.0_dp 181 dq2div6=dq**2/6.0_dp 182 cutoff=gsqcut*tolfix 183 id1=n1/2+2 184 id2=n2/2+2 185 id3=n3/2+2 186 187 ia1=1 188 do itypat=1,ntypat 189 ! ia1,ia2 sets range of loop over atoms: 190 ia2=ia1+nattyp(itypat)-1 191 ii=0 192 do i3=1,n3 193 ig3=i3-(i3/id3)*n3-1 194 do i2=1,n2 195 if (fftn2_distrib(i2)==me_fft) then 196 ig2=i2-(i2/id2)*n2-1 197 do i1=1,n1 198 ig1=i1-(i1/id1)*n1-1 199 200 ii=ii+1 201 ! Skip G=0: 202 if (ig1==0 .and. ig2==0 .and. ig3==0) cycle 203 204 ! Skip G**2 outside cutoff: 205 gsquar=gsq_elt(ig1,ig2,ig3) 206 if (gsquar<=cutoff) then 207 gmag=sqrt(gsquar) 208 209 ! Compute vion(G) for given type of atom 210 jj=1+int(gmag*dqm1) 211 diff=gmag-qgrid(jj) 212 213 ! Evaluate spline fit from q^2 V(q) to get V(q): 214 ! (p. 86 Numerical Recipes, Press et al; NOTE error in book for sign 215 ! of "aa" term in derivative; also see splfit routine). 216 bb = diff*dqm1 217 aa = 1.0_dp-bb 218 cc = aa*(aa**2-1.0_dp)*dq2div6 219 dd = bb*(bb**2-1.0_dp)*dq2div6 220 term1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +& 221 & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat)) 222 vion1=term1 / gsquar 223 224 ! Also get dV(q)/dq: 225 ! (note correction of Numerical Recipes sign error 226 ! before (3._dp*aa**2-1._dp) 227 ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat) 228 ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) & 229 & - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat) 230 term2 = ee*dqm1 + ff*dqdiv6 231 vion2 = term2/gsquar - 2._dp*term1/(gsquar*gmag) 232 233 ! Also get V''(q) 234 term3=aa*vlspl(jj,2,itypat)+bb*vlspl(jj+1,2,itypat) 235 vion3 = (term3 - 4.0_dp*term2/gmag + 6._dp*term1/gsquar)/gsquar 236 237 ! Assemble structure factor over all atoms of given type: 238 sfr=zero;sfi=zero 239 do ia=ia1,ia2 240 sfr=sfr+phre_elt(ig1,ig2,ig3,ia) 241 sfi=sfi-phimag_elt(ig1,ig2,ig3,ia) 242 end do 243 term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi) 244 245 ! Loop over 2nd strain index 246 do is2=1,6 247 dg2=0.5_dp*dgsqds_elt(ig1,ig2,ig3,is2)/gmag 248 ! Loop over 1st strain index, upper triangle only 249 do is1=1,is2 250 dg1=0.5_dp*dgsqds_elt(ig1,ig2,ig3,is1)/gmag 251 d2g=(0.25_dp*d2gsqds_elt(ig1,ig2,ig3,is1,is2)-dg1*dg2)/gmag 252 eltfrloc(is1,is2)=eltfrloc(is1,is2)+& 253 & term*(vion3*dg1*dg2+vion2*d2g) 254 if(is2<=3)& 255 & eltfrloc(is1,is2)=eltfrloc(is1,is2)-term*vion2*dg1 256 if(is1<=3)& 257 & eltfrloc(is1,is2)=eltfrloc(is1,is2)-term*vion2*dg2 258 if(is1<=3 .and. is2<=3)& 259 & eltfrloc(is1,is2)=eltfrloc(is1,is2)+term*vion1 260 end do !is1 261 262 ! Internal strain section - loop over current atoms 263 do ia=ia1,ia2 264 if(is2 <=3) then 265 term4=vion2*dg2-vion1 266 else 267 term4=vion2*dg2 268 end if 269 term5=-two_pi*(rhog(re,ii)*phimag_elt(ig1,ig2,ig3,ia)& 270 & +rhog(im,ii)*phre_elt(ig1,ig2,ig3,ia))*term4 271 eltfrloc(7+3*(ia-1),is2)=eltfrloc(7+3*(ia-1),is2)+term5*dble(ig1) 272 eltfrloc(8+3*(ia-1),is2)=eltfrloc(8+3*(ia-1),is2)+term5*dble(ig2) 273 eltfrloc(9+3*(ia-1),is2)=eltfrloc(9+3*(ia-1),is2)+term5*dble(ig3) 274 end do 275 276 end do !is2 277 278 ! End skip G**2 outside cutoff: 279 end if 280 281 ! End loop on n1, n2, n3. There is a "cycle" inside the loop 282 end do 283 end if 284 end do 285 end do 286 287 ! End loop on type of atoms 288 ia1=ia2+1 289 end do 290 !Init mpi_comm 291 call timab(48,1,tsec) 292 call xmpi_sum(eltfrloc,mpi_enreg%comm_fft,ierr) 293 call timab(48,2,tsec) 294 295 !Fill in lower triangle 296 do is2=2,6 297 do is1=1,is2-1 298 eltfrloc(is2,is1)=eltfrloc(is1,is2) 299 end do 300 end do 301 302 !The indexing array atindx is used to reestablish the correct 303 !order of atoms 304 ABI_ALLOCATE(elt_work,(6+3*natom,6)) 305 elt_work(1:6,1:6)=eltfrloc(1:6,1:6) 306 do ia=1,natom 307 ielt=7+3*(ia-1) 308 ieltx=7+3*(atindx(ia)-1) 309 elt_work(ielt:ielt+2,1:6)=eltfrloc(ieltx:ieltx+2,1:6) 310 end do 311 eltfrloc(:,:)=elt_work(:,:) 312 313 ABI_DEALLOCATE(d2gm) 314 ABI_DEALLOCATE(elt_work) 315 316 contains 317 318 !Real and imaginary parts of phase. 319 function phr_elt(x1,y1,x2,y2,x3,y3) 320 321 322 !This section has been created automatically by the script Abilint (TD). 323 !Do not modify the following lines by hand. 324 #undef ABI_FUNC 325 #define ABI_FUNC 'phr_elt' 326 !End of the abilint section 327 328 real(dp) :: phr_elt,x1,x2,x3,y1,y2,y3 329 phr_elt=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 330 end function phr_elt 331 332 function phi_elt(x1,y1,x2,y2,x3,y3) 333 334 335 !This section has been created automatically by the script Abilint (TD). 336 !Do not modify the following lines by hand. 337 #undef ABI_FUNC 338 #define ABI_FUNC 'phi_elt' 339 !End of the abilint section 340 341 real(dp):: phi_elt,x1,x2,x3,y1,y2,y3 342 phi_elt=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 343 end function phi_elt 344 345 function ph1_elt(nri,ig1,ia) 346 347 348 !This section has been created automatically by the script Abilint (TD). 349 !Do not modify the following lines by hand. 350 #undef ABI_FUNC 351 #define ABI_FUNC 'ph1_elt' 352 !End of the abilint section 353 354 real(dp):: ph1_elt 355 integer :: nri,ig1,ia 356 ph1_elt=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1)) 357 end function ph1_elt 358 359 function ph2_elt(nri,ig2,ia) 360 361 362 !This section has been created automatically by the script Abilint (TD). 363 !Do not modify the following lines by hand. 364 #undef ABI_FUNC 365 #define ABI_FUNC 'ph2_elt' 366 !End of the abilint section 367 368 real(dp):: ph2_elt 369 integer :: nri,ig2,ia 370 ph2_elt=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)) 371 end function ph2_elt 372 373 function ph3_elt(nri,ig3,ia) 374 375 376 !This section has been created automatically by the script Abilint (TD). 377 !Do not modify the following lines by hand. 378 #undef ABI_FUNC 379 #define ABI_FUNC 'ph3_elt' 380 !End of the abilint section 381 382 real(dp):: ph3_elt 383 integer :: nri,ig3,ia 384 ph3_elt=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 385 end function ph3_elt 386 387 function phre_elt(ig1,ig2,ig3,ia) 388 389 390 !This section has been created automatically by the script Abilint (TD). 391 !Do not modify the following lines by hand. 392 #undef ABI_FUNC 393 #define ABI_FUNC 'phre_elt' 394 !End of the abilint section 395 396 real(dp):: phre_elt 397 integer :: ig1,ig2,ig3,ia 398 phre_elt=phr_elt(ph1_elt(re,ig1,ia),ph1_elt(im,ig1,ia),& 399 & ph2_elt(re,ig2,ia),ph2_elt(im,ig2,ia),ph3_elt(re,ig3,ia),ph3_elt(im,ig3,ia)) 400 end function phre_elt 401 402 function phimag_elt(ig1,ig2,ig3,ia) 403 404 405 !This section has been created automatically by the script Abilint (TD). 406 !Do not modify the following lines by hand. 407 #undef ABI_FUNC 408 #define ABI_FUNC 'phimag_elt' 409 !End of the abilint section 410 411 real(dp) :: phimag_elt 412 integer :: ig1,ig2,ig3,ia 413 phimag_elt=phi_elt(ph1_elt(re,ig1,ia),ph1_elt(im,ig1,ia),& 414 & ph2_elt(re,ig2,ia),ph2_elt(im,ig2,ia),ph3_elt(re,ig3,ia),ph3_elt(im,ig3,ia)) 415 end function phimag_elt 416 417 function gsq_elt(i1,i2,i3) 418 419 420 !This section has been created automatically by the script Abilint (TD). 421 !Do not modify the following lines by hand. 422 #undef ABI_FUNC 423 #define ABI_FUNC 'gsq_elt' 424 !End of the abilint section 425 426 real(dp) :: gsq_elt 427 integer :: i1,i2,i3 428 !Define G^2 based on G space metric gmet. 429 gsq_elt=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 430 & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 431 & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 432 end function gsq_elt 433 434 function dgsqds_elt(i1,i2,i3,is) 435 436 437 !This section has been created automatically by the script Abilint (TD). 438 !Do not modify the following lines by hand. 439 #undef ABI_FUNC 440 #define ABI_FUNC 'dgsqds_elt' 441 !End of the abilint section 442 443 real(dp) :: dgsqds_elt 444 integer :: i1,i2,i3,is 445 !Define dG^2/ds based on G space metric derivative 446 dgsqds_elt=dble(i1*i1)*dgm(1,1,is)+dble(i2*i2)*dgm(2,2,is)+& 447 & dble(i3*i3)*dgm(3,3,is)+& 448 & dble(i1*i2)*(dgm(1,2,is)+dgm(2,1,is))+& 449 & dble(i1*i3)*(dgm(1,3,is)+dgm(3,1,is))+& 450 & dble(i2*i3)*(dgm(2,3,is)+dgm(3,2,is)) 451 end function dgsqds_elt 452 453 function d2gsqds_elt(i1,i2,i3,is1,is2) 454 455 456 !This section has been created automatically by the script Abilint (TD). 457 !Do not modify the following lines by hand. 458 #undef ABI_FUNC 459 #define ABI_FUNC 'd2gsqds_elt' 460 !End of the abilint section 461 462 real(dp) :: d2gsqds_elt 463 integer :: i1,i2,i3,is1,is2 464 !Define 2dG^2/ds1ds2 based on G space metric derivative 465 d2gsqds_elt=dble(i1*i1)*d2gm(1,1,is1,is2)+& 466 & dble(i2*i2)*d2gm(2,2,is1,is2)+dble(i3*i3)*d2gm(3,3,is1,is2)+& 467 & dble(i1*i2)*(d2gm(1,2,is1,is2)+d2gm(2,1,is1,is2))+& 468 & dble(i1*i3)*(d2gm(1,3,is1,is2)+d2gm(3,1,is1,is2))+& 469 & dble(i2*i3)*(d2gm(2,3,is1,is2)+d2gm(3,2,is1,is2)) 470 end function d2gsqds_elt 471 472 end subroutine dfpt_eltfrloc