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