1      subroutine argos_cafe_fpww(xw,xwm,fw,pw,pwp,idt,iwfrom,
2     + nwloc,lpbc,eww,
3     + vdw,chg,iwatm,iwq,lwwndx,lwwjpt,lwwin,lwwj,
4     + xi,xj,rwx,rwi1,rwi2,rwi6,rwc,
5     + f,fi,fj,facu,pl,pj)
6c
7c $Id$
8c
9      implicit none
10c
11#include "argos_cafe_common.fh"
12#include "argos_cafe_funcs_dec.fh"
13#include "bitops_decls.fh"
14c
15      real*8 xw(mwm,3,mwa),xwm(mwm,3),fw(mwm,3,mwa,2),eww(mpe,2)
16      integer idt(mwm)
17      integer iwfrom,nwloc
18      logical lpbc
19c
20      real*8 vdw(mat,mat,map,mset),chg(mqt,mqp,mset)
21      integer iwatm(mwa),iwq(mwa)
22c
23      real*8 xi(mscr,3,mwa),xj(mscr,3,mwa),rwx(mscr,3)
24      real*8 rwi1(mscr),rwi2(mscr),rwi6(mscr),rwc(mscr,3)
25      real*8 f(mscr),fi(mscr,3,mwa),fj(mscr,3,mwa)
26c
27      real*8 facu(mscr)
28c     real*8 rdf(mgl,mgr)
29c
30      integer lwwj(*)
31      integer lwwndx(0:mwm,2),lwwjpt(nwloc,2),lwwin(nwloc,2)
32c
33
34      real*8 pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2)
35      real*8 pl(mscr,3,mwa),pj(mscr,3,mwa)
36c      integer nax2,ipset
37      real*8 qai,qaj,pai,paj,pix,piy,piz,pjx,pjy,pjz
38      real*8 ri3,rmi,rmj,fri,fmi,fmj,rmm,qfaci
39      real*8 rx,ry,rz,ri1,ri2,ewwpsm,etermp
40      real*8 ewwqsm
41c
42      integer iwfr,ipww,number,iwm,iwpm,nax
43      integer iwmn,lwwptr,iwa,iax,jwa,iptr,jptr,iwpj
44      real*8 ewwl6,ewwl12,q
45      real*8 c64,c124,qi,qj,qi4,qj4,dercon
46      real*8 c6p,c12p,qp,ep2tmp,ep3tmp
47      real*8 c6,cf6,c12,cf12
48c
49      integer nwwlen(2)
50      real*8 eterml,etermq
51c
52#include "argos_cafe_funcs_sfn.fh"
53#include "bitops_funcs.fh"
54c
55cx new stuff begin
56c
57cx new stuff end
58c
59c     calculation of solvent-solvent intermolecular energies and forces
60c
61c     subtract 1 from first molecule index for use as offset
62c
63      iwfr=iwfrom-1
64c
65c     loop over short and long range parts
66c
67      do 1 ipww=1,lpww
68c
69c     Evaluate the outer index array
70c
71      nwwlen(ipww)=0
72      lwwndx(0,ipww)=0
73      number=0
74      do 2 iwm=1,nwloc
75      if(number+lwwin(iwm,ipww).gt.mscr) then
76      nwwlen(ipww)=nwwlen(ipww)+1
77      lwwndx(nwwlen(ipww),ipww)=iwm-1
78      number=0
79      endif
80      number=number+lwwin(iwm,ipww)
81    2 continue
82      if(number.gt.0) then
83      nwwlen(ipww)=nwwlen(ipww)+1
84      lwwndx(nwwlen(ipww),ipww)=nwloc
85      endif
86c
87c     loop over number of cycles to complete pairlist
88c
89      do 3 iwpm=1,nwwlen(ipww)
90      nax=0
91c
92c     collect coordinates into workarrays
93c
94      do 4 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww)
95      iwpj=lwwjpt(iwm,ipww)-1
96      do 5 iwmn=1,lwwin(iwm,ipww)
97      lwwptr=lwwj(iwpj+iwmn)
98      rwc(nax+iwmn,1)=xwm(iwfr+iwm,1)-xwm(lwwptr,1)
99      rwc(nax+iwmn,2)=xwm(iwfr+iwm,2)-xwm(lwwptr,2)
100      rwc(nax+iwmn,3)=xwm(iwfr+iwm,3)-xwm(lwwptr,3)
101      facu(nax+iwmn)=one
102c      if( (iand(idt(iwm),mdynam).eq.ldynam.and.
103c     + iand(idt(lwwptr),mdynam).ne.ldynam).or.
104c     + (iand(idt(iwm),mdynam).ne.ldynam.and.
105c     + iand(idt(lwwptr),mdynam).eq.ldynam) ) facu(nax+iwmn)=half
106      if(iand(idt(iwm),mdynam).ne.ldynam.and.
107     + iand(idt(lwwptr),mdynam).ne.ldynam) facu(nax+iwmn)=zero
108      if(includ.eq.1) facu(nax+iwmn)=one
109    5 continue
110c
111      do 6 iwa=1,mwa
112      do 7 iwmn=1,lwwin(iwm,ipww)
113      lwwptr=lwwj(iwpj+iwmn)
114      xi(nax+iwmn,1,iwa)=xw(iwfr+iwm,1,iwa)
115      xi(nax+iwmn,2,iwa)=xw(iwfr+iwm,2,iwa)
116      xi(nax+iwmn,3,iwa)=xw(iwfr+iwm,3,iwa)
117      xj(nax+iwmn,1,iwa)=xw(lwwptr,1,iwa)
118      xj(nax+iwmn,2,iwa)=xw(lwwptr,2,iwa)
119      xj(nax+iwmn,3,iwa)=xw(lwwptr,3,iwa)
120      pl(nax+iwmn,1,iwa)=pw(iwfr+iwm,1,iwa,1)
121      pl(nax+iwmn,2,iwa)=pw(iwfr+iwm,2,iwa,1)
122      pl(nax+iwmn,3,iwa)=pw(iwfr+iwm,3,iwa,1)
123      pj(nax+iwmn,1,iwa)=pw(lwwptr,1,iwa,1)
124      pj(nax+iwmn,2,iwa)=pw(lwwptr,2,iwa,1)
125      pj(nax+iwmn,3,iwa)=pw(lwwptr,3,iwa,1)
126    7 continue
127    6 continue
128      if(lpbc) then
129      call argos_cafe_pbc(0,rwc,mscr,rwx,mscr,nax,1,lwwin(iwm,ipww))
130      do 8 iwmn=1,lwwin(iwm,ipww)
131      rwc(nax+iwmn,1)=rwc(nax+iwmn,1)-rwx(iwmn,1)
132      rwc(nax+iwmn,2)=rwc(nax+iwmn,2)-rwx(iwmn,2)
133      rwc(nax+iwmn,3)=rwc(nax+iwmn,3)-rwx(iwmn,3)
134    8 continue
135      do 9 iwa=1,mwa
136      do 10 iwmn=1,lwwin(iwm,ipww)
137      lwwptr=lwwj(iwpj+iwmn)
138      xj(nax+iwmn,1,iwa)=xj(nax+iwmn,1,iwa)+rwx(iwmn,1)
139      xj(nax+iwmn,2,iwa)=xj(nax+iwmn,2,iwa)+rwx(iwmn,2)
140      xj(nax+iwmn,3,iwa)=xj(nax+iwmn,3,iwa)+rwx(iwmn,3)
141   10 continue
142    9 continue
143      endif
144c
145      nax=nax+lwwin(iwm,ipww)
146    4 continue
147c
148c     initializations
149c
150c      if(npener.ne.0) then
151c      do 12 iax=1,nax
152c      u(iax)=zero
153c   12 continue
154c      endif
155c
156c     loops over number of atoms in a solvent molecule
157c
158      qfaci=one/qfac
159      do 13 iwa=1,mwa
160      qi=chg(iwq(iwa),1,iset)
161      pai=chg(iwq(iwa),2,iset)
162      qai=qfaci*qi
163      do 14 jwa=1,mwa
164      qj=chg(iwq(jwa),1,iset)
165      q=qi*qj
166      paj=chg(iwq(jwa),2,iset)
167      qaj=qfaci*qj
168c
169      do 15 iax=1,nax
170      f(iax)=zero
171      rwx(iax,1)=xi(iax,1,iwa)-xj(iax,1,jwa)
172      rwx(iax,2)=xi(iax,2,iwa)-xj(iax,2,jwa)
173      rwx(iax,3)=xi(iax,3,iwa)-xj(iax,3,jwa)
174      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
175      rwi1(iax)=sqrt(rwi2(iax))
176   15 continue
177c
178c
179c     van der Waals contribution
180c     --------------------------
181c
182      iptr=iwatm(iwa)
183      jptr=iwatm(jwa)
184      c6=vdw(iptr,jptr,1,iset)
185      cf6=six*c6
186      c12=vdw(iptr,jptr,3,iset)
187      cf12=twelve*c12
188c
189      eterml=zero
190      if(c6.ne.zero.or.c12.ne.zero) then
191      ewwl6=zero
192      ewwl12=zero
193      do 20 iax=1,nax
194      rwi6(iax)=rwi2(iax)*rwi2(iax)*rwi2(iax)
195      ewwl6=ewwl6+facu(iax)*rwi6(iax)
196      ewwl12=ewwl12+facu(iax)*rwi6(iax)*rwi6(iax)
197      f(iax)=f(iax)+(cf12*rwi6(iax)-cf6)*rwi6(iax)*rwi2(iax)
198   20 continue
199      eterml=c12*ewwl12-c6*ewwl6
200      eww(7,ipww)=eww(7,ipww)+eterml
201      endif
202c
203c
204c     force vectors
205c     -------------
206c
207      if(iwa.eq.1) then
208      do 22 iax=1,nax
209      fj(iax,1,jwa)=(-f(iax))*rwx(iax,1)
210      fj(iax,2,jwa)=(-f(iax))*rwx(iax,2)
211      fj(iax,3,jwa)=(-f(iax))*rwx(iax,3)
212   22 continue
213      else
214      do 23 iax=1,nax
215      fj(iax,1,jwa)=fj(iax,1,jwa)-f(iax)*rwx(iax,1)
216      fj(iax,2,jwa)=fj(iax,2,jwa)-f(iax)*rwx(iax,2)
217      fj(iax,3,jwa)=fj(iax,3,jwa)-f(iax)*rwx(iax,3)
218   23 continue
219      endif
220c
221      if(jwa.eq.1) then
222      do 24 iax=1,nax
223      fi(iax,1,iwa)=f(iax)*rwx(iax,1)
224      fi(iax,2,iwa)=f(iax)*rwx(iax,2)
225      fi(iax,3,iwa)=f(iax)*rwx(iax,3)
226   24 continue
227      else
228      do 25 iax=1,nax
229      fi(iax,1,iwa)=fi(iax,1,iwa)+f(iax)*rwx(iax,1)
230      fi(iax,2,iwa)=fi(iax,2,iwa)+f(iax)*rwx(iax,2)
231      fi(iax,3,iwa)=fi(iax,3,iwa)+f(iax)*rwx(iax,3)
232   25 continue
233      endif
234      do 26 iax=1,nax
235      zw(1,1,ipww)=zw(1,1,ipww)-f(iax)*rwx(iax,1)*rwc(iax,1)
236      zw(2,1,ipww)=zw(2,1,ipww)-f(iax)*rwx(iax,1)*rwc(iax,2)
237      zw(3,1,ipww)=zw(3,1,ipww)-f(iax)*rwx(iax,1)*rwc(iax,3)
238      zw(1,2,ipww)=zw(1,2,ipww)-f(iax)*rwx(iax,2)*rwc(iax,1)
239      zw(2,2,ipww)=zw(2,2,ipww)-f(iax)*rwx(iax,2)*rwc(iax,2)
240      zw(3,2,ipww)=zw(3,2,ipww)-f(iax)*rwx(iax,2)*rwc(iax,3)
241      zw(1,3,ipww)=zw(1,3,ipww)-f(iax)*rwx(iax,3)*rwc(iax,1)
242      zw(2,3,ipww)=zw(2,3,ipww)-f(iax)*rwx(iax,3)*rwc(iax,2)
243      zw(3,3,ipww)=zw(3,3,ipww)-f(iax)*rwx(iax,3)*rwc(iax,3)
244   26 continue
245c
246c
247c     electrostatic and polarization contribution
248c     -------------------------------------------
249c
250      ewwqsm=zero
251      ewwpsm=zero
252      do 117 iax=1,nax
253      pix=pai*pl(iax,1,iwa)
254      piy=pai*pl(iax,2,iwa)
255      piz=pai*pl(iax,3,iwa)
256      pjx=paj*pj(iax,1,jwa)
257      pjy=paj*pj(iax,2,jwa)
258      pjz=paj*pj(iax,3,jwa)
259      rx=-rwx(iax,1)
260      ry=-rwx(iax,2)
261      rz=-rwx(iax,3)
262      ri1=rwi1(iax)
263      ri2=rwi2(iax)
264      ri3=qfac*qfac*ri1*ri2
265      rmi=three*(rx*pix+ry*piy+rz*piz)*ri2
266      rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2
267      if(ipolt.eq.1) then
268      fri=((-qai)*qaj+qai*rmj-qaj*rmi)*ri3
269      fmi=(qaj)*ri3
270      fmj=(-qai)*ri3
271      else
272      rmm=three*(pix*pjx+piy*pjy+piz*pjz)*ri2
273      fri=((-qai)*qaj+qai*rmj-qaj*rmi+5.0*rmi*rmj/three-rmm)*ri3
274      fmi=(qaj-rmj)*ri3
275      fmj=((-qai)-rmi)*ri3
276      endif
277      fi(iax,1,iwa)=fi(iax,1,iwa)+fri*rx+fmi*pix+fmj*pjx
278      fi(iax,2,iwa)=fi(iax,2,iwa)+fri*ry+fmi*piy+fmj*pjy
279      fi(iax,3,iwa)=fi(iax,3,iwa)+fri*rz+fmi*piz+fmj*pjz
280      fj(iax,1,jwa)=fj(iax,1,jwa)-(fri*rx+fmi*pix+fmj*pjx)
281      fj(iax,2,jwa)=fj(iax,2,jwa)-(fri*ry+fmi*piy+fmj*pjy)
282      fj(iax,3,jwa)=fj(iax,3,jwa)-(fri*rz+fmi*piz+fmj*pjz)
283      zw(1,1,ipww)=zw(1,1,ipww)-(fri*rx+fmi*pix+fmj*pjx)*rwc(iax,1)
284      zw(2,1,ipww)=zw(2,1,ipww)-(fri*rx+fmi*pix+fmj*pjx)*rwc(iax,2)
285      zw(3,1,ipww)=zw(3,1,ipww)-(fri*rx+fmi*pix+fmj*pjx)*rwc(iax,3)
286      zw(1,2,ipww)=zw(1,2,ipww)-(fri*ry+fmi*piy+fmj*pjy)*rwc(iax,1)
287      zw(2,2,ipww)=zw(2,2,ipww)-(fri*ry+fmi*piy+fmj*pjy)*rwc(iax,2)
288      zw(3,2,ipww)=zw(3,2,ipww)-(fri*ry+fmi*piy+fmj*pjy)*rwc(iax,3)
289      zw(1,3,ipww)=zw(1,3,ipww)-(fri*rz+fmi*piz+fmj*pjz)*rwc(iax,1)
290      zw(2,3,ipww)=zw(2,3,ipww)-(fri*rz+fmi*piz+fmj*pjz)*rwc(iax,2)
291      zw(3,3,ipww)=zw(3,3,ipww)-(fri*rz+fmi*piz+fmj*pjz)*rwc(iax,3)
292      ewwpsm=ewwpsm+facu(iax)*(qai*rmj-qaj*rmi)*ri1
293      ewwqsm=ewwqsm+facu(iax)*ri1
294  117 continue
295      etermp=-qfac*qfac*ewwpsm/three
296      eww(8,ipww)=eww(8,ipww)+etermp
297      etermq=q*ewwqsm
298      eww(8,ipww)=eww(8,ipww)+etermq
299c
300c     Radial distribution functions
301c
302c      if(ifstep-1.eq.((ifstep-1)/nfrdf)*nfrdf .and. ngrww.gt.0) then
303c      do 27 igc=1,ngc
304c      if(ngt(igc).eq.1) then
305c      if(iagc(igc).eq.iwa .and. jagc(igc).eq.jwa) then
306c      igr=igrc(igc)
307c      do 28 iax=1,nax
308c      indx=int(one/(rwi1(iax)*drdf))
309c      if(indx.le.ngl) rdf(indx,igr)=rdf(indx,igr)+rdfvol
310c   28 continue
311c      endif
312c      endif
313c   27 continue
314c      endif
315c
316c     Thermodynamic integration
317c
318      if(ithint) then
319      if(ith(2)) then
320      c64=vdw(iwatm(iwa),iwatm(jwa),1,4)
321      c124=vdw(iwatm(iwa),iwatm(jwa),3,4)
322      ewwl6=zero
323      ewwl12=zero
324      do 29 iax=1,nax
325      ewwl6=ewwl6+facu(iax)*rwi6(iax)
326      ewwl12=ewwl12+facu(iax)*rwi6(iax)*rwi6(iax)
327   29 continue
328      deriv(2,ipww)=deriv(2,ipww)+c124*ewwl12-c64*ewwl6
329      endif
330      if(ith(4)) then
331      qi=chg(iwq(iwa),1,iset)
332      qj=chg(iwq(jwa),1,iset)
333      qi4=chg(iwq(iwa),1,4)
334      qj4=chg(iwq(jwa),1,4)
335      dercon=zero
336      if(ipme.eq.0) then
337      do 30 iax=1,nax
338      dercon=dercon+rwi1(iax)
339   30 continue
340      else
341      do 130 iax=1,nax
342      dercon=dercon+rwi1(iax)
343  130 continue
344      endif
345      deriv(4,ipww)=deriv(4,ipww)+(qi*qj4+qj*qi4)*dercon
346      if(ireact.ne.0) then
347      dercon=zero
348      do 31 iax=1,nax
349      dercon=dercon+one/rwi2(iax)
350   31 continue
351      deriv(4,ipww)=deriv(4,ipww)+(qi*qj4+qj*qi4)*rffww*dercon
352      endif
353      endif
354      endif
355c
356c     Thermodynamic perturbation 1
357c
358      if(ipert2) then
359      if(ip2(2)) then
360      c6p=vdw(iwatm(iwa),iwatm(jwa),1,2)
361      c12p=vdw(iwatm(iwa),iwatm(jwa),3,2)
362      do 32 iax=1,nax
363      ep2(ipww)=ep2(ipww)+facu(iax)*(c12p*rwi6(iax)-c6p)*rwi6(iax)
364   32 continue
365      ep2(ipww)=ep2(ipww)-eterml
366      endif
367      if(ip2(4).or.ip2(5)) then
368      qp=chg(iwq(iwa),1,2)*chg(iwq(jwa),1,2)
369      ep2tmp=zero
370      do 33 iax=1,nax
371      rwx(iax,1)=xi(iax,1,iwa)-xj(iax,1,jwa)
372      rwx(iax,2)=xi(iax,2,iwa)-xj(iax,2,jwa)
373      rwx(iax,3)=xi(iax,3,iwa)-xj(iax,3,jwa)
374      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
375      rwi1(iax)=sqrt(rwi2(iax))
376      if(ipme.eq.0) then
377      ep2tmp=ep2tmp+facu(iax)*rwi1(iax)
378      else
379      ep2tmp=ep2tmp+facu(iax)*erfc(ealpha/rwi1(iax))*rwi1(iax)
380      endif
381   33 continue
382      ep2(ipww)=ep2(ipww)+qp*ep2tmp-etermq
383      if(ireact.ne.0) then
384      ep2tmp=zero
385      do 34 iax=1,nax
386      ep2tmp=ep2tmp+facu(iax)/rwi2(iax)
387   34 continue
388      ep2(ipww)=ep2(ipww)+qp*rffww*ep2tmp
389      endif
390      endif
391      endif
392c
393c     Thermodynamic perturbation 2
394c
395      if(ipert3) then
396      if(ip3(2)) then
397      c6p=vdw(iwatm(iwa),iwatm(jwa),1,3)
398      c12p=vdw(iwatm(iwa),iwatm(jwa),3,3)
399      do 35 iax=1,nax
400      ep3(ipww)=ep3(ipww)+facu(iax)*(c12p*rwi6(iax)-c6p)*rwi6(iax)
401   35 continue
402      ep3(ipww)=ep3(ipww)-eterml
403      endif
404      if(ip2(4).or.ip2(5)) then
405      qp=chg(iwatm(iwa),1,3)*chg(iwatm(jwa),1,3)
406      ep3tmp=zero
407      do 36 iax=1,nax
408      rwx(iax,1)=xi(iax,1,iwa)-xj(iax,1,jwa)
409      rwx(iax,2)=xi(iax,2,iwa)-xj(iax,2,jwa)
410      rwx(iax,3)=xi(iax,3,iwa)-xj(iax,3,jwa)
411      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
412      rwi1(iax)=sqrt(rwi2(iax))
413      if(ipme.eq.0) then
414      ep3tmp=ep3tmp+facu(iax)*rwi1(iax)
415      else
416      ep3tmp=ep3tmp+facu(iax)*erfc(ealpha/rwi1(iax))*rwi1(iax)
417      endif
418   36 continue
419      ep3(ipww)=ep3(ipww)+qp*ep3tmp-etermq
420      if(ireact.ne.0) then
421      ep3tmp=zero
422      do 37 iax=1,nax
423      ep3tmp=ep3tmp+facu(iax)/rwi2(iax)
424   37 continue
425      ep3(ipww)=ep3(ipww)+qp*rffww*ep3tmp
426      endif
427      endif
428      endif
429   14 continue
430   13 continue
431c
432c     Update force arrays
433c
434      iax=0
435      do 38 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww)
436      iwpj=lwwjpt(iwm,ipww)-1
437      do 39 iwa=1,mwa
438      do 40 iwmn=1,lwwin(iwm,ipww)
439      lwwptr=lwwj(iwpj+iwmn)
440      fw(iwfr+iwm,1,iwa,ipww)=fw(iwfr+iwm,1,iwa,ipww)+fi(iax+iwmn,1,iwa)
441      fw(iwfr+iwm,2,iwa,ipww)=fw(iwfr+iwm,2,iwa,ipww)+fi(iax+iwmn,2,iwa)
442      fw(iwfr+iwm,3,iwa,ipww)=fw(iwfr+iwm,3,iwa,ipww)+fi(iax+iwmn,3,iwa)
443      fw(lwwptr,1,iwa,ipww)=fw(lwwptr,1,iwa,ipww)+fj(iax+iwmn,1,iwa)
444      fw(lwwptr,2,iwa,ipww)=fw(lwwptr,2,iwa,ipww)+fj(iax+iwmn,2,iwa)
445      fw(lwwptr,3,iwa,ipww)=fw(lwwptr,3,iwa,ipww)+fj(iax+iwmn,3,iwa)
446   40 continue
447   39 continue
448c
449c     update energy arrays if appropriate print option was set
450c
451c      if(npener.ne.0) then
452c      do 41 iwmn=1,lwwin(iwm,ipww)
453c      lwwptr=lwwj(iwpj+iwmn)
454c      uwmw(iwfr+iwm)=uwmw(iwfr+iwm)+u(iax+iwmn)
455c      uwmw(lwwptr)=uwmw(lwwptr)+u(iax+iwmn)
456c   41 continue
457c      endif
458c
459      iax=iax+lwwin(iwm,ipww)
460   38 continue
461    3 continue
462c
463    1 continue
464c
465      return
466      end
467