Branch data 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 : : !!***
|
