1 subroutine cf_dww(xw,xwm,pw,pwp,iwfrom,nlocw,lpbc,chg,iwq, 2 + lwwndx,lwwjpt,lwwin,lwwj,rwc,xi,xj,rwx,pl,pj,fi,fj) 3c 4c $Id$ 5c 6 implicit none 7c 8#include "cf_common.fh" 9c 10 real*8 xw(mwm,3,mwa),xwm(mwm,3) 11 real*8 pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2) 12 integer iwfrom,nlocw 13 logical lpbc 14 real*8 chg(mqt,mqp,mset) 15 integer iwq(mwa) 16 integer lwwndx(0:mwm,2),lwwin(nlocw,2),lwwjpt(nlocw,2),lwwj(*) 17c 18 real*8 rwc(mscr,3),rwx(mscr,3) 19 real*8 xi(mscr,3,mwa),xj(mscr,3,mwa) 20 real*8 pl(mscr,3,mwa),pj(mscr,3,mwa) 21 real*8 fi(mscr,3,mwa),fj(mscr,3,mwa) 22c 23 integer ix,ipset,nax2,nwwlen(2) 24 integer iwfr,ipww,number,iwm,iwpm,nax,iwpj,iwa,jwa,iax,iwmn,lwwptr 25 real*8 pai,paj,qai,qaj 26 real*8 rx,ry,rz,ri2,ri3,rmi,rmj,pix,piy,piz,pjx,pjy,pjz 27c 28 real*8 qfaci 29c 30 qfaci=one/qfac 31c 32c calculation of solvent-solvent intermolecular energies and forces 33c 34 iwfr=iwfrom-1 35c 36c loop over both short and long range parts 37c 38 do 1 ipww=1,npww 39c 40c Evaluate the outer index array 41c 42 nwwlen(ipww)=0 43 lwwndx(0,ipww)=0 44 number=0 45 do 2 iwm=1,nlocw 46 if(number+lwwin(iwm,ipww).gt.mscr) then 47 nwwlen(ipww)=nwwlen(ipww)+1 48 lwwndx(nwwlen(ipww),ipww)=iwm-1 49 number=0 50 endif 51 number=number+lwwin(iwm,ipww) 52 2 continue 53 if(number.gt.0) then 54 nwwlen(ipww)=nwwlen(ipww)+1 55 lwwndx(nwwlen(ipww),ipww)=nlocw 56 endif 57c 58c loop over number of cycles to complete pairlist 59c 60 do 3 iwpm=1,nwwlen(ipww) 61 nax=0 62c 63c collect coordinates into workarrays 64c 65 do 4 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww) 66 iwpj=lwwjpt(iwm,ipww)-1 67 do 5 iwmn=1,lwwin(iwm,ipww) 68 lwwptr=lwwj(iwpj+iwmn) 69 rwc(nax+iwmn,1)=xwm(lwwptr,1)-xwm(iwfr+iwm,1) 70 rwc(nax+iwmn,2)=xwm(lwwptr,2)-xwm(iwfr+iwm,2) 71 rwc(nax+iwmn,3)=xwm(lwwptr,3)-xwm(iwfr+iwm,3) 72 5 continue 73c 74 if(.not.lpbc) then 75 do 6 iwa=1,mwa 76 do 7 iwmn=1,lwwin(iwm,ipww) 77 lwwptr=lwwj(iwpj+iwmn) 78 xi(nax+iwmn,1,iwa)=xw(iwfr+iwm,1,iwa) 79 xi(nax+iwmn,2,iwa)=xw(iwfr+iwm,2,iwa) 80 xi(nax+iwmn,3,iwa)=xw(iwfr+iwm,3,iwa) 81 xj(nax+iwmn,1,iwa)=xw(lwwptr,1,iwa) 82 xj(nax+iwmn,2,iwa)=xw(lwwptr,2,iwa) 83 xj(nax+iwmn,3,iwa)=xw(lwwptr,3,iwa) 84 pl(nax+iwmn,1,iwa)=pw(iwfr+iwm,1,iwa,2) 85 pl(nax+iwmn,2,iwa)=pw(iwfr+iwm,2,iwa,2) 86 pl(nax+iwmn,3,iwa)=pw(iwfr+iwm,3,iwa,2) 87 pj(nax+iwmn,1,iwa)=pw(lwwptr,1,iwa,2) 88 pj(nax+iwmn,2,iwa)=pw(lwwptr,2,iwa,2) 89 pj(nax+iwmn,3,iwa)=pw(lwwptr,3,iwa,2) 90 7 continue 91 6 continue 92 else 93 call cf_pbc(0,rwc,mscr,rwx,mscr,nax,1,lwwin(iwm,ipww)) 94 do 9 iwmn=1,lwwin(iwm,ipww) 95 lwwptr=lwwj(iwpj+iwmn) 96 rwc(nax+iwmn,1)=rwc(nax+iwmn,1)-rwx(iwmn,1) 97 rwc(nax+iwmn,2)=rwc(nax+iwmn,2)-rwx(iwmn,2) 98 rwc(nax+iwmn,3)=rwc(nax+iwmn,3)-rwx(iwmn,3) 99 9 continue 100c 101 do 10 iwa=1,mwa 102 do 11 iwmn=1,lwwin(iwm,ipww) 103 lwwptr=lwwj(iwpj+iwmn) 104 xi(nax+iwmn,1,iwa)=xw(iwfr+iwm,1,iwa) 105 xi(nax+iwmn,2,iwa)=xw(iwfr+iwm,2,iwa) 106 xi(nax+iwmn,3,iwa)=xw(iwfr+iwm,3,iwa) 107 xj(nax+iwmn,1,iwa)=xw(lwwptr,1,iwa)-rwx(iwmn,1) 108 xj(nax+iwmn,2,iwa)=xw(lwwptr,2,iwa)-rwx(iwmn,2) 109 xj(nax+iwmn,3,iwa)=xw(lwwptr,3,iwa)-rwx(iwmn,3) 110 pl(nax+iwmn,1,iwa)=pw(iwfr+iwm,1,iwa,2) 111 pl(nax+iwmn,2,iwa)=pw(iwfr+iwm,2,iwa,2) 112 pl(nax+iwmn,3,iwa)=pw(iwfr+iwm,3,iwa,2) 113 pj(nax+iwmn,1,iwa)=pw(lwwptr,1,iwa,2) 114 pj(nax+iwmn,2,iwa)=pw(lwwptr,2,iwa,2) 115 pj(nax+iwmn,3,iwa)=pw(lwwptr,3,iwa,2) 116 11 continue 117 10 continue 118 endif 119 nax=nax+lwwin(iwm,ipww) 120 4 continue 121c 122c zero temporary arrays fi and fj 123c 124 do 12 iwa=1,mwa 125 do 13 ix=1,3 126 do 14 iax=1,nax 127 fi(iax,ix,iwa)=zero 128 fj(iax,ix,iwa)=zero 129 14 continue 130 13 continue 131 12 continue 132c 133c loops over number of atoms in a solvent molecule 134c 135c calculated here is 4*pi*epsilon*field and not just the field 136c since the polarization is given in alpha/(4*pi*epsilon) in 137c stead of just alpha, the induced dipole is obtained by the 138c product of pwa and pw 139c 140 do 20 iwa=1,mwa 141 qai=qfaci*chg(iwq(iwa),1,iset) 142 pai=chg(iwq(iwa),2,iset) 143 do 21 jwa=1,mwa 144 qaj=qfaci*chg(iwq(jwa),1,iset) 145 paj=chg(iwq(jwa),2,iset) 146 do 22 iax=1,nax 147 rx=xj(iax,1,jwa)-xi(iax,1,iwa) 148 ry=xj(iax,2,jwa)-xi(iax,2,iwa) 149 rz=xj(iax,3,jwa)-xi(iax,3,iwa) 150 pix=pai*pl(iax,1,iwa) 151 piy=pai*pl(iax,2,iwa) 152 piz=pai*pl(iax,3,iwa) 153 pjx=paj*pj(iax,1,jwa) 154 pjy=paj*pj(iax,2,jwa) 155 pjz=paj*pj(iax,3,jwa) 156 ri2=one/(rx**2+ry**2+rz**2) 157 ri3=sqrt(ri2)*ri2 158 rmi=three*(rx*pix+ry*piy+rz*piz)*ri2 159 rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2 160 fi(iax,1,iwa)=fi(iax,1,iwa)+((rmj-qaj)*rx-pjx)*ri3 161 fi(iax,2,iwa)=fi(iax,2,iwa)+((rmj-qaj)*ry-pjy)*ri3 162 fi(iax,3,iwa)=fi(iax,3,iwa)+((rmj-qaj)*rz-pjz)*ri3 163 fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3 164 fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3 165 fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3 166 22 continue 167 21 continue 168 20 continue 169c 170c Update the electric field arrays 171c 172 iax=0 173 do 23 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww) 174 iwpj=lwwjpt(iwm,ipww)-1 175 do 24 iwa=1,mwa 176 do 25 iwmn=1,lwwin(iwm,ipww) 177 lwwptr=lwwj(iwpj+iwmn) 178 pw(iwfr+iwm,1,iwa,1)=pw(iwfr+iwm,1,iwa,1)+fi(iax+iwmn,1,iwa) 179 pw(iwfr+iwm,2,iwa,1)=pw(iwfr+iwm,2,iwa,1)+fi(iax+iwmn,2,iwa) 180 pw(iwfr+iwm,3,iwa,1)=pw(iwfr+iwm,3,iwa,1)+fi(iax+iwmn,3,iwa) 181 pw(lwwptr,1,iwa,1)=pw(lwwptr,1,iwa,1)+fj(iax+iwmn,1,iwa) 182 pw(lwwptr,2,iwa,1)=pw(lwwptr,2,iwa,1)+fj(iax+iwmn,2,iwa) 183 pw(lwwptr,3,iwa,1)=pw(lwwptr,3,iwa,1)+fj(iax+iwmn,3,iwa) 184 25 continue 185 24 continue 186 iax=iax+lwwin(iwm,ipww) 187 23 continue 188c 189c thermodynamic perturbation and integration 190c 191 do 30 ipset=2,3 192 if((ipset.eq.2.and.ipert2).or. 193 + (ipset.eq.3.and.ipert3)) then 194c 195 nax2=0 196 do 31 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww) 197 iwpj=lwwjpt(iwm,ipww)-1 198 do 32 iwa=1,mwa 199 do 33 iwmn=1,lwwin(iwm,ipww) 200 lwwptr=lwwj(iwpj+iwmn) 201 pl(nax2+iwmn,1,iwa)=pwp(iwfr+iwm,1,iwa,ipset-1,2) 202 pl(nax2+iwmn,2,iwa)=pwp(iwfr+iwm,2,iwa,ipset-1,2) 203 pl(nax2+iwmn,3,iwa)=pwp(iwfr+iwm,3,iwa,ipset-1,2) 204 pj(nax2+iwmn,1,iwa)=pwp(lwwptr,1,iwa,ipset-1,2) 205 pj(nax2+iwmn,2,iwa)=pwp(lwwptr,2,iwa,ipset-1,2) 206 pj(nax2+iwmn,3,iwa)=pwp(lwwptr,3,iwa,ipset-1,2) 207 33 continue 208 32 continue 209 nax2=nax2+lwwin(iwm,ipww) 210 31 continue 211c 212 if(nax.ne.nax2) call md_abort('Error in dipww',me) 213c 214c 215 do 40 iwa=1,mwa 216 do 41 ix=1,3 217 do 42 iax=1,nax 218 fi(iax,ix,iwa)=zero 219 fj(iax,ix,iwa)=zero 220 42 continue 221 41 continue 222 40 continue 223c 224 do 34 iwa=1,mwa 225 qai=qfaci*chg(iwq(iwa),1,ipset) 226 pai=chg(iwq(iwa),2,ipset) 227 do 35 jwa=1,mwa 228 qaj=qfaci*chg(iwq(jwa),1,ipset) 229 paj=chg(iwq(jwa),2,ipset) 230 do 36 iax=1,nax 231 rx=xj(iax,1,jwa)-xi(iax,1,iwa) 232 ry=xj(iax,2,jwa)-xi(iax,2,iwa) 233 rz=xj(iax,3,jwa)-xi(iax,3,iwa) 234 pix=pai*pl(iax,1,iwa) 235 piy=pai*pl(iax,2,iwa) 236 piz=pai*pl(iax,3,iwa) 237 pjx=paj*pj(iax,1,jwa) 238 pjy=paj*pj(iax,2,jwa) 239 pjz=paj*pj(iax,3,jwa) 240 ri2=one/(rx**2+ry**2+rz**2) 241 ri3=sqrt(ri2)*ri2 242 rmi=three*(rx*pix+ry*piy+rz*piz)*ri2 243 rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2 244 fi(iax,1,iwa)=fi(iax,1,iwa)+((rmj-qaj)*rx-pjx)*ri3 245 fi(iax,2,iwa)=fi(iax,2,iwa)+((rmj-qaj)*ry-pjy)*ri3 246 fi(iax,3,iwa)=fi(iax,3,iwa)+((rmj-qaj)*rz-pjz)*ri3 247 fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3 248 fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3 249 fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3 250 36 continue 251 35 continue 252 34 continue 253c 254c Update the electric field arrays 255c 256 iax=0 257 do 37 iwm=lwwndx(iwpm-1,ipww)+1,lwwndx(iwpm,ipww) 258 iwpj=lwwjpt(iwm,ipww)-1 259 do 38 iwa=1,mwa 260 do 39 iwmn=1,lwwin(iwm,ipww) 261 lwwptr=lwwj(iwpj+iwmn) 262 pwp(iwfr+iwm,1,iwa,ipset-1,1)=pwp(iwfr+iwm,1,iwa,ipset-1,1)+ 263 + fi(iax+iwmn,1,iwa) 264 pwp(iwfr+iwm,2,iwa,ipset-1,1)=pwp(iwfr+iwm,2,iwa,ipset-1,1)+ 265 + fi(iax+iwmn,2,iwa) 266 pwp(iwfr+iwm,3,iwa,ipset-1,1)=pwp(iwfr+iwm,3,iwa,ipset-1,1)+ 267 + fi(iax+iwmn,3,iwa) 268 pwp(lwwptr,1,iwa,ipset-1,1)=pwp(lwwptr,1,iwa,ipset-1,1)+ 269 + fj(iax+iwmn,1,iwa) 270 pwp(lwwptr,2,iwa,ipset-1,1)=pwp(lwwptr,2,iwa,ipset-1,1)+ 271 + fj(iax+iwmn,2,iwa) 272 pwp(lwwptr,3,iwa,ipset-1,1)=pwp(lwwptr,3,iwa,ipset-1,1)+ 273 + fj(iax+iwmn,3,iwa) 274 39 continue 275 38 continue 276 iax=iax+lwwin(iwm,ipww) 277 37 continue 278 endif 279 30 continue 280c 281c 282 3 continue 283 1 continue 284c 285 return 286 end 287 subroutine cf_dsw(xs,xsm,ps,psp,isdt,ismf,isml,isq1,isfrom,nums, 288 + xw,xwm,pw,pwp,iwq,lpbc,chg,lswndx,lswjpt,lswin,lswj, 289 + rwc,xi,xj,rwx,pl,pj,fi,fj,isal) 290c 291c $Id$ 292c 293 implicit none 294c 295#include "cf_common.fh" 296c 297 real*8 xs(msa,3),xsm(msm,3),ps(msa,3,2),psp(msa,3,2,2) 298 integer isdt(msa),ismf(msa),isml(msa),isq1(msa) 299 integer isfrom,nums 300 real*8 xw(mwm,3,mwa),xwm(mwm,3),pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2) 301 integer iwq(mwa) 302 logical lpbc 303 real*8 chg(mqt,mqp,mset) 304 integer lswndx(0:msa,2),lswjpt(nums,2),lswin(nums,2),lswj(*) 305 real*8 rwc(mscr,3),xi(mscr,3),xj(mscr,3,mwa),rwx(mscr,3) 306 real*8 pl(mscr,3),pj(mscr,3,mwa),fi(mscr,3),fj(mscr,3,mwa) 307 integer isal(mscr) 308c 309 integer ispj,ism,lswptr,ipset,nswlen(2) 310 integer isfr,ipsw,number,isa,jwa,ismn,ispm,iax,nax,nax2 311 real*8 qai,qaj,pai,paj,rx,ry,rz,pix,piy,piz,pjx,pjy,pjz 312 real*8 ri2,ri3,rmi,rmj 313#include "bitops.fh" 314c 315 real*8 qfaci 316c 317 qfaci=one/qfac 318c 319c this subroutine evaluates the solute-solvent forces for nums 320c solute atoms starting from isfrom. the interacting solvent 321c molecules are determined from the pairlist. 322c 323 isfr=isfrom-1 324c 325c loop over short and long range pairs 326c 327 do 1 ipsw=1,npsw 328c 329c evaluate outer index array 330c 331 nswlen(ipsw)=0 332 lswndx(0,ipsw)=0 333 number=0 334 do 2 isa=1,nums 335 if(number+lswin(isa,ipsw).gt.mscr .or. 336 + (ismf(isfr+isa).ne.ismf(isfr+isa-1).and.number.gt.0)) then 337 nswlen(ipsw)=nswlen(ipsw)+1 338 lswndx(nswlen(ipsw),ipsw)=isa-1 339 number=0 340 endif 341 number=number+lswin(isa,ipsw) 342 2 continue 343 if(number.gt.0) then 344 nswlen(ipsw)=nswlen(ipsw)+1 345 lswndx(nswlen(ipsw),ipsw)=nums 346 endif 347c 348c loop over number of cycles to complete pairlist 349c 350 do 3 ispm=1,nswlen(ipsw) 351 nax=0 352c 353c vacuo conditions 354c 355c if(npbtyp.eq.0) then 356c if(.not.lpbc) then 357 do 4 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw) 358 ispj=lswjpt(isa,ipsw)-1 359 ism=isml(isfr+isa) 360c 361c collect center of mass distance vectors 362c 363 if(lpbc.or.ism.eq.0) then 364 do 6 ismn=1,lswin(isa,ipsw) 365 lswptr=lswj(ispj+ismn) 366 rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1) 367 rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2) 368 rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3) 369 6 continue 370 if(lpbc) call cf_pbc(0,rwc,mscr,rwx,mscr,nax,1,lswin(isa,ipsw)) 371 else 372 do 5 ismn=1,lswin(isa,ipsw) 373 lswptr=lswj(ispj+ismn) 374 rwc(nax+ismn,1)=xsm(ism,1)-xwm(lswptr,1) 375 rwc(nax+ismn,2)=xsm(ism,2)-xwm(lswptr,2) 376 rwc(nax+ismn,3)=xsm(ism,3)-xwm(lswptr,3) 377 5 continue 378 endif 379c 380c collect solute coordinates and atomic polarization fields 381c 382c if(iand(isdt(isfr+isa),mdynam).eq.ldynam) then 383 do 7 ismn=1,lswin(isa,ipsw) 384 lswptr=lswj(ispj+ismn) 385 xi(nax+ismn,1)=xs(isfr+isa,1) 386 xi(nax+ismn,2)=xs(isfr+isa,2) 387 xi(nax+ismn,3)=xs(isfr+isa,3) 388 pl(nax+ismn,1)=ps(isfr+isa,1,2) 389 pl(nax+ismn,2)=ps(isfr+isa,2,2) 390 pl(nax+ismn,3)=ps(isfr+isa,3,2) 391 isal(nax+ismn)=isfr+isa 392 7 continue 393c else 394c do 8 ismn=1,lswin(isa,ipsw) 395c lswptr=lswj(ispj+ismn) 396c xi(nax+ismn,1)=xs(isfr+isa,1) 397c xi(nax+ismn,2)=xs(isfr+isa,2) 398c xi(nax+ismn,3)=xs(isfr+isa,3) 399c pl(nax+ismn,1)=ps(isfr+isa,1,2) 400c pl(nax+ismn,2)=ps(isfr+isa,2,2) 401c pl(nax+ismn,3)=ps(isfr+isa,3,2) 402c isal(nax+ismn)=isfr+isa 403c 8 continue 404c endif 405c 406c collect solvent coordinates and atomic polarization fields 407c 408 do 8 jwa=1,mwa 409 if(lpbc) then 410 do 9 ismn=1,lswin(isa,ipsw) 411 lswptr=lswj(ispj+ismn) 412 xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa)+rwx(ismn,1) 413 xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa)+rwx(ismn,2) 414 xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa)+rwx(ismn,3) 415 pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2) 416 pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2) 417 pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2) 418 9 continue 419 else 420 do 10 ismn=1,lswin(isa,ipsw) 421 lswptr=lswj(ispj+ismn) 422 xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa) 423 xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa) 424 xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa) 425 pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2) 426 pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2) 427 pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2) 428 10 continue 429 endif 430 8 continue 431 nax=nax+lswin(isa,ipsw) 432 4 continue 433c else 434cc 435cc periodic boundary conditions 436cc 437c do 11 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw) 438c ispj=lswjpt(isa,ipsw)-1 439c ism=isml(isfr+isa) 440cc 441cc collect center of mass distance vectors 442cc 443c do 12 ismn=1,lswin(isa,ipsw) 444c lswptr=lswj(ispj+ismn) 445c rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1) 446c rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2) 447c rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3) 448c rwx(ismn,1)=zero 449c rwx(ismn,2)=zero 450c rwx(ismn,3)=zero 451c if(abs(rwc(nax+ismn,1)).gt.boxh(1)) then 452c rwx(ismn,1)=sign(box(1),xs(isfr+isa,1)) 453c endif 454c if(abs(rwc(nax+ismn,2)).gt.boxh(2)) then 455c rwx(ismn,2)=sign(box(2),xs(isfr+isa,2)) 456c endif 457c if(npbtyp.eq.1) then 458c if(abs(rwc(nax+ismn,3)).gt.boxh(3)) then 459c rwx(ismn,3)=sign(box(3),xs(isfr+isa,3)) 460c endif 461c endif 462c if(ism.gt.0) then 463c rwc(nax+ismn,1)=xsm(ism,1)-xwm(lswptr,1)-rwx(ismn,1) 464c rwc(nax+ismn,2)=xsm(ism,2)-xwm(lswptr,2)-rwx(ismn,2) 465c rwc(nax+ismn,3)=xsm(ism,3)-xwm(lswptr,3)-rwx(ismn,3) 466c else 467c rwc(nax+ismn,1)=xs(isfr+isa,1)-xwm(lswptr,1)-rwx(ismn,1) 468c rwc(nax+ismn,2)=xs(isfr+isa,2)-xwm(lswptr,2)-rwx(ismn,2) 469c rwc(nax+ismn,3)=xs(isfr+isa,3)-xwm(lswptr,3)-rwx(ismn,3) 470c endif 471c 12 continue 472cc 473cc collect solute coordinates and atomic polarization fields 474cc 475c do 13 ismn=1,lswin(isa,ipsw) 476c lswptr=lswj(ispj+ismn) 477c xi(nax+ismn,1)=xs(isfr+isa,1) 478c xi(nax+ismn,2)=xs(isfr+isa,2) 479c xi(nax+ismn,3)=xs(isfr+isa,3) 480c pl(nax+ismn,1)=ps(isfr+isa,1,2) 481c pl(nax+ismn,2)=ps(isfr+isa,2,2) 482c pl(nax+ismn,3)=ps(isfr+isa,3,2) 483c isal(nax+ismn)=isfr+isa 484c 13 continue 485cc 486cc collect solvent coordinates and atomic polarization fields 487cc 488c do 14 jwa=1,mwa 489c do 15 ismn=1,lswin(isa,ipsw) 490c lswptr=lswj(ispj+ismn) 491c xj(nax+ismn,1,jwa)=xw(lswptr,1,jwa)+rwx(ismn,1) 492c xj(nax+ismn,2,jwa)=xw(lswptr,2,jwa)+rwx(ismn,2) 493c xj(nax+ismn,3,jwa)=xw(lswptr,3,jwa)+rwx(ismn,3) 494c pj(nax+ismn,1,jwa)=pw(lswptr,1,jwa,2) 495c pj(nax+ismn,2,jwa)=pw(lswptr,2,jwa,2) 496c pj(nax+ismn,3,jwa)=pw(lswptr,3,jwa,2) 497c 15 continue 498c 14 continue 499c nax=nax+lswin(isa,ipsw) 500c 11 continue 501c endif 502c 503c zero temparary arays fi and fj 504c 505 do 16 iax=1,nax 506 fi(iax,1)=zero 507 fi(iax,2)=zero 508 fi(iax,3)=zero 509 16 continue 510 do 17 jwa=1,mwa 511 do 18 iax=1,nax 512 fj(iax,1,jwa)=zero 513 fj(iax,2,jwa)=zero 514 fj(iax,3,jwa)=zero 515 18 continue 516 17 continue 517c 518c loop over the number of atoms in a solvent molecule 519c 520c calculated here is 4*pi*epsilon*field and not just the field 521c since the polarization is given in alpha/(4*pi*epsilon) in 522c stead of just alpha, the induced dipole is obtained by the 523c product of pwa and pw 524c 525 do 19 jwa=1,mwa 526 qaj=qfaci*chg(iwq(jwa),1,iset) 527 paj=chg(iwq(jwa),2,iset) 528 do 20 iax=1,nax 529 isa=isal(iax) 530 qai=qfaci*chg(isq1(isa),1,iset) 531 pai=chg(isq1(isa),2,iset) 532 rx=xj(iax,1,jwa)-xi(iax,1) 533 ry=xj(iax,2,jwa)-xi(iax,2) 534 rz=xj(iax,3,jwa)-xi(iax,3) 535 pix=pai*pl(iax,1) 536 piy=pai*pl(iax,2) 537 piz=pai*pl(iax,3) 538 pjx=paj*pj(iax,1,jwa) 539 pjy=paj*pj(iax,2,jwa) 540 pjz=paj*pj(iax,3,jwa) 541 ri2=one/(rx**2+ry**2+rz**2) 542 ri3=sqrt(ri2)*ri2 543 rmi=three*(rx*pix+ry*piy+rz*piz)*ri2 544 rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2 545 fi(iax,1)=fi(iax,1)+((rmj-qaj)*rx-pjx)*ri3 546 fi(iax,2)=fi(iax,2)+((rmj-qaj)*ry-pjy)*ri3 547 fi(iax,3)=fi(iax,3)+((rmj-qaj)*rz-pjz)*ri3 548 fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3 549 fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3 550 fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3 551 20 continue 552 19 continue 553c 554c update the electric field arrays 555c 556 iax=0 557 do 21 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw) 558 ispj=lswjpt(isa,ipsw)-1 559 do 22 ismn=1,lswin(isa,ipsw) 560 ps(isfr+isa,1,1)=ps(isfr+isa,1,1)+fi(iax+ismn,1) 561 ps(isfr+isa,2,1)=ps(isfr+isa,2,1)+fi(iax+ismn,2) 562 ps(isfr+isa,3,1)=ps(isfr+isa,3,1)+fi(iax+ismn,3) 563 22 continue 564 do 23 jwa=1,mwa 565 do 24 ismn=1,lswin(isa,ipsw) 566 lswptr=lswj(ispj+ismn) 567 pw(lswptr,1,jwa,1)=pw(lswptr,1,jwa,1)+fj(iax+ismn,1,jwa) 568 pw(lswptr,2,jwa,1)=pw(lswptr,2,jwa,1)+fj(iax+ismn,2,jwa) 569 pw(lswptr,3,jwa,1)=pw(lswptr,3,jwa,1)+fj(iax+ismn,3,jwa) 570 24 continue 571 23 continue 572 iax=iax+lswin(isa,ipsw) 573 21 continue 574c 575c thermodynamic integration and perturbation 576c 577 do 30 ipset=2,3 578 if((ipset.eq.2.and.ipert2).or. 579 + (ipset.eq.3.and.ipert3)) then 580c 581 nax2=0 582 do 31 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw) 583 ispj=lswjpt(isa,ipsw)-1 584 do 32 ismn=1,lswin(isa,ipsw) 585 lswptr=lswj(ispj+ismn) 586 pl(nax2+ismn,1)=psp(isfr+isa,1,2,ipset-1) 587 pl(nax2+ismn,2)=psp(isfr+isa,2,2,ipset-1) 588 pl(nax2+ismn,3)=psp(isfr+isa,3,2,ipset-1) 589 32 continue 590 do 33 jwa=1,mwa 591 do 34 ismn=1,lswin(isa,ipsw) 592 lswptr=lswj(ispj+ismn) 593 pj(nax2+ismn,1,jwa)=pwp(lswptr,1,jwa,2,ipset-1) 594 pj(nax2+ismn,2,jwa)=pwp(lswptr,2,jwa,2,ipset-1) 595 pj(nax2+ismn,3,jwa)=pwp(lswptr,3,jwa,2,ipset-1) 596 34 continue 597 33 continue 598 nax2=nax2+lswin(isa,ipsw) 599 31 continue 600c 601 if(nax.ne.nax2) call md_abort('Error in dipsw',me) 602c 603 do 41 iax=1,nax 604 fi(iax,1)=zero 605 fi(iax,2)=zero 606 fi(iax,3)=zero 607 41 continue 608 do 42 jwa=1,mwa 609 do 43 iax=1,nax 610 fj(iax,1,jwa)=zero 611 fj(iax,2,jwa)=zero 612 fj(iax,3,jwa)=zero 613 43 continue 614 42 continue 615c 616 do 35 jwa=1,mwa 617 qaj=qfaci*chg(iwq(jwa),1,ipset) 618 paj=chg(iwq(jwa),2,ipset) 619 do 36 iax=1,nax 620 isa=isal(iax) 621 qai=qfaci*chg(isq1(isa),1,ipset) 622 pai=chg(isq1(isa),2,ipset) 623 rx=xj(iax,1,jwa)-xi(iax,1) 624 ry=xj(iax,2,jwa)-xi(iax,2) 625 rz=xj(iax,3,jwa)-xi(iax,3) 626 pix=pai*pl(iax,1) 627 piy=pai*pl(iax,2) 628 piz=pai*pl(iax,3) 629 pjx=paj*pj(iax,1,jwa) 630 pjy=paj*pj(iax,2,jwa) 631 pjz=paj*pj(iax,3,jwa) 632 ri2=one/(rx**2+ry**2+rz**2) 633 ri3=sqrt(ri2)*ri2 634 rmi=three*(rx*pix+ry*piy+rz*piz)*ri2 635 rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2 636 fi(iax,1)=fi(iax,1)+((rmj-qaj)*rx-pjx)*ri3 637 fi(iax,2)=fi(iax,2)+((rmj-qaj)*ry-pjy)*ri3 638 fi(iax,3)=fi(iax,3)+((rmj-qaj)*rz-pjz)*ri3 639 fj(iax,1,jwa)=fj(iax,1,jwa)+((rmi+qai)*rx-pix)*ri3 640 fj(iax,2,jwa)=fj(iax,2,jwa)+((rmi+qai)*ry-piy)*ri3 641 fj(iax,3,jwa)=fj(iax,3,jwa)+((rmi+qai)*rz-piz)*ri3 642 36 continue 643 35 continue 644c 645c update the electric field arrays 646c 647 iax=0 648 do 37 isa=lswndx(ispm-1,ipsw)+1,lswndx(ispm,ipsw) 649 ispj=lswjpt(isa,ipsw)-1 650 do 38 ismn=1,lswin(isa,ipsw) 651 psp(isfr+isa,1,ipset-1,1)=psp(isfr+isa,1,ipset-1,1)+fi(iax+ismn,1) 652 psp(isfr+isa,2,ipset-1,1)=psp(isfr+isa,2,ipset-1,1)+fi(iax+ismn,2) 653 psp(isfr+isa,3,ipset-1,1)=psp(isfr+isa,3,ipset-1,1)+fi(iax+ismn,3) 654 38 continue 655 do 39 jwa=1,mwa 656 do 40 ismn=1,lswin(isa,ipsw) 657 lswptr=lswj(ispj+ismn) 658 pwp(lswptr,1,jwa,ipset-1,1)=pwp(lswptr,1,jwa,ipset-1,1)+ 659 + fj(iax+ismn,1,jwa) 660 pwp(lswptr,2,jwa,ipset-1,1)=pwp(lswptr,2,jwa,ipset-1,1)+ 661 + fj(iax+ismn,2,jwa) 662 pwp(lswptr,3,jwa,ipset-1,1)=pwp(lswptr,3,jwa,ipset-1,1)+ 663 + fj(iax+ismn,3,jwa) 664 40 continue 665 39 continue 666 iax=iax+lswin(isa,ipsw) 667 37 continue 668 endif 669 30 continue 670c 671 3 continue 672 1 continue 673c 674 return 675 end 676 subroutine cf_dss(xs,xsm,ps,psp,ismf,isml,isq2,isq3,isfrom,nums, 677 + lpbc,chg,lssndx,lssjpt,lssin,lssj,xi,xj,rwc,rwi1,rwi2,rwi6, 678 + rwx,rw,fi,fj,f,isal,jsal,jmal,jfal,pl,pj) 679c 680c $Id$ 681c 682 implicit none 683c 684#include "cf_common.fh" 685c 686 real*8 xs(msa,3),xsm(msm,3),ps(msa,3,2),psp(msa,3,2,2) 687 integer ismf(msa),isml(msa),isq2(msa),isq3(msa) 688 integer isfrom,nums 689 logical lpbc 690 real*8 chg(mqt,mqp,mset) 691 integer lssndx(0:msa,2),lssjpt(nums,2),lssin(nums,2) 692c 693 real*8 xi(mscr,3),xj(mscr,3),rwx(mscr,3),rwi1(mscr) 694 real*8 rwi2(mscr),rwi6(mscr),rwc(mscr,3),rw(mscr) 695 real*8 f(mscr),fi(mscr,3),fj(mscr,3),pl(mscr,3),pj(mscr,3) 696 integer isal(mscr),jsal(mscr),jmal(mscr),jfal(mscr) 697 integer lssj(*) 698c 699 integer isa,jsa,jsf,ipset 700 integer isfr,nsslen(2) 701 integer ipss,number,isslen,nax,nax2,jsaptr 702 integer jnum,lssptr,iax,ism,jsm 703c 704 real*8 ri1,ri2,ri3,rx,ry,rz,pix,piy,piz,pjx,pjy,pjz,rmi,rmj 705 real*8 qai,pai,qaj,paj 706c 707#include "bitops.fh" 708c 709 real*8 qfaci 710c 711 qfaci=one/qfac 712c 713c solute non-bonded pairs 714c ======================= 715c 716 isfr=isfrom-1 717c 718c loop over both short and long range pairlists 719c 720 do 1 ipss=1,npss 721c 722c evaluate outer index array 723c 724 nsslen(ipss)=0 725 lssndx(0,ipss)=0 726 number=0 727 do 2 isa=1,nums 728 if(number+lssin(isa,ipss).gt.mscr .or. 729 + (ismf(isfr+isa).ne.ismf(isfr+isa-1).and.number.gt.0)) then 730 nsslen(ipss)=nsslen(ipss)+1 731 lssndx(nsslen(ipss),ipss)=isa-1 732 number=0 733 endif 734 number=number+lssin(isa,ipss) 735 2 continue 736 if(number.gt.0) then 737 nsslen(ipss)=nsslen(ipss)+1 738 lssndx(nsslen(ipss),ipss)=nums 739 endif 740c 741c loop over number of cycles to complete pairlists 742c 743 do 3 isslen=1,nsslen(ipss) 744 nax=0 745 ism=isml(isfr+lssndx(isslen,ipss)) 746c 747c collect coordinates into workarrays 748c 749c if(npbtyp.eq.0) then 750 do 4 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss) 751 jsaptr=lssjpt(isa,ipss)-1 752 ism=isml(isfr+isa) 753 if(lpbc) then 754 do 151 jnum=1,lssin(isa,ipss) 755 lssptr=lssj(jsaptr+jnum) 756 rwc(nax+jnum,1)=xs(isfr+isa,1)-xs(lssptr,1) 757 rwc(nax+jnum,2)=xs(isfr+isa,2)-xs(lssptr,2) 758 rwc(nax+jnum,3)=xs(isfr+isa,3)-xs(lssptr,3) 759 151 continue 760 call cf_pbc(0,rwc,mscr,rwx,mscr,nax,1,lssin(isa,ipss)) 761 endif 762 do 5 jnum=1,lssin(isa,ipss) 763 lssptr=lssj(jsaptr+jnum) 764 jsf=ismf(lssptr) 765 isal(nax+jnum)=isfr+isa 766 jsal(nax+jnum)=lssptr 767 jmal(nax+jnum)=jsf 768 jsm=isml(lssptr) 769 jfal(nax+jnum)=2 770 if(ism.ne.jsm) jfal(nax+jnum)=3 771 if(ism.gt.0) then 772 if(jsm.gt.0) then 773 rwc(nax+jnum,1)=xsm(ism,1)-xsm(jsm,1) 774 rwc(nax+jnum,2)=xsm(ism,2)-xsm(jsm,2) 775 rwc(nax+jnum,3)=xsm(ism,3)-xsm(jsm,3) 776 else 777 rwc(nax+jnum,1)=xsm(ism,1)-xs(lssptr,1) 778 rwc(nax+jnum,2)=xsm(ism,2)-xs(lssptr,2) 779 rwc(nax+jnum,3)=xsm(ism,3)-xs(lssptr,3) 780 endif 781 else 782 if(jsm.gt.0) then 783 rwc(nax+jnum,1)=xs(isfr+isa,1)-xsm(jsm,1) 784 rwc(nax+jnum,2)=xs(isfr+isa,2)-xsm(jsm,2) 785 rwc(nax+jnum,3)=xs(isfr+isa,3)-xsm(jsm,3) 786 else 787 rwc(nax+jnum,1)=xs(isfr+isa,1)-xs(lssptr,1) 788 rwc(nax+jnum,2)=xs(isfr+isa,2)-xs(lssptr,2) 789 rwc(nax+jnum,3)=xs(isfr+isa,3)-xs(lssptr,3) 790 endif 791 endif 792 5 continue 793 if(.not.lpbc) then 794 do 7 jnum=1,lssin(isa,ipss) 795 lssptr=lssj(jsaptr+jnum) 796 xi(nax+jnum,1)=xs(isfr+isa,1) 797 xi(nax+jnum,2)=xs(isfr+isa,2) 798 xi(nax+jnum,3)=xs(isfr+isa,3) 799 xj(nax+jnum,1)=xs(lssptr,1) 800 xj(nax+jnum,2)=xs(lssptr,2) 801 xj(nax+jnum,3)=xs(lssptr,3) 802 pl(nax+jnum,1)=ps(isfr+isa,1,2) 803 pl(nax+jnum,2)=ps(isfr+isa,2,2) 804 pl(nax+jnum,3)=ps(isfr+isa,3,2) 805 pj(nax+jnum,1)=ps(lssptr,1,2) 806 pj(nax+jnum,2)=ps(lssptr,2,2) 807 pj(nax+jnum,3)=ps(lssptr,3,2) 808 isal(nax+jnum)=isfr+isa 809 jsal(nax+jnum)=lssptr 810 7 continue 811 else 812 do 8 jnum=1,lssin(isa,ipss) 813 rwc(nax+jnum,1)=rwc(nax+jnum,1)-rwx(jnum,1) 814 rwc(nax+jnum,2)=rwc(nax+jnum,2)-rwx(jnum,2) 815 rwc(nax+jnum,3)=rwc(nax+jnum,3)-rwx(jnum,3) 816 lssptr=lssj(jsaptr+jnum) 817 xi(nax+jnum,1)=xs(isfr+isa,1) 818 xi(nax+jnum,2)=xs(isfr+isa,2) 819 xi(nax+jnum,3)=xs(isfr+isa,3) 820 xj(nax+jnum,1)=xs(lssptr,1)+rwx(jnum,1) 821 xj(nax+jnum,2)=xs(lssptr,2)+rwx(jnum,2) 822 xj(nax+jnum,3)=xs(lssptr,3)+rwx(jnum,3) 823 pl(nax+jnum,1)=ps(isfr+isa,1,2) 824 pl(nax+jnum,2)=ps(isfr+isa,2,2) 825 pl(nax+jnum,3)=ps(isfr+isa,3,2) 826 pj(nax+jnum,1)=ps(lssptr,1,2) 827 pj(nax+jnum,2)=ps(lssptr,2,2) 828 pj(nax+jnum,3)=ps(lssptr,3,2) 829 isal(nax+jnum)=isfr+isa 830 jsal(nax+jnum)=lssptr 831 8 continue 832 endif 833 nax=nax+lssin(isa,ipss) 834 4 continue 835c else 836c do 8 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss) 837c jsaptr=lssjpt(isa,ipss)-1 838c ism=isl(isfr+isa,lsmol) 839c do 9 jnum=1,lssin(isa,ipss) 840c lssptr=lssj(jsaptr+jnum) 841c jsa=lssptr 842c rwx(jnum,1)=zero 843c rwx(jnum,2)=zero 844c rwx(jnum,3)=zero 845c if(abs(xs(isfr+isa,1)-xs(jsa,1)).gt.boxh(1)) then 846c rwx(jnum,1)=sign(box(1),xs(isfr+isa,1)) 847c endif 848c if(abs(xs(isfr+isa,2)-xs(jsa,2)).gt.boxh(2)) then 849c rwx(jnum,2)=sign(box(2),xs(isfr+isa,2)) 850c endif 851c if(npbtyp.eq.1) then 852c if(abs(xs(isfr+isa,3)-xs(jsa,3)).gt.boxh(3)) then 853c rwx(jnum,3)=sign(box(3),xs(isfr+isa,3)) 854c endif 855c endif 856c jsf=isl(lssptr,lsfrc) 857c isal(nax+jnum)=isfr+isa 858c jsal(nax+jnum)=lssptr 859c jmal(nax+jnum)=jsf 860c jsm=isl(lssptr,lsmol) 861c jfal(nax+jnum)=2 862c if(ism.ne.jsm) jfal(nax+jnum)=3 863c if(ism.gt.0) then 864c if(jsm.gt.0) then 865c rwc(nax+jnum,1)=xsm(ism,1)-xsm(jsm,1)-rwx(jnum,1) 866c rwc(nax+jnum,2)=xsm(ism,2)-xsm(jsm,2)-rwx(jnum,2) 867c rwc(nax+jnum,3)=xsm(ism,3)-xsm(jsm,3)-rwx(jnum,3) 868c else 869c rwc(nax+jnum,1)=xsm(ism,1)-xs(lssptr,1)-rwx(jnum,1) 870c rwc(nax+jnum,2)=xsm(ism,2)-xs(lssptr,2)-rwx(jnum,2) 871c rwc(nax+jnum,3)=xsm(ism,3)-xs(lssptr,3)-rwx(jnum,3) 872c endif 873c else 874c if(jsm.gt.0) then 875c rwc(nax+jnum,1)=xs(isfr+isa,1)-xsm(jsm,1)-rwx(jnum,1) 876c rwc(nax+jnum,2)=xs(isfr+isa,2)-xsm(jsm,2)-rwx(jnum,2) 877c rwc(nax+jnum,3)=xs(isfr+isa,3)-xsm(jsm,3)-rwx(jnum,3) 878c else 879c rwc(nax+jnum,1)=xs(isfr+isa,1)-xs(lssptr,1)-rwx(jnum,1) 880c rwc(nax+jnum,2)=xs(isfr+isa,2)-xs(lssptr,2)-rwx(jnum,2) 881c rwc(nax+jnum,3)=xs(isfr+isa,3)-xs(lssptr,3)-rwx(jnum,3) 882c endif 883c endif 884c 9 continue 885c do 11 jnum=1,lssin(isa,ipss) 886c lssptr=lssj(jsaptr+jnum) 887c xi(nax+jnum,1)=xs(isfr+isa,1) 888c xi(nax+jnum,2)=xs(isfr+isa,2) 889c xi(nax+jnum,3)=xs(isfr+isa,3) 890c xj(nax+jnum,1)=xs(lssptr,1)+rwx(jnum,1) 891c xj(nax+jnum,2)=xs(lssptr,2)+rwx(jnum,2) 892c xj(nax+jnum,3)=xs(lssptr,3)+rwx(jnum,3) 893c pl(nax+jnum,1)=ps(isfr+isa,1,2) 894c pl(nax+jnum,2)=ps(isfr+isa,2,2) 895c pl(nax+jnum,3)=ps(isfr+isa,3,2) 896c pj(nax+jnum,1)=ps(lssptr,1,2) 897c pj(nax+jnum,2)=ps(lssptr,2,2) 898c pj(nax+jnum,3)=ps(lssptr,3,2) 899c 11 continue 900c nax=nax+lssin(isa,ipss) 901c 8 continue 902c endif 903c 904 do 12 iax=1,nax 905 isa=isal(iax) 906 jsa=jsal(iax) 907 if(jfal(iax).eq.2) then 908 qai=qfaci*chg(isq2(isa),1,iset) 909 pai=chg(isq2(isa),2,iset) 910 qaj=qfaci*chg(isq2(jsa),1,iset) 911 paj=chg(isq2(jsa),2,iset) 912 else 913 qai=qfaci*chg(isq3(isa),1,iset) 914 pai=chg(isq3(isa),2,iset) 915 qaj=qfaci*chg(isq3(jsa),1,iset) 916 paj=chg(isq3(jsa),2,iset) 917 endif 918 rx=xj(iax,1)-xi(iax,1) 919 ry=xj(iax,2)-xi(iax,2) 920 rz=xj(iax,3)-xi(iax,3) 921 ri2=one/(rx*rx+ry*ry+rz*rz) 922 ri1=sqrt(ri2) 923 ri3=ri1*ri2 924 pix=pai*pl(iax,1) 925 piy=pai*pl(iax,2) 926 piz=pai*pl(iax,3) 927 pjx=paj*pj(iax,1) 928 pjy=paj*pj(iax,2) 929 pjz=paj*pj(iax,3) 930 rmi=three*(rx*pix+ry*piy+rz*piz)*ri2 931 rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2 932 fi(iax,1)=((rmj-qaj)*rx-pjx)*ri3 933 fi(iax,2)=((rmj-qaj)*ry-pjy)*ri3 934 fi(iax,3)=((rmj-qaj)*rz-pjz)*ri3 935 fj(iax,1)=((rmi+qai)*rx-pix)*ri3 936 fj(iax,2)=((rmi+qai)*ry-piy)*ri3 937 fj(iax,3)=((rmi+qai)*rz-piz)*ri3 938 12 continue 939c 940c accumulate fields into solute field arrays 941c 942 iax=0 943 do 13 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss) 944 jsaptr=lssjpt(isa,ipss)-1 945 do 14 jnum=1,lssin(isa,ipss) 946 lssptr=lssj(jsaptr+jnum) 947 ps(isfr+isa,1,1)=ps(isfr+isa,1,1)+fi(iax+jnum,1) 948 ps(isfr+isa,2,1)=ps(isfr+isa,2,1)+fi(iax+jnum,2) 949 ps(isfr+isa,3,1)=ps(isfr+isa,3,1)+fi(iax+jnum,3) 950 ps(lssptr,1,1)=ps(lssptr,1,1)+fj(iax+jnum,1) 951 ps(lssptr,2,1)=ps(lssptr,2,1)+fj(iax+jnum,2) 952 ps(lssptr,3,1)=ps(lssptr,3,1)+fj(iax+jnum,3) 953 14 continue 954 iax=iax+lssin(isa,ipss) 955 13 continue 956c 957c thermodynamic integration and perturbation 958c 959 do 15 ipset=2,3 960 if((ipset.eq.2.and.ipert2).or. 961 + (ipset.eq.3.and.ipert3)) then 962c 963 nax2=0 964 do 16 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss) 965 jsaptr=lssjpt(isa,ipss)-1 966 do 17 jnum=1,lssin(isa,ipss) 967 lssptr=lssj(jsaptr+jnum) 968 pl(nax2+jnum,1)=psp(isfr+isa,1,ipset-1,2) 969 pl(nax2+jnum,2)=psp(isfr+isa,2,ipset-1,2) 970 pl(nax2+jnum,3)=psp(isfr+isa,3,ipset-1,2) 971 pj(nax2+jnum,1)=psp(lssptr,1,ipset-1,2) 972 pj(nax2+jnum,2)=psp(lssptr,2,ipset-1,2) 973 pj(nax2+jnum,3)=psp(lssptr,3,ipset-1,2) 974 17 continue 975 nax2=nax2+lssin(isa,ipss) 976 16 continue 977c 978 if(nax2.ne.nax) call md_abort('Error in dips',me) 979c 980 do 18 iax=1,nax 981 isa=isal(iax) 982 jsa=jsal(iax) 983 if(jfal(iax).eq.2) then 984 qai=qfaci*chg(isq2(isa),1,ipset) 985 pai=chg(isq2(isa),2,ipset) 986 qaj=qfaci*chg(isq2(jsa),1,ipset) 987 paj=chg(isq2(jsa),2,ipset) 988 else 989 qai=qfaci*chg(isq3(isa),1,ipset) 990 pai=chg(isq3(isa),2,ipset) 991 qaj=qfaci*chg(isq3(jsa),1,ipset) 992 paj=chg(isq3(jsa),2,ipset) 993 endif 994 rx=xj(iax,1)-xi(iax,1) 995 ry=xj(iax,2)-xi(iax,2) 996 rz=xj(iax,3)-xi(iax,3) 997 ri2=one/(rx*rx+ry*ry+rz*rz) 998 ri1=sqrt(ri2) 999 ri3=ri1*ri2 1000 pix=pai*pl(iax,1) 1001 piy=pai*pl(iax,2) 1002 piz=pai*pl(iax,3) 1003 pjx=paj*pj(iax,1) 1004 pjy=paj*pj(iax,2) 1005 pjz=paj*pj(iax,3) 1006 rmi=three*(rx*pix+ry*piy+rz*piz)*ri2 1007 rmj=three*(rx*pjx+ry*pjy+rz*pjz)*ri2 1008 fi(iax,1)=((rmj-qaj)*rx-pjx)*ri3 1009 fi(iax,2)=((rmj-qaj)*ry-pjy)*ri3 1010 fi(iax,3)=((rmj-qaj)*rz-pjz)*ri3 1011 fj(iax,1)=((rmi+qai)*rx-pix)*ri3 1012 fj(iax,2)=((rmi+qai)*ry-piy)*ri3 1013 fj(iax,3)=((rmi+qai)*rz-piz)*ri3 1014 18 continue 1015c 1016c accumulate fields into solute field arrays 1017c 1018 iax=0 1019 do 19 isa=lssndx(isslen-1,ipss)+1,lssndx(isslen,ipss) 1020 jsaptr=lssjpt(isa,ipss)-1 1021 do 20 jnum=1,lssin(isa,ipss) 1022 lssptr=lssj(jsaptr+jnum) 1023 psp(isfr+isa,1,ipset-1,1)=psp(isfr+isa,1,ipset-1,1)+fi(iax+jnum,1) 1024 psp(isfr+isa,2,ipset-1,1)=psp(isfr+isa,2,ipset-1,1)+fi(iax+jnum,2) 1025 psp(isfr+isa,3,ipset-1,1)=psp(isfr+isa,3,ipset-1,1)+fi(iax+jnum,3) 1026 psp(lssptr,1,ipset-1,1)=psp(lssptr,1,ipset-1,1)+fj(iax+jnum,1) 1027 psp(lssptr,2,ipset-1,1)=psp(lssptr,2,ipset-1,1)+fj(iax+jnum,2) 1028 psp(lssptr,3,ipset-1,1)=psp(lssptr,3,ipset-1,1)+fj(iax+jnum,3) 1029 20 continue 1030 iax=iax+lssin(isa,ipss) 1031 19 continue 1032c 1033 endif 1034 15 continue 1035 3 continue 1036 1 continue 1037c 1038 return 1039 end 1040