1      subroutine argos_cafe_fw(iwfr,iwto,xw,fw,iwdt,iwatm,iwq,
2     + lpbc,eww,vdw,chg,
3     + mwb,nwb,nbp,ibnd,bnd,rbnd,mwh,nwh,nhp,iang,ang,rang,rub,
4     + mwd,nwd,ndp,idih,dih,rdih,mwo,nwo,nop,iimp,dimp,rimp,
5     + mwt,nwt,idwt,mwn,nwn,idwn)
6c
7c $Id$
8c
9      implicit none
10c
11#include "argos_cafe_common.fh"
12c
13      integer iwfr,iwto
14      integer mwb,mwh,mwd,mwo,nbp,nhp,ndp,nop,mwt,mwn
15      integer nwb,nwh,nwd,nwo,nwt,nwn
16      real*8 xw(mwm,3,mwa),fw(mwm,3,mwa,2)
17      integer iwdt(mwm),iwq(mwa),iwatm(mwa)
18      logical lpbc
19      real*8 vdw(mat,mat,map,mset),chg(mqt,mqp,mset)
20      integer ibnd(mwb,3),iang(mwh,4),idih(mwd,5),iimp(mwo,5)
21      real*8 bnd(mwb,nbp,mset),ang(mwh,nhp,mset)
22      real*8 dih(mwd,ndp,mset),dimp(mwo,nop,mset)
23      real*8 rbnd(mwb,2),rang(mwh,2),rub(mwh,2),rdih(mwd,2),rimp(mwo,2)
24c
25c      real*8 ca6(mat,mat,6),ca12(mat,mat,6)
26c      real*8 cb6(mat,mat,6),cb12(mat,mat,6)
27c      integer iwl(mwm,miw2),
28c
29       integer idwt(0:mwt,2),idwn(0:mwn,2)
30c
31c      real*8 cdwb(mwb,6),ddwb(mwb,6)
32c      integer iwbs(mwb),idwb(mwb),jdwb(mwb),iwatm(mwa)
33c      real*8 cdwh(mwh,6),ddwh(mwh,6)
34c      integer idwh(mwh),jdwh(mwh),kdwh(mwh)
35c      real*8 cdwd(mwd,6),ddwd(mwd,6),edwd(mwd,6)
36c      integer idwd(mwd),jdwd(mwd),kdwd(mwd),ldwd(mwd)
37c      real*8 cdwo(mwo,6),ddwo(mwo,6)
38c      integer idwo(mwo),jdwo(mwo),kdwo(mwo),ldwo(mwo)
39c      real*8 uwb(mwb),uwh(mwh),uwd(mwd),uwo(mwo)
40c
41      integer iwb,iwa,jwa,iwm,iwh,kwa,iwd,lwa,iwo,iwt,iwn
42      real*8 bond,for,rwx1,rwx2,rwx3,rww,rwwi,dbond,dfor,dfw1,dfw2,dfw3
43      real*8 angle,xwij1,xwij2,xwij3,xwkj1,xwkj2,xwkj3,rwij2,rwij2i
44      real*8 rwkj2,rwkj2i,cphi,phi,dangle,sphi,rmul
45      real*8 xwkl1,xwkl2,xwkl3,xwik1,xwik2,xwik3,xwjl1,xwjl2,xwjl3
46      real*8 xm1,xm2,xm3,xn1,xn2,xn3,rm2i,rn2i,rmni,s,rpa
47      real*8 xd1,xd2,xd3,xe1,xe2,xe3,dfwi1,dfwi2,dfwi3
48      real*8 dfwj1,dfwj2,dfwj3,dfwk1,dfwk2,dfwk3,dfwl1,dfwl2,dfwl3
49      real*8 danglep,c6p1,c12p1,c6p2,c12p2,qip1,qjp1,qip2,qjp2
50      real*8 c6,c12,c6t,c12t,qit,qjt,cf6,cf12,qi,qj,q,qp1,qp2
51      real*8 ep2l,ep3l,ep2q,ep3q,rxx,rxy,rxz,r2,r2i,r1i,r6i,dfw
52      real*8 rwikji,sphii,qij,rwi,ferfc,fderfc,eww(mpe,2)
53      real*8 etermq,eterml,eub
54c
55#include "argos_cafe_funcs_dec.fh"
56#include "bitops_decls.fh"
57#include "argos_cafe_funcs_sfn.fh"
58#include "bitops_funcs.fh"
59c
60      c6t=zero
61      c12t=zero
62      qit=zero
63      qjt=zero
64      qp1=zero
65      qp2=zero
66c
67      do 10 iwb=1,nwb
68      if(iand(ibnd(iwb,3),icnstr).eq.0) then
69      iwa=ibnd(iwb,1)
70      jwa=ibnd(iwb,2)
71      bond=bnd(iwb,1,iset)
72      for=bnd(iwb,2,iset)
73      rbnd(iwb,2)=zero
74      do 20 iwm=iwfr,iwto
75      rwx1=xw(iwm,1,iwa)-xw(iwm,1,jwa)
76      rwx2=xw(iwm,2,iwa)-xw(iwm,2,jwa)
77      rwx3=xw(iwm,3,iwa)-xw(iwm,3,jwa)
78      rww=sqrt(rwx1**2+rwx2**2+rwx3**2)
79      if(rww.lt.tiny) then
80      rwwi=one
81      else
82      rwwi=one/rww
83      endif
84      dbond=rww-bond
85      if(iand(iwdt(iwm),mdynam).eq.ldynam)
86     + rbnd(iwb,2)=rbnd(iwb,2)+half*for*(rww-bond)**2
87      dfor=for*dbond*rwwi
88      dfw1=dfor*rwx1
89      dfw2=dfor*rwx2
90      dfw3=dfor*rwx3
91      fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfw1
92      fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)+dfw1
93      fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfw2
94      fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)+dfw2
95      fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfw3
96      fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)+dfw3
97      if(ip2(6))
98     + ep2(1)=ep2(1)+half*bnd(iwb,2,2)*(rww-bnd(iwb,1,2))**2
99      if(ip3(6))
100     + ep3(1)=ep3(1)+half*bnd(iwb,2,3)*(rww-bnd(iwb,1,3))**2
101      if(ith(6)) then
102      deriv(6,1)=deriv(6,1)+
103     + dbond*(half*dbond*bnd(iwb,2,4)-for*bnd(iwb,1,4))
104      endif
105   20 continue
106      eww(1,1)=eww(1,1)+rbnd(iwb,2)
107      if(ip2(6)) ep2(1)=ep2(1)-rbnd(iwb,2)
108      if(ip3(6)) ep3(1)=ep3(1)-rbnd(iwb,2)
109      endif
110      if(ipme.ne.0) then
111      iwa=ibnd(iwb,1)
112      jwa=ibnd(iwb,2)
113      qij=chg(iwq(iwa),1,iset)*chg(iwq(jwa),1,iset)
114      do 21 iwm=iwfr,iwto
115      rwx1=xw(iwm,1,iwa)-xw(iwm,1,jwa)
116      rwx2=xw(iwm,2,iwa)-xw(iwm,2,jwa)
117      rwx3=xw(iwm,3,iwa)-xw(iwm,3,jwa)
118      rww=sqrt(rwx1**2+rwx2**2+rwx3**2)
119      rwi=one/rww
120      ferfc=one-erfc(ealpha*rww)
121      fderfc=-(ealpha*derfc(ealpha*rww))
122      epmecw=epmecw-ferfc*qij*rwi
123      eww(9,1)=eww(9,1)-ferfc*qij*rwi
124      dfor=-(qij*rwi*rwi*(ferfc*rwi-fderfc))
125      dfw1=dfor*rwx1
126      dfw2=dfor*rwx2
127      dfw3=dfor*rwx3
128      fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfw1
129      fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)+dfw1
130      fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfw2
131      fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)+dfw2
132      fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfw3
133      fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)+dfw3
134      vpmeb(1)=vpmeb(1)+dfw1*rwx1
135      vpmeb(2)=vpmeb(2)+dfw2*rwx1
136      vpmeb(3)=vpmeb(3)+dfw3*rwx1
137      vpmeb(4)=vpmeb(4)+dfw2*rwx2
138      vpmeb(5)=vpmeb(5)+dfw3*rwx2
139      vpmeb(6)=vpmeb(6)+dfw3*rwx3
140   21 continue
141      endif
142   10 continue
143      do 40 iwh=1,nwh
144      iwa=iang(iwh,1)
145      jwa=iang(iwh,2)
146      kwa=iang(iwh,3)
147      angle=ang(iwh,1,iset)
148      for=ang(iwh,2,iset)
149      rang(iwh,2)=zero
150      do 50 iwm=iwfr,iwto
151      xwij1=xw(iwm,1,iwa)-xw(iwm,1,jwa)
152      xwij2=xw(iwm,2,iwa)-xw(iwm,2,jwa)
153      xwij3=xw(iwm,3,iwa)-xw(iwm,3,jwa)
154      xwkj1=xw(iwm,1,kwa)-xw(iwm,1,jwa)
155      xwkj2=xw(iwm,2,kwa)-xw(iwm,2,jwa)
156      xwkj3=xw(iwm,3,kwa)-xw(iwm,3,jwa)
157      rwij2=xwij1**2+xwij2**2+xwij3**2
158      rwkj2=xwkj1**2+xwkj2**2+xwkj3**2
159      rwij2i=one/rwij2
160      rwkj2i=one/rwkj2
161      rwikji=one/sqrt(rwij2*rwkj2)
162      cphi=rwikji*(xwij1*xwkj1+xwij2*xwkj2+xwij3*xwkj3)
163      if(cphi.lt.-one) cphi=-one
164      if(cphi.gt. one) cphi= one
165      phi=acos(cphi)
166      dangle=phi-angle
167      if(iand(iwdt(iwm),mdynam).eq.ldynam)
168     + rang(iwh,2)=rang(iwh,2)+half*for*dangle*dangle
169      sphi=sin(phi)
170      if(abs(sphi).lt.small) sphi=small
171      dfor=for*dangle/sphi
172      dfw1=dfor*(xwkj1*rwikji-xwij1*rwij2i*cphi)
173      dfw2=dfor*(xwkj2*rwikji-xwij2*rwij2i*cphi)
174      dfw3=dfor*(xwkj3*rwikji-xwij3*rwij2i*cphi)
175      fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)+dfw1
176      fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfw1
177      fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)+dfw2
178      fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfw2
179      fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)+dfw3
180      fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfw3
181      dfw1=dfor*(xwij1*rwikji-xwkj1*rwkj2i*cphi)
182      dfw2=dfor*(xwij2*rwikji-xwkj2*rwkj2i*cphi)
183      dfw3=dfor*(xwij3*rwikji-xwkj3*rwkj2i*cphi)
184      fw(iwm,1,kwa,1)=fw(iwm,1,kwa,1)+dfw1
185      fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfw1
186      fw(iwm,2,kwa,1)=fw(iwm,2,kwa,1)+dfw2
187      fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfw2
188      fw(iwm,3,kwa,1)=fw(iwm,3,kwa,1)+dfw3
189      fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfw3
190      if(ip2(8))
191     + ep2(1)=ep2(1)+half*ang(iwh,2,2)*(phi-ang(iwh,1,2))**2
192      if(ip3(8))
193     + ep3(1)=ep3(1)+half*ang(iwh,2,3)*(phi-ang(iwh,1,3))**2
194      if(ith(8)) then
195      deriv(8,1)=deriv(8,1)+
196     + dangle*(half*dangle*ang(iwh,2,4)-for*ang(iwh,1,4))
197      endif
198   50 continue
199      eww(2,1)=eww(2,1)+rang(iwh,2)
200      if(ip2(8)) ep2(1)=ep2(1)-rang(iwh,2)
201      if(ip3(8)) ep3(1)=ep3(1)-rang(iwh,2)
202      if(ipme.ne.0) then
203      iwa=iang(iwh,1)
204      jwa=iang(iwh,3)
205      qij=chg(iwq(iwa),1,iset)*chg(iwq(jwa),1,iset)
206      do 41 iwm=iwfr,iwto
207      rwx1=xw(iwm,1,iwa)-xw(iwm,1,jwa)
208      rwx2=xw(iwm,2,iwa)-xw(iwm,2,jwa)
209      rwx3=xw(iwm,3,iwa)-xw(iwm,3,jwa)
210      rww=sqrt(rwx1**2+rwx2**2+rwx3**2)
211      rwi=one/rww
212      ferfc=one-erfc(ealpha*rww)
213      fderfc=-(ealpha*derfc(ealpha*rww))
214      epmecw=epmecw-ferfc*qij*rwi
215      eww(9,1)=eww(9,1)-ferfc*qij*rwi
216      dfor=-(qij*rwi*rwi*(ferfc*rwi-fderfc))
217      dfw1=dfor*rwx1
218      dfw2=dfor*rwx2
219      dfw3=dfor*rwx3
220      fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfw1
221      fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)+dfw1
222      fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfw2
223      fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)+dfw2
224      fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfw3
225      fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)+dfw3
226      vpmeb(1)=vpmeb(1)+dfw1*rwx1
227      vpmeb(2)=vpmeb(2)+dfw2*rwx1
228      vpmeb(3)=vpmeb(3)+dfw3*rwx1
229      vpmeb(4)=vpmeb(4)+dfw2*rwx2
230      vpmeb(5)=vpmeb(5)+dfw3*rwx2
231      vpmeb(6)=vpmeb(6)+dfw3*rwx3
232   41 continue
233      endif
234   40 continue
235      if(iffld.eq.2) then
236      do 1140 iwh=1,nwh
237      iwa=iang(iwh,1)
238      kwa=iang(iwh,3)
239      bond=ang(iwb,3,iset)
240      for=ang(iwb,4,iset)
241      eub=zero
242      do 150 iwm=iwfr,iwto
243      rwx1=xw(iwm,1,iwa)-xw(iwm,1,kwa)
244      rwx2=xw(iwm,2,iwa)-xw(iwm,2,kwa)
245      rwx3=xw(iwm,3,iwa)-xw(iwm,3,kwa)
246      rww=sqrt(rwx1**2+rwx2**2+rwx3**2)
247      if(rww.lt.tiny) then
248      rwwi=one
249      else
250      rwwi=one/rww
251      endif
252      dbond=rww-bond
253      if(iand(iwdt(iwm),mdynam).eq.ldynam)
254     + eub=eub+half*for*(rww-bond)**2
255      dfor=for*dbond*rwwi
256      dfw1=dfor*rwx1
257      dfw2=dfor*rwx2
258      dfw3=dfor*rwx3
259      fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfw1
260      fw(iwm,1,kwa,1)=fw(iwm,1,kwa,1)+dfw1
261      fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfw2
262      fw(iwm,2,kwa,1)=fw(iwm,2,kwa,1)+dfw2
263      fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfw3
264      fw(iwm,3,kwa,1)=fw(iwm,3,kwa,1)+dfw3
265      if(ip2(8))
266     + ep2(1)=ep2(1)+half*ang(iwh,4,2)*(rww-ang(iwh,3,2))**2
267      if(ip3(8))
268     + ep3(1)=ep3(1)+half*ang(iwh,4,3)*(rww-ang(iwh,3,3))**2
269      if(ith(8)) then
270      deriv(8,1)=deriv(8,1)+
271     + dbond*(half*dbond*ang(iwh,4,4)-for*ang(iwh,3,4))
272      endif
273  150 continue
274      rub(iwh,2)=rub(iwh,2)+eub
275      eww(13,1)=eww(13,1)+eub
276      if(ip2(8)) ep2(1)=ep2(1)-eub
277      if(ip3(8)) ep3(1)=ep3(1)-eub
278 1140 continue
279      endif
280      do 70 iwd=1,nwd
281      iwa=idih(iwd,1)
282      jwa=idih(iwd,2)
283      kwa=idih(iwd,3)
284      lwa=idih(iwd,4)
285      angle=dih(iwd,2,iset)
286      for=dih(iwd,3,iset)
287      rmul=dih(iwd,1,iset)
288      rdih(iwd,2)=zero
289      do 80 iwm=iwfr,iwto
290      xwij1=xw(iwm,1,iwa)-xw(iwm,1,jwa)
291      xwij2=xw(iwm,2,iwa)-xw(iwm,2,jwa)
292      xwij3=xw(iwm,3,iwa)-xw(iwm,3,jwa)
293      xwkj1=xw(iwm,1,kwa)-xw(iwm,1,jwa)
294      xwkj2=xw(iwm,2,kwa)-xw(iwm,2,jwa)
295      xwkj3=xw(iwm,3,kwa)-xw(iwm,3,jwa)
296      xwkl1=xw(iwm,1,kwa)-xw(iwm,1,lwa)
297      xwkl2=xw(iwm,2,kwa)-xw(iwm,2,lwa)
298      xwkl3=xw(iwm,3,kwa)-xw(iwm,3,lwa)
299      xwik1=xwij1-xwkj1
300      xwik2=xwij2-xwkj2
301      xwik3=xwij3-xwkj3
302      xwjl1=xwkl1-xwkj1
303      xwjl2=xwkl2-xwkj2
304      xwjl3=xwkl3-xwkj3
305      xm1=xwij2*xwkj3-xwkj2*xwij3
306      xm2=xwij3*xwkj1-xwkj3*xwij1
307      xm3=xwij1*xwkj2-xwkj1*xwij2
308      xn1=xwkj2*xwkl3-xwkl2*xwkj3
309      xn2=xwkj3*xwkl1-xwkl3*xwkj1
310      xn3=xwkj1*xwkl2-xwkl1*xwkj2
311      rm2i=one/(xm1**2+xm2**2+xm3**2)
312      rn2i=one/(xn1**2+xn2**2+xn3**2)
313      rmni=sqrt(rm2i*rn2i)
314      cphi=(xm1*xn1+xm2*xn2+xm3*xn3)*rmni
315      if(cphi.lt.-one) cphi=-one
316      if(cphi.gt. one) cphi= one
317      phi=acos(cphi)
318      s=xwkj1*(xm2*xn3-xm3*xn2) +xwkj2*(xm3*xn1-xm1*xn3)
319     + +xwkj3*(xm1*xn2-xm2*xn1)
320      if(s.lt.zero) phi=-phi
321      sphi=sin(phi)
322      rpa=rmul*phi-angle
323      if(iand(iwdt(iwm),mdynam).eq.ldynam)
324     + rdih(iwd,2)=rdih(iwd,2)+for*(one+cos(rpa))
325      dfor=(-for)*rmul*sin(rpa)
326      if(ip2(8)) ep2(1)=ep2(1)+
327     + dih(iwd,3,2)*(one+cos(dih(iwd,1,2)*phi-dih(iwd,2,2)))
328      if(ip3(8)) ep3(1)=ep3(1)+
329     + dih(iwd,3,3)*(one+cos(dih(iwd,1,3)*phi-dih(iwd,2,3)))
330      if(abs(sphi).lt.small) sphi=sign(small,sphi)
331      sphii=one/sphi
332      xd1=(-dfor)*sphii*(rmni*xn1-cphi*rm2i*xm1)
333      xe1=(-dfor)*sphii*(rmni*xm1-cphi*rn2i*xn1)
334      xd2=(-dfor)*sphii*(rmni*xn2-cphi*rm2i*xm2)
335      xe2=(-dfor)*sphii*(rmni*xm2-cphi*rn2i*xn2)
336      xd3=(-dfor)*sphii*(rmni*xn3-cphi*rm2i*xm3)
337      xe3=(-dfor)*sphii*(rmni*xm3-cphi*rn2i*xn3)
338      dfwi1=xwkj2*xd3-xwkj3*xd2
339      dfwi2=xwkj3*xd1-xwkj1*xd3
340      dfwi3=xwkj1*xd2-xwkj2*xd1
341      dfwj1=xwik2*xd3-xwik3*xd2-xwkl2*xe3+xwkl3*xe2
342      dfwj2=xwik3*xd1-xwik1*xd3-xwkl3*xe1+xwkl1*xe3
343      dfwj3=xwik1*xd2-xwik2*xd1-xwkl1*xe2+xwkl2*xe1
344      dfwk1=xwjl2*xe3-xwjl3*xe2-xwij2*xd3+xwij3*xd2
345      dfwk2=xwjl3*xe1-xwjl1*xe3-xwij3*xd1+xwij1*xd3
346      dfwk3=xwjl1*xe2-xwjl2*xe1-xwij1*xd2+xwij2*xd1
347      dfwl1=xwkj2*xe3-xwkj3*xe2
348      dfwl2=xwkj3*xe1-xwkj1*xe3
349      dfwl3=xwkj1*xe2-xwkj2*xe1
350      fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfwi1
351      fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfwi2
352      fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfwi3
353      fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfwj1
354      fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfwj2
355      fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfwj3
356      fw(iwm,1,kwa,1)=fw(iwm,1,kwa,1)-dfwk1
357      fw(iwm,2,kwa,1)=fw(iwm,2,kwa,1)-dfwk2
358      fw(iwm,3,kwa,1)=fw(iwm,3,kwa,1)-dfwk3
359      fw(iwm,1,lwa,1)=fw(iwm,1,lwa,1)-dfwl1
360      fw(iwm,2,lwa,1)=fw(iwm,2,lwa,1)-dfwl2
361      fw(iwm,3,lwa,1)=fw(iwm,3,lwa,1)-dfwl3
362      if(ith(9)) then
363      deriv(9,1)=deriv(9,1)+(one+cos(rpa))*dih(iwd,3,4)
364     + -for*sin(rpa)*(phi*dih(iwd,1,4)-dih(iwd,2,4))
365      endif
366   80 continue
367      eww(3,1)=eww(3,1)+rdih(iwd,2)
368      if(ip2(8)) ep2(1)=ep2(1)-rdih(iwd,2)
369      if(ip3(8)) ep3(1)=ep3(1)-rdih(iwd,2)
370   70 continue
371      do 90 iwo=1,nwo
372      iwa=iimp(iwo,1)
373      jwa=iimp(iwo,2)
374      kwa=iimp(iwo,3)
375      lwa=iimp(iwo,4)
376      angle=dimp(iwo,2,iset)
377      for=dimp(iwo,3,iset)
378      rimp(iwo,2)=zero
379      do 100 iwm=iwfr,iwto
380      xwij1=xw(iwm,1,iwa)-xw(iwm,1,jwa)
381      xwij2=xw(iwm,2,iwa)-xw(iwm,2,jwa)
382      xwij3=xw(iwm,3,iwa)-xw(iwm,3,jwa)
383      xwkj1=xw(iwm,1,kwa)-xw(iwm,1,jwa)
384      xwkj2=xw(iwm,2,kwa)-xw(iwm,2,jwa)
385      xwkj3=xw(iwm,3,kwa)-xw(iwm,3,jwa)
386      xwkl1=xw(iwm,1,kwa)-xw(iwm,1,lwa)
387      xwkl2=xw(iwm,2,kwa)-xw(iwm,2,lwa)
388      xwkl3=xw(iwm,3,kwa)-xw(iwm,3,lwa)
389      xwik1=xwij1-xwkj1
390      xwik2=xwij2-xwkj2
391      xwik3=xwij3-xwkj3
392      xwjl1=xwkl1-xwkj1
393      xwjl2=xwkl2-xwkj2
394      xwjl3=xwkl3-xwkj3
395      xm1=xwij2*xwkj3-xwkj2*xwij3
396      xm2=xwij3*xwkj1-xwkj3*xwij1
397      xm3=xwij1*xwkj2-xwkj1*xwij2
398      xn1=xwkj2*xwkl3-xwkl2*xwkj3
399      xn2=xwkj3*xwkl1-xwkl3*xwkj1
400      xn3=xwkj1*xwkl2-xwkl1*xwkj2
401      rm2i=one/(xm1**2+xm2**2+xm3**2)
402      rn2i=one/(xn1**2+xn2**2+xn3**2)
403      rmni=sqrt(rm2i*rn2i)
404      cphi=(xm1*xn1+xm2*xn2+xm3*xn3)
405      if(cphi.lt.-one) cphi=-one
406      if(cphi.gt. one) cphi= one
407      phi=acos(cphi)
408      s=xwkj1*(xm2*xn3-xm3*xn2) +xwkj2*(xm3*xn1-xm1*xn3)
409     + +xwkj3*(xm1*xn2-xm2*xn1)
410      if(s.lt.zero) phi=-phi
411      sphi=sin(phi)
412      dangle=(phi-angle)-nint((phi-angle)/twopi)*twopi
413      dfor=for*dangle
414      if(iand(iwdt(iwm),mdynam).eq.ldynam) rimp(iwo,2)=half*dfor*dangle
415      if(ip2(9)) then
416      danglep=(phi-dimp(iwo,2,2))-nint((phi-dimp(iwo,2,2))/twopi)*twopi
417      ep2(1)=ep2(1)+half*dimp(iwo,3,2)*danglep**2
418      endif
419      if(ip3(9)) then
420      danglep=(phi-dimp(iwo,2,3))-nint((phi-dimp(iwo,2,3))/twopi)*twopi
421      ep3(1)=ep3(1)+half*dimp(iwo,3,3)*danglep**2
422      endif
423      if(abs(sphi).lt.small) sphi=sign(small,sphi)
424      sphii=one/sphi
425      xd1=(-dfor)*sphii*(rmni*xn1-cphi*rm2i*xm1)
426      xe1=(-dfor)*sphii*(rmni*xm1-cphi*rn2i*xn1)
427      xd2=(-dfor)*sphii*(rmni*xn2-cphi*rm2i*xm2)
428      xe2=(-dfor)*sphii*(rmni*xm2-cphi*rn2i*xn2)
429      xd3=(-dfor)*sphii*(rmni*xn3-cphi*rm2i*xm3)
430      xe3=(-dfor)*sphii*(rmni*xm3-cphi*rn2i*xn3)
431      dfwi1=xwkj2*xd3-xwkj3*xd2
432      dfwi2=xwkj3*xd1-xwkj1*xd3
433      dfwi3=xwkj1*xd2-xwkj2*xd1
434      dfwj1=xwik2*xd3-xwik3*xd2-xwkl2*xe3+xwkl3*xe2
435      dfwj2=xwik3*xd1-xwik1*xd3-xwkl3*xe1+xwkl1*xe3
436      dfwj3=xwik1*xd2-xwik2*xd1-xwkl1*xe2+xwkl2*xe1
437      dfwk1=xwjl2*xe3-xwjl3*xe2-xwij2*xd3+xwij3*xd2
438      dfwk2=xwjl3*xe1-xwjl1*xe3-xwij3*xd1+xwij1*xd3
439      dfwk3=xwjl1*xe2-xwjl2*xe1-xwij1*xd2+xwij2*xd1
440      dfwl1=xwkj2*xe3-xwkj3*xe2
441      dfwl2=xwkj3*xe1-xwkj1*xe3
442      dfwl3=xwkj1*xe2-xwkj2*xe1
443      fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfwi1
444      fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfwi2
445      fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfwi3
446      fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfwj1
447      fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfwj2
448      fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfwj3
449      fw(iwm,1,kwa,1)=fw(iwm,1,kwa,1)-dfwk1
450      fw(iwm,2,kwa,1)=fw(iwm,2,kwa,1)-dfwk2
451      fw(iwm,3,kwa,1)=fw(iwm,3,kwa,1)-dfwk3
452      fw(iwm,1,lwa,1)=fw(iwm,1,lwa,1)-dfwl1
453      fw(iwm,2,lwa,1)=fw(iwm,2,lwa,1)-dfwl2
454      fw(iwm,3,lwa,1)=fw(iwm,3,lwa,1)-dfwl3
455      if(ith(10)) then
456      deriv(10,1)=deriv(10,1)+
457     + dangle*(half*dangle*dimp(iwo,3,4)-for*dimp(iwo,2,4))
458      endif
459  100 continue
460      eww(4,1)=eww(4,1)+rimp(iwo,2)
461      if(ip2(9)) ep2(1)=ep2(1)-rimp(iwo,2)
462      if(ip3(9)) ep3(1)=ep3(1)-rimp(iwo,2)
463   90 continue
464      c6p1=zero
465      c12p1=zero
466      c6p2=zero
467      c12p2=zero
468      qip1=zero
469      qjp1=zero
470      qip2=zero
471      qjp2=zero
472      do 110 iwt=1,nwt
473      iwa=idwt(iwt,1)
474      jwa=idwt(iwt,2)
475      c6=vdw(iwatm(iwa),iwatm(jwa),2,iset)
476      c12=vdw(iwatm(iwa),iwatm(jwa),4,iset)
477      if(ip2(2)) then
478      c6p1=vdw(iwatm(iwa),iwatm(jwa),2,2)
479      c12p1=vdw(iwatm(iwa),iwatm(jwa),4,2)
480      endif
481      if(ip3(2)) then
482      c6p2=vdw(iwatm(iwa),iwatm(jwa),2,3)
483      c12p2=vdw(iwatm(iwa),iwatm(jwa),4,3)
484      endif
485      if(ith(2).or.ith(4)) then
486      c6t=vdw(iwatm(iwa),iwatm(jwa),2,4)
487      c12t=vdw(iwatm(iwa),iwatm(jwa),4,4)
488      qit=chg(iwq(iwa),1,4)*q14fac
489      qjt=chg(iwq(jwa),1,4)
490      endif
491      cf6=six*c6
492      cf12=twelve*c12
493      qi=chg(iwq(iwa),1,iset)*q14fac
494      qj=chg(iwq(jwa),1,iset)
495      q=qi*qj
496      if(ip2(4)) then
497      qip1=chg(iwq(iwa),1,2)*q14fac
498      qjp1=chg(iwq(jwa),1,2)
499      qp1=qip1*qjp1
500      endif
501      if(ip3(4)) then
502      qip2=chg(iwq(iwa),1,3)*q14fac
503      qjp2=chg(iwq(jwa),1,3)
504      qp2=qip2*qjp2
505      endif
506      ep2l=zero
507      ep3l=zero
508      ep2q=zero
509      ep3q=zero
510      do 120 iwm=iwfr,iwto
511      rxx=xw(iwm,1,iwa)-xw(iwm,1,jwa)
512      rxy=xw(iwm,2,iwa)-xw(iwm,2,jwa)
513      rxz=xw(iwm,3,iwa)-xw(iwm,3,jwa)
514      r2=rxx*rxx+rxy*rxy+rxz*rxz
515      r2i=one/r2
516      r1i=sqrt(r2i)
517      r6i=r2i*r2i*r2i
518      eterml=(c12*r6i-c6)*r6i
519      etermq=q*r1i
520      if(iand(iwdt(iwm),mdynam).eq.ldynam) eww(5,1)=eww(5,1)+eterml
521      if(iand(iwdt(iwm),mdynam).eq.ldynam) eww(6,1)=eww(6,1)+etermq
522      if(ip2(2)) ep2l=ep2l-eterml+(c12p1*r6i-c6p1)*r6i
523      if(ip3(2)) ep3l=ep3l-eterml+(c12p2*r6i-c6p2)*r6i
524      if(ip2(4)) ep2q=ep2q-etermq+qp1*r1i
525      if(ip3(4)) ep3q=ep3q-etermq+qp2*r1i
526      dfw=((cf12*r6i-cf6)*r6i+q*r1i)*r2i
527      fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)+dfw*rxx
528      fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)+dfw*rxy
529      fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)+dfw*rxz
530      fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfw*rxx
531      fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfw*rxy
532      fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfw*rxz
533      if(ith(2)) then
534      deriv(2,1)=deriv(2,1)+(c12t*r6i-c6t)*r6i
535      endif
536      if(ith(4)) then
537      deriv(4,1)=deriv(4,1)+(qi*qjt+qj*qit)*r1i
538      endif
539  120 continue
540      ep2(1)=ep2(1)+ep2l+ep2q
541      ep3(1)=ep3(1)+ep3l+ep3q
542      if(ipme.ne.0) then
543      qij=(one-q14fac)*chg(iwq(iwa),1,iset)*chg(iwq(jwa),1,iset)
544      do 111 iwm=iwfr,iwto
545      rwx1=xw(iwm,1,iwa)-xw(iwm,1,jwa)
546      rwx2=xw(iwm,2,iwa)-xw(iwm,2,jwa)
547      rwx3=xw(iwm,3,iwa)-xw(iwm,3,jwa)
548      rww=sqrt(rwx1**2+rwx2**2+rwx3**2)
549      rwi=one/rww
550      ferfc=one-erfc(ealpha*rww)
551      fderfc=-(ealpha*derfc(ealpha*rww))
552      epmecw=epmecw-ferfc*qij*rwi
553      eww(6,1)=eww(6,1)-ferfc*qij*rwi
554      dfor=-(qij*rwi*rwi*(ferfc*rwi-fderfc))
555      dfw1=dfor*rwx1
556      dfw2=dfor*rwx2
557      dfw3=dfor*rwx3
558      fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfw1
559      fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)+dfw1
560      fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfw2
561      fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)+dfw2
562      fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfw3
563      fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)+dfw3
564      vpmeb(1)=vpmeb(1)+dfw1*rwx1
565      vpmeb(2)=vpmeb(2)+dfw2*rwx1
566      vpmeb(3)=vpmeb(3)+dfw3*rwx1
567      vpmeb(4)=vpmeb(4)+dfw2*rwx2
568      vpmeb(5)=vpmeb(5)+dfw3*rwx2
569      vpmeb(6)=vpmeb(6)+dfw3*rwx3
570  111 continue
571      endif
572  110 continue
573      do 130 iwn=1,nwn
574      iwa=idwn(iwn,1)
575      jwa=idwn(iwn,2)
576      c6=vdw(iwatm(iwa),iwatm(jwa),1,iset)
577      c12=vdw(iwatm(iwa),iwatm(jwa),3,iset)
578      if(ip2(2)) then
579      c6p1=vdw(iwatm(iwa),iwatm(jwa),1,2)
580      c12p1=vdw(iwatm(iwa),iwatm(jwa),3,2)
581      endif
582      if(ip3(2)) then
583      c6p2=vdw(iwatm(iwa),iwatm(jwa),1,3)
584      c12p2=vdw(iwatm(iwa),iwatm(jwa),3,3)
585      endif
586      if(ith(2).or.ith(4)) then
587      c6t=vdw(iwatm(iwa),iwatm(jwa),1,4)
588      c12t=vdw(iwatm(iwa),iwatm(jwa),3,4)
589      qit=chg(iwq(iwa),1,4)
590      qjt=chg(iwq(jwa),1,4)
591      endif
592      cf6=six*c6
593      cf12=twelve*c12
594      qi=chg(iwq(iwa),1,iset)
595      qj=chg(iwq(jwa),1,iset)
596      q=qi*qj
597      if(ip2(4)) then
598      qip1=chg(iwq(iwa),1,2)
599      qjp1=chg(iwq(jwa),1,2)
600      qp1=qip1*qjp1
601      endif
602      if(ip3(4)) then
603      qip2=chg(iwq(iwa),1,3)
604      qjp2=chg(iwq(jwa),1,3)
605      qp2=qip2*qjp2
606      endif
607      ep2l=zero
608      ep3l=zero
609      ep2q=zero
610      ep3q=zero
611      do 140 iwm=iwfr,iwto
612      rxx=xw(iwm,1,iwa)-xw(iwm,1,jwa)
613      rxy=xw(iwm,2,iwa)-xw(iwm,2,jwa)
614      rxz=xw(iwm,3,iwa)-xw(iwm,3,jwa)
615      r2=rxx*rxx+rxy*rxy+rxz*rxz
616      r2i=one/r2
617      r1i=sqrt(r2i)
618      r6i=r2i*r2i*r2i
619      ferfc=one
620      fderfc=zero
621      if(ipme.ne.0) then
622      ferfc=erfc(ealpha/r1i)
623      fderfc=ealpha+derfc(ealpha/r1i)
624      endif
625      eterml=(c12*r6i-c6)*r6i
626      etermq=ferfc*q*r1i
627      if(iand(iwdt(iwm),mdynam).eq.ldynam) then
628      eww(5,1)=eww(5,1)+eterml
629      eww(6,1)=eww(6,1)+etermq
630      endif
631      if(ip2(2)) ep2l=ep2l-eterml+(c12p1*r6i-c6p1)*r6i
632      if(ip3(2)) ep3l=ep3l-eterml+(c12p2*r6i-c6p2)*r6i
633      if(ip2(4)) ep2q=ep2q-etermq+qp1*r1i
634      if(ip3(4)) ep3q=ep3q-etermq+qp2*r1i
635      dfw=((cf12*r6i-cf6)*r6i+q*(ferfc*r1i-fderfc))*r2i
636      fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)+dfw*rxx
637      fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)+dfw*rxy
638      fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)+dfw*rxz
639      fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfw*rxx
640      fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfw*rxy
641      fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfw*rxz
642      if(ith(2)) deriv(2,1)=deriv(2,1)+(c12t*r6i-c6t)*r6i
643      if(ith(4)) deriv(4,1)=deriv(4,1)+(qi*qjt+qj*qit)*r1i
644  140 continue
645      ep2(1)=ep2(1)+ep2l+ep2q
646      ep3(1)=ep3(1)+ep3l+ep3q
647  130 continue
648c
649      return
650      end
651