1      subroutine argos_cafe_dsw(xs,xsm,ps,psp,isdt,ismf,isml,
2     + isq1,isfrom,nums,
3     + xw,xwm,pw,pwp,iwq,lpbc,chg,lswndx,lswjpt,lswin,lswj,
4     + rwc,xi,xj,rwx,pl,pj,fi,fj,isal)
5c
6c $Id$
7c
8      implicit none
9c
10#include "argos_cafe_common.fh"
11c
12      real*8 xs(msa,3),xsm(msm,3),ps(msa,3,2),psp(msa,3,2,2)
13      integer isdt(msa),ismf(msa),isml(msa),isq1(msa)
14      integer isfrom,nums
15      real*8 xw(mwm,3,mwa),xwm(mwm,3),pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2)
16      integer iwq(mwa)
17      logical lpbc
18      real*8 chg(mqt,mqp,mset)
19      integer lswndx(0:msa,2),lswjpt(nums,2),lswin(nums,2),lswj(*)
20      real*8 rwc(mscr,3),xi(mscr,3),xj(mscr,3,mwa),rwx(mscr,3)
21      real*8 pl(mscr,3),pj(mscr,3,mwa),fi(mscr,3),fj(mscr,3,mwa)
22      integer isal(mscr)
23c
24      integer ispj,ism,lswptr,ipset,nswlen(2)
25      integer isfr,ipsw,number,isa,jwa,ismn,ispm,iax,nax,nax2
26      real*8 qai,qaj,pai,paj,rx,ry,rz,pix,piy,piz,pjx,pjy,pjz
27      real*8 ri2,ri3,rmi,rmj
28#include "bitops.fh"
29c
30      real*8 qfaci
31c
32      qfaci=one/qfac
33c
34c     this subroutine evaluates the solute-solvent forces for nums
35c     solute atoms starting from isfrom. the interacting solvent
36c     molecules are determined from the pairlist.
37c
38      isfr=isfrom-1
39c
40c     loop over short and long range pairs
41c
42      do 1 ipsw=1,npsw
43c
44c     evaluate outer index array
45c
46      nswlen(ipsw)=0
47      lswndx(0,ipsw)=0
48      number=0
49      do 2 isa=1,nums
50      if(number+lswin(isa,ipsw).gt.mscr .or.
51     + (ismf(isfr+isa).ne.ismf(isfr+isa-1).and.number.gt.0)) then
52      nswlen(ipsw)=nswlen(ipsw)+1
53      lswndx(nswlen(ipsw),ipsw)=isa-1
54      number=0
55      endif
56      number=number+lswin(isa,ipsw)
57    2 continue
58      if(number.gt.0) then
59      nswlen(ipsw)=nswlen(ipsw)+1
60      lswndx(nswlen(ipsw),ipsw)=nums
61      endif
62c
63c     loop over number of cycles to complete pairlist
64c
65      do 3 ispm=1,nswlen(ipsw)
66      nax=0
67c
68c     vacuo conditions
69c
70c      if(npbtyp.eq.0) then
71c      if(.not.lpbc) then
72      do 4 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
73      ispj=lswjpt(isa,ipsw)-1
74      ism=isml(isfr+isa)
75c
76c     collect center of mass distance vectors
77c
78      if(lpbc.or.ism.eq.0) then
79      do 6 ismn=1,lswin(isa,ipsw)
80      lswptr=lswj(ispj+ismn)
81      rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1)
82      rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2)
83      rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3)
84    6 continue
85      if(lpbc)
86     + call argos_cafe_pbc(0,rwc,mscr,rwx,mscr,nax,1,lswin(isa,ipsw))
87      else
88      do 5 ismn=1,lswin(isa,ipsw)
89      lswptr=lswj(ispj+ismn)
90      rwc(nax+ismn,1)=xsm(ism,1)-xwm(lswptr,1)
91      rwc(nax+ismn,2)=xsm(ism,2)-xwm(lswptr,2)
92      rwc(nax+ismn,3)=xsm(ism,3)-xwm(lswptr,3)
93    5 continue
94      endif
95c
96c     collect solute coordinates and atomic polarization fields
97c
98c      if(iand(isdt(isfr+isa),mdynam).eq.ldynam) then
99      do 7 ismn=1,lswin(isa,ipsw)
100      lswptr=lswj(ispj+ismn)
101      xi(nax+ismn,1)=xs(isfr+isa,1)
102      xi(nax+ismn,2)=xs(isfr+isa,2)
103      xi(nax+ismn,3)=xs(isfr+isa,3)
104      pl(nax+ismn,1)=ps(isfr+isa,1,2)
105      pl(nax+ismn,2)=ps(isfr+isa,2,2)
106      pl(nax+ismn,3)=ps(isfr+isa,3,2)
107      isal(nax+ismn)=isfr+isa
108    7 continue
109c      else
110c      do 8 ismn=1,lswin(isa,ipsw)
111c      lswptr=lswj(ispj+ismn)
112c      xi(nax+ismn,1)=xs(isfr+isa,1)
113c      xi(nax+ismn,2)=xs(isfr+isa,2)
114c      xi(nax+ismn,3)=xs(isfr+isa,3)
115c      pl(nax+ismn,1)=ps(isfr+isa,1,2)
116c      pl(nax+ismn,2)=ps(isfr+isa,2,2)
117c      pl(nax+ismn,3)=ps(isfr+isa,3,2)
118c      isal(nax+ismn)=isfr+isa
119c    8 continue
120c      endif
121c
122c     collect solvent coordinates and atomic polarization fields
123c
124      do 8 jwa=1,mwa
125      if(lpbc) then
126      do 9 ismn=1,lswin(isa,ipsw)
127      lswptr=lswj(ispj+ismn)
128      xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa)+rwx(ismn,1)
129      xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa)+rwx(ismn,2)
130      xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa)+rwx(ismn,3)
131      pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2)
132      pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2)
133      pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2)
134    9 continue
135      else
136      do 10 ismn=1,lswin(isa,ipsw)
137      lswptr=lswj(ispj+ismn)
138      xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa)
139      xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa)
140      xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa)
141      pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2)
142      pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2)
143      pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2)
144   10 continue
145      endif
146    8 continue
147      nax=nax+lswin(isa,ipsw)
148    4 continue
149c      else
150cc
151cc     periodic boundary conditions
152cc
153c      do 11 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
154c      ispj=lswjpt(isa,ipsw)-1
155c      ism=isml(isfr+isa)
156cc
157cc     collect center of mass distance vectors
158cc
159c      do 12 ismn=1,lswin(isa,ipsw)
160c      lswptr=lswj(ispj+ismn)
161c      rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1)
162c      rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2)
163c      rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3)
164c      rwx(ismn,1)=zero
165c      rwx(ismn,2)=zero
166c      rwx(ismn,3)=zero
167c      if(abs(rwc(nax+ismn,1)).gt.boxh(1)) then
168c      rwx(ismn,1)=sign(box(1),xs(isfr+isa,1))
169c      endif
170c      if(abs(rwc(nax+ismn,2)).gt.boxh(2)) then
171c      rwx(ismn,2)=sign(box(2),xs(isfr+isa,2))
172c      endif
173c      if(npbtyp.eq.1) then
174c      if(abs(rwc(nax+ismn,3)).gt.boxh(3)) then
175c      rwx(ismn,3)=sign(box(3),xs(isfr+isa,3))
176c      endif
177c      endif
178c      if(ism.gt.0) then
179c      rwc(nax+ismn,1)=xsm(ism,1)-xwm(lswptr,1)-rwx(ismn,1)
180c      rwc(nax+ismn,2)=xsm(ism,2)-xwm(lswptr,2)-rwx(ismn,2)
181c      rwc(nax+ismn,3)=xsm(ism,3)-xwm(lswptr,3)-rwx(ismn,3)
182c      else
183c      rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1)-rwx(ismn,1)
184c      rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2)-rwx(ismn,2)
185c      rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3)-rwx(ismn,3)
186c      endif
187c   12 continue
188cc
189cc     collect solute coordinates and atomic polarization fields
190cc
191c      do 13 ismn=1,lswin(isa,ipsw)
192c      lswptr=lswj(ispj+ismn)
193c      xi(nax+ismn,1)=xs(isfr+isa,1)
194c      xi(nax+ismn,2)=xs(isfr+isa,2)
195c      xi(nax+ismn,3)=xs(isfr+isa,3)
196c      pl(nax+ismn,1)=ps(isfr+isa,1,2)
197c      pl(nax+ismn,2)=ps(isfr+isa,2,2)
198c      pl(nax+ismn,3)=ps(isfr+isa,3,2)
199c      isal(nax+ismn)=isfr+isa
200c   13 continue
201cc
202cc     collect solvent coordinates and atomic polarization fields
203cc
204c      do 14 jwa=1,mwa
205c      do 15 ismn=1,lswin(isa,ipsw)
206c      lswptr=lswj(ispj+ismn)
207c      xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa)+rwx(ismn,1)
208c      xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa)+rwx(ismn,2)
209c      xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa)+rwx(ismn,3)
210c      pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2)
211c      pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2)
212c      pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2)
213c   15 continue
214c   14 continue
215c      nax=nax+lswin(isa,ipsw)
216c   11 continue
217c      endif
218c
219c     zero temparary arays fi and fj
220c
221      do 16 iax=1,nax
222      fi(iax,1)=zero
223      fi(iax,2)=zero
224      fi(iax,3)=zero
225   16 continue
226      do 17 jwa=1,mwa
227      do 18 iax=1,nax
228      fj(iax,1,jwa)=zero
229      fj(iax,2,jwa)=zero
230      fj(iax,3,jwa)=zero
231   18 continue
232   17 continue
233c
234c     loop over the number of atoms in a solvent molecule
235c
236c     calculated here is 4*pi*epsilon*field and not just the field
237c     since the polarization is given in alpha/(4*pi*epsilon) in
238c     stead of just alpha, the induced dipole is obtained by the
239c     product of pwa and pw
240c
241      do 19 jwa=1,mwa
242      qaj=qfaci*chg(iwq(jwa),1,iset)
243      paj=chg(iwq(jwa),2,iset)
244      do 20 iax=1,nax
245      isa=isal(iax)
246      qai=qfaci*chg(isq1(isa),1,iset)
247      pai=chg(isq1(isa),2,iset)
248      rx=xj(iax,1,jwa)-xi(iax,1)
249      ry=xj(iax,2,jwa)-xi(iax,2)
250      rz=xj(iax,3,jwa)-xi(iax,3)
251      pix=pai*pl(iax,1)
252      piy=pai*pl(iax,2)
253      piz=pai*pl(iax,3)
254      pjx=paj*pj(iax,1,jwa)
255      pjy=paj*pj(iax,2,jwa)
256      pjz=paj*pj(iax,3,jwa)
257      ri2=one/(rx**2+ry**2+rz**2)
258      ri3=sqrt(ri2)*ri2
259      rmi=three*(rx*pix+ry*piy+rz*piz)*ri2
260      rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2
261      fi(iax,1)=fi(iax,1)+((rmj-qaj)*rx-pjx)*ri3
262      fi(iax,2)=fi(iax,2)+((rmj-qaj)*ry-pjy)*ri3
263      fi(iax,3)=fi(iax,3)+((rmj-qaj)*rz-pjz)*ri3
264      fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3
265      fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3
266      fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3
267   20 continue
268   19 continue
269c
270c     update the electric field arrays
271c
272      iax=0
273      do 21 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
274      ispj=lswjpt(isa,ipsw)-1
275      do 22 ismn=1,lswin(isa,ipsw)
276      ps(isfr+isa,1,1)=ps(isfr+isa,1,1)+fi(iax+ismn,1)
277      ps(isfr+isa,2,1)=ps(isfr+isa,2,1)+fi(iax+ismn,2)
278      ps(isfr+isa,3,1)=ps(isfr+isa,3,1)+fi(iax+ismn,3)
279   22 continue
280      do 23 jwa=1,mwa
281      do 24 ismn=1,lswin(isa,ipsw)
282      lswptr=lswj(ispj+ismn)
283      pw(lswptr,1,jwa,1)=pw(lswptr,1,jwa,1)+fj(iax+ismn,1,jwa)
284      pw(lswptr,2,jwa,1)=pw(lswptr,2,jwa,1)+fj(iax+ismn,2,jwa)
285      pw(lswptr,3,jwa,1)=pw(lswptr,3,jwa,1)+fj(iax+ismn,3,jwa)
286   24 continue
287   23 continue
288      iax=iax+lswin(isa,ipsw)
289   21 continue
290c
291c     thermodynamic integration and perturbation
292c
293      do 30 ipset=2,3
294      if((ipset.eq.2.and.ipert2).or.
295     + (ipset.eq.3.and.ipert3)) then
296c
297      nax2=0
298      do 31 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
299      ispj=lswjpt(isa,ipsw)-1
300      do 32 ismn=1,lswin(isa,ipsw)
301      lswptr=lswj(ispj+ismn)
302      pl(nax2+ismn,1)=psp(isfr+isa,1,2,ipset-1)
303      pl(nax2+ismn,2)=psp(isfr+isa,2,2,ipset-1)
304      pl(nax2+ismn,3)=psp(isfr+isa,3,2,ipset-1)
305   32 continue
306      do 33 jwa=1,mwa
307      do 34 ismn=1,lswin(isa,ipsw)
308      lswptr=lswj(ispj+ismn)
309      pj(nax2+ismn,1,jwa)=pwp(lswptr,1,jwa,2,ipset-1)
310      pj(nax2+ismn,2,jwa)=pwp(lswptr,2,jwa,2,ipset-1)
311      pj(nax2+ismn,3,jwa)=pwp(lswptr,3,jwa,2,ipset-1)
312   34 continue
313   33 continue
314      nax2=nax2+lswin(isa,ipsw)
315   31 continue
316c
317      if(nax.ne.nax2) call md_abort('Error in dipsw',me)
318c
319      do 41 iax=1,nax
320      fi(iax,1)=zero
321      fi(iax,2)=zero
322      fi(iax,3)=zero
323   41 continue
324      do 42 jwa=1,mwa
325      do 43 iax=1,nax
326      fj(iax,1,jwa)=zero
327      fj(iax,2,jwa)=zero
328      fj(iax,3,jwa)=zero
329   43 continue
330   42 continue
331c
332      do 35 jwa=1,mwa
333      qaj=qfaci*chg(iwq(jwa),1,ipset)
334      paj=chg(iwq(jwa),2,ipset)
335      do 36 iax=1,nax
336      isa=isal(iax)
337      qai=qfaci*chg(isq1(isa),1,ipset)
338      pai=chg(isq1(isa),2,ipset)
339      rx=xj(iax,1,jwa)-xi(iax,1)
340      ry=xj(iax,2,jwa)-xi(iax,2)
341      rz=xj(iax,3,jwa)-xi(iax,3)
342      pix=pai*pl(iax,1)
343      piy=pai*pl(iax,2)
344      piz=pai*pl(iax,3)
345      pjx=paj*pj(iax,1,jwa)
346      pjy=paj*pj(iax,2,jwa)
347      pjz=paj*pj(iax,3,jwa)
348      ri2=one/(rx**2+ry**2+rz**2)
349      ri3=sqrt(ri2)*ri2
350      rmi=three*(rx*pix+ry*piy+rz*piz)*ri2
351      rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2
352      fi(iax,1)=fi(iax,1)+((rmj-qaj)*rx-pjx)*ri3
353      fi(iax,2)=fi(iax,2)+((rmj-qaj)*ry-pjy)*ri3
354      fi(iax,3)=fi(iax,3)+((rmj-qaj)*rz-pjz)*ri3
355      fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3
356      fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3
357      fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3
358   36 continue
359   35 continue
360c
361c     update the electric field arrays
362c
363      iax=0
364      do 37 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
365      ispj=lswjpt(isa,ipsw)-1
366      do 38 ismn=1,lswin(isa,ipsw)
367      psp(isfr+isa,1,ipset-1,1)=psp(isfr+isa,1,ipset-1,1)+fi(iax+ismn,1)
368      psp(isfr+isa,2,ipset-1,1)=psp(isfr+isa,2,ipset-1,1)+fi(iax+ismn,2)
369      psp(isfr+isa,3,ipset-1,1)=psp(isfr+isa,3,ipset-1,1)+fi(iax+ismn,3)
370   38 continue
371      do 39 jwa=1,mwa
372      do 40 ismn=1,lswin(isa,ipsw)
373      lswptr=lswj(ispj+ismn)
374      pwp(lswptr,1,jwa,ipset-1,1)=pwp(lswptr,1,jwa,ipset-1,1)+
375     + fj(iax+ismn,1,jwa)
376      pwp(lswptr,2,jwa,ipset-1,1)=pwp(lswptr,2,jwa,ipset-1,1)+
377     + fj(iax+ismn,2,jwa)
378      pwp(lswptr,3,jwa,ipset-1,1)=pwp(lswptr,3,jwa,ipset-1,1)+
379     + fj(iax+ismn,3,jwa)
380   40 continue
381   39 continue
382      iax=iax+lswin(isa,ipsw)
383   37 continue
384      endif
385   30 continue
386c
387    3 continue
388    1 continue
389c
390      return
391      end
392