1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1990 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################## 9c ## ## 10c ## subroutine edipole1 -- dipole-dipole energy & derivs ## 11c ## ## 12c ############################################################## 13c 14c 15c "edipole1" calculates the dipole-dipole interaction energy 16c and first derivatives with respect to Cartesian coordinates 17c 18c 19 subroutine edipole1 20 use atoms 21 use bound 22 use cell 23 use chgpot 24 use deriv 25 use dipole 26 use energi 27 use group 28 use shunt 29 use units 30 use usage 31 use virial 32 implicit none 33 integer i,j,k 34 integer i1,i2,k1,k2 35 real*8 xi,yi,zi 36 real*8 xk,yk,zk 37 real*8 xq,yq,zq 38 real*8 xr,yr,zr 39 real*8 xq1,yq1,zq1 40 real*8 xq2,yq2,zq2 41 real*8 f,fi,fik,fgrp 42 real*8 e,r2,ri2,rk2,rirkr3 43 real*8 doti,dotk,dotp 44 real*8 si1,si2,sk1,sk2 45 real*8 de,dedr,dedrirk 46 real*8 deddoti,deddotk,deddotp 47 real*8 termx,termy,termz 48 real*8 dedrirkri2,dedrirkrk2 49 real*8 termxi,termyi,termzi 50 real*8 termxk,termyk,termzk 51 real*8 dedxi1,dedyi1,dedzi1 52 real*8 dedxi2,dedyi2,dedzi2 53 real*8 dedxk1,dedyk1,dedzk1 54 real*8 dedxk2,dedyk2,dedzk2 55 real*8 r,r3,r4,r5,taper,dtaper 56 real*8 dtaperx,dtapery,dtaperz 57 real*8 vxx,vyy,vzz 58 real*8 vyx,vzx,vzy 59 logical proceed 60 character*6 mode 61c 62c 63c zero out the overall dipole interaction energy and derivs, 64c then set up the constants for the calculation 65c 66 ed = 0.0d0 67 do i = 1, n 68 ded(1,i) = 0.0d0 69 ded(2,i) = 0.0d0 70 ded(3,i) = 0.0d0 71 end do 72 if (ndipole .eq. 0) return 73c 74c set conversion factor and switching function coefficients 75c 76 f = electric / (debye**2 * dielec) 77 mode = 'DIPOLE' 78 call switch (mode) 79c 80c compute the dipole interaction energy and first derivatives 81c 82 do i = 1, ndipole-1 83 i1 = idpl(1,i) 84 i2 = idpl(2,i) 85 si1 = 1.0d0 - sdpl(i) 86 si2 = sdpl(i) 87 xi = x(i2) - x(i1) 88 yi = y(i2) - y(i1) 89 zi = z(i2) - z(i1) 90 if (use_polymer) call imager (xi,yi,zi,-1) 91 ri2 = xi*xi + yi*yi + zi*zi 92 xq = x(i1) + xi*si2 93 yq = y(i1) + yi*si2 94 zq = z(i1) + zi*si2 95 fi = f * bdpl(i) 96c 97c decide whether to compute the current interaction 98c 99 do k = i+1, ndipole 100 k1 = idpl(1,k) 101 k2 = idpl(2,k) 102 proceed = .true. 103 if (use_group) call groups (proceed,fgrp,i1,i2,k1,k2,0,0) 104 if (proceed) proceed = (use(i1) .or. use(i2) .or. 105 & use(k1) .or. use(k2)) 106 if (proceed) proceed = (k1.ne.i1 .and. k1.ne.i2 .and. 107 & k2.ne.i1 .and. k2.ne.i2) 108c 109c compute the energy contribution for this interaction 110c 111 if (proceed) then 112 sk1 = 1.0d0 - sdpl(k) 113 sk2 = sdpl(k) 114 xk = x(k2) - x(k1) 115 yk = y(k2) - y(k1) 116 zk = z(k2) - z(k1) 117 if (use_polymer) call imager (xk,yk,zk,-1) 118 xr = xq - x(k1) - xk*sk2 119 yr = yq - y(k1) - yk*sk2 120 zr = zq - z(k1) - zk*sk2 121 call image (xr,yr,zr) 122 r2 = xr*xr + yr*yr + zr*zr 123 if (r2 .le. off2) then 124 rk2 = xk*xk + yk*yk + zk*zk 125 rirkr3 = sqrt(ri2*rk2*r2) * r2 126 dotp = xi*xk + yi*yk + zi*zk 127 doti = xi*xr + yi*yr + zi*zr 128 dotk = xk*xr + yk*yr + zk*zr 129 fik = fi * bdpl(k) 130c 131c form the energy and master chain rule term for derivatives 132c 133 e = fik * (dotp-3.0d0*doti*dotk/r2) / rirkr3 134 de = -fik / (rirkr3*r2) 135c 136c scale the interaction based on its group membership 137c 138 if (use_group) then 139 e = e * fgrp 140 de = de * fgrp 141 end if 142c 143c secondary chain rule terms for derivative expressions 144c 145 deddotp = -de * r2 146 deddoti = de * 3.0d0*dotk 147 deddotk = de * 3.0d0*doti 148 dedr = de * (3.0d0*dotp-15.0d0*doti*dotk/r2) 149 dedrirk = -e 150 dedrirkri2 = dedrirk / ri2 151 dedrirkrk2 = dedrirk / rk2 152c 153c more chain rule terms for derivative expressions 154c 155 termx = dedr*xr + deddoti*xi + deddotk*xk 156 termy = dedr*yr + deddoti*yi + deddotk*yk 157 termz = dedr*zr + deddoti*zi + deddotk*zk 158 termxi = dedrirkri2*xi + deddotp*xk + deddoti*xr 159 termyi = dedrirkri2*yi + deddotp*yk + deddoti*yr 160 termzi = dedrirkri2*zi + deddotp*zk + deddoti*zr 161 termxk = dedrirkrk2*xk + deddotp*xi + deddotk*xr 162 termyk = dedrirkrk2*yk + deddotp*yi + deddotk*yr 163 termzk = dedrirkrk2*zk + deddotp*zi + deddotk*zr 164c 165c finally, the individual first derivative components 166c 167 dedxi1 = si1*termx - termxi 168 dedyi1 = si1*termy - termyi 169 dedzi1 = si1*termz - termzi 170 dedxi2 = si2*termx + termxi 171 dedyi2 = si2*termy + termyi 172 dedzi2 = si2*termz + termzi 173 dedxk1 = -sk1*termx - termxk 174 dedyk1 = -sk1*termy - termyk 175 dedzk1 = -sk1*termz - termzk 176 dedxk2 = -sk2*termx + termxk 177 dedyk2 = -sk2*termy + termyk 178 dedzk2 = -sk2*termz + termzk 179c 180c use energy switching if near the cutoff distance 181c 182 if (r2 .gt. cut2) then 183 r = sqrt(r2) 184 r3 = r2 * r 185 r4 = r2 * r2 186 r5 = r2 * r3 187 taper = c5*r5 + c4*r4 + c3*r3 188 & + c2*r2 + c1*r + c0 189 dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 190 & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 191 dtaper = dtaper * e/r 192 dtaperx = xr * dtaper 193 dtapery = yr * dtaper 194 dtaperz = zr * dtaper 195 e = e * taper 196 dedxi1 = dedxi1*taper + si1*dtaperx 197 dedyi1 = dedyi1*taper + si1*dtapery 198 dedzi1 = dedzi1*taper + si1*dtaperz 199 dedxi2 = dedxi2*taper + si2*dtaperx 200 dedyi2 = dedyi2*taper + si2*dtapery 201 dedzi2 = dedzi2*taper + si2*dtaperz 202 dedxk1 = dedxk1*taper - sk1*dtaperx 203 dedyk1 = dedyk1*taper - sk1*dtapery 204 dedzk1 = dedzk1*taper - sk1*dtaperz 205 dedxk2 = dedxk2*taper - sk2*dtaperx 206 dedyk2 = dedyk2*taper - sk2*dtapery 207 dedzk2 = dedzk2*taper - sk2*dtaperz 208 end if 209c 210c increment the overall energy and derivative expressions 211c 212 ed = ed + e 213 ded(1,i1) = ded(1,i1) + dedxi1 214 ded(2,i1) = ded(2,i1) + dedyi1 215 ded(3,i1) = ded(3,i1) + dedzi1 216 ded(1,i2) = ded(1,i2) + dedxi2 217 ded(2,i2) = ded(2,i2) + dedyi2 218 ded(3,i2) = ded(3,i2) + dedzi2 219 ded(1,k1) = ded(1,k1) + dedxk1 220 ded(2,k1) = ded(2,k1) + dedyk1 221 ded(3,k1) = ded(3,k1) + dedzk1 222 ded(1,k2) = ded(1,k2) + dedxk2 223 ded(2,k2) = ded(2,k2) + dedyk2 224 ded(3,k2) = ded(3,k2) + dedzk2 225c 226c increment the internal virial tensor components 227c 228 xq1 = x(k1) - xq 229 yq1 = y(k1) - yq 230 zq1 = z(k1) - zq 231 xq2 = x(k2) - xq 232 yq2 = y(k2) - yq 233 zq2 = z(k2) - zq 234 vxx = xq1*dedxk1 + xq2*dedxk2 235 vyx = 0.5d0 * (yq1*dedxk1 + yq2*dedxk2 236 & + xq1*dedyk1 + xq2*dedyk2) 237 vzx = 0.5d0 * (zq1*dedxk1 + zq2*dedxk2 238 & + xq1*dedzk1 + xq2*dedzk2) 239 vyy = yq1*dedyk1 + yq2*dedyk2 240 vzy = 0.5d0 * (zq1*dedyk1 + zq2*dedyk2 241 & + yq1*dedzk1 + yq2*dedzk2) 242 vzz = zq1*dedzk1 + zq2*dedzk2 243 vir(1,1) = vir(1,1) + vxx 244 vir(2,1) = vir(2,1) + vyx 245 vir(3,1) = vir(3,1) + vzx 246 vir(1,2) = vir(1,2) + vyx 247 vir(2,2) = vir(2,2) + vyy 248 vir(3,2) = vir(3,2) + vzy 249 vir(1,3) = vir(1,3) + vzx 250 vir(2,3) = vir(2,3) + vzy 251 vir(3,3) = vir(3,3) + vzz 252 end if 253 end if 254 end do 255 end do 256c 257c for periodic boundary conditions with large cutoffs 258c neighbors must be found by the replicates method 259c 260 if (.not. use_replica) return 261c 262c calculate interaction energy with other unit cells 263c 264 do i = 1, ndipole 265 i1 = idpl(1,i) 266 i2 = idpl(2,i) 267 si1 = 1.0d0 - sdpl(i) 268 si2 = sdpl(i) 269 xi = x(i2) - x(i1) 270 yi = y(i2) - y(i1) 271 zi = z(i2) - z(i1) 272 if (use_polymer) call imager (xi,yi,zi,-1) 273 ri2 = xi*xi + yi*yi + zi*zi 274 xq = x(i1) + xi*si2 275 yq = y(i1) + yi*si2 276 zq = z(i1) + zi*si2 277 fi = f * bdpl(i) 278c 279c decide whether to compute the current interaction 280c 281 do k = i, ndipole 282 k1 = idpl(1,k) 283 k2 = idpl(2,k) 284 proceed = .true. 285 if (use_group) call groups (proceed,fgrp,i1,i2,k1,k2,0,0) 286 if (proceed) proceed = (use(i1) .or. use(i2) .or. 287 & use(k1) .or. use(k2)) 288c 289c compute the energy contribution for this interaction 290c 291 if (proceed) then 292 sk1 = 1.0d0 - sdpl(k) 293 sk2 = sdpl(k) 294 do j = 2, ncell 295 xk = x(k2) - x(k1) 296 yk = y(k2) - y(k1) 297 zk = z(k2) - z(k1) 298 if (use_polymer) call imager (xk,yk,zk,-1) 299 xr = xq - x(k1) - xk*sk2 300 yr = yq - y(k1) - yk*sk2 301 zr = zq - z(k1) - zk*sk2 302 call imager (xr,yr,zr,j) 303 r2 = xr*xr + yr*yr + zr*zr 304 if (r2 .le. off2) then 305 rk2 = xk*xk + yk*yk + zk*zk 306 rirkr3 = sqrt(ri2*rk2*r2) * r2 307 dotp = xi*xk + yi*yk + zi*zk 308 doti = xi*xr + yi*yr + zi*zr 309 dotk = xk*xr + yk*yr + zk*zr 310 fik = fi * bdpl(k) 311 if (use_polymer) then 312 if (r2 .lt. polycut2) then 313 if (k1.eq.i1 .or. k1.eq.i2 .or. 314 & k2.eq.i1 .or. k2.eq.i2) fik = 0.0d0 315 end if 316 end if 317c 318c form the energy and master chain rule term for derivatives 319c 320 e = fik * (dotp-3.0d0*doti*dotk/r2) / rirkr3 321 de = -fik / (rirkr3*r2) 322c 323c scale the interaction based on its group membership 324c 325 if (use_group) then 326 e = e * fgrp 327 de = de * fgrp 328 end if 329c 330c secondary chain rule terms for derivative expressions 331c 332 deddotp = -de * r2 333 deddoti = de * 3.0d0*dotk 334 deddotk = de * 3.0d0*doti 335 dedr = de * (3.0d0*dotp-15.0d0*doti*dotk/r2) 336 dedrirk = -e 337 dedrirkri2 = dedrirk / ri2 338 dedrirkrk2 = dedrirk / rk2 339c 340c more chain rule terms for derivative expressions 341c 342 termx = dedr*xr + deddoti*xi + deddotk*xk 343 termy = dedr*yr + deddoti*yi + deddotk*yk 344 termz = dedr*zr + deddoti*zi + deddotk*zk 345 termxi = dedrirkri2*xi + deddotp*xk + deddoti*xr 346 termyi = dedrirkri2*yi + deddotp*yk + deddoti*yr 347 termzi = dedrirkri2*zi + deddotp*zk + deddoti*zr 348 termxk = dedrirkrk2*xk + deddotp*xi + deddotk*xr 349 termyk = dedrirkrk2*yk + deddotp*yi + deddotk*yr 350 termzk = dedrirkrk2*zk + deddotp*zi + deddotk*zr 351c 352c finally, the individual first derivative components 353c 354 dedxi1 = si1*termx - termxi 355 dedyi1 = si1*termy - termyi 356 dedzi1 = si1*termz - termzi 357 dedxi2 = si2*termx + termxi 358 dedyi2 = si2*termy + termyi 359 dedzi2 = si2*termz + termzi 360 dedxk1 = -sk1*termx - termxk 361 dedyk1 = -sk1*termy - termyk 362 dedzk1 = -sk1*termz - termzk 363 dedxk2 = -sk2*termx + termxk 364 dedyk2 = -sk2*termy + termyk 365 dedzk2 = -sk2*termz + termzk 366c 367c use energy switching if near the cutoff distance 368c 369 if (r2 .gt. cut2) then 370 r = sqrt(r2) 371 r3 = r2 * r 372 r4 = r2 * r2 373 r5 = r2 * r3 374 taper = c5*r5 + c4*r4 + c3*r3 375 & + c2*r2 + c1*r + c0 376 dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 377 & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 378 dtaper = dtaper * e/r 379 dtaperx = xr * dtaper 380 dtapery = yr * dtaper 381 dtaperz = zr * dtaper 382 e = e * taper 383 dedxi1 = dedxi1*taper + si1*dtaperx 384 dedyi1 = dedyi1*taper + si1*dtapery 385 dedzi1 = dedzi1*taper + si1*dtaperz 386 dedxi2 = dedxi2*taper + si2*dtaperx 387 dedyi2 = dedyi2*taper + si2*dtapery 388 dedzi2 = dedzi2*taper + si2*dtaperz 389 dedxk1 = dedxk1*taper - sk1*dtaperx 390 dedyk1 = dedyk1*taper - sk1*dtapery 391 dedzk1 = dedzk1*taper - sk1*dtaperz 392 dedxk2 = dedxk2*taper - sk2*dtaperx 393 dedyk2 = dedyk2*taper - sk2*dtapery 394 dedzk2 = dedzk2*taper - sk2*dtaperz 395 end if 396c 397c increment the overall energy and derivative expressions 398c 399 if (i .eq. k) e = 0.5d0 * e 400 ed = ed + e 401 ded(1,i1) = ded(1,i1) + dedxi1 402 ded(2,i1) = ded(2,i1) + dedyi1 403 ded(3,i1) = ded(3,i1) + dedzi1 404 ded(1,i2) = ded(1,i2) + dedxi2 405 ded(2,i2) = ded(2,i2) + dedyi2 406 ded(3,i2) = ded(3,i2) + dedzi2 407 if (i .ne. k) then 408 ded(1,k1) = ded(1,k1) + dedxk1 409 ded(2,k1) = ded(2,k1) + dedyk1 410 ded(3,k1) = ded(3,k1) + dedzk1 411 ded(1,k2) = ded(1,k2) + dedxk2 412 ded(2,k2) = ded(2,k2) + dedyk2 413 ded(3,k2) = ded(3,k2) + dedzk2 414 end if 415c 416c increment the internal virial tensor components 417c 418 xq1 = x(k1) - xq 419 yq1 = y(k1) - yq 420 zq1 = z(k1) - zq 421 xq2 = x(k2) - xq 422 yq2 = y(k2) - yq 423 zq2 = z(k2) - zq 424 vxx = xq1*dedxk1 + xq2*dedxk2 425 vyx = 0.5d0 * (yq1*dedxk1 + yq2*dedxk2 426 & + xq1*dedyk1 + xq2*dedyk2) 427 vzx = 0.5d0 * (zq1*dedxk1 + zq2*dedxk2 428 & + xq1*dedzk1 + xq2*dedzk2) 429 vyy = yq1*dedyk1 + yq2*dedyk2 430 vzy = 0.5d0 * (zq1*dedyk1 + zq2*dedyk2 431 & + yq1*dedzk1 + yq2*dedzk2) 432 vzz = zq1*dedzk1 + zq2*dedzk2 433 vir(1,1) = vir(1,1) + vxx 434 vir(2,1) = vir(2,1) + vyx 435 vir(3,1) = vir(3,1) + vzx 436 vir(1,2) = vir(1,2) + vyx 437 vir(2,2) = vir(2,2) + vyy 438 vir(3,2) = vir(3,2) + vzy 439 vir(1,3) = vir(1,3) + vzx 440 vir(2,3) = vir(2,3) + vzy 441 vir(3,3) = vir(3,3) + vzz 442 end if 443 end do 444 end if 445 end do 446 end do 447 return 448 end 449