1      subroutine argos_cafe_fss(xs,xsm,fs,zs,
2     + isga,isat,isdt,ismf,isml,isss,isq2,isq3,isgm,
3     + isfrom,nums,lpbc,lpbcs,ess,fss,esa,
4     + vdw,chg,iass,
5     + lssndx,lssjpt,lssin,lssj,
6     + xi,xj,rwx,rwi1,rwi2,rwi6,rwc,f,fi,fj,facu,
7     + rw,isal,jsal,jmal,jfal,isrx,qsa2,qsa3,dera,lda,rda,uda,lseq)
8c
9c $Id$
10c
11      implicit none
12c
13#include "argos_cafe_common.fh"
14#include "mafdecls.fh"
15c
16      real*8 rtmp
17      real*8 xs(msa,3),xsm(msm,3),fs(msa,3,2)
18      real*8 zs(msf,3,3,2),ess(msf,msf,mpe,2)
19      real*8 fss(msf,msf,3,2)
20      real*8 esa(nsa)
21      integer isga(msa),isat(msa),isdt(msa),ismf(msa)
22      integer isml(msa),isss(msa),isq2(msa),isq3(msa)
23      integer isgm(msa),lseq(mseq)
24c
25      real*8 vdw(mat,mat,map,mset),chg(mqt,mqp,mset)
26      logical lpbc,lpbcs
27c
28      real*8 xi(mscr,3),xj(mscr,3),rwx(mscr,3),rwi1(mscr)
29      real*8 rwi2(mscr),rwi6(mscr),rwc(mscr,3),rw(mscr)
30      real*8 f(mscr),fi(mscr,3),fj(mscr,3)
31      real*8 qsa2(mscr),qsa3(mscr)
32      integer isal(mscr),jsal(mscr),jmal(mscr),jfal(mscr),isrx(mscr)
33      integer lssj(*)
34      real*8 facu(mscr)
35      integer nums
36      integer lssndx(0:msa,2),lssjpt(nums,2),lssin(nums,2)
37      integer iass(mat,mat),nsslen(2)
38c
39      real*8 dera(6,nsatot)
40      integer lda(16,*)
41      real*8 rda(11,*),uda(4,*)
42c
43      integer isa,jsa,i,isf,jsf,ix
44      integer isfr,isfrom,ism,jsm
45      integer ipss,number,isslen,nax,jsaptr
46      integer jnum,lssptr,iax
47      real*8 dercon
48c
49      real*8 c6,c12,cf6,cf12
50      real*8 c64,c124
51      real*8 q14,sumen1,sumen2,sumen3
52      real*8 etermq,eterml
53      integer istt,jstt
54      real*8 q,ferfc,fderfc
55c
56#include "argos_cafe_funcs_dec.fh"
57#include "bitops_decls.fh"
58#include "argos_cafe_funcs_sfn.fh"
59#include "bitops_funcs.fh"
60c
61c
62      if(nfhop.eq.0) then
63      do 112 i=1,msa
64      if(isq2(i).le.0.or.isq3(i).le.0.or.
65     + isq2(i).gt.mqt.or.isq3(i).gt.mqt) goto 113
66      qsa2(i)=chg(isq2(i),1,iset)
67      qsa3(i)=chg(isq3(i),1,iset)
68  112 continue
69  113 continue
70      else
71      do 1112 i=1,msa
72      if(isq2(i).le.0.or.isq3(i).le.0.or.
73     + isq2(i).gt.mqt.or.isq3(i).gt.mqt) goto 1113
74      qsa2(i)=chg(isq2(i),1,lseq(isgm(i)))
75      qsa3(i)=chg(isq3(i),1,lseq(isgm(i)))
76 1112 continue
77 1113 continue
78      endif
79c
80c     solute non-bonded interactions
81c     ==============================
82c
83      isfr=isfrom-1
84c
85c     loop over short and long range pairlists
86c
87      do 11 ipss=1,lpss
88c
89c     evaluate outer index array
90c
91      nsslen(ipss)=0
92      lssndx(0,ipss)=0
93      number=0
94      do 12 isa=1,nums
95      if(number+lssin(isa,ipss).gt.mscr.or.
96     + (ismf(isfr+isa).ne.ismf(isfr+isa-1).and.
97     + number.gt.0)) then
98      nsslen(ipss)=nsslen(ipss)+1
99      lssndx(nsslen(ipss),ipss)=isa-1
100      number=0
101      endif
102      number=number+lssin(isa,ipss)
103   12 continue
104      if(number.gt.0) then
105      nsslen(ipss)=nsslen(ipss)+1
106      lssndx(nsslen(ipss),ipss)=nums
107      endif
108c
109c     loop over number of cycles to complete pairlists
110c
111      do 13 isslen=1,nsslen(ipss)
112c
113      etermq=zero
114      eterml=zero
115c
116      nax=0
117      isf=ismf(isfr+lssndx(isslen,ipss))
118c
119c     collect coordinates into workarrays
120c
121      do 14 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss)
122      jsaptr=lssjpt(isa,ipss)-1
123      ism=isml(isfr+isa)
124      if(lpbc.or.lpbcs) then
125      if(ipbtyp.eq.1) then
126      do 15 jnum=1,lssin(isa,ipss)
127      lssptr=lssj(jsaptr+jnum)
128      rwc(nax+jnum,1)=xs(isfr+isa,1)-xs(lssptr,1)
129      rwc(nax+jnum,2)=xs(isfr+isa,2)-xs(lssptr,2)
130      rwc(nax+jnum,3)=xs(isfr+isa,3)-xs(lssptr,3)
131      isrx(nax+jnum)=0
132   15 continue
133      else
134      do 115 jnum=1,lssin(isa,ipss)
135      lssptr=lssj(jsaptr+jnum)
136      jsm=isml(lssptr)
137      rwc(nax+jnum,1)=xsm(ism,1)-xsm(jsm,1)
138      rwc(nax+jnum,2)=xsm(ism,2)-xsm(jsm,2)
139      rwc(nax+jnum,3)=xsm(ism,3)-xsm(jsm,3)
140      isrx(nax+jnum)=0
141  115 continue
142      endif
143      call argos_cafe_pbc(0,rwc,mscr,rwx,mscr,nax,1,lssin(isa,ipss))
144      endif
145      do 16 jnum=1,lssin(isa,ipss)
146      lssptr=lssj(jsaptr+jnum)
147      jsf=ismf(lssptr)
148      isal(nax+jnum)=isfr+isa
149      jsal(nax+jnum)=lssptr
150      jfal(nax+jnum)=jsf
151      jmal(nax+jnum)=0
152      jsm=isml(lssptr)
153      if(ism.ne.jsm) jmal(nax+jnum)=1
154      if(ism.gt.0) then
155      if(jsm.gt.0) then
156      rwc(nax+jnum,1)=xsm(ism,1)-xsm(jsm,1)
157      rwc(nax+jnum,2)=xsm(ism,2)-xsm(jsm,2)
158      rwc(nax+jnum,3)=xsm(ism,3)-xsm(jsm,3)
159      else
160      rwc(nax+jnum,1)=xsm(ism,1)-xs(lssptr,1)
161      rwc(nax+jnum,2)=xsm(ism,2)-xs(lssptr,2)
162      rwc(nax+jnum,3)=xsm(ism,3)-xs(lssptr,3)
163      endif
164      else
165      if(jsm.gt.0) then
166      rwc(nax+jnum,1)=xs(isfr+isa,1)-xsm(jsm,1)
167      rwc(nax+jnum,2)=xs(isfr+isa,2)-xsm(jsm,2)
168      rwc(nax+jnum,3)=xs(isfr+isa,3)-xsm(jsm,3)
169      else
170      rwc(nax+jnum,1)=xs(isfr+isa,1)-xs(lssptr,1)
171      rwc(nax+jnum,2)=xs(isfr+isa,2)-xs(lssptr,2)
172      rwc(nax+jnum,3)=xs(isfr+isa,3)-xs(lssptr,3)
173      endif
174      endif
175c
176      isrx(nax+jnum)=0
177c
178      if(lssscl) then
179c
180      istt=iand(isss(isfr+isa),48)
181      jstt=iand(isss(lssptr),48)
182      if(ism.ne.jsm) then
183      if(istt.eq.16.or.jstt.eq.16) isrx(nax+jnum)=-1
184      if(istt.eq.32.or.jstt.eq.32) isrx(nax+jnum)=1
185      endif
186c
187      istt=iand(isss(isfr+isa),384)
188      jstt=iand(isss(lssptr),384)
189      if(istt.eq.128.or.jstt.eq.128) isrx(nax+jnum)=-2
190      if(istt.eq.256.or.jstt.eq.256) isrx(nax+jnum)=2
191c
192      istt=iand(isss(isfr+isa),384)
193      jstt=iand(isss(lssptr),384)
194      if(istt.eq.128.and.jstt.eq.256) isrx(nax+jnum)=999
195      if(istt.eq.256.and.jstt.eq.128) isrx(nax+jnum)=999
196c
197c      write(*,'(5i5)')
198c     + isga(isfr+isa),isga(lssptr),istt,jstt,isrx(nax+jnum)
199c
200      endif
201c
202   16 continue
203c
204      do 17 jnum=1,lssin(isa,ipss)
205      lssptr=lssj(jsaptr+jnum)
206      facu(nax+jnum)=zero
207      if(iand(isdt(isfr+isa),mdynam).eq.ldynam.or.
208     + iand(isdt(lssptr),mdynam).eq.ldynam) facu(nax+jnum)=one
209c      if((iand(isdt(isfr+isa),mdynam).eq.ldynam.and.
210c     + iand(isdt(lssptr),mdynam).ne.ldynam) .or.
211c     + (iand(isdt(isfr+isa),mdynam).ne.ldynam.and.
212c     + iand(isdt(lssptr),mdynam).eq.ldynam)) facu(nax+jnum)=half
213      if(includ.eq.1) facu(nax+jnum)=one
214   17 continue
215c
216      if(.not.lpbc.and..not.lpbcs) then
217      do 18 jnum=1,lssin(isa,ipss)
218      lssptr=lssj(jsaptr+jnum)
219      xi(nax+jnum,1)=xs(isfr+isa,1)
220      xi(nax+jnum,2)=xs(isfr+isa,2)
221      xi(nax+jnum,3)=xs(isfr+isa,3)
222      xj(nax+jnum,1)=xs(lssptr,1)
223      xj(nax+jnum,2)=xs(lssptr,2)
224      xj(nax+jnum,3)=xs(lssptr,3)
225      isal(nax+jnum)=isfr+isa
226      jsal(nax+jnum)=lssptr
227   18 continue
228      else
229      do 19 jnum=1,lssin(isa,ipss)
230      rwc(nax+jnum,1)=rwc(nax+jnum,1)-rwx(jnum,1)
231      rwc(nax+jnum,2)=rwc(nax+jnum,2)-rwx(jnum,2)
232      rwc(nax+jnum,3)=rwc(nax+jnum,3)-rwx(jnum,3)
233      lssptr=lssj(jsaptr+jnum)
234      xi(nax+jnum,1)=xs(isfr+isa,1)
235      xi(nax+jnum,2)=xs(isfr+isa,2)
236      xi(nax+jnum,3)=xs(isfr+isa,3)
237      xj(nax+jnum,1)=xs(lssptr,1)+rwx(jnum,1)
238      xj(nax+jnum,2)=xs(lssptr,2)+rwx(jnum,2)
239      xj(nax+jnum,3)=xs(lssptr,3)+rwx(jnum,3)
240      isal(nax+jnum)=isfr+isa
241      jsal(nax+jnum)=lssptr
242   19 continue
243      endif
244c
245      nax=nax+lssin(isa,ipss)
246   14 continue
247c
248c
249c     evaluate electrostatic energies and forces
250c
251c      etermq=zero
252      if(.not.lssscl) then
253      if(ipme.eq.0.or.isolvo.ne.0) then
254      do 24 iax=1,nax
255      f(iax)=zero
256      rwx(iax,1)=xi(iax,1)-xj(iax,1)
257      rwx(iax,2)=xi(iax,2)-xj(iax,2)
258      rwx(iax,3)=xi(iax,3)-xj(iax,3)
259      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
260      rwi1(iax)=sqrt(rwi2(iax))
261      isa=isal(iax)
262      jsa=jsal(iax)
263c
264      if(jfal(iax).ne.isf) then
265      q=qsa2(isa)*qsa2(jsa)
266      else
267      q=qsa3(isa)*qsa3(jsa)
268      endif
269c
270      rw(iax)=facu(iax)*q*rwi1(iax)
271      f(iax)=f(iax)+q*rwi1(iax)*rwi2(iax)
272      if(ireact.ne.0) then
273      ess(isf,jfal(iax),6,ipss)=ess(isf,jfal(iax),6,ipss)+
274     + facu(iax)*q*rffss/rwi2(iax)
275      if(npener.ne.0) then
276      esa(isga(isa))=esa(isga(isa))+half*facu(iax)*q*rffss/rwi2(iax)
277      esa(isga(jsa))=esa(isga(jsa))+half*facu(iax)*q*rffss/rwi2(iax)
278      endif
279      f(iax)=f(iax)-two*q*rffss
280      endif
281cx      if(ihess.gt.0) then
282cx      h(iax)=three*q*rwi1(iax)*rwi2(iax)*rwi2(iax)
283cx      endif
284   24 continue
285      else
286      do 25 iax=1,nax
287      f(iax)=zero
288      rwx(iax,1)=xi(iax,1)-xj(iax,1)
289      rwx(iax,2)=xi(iax,2)-xj(iax,2)
290      rwx(iax,3)=xi(iax,3)-xj(iax,3)
291      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
292      if(isrx(iax).eq.999) then
293      rwi2(iax)=zero
294      isrx(iax)=0
295      endif
296      rwi1(iax)=sqrt(rwi2(iax))
297      isa=isal(iax)
298      jsa=jsal(iax)
299      if(jfal(iax).ne.isf) then
300      q=qsa2(isa)*qsa2(jsa)
301      else
302      q=qsa3(isa)*qsa3(jsa)
303      endif
304c
305      ferfc=erfc(ealpha/rwi1(iax))
306      fderfc=ealpha*derfc(ealpha/rwi1(iax))
307      rw(iax)=ferfc*facu(iax)*q*rwi1(iax)
308      f(iax)=f(iax)+q*rwi2(iax)*(ferfc*rwi1(iax)-fderfc)
309   25 continue
310      endif
311      else
312      if(ipme.eq.0.or.isolvo.ne.0) then
313      do 26 iax=1,nax
314      f(iax)=zero
315      rwx(iax,1)=xi(iax,1)-xj(iax,1)
316      rwx(iax,2)=xi(iax,2)-xj(iax,2)
317      rwx(iax,3)=xi(iax,3)-xj(iax,3)
318      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
319      if(isrx(iax).eq.999) then
320      rwi2(iax)=zero
321      isrx(iax)=0
322      endif
323      rtmp=rwi2(iax)
324      if(isrx(iax).gt.0) then
325      rwi2(iax)=one/(one/rwi2(iax)+shift0(1))
326      elseif(isrx(iax).lt.0) then
327      rwi2(iax)=one/(one/rwi2(iax)+shift1(1))
328      endif
329      rwi1(iax)=sqrt(rwi2(iax))
330      isa=isal(iax)
331      jsa=jsal(iax)
332c      write(*,'(3i5,4f12.6)')
333c     + isga(isa),isga(jsa),isrx(iax),shift0(1),shift1(1),rtmp,rwi2(iax)
334      if(jfal(iax).ne.isf) then
335      q=qsa2(isa)*qsa2(jsa)
336      else
337      q=qsa3(isa)*qsa3(jsa)
338      endif
339      rw(iax)=facu(iax)*q*rwi1(iax)
340      f(iax)=f(iax)+q*rwi1(iax)*rwi2(iax)
341      if(ireact.ne.0) then
342      rw(iax)=rw(iax)+facu(iax)*q*rffss/rwi2(iax)
343      f(iax)=f(iax)-two*q*rffss
344      endif
345   26 continue
346      else
347      do 126 iax=1,nax
348      f(iax)=zero
349      rwx(iax,1)=xi(iax,1)-xj(iax,1)
350      rwx(iax,2)=xi(iax,2)-xj(iax,2)
351      rwx(iax,3)=xi(iax,3)-xj(iax,3)
352      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
353      if(isrx(iax).eq.999) then
354      rwi2(iax)=zero
355      isrx(iax)=0
356      endif
357      if(isrx(iax).gt.0) then
358      rwi2(iax)=one/(one/rwi2(iax)+shift0(1))
359      elseif(isrx(iax).lt.0) then
360      rwi2(iax)=one/(one/rwi2(iax)+shift1(1))
361      endif
362      rwi1(iax)=sqrt(rwi2(iax))
363      isa=isal(iax)
364      jsa=jsal(iax)
365      if(jfal(iax).ne.isf) then
366      q=qsa2(isa)*qsa2(jsa)
367      else
368      q=qsa3(isa)*qsa3(jsa)
369      endif
370      ferfc=erfc(ealpha/rwi1(iax))
371      fderfc=ealpha*derfc(ealpha/rwi1(iax))
372      rw(iax)=ferfc*facu(iax)*q*rwi1(iax)
373      f(iax)=f(iax)+q*rwi2(iax)*(ferfc*rwi1(iax)-fderfc)
374      if(ireact.ne.0) then
375      rw(iax)=rw(iax)+facu(iax)*q*rffss/rwi2(iax)
376      f(iax)=f(iax)-two*q*rffss
377      endif
378  126 continue
379      endif
380      endif
381c
382c     accumulate electrostatic energies per solute molecule
383c
384c      etermq=zero
385      do 27 iax=1,nax
386      if(npener.ne.0) then
387      esa(isga(isal(iax)))=esa(isga(isal(iax)))+half*rw(iax)
388      esa(isga(jsal(iax)))=esa(isga(jsal(iax)))+half*rw(iax)
389      endif
390      ess(isf,jfal(iax),6,ipss)=ess(isf,jfal(iax),6,ipss)+rw(iax)
391      etermq=etermq+rw(iax)
392   27 continue
393c
394c      do 27 jsf=1,msf
395c      sumen=zero
396c      do 28 iax=1,nax
397c      if(jfal(iax).eq.jsf) sumen=sumen+rw(iax)
398c      if(npener.ne.0) then
399c      esa(isga(isal(iax)))=esa(isga(isal(iax)))+half*rw(iax)
400c      esa(isga(jsal(iax)))=esa(isga(jsal(iax)))+half*rw(iax)
401c      endif
402c   28 continue
403c      ess(isf,jsf,6,ipss)=ess(isf,jsf,6,ipss)+sumen
404c      etermq=etermq+sumen
405c   27 continue
406c
407c
408c
409c     Lennard-Jones energies and forces
410c     =================================
411c
412      do 29 iax=1,nax
413      isa=isal(iax)
414      jsa=jsal(iax)
415      rwi6(iax)=rwi2(iax)*rwi2(iax)*rwi2(iax)
416      c6=vdw(isat(isa),isat(jsa),1,iset)
417      c12=vdw(isat(isa),isat(jsa),3,iset)
418      cf6=six*c6
419      cf12=twelve*c12
420      rw(iax)=facu(iax)*(c12*rwi6(iax)-c6)*rwi6(iax)
421      f(iax)=f(iax)+(cf12*rwi6(iax)-cf6)*rwi6(iax)*rwi2(iax)
422cx      if(ihess.gt.0) then
423cx      h(iax)=h(iax)+(forten*cf12*rwi6(iax)-eight*cf6)*rwi6(iax)*
424cx     + rwi2(iax)*rwi2(iax)
425cx      endif
426   29 continue
427c
428c     accumulate Lennard-Jones energies per solute molecule
429c
430c      eterml=zero
431c      etermq=zero
432      do 30 iax=1,nax
433      if(npener.ne.0) then
434      esa(isga(isal(iax)))=esa(isga(isal(iax)))+half*rw(iax)
435      esa(isga(jsal(iax)))=esa(isga(jsal(iax)))+half*rw(iax)
436      endif
437      ess(isf,jfal(iax),5,ipss)=ess(isf,jfal(iax),5,ipss)+rw(iax)
438      eterml=eterml+rw(iax)
439   30 continue
440c
441c      do 30 jsf=1,msf
442c      sumen=zero
443c      do 31 iax=1,nax
444c      if(jfal(iax).eq.jsf) sumen=sumen+rw(iax)
445c      if(npener.ne.0) then
446c      esa(isga(isal(iax)))=esa(isga(isal(iax)))+half*rw(iax)
447c      esa(isga(jsal(iax)))=esa(isga(jsal(iax)))+half*rw(iax)
448c      endif
449c   31 continue
450c      ess(isf,jsf,5,ipss)=ess(isf,jsf,5,ipss)+sumen
451c      eterml=eterml+sumen
452c   30 continue
453c
454c     evaluate and accumulate the solute-solute virial contributions
455c     allow virial contributions from interactions between a solute
456c     molecule and its own image
457c
458      do 132 ix=1,3
459      do 32 jsf=1,msf
460      sumen1=zero
461      sumen2=zero
462      sumen3=zero
463      do 33 iax=1,nax
464cx      if(jfal(iax).eq.jsf.and.jmal(iax).eq.1) then
465      if(jfal(iax).eq.jsf) then
466      sumen1=sumen1-half*f(iax)*rwx(iax,1)*rwc(iax,ix)
467      sumen2=sumen2-half*f(iax)*rwx(iax,2)*rwc(iax,ix)
468      sumen3=sumen3-half*f(iax)*rwx(iax,3)*rwc(iax,ix)
469      endif
470   33 continue
471      zs(isf,ix,1,ipss)=zs(isf,ix,1,ipss)+sumen1
472      zs(jsf,ix,1,ipss)=zs(jsf,ix,1,ipss)+sumen1
473      zs(isf,ix,2,ipss)=zs(isf,ix,2,ipss)+sumen2
474      zs(jsf,ix,2,ipss)=zs(jsf,ix,2,ipss)+sumen2
475      zs(isf,ix,3,ipss)=zs(isf,ix,3,ipss)+sumen3
476      zs(jsf,ix,3,ipss)=zs(jsf,ix,3,ipss)+sumen3
477   32 continue
478  132 continue
479c
480c     generate radial distribution functions
481c
482c      if(ifstep-1.eq.((ifstep-1)/nfrdf)*nfrdf.and.ngrss.gt.0) then
483c      do 34 iax=1,nax
484c      isa=isal(iax)
485c      jsa=jsal(iax)
486c      do 35 igc=1,ngc
487c      if(ngt(igc).eq.3) then
488c      if((isga(isa).eq.iagc(igc).and.
489c     + isga(jsa).eq.jagc(igc)).or.
490c     + (isga(isa).eq.jagc(igc).and.
491c     + isga(jsa).eq.iagc(igc))) then
492c      igr=igrc(igc)
493c      indx=int(one/(rwi1(iax)*drdf))
494c      if(indx.gt.ngl) indx=ngl
495c      rdf(indx,igr)=rdf(indx,igr)+rdfvol
496c      endif
497c      endif
498c   35 continue
499c   34 continue
500c      endif
501c
502c     accumulate forces into solute force arrays
503c
504      nax=0
505      do 36 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss)
506      jsaptr=lssjpt(isa,ipss)-1
507      do 37 jnum=1,lssin(isa,ipss)
508      lssptr=lssj(jsaptr+jnum)
509      fs(isfr+isa,1,ipss)=fs(isfr+isa,1,ipss)+
510     + f(nax+jnum)*rwx(nax+jnum,1)
511      fs(isfr+isa,2,ipss)=fs(isfr+isa,2,ipss)+
512     + f(nax+jnum)*rwx(nax+jnum,2)
513      fs(isfr+isa,3,ipss)=fs(isfr+isa,3,ipss)+
514     + f(nax+jnum)*rwx(nax+jnum,3)
515      fs(lssptr,1,ipss)=fs(lssptr,1,ipss)-f(nax+jnum)*rwx(nax+jnum,1)
516      fs(lssptr,2,ipss)=fs(lssptr,2,ipss)-f(nax+jnum)*rwx(nax+jnum,2)
517      fs(lssptr,3,ipss)=fs(lssptr,3,ipss)-f(nax+jnum)*rwx(nax+jnum,3)
518      isf=ismf(isfr+isa)
519      jsf=ismf(lssptr)
520      fss(isf,jsf,1,ipss)=fss(isf,jsf,1,ipss)+
521     + f(nax+jnum)*rwx(nax+jnum,1)
522      fss(isf,jsf,2,ipss)=fss(isf,jsf,2,ipss)+
523     + f(nax+jnum)*rwx(nax+jnum,2)
524      fss(isf,jsf,3,ipss)=fss(isf,jsf,3,ipss)+
525     + f(nax+jnum)*rwx(nax+jnum,3)
526   37 continue
527cx      if(ihess.gt.0) then
528cx      do 137 jnum=1,lssin(isa,ipss)
529cx      lssptr=lssj(jsaptr+jnum)
530cx      hs(isfr+isa,1,ipss)=hs(isfr+isa,1,ipss)-f(nax+jnum)+
531cx     + h(nax+jnum)**rwx(nax+jnum,1)*rwx(nax+jnum,1)
532cx      hs(isfr+isa,2,ipss)=hs(isfr+isa,2,ipss)+
533cx     + h(nax+jnum)*rwx(nax+jnum,1)*rwx(nax+jnum,2)
534cx      hs(isfr+isa,3,ipss)=hs(isfr+isa,3,ipss)+
535cx     + h(nax+jnum)*rwx(nax+jnum,1)*rwx(nax+jnum,3)
536cx      hs(isfr+isa,4,ipss)=hs(isfr+isa,4,ipss)-f(nax+jnum)+
537cx     + h(nax+jnum)*rwx(nax+jnum,2)*rwx(nax+jnum,2)
538cx      hs(isfr+isa,5,ipss)=hs(isfr+isa,5,ipss)+
539cx     + h(nax+jnum)*rwx(nax+jnum,2)*rwx(nax+jnum,3)
540cx      hs(isfr+isa,6,ipss)=hs(isfr+isa,6,ipss)-f(nax+jnum)+
541cx     + h(nax+jnum)*rwx(nax+jnum,3)*rwx(nax+jnum,3)
542cx      hs(lssptr,1,ipss)=hs(lssptr,1,ipss)+f(nax+jnum)-
543cx     + h(nax+jnum)*rwx(nax+jnum,1)*rwx(nax+jnum,1)
544cx      hs(lssptr,2,ipss)=hs(lssptr,2,ipss)-
545cx     + h(nax+jnum)*rwx(nax+jnum,1)*rwx(nax+jnum,2)
546cx      hs(lssptr,3,ipss)=hs(lssptr,3,ipss)-
547cx     + h(nax+jnum)*rwx(nax+jnum,1)*rwx(nax+jnum,3)
548cx      hs(lssptr,4,ipss)=hs(lssptr,4,ipss)+f(nax+jnum)-
549cx     + h(nax+jnum)*rwx(nax+jnum,2)*rwx(nax+jnum,2)
550cx      hs(lssptr,5,ipss)=hs(lssptr,5,ipss)-
551cx     + h(nax+jnum)*rwx(nax+jnum,2)*rwx(nax+jnum,3)
552cx      hs(lssptr,6,ipss)=hs(lssptr,6,ipss)+f(nax+jnum)-
553cx     + h(nax+jnum)*rwx(nax+jnum,3)*rwx(nax+jnum,3)
554cx  137 continue
555cx      endif
556      nax=nax+lssin(isa,ipss)
557   36 continue
558c
559c     thermodynamic integration
560c
561      if(ithint) then
562      if(ith(14)) then
563      nax=0
564      do 38 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss)
565      jsaptr=lssjpt(isa,ipss)-1
566c
567      if(.not.lssscl) then
568      do 39 jnum=1,lssin(isa,ipss)
569      jsa=lssj(jsaptr+jnum)
570      dercon=(vdw(isat(isfr+isa),isat(jsa),3,4)*rwi6(nax+jnum)
571     + -vdw(isat(isfr+isa),isat(jsa),1,4))*rwi6(nax+jnum)
572      deriv(15,ipss)=deriv(15,ipss)+dercon
573      if(npgdec.gt.1) then
574      dera(3,isga(isa))=dera(3,isga(isa))+half*dercon
575      dera(3,isga(jsa))=dera(3,isga(jsa))+half*dercon
576      endif
577   39 continue
578      else
579      do 40 jnum=1,lssin(isa,ipss)
580      jsa=lssj(jsaptr+jnum)
581      dercon=(vdw(isat(isfr+isa),isat(jsa),3,4)*rwi6(nax+jnum)
582     + -vdw(isat(isfr+isa),isat(jsa),1,4))*rwi6(nax+jnum)
583      if(isrx(nax+jnum).gt.0) then
584      c64=three*vdw(isat(isfr+isa),isat(jsa),1,iset)
585      c124=six*vdw(isat(isfr+isa),isat(jsa),3,iset)
586      dercon=dercon+shift0(4)*
587     + rwi2(nax+jnum)*rwi6(nax+jnum)*(c64-c124*rwi6(nax+jnum))
588      elseif(isrx(nax+jnum).lt.0) then
589      c64=three*vdw(isat(isfr+isa),isat(jsa),1,iset)
590      c124=six*vdw(isat(isfr+isa),isat(jsa),3,iset)
591      dercon=dercon+shift1(4)*
592     + rwi2(nax+jnum)*rwi6(nax+jnum)*(c64-c124*rwi6(nax+jnum))
593      endif
594      deriv(15,ipss)=deriv(15,ipss)+dercon
595c      write(*,'(a,3i5,4f12.6)') 'gv ',
596c     + isga(isfr+isa),isga(jsa),isrx(nax+jnum),shift0(4),shift1(4),
597c     + dercon,deriv(15,ipss)
598      if(npgdec.gt.1) then
599      dera(3,isga(isfr+isa))=dera(3,isga(isfr+isa))+half*dercon
600      dera(3,isga(jsa))=dera(3,isga(jsa))+half*dercon
601      endif
602   40 continue
603      endif
604c
605      nax=nax+lssin(isa,ipss)
606   38 continue
607      endif
608c
609      if(ith(16)) then
610      nax=0
611      do 41 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss)
612      jsaptr=lssjpt(isa,ipss)-1
613      ism=isml(isfr+isa)
614      if(ipme.eq.0) then
615      if(.not.lssscl) then
616      do 42 jnum=1,lssin(isa,ipss)
617      jsa=lssj(jsaptr+jnum)
618      if(isml(jsa).ne.ism) then
619      dercon=(chg(isq2(isfr+isa),1,iset)*chg(isq2(jsa),1,4)
620     + +chg(isq2(isfr+isa),1,4)*chg(isq2(jsa),1,iset))
621      else
622      dercon=(chg(isq3(isfr+isa),1,iset)*chg(isq3(jsa),1,4)
623     + +chg(isq3(isfr+isa),1,4)*chg(isq3(jsa),1,iset))
624      endif
625      deriv(17,ipss)=deriv(17,ipss)+dercon*rwi1(nax+jnum)
626      if(ireact.ne.0) then
627      deriv(17,ipss)=deriv(17,ipss)+dercon*rffss/rwi2(nax+jnum)
628      endif
629c      write(*,'(a,3i5,4f12.6)') 'gq ',
630c     + isga(isfr+isa),isga(jsa),isrx(nax+jnum),shift0(4),shift1(4),
631c     + dercon,deriv(17,ipss)
632      if(npgdec.gt.1) then
633      dera(4,isga(isfr+isa))=dera(4,isga(isfr+isa))+
634     + half*dercon*rwi1(nax+jnum)
635      dera(4,isga(jsa))=dera(4,isga(jsa))+half*dercon*rwi1(nax+jnum)
636      endif
637   42 continue
638      else
639      do 43 jnum=1,lssin(isa,ipss)
640      jsa=lssj(jsaptr+jnum)
641      if(isml(jsa).ne.ism) then
642      dercon=(chg(isq2(isfr+isa),1,iset)*chg(isq2(jsa),1,4)
643     + +chg(isq2(isfr+isa),1,4)*chg(isq2(jsa),1,iset))
644      if(isrx(nax+jnum).gt.0) then
645      dercon=dercon-half*shift0(4)*
646     + chg(isq2(isfr+isa),1,iset)*chg(isq2(jsa),1,iset)*rwi2(nax+jnum)
647      elseif(isrx(nax+jnum).lt.0) then
648      dercon=dercon-half*shift1(4)*
649     + chg(isq2(isfr+isa),1,iset)*chg(isq2(jsa),1,iset)*rwi2(nax+jnum)
650      endif
651      else
652      dercon=(chg(isq3(isfr+isa),1,iset)*chg(isq3(jsa),1,4)
653     + +chg(isq3(isfr+isa),1,4)*chg(isq3(jsa),1,iset))
654      if(isrx(nax+jnum).gt.1) then
655      dercon=dercon-half*shift0(4)*
656     + chg(isq3(isfr+isa),1,iset)*chg(isq3(jsa),1,iset)*rwi2(nax+jnum)
657      elseif(isrx(nax+jnum).lt.-1) then
658      dercon=dercon-half*shift1(4)*
659     + chg(isq3(isfr+isa),1,iset)*chg(isq3(jsa),1,iset)*rwi2(nax+jnum)
660      endif
661      endif
662      deriv(17,ipss)=deriv(17,ipss)+dercon*rwi1(nax+jnum)
663      if(ireact.ne.0) then
664      deriv(17,ipss)=deriv(17,ipss)+dercon*rffss/rwi2(nax+jnum)
665      endif
666      if(npgdec.gt.1) then
667      dera(4,isga(isfr+isa))=dera(4,isga(isfr+isa))+
668     + half*dercon*rwi1(nax+jnum)
669      dera(4,isga(jsa))=dera(4,isga(jsa))+half*dercon*rwi1(nax+jnum)
670      endif
671   43 continue
672      endif
673      else
674      if(.not.lssscl) then
675      do 142 jnum=1,lssin(isa,ipss)
676      jsa=lssj(jsaptr+jnum)
677      if(isml(jsa).ne.ism) then
678      dercon=(chg(isq2(isfr+isa),1,iset)*chg(isq2(jsa),1,4)
679     + +chg(isq2(isfr+isa),1,4)*chg(isq2(jsa),1,iset))
680      else
681      dercon=(chg(isq3(isfr+isa),1,iset)*chg(isq3(jsa),1,4)
682     + +chg(isq3(isfr+isa),1,4)*chg(isq3(jsa),1,iset))
683      endif
684      deriv(17,ipss)=deriv(17,ipss)+dercon*rwi1(nax+jnum)
685      if(ireact.ne.0) then
686      deriv(17,ipss)=deriv(17,ipss)+dercon*rffss/rwi2(nax+jnum)
687      endif
688      if(npgdec.gt.1) then
689      dera(4,isga(isfr+isa))=dera(4,isga(isfr+isa))+
690     + half*dercon*rwi1(nax+jnum)
691      dera(4,isga(jsa))=dera(4,isga(jsa))+half*dercon*rwi1(nax+jnum)
692      endif
693  142 continue
694      else
695      do 143 jnum=1,lssin(isa,ipss)
696      jsa=lssj(jsaptr+jnum)
697      if(isml(jsa).ne.ism) then
698      dercon=(chg(isq2(isfr+isa),1,iset)*chg(isq2(jsa),1,4)
699     + +chg(isq2(isfr+isa),1,4)*chg(isq2(jsa),1,iset))
700      if(isrx(nax+jnum).gt.0) then
701      dercon=dercon-half*shift0(4)*
702     + chg(isq2(isfr+isa),1,iset)*chg(isq2(jsa),1,iset)*rwi2(nax+jnum)
703      elseif(isrx(nax+jnum).lt.0) then
704      dercon=dercon-half*shift1(4)*
705     + chg(isq2(isfr+isa),1,iset)*chg(isq2(jsa),1,iset)*rwi2(nax+jnum)
706      endif
707      else
708      dercon=(chg(isq3(isfr+isa),1,iset)*chg(isq3(jsa),1,4)
709     + +chg(isq3(isfr+isa),1,4)*chg(isq3(jsa),1,iset))
710      if(isrx(nax+jnum).gt.1) dercon=dercon-half*shift0(4)*
711     + chg(isq3(isfr+isa),1,iset)*chg(isq3(jsa),1,iset)*rwi2(nax+jnum)
712      if(isrx(nax+jnum).lt.-1) dercon=dercon-half*shift1(4)*
713     + chg(isq3(isfr+isa),1,iset)*chg(isq3(jsa),1,iset)*rwi2(nax+jnum)
714      endif
715      deriv(17,ipss)=deriv(17,ipss)+dercon*rwi1(nax+jnum)
716      if(ireact.ne.0) then
717      deriv(17,ipss)=deriv(17,ipss)+dercon*rffss/rwi2(nax+jnum)
718      endif
719      if(npgdec.gt.1) then
720      dera(4,isga(isfr+isa))=dera(4,isga(isfr+isa))+
721     + half*dercon*rwi1(nax+jnum)
722      dera(4,isga(jsa))=dera(4,isga(jsa))+half*dercon*rwi1(nax+jnum)
723      endif
724  143 continue
725      endif
726      endif
727      nax=nax+lssin(isa,ipss)
728   41 continue
729      endif
730      endif
731c
732c     thermodynamic perturbation 1
733c
734      if(ipert2) then
735      if(ip2(14)) then
736      if(.not.lssscl) then
737      do 44 iax=1,nax
738      isa=isal(iax)
739      jsa=jsal(iax)
740      ep2(ipss)=ep2(ipss)
741     + +facu(iax)*(vdw(isat(isa),isat(jsa),3,2)*rwi6(iax)
742     + -vdw(isat(isa),isat(jsa),1,2))*rwi6(iax)
743   44 continue
744      else
745      do 45 iax=1,nax
746      isa=isal(iax)
747      jsa=jsal(iax)
748      rwi6(iax)=rwi2(iax)**3
749      if(isrx(iax).gt.0) then
750      rwi6(iax)=(one/(one/rwi2(iax)-shift0(1)+shift0(2)))**3
751      elseif(isrx(iax).lt.0) then
752      rwi6(iax)=(one/(one/rwi2(iax)-shift1(1)+shift1(2)))**3
753      endif
754      ep2(ipss)=ep2(ipss)
755     + +facu(iax)*(vdw(isat(isa),isat(jsa),3,2)*rwi6(iax)
756     + -vdw(isat(isa),isat(jsa),1,2))*rwi6(iax)
757   45 continue
758      endif
759      ep2(ipss)=ep2(ipss)-eterml
760      endif
761      if(ip2(16).or.ip2(17)) then
762      if(ipme.eq.0) then
763      if(.not.lssscl) then
764      do 46 iax=1,nax
765      isa=isal(iax)
766      jsa=jsal(iax)
767      if(jmal(iax).ne.0) then
768      q14=chg(isq2(isa),1,2)*chg(isq2(jsa),1,2)
769      else
770      q14=chg(isq3(isa),1,2)*chg(isq3(jsa),1,2)
771      endif
772      rwx(iax,1)=xi(iax,1)-xj(iax,1)
773      rwx(iax,2)=xi(iax,2)-xj(iax,2)
774      rwx(iax,3)=xi(iax,3)-xj(iax,3)
775      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
776      rwi1(iax)=sqrt(rwi2(iax))
777      ep2(ipss)=ep2(ipss)+facu(iax)*q14*rwi1(iax)
778      if(ireact.ne.0) then
779      ep2(ipss)=ep2(ipss)+facu(iax)*q14*rffss/rwi2(iax)
780      endif
781   46 continue
782      else
783      do 47 iax=1,nax
784      isa=isal(iax)
785      jsa=jsal(iax)
786      if(jmal(iax).ne.0) then
787      q14=chg(isq2(isa),1,2)*chg(isq2(jsa),1,2)
788      istt=0
789      else
790      q14=chg(isq3(isa),1,2)*chg(isq3(jsa),1,2)
791      istt=1
792      endif
793      rwx(iax,1)=xi(iax,1)-xj(iax,1)
794      rwx(iax,2)=xi(iax,2)-xj(iax,2)
795      rwx(iax,3)=xi(iax,3)-xj(iax,3)
796      rwi6(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
797      if(isrx(iax).gt.istt) then
798      rwi6(iax)=one/(one/rwi6(iax)+shift0(2))
799      elseif(isrx(iax).lt.-istt) then
800      rwi6(iax)=one/(one/rwi6(iax)+shift1(2))
801      endif
802      rwi1(iax)=sqrt(rwi6(iax))
803      ep2(ipss)=ep2(ipss)+facu(iax)*q14*rwi1(iax)
804      if(ireact.ne.0) then
805      ep2(ipss)=ep2(ipss)+facu(iax)*q14*rffss/rwi2(iax)
806      endif
807   47 continue
808      endif
809      else
810      if(.not.lssscl) then
811      do 146 iax=1,nax
812      isa=isal(iax)
813      jsa=jsal(iax)
814      if(jmal(iax).ne.0) then
815      q14=chg(isq2(isa),1,2)*chg(isq2(jsa),1,2)*
816     + erfc(ealpha/rwi1(iax))
817      else
818      q14=chg(isq3(isa),1,2)*chg(isq3(jsa),1,2)*
819     + erfc(ealpha/rwi1(iax))
820      endif
821      rwx(iax,1)=xi(iax,1)-xj(iax,1)
822      rwx(iax,2)=xi(iax,2)-xj(iax,2)
823      rwx(iax,3)=xi(iax,3)-xj(iax,3)
824      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
825      rwi1(iax)=sqrt(rwi2(iax))
826      ep2(ipss)=ep2(ipss)+facu(iax)*q14*rwi1(iax)
827      if(ireact.ne.0) then
828      ep2(ipss)=ep2(ipss)+facu(iax)*q14*rffss/rwi2(iax)
829      endif
830  146 continue
831      else
832      do 147 iax=1,nax
833      isa=isal(iax)
834      jsa=jsal(iax)
835      if(jmal(iax).ne.0) then
836      q14=chg(isq2(isa),1,2)*chg(isq2(jsa),1,2)*
837     + erfc(ealpha/rwi1(iax))
838      istt=0
839      else
840      q14=chg(isq3(isa),1,2)*chg(isq3(jsa),1,2)*
841     + erfc(ealpha/rwi1(iax))
842      istt=1
843      endif
844      rwx(iax,1)=xi(iax,1)-xj(iax,1)
845      rwx(iax,2)=xi(iax,2)-xj(iax,2)
846      rwx(iax,3)=xi(iax,3)-xj(iax,3)
847      rwi6(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
848      if(isrx(iax).gt.istt) then
849      rwi6(iax)=one/(one/rwi6(iax)+shift0(2))
850      elseif(isrx(iax).lt.-istt) then
851      rwi6(iax)=one/(one/rwi6(iax)+shift1(2))
852      endif
853      rwi1(iax)=sqrt(rwi6(iax))
854      ep2(ipss)=ep2(ipss)+facu(iax)*q14*rwi1(iax)
855      if(ireact.ne.0) then
856      ep2(ipss)=ep2(ipss)+facu(iax)*q14*rffss/rwi2(iax)
857      endif
858  147 continue
859      endif
860      endif
861      ep2(ipss)=ep2(ipss)-etermq
862      endif
863      endif
864c
865c     thermodynamic perturbation 2
866c
867      if(ipert3) then
868      if(ip3(14)) then
869      if(.not.lssscl) then
870      do 48 iax=1,nax
871      isa=isal(iax)
872      jsa=jsal(iax)
873      ep3(ipss)=ep3(ipss)
874     + +facu(iax)*(vdw(isat(isa),isat(jsa),3,3)*rwi6(iax)
875     + -vdw(isat(isa),isat(jsa),1,3))*rwi6(iax)
876   48 continue
877      else
878      do 49 iax=1,nax
879      isa=isal(iax)
880      jsa=jsal(iax)
881      rwi6(iax)=rwi2(iax)**3
882      if(isrx(iax).gt.0) then
883      rwi6(iax)=(one/(one/rwi2(iax)-shift0(1)+shift0(3)))**3
884      elseif(isrx(iax).lt.0) then
885      rwi6(iax)=(one/(one/rwi2(iax)-shift1(1)+shift1(3)))**3
886      endif
887      ep3(ipss)=ep3(ipss)
888     + +facu(iax)*(vdw(isat(isa),isat(jsa),3,3)*rwi6(iax)
889     + -vdw(isat(isa),isat(jsa),1,3))*rwi6(iax)
890   49 continue
891      endif
892      ep3(ipss)=ep3(ipss)-eterml
893      endif
894      if(ip2(16).or.ip2(17)) then
895      if(ipme.eq.0) then
896      if(.not.lssscl) then
897      do 50 iax=1,nax
898      isa=isal(iax)
899      jsa=jsal(iax)
900      if(jmal(iax).ne.0) then
901      q14=chg(isq2(isa),1,3)*chg(isq2(jsa),1,3)
902      else
903      q14=chg(isq3(isa),1,3)*chg(isq3(jsa),1,3)
904      endif
905      rwx(iax,1)=xi(iax,1)-xj(iax,1)
906      rwx(iax,2)=xi(iax,2)-xj(iax,2)
907      rwx(iax,3)=xi(iax,3)-xj(iax,3)
908      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
909      rwi1(iax)=sqrt(rwi2(iax))
910      ep3(ipss)=ep3(ipss)+facu(iax)*q14*rwi1(iax)
911      if(ireact.ne.0) then
912      ep3(ipss)=ep3(ipss)+facu(iax)*q14*rffss/rwi2(iax)
913      endif
914   50 continue
915      else
916      do 51 iax=1,nax
917      isa=isal(iax)
918      jsa=jsal(iax)
919      if(jmal(iax).ne.0) then
920      q14=chg(isq2(isa),1,3)*chg(isq2(jsa),1,3)
921      istt=0
922      else
923      q14=chg(isq3(isa),1,3)*chg(isq3(jsa),1,3)
924      istt=1
925      endif
926      rwx(iax,1)=xi(iax,1)-xj(iax,1)
927      rwx(iax,2)=xi(iax,2)-xj(iax,2)
928      rwx(iax,3)=xi(iax,3)-xj(iax,3)
929      rwi6(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
930      if(isrx(iax).gt.istt) then
931      rwi6(iax)=one/(one/rwi6(iax)+shift0(3))
932      elseif(isrx(iax).lt.-istt) then
933      rwi6(iax)=one/(one/rwi6(iax)+shift1(3))
934      endif
935      rwi1(iax)=sqrt(rwi6(iax))
936      ep3(ipss)=ep3(ipss)+facu(iax)*q14*rwi1(iax)
937      if(ireact.ne.0) then
938      ep3(ipss)=ep3(ipss)+facu(iax)*q14*rffss/rwi2(iax)
939      endif
940   51 continue
941      endif
942      else
943      if(.not.lssscl) then
944      do 150 iax=1,nax
945      isa=isal(iax)
946      jsa=jsal(iax)
947      if(jmal(iax).ne.0) then
948      q14=chg(isq2(isa),1,3)*chg(isq2(jsa),1,3)*
949     + erfc(ealpha/rwi1(iax))
950      else
951      q14=chg(isq3(isa),1,3)*chg(isq3(jsa),1,3)*
952     + erfc(ealpha/rwi1(iax))
953      endif
954      rwx(iax,1)=xi(iax,1)-xj(iax,1)
955      rwx(iax,2)=xi(iax,2)-xj(iax,2)
956      rwx(iax,3)=xi(iax,3)-xj(iax,3)
957      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
958      rwi1(iax)=sqrt(rwi2(iax))
959      ep3(ipss)=ep3(ipss)+facu(iax)*q14*rwi1(iax)
960      if(ireact.ne.0) then
961      ep3(ipss)=ep3(ipss)+facu(iax)*q14*rffss/rwi2(iax)
962      endif
963  150 continue
964      else
965      do 151 iax=1,nax
966      isa=isal(iax)
967      jsa=jsal(iax)
968      if(jmal(iax).ne.0) then
969      q14=chg(isq2(isa),1,3)*chg(isq2(jsa),1,3)*
970     + erfc(ealpha/rwi1(iax))
971      istt=0
972      else
973      q14=chg(isq3(isa),1,3)*chg(isq3(jsa),1,3)*
974     + erfc(ealpha/rwi1(iax))
975      istt=1
976      endif
977      rwx(iax,1)=xi(iax,1)-xj(iax,1)
978      rwx(iax,2)=xi(iax,2)-xj(iax,2)
979      rwx(iax,3)=xi(iax,3)-xj(iax,3)
980      rwi6(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
981      if(isrx(iax).gt.istt) then
982      rwi6(iax)=one/(one/rwi6(iax)+shift0(3))
983      elseif(isrx(iax).lt.-istt) then
984      rwi6(iax)=one/(one/rwi6(iax)+shift1(3))
985      endif
986      rwi1(iax)=sqrt(rwi6(iax))
987      ep3(ipss)=ep3(ipss)+facu(iax)*q14*rwi1(iax)
988      if(ireact.ne.0) then
989      ep3(ipss)=ep3(ipss)+facu(iax)*q14*rffss/rwi2(iax)
990      endif
991  151 continue
992      endif
993      endif
994      ep3(ipss)=ep3(ipss)-etermq
995      endif
996      endif
997   13 continue
998   11 continue
999c
1000c     accumulate radial distribution function contributions from
1001c     the excluded pairlist
1002c
1003c      if(ifstep-1.eq.((ifstep-1)/nfrdf)*nfrdf.and.ngrss.gt.0) then
1004c      do 52 isx=1,nsx
1005c      isa=idsx(isx)
1006c      jsa=jdsx(isx)
1007c      do 53 igc=1,ngc
1008c      if(ngt(igc).eq.3) then
1009c      if((isa.eq.iagc(igc).and.jsa.eq.jagc(igc)).or.
1010c     + (isa.eq.iagc(igc).and.jsa.eq.jagc(igc))) then
1011c      igr=igrc(igc)
1012c      indx=int(sqrt((xs(isa,1)-xs(jsa,1))**2+(xs(isa,2)-xs(jsa,2))**2+
1013c     + (xs(isa,3)-xs(jsa,3))**2)/drdf)
1014c      if(indx.gt.ngl) indx=ngl
1015c      rdf(indx,igr)=rdf(indx,igr)+rdfvol
1016c      endif
1017c      endif
1018c   53 continue
1019c   52 continue
1020c      endif
1021c
1022c
1023      return
1024      end
1025