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