1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1995 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ######################################################### 9c ## ## 10c ## subroutine epolar -- dipole polarization energy ## 11c ## ## 12c ######################################################### 13c 14c 15c "epolar" calculates the induced dipole polarization energy 16c 17c 18 subroutine epolar 19 implicit none 20 include 'sizes.i' 21 include 'atoms.i' 22 include 'bound.i' 23 include 'couple.i' 24 include 'energi.i' 25 include 'mpole.i' 26 include 'polar.i' 27 include 'shunt.i' 28 include 'usage.i' 29 integer i,j,k,skip(maxatm) 30 integer ii,kk,iz,ix 31 real*8 e,epolik,a(3,3) 32 real*8 xi,yi,zi,xr,yr,zr 33 real*8 r,r2,r3,r4,r5,taper 34 real*8 rpi(13),indk(3) 35 logical iuse,kuse 36c 37c 38c zero out the induced dipole polarization energy 39c 40 ep = 0.0d0 41c 42c zero out the list of atoms to be skipped 43c 44 do i = 1, n 45 skip(i) = 0 46 end do 47c 48c set the switching function coefficients 49c 50 call switch ('CHGDPL') 51c 52c rotate the multipole components into the global frame 53c 54 call rotpole 55c 56c compute the induced dipoles at each atom 57c 58 call induce 59c 60c calculate the induced dipole polarization energy term 61c 62 do ii = 1, npole 63 i = ipole(ii) 64 iz = zaxis(ii) 65 ix = xaxis(ii) 66 iuse = (use(i) .or. use(iz) .or. use(ix)) 67 xi = x(i) 68 yi = y(i) 69 zi = z(i) 70 skip(i) = i 71 do j = 1, n12(i) 72 skip(i12(j,i)) = i 73 end do 74 do j = 1, n13(i) 75 skip(i13(j,i)) = i 76 end do 77 do j = 1, mdqsiz 78 rpi(j) = rpole(j,ii) 79 end do 80 do kk = 1, npolar 81 k = ipolar(kk) 82 kuse = use(k) 83 if (iuse .or. kuse) then 84 if (skip(k) .ne. i) then 85 xr = x(k) - xi 86 yr = y(k) - yi 87 zr = z(k) - zi 88 if (use_image) call image (xr,yr,zr,0) 89 r2 = xr*xr + yr*yr + zr*zr 90 if (r2 .le. off2) then 91 indk(1) = uind(1,k) 92 indk(2) = uind(2,k) 93 indk(3) = uind(3,k) 94 r = sqrt(r2) 95 e = epolik (r,xr,yr,zr,rpi,indk) 96 if (r2 .gt. cut2) then 97 r3 = r2 * r 98 r4 = r2 * r2 99 r5 = r2 * r3 100 taper = c5*r5 + c4*r4 + c3*r3 101 & + c2*r2 + c1*r + c0 102 e = e * taper 103 end if 104 ep = ep + e 105 end if 106 end if 107 end if 108 end do 109 end do 110 return 111 end 112c 113c 114c ####################### 115c ## ## 116c ## function epolik ## 117c ## ## 118c ####################### 119c 120c 121 function epolik (r,xr,yr,zr,rpi,indk) 122 implicit none 123 include 'sizes.i' 124 include 'atoms.i' 125 include 'electr.i' 126 include 'mpole.i' 127 include 'units.i' 128 integer i,j 129 real*8 epolik,t2matrix 130 real*8 r,r2,r3,r4 131 real*8 xr,yr,zr,zrr 132 real*8 xrr,xrr2,xrr3 133 real*8 yrr,yrr2,yrr3 134 real*8 rpi(13),indk(3) 135 real*8 m2t2(13),t2(13,13) 136 real*8 bx(0:5,0:5),by(11,0:5,0:1),bz(11,0:1) 137c 138c 139c zeroth order coefficients to speed the T2 calculation 140c 141 if (mdqsiz .ge. 1) then 142 bz(1,0) = c(1,0,0) 143 by(1,0,0) = c(1,0,0) * bz(1,0) 144 bx(0,0) = c(1,0,0) 145c 146c first order coefficients to speed the T2 calculation 147c 148 r2 = r * r 149 zrr = zr / r 150 bz(3,0) = c(3,0,0) 151 bz(1,1) = c(1,1,1) * zrr 152 yrr = yr / r 153 by(3,0,0) = c(3,0,0) * bz(3,0) 154 by(1,0,1) = c(1,0,0) * bz(1,1) 155 by(1,1,0) = c(1,1,1) * bz(3,0) * yrr 156 xrr = xr / r 157 bx(1,1) = c(1,1,1) * xrr 158 end if 159c 160c second order coefficients to speed the T2 calculation 161c 162 if (mdqsiz .ge. 4) then 163 r3 = r2 * r 164 bz(5,0) = c(5,0,0) 165 bz(3,1) = c(3,1,1) * zrr 166 yrr2 = yrr * yrr 167 by(5,0,0) = c(5,0,0) * bz(5,0) 168 by(3,0,1) = c(3,0,0) * bz(3,1) 169 by(3,1,0) = c(3,1,1) * bz(5,0) * yrr 170 by(1,1,1) = c(1,1,1) * bz(3,1) * yrr 171 by(1,2,0) = c(1,2,0)*bz(3,0) + c(1,2,2)*bz(5,0)*yrr2 172 xrr2 = xrr * xrr 173 bx(2,0) = c(1,2,0) 174 bx(2,2) = c(1,2,2) * xrr2 175 end if 176c 177c third order coefficients to speed the T2 calculation 178c 179 if (mdqsiz .ge. 13) then 180 r4 = r2 * r2 181 bz(7,0) = c(7,0,0) 182 bz(5,1) = c(5,1,1) * zrr 183 yrr3 = yrr2 * yrr 184 by(7,0,0) = c(7,0,0)*bz(7,0) 185 by(5,0,1) = c(5,0,0)*bz(5,1) 186 by(5,1,0) = c(5,1,1)*bz(7,0)*yrr 187 by(3,1,1) = c(3,1,1)*bz(5,1)*yrr 188 by(3,2,0) = c(3,2,0)*bz(5,0) + c(3,2,2)*bz(7,0)*yrr2 189 by(1,2,1) = c(1,2,0)*bz(3,1) + c(1,2,2)*bz(5,1)*yrr2 190 by(1,3,0) = c(1,3,1)*bz(5,0)*yrr + c(1,3,3)*bz(7,0)*yrr3 191 xrr3 = xrr2 * xrr 192 bx(3,1) = c(1,3,1) * xrr 193 bx(3,3) = c(1,3,3) * xrr3 194 end if 195c 196c first order T2 matrix elements 197c 198 if (mdqsiz .ge. 1) then 199 t2(2,1) = t2matrix (0,0,1,r2,bx,by) 200 t2(3,1) = t2matrix (0,1,0,r2,bx,by) 201 t2(4,1) = t2matrix (1,0,0,r2,bx,by) 202 end if 203c 204c second order T2 matrix elements 205c 206 if (mdqsiz .ge. 4) then 207 t2(2,2) = -t2matrix (0,0,2,r3,bx,by) 208 t2(3,2) = -t2matrix (0,1,1,r3,bx,by) 209 t2(4,2) = -t2matrix (1,0,1,r3,bx,by) 210 t2(3,3) = -t2matrix (0,2,0,r3,bx,by) 211 t2(4,3) = -t2matrix (1,1,0,r3,bx,by) 212 t2(4,4) = -t2(2,2) - t2(3,3) 213 t2(2,3) = t2(3,2) 214 t2(2,4) = t2(4,2) 215 t2(3,4) = t2(4,3) 216 end if 217c 218c third order T2 matrix elements 219c 220 if (mdqsiz .ge. 13) then 221 t2(2,5) = t2matrix (0,0,3,r4,bx,by) 222 t2(3,5) = t2matrix (0,1,2,r4,bx,by) 223 t2(4,5) = t2matrix (1,0,2,r4,bx,by) 224 t2(3,6) = t2matrix (0,2,1,r4,bx,by) 225 t2(4,6) = t2matrix (1,1,1,r4,bx,by) 226 t2(4,7) = -t2(2,5) - t2(3,6) 227 t2(2,6) = t2(3,5) 228 t2(2,7) = t2(4,5) 229 t2(3,7) = t2(4,6) 230 t2(2,8) = t2(2,6) 231 t2(2,9) = t2(3,6) 232 t2(2,10) = t2(4,6) 233 t2(3,9) = t2matrix (0,3,0,r4,bx,by) 234 t2(3,10) = t2matrix (1,2,0,r4,bx,by) 235 t2(4,10) = -t2(2,8) - t2(3,9) 236 t2(3,8) = t2(2,9) 237 t2(4,8) = t2(2,10) 238 t2(4,9) = t2(3,10) 239 t2(2,11) = t2(2,7) 240 t2(2,12) = t2(3,7) 241 t2(2,13) = t2(4,7) 242 t2(3,12) = t2(3,10) 243 t2(3,13) = t2(4,10) 244 t2(4,13) = -t2(2,11) - t2(3,12) 245 t2(3,11) = t2(2,12) 246 t2(4,11) = t2(2,13) 247 t2(4,12) = t2(3,13) 248 end if 249c 250c compute energy of the induced dipole and multipole sites 251c 252 do i = 1, mdqsiz 253 m2t2(i) = 0.0d0 254 do j = 2, 4 255 m2t2(i) = m2t2(i) + indk(j-1)*t2(j,i) 256 end do 257 end do 258 epolik = 0.0d0 259 do i = 1, mdqsiz 260 epolik = epolik + m2t2(i)*rpi(i) 261 end do 262 epolik = 0.5d0 * (electric / dielec) * epolik 263 return 264 end 265c 266c 267c ################################################### 268c ## COPYRIGHT (C) 1995 by Jay William Ponder ## 269c ################################################### 270c 271c ################################################################# 272c ## ## 273c ## subroutine epolar1 -- polarization energy & derivatives ## 274c ## ## 275c ################################################################# 276c 277c 278c "epolar1" calculates the induced dipole polarization energy 279c and first derivatives with respect to Cartesian coordinates 280c 281c 282 subroutine epolar1 283 implicit none 284 include 'sizes.i' 285 include 'atoms.i' 286 include 'bath.i' 287 include 'bound.i' 288 include 'couple.i' 289 include 'deriv.i' 290 include 'energi.i' 291 include 'inter.i' 292 include 'molcul.i' 293 include 'mpole.i' 294 include 'polar.i' 295 include 'polpot.i' 296 include 'shunt.i' 297 include 'usage.i' 298 include 'virial.i' 299 integer i,j,k,skip(maxatm) 300 integer ii,kk,iz,ix 301 real*8 e,de,epolik1,dutu 302 real*8 xi,yi,zi,xr,yr,zr 303 real*8 xiz,yiz,ziz,xix,yix,zix 304 real*8 taper,dtaper 305 real*8 r,r2,r3,r4,r5 306 real*8 a(3,3),d(3,3,3,3) 307 real*8 rpi(13),indi(3),indk(3) 308 real*8 d1(3),d2(3,3),d3(3),utu 309 logical iuse,kuse 310c 311c 312c zero out the induced dipole energy and derivatives 313c 314 ep = 0.0d0 315 do i = 1, n 316 do j = 1, 3 317 dep(j,i) = 0.0d0 318 end do 319 end do 320c 321c zero out the list of atoms to be skipped 322c 323 do i = 1, n 324 skip(i) = 0 325 end do 326c 327c set the switching function coefficients 328c 329 call switch ('CHGDPL') 330c 331c rotate multipole components and get rotation derivatives 332c 333 do i = 1, npole 334 call rotmat (i,a) 335 call rotsite (i,a) 336 do j = 1, 3 337 do k = 1, 3 338 call drotmat (i,d) 339 call drotpole (i,a,d,j,k) 340 end do 341 end do 342 end do 343c 344c compute the induced dipoles at each atom 345c 346 call induce 347c 348c calculate the induced dipole energy and derivatives 349c 350 do ii = 1, npole 351 i = ipole(ii) 352 iz = zaxis(ii) 353 ix = xaxis(ii) 354 iuse = (use(i) .or. use(iz) .or. use(ix)) 355 xi = x(i) 356 yi = y(i) 357 zi = z(i) 358 skip(i) = i 359 do j = 1, n12(i) 360 skip(i12(j,i)) = i 361 end do 362 do j = 1, n13(i) 363 skip(i13(j,i)) = i 364 end do 365 do j = 1, mdqsiz 366 rpi(j) = rpole(j,ii) 367 end do 368 indi(1) = uind(1,i) 369 indi(2) = uind(2,i) 370 indi(3) = uind(3,i) 371 do kk = 1, npolar 372 k = ipolar(kk) 373 kuse = use(k) 374 if (iuse .or. kuse) then 375 if (skip(k) .ne. i) then 376 xr = x(k) - xi 377 yr = y(k) - yi 378 zr = z(k) - zi 379 if (use_image) call image (xr,yr,zr,0) 380 r2 = xr*xr + yr*yr + zr*zr 381 if (r2 .le. off2) then 382 indk(1) = uind(1,k) 383 indk(2) = uind(2,k) 384 indk(3) = uind(3,k) 385 r = sqrt(r2) 386 e = epolik1 (ii,r,xr,yr,zr,rpi,indi,indk, 387 & d1,d2,d3,utu) 388 if (r2 .gt. cut2) then 389 r3 = r2 * r 390 r4 = r2 * r2 391 r5 = r2 * r3 392 taper = c5*r5 + c4*r4 + c3*r3 393 & + c2*r2 + c1*r + c0 394 dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 395 & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 396 de = 2.0d0 * e * dtaper 397 dutu = utu * dtaper 398 e = e * taper 399 d1(1) = d1(1)*taper + de*(xr/r) 400 d1(2) = d1(2)*taper + de*(yr/r) 401 d1(3) = d1(3)*taper + de*(zr/r) 402 d2(1,1) = d2(1,1) * taper 403 d2(1,2) = d2(1,2) * taper 404 d2(1,3) = d2(1,3) * taper 405 d2(2,1) = d2(2,1) * taper 406 d2(2,2) = d2(2,2) * taper 407 d2(2,3) = d2(2,3) * taper 408 d2(3,1) = d2(3,1) * taper 409 d2(3,2) = d2(3,2) * taper 410 d2(3,3) = d2(3,3) * taper 411 d3(1) = d3(1)*taper + dutu*(xr/r) 412 d3(2) = d3(2)*taper + dutu*(yr/r) 413 d3(3) = d3(3)*taper + dutu*(zr/r) 414 end if 415 ep = ep + e 416 dep(1,i) = dep(1,i) - d1(1) - d2(1,1) 417 dep(2,i) = dep(2,i) - d1(2) - d2(1,2) 418 dep(3,i) = dep(3,i) - d1(3) - d2(1,3) 419 if (poltyp .ne. 'DIRECT') then 420 dep(1,i) = dep(1,i) - d3(1) 421 dep(2,i) = dep(2,i) - d3(2) 422 dep(3,i) = dep(3,i) - d3(3) 423 end if 424 dep(1,iz) = dep(1,iz) - d2(2,1) 425 dep(2,iz) = dep(2,iz) - d2(2,2) 426 dep(3,iz) = dep(3,iz) - d2(2,3) 427 dep(1,ix) = dep(1,ix) - d2(3,1) 428 dep(2,ix) = dep(2,ix) - d2(3,2) 429 dep(3,ix) = dep(3,ix) - d2(3,3) 430 dep(1,k) = dep(1,k) + d1(1) 431 dep(2,k) = dep(2,k) + d1(2) 432 dep(3,k) = dep(3,k) + d1(3) 433c 434c increment the total intermolecular energy 435c 436 if (molcule(i) .ne. molcule(k)) then 437 einter = einter + e 438 end if 439c 440c increment the virial for use in pressure computation 441c 442 if (isobaric) then 443 xiz = x(i) - x(iz) 444 yiz = y(i) - y(iz) 445 ziz = z(i) - z(iz) 446 xix = x(i) - x(ix) 447 yix = y(i) - y(ix) 448 zix = z(i) - z(ix) 449 virx = virx + xr*d1(1) + xiz*d2(2,1) 450 & + xix*d2(3,1) 451 viry = viry + yr*d1(2) + yiz*d2(2,2) 452 & + yix*d2(3,2) 453 virz = virz + zr*d1(3) + ziz*d2(2,3) 454 & + zix*d2(3,3) 455 if (poltyp .ne. 'DIRECT') then 456 virx = virx + 0.5d0*xr*d3(1) 457 viry = viry + 0.5d0*yr*d3(2) 458 virz = virz + 0.5d0*zr*d3(3) 459 end if 460 end if 461 end if 462 end if 463 end if 464 end do 465 end do 466 return 467 end 468c 469c 470c ######################## 471c ## ## 472c ## function epolik1 ## 473c ## ## 474c ######################## 475c 476c 477 function epolik1 (ii,r,xr,yr,zr,rpi,indi,indk,d1,d2,d3,utu) 478 implicit none 479 include 'sizes.i' 480 include 'atoms.i' 481 include 'electr.i' 482 include 'mpole.i' 483 include 'units.i' 484 integer i,j,k,m,ii 485 real*8 epolik1,t2matrix 486 real*8 r,r2,r3,r4,r5 487 real*8 xr,yr,zr,zrr,factor 488 real*8 xrr,xrr2,xrr3,xrr4 489 real*8 yrr,yrr2,yrr3,yrr4 490 real*8 t3x(3,13),t3y(3,13),t3z(3,13),tt2(3,13) 491 real*8 rpi(13),indi(3),indk(3) 492 real*8 m2t2(13),t2(13,13) 493 real*8 interx(3),intery(3),interz(3),m1t(13) 494 real*8 bx(0:5,0:5),by(11,0:5,0:1),bz(11,0:1) 495 real*8 d1(3),d2(3,3),d3(3),utu 496c 497c 498c zeroth order coefficients to speed the T2 calculation 499c 500 if (mdqsiz .ge. 1) then 501 bz(1,0) = c(1,0,0) 502 by(1,0,0) = c(1,0,0) * bz(1,0) 503 bx(0,0) = c(1,0,0) 504c 505c first order coefficients to speed the T2 calculation 506c 507 r2 = r * r 508 zrr = zr / r 509 bz(3,0) = c(3,0,0) 510 bz(1,1) = c(1,1,1) * zrr 511 yrr = yr / r 512 by(3,0,0) = c(3,0,0) * bz(3,0) 513 by(1,0,1) = c(1,0,0) * bz(1,1) 514 by(1,1,0) = c(1,1,1) * bz(3,0) * yrr 515 xrr = xr / r 516 bx(1,1) = c(1,1,1) * xrr 517c 518c second order coefficients to speed the T2 calculation 519c 520 r3 = r2 * r 521 bz(5,0) = c(5,0,0) 522 bz(3,1) = c(3,1,1) * zrr 523 yrr2 = yrr * yrr 524 by(5,0,0) = c(5,0,0) * bz(5,0) 525 by(3,0,1) = c(3,0,0) * bz(3,1) 526 by(3,1,0) = c(3,1,1) * bz(5,0) * yrr 527 by(1,1,1) = c(1,1,1) * bz(3,1) * yrr 528 by(1,2,0) = c(1,2,0)*bz(3,0) + c(1,2,2)*bz(5,0)*yrr2 529 xrr2 = xrr * xrr 530 bx(2,0) = c(1,2,0) 531 bx(2,2) = c(1,2,2) * xrr2 532 end if 533c 534c third order coefficients to speed the T2 calculation 535c 536 if (mdqsiz .ge. 4) then 537 r4 = r2 * r2 538 bz(7,0) = c(7,0,0) 539 bz(5,1) = c(5,1,1) * zrr 540 yrr3 = yrr2 * yrr 541 by(7,0,0) = c(7,0,0)*bz(7,0) 542 by(5,0,1) = c(5,0,0)*bz(5,1) 543 by(5,1,0) = c(5,1,1)*bz(7,0)*yrr 544 by(3,1,1) = c(3,1,1)*bz(5,1)*yrr 545 by(3,2,0) = c(3,2,0)*bz(5,0) + c(3,2,2)*bz(7,0)*yrr2 546 by(1,2,1) = c(1,2,0)*bz(3,1) + c(1,2,2)*bz(5,1)*yrr2 547 by(1,3,0) = c(1,3,1)*bz(5,0)*yrr + c(1,3,3)*bz(7,0)*yrr3 548 xrr3 = xrr2 * xrr 549 bx(3,1) = c(1,3,1) * xrr 550 bx(3,3) = c(1,3,3) * xrr3 551 end if 552c 553c fourth order coefficients to speed the T2 calculation 554c 555 if (mdqsiz .ge. 13) then 556 r5 = r2 * r3 557 bz(9,0) = c(9,0,0) 558 bz(7,1) = c(7,1,1) * zrr 559 yrr4 = yrr2 * yrr2 560 by(9,0,0) = c(9,0,0)*bz(9,0) 561 by(7,0,1) = c(7,0,0)*bz(7,1) 562 by(7,1,0) = c(7,1,1)*bz(9,0)*yrr 563 by(5,1,1) = c(5,1,1)*bz(7,1)*yrr 564 by(5,2,0) = c(5,2,0)*bz(7,0) + c(5,2,2)*bz(9,0)*yrr2 565 by(3,2,1) = c(3,2,0)*bz(5,1) + c(3,2,2)*bz(7,1)*yrr2 566 by(3,3,0) = c(3,3,1)*bz(7,0)*yrr + c(3,3,3)*bz(9,0)*yrr3 567 by(1,3,1) = c(1,3,1)*bz(5,1)*yrr + c(1,3,3)*bz(7,1)*yrr3 568 by(1,4,0) = c(1,4,0)*bz(5,0) + c(1,4,2)*bz(7,0)*yrr2 569 & + c(1,4,4)*bz(9,0)*yrr4 570 xrr4 = xrr2 * xrr2 571 bx(4,0) = c(1,4,0) 572 bx(4,2) = c(1,4,2) * xrr2 573 bx(4,4) = c(1,4,4) * xrr4 574 end if 575c 576c zeroth order T2 matrix elements 577c 578 if (mdqsiz .ge. 1) then 579 t2(1,1) = t2matrix (0,0,0,r,bx,by) 580c 581c first order T2 matrix elements 582c 583 t2(2,1) = t2matrix (0,0,1,r2,bx,by) 584 t2(3,1) = t2matrix (0,1,0,r2,bx,by) 585 t2(4,1) = t2matrix (1,0,0,r2,bx,by) 586 t2(1,2) = -t2(2,1) 587 t2(1,3) = -t2(3,1) 588 t2(1,4) = -t2(4,1) 589c 590c second order T2 matrix elements 591c 592 t2(2,2) = -t2matrix (0,0,2,r3,bx,by) 593 t2(3,2) = -t2matrix (0,1,1,r3,bx,by) 594 t2(4,2) = -t2matrix (1,0,1,r3,bx,by) 595 t2(3,3) = -t2matrix (0,2,0,r3,bx,by) 596 t2(4,3) = -t2matrix (1,1,0,r3,bx,by) 597 t2(4,4) = -t2(2,2) - t2(3,3) 598 t2(2,3) = t2(3,2) 599 t2(2,4) = t2(4,2) 600 t2(3,4) = t2(4,3) 601 k = 5 602 do i = 2, 4 603 do j = 2, 4 604 t2(1,k) = -t2(i,j) 605 t2(k,1) = t2(1,k) 606 k = k + 1 607 end do 608 end do 609 end if 610c 611c third order T2 matrix elements 612c 613 if (mdqsiz .ge. 4) then 614 t2(2,5) = t2matrix (0,0,3,r4,bx,by) 615 t2(3,5) = t2matrix (0,1,2,r4,bx,by) 616 t2(4,5) = t2matrix (1,0,2,r4,bx,by) 617 t2(3,6) = t2matrix (0,2,1,r4,bx,by) 618 t2(4,6) = t2matrix (1,1,1,r4,bx,by) 619 t2(4,7) = -t2(2,5) - t2(3,6) 620 t2(2,6) = t2(3,5) 621 t2(2,7) = t2(4,5) 622 t2(3,7) = t2(4,6) 623 t2(2,8) = t2(2,6) 624 t2(2,9) = t2(3,6) 625 t2(2,10) = t2(4,6) 626 t2(3,9) = t2matrix (0,3,0,r4,bx,by) 627 t2(3,10) = t2matrix (1,2,0,r4,bx,by) 628 t2(4,10) = -t2(2,8) - t2(3,9) 629 t2(3,8) = t2(2,9) 630 t2(4,8) = t2(2,10) 631 t2(4,9) = t2(3,10) 632 t2(2,11) = t2(2,7) 633 t2(2,12) = t2(3,7) 634 t2(2,13) = t2(4,7) 635 t2(3,12) = t2(3,10) 636 t2(3,13) = t2(4,10) 637 t2(4,13) = -t2(2,11) - t2(3,12) 638 t2(3,11) = t2(2,12) 639 t2(4,11) = t2(2,13) 640 t2(4,12) = t2(3,13) 641 do i = 5, 13 642 do j = 2, 4 643 t2(i,j) = -t2(j,i) 644 end do 645 end do 646 end if 647c 648c fourth order T2 matrix elements 649c 650 if (mdqsiz .ge.13) then 651 t2(5,5) = t2matrix (0,0,4,r5,bx,by) 652 t2(6,5) = t2matrix (0,1,3,r5,bx,by) 653 t2(7,5) = t2matrix (1,0,3,r5,bx,by) 654 t2(6,6) = t2matrix (0,2,2,r5,bx,by) 655 t2(7,6) = t2matrix (1,1,2,r5,bx,by) 656 t2(7,7) = -t2(5,5) - t2(6,6) 657 t2(8,5) = t2(6,5) 658 t2(9,5) = t2(6,6) 659 t2(10,5) = t2(7,6) 660 t2(9,6) = t2matrix (0,3,1,r5,bx,by) 661 t2(10,6) = t2matrix (1,2,1,r5,bx,by) 662 t2(10,7) = -t2(8,5) - t2(9,6) 663 t2(8,6) = t2(9,5) 664 t2(8,7) = t2(10,5) 665 t2(9,7) = t2(10,6) 666 t2(11,5) = t2(7,5) 667 t2(12,5) = t2(7,6) 668 t2(13,5) = t2(7,7) 669 t2(12,6) = t2(10,6) 670 t2(13,6) = t2(10,7) 671 t2(13,7) = -t2(11,5) - t2(12,6) 672 t2(11,6) = t2(12,5) 673 t2(11,7) = t2(13,5) 674 t2(12,7) = t2(13,6) 675 t2(8,8) = t2(8,6) 676 t2(9,8) = t2(9,6) 677 t2(10,8) = t2(10,6) 678 t2(9,9) = t2matrix (0,4,0,r5,bx,by) 679 t2(10,9) = t2matrix (1,3,0,r5,bx,by) 680 t2(10,10) = -t2(8,8) - t2(9,9) 681 t2(11,8) = t2(11,6) 682 t2(12,8) = t2(12,6) 683 t2(13,8) = t2(13,6) 684 t2(12,9) = t2(10,9) 685 t2(13,9) = t2(10,10) 686 t2(13,10) = -t2(11,8) - t2(12,9) 687 t2(11,9) = t2(12,8) 688 t2(11,10) = t2(13,8) 689 t2(12,10) = t2(13,9) 690 t2(11,11) = t2(11,7) 691 t2(12,11) = t2(12,7) 692 t2(13,11) = t2(13,7) 693 t2(12,12) = t2(12,10) 694 t2(13,12) = t2(13,10) 695 t2(13,13) = -t2(11,11) - t2(12,12) 696 do i = 5, 12 697 do j = i+1, 13 698 t2(i,j) = t2(j,i) 699 end do 700 end do 701 end if 702c 703c set some additional T matrix elements 704c 705 k = 1 706 m = 1 707 do i = 5, 13 708 do j = 1, 13 709 if (k .eq. 1) then 710 t3x(m,j) = t2(i,j) 711 else if (k .eq. 2) then 712 t3y(m,j) = t2(i,j) 713 else if (k .eq. 3) then 714 t3z(m,j) = t2(i,j) 715 end if 716 end do 717 k = k + 1 718 if (k .eq. 4) then 719 k = 1 720 m = m + 1 721 end if 722 end do 723 do i = 1, 13 724 do j = 2, 4 725 tt2(j-1,i) = t2(i,j) 726 end do 727 end do 728 do i = 1, 3 729 do j = 2, 4 730 tt2(i,j) = -tt2(i,j) 731 end do 732 end do 733c 734c compute the induced dipole polarization energy 735c 736 factor = electric / dielec 737 do i = 1, mdqsiz 738 m2t2(i) = 0.0d0 739 do j = 2, 4 740 m2t2(i) = m2t2(i) + indk(j-1)*t2(j,i) 741 end do 742 end do 743 epolik1 = 0.0d0 744 do i = 1, mdqsiz 745 epolik1 = epolik1 + m2t2(i)*rpi(i) 746 end do 747 epolik1 = 0.5d0 * factor * epolik1 748c 749c compute the d1 components for the induced dipole gradient 750c 751 do i = 1, 3 752 interx(i) = 0.0d0 753 intery(i) = 0.0d0 754 interz(i) = 0.0d0 755 do j = 1, mdqsiz 756 interx(i) = interx(i) + rpi(j)*t3x(i,j) 757 intery(i) = intery(i) + rpi(j)*t3y(i,j) 758 interz(i) = interz(i) + rpi(j)*t3z(i,j) 759 end do 760 end do 761 do i = 1, 3 762 d1(i) = 0.0d0 763 end do 764 do i = 1, 3 765 d1(1) = d1(1) + interx(i)*indk(i) 766 d1(2) = d1(2) + intery(i)*indk(i) 767 d1(3) = d1(3) + interz(i)*indk(i) 768 end do 769 do i = 1, 3 770 d1(i) = factor * d1(i) 771 end do 772c 773c compute the d2 components for the induced dipole gradient 774c 775 do i = 1, mdqsiz 776 m1t(i) = 0.0d0 777 do j = 1, 3 778 m1t(i) = m1t(i) + indk(j)*tt2(j,i) 779 end do 780 end do 781 do i = 1, 3 782 do j = 1, 3 783 d2(i,j) = 0.0d0 784 do k = 1, mdqsiz 785 d2(i,j) = d2(i,j) + m1t(k)*dpole(k,i,j,ii) 786 end do 787 d2(i,j) = factor * d2(i,j) 788 end do 789 end do 790c 791c compute the d3 components for the induced dipole gradient 792c 793 do i = 1, 3 794 interx(i) = 0.0d0 795 intery(i) = 0.0d0 796 interz(i) = 0.0d0 797 do j = 2, 4 798 interx(i) = interx(i) + indk(j-1)*t3x(i,j) 799 intery(i) = intery(i) + indk(j-1)*t3y(i,j) 800 interz(i) = interz(i) + indk(j-1)*t3z(i,j) 801 end do 802 end do 803 do i = 1, 3 804 d3(i) = 0.0d0 805 end do 806 do i = 1, 3 807 d3(1) = d3(1) + interx(i)*indi(i) 808 d3(2) = d3(2) + intery(i)*indi(i) 809 d3(3) = d3(3) + interz(i)*indi(i) 810 end do 811 do i = 1, 3 812 d3(i) = factor * d3(i) 813 end do 814c 815c compute the utu component for the induced dipole gradient 816c 817 utu = 0.0d0 818 do i = 1, 3 819 utu = utu + m2t2(i+1)*indi(i) 820 end do 821 utu = factor * utu 822 return 823 end 824c 825c 826c ################################################### 827c ## COPYRIGHT (C) 1995 by Jay William Ponder ## 828c ################################################### 829c 830c ################################################################# 831c ## ## 832c ## subroutine epolar2 -- atom-by-atom polarization Hessian ## 833c ## ## 834c ################################################################# 835c 836c 837c "epolar2" calculates second derivatives of the induced 838c dipole polarization energy for a single atom at a time 839c 840c 841 subroutine epolar2 (i) 842 implicit none 843 integer i 844c 845c 846 return 847 end 848c 849c 850c ################################################### 851c ## COPYRIGHT (C) 1995 by Jay William Ponder ## 852c ################################################### 853c 854c ################################################################ 855c ## ## 856c ## subroutine epolar3 -- polarization energy and analysis ## 857c ## ## 858c ################################################################ 859c 860c 861 subroutine epolar3 862 implicit none 863 include 'sizes.i' 864 include 'analyz.i' 865 include 'atmtyp.i' 866 include 'atoms.i' 867 include 'bound.i' 868 include 'couple.i' 869 include 'electr.i' 870 include 'energi.i' 871 include 'inform.i' 872 include 'iounit.i' 873 include 'moment.i' 874 include 'mpole.i' 875 include 'polar.i' 876 include 'shunt.i' 877 include 'units.i' 878 include 'usage.i' 879 integer i,j,k,skip(maxatm) 880 integer ii,kk,iz,ix 881 real*8 e,epolik,a(3,3) 882 real*8 xi,yi,zi,xr,yr,zr 883 real*8 r,r2,r3,r4,r5,taper 884 real*8 rpi(13),indk(3) 885 real*8 utotal,xsum,ysum,zsum 886 logical iuse,kuse,huge,header 887c 888c 889c zero out the induced dipole energy and partitioning 890c 891 ep = 0.0d0 892 do i = 1, n 893 aep(i) = 0.0d0 894 end do 895 nplr = 0 896 header = .true. 897c 898c zero out the list of atoms to be skipped 899c 900 do i = 1, n 901 skip(i) = 0 902 end do 903c 904c set the switching function coefficients 905c 906 call switch ('CHGDPL') 907c 908c rotate the multipole components into the global frame 909c 910 call rotpole 911c 912c compute the induced dipoles at each atom 913c 914 call induce 915c 916c add induced dipole components to the permanent components 917c 918 xsum = 0.0d0 919 ysum = 0.0d0 920 zsum = 0.0d0 921 do ii = 1, npolar 922 i = ipolar(ii) 923 xsum = xsum + uind(1,i) 924 ysum = ysum + uind(2,i) 925 zsum = zsum + uind(3,i) 926 end do 927 xdipole = xdipole + debye*xsum 928 ydipole = ydipole + debye*ysum 929 zdipole = zdipole + debye*zsum 930c 931c compute and partition the induced dipole energy 932c 933 do ii = 1, npole 934 i = ipole(ii) 935 iz = zaxis(ii) 936 ix = xaxis(ii) 937 iuse = (use(i) .or. use(iz) .or. use(ix)) 938 xi = x(i) 939 yi = y(i) 940 zi = z(i) 941 skip(i) = i 942 do j = 1, n12(i) 943 skip(i12(j,i)) = i 944 end do 945 do j = 1, n13(i) 946 skip(i13(j,i)) = i 947 end do 948 do j = 1, mdqsiz 949 rpi(j) = rpole(j,ii) 950 end do 951 do kk = 1, npolar 952 k = ipolar(kk) 953 kuse = use(k) 954 if (iuse .or. kuse) then 955 if (skip(k) .ne. i) then 956 xr = x(k) - xi 957 yr = y(k) - yi 958 zr = z(k) - zi 959 if (use_image) call image (xr,yr,zr,0) 960 r2 = xr*xr + yr*yr + zr*zr 961 if (r2 .le. off2) then 962 indk(1) = uind(1,k) 963 indk(2) = uind(2,k) 964 indk(3) = uind(3,k) 965 r = sqrt(r2) 966 e = epolik (r,xr,yr,zr,rpi,indk) 967 if (r2 .gt. cut2) then 968 r3 = r2 * r 969 r4 = r2 * r2 970 r5 = r2 * r3 971 taper = c5*r5 + c4*r4 + c3*r3 972 & + c2*r2 + c1*r + c0 973 e = e * taper 974 end if 975 nplr = nplr + 1 976 ep = ep + e 977 aep(i) = aep(i) + 0.5d0*e 978 aep(k) = aep(k) + 0.5d0*e 979c 980c print a warning if the energy of dipole interaction is large 981c 982 huge = (abs(e) .gt. 1.0d0) 983 if (debug .or. (verbose.and.huge)) then 984 if (header) then 985 header = .false. 986 write (iout,10) 987 10 format (/,' Individual Induced Dipole', 988 & ' Interactions :', 989 & //,' Type',11x,'Atom Names', 990 & 13x,'IndDpl - Chg',3x,'Distance', 991 & 5x,'Energy',/) 992 end if 993 utotal = sqrt(uind(1,k)**2 + uind(2,k)**2 994 & + uind(3,k)**2) * debye 995 write (iout,20) k,name(k),i,name(i), 996 & utotal,pole(1,ii),r,e 997 20 format (' Polarize ',i5,'-',a3,1x,i5,'-', 998 & a3,8x,2f7.2,f10.4,f12.4) 999 end if 1000 end if 1001 end if 1002 end if 1003 end do 1004 end do 1005 return 1006 end 1007