← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/nwxc/nwxc_vdw_der.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwxc/nwxc_vdw_der.F
 
1
c
 
2
c     empirical dispersion: derivatives
 
3
c
 
4
      subroutine nwxc_vdw_der(s6,s8,sr6,sr8,n,x,z,force)
 
5
c
 
6
c     S. Grimme J Comp Chem 25, 1463 (2004)
 
7
c     U. Zimmerli, M Parrinello and P. Koumoutsakos, JCP. 120, 2693 (2004)
 
8
c     Q. Wu and W. Yang, JCP. 116, 515 (2002)
 
9
c
 
10
      implicit none
 
11
c
 
12
#include "mafdecls.fh"
 
13
#include "errquit.fh"
 
14
#include "nwxc_vdw.fh"
 
15
#include "global.fh"
 
16
c
 
17
      double precision s6,s8,sr6,sr8
 
18
      integer n
 
19
      double precision x(3,n),force(3,n)
 
20
      integer z(n)
 
21
c
 
22
      integer i,j,k,A,l_cnij,k_cnij,l_cnijk,k_cnijk
 
23
      double precision c6ij_sk
 
24
      external c6ij_sk
 
25
      double precision drajdxa
 
26
      double precision ff1,rr,ff
 
27
      double precision fdmp,f1dmp,cnA,cnj
 
28
      external c6cn,crd_nr
 
29
      double precision c6cn,crd_nr
 
30
      double precision fac6,fac8,fdmp6,fdmp6a,fdmp8,fdmp8a,Qfac
 
31
      double precision rAj,rAk,rjk,r0aj,r0ak,r0jk,c6Aj,grad_c6(3)
 
32
      double precision dxAj,dyAj,dzAj,dxAk,dyAk,dzAk,dxjk,dyjk,dzjk
 
33
      double precision tmp6,tmp6a,tmp8,tmp8a
 
34
c
 
35
c     Derivatives of Grimme dispersion term
 
36
c
 
37
c  DFT-D1 / DFT-D2
 
38
c
 
39
      if (ivdw.le.2) then
 
40
         do A=1,n
 
41
            force(1,A)=0d0
 
42
            force(2,A)=0d0
 
43
            force(3,A)=0d0
 
44
            if (Z(A).ne.0) then
 
45
              do j=1,n
 
46
                 if(A.ne.j) then
 
47
                    rAj=sqrt(
 
48
     +                 (x(1,A)-x(1,j))**2 +
 
49
     +                 (x(2,A)-x(2,j))**2 +
 
50
     +                 (x(3,A)-x(3,j))**2)
 
51
                    r0aj=r0(z(A))+r0(z(j))
 
52
                    ff= fdmp(rAj,r0aj)
 
53
                    ff1= f1dmp(rAj,r0aj,ff)
 
54
                    rr=c6ij_sk(A,j,z)/(rAj**6)*
 
55
     *               ((-6d0*ff/rAj)+ff1)
 
56
                    do i=1,3
 
57
                       drAjdxa=(x(i,A)-x(i,j))/rAj
 
58
                       force(i,A)=force(i,A)-rr*drAjdxa
 
59
                    enddo
 
60
                 endif
 
61
              enddo
 
62
            endif
 
63
         enddo
 
64
         if(abs(s6-1d0).gt.1d-9) 
 
65
     F        call dscal(3*n,s6,force,1)
 
66
c
 
67
c DFT-D3
 
68
c
 
69
      else if (ivdw.eq.3) then
 
70
c
 
71
c        Precompute coordinate derivatives C6 dependency
 
72
c
 
73
         if (.not.ma_push_get(mt_dbl,3*n,'xcvdw cnij',l_cnij,k_cnij))
 
74
     &      call errquit('xcvdw cnij: cannot allocate cnij',0, MA_ERR)
 
75
         if (.not.ma_push_get(mt_dbl,3*n*n,'vdw cnijk',l_cnijk,k_cnijk))
 
76
     &      call errquit('vdw cnijk: cannot allocate cnijk',0, MA_ERR)
 
77
c
 
78
         call crd_nr_der(n,x,z,dbl_mb(k_cnij),dbl_mb(k_cnijk))
 
79
c
 
80
         do A=1,n
 
81
           force(1,A)=0.0d0
 
82
           force(2,A)=0.0d0
 
83
           force(3,A)=0.0d0
 
84
           if (Z(A).ne.0) then
 
85
             do j=1,n
 
86
               if(A.ne.j) then      
 
87
                  dxAj=x(1,A)-x(1,j)
 
88
                  dyAj=x(2,A)-x(2,j)
 
89
                  dzAj=x(3,A)-x(3,j)
 
90
                  rAj=dxAj**2+dyAj**2+dzAj**2
 
91
c
 
92
c                 Two center derivatives. Grimme uses screening to reduce 
 
93
c                 computational work
 
94
c
 
95
c                 Screening r^2 distance vs threshold of 20000.0
 
96
c
 
97
                  if (rAj.gt.20000.d0) goto 901
 
98
c
 
99
c                 Factors
 
100
c
 
101
                  r0aj=r0AB(z(A),z(j))
 
102
                  Qfac=Qatom(z(A))*Qatom(z(j))
 
103
                  fac6=(dsqrt(rAj)/(sr6*r0aj))**(-alpha)
 
104
                  fac8=(dsqrt(rAj)/(sr8*r0aj))**(-(alpha+2.0d0))
 
105
                  fdmp6=1.0d0/(1.0d0+6.0d0*fac6)
 
106
                  fdmp8=1.0d0/(1.0d0+6.0d0*fac8)
 
107
c
 
108
c                 Coordination dependent C6_AB value
 
109
c
 
110
                  cnA=crd_nr(A,n,x,z)
 
111
                  cnj=crd_nr(j,n,x,z)
 
112
                  c6Aj=c6cn(z(A),z(j),cnA,cnj)
 
113
c
 
114
c                 Get gradient for coordination number dependent C6
 
115
c
 
116
                  call c6_grad(grad_c6,A,j,A,x,z,n,
 
117
     &                         dbl_mb(k_cnij),dbl_mb(k_cnijk))
 
118
c
 
119
                  tmp6=6.0d0*fdmp6*s6*c6Aj/(rAj**4.0d0)
 
120
                  tmp8=6.0d0*fdmp8*s8*c6Aj*Qfac/(rAj**5.0d0)
 
121
c
 
122
c                 dx contribution to A
 
123
c
 
124
                  tmp6a=tmp6*dxAj
 
125
                  tmp8a=tmp8*dxAj
 
126
                  force(1,A)=force(1,A)
 
127
     $              +(1.0d0-fdmp6*fac6*alpha)*tmp6a
 
128
     $              -fdmp6*s6*grad_c6(1)/(rAj**3.0d0)
 
129
     $              +(4.0d0-3.0d0*fdmp8*fac8*(alpha+2.0d0))*tmp8a
 
130
     $              -3.0d0*fdmp8*s8*grad_c6(1)*Qfac/(rAj**4.0d0)
 
131
c
 
132
c                 dy contribution to A
 
133
c
 
134
                  tmp6a=tmp6*dyAj
 
135
                  tmp8a=tmp8*dyAj
 
136
                  force(2,A)=force(2,A)
 
137
     $              +(1.0d0-fdmp6*fac6*alpha)*tmp6a
 
138
     $              -fdmp6*s6*grad_c6(2)/(rAj**3.0d0)
 
139
     $              +(4.0d0-3.0d0*fdmp8*fac8*(alpha+2.0d0))*tmp8a
 
140
     $              -3.0d0*fdmp8*s8*grad_c6(2)*Qfac/(rAj**4.0d0)
 
141
c
 
142
c                 dz contribution to A
 
143
c
 
144
                  tmp6a=tmp6*dzAj
 
145
                  tmp8a=tmp8*dzAj
 
146
                  force(3,A)=force(3,A)
 
147
     $              +(1.0d0-fdmp6*fac6*alpha)*tmp6a
 
148
     $              -fdmp6*s6*grad_c6(3)/(rAj**3.0d0)
 
149
     $              +(4.0d0-3.0d0*fdmp8*fac8*(alpha+2.0d0))*tmp8a
 
150
     $              -3.0d0*fdmp8*s8*grad_c6(3)*Qfac/(rAj**4.0d0)
 
151
 901              continue
 
152
               endif
 
153
             enddo
 
154
c
 
155
c            Three center derivatives. Grimme uses aggressive screening
 
156
c            to get this N^3 contribution back to N^2
 
157
c
 
158
             do j=2,n
 
159
               if(A.ne.j) then      
 
160
                  rAj=sqrt(
 
161
     +                 (x(1,A)-x(1,j))**2 +
 
162
     +                 (x(2,A)-x(2,j))**2 +
 
163
     +                 (x(3,A)-x(3,j))**2)
 
164
                  r0aj=r0AB(z(A),z(j))
 
165
c
 
166
c                 Screening per Grimme
 
167
c
 
168
                  if (rAj.gt.1600d0*r0aj/r0AB(1,1)) goto 910
 
169
c
 
170
c                 Third center involved
 
171
c
 
172
                  do k=1,j-1
 
173
                     if(A.ne.k) then      
 
174
                       dxAk=x(1,A)-x(1,k)
 
175
                       dyAk=x(2,A)-x(2,k)
 
176
                       dzAk=x(3,A)-x(3,k)
 
177
                       rAk=dxAk**2+dyAk**2+dzAk**2
 
178
                       r0ak=r0AB(z(A),z(k))
 
179
                       dxjk=x(1,j)-x(1,k)
 
180
                       dyjk=x(2,j)-x(2,k)
 
181
                       dzjk=x(3,j)-x(3,k)
 
182
                       rjk=dxjk**2+dyjk**2+dzjk**2
 
183
                       r0jk=r0AB(z(j),z(k))
 
184
c
 
185
c                      Screening r^2 distance vs threshold of 1600.0*(radii Ak)
 
186
c
 
187
                       if ((rAk.gt.1600.0d0*r0ak/r0AB(1,1)).or.
 
188
     $                     (rjk.gt.1600.0d0*r0jk/r0AB(1,1))) goto 911
 
189
c
 
190
c                      Get gradient for coordination number dependent C6 for three centers
 
191
c
 
192
                       call c6_grad(grad_c6,j,k,A,x,z,n,
 
193
     &                              dbl_mb(k_cnij),dbl_mb(k_cnijk))
 
194
                       fac6=(sr6*r0jk/dsqrt(rjk))**(alpha)
 
195
                       fac8=(sr8*r0jk/dsqrt(rjk))**(alpha+2.0d0)
 
196
                       fdmp6=1.0d0/(1.0d0+6.0d0*fac6)
 
197
                       fdmp8=1.0d0/(1.0d0+6.0d0*fac8)
 
198
c
 
199
c                      dx, dy, and dz contribution to A
 
200
c
 
201
                       Qfac=Qatom(z(j))*Qatom(z(k))
 
202
                       force(1,A)=force(1,A)
 
203
     $                      -fdmp6*s6*grad_c6(1)/(rjk**3.0d0)
 
204
     $                      -3.0d0*fdmp8*s8*grad_c6(1)*Qfac/(rjk**4.0d0)
 
205
                       force(2,A)=force(2,A)
 
206
     $                      -fdmp6*s6*grad_c6(2)/(rjk**3.0d0)
 
207
     $                      -3.0d0*fdmp8*s8*grad_c6(2)*Qfac/(rjk**4.0d0)
 
208
                       force(3,A)=force(3,A)
 
209
     $                      -fdmp6*s6*grad_c6(3)/(rjk**3.0d0)
 
210
     $                      -3.0d0*fdmp8*s8*grad_c6(3)*Qfac/(rjk**4.0d0)
 
211
                     endif
 
212
 911              continue
 
213
                  enddo
 
214
 910           continue
 
215
               endif
 
216
             enddo
 
217
           endif
 
218
         enddo
 
219
         if (.not.ma_pop_stack(l_cnijk))
 
220
     $      call errquit('xcvdw cnijk: cannot pop cnijk',4, MA_ERR)
 
221
         if (.not.ma_pop_stack(l_cnij))
 
222
     $      call errquit('xcvdw cnij: cannot pop cnij',4, MA_ERR)
 
223
      endif
 
224
c
 
225
#ifdef DEBUG
 
226
      write(6,*) ' gradient vdw called'
 
227
#endif
 
228
      return
 
229
      end