TABLE OF CONTENTS
ABINIT/vdw_dftd2 [ Functions ]
NAME
vdw_dftd2
FUNCTION
Compute energy and derivatives with respect to dimensionless reduced atom coordinates due to Van der Waals interaction. The formalism here follows the DFT-D2 approach of Grimme which consists in adding a semi-empirical dispersion potential (pair-wise force field) to the conventional Kohn-Sham DFT energy.
COPYRIGHT
Copyright (C) 2012-2018 ABINIT group (MT) 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
ixc=choice of exchange-correlation functional natom=number of atoms ntypat=number of atom types prtvol=printing volume (if >0, print computation parameters) typat(natom)=type integer for each atom in cell rprimd(3,3)=real space primitive translations vdw_tol=tolerance use to converge the potential (a pair of atoms is included in potential if its contribution is larger than vdw_tol) vdw_tol<0 takes default value (10^-10) xred(3,natom)=reduced atomic coordinates znucl(ntypat)=atomic number of atom type === optional inputs === [qphon(3)]=wavevector of the phonon; used only for dynamical matrix computation
OUTPUT
e_vdw_dftd2=contribution to energy from DFT-D2 dispersion potential === optional outputs === [dyn_vdw_dftd2(2,3,natom,3,natom)]=contribution to dynamical matrix from DFT-D2 dispersion potential [elt_vdw_dftd2(6+3*natom,6)]=contribution to elastic tensor and internal strains from DFT-D2 disp. pot. [fred_vdw_dftd2(3,natom)]=contribution to forces from DFT-D2 dispersion potential [str_vdw_dftd2(6)]=contribution to stress tensor from DFT-D2 dispersion potential
NOTES
Ref.: S. Grimme, Semiempirical GGA-type density functional
constructed with a long-range dispersion correction,
J. Comp. Chem. 27, 1787 (2006)
PARENTS
respfn,setvtr,stress
CHILDREN
SOURCE
57 #if defined HAVE_CONFIG_H 58 #include "config.h" 59 #endif 60 61 #include "abi_common.h" 62 63 subroutine vdw_dftd2(e_vdw_dftd2,ixc,natom,ntypat,prtvol,typat,rprimd,vdw_tol,xred,znucl,& 64 & dyn_vdw_dftd2,elt_vdw_dftd2,fred_vdw_dftd2,str_vdw_dftd2,qphon) ! Optionals 65 66 use defs_basis 67 use m_profiling_abi 68 use m_errors 69 use m_atomdata 70 71 use m_geometry, only : metric 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 'vdw_dftd2' 77 use interfaces_14_hidewrite 78 !End of the abilint section 79 80 implicit none 81 82 !Arguments ------------------------------------ 83 !scalars 84 integer,intent(in) :: ixc,natom,ntypat,prtvol 85 real(dp),intent(in) :: vdw_tol 86 real(dp),intent(out) :: e_vdw_dftd2 87 !arrays 88 integer,intent(in) :: typat(natom) 89 real(dp),intent(in) :: rprimd(3,3),xred(3,natom),znucl(ntypat) 90 real(dp),intent(in),optional :: qphon(3) 91 real(dp),intent(out),optional :: dyn_vdw_dftd2(2,3,natom,3,natom) 92 real(dp),intent(out),optional :: elt_vdw_dftd2(6+3*natom,6) 93 real(dp),intent(out),optional :: fred_vdw_dftd2(3,natom) 94 real(dp),intent(out),optional :: str_vdw_dftd2(6) 95 96 !Local variables------------------------------- 97 !scalars 98 integer,parameter :: vdw_nspecies=55 99 integer :: ia,ia1,ii,is1,is2,is3,itypat,ja,ja1,jj,jtypat,kk,ll,mu,npairs,nshell,nu 100 logical :: need_dynmat,need_elast,need_forces,need_intstr,need_stress 101 logical :: need_gradient,need_gradient2,newshell,qeq0=.true. 102 real(dp),parameter :: e_conv=(10/Bohr_Ang)**6/Ha_J/Avogadro ! 1 J.nm^6.mol^-1 in Ha.Bohr^6 103 real(dp),parameter :: vdw_d=20._dp,vdw_tol_default=tol10 104 real(dp),parameter :: vdw_s_pbe=0.75_dp, vdw_s_blyp=1.2_dp, vdw_s_b3lyp=1.05_dp 105 real(dp),parameter :: vdw_s_bp86=1.05_dp, vdw_s_tpss=1.0_dp, vdw_s_b97d=1.25_dp 106 real(dp) :: c6,c6r6,ex,fr,gr,gr2,grad,grad2,ph,ph1r,ph1i 107 real(dp) :: r0,r1,r2,r3,rcut,rcut2,rsq,rr,sfact,ucvol,vdw_s 108 character(len=500) :: msg 109 type(atomdata_t) :: atom 110 !arrays 111 integer,allocatable :: ivdw(:) 112 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/) 113 real(dp) :: gmet(3,3),gprimd(3,3),mat(3,3),rcart(3),rmet(3,3),vec(3) 114 real(dp),allocatable :: vdw_c6(:,:),vdw_r0(:,:),xred01(:,:) 115 real(dp),parameter :: vdw_c6_dftd2(vdw_nspecies)= & 116 & (/ 0.14, 0.08, 1.61, 1.61, 3.13, 1.75, 1.23, 0.70, 0.75, 0.63,& 117 & 5.71, 5.71,10.79, 9.23, 7.84, 5.57, 5.07, 4.61,10.80,10.80,& 118 & 10.80,10.80,10.80,10.80,10.80,10.80,10.80,10.80,10.80,10.80,& 119 & 16.99,17.10,16.37,12.64,12.47,12.01,24.67,24.67,24.67,24.67,& 120 & 24.67,24.67,24.67,24.67,24.67,24.67,24.67,24.67,37.32,38.71,& 121 & 38.44,31.74,31.50,29.99, 0.00/) 122 real(dp),parameter :: vdw_r0_dftd2(vdw_nspecies)= & 123 & (/1.001,1.012,0.825,1.408,1.485,1.452,1.397,1.342,1.287,1.243,& 124 & 1.144,1.364,1.639,1.716,1.705,1.683,1.639,1.595,1.485,1.474,& 125 & 1.562,1.562,1.562,1.562,1.562,1.562,1.562,1.562,1.562,1.562,& 126 & 1.650,1.727,1.760,1.771,1.749,1.727,1.628,1.606,1.639,1.639,& 127 & 1.639,1.639,1.639,1.639,1.639,1.639,1.639,1.639,1.672,1.804,& 128 & 1.881,1.892,1.892,1.881,1.000/) 129 character(len=2),parameter :: vdw_symb(vdw_nspecies)= & 130 & (/' H','He','Li','Be',' B',' C',' N',' O',' F','Ne',& 131 & 'Na','Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca',& 132 & 'Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn',& 133 & 'Ga','Ge','As','Se','Br','Kr','Rb','Sr',' Y','Zr',& 134 & 'Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn',& 135 & 'Sb','Te',' I','Xe','no'/) 136 137 ! ************************************************************************* 138 139 DBG_ENTER("COLL") 140 141 !Extract options 142 need_forces=present(fred_vdw_dftd2) 143 need_stress=present(str_vdw_dftd2) 144 need_dynmat=present(dyn_vdw_dftd2) 145 need_elast=present(elt_vdw_dftd2) 146 need_intstr=present(elt_vdw_dftd2) 147 need_gradient=(need_forces.or.need_stress) 148 need_gradient2=(need_dynmat.or.need_elast.or.need_intstr) 149 if (need_dynmat) then 150 if (.not.present(qphon)) then 151 msg='Dynamical matrix required without a q-vector' 152 MSG_BUG(msg) 153 end if 154 qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) 155 end if 156 157 !Identify type(s) of atoms 158 ABI_ALLOCATE(ivdw,(ntypat)) 159 do itypat=1,ntypat 160 ivdw(itypat)=-1;jtypat=0 161 call atomdata_from_znucl(atom,znucl(itypat)) 162 do while ((ivdw(itypat)<0).and.(jtypat<vdw_nspecies)) 163 jtypat=jtypat+1;if (vdw_symb(jtypat)==atom%symbol) ivdw(itypat)=jtypat 164 end do 165 if (ivdw(itypat)<0) then 166 write(msg,'(3a)') & 167 & 'Van der Waals DFT-D2 correction not available for atom type: ',atom%symbol,' !' 168 MSG_ERROR(msg) 169 end if 170 end do 171 172 !Select DFT-D2 VdW parameters according to system data 173 vdw_s=e_conv 174 if (ixc==11.or.ixc==-101130.or.ixc==-130101) then 175 vdw_s=vdw_s*vdw_s_pbe 176 else if (ixc==18.or.ixc==-106131.or.ixc==-131106) then 177 vdw_s=vdw_s*vdw_s_blyp 178 else if (ixc==19.or.ixc==-106132.or.ixc==-132106) then 179 vdw_s=vdw_s*vdw_s_bp86 180 else if (ixc==-202231.or.ixc==-231202) then 181 vdw_s=vdw_s*vdw_s_tpss 182 else 183 write(msg,'(a,i8,a)')' Van der Waals DFT-D2 correction not compatible with ixc=',ixc,' !' 184 MSG_ERROR(msg) 185 end if 186 ABI_ALLOCATE(vdw_c6,(ntypat,ntypat)) 187 ABI_ALLOCATE(vdw_r0,(ntypat,ntypat)) 188 do itypat=1,ntypat 189 do jtypat=1,ntypat 190 vdw_c6(itypat,jtypat)=sqrt(vdw_c6_dftd2(ivdw(itypat))*vdw_c6_dftd2(ivdw(jtypat))) 191 vdw_r0(itypat,jtypat)=(vdw_r0_dftd2(ivdw(itypat))+vdw_r0_dftd2(ivdw(jtypat)))/Bohr_Ang 192 end do 193 end do 194 195 !Computation of cut-off radius according to tolerance 196 !We take: r_cut=(s6*max(C6)/tol)**(1/6) and rcut<=75 bohr 197 if (vdw_tol<zero) then 198 rcut=(vdw_s/vdw_tol_default*maxval(vdw_c6))**sixth 199 else 200 rcut=(vdw_s/vdw_tol*maxval(vdw_c6))**sixth 201 end if 202 !rcut=min(rcut,100._dp) 203 rcut2=rcut*rcut 204 205 !Retrieve cell geometry data 206 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 207 208 !Map reduced coordinates into [0,1] 209 ABI_ALLOCATE(xred01,(3,natom)) 210 do ia=1,natom 211 xred01(1,ia)=xred(1,ia)-aint(xred(1,ia))+half-sign(half,xred(1,ia)) 212 xred01(2,ia)=xred(2,ia)-aint(xred(2,ia))+half-sign(half,xred(2,ia)) 213 xred01(3,ia)=xred(3,ia)-aint(xred(3,ia))+half-sign(half,xred(3,ia)) 214 end do 215 216 !Set accumulated quantities to zero 217 npairs=0 218 e_vdw_dftd2=zero 219 if (need_forces) fred_vdw_dftd2=zero 220 if (need_stress) str_vdw_dftd2=zero 221 if (need_dynmat) dyn_vdw_dftd2=zero 222 if (need_elast) elt_vdw_dftd2(1:6,1:6)=zero 223 if (need_intstr) elt_vdw_dftd2(7:6+3*natom,1:6)=zero 224 225 !Loop over shells of cell replicas 226 nshell=0 227 do 228 newshell=.false.;nshell=nshell+1 229 230 ! Loop over cell replicas in the shell 231 ! ns1=1+int(rcut*sqrt(SUM(gprimd(:,1)**2)) 232 ! ns2=1+int(rcut*sqrt(SUM(gprimd(:,2)**2)) 233 ! ns3=1+int(rcut*sqrt(SUM(gprimd(:,3)**2)) 234 do is3=-nshell,nshell 235 do is2=-nshell,nshell 236 do is1=-nshell,nshell 237 if (nshell==1.or. & 238 & abs(is3)==nshell.or.abs(is2)==nshell.or.abs(is1)==nshell) then 239 240 ! Phase for dynamical matrix 241 if (need_dynmat) then 242 ph1r=one;ph1i=zero !ph1=exp(-iqL) 243 if (.not.qeq0) then 244 ph=-two_pi*(qphon(1)*is1+qphon(2)*is2+qphon(3)*is3) 245 ph1r=cos(ph);ph1i=sin(ph) 246 end if 247 end if 248 249 ! Loops over atoms a and b 250 do ia=1,natom 251 itypat=typat(ia) 252 do ja=1,ia 253 jtypat=typat(ja) 254 r1=xred01(1,ia)-xred01(1,ja)-dble(is1) 255 r2=xred01(2,ia)-xred01(2,ja)-dble(is2) 256 r3=xred01(3,ia)-xred01(3,ja)-dble(is3) 257 rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3 & 258 & +two*(rmet(2,1)*r2*r1+rmet(3,2)*r3*r2+rmet(3,1)*r1*r3) 259 260 ! Select atomic pairs (a,b) and avoid atom_a=atom_b 261 if (rsq>=tol16.and.rsq<=rcut2) then 262 263 ! Data for the selected pair 264 npairs=npairs+1;newshell=.true. 265 sfact=vdw_s;if (ia==ja) sfact=half*sfact 266 rr=sqrt(rsq) 267 c6=vdw_c6(itypat,jtypat) 268 r0=vdw_r0(itypat,jtypat) 269 270 ! Computation of pair-wise potential 271 ex=exp(-vdw_d*(rr/r0-one)) 272 fr=one/(one+ex) 273 c6r6=c6/rr**6 274 275 ! Contribution to energy 276 e_vdw_dftd2=e_vdw_dftd2-sfact*fr*c6r6 277 278 if (need_gradient.or.need_gradient2) then 279 gr=(vdw_d/r0)*(fr**2)*ex 280 grad=-sfact*(gr-six*fr/rr)*c6r6/rr 281 rcart(1)=rprimd(1,1)*r1+rprimd(1,2)*r2+rprimd(1,3)*r3 282 rcart(2)=rprimd(2,1)*r1+rprimd(2,2)*r2+rprimd(2,3)*r3 283 rcart(3)=rprimd(3,1)*r1+rprimd(3,2)*r2+rprimd(3,3)*r3 284 285 ! Contribution to forces 286 if (need_forces.and.ia/=ja) then 287 vec(1:3)=grad*rcart(1:3) 288 fred_vdw_dftd2(1:3,ia)=fred_vdw_dftd2(1:3,ia)+vec(1:3) 289 fred_vdw_dftd2(1:3,ja)=fred_vdw_dftd2(1:3,ja)-vec(1:3) 290 end if 291 292 ! Contribution to stress tensor 293 if (need_stress) then 294 do mu=1,6 295 ii=alpha(mu);jj=beta(mu) 296 str_vdw_dftd2(mu)=str_vdw_dftd2(mu)+grad*rcart(ii)*rcart(jj) 297 end do 298 end if 299 300 if (need_gradient2) then 301 gr2=(vdw_d/r0)*gr*(2*fr*ex-one) 302 grad2=-sfact*(gr2-13._dp*gr/rr+48._dp*fr/rr**2)*c6r6/rr**2 303 304 ! Contribution to dynamical matrix (phase factors are subtle!) 305 if (need_dynmat) then 306 mat(1:3,1)=grad2*rcart(1:3)*rcart(1) ; mat(1,1)=mat(1,1)+grad 307 mat(1:3,2)=grad2*rcart(1:3)*rcart(2) ; mat(2,2)=mat(2,2)+grad 308 mat(1:3,3)=grad2*rcart(1:3)*rcart(3) ; mat(3,3)=mat(3,3)+grad 309 if (ia/=ja) then 310 do ii=1,3 311 dyn_vdw_dftd2(1,1:3,ia,ii,ia)=dyn_vdw_dftd2(1,1:3,ia,ii,ia)+mat(1:3,ii) 312 dyn_vdw_dftd2(1,1:3,ja,ii,ja)=dyn_vdw_dftd2(1,1:3,ja,ii,ja)+mat(1:3,ii) 313 dyn_vdw_dftd2(1,1:3,ia,ii,ja)=dyn_vdw_dftd2(1,1:3,ia,ii,ja)-mat(1:3,ii)*ph1r 314 dyn_vdw_dftd2(2,1:3,ia,ii,ja)=dyn_vdw_dftd2(2,1:3,ia,ii,ja)-mat(1:3,ii)*ph1i 315 dyn_vdw_dftd2(1,1:3,ja,ii,ia)=dyn_vdw_dftd2(1,1:3,ja,ii,ia)-mat(1:3,ii)*ph1r 316 dyn_vdw_dftd2(2,1:3,ja,ii,ia)=dyn_vdw_dftd2(2,1:3,ja,ii,ia)+mat(1:3,ii)*ph1i 317 end do 318 else if (.not.qeq0) then 319 do ii=1,3 320 dyn_vdw_dftd2(1,1:3,ia,ii,ia)=dyn_vdw_dftd2(1,1:3,ia,ii,ia) & 321 & +two*mat(1:3,ii)*(one-ph1r) 322 end do 323 end if 324 end if 325 326 ! Contribution to elastic tensor 327 if (need_elast) then 328 do mu=1,6 329 ii=alpha(mu);jj=beta(mu) 330 do nu=1,6 331 kk=alpha(nu);ll=beta(nu) 332 elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) & 333 & +grad2*rcart(ii)*rcart(jj)*rcart(kk)*rcart(ll) 334 if (ii==kk) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) & 335 & +half*grad*rcart(jj)*rcart(ll) 336 if (ii==ll) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) & 337 & +half*grad*rcart(jj)*rcart(kk) 338 if (jj==kk) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) & 339 & +half*grad*rcart(ii)*rcart(ll) 340 if (jj==ll) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) & 341 & +half*grad*rcart(ii)*rcart(kk) 342 end do 343 end do 344 end if 345 346 ! Contribution to internal strains 347 if (need_intstr.and.ia/=ja) then 348 ia1=6+3*(ia-1);ja1=6+3*(ja-1) 349 do mu=1,6 350 ii=alpha(mu);jj=beta(mu) 351 vec(1:3)=grad2*rcart(ii)*rcart(jj)*rcart(1:3) 352 vec(ii)=vec(ii)+half*grad*rcart(jj) 353 vec(jj)=vec(jj)+half*grad*rcart(ii) 354 elt_vdw_dftd2(ia1+1:ia1+3,mu)=elt_vdw_dftd2(ia1+1:ia1+3,mu)+vec(1:3) 355 elt_vdw_dftd2(ja1+1:ja1+3,mu)=elt_vdw_dftd2(ja1+1:ja1+3,mu)-vec(1:3) 356 end do 357 end if 358 359 end if ! Computation of 2nd gradient 360 end if ! Computation of gradient 361 end if ! Pairs selection 362 end do ! Loop over atom b 363 end do ! Loop over atom a 364 end if ! Triple loop over cell replicas in shell 365 end do 366 end do 367 end do 368 if(.not.newshell) exit ! Check if new shell must be calculated 369 end do ! Loop over shells 370 371 !Gradients: convert them from cartesian to reduced coordinates 372 if (need_forces) then 373 do ia=1,natom 374 call grad_cart2red(fred_vdw_dftd2(:,ia)) 375 end do 376 end if 377 if (need_dynmat) then 378 do ja=1,natom 379 do ia=1,natom 380 do kk=1,merge(2,1,qeq0) 381 do ii=1,3 382 vec(1:3)=dyn_vdw_dftd2(kk,1:3,ia,ii,ja) 383 call grad_cart2red(vec) 384 dyn_vdw_dftd2(kk,1:3,ia,ii,ja)=vec(1:3) 385 end do 386 do ii=1,3 387 vec(1:3)=dyn_vdw_dftd2(kk,ii,ia,1:3,ja) 388 call grad_cart2red(vec) 389 dyn_vdw_dftd2(kk,ii,ia,1:3,ja)=vec(1:3) 390 end do 391 end do 392 end do 393 end do 394 end if 395 if (need_intstr) then 396 do mu=1,6 397 ia1=6 398 do ia=1,natom 399 call grad_cart2red(elt_vdw_dftd2(ia1+1:ia1+3,mu)) 400 ia1=ia1+3 401 end do 402 end do 403 end if 404 405 !DEBUG 406 !write(77,*) "---------------" 407 !write(77,*) "E=",e_vdw_dftd2 408 !if (need_forces) then 409 ! do ia=1,natom 410 ! write(77,*) "F=",ia,fred_vdw_dftd2(:,ia) 411 ! end do 412 !end if 413 !if (need_stress) write(77,*) "S=",str_vdw_dftd2(:) 414 !if (need_dynmat) then 415 ! do ia=1,natom 416 ! do ii=1,3 417 ! do ja=1,natom 418 ! write(77,*) "D=",ia,ii,ja,dyn_vdw_dftd2(:,:,ja,ii,ia) 419 ! end do 420 ! end do 421 ! end do 422 !end if 423 !if (need_elast) then 424 ! do ii=1,6 425 ! write(77,*) "e=",ii,elt_vdw_dftd2(1:6,ii) 426 ! end do 427 !end if 428 !if (need_intstr) then 429 ! do ii=1,6 430 ! do ia=1,natom 431 ! write(77,*) "I=",ii,ia,elt_vdw_dftd2(7+3*(ia-1):9+3*(ia-1),ii) 432 ! end do 433 ! end do 434 !end if 435 !flush(77) 436 !DEBUG 437 438 !Stress tensor: divide by volume 439 if (need_stress) str_vdw_dftd2=str_vdw_dftd2/ucvol 440 441 !Printing 442 if (prtvol>0) then 443 write(msg,'(10a)') ch10,& 444 & ' --------------------------------------------------------------',ch10,& 445 & ' Van der Waals DFT-D2 semi-empirical dispersion potential added',ch10,& 446 & ' with following parameters:',ch10,& 447 & ' Specie C6 (J.nm^6.mol^-1) R0 (Ang)',ch10,& 448 & ' ------------------------------------' 449 call wrtout(std_out,msg,'COLL') 450 do itypat=1,ntypat 451 write(msg,'(9X,a2,11X,f5.2,8X,f6.3)') & 452 & vdw_symb(ivdw(itypat)),vdw_c6_dftd2(ivdw(itypat)),vdw_r0_dftd2(ivdw(itypat)) 453 call wrtout(std_out,msg,'COLL') 454 end do 455 write(msg,'(2a,f6.2,2a,f6.2,2a,f6.2,a)') ch10,& 456 & ' Scaling factor = ',vdw_s/e_conv,ch10,& 457 & ' Damping parameter= ',vdw_d,ch10,& 458 & ' Cut-off radius = ',rcut,' bohr' 459 call wrtout(std_out,msg,'COLL') 460 write(msg,'(2a,i8,2a,es14.5,4a)') ch10,& 461 & ' Number of pairs contributing = ',npairs,ch10,& 462 & ' DFT-D2 energy contribution = ',e_vdw_dftd2,' Ha',ch10,& 463 & ' --------------------------------------------------------------',ch10 464 call wrtout(std_out,msg,'COLL') 465 end if 466 467 ABI_DEALLOCATE(ivdw) 468 ABI_DEALLOCATE(vdw_c6) 469 ABI_DEALLOCATE(vdw_r0) 470 ABI_DEALLOCATE(xred01) 471 472 DBG_EXIT("COLL") 473 474 contains
vdw_dftd2/grad_cart2red [ Functions ]
[ Top ] [ vdw_dftd2 ] [ Functions ]
NAME
grad_cart2red
FUNCTION
Convert gradients from cartesian to reduced coordinates
PARENTS
vdw_dftd2
CHILDREN
SOURCE
492 subroutine grad_cart2red(grad) 493 494 495 !This section has been created automatically by the script Abilint (TD). 496 !Do not modify the following lines by hand. 497 #undef ABI_FUNC 498 #define ABI_FUNC 'grad_cart2red' 499 !End of the abilint section 500 501 implicit none 502 503 !Arguments ------------------------------------ 504 real(dp),intent(inout) :: grad(3) 505 !Local variables------------------------------- 506 real(dp) :: tmp(3) 507 508 ! ********************************************************************* 509 510 tmp(1)=rprimd(1,1)*grad(1)+rprimd(2,1)*grad(2)+rprimd(3,1)*grad(3) 511 tmp(2)=rprimd(1,2)*grad(1)+rprimd(2,2)*grad(2)+rprimd(3,2)*grad(3) 512 tmp(3)=rprimd(1,3)*grad(1)+rprimd(2,3)*grad(2)+rprimd(3,3)*grad(3) 513 grad(1:3)=tmp(1:3) 514 515 end subroutine grad_cart2red