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