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