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