1 subroutine argos_cafe_dsw(xs,xsm,ps,psp,isdt,ismf,isml, 2 + isq1,isfrom,nums, 3 + xw,xwm,pw,pwp,iwq,lpbc,chg,lswndx,lswjpt,lswin,lswj, 4 + rwc,xi,xj,rwx,pl,pj,fi,fj,isal) 5c 6c $Id$ 7c 8 implicit none 9c 10#include "argos_cafe_common.fh" 11c 12 real*8 xs(msa,3),xsm(msm,3),ps(msa,3,2),psp(msa,3,2,2) 13 integer isdt(msa),ismf(msa),isml(msa),isq1(msa) 14 integer isfrom,nums 15 real*8 xw(mwm,3,mwa),xwm(mwm,3),pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2) 16 integer iwq(mwa) 17 logical lpbc 18 real*8 chg(mqt,mqp,mset) 19 integer lswndx(0:msa,2),lswjpt(nums,2),lswin(nums,2),lswj(*) 20 real*8 rwc(mscr,3),xi(mscr,3),xj(mscr,3,mwa),rwx(mscr,3) 21 real*8 pl(mscr,3),pj(mscr,3,mwa),fi(mscr,3),fj(mscr,3,mwa) 22 integer isal(mscr) 23c 24 integer ispj,ism,lswptr,ipset,nswlen(2) 25 integer isfr,ipsw,number,isa,jwa,ismn,ispm,iax,nax,nax2 26 real*8 qai,qaj,pai,paj,rx,ry,rz,pix,piy,piz,pjx,pjy,pjz 27 real*8 ri2,ri3,rmi,rmj 28#include "bitops.fh" 29c 30 real*8 qfaci 31c 32 qfaci=one/qfac 33c 34c this subroutine evaluates the solute-solvent forces for nums 35c solute atoms starting from isfrom. the interacting solvent 36c molecules are determined from the pairlist. 37c 38 isfr=isfrom-1 39c 40c loop over short and long range pairs 41c 42 do 1 ipsw=1,npsw 43c 44c evaluate outer index array 45c 46 nswlen(ipsw)=0 47 lswndx(0,ipsw)=0 48 number=0 49 do 2 isa=1,nums 50 if(number+lswin(isa,ipsw).gt.mscr .or. 51 + (ismf(isfr+isa).ne.ismf(isfr+isa-1).and.number.gt.0)) then 52 nswlen(ipsw)=nswlen(ipsw)+1 53 lswndx(nswlen(ipsw),ipsw)=isa-1 54 number=0 55 endif 56 number=number+lswin(isa,ipsw) 57 2 continue 58 if(number.gt.0) then 59 nswlen(ipsw)=nswlen(ipsw)+1 60 lswndx(nswlen(ipsw),ipsw)=nums 61 endif 62c 63c loop over number of cycles to complete pairlist 64c 65 do 3 ispm=1,nswlen(ipsw) 66 nax=0 67c 68c vacuo conditions 69c 70c if(npbtyp.eq.0) then 71c if(.not.lpbc) then 72 do 4 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw) 73 ispj=lswjpt(isa,ipsw)-1 74 ism=isml(isfr+isa) 75c 76c collect center of mass distance vectors 77c 78 if(lpbc.or.ism.eq.0) then 79 do 6 ismn=1,lswin(isa,ipsw) 80 lswptr=lswj(ispj+ismn) 81 rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1) 82 rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2) 83 rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3) 84 6 continue 85 if(lpbc) 86 + call argos_cafe_pbc(0,rwc,mscr,rwx,mscr,nax,1,lswin(isa,ipsw)) 87 else 88 do 5 ismn=1,lswin(isa,ipsw) 89 lswptr=lswj(ispj+ismn) 90 rwc(nax+ismn,1)=xsm(ism,1)-xwm(lswptr,1) 91 rwc(nax+ismn,2)=xsm(ism,2)-xwm(lswptr,2) 92 rwc(nax+ismn,3)=xsm(ism,3)-xwm(lswptr,3) 93 5 continue 94 endif 95c 96c collect solute coordinates and atomic polarization fields 97c 98c if(iand(isdt(isfr+isa),mdynam).eq.ldynam) then 99 do 7 ismn=1,lswin(isa,ipsw) 100 lswptr=lswj(ispj+ismn) 101 xi(nax+ismn,1)=xs(isfr+isa,1) 102 xi(nax+ismn,2)=xs(isfr+isa,2) 103 xi(nax+ismn,3)=xs(isfr+isa,3) 104 pl(nax+ismn,1)=ps(isfr+isa,1,2) 105 pl(nax+ismn,2)=ps(isfr+isa,2,2) 106 pl(nax+ismn,3)=ps(isfr+isa,3,2) 107 isal(nax+ismn)=isfr+isa 108 7 continue 109c else 110c do 8 ismn=1,lswin(isa,ipsw) 111c lswptr=lswj(ispj+ismn) 112c xi(nax+ismn,1)=xs(isfr+isa,1) 113c xi(nax+ismn,2)=xs(isfr+isa,2) 114c xi(nax+ismn,3)=xs(isfr+isa,3) 115c pl(nax+ismn,1)=ps(isfr+isa,1,2) 116c pl(nax+ismn,2)=ps(isfr+isa,2,2) 117c pl(nax+ismn,3)=ps(isfr+isa,3,2) 118c isal(nax+ismn)=isfr+isa 119c 8 continue 120c endif 121c 122c collect solvent coordinates and atomic polarization fields 123c 124 do 8 jwa=1,mwa 125 if(lpbc) then 126 do 9 ismn=1,lswin(isa,ipsw) 127 lswptr=lswj(ispj+ismn) 128 xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa)+rwx(ismn,1) 129 xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa)+rwx(ismn,2) 130 xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa)+rwx(ismn,3) 131 pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2) 132 pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2) 133 pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2) 134 9 continue 135 else 136 do 10 ismn=1,lswin(isa,ipsw) 137 lswptr=lswj(ispj+ismn) 138 xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa) 139 xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa) 140 xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa) 141 pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2) 142 pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2) 143 pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2) 144 10 continue 145 endif 146 8 continue 147 nax=nax+lswin(isa,ipsw) 148 4 continue 149c else 150cc 151cc periodic boundary conditions 152cc 153c do 11 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw) 154c ispj=lswjpt(isa,ipsw)-1 155c ism=isml(isfr+isa) 156cc 157cc collect center of mass distance vectors 158cc 159c do 12 ismn=1,lswin(isa,ipsw) 160c lswptr=lswj(ispj+ismn) 161c rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1) 162c rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2) 163c rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3) 164c rwx(ismn,1)=zero 165c rwx(ismn,2)=zero 166c rwx(ismn,3)=zero 167c if(abs(rwc(nax+ismn,1)).gt.boxh(1)) then 168c rwx(ismn,1)=sign(box(1),xs(isfr+isa,1)) 169c endif 170c if(abs(rwc(nax+ismn,2)).gt.boxh(2)) then 171c rwx(ismn,2)=sign(box(2),xs(isfr+isa,2)) 172c endif 173c if(npbtyp.eq.1) then 174c if(abs(rwc(nax+ismn,3)).gt.boxh(3)) then 175c rwx(ismn,3)=sign(box(3),xs(isfr+isa,3)) 176c endif 177c endif 178c if(ism.gt.0) then 179c rwc(nax+ismn,1)=xsm(ism,1)-xwm(lswptr,1)-rwx(ismn,1) 180c rwc(nax+ismn,2)=xsm(ism,2)-xwm(lswptr,2)-rwx(ismn,2) 181c rwc(nax+ismn,3)=xsm(ism,3)-xwm(lswptr,3)-rwx(ismn,3) 182c else 183c rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1)-rwx(ismn,1) 184c rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2)-rwx(ismn,2) 185c rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3)-rwx(ismn,3) 186c endif 187c 12 continue 188cc 189cc collect solute coordinates and atomic polarization fields 190cc 191c do 13 ismn=1,lswin(isa,ipsw) 192c lswptr=lswj(ispj+ismn) 193c xi(nax+ismn,1)=xs(isfr+isa,1) 194c xi(nax+ismn,2)=xs(isfr+isa,2) 195c xi(nax+ismn,3)=xs(isfr+isa,3) 196c pl(nax+ismn,1)=ps(isfr+isa,1,2) 197c pl(nax+ismn,2)=ps(isfr+isa,2,2) 198c pl(nax+ismn,3)=ps(isfr+isa,3,2) 199c isal(nax+ismn)=isfr+isa 200c 13 continue 201cc 202cc collect solvent coordinates and atomic polarization fields 203cc 204c do 14 jwa=1,mwa 205c do 15 ismn=1,lswin(isa,ipsw) 206c lswptr=lswj(ispj+ismn) 207c xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa)+rwx(ismn,1) 208c xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa)+rwx(ismn,2) 209c xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa)+rwx(ismn,3) 210c pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2) 211c pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2) 212c pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2) 213c 15 continue 214c 14 continue 215c nax=nax+lswin(isa,ipsw) 216c 11 continue 217c endif 218c 219c zero temparary arays fi and fj 220c 221 do 16 iax=1,nax 222 fi(iax,1)=zero 223 fi(iax,2)=zero 224 fi(iax,3)=zero 225 16 continue 226 do 17 jwa=1,mwa 227 do 18 iax=1,nax 228 fj(iax,1,jwa)=zero 229 fj(iax,2,jwa)=zero 230 fj(iax,3,jwa)=zero 231 18 continue 232 17 continue 233c 234c loop over the number of atoms in a solvent molecule 235c 236c calculated here is 4*pi*epsilon*field and not just the field 237c since the polarization is given in alpha/(4*pi*epsilon) in 238c stead of just alpha, the induced dipole is obtained by the 239c product of pwa and pw 240c 241 do 19 jwa=1,mwa 242 qaj=qfaci*chg(iwq(jwa),1,iset) 243 paj=chg(iwq(jwa),2,iset) 244 do 20 iax=1,nax 245 isa=isal(iax) 246 qai=qfaci*chg(isq1(isa),1,iset) 247 pai=chg(isq1(isa),2,iset) 248 rx=xj(iax,1,jwa)-xi(iax,1) 249 ry=xj(iax,2,jwa)-xi(iax,2) 250 rz=xj(iax,3,jwa)-xi(iax,3) 251 pix=pai*pl(iax,1) 252 piy=pai*pl(iax,2) 253 piz=pai*pl(iax,3) 254 pjx=paj*pj(iax,1,jwa) 255 pjy=paj*pj(iax,2,jwa) 256 pjz=paj*pj(iax,3,jwa) 257 ri2=one/(rx**2+ry**2+rz**2) 258 ri3=sqrt(ri2)*ri2 259 rmi=three*(rx*pix+ry*piy+rz*piz)*ri2 260 rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2 261 fi(iax,1)=fi(iax,1)+((rmj-qaj)*rx-pjx)*ri3 262 fi(iax,2)=fi(iax,2)+((rmj-qaj)*ry-pjy)*ri3 263 fi(iax,3)=fi(iax,3)+((rmj-qaj)*rz-pjz)*ri3 264 fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3 265 fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3 266 fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3 267 20 continue 268 19 continue 269c 270c update the electric field arrays 271c 272 iax=0 273 do 21 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw) 274 ispj=lswjpt(isa,ipsw)-1 275 do 22 ismn=1,lswin(isa,ipsw) 276 ps(isfr+isa,1,1)=ps(isfr+isa,1,1)+fi(iax+ismn,1) 277 ps(isfr+isa,2,1)=ps(isfr+isa,2,1)+fi(iax+ismn,2) 278 ps(isfr+isa,3,1)=ps(isfr+isa,3,1)+fi(iax+ismn,3) 279 22 continue 280 do 23 jwa=1,mwa 281 do 24 ismn=1,lswin(isa,ipsw) 282 lswptr=lswj(ispj+ismn) 283 pw(lswptr,1,jwa,1)=pw(lswptr,1,jwa,1)+fj(iax+ismn,1,jwa) 284 pw(lswptr,2,jwa,1)=pw(lswptr,2,jwa,1)+fj(iax+ismn,2,jwa) 285 pw(lswptr,3,jwa,1)=pw(lswptr,3,jwa,1)+fj(iax+ismn,3,jwa) 286 24 continue 287 23 continue 288 iax=iax+lswin(isa,ipsw) 289 21 continue 290c 291c thermodynamic integration and perturbation 292c 293 do 30 ipset=2,3 294 if((ipset.eq.2.and.ipert2).or. 295 + (ipset.eq.3.and.ipert3)) then 296c 297 nax2=0 298 do 31 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw) 299 ispj=lswjpt(isa,ipsw)-1 300 do 32 ismn=1,lswin(isa,ipsw) 301 lswptr=lswj(ispj+ismn) 302 pl(nax2+ismn,1)=psp(isfr+isa,1,2,ipset-1) 303 pl(nax2+ismn,2)=psp(isfr+isa,2,2,ipset-1) 304 pl(nax2+ismn,3)=psp(isfr+isa,3,2,ipset-1) 305 32 continue 306 do 33 jwa=1,mwa 307 do 34 ismn=1,lswin(isa,ipsw) 308 lswptr=lswj(ispj+ismn) 309 pj(nax2+ismn,1,jwa)=pwp(lswptr,1,jwa,2,ipset-1) 310 pj(nax2+ismn,2,jwa)=pwp(lswptr,2,jwa,2,ipset-1) 311 pj(nax2+ismn,3,jwa)=pwp(lswptr,3,jwa,2,ipset-1) 312 34 continue 313 33 continue 314 nax2=nax2+lswin(isa,ipsw) 315 31 continue 316c 317 if(nax.ne.nax2) call md_abort('Error in dipsw',me) 318c 319 do 41 iax=1,nax 320 fi(iax,1)=zero 321 fi(iax,2)=zero 322 fi(iax,3)=zero 323 41 continue 324 do 42 jwa=1,mwa 325 do 43 iax=1,nax 326 fj(iax,1,jwa)=zero 327 fj(iax,2,jwa)=zero 328 fj(iax,3,jwa)=zero 329 43 continue 330 42 continue 331c 332 do 35 jwa=1,mwa 333 qaj=qfaci*chg(iwq(jwa),1,ipset) 334 paj=chg(iwq(jwa),2,ipset) 335 do 36 iax=1,nax 336 isa=isal(iax) 337 qai=qfaci*chg(isq1(isa),1,ipset) 338 pai=chg(isq1(isa),2,ipset) 339 rx=xj(iax,1,jwa)-xi(iax,1) 340 ry=xj(iax,2,jwa)-xi(iax,2) 341 rz=xj(iax,3,jwa)-xi(iax,3) 342 pix=pai*pl(iax,1) 343 piy=pai*pl(iax,2) 344 piz=pai*pl(iax,3) 345 pjx=paj*pj(iax,1,jwa) 346 pjy=paj*pj(iax,2,jwa) 347 pjz=paj*pj(iax,3,jwa) 348 ri2=one/(rx**2+ry**2+rz**2) 349 ri3=sqrt(ri2)*ri2 350 rmi=three*(rx*pix+ry*piy+rz*piz)*ri2 351 rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2 352 fi(iax,1)=fi(iax,1)+((rmj-qaj)*rx-pjx)*ri3 353 fi(iax,2)=fi(iax,2)+((rmj-qaj)*ry-pjy)*ri3 354 fi(iax,3)=fi(iax,3)+((rmj-qaj)*rz-pjz)*ri3 355 fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3 356 fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3 357 fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3 358 36 continue 359 35 continue 360c 361c update the electric field arrays 362c 363 iax=0 364 do 37 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw) 365 ispj=lswjpt(isa,ipsw)-1 366 do 38 ismn=1,lswin(isa,ipsw) 367 psp(isfr+isa,1,ipset-1,1)=psp(isfr+isa,1,ipset-1,1)+fi(iax+ismn,1) 368 psp(isfr+isa,2,ipset-1,1)=psp(isfr+isa,2,ipset-1,1)+fi(iax+ismn,2) 369 psp(isfr+isa,3,ipset-1,1)=psp(isfr+isa,3,ipset-1,1)+fi(iax+ismn,3) 370 38 continue 371 do 39 jwa=1,mwa 372 do 40 ismn=1,lswin(isa,ipsw) 373 lswptr=lswj(ispj+ismn) 374 pwp(lswptr,1,jwa,ipset-1,1)=pwp(lswptr,1,jwa,ipset-1,1)+ 375 + fj(iax+ismn,1,jwa) 376 pwp(lswptr,2,jwa,ipset-1,1)=pwp(lswptr,2,jwa,ipset-1,1)+ 377 + fj(iax+ismn,2,jwa) 378 pwp(lswptr,3,jwa,ipset-1,1)=pwp(lswptr,3,jwa,ipset-1,1)+ 379 + fj(iax+ismn,3,jwa) 380 40 continue 381 39 continue 382 iax=iax+lswin(isa,ipsw) 383 37 continue 384 endif 385 30 continue 386c 387 3 continue 388 1 continue 389c 390 return 391 end 392