TABLE OF CONTENTS
ABINIT/vlocalstr [ Functions ]
NAME
vlocalstr
FUNCTION
Compute strain derivatives of local ionic potential
second derivative of E wrt xred
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DRH, DCA, XG, GMR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$). gprimd(3,3)=reciprocal space dimensional primitive translations gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere). istr=1,...6 specifies cartesian strain component 11,22,33,32,31,21 mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization mqgrid=number of grid pts in q array for f(q) spline. natom=number of atoms in unit cell. nattyp(ntypat)=number of atoms of each type in cell. 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)=1-dim structure factor phase information. qgrid(mqgrid)=q grid for spline from 0 to qmax. ucvol=unit cell volume ($\textrm{Bohr}^3$). vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
OUTPUT
vpsp1(nfft)=first-order local crystal pseudopotential in real space.
SIDE EFFECTS
NOTES
* Note that the present routine is tightly connected to the dfpt_vlocal.f routine, that compute the derivative of the local ionic potential with respect to one atomic displacement. The argument list and the internal loops to be considered were sufficiently different as to make the two routines different. * The routine was adapted from mklocl.F90
PARENTS
dfpt_looppert,dfpt_nselt,dfpt_nstpaw
CHILDREN
fourdp,ptabs_fourdp
SOURCE
56 #if defined HAVE_CONFIG_H 57 #include "config.h" 58 #endif 59 60 #include "abi_common.h" 61 62 63 subroutine vlocalstr(gmet,gprimd,gsqcut,istr,mgfft,mpi_enreg,& 64 & mqgrid,natom,nattyp,nfft,ngfft,ntypat,paral_kgb,ph1d,qgrid,& 65 & ucvol,vlspl,vpsp1) 66 67 use defs_basis 68 use defs_abitypes 69 use m_errors 70 use m_profiling_abi 71 72 use m_mpinfo, only : ptabs_fourdp 73 74 !This section has been created automatically by the script Abilint (TD). 75 !Do not modify the following lines by hand. 76 #undef ABI_FUNC 77 #define ABI_FUNC 'vlocalstr' 78 use interfaces_53_ffts 79 !End of the abilint section 80 81 implicit none 82 83 !Arguments ------------------------------------ 84 !scalars 85 integer,intent(in) :: istr,mgfft,mqgrid,natom,nfft,ntypat,paral_kgb 86 real(dp),intent(in) :: gsqcut,ucvol 87 type(MPI_type),intent(in) :: mpi_enreg 88 !arrays 89 integer,intent(in) :: nattyp(ntypat),ngfft(18) 90 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom) 91 real(dp),intent(in) :: qgrid(mqgrid),vlspl(mqgrid,2,ntypat) 92 real(dp),intent(out) :: vpsp1(nfft) 93 94 !Local variables------------------------------- 95 !scalars 96 integer,parameter :: im=2,re=1 97 integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ig1,ig2,ig3,ii,itypat,jj 98 integer :: ka,kb,n1,n2,n3 99 real(dp),parameter :: tolfix=1.0000001_dp 100 real(dp) :: aa,bb,cc,cutoff,dd,dgsquards,diff 101 real(dp) :: dq,dq2div6,dqdiv6,dqm1,ee,ff,gmag,gsquar 102 real(dp) :: sfi,sfr,term,vion1,vion2 103 real(dp) :: xnorm 104 character(len=500) :: message 105 !arrays 106 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 107 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 108 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 109 real(dp) :: dgmetds(3,3) 110 real(dp),allocatable :: work1(:,:) 111 112 ! ************************************************************************* 113 114 !Define G^2 based on G space metric gmet. 115 ! gsq_vl(i1,i2,i3)=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 116 !& dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 117 !& dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 118 119 !Define dG^2/ds based on G space metric derivative dgmetds. 120 ! dgsqds_vl(i1,i2,i3)=dble(i1*i1)*dgmetds(1,1)+dble(i2*i2)*dgmetds(2,2)+& 121 !& dble(i3*i3)*dgmetds(3,3)+& 122 !& dble(i1*i2)*(dgmetds(1,2)+dgmetds(2,1))+& 123 !& dble(i1*i3)*(dgmetds(1,3)+dgmetds(3,1))+& 124 !& dble(i2*i3)*(dgmetds(2,3)+dgmetds(3,2)) 125 126 !Real and imaginary parts of phase--statment functions: 127 ! phr_vl(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 128 ! phi_vl(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 129 ! ph1_vl(nri,i1,ia)=ph1d(nri,i1+1+n1+(ia-1)*(2*n1+1)) 130 ! ph2_vl(nri,i2,ia)=ph1d(nri,i2+1+n2+(ia-1)*(2*n2+1)+& 131 !& natom*(2*n1+1)) 132 ! ph3_vl(nri,i3,ia)=ph1d(nri,i3+1+n3+(ia-1)*(2*n3+1)+& 133 !& natom*(2*n1+1+2*n2+1)) 134 ! phre_vl(i1,i2,i3,ia)=phr_vl(ph1_vl(re,i1,ia),ph1_vl(im,i1,ia),ph2_vl(re,i2,ia),& 135 !& ph2_vl(im,i2,ia),ph3_vl(re,i3,ia),ph3_vl(im,i3,ia)) 136 ! phimag_vl(i1,i2,i3,ia)=phi_vl(ph1_vl(re,i1,ia),ph1_vl(im,i1,ia),ph2_vl(re,i2,ia),& 137 !& ph2_vl(im,i2,ia),ph3_vl(re,i3,ia),ph3_vl(im,i3,ia)) 138 139 !----- 140 !Compute derivative of metric tensor wrt strain component istr 141 if(istr<1 .or. istr>6)then 142 write(message, '(a,i10,a,a,a)' )& 143 & ' Input istr=',istr,' not allowed.',ch10,& 144 & ' Possible values are 1,2,3,4,5,6 only.' 145 MSG_BUG(message) 146 end if 147 148 ka=idx(2*istr-1);kb=idx(2*istr) 149 do ii = 1,3 150 dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii)) 151 end do 152 !For historical reasons: 153 dgmetds(:,:)=0.5_dp*dgmetds(:,:) 154 155 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 156 157 !Get the distrib associated with this fft_grid 158 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 159 160 !Zero out array to permit accumulation over atom types below: 161 ABI_ALLOCATE(work1,(2,nfft)) 162 work1(:,:)=0.0_dp 163 ! 164 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 165 dqm1=1.0_dp/dq 166 dqdiv6=dq/6.0_dp 167 dq2div6=dq**2/6.0_dp 168 cutoff=gsqcut*tolfix 169 id1=n1/2+2 170 id2=n2/2+2 171 id3=n3/2+2 172 173 ia1=1 174 do itypat=1,ntypat 175 ! ia1,ia2 sets range of loop over atoms: 176 ia2=ia1+nattyp(itypat)-1 177 178 ii=0 179 do i3=1,n3 180 ig3=i3-(i3/id3)*n3-1 181 do i2=1,n2 182 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 183 ig2=i2-(i2/id2)*n2-1 184 do i1=1,n1 185 ig1=i1-(i1/id1)*n1-1 186 ii=ii+1 187 ! *** GET RID OF THIS THESE IF STATEMENTS (if they slow code) 188 ! Skip G=0: 189 ! if (ii==1) cycle 190 if (ig1==0 .and. ig2==0 .and. ig3==0) cycle 191 gsquar=gsq_vl(ig1,ig2,ig3) 192 193 ! Skip G**2 outside cutoff: 194 if (gsquar<=cutoff) then 195 gmag=sqrt(gsquar) 196 dgsquards=dgsqds_vl(ig1,ig2,ig3) 197 ! Compute vion(G) for given type of atom 198 jj=1+int(gmag*dqm1) 199 diff=gmag-qgrid(jj) 200 201 ! Evaluate spline fit from q^2 V(q) to get V(q): 202 ! (p. 86 Numerical Recipes, Press et al; 203 ! NOTE error in book for sign 204 ! of "aa" term in derivative; also see splfit routine). 205 206 bb = diff*dqm1 207 aa = 1.0_dp-bb 208 cc = aa*(aa**2-1.0_dp)*dq2div6 209 dd = bb*(bb**2-1.0_dp)*dq2div6 210 211 vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +& 212 & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) & 213 & / gsquar 214 215 ! Also get (dV(q)/dq)/q: 216 ! (note correction of Numerical Recipes sign error 217 ! before (3._dp*aa**2-1._dp) 218 ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat) 219 ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) & 220 & - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat) 221 vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag& 222 & - 2.0_dp*vion1 ) / gsquar 223 224 225 ! Assemble structure factor over all atoms of given type: 226 sfr=0.0_dp 227 sfi=0.0_dp 228 do ia=ia1,ia2 229 sfr=sfr+phre_vl(ig1,ig2,ig3,ia) 230 sfi=sfi-phimag_vl(ig1,ig2,ig3,ia) 231 end do 232 233 term=dgsquards*vion2 234 ! Add potential for diagonal strain components 235 if(istr <=3) then 236 term=term-vion1 237 end if 238 239 ! Multiply structure factor times vion derivatives: 240 work1(re,ii)=work1(re,ii)+sfr*term 241 work1(im,ii)=work1(im,ii)+sfi*term 242 243 ! End skip G**2 outside cutoff: 244 end if 245 ! End loop on n1, n2, n3. There is a "cycle" inside the loop 246 end do 247 end if 248 end do 249 end do 250 251 ia1=ia2+1 252 253 ! End loop on type of atoms 254 end do 255 256 257 !Set Vloc(G=0)=0: 258 work1(re,1)=0.0_dp 259 work1(im,1)=0.0_dp 260 261 262 !Transform back to real space 263 call fourdp(1,work1,vpsp1,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 264 265 !Divide by unit cell volume 266 xnorm=1.0_dp/ucvol 267 vpsp1(:)=vpsp1(:)*xnorm 268 269 ABI_DEALLOCATE(work1) 270 271 contains 272 273 !Real and imaginary parts of phase. 274 function phr_vl(x1,y1,x2,y2,x3,y3) 275 276 277 !This section has been created automatically by the script Abilint (TD). 278 !Do not modify the following lines by hand. 279 #undef ABI_FUNC 280 #define ABI_FUNC 'phr_vl' 281 !End of the abilint section 282 283 real(dp) :: phr_vl,x1,x2,x3,y1,y2,y3 284 phr_vl=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 285 end function phr_vl 286 287 function phi_vl(x1,y1,x2,y2,x3,y3) 288 289 290 !This section has been created automatically by the script Abilint (TD). 291 !Do not modify the following lines by hand. 292 #undef ABI_FUNC 293 #define ABI_FUNC 'phi_vl' 294 !End of the abilint section 295 296 real(dp):: phi_vl,x1,x2,x3,y1,y2,y3 297 phi_vl=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 298 end function phi_vl 299 300 function ph1_vl(nri,ig1,ia) 301 302 303 !This section has been created automatically by the script Abilint (TD). 304 !Do not modify the following lines by hand. 305 #undef ABI_FUNC 306 #define ABI_FUNC 'ph1_vl' 307 !End of the abilint section 308 309 real(dp):: ph1_vl 310 integer :: nri,ig1,ia 311 ph1_vl=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1)) 312 end function ph1_vl 313 314 function ph2_vl(nri,ig2,ia) 315 316 317 !This section has been created automatically by the script Abilint (TD). 318 !Do not modify the following lines by hand. 319 #undef ABI_FUNC 320 #define ABI_FUNC 'ph2_vl' 321 !End of the abilint section 322 323 real(dp):: ph2_vl 324 integer :: nri,ig2,ia 325 ph2_vl=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)) 326 end function ph2_vl 327 328 function ph3_vl(nri,ig3,ia) 329 330 331 !This section has been created automatically by the script Abilint (TD). 332 !Do not modify the following lines by hand. 333 #undef ABI_FUNC 334 #define ABI_FUNC 'ph3_vl' 335 !End of the abilint section 336 337 real(dp):: ph3_vl 338 integer :: nri,ig3,ia 339 ph3_vl=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 340 end function ph3_vl 341 342 function phre_vl(ig1,ig2,ig3,ia) 343 344 345 !This section has been created automatically by the script Abilint (TD). 346 !Do not modify the following lines by hand. 347 #undef ABI_FUNC 348 #define ABI_FUNC 'phre_vl' 349 !End of the abilint section 350 351 real(dp):: phre_vl 352 integer :: ig1,ig2,ig3,ia 353 phre_vl=phr_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),& 354 & ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia)) 355 end function phre_vl 356 357 function phimag_vl(ig1,ig2,ig3,ia) 358 359 360 !This section has been created automatically by the script Abilint (TD). 361 !Do not modify the following lines by hand. 362 #undef ABI_FUNC 363 #define ABI_FUNC 'phimag_vl' 364 !End of the abilint section 365 366 real(dp) :: phimag_vl 367 integer :: ig1,ig2,ig3,ia 368 phimag_vl=phi_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),& 369 & ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia)) 370 end function phimag_vl 371 372 function gsq_vl(i1,i2,i3) 373 374 375 !This section has been created automatically by the script Abilint (TD). 376 !Do not modify the following lines by hand. 377 #undef ABI_FUNC 378 #define ABI_FUNC 'gsq_vl' 379 !End of the abilint section 380 381 real(dp) :: gsq_vl 382 integer :: i1,i2,i3 383 !Define G^2 based on G space metric gmet. 384 gsq_vl=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 385 & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 386 & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 387 end function gsq_vl 388 389 function dgsqds_vl(i1,i2,i3) 390 391 392 !This section has been created automatically by the script Abilint (TD). 393 !Do not modify the following lines by hand. 394 #undef ABI_FUNC 395 #define ABI_FUNC 'dgsqds_vl' 396 !End of the abilint section 397 398 real(dp) :: dgsqds_vl 399 integer :: i1,i2,i3 400 !Define dG^2/ds based on G space metric derivative dgmetds. 401 dgsqds_vl=dble(i1*i1)*dgmetds(1,1)+dble(i2*i2)*dgmetds(2,2)+& 402 & dble(i3*i3)*dgmetds(3,3)+& 403 & dble(i1*i2)*(dgmetds(1,2)+dgmetds(2,1))+& 404 & dble(i1*i3)*(dgmetds(1,3)+dgmetds(3,1))+& 405 & dble(i2*i3)*(dgmetds(2,3)+dgmetds(3,2)) 406 end function dgsqds_vl 407 408 end subroutine vlocalstr