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 141 ! s6 parameters (BJ case) 142 real(dp),parameter :: vdwbj_s6_b2gpplyp=0.560_dp, vdwbj_s6_ptpss=0.750_dp 143 real(dp),parameter :: vdwbj_s6_b2plyp=0.640_dp, vdwbj_s6_dsdblyp=0.500_dp 144 real(dp),parameter :: vdwbj_s6_pwpb95=0.820_dp 145 ! s8 parameters (BJ case) 146 real(dp),parameter :: vdwbj_s8_b1b95=1.4507_dp, vdwbj_s8_b2gpplyp=0.2597_dp 147 real(dp),parameter :: vdwbj_s8_b3pw91=2.8524_dp, vdwbj_s8_bhlyp=1.0354_dp 148 real(dp),parameter :: vdwbj_s8_bmk=2.0860_dp, vdwbj_s8_bop=3.295_dp 149 real(dp),parameter :: vdwbj_s8_bpbe=4.0728_dp, vdwbj_s8_camb3lyp=2.0674_dp 150 real(dp),parameter :: vdwbj_s8_lcwpbe=1.8541_dp, vdwbj_s8_mpw1b95=1.0508_dp 151 real(dp),parameter :: vdwbj_s8_mpwb1k=0.9499_dp, vdwbj_s8_mpwlyp=2.0077_dp 152 real(dp),parameter :: vdwbj_s8_olyp=2.6205_dp, vdwbj_s8_opbe=3.3816_dp 153 real(dp),parameter :: vdwbj_s8_otpss=2.7495_dp, vdwbj_s8_pbe38=1.4623_dp 154 real(dp),parameter :: vdwbj_s8_pbesol=2.9491_dp, vdwbj_s8_ptpss=0.2804_dp 155 real(dp),parameter :: vdwbj_s8_pwb6k=0.9383_dp, vdwbj_s8_revssb=0.4389_dp 156 real(dp),parameter :: vdwbj_s8_ssb=-0.1744_dp, vdwbj_s8_tpssh=2.2382_dp 157 real(dp),parameter :: vdwbj_s8_hcth120=1.0821_dp, vdwbj_s8_b2plyp=0.9147_dp 158 real(dp),parameter :: vdwbj_s8_b3lyp=1.9889_dp, vdwbj_s8_b97d=2.2609_dp 159 real(dp),parameter :: vdwbj_s8_blyp=2.6996_dp, vdwbj_s8_bp86=3.2822_dp 160 real(dp),parameter :: vdwbj_s8_dsdblyp=0.2130_dp, vdwbj_s8_pbe0=1.2177_dp 161 real(dp),parameter :: vdwbj_s8_pbe=0.7875_dp, vdwbj_s8_pw6b95=0.7257_dp 162 real(dp),parameter :: vdwbj_s8_pwpb95=0.2904_dp, vdwbj_s8_revpbe0=1.7588_dp 163 real(dp),parameter :: vdwbj_s8_revpbe38=1.4760_dp, vdwbj_s8_revpbe=2.3550_dp 164 real(dp),parameter :: vdwbj_s8_rpw86pbe=1.3845_dp, vdwbj_s8_tpss0=1.2576_dp 165 real(dp),parameter :: vdwbj_s8_tpss=1.9435_dp 166 ! a1 parameters (BJ only) 167 real(dp),parameter :: vdwbj_a1_b1b95=0.2092_dp, vdwbj_a1_b2gpplyp=0.0000_dp 168 real(dp),parameter :: vdwbj_a1_b3pw91=0.4312_dp, vdwbj_a1_bhlyp=0.2793_dp 169 real(dp),parameter :: vdwbj_a1_bmk=0.1940_dp, vdwbj_a1_bop=0.4870_dp 170 real(dp),parameter :: vdwbj_a1_bpbe=0.4567_dp, vdwbj_a1_camb3lyp=0.3708_dp 171 real(dp),parameter :: vdwbj_a1_lcwpbe=0.3919_dp, vdwbj_a1_mpw1b95=0.1955_dp 172 real(dp),parameter :: vdwbj_a1_mpwb1k=0.1474_dp, vdwbj_a1_mpwlyp=0.4831_dp 173 real(dp),parameter :: vdwbj_a1_olyp=0.5299_dp, vdwbj_a1_opbe=0.5512_dp 174 real(dp),parameter :: vdwbj_a1_otpss=0.4634_dp, vdwbj_a1_pbe38=0.3995_dp 175 real(dp),parameter :: vdwbj_a1_pbesol=0.4466_dp, vdwbj_a1_ptpss=0.000_dp 176 real(dp),parameter :: vdwbj_a1_pwb6k=0.1805_dp, vdwbj_a1_revssb=0.4720_dp 177 real(dp),parameter :: vdwbj_a1_ssb=-0.0952_dp, vdwbj_a1_tpssh=0.4529_dp 178 real(dp),parameter :: vdwbj_a1_hcth120=0.3563_dp, vdwbj_a1_b2plyp=0.3065_dp 179 real(dp),parameter :: vdwbj_a1_b3lyp=0.3981_dp, vdwbj_a1_b97d=0.5545_dp 180 real(dp),parameter :: vdwbj_a1_blyp=0.4298_dp, vdwbj_a1_bp86=0.3946_dp 181 real(dp),parameter :: vdwbj_a1_dsdblyp=0.000_dp, vdwbj_a1_pbe0=0.4145_dp 182 real(dp),parameter :: vdwbj_a1_pbe=0.4289_dp, vdwbj_a1_pw6b95=0.2076_dp 183 real(dp),parameter :: vdwbj_a1_pwpb95=0.0000_dp, vdwbj_a1_revpbe0=0.4679_dp 184 real(dp),parameter :: vdwbj_a1_revpbe38=0.4309_dp, vdwbj_a1_revpbe=0.5238_dp 185 real(dp),parameter :: vdwbj_a1_rpw86pbe=0.4613_dp, vdwbj_a1_tpss0=0.3768_dp 186 real(dp),parameter :: vdwbj_a1_tpss=0.4535_dp 187 ! a2 parameters (BJ only) 188 real(dp),parameter :: vdwbj_a2_b1b95=5.5545_dp, vdwbj_a2_b2gpplyp=6.3332_dp 189 real(dp),parameter :: vdwbj_a2_b3pw91=4.4693_dp, vdwbj_a2_bhlyp=4.9615_dp 190 real(dp),parameter :: vdwbj_a2_bmk=5.9197_dp, vdwbj_a2_bop=3.5043_dp 191 real(dp),parameter :: vdwbj_a2_bpbe=4.3908_dp, vdwbj_a2_camb3lyp=5.4743_dp 192 real(dp),parameter :: vdwbj_a2_lcwpbe=5.0897_dp, vdwbj_a2_mpw1b95=6.4177_dp 193 real(dp),parameter :: vdwbj_a2_mpwb1k=6.6223_dp, vdwbj_a2_mpwlyp=4.5323_dp 194 real(dp),parameter :: vdwbj_a2_olyp=2.8065_dp, vdwbj_a2_opbe=2.9444_dp 195 real(dp),parameter :: vdwbj_a2_otpss=4.3153_dp, vdwbj_a2_pbe38=5.1405_dp 196 real(dp),parameter :: vdwbj_a2_pbesol=6.1742_dp, vdwbj_a2_ptpss=6.5745_dp 197 real(dp),parameter :: vdwbj_a2_pwb6k=7.7627_dp, vdwbj_a2_revssb=4.0986_dp 198 real(dp),parameter :: vdwbj_a2_ssb=5.2170_dp, vdwbj_a2_tpssh=4.6550_dp 199 real(dp),parameter :: vdwbj_a2_hcth120=4.3359_dp, vdwbj_a2_b2plyp=5.0570_dp 200 real(dp),parameter :: vdwbj_a2_b3lyp=4.4211_dp, vdwbj_a2_b97d=3.2297_dp 201 real(dp),parameter :: vdwbj_a2_blyp=4.2359_dp, vdwbj_a2_bp86=4.8516_dp 202 real(dp),parameter :: vdwbj_a2_dsdblyp=6.0519_dp, vdwbj_a2_pbe0=4.8593_dp 203 real(dp),parameter :: vdwbj_a2_pbe=4.4407_dp, vdwbj_a2_pw6b95=6.3750_dp 204 real(dp),parameter :: vdwbj_a2_pwpb95=7.3141_dp, vdwbj_a2_revpbe0=3.7619_dp 205 real(dp),parameter :: vdwbj_a2_revpbe38=3.9446_dp, vdwbj_a2_revpbe=3.5016_dp 206 real(dp),parameter :: vdwbj_a2_rpw86pbe=4.5062_dp, vdwbj_a2_tpss0=4.5865_dp 207 real(dp),parameter :: vdwbj_a2_tpss=4.4752_dp 208 ! s6 parameters (zero damping) 209 real(dp),parameter :: vdw_s6_b2gpplyp=0.56_dp, vdw_s6_b2plyp=0.64_dp 210 real(dp),parameter :: vdw_s6_dsdblyp=0.50_dp, vdw_s6_ptpss=0.75_dp 211 real(dp),parameter :: vdw_s6_pwpb95=0.82_dp 212 ! s8 parameters (zero damping) 213 real(dp),parameter :: vdw_s8_b1b95=1.868_dp, vdw_s8_b2gpplyp=0.760_dp 214 real(dp),parameter :: vdw_s8_b3lyp=1.703_dp, vdw_s8_b97d=0.909_dp 215 real(dp),parameter :: vdw_s8_bhlyp=1.442_dp, vdw_s8_blyp=1.682_dp 216 real(dp),parameter :: vdw_s8_bp86=1.683_dp, vdw_s8_bpbe=2.033_dp 217 real(dp),parameter :: vdw_s8_mpwlyp=1.098_dp, vdw_s8_pbe=0.722_dp 218 real(dp),parameter :: vdw_s8_pbe0=0.928_dp, vdw_s8_pw6b95=0.862_dp 219 real(dp),parameter :: vdw_s8_pwb6k=0.550_dp, vdw_s8_revpbe=1.010_dp 220 real(dp),parameter :: vdw_s8_tpss=1.105_dp, vdw_s8_tpss0=1.242_dp 221 real(dp),parameter :: vdw_s8_tpssh=1.219_dp, vdw_s8_bop=1.975_dp 222 real(dp),parameter :: vdw_s8_mpw1b95=1.118_dp, vdw_s8_mpwb1k=1.061_dp 223 real(dp),parameter :: vdw_s8_olyp=1.764_dp, vdw_s8_opbe=2.055_dp 224 real(dp),parameter :: vdw_s8_otpss=1.494_dp, vdw_s8_pbe38=0.998_dp 225 real(dp),parameter :: vdw_s8_pbesol=0.612_dp, vdw_s8_revssb=0.560_dp 226 real(dp),parameter :: vdw_s8_ssb=0.663_dp, vdw_s8_b3pw91=1.775_dp 227 real(dp),parameter :: vdw_s8_bmk=2.168_dp, vdw_s8_camb3lyp=1.217_dp 228 real(dp),parameter :: vdw_s8_lcwpbe=1.279_dp, vdw_s8_m052x=0.00_dp 229 real(dp),parameter :: vdw_s8_m05=0.595_dp, vdw_s8_m062x=0.00_dp 230 real(dp),parameter :: vdw_s8_m06hf=0.00_dp, vdw_s8_m06l=0.00_dp 231 real(dp),parameter :: vdw_s8_m06=0.00_dp, vdw_s8_hcth120=1.206_dp 232 real(dp),parameter :: vdw_s8_b2plyp=1.022_dp, vdw_s8_dsdblyp=0.705_dp 233 real(dp),parameter :: vdw_s8_ptpss=0.879_dp, vdw_s8_pwpb95=0.705_dp 234 real(dp),parameter :: vdw_s8_revpbe0=0.792_dp, vdw_s8_revpbe38=0.862_dp 235 real(dp),parameter :: vdw_s8_rpw86pbe=0.901_dp 236 ! sr6 parameters (zero damping) 237 real(dp),parameter :: vdw_sr6_b1b95=1.613_dp, vdw_sr6_b2gpplyp=1.586_dp 238 real(dp),parameter :: vdw_sr6_b3lyp=1.261_dp, vdw_sr6_b97d=0.892_dp 239 real(dp),parameter :: vdw_sr6_bhlyp=1.370_dp, vdw_sr6_blyp=1.094_dp 240 real(dp),parameter :: vdw_sr6_bp86=1.139_dp, vdw_sr6_bpbe=1.087_dp 241 real(dp),parameter :: vdw_sr6_mpwlyp=1.239_dp, vdw_sr6_pbe=1.217_dp 242 real(dp),parameter :: vdw_sr6_pbe0=1.287_dp, vdw_sr6_pw6b95=1.532_dp 243 real(dp),parameter :: vdw_sr6_pwb6k=1.660_dp, vdw_sr6_revpbe=0.923_dp 244 real(dp),parameter :: vdw_sr6_tpss=1.166_dp, vdw_sr6_tpss0=1.252_dp 245 real(dp),parameter :: vdw_sr6_tpssh=1.223_dp, vdw_sr6_bop=0.929_dp 246 real(dp),parameter :: vdw_sr6_mpw1b95=1.605_dp, vdw_sr6_mpwb1k=1.671_dp 247 real(dp),parameter :: vdw_sr6_olyp=0.806_dp, vdw_sr6_opbe=0.837_dp 248 real(dp),parameter :: vdw_sr6_otpss=1.128_dp, vdw_sr6_pbe38=1.333_dp 249 real(dp),parameter :: vdw_sr6_pbesol=1.345_dp, vdw_sr6_revssb=1.221_dp 250 real(dp),parameter :: vdw_sr6_ssb=1.215_dp, vdw_sr6_b3pw91=1.176_dp 251 real(dp),parameter :: vdw_sr6_bmk=1.931_dp, vdw_sr6_camb3lyp=1.378_dp 252 real(dp),parameter :: vdw_sr6_lcwpbe=1.355_dp, vdw_sr6_m052x=1.417_dp 253 real(dp),parameter :: vdw_sr6_m05=1.373_dp, vdw_sr6_m062x=1.619_dp 254 real(dp),parameter :: vdw_sr6_m06hf=1.446_dp, vdw_sr6_m06l=1.581_dp 255 real(dp),parameter :: vdw_sr6_m06=1.325_dp, vdw_sr6_hcth120=1.221_dp 256 real(dp),parameter :: vdw_sr6_b2plyp=1.427_dp, vdw_sr6_dsdblyp=1.569_dp 257 real(dp),parameter :: vdw_sr6_ptpss=1.541_dp, vdw_sr6_pwpb95=1.557_dp 258 real(dp),parameter :: vdw_sr6_revpbe0=0.949_dp, vdw_sr6_revpbe38=1.021_dp 259 real(dp),parameter :: vdw_sr6_rpw86pbe=1.224_dp 260 ! sr8 parameters (zero damping) = 1.000_dp 261 262 real(dp),parameter :: vdw_sr9=3.0/4.0 263 real(dp),parameter :: vdw_tol_default=tol10 264 real(dp) :: ang,arg,cn_dmp,cosa,cosb,cosc,c6,c8 265 real(dp) :: dcosa_r3drij,dcosa_r3drjk,dcosa_r3drki 266 real(dp) :: dcosb_r3drij,dcosb_r3drjk,dcosb_r3drki 267 real(dp) :: dcosc_r3drij,dcosc_r3drjk,dcosc_r3drki 268 real(dp) :: dcn_dmp,dexp_cn,dfdmp,dfdmp_drij 269 real(dp) :: dfdmp_drjk,dfdmp_drki,dlri,dlrj,dmp,dmp6,dmp8,dmp9,dr,d2lri,d2lrj,d2lrirj 270 real(dp) :: dsysref,dsysref_a, dsysref_b 271 real(dp) :: d1_r3drij,d1_r3drjk,d1_r3drki,d2cn_dmp,d2cn_exp,d2frac_cn 272 real(dp) :: d_drij,d_drjk,d_drki 273 real(dp) :: exp_cn,e_no_c6,e_no_c8,e_3bt,fdmp6,fdmp8,fdmp9,frac_cn 274 real(dp) :: grad,grad_no_c,grad6,grad6_no_c6,grad8,grad8_no_c8,gr6,gr8 275 real(dp) :: hess,hessij, hess6, hess8,im_arg,l,ltot 276 real(dp) :: max_vdw_c6,min_dsys,re_arg,rcovij,rcut,rcutcn,rcut2,rcut9 277 real(dp) :: rsq,rsqij,rsqjk,rsqki,rmean,rr,rrij,rrjk,rrki,rijk,r0,r6,r8 278 real(dp) :: sfact6,sfact8,sfact9,sum_dlri,sum_dlrj,sum_dlc6ri,sum_dlc6rj 279 real(dp) :: sum_d2lri,sum_d2lrj,sum_d2lrirj,sum_d2lc6ri,sum_d2lc6rj,sum_d2lc6rirj 280 real(dp) :: temp,temp2 281 real(dp) :: ucvol,vdw_s6,vdw_s8,vdw_sr6,vdw_sr8,vdw_a1,vdw_a2,vdw_q 282 character(len=500) :: msg 283 type(atomdata_t) :: atom1,atom2 284 285 !arrays 286 287 ! Covalence radius of the different species for CN (coordination number) 288 real(dp),parameter:: rcov(vdw_nspecies)=& 289 & (/0.80628308, 1.15903197, 3.02356173, 2.36845659, 1.94011865, & 290 & 1.88972601, 1.78894056, 1.58736983, 1.61256616, 1.68815527, & 291 & 3.52748848, 3.14954334, 2.84718717, 2.62041997, 2.77159820, & 292 & 2.57002732, 2.49443835, 2.41884923, 4.43455700, 3.88023730, & 293 & 3.35111422, 3.07395437, 3.04875805, 2.77159820, 2.69600923, & 294 & 2.62041997, 2.51963467, 2.49443835, 2.54483100, 2.74640188, & 295 & 2.82199085, 2.74640188, 2.89757982, 2.77159820, 2.87238349, & 296 & 2.94797246, 4.76210950, 4.20778980, 3.70386304, 3.50229216, & 297 & 3.32591790, 3.12434702, 2.89757982, 2.84718717, 2.84718717, & 298 & 2.72120556, 2.89757982, 3.09915070, 3.22513231, 3.17473967, & 299 & 3.17473967, 3.09915070, 3.32591790, 3.30072128, 5.26603625, & 300 & 4.43455700, 4.08180818, 3.70386304, 3.98102289, 3.95582657, & 301 & 3.93062995, 3.90543362, 3.80464833, 3.82984466, 3.80464833, & 302 & 3.77945201, 3.75425569, 3.75425569, 3.72905937, 3.85504098, & 303 & 3.67866672, 3.45189952, 3.30072128, 3.09915070, 2.97316878, & 304 & 2.92277614, 2.79679452, 2.82199085, 2.84718717, 3.32591790, & 305 & 3.27552496, 3.27552496, 3.42670319, 3.30072128, 3.47709584, & 306 & 3.57788113, 5.06446567, 4.56053862, 4.20778980, 3.98102289, & 307 & 3.82984466, 3.85504098, 3.88023730, 3.90543362 /) 308 309 ! q = arrays of vdw_species elements containing the link between C6ij and C8ij: 310 ! C8ij = 3sqrt(qi)sqrt(qj)C6ij 311 real(dp),parameter :: vdw_q_dftd3(vdw_nspecies)= & 312 & (/2.00734898, 1.56637132, 5.01986934, 3.85379032, 3.64446594, & 313 & 3.10492822, 2.71175247, 2.59361680, 2.38825250, 2.21522516, & 314 & 6.58585536, 5.46295967, 5.65216669, 4.88284902, 4.29727576, & 315 & 4.04108902, 3.72932356, 3.44677275, 7.97762753, 7.07623947, & 316 & 6.60844053, 6.28791364, 6.07728703, 5.54643096, 5.80491167, & 317 & 5.58415602, 5.41374528, 5.28497229, 5.22592821, 5.09817141, & 318 & 6.12149689, 5.54083734, 5.06696878, 4.87005108, 4.59089647, & 319 & 4.31176304, 9.55461698, 8.67396077, 7.97210197, 7.43439917, & 320 & 6.58711862, 6.19536215, 6.01517290, 5.81623410, 5.65710424, & 321 & 5.52640661, 5.44263305, 5.58285373, 7.02081898, 6.46815523, & 322 & 5.98089120, 5.81686657, 5.53321815, 5.25477007, 11.02204549, & 323 & 10.15679528, 9.35167836, 9.06926079, 8.97241155, 8.90092807, & 324 & 8.85984840, 8.81736827, 8.79317710, 7.89969626, 8.80588454, & 325 & 8.42439218, 8.54289262, 8.47583370, 8.45090888, 8.47339339, & 326 & 7.83525634, 8.20702843, 7.70559063, 7.32755997, 7.03887381, & 327 & 6.68978720, 6.05450052, 5.88752022, 5.70661499, 5.78450695, & 328 & 7.79780729, 7.26443867, 6.78151984, 6.67883169, 6.39024318, & 329 & 6.09527958, 11.79156076, 11.10997644, 9.51377795, 8.67197068, & 330 & 8.77140725, 8.65402716, 8.53923501, 8.85024712 /) 331 332 integer :: is(3), nshell_3bt(3) 333 integer,allocatable :: ivdw(:) 334 integer :: jmin(3), jmax(3), js(3) 335 integer,parameter :: voigt1(6)=(/1,2,3,2,1,1/),voigt2(6)=(/1,2,3,3,3,2/) 336 real(dp),allocatable :: cn(:),cfgrad_no_c(:,:,:,:) 337 real(dp),allocatable :: dcn(:,:,:),dcn_cart(:,:,:),dc6ri(:,:),dc6rj(:,:),dc9ijri(:,:),dc9ijrj(:,:) 338 !real(dp),allocatable:: d2cn(:,:,:,:,:,:) 339 real(dp),allocatable:: d2cn_iii(:,:,:,:) 340 real(dp),allocatable:: d2cn_jji(:,:,:,:,:) 341 real(dp),allocatable:: d2cn_iji(:,:,:,:,:) 342 real(dp),allocatable:: d2cn_jii(:,:,:,:,:) 343 real(dp),allocatable:: d2cn_tmp(:) 344 real(dp),allocatable :: d2c6ri(:,:),d2c6rj(:,:),d2c6rirj(:,:) 345 real(dp),allocatable:: elt_cn(:,:,:),e3bt_ij(:,:),e3bt_jk(:,:),e3bt_ki(:,:),e_no_c(:,:) 346 real(dp),allocatable:: e_alpha1(:),e_alpha2(:),e_alpha3(:),e_alpha4(:) 347 real(dp) :: gred(3),gredij(3),gredjk(3),gredki(3) 348 real(dp),allocatable:: fe_no_c(:,:,:),cfdcn(:,:,:,:),fdcn(:,:,:,:),fgrad_no_c(:,:,:,:),gred_vdw_3bt(:,:) 349 real(dp) :: gmet(3,3),gprimd(3,3) 350 real(dp),allocatable:: grad_no_cij(:,:,:) 351 real(dp) :: mcart(3,3) 352 real(dp) :: r(3),rcart(3),rcart2(3,3),rcartij(3),rcartjk(3),rcartki(3) 353 real(dp):: rij(3), rjk(3), rki(3),rmet(3,3),rred(3) 354 real(dp),allocatable:: r0ijk(:,:,:) 355 real(dp),allocatable :: str_alpha1(:,:),str_alpha2(:,:),str_dcn(:,:),str_no_c(:,:,:) 356 real(dp) :: str_3bt(6) 357 real(dp) :: temp_comp(2),temp_comp2(2) 358 real(dp),allocatable:: temp_prod(:,:) 359 real(dp) :: vec(6),vecij(6), vecjk(6),vecki(6) 360 real(dp),allocatable :: vdw_cnrefi(:,:,:,:),vdw_cnrefj(:,:,:,:) 361 real(dp),allocatable :: vdw_c6(:,:),vdw_c6ref(:,:,:,:),vdw_c8(:,:),vdw_c9(:,:,:),vdw_r0(:,:) 362 real(dp):: vdw_dftd3_r0(4465) 363 real(dp):: vdw_dftd3_c6(32385) 364 integer:: index_c6(254) 365 real(dp):: vdw_dftd3_cni(27884) 366 integer:: index_cni(27884) 367 real(dp):: vdw_dftd3_cnj(13171) 368 integer:: index_cnj(13171) 369 real(dp),allocatable :: xred01(:,:) 370 371 ! ************************************************************************* 372 373 DBG_ENTER("COLL") 374 375 write(msg,'(1a)')& 376 & '====> STARTING DFT-D3 computation' 377 call wrtout(std_out,msg,'COLL') 378 379 call vdw_dftd3_data(vdw_dftd3_r0,vdw_dftd3_c6,index_c6,vdw_dftd3_cni,index_cni,& 380 & vdw_dftd3_cnj,index_cnj) 381 ! Determine the properties which have to be studied 382 bol_3bt = (vdw_tol_3bt>0) 383 need_forces = present(gred_vdw_dftd3) 384 need_stress= present(str_vdw_dftd3) 385 need_dynmat= present(dyn_vdw_dftd3) 386 need_elast= present(elt_vdw_dftd3) 387 need_grad=(need_forces.or.need_stress.or.need_dynmat.or.need_elast) 388 need_hess=(need_dynmat.or.need_elast) 389 if (need_dynmat) then 390 if (.not.present(qphon)) then 391 msg='Dynamical matrix required without a q-vector' 392 ABI_BUG(msg) 393 end if 394 dyn_vdw_dftd3=zero 395 end if 396 e_vdw_dftd3 = zero 397 if (need_forces) gred_vdw_dftd3=zero 398 if (need_stress) str_vdw_dftd3=zero 399 if (need_elast) elt_vdw_dftd3=zero 400 401 !Identify type(s) of atoms 402 ABI_MALLOC(ivdw,(ntypat)) 403 do itypat=1,ntypat 404 call atomdata_from_znucl(atom1,znucl(itypat)) 405 if (znucl(itypat).gt.94.0_dp) then 406 write(msg,'(3a,es14.2)') & 407 & 'Van der Waals DFT-D3 correction not available for atom type: ',znucl(itypat),' !' 408 ABI_ERROR(msg) 409 else 410 ivdw(itypat) = znucl(itypat) 411 end if 412 end do 413 414 ! Determination of coefficients that depend of the 415 ! exchange-correlation functional 416 vdw_s6 =one; vdw_s8 =one 417 vdw_sr6=one; vdw_sr8=one 418 vdw_a1 =one; vdw_a2 =one 419 ! Case one : DFT-D3 420 if (vdw_xc == 6) then 421 select case (ixc) 422 case(11, -130101, -101130) 423 vdw_sr6=vdw_sr6_pbe; vdw_s8=vdw_s8_pbe 424 case(14, -130102, -102130) 425 vdw_sr6=vdw_sr6_revpbe; vdw_s8=vdw_s8_revpbe 426 case(18, -131106, -106131) 427 vdw_sr6=vdw_sr6_blyp; vdw_s8=vdw_s8_blyp 428 case(19, -132106, -106132) 429 vdw_sr6=vdw_sr6_bp86; vdw_s8=vdw_s8_bp86 430 case(41, -406) 431 vdw_sr6=vdw_sr6_pbe0; vdw_s8=vdw_s8_pbe0 432 case(-440) 433 vdw_sr6=vdw_sr6_b1b95; vdw_s8=vdw_s8_b1b95 434 case(-402) 435 vdw_sr6=vdw_sr6_b3lyp; vdw_s8=vdw_s8_b3lyp 436 case(-170) 437 vdw_sr6=vdw_sr6_b97d; vdw_s8=vdw_s8_b97d 438 case(-435) 439 vdw_sr6=vdw_sr6_bhlyp; vdw_s8=vdw_s8_bhlyp 440 case(-130106, -106130) 441 vdw_sr6=vdw_sr6_bpbe; vdw_s8=vdw_s8_bpbe 442 case(-174) 443 vdw_sr6=vdw_sr6_mpwlyp; vdw_s8=vdw_s8_mpwlyp 444 case(-451) 445 vdw_sr6=vdw_sr6_pw6b95; vdw_s8=vdw_s8_pw6b95 446 case(-452) 447 vdw_sr6=vdw_sr6_pwb6k; vdw_s8=vdw_s8_pwb6k 448 case(-231202, -202231) 449 vdw_sr6=vdw_sr6_tpss; vdw_s8=vdw_s8_tpss 450 case(-396) 451 vdw_sr6=vdw_sr6_tpss0; vdw_s8=vdw_s8_tpss0 452 case(-457) 453 vdw_sr6=vdw_sr6_tpssh; vdw_s8=vdw_s8_tpssh 454 case(-636) 455 vdw_sr6=vdw_sr6_bop; vdw_s8=vdw_s8_bop 456 case(-445) 457 vdw_sr6=vdw_sr6_mpw1b95; vdw_s8=vdw_s8_mpw1b95 458 case(-446) 459 vdw_sr6=vdw_sr6_mpwb1k; vdw_s8=vdw_s8_mpwb1k 460 case(-67) 461 vdw_sr6=vdw_sr6_olyp; vdw_s8=vdw_s8_olyp 462 case(-65) 463 vdw_sr6=vdw_sr6_opbe; vdw_s8=vdw_s8_opbe 464 case(-64) 465 vdw_sr6=vdw_sr6_otpss; vdw_s8=vdw_s8_otpss 466 case(-393) 467 vdw_sr6=vdw_sr6_pbe38; vdw_s8=vdw_s8_pbe38 468 case(-133116, -116133) 469 vdw_sr6=vdw_sr6_pbesol; vdw_s8=vdw_s8_pbesol 470 case(-312089, -89312) 471 vdw_sr6=vdw_sr6_revssb; vdw_s8=vdw_s8_revssb 472 case(-91089, -89091) 473 vdw_sr6=vdw_sr6_ssb; vdw_s8=vdw_s8_ssb 474 case(-401) 475 vdw_sr6=vdw_sr6_b3pw91; vdw_s8=vdw_s8_b3pw91 476 case(-280279, -279280) 477 vdw_sr6=vdw_sr6_bmk; vdw_s8=vdw_s8_bmk 478 case(-433) 479 vdw_sr6=vdw_sr6_camb3lyp; vdw_s8=vdw_s8_camb3lyp 480 case(-478) 481 vdw_sr6=vdw_sr6_lcwpbe; vdw_s8=vdw_s8_lcwpbe 482 case(-439238, -238439) 483 vdw_sr6=vdw_sr6_m052x; vdw_s8=vdw_s8_m052x 484 case(-438237, -237438) 485 vdw_sr6=vdw_sr6_m05; vdw_s8=vdw_s8_m05 486 case(-450236, -236450) 487 vdw_sr6=vdw_sr6_m062x; vdw_s8=vdw_s8_m062x 488 case(-444234, -234444) 489 vdw_sr6=vdw_sr6_m06hf; vdw_s8=vdw_s8_m06hf 490 case(-203233, -233203) 491 vdw_sr6=vdw_sr6_m06l; vdw_s8=vdw_s8_m06l 492 case(-449235, -235449) 493 vdw_sr6=vdw_sr6_m06; vdw_s8=vdw_s8_m06 494 case(-162) 495 vdw_sr6=vdw_sr6_hcth120; vdw_s8=vdw_s8_hcth120 496 case(-456) 497 vdw_sr6=vdw_sr6_revpbe0; vdw_s8=vdw_s8_revpbe0 498 case(-30108, -108030) 499 vdw_sr6=vdw_sr6_rpw86pbe; vdw_s8=vdw_s8_rpw86pbe 500 case default 501 write(msg,'(a,i8,a)')' Van der Waals DFT-D3 correction not compatible with ixc=',ixc,' !' 502 ABI_ERROR(msg) 503 end select 504 ! Case DFT-D3(BJ) 505 elseif (vdw_xc == 7) then 506 select case (ixc) 507 case(11, -130101, -101130) 508 vdw_s8=vdwbj_s8_pbe; vdw_a1=vdwbj_a1_pbe; vdw_a2=vdwbj_a2_pbe 509 case(14, -130102, -102130) 510 vdw_s8=vdwbj_s8_revpbe; vdw_a1=vdwbj_a1_revpbe; vdw_a2=vdwbj_a2_revpbe 511 case(18, -131106, -106131) 512 vdw_s8=vdwbj_s8_blyp; vdw_a1=vdwbj_a1_blyp; vdw_a2=vdwbj_a2_blyp 513 case(19, -132106, -106132) 514 vdw_s8=vdwbj_s8_bp86; vdw_a1=vdwbj_a1_bp86; vdw_a2=vdwbj_a2_bp86 515 case(41, -406) 516 vdw_s8=vdwbj_s8_pbe0; vdw_a1=vdwbj_a1_pbe0; vdw_a2=vdwbj_a2_pbe0 517 case(-440) 518 vdw_s8=vdwbj_s8_b1b95; vdw_a1=vdwbj_a1_b1b95; vdw_a2=vdwbj_a2_b1b95 519 case(-401) 520 vdw_s8=vdwbj_s8_b3pw91; vdw_a1=vdwbj_a1_b3pw91; vdw_a2=vdwbj_a2_b3pw91 521 case(-435) 522 vdw_s8=vdwbj_s8_bhlyp; vdw_a1=vdwbj_a1_bhlyp; vdw_a2=vdwbj_a2_bhlyp 523 case(-280279, -279280) 524 vdw_s8=vdwbj_s8_bmk; vdw_a1=vdwbj_a1_bmk; vdw_a2=vdwbj_a2_bmk 525 case(-636) 526 vdw_s8=vdwbj_s8_bop; vdw_a1=vdwbj_a1_bop; vdw_a2=vdwbj_a2_bop 527 case(-130106, -106130) 528 vdw_s8=vdwbj_s8_bpbe; vdw_a1=vdwbj_a1_bpbe; vdw_a2=vdwbj_a2_bpbe 529 case(-433) 530 vdw_s8=vdwbj_s8_camb3lyp; vdw_a1=vdwbj_a1_camb3lyp; vdw_a2=vdwbj_a2_camb3lyp 531 case(-478) 532 vdw_s8=vdwbj_s8_lcwpbe; vdw_a1=vdwbj_a1_lcwpbe; vdw_a2=vdwbj_a2_lcwpbe 533 case(-445) 534 vdw_s8=vdwbj_s8_mpw1b95; vdw_a1=vdwbj_a1_mpw1b95; vdw_a2=vdwbj_a2_mpw1b95 535 case(-446) 536 vdw_s8=vdwbj_s8_mpwb1k; vdw_a1=vdwbj_a1_mpwb1k; vdw_a2=vdwbj_a2_mpwb1k 537 case(-174) 538 vdw_s8=vdwbj_s8_mpwlyp; vdw_a1=vdwbj_a1_mpwlyp; vdw_a2=vdwbj_a2_mpwlyp 539 case(-67) 540 vdw_s8=vdwbj_s8_olyp; vdw_a1=vdwbj_a1_olyp; vdw_a2=vdwbj_a2_olyp 541 case(-65) 542 vdw_s8=vdwbj_s8_opbe; vdw_a1=vdwbj_a1_opbe; vdw_a2=vdwbj_a2_opbe 543 case(-64) 544 vdw_s8=vdwbj_s8_otpss; vdw_a1=vdwbj_a1_otpss; vdw_a2=vdwbj_a2_otpss 545 case(-393) 546 vdw_s8=vdwbj_s8_pbe38; vdw_a1=vdwbj_a1_pbe38; vdw_a2=vdwbj_a2_pbe38 547 case(-133116, -116133) 548 vdw_s8=vdwbj_s8_pbesol; vdw_a1=vdwbj_a1_pbesol; vdw_a2=vdwbj_a2_pbesol 549 case(-452) 550 vdw_s8=vdwbj_s8_pwb6k; vdw_a1=vdwbj_a1_pwb6k; vdw_a2=vdwbj_a2_pwb6k 551 case(-312089, -89312) 552 vdw_s8=vdwbj_s8_revssb; vdw_a1=vdwbj_a1_revssb; vdw_a2=vdwbj_a2_revssb 553 case(-91089, -89091) 554 vdw_s8=vdwbj_s8_ssb; vdw_a1=vdwbj_a1_ssb; vdw_a2=vdwbj_a2_ssb 555 case(-457) 556 vdw_s8=vdwbj_s8_tpssh; vdw_a1=vdwbj_a1_tpssh; vdw_a2=vdwbj_a2_tpssh 557 case(-162) 558 vdw_s8=vdwbj_s8_hcth120; vdw_a1=vdwbj_a1_hcth120; vdw_a2=vdwbj_a2_hcth120 559 case(-402) 560 vdw_s8=vdwbj_s8_b3lyp; vdw_a1=vdwbj_a1_b3lyp; vdw_a2=vdwbj_a2_b3lyp 561 case(-170) 562 vdw_s8=vdwbj_s8_b97d; vdw_a1=vdwbj_a1_b97d; vdw_a2=vdwbj_a2_b97d 563 case(-451) 564 vdw_s8=vdwbj_s8_pw6b95; vdw_a1=vdwbj_a1_pw6b95; vdw_a2=vdwbj_a2_pw6b95 565 case(-30108, -108030) 566 vdw_s8=vdwbj_s8_rpw86pbe; vdw_a1=vdwbj_a1_rpw86pbe; vdw_a2=vdwbj_a2_rpw86pbe 567 case(-396) 568 vdw_s8=vdwbj_s8_tpss0; vdw_a1=vdwbj_a1_tpss0; vdw_a2=vdwbj_a2_tpss0 569 case(-231202, -202231) 570 vdw_s8=vdwbj_s8_tpss; vdw_a1=vdwbj_a1_tpss; vdw_a2=vdwbj_a2_tpss 571 case default 572 write(msg,'(a,i8,a)')' Van der Waals DFT-D3(BJ) correction not compatible with ixc=',ixc,' !' 573 ABI_ERROR(msg) 574 end select 575 end if 576 577 ! -------------------------------------------------------------- 578 ! Retrieve the data for the referenced c6, cn and r0 coefficients 579 !--------------------------------------------------------------- 580 refmax = 5 581 582 ABI_MALLOC(vdw_c6ref,(ntypat,ntypat,refmax,refmax)) 583 ABI_MALLOC(vdw_cnrefi,(ntypat,ntypat,refmax,refmax)) 584 ABI_MALLOC(vdw_cnrefj,(ntypat,ntypat,refmax,refmax)) 585 ABI_MALLOC(vdw_r0,(ntypat,ntypat)) 586 if (bol_3bt) then 587 ABI_MALLOC(r0ijk,(ntypat,ntypat,ntypat)) 588 end if 589 590 vdw_c6ref = zero 591 vdw_cnrefi = 100 ; vdw_cnrefj = 100 ; 592 593 do refi=1,refmax 594 do refj=1,refi 595 do itypat=1,ntypat 596 do jtypat=1,ntypat 597 indi = ivdw(itypat)+100*(refi-1) 598 indj = ivdw(jtypat)+100*(refj-1) 599 found = .false. 600 do ia=1,size(index_c6) 601 do ja=1,size(index_c6) 602 if (index_c6(ia)==indi.and.index_c6(ja)==indj) then 603 if (ia>=ja)then 604 nline = ia*(ia-1)/2 + ja 605 else 606 nline = ja*(ja-1)/2 + ia 607 endif 608 vdw_c6ref(itypat,jtypat,refi,refj) = vdw_dftd3_c6(nline) 609 vdw_c6ref(jtypat,itypat,refj,refi) = vdw_dftd3_c6(nline) 610 found = .false. 611 do la=1,size(index_cni) 612 if (index_cni(la)==nline) then 613 found=.true. 614 vdw_cnrefi(itypat,jtypat,refi,refj)= vdw_dftd3_cni(la) 615 vdw_cnrefj(jtypat,itypat,refj,refi)= vdw_dftd3_cni(la) 616 else 617 vdw_cnrefi(itypat,jtypat,refi,refj) = zero 618 vdw_cnrefj(jtypat,itypat,refj,refi) = zero 619 end if 620 if (found) exit 621 end do 622 found = .false. 623 do la=1,size(index_cnj) 624 if (index_cnj(la)==nline) then 625 found=.true. 626 vdw_cnrefj(itypat,jtypat,refi,refj)= vdw_dftd3_cnj(la) 627 vdw_cnrefi(jtypat,itypat,refj,refi)= vdw_dftd3_cnj(la) 628 else 629 vdw_cnrefj(itypat,jtypat,refi,refj) = zero 630 vdw_cnrefi(jtypat,itypat,refj,refi) = zero 631 end if 632 if (found) exit 633 end do 634 found = .true. 635 end if 636 if (found) exit 637 end do 638 if (found) exit 639 end do 640 if (refi.eq.1.and.refj.eq.1) then 641 nline = ia*(ia-1)/2 + ja 642 vdw_r0(itypat,jtypat)=vdw_dftd3_r0(nline)/Bohr_Ang 643 if (bol_3bt) then 644 do ktypat=1,ntypat 645 r0ijk(itypat,jtypat,ktypat)=one/(vdw_r0(itypat,jtypat)*vdw_r0(jtypat,ktypat)*vdw_r0(ktypat,itypat))**third 646 end do ! ka atom 647 end if ! Only if 3bt required 648 end if ! Only for the first set of references 649 end do ! Loop on references j 650 end do ! Loop on references i 651 end do ! Loop on atom j 652 end do ! Loop on atom i 653 654 !if (vdw_d3_cov==1) then 655 ! vdw_cnrefi(:,:,:,refmax) =vdw_cnrefi(:,:,:,refmax-1) 656 ! vdw_cnrefi(:,:,refmax,:) = 14.0_dp 657 ! vdw_cnrefj(:,:,refmax,:) =vdw_cnrefj(:,:,refmax-1,:) 658 ! vdw_cnrefj(:,:,:,refmax) = 14.0_dp 659 !end if 660 !Retrieve cell geometry data 661 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 662 663 !Map reduced coordinates into [0,1[ 664 ABI_MALLOC(xred01,(3,natom)) 665 do ia=1,natom 666 xred01(:,ia)=xred(:,ia)-aint(xred(:,ia)) ! Map into ]-1,1[ 667 do alpha=1,3 668 if (abs(xred01(alpha,ia)).ge.tol8) xred01(alpha,ia) = xred01(alpha,ia)+half-sign(half,xred(alpha,ia)) 669 end do 670 end do 671 672 ! ------------------------------------------------------------------- 673 ! Computation of the coordination number (CN) for the different atoms 674 ! ------------------------------------------------------------------- 675 676 write(msg,'(3a)')& 677 & ' Begin the computation of the Coordination Numbers (CN)',ch10,& 678 & ' required for DFT-D3 energy corrections...' 679 call wrtout(std_out,msg,'COLL') 680 681 ! Allocation of the CN coefficients and derivatives 682 ABI_MALLOC(cn,(natom)) 683 ABI_MALLOC(dcn,(3,natom,natom)) 684 ABI_MALLOC(dcn_cart,(3,natom,natom)) 685 ABI_MALLOC(str_dcn,(6,natom)) 686 ! ABI_MALLOC_OR_DIE(d2cn, (2,3,natom,3,natom,natom), ierr) 687 ABI_MALLOC_OR_DIE(d2cn_iii, (2,3,3,natom), ierr) 688 ABI_MALLOC_OR_DIE(d2cn_jji, (2,3,3,natom,natom), ierr) 689 ABI_MALLOC_OR_DIE(d2cn_iji, (2,3,3,natom,natom), ierr) 690 ABI_MALLOC_OR_DIE(d2cn_jii, (2,3,3,natom,natom), ierr) 691 ABI_MALLOC(d2cn_tmp, (2)) 692 ABI_MALLOC(fdcn,(2,3,natom,natom)) 693 ABI_MALLOC(cfdcn,(2,3,natom,natom)) 694 ABI_MALLOC(elt_cn,(6+3*natom,6,natom)) 695 696 ! Initializing the computed quantities to zero 697 nshell = 0 698 cn = zero 699 ! Initializing the derivative of the computed quantities to zero (if required) 700 dcn = zero ; str_dcn = zero 701 d2cn_iii = zero 702 d2cn_jji = zero 703 d2cn_iji = zero 704 d2cn_jii = zero 705 fdcn = zero; cfdcn = zero 706 elt_cn = zero ; dcn_cart = zero 707 708 re_arg = zero ; im_arg = zero 709 if (need_hess) then 710 mcart = zero 711 do alpha=1,3 712 mcart(alpha,alpha) = one 713 end do 714 end if 715 rcutcn = 200**2 ! Bohr 716 !Loop over shells of cell replicas 717 do 718 newshell=.false.;nshell=nshell+1 719 ! Loop over cell replicas in the shell 720 do is3=-nshell,nshell 721 do is2=-nshell,nshell 722 do is1=-nshell,nshell 723 if (nshell==1.or. & 724 & abs(is3)==nshell.or.abs(is2)==nshell.or.abs(is1)==nshell) then 725 is(3) = is3 ; is(2) = is2 ; is(1) = is1 726 ! Computation of phase factor for discrete Fourier transform 727 ! if phonon at qphon is required 728 if (need_dynmat) then 729 arg=two_pi*dot_product(qphon,is) 730 re_arg=cos(arg) ; im_arg=sin(arg) 731 end if 732 ! Loop over atoms ia and ja 733 do ia=1,natom 734 itypat=typat(ia) 735 do ja=1,natom 736 jtypat=typat(ja) 737 r(:)=xred01(:,ia)-xred01(:,ja)-dble(is(:)) 738 rsq=dot_product(r,matmul(rmet,r)) 739 ! atom i =/= j 740 if (rsq.ge.tol16.and.rsq<rcutcn) then 741 newshell=.true. 742 rr = sqrt(rsq) 743 rcovij = rcov(ivdw(itypat))+rcov(ivdw(jtypat)) 744 745 ! Computation of partial contribution to cn coefficients 746 exp_cn = exp(-k1*(rcovij/rr-one)) 747 frac_cn= one/(one+exp_cn) 748 ! Introduction of a damping function for the coordination 749 ! number because of the divergence with increasing 750 ! number of cells of this quantity in periodic systems 751 ! See Reckien et al., J. Comp. Chem. 33, 2023 (2012) [[cite:Reckien2012]] 752 dr = rr-k2*rcovij 753 cn_dmp = half*abi_derfc(dr) 754 cn(ia) = cn(ia)+frac_cn*cn_dmp 755 756 ! If force, stress, IFC or Elastic constants are required, 757 ! computation of the first derivative of CN 758 if (need_grad) then 759 rcart=matmul(rprimd,r) 760 dexp_cn= k1*rcovij*exp_cn/rsq 761 dcn_dmp = -one/sqrt(pi)*exp(-dr*dr) 762 grad=(-frac_cn*frac_cn*cn_dmp*dexp_cn+dcn_dmp*frac_cn)/rr 763 if (need_forces.and.ia/=ja) then 764 ! Variation of CN(ia) w.r. to displacement of atom 765 ! ja. If ia==ka then all the other atoms contribute 766 ! to the derivative. Required for the computation of 767 ! the forces applied on atom k 768 rred = matmul(transpose(rprimd),rcart) 769 dcn(:,ia,ia) = dcn(:,ia,ia)+grad*rred(:) 770 dcn(:,ia,ja) = dcn(:,ia,ja)-grad*rred(:) 771 elseif (need_stress.or.need_elast) then 772 ! The following quantity (str_dcn) is used for the computation 773 ! of the DFT-D3 contribution to stress and elastic constants 774 vec(1:3)=rcart(1:3)*rcart(1:3); vec(4)=rcart(2)*rcart(3) 775 vec(5)=rcart(1)*rcart(3); vec(6)=rcart(1)*rcart(2) 776 str_dcn(:,ia)=str_dcn(:,ia)+grad*vec(:) 777 end if 778 ! If dynamical matrix or elastic constants are required, compute 779 ! the second derivative 780 if (need_hess) then 781 d2cn_dmp = two*(rr-k2*rcovij)/sqrt(pi)*exp(-(rr-k2*rcovij)**two) 782 d2cn_exp = dexp_cn*(k1*rcovij/rsq-two/rr) 783 d2frac_cn =frac_cn**two*(two*frac_cn*dexp_cn**two-d2cn_exp) 784 hess = (d2frac_cn*cn_dmp+d2cn_dmp*frac_cn-& 785 & two*dcn_dmp*frac_cn**two*dexp_cn-grad)/rsq 786 if (need_dynmat) then 787 ! Discrete Fourier Transform of dCN/drk in cartesian 788 ! coordinates. Note that the phase factor is different 789 ! if (ka=ia or ka/=ia). 790 ! This Fourier transform is summed over cell replica 791 ! See NOTE: add reference for more information 792 fdcn(1,:,ia,ia) = fdcn(1,:,ia,ia)+grad*rcart(:) 793 fdcn(1,:,ia,ja) = fdcn(1,:,ia,ja)-grad*rcart(:)*re_arg 794 fdcn(2,:,ia,ja) = fdcn(2,:,ia,ja)-grad*rcart(:)*im_arg 795 ! Conjugate of fdcn 796 cfdcn(1,:,ia,ia) = cfdcn(1,:,ia,ia)+grad*rcart(:) 797 cfdcn(1,:,ia,ja) = cfdcn(1,:,ia,ja)-grad*rcart(:)*re_arg 798 cfdcn(2,:,ia,ja) = cfdcn(2,:,ia,ja)+grad*rcart(:)*im_arg 799 ! Computation of second derivative of CN required for the 800 ! interatomic force constants in reciprocal space 801 do alpha=1,3 802 rcart2(alpha,:) = rcart(alpha)*rcart(:) 803 end do 804 ! Computation of second derivative of CN required for the 805 ! interatomic force constants in reciprocal space 806 ! This Fourier transform is summed over cell replica 807 ! as it appears in the theory 808 do alpha=1,3 809 if (ia/=ja) then 810 ! ka = ia ; la = ia 811 d2cn_iii(1,alpha,:,ia) = d2cn_iii(1,alpha,:,ia)+& 812 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:)) 813 ! ka = ja ; la = ja 814 d2cn_jji(1,alpha,:,ja,ia) = d2cn_jji(1,alpha,:,ja,ia)+& 815 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:)) 816 if (abs(re_arg)>tol12) then 817 ! ka = ia ; la = ja 818 d2cn_iji(1,alpha,:,ja,ia) = d2cn_iji(1,alpha,:,ja,ia)-& 819 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:))*re_arg 820 ! ka = ja ; la = ia 821 d2cn_jii(1,alpha,:,ja,ia) = d2cn_jii(1,alpha,:,ja,ia)-& 822 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:))*re_arg 823 end if 824 if (abs(im_arg)>tol12) then 825 ! ka = ia ; la = ja 826 d2cn_iji(2,alpha,:,ja,ia) = d2cn_iji(2,alpha,:,ja,ia)-& 827 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:))*im_arg 828 ! ka = ja ; la = ia 829 d2cn_jii(2,alpha,:,ja,ia) = d2cn_jii(2,alpha,:,ja,ia)+& 830 & (hess*rcart2(alpha,:)+grad*mcart(alpha,:))*im_arg 831 end if 832 else 833 if (abs(re_arg-one)>tol12) then 834 d2cn_iji(1,alpha,:,ja,ia) = d2cn_iji(1,alpha,:,ja,ia)+& 835 & two*(hess*rcart2(alpha,:)+grad*mcart(alpha,:))*(one-re_arg) 836 end if 837 end if 838 end do 839 end if ! Boolean Need_dynmat 840 if (need_elast) then 841 ! Derivative of str_dcn w.r. to strain for elastic tensor 842 vec(1:3)=rcart(1:3)*rcart(1:3); vec(4)=rcart(2)*rcart(3) 843 vec(5)=rcart(1)*rcart(3); vec(6)=rcart(1)*rcart(2) 844 do alpha=1,6 845 ii = voigt1(alpha) ; jj=voigt2(alpha) 846 do beta=1,6 847 kk = voigt1(beta) ; ll=voigt2(beta) 848 elt_cn(alpha,beta,ia) = elt_cn(alpha,beta,ia)+hess*rcart(ii)*rcart(jj)*rcart(kk)*rcart(ll) 849 if (ii==kk) elt_cn(alpha,beta,ia) = elt_cn(alpha,beta,ia)+half*grad*rcart(jj)*rcart(ll) 850 if (jj==kk) elt_cn(alpha,beta,ia) = elt_cn(alpha,beta,ia)+half*grad*rcart(ii)*rcart(ll) 851 if (ii==ll) elt_cn(alpha,beta,ia) = elt_cn(alpha,beta,ia)+half*grad*rcart(jj)*rcart(kk) 852 if (jj==ll) elt_cn(alpha,beta,ia) = elt_cn(alpha,beta,ia)+half*grad*rcart(ii)*rcart(kk) 853 end do 854 end do 855 ! Derivative of str_dcn w.r. to atomic displacement 856 ! for internal strains 857 dcn_cart(:,ia,ia) = dcn_cart(:,ia,ia)+grad*rcart(:) 858 dcn_cart(:,ia,ja) = dcn_cart(:,ia,ja)-grad*rcart(:) 859 if (ia/=ja) then 860 index_ia = 6+3*(ia-1) 861 index_ja = 6+3*(ja-1) 862 do alpha=1,6 863 do beta=1,3 864 ii = voigt1(alpha) ; jj=voigt2(alpha) 865 elt_cn(index_ia+beta,alpha,ia)=elt_cn(index_ia+beta,alpha,ia)& 866 & +hess*vec(alpha)*rcart(beta) 867 elt_cn(index_ja+beta,alpha,ia)=elt_cn(index_ja+beta,alpha,ia)& 868 & -hess*vec(alpha)*rcart(beta) 869 if (ii==beta) then 870 elt_cn(index_ia+beta,alpha,ia)=elt_cn(index_ia+beta,alpha,ia)+grad*rcart(jj) 871 elt_cn(index_ja+beta,alpha,ia)=elt_cn(index_ja+beta,alpha,ia)-grad*rcart(jj) 872 end if 873 if (jj==beta) then 874 elt_cn(index_ia+beta,alpha,ia)=elt_cn(index_ia+beta,alpha,ia)+grad*rcart(ii) 875 elt_cn(index_ja+beta,alpha,ia)=elt_cn(index_ja+beta,alpha,ia)-grad*rcart(ii) 876 end if 877 end do 878 end do 879 end if ! ia/=ja 880 end if ! Need strain derivative 881 end if ! Boolean second derivative 882 end if ! Boolean first derivative 883 end if ! Tolerence 884 end do ! Loop over ia atom 885 end do ! Loop over ja atom 886 end if ! Bondary Condition 887 end do ! Loop over is1 888 end do ! Loop over is2 889 end do ! Loop over is3 890 if(.not.newshell) exit ! Check if a new shell must be considered 891 end do ! Loop over shell 892 write(msg,'(3a,f8.5,1a,i3,1a,f8.5,1a,i3,1a)')& 893 & ' ... Done.',ch10,& 894 & ' max(CN) =', maxval(cn), ' (atom ',maxloc(cn),') ; min(CN) =', minval(cn), ' (atom ', minloc(cn),')' 895 call wrtout(std_out,msg,'COLL') 896 897 !---------------------------------------------------------------- 898 ! Computation of the C6 coefficient 899 ! --------------------------------------------------------------- 900 901 write(msg,'(3a)')& 902 & ' Begin the computation of the C6(CN)',ch10,& 903 & ' required for DFT-D3 energy corrections...' 904 call wrtout(std_out,msg,'COLL') 905 ! Allocation 906 ABI_MALLOC(vdw_c6,(natom,natom)) 907 ABI_MALLOC(vdw_c8,(natom,natom)) 908 ABI_MALLOC(dc6ri,(natom,natom)) 909 ABI_MALLOC(dc6rj,(natom,natom)) 910 ABI_MALLOC(d2c6ri,(natom,natom)) 911 ABI_MALLOC(d2c6rj,(natom,natom)) 912 ABI_MALLOC(d2c6rirj,(natom,natom)) 913 if (bol_3bt) then 914 ABI_MALLOC(vdw_c9,(natom,natom,natom)) 915 ABI_MALLOC(dc9ijri,(natom,natom)) 916 ABI_MALLOC(dc9ijrj,(natom,natom)) 917 end if 918 ! Set accumulating quantities to zero 919 vdw_c6 = zero ; vdw_c8 = zero 920 dc6ri = zero ; dc6rj = zero 921 d2c6ri = zero ; d2c6rj = zero 922 d2c6rirj = zero 923 if (bol_3bt) then 924 dc9ijri = zero; dc9ijrj = zero 925 end if 926 ! C6 coefficients are interpolated from tabulated 927 ! ab initio C6 values (following loop). 928 ! C8 coefficients are obtained by: 929 ! C8 = vdw_dftd3_q(itypat)*vdw_dftd3_q(jtypat)*C6 930 931 do ia=1,natom 932 itypat=typat(ia) 933 do ja=1,natom 934 jtypat=typat(ja) 935 ! Set accumulating quantities to zero 936 ltot=zero 937 sum_dlri = zero ; sum_dlc6ri= zero 938 sum_dlrj = zero ; sum_dlc6rj= zero 939 sum_d2lri = zero ; sum_d2lc6ri = zero 940 sum_d2lrj = zero ; sum_d2lc6rj = zero 941 sum_d2lrirj = zero ; sum_d2lc6rirj = zero 942 min_dsys = 10000 943 max_vdw_c6 = zero 944 ! Loop over references 945 do refi=1,refmax 946 do refj=1,refmax 947 dsysref_a = cn(ia)-vdw_cnrefi(itypat,jtypat,refi,refj) 948 dsysref_b = cn(ja)-vdw_cnrefj(itypat,jtypat,refi,refj) 949 dsysref=(dsysref_a)**two+(dsysref_b)**two 950 if (dsysref<min_dsys) then 951 ! Keep in memory the smallest value of dsysref 952 ! And the associated tabulated C6 value 953 min_dsys = dsysref 954 max_vdw_c6 = vdw_c6ref(itypat,jtypat,refi,refj) 955 end if 956 l = dexp(-k3*dsysref) 957 ltot = ltot+l 958 vdw_c6(ia,ja)=vdw_c6(ia,ja)+vdw_c6ref(itypat,jtypat,refi,refj)*l 959 960 if (need_grad) then 961 ! Derivative of l(ia,ja) with respect to the displacement 962 ! of atom ka in reduced coordinates. 963 ! This factor is identical in the case of stress. 964 ! In purpose of speed up this routine, the prefactor of 965 ! dCNi/drk and dCNj/drk are separated 966 ! See NOTE: article to be added 967 dlri=-k3*l*two*dsysref_a ;dlrj=-k3*l*two*dsysref_b 968 sum_dlri=sum_dlri+dlri ; sum_dlrj=sum_dlrj+dlrj 969 sum_dlc6ri=sum_dlc6ri+dlri*vdw_c6ref(itypat,jtypat,refi,refj) 970 sum_dlc6rj=sum_dlc6rj+dlrj*vdw_c6ref(itypat,jtypat,refi,refj) 971 if (need_hess) then 972 ! Second derivative of l(ia,ja). Once again, it is separately in 973 ! different contributions: 974 ! d2lri: prefactor of dCNi/drk*dCNi/drl 975 ! d2lrj: prefactor of dCNj/drk*dCNj/drl 976 ! d2lrirj: prefacto of dCNi/drk*dCNj/drl 977 ! The prefactor for d2CNi/drkdrl is dlri; for d2CNj/drkdrl is dlrj 978 d2lri = -two*k3*l*(one-two*k3*dsysref_a**two) 979 d2lrj = -two*k3*l*(one-two*k3*dsysref_b**two) 980 d2lrirj = four*k3*k3*l*(dsysref_a*dsysref_b) 981 sum_d2lri=sum_d2lri+d2lri ; sum_d2lrj=sum_d2lrj+d2lrj 982 sum_d2lrirj = sum_d2lrirj+d2lrirj 983 sum_d2lc6ri=sum_d2lc6ri+d2lri*vdw_c6ref(itypat,jtypat,refi,refj) 984 sum_d2lc6rj=sum_d2lc6rj+d2lrj*vdw_c6ref(itypat,jtypat,refi,refj) 985 sum_d2lc6rirj = sum_d2lc6rirj + d2lrirj*vdw_c6ref(itypat,jtypat,refi,refj) 986 end if ! Boolean second derivative 987 end if ! Boolean gradient 988 end do ! Loop over references 989 end do ! Loop over references 990 ! In some specific case (really covalently bound compounds) ltot -> 0 991 ! which may cause numerical problems for all quantities related to dispersion coefficient. 992 ! To be consistent with VASP implementation, the c6 value is taken as the last 993 ! referenced value of the dispersion coefficient. 994 if (ltot>tol12) then 995 vdw_c6(ia,ja)=vdw_c6(ia,ja)/ltot 996 vdw_c8(ia,ja)=three*vdw_q_dftd3(ivdw(itypat))*vdw_q_dftd3(ivdw(jtypat))*vdw_c6(ia,ja) 997 ! If Force of Stress is required 998 if (need_grad) then 999 ! Computation of the derivative of C6 w.r.to the displacement 1000 ! of atom ka, in reduced coordinates (separated for dCNi/drk and dCNj/drk) 1001 ! This is the crucial step to reduce the scaling from O(N^3) to O(N^2) for 1002 ! the gradients 1003 dc6ri(ia,ja)=(sum_dlc6ri-vdw_c6(ia,ja)*sum_dlri)/ltot 1004 dc6rj(ia,ja)=(sum_dlc6rj-vdw_c6(ia,ja)*sum_dlrj)/ltot 1005 if (need_hess) then 1006 ! Computation of the second derivative of C6 w.r.to the displacement of atom ka 1007 ! and atom la 1008 d2c6ri(ia,ja)=(sum_d2lc6ri-vdw_c6(ia,ja)*sum_d2lri-two*dc6ri(ia,ja)*sum_dlri)/ltot 1009 d2c6rj(ia,ja)=(sum_d2lc6rj-vdw_c6(ia,ja)*sum_d2lrj-two*dc6rj(ia,ja)*sum_dlrj)/ltot 1010 d2c6rirj(ia,ja) = (sum_d2lc6rirj-vdw_c6(ia,ja)*sum_d2lrirj-dc6ri(ia,ja)*& 1011 & sum_dlrj-dc6rj(ia,ja)*sum_dlri)/ltot 1012 end if ! Boolean second derivative 1013 end if ! Boolean gradient 1014 else 1015 vdw_c6(ia,ja)= max_vdw_c6 1016 vdw_c8(ia,ja)=three*vdw_q_dftd3(ivdw(itypat))*vdw_q_dftd3(ivdw(jtypat))*vdw_c6(ia,ja) 1017 end if 1018 end do 1019 end do 1020 ! Computation of the three-body term dispersion coefficient 1021 if (bol_3bt) then 1022 do ia=1,natom 1023 do ja=1,natom 1024 do ka=1,natom 1025 vdw_c9(ia,ja,ka) =-sqrt(vdw_c6(ia,ja)*vdw_c6(ja,ka)*vdw_c6(ka,ia)) 1026 end do 1027 if (need_grad) then 1028 dc9ijri(ia,ja) = half/vdw_c6(ia,ja)*dc6ri(ia,ja) 1029 dc9ijrj(ia,ja) = half/vdw_c6(ia,ja)*dc6rj(ia,ja) 1030 end if 1031 end do 1032 end do 1033 end if 1034 1035 write(msg,'(3a,f10.5,1a,f10.5)')& 1036 & ' ... Done.',ch10,& 1037 & ' max(C6) =', maxval(vdw_c6),' ; min(C6) =', minval(vdw_c6) 1038 call wrtout(std_out,msg,'COLL') 1039 1040 ! Deallocation of used variables not needed anymore 1041 ABI_FREE(vdw_c6ref) 1042 ABI_FREE(vdw_cnrefi) 1043 ABI_FREE(vdw_cnrefj) 1044 1045 !---------------------------------------------------- 1046 ! Computation of cut-off radii according to tolerance 1047 !---------------------------------------------------- 1048 1049 ! Cut-off radius for pair-wise term 1050 if (vdw_tol<zero) then 1051 rcut=max((vdw_s6/vdw_tol_default*maxval(vdw_c6))**sixth, & 1052 & (vdw_s8/vdw_tol_default*maxval(vdw_c8))**(one/eight)) 1053 else 1054 rcut=max((vdw_s6/vdw_tol*maxval(vdw_c6))**sixth,& 1055 & (vdw_s8/vdw_tol*maxval(vdw_c8))**(one/eight)) 1056 end if 1057 ! Cut-off radius for three-body term 1058 rcut9 = zero 1059 if (bol_3bt) then 1060 rcut9=(128.0_dp*vdw_s6/(vdw_tol_3bt)*maxval(vdw_c6)**(3.0/2.0))**(1.0/9.0) 1061 end if 1062 rcut2=rcut*rcut 1063 1064 !-------------------------------------------------------------------- 1065 ! Computation of the two bodies contribution to the dispersion energy 1066 !-------------------------------------------------------------------- 1067 1068 write(msg,'(3a)')& 1069 & ' Begin the computation of pair-wise term',ch10,& 1070 & ' of DFT-D3 energy contribution...' 1071 call wrtout(std_out,msg,'COLL') 1072 nshell=0 1073 npairs=0 1074 ABI_MALLOC(e_alpha1,(natom)) 1075 ABI_MALLOC(e_alpha2,(natom)) 1076 ABI_MALLOC(e_alpha3,(natom)) 1077 ABI_MALLOC(e_alpha4,(natom)) 1078 ABI_MALLOC(e_no_c,(natom,natom)) 1079 ABI_MALLOC(fe_no_c,(2,natom,natom)) 1080 e_alpha1 =zero ; e_alpha2 = zero 1081 e_alpha3 =zero ; e_alpha4 = zero 1082 e_no_c=zero ; fe_no_c = zero 1083 ABI_MALLOC(grad_no_cij,(3,natom,natom)) 1084 ABI_MALLOC(fgrad_no_c,(2,3,natom,natom)) 1085 ABI_MALLOC(cfgrad_no_c,(2,3,natom,natom)) 1086 ABI_MALLOC(str_no_c,(6,natom,natom)) 1087 ABI_MALLOC(str_alpha1,(6,natom)) 1088 ABI_MALLOC(str_alpha2,(6,natom)) 1089 grad_no_cij=zero ; str_no_c=zero 1090 fgrad_no_c = zero ; cfgrad_no_c = zero 1091 str_alpha1 = zero ; str_alpha2 = zero 1092 1093 re_arg = zero ; im_arg = zero 1094 dmp6 = zero ; dmp8 = zero 1095 e_no_c6 = zero ; e_no_c8 = zero 1096 fdmp6 = zero ; fdmp8 = zero 1097 grad6 = zero ; grad8 = zero 1098 grad6_no_c6 = zero ; grad8_no_c8 = zero 1099 hess6 = zero ; hess8 = zero 1100 do 1101 newshell=.false.;nshell=nshell+1 1102 do is3=-nshell,nshell 1103 do is2=-nshell,nshell 1104 do is1=-nshell,nshell 1105 if (nshell==1.or.abs(is3)==nshell.or.abs(is2)==nshell.or.abs(is1)==nshell) then 1106 is(1) = is1 ; is(2) = is2 ; is(3) = is3 1107 ! Computation of phase factor for discrete Fourier transform 1108 ! if phonon at qphon is required 1109 if (need_dynmat) then 1110 arg= two_pi*dot_product(qphon,is) 1111 re_arg=cos(arg) 1112 im_arg=sin(arg) 1113 end if 1114 do ia=1,natom 1115 itypat=typat(ia) 1116 do ja=1,natom 1117 jtypat=typat(ja) 1118 r(:)=xred01(:,ia)-xred01(:,ja)-dble(is(:)) 1119 rsq=dot_product(r,matmul(rmet,r)) 1120 if (rsq>=tol16.and.rsq<rcut2) then 1121 npairs=npairs+1;newshell=.true. 1122 sfact6 = half*vdw_s6 ; sfact8 = half*vdw_s8 1123 rr=sqrt(rsq); r6 = rr**six ; r8 = rr**eight 1124 c6=vdw_c6(ia,ja) ; c8=vdw_c8(ia,ja) 1125 vdw_q = three*vdw_q_dftd3(ivdw(itypat))*vdw_q_dftd3(ivdw(jtypat)) 1126 r0=vdw_r0(itypat,jtypat) 1127 ! Computation of e_vdw_dftd3 (case DFT+D3) 1128 if (vdw_xc == 6) then 1129 dmp6=six*(rr/(vdw_sr6*r0))**(-alpha6) 1130 fdmp6=one/(one+dmp6) 1131 dmp8=six*(rr/(vdw_sr8*r0))**(-alpha8) 1132 fdmp8=one/(one+dmp8) 1133 ! Contribution to energy 1134 e_no_c6 = -sfact6*fdmp6/r6 ; e_no_c8 = -sfact8*fdmp8/r8 1135 e_vdw_dftd3=e_vdw_dftd3+e_no_c6*c6 +e_no_c8*c8 1136 ! Computation of e_vdw_dftd3 (case DFT+D3-BJ) 1137 elseif (vdw_xc == 7) then 1138 dmp = (vdw_a1*sqrt(vdw_q)+vdw_a2) 1139 fdmp6 = one/(dmp**six+rr**six) 1140 fdmp8 = one/(dmp**eight+rr**eight) 1141 e_no_c6 = -sfact6*fdmp6 ; e_no_c8 = -sfact8*fdmp8 1142 e_vdw_dftd3=e_vdw_dftd3-sfact6*c6*fdmp6-sfact8*c8*fdmp8 1143 end if 1144 ! Computation of the gradients (if required) 1145 if (need_grad) then 1146 if (vdw_xc == 6) then 1147 gr6 = alpha6*dmp6*fdmp6**two 1148 grad6_no_c6 = sfact6*(gr6-six*fdmp6)/r8 1149 grad6 = grad6_no_c6*c6 1150 gr8 = alpha8*dmp8*fdmp8**two 1151 grad8_no_c8 = sfact8*(gr8-eight*fdmp8)/r8/rsq 1152 grad8 = grad8_no_c8*c8 1153 elseif (vdw_xc == 7) then 1154 grad6_no_c6 = -sfact6*six*(fdmp6*rsq)**two 1155 grad6 = grad6_no_c6*c6 1156 grad8_no_c8 = -sfact8*eight*(fdmp8)**two*rsq**three 1157 grad8 = grad8_no_c8*c8 1158 end if 1159 grad =grad6+grad8 1160 grad_no_c = grad6_no_c6+grad8_no_c8*vdw_q 1161 rcart=matmul(rprimd,r) 1162 rred= matmul(transpose(rprimd),rcart) 1163 ! Additional contribution due to c6(cn(r)) 1164 ! Not yet multiply by dCN/drk and summed to reduce 1165 ! computational time 1166 e_no_c(ia,ja) = e_no_c(ia,ja)+(e_no_c6+vdw_q*e_no_c8) 1167 ! Part related to alpha1ij/alpha2ij 1168 e_alpha1(ia) = e_alpha1(ia)+(e_no_c6+vdw_q*e_no_c8)*dc6ri(ia,ja) 1169 e_alpha2(ja) = e_alpha2(ja)+(e_no_c6+vdw_q*e_no_c8)*dc6rj(ia,ja) 1170 ! Contribution to gradients wr to atomic displacement 1171 ! (forces) 1172 if (need_forces.and.ia/=ja) then 1173 gred(:)=grad*rred(:) 1174 do alpha=1,3 1175 gred_vdw_dftd3(alpha,ia)=gred_vdw_dftd3(alpha,ia)-gred(alpha) 1176 gred_vdw_dftd3(alpha,ja)=gred_vdw_dftd3(alpha,ja)+gred(alpha) 1177 end do 1178 elseif (need_stress) then 1179 ! Computation of the DFT-D3 contribution to stress 1180 vec(1:3)=rcart(1:3)*rcart(1:3); vec(4)=rcart(2)*rcart(3) 1181 vec(5)=rcart(1)*rcart(3); vec(6)=rcart(1)*rcart(2) 1182 do alpha=1,6 1183 str_vdw_dftd3(alpha)=str_vdw_dftd3(alpha)-grad*vec(alpha) 1184 end do 1185 end if 1186 ! Second derivative (if required) 1187 if (need_hess) then 1188 if (vdw_xc==6) then 1189 hess6 = (grad6*(alpha6*fdmp6*dmp6-8.0_dp)+& 1190 & sfact6*c6/r6*dmp6*((alpha6*fdmp6)**two)*& 1191 & (fdmp6*dmp6-one)/rsq)/rsq 1192 hess8 = (grad8*(alpha8*fdmp8*dmp8-10.0_dp)+& 1193 & sfact8*c8/r8*dmp8*((alpha8*fdmp8)**two)*& 1194 & (fdmp8*dmp8-one)/rsq)/rsq 1195 elseif (vdw_xc==7) then 1196 hess6 = -four*grad6*(three*rsq**two*fdmp6-one/rsq) 1197 hess8 = -two*grad8*(eight*rsq**three*fdmp8-three/rsq) 1198 end if 1199 ! Contribution of d2C6 to the interatomic force constants 1200 ! Not yet multiply by CN derivative and summed to reduce the scaling from O(N^3) to O(N^2) 1201 hessij = hess6+hess8 1202 ! Contribution of cross-derivative dC6 and grad 1203 do alpha=1,3 1204 grad_no_cij(alpha,ia,ja) = grad_no_cij(alpha,ia,ja) - grad_no_c*rcart(alpha) 1205 end do 1206 e_alpha3(ia) = e_alpha3(ia)+(e_no_c6+vdw_q*e_no_c8)*d2c6ri(ia,ja) 1207 e_alpha4(ja) = e_alpha4(ja)+(e_no_c6+vdw_q*e_no_c8)*d2c6rj(ia,ja) 1208 if (need_dynmat) then 1209 ! Fourier transform of the partial contribution to the dispersion potential 1210 fe_no_c(1,ia,ja) = fe_no_c(1,ia,ja)+(e_no_c6+vdw_q*e_no_c8)*re_arg 1211 fe_no_c(2,ia,ja) = fe_no_c(2,ia,ja)+(e_no_c6+vdw_q*e_no_c8)*im_arg 1212 do alpha=1,3 1213 ! Fourier transform of the gradient (required for the IFCs) 1214 fgrad_no_c(1,alpha,ia,ja) = fgrad_no_c(1,alpha,ia,ja)-grad_no_c*rcart(alpha)*re_arg 1215 fgrad_no_c(2,alpha,ia,ja) = fgrad_no_c(2,alpha,ia,ja)-grad_no_c*rcart(alpha)*im_arg 1216 ! Complex conjugated of the Fourier transform of the gradient 1217 cfgrad_no_c(1,alpha,ia,ja) = cfgrad_no_c(1,alpha,ia,ja)-grad_no_c*rcart(alpha)*re_arg 1218 cfgrad_no_c(2,alpha,ia,ja) = cfgrad_no_c(2,alpha,ia,ja)+grad_no_c*rcart(alpha)*im_arg 1219 end do 1220 ! Contribution to the IFCs (reciprocal space) of the 2nd derivative of e_no_c part 1221 do alpha=1,3 1222 do beta=1,3 1223 rcart2(alpha,beta) = rcart(alpha)*rcart(beta) 1224 end do 1225 end do 1226 if (ia/=ja) then 1227 do alpha=1,3 1228 dyn_vdw_dftd3(1,alpha,ja,:,ja) = dyn_vdw_dftd3(1,alpha,ja,:,ja) -& 1229 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:)) 1230 dyn_vdw_dftd3(1,alpha,ia,:,ia) = dyn_vdw_dftd3(1,alpha,ia,:,ia) -& 1231 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:)) 1232 if (abs(re_arg)>tol12) then 1233 dyn_vdw_dftd3(1,alpha,ia,:,ja) = dyn_vdw_dftd3(1,alpha,ia,:,ja) +& 1234 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:))*re_arg 1235 dyn_vdw_dftd3(1,alpha,ja,:,ia) = dyn_vdw_dftd3(1,alpha,ja,:,ia) +& 1236 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:))*re_arg 1237 end if 1238 if (abs(im_arg)>tol12) then 1239 dyn_vdw_dftd3(2,alpha,ia,:,ja) = dyn_vdw_dftd3(2,alpha,ia,:,ja) +& 1240 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:))*im_arg 1241 dyn_vdw_dftd3(2,alpha,ja,:,ia) = dyn_vdw_dftd3(2,alpha,ja,:,ia) -& 1242 & (hessij*rcart2(alpha,:)+grad*mcart(alpha,:))*im_arg 1243 end if 1244 end do 1245 else ! ia==ja 1246 do alpha=1,3 1247 if (abs(re_arg-one)>tol12) then 1248 dyn_vdw_dftd3(1,alpha,ia,:,ia) = dyn_vdw_dftd3(1,alpha,ia,:,ia) -& 1249 & two*(hessij*rcart2(alpha,:)+grad*mcart(alpha,:))*(one-re_arg) 1250 end if 1251 end do 1252 end if 1253 end if 1254 ! Now compute the contribution to the elastic constants !!! Still under development 1255 if (need_elast) then 1256 vec(1:3)=rcart(1:3)*rcart(1:3); vec(4)=rcart(2)*rcart(3) 1257 vec(5)=rcart(1)*rcart(3); vec(6)=rcart(1)*rcart(2) 1258 str_no_c(:,ia,ja)=str_no_c(:,ia,ja)-grad_no_c*vec(:) 1259 str_alpha1(:,ia)=str_alpha1(:,ia)-dc6ri(ia,ja)*grad_no_c*vec(:) 1260 str_alpha2(:,ja)=str_alpha2(:,ja)-dc6rj(ia,ja)*grad_no_c*vec(:) 1261 ! Contribution to elastic constants of DFT-D3 dispersion potential (no C6 derivative) 1262 do alpha=1,6 1263 ii = voigt1(alpha) ; jj=voigt2(alpha) 1264 do beta=1,6 1265 kk = voigt1(beta) ; ll=voigt2(beta) 1266 elt_vdw_dftd3(alpha,beta) = elt_vdw_dftd3(alpha,beta)-hessij*vec(alpha)*vec(beta) 1267 if (ii==kk) elt_vdw_dftd3(alpha,beta) = elt_vdw_dftd3(alpha,beta)-half*grad*rcart(jj)*rcart(ll) 1268 if (jj==kk) elt_vdw_dftd3(alpha,beta) = elt_vdw_dftd3(alpha,beta)-half*grad*rcart(ii)*rcart(ll) 1269 if (ii==ll) elt_vdw_dftd3(alpha,beta) = elt_vdw_dftd3(alpha,beta)-half*grad*rcart(jj)*rcart(kk) 1270 if (jj==ll) elt_vdw_dftd3(alpha,beta) = elt_vdw_dftd3(alpha,beta)-half*grad*rcart(ii)*rcart(kk) 1271 end do 1272 end do 1273 ! Contribution to internal strain of DFT-D3 dispersion potential (no C6 derivative) 1274 if (ia/=ja) then 1275 index_ia = 6+3*(ia-1) 1276 index_ja = 6+3*(ja-1) 1277 do alpha=1,6 1278 do beta=1,3 1279 ii = voigt1(alpha) ; jj=voigt2(alpha) 1280 elt_vdw_dftd3(index_ia+beta,alpha)=elt_vdw_dftd3(index_ia+beta,alpha)-& 1281 & hessij*vec(alpha)*rcart(beta) 1282 elt_vdw_dftd3(index_ja+beta,alpha)=elt_vdw_dftd3(index_ja+beta,alpha)+& 1283 & hessij*vec(alpha)*rcart(beta) 1284 if (ii==beta) then 1285 elt_vdw_dftd3(index_ia+beta,alpha)=elt_vdw_dftd3(index_ia+beta,alpha)-grad*rcart(jj) 1286 elt_vdw_dftd3(index_ja+beta,alpha)=elt_vdw_dftd3(index_ja+beta,alpha)+grad*rcart(jj) 1287 end if 1288 if (jj==beta) then 1289 elt_vdw_dftd3(index_ia+beta,alpha)=elt_vdw_dftd3(index_ia+beta,alpha)-grad*rcart(ii) 1290 elt_vdw_dftd3(index_ja+beta,alpha)=elt_vdw_dftd3(index_ja+beta,alpha)+grad*rcart(ii) 1291 end if 1292 end do ! Direction beta 1293 end do ! Strain alpha 1294 end if ! ia/=ja 1295 end if ! Need elastic constant 1296 end if ! Need hessian 1297 end if ! Need gradient 1298 end if ! Tolerance 1299 end do ! Loop over atom j 1300 end do ! Loop over atom i 1301 end if ! Triple loop over cell replicas in shell 1302 end do ! Is1 1303 end do ! Is2 1304 end do ! Is3 1305 if(.not.newshell) exit ! Check if new shell must be calculated 1306 end do ! Loop over shell 1307 ABI_MALLOC(temp_prod,(2,natom)) 1308 if (need_grad) then 1309 ! Additional contribution to force due dc6_drk 1310 if (need_forces) then 1311 do ka=1,natom 1312 do ia=1,natom 1313 !do ja=1,natom 1314 !gred_vdw_dftd3(:,ka) = gred_vdw_dftd3(:,ka)+e_no_c(ia,ja)*(& 1315 !& !dcn(:,ia,ka)*dc6ri(ia,ja)+dcn(:,ja,ka)*dc6rj(ia,ja)) 1316 !end do 1317 gred_vdw_dftd3(:,ka) = gred_vdw_dftd3(:,ka)+e_alpha1(ia)*dcn(:,ia,ka)+& 1318 & e_alpha2(ia)*dcn(:,ia,ka) 1319 end do 1320 end do 1321 elseif (need_stress) then 1322 do ia=1,natom 1323 !do ja=1,natom 1324 ! str_vdw_dftd3(:) = str_vdw_dftd3(:)+e_no_c(ia,ja)*(str_dcn(:,ia)*& 1325 !& ! dc6ri(ia,ja)+str_dcn(:,ja)*dc6rj(ia,ja)) 1326 !end do 1327 str_vdw_dftd3(:) = str_vdw_dftd3(:)+e_alpha1(ia)*str_dcn(:,ia)+& 1328 & e_alpha2(ia)*str_dcn(:,ia) 1329 end do 1330 end if ! Optimization 1331 ! If dynmat is required, add all the terms related to dc6, d2c6, ... 1332 if (need_hess) then 1333 if (need_dynmat) then 1334 do ka=1,natom 1335 do la =1,natom 1336 do alpha=1,3 1337 do beta=1,3 1338 do ia=1,natom 1339 !TODO: avoid stupid if clauses inside the loops 1340 d2cn_tmp = zero 1341 if (ia==la) then 1342 if (ia==ka) then ! iii 1343 d2cn_tmp(:) = d2cn_iii(:,alpha,beta,ia) 1344 else ! jii 1345 d2cn_tmp(:) = d2cn_jii(:,alpha,beta,ka,ia) 1346 end if 1347 else if (ia==ka) then !iji 1348 d2cn_tmp(:) = d2cn_iji(:,alpha,beta,la,ia) 1349 else if (ka==la) then ! jji 1350 d2cn_tmp(:) = d2cn_jji(:,alpha,beta,ka,ia) 1351 end if 1352 ! Add the second derivative of C6 contribution to the dynamical matrix 1353 ! First, add the second derivative of CN-related term 1354 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1355 & (e_alpha1(ia)+e_alpha2(ia))*d2cn_tmp(:) 1356 !OLDVERSION dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1357 !& (e_alpha1(ia)+e_alpha2(ia))*d2cn(:,alpha,ka,beta,la,ia) 1358 1359 1360 1361 ! Then the term related to dCNi/dr*dCNi/dr and dCNj/dr*dCNj/dr 1362 call comp_prod(cfdcn(:,alpha,ia,ka),fdcn(:,beta,ia,la),temp_comp) 1363 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1364 & (e_alpha3(ia)+e_alpha4(ia))*temp_comp(:) 1365 ! Add the cross derivative of fdmp/rij**6 and C6 contribution to the dynamical matrix 1366 ! !!!! The products are kind of tricky... 1367 ! First, add the dCNk/drl gradik and dCNk/drl gradjk terms... 1368 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1369 & cfdcn(:,alpha,la,ka)*grad_no_cij(beta,la,ia)*(dc6ri(la,ia)+dc6rj(ia,la)) 1370 ! Then the dCNk/drl gradjk and dCNk/drl gradik terms... 1371 call comp_prod(cfdcn(:,alpha,ia,ka),cfgrad_no_c(:,beta,la,ia),temp_comp) 1372 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1373 & temp_comp(:)*(dc6ri(ia,la)+dc6rj(la,ia)) 1374 ! Here the symmetrical term (for dCNl/drk) are added... 1375 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1376 & fdcn(:,beta,ka,la)*grad_no_cij(alpha,ka,ia)*(dc6ri(ka,ia)+dc6rj(ia,ka)) 1377 call comp_prod(fdcn(:,beta,ia,la),fgrad_no_c(:,alpha,ka,ia),temp_comp) 1378 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1379 & temp_comp(:)*(dc6ri(ia,ka)+dc6rj(ka,ia)) 1380 end do ! ia 1381 end do ! alpha 1382 end do ! beta 1383 end do ! la 1384 do alpha=1,3 1385 temp_prod(:,:) = zero 1386 do ja=1,natom 1387 do ia=1,natom 1388 ! Finally the cross derivative dCNi/dr dCNj/dr 1389 temp_comp2(:) = d2c6rirj(ia,ja)*fe_no_c(:,ia,ja) 1390 call comp_prod(cfdcn(:,alpha,ia,ka),temp_comp2,temp_comp) 1391 temp_prod(:,ja) = temp_prod(:,ja)+temp_comp(:) 1392 end do 1393 do la = 1,natom 1394 do beta=1,3 1395 call comp_prod(fdcn(:,beta,ja,la),temp_prod(:,ja),temp_comp2) 1396 dyn_vdw_dftd3(:,alpha,ka,beta,la)=dyn_vdw_dftd3(:,alpha,ka,beta,la)+& 1397 & two*temp_comp2 1398 end do ! beta 1399 end do ! la 1400 end do ! ja 1401 end do ! alpha 1402 end do !ka 1403 ! Transformation from cartesian coordinates to reduced coordinates 1404 do ka=1,natom 1405 do la=1,natom 1406 do kk=1,2 1407 do alpha=1,3 1408 vec(1:3)=dyn_vdw_dftd3(kk,1:3,ka,alpha,la) 1409 call d3_cart2red(vec) 1410 dyn_vdw_dftd3(kk,1:3,ka,alpha,la)=vec(1:3) 1411 end do 1412 do alpha=1,3 1413 vec(1:3)=dyn_vdw_dftd3(kk,alpha,ka,1:3,la) 1414 call d3_cart2red(vec) 1415 dyn_vdw_dftd3(kk,alpha,ka,1:3,la)=vec(1:3) 1416 end do ! alpha 1417 end do ! real/im 1418 end do ! Atom la 1419 end do ! Atom ka 1420 end if ! Boolean dynamical matrix 1421 if (need_elast) then 1422 do ia=1,natom 1423 index_ia = 6+3*(ia-1) 1424 do alpha=1,6 1425 ! Add the second derivative of C6 contribution to the elastic tensor 1426 ! First, the second derivative of CN with strain 1427 elt_vdw_dftd3(alpha,:) = elt_vdw_dftd3(alpha,:)+(& 1428 & e_alpha1(ia)+e_alpha2(ia))*elt_cn(alpha,:,ia) 1429 ! Then, the derivative product dCNi/deta dCNi/deta 1430 elt_vdw_dftd3(alpha,:) = elt_vdw_dftd3(alpha,:)+(& 1431 & e_alpha3(ia)+e_alpha4(ia))*str_dcn(alpha,ia)*str_dcn(:,ia) 1432 ! Then, the dCNi/deta dCNj/deta 1433 do ja=1,natom 1434 elt_vdw_dftd3(alpha,:)=elt_vdw_dftd3(alpha,:)+two*e_no_c(ia,ja)*& 1435 & d2c6rirj(ia,ja)*str_dcn(alpha,ia)*str_dcn(:,ja) 1436 end do 1437 ! Add the cross derivative of fij and C6 contribution to the elastic tensor 1438 elt_vdw_dftd3(alpha,:)=elt_vdw_dftd3(alpha,:)+str_alpha1(:,ia)*str_dcn(alpha,ia)+str_alpha1(alpha,ia)*str_dcn(:,ia)& 1439 & +str_dcn(alpha,ia)*str_alpha2(:,ia)+str_dcn(:,ia)*str_alpha2(alpha,ia) 1440 end do 1441 do alpha=1,6 1442 ii = voigt1(alpha) ; jj=voigt2(alpha) 1443 do ka=1,natom 1444 index_ka = 6+3*(ka-1) 1445 do beta=1,3 1446 ! Add the second derivative of C6 contribution to the internal strains 1447 ii = voigt1(alpha) ; jj=voigt2(alpha) 1448 ! Second derivative of CN 1449 elt_vdw_dftd3(index_ka+beta,alpha)=elt_vdw_dftd3(index_ka+beta,alpha)+(e_alpha1(ia)+e_alpha2(ia))*& 1450 & elt_cn(index_ka+beta,alpha,ia) 1451 ! Cross-derivatives of CN 1452 elt_vdw_dftd3(index_ka+beta,alpha)=elt_vdw_dftd3(index_ka+beta,alpha)+(e_alpha3(ia)+e_alpha4(ia))*& 1453 & dcn_cart(beta,ia,ka)*str_dcn(alpha,ia) !OK 1454 do ja=1,natom 1455 index_ja = 6+3*(ja-1) 1456 elt_vdw_dftd3(index_ka+beta,alpha)=elt_vdw_dftd3(index_ka+beta,alpha)+two*d2c6rirj(ia,ja)*e_no_c(ia,ja)*& 1457 & dcn_cart(beta,ia,ka)*str_dcn(alpha,ja) 1458 end do 1459 elt_vdw_dftd3(index_ka+beta,alpha)=elt_vdw_dftd3(index_ka+beta,alpha)+(str_alpha1(alpha,ia)+str_alpha2(alpha,ia))*& 1460 & dcn_cart(beta,ia,ka) 1461 elt_vdw_dftd3(index_ka+beta,alpha)=elt_vdw_dftd3(index_ka+beta,alpha)+grad_no_cij(beta,ka,ia)*& 1462 & (dc6ri(ia,ka)*str_dcn(alpha,ia)+dc6rj(ia,ka)*str_dcn(alpha,ka))-grad_no_cij(beta,ia,ka)*& 1463 & (dc6ri(ka,ia)*str_dcn(alpha,ka)+dc6rj(ka,ia)*str_dcn(alpha,ia)) 1464 end do ! beta 1465 end do ! ka 1466 end do ! alpha 1467 end do ! ia 1468 do alpha=1,6 1469 index_ia=6 1470 do ia=1,natom 1471 elt_vdw_dftd3(index_ia+1:index_ia+3,alpha)=matmul(transpose(rprimd),elt_vdw_dftd3(index_ia+1:index_ia+3,alpha)) 1472 index_ia=index_ia+3 1473 end do ! Atom ia 1474 end do ! Strain alpha 1475 end if ! Boolean elastic tensor 1476 end if ! Boolean hessian 1477 end if ! Boolean need_gradient 1478 ABI_FREE(temp_prod) 1479 write(msg,'(3a)')& 1480 & ' ...Done.' 1481 call wrtout(std_out,msg,'COLL') 1482 1483 !print *, 'Evdw', e_vdw_dftd3 1484 !if (need_forces) print *, 'fvdw', gred_vdw_dftd3(3,:) 1485 !if (need_stress) print *, 'strvdw', str_vdw_dftd3(6) 1486 !if (need_elast) print *, 'Elast(3,3,3,3)', elt_vdw_dftd3(3,3) 1487 !if (need_elast) print *, 'Elast(3,3,3,3)', elt_vdw_dftd3(6,6) 1488 !if (need_elast) print *, 'Elast(3,3,3,3)', elt_vdw_dftd3(3,6) 1489 !if (need_elast) print *, 'Internal(1,3,3)', elt_vdw_dftd3(9,3) 1490 !if (need_elast) print *, 'Internal(1,3,6)', elt_vdw_dftd3(9,6) 1491 !if (need_dynmat) print *, 'Dynmat(1,3,:,3,:)', dyn_vdw_dftd3(1,3,:,3,:) 1492 !if (need_dynmat) print *, 'Dynmat(1,3,:,3,:)', dyn_vdw_dftd3(1,2,:,1,:) 1493 1494 !--------------------------------------------- 1495 ! Computation of the 3 body term (if required) 1496 !--------------------------------------------- 1497 1498 e_3bt=zero 1499 if (bol_3bt) then 1500 if (need_grad) then 1501 ABI_MALLOC(e3bt_ij,(natom,natom)) 1502 ABI_MALLOC(e3bt_jk,(natom,natom)) 1503 ABI_MALLOC(e3bt_ki,(natom,natom)) 1504 ABI_MALLOC(gred_vdw_3bt,(3,natom)) 1505 e3bt_ij=zero; e3bt_jk=zero; e3bt_ki=zero 1506 if (need_forces) then 1507 gred_vdw_3bt = zero 1508 elseif (need_stress) then 1509 str_3bt=zero 1510 end if 1511 end if 1512 nshell_3bt(1) = int(0.5+rcut9/sqrt(rmet(1,1)+rmet(2,1)+rmet(3,1))) 1513 nshell_3bt(2) = int(0.5+rcut9/sqrt(rmet(1,2)+rmet(2,2)+rmet(3,2))) 1514 nshell_3bt(3) = int(0.5+rcut9/sqrt(rmet(1,3)+rmet(2,3)+rmet(3,3))) 1515 1516 do is3 = -nshell_3bt(3),nshell_3bt(3) 1517 do is2 = -nshell_3bt(2),nshell_3bt(2) 1518 do is1 = -nshell_3bt(1),nshell_3bt(1) 1519 is(1) = is1 ; is(2)=is2 ; is(3) = is3 1520 do alpha=1,3 1521 jmin(alpha) = max(-nshell_3bt(alpha), -nshell_3bt(alpha)+is(alpha)) 1522 jmax(alpha) = min(nshell_3bt(alpha), nshell_3bt(alpha)+is(alpha)) 1523 end do 1524 do js3=jmin(3),jmax(3) 1525 do js2=jmin(2),jmax(2) 1526 do js1=jmin(1),jmax(1) 1527 js(1) = js1 ; js(2)=js2 ; js(3) = js3 1528 do ia=1,natom 1529 itypat=typat(ia) 1530 do ja=1,ia 1531 jtypat=typat(ja) 1532 do ka=1,ja 1533 ktypat=typat(ka) 1534 rij(:) = xred01(:,ia)-xred01(:,ja)-dble(is(:)) 1535 rsqij = dot_product(rij(:),matmul(rmet,rij(:))) 1536 rrij = dsqrt(rsqij) 1537 rjk(:) = xred01(:,ja)-xred01(:,ka)+dble(is(:))-dble(js(:)) 1538 rsqjk = dot_product(rjk(:),matmul(rmet,rjk(:))) 1539 rrjk = dsqrt(rsqjk) 1540 rki(:) = xred01(:,ka)-xred01(:,ia)+dble(js(:)) 1541 rsqki = dot_product(rki(:),matmul(rmet,rki(:))) 1542 rrki = dsqrt(rsqki) 1543 if (rsqij>=tol16.and.rsqjk>=tol16.and.rsqki>=tol16) then 1544 rmean = (rrij*rrjk*rrki)**third 1545 if (rrij>rcut9.or.rrjk>rcut9.or.rrki>rcut9) cycle 1546 sfact9=vdw_s6 1547 if (ia==ja.and.ja==ka) sfact9 = sixth*sfact9 1548 if (ia==ja.and.ja/=ka) sfact9 = half*sfact9 1549 if (ia/=ja.and.ja==ka) sfact9 = half*sfact9 1550 rijk = one/(rrij*rrjk*rrki) 1551 dmp9 = six*(rmean*vdw_sr9*r0ijk(itypat,jtypat,ktypat))**(-alpha8) 1552 fdmp9 = one/(one+dmp9) 1553 cosa = half*rrjk*(rsqij+rsqki-rsqjk)*rijk 1554 cosb = half*rrki*(rsqij+rsqjk-rsqki)*rijk 1555 cosc = half*rrij*(rsqjk+rsqki-rsqij)*rijk 1556 ang = one+three*cosa*cosb*cosc 1557 temp = sfact9*rijk*rijk*rijk 1558 temp2 = temp*fdmp9*ang 1559 ! Contribution to energy 1560 e_3bt = e_3bt-temp2*vdw_c9(ia,ja,ka) !*temp2 1561 e3bt_ij(ia,ja) = e3bt_ij(ia,ja)-temp2*vdw_c9(ia,ja,ka) 1562 e3bt_jk(ja,ka) = e3bt_jk(ja,ka)-temp2*vdw_c9(ia,ja,ka) 1563 e3bt_ki(ka,ia) = e3bt_ki(ka,ia)-temp2*vdw_c9(ia,ja,ka) 1564 if (need_grad) then 1565 dfdmp = third*alpha8*fdmp9*fdmp9*dmp9 1566 if (ia/=ja.or.need_stress) then 1567 d1_r3drij = -three*rrki*rrjk 1568 dcosa_r3drij = (rrij-two*cosa*rrki)*rrjk 1569 dcosb_r3drij = (rrij-two*cosb*rrjk)*rrki 1570 dcosc_r3drij =-(rsqij+cosc*rrjk*rrki) 1571 dfdmp_drij = dfdmp*rrki*rrjk 1572 d_drij = vdw_c9(ia,ja,ka)*temp*rijk*((d1_r3drij+three*& 1573 & (cosb*cosc*dcosa_r3drij+cosa*cosc*dcosb_r3drij+cosa*cosb*dcosc_r3drij))*& 1574 & fdmp9+dfdmp_drij*ang)*rijk*rrjk*rrki 1575 rcartij=matmul(rprimd,rij) 1576 end if 1577 if (ja/=ka.or.need_stress) then 1578 d1_r3drjk = -three*rrij*rrki 1579 dcosa_r3drjk =-(rsqjk+cosa*rrij*rrki) 1580 dcosb_r3drjk = (rrjk-two*cosb*rrij)*rrki 1581 dcosc_r3drjk = (rrjk-two*cosc*rrki)*rrij 1582 dfdmp_drjk = dfdmp*rrij*rrki 1583 d_drjk = vdw_c9(ia,ja,ka)*temp*rijk*((d1_r3drjk+three*& 1584 & (cosb*cosc*dcosa_r3drjk+cosa*cosc*dcosb_r3drjk+cosa*cosb*dcosc_r3drjk))*& 1585 & fdmp9+dfdmp_drjk*ang)*rijk*rrij*rrki 1586 rcartjk=matmul(rprimd,rjk) 1587 end if 1588 if (ka/=ia.or.need_stress) then 1589 d1_r3drki = -three*rrjk*rrij 1590 dcosa_r3drki = (rrki-two*cosa*rrij)*rrjk 1591 dcosb_r3drki =-(rsqki+cosb*rrij*rrjk) 1592 dcosc_r3drki = (rrki-two*cosc*rrjk)*rrij 1593 dfdmp_drki = dfdmp*rrij*rrjk 1594 d_drki = vdw_c9(ia,ja,ka)*temp*rijk*((d1_r3drki+three*& 1595 & (cosb*cosc*dcosa_r3drki+cosa*cosc*dcosb_r3drki+cosa*cosb*dcosc_r3drki))*& 1596 & fdmp9+dfdmp_drki*ang)*rijk*rrij*rrjk 1597 rcartki=matmul(rprimd,rki) 1598 end if 1599 ! Contribution to gradients wr to atomic displacement 1600 ! (forces) 1601 if (need_forces) then 1602 if (ia/=ja) gredij=d_drij*matmul(transpose(rprimd),rcartij) 1603 if (ja/=ka) gredjk=d_drjk*matmul(transpose(rprimd),rcartjk) 1604 if (ka/=ia) gredki=d_drki*matmul(transpose(rprimd),rcartki) 1605 if (ia/=ja.and.ka/=ia) then 1606 gred_vdw_3bt(:,ia)=gred_vdw_3bt(:,ia)-gredij(:)+gredki(:) 1607 gred_vdw_3bt(:,ja)=gred_vdw_3bt(:,ja)+gredij(:)-gredjk(:) 1608 gred_vdw_3bt(:,ka)=gred_vdw_3bt(:,ka)-gredki(:)+gredjk(:) 1609 else if (ia==ja.and.ia/=ka) then 1610 gred_vdw_3bt(:,ia)=gred_vdw_3bt(:,ia)+gredki(:) 1611 gred_vdw_3bt(:,ja)=gred_vdw_3bt(:,ja)-gredjk(:) 1612 gred_vdw_3bt(:,ka)=gred_vdw_3bt(:,ka)-gredki(:)+gredjk(:) 1613 elseif (ia==ka.and.ia/=ja) then 1614 gred_vdw_3bt(:,ia)=gred_vdw_3bt(:,ia)-gredij(:) 1615 gred_vdw_3bt(:,ja)=gred_vdw_3bt(:,ja)+gredij(:)-gredjk(:) 1616 gred_vdw_3bt(:,ka)=gred_vdw_3bt(:,ka)+gredjk(:) 1617 elseif (ja==ka.and.ia/=ja) then 1618 gred_vdw_3bt(:,ia)=gred_vdw_3bt(:,ia)-gredij(:)+gredki(:) 1619 gred_vdw_3bt(:,ja)=gred_vdw_3bt(:,ja)+gredij(:) 1620 gred_vdw_3bt(:,ka)=gred_vdw_3bt(:,ka)-gredki(:) 1621 end if 1622 end if 1623 ! Contribution to stress tensor 1624 if (need_stress) then 1625 vecij(1:3)=rcartij(1:3)*rcartij(1:3); vecij(4)=rcartij(2)*rcartij(3) 1626 vecij(5)=rcartij(1)*rcartij(3); vecij(6)=rcartij(1)*rcartij(2) 1627 vecjk(1:3)=rcartjk(1:3)*rcartjk(1:3); vecjk(4)=rcartjk(2)*rcartjk(3) 1628 vecjk(5)=rcartjk(1)*rcartjk(3); vecjk(6)=rcartjk(1)*rcartjk(2) 1629 vecki(1:3)=rcartki(1:3)*rcartki(1:3); vecki(4)=rcartki(2)*rcartki(3) 1630 vecki(5)=rcartki(1)*rcartki(3); vecki(6)=rcartki(1)*rcartki(2) 1631 str_3bt(:)=str_3bt(:)-d_drij*vecij(:)-d_drjk*vecjk(:)-d_drki*vecki(:) !-str_3bt_dcn9(:) 1632 end if 1633 end if ! Optimization 1634 end if ! Tolerance 1635 end do ! Loop over atom k 1636 end do ! Loop over atom j 1637 end do ! Loop over atom i 1638 end do ! j3 1639 end do ! j2 1640 end do ! j1 1641 end do 1642 end do 1643 end do 1644 if (need_forces) then 1645 do ia=1,natom 1646 do ja=1,natom 1647 do la=1,natom 1648 gred_vdw_3bt(:,la) = gred_vdw_3bt(:,la)+e3bt_ij(ia,ja)*(dc9ijri(ia,ja)*dcn(:,ia,la)+dc9ijrj(ia,ja)*dcn(:,ja,la)) 1649 gred_vdw_3bt(:,la) = gred_vdw_3bt(:,la)+e3bt_jk(ia,ja)*(dc9ijri(ia,ja)*dcn(:,ia,la)+dc9ijrj(ia,ja)*dcn(:,ja,la)) 1650 gred_vdw_3bt(:,la) = gred_vdw_3bt(:,la)+e3bt_ki(ia,ja)*(dc9ijri(ia,ja)*dcn(:,ia,la)+dc9ijrj(ia,ja)*dcn(:,ja,la)) 1651 end do 1652 end do 1653 end do 1654 elseif (need_stress) then 1655 do ia=1,natom 1656 do ja=1,natom 1657 str_3bt(:) = str_3bt(:)+e3bt_ij(ia,ja)*(dc9ijri(ia,ja)*str_dcn(:,ia)+dc9ijrj(ia,ja)*str_dcn(:,ja)) 1658 str_3bt(:) = str_3bt(:)+e3bt_jk(ia,ja)*(dc9ijri(ia,ja)*str_dcn(:,ia)+dc9ijrj(ia,ja)*str_dcn(:,ja)) 1659 str_3bt(:) = str_3bt(:)+e3bt_ki(ia,ja)*(dc9ijri(ia,ja)*str_dcn(:,ia)+dc9ijrj(ia,ja)*str_dcn(:,ja)) 1660 end do 1661 end do 1662 end if 1663 e_vdw_dftd3 = e_vdw_dftd3+e_3bt 1664 if (need_forces) gred_vdw_dftd3= gred_vdw_dftd3+gred_vdw_3bt 1665 if (need_stress) str_vdw_dftd3 = str_vdw_dftd3+str_3bt 1666 ABI_FREE(dc9ijri) 1667 ABI_FREE(dc9ijrj) 1668 ABI_FREE(e3bt_ij) 1669 ABI_FREE(e3bt_jk) 1670 ABI_FREE(e3bt_ki) 1671 ABI_FREE(vdw_c9) 1672 ABI_FREE(r0ijk) 1673 ABI_FREE(gred_vdw_3bt) 1674 end if 1675 if (need_stress) str_vdw_dftd3=str_vdw_dftd3/ucvol 1676 1677 !Printing 1678 if (prtvol>0) then 1679 write(msg,'(2a)') ch10,& 1680 & ' --------------------------------------------------------------' 1681 call wrtout(std_out,msg,'COLL') 1682 if (vdw_xc==6) then 1683 write(msg,'(3a)') & 1684 & ' Van der Waals DFT-D3 semi-empirical dispersion potential as',ch10,& 1685 & ' proposed by Grimme et al., J. Chem. Phys. 132, 154104 (2010)' ! [[cite:Grimme2010]] 1686 call wrtout(std_out,msg,'COLL') 1687 elseif (vdw_xc==7) then 1688 write(msg,'(5a)') & 1689 & ' Van der Waals DFT-D3 semi-empirical dispersion potential ' ,ch10,& 1690 & ' with Becke-Jonhson (BJ) refined by Grimme et al. J. ',ch10,& 1691 & ' Comput. Chem. 32, 1456 (2011) ' ! [[cite:Grimme2011]] 1692 call wrtout(std_out,msg,'COLL') 1693 end if 1694 if (natom<5) then 1695 write(msg,'(3a)')& 1696 & ' Pair C6 (a.u.) C8 (a.u.) R0 (Ang) ',ch10,& 1697 & ' ---------------------------------------------------------------' 1698 call wrtout(std_out,msg,'COLL') 1699 do ia=1,natom 1700 do ja=1,ia 1701 itypat = typat(ia) ; jtypat = typat(ja) 1702 call atomdata_from_znucl(atom1,znucl(itypat)) 1703 call atomdata_from_znucl(atom2,znucl(jtypat)) 1704 write(msg,'(4X,2a,i2,3a,i2,1a,1X,es12.4,4X,es12.4,4X,es12.4,1X)') & 1705 atom1%symbol,'(',ia,')-',atom2%symbol,'(',ja,')', & 1706 vdw_c6(ia,ja), vdw_c8(ia,ja),& 1707 vdw_r0(itypat,jtypat) 1708 call wrtout(std_out,msg,'COLL') 1709 end do 1710 end do 1711 end if 1712 write(msg, '(3a,f6.3,a,f6.3)') & 1713 & ' ---------------------------------------------------------------',ch10,& 1714 & ' Scaling factors: s6 = ', vdw_s6,', s8 = ',vdw_s8 1715 call wrtout(std_out,msg,'COLL') 1716 if (vdw_xc==6) then 1717 write(msg,'(a,f6.3,a,f6.3)') & 1718 & ' Damping parameters: sr6 = ', vdw_sr6,', sr8 = ',vdw_sr8 1719 call wrtout(std_out,msg,'COLL') 1720 elseif (vdw_xc==7) then 1721 write(msg,'(a,f6.3,a,f6.3)') & 1722 & ' Damping parameters: a1 = ', vdw_a1, ', a2 = ', vdw_a2 1723 call wrtout(std_out,msg,'COLL') 1724 end if 1725 write(msg,'(a,es12.5,3a,i14,2a,es12.5,1a)') & 1726 & ' Cut-off radius = ',rcut,' Bohr',ch10,& 1727 & ' Number of pairs contributing = ',npairs,ch10,& 1728 & ' DFT-D3 (no 3-body) energy contribution = ',e_vdw_dftd3-e_3bt,' Ha' 1729 call wrtout(std_out,msg,'COLL') 1730 if (bol_3bt) then 1731 write(msg,'(6a,i5,2a,es20.11,3a,es20.11,1a)')ch10,& 1732 & ' ---------------------------------------------------------------',ch10,& 1733 & ' 3-Body Term Contribution:', ch10,& 1734 & ' Number of shells considered = ', nshell, ch10,& 1735 & ' Additional 3-body contribution = ', e_3bt, ' Ha',ch10,& 1736 & ' Total E (2-body and 3-body) = ', e_vdw_dftd3, 'Ha' 1737 call wrtout(std_out,msg,'COLL') 1738 end if 1739 write(msg,'(2a)')& 1740 & ' ----------------------------------------------------------------',ch10 1741 call wrtout(std_out,msg,'COLL') 1742 end if 1743 ABI_FREE(ivdw) 1744 ABI_FREE(xred01) 1745 ABI_FREE(vdw_r0) 1746 ABI_FREE(fe_no_c) 1747 ABI_FREE(e_no_c) 1748 ABI_FREE(e_alpha1) 1749 ABI_FREE(e_alpha2) 1750 ABI_FREE(e_alpha3) 1751 ABI_FREE(e_alpha4) 1752 ABI_FREE(grad_no_cij) 1753 ABI_FREE(fgrad_no_c) 1754 ABI_FREE(cfgrad_no_c) 1755 ABI_FREE(vdw_c6) 1756 ABI_FREE(vdw_c8) 1757 ABI_FREE(dc6ri) 1758 ABI_FREE(dc6rj) 1759 ABI_FREE(d2c6ri) 1760 ABI_FREE(d2c6rj) 1761 ABI_FREE(d2c6rirj) 1762 ABI_FREE(cn) 1763 ABI_FREE(d2cn_iii) 1764 ABI_FREE(d2cn_jji) 1765 ABI_FREE(d2cn_iji) 1766 ABI_FREE(d2cn_jii) 1767 ABI_FREE(d2cn_tmp) 1768 ABI_FREE(dcn) 1769 ABI_FREE(fdcn) 1770 ABI_FREE(cfdcn) 1771 ABI_FREE(str_dcn) 1772 ABI_FREE(elt_cn) 1773 ABI_FREE(str_no_c) 1774 ABI_FREE(str_alpha1) 1775 ABI_FREE(str_alpha2) 1776 ABI_FREE(dcn_cart) 1777 DBG_EXIT("COLL") 1778 1779 contains 1780 !! ***