1      subroutine vdwlj(ev,nac,iac,nad,iad,coo,iconn,ityp,forces,potscl,
2     &                 iopt)
3      implicit real (a-h,o-z), integer (i-n)
4      parameter (mxcon=10)
5      parameter (mxac=3*mxcon)
6      parameter (mxad=9*mxcon)
7      parameter (mxamb=1590)
8      integer ambcls
9      common /typpar/ ambchg(mxamb),ambcls(mxamb)
10      parameter (mxcls=50)
11      common /clspar/ ambvdwr(mxcls),ambvdwe(mxcls),mapagf(mxcls)
12      parameter (mxgff=72)
13      integer amb2gf
14      common /gftyp/  gfvdw(2,mxgff),amb2gf(mxamb)
15      logical usecut,usesw
16      common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw
17      common /athlp/  iatoms, mxnat
18      logical box,cell,fast
19      common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast
20      integer*2 ityp
21      dimension ded(3),vr(3)
22      dimension coo(3,*),ityp(*),iconn(mxcon+1,*),
23     &          forces(3,*),potscl(*),iopt(*),
24     &          nac(*),nad(*),iac(mxac,*),iad(mxad,*)
25
26      ev = 0.0e0
27      v14sc = 0.5e0
28
29      do i=1,iatoms
30
31         i1 = int(ityp(i))
32         if (i1.gt.0) then
33            il = ambcls(i1)
34            vdwr1 = ambvdwr(il)
35            vdwe1 = ambvdwe(il)
36         elseif (i1.le.0) then
37            i1 = iabs(i1)
38            vdwr1 = gfvdw(1,i1)
39            vdwe1 = gfvdw(2,i1)
40         endif
41
42         do j=i+1,iatoms
43            potscl(j) = 1.0e0
44         end do
45
46         do j=1,iconn(1,i)
47            jj = iconn(j+1,i)
48            if (jj.gt.0) then
49               potscl(jj) = 0.0e0
50            endif
51         end do
52
53         do j=1,nac(i)
54            potscl(iac(j,i)) = 0.0e0
55         end do
56
57         do j=1,nad(i)
58            potscl(iad(j,i)) = v14sc
59         end do
60
61         if (vdwe1.ne.0.0e0) then
62
63           do k=i+1,iatoms
64
65            if (potscl(k).ne.0.0e0.and.
66     &          (iopt(i).eq.1.or.iopt(k).eq.1) ) then
67
68               i2 = int(ityp(k))
69               if (i2.gt.0) then
70                  kl = ambcls(i2)
71                  vdwr2 = ambvdwr(kl)
72                  vdwe2 = ambvdwe(kl)
73               elseif (i2.le.0) then
74                  i2 = iabs(i2)
75                  vdwr2 = gfvdw(1,i2)
76                  vdwe2 = gfvdw(2,i2)
77               endif
78
79               if (vdwe2.ne.0.0e0) then
80
81                  do j=1,3
82                     vr(j) = coo(j,i) - coo(j,k)
83                  end do
84
85                  if (box) call reddis(vr)
86
87                  rv2 = vr(1)*vr(1) + vr(2)*vr(2) + vr(3)*vr(3)
88
89                  if (rv2.le.cutof2) then
90
91c          [ (Rmin)**12       (Rmin)**6 ]
92c e = eps  [ (----)     - 2.0 (----)    ]
93c          [ ( r  )           ( r  )    ]
94
95
96c alternatively we could precalculate these vdwr(il,kl)
97c                                           vdwe(il,kl)
98                     rsum = vdwr1 + vdwr2
99                     epsm = sqrt(vdwe1 * vdwe2)
100
101                     epsm = epsm * potscl(k)
102                     rv   = sqrt(rv2)
103		     rs2  = rsum*rsum
104		     rs3  = rs2*rsum
105                     p6   = (rs3*rs3) / (rv2*rv2*rv2)
106                     p12  = p6 * p6
107
108                     e    = epsm * (p12 - 2.0e0*p6)
109                     de   = epsm * (p12 - p6) * (-12.0e0/rv)
110
111                     de   = de / rv
112
113                     do j=1,3
114                        ded(j) = de * vr(j)
115                     end do
116
117                     ev   = ev + e
118
119                     do j=1,3
120                        forces(j,i) = forces(j,i) + ded(j)
121                        forces(j,k) = forces(j,k) - ded(j)
122                     end do
123
124                  end if
125
126               end if
127            end if
128           end do
129         end if
130      end do
131
132      return
133      end
134
135