1      subroutine argos_cafe_fpsw(xs,xsm,fs,zs,ps,psp,
2     + isga,isat,isdt,ismf,isml,isss,isq1,
3     + isfrom,nums,xw,xwm,fw,pw,pwp,rtos,iwdt,lpbc,lpbcs,esw,esa,
4     + vdw,chg,iwatm,iwq,iass,lswndx,lswjpt,lswin,lswj,
5     + xi,xj,rwx,rwi1,rwi2,rwi6,rwc,f,fi,fj,facu,
6     + rw,isal,isrx,list,pl,pj)
7c
8c $Id$
9c
10      implicit none
11c
12#include "argos_cafe_common.fh"
13#include "argos_cafe_funcs_dec.fh"
14#include "bitops_decls.fh"
15c
16      real*8 xs(msa,3),xsm(msm,3),fs(msa,3,2)
17      real*8 zs(msf,3,3,2),esw(msf,mpe,2)
18      integer isga(msa),isat(msa),isdt(msa),ismf(msa)
19      integer isml(msa),isss(msa),isq1(msa)
20      real*8 xw(mwm,3,mwa),xwm(mwm,3),fw(mwm,3,mwa,2),rtos(mwm)
21      real*8 esa(nsa)
22      integer iwdt(mwm)
23      integer isfrom
24      logical lpbc,lpbcs
25c
26      real*8 vdw(mat,mat,map,mset),chg(mqt,mqp,mset)
27      integer iass(mat,mat),iwatm(mwa),iwq(mwa)
28c
29      real*8 xi(mscr,3),xj(mscr,3,mwa),rwx(mscr,3)
30      real*8 rwi1(mscr),rwi2(mscr),rwi6(mscr),rw(mscr),rwc(mscr,3)
31      real*8 f(mscr),fi(mscr,3,mwa),fj(mscr,3,mwa),facu(mscr)
32      integer isal(mscr),isrx(mscr)
33c
34      integer lswj(*)
35      integer nums,i
36      integer lswndx(0:msa,2),lswjpt(nums,2),lswin(nums,2)
37      integer list(0:msa)
38c
39      real*8 ps(msa,3,2),psp(msa,3,2,2)
40      real*8 pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2)
41      real*8 pl(mscr,3),pj(mscr,3,mwa)
42c
43      integer isatm,nswlen(2)
44      integer isfr,iwm,ipsw,number,isa,ispm,isf,nax,ism
45      integer ispj,ismn,lswptr,iwa,iax,iwatmi,ix,iy
46      integer iwatyp
47      real*8 c6,cf6,c12,cf12,sumen
48      real*8 c64,c124,dercon,qj,qj4,derco1,derco2
49      real*8 drvco1,drvco2,derco3,drvco3,c6p,c12p,etermq,eterml
50      real*8 qi,qai,qaj,pai,paj,pix,piy,piz,pjx,pjy,pjz
51      real*8 rx,ry,rz,ri1,ri2,ri3,rmi,rmj,fri,fmi,fmj,rmm
52      real*8 zxx,zxy,zxz,zyx,zyy,zyz,zzx,zzy,zzz
53      real*8 eswqsm,eswpsm,qfaci
54      real*8 rtmp
55c
56#include "argos_cafe_funcs_sfn.fh"
57#include "bitops_funcs.fh"
58c
59      etermq=zero
60c
61c
62cx new stuff end
63c     this subroutine evaluates the solute-solvent forces for nums
64c     solute atoms starting from isfrom. the interacting solvent
65c     molecules are determined from the pairlist.
66c
67      isfr=isfrom-1
68c
69      if(nrwrec.gt.0) then
70      do 1 iwm=1,mwm
71      rtos(iwm)=zero
72    1 continue
73      endif
74c
75      qfaci=one/qfac
76c
77      do 2 ipsw=1,lpsw
78c
79c     evaluate outer index array
80c
81      nswlen(ipsw)=0
82      lswndx(0,ipsw)=0
83      number=0
84      do 3 isa=1,nums
85      if(number+lswin(isa,ipsw).gt.mscr .or.
86     + (ismf(isfr+isa).ne.ismf(isfr+isa-1).and.
87     + number.gt.0)) then
88      nswlen(ipsw)=nswlen(ipsw)+1
89      lswndx(nswlen(ipsw),ipsw)=isa-1
90      number=0
91      endif
92      number=number+lswin(isa,ipsw)
93    3 continue
94      if(number.gt.0) then
95      nswlen(ipsw)=nswlen(ipsw)+1
96      lswndx(nswlen(ipsw),ipsw)=nums
97      endif
98c
99      do 4 ispm=1,nswlen(ipsw)
100      isf=ismf(isfr+lswndx(ispm,ipsw))
101      do 5 isa=0,nums
102      list(isa)=0
103    5 continue
104      nax=0
105c
106      do 6 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
107      ispj=lswjpt(isa,ipsw)-1
108      ism=isml(isfr+isa)
109      if(lpbc.or.lpbcs.or.ism.eq.0) then
110      do 7 ismn=1,lswin(isa,ipsw)
111      lswptr=lswj(ispj+ismn)
112      rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1)
113      rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2)
114      rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3)
115      isrx(nax+ismn)=0
116    7 continue
117      if(lpbc.or.lpbcs)
118     + call argos_cafe_pbc(0,rwc,mscr,rwx,mscr,nax,1,lswin(isa,ipsw))
119      endif
120      if(ism.gt.0) then
121      do 8 ismn=1,lswin(isa,ipsw)
122      lswptr=lswj(ispj+ismn)
123      rwc(nax+ismn,1)=xsm(ism,1)-xwm(lswptr,1)
124      rwc(nax+ismn,2)=xsm(ism,2)-xwm(lswptr,2)
125      rwc(nax+ismn,3)=xsm(ism,3)-xwm(lswptr,3)
126    8 continue
127      endif
128c
129c      if(lssscl) then
130c      isrst=iand(isss(isfr+isa),3)
131c      isatm=isat(isfr+isa)
132c      do 9 iwa=1,mwa
133c      iasst=iass(isatm,iwatm(iwa))
134c      if(iasst.le.0.or.iasst.ge.3.or.isrst.ne.iasst) isrst=0
135c    9 continue
136c      do 10 ismn=1,lswin(isa,ipsw)
137c      isrx(nax+ismn)=isrst
138c   10 continue
139c      endif
140c
141c      write(*,'(4i5,2f12.6)')
142c     + lssscl,isga(isa),isss(isfr+isa),iand(isss(isfr+isa),6),
143c     + shift0(1),shift1(1)
144      if(lssscl) then
145      do 10 ismn=1,lswin(isa,ipsw)
146c      isrx(nax+ismn)=isss(isfr+isa)
147      if(iand(isss(isfr+isa),6).eq.2) isrx(nax+ismn)=-1
148      if(iand(isss(isfr+isa),6).eq.4) isrx(nax+ismn)=1
149   10 continue
150      endif
151c
152      if(iand(isdt(isfr+isa),mdynam).eq.ldynam) then
153      do 11 ismn=1,lswin(isa,ipsw)
154      lswptr=lswj(ispj+ismn)
155      xi(nax+ismn,1)=xs(isfr+isa,1)
156      xi(nax+ismn,2)=xs(isfr+isa,2)
157      xi(nax+ismn,3)=xs(isfr+isa,3)
158      pl(nax+ismn,1)=ps(isfr+isa,1,1)
159      pl(nax+ismn,2)=ps(isfr+isa,2,1)
160      pl(nax+ismn,3)=ps(isfr+isa,3,1)
161      isal(nax+ismn)=isfr+isa
162c      if(iand(iwdt(lswptr),mdynam).ne.ldynam) then
163c      facu(nax+ismn)=half
164c      else
165      facu(nax+ismn)=one
166c      endif
167c      if(includ.eq.1) facu(nax+ismn)=one
168   11 continue
169      else
170      do 12 ismn=1,lswin(isa,ipsw)
171      lswptr=lswj(ispj+ismn)
172      xi(nax+ismn,1)=xs(isfr+isa,1)
173      xi(nax+ismn,2)=xs(isfr+isa,2)
174      xi(nax+ismn,3)=xs(isfr+isa,3)
175      pl(nax+ismn,1)=ps(isfr+isa,1,1)
176      pl(nax+ismn,2)=ps(isfr+isa,2,1)
177      pl(nax+ismn,3)=ps(isfr+isa,3,1)
178      isal(nax+ismn)=isfr+isa
179      if(iand(iwdt(lswptr),mdynam).eq.ldynam) then
180      facu(nax+ismn)=one
181      else
182      facu(nax+ismn)=zero
183      endif
184      if(includ.eq.1) facu(nax+ismn)=one
185   12 continue
186      endif
187c
188      if(.not.lpbc.and..not.lpbcs) then
189      do 13 iwa=1,mwa
190      do 14 ismn=1,lswin(isa,ipsw)
191      lswptr=lswj(ispj+ismn)
192      xj(nax+ismn,1,iwa)=xw(lswptr,1,iwa)
193      xj(nax+ismn,2,iwa)=xw(lswptr,2,iwa)
194      xj(nax+ismn,3,iwa)=xw(lswptr,3,iwa)
195      pj(nax+ismn,1,iwa)=pw(lswptr,1,iwa,1)
196      pj(nax+ismn,2,iwa)=pw(lswptr,2,iwa,1)
197      pj(nax+ismn,3,iwa)=pw(lswptr,3,iwa,1)
198   14 continue
199   13 continue
200      else
201      do 15 ismn=1,lswin(isa,ipsw)
202      rwc(nax+ismn,1)=rwc(nax+ismn,1)-rwx(ismn,1)
203      rwc(nax+ismn,2)=rwc(nax+ismn,2)-rwx(ismn,2)
204      rwc(nax+ismn,3)=rwc(nax+ismn,3)-rwx(ismn,3)
205   15 continue
206      do 16 iwa=1,mwa
207      do 17 ismn=1,lswin(isa,ipsw)
208      lswptr=lswj(ispj+ismn)
209      xj(nax+ismn,1,iwa)=xw(lswptr,1,iwa)+rwx(ismn,1)
210      xj(nax+ismn,2,iwa)=xw(lswptr,2,iwa)+rwx(ismn,2)
211      xj(nax+ismn,3,iwa)=xw(lswptr,3,iwa)+rwx(ismn,3)
212      pj(nax+ismn,1,iwa)=pw(lswptr,1,iwa,1)
213      pj(nax+ismn,2,iwa)=pw(lswptr,2,iwa,1)
214      pj(nax+ismn,3,iwa)=pw(lswptr,3,iwa,1)
215   17 continue
216   16 continue
217      endif
218c
219      nax=nax+lswin(isa,ipsw)
220      list(isa)=nax
221    6 continue
222c
223      do 22 iax=1,nax
224      fi(iax,1,1)=zero
225      fi(iax,2,1)=zero
226      fi(iax,3,1)=zero
227   22 continue
228      do 23 iwa=1,mwa
229      do 24 iax=1,nax
230      fj(iax,1,iwa)=zero
231      fj(iax,2,iwa)=zero
232      fj(iax,3,iwa)=zero
233   24 continue
234   23 continue
235c      if(npener.ne.0) then
236c      do 25 iax=1,nax
237c      u(iax)=zero
238c   25 continue
239c      endif
240      do 26 iwa=1,mwa
241      do 27 iax=1,nax
242      f(iax)=zero
243      rwx(iax,1)=xi(iax,1)-xj(iax,1,iwa)
244      rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
245      rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
246      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
247      rtmp=rwi2(iax)
248      if(isrx(iax).gt.0) rwi2(iax)=one/(one/rwi2(iax)+shift0(1))
249      if(isrx(iax).lt.0) rwi2(iax)=one/(one/rwi2(iax)+shift1(1))
250c      write(*,'(3i5,2f12.6)')
251c     + isga(isal(iax)),isal(iax),isrx(iax),rtmp,rwi2(iax)
252   27 continue
253c
254c     Lennard-Jones interactions
255c
256      iwatmi=iwatm(iwa)
257      eterml=zero
258      do 28 iax=1,nax
259      isa=isal(iax)
260      isatm=isat(isa)
261      c6=vdw(isatm,iwatmi,1,iset)
262      cf6=six*c6
263      c12=vdw(isatm,iwatmi,3,iset)
264      cf12=twelve*c12
265      rwi6(iax)=rwi2(iax)*rwi2(iax)*rwi2(iax)
266      rw(iax)=facu(iax)*(c12*rwi6(iax)-c6)*rwi6(iax)
267      eterml=eterml+rw(iax)
268      if(npener.ne.0) then
269      esa(isga(isa))=esa(isga(isa))+half*rw(iax)
270      endif
271      f(iax)=f(iax)+(cf12*rwi6(iax)-cf6)*rwi6(iax)*rwi2(iax)
272   28 continue
273      esw(isf,5,ipsw)=esw(isf,5,ipsw)+eterml
274c
275      qj=chg(iwq(iwa),1,iset)
276      qaj=qfaci*qj
277c      dqj=qwa(iwa,4)
278c      dqaj=qfaci*dqj
279      paj=chg(iwq(iwa),2,iset)
280c      dpaj=pwa(iwa,4)
281      eswqsm=zero
282      eswpsm=zero
283c      dswqsm=zero
284c      dswqws=zero
285c      dswqps=zero
286c      dswpss=zero
287c      dswpws=zero
288      do 21 iax=1,nax
289      isa=isal(iax)
290      qi=chg(isq1(isa),1,iset)
291c      dqi=qsa(isa,4,1)
292      qai=qfaci*qi
293c      dqai=qfaci*dqi
294      pai=chg(isq1(isa),2,iset)
295c      dpai=psa(isa,4)
296      pix=pai*pl(iax,1)
297      piy=pai*pl(iax,2)
298      piz=pai*pl(iax,3)
299      pjx=paj*pj(iax,1,iwa)
300      pjy=paj*pj(iax,2,iwa)
301      pjz=paj*pj(iax,3,iwa)
302      rwx(iax,1)=xi(iax,1)-xj(iax,1,iwa)
303      rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
304      rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
305      rx=-rwx(iax,1)
306      ry=-rwx(iax,2)
307      rz=-rwx(iax,3)
308      rwi2(iax)=one/(rx**2+ry**2+rz**2)
309      rwi1(iax)=sqrt(rwi2(iax))
310      ri1=rwi1(iax)
311      ri2=rwi2(iax)
312      ri3=qfac*qfac*ri1*ri2
313      rmi=three*(rx*pix+ry*piy+rz*piz)*ri2
314      rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2
315      if(ipolt.eq.1) then
316      fri=((-qai)*qaj+qai*rmj-qaj*rmi)*ri3
317      fmi=(qaj)*ri3
318      fmj=(-qai)*ri3
319      else
320      rmm=three*(pix*pjx+piy*pjy+piz*pjz)*ri2
321      fri=((-qai)*qaj+qai*rmj-qaj*rmi+5.0*rmi*rmj/three-rmm)*ri3
322      fmi=(qaj-rmj)*ri3
323      fmj=((-qai)-rmi)*ri3
324      endif
325      fi(iax,1,1)=fi(iax,1,1)+fri*rx+fmi*pix+fmj*pjx
326      fi(iax,2,1)=fi(iax,2,1)+fri*ry+fmi*piy+fmj*pjy
327      fi(iax,3,1)=fi(iax,3,1)+fri*rz+fmi*piz+fmj*pjz
328      fj(iax,1,iwa)=fj(iax,1,iwa)-(fri*rx+fmi*pix+fmj*pjx)
329      fj(iax,2,iwa)=fj(iax,2,iwa)-(fri*ry+fmi*piy+fmj*pjy)
330      fj(iax,3,iwa)=fj(iax,3,iwa)-(fri*rz+fmi*piz+fmj*pjz)
331      eswqsm=eswqsm+qi*facu(iax)*ri1
332      eswpsm=eswpsm+facu(iax)*(qai*rmj-qaj*rmi)*ri1
333c      if(ithint.ne.0) then
334c      dpix=dpai*pl(iax,1)
335c      dpiy=dpai*pl(iax,2)
336c      dpiz=dpai*pl(iax,3)
337c      dpjx=dpaj*pj(iax,1,jwa)
338c      dpjy=dpaj*pj(iax,2,jwa)
339c      dpjz=dpaj*pj(iax,3,jwa)
340c      drmi=three*(rx*dpix+ry*dpiy+rz*dpiz)*ri2
341c      drmj=three*(rx*dpjx+ry*dpjy+rz*dpjz)*ri2
342c      dswqsm=dswqsm+dqi*facu(iax)*ri1
343c      dswqws=dswqws+drmj*ri1
344c      dswqps=dswqps+dqai*rmj*ri1
345c      dswpss=dswpss-drmi*ri1
346c      dswpws=dswpws+qai*drmj*ri1
347c      endif
348      zxx=(-half)*(fri*rx+fmi*pix+fmj*pjx)*rwc(iax,1)
349      zyx=(-half)*(fri*rx+fmi*pix+fmj*pjx)*rwc(iax,2)
350      zzx=(-half)*(fri*rx+fmi*pix+fmj*pjx)*rwc(iax,3)
351      zxy=(-half)*(fri*ry+fmi*piy+fmj*pjy)*rwc(iax,1)
352      zyy=(-half)*(fri*ry+fmi*piy+fmj*pjy)*rwc(iax,2)
353      zzy=(-half)*(fri*ry+fmi*piy+fmj*pjy)*rwc(iax,3)
354      zxz=(-half)*(fri*rz+fmi*piz+fmj*pjz)*rwc(iax,1)
355      zyz=(-half)*(fri*rz+fmi*piz+fmj*pjz)*rwc(iax,2)
356      zzz=(-half)*(fri*rz+fmi*piz+fmj*pjz)*rwc(iax,3)
357      zw(1,1,ipsw)=zw(1,1,ipsw)+zxx
358      zw(2,1,ipsw)=zw(2,1,ipsw)+zyx
359      zw(3,1,ipsw)=zw(3,1,ipsw)+zzx
360      zw(1,2,ipsw)=zw(1,2,ipsw)+zxy
361      zw(2,2,ipsw)=zw(2,2,ipsw)+zyy
362      zw(3,2,ipsw)=zw(3,2,ipsw)+zzy
363      zw(1,3,ipsw)=zw(1,3,ipsw)+zxz
364      zw(2,3,ipsw)=zw(2,3,ipsw)+zyz
365      zw(3,3,ipsw)=zw(3,3,ipsw)+zzz
366      zs(isf,1,1,ipsw)=zs(isf,1,1,ipsw)+zxx
367      zs(isf,2,1,ipsw)=zs(isf,2,1,ipsw)+zyx
368      zs(isf,3,1,ipsw)=zs(isf,3,1,ipsw)+zzx
369      zs(isf,1,2,ipsw)=zs(isf,1,2,ipsw)+zxy
370      zs(isf,2,2,ipsw)=zs(isf,2,2,ipsw)+zyy
371      zs(isf,3,2,ipsw)=zs(isf,3,2,ipsw)+zzy
372      zs(isf,1,3,ipsw)=zs(isf,1,3,ipsw)+zxz
373      zs(isf,2,3,ipsw)=zs(isf,2,3,ipsw)+zyz
374      zs(isf,3,3,ipsw)=zs(isf,3,3,ipsw)+zzz
375   21 continue
376c
377      do 35 iax=1,nax
378      fi(iax,1,1)=fi(iax,1,1)+f(iax)*rwx(iax,1)
379      fi(iax,2,1)=fi(iax,2,1)+f(iax)*rwx(iax,2)
380      fi(iax,3,1)=fi(iax,3,1)+f(iax)*rwx(iax,3)
381      fj(iax,1,iwa)=fj(iax,1,iwa)-f(iax)*rwx(iax,1)
382      fj(iax,2,iwa)=fj(iax,2,iwa)-f(iax)*rwx(iax,2)
383      fj(iax,3,iwa)=fj(iax,3,iwa)-f(iax)*rwx(iax,3)
384   35 continue
385      do 136 iy=1,3
386      do 36 ix=1,3
387      sumen=zero
388      do 37 iax=1,nax
389      sumen=sumen-half*f(iax)*rwx(iax,iy)*rwc(iax,ix)
390   37 continue
391      zs(isf,ix,iy,ipsw)=zs(isf,ix,iy,ipsw)+sumen
392      zw(ix,iy,ipsw)=zw(ix,iy,ipsw)+sumen
393   36 continue
394  136 continue
395c
396c     Radial distribution functions
397c
398c      if(ifstep-1.eq.((ifstep-1)/nfrdf)*nfrdf .and. ngrsw.gt.0) then
399c      do 38 igc=1,ngc
400c      if(ngt(igc).eq.2) then
401c      if(iagc(igc).eq.iwa) then
402c      igr=igrc(igc)
403c      do 39 iax=1,nax
404c      if(isga(isal(iax)).eq.jagc(igc)) then
405c      indx=int(one/(rwi1(iax)*drdf))
406c      if(indx.gt.ngl) indx=ngl
407c      rdf(indx,igr)=rdf(indx,igr)+rdfvol
408c      endif
409c   39 continue
410c      endif
411c      endif
412c   38 continue
413c      endif
414c
415c     Thermodynamic integration
416c
417      if(ithint) then
418      if(ith(2).or.ith(14)) then
419      if(.not.lssscl) then
420      do 40 iax=1,nax
421      isa=isal(iax)
422      isatm=isat(isa)
423      c64=vdw(isatm,iwatm(iwa),1,4)
424      c124=vdw(isatm,iwatm(iwa),3,4)
425      dercon=half*(c124*rwi6(iax)-c64)*rwi6(iax)
426      deriv(3,ipsw)=deriv(3,ipsw)+dercon
427      deriv(14,ipsw)=deriv(14,ipsw)+dercon
428   40 continue
429      else
430      do 41 iax=1,nax
431      isa=isal(iax)
432      isatm=isat(isa)
433      c64=vdw(isatm,iwatm(iwa),1,4)
434      c124=vdw(isatm,iwatm(iwa),3,4)
435      dercon=half*(c124*rwi6(iax)-c64)*rwi6(iax)
436      if(isrx(iax).gt.0) then
437      c64=half*three*vdw(isatm,iwatm(iwa),1,iset)
438      c124=three*vdw(isatm,iwatm(iwa),3,iset)
439      dercon=dercon+shift0(4)*rwi2(iax)*rwi6(iax)*(c64-c124*rwi6(iax))
440      elseif(isrx(iax).lt.0) then
441      c64=half*three*vdw(isatm,iwatm(iwa),1,iset)
442      c124=three*vdw(isatm,iwatm(iwa),3,iset)
443      dercon=dercon+shift1(4)*rwi2(iax)*rwi6(iax)*(c64-c124*rwi6(iax))
444      else
445      c64=vdw(isatm,iwatm(iwa),1,4)
446      c124=vdw(isatm,iwatm(iwa),3,4)
447      dercon=half*(c124*rwi6(iax)-c64)*rwi6(iax)
448      endif
449      deriv(3,ipsw)=deriv(3,ipsw)+dercon
450      deriv(14,ipsw)=deriv(14,ipsw)+dercon
451   41 continue
452      endif
453      endif
454      if(ith(4).or.ith(16)) then
455      qj=chg(iwq(iwa),1,iset)
456      qj4=chg(iwq(iwa),1,4)
457      derco1=zero
458      derco2=zero
459      if(ipme.eq.0) then
460      if(.not.lssscl) then
461      do 42 iax=1,nax
462      isa=isal(iax)
463      drvco1=qj*chg(isq1(isa),1,4)*rwi1(iax)
464      derco1=derco1+drvco1
465      drvco2=chg(isq1(isa),1,iset)*qj4*rwi1(iax)
466      derco2=derco2+drvco2
467   42 continue
468      deriv(5,ipsw)=deriv(5,ipsw)+derco1
469      deriv(16,ipsw)=deriv(16,ipsw)+derco2
470      else
471      derco3=zero
472      do 43 iax=1,nax
473      isa=isal(iax)
474      drvco1=qj*chg(isq1(isa),1,4)*rwi1(iax)
475      derco1=derco1+drvco1
476      drvco2=chg(isq1(isa),1,iset)*qj4*rwi1(iax)
477      derco2=derco2+drvco2
478      drvco3=zero
479      if(isrx(iax).gt.0) then
480      drvco3=(-half)*shift0(4)*chg(isq1(isa),1,iset)*
481     + qj*rwi1(iax)*rwi2(iax)
482      elseif(isrx(iax).lt.0) then
483      drvco3=(-half)*shift1(4)*chg(isq1(isa),1,iset)*
484     + qj*rwi1(iax)*rwi2(iax)
485      endif
486      derco3=derco3+drvco3
487   43 continue
488      deriv(5,ipsw)=deriv(5,ipsw)+derco1+half*derco3
489      deriv(16,ipsw)=deriv(16,ipsw)+derco2+half*derco3
490      endif
491      else
492      if(.not.lssscl) then
493      do 142 iax=1,nax
494      isa=isal(iax)
495      drvco1=qj*chg(isq1(isa),1,4)*rwi1(iax)
496      derco1=derco1+drvco1
497      drvco2=chg(isq1(isa),1,iset)*qj4*rwi1(iax)
498      derco2=derco2+drvco2
499  142 continue
500      deriv(5,ipsw)=deriv(5,ipsw)+derco1
501      deriv(16,ipsw)=deriv(16,ipsw)+derco2
502      else
503      derco3=zero
504      do 143 iax=1,nax
505      isa=isal(iax)
506      drvco1=qj*chg(isq1(isa),1,4)*rwi1(iax)
507      derco1=derco1+drvco1
508      drvco2=chg(isq1(isa),1,iset)*qj4*rwi1(iax)
509      derco2=derco2+drvco2
510      drvco3=zero
511      if(isrx(iax).gt.0) then
512      drvco3=(-half)*shift0(4)*chg(isq1(isa),1,iset)*
513     + qj*rwi1(iax)*rwi2(iax)
514      elseif(isrx(iax).lt.0) then
515      drvco3=(-half)*shift1(4)*chg(isq1(isa),1,iset)*
516     + qj*rwi1(iax)*rwi2(iax)
517      endif
518      derco3=derco3+drvco3
519  143 continue
520      deriv(5,ipsw)=deriv(5,ipsw)+derco1+half*derco3
521      deriv(16,ipsw)=deriv(16,ipsw)+derco2+half*derco3
522      endif
523      endif
524      endif
525      endif
526c
527c     Thermodynamic perturbation 1
528c
529      if(ipert2) then
530      if(ip2(2).or.ip2(14)) then
531      iwatyp=iwatm(iwa)
532      if(.not.lssscl) then
533      do 44 iax=1,nax
534      isa=isal(iax)
535      c6p=vdw(isat(isa),iwatyp,1,2)
536      c12p=vdw(isat(isa),iwatyp,3,2)
537      ep2(ipsw)=ep2(ipsw)+facu(iax)*(c12p*rwi6(iax)-c6p)*rwi6(iax)
538   44 continue
539      else
540      do 45 iax=1,nax
541      isa=isal(iax)
542      c6p=vdw(isat(isa),iwatyp,1,2)
543      c12p=vdw(isat(isa),iwatyp,3,2)
544      if(isrx(iax).gt.0) then
545      rwi6(iax)=(one/(one/rwi2(iax)-shift0(1)+shift0(2)))**3
546      elseif(isrx(iax).lt.0) then
547      rwi6(iax)=(one/(one/rwi2(iax)-shift1(1)+shift1(2)))**3
548      else
549      rwi6(iax)=rwi2(iax)**3
550      endif
551      ep2(ipsw)=ep2(ipsw)+facu(iax)*(c12p*rwi6(iax)-c6p)*rwi6(iax)
552   45 continue
553      endif
554      ep2(ipsw)=ep2(ipsw)-eterml
555      endif
556      if(ip2(4).or.ip2(5).or.ip2(16).or.ip2(17)) then
557      qj=chg(iwq(iwa),1,2)
558      if(ipme.eq.0) then
559      if(.not.lssscl) then
560      do 46 iax=1,nax
561      isa=isal(iax)
562      rwx(iax,1)=xi(iax,1)-xj(iax,1,iwa)
563      rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
564      rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
565      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
566      rwi1(iax)=sqrt(rwi2(iax))
567      ep2(ipsw)=ep2(ipsw)+facu(iax)*chg(isq1(isa),1,2)*qj*rwi1(iax)
568   46 continue
569      else
570      do 47 iax=1,nax
571      isa=isal(iax)
572      rwx(iax,1)=xi(iax,1)-xj(iax,1,iwa)
573      rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
574      rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
575      rwi6(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
576      if(isrx(iax).gt.0) then
577      rwi6(iax)=one/(one/rwi6(iax)+shift0(2))
578      elseif(isrx(iax).lt.0) then
579      rwi6(iax)=one/(one/rwi6(iax)+shift1(2))
580      endif
581      rwi1(iax)=sqrt(rwi6(iax))
582      ep2(ipsw)=ep2(ipsw)+facu(iax)*chg(isq1(isa),1,2)*qj*rwi1(iax)
583   47 continue
584      endif
585      else
586      if(.not.lssscl) then
587      do 146 iax=1,nax
588      isa=isal(iax)
589      rwx(iax,1)=xi(iax,1)-xj(iax,1,iwa)
590      rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
591      rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
592      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
593      rwi1(iax)=sqrt(rwi2(iax))
594      ep2(ipsw)=ep2(ipsw)+facu(iax)*erfc(ealpha/rwi1(iax))*
595     + chg(isq1(isa),1,2)*qj*rwi1(iax)
596  146 continue
597      else
598      do 147 iax=1,nax
599      isa=isal(iax)
600      rwx(iax,1)=xi(iax,1)-xj(iax,1,iwa)
601      rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
602      rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
603      rwi6(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
604      if(isrx(iax).gt.0) then
605      rwi6(iax)=one/(one/rwi6(iax)+shift0(2))
606      elseif(isrx(iax).lt.0) then
607      rwi6(iax)=one/(one/rwi6(iax)+shift1(2))
608      endif
609      rwi1(iax)=sqrt(rwi6(iax))
610      ep2(ipsw)=ep2(ipsw)+facu(iax)*erfc(ealpha/rwi1(iax))*
611     + chg(isq1(isa),1,2)*qj*rwi1(iax)
612  147 continue
613      endif
614      endif
615      ep2(ipsw)=ep2(ipsw)-etermq
616      endif
617      endif
618c
619c     Thermodynamic perturbation 2
620c
621      if(ipert3) then
622      if(ip3(2).or.ip3(14)) then
623      iwatyp=iwatm(iwa)
624      if(.not.lssscl) then
625      do 48 iax=1,nax
626      isa=isal(iax)
627      c6p=vdw(isat(isa),iwatyp,1,3)
628      c12p=vdw(isat(isa),iwatyp,3,3)
629      ep3(ipsw)=ep3(ipsw)+facu(iax)*(c12p*rwi6(iax)-c6p)*rwi6(iax)
630   48 continue
631      else
632      do 49 iax=1,nax
633      isa=isal(iax)
634      c6p=vdw(isat(isa),iwatyp,1,3)
635      c12p=vdw(isat(isa),iwatyp,3,3)
636      if(isrx(iax).gt.0) then
637      rwi6(iax)=(one/(one/rwi2(iax)-shift0(1)+shift0(3)))**3
638      elseif(isrx(iax).lt.0) then
639      rwi6(iax)=(one/(one/rwi2(iax)-shift1(1)+shift1(3)))**3
640      else
641      rwi6(iax)=rwi2(iax)**3
642      endif
643      ep3(ipsw)=ep3(ipsw)+facu(iax)*(c12p*rwi6(iax)-c6p)*rwi6(iax)
644   49 continue
645      endif
646      ep3(ipsw)=ep3(ipsw)-eterml
647      endif
648      if(ip2(4).or.ip2(5).or.ip2(16).or.ip2(17)) then
649      qj=chg(iwq(iwa),1,3)
650      if(ipme.eq.0) then
651      if(.not.lssscl) then
652      do 50 iax=1,nax
653      isa=isal(iax)
654      rwx(iax,1)=xi(iax,1)-xj(iax,1,iwa)
655      rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
656      rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
657      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
658      rwi1(iax)=sqrt(rwi2(iax))
659      ep3(ipsw)=ep3(ipsw)+facu(iax)*chg(isq1(isa),1,3)*qj*rwi1(iax)
660   50 continue
661      else
662      do 51 iax=1,nax
663      isa=isal(iax)
664      rwx(iax,1)=xi(iax,1)-xj(iax,1,iwa)
665      rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
666      rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
667      rwi6(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
668      if(isrx(iax).gt.0) then
669      rwi6(iax)=one/(one/rwi6(iax)+shift0(3))
670      elseif(isrx(iax).lt.0) then
671      rwi6(iax)=one/(one/rwi6(iax)+shift1(3))
672      endif
673      rwi1(iax)=sqrt(rwi6(iax))
674      ep3(ipsw)=ep3(ipsw)+facu(iax)*chg(isq1(isa),1,3)*qj*rwi1(iax)
675   51 continue
676      endif
677      else
678      if(.not.lssscl) then
679      do 150 iax=1,nax
680      isa=isal(iax)
681      rwx(iax,1)=xi(iax,1)-xj(iax,1,iwa)
682      rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
683      rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
684      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
685      rwi1(iax)=sqrt(rwi2(iax))
686      ep3(ipsw)=ep3(ipsw)+facu(iax)*erfc(ealpha/rwi1(iax))*
687     + chg(isq1(isa),1,3)*qj*rwi1(iax)
688  150 continue
689      else
690      do 151 iax=1,nax
691      isa=isal(iax)
692      rwx(iax,1)=xi(iax,1)-xj(iax,1,iwa)
693      rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
694      rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
695      rwi6(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
696      if(isrx(iax).gt.0) then
697      rwi6(iax)=one/(one/rwi6(iax)+shift0(3))
698      elseif(isrx(iax).lt.0) then
699      rwi6(iax)=one/(one/rwi6(iax)+shift1(3))
700      endif
701      rwi1(iax)=sqrt(rwi6(iax))
702      ep3(ipsw)=ep3(ipsw)+facu(iax)*erfc(ealpha/rwi1(iax))*
703     + chg(isq1(isa),1,3)*qj*rwi1(iax)
704  151 continue
705      endif
706      endif
707      ep3(ipsw)=ep3(ipsw)-etermq
708      endif
709      endif
710   26 continue
711c
712      iax=0
713      do 52 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw)
714      ispj=lswjpt(isa,ipsw)-1
715      do 53 ismn=1,lswin(isa,ipsw)
716      fs(isfr+isa,1,ipsw)=fs(isfr+isa,1,ipsw)+fi(iax+ismn,1,1)
717      fs(isfr+isa,2,ipsw)=fs(isfr+isa,2,ipsw)+fi(iax+ismn,2,1)
718      fs(isfr+isa,3,ipsw)=fs(isfr+isa,3,ipsw)+fi(iax+ismn,3,1)
719   53 continue
720      do 54 iwa=1,mwa
721      do 55 ismn=1,lswin(isa,ipsw)
722      lswptr=lswj(ispj+ismn)
723      fw(lswptr,1,iwa,ipsw)=fw(lswptr,1,iwa,ipsw)+fj(iax+ismn,1,iwa)
724      fw(lswptr,2,iwa,ipsw)=fw(lswptr,2,iwa,ipsw)+fj(iax+ismn,2,iwa)
725      fw(lswptr,3,iwa,ipsw)=fw(lswptr,3,iwa,ipsw)+fj(iax+ismn,3,iwa)
726   55 continue
727c
728      if(nrwrec.gt.0) then
729      do 56 ismn=1,lswin(isa,ipsw)
730      lswptr=lswj(ispj+ismn)
731      if(rtos(lswptr).lt.rwi2(iax+ismn)) rtos(lswptr)=rwi2(iax+ismn)
732   56 continue
733      endif
734   54 continue
735c
736c      if(npener.ne.0) then
737c      do 57 ismn=1,lswin(isa,ipsw)
738c      lswptr=lswj(ispj+ismn)
739c      uwms(lswptr)=uwms(lswptr)+u(iax+ismn)
740c   57 continue
741c      endif
742c
743      iax=iax+lswin(isa,ipsw)
744   52 continue
745    4 continue
746    2 continue
747c
748      return
749      end
750