Branch data Line data Source code
1 : : !!****m* ABINIT/m_vdw_dftd2
2 : : !! NAME
3 : : !! m_vdw_dftd2
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2008-2026 ABINIT group (MT)
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_dftd2
22 : :
23 : : use defs_basis
24 : : use m_abicore
25 : : use m_errors
26 : : use m_atomdata
27 : :
28 : : use m_geometry, only : metric
29 : :
30 : : implicit none
31 : :
32 : : private
33 : : !!***
34 : :
35 : : public :: vdw_dftd2
36 : : !!***
37 : :
38 : : contains
39 : : !!***
40 : :
41 : : !!****f* ABINIT/vdw_dftd2
42 : : !!
43 : : !! NAME
44 : : !! vdw_dftd2
45 : : !!
46 : : !! FUNCTION
47 : : !! Compute energy and derivatives with respect to dimensionless
48 : : !! reduced atom coordinates due to Van der Waals interaction.
49 : : !! The formalism here follows the DFT-D2 approach of Grimme
50 : : !! which consists in adding a semi-empirical dispersion potential
51 : : !! (pair-wise force field) to the conventional Kohn-Sham DFT energy.
52 : : !!
53 : : !! INPUTS
54 : : !! ixc=choice of exchange-correlation functional
55 : : !! natom=number of atoms
56 : : !! ntypat=number of atom types
57 : : !! prtvol=printing volume (if >0, print computation parameters)
58 : : !! typat(natom)=type integer for each atom in cell
59 : : !! rprimd(3,3)=real space primitive translations
60 : : !! vdw_tol=tolerance use to converge the potential (a pair of atoms is included
61 : : !! in potential if its contribution is larger than vdw_tol)
62 : : !! vdw_tol<0 takes default value (10^-10)
63 : : !! xred(3,natom)=reduced atomic coordinates
64 : : !! znucl(ntypat)=atomic number of atom type
65 : : !! === optional inputs ===
66 : : !! [qphon(3)]=wavevector of the phonon;
67 : : !! used only for dynamical matrix computation
68 : : !!
69 : : !! OUTPUT
70 : : !! e_vdw_dftd2=contribution to energy from DFT-D2 dispersion potential
71 : : !! === optional outputs ===
72 : : !! [dyn_vdw_dftd2(2,3,natom,3,natom)]=contribution to dynamical matrix from DFT-D2 dispersion potential
73 : : !! [elt_vdw_dftd2(6+3*natom,6)]=contribution to elastic tensor and internal strains from DFT-D2 disp. pot.
74 : : !! [gred_vdw_dftd2(3,natom)]=contribution to gradients wrt nuclear positions from DFT-D2 dispersion potential
75 : : !! [str_vdw_dftd2(6)]=contribution to stress tensor from DFT-D2 dispersion potential
76 : : !!
77 : : !! NOTES
78 : : !! Ref.: S. Grimme, Semiempirical GGA-type density functional
79 : : !! constructed with a long-range dispersion correction,
80 : : !! J. Comp. Chem. 27, 1787 (2006) [[cite:Grimme2006]]
81 : : !!
82 : : !! SOURCE
83 : :
84 : 40 : subroutine vdw_dftd2(e_vdw_dftd2,ixc,natom,ntypat,prtvol,typat,rprimd,vdw_tol,xred,znucl,&
85 : 10 : & dyn_vdw_dftd2,elt_vdw_dftd2,gred_vdw_dftd2,str_vdw_dftd2,qphon) ! Optionals
86 : :
87 : : !Arguments ------------------------------------
88 : : !scalars
89 : : integer,intent(in) :: ixc,natom,ntypat,prtvol
90 : : real(dp),intent(in) :: vdw_tol
91 : : real(dp),intent(out) :: e_vdw_dftd2
92 : : !arrays
93 : : integer,intent(in) :: typat(natom)
94 : : real(dp),intent(in) :: rprimd(3,3),xred(3,natom),znucl(ntypat)
95 : : real(dp),intent(in),optional :: qphon(3)
96 : : real(dp),intent(out),optional :: dyn_vdw_dftd2(2,3,natom,3,natom)
97 : : real(dp),intent(out),optional :: elt_vdw_dftd2(6+3*natom,6)
98 : : real(dp),intent(out),optional :: gred_vdw_dftd2(3,natom)
99 : : real(dp),intent(out),optional :: str_vdw_dftd2(6)
100 : :
101 : : !Local variables-------------------------------
102 : : !scalars
103 : : integer,parameter :: vdw_nspecies=55
104 : : integer :: ia,ia1,ii,is1,is2,is3,itypat,ja,ja1,jj,jtypat,kk,ll,mu,npairs,nshell,nu
105 : : logical :: need_dynmat,need_elast,need_forces,need_intstr,need_stress
106 : : logical :: need_gradient,need_gradient2,newshell,qeq0=.true.
107 : : real(dp),parameter :: e_conv=(10/Bohr_Ang)**6/Ha_J/Avogadro ! 1 J.nm^6.mol^-1 in Ha.Bohr^6
108 : : real(dp),parameter :: vdw_d=20._dp,vdw_tol_default=tol10
109 : : real(dp),parameter :: vdw_s_pbe=0.75_dp, vdw_s_blyp=1.2_dp, vdw_s_b3lyp=1.05_dp
110 : : real(dp),parameter :: vdw_s_bp86=1.05_dp, vdw_s_tpss=1.0_dp, vdw_s_b97d=1.25_dp
111 : : real(dp) :: c6,c6r6,ex,fr,gr,gr2,grad,grad2,ph,ph1r,ph1i
112 : : real(dp) :: r0,r1,r2,r3,rcut,rcut2,rsq,rr,sfact,ucvol,vdw_s
113 : : character(len=500) :: msg
114 : : type(atomdata_t) :: atom
115 : : !arrays
116 : 20 : integer,allocatable :: ivdw(:)
117 : : integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
118 : : real(dp) :: gmet(3,3),gprimd(3,3),mat(3,3),rcart(3),rmet(3,3),vec(3)
119 : 20 : real(dp),allocatable :: vdw_c6(:,:),vdw_r0(:,:),xred01(:,:)
120 : : real(dp),parameter :: vdw_c6_dftd2(vdw_nspecies)= &
121 : : & (/ 0.14, 0.08, 1.61, 1.61, 3.13, 1.75, 1.23, 0.70, 0.75, 0.63,&
122 : : & 5.71, 5.71,10.79, 9.23, 7.84, 5.57, 5.07, 4.61,10.80,10.80,&
123 : : & 10.80,10.80,10.80,10.80,10.80,10.80,10.80,10.80,10.80,10.80,&
124 : : & 16.99,17.10,16.37,12.64,12.47,12.01,24.67,24.67,24.67,24.67,&
125 : : & 24.67,24.67,24.67,24.67,24.67,24.67,24.67,24.67,37.32,38.71,&
126 : : & 38.44,31.74,31.50,29.99, 0.00/)
127 : : real(dp),parameter :: vdw_r0_dftd2(vdw_nspecies)= &
128 : : & (/1.001,1.012,0.825,1.408,1.485,1.452,1.397,1.342,1.287,1.243,&
129 : : & 1.144,1.364,1.639,1.716,1.705,1.683,1.639,1.595,1.485,1.474,&
130 : : & 1.562,1.562,1.562,1.562,1.562,1.562,1.562,1.562,1.562,1.562,&
131 : : & 1.650,1.727,1.760,1.771,1.749,1.727,1.628,1.606,1.639,1.639,&
132 : : & 1.639,1.639,1.639,1.639,1.639,1.639,1.639,1.639,1.672,1.804,&
133 : : & 1.881,1.892,1.892,1.881,1.000/)
134 : : character(len=2),parameter :: vdw_symb(vdw_nspecies)= &
135 : : & (/' H','He','Li','Be',' B',' C',' N',' O',' F','Ne',&
136 : : & 'Na','Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca',&
137 : : & 'Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn',&
138 : : & 'Ga','Ge','As','Se','Br','Kr','Rb','Sr',' Y','Zr',&
139 : : & 'Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn',&
140 : : & 'Sb','Te',' I','Xe','no'/)
141 : :
142 : : ! *************************************************************************
143 : :
144 : : DBG_ENTER("COLL")
145 : :
146 : : !Extract options
147 : 20 : need_forces=present(gred_vdw_dftd2)
148 : 20 : need_stress=present(str_vdw_dftd2)
149 : 20 : need_dynmat=present(dyn_vdw_dftd2)
150 : 20 : need_elast=present(elt_vdw_dftd2)
151 : 20 : need_intstr=present(elt_vdw_dftd2)
152 : 20 : need_gradient=(need_forces.or.need_stress)
153 [ + - ]: 20 : need_gradient2=(need_dynmat.or.need_elast.or.need_intstr)
154 [ - + ]: 20 : if (need_dynmat) then
155 [ # # ]: 0 : if (.not.present(qphon)) then
156 : 0 : msg='Dynamical matrix required without a q-vector'
157 : 0 : ABI_BUG(msg)
158 : : end if
159 : 0 : qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
160 : : end if
161 : :
162 : : !Identify type(s) of atoms
163 [ - + + - : 60 : ABI_MALLOC(ivdw,(ntypat))
+ - + - -
+ - + ]
164 [ + + ]: 40 : do itypat=1,ntypat
165 : 20 : ivdw(itypat)=-1;jtypat=0
166 : 20 : call atomdata_from_znucl(atom,znucl(itypat))
167 [ + + - + ]: 380 : do while ((ivdw(itypat)<0).and.(jtypat<vdw_nspecies))
168 [ + + ]: 380 : jtypat=jtypat+1;if (vdw_symb(jtypat)==atom%symbol) ivdw(itypat)=jtypat
169 : : end do
170 [ - + ]: 60 : if (ivdw(itypat)<0) then
171 : : write(msg,'(3a)') &
172 : 0 : & 'Van der Waals DFT-D2 correction not available for atom type: ',atom%symbol,' !'
173 : 0 : ABI_ERROR(msg)
174 : : end if
175 : : end do
176 : :
177 : : !Select DFT-D2 VdW parameters according to system data
178 : 20 : vdw_s=e_conv
179 [ - - - - : 20 : if (ixc==11.or.ixc==-101130.or.ixc==-130101) then
- - + ]
180 : : vdw_s=vdw_s*vdw_s_pbe
181 : : else if (ixc==18.or.ixc==-106131.or.ixc==-131106) then
182 : : vdw_s=vdw_s*vdw_s_blyp
183 : : else if (ixc==19.or.ixc==-106132.or.ixc==-132106) then
184 : : vdw_s=vdw_s*vdw_s_bp86
185 : : else if (ixc==-202231.or.ixc==-231202) then
186 : : vdw_s=vdw_s*vdw_s_tpss
187 : : else
188 : 0 : write(msg,'(a,i8,a)')' Van der Waals DFT-D2 correction not compatible with ixc=',ixc,' !'
189 : 0 : ABI_ERROR(msg)
190 : : end if
191 [ - + + - : 80 : ABI_MALLOC(vdw_c6,(ntypat,ntypat))
- + + - +
- + - - +
- + ]
192 [ - + + - : 60 : ABI_MALLOC(vdw_r0,(ntypat,ntypat))
- + + - -
+ - + ]
193 [ + + ]: 40 : do itypat=1,ntypat
194 [ + + ]: 60 : do jtypat=1,ntypat
195 : 20 : vdw_c6(itypat,jtypat)=sqrt(vdw_c6_dftd2(ivdw(itypat))*vdw_c6_dftd2(ivdw(jtypat)))
196 : 40 : vdw_r0(itypat,jtypat)=(vdw_r0_dftd2(ivdw(itypat))+vdw_r0_dftd2(ivdw(jtypat)))/Bohr_Ang
197 : : end do
198 : : end do
199 : :
200 : : !Computation of cut-off radius according to tolerance
201 : : !We take: r_cut=(s6*max(C6)/tol)**(1/6) and rcut<=75 bohr
202 [ - + ]: 20 : if (vdw_tol<zero) then
203 [ # # # # : 0 : rcut=(vdw_s/vdw_tol_default*maxval(vdw_c6))**sixth
# # # # #
# # # ]
204 : : else
205 [ + + + + : 60 : rcut=(vdw_s/vdw_tol*maxval(vdw_c6))**sixth
- + + - -
+ - - ]
206 : : end if
207 : : !rcut=min(rcut,100._dp)
208 : 20 : rcut2=rcut*rcut
209 : :
210 : : !Retrieve cell geometry data
211 : 20 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
212 : :
213 : : !Map reduced coordinates into [0,1]
214 [ - + + - : 60 : ABI_MALLOC(xred01,(3,natom))
+ - + - -
+ - + ]
215 [ + + ]: 40 : do ia=1,natom
216 : 20 : xred01(1,ia)=xred(1,ia)-aint(xred(1,ia))+half-sign(half,xred(1,ia))
217 : 20 : xred01(2,ia)=xred(2,ia)-aint(xred(2,ia))+half-sign(half,xred(2,ia))
218 : 40 : xred01(3,ia)=xred(3,ia)-aint(xred(3,ia))+half-sign(half,xred(3,ia))
219 : : end do
220 : :
221 : : !Set accumulated quantities to zero
222 : 20 : npairs=0
223 : 20 : e_vdw_dftd2=zero
224 [ + + + + : 60 : if (need_forces) gred_vdw_dftd2=zero
+ + ]
225 [ + + ]: 20 : if (need_stress) str_vdw_dftd2=zero
226 [ - + - - : 20 : if (need_dynmat) dyn_vdw_dftd2=zero
- - - - -
- - - ]
227 [ - + - - : 20 : if (need_elast) elt_vdw_dftd2(1:6,1:6)=zero
- - ]
228 [ # # # # ]: 0 : if (need_intstr) elt_vdw_dftd2(7:6+3*natom,1:6)=zero
229 : :
230 : : !Loop over shells of cell replicas
231 : : nshell=0
232 : : do
233 : 140 : newshell=.false.;nshell=nshell+1
234 : :
235 : : ! Loop over cell replicas in the shell
236 : : ! ns1=1+int(rcut*sqrt(SUM(gprimd(:,1)**2))
237 : : ! ns2=1+int(rcut*sqrt(SUM(gprimd(:,2)**2))
238 : : ! ns3=1+int(rcut*sqrt(SUM(gprimd(:,3)**2))
239 [ + + ]: 1400 : do is3=-nshell,nshell
240 [ + + ]: 14980 : do is2=-nshell,nshell
241 [ + + ]: 177380 : do is1=-nshell,nshell
242 : : if (nshell==1.or. &
243 [ + + + + : 176120 : & abs(is3)==nshell.or.abs(is2)==nshell.or.abs(is1)==nshell) then
+ + + + ]
244 : :
245 : : ! Phase for dynamical matrix
246 [ - + ]: 67500 : if (need_dynmat) then
247 : 0 : ph1r=one;ph1i=zero !ph1=exp(-iqL)
248 [ # # ]: 0 : if (.not.qeq0) then
249 : 0 : ph=-two_pi*(qphon(1)*is1+qphon(2)*is2+qphon(3)*is3)
250 : 0 : ph1r=cos(ph);ph1i=sin(ph)
251 : : end if
252 : : end if
253 : :
254 : : ! Loops over atoms a and b
255 [ + + ]: 135000 : do ia=1,natom
256 : 67500 : itypat=typat(ia)
257 [ + + ]: 297540 : do ja=1,ia
258 : 67500 : jtypat=typat(ja)
259 : 67500 : r1=xred01(1,ia)-xred01(1,ja)-dble(is1)
260 : 67500 : r2=xred01(2,ia)-xred01(2,ja)-dble(is2)
261 : 67500 : r3=xred01(3,ia)-xred01(3,ja)-dble(is3)
262 : : rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3 &
263 : 67500 : & +two*(rmet(2,1)*r2*r1+rmet(3,2)*r3*r2+rmet(3,1)*r1*r3)
264 : :
265 : : ! Select atomic pairs (a,b) and avoid atom_a=atom_b
266 [ + + + + ]: 135000 : if (rsq>=tol16.and.rsq<=rcut2) then
267 : :
268 : : ! Data for the selected pair
269 : 15416 : npairs=npairs+1;newshell=.true.
270 [ + - ]: 15416 : sfact=vdw_s;if (ia==ja) sfact=half*sfact
271 : 15416 : rr=sqrt(rsq)
272 : 15416 : c6=vdw_c6(itypat,jtypat)
273 : 15416 : r0=vdw_r0(itypat,jtypat)
274 : :
275 : : ! Computation of pair-wise potential
276 : 15416 : ex=exp(-vdw_d*(rr/r0-one))
277 : 15416 : fr=one/(one+ex)
278 : 15416 : c6r6=c6/rr**6
279 : :
280 : : ! Contribution to energy
281 : 15416 : e_vdw_dftd2=e_vdw_dftd2-sfact*fr*c6r6
282 : :
283 [ + - ]: 15416 : if (need_gradient.or.need_gradient2) then
284 : 15416 : gr=(vdw_d/r0)*(fr**2)*ex
285 : 15416 : grad=-sfact*(gr-six*fr/rr)*c6r6/rr
286 : 15416 : rcart(1)=rprimd(1,1)*r1+rprimd(1,2)*r2+rprimd(1,3)*r3
287 : 15416 : rcart(2)=rprimd(2,1)*r1+rprimd(2,2)*r2+rprimd(2,3)*r3
288 : 15416 : rcart(3)=rprimd(3,1)*r1+rprimd(3,2)*r2+rprimd(3,3)*r3
289 : :
290 : : ! Contribution to gradients wrt nuclear positions
291 [ - + ]: 15416 : if (need_forces.and.ia/=ja) then
292 [ # # ]: 0 : vec(1:3)=grad*rcart(1:3)
293 [ # # ]: 0 : gred_vdw_dftd2(1:3,ia)=gred_vdw_dftd2(1:3,ia)+vec(1:3)
294 [ # # ]: 0 : gred_vdw_dftd2(1:3,ja)=gred_vdw_dftd2(1:3,ja)-vec(1:3)
295 : : end if
296 : :
297 : : ! Contribution to stress tensor
298 [ + + ]: 15416 : if (need_stress) then
299 [ + + ]: 53956 : do mu=1,6
300 : 46248 : ii=alpha(mu);jj=beta(mu)
301 : 53956 : str_vdw_dftd2(mu)=str_vdw_dftd2(mu)+grad*rcart(ii)*rcart(jj)
302 : : end do
303 : : end if
304 : :
305 [ - + ]: 15416 : if (need_gradient2) then
306 : 0 : gr2=(vdw_d/r0)*gr*(2*fr*ex-one)
307 : 0 : grad2=-sfact*(gr2-13._dp*gr/rr+48._dp*fr/rr**2)*c6r6/rr**2
308 : :
309 : : ! Contribution to dynamical matrix (phase factors are subtle!)
310 [ # # ]: 0 : if (need_dynmat) then
311 [ # # ]: 0 : mat(1:3,1)=grad2*rcart(1:3)*rcart(1) ; mat(1,1)=mat(1,1)+grad
312 [ # # ]: 0 : mat(1:3,2)=grad2*rcart(1:3)*rcart(2) ; mat(2,2)=mat(2,2)+grad
313 [ # # ]: 0 : mat(1:3,3)=grad2*rcart(1:3)*rcart(3) ; mat(3,3)=mat(3,3)+grad
314 [ # # ]: 0 : if (ia/=ja) then
315 [ # # ]: 0 : do ii=1,3
316 [ # # ]: 0 : dyn_vdw_dftd2(1,1:3,ia,ii,ia)=dyn_vdw_dftd2(1,1:3,ia,ii,ia)+mat(1:3,ii)
317 [ # # ]: 0 : dyn_vdw_dftd2(1,1:3,ja,ii,ja)=dyn_vdw_dftd2(1,1:3,ja,ii,ja)+mat(1:3,ii)
318 [ # # ]: 0 : dyn_vdw_dftd2(1,1:3,ia,ii,ja)=dyn_vdw_dftd2(1,1:3,ia,ii,ja)-mat(1:3,ii)*ph1r
319 [ # # ]: 0 : dyn_vdw_dftd2(2,1:3,ia,ii,ja)=dyn_vdw_dftd2(2,1:3,ia,ii,ja)-mat(1:3,ii)*ph1i
320 [ # # ]: 0 : dyn_vdw_dftd2(1,1:3,ja,ii,ia)=dyn_vdw_dftd2(1,1:3,ja,ii,ia)-mat(1:3,ii)*ph1r
321 [ # # ]: 0 : dyn_vdw_dftd2(2,1:3,ja,ii,ia)=dyn_vdw_dftd2(2,1:3,ja,ii,ia)+mat(1:3,ii)*ph1i
322 : : end do
323 [ # # ]: 0 : else if (.not.qeq0) then
324 [ # # ]: 0 : do ii=1,3
325 : : dyn_vdw_dftd2(1,1:3,ia,ii,ia)=dyn_vdw_dftd2(1,1:3,ia,ii,ia) &
326 [ # # ]: 0 : & +two*mat(1:3,ii)*(one-ph1r)
327 : : end do
328 : : end if
329 : : end if
330 : :
331 : : ! Contribution to elastic tensor
332 [ # # ]: 0 : if (need_elast) then
333 [ # # ]: 0 : do mu=1,6
334 : 0 : ii=alpha(mu);jj=beta(mu)
335 [ # # ]: 0 : do nu=1,6
336 : 0 : kk=alpha(nu);ll=beta(nu)
337 : : elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) &
338 : 0 : & +grad2*rcart(ii)*rcart(jj)*rcart(kk)*rcart(ll)
339 [ # # ]: 0 : if (ii==kk) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) &
340 : 0 : & +half*grad*rcart(jj)*rcart(ll)
341 [ # # ]: 0 : if (ii==ll) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) &
342 : 0 : & +half*grad*rcart(jj)*rcart(kk)
343 [ # # ]: 0 : if (jj==kk) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) &
344 : 0 : & +half*grad*rcart(ii)*rcart(ll)
345 [ # # ]: 0 : if (jj==ll) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) &
346 : 0 : & +half*grad*rcart(ii)*rcart(kk)
347 : : end do
348 : : end do
349 : : end if
350 : :
351 : : ! Contribution to internal strains
352 [ # # ]: 0 : if (need_intstr.and.ia/=ja) then
353 : 0 : ia1=6+3*(ia-1);ja1=6+3*(ja-1)
354 [ # # ]: 0 : do mu=1,6
355 : 0 : ii=alpha(mu);jj=beta(mu)
356 [ # # ]: 0 : vec(1:3)=grad2*rcart(ii)*rcart(jj)*rcart(1:3)
357 : 0 : vec(ii)=vec(ii)+half*grad*rcart(jj)
358 : 0 : vec(jj)=vec(jj)+half*grad*rcart(ii)
359 [ # # ]: 0 : elt_vdw_dftd2(ia1+1:ia1+3,mu)=elt_vdw_dftd2(ia1+1:ia1+3,mu)+vec(1:3)
360 [ # # ]: 0 : elt_vdw_dftd2(ja1+1:ja1+3,mu)=elt_vdw_dftd2(ja1+1:ja1+3,mu)-vec(1:3)
361 : : end do
362 : : end if
363 : :
364 : : end if ! Computation of 2nd gradient
365 : : end if ! Computation of gradient
366 : : end if ! Pairs selection
367 : : end do ! Loop over atom b
368 : : end do ! Loop over atom a
369 : : end if ! Triple loop over cell replicas in shell
370 : : end do
371 : : end do
372 : : end do
373 [ + + ]: 140 : if(.not.newshell) exit ! Check if new shell must be calculated
374 : : end do ! Loop over shells
375 : :
376 : : !Gradients: convert them from cartesian to reduced coordinates
377 [ + + ]: 20 : if (need_forces) then
378 [ + + ]: 20 : do ia=1,natom
379 : 20 : call grad_cart2red(gred_vdw_dftd2(:,ia))
380 : : end do
381 : : end if
382 [ - + ]: 20 : if (need_dynmat) then
383 [ # # ]: 0 : do ja=1,natom
384 [ # # ]: 0 : do ia=1,natom
385 [ # # # # ]: 0 : do kk=1,merge(2,1,qeq0)
386 [ # # ]: 0 : do ii=1,3
387 [ # # ]: 0 : vec(1:3)=dyn_vdw_dftd2(kk,1:3,ia,ii,ja)
388 : 0 : call grad_cart2red(vec)
389 [ # # ]: 0 : dyn_vdw_dftd2(kk,1:3,ia,ii,ja)=vec(1:3)
390 : : end do
391 [ # # ]: 0 : do ii=1,3
392 [ # # ]: 0 : vec(1:3)=dyn_vdw_dftd2(kk,ii,ia,1:3,ja)
393 : 0 : call grad_cart2red(vec)
394 [ # # ]: 0 : dyn_vdw_dftd2(kk,ii,ia,1:3,ja)=vec(1:3)
395 : : end do
396 : : end do
397 : : end do
398 : : end do
399 : : end if
400 [ - + ]: 20 : if (need_intstr) then
401 [ # # ]: 0 : do mu=1,6
402 : : ia1=6
403 [ # # ]: 0 : do ia=1,natom
404 : 0 : call grad_cart2red(elt_vdw_dftd2(ia1+1:ia1+3,mu))
405 : 0 : ia1=ia1+3
406 : : end do
407 : : end do
408 : : end if
409 : :
410 : : !DEBUG
411 : : !write(77,*) "---------------"
412 : : !write(77,*) "E=",e_vdw_dftd2
413 : : !if (need_forces) then
414 : : ! do ia=1,natom
415 : : ! write(77,*) "F=",ia,gred_vdw_dftd2(:,ia)
416 : : ! end do
417 : : !end if
418 : : !if (need_stress) write(77,*) "S=",str_vdw_dftd2(:)
419 : : !if (need_dynmat) then
420 : : ! do ia=1,natom
421 : : ! do ii=1,3
422 : : ! do ja=1,natom
423 : : ! write(77,*) "D=",ia,ii,ja,dyn_vdw_dftd2(:,:,ja,ii,ia)
424 : : ! end do
425 : : ! end do
426 : : ! end do
427 : : !end if
428 : : !if (need_elast) then
429 : : ! do ii=1,6
430 : : ! write(77,*) "e=",ii,elt_vdw_dftd2(1:6,ii)
431 : : ! end do
432 : : !end if
433 : : !if (need_intstr) then
434 : : ! do ii=1,6
435 : : ! do ia=1,natom
436 : : ! write(77,*) "I=",ii,ia,elt_vdw_dftd2(7+3*(ia-1):9+3*(ia-1),ii)
437 : : ! end do
438 : : ! end do
439 : : !end if
440 : : !flush(77)
441 : : !DEBUG
442 : :
443 : : !Stress tensor: divide by volume
444 [ + + + + ]: 80 : if (need_stress) str_vdw_dftd2=str_vdw_dftd2/ucvol
445 : :
446 : : !Printing
447 [ + + ]: 20 : if (prtvol>0) then
448 : 10 : write(msg,'(10a)') ch10,&
449 : 10 : & ' --------------------------------------------------------------',ch10,&
450 : 10 : & ' Van der Waals DFT-D2 semi-empirical dispersion potential added',ch10,&
451 : 10 : & ' with following parameters:',ch10,&
452 : 10 : & ' Specie C6 (J.nm^6.mol^-1) R0 (Ang)',ch10,&
453 : 20 : & ' ------------------------------------'
454 : 10 : call wrtout(std_out,msg,'COLL')
455 [ + + ]: 20 : do itypat=1,ntypat
456 : : write(msg,'(9X,a2,11X,f5.2,8X,f6.3)') &
457 : 10 : & vdw_symb(ivdw(itypat)),vdw_c6_dftd2(ivdw(itypat)),vdw_r0_dftd2(ivdw(itypat))
458 : 20 : call wrtout(std_out,msg,'COLL')
459 : : end do
460 : 10 : write(msg,'(2a,f6.2,2a,f6.2,2a,f6.2,a)') ch10,&
461 : 10 : & ' Scaling factor = ',vdw_s/e_conv,ch10,&
462 : 10 : & ' Damping parameter= ',vdw_d,ch10,&
463 : 20 : & ' Cut-off radius = ',rcut,' bohr'
464 : 10 : call wrtout(std_out,msg,'COLL')
465 : 10 : write(msg,'(2a,i14,2a,es14.5,4a)') ch10,&
466 : 10 : & ' Number of pairs contributing = ',npairs,ch10,&
467 : 10 : & ' DFT-D2 energy contribution = ',e_vdw_dftd2,' Ha',ch10,&
468 : 20 : & ' --------------------------------------------------------------',ch10
469 : 10 : call wrtout(std_out,msg,'COLL')
470 : : end if
471 : :
472 [ - + ]: 20 : ABI_FREE(ivdw)
473 [ - + ]: 20 : ABI_FREE(vdw_c6)
474 [ - + ]: 20 : ABI_FREE(vdw_r0)
475 [ + + - + : 50 : ABI_FREE(xred01)
- + - + ]
476 : :
477 : : DBG_EXIT("COLL")
478 : :
479 : : contains
480 : : !!***
481 : :
482 : : !!****f* vdw_dftd2/grad_cart2red
483 : : !!
484 : : !! NAME
485 : : !! grad_cart2red
486 : : !!
487 : : !! FUNCTION
488 : : !! Convert gradients from cartesian to reduced coordinates
489 : : !!
490 : : !! SOURCE
491 : :
492 : 10 : subroutine grad_cart2red(grad)
493 : :
494 : : !Arguments ------------------------------------
495 : : real(dp),intent(inout) :: grad(3)
496 : : !Local variables-------------------------------
497 : : real(dp) :: tmp(3)
498 : :
499 : : ! *********************************************************************
500 : :
501 : 10 : tmp(1)=rprimd(1,1)*grad(1)+rprimd(2,1)*grad(2)+rprimd(3,1)*grad(3)
502 : 10 : tmp(2)=rprimd(1,2)*grad(1)+rprimd(2,2)*grad(2)+rprimd(3,2)*grad(3)
503 : 10 : tmp(3)=rprimd(1,3)*grad(1)+rprimd(2,3)*grad(2)+rprimd(3,3)*grad(3)
504 : 10 : grad(1:3)=tmp(1:3)
505 : :
506 : 10 : end subroutine grad_cart2red
507 : : !!***
508 : :
509 : : end subroutine vdw_dftd2
510 : : !!***
511 : :
512 : : end module m_vdw_dftd2
513 : : !!***
|
