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