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