TABLE OF CONTENTS


ABINIT/m_vdw_dftd2 [ Modules ]

[ Top ] [ Modules ]

NAME

  m_vdw_dftd2

FUNCTION

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

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

ABINIT/vdw_dftd2 [ Functions ]

[ Top ] [ Functions ]

NAME

 vdw_dftd2

FUNCTION

 Compute energy and derivatives with respect to dimensionless
 reduced atom coordinates due to Van der Waals interaction.
 The formalism here follows the DFT-D2 approach of Grimme
 which consists in adding a semi-empirical dispersion potential
 (pair-wise force field) to the conventional Kohn-Sham DFT energy.

INPUTS

  ixc=choice of exchange-correlation functional
  natom=number of atoms
  ntypat=number of atom types
  prtvol=printing volume (if >0, print computation parameters)
  typat(natom)=type integer for each atom in cell
  rprimd(3,3)=real space primitive translations
  vdw_tol=tolerance use to converge the potential (a pair of atoms is included
          in potential if its contribution is larger than vdw_tol)
          vdw_tol<0 takes default value (10^-10)
  xred(3,natom)=reduced atomic coordinates
  znucl(ntypat)=atomic number of atom type
  === optional inputs ===
  [qphon(3)]=wavevector of the phonon;
             used only for dynamical matrix computation

OUTPUT

  e_vdw_dftd2=contribution to energy from DFT-D2 dispersion potential
  === optional outputs ===
  [dyn_vdw_dftd2(2,3,natom,3,natom)]=contribution to dynamical matrix from DFT-D2 dispersion potential
  [elt_vdw_dftd2(6+3*natom,6)]=contribution to elastic tensor and internal strains from DFT-D2 disp. pot.
  [fred_vdw_dftd2(3,natom)]=contribution to forces from DFT-D2 dispersion potential
  [str_vdw_dftd2(6)]=contribution to stress tensor from DFT-D2 dispersion potential

NOTES

  Ref.: S. Grimme, Semiempirical GGA-type density functional
        constructed with a long-range dispersion correction,
        J. Comp. Chem. 27, 1787 (2006) [[cite:Grimme2006]]

PARENTS

      respfn,setvtr,stress

CHILDREN

SOURCE

 94 subroutine vdw_dftd2(e_vdw_dftd2,ixc,natom,ntypat,prtvol,typat,rprimd,vdw_tol,xred,znucl,&
 95 &          dyn_vdw_dftd2,elt_vdw_dftd2,fred_vdw_dftd2,str_vdw_dftd2,qphon) ! Optionals
 96 
 97 
 98 !This section has been created automatically by the script Abilint (TD).
 99 !Do not modify the following lines by hand.
100 #undef ABI_FUNC
101 #define ABI_FUNC 'vdw_dftd2'
102 !End of the abilint section
103 
104  implicit none
105 
106 !Arguments ------------------------------------
107 !scalars
108  integer,intent(in) :: ixc,natom,ntypat,prtvol
109  real(dp),intent(in) :: vdw_tol
110  real(dp),intent(out) :: e_vdw_dftd2
111 !arrays
112  integer,intent(in) :: typat(natom)
113  real(dp),intent(in) :: rprimd(3,3),xred(3,natom),znucl(ntypat)
114  real(dp),intent(in),optional :: qphon(3)
115  real(dp),intent(out),optional :: dyn_vdw_dftd2(2,3,natom,3,natom)
116  real(dp),intent(out),optional :: elt_vdw_dftd2(6+3*natom,6)
117  real(dp),intent(out),optional :: fred_vdw_dftd2(3,natom)
118  real(dp),intent(out),optional :: str_vdw_dftd2(6)
119 
120 !Local variables-------------------------------
121 !scalars
122  integer,parameter :: vdw_nspecies=55
123  integer :: ia,ia1,ii,is1,is2,is3,itypat,ja,ja1,jj,jtypat,kk,ll,mu,npairs,nshell,nu
124  logical :: need_dynmat,need_elast,need_forces,need_intstr,need_stress
125  logical :: need_gradient,need_gradient2,newshell,qeq0=.true.
126  real(dp),parameter :: e_conv=(10/Bohr_Ang)**6/Ha_J/Avogadro ! 1 J.nm^6.mol^-1 in Ha.Bohr^6
127  real(dp),parameter :: vdw_d=20._dp,vdw_tol_default=tol10
128  real(dp),parameter :: vdw_s_pbe=0.75_dp, vdw_s_blyp=1.2_dp, vdw_s_b3lyp=1.05_dp
129  real(dp),parameter :: vdw_s_bp86=1.05_dp, vdw_s_tpss=1.0_dp, vdw_s_b97d=1.25_dp
130  real(dp) :: c6,c6r6,ex,fr,gr,gr2,grad,grad2,ph,ph1r,ph1i
131  real(dp) :: r0,r1,r2,r3,rcut,rcut2,rsq,rr,sfact,ucvol,vdw_s
132  character(len=500) :: msg
133  type(atomdata_t) :: atom
134 !arrays
135  integer,allocatable :: ivdw(:)
136  integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
137  real(dp) :: gmet(3,3),gprimd(3,3),mat(3,3),rcart(3),rmet(3,3),vec(3)
138  real(dp),allocatable :: vdw_c6(:,:),vdw_r0(:,:),xred01(:,:)
139  real(dp),parameter :: vdw_c6_dftd2(vdw_nspecies)= &
140 &      (/ 0.14, 0.08, 1.61, 1.61, 3.13, 1.75, 1.23, 0.70, 0.75, 0.63,&
141 &         5.71, 5.71,10.79, 9.23, 7.84, 5.57, 5.07, 4.61,10.80,10.80,&
142 &        10.80,10.80,10.80,10.80,10.80,10.80,10.80,10.80,10.80,10.80,&
143 &        16.99,17.10,16.37,12.64,12.47,12.01,24.67,24.67,24.67,24.67,&
144 &        24.67,24.67,24.67,24.67,24.67,24.67,24.67,24.67,37.32,38.71,&
145 &        38.44,31.74,31.50,29.99, 0.00/)
146  real(dp),parameter :: vdw_r0_dftd2(vdw_nspecies)= &
147 &      (/1.001,1.012,0.825,1.408,1.485,1.452,1.397,1.342,1.287,1.243,&
148 &        1.144,1.364,1.639,1.716,1.705,1.683,1.639,1.595,1.485,1.474,&
149 &        1.562,1.562,1.562,1.562,1.562,1.562,1.562,1.562,1.562,1.562,&
150 &        1.650,1.727,1.760,1.771,1.749,1.727,1.628,1.606,1.639,1.639,&
151 &        1.639,1.639,1.639,1.639,1.639,1.639,1.639,1.639,1.672,1.804,&
152 &        1.881,1.892,1.892,1.881,1.000/)
153  character(len=2),parameter :: vdw_symb(vdw_nspecies)= &
154 &      (/' H','He','Li','Be',' B',' C',' N',' O',' F','Ne',&
155 &        'Na','Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca',&
156 &        'Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn',&
157 &        'Ga','Ge','As','Se','Br','Kr','Rb','Sr',' Y','Zr',&
158 &        'Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn',&
159 &        'Sb','Te',' I','Xe','no'/)
160 
161 ! *************************************************************************
162 
163  DBG_ENTER("COLL")
164 
165 !Extract options
166  need_forces=present(fred_vdw_dftd2)
167  need_stress=present(str_vdw_dftd2)
168  need_dynmat=present(dyn_vdw_dftd2)
169  need_elast=present(elt_vdw_dftd2)
170  need_intstr=present(elt_vdw_dftd2)
171  need_gradient=(need_forces.or.need_stress)
172  need_gradient2=(need_dynmat.or.need_elast.or.need_intstr)
173  if (need_dynmat) then
174    if (.not.present(qphon)) then
175      msg='Dynamical matrix required without a q-vector'
176      MSG_BUG(msg)
177    end if
178    qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
179  end if
180 
181 !Identify type(s) of atoms
182  ABI_ALLOCATE(ivdw,(ntypat))
183  do itypat=1,ntypat
184    ivdw(itypat)=-1;jtypat=0
185    call atomdata_from_znucl(atom,znucl(itypat))
186    do while ((ivdw(itypat)<0).and.(jtypat<vdw_nspecies))
187      jtypat=jtypat+1;if (vdw_symb(jtypat)==atom%symbol) ivdw(itypat)=jtypat
188    end do
189    if (ivdw(itypat)<0) then
190      write(msg,'(3a)') &
191 &     'Van der Waals DFT-D2 correction not available for atom type: ',atom%symbol,' !'
192      MSG_ERROR(msg)
193    end if
194  end do
195 
196 !Select DFT-D2 VdW parameters according to system data
197  vdw_s=e_conv
198  if (ixc==11.or.ixc==-101130.or.ixc==-130101) then
199    vdw_s=vdw_s*vdw_s_pbe
200  else if (ixc==18.or.ixc==-106131.or.ixc==-131106) then
201    vdw_s=vdw_s*vdw_s_blyp
202  else if (ixc==19.or.ixc==-106132.or.ixc==-132106) then
203    vdw_s=vdw_s*vdw_s_bp86
204  else if (ixc==-202231.or.ixc==-231202) then
205    vdw_s=vdw_s*vdw_s_tpss
206  else
207    write(msg,'(a,i8,a)')'  Van der Waals DFT-D2 correction not compatible with ixc=',ixc,' !'
208    MSG_ERROR(msg)
209  end if
210  ABI_ALLOCATE(vdw_c6,(ntypat,ntypat))
211  ABI_ALLOCATE(vdw_r0,(ntypat,ntypat))
212  do itypat=1,ntypat
213    do jtypat=1,ntypat
214      vdw_c6(itypat,jtypat)=sqrt(vdw_c6_dftd2(ivdw(itypat))*vdw_c6_dftd2(ivdw(jtypat)))
215      vdw_r0(itypat,jtypat)=(vdw_r0_dftd2(ivdw(itypat))+vdw_r0_dftd2(ivdw(jtypat)))/Bohr_Ang
216    end do
217  end do
218 
219 !Computation of cut-off radius according to tolerance
220 !We take: r_cut=(s6*max(C6)/tol)**(1/6) and rcut<=75 bohr
221  if (vdw_tol<zero) then
222    rcut=(vdw_s/vdw_tol_default*maxval(vdw_c6))**sixth
223  else
224    rcut=(vdw_s/vdw_tol*maxval(vdw_c6))**sixth
225  end if
226 !rcut=min(rcut,100._dp)
227  rcut2=rcut*rcut
228 
229 !Retrieve cell geometry data
230  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
231 
232 !Map reduced coordinates into [0,1]
233  ABI_ALLOCATE(xred01,(3,natom))
234  do ia=1,natom
235    xred01(1,ia)=xred(1,ia)-aint(xred(1,ia))+half-sign(half,xred(1,ia))
236    xred01(2,ia)=xred(2,ia)-aint(xred(2,ia))+half-sign(half,xred(2,ia))
237    xred01(3,ia)=xred(3,ia)-aint(xred(3,ia))+half-sign(half,xred(3,ia))
238  end do
239 
240 !Set accumulated quantities to zero
241  npairs=0
242  e_vdw_dftd2=zero
243  if (need_forces) fred_vdw_dftd2=zero
244  if (need_stress) str_vdw_dftd2=zero
245  if (need_dynmat) dyn_vdw_dftd2=zero
246  if (need_elast)  elt_vdw_dftd2(1:6,1:6)=zero
247  if (need_intstr) elt_vdw_dftd2(7:6+3*natom,1:6)=zero
248 
249 !Loop over shells of cell replicas
250  nshell=0
251  do
252    newshell=.false.;nshell=nshell+1
253 
254 !  Loop over cell replicas in the shell
255 !  ns1=1+int(rcut*sqrt(SUM(gprimd(:,1)**2))
256 !  ns2=1+int(rcut*sqrt(SUM(gprimd(:,2)**2))
257 !  ns3=1+int(rcut*sqrt(SUM(gprimd(:,3)**2))
258    do is3=-nshell,nshell
259      do is2=-nshell,nshell
260        do is1=-nshell,nshell
261          if (nshell==1.or. &
262 &         abs(is3)==nshell.or.abs(is2)==nshell.or.abs(is1)==nshell) then
263 
264 !          Phase for dynamical matrix
265            if (need_dynmat) then
266              ph1r=one;ph1i=zero  !ph1=exp(-iqL)
267              if (.not.qeq0) then
268                ph=-two_pi*(qphon(1)*is1+qphon(2)*is2+qphon(3)*is3)
269                ph1r=cos(ph);ph1i=sin(ph)
270              end if
271            end if
272 
273 !          Loops over atoms a and b
274            do ia=1,natom
275              itypat=typat(ia)
276              do ja=1,ia
277                jtypat=typat(ja)
278                r1=xred01(1,ia)-xred01(1,ja)-dble(is1)
279                r2=xred01(2,ia)-xred01(2,ja)-dble(is2)
280                r3=xred01(3,ia)-xred01(3,ja)-dble(is3)
281                rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3 &
282 &               +two*(rmet(2,1)*r2*r1+rmet(3,2)*r3*r2+rmet(3,1)*r1*r3)
283 
284 !              Select atomic pairs (a,b) and avoid atom_a=atom_b
285                if (rsq>=tol16.and.rsq<=rcut2) then
286 
287 !                Data for the selected pair
288                  npairs=npairs+1;newshell=.true.
289                  sfact=vdw_s;if (ia==ja) sfact=half*sfact
290                  rr=sqrt(rsq)
291                  c6=vdw_c6(itypat,jtypat)
292                  r0=vdw_r0(itypat,jtypat)
293 
294 !                Computation of pair-wise potential
295                  ex=exp(-vdw_d*(rr/r0-one))
296                  fr=one/(one+ex)
297                  c6r6=c6/rr**6
298 
299 !                Contribution to energy
300                  e_vdw_dftd2=e_vdw_dftd2-sfact*fr*c6r6
301 
302                  if (need_gradient.or.need_gradient2) then
303                    gr=(vdw_d/r0)*(fr**2)*ex
304                    grad=-sfact*(gr-six*fr/rr)*c6r6/rr
305                    rcart(1)=rprimd(1,1)*r1+rprimd(1,2)*r2+rprimd(1,3)*r3
306                    rcart(2)=rprimd(2,1)*r1+rprimd(2,2)*r2+rprimd(2,3)*r3
307                    rcart(3)=rprimd(3,1)*r1+rprimd(3,2)*r2+rprimd(3,3)*r3
308 
309 !                  Contribution to forces
310                    if (need_forces.and.ia/=ja) then
311                      vec(1:3)=grad*rcart(1:3)
312                      fred_vdw_dftd2(1:3,ia)=fred_vdw_dftd2(1:3,ia)+vec(1:3)
313                      fred_vdw_dftd2(1:3,ja)=fred_vdw_dftd2(1:3,ja)-vec(1:3)
314                    end if
315 
316 !                  Contribution to stress tensor
317                    if (need_stress) then
318                      do mu=1,6
319                        ii=alpha(mu);jj=beta(mu)
320                        str_vdw_dftd2(mu)=str_vdw_dftd2(mu)+grad*rcart(ii)*rcart(jj)
321                      end do
322                    end if
323 
324                    if (need_gradient2) then
325                      gr2=(vdw_d/r0)*gr*(2*fr*ex-one)
326                      grad2=-sfact*(gr2-13._dp*gr/rr+48._dp*fr/rr**2)*c6r6/rr**2
327 
328 !                    Contribution to dynamical matrix (phase factors are subtle!)
329                      if (need_dynmat) then
330                        mat(1:3,1)=grad2*rcart(1:3)*rcart(1) ; mat(1,1)=mat(1,1)+grad
331                        mat(1:3,2)=grad2*rcart(1:3)*rcart(2) ; mat(2,2)=mat(2,2)+grad
332                        mat(1:3,3)=grad2*rcart(1:3)*rcart(3) ; mat(3,3)=mat(3,3)+grad
333                        if (ia/=ja) then
334                          do ii=1,3
335                            dyn_vdw_dftd2(1,1:3,ia,ii,ia)=dyn_vdw_dftd2(1,1:3,ia,ii,ia)+mat(1:3,ii)
336                            dyn_vdw_dftd2(1,1:3,ja,ii,ja)=dyn_vdw_dftd2(1,1:3,ja,ii,ja)+mat(1:3,ii)
337                            dyn_vdw_dftd2(1,1:3,ia,ii,ja)=dyn_vdw_dftd2(1,1:3,ia,ii,ja)-mat(1:3,ii)*ph1r
338                            dyn_vdw_dftd2(2,1:3,ia,ii,ja)=dyn_vdw_dftd2(2,1:3,ia,ii,ja)-mat(1:3,ii)*ph1i
339                            dyn_vdw_dftd2(1,1:3,ja,ii,ia)=dyn_vdw_dftd2(1,1:3,ja,ii,ia)-mat(1:3,ii)*ph1r
340                            dyn_vdw_dftd2(2,1:3,ja,ii,ia)=dyn_vdw_dftd2(2,1:3,ja,ii,ia)+mat(1:3,ii)*ph1i
341                          end do
342                        else if (.not.qeq0) then
343                          do ii=1,3
344                            dyn_vdw_dftd2(1,1:3,ia,ii,ia)=dyn_vdw_dftd2(1,1:3,ia,ii,ia) &
345 &                           +two*mat(1:3,ii)*(one-ph1r)
346                          end do
347                        end if
348                      end if
349 
350 !                    Contribution to elastic tensor
351                      if (need_elast) then
352                        do mu=1,6
353                          ii=alpha(mu);jj=beta(mu)
354                          do nu=1,6
355                            kk=alpha(nu);ll=beta(nu)
356                            elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) &
357 &                           +grad2*rcart(ii)*rcart(jj)*rcart(kk)*rcart(ll)
358                            if (ii==kk) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) &
359 &                           +half*grad*rcart(jj)*rcart(ll)
360                            if (ii==ll) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) &
361 &                           +half*grad*rcart(jj)*rcart(kk)
362                            if (jj==kk) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) &
363 &                           +half*grad*rcart(ii)*rcart(ll)
364                            if (jj==ll) elt_vdw_dftd2(mu,nu)=elt_vdw_dftd2(mu,nu) &
365 &                           +half*grad*rcart(ii)*rcart(kk)
366                          end do
367                        end do
368                      end if
369 
370 !                    Contribution to internal strains
371                      if (need_intstr.and.ia/=ja) then
372                        ia1=6+3*(ia-1);ja1=6+3*(ja-1)
373                        do mu=1,6
374                          ii=alpha(mu);jj=beta(mu)
375                          vec(1:3)=grad2*rcart(ii)*rcart(jj)*rcart(1:3)
376                          vec(ii)=vec(ii)+half*grad*rcart(jj)
377                          vec(jj)=vec(jj)+half*grad*rcart(ii)
378                          elt_vdw_dftd2(ia1+1:ia1+3,mu)=elt_vdw_dftd2(ia1+1:ia1+3,mu)+vec(1:3)
379                          elt_vdw_dftd2(ja1+1:ja1+3,mu)=elt_vdw_dftd2(ja1+1:ja1+3,mu)-vec(1:3)
380                        end do
381                      end if
382 
383                    end if ! Computation of 2nd gradient
384                  end if ! Computation of gradient
385                end if   ! Pairs selection
386              end do     ! Loop over atom b
387            end do       ! Loop over atom a
388          end if         ! Triple loop over cell replicas in shell
389        end do
390      end do
391    end do
392    if(.not.newshell) exit ! Check if new shell must be calculated
393  end do ! Loop over shells
394 
395 !Gradients: convert them from cartesian to reduced coordinates
396  if (need_forces) then
397    do ia=1,natom
398      call grad_cart2red(fred_vdw_dftd2(:,ia))
399    end do
400  end if
401  if (need_dynmat) then
402    do ja=1,natom
403      do ia=1,natom
404        do kk=1,merge(2,1,qeq0)
405          do ii=1,3
406            vec(1:3)=dyn_vdw_dftd2(kk,1:3,ia,ii,ja)
407            call grad_cart2red(vec)
408            dyn_vdw_dftd2(kk,1:3,ia,ii,ja)=vec(1:3)
409          end do
410          do ii=1,3
411            vec(1:3)=dyn_vdw_dftd2(kk,ii,ia,1:3,ja)
412            call grad_cart2red(vec)
413            dyn_vdw_dftd2(kk,ii,ia,1:3,ja)=vec(1:3)
414          end do
415        end do
416      end do
417    end do
418  end if
419  if (need_intstr) then
420    do mu=1,6
421      ia1=6
422      do ia=1,natom
423        call grad_cart2red(elt_vdw_dftd2(ia1+1:ia1+3,mu))
424        ia1=ia1+3
425      end do
426    end do
427  end if
428 
429 !DEBUG
430 !write(77,*) "---------------"
431 !write(77,*) "E=",e_vdw_dftd2
432 !if (need_forces) then
433 ! do ia=1,natom
434 !  write(77,*) "F=",ia,fred_vdw_dftd2(:,ia)
435 ! end do
436 !end if
437 !if (need_stress) write(77,*) "S=",str_vdw_dftd2(:)
438 !if (need_dynmat) then
439 ! do ia=1,natom
440 !  do ii=1,3
441 !   do ja=1,natom
442 !    write(77,*) "D=",ia,ii,ja,dyn_vdw_dftd2(:,:,ja,ii,ia)
443 !   end do
444 !  end do
445 ! end do
446 !end if
447 !if (need_elast) then
448 ! do ii=1,6
449 !  write(77,*) "e=",ii,elt_vdw_dftd2(1:6,ii)
450 ! end do
451 !end if
452 !if (need_intstr) then
453 ! do ii=1,6
454 !  do ia=1,natom
455 !   write(77,*) "I=",ii,ia,elt_vdw_dftd2(7+3*(ia-1):9+3*(ia-1),ii)
456 !  end do
457 ! end do
458 !end if
459 !flush(77)
460 !DEBUG
461 
462 !Stress tensor: divide by volume
463  if (need_stress) str_vdw_dftd2=str_vdw_dftd2/ucvol
464 
465 !Printing
466  if (prtvol>0) then
467    write(msg,'(10a)') ch10,&
468 &   '  --------------------------------------------------------------',ch10,&
469 &   '  Van der Waals DFT-D2 semi-empirical dispersion potential added',ch10,&
470 &   '      with following parameters:',ch10,&
471 &   '      Specie  C6 (J.nm^6.mol^-1)  R0 (Ang)',ch10,&
472 &   '      ------------------------------------'
473    call wrtout(std_out,msg,'COLL')
474    do itypat=1,ntypat
475      write(msg,'(9X,a2,11X,f5.2,8X,f6.3)') &
476 &     vdw_symb(ivdw(itypat)),vdw_c6_dftd2(ivdw(itypat)),vdw_r0_dftd2(ivdw(itypat))
477      call wrtout(std_out,msg,'COLL')
478    end do
479    write(msg,'(2a,f6.2,2a,f6.2,2a,f6.2,a)') ch10,&
480 &   '      Scaling factor   = ',vdw_s/e_conv,ch10,&
481 &   '      Damping parameter= ',vdw_d,ch10,&
482 &   '      Cut-off radius   = ',rcut,' bohr'
483    call wrtout(std_out,msg,'COLL')
484    write(msg,'(2a,i8,2a,es14.5,4a)') ch10,&
485 &   '      Number of pairs contributing = ',npairs,ch10,&
486 &   '      DFT-D2 energy contribution   = ',e_vdw_dftd2,' Ha',ch10,&
487 &   '  --------------------------------------------------------------',ch10
488    call wrtout(std_out,msg,'COLL')
489  end if
490 
491  ABI_DEALLOCATE(ivdw)
492  ABI_DEALLOCATE(vdw_c6)
493  ABI_DEALLOCATE(vdw_r0)
494  ABI_DEALLOCATE(xred01)
495 
496  DBG_EXIT("COLL")
497 
498  contains

vdw_dftd2/grad_cart2red [ Functions ]

[ Top ] [ vdw_dftd2 ] [ Functions ]

NAME

 grad_cart2red

FUNCTION

 Convert gradients from cartesian to reduced coordinates

PARENTS

      vdw_dftd2

CHILDREN

SOURCE

516 subroutine grad_cart2red(grad)
517 
518 
519 !This section has been created automatically by the script Abilint (TD).
520 !Do not modify the following lines by hand.
521 #undef ABI_FUNC
522 #define ABI_FUNC 'grad_cart2red'
523 !End of the abilint section
524 
525 implicit none
526 
527 !Arguments ------------------------------------
528  real(dp),intent(inout) :: grad(3)
529 !Local variables-------------------------------
530  real(dp) :: tmp(3)
531 
532 ! *********************************************************************
533 
534    tmp(1)=rprimd(1,1)*grad(1)+rprimd(2,1)*grad(2)+rprimd(3,1)*grad(3)
535    tmp(2)=rprimd(1,2)*grad(1)+rprimd(2,2)*grad(2)+rprimd(3,2)*grad(3)
536    tmp(3)=rprimd(1,3)*grad(1)+rprimd(2,3)*grad(2)+rprimd(3,3)*grad(3)
537    grad(1:3)=tmp(1:3)
538 
539  end subroutine grad_cart2red