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