1 subroutine argos_cafe_fpww(xw,xwm,fw,pw,pwp,idt,iwfrom, 2 + nwloc,lpbc,eww, 3 + vdw,chg,iwatm,iwq,lwwndx,lwwjpt,lwwin,lwwj, 4 + xi,xj,rwx,rwi1,rwi2,rwi6,rwc, 5 + f,fi,fj,facu,pl,pj) 6c 7c $Id$ 8c 9 implicit none 10c 11#include "argos_cafe_common.fh" 12#include "argos_cafe_funcs_dec.fh" 13#include "bitops_decls.fh" 14c 15 real*8 xw(mwm,3,mwa),xwm(mwm,3),fw(mwm,3,mwa,2),eww(mpe,2) 16 integer idt(mwm) 17 integer iwfrom,nwloc 18 logical lpbc 19c 20 real*8 vdw(mat,mat,map,mset),chg(mqt,mqp,mset) 21 integer iwatm(mwa),iwq(mwa) 22c 23 real*8 xi(mscr,3,mwa),xj(mscr,3,mwa),rwx(mscr,3) 24 real*8 rwi1(mscr),rwi2(mscr),rwi6(mscr),rwc(mscr,3) 25 real*8 f(mscr),fi(mscr,3,mwa),fj(mscr,3,mwa) 26c 27 real*8 facu(mscr) 28c real*8 rdf(mgl,mgr) 29c 30 integer lwwj(*) 31 integer lwwndx(0:mwm,2),lwwjpt(nwloc,2),lwwin(nwloc,2) 32c 33 34 real*8 pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2) 35 real*8 pl(mscr,3,mwa),pj(mscr,3,mwa) 36c integer nax2,ipset 37 real*8 qai,qaj,pai,paj,pix,piy,piz,pjx,pjy,pjz 38 real*8 ri3,rmi,rmj,fri,fmi,fmj,rmm,qfaci 39 real*8 rx,ry,rz,ri1,ri2,ewwpsm,etermp 40 real*8 ewwqsm 41c 42 integer iwfr,ipww,number,iwm,iwpm,nax 43 integer iwmn,lwwptr,iwa,iax,jwa,iptr,jptr,iwpj 44 real*8 ewwl6,ewwl12,q 45 real*8 c64,c124,qi,qj,qi4,qj4,dercon 46 real*8 c6p,c12p,qp,ep2tmp,ep3tmp 47 real*8 c6,cf6,c12,cf12 48c 49 integer nwwlen(2) 50 real*8 eterml,etermq 51c 52#include "argos_cafe_funcs_sfn.fh" 53#include "bitops_funcs.fh" 54c 55cx new stuff begin 56c 57cx new stuff end 58c 59c calculation of solvent-solvent intermolecular energies and forces 60c 61c subtract 1 from first molecule index for use as offset 62c 63 iwfr=iwfrom-1 64c 65c loop over short and long range parts 66c 67 do 1 ipww=1,lpww 68c 69c Evaluate the outer index array 70c 71 nwwlen(ipww)=0 72 lwwndx(0,ipww)=0 73 number=0 74 do 2 iwm=1,nwloc 75 if(number+lwwin(iwm,ipww).gt.mscr) then 76 nwwlen(ipww)=nwwlen(ipww)+1 77 lwwndx(nwwlen(ipww),ipww)=iwm-1 78 number=0 79 endif 80 number=number+lwwin(iwm,ipww) 81 2 continue 82 if(number.gt.0) then 83 nwwlen(ipww)=nwwlen(ipww)+1 84 lwwndx(nwwlen(ipww),ipww)=nwloc 85 endif 86c 87c loop over number of cycles to complete pairlist 88c 89 do 3 iwpm=1,nwwlen(ipww) 90 nax=0 91c 92c collect coordinates into workarrays 93c 94 do 4 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww) 95 iwpj=lwwjpt(iwm,ipww)-1 96 do 5 iwmn=1,lwwin(iwm,ipww) 97 lwwptr=lwwj(iwpj+iwmn) 98 rwc(nax+iwmn,1)=xwm(iwfr+iwm,1)-xwm(lwwptr,1) 99 rwc(nax+iwmn,2)=xwm(iwfr+iwm,2)-xwm(lwwptr,2) 100 rwc(nax+iwmn,3)=xwm(iwfr+iwm,3)-xwm(lwwptr,3) 101 facu(nax+iwmn)=one 102c if( (iand(idt(iwm),mdynam).eq.ldynam.and. 103c + iand(idt(lwwptr),mdynam).ne.ldynam).or. 104c + (iand(idt(iwm),mdynam).ne.ldynam.and. 105c + iand(idt(lwwptr),mdynam).eq.ldynam) ) facu(nax+iwmn)=half 106 if(iand(idt(iwm),mdynam).ne.ldynam.and. 107 + iand(idt(lwwptr),mdynam).ne.ldynam) facu(nax+iwmn)=zero 108 if(includ.eq.1) facu(nax+iwmn)=one 109 5 continue 110c 111 do 6 iwa=1,mwa 112 do 7 iwmn=1,lwwin(iwm,ipww) 113 lwwptr=lwwj(iwpj+iwmn) 114 xi(nax+iwmn,1,iwa)=xw(iwfr+iwm,1,iwa) 115 xi(nax+iwmn,2,iwa)=xw(iwfr+iwm,2,iwa) 116 xi(nax+iwmn,3,iwa)=xw(iwfr+iwm,3,iwa) 117 xj(nax+iwmn,1,iwa)=xw(lwwptr,1,iwa) 118 xj(nax+iwmn,2,iwa)=xw(lwwptr,2,iwa) 119 xj(nax+iwmn,3,iwa)=xw(lwwptr,3,iwa) 120 pl(nax+iwmn,1,iwa)=pw(iwfr+iwm,1,iwa,1) 121 pl(nax+iwmn,2,iwa)=pw(iwfr+iwm,2,iwa,1) 122 pl(nax+iwmn,3,iwa)=pw(iwfr+iwm,3,iwa,1) 123 pj(nax+iwmn,1,iwa)=pw(lwwptr,1,iwa,1) 124 pj(nax+iwmn,2,iwa)=pw(lwwptr,2,iwa,1) 125 pj(nax+iwmn,3,iwa)=pw(lwwptr,3,iwa,1) 126 7 continue 127 6 continue 128 if(lpbc) then 129 call argos_cafe_pbc(0,rwc,mscr,rwx,mscr,nax,1,lwwin(iwm,ipww)) 130 do 8 iwmn=1,lwwin(iwm,ipww) 131 rwc(nax+iwmn,1)=rwc(nax+iwmn,1)-rwx(iwmn,1) 132 rwc(nax+iwmn,2)=rwc(nax+iwmn,2)-rwx(iwmn,2) 133 rwc(nax+iwmn,3)=rwc(nax+iwmn,3)-rwx(iwmn,3) 134 8 continue 135 do 9 iwa=1,mwa 136 do 10 iwmn=1,lwwin(iwm,ipww) 137 lwwptr=lwwj(iwpj+iwmn) 138 xj(nax+iwmn,1,iwa)=xj(nax+iwmn,1,iwa)+rwx(iwmn,1) 139 xj(nax+iwmn,2,iwa)=xj(nax+iwmn,2,iwa)+rwx(iwmn,2) 140 xj(nax+iwmn,3,iwa)=xj(nax+iwmn,3,iwa)+rwx(iwmn,3) 141 10 continue 142 9 continue 143 endif 144c 145 nax=nax+lwwin(iwm,ipww) 146 4 continue 147c 148c initializations 149c 150c if(npener.ne.0) then 151c do 12 iax=1,nax 152c u(iax)=zero 153c 12 continue 154c endif 155c 156c loops over number of atoms in a solvent molecule 157c 158 qfaci=one/qfac 159 do 13 iwa=1,mwa 160 qi=chg(iwq(iwa),1,iset) 161 pai=chg(iwq(iwa),2,iset) 162 qai=qfaci*qi 163 do 14 jwa=1,mwa 164 qj=chg(iwq(jwa),1,iset) 165 q=qi*qj 166 paj=chg(iwq(jwa),2,iset) 167 qaj=qfaci*qj 168c 169 do 15 iax=1,nax 170 f(iax)=zero 171 rwx(iax,1)=xi(iax,1,iwa)-xj(iax,1,jwa) 172 rwx(iax,2)=xi(iax,2,iwa)-xj(iax,2,jwa) 173 rwx(iax,3)=xi(iax,3,iwa)-xj(iax,3,jwa) 174 rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2) 175 rwi1(iax)=sqrt(rwi2(iax)) 176 15 continue 177c 178c 179c van der Waals contribution 180c -------------------------- 181c 182 iptr=iwatm(iwa) 183 jptr=iwatm(jwa) 184 c6=vdw(iptr,jptr,1,iset) 185 cf6=six*c6 186 c12=vdw(iptr,jptr,3,iset) 187 cf12=twelve*c12 188c 189 eterml=zero 190 if(c6.ne.zero.or.c12.ne.zero) then 191 ewwl6=zero 192 ewwl12=zero 193 do 20 iax=1,nax 194 rwi6(iax)=rwi2(iax)*rwi2(iax)*rwi2(iax) 195 ewwl6=ewwl6+facu(iax)*rwi6(iax) 196 ewwl12=ewwl12+facu(iax)*rwi6(iax)*rwi6(iax) 197 f(iax)=f(iax)+(cf12*rwi6(iax)-cf6)*rwi6(iax)*rwi2(iax) 198 20 continue 199 eterml=c12*ewwl12-c6*ewwl6 200 eww(7,ipww)=eww(7,ipww)+eterml 201 endif 202c 203c 204c force vectors 205c ------------- 206c 207 if(iwa.eq.1) then 208 do 22 iax=1,nax 209 fj(iax,1,jwa)=(-f(iax))*rwx(iax,1) 210 fj(iax,2,jwa)=(-f(iax))*rwx(iax,2) 211 fj(iax,3,jwa)=(-f(iax))*rwx(iax,3) 212 22 continue 213 else 214 do 23 iax=1,nax 215 fj(iax,1,jwa)=fj(iax,1,jwa)-f(iax)*rwx(iax,1) 216 fj(iax,2,jwa)=fj(iax,2,jwa)-f(iax)*rwx(iax,2) 217 fj(iax,3,jwa)=fj(iax,3,jwa)-f(iax)*rwx(iax,3) 218 23 continue 219 endif 220c 221 if(jwa.eq.1) then 222 do 24 iax=1,nax 223 fi(iax,1,iwa)=f(iax)*rwx(iax,1) 224 fi(iax,2,iwa)=f(iax)*rwx(iax,2) 225 fi(iax,3,iwa)=f(iax)*rwx(iax,3) 226 24 continue 227 else 228 do 25 iax=1,nax 229 fi(iax,1,iwa)=fi(iax,1,iwa)+f(iax)*rwx(iax,1) 230 fi(iax,2,iwa)=fi(iax,2,iwa)+f(iax)*rwx(iax,2) 231 fi(iax,3,iwa)=fi(iax,3,iwa)+f(iax)*rwx(iax,3) 232 25 continue 233 endif 234 do 26 iax=1,nax 235 zw(1,1,ipww)=zw(1,1,ipww)-f(iax)*rwx(iax,1)*rwc(iax,1) 236 zw(2,1,ipww)=zw(2,1,ipww)-f(iax)*rwx(iax,1)*rwc(iax,2) 237 zw(3,1,ipww)=zw(3,1,ipww)-f(iax)*rwx(iax,1)*rwc(iax,3) 238 zw(1,2,ipww)=zw(1,2,ipww)-f(iax)*rwx(iax,2)*rwc(iax,1) 239 zw(2,2,ipww)=zw(2,2,ipww)-f(iax)*rwx(iax,2)*rwc(iax,2) 240 zw(3,2,ipww)=zw(3,2,ipww)-f(iax)*rwx(iax,2)*rwc(iax,3) 241 zw(1,3,ipww)=zw(1,3,ipww)-f(iax)*rwx(iax,3)*rwc(iax,1) 242 zw(2,3,ipww)=zw(2,3,ipww)-f(iax)*rwx(iax,3)*rwc(iax,2) 243 zw(3,3,ipww)=zw(3,3,ipww)-f(iax)*rwx(iax,3)*rwc(iax,3) 244 26 continue 245c 246c 247c electrostatic and polarization contribution 248c ------------------------------------------- 249c 250 ewwqsm=zero 251 ewwpsm=zero 252 do 117 iax=1,nax 253 pix=pai*pl(iax,1,iwa) 254 piy=pai*pl(iax,2,iwa) 255 piz=pai*pl(iax,3,iwa) 256 pjx=paj*pj(iax,1,jwa) 257 pjy=paj*pj(iax,2,jwa) 258 pjz=paj*pj(iax,3,jwa) 259 rx=-rwx(iax,1) 260 ry=-rwx(iax,2) 261 rz=-rwx(iax,3) 262 ri1=rwi1(iax) 263 ri2=rwi2(iax) 264 ri3=qfac*qfac*ri1*ri2 265 rmi=three*(rx*pix+ry*piy+rz*piz)*ri2 266 rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2 267 if(ipolt.eq.1) then 268 fri=((-qai)*qaj+qai*rmj-qaj*rmi)*ri3 269 fmi=(qaj)*ri3 270 fmj=(-qai)*ri3 271 else 272 rmm=three*(pix*pjx+piy*pjy+piz*pjz)*ri2 273 fri=((-qai)*qaj+qai*rmj-qaj*rmi+5.0*rmi*rmj/three-rmm)*ri3 274 fmi=(qaj-rmj)*ri3 275 fmj=((-qai)-rmi)*ri3 276 endif 277 fi(iax,1,iwa)=fi(iax,1,iwa)+fri*rx+fmi*pix+fmj*pjx 278 fi(iax,2,iwa)=fi(iax,2,iwa)+fri*ry+fmi*piy+fmj*pjy 279 fi(iax,3,iwa)=fi(iax,3,iwa)+fri*rz+fmi*piz+fmj*pjz 280 fj(iax,1,jwa)=fj(iax,1,jwa)-(fri*rx+fmi*pix+fmj*pjx) 281 fj(iax,2,jwa)=fj(iax,2,jwa)-(fri*ry+fmi*piy+fmj*pjy) 282 fj(iax,3,jwa)=fj(iax,3,jwa)-(fri*rz+fmi*piz+fmj*pjz) 283 zw(1,1,ipww)=zw(1,1,ipww)-(fri*rx+fmi*pix+fmj*pjx)*rwc(iax,1) 284 zw(2,1,ipww)=zw(2,1,ipww)-(fri*rx+fmi*pix+fmj*pjx)*rwc(iax,2) 285 zw(3,1,ipww)=zw(3,1,ipww)-(fri*rx+fmi*pix+fmj*pjx)*rwc(iax,3) 286 zw(1,2,ipww)=zw(1,2,ipww)-(fri*ry+fmi*piy+fmj*pjy)*rwc(iax,1) 287 zw(2,2,ipww)=zw(2,2,ipww)-(fri*ry+fmi*piy+fmj*pjy)*rwc(iax,2) 288 zw(3,2,ipww)=zw(3,2,ipww)-(fri*ry+fmi*piy+fmj*pjy)*rwc(iax,3) 289 zw(1,3,ipww)=zw(1,3,ipww)-(fri*rz+fmi*piz+fmj*pjz)*rwc(iax,1) 290 zw(2,3,ipww)=zw(2,3,ipww)-(fri*rz+fmi*piz+fmj*pjz)*rwc(iax,2) 291 zw(3,3,ipww)=zw(3,3,ipww)-(fri*rz+fmi*piz+fmj*pjz)*rwc(iax,3) 292 ewwpsm=ewwpsm+facu(iax)*(qai*rmj-qaj*rmi)*ri1 293 ewwqsm=ewwqsm+facu(iax)*ri1 294 117 continue 295 etermp=-qfac*qfac*ewwpsm/three 296 eww(8,ipww)=eww(8,ipww)+etermp 297 etermq=q*ewwqsm 298 eww(8,ipww)=eww(8,ipww)+etermq 299c 300c Radial distribution functions 301c 302c if(ifstep-1.eq.((ifstep-1)/nfrdf)*nfrdf .and. ngrww.gt.0) then 303c do 27 igc=1,ngc 304c if(ngt(igc).eq.1) then 305c if(iagc(igc).eq.iwa .and. jagc(igc).eq.jwa) then 306c igr=igrc(igc) 307c do 28 iax=1,nax 308c indx=int(one/(rwi1(iax)*drdf)) 309c if(indx.le.ngl) rdf(indx,igr)=rdf(indx,igr)+rdfvol 310c 28 continue 311c endif 312c endif 313c 27 continue 314c endif 315c 316c Thermodynamic integration 317c 318 if(ithint) then 319 if(ith(2)) then 320 c64=vdw(iwatm(iwa),iwatm(jwa),1,4) 321 c124=vdw(iwatm(iwa),iwatm(jwa),3,4) 322 ewwl6=zero 323 ewwl12=zero 324 do 29 iax=1,nax 325 ewwl6=ewwl6+facu(iax)*rwi6(iax) 326 ewwl12=ewwl12+facu(iax)*rwi6(iax)*rwi6(iax) 327 29 continue 328 deriv(2,ipww)=deriv(2,ipww)+c124*ewwl12-c64*ewwl6 329 endif 330 if(ith(4)) then 331 qi=chg(iwq(iwa),1,iset) 332 qj=chg(iwq(jwa),1,iset) 333 qi4=chg(iwq(iwa),1,4) 334 qj4=chg(iwq(jwa),1,4) 335 dercon=zero 336 if(ipme.eq.0) then 337 do 30 iax=1,nax 338 dercon=dercon+rwi1(iax) 339 30 continue 340 else 341 do 130 iax=1,nax 342 dercon=dercon+rwi1(iax) 343 130 continue 344 endif 345 deriv(4,ipww)=deriv(4,ipww)+(qi*qj4+qj*qi4)*dercon 346 if(ireact.ne.0) then 347 dercon=zero 348 do 31 iax=1,nax 349 dercon=dercon+one/rwi2(iax) 350 31 continue 351 deriv(4,ipww)=deriv(4,ipww)+(qi*qj4+qj*qi4)*rffww*dercon 352 endif 353 endif 354 endif 355c 356c Thermodynamic perturbation 1 357c 358 if(ipert2) then 359 if(ip2(2)) then 360 c6p=vdw(iwatm(iwa),iwatm(jwa),1,2) 361 c12p=vdw(iwatm(iwa),iwatm(jwa),3,2) 362 do 32 iax=1,nax 363 ep2(ipww)=ep2(ipww)+facu(iax)*(c12p*rwi6(iax)-c6p)*rwi6(iax) 364 32 continue 365 ep2(ipww)=ep2(ipww)-eterml 366 endif 367 if(ip2(4).or.ip2(5)) then 368 qp=chg(iwq(iwa),1,2)*chg(iwq(jwa),1,2) 369 ep2tmp=zero 370 do 33 iax=1,nax 371 rwx(iax,1)=xi(iax,1,iwa)-xj(iax,1,jwa) 372 rwx(iax,2)=xi(iax,2,iwa)-xj(iax,2,jwa) 373 rwx(iax,3)=xi(iax,3,iwa)-xj(iax,3,jwa) 374 rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2) 375 rwi1(iax)=sqrt(rwi2(iax)) 376 if(ipme.eq.0) then 377 ep2tmp=ep2tmp+facu(iax)*rwi1(iax) 378 else 379 ep2tmp=ep2tmp+facu(iax)*erfc(ealpha/rwi1(iax))*rwi1(iax) 380 endif 381 33 continue 382 ep2(ipww)=ep2(ipww)+qp*ep2tmp-etermq 383 if(ireact.ne.0) then 384 ep2tmp=zero 385 do 34 iax=1,nax 386 ep2tmp=ep2tmp+facu(iax)/rwi2(iax) 387 34 continue 388 ep2(ipww)=ep2(ipww)+qp*rffww*ep2tmp 389 endif 390 endif 391 endif 392c 393c Thermodynamic perturbation 2 394c 395 if(ipert3) then 396 if(ip3(2)) then 397 c6p=vdw(iwatm(iwa),iwatm(jwa),1,3) 398 c12p=vdw(iwatm(iwa),iwatm(jwa),3,3) 399 do 35 iax=1,nax 400 ep3(ipww)=ep3(ipww)+facu(iax)*(c12p*rwi6(iax)-c6p)*rwi6(iax) 401 35 continue 402 ep3(ipww)=ep3(ipww)-eterml 403 endif 404 if(ip2(4).or.ip2(5)) then 405 qp=chg(iwatm(iwa),1,3)*chg(iwatm(jwa),1,3) 406 ep3tmp=zero 407 do 36 iax=1,nax 408 rwx(iax,1)=xi(iax,1,iwa)-xj(iax,1,jwa) 409 rwx(iax,2)=xi(iax,2,iwa)-xj(iax,2,jwa) 410 rwx(iax,3)=xi(iax,3,iwa)-xj(iax,3,jwa) 411 rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2) 412 rwi1(iax)=sqrt(rwi2(iax)) 413 if(ipme.eq.0) then 414 ep3tmp=ep3tmp+facu(iax)*rwi1(iax) 415 else 416 ep3tmp=ep3tmp+facu(iax)*erfc(ealpha/rwi1(iax))*rwi1(iax) 417 endif 418 36 continue 419 ep3(ipww)=ep3(ipww)+qp*ep3tmp-etermq 420 if(ireact.ne.0) then 421 ep3tmp=zero 422 do 37 iax=1,nax 423 ep3tmp=ep3tmp+facu(iax)/rwi2(iax) 424 37 continue 425 ep3(ipww)=ep3(ipww)+qp*rffww*ep3tmp 426 endif 427 endif 428 endif 429 14 continue 430 13 continue 431c 432c Update force arrays 433c 434 iax=0 435 do 38 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww) 436 iwpj=lwwjpt(iwm,ipww)-1 437 do 39 iwa=1,mwa 438 do 40 iwmn=1,lwwin(iwm,ipww) 439 lwwptr=lwwj(iwpj+iwmn) 440 fw(iwfr+iwm,1,iwa,ipww)=fw(iwfr+iwm,1,iwa,ipww)+fi(iax+iwmn,1,iwa) 441 fw(iwfr+iwm,2,iwa,ipww)=fw(iwfr+iwm,2,iwa,ipww)+fi(iax+iwmn,2,iwa) 442 fw(iwfr+iwm,3,iwa,ipww)=fw(iwfr+iwm,3,iwa,ipww)+fi(iax+iwmn,3,iwa) 443 fw(lwwptr,1,iwa,ipww)=fw(lwwptr,1,iwa,ipww)+fj(iax+iwmn,1,iwa) 444 fw(lwwptr,2,iwa,ipww)=fw(lwwptr,2,iwa,ipww)+fj(iax+iwmn,2,iwa) 445 fw(lwwptr,3,iwa,ipww)=fw(lwwptr,3,iwa,ipww)+fj(iax+iwmn,3,iwa) 446 40 continue 447 39 continue 448c 449c update energy arrays if appropriate print option was set 450c 451c if(npener.ne.0) then 452c do 41 iwmn=1,lwwin(iwm,ipww) 453c lwwptr=lwwj(iwpj+iwmn) 454c uwmw(iwfr+iwm)=uwmw(iwfr+iwm)+u(iax+iwmn) 455c uwmw(lwwptr)=uwmw(lwwptr)+u(iax+iwmn) 456c 41 continue 457c endif 458c 459 iax=iax+lwwin(iwm,ipww) 460 38 continue 461 3 continue 462c 463 1 continue 464c 465 return 466 end 467