1      subroutine cf_dww(xw,xwm,pw,pwp,iwfrom,nlocw,lpbc,chg,iwq,
2     + lwwndx,lwwjpt,lwwin,lwwj,rwc,xi,xj,rwx,pl,pj,fi,fj)
3c
4c $Id$
5c
6      implicit none
7c
8#include "cf_common.fh"
9c
10      real*8 xw(mwm,3,mwa),xwm(mwm,3)
11      real*8 pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2)
12      integer iwfrom,nlocw
13      logical lpbc
14      real*8 chg(mqt,mqp,mset)
15      integer iwq(mwa)
16      integer lwwndx(0:mwm,2),lwwin(nlocw,2),lwwjpt(nlocw,2),lwwj(*)
17c
18      real*8 rwc(mscr,3),rwx(mscr,3)
19      real*8 xi(mscr,3,mwa),xj(mscr,3,mwa)
20      real*8 pl(mscr,3,mwa),pj(mscr,3,mwa)
21      real*8 fi(mscr,3,mwa),fj(mscr,3,mwa)
22c
23      integer ix,ipset,nax2,nwwlen(2)
24      integer iwfr,ipww,number,iwm,iwpm,nax,iwpj,iwa,jwa,iax,iwmn,lwwptr
25      real*8 pai,paj,qai,qaj
26      real*8 rx,ry,rz,ri2,ri3,rmi,rmj,pix,piy,piz,pjx,pjy,pjz
27c
28      real*8 qfaci
29c
30      qfaci=one/qfac
31c
32c     calculation of solvent-solvent intermolecular energies and forces
33c
34      iwfr=iwfrom-1
35c
36c     loop over both short and long range parts
37c
38      do 1 ipww=1,npww
39c
40c     Evaluate the outer index array
41c
42      nwwlen(ipww)=0
43      lwwndx(0,ipww)=0
44      number=0
45      do 2 iwm=1,nlocw
46      if(number+lwwin(iwm,ipww).gt.mscr) then
47      nwwlen(ipww)=nwwlen(ipww)+1
48      lwwndx(nwwlen(ipww),ipww)=iwm-1
49      number=0
50      endif
51      number=number+lwwin(iwm,ipww)
52    2 continue
53      if(number.gt.0) then
54      nwwlen(ipww)=nwwlen(ipww)+1
55      lwwndx(nwwlen(ipww),ipww)=nlocw
56      endif
57c
58c     loop over number of cycles to complete pairlist
59c
60      do 3 iwpm=1,nwwlen(ipww)
61      nax=0
62c
63c     collect coordinates into workarrays
64c
65      do 4 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww)
66      iwpj=lwwjpt(iwm,ipww)-1
67      do 5 iwmn=1,lwwin(iwm,ipww)
68      lwwptr=lwwj(iwpj+iwmn)
69      rwc(nax+iwmn,1)=xwm(lwwptr,1)-xwm(iwfr+iwm,1)
70      rwc(nax+iwmn,2)=xwm(lwwptr,2)-xwm(iwfr+iwm,2)
71      rwc(nax+iwmn,3)=xwm(lwwptr,3)-xwm(iwfr+iwm,3)
72    5 continue
73c
74      if(.not.lpbc) then
75      do 6 iwa=1,mwa
76      do 7 iwmn=1,lwwin(iwm,ipww)
77      lwwptr=lwwj(iwpj+iwmn)
78      xi(nax+iwmn,1,iwa)=xw(iwfr+iwm,1,iwa)
79      xi(nax+iwmn,2,iwa)=xw(iwfr+iwm,2,iwa)
80      xi(nax+iwmn,3,iwa)=xw(iwfr+iwm,3,iwa)
81      xj(nax+iwmn,1,iwa)=xw(lwwptr,1,iwa)
82      xj(nax+iwmn,2,iwa)=xw(lwwptr,2,iwa)
83      xj(nax+iwmn,3,iwa)=xw(lwwptr,3,iwa)
84      pl(nax+iwmn,1,iwa)=pw(iwfr+iwm,1,iwa,2)
85      pl(nax+iwmn,2,iwa)=pw(iwfr+iwm,2,iwa,2)
86      pl(nax+iwmn,3,iwa)=pw(iwfr+iwm,3,iwa,2)
87      pj(nax+iwmn,1,iwa)=pw(lwwptr,1,iwa,2)
88      pj(nax+iwmn,2,iwa)=pw(lwwptr,2,iwa,2)
89      pj(nax+iwmn,3,iwa)=pw(lwwptr,3,iwa,2)
90    7 continue
91    6 continue
92      else
93      call cf_pbc(0,rwc,mscr,rwx,mscr,nax,1,lwwin(iwm,ipww))
94      do 9 iwmn=1,lwwin(iwm,ipww)
95      lwwptr=lwwj(iwpj+iwmn)
96      rwc(nax+iwmn,1)=rwc(nax+iwmn,1)-rwx(iwmn,1)
97      rwc(nax+iwmn,2)=rwc(nax+iwmn,2)-rwx(iwmn,2)
98      rwc(nax+iwmn,3)=rwc(nax+iwmn,3)-rwx(iwmn,3)
99    9 continue
100c
101      do 10 iwa=1,mwa
102      do 11 iwmn=1,lwwin(iwm,ipww)
103      lwwptr=lwwj(iwpj+iwmn)
104      xi(nax+iwmn,1,iwa)=xw(iwfr+iwm,1,iwa)
105      xi(nax+iwmn,2,iwa)=xw(iwfr+iwm,2,iwa)
106      xi(nax+iwmn,3,iwa)=xw(iwfr+iwm,3,iwa)
107      xj(nax+iwmn,1,iwa)=xw(lwwptr,1,iwa)-rwx(iwmn,1)
108      xj(nax+iwmn,2,iwa)=xw(lwwptr,2,iwa)-rwx(iwmn,2)
109      xj(nax+iwmn,3,iwa)=xw(lwwptr,3,iwa)-rwx(iwmn,3)
110      pl(nax+iwmn,1,iwa)=pw(iwfr+iwm,1,iwa,2)
111      pl(nax+iwmn,2,iwa)=pw(iwfr+iwm,2,iwa,2)
112      pl(nax+iwmn,3,iwa)=pw(iwfr+iwm,3,iwa,2)
113      pj(nax+iwmn,1,iwa)=pw(lwwptr,1,iwa,2)
114      pj(nax+iwmn,2,iwa)=pw(lwwptr,2,iwa,2)
115      pj(nax+iwmn,3,iwa)=pw(lwwptr,3,iwa,2)
116   11 continue
117   10 continue
118      endif
119      nax=nax+lwwin(iwm,ipww)
120    4 continue
121c
122c     zero temporary arrays fi and fj
123c
124      do 12 iwa=1,mwa
125      do 13 ix=1,3
126      do 14 iax=1,nax
127      fi(iax,ix,iwa)=zero
128      fj(iax,ix,iwa)=zero
129   14 continue
130   13 continue
131   12 continue
132c
133c     loops over number of atoms in a solvent molecule
134c
135c     calculated here is 4*pi*epsilon*field and not just the field
136c     since the polarization is given in alpha/(4*pi*epsilon) in
137c     stead of just alpha, the induced dipole is obtained by the
138c     product of pwa and pw
139c
140      do 20 iwa=1,mwa
141      qai=qfaci*chg(iwq(iwa),1,iset)
142      pai=chg(iwq(iwa),2,iset)
143      do 21 jwa=1,mwa
144      qaj=qfaci*chg(iwq(jwa),1,iset)
145      paj=chg(iwq(jwa),2,iset)
146      do 22 iax=1,nax
147      rx=xj(iax,1,jwa)-xi(iax,1,iwa)
148      ry=xj(iax,2,jwa)-xi(iax,2,iwa)
149      rz=xj(iax,3,jwa)-xi(iax,3,iwa)
150      pix=pai*pl(iax,1,iwa)
151      piy=pai*pl(iax,2,iwa)
152      piz=pai*pl(iax,3,iwa)
153      pjx=paj*pj(iax,1,jwa)
154      pjy=paj*pj(iax,2,jwa)
155      pjz=paj*pj(iax,3,jwa)
156      ri2=one/(rx**2+ry**2+rz**2)
157      ri3=sqrt(ri2)*ri2
158      rmi=three*(rx*pix+ry*piy+rz*piz)*ri2
159      rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2
160      fi(iax,1,iwa)=fi(iax,1,iwa)+((rmj-qaj)*rx-pjx)*ri3
161      fi(iax,2,iwa)=fi(iax,2,iwa)+((rmj-qaj)*ry-pjy)*ri3
162      fi(iax,3,iwa)=fi(iax,3,iwa)+((rmj-qaj)*rz-pjz)*ri3
163      fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3
164      fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3
165      fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3
166   22 continue
167   21 continue
168   20 continue
169c
170c     Update the electric field arrays
171c
172      iax=0
173      do 23 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww)
174      iwpj=lwwjpt(iwm,ipww)-1
175      do 24 iwa=1,mwa
176      do 25 iwmn=1,lwwin(iwm,ipww)
177      lwwptr=lwwj(iwpj+iwmn)
178      pw(iwfr+iwm,1,iwa,1)=pw(iwfr+iwm,1,iwa,1)+fi(iax+iwmn,1,iwa)
179      pw(iwfr+iwm,2,iwa,1)=pw(iwfr+iwm,2,iwa,1)+fi(iax+iwmn,2,iwa)
180      pw(iwfr+iwm,3,iwa,1)=pw(iwfr+iwm,3,iwa,1)+fi(iax+iwmn,3,iwa)
181      pw(lwwptr,1,iwa,1)=pw(lwwptr,1,iwa,1)+fj(iax+iwmn,1,iwa)
182      pw(lwwptr,2,iwa,1)=pw(lwwptr,2,iwa,1)+fj(iax+iwmn,2,iwa)
183      pw(lwwptr,3,iwa,1)=pw(lwwptr,3,iwa,1)+fj(iax+iwmn,3,iwa)
184   25 continue
185   24 continue
186      iax=iax+lwwin(iwm,ipww)
187   23 continue
188c
189c     thermodynamic perturbation and integration
190c
191      do 30 ipset=2,3
192      if((ipset.eq.2.and.ipert2).or.
193     + (ipset.eq.3.and.ipert3)) then
194c
195      nax2=0
196      do 31 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww)
197      iwpj=lwwjpt(iwm,ipww)-1
198      do 32 iwa=1,mwa
199      do 33 iwmn=1,lwwin(iwm,ipww)
200      lwwptr=lwwj(iwpj+iwmn)
201      pl(nax2+iwmn,1,iwa)=pwp(iwfr+iwm,1,iwa,ipset-1,2)
202      pl(nax2+iwmn,2,iwa)=pwp(iwfr+iwm,2,iwa,ipset-1,2)
203      pl(nax2+iwmn,3,iwa)=pwp(iwfr+iwm,3,iwa,ipset-1,2)
204      pj(nax2+iwmn,1,iwa)=pwp(lwwptr,1,iwa,ipset-1,2)
205      pj(nax2+iwmn,2,iwa)=pwp(lwwptr,2,iwa,ipset-1,2)
206      pj(nax2+iwmn,3,iwa)=pwp(lwwptr,3,iwa,ipset-1,2)
207   33 continue
208   32 continue
209      nax2=nax2+lwwin(iwm,ipww)
210   31 continue
211c
212      if(nax.ne.nax2) call md_abort('Error in dipww',me)
213c
214c
215      do 40 iwa=1,mwa
216      do 41 ix=1,3
217      do 42 iax=1,nax
218      fi(iax,ix,iwa)=zero
219      fj(iax,ix,iwa)=zero
220   42 continue
221   41 continue
222   40 continue
223c
224      do 34 iwa=1,mwa
225      qai=qfaci*chg(iwq(iwa),1,ipset)
226      pai=chg(iwq(iwa),2,ipset)
227      do 35 jwa=1,mwa
228      qaj=qfaci*chg(iwq(jwa),1,ipset)
229      paj=chg(iwq(jwa),2,ipset)
230      do 36 iax=1,nax
231      rx=xj(iax,1,jwa)-xi(iax,1,iwa)
232      ry=xj(iax,2,jwa)-xi(iax,2,iwa)
233      rz=xj(iax,3,jwa)-xi(iax,3,iwa)
234      pix=pai*pl(iax,1,iwa)
235      piy=pai*pl(iax,2,iwa)
236      piz=pai*pl(iax,3,iwa)
237      pjx=paj*pj(iax,1,jwa)
238      pjy=paj*pj(iax,2,jwa)
239      pjz=paj*pj(iax,3,jwa)
240      ri2=one/(rx**2+ry**2+rz**2)
241      ri3=sqrt(ri2)*ri2
242      rmi=three*(rx*pix+ry*piy+rz*piz)*ri2
243      rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2
244      fi(iax,1,iwa)=fi(iax,1,iwa)+((rmj-qaj)*rx-pjx)*ri3
245      fi(iax,2,iwa)=fi(iax,2,iwa)+((rmj-qaj)*ry-pjy)*ri3
246      fi(iax,3,iwa)=fi(iax,3,iwa)+((rmj-qaj)*rz-pjz)*ri3
247      fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3
248      fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3
249      fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3
250   36 continue
251   35 continue
252   34 continue
253c
254c     Update the electric field arrays
255c
256      iax=0
257      do 37 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww)
258      iwpj=lwwjpt(iwm,ipww)-1
259      do 38 iwa=1,mwa
260      do 39 iwmn=1,lwwin(iwm,ipww)
261      lwwptr=lwwj(iwpj+iwmn)
262      pwp(iwfr+iwm,1,iwa,ipset-1,1)=pwp(iwfr+iwm,1,iwa,ipset-1,1)+
263     + fi(iax+iwmn,1,iwa)
264      pwp(iwfr+iwm,2,iwa,ipset-1,1)=pwp(iwfr+iwm,2,iwa,ipset-1,1)+
265     + fi(iax+iwmn,2,iwa)
266      pwp(iwfr+iwm,3,iwa,ipset-1,1)=pwp(iwfr+iwm,3,iwa,ipset-1,1)+
267     + fi(iax+iwmn,3,iwa)
268      pwp(lwwptr,1,iwa,ipset-1,1)=pwp(lwwptr,1,iwa,ipset-1,1)+
269     + fj(iax+iwmn,1,iwa)
270      pwp(lwwptr,2,iwa,ipset-1,1)=pwp(lwwptr,2,iwa,ipset-1,1)+
271     + fj(iax+iwmn,2,iwa)
272      pwp(lwwptr,3,iwa,ipset-1,1)=pwp(lwwptr,3,iwa,ipset-1,1)+
273     + fj(iax+iwmn,3,iwa)
274   39 continue
275   38 continue
276      iax=iax+lwwin(iwm,ipww)
277   37 continue
278      endif
279   30 continue
280c
281c
282    3 continue
283    1 continue
284c
285      return
286      end
287      subroutine cf_dsw(xs,xsm,ps,psp,isdt,ismf,isml,isq1,isfrom,nums,
288     + xw,xwm,pw,pwp,iwq,lpbc,chg,lswndx,lswjpt,lswin,lswj,
289     + rwc,xi,xj,rwx,pl,pj,fi,fj,isal)
290c
291c $Id$
292c
293      implicit none
294c
295#include "cf_common.fh"
296c
297      real*8 xs(msa,3),xsm(msm,3),ps(msa,3,2),psp(msa,3,2,2)
298      integer isdt(msa),ismf(msa),isml(msa),isq1(msa)
299      integer isfrom,nums
300      real*8 xw(mwm,3,mwa),xwm(mwm,3),pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2)
301      integer iwq(mwa)
302      logical lpbc
303      real*8 chg(mqt,mqp,mset)
304      integer lswndx(0:msa,2),lswjpt(nums,2),lswin(nums,2),lswj(*)
305      real*8 rwc(mscr,3),xi(mscr,3),xj(mscr,3,mwa),rwx(mscr,3)
306      real*8 pl(mscr,3),pj(mscr,3,mwa),fi(mscr,3),fj(mscr,3,mwa)
307      integer isal(mscr)
308c
309      integer ispj,ism,lswptr,ipset,nswlen(2)
310      integer isfr,ipsw,number,isa,jwa,ismn,ispm,iax,nax,nax2
311      real*8 qai,qaj,pai,paj,rx,ry,rz,pix,piy,piz,pjx,pjy,pjz
312      real*8 ri2,ri3,rmi,rmj
313#include "bitops.fh"
314c
315      real*8 qfaci
316c
317      qfaci=one/qfac
318c
319c     this subroutine evaluates the solute-solvent forces for nums
320c     solute atoms starting from isfrom. the interacting solvent
321c     molecules are determined from the pairlist.
322c
323      isfr=isfrom-1
324c
325c     loop over short and long range pairs
326c
327      do 1 ipsw=1,npsw
328c
329c     evaluate outer index array
330c
331      nswlen(ipsw)=0
332      lswndx(0,ipsw)=0
333      number=0
334      do 2 isa=1,nums
335      if(number+lswin(isa,ipsw).gt.mscr .or.
336     + (ismf(isfr+isa).ne.ismf(isfr+isa-1).and.number.gt.0)) then
337      nswlen(ipsw)=nswlen(ipsw)+1
338      lswndx(nswlen(ipsw),ipsw)=isa-1
339      number=0
340      endif
341      number=number+lswin(isa,ipsw)
342    2 continue
343      if(number.gt.0) then
344      nswlen(ipsw)=nswlen(ipsw)+1
345      lswndx(nswlen(ipsw),ipsw)=nums
346      endif
347c
348c     loop over number of cycles to complete pairlist
349c
350      do 3 ispm=1,nswlen(ipsw)
351      nax=0
352c
353c     vacuo conditions
354c
355c      if(npbtyp.eq.0) then
356c      if(.not.lpbc) then
357      do 4 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
358      ispj=lswjpt(isa,ipsw)-1
359      ism=isml(isfr+isa)
360c
361c     collect center of mass distance vectors
362c
363      if(lpbc.or.ism.eq.0) then
364      do 6 ismn=1,lswin(isa,ipsw)
365      lswptr=lswj(ispj+ismn)
366      rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1)
367      rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2)
368      rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3)
369    6 continue
370      if(lpbc) call cf_pbc(0,rwc,mscr,rwx,mscr,nax,1,lswin(isa,ipsw))
371      else
372      do 5 ismn=1,lswin(isa,ipsw)
373      lswptr=lswj(ispj+ismn)
374      rwc(nax+ismn,1)=xsm(ism,1)-xwm(lswptr,1)
375      rwc(nax+ismn,2)=xsm(ism,2)-xwm(lswptr,2)
376      rwc(nax+ismn,3)=xsm(ism,3)-xwm(lswptr,3)
377    5 continue
378      endif
379c
380c     collect solute coordinates and atomic polarization fields
381c
382c      if(iand(isdt(isfr+isa),mdynam).eq.ldynam) then
383      do 7 ismn=1,lswin(isa,ipsw)
384      lswptr=lswj(ispj+ismn)
385      xi(nax+ismn,1)=xs(isfr+isa,1)
386      xi(nax+ismn,2)=xs(isfr+isa,2)
387      xi(nax+ismn,3)=xs(isfr+isa,3)
388      pl(nax+ismn,1)=ps(isfr+isa,1,2)
389      pl(nax+ismn,2)=ps(isfr+isa,2,2)
390      pl(nax+ismn,3)=ps(isfr+isa,3,2)
391      isal(nax+ismn)=isfr+isa
392    7 continue
393c      else
394c      do 8 ismn=1,lswin(isa,ipsw)
395c      lswptr=lswj(ispj+ismn)
396c      xi(nax+ismn,1)=xs(isfr+isa,1)
397c      xi(nax+ismn,2)=xs(isfr+isa,2)
398c      xi(nax+ismn,3)=xs(isfr+isa,3)
399c      pl(nax+ismn,1)=ps(isfr+isa,1,2)
400c      pl(nax+ismn,2)=ps(isfr+isa,2,2)
401c      pl(nax+ismn,3)=ps(isfr+isa,3,2)
402c      isal(nax+ismn)=isfr+isa
403c    8 continue
404c      endif
405c
406c     collect solvent coordinates and atomic polarization fields
407c
408      do 8 jwa=1,mwa
409      if(lpbc) then
410      do 9 ismn=1,lswin(isa,ipsw)
411      lswptr=lswj(ispj+ismn)
412      xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa)+rwx(ismn,1)
413      xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa)+rwx(ismn,2)
414      xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa)+rwx(ismn,3)
415      pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2)
416      pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2)
417      pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2)
418    9 continue
419      else
420      do 10 ismn=1,lswin(isa,ipsw)
421      lswptr=lswj(ispj+ismn)
422      xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa)
423      xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa)
424      xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa)
425      pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2)
426      pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2)
427      pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2)
428   10 continue
429      endif
430    8 continue
431      nax=nax+lswin(isa,ipsw)
432    4 continue
433c      else
434cc
435cc     periodic boundary conditions
436cc
437c      do 11 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
438c      ispj=lswjpt(isa,ipsw)-1
439c      ism=isml(isfr+isa)
440cc
441cc     collect center of mass distance vectors
442cc
443c      do 12 ismn=1,lswin(isa,ipsw)
444c      lswptr=lswj(ispj+ismn)
445c      rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1)
446c      rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2)
447c      rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3)
448c      rwx(ismn,1)=zero
449c      rwx(ismn,2)=zero
450c      rwx(ismn,3)=zero
451c      if(abs(rwc(nax+ismn,1)).gt.boxh(1)) then
452c      rwx(ismn,1)=sign(box(1),xs(isfr+isa,1))
453c      endif
454c      if(abs(rwc(nax+ismn,2)).gt.boxh(2)) then
455c      rwx(ismn,2)=sign(box(2),xs(isfr+isa,2))
456c      endif
457c      if(npbtyp.eq.1) then
458c      if(abs(rwc(nax+ismn,3)).gt.boxh(3)) then
459c      rwx(ismn,3)=sign(box(3),xs(isfr+isa,3))
460c      endif
461c      endif
462c      if(ism.gt.0) then
463c      rwc(nax+ismn,1)=xsm(ism,1)-xwm(lswptr,1)-rwx(ismn,1)
464c      rwc(nax+ismn,2)=xsm(ism,2)-xwm(lswptr,2)-rwx(ismn,2)
465c      rwc(nax+ismn,3)=xsm(ism,3)-xwm(lswptr,3)-rwx(ismn,3)
466c      else
467c      rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1)-rwx(ismn,1)
468c      rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2)-rwx(ismn,2)
469c      rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3)-rwx(ismn,3)
470c      endif
471c   12 continue
472cc
473cc     collect solute coordinates and atomic polarization fields
474cc
475c      do 13 ismn=1,lswin(isa,ipsw)
476c      lswptr=lswj(ispj+ismn)
477c      xi(nax+ismn,1)=xs(isfr+isa,1)
478c      xi(nax+ismn,2)=xs(isfr+isa,2)
479c      xi(nax+ismn,3)=xs(isfr+isa,3)
480c      pl(nax+ismn,1)=ps(isfr+isa,1,2)
481c      pl(nax+ismn,2)=ps(isfr+isa,2,2)
482c      pl(nax+ismn,3)=ps(isfr+isa,3,2)
483c      isal(nax+ismn)=isfr+isa
484c   13 continue
485cc
486cc     collect solvent coordinates and atomic polarization fields
487cc
488c      do 14 jwa=1,mwa
489c      do 15 ismn=1,lswin(isa,ipsw)
490c      lswptr=lswj(ispj+ismn)
491c      xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa)+rwx(ismn,1)
492c      xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa)+rwx(ismn,2)
493c      xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa)+rwx(ismn,3)
494c      pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2)
495c      pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2)
496c      pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2)
497c   15 continue
498c   14 continue
499c      nax=nax+lswin(isa,ipsw)
500c   11 continue
501c      endif
502c
503c     zero temparary arays fi and fj
504c
505      do 16 iax=1,nax
506      fi(iax,1)=zero
507      fi(iax,2)=zero
508      fi(iax,3)=zero
509   16 continue
510      do 17 jwa=1,mwa
511      do 18 iax=1,nax
512      fj(iax,1,jwa)=zero
513      fj(iax,2,jwa)=zero
514      fj(iax,3,jwa)=zero
515   18 continue
516   17 continue
517c
518c     loop over the number of atoms in a solvent molecule
519c
520c     calculated here is 4*pi*epsilon*field and not just the field
521c     since the polarization is given in alpha/(4*pi*epsilon) in
522c     stead of just alpha, the induced dipole is obtained by the
523c     product of pwa and pw
524c
525      do 19 jwa=1,mwa
526      qaj=qfaci*chg(iwq(jwa),1,iset)
527      paj=chg(iwq(jwa),2,iset)
528      do 20 iax=1,nax
529      isa=isal(iax)
530      qai=qfaci*chg(isq1(isa),1,iset)
531      pai=chg(isq1(isa),2,iset)
532      rx=xj(iax,1,jwa)-xi(iax,1)
533      ry=xj(iax,2,jwa)-xi(iax,2)
534      rz=xj(iax,3,jwa)-xi(iax,3)
535      pix=pai*pl(iax,1)
536      piy=pai*pl(iax,2)
537      piz=pai*pl(iax,3)
538      pjx=paj*pj(iax,1,jwa)
539      pjy=paj*pj(iax,2,jwa)
540      pjz=paj*pj(iax,3,jwa)
541      ri2=one/(rx**2+ry**2+rz**2)
542      ri3=sqrt(ri2)*ri2
543      rmi=three*(rx*pix+ry*piy+rz*piz)*ri2
544      rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2
545      fi(iax,1)=fi(iax,1)+((rmj-qaj)*rx-pjx)*ri3
546      fi(iax,2)=fi(iax,2)+((rmj-qaj)*ry-pjy)*ri3
547      fi(iax,3)=fi(iax,3)+((rmj-qaj)*rz-pjz)*ri3
548      fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3
549      fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3
550      fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3
551   20 continue
552   19 continue
553c
554c     update the electric field arrays
555c
556      iax=0
557      do 21 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
558      ispj=lswjpt(isa,ipsw)-1
559      do 22 ismn=1,lswin(isa,ipsw)
560      ps(isfr+isa,1,1)=ps(isfr+isa,1,1)+fi(iax+ismn,1)
561      ps(isfr+isa,2,1)=ps(isfr+isa,2,1)+fi(iax+ismn,2)
562      ps(isfr+isa,3,1)=ps(isfr+isa,3,1)+fi(iax+ismn,3)
563   22 continue
564      do 23 jwa=1,mwa
565      do 24 ismn=1,lswin(isa,ipsw)
566      lswptr=lswj(ispj+ismn)
567      pw(lswptr,1,jwa,1)=pw(lswptr,1,jwa,1)+fj(iax+ismn,1,jwa)
568      pw(lswptr,2,jwa,1)=pw(lswptr,2,jwa,1)+fj(iax+ismn,2,jwa)
569      pw(lswptr,3,jwa,1)=pw(lswptr,3,jwa,1)+fj(iax+ismn,3,jwa)
570   24 continue
571   23 continue
572      iax=iax+lswin(isa,ipsw)
573   21 continue
574c
575c     thermodynamic integration and perturbation
576c
577      do 30 ipset=2,3
578      if((ipset.eq.2.and.ipert2).or.
579     + (ipset.eq.3.and.ipert3)) then
580c
581      nax2=0
582      do 31 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
583      ispj=lswjpt(isa,ipsw)-1
584      do 32 ismn=1,lswin(isa,ipsw)
585      lswptr=lswj(ispj+ismn)
586      pl(nax2+ismn,1)=psp(isfr+isa,1,2,ipset-1)
587      pl(nax2+ismn,2)=psp(isfr+isa,2,2,ipset-1)
588      pl(nax2+ismn,3)=psp(isfr+isa,3,2,ipset-1)
589   32 continue
590      do 33 jwa=1,mwa
591      do 34 ismn=1,lswin(isa,ipsw)
592      lswptr=lswj(ispj+ismn)
593      pj(nax2+ismn,1,jwa)=pwp(lswptr,1,jwa,2,ipset-1)
594      pj(nax2+ismn,2,jwa)=pwp(lswptr,2,jwa,2,ipset-1)
595      pj(nax2+ismn,3,jwa)=pwp(lswptr,3,jwa,2,ipset-1)
596   34 continue
597   33 continue
598      nax2=nax2+lswin(isa,ipsw)
599   31 continue
600c
601      if(nax.ne.nax2) call md_abort('Error in dipsw',me)
602c
603      do 41 iax=1,nax
604      fi(iax,1)=zero
605      fi(iax,2)=zero
606      fi(iax,3)=zero
607   41 continue
608      do 42 jwa=1,mwa
609      do 43 iax=1,nax
610      fj(iax,1,jwa)=zero
611      fj(iax,2,jwa)=zero
612      fj(iax,3,jwa)=zero
613   43 continue
614   42 continue
615c
616      do 35 jwa=1,mwa
617      qaj=qfaci*chg(iwq(jwa),1,ipset)
618      paj=chg(iwq(jwa),2,ipset)
619      do 36 iax=1,nax
620      isa=isal(iax)
621      qai=qfaci*chg(isq1(isa),1,ipset)
622      pai=chg(isq1(isa),2,ipset)
623      rx=xj(iax,1,jwa)-xi(iax,1)
624      ry=xj(iax,2,jwa)-xi(iax,2)
625      rz=xj(iax,3,jwa)-xi(iax,3)
626      pix=pai*pl(iax,1)
627      piy=pai*pl(iax,2)
628      piz=pai*pl(iax,3)
629      pjx=paj*pj(iax,1,jwa)
630      pjy=paj*pj(iax,2,jwa)
631      pjz=paj*pj(iax,3,jwa)
632      ri2=one/(rx**2+ry**2+rz**2)
633      ri3=sqrt(ri2)*ri2
634      rmi=three*(rx*pix+ry*piy+rz*piz)*ri2
635      rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2
636      fi(iax,1)=fi(iax,1)+((rmj-qaj)*rx-pjx)*ri3
637      fi(iax,2)=fi(iax,2)+((rmj-qaj)*ry-pjy)*ri3
638      fi(iax,3)=fi(iax,3)+((rmj-qaj)*rz-pjz)*ri3
639      fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3
640      fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3
641      fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3
642   36 continue
643   35 continue
644c
645c     update the electric field arrays
646c
647      iax=0
648      do 37 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
649      ispj=lswjpt(isa,ipsw)-1
650      do 38 ismn=1,lswin(isa,ipsw)
651      psp(isfr+isa,1,ipset-1,1)=psp(isfr+isa,1,ipset-1,1)+fi(iax+ismn,1)
652      psp(isfr+isa,2,ipset-1,1)=psp(isfr+isa,2,ipset-1,1)+fi(iax+ismn,2)
653      psp(isfr+isa,3,ipset-1,1)=psp(isfr+isa,3,ipset-1,1)+fi(iax+ismn,3)
654   38 continue
655      do 39 jwa=1,mwa
656      do 40 ismn=1,lswin(isa,ipsw)
657      lswptr=lswj(ispj+ismn)
658      pwp(lswptr,1,jwa,ipset-1,1)=pwp(lswptr,1,jwa,ipset-1,1)+
659     + fj(iax+ismn,1,jwa)
660      pwp(lswptr,2,jwa,ipset-1,1)=pwp(lswptr,2,jwa,ipset-1,1)+
661     + fj(iax+ismn,2,jwa)
662      pwp(lswptr,3,jwa,ipset-1,1)=pwp(lswptr,3,jwa,ipset-1,1)+
663     + fj(iax+ismn,3,jwa)
664   40 continue
665   39 continue
666      iax=iax+lswin(isa,ipsw)
667   37 continue
668      endif
669   30 continue
670c
671    3 continue
672    1 continue
673c
674      return
675      end
676      subroutine cf_dss(xs,xsm,ps,psp,ismf,isml,isq2,isq3,isfrom,nums,
677     + lpbc,chg,lssndx,lssjpt,lssin,lssj,xi,xj,rwc,rwi1,rwi2,rwi6,
678     + rwx,rw,fi,fj,f,isal,jsal,jmal,jfal,pl,pj)
679c
680c $Id$
681c
682      implicit none
683c
684#include "cf_common.fh"
685c
686      real*8 xs(msa,3),xsm(msm,3),ps(msa,3,2),psp(msa,3,2,2)
687      integer ismf(msa),isml(msa),isq2(msa),isq3(msa)
688      integer isfrom,nums
689      logical lpbc
690      real*8 chg(mqt,mqp,mset)
691      integer lssndx(0:msa,2),lssjpt(nums,2),lssin(nums,2)
692c
693      real*8 xi(mscr,3),xj(mscr,3),rwx(mscr,3),rwi1(mscr)
694      real*8 rwi2(mscr),rwi6(mscr),rwc(mscr,3),rw(mscr)
695      real*8 f(mscr),fi(mscr,3),fj(mscr,3),pl(mscr,3),pj(mscr,3)
696      integer isal(mscr),jsal(mscr),jmal(mscr),jfal(mscr)
697      integer lssj(*)
698c
699      integer isa,jsa,jsf,ipset
700      integer isfr,nsslen(2)
701      integer ipss,number,isslen,nax,nax2,jsaptr
702      integer jnum,lssptr,iax,ism,jsm
703c
704      real*8 ri1,ri2,ri3,rx,ry,rz,pix,piy,piz,pjx,pjy,pjz,rmi,rmj
705      real*8 qai,pai,qaj,paj
706c
707#include "bitops.fh"
708c
709      real*8 qfaci
710c
711      qfaci=one/qfac
712c
713c     solute non-bonded pairs
714c     =======================
715c
716      isfr=isfrom-1
717c
718c     loop over both short and long range pairlists
719c
720      do 1 ipss=1,npss
721c
722c     evaluate outer index array
723c
724      nsslen(ipss)=0
725      lssndx(0,ipss)=0
726      number=0
727      do 2 isa=1,nums
728      if(number+lssin(isa,ipss).gt.mscr .or.
729     + (ismf(isfr+isa).ne.ismf(isfr+isa-1).and.number.gt.0)) then
730      nsslen(ipss)=nsslen(ipss)+1
731      lssndx(nsslen(ipss),ipss)=isa-1
732      number=0
733      endif
734      number=number+lssin(isa,ipss)
735    2 continue
736      if(number.gt.0) then
737      nsslen(ipss)=nsslen(ipss)+1
738      lssndx(nsslen(ipss),ipss)=nums
739      endif
740c
741c     loop over number of cycles to complete pairlists
742c
743      do 3 isslen=1,nsslen(ipss)
744      nax=0
745      ism=isml(isfr+lssndx(isslen,ipss))
746c
747c     collect coordinates into workarrays
748c
749c      if(npbtyp.eq.0) then
750      do 4 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss)
751      jsaptr=lssjpt(isa,ipss)-1
752      ism=isml(isfr+isa)
753      if(lpbc) then
754      do 151 jnum=1,lssin(isa,ipss)
755      lssptr=lssj(jsaptr+jnum)
756      rwc(nax+jnum,1)=xs(isfr+isa,1)-xs(lssptr,1)
757      rwc(nax+jnum,2)=xs(isfr+isa,2)-xs(lssptr,2)
758      rwc(nax+jnum,3)=xs(isfr+isa,3)-xs(lssptr,3)
759  151 continue
760      call cf_pbc(0,rwc,mscr,rwx,mscr,nax,1,lssin(isa,ipss))
761      endif
762      do 5 jnum=1,lssin(isa,ipss)
763      lssptr=lssj(jsaptr+jnum)
764      jsf=ismf(lssptr)
765      isal(nax+jnum)=isfr+isa
766      jsal(nax+jnum)=lssptr
767      jmal(nax+jnum)=jsf
768      jsm=isml(lssptr)
769      jfal(nax+jnum)=2
770      if(ism.ne.jsm) jfal(nax+jnum)=3
771      if(ism.gt.0) then
772      if(jsm.gt.0) then
773      rwc(nax+jnum,1)=xsm(ism,1)-xsm(jsm,1)
774      rwc(nax+jnum,2)=xsm(ism,2)-xsm(jsm,2)
775      rwc(nax+jnum,3)=xsm(ism,3)-xsm(jsm,3)
776      else
777      rwc(nax+jnum,1)=xsm(ism,1)-xs(lssptr,1)
778      rwc(nax+jnum,2)=xsm(ism,2)-xs(lssptr,2)
779      rwc(nax+jnum,3)=xsm(ism,3)-xs(lssptr,3)
780      endif
781      else
782      if(jsm.gt.0) then
783      rwc(nax+jnum,1)=xs(isfr+isa,1)-xsm(jsm,1)
784      rwc(nax+jnum,2)=xs(isfr+isa,2)-xsm(jsm,2)
785      rwc(nax+jnum,3)=xs(isfr+isa,3)-xsm(jsm,3)
786      else
787      rwc(nax+jnum,1)=xs(isfr+isa,1)-xs(lssptr,1)
788      rwc(nax+jnum,2)=xs(isfr+isa,2)-xs(lssptr,2)
789      rwc(nax+jnum,3)=xs(isfr+isa,3)-xs(lssptr,3)
790      endif
791      endif
792    5 continue
793      if(.not.lpbc) then
794      do 7 jnum=1,lssin(isa,ipss)
795      lssptr=lssj(jsaptr+jnum)
796      xi(nax+jnum,1)=xs(isfr+isa,1)
797      xi(nax+jnum,2)=xs(isfr+isa,2)
798      xi(nax+jnum,3)=xs(isfr+isa,3)
799      xj(nax+jnum,1)=xs(lssptr,1)
800      xj(nax+jnum,2)=xs(lssptr,2)
801      xj(nax+jnum,3)=xs(lssptr,3)
802      pl(nax+jnum,1)=ps(isfr+isa,1,2)
803      pl(nax+jnum,2)=ps(isfr+isa,2,2)
804      pl(nax+jnum,3)=ps(isfr+isa,3,2)
805      pj(nax+jnum,1)=ps(lssptr,1,2)
806      pj(nax+jnum,2)=ps(lssptr,2,2)
807      pj(nax+jnum,3)=ps(lssptr,3,2)
808      isal(nax+jnum)=isfr+isa
809      jsal(nax+jnum)=lssptr
810    7 continue
811      else
812      do 8 jnum=1,lssin(isa,ipss)
813      rwc(nax+jnum,1)=rwc(nax+jnum,1)-rwx(jnum,1)
814      rwc(nax+jnum,2)=rwc(nax+jnum,2)-rwx(jnum,2)
815      rwc(nax+jnum,3)=rwc(nax+jnum,3)-rwx(jnum,3)
816      lssptr=lssj(jsaptr+jnum)
817      xi(nax+jnum,1)=xs(isfr+isa,1)
818      xi(nax+jnum,2)=xs(isfr+isa,2)
819      xi(nax+jnum,3)=xs(isfr+isa,3)
820      xj(nax+jnum,1)=xs(lssptr,1)+rwx(jnum,1)
821      xj(nax+jnum,2)=xs(lssptr,2)+rwx(jnum,2)
822      xj(nax+jnum,3)=xs(lssptr,3)+rwx(jnum,3)
823      pl(nax+jnum,1)=ps(isfr+isa,1,2)
824      pl(nax+jnum,2)=ps(isfr+isa,2,2)
825      pl(nax+jnum,3)=ps(isfr+isa,3,2)
826      pj(nax+jnum,1)=ps(lssptr,1,2)
827      pj(nax+jnum,2)=ps(lssptr,2,2)
828      pj(nax+jnum,3)=ps(lssptr,3,2)
829      isal(nax+jnum)=isfr+isa
830      jsal(nax+jnum)=lssptr
831    8 continue
832      endif
833      nax=nax+lssin(isa,ipss)
834    4 continue
835c      else
836c      do 8 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss)
837c      jsaptr=lssjpt(isa,ipss)-1
838c      ism=isl(isfr+isa,lsmol)
839c      do 9 jnum=1,lssin(isa,ipss)
840c      lssptr=lssj(jsaptr+jnum)
841c      jsa=lssptr
842c      rwx(jnum,1)=zero
843c      rwx(jnum,2)=zero
844c      rwx(jnum,3)=zero
845c      if(abs(xs(isfr+isa,1)-xs(jsa,1)).gt.boxh(1)) then
846c      rwx(jnum,1)=sign(box(1),xs(isfr+isa,1))
847c      endif
848c      if(abs(xs(isfr+isa,2)-xs(jsa,2)).gt.boxh(2)) then
849c      rwx(jnum,2)=sign(box(2),xs(isfr+isa,2))
850c      endif
851c      if(npbtyp.eq.1) then
852c      if(abs(xs(isfr+isa,3)-xs(jsa,3)).gt.boxh(3)) then
853c      rwx(jnum,3)=sign(box(3),xs(isfr+isa,3))
854c      endif
855c      endif
856c      jsf=isl(lssptr,lsfrc)
857c      isal(nax+jnum)=isfr+isa
858c      jsal(nax+jnum)=lssptr
859c      jmal(nax+jnum)=jsf
860c      jsm=isl(lssptr,lsmol)
861c      jfal(nax+jnum)=2
862c      if(ism.ne.jsm) jfal(nax+jnum)=3
863c      if(ism.gt.0) then
864c      if(jsm.gt.0) then
865c      rwc(nax+jnum,1)=xsm(ism,1)-xsm(jsm,1)-rwx(jnum,1)
866c      rwc(nax+jnum,2)=xsm(ism,2)-xsm(jsm,2)-rwx(jnum,2)
867c      rwc(nax+jnum,3)=xsm(ism,3)-xsm(jsm,3)-rwx(jnum,3)
868c      else
869c      rwc(nax+jnum,1)=xsm(ism,1)-xs(lssptr,1)-rwx(jnum,1)
870c      rwc(nax+jnum,2)=xsm(ism,2)-xs(lssptr,2)-rwx(jnum,2)
871c      rwc(nax+jnum,3)=xsm(ism,3)-xs(lssptr,3)-rwx(jnum,3)
872c      endif
873c      else
874c      if(jsm.gt.0) then
875c      rwc(nax+jnum,1)=xs(isfr+isa,1)-xsm(jsm,1)-rwx(jnum,1)
876c      rwc(nax+jnum,2)=xs(isfr+isa,2)-xsm(jsm,2)-rwx(jnum,2)
877c      rwc(nax+jnum,3)=xs(isfr+isa,3)-xsm(jsm,3)-rwx(jnum,3)
878c      else
879c      rwc(nax+jnum,1)=xs(isfr+isa,1)-xs(lssptr,1)-rwx(jnum,1)
880c      rwc(nax+jnum,2)=xs(isfr+isa,2)-xs(lssptr,2)-rwx(jnum,2)
881c      rwc(nax+jnum,3)=xs(isfr+isa,3)-xs(lssptr,3)-rwx(jnum,3)
882c      endif
883c      endif
884c    9 continue
885c      do 11 jnum=1,lssin(isa,ipss)
886c      lssptr=lssj(jsaptr+jnum)
887c      xi(nax+jnum,1)=xs(isfr+isa,1)
888c      xi(nax+jnum,2)=xs(isfr+isa,2)
889c      xi(nax+jnum,3)=xs(isfr+isa,3)
890c      xj(nax+jnum,1)=xs(lssptr,1)+rwx(jnum,1)
891c      xj(nax+jnum,2)=xs(lssptr,2)+rwx(jnum,2)
892c      xj(nax+jnum,3)=xs(lssptr,3)+rwx(jnum,3)
893c      pl(nax+jnum,1)=ps(isfr+isa,1,2)
894c      pl(nax+jnum,2)=ps(isfr+isa,2,2)
895c      pl(nax+jnum,3)=ps(isfr+isa,3,2)
896c      pj(nax+jnum,1)=ps(lssptr,1,2)
897c      pj(nax+jnum,2)=ps(lssptr,2,2)
898c      pj(nax+jnum,3)=ps(lssptr,3,2)
899c   11 continue
900c      nax=nax+lssin(isa,ipss)
901c    8 continue
902c      endif
903c
904      do 12 iax=1,nax
905      isa=isal(iax)
906      jsa=jsal(iax)
907      if(jfal(iax).eq.2) then
908      qai=qfaci*chg(isq2(isa),1,iset)
909      pai=chg(isq2(isa),2,iset)
910      qaj=qfaci*chg(isq2(jsa),1,iset)
911      paj=chg(isq2(jsa),2,iset)
912      else
913      qai=qfaci*chg(isq3(isa),1,iset)
914      pai=chg(isq3(isa),2,iset)
915      qaj=qfaci*chg(isq3(jsa),1,iset)
916      paj=chg(isq3(jsa),2,iset)
917      endif
918      rx=xj(iax,1)-xi(iax,1)
919      ry=xj(iax,2)-xi(iax,2)
920      rz=xj(iax,3)-xi(iax,3)
921      ri2=one/(rx*rx+ry*ry+rz*rz)
922      ri1=sqrt(ri2)
923      ri3=ri1*ri2
924      pix=pai*pl(iax,1)
925      piy=pai*pl(iax,2)
926      piz=pai*pl(iax,3)
927      pjx=paj*pj(iax,1)
928      pjy=paj*pj(iax,2)
929      pjz=paj*pj(iax,3)
930      rmi=three*(rx*pix+ry*piy+rz*piz)*ri2
931      rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2
932      fi(iax,1)=((rmj-qaj)*rx-pjx)*ri3
933      fi(iax,2)=((rmj-qaj)*ry-pjy)*ri3
934      fi(iax,3)=((rmj-qaj)*rz-pjz)*ri3
935      fj(iax,1)=((rmi+qai)*rx-pix)*ri3
936      fj(iax,2)=((rmi+qai)*ry-piy)*ri3
937      fj(iax,3)=((rmi+qai)*rz-piz)*ri3
938   12 continue
939c
940c     accumulate fields into solute field arrays
941c
942      iax=0
943      do 13 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss)
944      jsaptr=lssjpt(isa,ipss)-1
945      do 14 jnum=1,lssin(isa,ipss)
946      lssptr=lssj(jsaptr+jnum)
947      ps(isfr+isa,1,1)=ps(isfr+isa,1,1)+fi(iax+jnum,1)
948      ps(isfr+isa,2,1)=ps(isfr+isa,2,1)+fi(iax+jnum,2)
949      ps(isfr+isa,3,1)=ps(isfr+isa,3,1)+fi(iax+jnum,3)
950      ps(lssptr,1,1)=ps(lssptr,1,1)+fj(iax+jnum,1)
951      ps(lssptr,2,1)=ps(lssptr,2,1)+fj(iax+jnum,2)
952      ps(lssptr,3,1)=ps(lssptr,3,1)+fj(iax+jnum,3)
953   14 continue
954      iax=iax+lssin(isa,ipss)
955   13 continue
956c
957c     thermodynamic integration and perturbation
958c
959      do 15 ipset=2,3
960      if((ipset.eq.2.and.ipert2).or.
961     + (ipset.eq.3.and.ipert3)) then
962c
963      nax2=0
964      do 16 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss)
965      jsaptr=lssjpt(isa,ipss)-1
966      do 17 jnum=1,lssin(isa,ipss)
967      lssptr=lssj(jsaptr+jnum)
968      pl(nax2+jnum,1)=psp(isfr+isa,1,ipset-1,2)
969      pl(nax2+jnum,2)=psp(isfr+isa,2,ipset-1,2)
970      pl(nax2+jnum,3)=psp(isfr+isa,3,ipset-1,2)
971      pj(nax2+jnum,1)=psp(lssptr,1,ipset-1,2)
972      pj(nax2+jnum,2)=psp(lssptr,2,ipset-1,2)
973      pj(nax2+jnum,3)=psp(lssptr,3,ipset-1,2)
974   17 continue
975      nax2=nax2+lssin(isa,ipss)
976   16 continue
977c
978      if(nax2.ne.nax) call md_abort('Error in dips',me)
979c
980      do 18 iax=1,nax
981      isa=isal(iax)
982      jsa=jsal(iax)
983      if(jfal(iax).eq.2) then
984      qai=qfaci*chg(isq2(isa),1,ipset)
985      pai=chg(isq2(isa),2,ipset)
986      qaj=qfaci*chg(isq2(jsa),1,ipset)
987      paj=chg(isq2(jsa),2,ipset)
988      else
989      qai=qfaci*chg(isq3(isa),1,ipset)
990      pai=chg(isq3(isa),2,ipset)
991      qaj=qfaci*chg(isq3(jsa),1,ipset)
992      paj=chg(isq3(jsa),2,ipset)
993      endif
994      rx=xj(iax,1)-xi(iax,1)
995      ry=xj(iax,2)-xi(iax,2)
996      rz=xj(iax,3)-xi(iax,3)
997      ri2=one/(rx*rx+ry*ry+rz*rz)
998      ri1=sqrt(ri2)
999      ri3=ri1*ri2
1000      pix=pai*pl(iax,1)
1001      piy=pai*pl(iax,2)
1002      piz=pai*pl(iax,3)
1003      pjx=paj*pj(iax,1)
1004      pjy=paj*pj(iax,2)
1005      pjz=paj*pj(iax,3)
1006      rmi=three*(rx*pix+ry*piy+rz*piz)*ri2
1007      rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2
1008      fi(iax,1)=((rmj-qaj)*rx-pjx)*ri3
1009      fi(iax,2)=((rmj-qaj)*ry-pjy)*ri3
1010      fi(iax,3)=((rmj-qaj)*rz-pjz)*ri3
1011      fj(iax,1)=((rmi+qai)*rx-pix)*ri3
1012      fj(iax,2)=((rmi+qai)*ry-piy)*ri3
1013      fj(iax,3)=((rmi+qai)*rz-piz)*ri3
1014   18 continue
1015c
1016c     accumulate fields into solute field arrays
1017c
1018      iax=0
1019      do 19 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss)
1020      jsaptr=lssjpt(isa,ipss)-1
1021      do 20 jnum=1,lssin(isa,ipss)
1022      lssptr=lssj(jsaptr+jnum)
1023      psp(isfr+isa,1,ipset-1,1)=psp(isfr+isa,1,ipset-1,1)+fi(iax+jnum,1)
1024      psp(isfr+isa,2,ipset-1,1)=psp(isfr+isa,2,ipset-1,1)+fi(iax+jnum,2)
1025      psp(isfr+isa,3,ipset-1,1)=psp(isfr+isa,3,ipset-1,1)+fi(iax+jnum,3)
1026      psp(lssptr,1,ipset-1,1)=psp(lssptr,1,ipset-1,1)+fj(iax+jnum,1)
1027      psp(lssptr,2,ipset-1,1)=psp(lssptr,2,ipset-1,1)+fj(iax+jnum,2)
1028      psp(lssptr,3,ipset-1,1)=psp(lssptr,3,ipset-1,1)+fj(iax+jnum,3)
1029   20 continue
1030      iax=iax+lssin(isa,ipss)
1031   19 continue
1032c
1033      endif
1034   15 continue
1035    3 continue
1036    1 continue
1037c
1038      return
1039      end
1040