TABLE OF CONTENTS
ABINIT/dfpt_vlocal [ Functions ]
NAME
dfpt_vlocal
FUNCTION
Compute local part of 1st-order potential from the appropriate atomic pseudopotential with structure and derivative factor. In case of derivative with respect to k or electric (magnetic Zeeman) field perturbation, the 1st-order local potential vanishes.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (XG,MM) 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) cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX gmet(3,3)=reciprocal space metric (Bohr**-2) gsqcut=cutoff G**2 for included G s in fft box. idir=direction of atomic displacement (=1,2 or 3 : displacement of atom ipert along the 1st, 2nd or 3rd axis). ipert=number of the atom being displaced in the frozen-phonon mpi_enreg=information about MPI parallelization mqgrid=dimension of q grid for pseudopotentials natom=number of atoms in 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 in cell. n1,n2,n3=fft grid. ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information. qgrid(mqgrid)=grid of q points from 0 to qmax. qphon(3)=wavevector of the phonon ucvol=unit cell volume (Bohr**3). vlspl(mqgrid,2,ntypat)=spline fit of q^2 V(q) for each type of atom. xred(3,natom)=reduced atomic coordinates
OUTPUT
vpsp1(cplex*nfft)=first-order local crystal pseudopotential in real space (including the minus sign, forgotten in the paper non-linear..
PARENTS
dfpt_looppert,dfpt_nstdy,dfpt_nstpaw,dfptnl_loop
CHILDREN
fourdp,ptabs_fourdp
SOURCE
55 #if defined HAVE_CONFIG_H 56 #include "config.h" 57 #endif 58 59 #include "abi_common.h" 60 61 62 subroutine dfpt_vlocal(atindx,cplex,gmet,gsqcut,idir,ipert,& 63 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,& 64 & ntypat,n1,n2,n3,paral_kgb,ph1d,qgrid,qphon,ucvol,vlspl,vpsp1,xred) 65 66 use defs_basis 67 use defs_abitypes 68 use m_errors 69 use m_profiling_abi 70 71 use m_mpinfo, only : ptabs_fourdp 72 73 !This section has been created automatically by the script Abilint (TD). 74 !Do not modify the following lines by hand. 75 #undef ABI_FUNC 76 #define ABI_FUNC 'dfpt_vlocal' 77 use interfaces_53_ffts 78 !End of the abilint section 79 80 implicit none 81 82 !Arguments ------------------------------- 83 !scalars 84 integer,intent(in) :: cplex,idir,ipert,mqgrid,n1,n2,n3,natom,nfft,ntypat 85 integer,intent(in) :: paral_kgb 86 real(dp),intent(in) :: gsqcut,ucvol 87 type(MPI_type),intent(in) :: mpi_enreg 88 !arrays 89 integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18) 90 real(dp),intent(in) :: gmet(3,3),ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom) 91 real(dp),intent(in) :: qgrid(mqgrid),qphon(3),vlspl(mqgrid,2,ntypat) 92 real(dp),intent(in) :: xred(3,natom) 93 real(dp),intent(out) :: vpsp1(cplex*nfft) 94 95 !Local variables ------------------------- 96 !scalars 97 integer :: i1,i2,i3,ia1,iatom,id1,id2,id3,ig1,ig2,ig3,ii,ii1,im=2 98 integer :: itypat,jj,re=1 99 real(dp),parameter :: tolfix=1.000000001_dp 100 real(dp) :: aa,bb,cc,cutoff,dd,diff,dq,dq2div6,dqdiv6,dqm1,gmag,gq1 101 real(dp) :: gq2,gq3,gsquar,phqim,phqre 102 real(dp) :: qxred2pi,sfi,sfr,vion1,xnorm 103 logical :: qeq0 104 !arrays 105 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 106 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 107 real(dp) :: gq(3) 108 real(dp),allocatable :: work1(:,:) 109 110 ! ********************************************************************* 111 112 iatom=ipert 113 114 if(iatom==natom+1 .or. iatom==natom+2 .or. iatom==natom+10 .or. iatom==natom+11 .or. iatom==natom+5)then 115 116 ! (In case of d/dk or an electric field, or magnetic (Zeeman) field->[natom+5] SPr deb ) 117 vpsp1(1:cplex*nfft)=zero 118 119 else 120 121 ! (In case of a phonon perturbation) 122 ABI_ALLOCATE(work1,(2,nfft)) 123 work1(1:2,1:nfft)=0.0_dp 124 125 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 126 dqm1=1.0_dp/dq 127 dqdiv6=dq/6.0_dp 128 dq2div6=dq**2/6.0_dp 129 cutoff=gsqcut*tolfix 130 id1=n1/2+2 131 id2=n2/2+2 132 id3=n3/2+2 133 134 ! Get the distrib associated with this fft_grid 135 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 136 137 ! This is to allow q=0 138 qeq0=.false. 139 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)qeq0=.true. 140 141 ! Determination of the atom type 142 ia1=0 143 itypat=0 144 do ii=1,ntypat 145 ia1=ia1+nattyp(ii) 146 if(atindx(iatom)<=ia1.and.itypat==0)itypat=ii 147 end do 148 149 ! Determination of phase qxred* 150 qxred2pi=2.0_dp*pi*(qphon(1)*xred(1,iatom)+ & 151 & qphon(2)*xred(2,iatom)+ & 152 & qphon(3)*xred(3,iatom) ) 153 phqre=cos(qxred2pi) 154 phqim=sin(qxred2pi) 155 ii=0 156 157 do i3=1,n3 158 ig3=i3-(i3/id3)*n3-1 159 gq3=dble(ig3)+qphon(3) 160 gq(3)=gq3 161 do i2=1,n2 162 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 163 ig2=i2-(i2/id2)*n2-1 164 gq2=dble(ig2)+qphon(2) 165 gq(2)=gq2 166 167 ! Note the lower limit of the next loop 168 ii1=1 169 if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then 170 ii1=2 171 ii=ii+1 172 end if 173 do i1=ii1,n1 174 ig1=i1-(i1/id1)*n1-1 175 gq1=dble(ig1)+qphon(1) 176 gq(1)=gq1 177 ii=ii+1 178 gsquar=gsq_vl3(gq1,gq2,gq3) 179 ! Skip G**2 outside cutoff: 180 if (gsquar<=cutoff) then 181 gmag=sqrt(gsquar) 182 183 ! Compute vion(G) for given type of atom 184 jj=1+int(gmag*dqm1) 185 diff=gmag-qgrid(jj) 186 187 ! Evaluate spline fit from q^2 V(q) to get V(q): 188 ! (p. 86 Numerical Recipes, Press et al; NOTE error in book for sign 189 ! of "aa" term in derivative; also see splfit routine. 190 ! This bug fixed here 27 Jan 1992.) 191 192 bb = diff*dqm1 193 aa = 1.0_dp-bb 194 cc = aa*(aa**2-1.0_dp)*dq2div6 195 dd = bb*(bb**2-1.0_dp)*dq2div6 196 vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) + & 197 & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) & 198 & / gsquar 199 200 ! Phase G*xred (complex conjugate) * -i *2pi*(g+q)*vion 201 sfr=-phimag_vl3(ig1,ig2,ig3,iatom)*2.0_dp*pi*gq(idir)*vion1 202 sfi=-phre_vl3(ig1,ig2,ig3,iatom)*2.0_dp*pi*gq(idir)*vion1 203 204 ! Phase q*xred (complex conjugate) 205 work1(re,ii)=sfr*phqre+sfi*phqim 206 work1(im,ii)=-sfr*phqim+sfi*phqre 207 end if 208 209 end do 210 end if 211 end do 212 end do 213 214 ! Transform back to real space 215 call fourdp(cplex,work1,vpsp1,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 216 217 xnorm=1.0_dp/ucvol 218 vpsp1(1:cplex*nfft)=vpsp1(1:cplex*nfft)*xnorm 219 220 ABI_DEALLOCATE(work1) 221 222 ! End the condition of non-electric-field 223 end if 224 225 contains 226 227 !Real and imaginary parts of phase. 228 function phr_vl3(x1,y1,x2,y2,x3,y3) 229 230 231 !This section has been created automatically by the script Abilint (TD). 232 !Do not modify the following lines by hand. 233 #undef ABI_FUNC 234 #define ABI_FUNC 'phr_vl3' 235 !End of the abilint section 236 237 real(dp) :: phr_vl3 238 real(dp),intent(in) :: x1,x2,x3,y1,y2,y3 239 phr_vl3=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 240 end function phr_vl3 241 242 function phi_vl3(x1,y1,x2,y2,x3,y3) 243 244 245 !This section has been created automatically by the script Abilint (TD). 246 !Do not modify the following lines by hand. 247 #undef ABI_FUNC 248 #define ABI_FUNC 'phi_vl3' 249 !End of the abilint section 250 251 real(dp) :: phi_vl3 252 real(dp),intent(in) :: x1,x2,x3,y1,y2,y3 253 phi_vl3=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 254 end function phi_vl3 255 256 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 257 function ph1_vl3(nri,ig1,ia) 258 259 260 !This section has been created automatically by the script Abilint (TD). 261 !Do not modify the following lines by hand. 262 #undef ABI_FUNC 263 #define ABI_FUNC 'ph1_vl3' 264 !End of the abilint section 265 266 real(dp) :: ph1_vl3 267 integer,intent(in) :: nri,ig1,ia 268 ph1_vl3=ph1d(nri,ig1+1+n1+(atindx(ia)-1)*(2*n1+1)) 269 end function ph1_vl3 270 271 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 272 function ph2_vl3(nri,ig2,ia) 273 274 275 !This section has been created automatically by the script Abilint (TD). 276 !Do not modify the following lines by hand. 277 #undef ABI_FUNC 278 #define ABI_FUNC 'ph2_vl3' 279 !End of the abilint section 280 281 real(dp) :: ph2_vl3 282 integer,intent(in) :: nri,ig2,ia 283 ph2_vl3=ph1d(nri,ig2+1+n2+(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1)) 284 end function ph2_vl3 285 286 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 287 function ph3_vl3(nri,ig3,ia) 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 'ph3_vl3' 294 !End of the abilint section 295 296 real(dp) :: ph3_vl3 297 integer,intent(in) :: nri,ig3,ia 298 ph3_vl3=ph1d(nri,ig3+1+n3+(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 299 end function ph3_vl3 300 301 function phre_vl3(ig1,ig2,ig3,ia) 302 303 304 !This section has been created automatically by the script Abilint (TD). 305 !Do not modify the following lines by hand. 306 #undef ABI_FUNC 307 #define ABI_FUNC 'phre_vl3' 308 !End of the abilint section 309 310 real(dp) :: phre_vl3 311 integer,intent(in) :: ig1,ig2,ig3,ia 312 phre_vl3=phr_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),& 313 & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia)) 314 end function phre_vl3 315 316 function phimag_vl3(ig1,ig2,ig3,ia) 317 318 319 !This section has been created automatically by the script Abilint (TD). 320 !Do not modify the following lines by hand. 321 #undef ABI_FUNC 322 #define ABI_FUNC 'phimag_vl3' 323 !End of the abilint section 324 325 real(dp) :: phimag_vl3 326 integer,intent(in) :: ig1,ig2,ig3,ia 327 phimag_vl3=phi_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),& 328 & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia)) 329 end function phimag_vl3 330 331 function gsq_vl3(g1,g2,g3) 332 333 334 !This section has been created automatically by the script Abilint (TD). 335 !Do not modify the following lines by hand. 336 #undef ABI_FUNC 337 #define ABI_FUNC 'gsq_vl3' 338 !End of the abilint section 339 340 real(dp) :: gsq_vl3 341 real(dp),intent(in) :: g1,g2,g3 ! Note that they are real, unlike in other similar function definitions 342 !Define G^2 based on G space metric gmet. 343 gsq_vl3=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+& 344 & g3*g3*gmet(3,3)+2.0_dp*g1*g2*gmet(1,2)+& 345 & 2.0_dp*g2*g3*gmet(2,3)+2.0_dp*g3*g1*gmet(3,1) 346 end function gsq_vl3 347 348 end subroutine dfpt_vlocal