TABLE OF CONTENTS


ABINIT/m_vdw_dftd3 [ Modules ]

[ Top ] [ Modules ]

NAME

  m_vdw_dftd3

FUNCTION

COPYRIGHT

  Copyright (C) 2015-2018 ABINIT group (BVT)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_vdw_dftd3
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32  use m_atomdata
33 
34  use m_special_funcs,  only : abi_derfc
35  use m_geometry,       only : metric
36  use m_vdw_dftd3_data, only : vdw_dftd3_data
37 
38  implicit none
39 
40  private

ABINIT/vdw_dftd3 [ Functions ]

[ Top ] [ 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
  [fred_vdw_dftd3(3,natom)]=contribution to gradient w.r.to atomic displ.
  from DFT-D3 dispersion potential
  [str_vdw_dftd3(6)]=contribution to stress tensor from DFT-D3 dispersion potential
  [dyn_vdw_dftd3(2,3,natom,3,natom)]= contribution to the interatomic force
  constants (in reciprocal space) at given input q-vector
  from DFT-D3 dispersion potential

NOTES

  Ref.:
  DFT-D3: S. Grimme, J. Antony, S. Ehrlich, and H. Krieg
  A consistent and accurate ab initio parametrization of density functional
  dispersion correction (DFT-D) for the 94 elements H-Pu
  J. Chem. Phys. 132, 154104 (2010) [[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]]

PARENTS

      respfn,setvtr,stress

CHILDREN

SOURCE

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