TABLE OF CONTENTS
ABINIT/vdw_dftd3 [ Functions ]
NAME
vdw_dftd3
FUNCTION
Compute energy, forces, stress, interatomic force constant and elastic contribution due to dispersion interaction as formulated by Grimme in the DFT-D3 approach. The last cited adds a dispersion potential (pair-wise force field, rij^6 and rij^8) to Kohn-Sham DFT energy. It is also possible to include a three-body term and molecular dispersion (using vdw_tol_3bt>0). DFT-D3(Becke and Johnson), another formulation which avoids the use of a damping function to remove the undesired short-range behaviour is also activable using vdw_xc=7
COPYRIGHT
Copyright (C) 2015-2018 ABINIT group (BVT) 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 vdw_xc= select van-der-Waals correction if =6: DFT-D3 as in Grimme, J. Chem. Phys. 132, 154104 (2010) if =7: DFT-D3(BJ) as in Grimme, Comput. Chem. 32, 1456 (2011) Only the use of R0 = a1 C8/C6 + a2 is available here vdw_tol=tolerance use to converge the pair-wise 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) vdw_tol_3bt= tolerance use to converge three body terms (only for vdw_xc=6) a triplet of atom contributes to the correction if its contribution is larger than vdw_tol_3bt xred(3,natom)=reduced atomic coordinates znucl(ntypat)=atomic number of atom type === optional input === [qphon(3)]= reduced q-vector along which is computed the DFT-D3 contribution to the IFCs in reciprocal space
OUTPUT
e_vdw_dftd3=contribution to energy from DFT-D3 dispersion potential === optional outputs === [elt_vdw_dftd3(6+3*natom,6)]= contribution to elastic constant and internal strains from DFT-D3 dispersion potential [fred_vdw_dftd3(3,natom)]=contribution to gradient w.r.to atomic displ. from DFT-D3 dispersion potential [str_vdw_dftd3(6)]=contribution to stress tensor from DFT-D3 dispersion potential [dyn_vdw_dftd3(2,3,natom,3,natom)]= contribution to the interatomic force constants (in reciprocal space) at given input q-vector from DFT-D3 dispersion potential
NOTES
Ref.: DFT-D3: S. Grimme, J. Antony, S. Ehrlich, and H. Krieg A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu J. Chem. Phys. 132, 154104 (2010) DFT-D3(BJ) S. Grimme, S. Ehrlich and L. Goerigk Effect of the damping function in dispersion corrected density functional theory Comput. Chem. 32, 1456 (2011)
PARENTS
respfn,setvtr,stress
CHILDREN
SOURCE
76 #if defined HAVE_CONFIG_H 77 #include "config.h" 78 #endif 79 80 #include "abi_common.h" 81 82 subroutine vdw_dftd3(e_vdw_dftd3,ixc,natom,ntypat,prtvol,typat,rprimd,vdw_xc,& 83 & vdw_tol,vdw_tol_3bt,xred,znucl,dyn_vdw_dftd3,elt_vdw_dftd3,& 84 & fred_vdw_dftd3,str_vdw_dftd3,qphon) 85 86 use defs_basis 87 use m_profiling_abi 88 use m_errors 89 use m_atomdata 90 91 use m_special_funcs, only : abi_derfc 92 use m_geometry, only : metric 93 94 !This section has been created automatically by the script Abilint (TD). 95 !Do not modify the following lines by hand. 96 #undef ABI_FUNC 97 #define ABI_FUNC 'vdw_dftd3' 98 use interfaces_14_hidewrite 99 use interfaces_20_datashare 100 !End of the abilint section 101 102 implicit none 103 104 !Arguments ------------------------------------ 105 !scalars 106 integer,intent(in) :: ixc,natom,ntypat,prtvol,vdw_xc 107 real(dp),intent(in) :: vdw_tol,vdw_tol_3bt 108 real(dp),intent(out) :: e_vdw_dftd3 109 !arrays 110 integer,intent(in) :: typat(natom) 111 real(dp),intent(in) :: rprimd(3,3),xred(3,natom),znucl(ntypat) 112 real(dp),intent(in),optional :: qphon(3) 113 real(dp),intent(out),optional :: dyn_vdw_dftd3(2,3,natom,3,natom) 114 real(dp),intent(out),optional :: elt_vdw_dftd3(6+3*natom,6) 115 real(dp),intent(out),optional :: fred_vdw_dftd3(3,natom) 116 real(dp),intent(out),optional :: str_vdw_dftd3(6) 117 118 !Local variables------------------------------- 119 !scalars 120 ! The maximal number of reference systems for c6 is 5 (for C) 121 integer,parameter :: vdw_nspecies=94 122 integer:: alpha,beta,ia,ii,indi,indj,index_ia,index_ja,index_ka 123 integer :: is1,is2,is3,itypat,ja,jj,js1,js2,js3 124 integer :: jtypat,ka,kk,ktypat,la,ll 125 integer :: nline,npairs,nshell 126 integer :: refi,refj,refmax 127 logical :: bol_3bt,found 128 logical :: need_dynmat,need_elast,need_forces,need_grad,need_hess,need_stress,newshell 129 real(dp),parameter :: alpha6=14.0_dp, alpha8=16.0_dp 130 real(dp),parameter:: k1=16.0_dp, k2=15.0_dp, k3=4.0_dp 131 ! s8 parameters (BJ case) 132 real(dp),parameter :: vdwbj_s8_b2plyp=1.0860_dp, vdwbj_s8_pw6b95=0.7257_dp 133 real(dp),parameter :: vdwbj_s8_b97d=2.2609_dp, vdwbj_s8_revpbe=2.3550_dp 134 real(dp),parameter :: vdwbj_s8_b3lyp=1.9889_dp, vdwbj_s8_blyp=2.6996_dp 135 real(dp),parameter :: vdwbj_s8_tpss0=1.2576_dp, vdwbj_s8_pbe0=1.2177_dp 136 real(dp),parameter :: vdwbj_s8_tpss=1.9435_dp, vdwbj_s8_pbe=0.7875_dp 137 real(dp),parameter :: vdwbj_s8_bp86=3.2822_dp 138 ! a1 parameters (BJ only) 139 real(dp),parameter :: vdw_a1_b2lyp = 0.3451_dp, vdw_a1_pw6b95= 0.2076_dp 140 real(dp),parameter :: vdw_a1_b97d= 0.5545_dp, vdw_a1_revpbe= 0.5238_dp 141 real(dp),parameter :: vdw_a1_b3lyp = 0.3981_dp, vdw_a1_blyp = 0.4298_dp 142 real(dp),parameter :: vdw_a1_tpss0 = 0.3768_dp, vdw_a1_pbe0 = 0.4145_dp 143 real(dp),parameter :: vdw_a1_tpss = 0.4535_dp, vdw_a1_pbe = 0.4289_dp 144 real(dp),parameter :: vdw_a1_bp86 = 0.3946_dp 145 ! a2 parameters (BJ only) 146 real(dp),parameter :: vdw_a2_b2lyp = 4.7735_dp, vdw_a2_pw6b95= 6.3750_dp 147 real(dp),parameter :: vdw_a2_b97d= 3.2287_dp, vdw_a2_revpbe= 3.5016_dp 148 real(dp),parameter :: vdw_a2_b3lyp = 4.4211_dp, vdw_a2_blyp = 4.2359_dp 149 real(dp),parameter :: vdw_a2_tpss0 = 4.5865_dp, vdw_a2_pbe0 = 4.8593_dp 150 real(dp),parameter :: vdw_a2_tpss = 4.4752_dp, vdw_a2_pbe = 4.4407_dp 151 real(dp),parameter :: vdw_a2_bp86 = 4.8516_dp 152 ! s6=1 except for double hybrid functionals 153 real(dp),parameter :: vdw_s6_b2plyp=0.5_dp 154 ! s8 parameters 155 real(dp),parameter :: vdw_s8_b2plyp=1.000_dp, vdw_s8_pw6b95=0.862_dp 156 real(dp),parameter :: vdw_s8_b97d=0.909_dp, vdw_s8_revpbe=1.010_dp 157 real(dp),parameter :: vdw_s8_b3lyp=1.703_dp, vdw_s8_blyp=1.682_dp 158 real(dp),parameter :: vdw_s8_tpss0=1.242_dp, vdw_s8_pbe0=0.928_dp 159 real(dp),parameter :: vdw_s8_tpss=1.105_dp, vdw_s8_pbe=0.722_dp 160 real(dp),parameter :: vdw_s8_bp86=1.683_dp 161 ! sr6 parameters 162 real(dp),parameter :: vdw_sr6_b2plyp=1.332_dp, vdw_sr6_pw6b95=1.532_dp 163 real(dp),parameter :: vdw_sr6_b97d=0.892_dp, vdw_sr6_revpbe=0.923_dp 164 real(dp),parameter :: vdw_sr6_b3lyp=1.261_dp, vdw_sr6_blyp=1.094_dp 165 real(dp),parameter :: vdw_sr6_tpss0=1.252_dp, vdw_sr6_pbe0=1.287_dp 166 real(dp),parameter :: vdw_sr6_tpss=1.166_dp, vdw_sr6_pbe=1.217_dp 167 real(dp),parameter :: vdw_sr6_bp86=1.139_dp 168 ! sr8 parameters 169 real(dp),parameter :: vdw_sr8=one, vdw_sr9=3.0/4.0 170 real(dp),parameter :: vdw_tol_default=tol10 171 real(dp) :: ang,arg,cn_dmp,cosa,cosb,cosc,c6,c8 172 real(dp) :: dcosa_r3drij,dcosa_r3drjk,dcosa_r3drki 173 real(dp) :: dcosb_r3drij,dcosb_r3drjk,dcosb_r3drki 174 real(dp) :: dcosc_r3drij,dcosc_r3drjk,dcosc_r3drki 175 real(dp) :: dcn_dmp,dexp_cn,dfdmp,dfdmp_drij 176 real(dp) :: dfdmp_drjk,dfdmp_drki,dlri,dlrj,dmp,dmp6,dmp8,dmp9,dr,d2lri,d2lrj,d2lrirj 177 real(dp) :: dsysref,dsysref_a, dsysref_b 178 real(dp) :: d1_r3drij,d1_r3drjk,d1_r3drki,d2cn_dmp,d2cn_exp,d2frac_cn 179 real(dp) :: d_drij,d_drjk,d_drki 180 real(dp) :: exp_cn,e_no_c6,e_no_c8,e_3bt,fdmp6,fdmp8,fdmp9,frac_cn 181 real(dp) :: grad,grad_no_c,grad6,grad6_no_c6,grad8,grad8_no_c8,gr6,gr8 182 real(dp) :: hess,hessij, hess6, hess8,im_arg,l,ltot 183 real(dp) :: max_vdw_c6,min_dsys,re_arg,rcovij,rcut,rcutcn,rcut2,rcut9 184 real(dp) :: rsq,rsqij,rsqjk,rsqki,rmean,rr,rrij,rrjk,rrki,rijk,r0,r6,r8 185 real(dp) :: sfact6,sfact8,sfact9,sum_dlri,sum_dlrj,sum_dlc6ri,sum_dlc6rj 186 real(dp) :: sum_d2lri,sum_d2lrj,sum_d2lrirj,sum_d2lc6ri,sum_d2lc6rj,sum_d2lc6rirj 187 real(dp) :: temp,temp2 188 real(dp) :: ucvol,vdw_s6,vdw_s8,vdw_sr6,vdw_a1,vdw_a2,vdw_q 189 character(len=500) :: msg 190 type(atomdata_t) :: atom1,atom2 191 192 !arrays 193 194 ! Covalence radius of the different species for CN (coordination number) 195 real(dp),parameter:: rcov(vdw_nspecies)=& 196 & (/0.80628308, 1.15903197, 3.02356173, 2.36845659, 1.94011865, & 197 & 1.88972601, 1.78894056, 1.58736983, 1.61256616, 1.68815527, & 198 & 3.52748848, 3.14954334, 2.84718717, 2.62041997, 2.77159820, & 199 & 2.57002732, 2.49443835, 2.41884923, 4.43455700, 3.88023730, & 200 & 3.35111422, 3.07395437, 3.04875805, 2.77159820, 2.69600923, & 201 & 2.62041997, 2.51963467, 2.49443835, 2.54483100, 2.74640188, & 202 & 2.82199085, 2.74640188, 2.89757982, 2.77159820, 2.87238349, & 203 & 2.94797246, 4.76210950, 4.20778980, 3.70386304, 3.50229216, & 204 & 3.32591790, 3.12434702, 2.89757982, 2.84718717, 2.84718717, & 205 & 2.72120556, 2.89757982, 3.09915070, 3.22513231, 3.17473967, & 206 & 3.17473967, 3.09915070, 3.32591790, 3.30072128, 5.26603625, & 207 & 4.43455700, 4.08180818, 3.70386304, 3.98102289, 3.95582657, & 208 & 3.93062995, 3.90543362, 3.80464833, 3.82984466, 3.80464833, & 209 & 3.77945201, 3.75425569, 3.75425569, 3.72905937, 3.85504098, & 210 & 3.67866672, 3.45189952, 3.30072128, 3.09915070, 2.97316878, & 211 & 2.92277614, 2.79679452, 2.82199085, 2.84718717, 3.32591790, & 212 & 3.27552496, 3.27552496, 3.42670319, 3.30072128, 3.47709584, & 213 & 3.57788113, 5.06446567, 4.56053862, 4.20778980, 3.98102289, & 214 & 3.82984466, 3.85504098, 3.88023730, 3.90543362 /) 215 216 ! q = arrays of vdw_species elements containing the link between C6ij and C8ij: 217 ! C8ij = 3sqrt(qi)sqrt(qj)C6ij 218 real(dp),parameter :: vdw_q_dftd3(vdw_nspecies)= & 219 & (/2.00734898, 1.56637132, 5.01986934, 3.85379032, 3.64446594, & 220 & 3.10492822, 2.71175247, 2.59361680, 2.38825250, 2.21522516, & 221 & 6.58585536, 5.46295967, 5.65216669, 4.88284902, 4.29727576, & 222 & 4.04108902, 3.72932356, 3.44677275, 7.97762753, 7.07623947, & 223 & 6.60844053, 6.28791364, 6.07728703, 5.54643096, 5.80491167, & 224 & 5.58415602, 5.41374528, 5.28497229, 5.22592821, 5.09817141, & 225 & 6.12149689, 5.54083734, 5.06696878, 4.87005108, 4.59089647, & 226 & 4.31176304, 9.55461698, 8.67396077, 7.97210197, 7.43439917, & 227 & 6.58711862, 6.19536215, 6.01517290, 5.81623410, 5.65710424, & 228 & 5.52640661, 5.44263305, 5.58285373, 7.02081898, 6.46815523, & 229 & 5.98089120, 5.81686657, 5.53321815, 5.25477007, 11.02204549, & 230 & 10.15679528, 9.35167836, 9.06926079, 8.97241155, 8.90092807, & 231 & 8.85984840, 8.81736827, 8.79317710, 7.89969626, 8.80588454, & 232 & 8.42439218, 8.54289262, 8.47583370, 8.45090888, 8.47339339, & 233 & 7.83525634, 8.20702843, 7.70559063, 7.32755997, 7.03887381, & 234 & 6.68978720, 6.05450052, 5.88752022, 5.70661499, 5.78450695, & 235 & 7.79780729, 7.26443867, 6.78151984, 6.67883169, 6.39024318, & 236 & 6.09527958, 11.79156076, 11.10997644, 9.51377795, 8.67197068, & 237 & 8.77140725, 8.65402716, 8.53923501, 8.85024712 /) 238 239 integer :: is(3), nshell_3bt(3) 240 integer,allocatable :: ivdw(:) 241 integer :: jmin(3), jmax(3), js(3) 242 integer,parameter :: voigt1(6)=(/1,2,3,2,1,1/),voigt2(6)=(/1,2,3,3,3,2/) 243 real(dp),allocatable :: cn(:),cfgrad_no_c(:,:,:,:) 244 real(dp),allocatable :: dcn(:,:,:),dcn_cart(:,:,:),dc6ri(:,:),dc6rj(:,:),dc9ijri(:,:),dc9ijrj(:,:) 245 real(dp),allocatable:: d2cn(:,:,:,:,:,:) 246 real(dp),allocatable :: d2c6ri(:,:),d2c6rj(:,:),d2c6rirj(:,:) 247 real(dp),allocatable:: elt_cn(:,:,:),e3bt_ij(:,:),e3bt_jk(:,:),e3bt_ki(:,:),e_no_c(:,:) 248 real(dp),allocatable:: e_alpha1(:),e_alpha2(:),e_alpha3(:),e_alpha4(:) 249 real(dp) :: fred(3),fredij(3),fredjk(3),fredki(3) 250 real(dp),allocatable:: fe_no_c(:,:,:),cfdcn(:,:,:,:),fdcn(:,:,:,:),fgrad_no_c(:,:,:,:),fred_vdw_3bt(:,:) 251 real(dp) :: gmet(3,3),gprimd(3,3) 252 real(dp),allocatable:: grad_no_cij(:,:,:) 253 real(dp) :: mcart(3,3) 254 real(dp) :: r(3),rcart(3),rcart2(3,3),rcartij(3),rcartjk(3),rcartki(3) 255 real(dp):: rij(3), rjk(3), rki(3),rmet(3,3),rred(3) 256 real(dp),allocatable:: r0ijk(:,:,:) 257 real(dp),allocatable :: str_alpha1(:,:),str_alpha2(:,:),str_dcn(:,:),str_no_c(:,:,:) 258 real(dp) :: str_3bt(6) 259 real(dp) :: temp_comp(2),temp_comp2(2) 260 real(dp),allocatable:: temp_prod(:,:) 261 real(dp) :: vec(6),vecij(6), vecjk(6),vecki(6) 262 real(dp),allocatable :: vdw_cnrefi(:,:,:,:),vdw_cnrefj(:,:,:,:) 263 real(dp),allocatable :: vdw_c6(:,:),vdw_c6ref(:,:,:,:),vdw_c8(:,:),vdw_c9(:,:,:),vdw_r0(:,:) 264 real(dp):: vdw_dftd3_r0(4465) 265 real(dp):: vdw_dftd3_c6(32385) 266 integer:: index_c6(254) 267 real(dp):: vdw_dftd3_cni(27884) 268 integer:: index_cni(27884) 269 real(dp):: vdw_dftd3_cnj(13171) 270 integer:: index_cnj(13171) 271 real(dp),allocatable :: xred01(:,:) 272 273 ! ************************************************************************* 274 275 DBG_ENTER("COLL") 276 277 write(msg,'(1a)')& 278 & '====> STARTING DFT-D3 computation' 279 call wrtout(std_out,msg,'COLL') 280 281 call vdw_dftd3_data(vdw_dftd3_r0,vdw_dftd3_c6,index_c6,vdw_dftd3_cni,index_cni,& 282 & vdw_dftd3_cnj,index_cnj) 283 ! Determine the properties which have to be studied 284 bol_3bt = (vdw_tol_3bt>0) 285 need_forces = present(fred_vdw_dftd3) 286 need_stress= present(str_vdw_dftd3) 287 need_dynmat= present(dyn_vdw_dftd3) 288 need_elast= present(elt_vdw_dftd3) 289 need_grad=(need_forces.or.need_stress.or.need_dynmat.or.need_elast) 290 need_hess=(need_dynmat.or.need_elast) 291 if (need_dynmat) then 292 if (.not.present(qphon)) then 293 msg='Dynamical matrix required without a q-vector' 294 MSG_BUG(msg) 295 end if 296 dyn_vdw_dftd3=zero 297 end if 298 e_vdw_dftd3 = zero 299 if (need_forces) fred_vdw_dftd3=zero 300 if (need_stress) str_vdw_dftd3=zero 301 if (need_elast) elt_vdw_dftd3=zero 302 303 !Identify type(s) of atoms 304 ABI_ALLOCATE(ivdw,(ntypat)) 305 do itypat=1,ntypat 306 call atomdata_from_znucl(atom1,znucl(itypat)) 307 if (znucl(itypat).gt.94.0_dp) then 308 write(msg,'(3a,es14.2)') & 309 & 'Van der Waals DFT-D3 correction not available for atom type: ',znucl(itypat),' !' 310 MSG_ERROR(msg) 311 else 312 ivdw(itypat) = znucl(itypat) 313 end if 314 end do 315 316 ! Determination of coefficients that depend of the 317 ! exchange-correlation functional 318 vdw_s6=one 319 vdw_s8=one 320 ! Case one : DFT-D3 321 if (vdw_xc == 6) then 322 select case (ixc) 323 case (11,-101130,-130101) 324 vdw_sr6 = vdw_sr6_pbe ; vdw_s8 = vdw_s8*vdw_s8_pbe 325 case (18,-106131,-131106) 326 vdw_sr6=vdw_sr6*vdw_sr6_blyp ; vdw_s8 =vdw_s8*vdw_s8_blyp 327 case (19,-106132,-132106) 328 vdw_sr6=vdw_sr6*vdw_sr6_bp86 ; vdw_s8=vdw_s8*vdw_s8_bp86 329 case (-202231,-231202) 330 vdw_sr6=vdw_sr6*vdw_sr6_tpss ; vdw_s8=vdw_s8*vdw_s8_tpss 331 case (14,-102130,-130102) 332 vdw_sr6=vdw_sr6*vdw_sr6_revpbe ; vdw_s8=vdw_s8*vdw_s8_revpbe 333 case (-170) 334 vdw_sr6=vdw_sr6*vdw_sr6_b97d ; vdw_s8=vdw_s8*vdw_s8_b97d 335 case (41,-406) 336 vdw_sr6=vdw_sr6*vdw_sr6_pbe0 ; vdw_s8=vdw_s8*vdw_s8_pbe0 337 case default 338 write(msg,'(a,i8,a)')' Van der Waals DFT-D3 correction not compatible with ixc=',ixc,' !' 339 MSG_ERROR(msg) 340 end select 341 ! Case DFT-D3(BJ) 342 elseif (vdw_xc == 7) then 343 select case (ixc) 344 case (11,-101130,-130101) 345 vdw_a1 = vdw_a1_pbe ; vdw_a2 = vdw_a2_pbe ; vdw_s8= vdwbj_s8_pbe 346 case (18,-106131,-131106) 347 vdw_a1=vdw_a1_blyp ; vdw_a2=vdw_a2_blyp ; vdw_s8=vdwbj_s8_blyp 348 case (19,-106132,-132106) 349 vdw_a1=vdw_a1_bp86 ; vdw_a2=vdw_a2_bp86 ; vdw_s8=vdwbj_s8_bp86 350 case (-202231,-231202) 351 vdw_a1=vdw_a1_tpss ; vdw_a2=vdw_a2_tpss ; vdw_s8=vdwbj_s8_tpss 352 case (14,-102130,-130102) 353 vdw_a1=vdw_a1_revpbe ; vdw_a2=vdw_a2_revpbe ; vdw_s8=vdwbj_s8_revpbe 354 case (-170) 355 vdw_a1=vdw_a1_b97d ; vdw_a2=vdw_a2_b97d ; vdw_s8=vdwbj_s8_b97d 356 case (41,-406) 357 vdw_a1=vdw_a1_pbe0 ; vdw_a2=vdw_a2_pbe0 ; vdw_s8=vdwbj_s8_pbe0 358 end select 359 end if 360 ! -------------------------------------------------------------- 361 ! Retrieve the data for the referenced c6, cn and r0 coefficients 362 !--------------------------------------------------------------- 363 refmax = 5 364 365 ABI_ALLOCATE(vdw_c6ref,(ntypat,ntypat,refmax,refmax)) 366 ABI_ALLOCATE(vdw_cnrefi,(ntypat,ntypat,refmax,refmax)) 367 ABI_ALLOCATE(vdw_cnrefj,(ntypat,ntypat,refmax,refmax)) 368 ABI_ALLOCATE(vdw_r0,(ntypat,ntypat)) 369 if (bol_3bt) then 370 ABI_ALLOCATE(r0ijk,(ntypat,ntypat,ntypat)) 371 end if 372 373 vdw_c6ref = zero 374 vdw_cnrefi = 100 ; vdw_cnrefj = 100 ; 375 376 do refi=1,refmax 377 do refj=1,refi 378 do itypat=1,ntypat 379 do jtypat=1,ntypat 380 indi = ivdw(itypat)+100*(refi-1) 381 indj = ivdw(jtypat)+100*(refj-1) 382 found = .false. 383 do ia=1,size(index_c6) 384 do ja=1,size(index_c6) 385 if (index_c6(ia)==indi.and.index_c6(ja)==indj) then 386 nline = ia*(ia-1)/2 + ja 387 vdw_c6ref(itypat,jtypat,refi,refj) = vdw_dftd3_c6(nline) 388 vdw_c6ref(jtypat,itypat,refj,refi) = vdw_dftd3_c6(nline) 389 found = .false. 390 do la=1,size(index_cni) 391 if (index_cni(la)==nline) then 392 found=.true. 393 vdw_cnrefi(itypat,jtypat,refi,refj)= vdw_dftd3_cni(la) 394 vdw_cnrefj(jtypat,itypat,refj,refi)= vdw_dftd3_cni(la) 395 else 396 vdw_cnrefi(itypat,jtypat,refi,refj) = zero 397 vdw_cnrefj(jtypat,itypat,refj,refi) = zero 398 end if 399 if (found) exit 400 end do 401 found = .false. 402 do la=1,size(index_cnj) 403 if (index_cnj(la)==nline) then 404 found=.true. 405 vdw_cnrefj(itypat,jtypat,refi,refj)= vdw_dftd3_cnj(la) 406 vdw_cnrefi(jtypat,itypat,refj,refi)= vdw_dftd3_cnj(la) 407 else 408 vdw_cnrefj(itypat,jtypat,refi,refj) = zero 409 vdw_cnrefi(jtypat,itypat,refj,refi) = zero 410 end if 411 if (found) exit 412 end do 413 found = .true. 414 end if 415 if (found) exit 416 end do 417 if (found) exit 418 end do 419 if (refi.eq.1.and.refj.eq.1) then 420 nline = ia*(ia-1)/2 + ja 421 vdw_r0(itypat,jtypat)=vdw_dftd3_r0(nline)/Bohr_Ang 422 if (bol_3bt) then 423 do ktypat=1,ntypat 424 r0ijk(itypat,jtypat,ktypat)=one/(vdw_r0(itypat,jtypat)*vdw_r0(jtypat,ktypat)*vdw_r0(ktypat,itypat))**third 425 end do ! ka atom 426 end if ! Only if 3bt required 427 end if ! Only for the first set of references 428 end do ! Loop on references j 429 end do ! Loop on references i 430 end do ! Loop on atom j 431 end do ! Loop on atom i 432 433 !if (vdw_d3_cov==1) then 434 ! vdw_cnrefi(:,:,:,refmax) =vdw_cnrefi(:,:,:,refmax-1) 435 ! vdw_cnrefi(:,:,refmax,:) = 14.0_dp 436 ! vdw_cnrefj(:,:,refmax,:) =vdw_cnrefj(:,:,refmax-1,:) 437 ! vdw_cnrefj(:,:,:,refmax) = 14.0_dp 438 !end if 439 !Retrieve cell geometry data 440 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 441 442 !Map reduced coordinates into [0,1[ 443 ABI_ALLOCATE(xred01,(3,natom)) 444 do ia=1,natom 445 xred01(:,ia)=xred(:,ia)-aint(xred(:,ia)) ! Map into ]-1,1[ 446 do alpha=1,3 447 if (abs(xred01(alpha,ia)).ge.tol8) xred01(alpha,ia) = xred01(alpha,ia)+half-sign(half,xred(alpha,ia)) 448 end do 449 end do 450 451 ! ------------------------------------------------------------------- 452 ! Computation of the coordination number (CN) for the different atoms 453 ! ------------------------------------------------------------------- 454 455 write(msg,'(3a)')& 456 & ' Begin the computation of the Coordination Numbers (CN)',ch10,& 457 & ' required for DFT-D3 energy corrections...' 458 call wrtout(std_out,msg,'COLL') 459 460 ! Allocation of the CN coefficients and derivatives 461 ABI_ALLOCATE(cn,(natom)) 462 ABI_ALLOCATE(dcn,(3,natom,natom)) 463 ABI_ALLOCATE(dcn_cart,(3,natom,natom)) 464 ABI_ALLOCATE(str_dcn,(6,natom)) 465 ABI_ALLOCATE(d2cn,(2,3,natom,3,natom,natom)) 466 ABI_ALLOCATE(fdcn,(2,3,natom,natom)) 467 ABI_ALLOCATE(cfdcn,(2,3,natom,natom)) 468 ABI_ALLOCATE(elt_cn,(6+3*natom,6,natom)) 469 470 ! Initializing the computed quantities to zero 471 nshell = 0 472 cn = zero 473 ! Initializing the derivative of the computed quantities to zero (if required) 474 dcn = zero ; str_dcn = zero 475 d2cn = zero ;fdcn = zero; cfdcn = zero 476 elt_cn = zero ; dcn_cart = zero 477 478 re_arg = zero ; im_arg = zero 479 if (need_hess) then 480 mcart = zero 481 do alpha=1,3 482 mcart(alpha,alpha) = one 483 end do 484 end if 485 rcutcn = 200**2 ! Bohr 486 !Loop over shells of cell replicas 487 do 488 newshell=.false.;nshell=nshell+1 489 ! Loop over cell replicas in the shell 490 do is3=-nshell,nshell 491 do is2=-nshell,nshell 492 do is1=-nshell,nshell 493 if (nshell==1.or. & 494 & abs(is3)==nshell.or.abs(is2)==nshell.or.abs(is1)==nshell) then 495 is(3) = is3 ; is(2) = is2 ; is(1) = is1 496 ! Computation of phase factor for discrete Fourier transform 497 ! if phonon at qphon is required 498 if (need_dynmat) then 499 arg=two_pi*dot_product(qphon,is) 500 re_arg=cos(arg) ; im_arg=sin(arg) 501 end if 502 ! Loop over atoms ia and ja 503 do ia=1,natom 504 itypat=typat(ia) 505 do ja=1,natom 506 jtypat=typat(ja) 507 r(:)=xred01(:,ia)-xred01(:,ja)-dble(is(:)) 508 rsq=dot_product(r,matmul(rmet,r)) 509 ! atom i =/= j 510 if (rsq.ge.tol16.and.rsq<rcutcn) then 511 newshell=.true. 512 rr = sqrt(rsq) 513 rcovij = rcov(ivdw(itypat))+rcov(ivdw(jtypat)) 514 515 ! Computation of partial contribution to cn coefficients 516 exp_cn = exp(-k1*(rcovij/rr-one)) 517 frac_cn= one/(one+exp_cn) 518 ! Introduction of a damping function for the coordination 519 ! number because of the divergence with increasing 520 ! number of cells of this quantity in periodic systems 521 ! See Reckien et al., J. Chem. Phys. 132, 154104 (2010) 522 dr = rr-k2*rcovij 523 cn_dmp = half*abi_derfc(dr) 524 cn(ia) = cn(ia)+frac_cn*cn_dmp 525 526 ! If force, stress, IFC or Elastic constants are required, 527 ! computation of the first derivative of CN 528 if (need_grad) then 529 rcart=matmul(rprimd,r) 530 dexp_cn= k1*rcovij*exp_cn/rsq 531 dcn_dmp = -one/sqrt(pi)*exp(-dr*dr) 532 grad=(-frac_cn*frac_cn*cn_dmp*dexp_cn+dcn_dmp*frac_cn)/rr 533 if (need_forces.and.ia/=ja) then 534 ! Variation of CN(ia) w.r. to displacement of atom 535 ! ja. If ia==ka then all the other atoms contribute 536 ! to the derivative. Required for the computation of 537 ! the forces applied on atom k 538 rred = matmul(transpose(rprimd),rcart) 539 dcn(:,ia,ia) = dcn(:,ia,ia)+grad*rred(:) 540 dcn(:,ia,ja) = dcn(:,ia,ja)-grad*rred(:) 541 elseif (need_stress.or.need_elast) then 542 ! The following quantity (str_dcn) is used for the computation 543 ! of the DFT-D3 contribution to stress and elastic constants 544 vec(1:3)=rcart(1:3)*rcart(1:3); vec(4)=rcart(2)*rcart(3) 545 vec(5)=rcart(1)*rcart(3); vec(6)=rcart(1)*rcart(2) 546 str_dcn(:,ia)=str_dcn(:,ia)+grad*vec(:) 547 end if 548 ! If dynamical matrix or elastic constants are required, compute 549 ! the second derivative 550 if (need_hess) then 551 d2cn_dmp = two*(rr-k2*rcovij)/sqrt(pi)*exp(-(rr-k2*rcovij)**two) 552 d2cn_exp = dexp_cn*(k1*rcovij/rsq-two/rr) 553 d2frac_cn =frac_cn**two*(two*frac_cn*dexp_cn**two-d2cn_exp) 554 hess = (d2frac_cn*cn_dmp+d2cn_dmp*frac_cn-& 555 & two*dcn_dmp*frac_cn**two*dexp_cn-grad)/rsq 556 if (need_dynmat) then 557 ! Discrete Fourier Transform of dCN/drk in cartesian 558 ! coordinates. Note that the phase factor is different 559 ! if (ka=ia or ka/=ia). 560 ! This Fourier transform is summed over cell replica 561 ! See NOTE: add reference for more informations 562 fdcn(1,:,ia,ia) = fdcn(1,:,ia,ia)+grad*rcart(:) 563 fdcn(1,:,ia,ja) = fdcn(1,:,ia,ja)-grad*rcart(:)*re_arg 564 fdcn(2,:,ia,ja) = fdcn(2,:,ia,ja)-grad*rcart(:)*im_arg 565 ! Conjugate of fdcn 566 cfdcn(1,:,ia,ia) = cfdcn(1,:,ia,ia)+grad*rcart(:) 567 cfdcn(1,:,ia,ja) = cfdcn(1,:,ia,ja)-grad*rcart(:)*re_arg 568 cfdcn(2,:,ia,ja) = cfdcn(2,:,ia,ja)+grad*rcart(:)*im_arg 569 ! Computation of second derivative of CN required for the 570 ! interatomic force constants in reciprocal space 571 do alpha=1,3 572 rcart2(alpha,:) = rcart(alpha)*rcart(:) 573 end do 574 ! Computation of second derivative of CN required for the 575 ! interatomic force constants in reciprocal space 576 ! This Fourier transform is summed over cell replica 577 ! as it appears in the theory 578 do alpha=1,3 579 if (ia/=ja) then 580 ! ka = ia ; la = ia 581 d2cn(1,alpha,ia,:,ia,ia) = d2cn(1,alpha,ia,:,ia,ia)+& 582 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:)) 583 ! ka = ja ; la = ja 584 d2cn(1,alpha,ja,:,ja,ia) = d2cn(1,alpha,ja,:,ja,ia)+& 585 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:)) 586 if (abs(re_arg)>tol12) then 587 ! ka = ia ; la = ja 588 d2cn(1,alpha,ia,:,ja,ia) = d2cn(1,alpha,ia,:,ja,ia)-& 589 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:))*re_arg 590 ! ka = ja ; la = ia 591 d2cn(1,alpha,ja,:,ia,ia) = d2cn(1,alpha,ja,:,ia,ia)-& 592 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:))*re_arg 593 end if 594 if (abs(im_arg)>tol12) then 595 ! ka = ia ; la = ja 596 d2cn(2,alpha,ia,:,ja,ia) = d2cn(2,alpha,ia,:,ja,ia)-& 597 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:))*im_arg 598 ! ka = ja ; la = ia 599 d2cn(2,alpha,ja,:,ia,ia) = d2cn(2,alpha,ja,:,ia,ia)+& 600 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:))*im_arg 601 end if 602 else 603 if (abs(re_arg-one)>tol12) then 604 d2cn(1,alpha,ia,:,ja,ia) = d2cn(1,alpha,ia,:,ja,ia)+& 605 & two*(hess*rcart2(alpha,:)+grad*mcart(alpha,:))*(one-re_arg) 606 end if 607 end if 608 end do 609 end if ! Boolean Need_dynmat 610 if (need_elast) then 611 ! Derivative of str_dcn w.r. to strain for elastic tensor 612 vec(1:3)=rcart(1:3)*rcart(1:3); vec(4)=rcart(2)*rcart(3) 613 vec(5)=rcart(1)*rcart(3); vec(6)=rcart(1)*rcart(2) 614 do alpha=1,6 615 ii = voigt1(alpha) ; jj=voigt2(alpha) 616 do beta=1,6 617 kk = voigt1(beta) ; ll=voigt2(beta) 618 elt_cn(alpha,beta,ia) = elt_cn(alpha,beta,ia)+hess*rcart(ii)*rcart(jj)*rcart(kk)*rcart(ll) 619 if (ii==kk) elt_cn(alpha,beta,ia) = elt_cn(alpha,beta,ia)+half*grad*rcart(jj)*rcart(ll) 620 if (jj==kk) elt_cn(alpha,beta,ia) = elt_cn(alpha,beta,ia)+half*grad*rcart(ii)*rcart(ll) 621 if (ii==ll) elt_cn(alpha,beta,ia) = elt_cn(alpha,beta,ia)+half*grad*rcart(jj)*rcart(kk) 622 if (jj==ll) elt_cn(alpha,beta,ia) = elt_cn(alpha,beta,ia)+half*grad*rcart(ii)*rcart(kk) 623 end do 624 end do 625 ! Derivative of str_dcn w.r. to atomic displacement 626 ! for internal strains 627 dcn_cart(:,ia,ia) = dcn_cart(:,ia,ia)+grad*rcart(:) 628 dcn_cart(:,ia,ja) = dcn_cart(:,ia,ja)-grad*rcart(:) 629 if (ia/=ja) then 630 index_ia = 6+3*(ia-1) 631 index_ja = 6+3*(ja-1) 632 do alpha=1,6 633 do beta=1,3 634 ii = voigt1(alpha) ; jj=voigt2(alpha) 635 elt_cn(index_ia+beta,alpha,ia)=elt_cn(index_ia+beta,alpha,ia)& 636 & +hess*vec(alpha)*rcart(beta) 637 elt_cn(index_ja+beta,alpha,ia)=elt_cn(index_ja+beta,alpha,ia)& 638 & -hess*vec(alpha)*rcart(beta) 639 if (ii==beta) then 640 elt_cn(index_ia+beta,alpha,ia)=elt_cn(index_ia+beta,alpha,ia)+grad*rcart(jj) 641 elt_cn(index_ja+beta,alpha,ia)=elt_cn(index_ja+beta,alpha,ia)-grad*rcart(jj) 642 end if 643 if (jj==beta) then 644 elt_cn(index_ia+beta,alpha,ia)=elt_cn(index_ia+beta,alpha,ia)+grad*rcart(ii) 645 elt_cn(index_ja+beta,alpha,ia)=elt_cn(index_ja+beta,alpha,ia)-grad*rcart(ii) 646 end if 647 end do 648 end do 649 end if ! ia/=ja 650 end if ! Need strain derivative 651 end if ! Boolean second derivative 652 end if ! Boolean first derivative 653 end if ! Tolerence 654 end do ! Loop over ia atom 655 end do ! Loop over ja atom 656 end if ! Bondary Condition 657 end do ! Loop over is1 658 end do ! Loop over is2 659 end do ! Loop over is3 660 if(.not.newshell) exit ! Check if a new shell must be considered 661 end do ! Loop over shell 662 write(msg,'(3a,f8.5,1a,i3,1a,f8.5,1a,i3,1a)')& 663 & ' ... Done.',ch10,& 664 & ' max(CN) =', maxval(cn), ' (atom ',maxloc(cn),') ; min(CN) =', minval(cn), ' (atom ', minloc(cn),')' 665 call wrtout(std_out,msg,'COLL') 666 667 !---------------------------------------------------------------- 668 ! Computation of the C6 coefficient 669 ! --------------------------------------------------------------- 670 671 write(msg,'(3a)')& 672 & ' Begin the computation of the C6(CN)',ch10,& 673 & ' required for DFT-D3 energy corrections...' 674 call wrtout(std_out,msg,'COLL') 675 ! Allocation 676 ABI_ALLOCATE(vdw_c6,(natom,natom)) 677 ABI_ALLOCATE(vdw_c8,(natom,natom)) 678 ABI_ALLOCATE(dc6ri,(natom,natom)) 679 ABI_ALLOCATE(dc6rj,(natom,natom)) 680 ABI_ALLOCATE(d2c6ri,(natom,natom)) 681 ABI_ALLOCATE(d2c6rj,(natom,natom)) 682 ABI_ALLOCATE(d2c6rirj,(natom,natom)) 683 if (bol_3bt) then 684 ABI_ALLOCATE(vdw_c9,(natom,natom,natom)) 685 ABI_ALLOCATE(dc9ijri,(natom,natom)) 686 ABI_ALLOCATE(dc9ijrj,(natom,natom)) 687 end if 688 ! Set accumulating quantities to zero 689 vdw_c6 = zero ; vdw_c8 = zero 690 dc6ri = zero ; dc6rj = zero 691 d2c6ri = zero ; d2c6rj = zero 692 d2c6rirj = zero 693 if (bol_3bt) then 694 dc9ijri = zero; dc9ijrj = zero 695 end if 696 ! C6 coefficients are interpolated from tabulated 697 ! ab initio C6 values (following loop). 698 ! C8 coefficients are obtained by: 699 ! C8 = vdw_dftd3_q(itypat)*vdw_dftd3_q(jtypat)*C6 700 701 do ia=1,natom 702 itypat=typat(ia) 703 do ja=1,natom 704 jtypat=typat(ja) 705 ! Set accumulating quantities to zero 706 ltot=zero 707 sum_dlri = zero ; sum_dlc6ri= zero 708 sum_dlrj = zero ; sum_dlc6rj= zero 709 sum_d2lri = zero ; sum_d2lc6ri = zero 710 sum_d2lrj = zero ; sum_d2lc6rj = zero 711 sum_d2lrirj = zero ; sum_d2lc6rirj = zero 712 min_dsys = 10000 713 max_vdw_c6 = zero 714 ! Loop over references 715 do refi=1,refmax 716 do refj=1,refmax 717 dsysref_a = cn(ia)-vdw_cnrefi(itypat,jtypat,refi,refj) 718 dsysref_b = cn(ja)-vdw_cnrefj(itypat,jtypat,refi,refj) 719 dsysref=(dsysref_a)**two+(dsysref_b)**two 720 if (dsysref<min_dsys) then 721 ! Keep in memory the smallest value of dsysref 722 ! And the associated tabulated C6 value 723 min_dsys = dsysref 724 max_vdw_c6 = vdw_c6ref(itypat,jtypat,refi,refj) 725 end if 726 l = dexp(-k3*dsysref) 727 ltot = ltot+l 728 vdw_c6(ia,ja)=vdw_c6(ia,ja)+vdw_c6ref(itypat,jtypat,refi,refj)*l 729 730 if (need_grad) then 731 ! Derivative of l(ia,ja) with respect to the displacement 732 ! of atom ka in reduced coordinates. 733 ! This factor is identical in the case of stress. 734 ! In purpose of speed up this routine, the prefactor of 735 ! dCNi/drk and dCNj/drk are separated 736 ! See NOTE: article to be added 737 dlri=-k3*l*two*dsysref_a ;dlrj=-k3*l*two*dsysref_b 738 sum_dlri=sum_dlri+dlri ; sum_dlrj=sum_dlrj+dlrj 739 sum_dlc6ri=sum_dlc6ri+dlri*vdw_c6ref(itypat,jtypat,refi,refj) 740 sum_dlc6rj=sum_dlc6rj+dlrj*vdw_c6ref(itypat,jtypat,refi,refj) 741 if (need_hess) then 742 ! Second derivative of l(ia,ja). Once again, it is separately in 743 ! different contributions: 744 ! d2lri: prefactor of dCNi/drk*dCNi/drl 745 ! d2lrj: prefactor of dCNj/drk*dCNj/drl 746 ! d2lrirj: prefacto of dCNi/drk*dCNj/drl 747 ! The prefactor for d2CNi/drkdrl is dlri; for d2CNj/drkdrl is dlrj 748 d2lri = -two*k3*l*(one-two*k3*dsysref_a**two) 749 d2lrj = -two*k3*l*(one-two*k3*dsysref_b**two) 750 d2lrirj = four*k3*k3*l*(dsysref_a*dsysref_b) 751 sum_d2lri=sum_d2lri+d2lri ; sum_d2lrj=sum_d2lrj+d2lrj 752 sum_d2lrirj = sum_d2lrirj+d2lrirj 753 sum_d2lc6ri=sum_d2lc6ri+d2lri*vdw_c6ref(itypat,jtypat,refi,refj) 754 sum_d2lc6rj=sum_d2lc6rj+d2lrj*vdw_c6ref(itypat,jtypat,refi,refj) 755 sum_d2lc6rirj = sum_d2lc6rirj + d2lrirj*vdw_c6ref(itypat,jtypat,refi,refj) 756 end if ! Boolean second derivative 757 end if ! Boolean gradient 758 end do ! Loop over references 759 end do ! Loop over references 760 ! In some specific case (really covalently bound compounds) ltot -> 0 761 ! which may cause numerical problems for all quantities related to dispersion coefficient. 762 ! To be consistent with VASP implementation, the c6 value is taken as the last 763 ! referenced value of the dispersion coefficient. 764 if (ltot>tol12) then 765 vdw_c6(ia,ja)=vdw_c6(ia,ja)/ltot 766 vdw_c8(ia,ja)=three*vdw_q_dftd3(ivdw(itypat))*vdw_q_dftd3(ivdw(jtypat))*vdw_c6(ia,ja) 767 ! If Force of Stress is required 768 if (need_grad) then 769 ! Computation of the derivative of C6 w.r.to the displacement 770 ! of atom ka, in reduced coordinates (separated for dCNi/drk and dCNj/drk) 771 ! This is the crucial step to reduce the scaling from O(N^3) to O(N^2) for 772 ! the gradients 773 dc6ri(ia,ja)=(sum_dlc6ri-vdw_c6(ia,ja)*sum_dlri)/ltot 774 dc6rj(ia,ja)=(sum_dlc6rj-vdw_c6(ia,ja)*sum_dlrj)/ltot 775 if (need_hess) then 776 ! Computation of the second derivative of C6 w.r.to the displacement of atom ka 777 ! and atom la 778 d2c6ri(ia,ja)=(sum_d2lc6ri-vdw_c6(ia,ja)*sum_d2lri-two*dc6ri(ia,ja)*sum_dlri)/ltot 779 d2c6rj(ia,ja)=(sum_d2lc6rj-vdw_c6(ia,ja)*sum_d2lrj-two*dc6rj(ia,ja)*sum_dlrj)/ltot 780 d2c6rirj(ia,ja) = (sum_d2lc6rirj-vdw_c6(ia,ja)*sum_d2lrirj-dc6ri(ia,ja)*& 781 & sum_dlrj-dc6rj(ia,ja)*sum_dlri)/ltot 782 end if ! Boolean second derivative 783 end if ! Boolean gradient 784 else 785 vdw_c6(ia,ja)= max_vdw_c6 786 vdw_c8(ia,ja)=three*vdw_q_dftd3(ivdw(itypat))*vdw_q_dftd3(ivdw(jtypat))*vdw_c6(ia,ja) 787 end if 788 end do 789 end do 790 ! Computation of the three-body term dispersion coefficient 791 if (bol_3bt) then 792 do ia=1,natom 793 do ja=1,natom 794 do ka=1,natom 795 vdw_c9(ia,ja,ka) =-sqrt(vdw_c6(ia,ja)*vdw_c6(ja,ka)*vdw_c6(ka,ia)) 796 end do 797 if (need_grad) then 798 dc9ijri(ia,ja) = half/vdw_c6(ia,ja)*dc6ri(ia,ja) 799 dc9ijrj(ia,ja) = half/vdw_c6(ia,ja)*dc6rj(ia,ja) 800 end if 801 end do 802 end do 803 end if 804 805 write(msg,'(3a,f8.5,1a,f8.5)')& 806 & ' ... Done.',ch10,& 807 & ' max(C6) =', maxval(vdw_c6),' ; min(C6) =', minval(vdw_c6) 808 call wrtout(std_out,msg,'COLL') 809 810 ! Deallocation of used variables not needed anymore 811 ABI_DEALLOCATE(vdw_c6ref) 812 ABI_DEALLOCATE(vdw_cnrefi) 813 ABI_DEALLOCATE(vdw_cnrefj) 814 815 !---------------------------------------------------- 816 ! Computation of cut-off radii according to tolerance 817 !---------------------------------------------------- 818 819 ! Cut-off radius for pair-wise term 820 if (vdw_tol<zero) then 821 rcut=max((vdw_s6/vdw_tol_default*maxval(vdw_c6))**sixth, & 822 & (vdw_s8/vdw_tol_default*maxval(vdw_c8))**(one/eight)) 823 else 824 rcut=max((vdw_s6/vdw_tol*maxval(vdw_c6))**sixth,& 825 & (vdw_s8/vdw_tol*maxval(vdw_c8))**(one/eight)) 826 end if 827 ! Cut-off radius for three-body term 828 rcut9 = zero 829 if (bol_3bt) then 830 rcut9=(128.0_dp*vdw_s6/(vdw_tol_3bt)*maxval(vdw_c6)**(3.0/2.0))**(1.0/9.0) 831 end if 832 rcut2=rcut*rcut 833 834 !-------------------------------------------------------------------- 835 ! Computation of the two bodies contribution to the dispersion energy 836 !-------------------------------------------------------------------- 837 838 write(msg,'(3a)')& 839 & ' Begin the computation of pair-wise term',ch10,& 840 & ' of DFT-D3 energy contribution...' 841 call wrtout(std_out,msg,'COLL') 842 nshell=0 843 npairs=0 844 ABI_ALLOCATE(e_alpha1,(natom)) 845 ABI_ALLOCATE(e_alpha2,(natom)) 846 ABI_ALLOCATE(e_alpha3,(natom)) 847 ABI_ALLOCATE(e_alpha4,(natom)) 848 ABI_ALLOCATE(e_no_c,(natom,natom)) 849 ABI_ALLOCATE(fe_no_c,(2,natom,natom)) 850 e_alpha1 =zero ; e_alpha2 = zero 851 e_alpha3 =zero ; e_alpha4 = zero 852 e_no_c=zero ; fe_no_c = zero 853 ABI_ALLOCATE(grad_no_cij,(3,natom,natom)) 854 ABI_ALLOCATE(fgrad_no_c,(2,3,natom,natom)) 855 ABI_ALLOCATE(cfgrad_no_c,(2,3,natom,natom)) 856 ABI_ALLOCATE(str_no_c,(6,natom,natom)) 857 ABI_ALLOCATE(str_alpha1,(6,natom)) 858 ABI_ALLOCATE(str_alpha2,(6,natom)) 859 grad_no_cij=zero ; str_no_c=zero 860 fgrad_no_c = zero ; cfgrad_no_c = zero 861 str_alpha1 = zero ; str_alpha2 = zero 862 863 re_arg = zero ; im_arg = zero 864 dmp6 = zero ; dmp8 = zero 865 e_no_c6 = zero ; e_no_c8 = zero 866 fdmp6 = zero ; fdmp8 = zero 867 grad6 = zero ; grad8 = zero 868 grad6_no_c6 = zero ; grad8_no_c8 = zero 869 hess6 = zero ; hess8 = zero 870 do 871 newshell=.false.;nshell=nshell+1 872 do is3=-nshell,nshell 873 do is2=-nshell,nshell 874 do is1=-nshell,nshell 875 if (nshell==1.or.abs(is3)==nshell.or.abs(is2)==nshell.or.abs(is1)==nshell) then 876 is(1) = is1 ; is(2) = is2 ; is(3) = is3 877 ! Computation of phase factor for discrete Fourier transform 878 ! if phonon at qphon is required 879 if (need_dynmat) then 880 arg= two_pi*dot_product(qphon,is) 881 re_arg=cos(arg) 882 im_arg=sin(arg) 883 end if 884 do ia=1,natom 885 itypat=typat(ia) 886 do ja=1,natom 887 jtypat=typat(ja) 888 r(:)=xred01(:,ia)-xred01(:,ja)-dble(is(:)) 889 rsq=dot_product(r,matmul(rmet,r)) 890 if (rsq>=tol16.and.rsq<rcut2) then 891 npairs=npairs+1;newshell=.true. 892 sfact6 = half*vdw_s6 ; sfact8 = half*vdw_s8 893 rr=sqrt(rsq); r6 = rr**six ; r8 = rr**eight 894 c6=vdw_c6(ia,ja) ; c8=vdw_c8(ia,ja) 895 vdw_q = three*vdw_q_dftd3(ivdw(itypat))*vdw_q_dftd3(ivdw(jtypat)) 896 r0=vdw_r0(itypat,jtypat) 897 ! Computation of e_vdw_dftd3 (case DFT+D3) 898 if (vdw_xc == 6) then 899 dmp6=six*(rr/(vdw_sr6*r0))**(-alpha6) 900 fdmp6=one/(one+dmp6) 901 dmp8=six*(rr/(vdw_sr8*r0))**(-alpha8) 902 fdmp8=one/(one+dmp8) 903 ! Contribution to energy 904 e_no_c6 = -sfact6*fdmp6/r6 ; e_no_c8 = -sfact8*fdmp8/r8 905 e_vdw_dftd3=e_vdw_dftd3+e_no_c6*c6 +e_no_c8*c8 906 ! Computation of e_vdw_dftd3 (case DFT+D3-BJ) 907 elseif (vdw_xc == 7) then 908 dmp = (vdw_a1*sqrt(vdw_q)+vdw_a2) 909 fdmp6 = one/(dmp**six+rr**six) 910 fdmp8 = one/(dmp**eight+rr**eight) 911 e_no_c6 = -sfact6*fdmp6 ; e_no_c8 = -sfact8*fdmp8 912 e_vdw_dftd3=e_vdw_dftd3-sfact6*c6*fdmp6-sfact8*c8*fdmp8 913 end if 914 ! Computation of the gradients (if required) 915 if (need_grad) then 916 if (vdw_xc == 6) then 917 gr6 = alpha6*dmp6*fdmp6**two 918 grad6_no_c6 = sfact6*(gr6-six*fdmp6)/r8 919 grad6 = grad6_no_c6*c6 920 gr8 = alpha8*dmp8*fdmp8**two 921 grad8_no_c8 = sfact8*(gr8-eight*fdmp8)/r8/rsq 922 grad8 = grad8_no_c8*c8 923 elseif (vdw_xc == 7) then 924 grad6_no_c6 = -sfact6*six*(fdmp6*rsq)**two 925 grad6 = grad6_no_c6*c6 926 grad8_no_c8 = -sfact8*eight*(fdmp8)**two*rsq**three 927 grad8 = grad8_no_c8*c8 928 end if 929 grad =grad6+grad8 930 grad_no_c = grad6_no_c6+grad8_no_c8*vdw_q 931 rcart=matmul(rprimd,r) 932 rred= matmul(transpose(rprimd),rcart) 933 ! Additional contribution due to c6(cn(r)) 934 ! Not yet multiply by dCN/drk and summed to reduce 935 ! computational time 936 e_no_c(ia,ja) = e_no_c(ia,ja)+(e_no_c6+vdw_q*e_no_c8) 937 ! Part related to alpha1ij/alpha2ij 938 e_alpha1(ia) = e_alpha1(ia)+(e_no_c6+vdw_q*e_no_c8)*dc6ri(ia,ja) 939 e_alpha2(ja) = e_alpha2(ja)+(e_no_c6+vdw_q*e_no_c8)*dc6rj(ia,ja) 940 ! Contribution to gradients wr to atomic displacement 941 ! (forces) 942 if (need_forces.and.ia/=ja) then 943 fred(:)=grad*rred(:) 944 do alpha=1,3 945 fred_vdw_dftd3(alpha,ia)=fred_vdw_dftd3(alpha,ia)-fred(alpha) 946 fred_vdw_dftd3(alpha,ja)=fred_vdw_dftd3(alpha,ja)+fred(alpha) 947 end do 948 elseif (need_stress) then 949 ! Computation of the DFT-D3 contribution to stress 950 vec(1:3)=rcart(1:3)*rcart(1:3); vec(4)=rcart(2)*rcart(3) 951 vec(5)=rcart(1)*rcart(3); vec(6)=rcart(1)*rcart(2) 952 do alpha=1,6 953 str_vdw_dftd3(alpha)=str_vdw_dftd3(alpha)-grad*vec(alpha) 954 end do 955 end if 956 ! Second derivative (if required) 957 if (need_hess) then 958 if (vdw_xc==6) then 959 hess6 = (grad6*(alpha6*fdmp6*dmp6-8.0_dp)+& 960 & sfact6*c6/r6*dmp6*((alpha6*fdmp6)**two)*& 961 & (fdmp6*dmp6-one)/rsq)/rsq 962 hess8 = (grad8*(alpha8*fdmp8*dmp8-10.0_dp)+& 963 & sfact8*c8/r8*dmp8*((alpha8*fdmp8)**two)*& 964 & (fdmp8*dmp8-one)/rsq)/rsq 965 elseif (vdw_xc==7) then 966 hess6 = -four*grad6*(three*rsq**two*fdmp6-one/rsq) 967 hess8 = -two*grad8*(eight*rsq**three*fdmp8-three/rsq) 968 end if 969 ! Contribution of d2C6 to the interatomic force constants 970 ! Not yet multiply by CN derivative and summed to reduce the scaling from O(N^3) to O(N^2) 971 hessij = hess6+hess8 972 ! Contribution of cross-derivative dC6 and grad 973 do alpha=1,3 974 grad_no_cij(alpha,ia,ja) = grad_no_cij(alpha,ia,ja) - grad_no_c*rcart(alpha) 975 end do 976 e_alpha3(ia) = e_alpha3(ia)+(e_no_c6+vdw_q*e_no_c8)*d2c6ri(ia,ja) 977 e_alpha4(ja) = e_alpha4(ja)+(e_no_c6+vdw_q*e_no_c8)*d2c6rj(ia,ja) 978 if (need_dynmat) then 979 ! Fourier transform of the partial contribution to the dispersion potential 980 fe_no_c(1,ia,ja) = fe_no_c(1,ia,ja)+(e_no_c6+vdw_q*e_no_c8)*re_arg 981 fe_no_c(2,ia,ja) = fe_no_c(2,ia,ja)+(e_no_c6+vdw_q*e_no_c8)*im_arg 982 do alpha=1,3 983 ! Fourier transform of the gradient (required for the IFCs) 984 fgrad_no_c(1,alpha,ia,ja) = fgrad_no_c(1,alpha,ia,ja)-grad_no_c*rcart(alpha)*re_arg 985 fgrad_no_c(2,alpha,ia,ja) = fgrad_no_c(2,alpha,ia,ja)-grad_no_c*rcart(alpha)*im_arg 986 ! Complex conjugated of the Fourier transform of the gradient 987 cfgrad_no_c(1,alpha,ia,ja) = cfgrad_no_c(1,alpha,ia,ja)-grad_no_c*rcart(alpha)*re_arg 988 cfgrad_no_c(2,alpha,ia,ja) = cfgrad_no_c(2,alpha,ia,ja)+grad_no_c*rcart(alpha)*im_arg 989 end do 990 ! Contribution to the IFCs (reciprocal space) of the 2nd derivative of e_no_c part 991 do alpha=1,3 992 do beta=1,3 993 rcart2(alpha,beta) = rcart(alpha)*rcart(beta) 994 end do 995 end do 996 if (ia/=ja) then 997 do alpha=1,3 998 dyn_vdw_dftd3(1,alpha,ja,:,ja) = dyn_vdw_dftd3(1,alpha,ja,:,ja) -& 999 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:)) 1000 dyn_vdw_dftd3(1,alpha,ia,:,ia) = dyn_vdw_dftd3(1,alpha,ia,:,ia) -& 1001 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:)) 1002 if (abs(re_arg)>tol12) then 1003 dyn_vdw_dftd3(1,alpha,ia,:,ja) = dyn_vdw_dftd3(1,alpha,ia,:,ja) +& 1004 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:))*re_arg 1005 dyn_vdw_dftd3(1,alpha,ja,:,ia) = dyn_vdw_dftd3(1,alpha,ja,:,ia) +& 1006 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:))*re_arg 1007 end if 1008 if (abs(im_arg)>tol12) then 1009 dyn_vdw_dftd3(2,alpha,ia,:,ja) = dyn_vdw_dftd3(2,alpha,ia,:,ja) +& 1010 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:))*im_arg 1011 dyn_vdw_dftd3(2,alpha,ja,:,ia) = dyn_vdw_dftd3(2,alpha,ja,:,ia) -& 1012 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:))*im_arg 1013 end if 1014 end do 1015 else ! ia==ja 1016 do alpha=1,3 1017 if (abs(re_arg-one)>tol12) then 1018 dyn_vdw_dftd3(1,alpha,ia,:,ia) = dyn_vdw_dftd3(1,alpha,ia,:,ia) -& 1019 & two*(hessij*rcart2(alpha,:)+grad*mcart(alpha,:))*(one-re_arg) 1020 end if 1021 end do 1022 end if 1023 end if 1024 ! Now compute the contribution to the elastic constants !!! Still under development 1025 if (need_elast) then 1026 vec(1:3)=rcart(1:3)*rcart(1:3); vec(4)=rcart(2)*rcart(3) 1027 vec(5)=rcart(1)*rcart(3); vec(6)=rcart(1)*rcart(2) 1028 str_no_c(:,ia,ja)=str_no_c(:,ia,ja)-grad_no_c*vec(:) 1029 str_alpha1(:,ia)=str_alpha1(:,ia)-dc6ri(ia,ja)*grad_no_c*vec(:) 1030 str_alpha2(:,ja)=str_alpha2(:,ja)-dc6rj(ia,ja)*grad_no_c*vec(:) 1031 ! Contribution to elastic constants of DFT-D3 dispersion potential (no C6 derivative) 1032 do alpha=1,6 1033 ii = voigt1(alpha) ; jj=voigt2(alpha) 1034 do beta=1,6 1035 kk = voigt1(beta) ; ll=voigt2(beta) 1036 elt_vdw_dftd3(alpha,beta) = elt_vdw_dftd3(alpha,beta)-hessij*vec(alpha)*vec(beta) 1037 if (ii==kk) elt_vdw_dftd3(alpha,beta) = elt_vdw_dftd3(alpha,beta)-half*grad*rcart(jj)*rcart(ll) 1038 if (jj==kk) elt_vdw_dftd3(alpha,beta) = elt_vdw_dftd3(alpha,beta)-half*grad*rcart(ii)*rcart(ll) 1039 if (ii==ll) elt_vdw_dftd3(alpha,beta) = elt_vdw_dftd3(alpha,beta)-half*grad*rcart(jj)*rcart(kk) 1040 if (jj==ll) elt_vdw_dftd3(alpha,beta) = elt_vdw_dftd3(alpha,beta)-half*grad*rcart(ii)*rcart(kk) 1041 end do 1042 end do 1043 ! Contribution to internal strain of DFT-D3 dispersion potential (no C6 derivative) 1044 if (ia/=ja) then 1045 index_ia = 6+3*(ia-1) 1046 index_ja = 6+3*(ja-1) 1047 do alpha=1,6 1048 do beta=1,3 1049 ii = voigt1(alpha) ; jj=voigt2(alpha) 1050 elt_vdw_dftd3(index_ia+beta,alpha)=elt_vdw_dftd3(index_ia+beta,alpha)-& 1051 & hessij*vec(alpha)*rcart(beta) 1052 elt_vdw_dftd3(index_ja+beta,alpha)=elt_vdw_dftd3(index_ja+beta,alpha)+& 1053 & hessij*vec(alpha)*rcart(beta) 1054 if (ii==beta) then 1055 elt_vdw_dftd3(index_ia+beta,alpha)=elt_vdw_dftd3(index_ia+beta,alpha)-grad*rcart(jj) 1056 elt_vdw_dftd3(index_ja+beta,alpha)=elt_vdw_dftd3(index_ja+beta,alpha)+grad*rcart(jj) 1057 end if 1058 if (jj==beta) then 1059 elt_vdw_dftd3(index_ia+beta,alpha)=elt_vdw_dftd3(index_ia+beta,alpha)-grad*rcart(ii) 1060 elt_vdw_dftd3(index_ja+beta,alpha)=elt_vdw_dftd3(index_ja+beta,alpha)+grad*rcart(ii) 1061 end if 1062 end do ! Direction beta 1063 end do ! Strain alpha 1064 end if ! ia/=ja 1065 end if ! Need elastic constant 1066 end if ! Need hessian 1067 end if ! Need gradient 1068 end if ! Tolerance 1069 end do ! Loop over atom j 1070 end do ! Loop over atom i 1071 end if ! Triple loop over cell replicas in shell 1072 end do ! Is1 1073 end do ! Is2 1074 end do ! Is3 1075 if(.not.newshell) exit ! Check if new shell must be calculated 1076 end do ! Loop over shell 1077 ABI_ALLOCATE(temp_prod,(2,natom)) 1078 if (need_grad) then 1079 ! Additional contribution to force due dc6_drk 1080 if (need_forces) then 1081 do ka=1,natom 1082 do ia=1,natom 1083 !do ja=1,natom 1084 !fred_vdw_dftd3(:,ka) = fred_vdw_dftd3(:,ka)+e_no_c(ia,ja)*(& 1085 !& !dcn(:,ia,ka)*dc6ri(ia,ja)+dcn(:,ja,ka)*dc6rj(ia,ja)) 1086 !end do 1087 fred_vdw_dftd3(:,ka) = fred_vdw_dftd3(:,ka)+e_alpha1(ia)*dcn(:,ia,ka)+& 1088 & e_alpha2(ia)*dcn(:,ia,ka) 1089 end do 1090 end do 1091 elseif (need_stress) then 1092 do ia=1,natom 1093 !do ja=1,natom 1094 ! str_vdw_dftd3(:) = str_vdw_dftd3(:)+e_no_c(ia,ja)*(str_dcn(:,ia)*& 1095 !& ! dc6ri(ia,ja)+str_dcn(:,ja)*dc6rj(ia,ja)) 1096 !end do 1097 str_vdw_dftd3(:) = str_vdw_dftd3(:)+e_alpha1(ia)*str_dcn(:,ia)+& 1098 & e_alpha2(ia)*str_dcn(:,ia) 1099 end do 1100 end if ! Optimization 1101 ! If dynmat is required, add all the terms related to dc6, d2c6, ... 1102 if (need_hess) then 1103 if (need_dynmat) then 1104 do ka=1,natom 1105 do la =1,natom 1106 do alpha=1,3 1107 do beta=1,3 1108 do ia=1,natom 1109 ! Add the second derivative of C6 contribution to the dynamical matrix 1110 ! First, add the second derivative of CN-related term 1111 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1112 & (e_alpha1(ia)+e_alpha2(ia))*d2cn(:,alpha,ka,beta,la,ia) 1113 ! Then the term related to dCNi/dr*dCNi/dr and dCNj/dr*dCNj/dr 1114 call comp_prod(cfdcn(:,alpha,ia,ka),fdcn(:,beta,ia,la),temp_comp) 1115 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1116 & (e_alpha3(ia)+e_alpha4(ia))*temp_comp(:) 1117 ! Add the cross derivative of fdmp/rij**6 and C6 contribution to the dynamical matrix 1118 ! !!!! The products are kind of tricky... 1119 ! First, add the dCNk/drl gradik and dCNk/drl gradjk terms... 1120 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1121 & cfdcn(:,alpha,la,ka)*grad_no_cij(beta,la,ia)*(dc6ri(la,ia)+dc6rj(ia,la)) 1122 ! Then the dCNk/drl gradjk and dCNk/drl gradik terms... 1123 call comp_prod(cfdcn(:,alpha,ia,ka),cfgrad_no_c(:,beta,la,ia),temp_comp) 1124 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1125 & temp_comp(:)*(dc6ri(ia,la)+dc6rj(la,ia)) 1126 ! Here the symmetrical term (for dCNl/drk) are added... 1127 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1128 & fdcn(:,beta,ka,la)*grad_no_cij(alpha,ka,ia)*(dc6ri(ka,ia)+dc6rj(ia,ka)) 1129 call comp_prod(fdcn(:,beta,ia,la),fgrad_no_c(:,alpha,ka,ia),temp_comp) 1130 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1131 & temp_comp(:)*(dc6ri(ia,ka)+dc6rj(ka,ia)) 1132 end do ! ia 1133 end do ! alpha 1134 end do ! beta 1135 end do ! la 1136 do alpha=1,3 1137 temp_prod(:,:) = zero 1138 do ja=1,natom 1139 do ia=1,natom 1140 ! Finally the cross derivative dCNi/dr dCNj/dr 1141 temp_comp2(:) = d2c6rirj(ia,ja)*fe_no_c(:,ia,ja) 1142 call comp_prod(cfdcn(:,alpha,ia,ka),temp_comp2,temp_comp) 1143 temp_prod(:,ja) = temp_prod(:,ja)+temp_comp(:) 1144 end do 1145 do la = 1,natom 1146 do beta=1,3 1147 call comp_prod(fdcn(:,beta,ja,la),temp_prod(:,ja),temp_comp2) 1148 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1149 & two*temp_comp2 1150 end do ! beta 1151 end do ! la 1152 end do ! ja 1153 end do ! alpha 1154 end do !ka 1155 ! Transformation from cartesian coordinates to reduced coordinates 1156 do ka=1,natom 1157 do la=1,natom 1158 do kk=1,2 1159 do alpha=1,3 1160 vec(1:3)=dyn_vdw_dftd3(kk,1:3,ka,alpha,la) 1161 call d3_cart2red(vec) 1162 dyn_vdw_dftd3(kk,1:3,ka,alpha,la)=vec(1:3) 1163 end do 1164 do alpha=1,3 1165 vec(1:3)=dyn_vdw_dftd3(kk,alpha,ka,1:3,la) 1166 call d3_cart2red(vec) 1167 dyn_vdw_dftd3(kk,alpha,ka,1:3,la)=vec(1:3) 1168 end do ! alpha 1169 end do ! real/im 1170 end do ! Atom la 1171 end do ! Atom ka 1172 end if ! Boolean dynamical matrix 1173 if (need_elast) then 1174 do ia=1,natom 1175 index_ia = 6+3*(ia-1) 1176 do alpha=1,6 1177 ! Add the second derivative of C6 contribution to the elastic tensor 1178 ! First, the second derivative of CN with strain 1179 elt_vdw_dftd3(alpha,:) = elt_vdw_dftd3(alpha,:)+(& 1180 & e_alpha1(ia)+e_alpha2(ia))*elt_cn(alpha,:,ia) 1181 ! Then, the derivative product dCNi/deta dCNi/deta 1182 elt_vdw_dftd3(alpha,:) = elt_vdw_dftd3(alpha,:)+(& 1183 & e_alpha3(ia)+e_alpha4(ia))*str_dcn(alpha,ia)*str_dcn(:,ia) 1184 ! Then, the dCNi/deta dCNj/deta 1185 do ja=1,natom 1186 elt_vdw_dftd3(alpha,:)=elt_vdw_dftd3(alpha,:)+two*e_no_c(ia,ja)*& 1187 & d2c6rirj(ia,ja)*str_dcn(alpha,ia)*str_dcn(:,ja) 1188 end do 1189 ! Add the cross derivative of fij and C6 contribution to the elastic tensor 1190 elt_vdw_dftd3(alpha,:)=elt_vdw_dftd3(alpha,:)+str_alpha1(:,ia)*str_dcn(alpha,ia)+str_alpha1(alpha,ia)*str_dcn(:,ia)& 1191 & +str_dcn(alpha,ia)*str_alpha2(:,ia)+str_dcn(:,ia)*str_alpha2(alpha,ia) 1192 end do 1193 do alpha=1,6 1194 ii = voigt1(alpha) ; jj=voigt2(alpha) 1195 do ka=1,natom 1196 index_ka = 6+3*(ka-1) 1197 do beta=1,3 1198 ! Add the second derivative of C6 contribution to the internal strains 1199 ii = voigt1(alpha) ; jj=voigt2(alpha) 1200 ! Second derivative of CN 1201 elt_vdw_dftd3(index_ka+beta,alpha)=elt_vdw_dftd3(index_ka+beta,alpha)+(e_alpha1(ia)+e_alpha2(ia))*& 1202 & elt_cn(index_ka+beta,alpha,ia) 1203 ! Cross-derivatives of CN 1204 elt_vdw_dftd3(index_ka+beta,alpha)=elt_vdw_dftd3(index_ka+beta,alpha)+(e_alpha3(ia)+e_alpha4(ia))*& 1205 & dcn_cart(beta,ia,ka)*str_dcn(alpha,ia) !OK 1206 do ja=1,natom 1207 index_ja = 6+3*(ja-1) 1208 elt_vdw_dftd3(index_ka+beta,alpha)=elt_vdw_dftd3(index_ka+beta,alpha)+two*d2c6rirj(ia,ja)*e_no_c(ia,ja)*& 1209 & dcn_cart(beta,ia,ka)*str_dcn(alpha,ja) 1210 end do 1211 elt_vdw_dftd3(index_ka+beta,alpha)=elt_vdw_dftd3(index_ka+beta,alpha)+(str_alpha1(alpha,ia)+str_alpha2(alpha,ia))*& 1212 & dcn_cart(beta,ia,ka) 1213 elt_vdw_dftd3(index_ka+beta,alpha)=elt_vdw_dftd3(index_ka+beta,alpha)+grad_no_cij(beta,ka,ia)*& 1214 & (dc6ri(ia,ka)*str_dcn(alpha,ia)+dc6rj(ia,ka)*str_dcn(alpha,ka))-grad_no_cij(beta,ia,ka)*& 1215 & (dc6ri(ka,ia)*str_dcn(alpha,ka)+dc6rj(ka,ia)*str_dcn(alpha,ia)) 1216 end do ! beta 1217 end do ! ka 1218 end do ! alpha 1219 end do ! ia 1220 do alpha=1,6 1221 index_ia=6 1222 do ia=1,natom 1223 elt_vdw_dftd3(index_ia+1:index_ia+3,alpha)=matmul(transpose(rprimd),elt_vdw_dftd3(index_ia+1:index_ia+3,alpha)) 1224 index_ia=index_ia+3 1225 end do ! Atom ia 1226 end do ! Strain alpha 1227 end if ! Boolean elastic tensor 1228 end if ! Boolean hessian 1229 end if ! Boolean need_gradient 1230 ABI_DEALLOCATE(temp_prod) 1231 write(msg,'(3a)')& 1232 & ' ...Done.' 1233 call wrtout(std_out,msg,'COLL') 1234 1235 !print *, 'Evdw', e_vdw_dftd3 1236 !if (need_forces) print *, 'fvdw', fred_vdw_dftd3(3,:) 1237 !if (need_stress) print *, 'strvdw', str_vdw_dftd3(6) 1238 !if (need_elast) print *, 'Elast(3,3,3,3)', elt_vdw_dftd3(3,3) 1239 !if (need_elast) print *, 'Elast(3,3,3,3)', elt_vdw_dftd3(6,6) 1240 !if (need_elast) print *, 'Elast(3,3,3,3)', elt_vdw_dftd3(3,6) 1241 !if (need_elast) print *, 'Internal(1,3,3)', elt_vdw_dftd3(9,3) 1242 !if (need_elast) print *, 'Internal(1,3,6)', elt_vdw_dftd3(9,6) 1243 !if (need_dynmat) print *, 'Dynmat(1,3,:,3,:)', dyn_vdw_dftd3(1,3,:,3,:) 1244 !if (need_dynmat) print *, 'Dynmat(1,3,:,3,:)', dyn_vdw_dftd3(1,2,:,1,:) 1245 1246 !--------------------------------------------- 1247 ! Computation of the 3 body term (if required) 1248 !--------------------------------------------- 1249 1250 e_3bt=zero 1251 if (bol_3bt) then 1252 if (need_grad) then 1253 ABI_ALLOCATE(e3bt_ij,(natom,natom)) 1254 ABI_ALLOCATE(e3bt_jk,(natom,natom)) 1255 ABI_ALLOCATE(e3bt_ki,(natom,natom)) 1256 ABI_ALLOCATE(fred_vdw_3bt,(3,natom)) 1257 e3bt_ij=zero; e3bt_jk=zero; e3bt_ki=zero 1258 if (need_forces) then 1259 fred_vdw_3bt = zero 1260 elseif (need_stress) then 1261 str_3bt=zero 1262 end if 1263 end if 1264 nshell_3bt(1) = int(0.5+rcut9/sqrt(rmet(1,1)+rmet(2,1)+rmet(3,1))) 1265 nshell_3bt(2) = int(0.5+rcut9/sqrt(rmet(1,2)+rmet(2,2)+rmet(3,2))) 1266 nshell_3bt(3) = int(0.5+rcut9/sqrt(rmet(1,3)+rmet(2,3)+rmet(3,3))) 1267 1268 do is3 = -nshell_3bt(3),nshell_3bt(3) 1269 do is2 = -nshell_3bt(2),nshell_3bt(2) 1270 do is1 = -nshell_3bt(1),nshell_3bt(1) 1271 is(1) = is1 ; is(2)=is2 ; is(3) = is3 1272 do alpha=1,3 1273 jmin(alpha) = max(-nshell_3bt(alpha), -nshell_3bt(alpha)+is(alpha)) 1274 jmax(alpha) = min(nshell_3bt(alpha), nshell_3bt(alpha)+is(alpha)) 1275 end do 1276 do js3=jmin(3),jmax(3) 1277 do js2=jmin(2),jmax(2) 1278 do js1=jmin(1),jmax(1) 1279 js(1) = js1 ; js(2)=js2 ; js(3) = js3 1280 do ia=1,natom 1281 itypat=typat(ia) 1282 do ja=1,ia 1283 jtypat=typat(ja) 1284 do ka=1,ja 1285 ktypat=typat(ka) 1286 rij(:) = xred01(:,ia)-xred01(:,ja)-dble(is(:)) 1287 rsqij = dot_product(rij(:),matmul(rmet,rij(:))) 1288 rrij = dsqrt(rsqij) 1289 rjk(:) = xred01(:,ja)-xred01(:,ka)+dble(is(:))-dble(js(:)) 1290 rsqjk = dot_product(rjk(:),matmul(rmet,rjk(:))) 1291 rrjk = dsqrt(rsqjk) 1292 rki(:) = xred01(:,ka)-xred01(:,ia)+dble(js(:)) 1293 rsqki = dot_product(rki(:),matmul(rmet,rki(:))) 1294 rrki = dsqrt(rsqki) 1295 if (rsqij>=tol16.and.rsqjk>=tol16.and.rsqki>=tol16) then 1296 rmean = (rrij*rrjk*rrki)**third 1297 if (rrij>rcut9.or.rrjk>rcut9.or.rrki>rcut9) cycle 1298 sfact9=vdw_s6 1299 if (ia==ja.and.ja==ka) sfact9 = sixth*sfact9 1300 if (ia==ja.and.ja/=ka) sfact9 = half*sfact9 1301 if (ia/=ja.and.ja==ka) sfact9 = half*sfact9 1302 rijk = one/(rrij*rrjk*rrki) 1303 dmp9 = six*(rmean*vdw_sr9*r0ijk(itypat,jtypat,ktypat))**(-alpha8) 1304 fdmp9 = one/(one+dmp9) 1305 cosa = half*rrjk*(rsqij+rsqki-rsqjk)*rijk 1306 cosb = half*rrki*(rsqij+rsqjk-rsqki)*rijk 1307 cosc = half*rrij*(rsqjk+rsqki-rsqij)*rijk 1308 ang = one+three*cosa*cosb*cosc 1309 temp = sfact9*rijk*rijk*rijk 1310 temp2 = temp*fdmp9*ang 1311 ! Contribution to energy 1312 e_3bt = e_3bt-temp2*vdw_c9(ia,ja,ka) !*temp2 1313 e3bt_ij(ia,ja) = e3bt_ij(ia,ja)-temp2*vdw_c9(ia,ja,ka) 1314 e3bt_jk(ja,ka) = e3bt_jk(ja,ka)-temp2*vdw_c9(ia,ja,ka) 1315 e3bt_ki(ka,ia) = e3bt_ki(ka,ia)-temp2*vdw_c9(ia,ja,ka) 1316 if (need_grad) then 1317 dfdmp = third*alpha8*fdmp9*fdmp9*dmp9 1318 if (ia/=ja.or.need_stress) then 1319 d1_r3drij = -three*rrki*rrjk 1320 dcosa_r3drij = (rrij-two*cosa*rrki)*rrjk 1321 dcosb_r3drij = (rrij-two*cosb*rrjk)*rrki 1322 dcosc_r3drij =-(rsqij+cosc*rrjk*rrki) 1323 dfdmp_drij = dfdmp*rrki*rrjk 1324 d_drij = vdw_c9(ia,ja,ka)*temp*rijk*((d1_r3drij+three*& 1325 & (cosb*cosc*dcosa_r3drij+cosa*cosc*dcosb_r3drij+cosa*cosb*dcosc_r3drij))*& 1326 & fdmp9+dfdmp_drij*ang)*rijk*rrjk*rrki 1327 rcartij=matmul(rprimd,rij) 1328 end if 1329 if (ja/=ka.or.need_stress) then 1330 d1_r3drjk = -three*rrij*rrki 1331 dcosa_r3drjk =-(rsqjk+cosa*rrij*rrki) 1332 dcosb_r3drjk = (rrjk-two*cosb*rrij)*rrki 1333 dcosc_r3drjk = (rrjk-two*cosc*rrki)*rrij 1334 dfdmp_drjk = dfdmp*rrij*rrki 1335 d_drjk = vdw_c9(ia,ja,ka)*temp*rijk*((d1_r3drjk+three*& 1336 & (cosb*cosc*dcosa_r3drjk+cosa*cosc*dcosb_r3drjk+cosa*cosb*dcosc_r3drjk))*& 1337 & fdmp9+dfdmp_drjk*ang)*rijk*rrij*rrki 1338 rcartjk=matmul(rprimd,rjk) 1339 end if 1340 if (ka/=ia.or.need_stress) then 1341 d1_r3drki = -three*rrjk*rrij 1342 dcosa_r3drki = (rrki-two*cosa*rrij)*rrjk 1343 dcosb_r3drki =-(rsqki+cosb*rrij*rrjk) 1344 dcosc_r3drki = (rrki-two*cosc*rrjk)*rrij 1345 dfdmp_drki = dfdmp*rrij*rrjk 1346 d_drki = vdw_c9(ia,ja,ka)*temp*rijk*((d1_r3drki+three*& 1347 & (cosb*cosc*dcosa_r3drki+cosa*cosc*dcosb_r3drki+cosa*cosb*dcosc_r3drki))*& 1348 & fdmp9+dfdmp_drki*ang)*rijk*rrij*rrjk 1349 rcartki=matmul(rprimd,rki) 1350 end if 1351 ! Contribution to gradients wr to atomic displacement 1352 ! (forces) 1353 if (need_forces) then 1354 if (ia/=ja) fredij=d_drij*matmul(transpose(rprimd),rcartij) 1355 if (ja/=ka) fredjk=d_drjk*matmul(transpose(rprimd),rcartjk) 1356 if (ka/=ia) fredki=d_drki*matmul(transpose(rprimd),rcartki) 1357 if (ia/=ja.and.ka/=ia) then 1358 fred_vdw_3bt(:,ia)=fred_vdw_3bt(:,ia)-fredij(:)+fredki(:) 1359 fred_vdw_3bt(:,ja)=fred_vdw_3bt(:,ja)+fredij(:)-fredjk(:) 1360 fred_vdw_3bt(:,ka)=fred_vdw_3bt(:,ka)-fredki(:)+fredjk(:) 1361 else if (ia==ja.and.ia/=ka) then 1362 fred_vdw_3bt(:,ia)=fred_vdw_3bt(:,ia)+fredki(:) 1363 fred_vdw_3bt(:,ja)=fred_vdw_3bt(:,ja)-fredjk(:) 1364 fred_vdw_3bt(:,ka)=fred_vdw_3bt(:,ka)-fredki(:)+fredjk(:) 1365 elseif (ia==ka.and.ia/=ja) then 1366 fred_vdw_3bt(:,ia)=fred_vdw_3bt(:,ia)-fredij(:) 1367 fred_vdw_3bt(:,ja)=fred_vdw_3bt(:,ja)+fredij(:)-fredjk(:) 1368 fred_vdw_3bt(:,ka)=fred_vdw_3bt(:,ka)+fredjk(:) 1369 elseif (ja==ka.and.ia/=ja) then 1370 fred_vdw_3bt(:,ia)=fred_vdw_3bt(:,ia)-fredij(:)+fredki(:) 1371 fred_vdw_3bt(:,ja)=fred_vdw_3bt(:,ja)+fredij(:) 1372 fred_vdw_3bt(:,ka)=fred_vdw_3bt(:,ka)-fredki(:) 1373 end if 1374 end if 1375 ! Contribution to stress tensor 1376 if (need_stress) then 1377 vecij(1:3)=rcartij(1:3)*rcartij(1:3); vecij(4)=rcartij(2)*rcartij(3) 1378 vecij(5)=rcartij(1)*rcartij(3); vecij(6)=rcartij(1)*rcartij(2) 1379 vecjk(1:3)=rcartjk(1:3)*rcartjk(1:3); vecjk(4)=rcartjk(2)*rcartjk(3) 1380 vecjk(5)=rcartjk(1)*rcartjk(3); vecjk(6)=rcartjk(1)*rcartjk(2) 1381 vecki(1:3)=rcartki(1:3)*rcartki(1:3); vecki(4)=rcartki(2)*rcartki(3) 1382 vecki(5)=rcartki(1)*rcartki(3); vecki(6)=rcartki(1)*rcartki(2) 1383 str_3bt(:)=str_3bt(:)-d_drij*vecij(:)-d_drjk*vecjk(:)-d_drki*vecki(:) !-str_3bt_dcn9(:) 1384 end if 1385 end if ! Optimization 1386 end if ! Tolerance 1387 end do ! Loop over atom k 1388 end do ! Loop over atom j 1389 end do ! Loop over atom i 1390 end do ! j3 1391 end do ! j2 1392 end do ! j1 1393 end do 1394 end do 1395 end do 1396 if (need_forces) then 1397 do ia=1,natom 1398 do ja=1,natom 1399 do la=1,natom 1400 fred_vdw_3bt(:,la) = fred_vdw_3bt(:,la)+e3bt_ij(ia,ja)*(dc9ijri(ia,ja)*dcn(:,ia,la)+dc9ijrj(ia,ja)*dcn(:,ja,la)) 1401 fred_vdw_3bt(:,la) = fred_vdw_3bt(:,la)+e3bt_jk(ia,ja)*(dc9ijri(ia,ja)*dcn(:,ia,la)+dc9ijrj(ia,ja)*dcn(:,ja,la)) 1402 fred_vdw_3bt(:,la) = fred_vdw_3bt(:,la)+e3bt_ki(ia,ja)*(dc9ijri(ia,ja)*dcn(:,ia,la)+dc9ijrj(ia,ja)*dcn(:,ja,la)) 1403 end do 1404 end do 1405 end do 1406 elseif (need_stress) then 1407 do ia=1,natom 1408 do ja=1,natom 1409 str_3bt(:) = str_3bt(:)+e3bt_ij(ia,ja)*(dc9ijri(ia,ja)*str_dcn(:,ia)+dc9ijrj(ia,ja)*str_dcn(:,ja)) 1410 str_3bt(:) = str_3bt(:)+e3bt_jk(ia,ja)*(dc9ijri(ia,ja)*str_dcn(:,ia)+dc9ijrj(ia,ja)*str_dcn(:,ja)) 1411 str_3bt(:) = str_3bt(:)+e3bt_ki(ia,ja)*(dc9ijri(ia,ja)*str_dcn(:,ia)+dc9ijrj(ia,ja)*str_dcn(:,ja)) 1412 end do 1413 end do 1414 end if 1415 e_vdw_dftd3 = e_vdw_dftd3+e_3bt 1416 if (need_forces) fred_vdw_dftd3= fred_vdw_dftd3+fred_vdw_3bt 1417 if (need_stress) str_vdw_dftd3 = str_vdw_dftd3+str_3bt 1418 ABI_DEALLOCATE(dc9ijri) 1419 ABI_DEALLOCATE(dc9ijrj) 1420 ABI_DEALLOCATE(e3bt_ij) 1421 ABI_DEALLOCATE(e3bt_jk) 1422 ABI_DEALLOCATE(e3bt_ki) 1423 ABI_DEALLOCATE(vdw_c9) 1424 ABI_DEALLOCATE(r0ijk) 1425 ABI_DEALLOCATE(fred_vdw_3bt) 1426 end if 1427 if (need_stress) str_vdw_dftd3=str_vdw_dftd3/ucvol 1428 1429 !Printing 1430 if (prtvol>0) then 1431 write(msg,'(2a)') ch10,& 1432 & ' --------------------------------------------------------------' 1433 call wrtout(std_out,msg,'COLL') 1434 if (vdw_xc==6) then 1435 write(msg,'(3a)') & 1436 & ' Van der Waals DFT-D3 semi-empirical dispersion potential as',ch10,& 1437 & ' proposed by Grimme et al., J. Chem. Phys. 132, 154104 (2010)' 1438 call wrtout(std_out,msg,'COLL') 1439 elseif (vdw_xc==7) then 1440 write(msg,'(5a)') & 1441 & ' Van der Waals DFT-D3 semi-empirical dispersion potential ' ,ch10,& 1442 & ' with Becke-Jonhson (BJ) refined by Grimme et al. J. ',ch10,& 1443 & ' Comput. Chem. 32, 1456 (2011) ' 1444 call wrtout(std_out,msg,'COLL') 1445 end if 1446 if (natom<5) then 1447 write(msg,'(3a)')& 1448 & ' Pair C6 (a.u.) C8 (a.u.) R0 (Ang) ',ch10,& 1449 & ' ---------------------------------------------------------------' 1450 call wrtout(std_out,msg,'COLL') 1451 do ia=1,natom 1452 do ja=1,ia 1453 itypat = typat(ia) ; jtypat = typat(ja) 1454 call atomdata_from_znucl(atom1,znucl(itypat)) 1455 call atomdata_from_znucl(atom2,znucl(jtypat)) 1456 write(msg,'(4X,2a,i2,3a,i2,1a,1X,es12.4,4X,es12.4,4X,es12.4,1X)') & 1457 atom1%symbol,'(',ia,')-',atom2%symbol,'(',ja,')', & 1458 vdw_c6(ia,ja), vdw_c8(ia,ja),& 1459 vdw_r0(itypat,jtypat) 1460 call wrtout(std_out,msg,'COLL') 1461 end do 1462 end do 1463 end if 1464 write(msg, '(3a,f6.3,a,f6.3)') & 1465 & ' ---------------------------------------------------------------',ch10,& 1466 & ' Scaling factors: s6 = ', vdw_s6,', s8 = ',vdw_s8 1467 call wrtout(std_out,msg,'COLL') 1468 if (vdw_xc==6) then 1469 write(msg,'(a,f6.3,a,f6.3)') & 1470 & ' Damping parameters: sr6 = ', vdw_sr6,', sr8 = ',vdw_sr8 1471 call wrtout(std_out,msg,'COLL') 1472 elseif (vdw_xc==7) then 1473 write(msg,'(a,f6.3,a,f6.3)') & 1474 & ' Damping parameters: a1 = ', vdw_a1, ', a2 = ', vdw_a2 1475 call wrtout(std_out,msg,'COLL') 1476 end if 1477 write(msg,'(a,es12.5,3a,i8,2a,es12.5,1a)') & 1478 & ' Cut-off radius = ',rcut,' Bohr',ch10,& 1479 & ' Number of pairs contributing = ',npairs,ch10,& 1480 & ' DFT-D3 (no 3-body) energy contribution = ',e_vdw_dftd3-e_3bt,' Ha' 1481 call wrtout(std_out,msg,'COLL') 1482 if (bol_3bt) then 1483 write(msg,'(6a,i5,2a,es20.11,3a,es20.11,1a)')ch10,& 1484 & ' ---------------------------------------------------------------',ch10,& 1485 & ' 3-Body Term Contribution:', ch10,& 1486 & ' Number of shells considered = ', nshell, ch10,& 1487 & ' Additional 3-body contribution = ', e_3bt, ' Ha',ch10,& 1488 & ' Total E (2-body and 3-body) = ', e_vdw_dftd3, 'Ha' 1489 call wrtout(std_out,msg,'COLL') 1490 end if 1491 write(msg,'(2a)')& 1492 & ' ----------------------------------------------------------------',ch10 1493 call wrtout(std_out,msg,'COLL') 1494 end if 1495 ABI_DEALLOCATE(ivdw) 1496 ABI_DEALLOCATE(xred01) 1497 ABI_DEALLOCATE(vdw_r0) 1498 ABI_DEALLOCATE(fe_no_c) 1499 ABI_DEALLOCATE(e_no_c) 1500 ABI_DEALLOCATE(e_alpha1) 1501 ABI_DEALLOCATE(e_alpha2) 1502 ABI_DEALLOCATE(e_alpha3) 1503 ABI_DEALLOCATE(e_alpha4) 1504 ABI_DEALLOCATE(grad_no_cij) 1505 ABI_DEALLOCATE(fgrad_no_c) 1506 ABI_DEALLOCATE(cfgrad_no_c) 1507 ABI_DEALLOCATE(vdw_c6) 1508 ABI_DEALLOCATE(vdw_c8) 1509 ABI_DEALLOCATE(dc6ri) 1510 ABI_DEALLOCATE(dc6rj) 1511 ABI_DEALLOCATE(d2c6ri) 1512 ABI_DEALLOCATE(d2c6rj) 1513 ABI_DEALLOCATE(d2c6rirj) 1514 ABI_DEALLOCATE(cn) 1515 ABI_DEALLOCATE(d2cn) 1516 ABI_DEALLOCATE(dcn) 1517 ABI_DEALLOCATE(fdcn) 1518 ABI_DEALLOCATE(cfdcn) 1519 ABI_DEALLOCATE(str_dcn) 1520 ABI_DEALLOCATE(elt_cn) 1521 ABI_DEALLOCATE(str_no_c) 1522 ABI_DEALLOCATE(str_alpha1) 1523 ABI_DEALLOCATE(str_alpha2) 1524 ABI_DEALLOCATE(dcn_cart) 1525 DBG_EXIT("COLL") 1526 1527 contains 1528 !! ***