1c 2c 3c ############################################################# 4c ## COPYRIGHT (C) 1995 by Yong Kong & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ############################################################# 7c 8c ################################################################ 9c ## ## 10c ## subroutine induce0a -- induced dipoles via double loop ## 11c ## ## 12c ################################################################ 13c 14c 15c "induce0a" computes the induced dipole moment at each 16c polarizable site using a pairwise double loop 17c 18c 19 subroutine induce0a 20 implicit none 21 include 'sizes.i' 22 include 'atoms.i' 23 include 'bound.i' 24 include 'cell.i' 25 include 'chgpot.i' 26 include 'couple.i' 27 include 'inform.i' 28 include 'iounit.i' 29 include 'mpole.i' 30 include 'polar.i' 31 include 'polpot.i' 32 include 'potent.i' 33 include 'shunt.i' 34 include 'units.i' 35 include 'usage.i' 36 integer i,j,k,m,skip(maxatm) 37 integer ii,iz,ix,kk,kz,kx 38 integer iter,maxiter 39 real*8 xr,yr,zr,eps,taper 40 real*8 r,r2,r3,r4,r5 41 real*8 rpi(13),rpk(13) 42 real*8 udir(3,maxatm) 43 real*8 uold(3,maxatm) 44 real*8 fieldi(3),fieldk(3) 45 logical iuse,kuse 46c 47c 48c zero out induced dipoles and count the polarizable sites 49c 50 do i = 1, npole 51 do j = 1, 3 52 uind(j,i) = 0.0d0 53 end do 54 end do 55 if (.not. use_polar) return 56c 57c zero out the list of atoms to be skipped 58c 59 do i = 1, n 60 skip(i) = 0 61 end do 62c 63c set the switching function coefficients 64c 65 call switch ('MPOLE') 66c 67c compute the direct induced dipole moment at each atom 68c 69 do ii = 1, npole-1 70 i = ipole(ii) 71 iz = zaxis(ii) 72 ix = xaxis(ii) 73 iuse = (use(i) .or. use(iz) .or. use(ix)) 74 do j = 1, n12(i) 75 skip(i12(j,i)) = i * chg12use 76 end do 77 do j = 1, n13(i) 78 skip(i13(j,i)) = i * chg13use 79 end do 80 do j = 1, n14(i) 81 skip(i14(j,i)) = i * chg14use 82 end do 83 do j = 1, maxpole 84 rpi(j) = rpole(j,ii) 85 end do 86 do kk = ii+1, npole 87 k = ipole(kk) 88 kz = zaxis(kk) 89 kx = xaxis(kk) 90 kuse = (use(k) .or. use(kz) .or. use(kx)) 91 if (iuse .or. kuse) then 92 if (skip(k) .ne. i) then 93 xr = x(k) - x(i) 94 yr = y(k) - y(i) 95 zr = z(k) - z(i) 96 if (use_image) call image (xr,yr,zr,0) 97 r2 = xr*xr + yr*yr + zr*zr 98 if (r2 .le. off2) then 99 r = sqrt(r2) 100 do j = 1, maxpole 101 rpk(j) = rpole(j,kk) 102 end do 103 call udirect (ii,kk,xr,yr,zr,r,r2,rpi, 104 & rpk,fieldi,fieldk) 105 if (skip(k) .eq. -i) then 106 do j = 1, 3 107 fieldi(j) = fieldi(j) / chgscale 108 fieldk(j) = fieldk(j) / chgscale 109 end do 110 end if 111 if (r2 .gt. cut2) then 112 r3 = r2 * r 113 r4 = r2 * r2 114 r5 = r2 * r3 115 taper = c5*r5 + c4*r4 + c3*r3 116 & + c2*r2 + c1*r + c0 117 do j = 1, 3 118 fieldi(j) = fieldi(j) * taper 119 fieldk(j) = fieldk(j) * taper 120 end do 121 end if 122 do j = 1, 3 123 uind(j,ii) = uind(j,ii) + fieldi(j) 124 uind(j,kk) = uind(j,kk) + fieldk(j) 125 end do 126 end if 127 end if 128 end if 129 end do 130 end do 131c 132c periodic boundary for large cutoffs via replicates method 133c 134 if (use_replica) then 135 do ii = 1, npole 136 i = ipole(ii) 137 iz = zaxis(ii) 138 ix = xaxis(ii) 139 iuse = (use(i) .or. use(iz) .or. use(ix)) 140 do j = 1, maxpole 141 rpi(j) = rpole(j,ii) 142 end do 143 do kk = ii, npole 144 k = ipole(kk) 145 kz = zaxis(kk) 146 kx = xaxis(kk) 147 kuse = (use(k) .or. use(kz) .or. use(kx)) 148 if (iuse .or. kuse) then 149 do m = 2, ncell 150 xr = x(k) - x(i) 151 yr = y(k) - y(i) 152 zr = z(k) - z(i) 153 call image (xr,yr,zr,m) 154 r2 = xr*xr + yr*yr + zr*zr 155 if (r2 .le. off2) then 156 r = sqrt(r2) 157 do j = 1, maxpole 158 rpk(j) = rpole(j,kk) 159 end do 160 call udirect (ii,kk,xr,yr,zr,r,r2,rpi, 161 & rpk,fieldi,fieldk) 162 if (r2 .gt. cut2) then 163 r3 = r2 * r 164 r4 = r2 * r2 165 r5 = r2 * r3 166 taper = c5*r5 + c4*r4 + c3*r3 167 & + c2*r2 + c1*r + c0 168 do j = 1, 3 169 fieldi(j) = fieldi(j) * taper 170 fieldk(j) = fieldk(j) * taper 171 end do 172 end if 173 do j = 1, 3 174 uind(j,ii) = uind(j,ii) + fieldi(j) 175 if (ii .ne. kk) 176 & uind(j,kk) = uind(j,kk) + fieldk(j) 177 end do 178 end if 179 end do 180 end if 181 end do 182 end do 183 end if 184c 185c convert total field at each atom to induced dipole moment 186c 187 do i = 1, npole 188 do j = 1, 3 189 uind(j,i) = polarize(i) * uind(j,i) 190 end do 191 end do 192c 193c set tolerances for computation of mutual induced dipoles 194c 195 if (poltyp .eq. 'MUTUAL') then 196 maxiter = 100 197 iter = 0 198 eps = 1.0d0 199 do i = 1, npole 200 do j = 1, 3 201 udir(j,i) = uind(j,i) 202 end do 203 end do 204c 205c compute mutual induced dipole moments by an iterative method 206c 207 dowhile (iter.lt.maxiter .and. eps.gt.poleps) 208 do ii = 1, npole 209 do j = 1, 3 210 uold(j,ii) = uind(j,ii) 211 uind(j,ii) = 0.0d0 212 end do 213 end do 214 do ii = 1, npole-1 215 i = ipole(ii) 216 iz = zaxis(ii) 217 ix = xaxis(ii) 218 iuse = (use(i) .or. use(iz) .or. use(ix)) 219 do j = 1, n12(i) 220 skip(i12(j,i)) = i * chg12use 221 end do 222 do j = 1, n13(i) 223 skip(i13(j,i)) = i * chg13use 224 end do 225 do j = 1, n14(i) 226 skip(i14(j,i)) = i * chg14use 227 end do 228 do kk = ii+1, npole 229 k = ipole(kk) 230 kz = zaxis(kk) 231 kx = xaxis(kk) 232 kuse = (use(k) .or. use(kz) .or. use(kx)) 233 if (iuse .or. kuse) then 234 if (skip(k) .ne. i) then 235 xr = x(k) - x(i) 236 yr = y(k) - y(i) 237 zr = z(k) - z(i) 238 if (use_image) call image (xr,yr,zr,0) 239 r2 = xr*xr + yr*yr + zr*zr 240 if (r2 .le. off2) then 241 r = sqrt(r2) 242 call umutual (ii,kk,xr,yr,zr,r,r2,uold(1,ii), 243 & uold(1,kk),fieldi,fieldk) 244 if (skip(k) .eq. -i) then 245 do j = 1, 3 246 fieldi(j) = fieldi(j) / chgscale 247 fieldk(j) = fieldk(j) / chgscale 248 end do 249 end if 250 if (r2 .gt. cut2) then 251 r3 = r2 * r 252 r4 = r2 * r2 253 r5 = r2 * r3 254 taper = c5*r5 + c4*r4 + c3*r3 255 & + c2*r2 + c1*r + c0 256 do j = 1, 3 257 fieldi(j) = fieldi(j) * taper 258 fieldk(j) = fieldk(j) * taper 259 end do 260 end if 261 do j = 1, 3 262 uind(j,ii) = uind(j,ii) + fieldi(j) 263 uind(j,kk) = uind(j,kk) + fieldk(j) 264 end do 265 end if 266 end if 267 end if 268 end do 269 end do 270c 271c periodic boundary for large cutoffs via replicates method 272c 273 if (use_replica) then 274 do ii = 1, npole 275 i = ipole(ii) 276 iz = zaxis(ii) 277 ix = xaxis(ii) 278 iuse = (use(i) .or. use(iz) .or. use(ix)) 279 do kk = ii, npole 280 k = ipole(kk) 281 kz = zaxis(kk) 282 kx = xaxis(kk) 283 kuse = (use(k) .or. use(kz) .or. use(kx)) 284 if (iuse .or. kuse) then 285 do m = 2, ncell 286 xr = x(k) - x(i) 287 yr = y(k) - y(i) 288 zr = z(k) - z(i) 289 call image (xr,yr,zr,m) 290 r2 = xr*xr + yr*yr + zr*zr 291 if (r2 .le. off2) then 292 r = sqrt(r2) 293 call umutual (ii,kk,xr,yr,zr,r,r2, 294 & uold(1,ii),uold(1,kk), 295 & fieldi,fieldk) 296 if (r2 .gt. cut2) then 297 r3 = r2 * r 298 r4 = r2 * r2 299 r5 = r2 * r3 300 taper = c5*r5 + c4*r4 + c3*r3 301 & + c2*r2 + c1*r + c0 302 do j = 1, 3 303 fieldi(j) = fieldi(j) * taper 304 fieldk(j) = fieldk(j) * taper 305 end do 306 end if 307 do j = 1, 3 308 uind(j,ii) = uind(j,ii) + fieldi(j) 309 if (ii .ne. kk) 310 & uind(j,kk) = uind(j,kk) + fieldk(j) 311 end do 312 end if 313 end do 314 end if 315 end do 316 end do 317 end if 318c 319c check to see if the mutual induced dipoles have converged 320c 321 iter = iter + 1 322 eps = 0.0d0 323 do i = 1, npole 324 do j = 1, 3 325 uind(j,i) = udir(j,i) + polarize(i)*uind(j,i) 326 eps = eps + (uind(j,i)-uold(j,i))**2 327 end do 328 end do 329 eps = debye * sqrt(eps/dble(npolar)) 330 if (debug) then 331 if (iter .eq. 1) then 332 write (iout,10) 333 10 format (/,' Determination of Induced Dipole', 334 & ' Moments :', 335 & //,4x,'Iter',8x,'RMS Change (Debyes)',/) 336 end if 337 write (iout,20) iter,eps 338 20 format (i8,7x,f16.10) 339 end if 340 end do 341c 342c terminate the calculation if dipoles failed to converge 343c 344 if (eps .gt. poleps) then 345 write (iout,30) 346 30 format (/,' INDUCE -- Warning, Induced Dipoles', 347 & ' are not Converged') 348 call prterr 349 call fatal 350 end if 351 end if 352 return 353 end 354c 355c 356c ################################################################ 357c ## ## 358c ## subroutine udirect -- field for direct induced dipoles ## 359c ## ## 360c ################################################################ 361c 362c 363c "udirect" evaluates the electric field at a polarizable atom 364c due to permanent atomic multipoles at a second atom, and vice 365c versa, for use in computation of direct induced dipole moments 366c 367c 368 subroutine udirect (ii,kk,xr,yr,zr,r,r2,rpi,rpk,fieldi,fieldk) 369 implicit none 370 include 'sizes.i' 371 include 'mpole.i' 372 include 'polar.i' 373 include 'polpot.i' 374 integer i,ii,kk 375 real*8 xr,yr,zr 376 real*8 r,r2,damp 377 real*8 rr3,rr5,rr7 378 real*8 xr5,yr5,xyz 379 real*8 xr7,yr7,rr53 380 real*8 t2,t3,t4,t5,t6,t7,t8,t9,t10,t11 381 real*8 t12,t13,t14,t15,t16,t17,t18,t19 382 real*8 i0,i1,i2,i3,i4,i5,i6,i7,i8,i9 383 real*8 k0,k1,k2,k3,k4,k5,k6,k7,k8,k9 384 real*8 i9i4,i9i7,k4k9,k7k9 385 real*8 rpi(13),rpk(13) 386 real*8 fieldi(3),fieldk(3) 387 logical iquad,kquad 388c 389c 390c set extent of the multipole expansion at each site 391c 392 iquad = (polsiz(ii) .ge. 13) 393 kquad = (polsiz(kk) .ge. 13) 394c 395c compute the first order T2 matrix elements 396c 397 rr3 = -1.0d0 / (r2*r) 398 t2 = xr * rr3 399 t3 = yr * rr3 400 t4 = zr * rr3 401c 402c compute the second order T2 matrix elements 403c 404 rr5 = -3.0d0 * rr3 / r2 405 xr5 = xr * rr5 406 yr5 = yr * rr5 407 t5 = rr3 + xr5*xr 408 t6 = xr5 * yr 409 t7 = xr5 * zr 410 t8 = rr3 + yr5*yr 411 t9 = yr5 * zr 412 t10 = -t5 - t8 413c 414c compute the third order T2 matrix elements 415c 416 if (iquad .or. kquad) then 417 rr7 = -5.0d0 * rr5 / r2 418 xyz = xr * yr * zr 419 xr7 = xr * xr * rr7 420 yr7 = yr * yr * rr7 421 rr53 = 3.0d0 * rr5 422 t11 = xr * (xr7+rr53) 423 t12 = yr * (xr7+rr5) 424 t13 = zr * (xr7+rr5) 425 t14 = xr * (yr7+rr5) 426 t15 = xyz * rr7 427 t16 = -t11 - t14 428 t17 = yr * (yr7+rr53) 429 t18 = zr * (yr7+rr5) 430 t19 = -t12 - t17 431 end if 432c 433c compute the field at site k due to multipoles at site i 434c 435 i0 = rpi(1) 436 i1 = rpi(2) 437 i2 = rpi(3) 438 i3 = rpi(4) 439 fieldk(1) = -i0*t2 + i1*t5 + i2*t6 + i3*t7 440 fieldk(2) = -i0*t3 + i1*t6 + i2*t8 + i3*t9 441 fieldk(3) = -i0*t4 + i1*t7 + i2*t9 + i3*t10 442 if (iquad) then 443 i4 = rpi(5) 444 i5 = rpi(6) 445 i6 = rpi(7) 446 i7 = rpi(9) 447 i8 = rpi(10) 448 i9 = rpi(13) 449 i9i4 = i9 - i4 450 i9i7 = i9 - i7 451 fieldk(1) = fieldk(1) + i9i4*t11 + i9i7*t14 452 & - 2.0d0*(i5*t12 + i6*t13 + i8*t15) 453 fieldk(2) = fieldk(2) + i9i4*t12 + i9i7*t17 454 & - 2.0d0*(i5*t14 + i6*t15 + i8*t18) 455 fieldk(3) = fieldk(3) + i9i4*t13 + i9i7*t18 456 & - 2.0d0*(i5*t15 + i6*t16 + i8*t19) 457 end if 458c 459c compute the field at site i due to multipoles at site k 460c 461 k0 = rpk(1) 462 k1 = rpk(2) 463 k2 = rpk(3) 464 k3 = rpk(4) 465 fieldi(1) = k0*t2 + k1*t5 + k2*t6 + k3*t7 466 fieldi(2) = k0*t3 + k1*t6 + k2*t8 + k3*t9 467 fieldi(3) = k0*t4 + k1*t7 + k2*t9 + k3*t10 468 if (kquad) then 469 k4 = rpk(5) 470 k5 = rpk(6) 471 k6 = rpk(7) 472 k7 = rpk(9) 473 k8 = rpk(10) 474 k9 = rpk(13) 475 k4k9 = k4 - k9 476 k7k9 = k7 - k9 477 fieldi(1) = fieldi(1) + k4k9*t11 + k7k9*t14 478 & + 2.0d0*(k5*t12 + k6*t13 + k8*t15) 479 fieldi(2) = fieldi(2) + k4k9*t12 + k7k9*t17 480 & + 2.0d0*(k5*t14 + k6*t15 + k8*t18) 481 fieldi(3) = fieldi(3) + k4k9*t13 + k7k9*t18 482 & + 2.0d0*(k5*t15 + k6*t16 + k8*t19) 483 end if 484c 485c apply a damping factor to reduce the field at short range 486c 487 damp = pdamp(ii) * pdamp(kk) 488 if (damp .ne. 0.0d0) then 489 damp = -pgamma * (r/damp)**3 490 if (damp .gt. -50.0d0) then 491 damp = 1.0d0 - exp(damp) 492 do i = 1, 3 493 fieldi(i) = fieldi(i) * damp 494 fieldk(i) = fieldk(i) * damp 495 end do 496 end if 497 end if 498 return 499 end 500c 501c 502c ################################################################ 503c ## ## 504c ## subroutine umutual -- field for mutual induced dipoles ## 505c ## ## 506c ################################################################ 507c 508c 509c "umutual" evaluates the electric field at a polarizable atom 510c due to the induced atomic dipoles at a second atom, and vice 511c versa, for use in computation of mutual induced dipole moments 512c 513c 514 subroutine umutual (ii,kk,xr,yr,zr,r,r2,upi,upk,fieldi,fieldk) 515 implicit none 516 include 'sizes.i' 517 include 'mpole.i' 518 include 'polar.i' 519 include 'polpot.i' 520 integer i,ii,kk 521 real*8 xr,yr,zr 522 real*8 r,r2,damp 523 real*8 rr3,rr5,xr5,yr5 524 real*8 upi(3),upk(3) 525 real*8 fieldi(3),fieldk(3) 526 real*8 u1,u2,u3,v1,v2,v3 527 real*8 t5,t6,t7,t8,t9,t10 528c 529c 530c set the current values of the induced dipoles 531c 532 u1 = upi(1) 533 u2 = upi(2) 534 u3 = upi(3) 535 v1 = upk(1) 536 v2 = upk(2) 537 v3 = upk(3) 538c 539c compute the second order T2 matrix elements 540c 541 rr3 = -1.0d0 / (r2*r) 542 rr5 = -3.0d0 * rr3 / r2 543 xr5 = xr * rr5 544 yr5 = yr * rr5 545 t5 = rr3 + xr5*xr 546 t6 = xr5 * yr 547 t7 = xr5 * zr 548 t8 = rr3 + yr5*yr 549 t9 = yr5 * zr 550 t10 = -t5 - t8 551c 552c compute the field at site i due to site k and vice versa 553c 554 fieldi(1) = v1*t5 + v2*t6 + v3*t7 555 fieldi(2) = v1*t6 + v2*t8 + v3*t9 556 fieldi(3) = v1*t7 + v2*t9 + v3*t10 557 fieldk(1) = u1*t5 + u2*t6 + u3*t7 558 fieldk(2) = u1*t6 + u2*t8 + u3*t9 559 fieldk(3) = u1*t7 + u2*t9 + u3*t10 560c 561c apply a damping factor to reduce the field at short range 562c 563 damp = pdamp(ii) * pdamp(kk) 564 if (damp .ne. 0.0d0) then 565 damp = -pgamma * (r/damp)**3 566 if (damp .gt. -50.0d0) then 567 damp = 1.0d0 - exp(damp) 568 do i = 1, 3 569 fieldi(i) = fieldi(i) * damp 570 fieldk(i) = fieldk(i) * damp 571 end do 572 end if 573 end if 574 return 575 end 576c 577c 578c ############################################################# 579c ## ## 580c ## subroutine empole0a -- double loop multipole energy ## 581c ## ## 582c ############################################################# 583c 584c 585c "empole0a" calculates the electrostatic energy due to 586c atomic multipole interactions and dipole polarizability 587c using a pairwise double loop 588c 589c 590 subroutine empole0a 591 implicit none 592 include 'sizes.i' 593 include 'atoms.i' 594 include 'bound.i' 595 include 'cell.i' 596 include 'chgpot.i' 597 include 'couple.i' 598 include 'energi.i' 599 include 'group.i' 600 include 'mpole.i' 601 include 'polar.i' 602 include 'shunt.i' 603 include 'units.i' 604 include 'usage.i' 605 integer i,j,k,m 606 integer ii,iz,ix 607 integer kk,kz,kx 608 integer skip(maxatm) 609 real*8 eik,ei,ek,a(3,3) 610 real*8 shift,taper,trans 611 real*8 f,fik,fgrp 612 real*8 xr,yr,zr,r 613 real*8 r2,r3,r4,r5,r6,r7 614 real*8 rpi(13),rpk(13) 615 real*8 indi(3),indk(3) 616 logical proceed,iuse,kuse 617c 618c 619c zero out the atomic multipole interaction energy 620c 621 em = 0.0d0 622 ep = 0.0d0 623c 624c zero out the list of atoms to be skipped 625c 626 do i = 1, n 627 skip(i) = 0 628 end do 629c 630c set conversion factor and cutoff values 631c 632 f = electric / dielec 633 call switch ('MPOLE') 634c 635c rotate the multipole components into the global frame 636c 637 call rotpole 638c 639c compute the induced dipoles at each atom 640c 641 call induce 642c 643c calculate the multipole interaction energy term 644c 645 do ii = 1, npole-1 646 i = ipole(ii) 647 iz = zaxis(ii) 648 ix = xaxis(ii) 649 iuse = (use(i) .or. use(iz) .or. use(ix)) 650 do j = 1, n12(i) 651 skip(i12(j,i)) = i * chg12use 652 end do 653 do j = 1, n13(i) 654 skip(i13(j,i)) = i * chg13use 655 end do 656 do j = 1, n14(i) 657 skip(i14(j,i)) = i * chg14use 658 end do 659 do j = 1, maxpole 660 rpi(j) = rpole(j,ii) 661 end do 662 do j = 1, 3 663 indi(j) = uind(j,ii) 664 end do 665 do kk = ii+1, npole 666 k = ipole(kk) 667 kz = zaxis(kk) 668 kx = xaxis(kk) 669 kuse = (use(k) .or. use(kz) .or. use(kx)) 670c 671c decide whether to compute the current interaction 672c 673 proceed = .true. 674 if (use_group) call groups (proceed,fgrp,2,i,k,0,0,0) 675 if (proceed) proceed = (iuse .or. kuse) 676 if (proceed) proceed = (skip(k) .ne. i) 677c 678c compute the energy contribution for this interaction 679c 680 if (proceed) then 681 xr = x(k) - x(i) 682 yr = y(k) - y(i) 683 zr = z(k) - z(i) 684 if (use_image) call image (xr,yr,zr,0) 685 r2 = xr*xr + yr*yr + zr*zr 686 if (r2 .le. off2) then 687 do j = 1, maxpole 688 rpk(j) = rpole(j,kk) 689 end do 690 do j = 1, 3 691 indk(j) = uind(j,kk) 692 end do 693 r = sqrt(r2) 694 call empik (ii,kk,xr,yr,zr,r,rpi,rpk, 695 & indi,indk,eik,ei,ek) 696 fik = f 697 if (skip(k) .eq. -i) fik = fik / chgscale 698 eik = fik * eik 699 ei = fik * ei 700 ek = fik * ek 701c 702c use shifted energy switching if near the cutoff distance 703c 704 fik = fik * rpi(1) * rpk(1) 705 shift = fik / (0.5d0*(off+cut)) 706 eik = eik - shift 707 if (r2 .gt. cut2) then 708 r3 = r2 * r 709 r4 = r2 * r2 710 r5 = r2 * r3 711 r6 = r3 * r3 712 r7 = r3 * r4 713 taper = c5*r5 + c4*r4 + c3*r3 714 & + c2*r2 + c1*r + c0 715 trans = fik * (f7*r7 + f6*r6 + f5*r5 + f4*r4 716 & + f3*r3 + f2*r2 + f1*r + f0) 717 eik = eik * taper + trans 718 ei = ei * taper 719 ek = ek * taper 720 end if 721c 722c scale the interaction based on its group membership 723c 724 if (use_group) then 725 eik = eik * fgrp 726 ei = ei * fgrp 727 ek = ek * fgrp 728 end if 729c 730c increment the overall multipole and polarization energies 731c 732 em = em + eik 733 ep = ep + ei + ek 734 end if 735 end if 736 end do 737 end do 738c 739c for periodic boundary conditions with large cutoffs 740c neighbors must be found by the replicates method 741c 742 if (.not. use_replica) return 743c 744c calculate interaction energy with other unit cells 745c 746 do ii = 1, npole 747 i = ipole(ii) 748 iz = zaxis(ii) 749 ix = xaxis(ii) 750 iuse = (use(i) .or. use(iz) .or. use(ix)) 751 do j = 1, maxpole 752 rpi(j) = rpole(j,ii) 753 end do 754 do j = 1, 3 755 indi(j) = uind(j,ii) 756 end do 757 do kk = ii, npole 758 k = ipole(kk) 759 kz = zaxis(kk) 760 kx = xaxis(kk) 761 kuse = (use(k) .or. use(kz) .or. use(kx)) 762c 763c decide whether to compute the current interaction 764c 765 proceed = .true. 766 if (use_group) call groups (proceed,fgrp,2,i,k,0,0,0) 767 if (proceed) proceed = (iuse .or. kuse) 768c 769c compute the energy contribution for this interaction 770c 771 if (proceed) then 772 do m = 2, ncell 773 xr = x(k) - x(i) 774 yr = y(k) - y(i) 775 zr = z(k) - z(i) 776 call image (xr,yr,zr,m) 777 r2 = xr*xr + yr*yr + zr*zr 778 if (r2 .le. off2) then 779 do j = 1, maxpole 780 rpk(j) = rpole(j,kk) 781 end do 782 do j = 1, 3 783 indk(j) = uind(j,kk) 784 end do 785 r = sqrt(r2) 786 call empik (ii,kk,xr,yr,zr,r,rpi,rpk, 787 & indi,indk,eik,ei,ek) 788 fik = f 789 eik = fik * eik 790 ei = fik * ei 791 ek = fik * ek 792c 793c use shifted energy switching if near the cutoff distance 794c 795 fik = fik * rpi(1) * rpk(1) 796 shift = fik / (0.5d0*(off+cut)) 797 eik = eik - shift 798 if (r2 .gt. cut2) then 799 r3 = r2 * r 800 r4 = r2 * r2 801 r5 = r2 * r3 802 r6 = r3 * r3 803 r7 = r3 * r4 804 taper = c5*r5 + c4*r4 + c3*r3 805 & + c2*r2 + c1*r + c0 806 trans = fik * (f7*r7 + f6*r6 + f5*r5 + f4*r4 807 & + f3*r3 + f2*r2 + f1*r + f0) 808 eik = eik * taper + trans 809 ei = ei * taper 810 ek = ek * taper 811 end if 812c 813c scale the interaction based on its group membership 814c 815 if (use_group) then 816 eik = eik * fgrp 817 ei = ei * fgrp 818 ek = ek * fgrp 819 end if 820c 821c increment the overall multipole and polarization energies 822c 823 if (i .eq. k) then 824 eik = 0.5d0 * eik 825 ei = 0.5d0 * ei 826 ek = 0.5d0 * ek 827 end if 828 em = em + eik 829 ep = ep + ei + ek 830 end if 831 end do 832 end if 833 end do 834 end do 835 return 836 end 837c 838c 839c ################################################################## 840c ## ## 841c ## subroutine empik -- multipole & polarization pair energy ## 842c ## ## 843c ################################################################## 844c 845c 846c "empik" computes the permanent multipole and induced dipole 847c energies between a specified pair of atomic multipole sites 848c 849c 850 subroutine empik (ii,kk,xr,yr,zr,r,rpi,rpk,indi,indk,eik,ei,ek) 851 implicit none 852 include 'sizes.i' 853 include 'atoms.i' 854 include 'mpole.i' 855 include 'polar.i' 856 include 'polpot.i' 857 include 'units.i' 858 integer ii,kk 859 real*8 eik,ei,ek 860 real*8 xr,yr,zr,r,damp 861 real*8 rr1,rr2,rr3,rr5,rr7,rr9 862 real*8 xx,yy,xyz,xry,yrx,zrx,zry 863 real*8 xr5,yr5,xr7,yr7 864 real*8 xr9,xyr9,yr9,xr9x7,yr9y7 865 real*8 rr53,xr73,yr73 866 real*8 t1,t2,t3,t4,t5,t6,t7,t8,t9,t11 867 real*8 t12,t13,t14,t15,t17,t18,t21,t22 868 real*8 t23,t24,t25,t27,t28,t31,t32 869 real*8 i0,i1,i2,i3,i4,i5,i6,i7,i8,i9 870 real*8 k0,k1,k2,k3,k4,k5,k6,k7,k8,k9 871 real*8 u1,u2,u3,v1,v2,v3 872 real*8 i3k3,i4i9,i7i9,k4k9,k7k9,i3k6 873 real*8 i3k8,i6k3,i8k3,ik49,ik59,ik69 874 real*8 ik79,ik89,ik68,i6k6,i8k8,i9k9 875 real*8 i3v3,i6v3,i8v3,k3u3,k6u3,k8u3 876 real*8 rpi(13),rpk(13) 877 real*8 indi(3),indk(3) 878 logical iquad,kquad 879c 880c 881c zero out multipole and polarization energy components 882c 883 eik = 0.0d0 884 ei = 0.0d0 885 ek = 0.0d0 886c 887c check for presence of quadrupole components at either site 888c 889 iquad = (polsiz(ii) .ge. 13) 890 kquad = (polsiz(kk) .ge. 13) 891c 892c set permanent and induced multipoles for first site 893c 894 i0 = rpi(1) 895 i1 = rpi(2) 896 i2 = rpi(3) 897 i3 = rpi(4) 898 i4 = rpi(5) 899 i5 = rpi(6) 900 i6 = rpi(7) 901 i7 = rpi(9) 902 i8 = rpi(10) 903 i9 = rpi(13) 904 u1 = indi(1) 905 u2 = indi(2) 906 u3 = indi(3) 907c 908c set permanent and induced multipoles for second site 909c 910 k0 = rpk(1) 911 k1 = rpk(2) 912 k2 = rpk(3) 913 k3 = rpk(4) 914 k4 = rpk(5) 915 k5 = rpk(6) 916 k6 = rpk(7) 917 k7 = rpk(9) 918 k8 = rpk(10) 919 k9 = rpk(13) 920 v1 = indk(1) 921 v2 = indk(2) 922 v3 = indk(3) 923c 924c compute the zeroth order T2 matrix element 925c 926 rr1 = 1.0d0 / r 927 t1 = rr1 928c 929c compute the first order T2 matrix elements 930c 931 rr2 = rr1 * rr1 932 rr3 = rr1 * rr2 933 t2 = -xr * rr3 934 t3 = -yr * rr3 935 t4 = -zr * rr3 936c 937c compute the second order T2 matrix elements 938c 939 rr5 = 3.0d0 * rr3 * rr2 940 xr5 = xr * rr5 941 yr5 = yr * rr5 942 t5 = -rr3 + xr5*xr 943 t6 = xr5 * yr 944 t7 = xr5 * zr 945 t8 = -rr3 + yr5*yr 946 t9 = yr5 * zr 947c 948c compute the third order T2 matrix elements 949c 950 if (iquad .or. kquad) then 951 rr7 = 5.0d0 * rr5 * rr2 952 xx = xr * xr 953 yy = yr * yr 954 xyz = xr * yr * zr 955 xr7 = xx * rr7 956 yr7 = yy * rr7 957 rr53 = 3.0d0 * rr5 958 t11 = xr * (rr53-xr7) 959 t12 = yr * (rr5-xr7) 960 t13 = zr * (rr5-xr7) 961 t14 = xr * (rr5-yr7) 962 t15 = -xyz * rr7 963 t17 = yr * (rr53-yr7) 964 t18 = zr * (rr5-yr7) 965c 966c compute the fourth order T2 matrix elements 967c 968 rr9 = 7.0d0 * rr7 * rr2 969 if (xr .eq. 0.0d0) then 970 yrx = 0.0d0 971 zrx = 0.0d0 972 else 973 yrx = yr / xr 974 zrx = zr / xr 975 end if 976 if (yr .eq. 0.0d0) then 977 xry = 0.0d0 978 zry = 0.0d0 979 else 980 xry = xr / yr 981 zry = zr / yr 982 end if 983 xr9 = xx * xx * rr9 984 xyr9 = xx * yy * rr9 985 yr9 = yy * yy * rr9 986 xr73 = 3.0d0 * xr7 987 yr73 = 3.0d0 * yr7 988 xr9x7 = xr9 - xr73 989 yr9y7 = yr9 - yr73 990 t21 = xr9x7 - xr73 + rr53 991 t22 = yrx * xr9x7 992 t23 = zrx * xr9x7 993 t24 = xyr9 - xr7 - yr7 + rr5 994 t25 = zry * (xyr9-yr7) 995 t27 = xry * yr9y7 996 t28 = zrx * (xyr9-xr7) 997 t31 = yr9y7 - yr73 + rr53 998 t32 = zry * yr9y7 999 end if 1000c 1001c get the M-M, M-D and D-D parts of the multipole energy 1002c 1003 i3k3 = i3 * k3 1004 eik = i0*k0*t1 + (i0*k1-i1*k0)*t2 + (i0*k2-i2*k0)*t3 1005 & + (i0*k3-i3*k0)*t4 + (i3k3-i1*k1)*t5 - (i1*k2+i2*k1)*t6 1006 & - (i1*k3+i3*k1)*t7 + (i3k3-i2*k2)*t8 - (i2*k3+i3*k2)*t9 1007c 1008c get the M-indD and D-indD parts of the polarization energy 1009c 1010 i3v3 = i3 * v3 1011 k3u3 = k3 * u3 1012 ei = -k0*(u1*t2+u2*t3+u3*t4) + (k3u3-k1*u1)*t5 - (k1*u2+k2*u1)*t6 1013 & - (k3*u1+k1*u3)*t7 + (k3u3-k2*u2)*t8 - (k2*u3+k3*u2)*t9 1014 ek = i0*(v1*t2+v2*t3+v3*t4) + (i3v3-i1*v1)*t5 - (i1*v2+i2*v1)*t6 1015 & - (i3*v1+i1*v3)*t7 + (i3v3-i2*v2)*t8 - (i2*v3+i3*v2)*t9 1016c 1017c get the M-Q, D-Q and Q-Dind energies, if necessary 1018c 1019 if (kquad) then 1020 k4k9 = k4 - k9 1021 k7k9 = k7 - k9 1022 i3k6 = 2.0d0 * i3 * k6 1023 i3k8 = 2.0d0 * i3 * k8 1024 eik = eik + i0*k4k9*t5 + 2.0d0*i0*k5*t6 + 2.0d0*i0*k6*t7 1025 & + i0*k7k9*t8 + 2.0d0*i0*k8*t9 + (i3k6-i1*k4k9)*t11 1026 & + (i3k8-i2*k4k9-2.0d0*i1*k5)*t12 1027 & - (2.0d0*i1*k6+i3*k4k9)*t13 1028 & + (i3k6-i1*k7k9-2.0d0*i2*k5)*t14 1029 & - 2.0d0*(i1*k8+i2*k6+i3*k5)*t15 1030 & + (i3k8-i2*k7k9)*t17 1031 & - (2.0d0*i2*k8+i3*k7k9)*t18 1032 k6u3 = k6 * u3 1033 k8u3 = k8 * u3 1034 ei = ei + (-k4k9*u1+2.0d0*k6u3)*t11 1035 & + (-k4k9*u2+2.0d0*(k8u3-k5*u1))*t12 1036 & + (-k4k9*u3-2.0d0*k6*u1)*t13 1037 & + (-k7k9*u1+2.0d0*(k6u3-k5*u2))*t14 1038 & - 2.0d0*(k5*u3+k6*u2+k8*u1)*t15 1039 & + (-k7k9*u2+2.0d0*k8u3)*t17 1040 & + (-k7k9*u3-2.0d0*k8*u2)*t18 1041 end if 1042 if (iquad) then 1043 i4i9 = i4 - i9 1044 i7i9 = i7 - i9 1045 i6k3 = -2.0d0 * i6 * k3 1046 i8k3 = -2.0d0 * i8 * k3 1047 eik = eik + i4i9*k0*t5 + 2.0d0*i5*k0*t6 + 2.0d0*i6*k0*t7 1048 & + i7i9*k0*t8 + 2.0d0*i8*k0*t9 + (i4i9*k1+i6k3)*t11 1049 & + (i8k3+i4i9*k2+2.0d0*i5*k1)*t12 1050 & + (2.0d0*i6*k1+i4i9*k3)*t13 1051 & + (i6k3+i7i9*k1+2.0d0*i5*k2)*t14 1052 & + 2.0d0*(i8*k1+i6*k2+i5*k3)*t15 1053 & + (i8k3+i7i9*k2)*t17 1054 & + (2.0d0*i8*k2+i7i9*k3)*t18 1055 i6v3 = i6 * v3 1056 i8v3 = i8 * v3 1057 ek = ek + (i4i9*v1-2.0d0*i6v3)*t11 1058 & + (i4i9*v2+2.0d0*(i5*v1-i8v3))*t12 1059 & + (i4i9*v3+2.0d0*i6*v1)*t13 1060 & + (i7i9*v1+2.0d0*(i5*v2-i6v3))*t14 1061 & + 2.0d0*(i5*v3+i6*v2+i8*v1)*t15 1062 & + (i7i9*v2-2.0d0*i8v3)*t17 1063 & + (i7i9*v3+2.0d0*i8*v2)*t18 1064 end if 1065c 1066c get the Q-Q part of the multipole interaction energy 1067c 1068 if (iquad .and. kquad) then 1069 ik49 = i4*k9 + i9*k4 1070 ik59 = i5*k9 + i9*k5 1071 ik69 = i6*k9 + i9*k6 1072 ik79 = i7*k9 + i9*k7 1073 ik89 = i8*k9 + i9*k8 1074 ik68 = 2.0d0 * (i6*k8 + i8*k6) 1075 i6k6 = i6 * k6 1076 i8k8 = i8 * k8 1077 i9k9 = i9 * k9 1078 eik = eik + (i4*k4-ik49+i9k9-4.0d0*i6k6)*t21 1079 & + 2.0d0*(i5*k4+i4*k5-ik68-ik59)*t22 1080 & + 2.0d0*(i4*k6+i6*k4-ik69)*t23 1081 & + (i4*k7+i7*k4-ik49-ik79+2.0d0*i9k9 1082 & +4.0d0*(i5*k5-i6k6-i8k8))*t24 1083 & + 2.0d0*(i4*k8+i8*k4+2.0d0*(i5*k6+i6*k5)-ik89)*t25 1084 & + 2.0d0*(i5*k7+i7*k5-ik68-ik59)*t27 1085 & + 2.0d0*(i6*k7+i7*k6+2.0d0*(i5*k8+i8*k5)-ik69)*t28 1086 & + (i7*k7-ik79+i9k9-4.0d0*i8k8)*t31 1087 & + 2.0d0*(i7*k8+i8*k7-ik89)*t32 1088 end if 1089c 1090c apply exponential damping to the polarization term 1091c 1092 damp = pdamp(ii) * pdamp(kk) 1093 if (damp .ne. 0.0d0) then 1094 damp = -pgamma * (r/damp)**3 1095 if (damp .gt. -50.0d0) then 1096 damp = 1.0d0 - exp(damp) 1097 ei = ei * damp 1098 ek = ek * damp 1099 end if 1100 end if 1101c 1102c final polarization energy is half of value computed above 1103c 1104 ei = 0.5d0 * ei 1105 ek = 0.5d0 * ek 1106 return 1107 end 1108c 1109c 1110c ################################################################# 1111c ## ## 1112c ## subroutine empik1 -- mpole & polarization pair gradient ## 1113c ## ## 1114c ################################################################# 1115c 1116c 1117c "empik1" computes the permanent multipole and induced dipole 1118c energies and derivatives between a pair of multipole sites 1119c 1120c 1121 subroutine empik1 (ii,kk,xr,yr,zr,r,r2,rpi,rpk,indi,indk, 1122 & eik,ei,ek,dm,dmi,dmk,dp,dpi,dpk,utu) 1123 implicit none 1124 include 'sizes.i' 1125 include 'atoms.i' 1126 include 'mpole.i' 1127 include 'polar.i' 1128 include 'polpot.i' 1129 include 'units.i' 1130 integer i,j,ii,kk 1131 real*8 eik,ei,ek 1132 real*8 damp,ddamp,de,term 1133 real*8 xr,yr,zr,r,r2,r4 1134 real*8 rpi(13),rpk(13) 1135 real*8 indi(3),indk(3) 1136 real*8 dm(3),dp(3),utu 1137 real*8 dmi(3,3),dmk(3,3) 1138 real*8 dpi(3,3),dpk(3,3) 1139 real*8 rr1,rr2,rr3,rr5,rr7,rr9,rr11 1140 real*8 xx,yy,xyz,xry,yrx,zrx,zry 1141 real*8 xr5,yr5,xr7,yr7 1142 real*8 xr9,xyr9,yr9,xr9x7,yr9y7 1143 real*8 rr53,xr73,yr73 1144 real*8 x2,y2,z2,x4,y4,x2y2,x2r2,y2r2 1145 real*8 xrr11,yrr11,zrr11,xyzrr11 1146 real*8 r445,r4225,nnn1,nnn2 1147 real*8 n945x4,n945y4,n945x2y2 1148 real*8 i3k3,i6k6,i8k8,i9k9 1149 real*8 ik49,ik59,ik69,ik79,ik89,ik68 1150 real*8 k4k9,k7k9,i4i9,i7i9 1151 real*8 k3u3,i3k6,i3k8,k6u3,k8u3 1152 real*8 i3v3,i6k3,i8k3,i6v3,i8v3 1153 real*8 i0,i1,i2,i3,i4,i5,i6,i7,i8,i9 1154 real*8 k0,k1,k2,k3,k4,k5,k6,k7,k8,k9 1155 real*8 u1,u2,u3,v1,v2,v3 1156 real*8 di1,di2,di3,di4,di5,di6,di7,di8,di9 1157 real*8 dk1,dk2,dk3,dk4,dk5,dk6,dk7,dk8,dk9 1158 real*8 w1,w2,w3,w4,w5,w6,w7,w8,w9,w10,w11,w12 1159 real*8 t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13 1160 real*8 t14,t15,t16,t17,t18,t19,t21,t22,t23,t24,t25 1161 real*8 t26,t27,t28,t29,t31,t32,t33,t36,t37,t38,t39 1162 real*8 t40,t41,t42,t43,t44,t46,t47,t48,t51,t52,t53 1163 logical iquad,kquad 1164c 1165c 1166c zero out the energy and first derivative components 1167c 1168 eik = 0.0d0 1169 ei = 0.0d0 1170 ek = 0.0d0 1171 utu = 0.0d0 1172 do i = 1, 3 1173 dm(i) = 0.0d0 1174 dp(i) = 0.0d0 1175 do j = 1, 3 1176 dmi(j,i) = 0.0d0 1177 dmk(j,i) = 0.0d0 1178 dpi(j,i) = 0.0d0 1179 dpk(j,i) = 0.0d0 1180 end do 1181 end do 1182c 1183c check for presence of quadrupole components at either site 1184c 1185 iquad = (polsiz(ii) .ge. 13) 1186 kquad = (polsiz(kk) .ge. 13) 1187c 1188c set permanent and induced multipoles for first site 1189c 1190 i0 = rpi(1) 1191 i1 = rpi(2) 1192 i2 = rpi(3) 1193 i3 = rpi(4) 1194 i4 = rpi(5) 1195 i5 = rpi(6) 1196 i6 = rpi(7) 1197 i7 = rpi(9) 1198 i8 = rpi(10) 1199 i9 = rpi(13) 1200 u1 = indi(1) 1201 u2 = indi(2) 1202 u3 = indi(3) 1203c 1204c set permanent and induced multipoles for second site 1205c 1206 k0 = rpk(1) 1207 k1 = rpk(2) 1208 k2 = rpk(3) 1209 k3 = rpk(4) 1210 k4 = rpk(5) 1211 k5 = rpk(6) 1212 k6 = rpk(7) 1213 k7 = rpk(9) 1214 k8 = rpk(10) 1215 k9 = rpk(13) 1216 v1 = indk(1) 1217 v2 = indk(2) 1218 v3 = indk(3) 1219c 1220c compute the zeroth order T2 matrix element 1221c 1222 rr1 = 1.0d0 / r 1223 t1 = rr1 1224c 1225c compute the first order T2 matrix elements 1226c 1227 rr2 = rr1 * rr1 1228 rr3 = rr1 * rr2 1229 t2 = -xr * rr3 1230 t3 = -yr * rr3 1231 t4 = -zr * rr3 1232c 1233c compute the second order T2 matrix elements 1234c 1235 rr5 = 3.0d0 * rr3 * rr2 1236 xr5 = xr * rr5 1237 yr5 = yr * rr5 1238 t5 = -rr3 + xr5*xr 1239 t6 = xr5 * yr 1240 t7 = xr5 * zr 1241 t8 = -rr3 + yr5*yr 1242 t9 = yr5 * zr 1243 t10 = -t5 - t8 1244c 1245c compute the third order T2 matrix elements 1246c 1247 rr7 = 5.0d0 * rr5 * rr2 1248 xx = xr * xr 1249 yy = yr * yr 1250 xyz = xr * yr * zr 1251 xr7 = xx * rr7 1252 yr7 = yy * rr7 1253 rr53 = 3.0d0 * rr5 1254 t11 = xr * (rr53-xr7) 1255 t12 = yr * (rr5-xr7) 1256 t13 = zr * (rr5-xr7) 1257 t14 = xr * (rr5-yr7) 1258 t15 = -xyz * rr7 1259 t16 = -t11 - t14 1260 t17 = yr * (rr53-yr7) 1261 t18 = zr * (rr5-yr7) 1262 t19 = -t12 - t17 1263c 1264c compute the fourth order T2 matrix elements 1265c 1266 if (iquad .or. kquad) then 1267 rr9 = 7.0d0 * rr7 * rr2 1268 if (xr .eq. 0.0d0) then 1269 yrx = 0.0d0 1270 zrx = 0.0d0 1271 else 1272 yrx = yr / xr 1273 zrx = zr / xr 1274 end if 1275 if (yr .eq. 0.0d0) then 1276 xry = 0.0d0 1277 zry = 0.0d0 1278 else 1279 xry = xr / yr 1280 zry = zr / yr 1281 end if 1282 xr9 = xx * xx * rr9 1283 xyr9 = xx * yy * rr9 1284 yr9 = yy * yy * rr9 1285 xr73 = 3.0d0 * xr7 1286 yr73 = 3.0d0 * yr7 1287 xr9x7 = xr9 - xr73 1288 yr9y7 = yr9 - yr73 1289 t21 = xr9x7 - xr73 + rr53 1290 t22 = yrx * xr9x7 1291 t23 = zrx * xr9x7 1292 t24 = xyr9 - xr7 - yr7 + rr5 1293 t25 = zry * (xyr9-yr7) 1294 t26 = -t21 - t24 1295 t27 = xry * yr9y7 1296 t28 = zrx * (xyr9-xr7) 1297 t29 = -t22 - t27 1298 t31 = yr9y7 - yr73 + rr53 1299 t32 = zry * yr9y7 1300 t33 = -t24 - t31 1301c 1302c compute the fifth order T2 matrix elements 1303c 1304 r4 = r2 * r2 1305 rr5 = rr2 * rr3 1306 rr7 = rr2 * rr5 1307 rr9 = rr2 * rr7 1308 rr11 = rr2 * rr9 1309 x2 = xr * xr 1310 y2 = yr * yr 1311 z2 = zr * zr 1312 x4 = x2 * x2 1313 y4 = y2 * y2 1314 x2r2 = x2 * r2 1315 y2r2 = y2 * r2 1316 x2y2 = x2 * y2 1317 xrr11 = xr * rr11 1318 yrr11 = yr * rr11 1319 zrr11 = zr * rr11 1320 xyzrr11 = 315.0d0 * xyz * rr11 1321 r445 = 45.0d0 * r4 1322 n945x4 = -945.0d0 * x4 1323 n945y4 = -945.0d0 * y4 1324 n945x2y2 = -945.0d0 * x2y2 1325 nnn1 = n945x4 + 630.0d0*x2r2 - r445 1326 nnn2 = n945y4 + 630.0d0*y2r2 - r445 1327 r4225 = 225.0d0 * r4 1328 t36 = (n945x4+1050.0d0*x2r2-r4225) * xrr11 1329 t37 = nnn1 * yrr11 1330 t38 = nnn1 * zrr11 1331 t39 = (n945x2y2+105.0d0*x2r2+315.0d0*y2r2-r445) * xrr11 1332 t40 = (r2-3.0d0*x2) * xyzrr11 1333 t41 = -t36 - t39 1334 t42 = (n945x2y2+105.0d0*y2r2+315.0d0*x2r2-r445) * yrr11 1335 t43 = (n945x2y2+105.0d0*(x2r2+y2r2)-15.0d0*r4) * zrr11 1336 t44 = -t37 - t42 1337 t46 = nnn2 * xrr11 1338 t47 = (r2-3.0d0*y2) * xyzrr11 1339 t48 = -t39 - t46 1340 t51 = (n945y4+1050.0d0*y2r2-r4225) * yrr11 1341 t52 = nnn2 * zrr11 1342 t53 = -t42 - t51 1343 end if 1344c 1345c get the M-M, M-D and D-D parts of the multipole energy 1346c 1347 i3k3 = i3 * k3 1348 w1 = i0 * k0 1349 w2 = i0*k1 - i1*k0 1350 w3 = i0*k2 - i2*k0 1351 w4 = i0*k3 - i3*k0 1352 w5 = i3k3 - i1*k1 1353 w6 = -i1*k2 - i2*k1 1354 w7 = -i1*k3 - i3*k1 1355 w8 = i3k3 - i2*k2 1356 w9 = -i2*k3 - i3*k2 1357 eik = w1*t1 + w2*t2 + w3*t3 + w4*t4 + w5*t5 1358 & + w6*t6 + w7*t7 + w8*t8 + w9*t9 1359 dm(1) = w1*t2 + w2*t5 + w3*t6 + w4*t7 + w5*t11 1360 & + w6*t12 + w7*t13 + w8*t14 + w9*t15 1361 dm(2) = w1*t3 + w2*t6 + w3*t8 + w4*t9 + w5*t12 1362 & + w6*t14 + w7*t15 + w8*t17 + w9*t18 1363 dm(3) = w1*t4 + w2*t7 + w3*t9 + w4*t10 + w5*t13 1364 & + w6*t15 + w7*t16 + w8*t18 + w9*t19 1365c 1366c get the M-Q and D-Q parts of the multipole energy 1367c 1368 if (kquad) then 1369 k4k9 = k4 - k9 1370 k7k9 = k7 - k9 1371 i3k6 = 2.0d0 * i3 * k6 1372 i3k8 = 2.0d0 * i3 * k8 1373 w1 = i0 * k4k9 1374 w2 = 2.0d0 * i0 * k5 1375 w3 = 2.0d0 * i0 * k6 1376 w4 = i0 * k7k9 1377 w5 = 2.0d0 * i0 * k8 1378 w6 = i3k6 - i1*k4k9 1379 w7 = i3k8 - i2*k4k9 - 2.0d0*i1*k5 1380 w8 = -2.0d0*i1*k6 - i3*k4k9 1381 w9 = i3k6 - i1*k7k9 - 2.0d0*i2*k5 1382 w10 = -2.0d0 * (i1*k8 + i2*k6 + i3*k5) 1383 w11 = i3k8 - i2*k7k9 1384 w12 = -2.0d0*i2*k8 - i3*k7k9 1385 eik = eik + w1*t5 + w2*t6 + w3*t7 + w4*t8 1386 & + w5*t9 + w6*t11 + w7*t12 + w8*t13 1387 & + w9*t14 + w10*t15 + w11*t17 + w12*t18 1388 dm(1) = dm(1) + w1*t11 + w2*t12 + w3*t13 + w4*t14 1389 & + w5*t15 + w6*t21 + w7*t22 + w8*t23 1390 & + w9*t24 + w10*t25 + w11*t27 + w12*t28 1391 dm(2) = dm(2) + w1*t12 + w2*t14 + w3*t15 + w4*t17 1392 & + w5*t18 + w6*t22 + w7*t24 + w8*t25 1393 & + w9*t27 + w10*t28 + w11*t31 + w12*t32 1394 dm(3) = dm(3) + w1*t13 + w2*t15 + w3*t16 + w4*t18 1395 & + w5*t19 + w6*t23 + w7*t25 + w8*t26 1396 & + w9*t28 + w10*t29 + w11*t32 + w12*t33 1397 end if 1398c 1399c get the M-Q and D-Q parts of the multipole energy 1400c 1401 if (iquad) then 1402 i4i9 = i4 - i9 1403 i7i9 = i7 - i9 1404 i6k3 = -2.0d0 * i6 * k3 1405 i8k3 = -2.0d0 * i8 * k3 1406 w1 = i4i9 * k0 1407 w2 = 2.0d0 * i5 * k0 1408 w3 = 2.0d0 * i6 * k0 1409 w4 = i7i9 * k0 1410 w5 = 2.0d0 * i8 * k0 1411 w6 = i4i9*k1 + i6k3 1412 w7 = i8k3 + i4i9*k2 + 2.0d0*i5*k1 1413 w8 = 2.0d0*i6*k1 + i4i9*k3 1414 w9 = i6k3 + i7i9*k1 + 2.0d0*i5*k2 1415 w10 = 2.0d0 * (i8*k1 + i6*k2 + i5*k3) 1416 w11 = i8k3 + i7i9*k2 1417 w12 = 2.0d0*i8*k2 + i7i9*k3 1418 eik = eik + w1*t5 + w2*t6 + w3*t7 + w4*t8 1419 & + w5*t9 + w6*t11 + w7*t12 + w8*t13 1420 & + w9*t14 + w10*t15 + w11*t17 + w12*t18 1421 dm(1) = dm(1) + w1*t11 + w2*t12 + w3*t13 + w4*t14 1422 & + w5*t15 + w6*t21 + w7*t22 + w8*t23 1423 & + w9*t24 + w10*t25 + w11*t27 + w12*t28 1424 dm(2) = dm(2) + w1*t12 + w2*t14 + w3*t15 + w4*t17 1425 & + w5*t18 + w6*t22 + w7*t24 + w8*t25 1426 & + w9*t27 + w10*t28 + w11*t31 + w12*t32 1427 dm(3) = dm(3) + w1*t13 + w2*t15 + w3*t16 + w4*t18 1428 & + w5*t19 + w6*t23 + w7*t25 + w8*t26 1429 & + w9*t28 + w10*t29 + w11*t32 + w12*t33 1430 end if 1431c 1432c get the Q-Q part of the multipole interaction energy 1433c 1434 if (iquad .and. kquad) then 1435 ik49 = i4*k9 + i9*k4 1436 ik59 = i5*k9 + i9*k5 1437 ik69 = i6*k9 + i9*k6 1438 ik79 = i7*k9 + i9*k7 1439 ik89 = i8*k9 + i9*k8 1440 ik68 = 2.0d0 * (i6*k8 + i8*k6) 1441 i6k6 = i6 * k6 1442 i8k8 = i8 * k8 1443 i9k9 = i9 * k9 1444 w1 = i4*k4 - ik49 + i9k9 - 4.0d0*i6k6 1445 w2 = 2.0d0 * (i5*k4 + i4*k5 - ik68 - ik59) 1446 w3 = 2.0d0 * (i4*k6 + i6*k4 - ik69) 1447 w4 = i4*k7 + i7*k4 - ik49 - ik79 + 2.0d0*i9k9 1448 & +4.0d0*(i5*k5-i6k6-i8k8) 1449 w5 = 2.0d0 * (i4*k8 + i8*k4 + 2.0d0*(i5*k6+i6*k5) - ik89) 1450 w6 = 2.0d0 * (i5*k7 + i7*k5 - ik68 - ik59) 1451 w7 = 2.0d0 * (i6*k7 + i7*k6 + 2.0d0*(i5*k8+i8*k5) - ik69) 1452 w8 = i7*k7 - ik79 + i9k9 - 4.0d0*i8k8 1453 w9 = 2.0d0 * (i7*k8 + i8*k7 - ik89) 1454 eik = eik + w1*t21 + w2*t22 + w3*t23 + w4*t24 1455 & + w5*t25 + w6*t27 + w7*t28 + w8*t31 + w9*t32 1456 dm(1) = dm(1) + w1*t36 + w2*t37 + w3*t38 + w4*t39 1457 & + w5*t40 + w6*t42 + w7*t43 + w8*t46 + w9*t47 1458 dm(2) = dm(2) + w1*t37 + w2*t39 + w3*t40 + w4*t42 1459 & + w5*t43 + w6*t46 + w7*t47 + w8*t51 + w9*t52 1460 dm(3) = dm(3) + w1*t38 + w2*t40 + w3*t41 + w4*t43 1461 & + w5*t44 + w6*t47 + w7*t48 + w8*t52 + w9*t53 1462 end if 1463c 1464c get the (dM2/dx)*T*M1 terms for dipoles at both sites 1465c 1466 do i = 1, 3 1467 do j = 1, 3 1468 di1 = dpole(2,j,i,ii) 1469 di2 = dpole(3,j,i,ii) 1470 di3 = dpole(4,j,i,ii) 1471 dk1 = dpole(2,j,i,kk) 1472 dk2 = dpole(3,j,i,kk) 1473 dk3 = dpole(4,j,i,kk) 1474 w1 = k3 * di3 1475 dmi(j,i) = -k0*di1*t2 - k0*di2*t3 - k0*di3*t4 1476 & + (w1-k1*di1)*t5 - (k2*di1+k1*di2)*t6 1477 & - (k1*di3+k3*di1)*t7 + (w1-k2*di2)*t8 1478 & - (k2*di3+k3*di2)*t9 1479 w1 = dk3 * i3 1480 dmk(j,i) = dk1*i0*t2 + dk2*i0*t3 + dk3*i0*t4 1481 & + (w1-dk1*i1)*t5 - (dk2*i1+dk1*i2)*t6 1482 & - (dk1*i3+dk3*i1)*t7 + (w1-dk2*i2)*t8 1483 & - (dk2*i3+dk3*i2)*t9 1484 w1 = v3 * di3 1485 dpi(j,i) = (w1-v1*di1)*t5 - (v2*di1+v1*di2)*t6 1486 & - (v1*di3+v3*di1)*t7 + (w1-v2*di2)*t8 1487 & - (v3*di2+v2*di3)*t9 1488 w1 = u3 * dk3 1489 dpk(j,i) = (w1-u1*dk1)*t5 - (u2*dk1+u1*dk2)*t6 1490 & - (u1*dk3+u3*dk1)*t7 + (w1-u2*dk2)*t8 1491 & - (u3*dk2+u2*dk3)*t9 1492c 1493c get the (dM2/dx)*T*M1 terms for quadrupole at first site 1494c 1495 if (iquad) then 1496 di4 = dpole(5,j,i,ii) 1497 di5 = dpole(6,j,i,ii) 1498 di6 = dpole(7,j,i,ii) 1499 di7 = dpole(9,j,i,ii) 1500 di8 = dpole(10,j,i,ii) 1501 di9 = dpole(13,j,i,ii) 1502 w1 = i9 - i4 1503 w2 = i9 - i7 1504 w3 = i6 * dk3 1505 w4 = i8 * dk3 1506 dmk(j,i) = dmk(j,i) - (w1*dk1+2.0d0*w3)*t11 1507 & - (w1*dk2+2.0d0*(w4-dk1*i5))*t12 1508 & - (w1*dk3-2.0d0*dk1*i6)*t13 1509 & - (w2*dk1+2.0d0*(w3-dk2*i5))*t14 1510 & + 2.0d0*(dk3*i5+dk2*i6+dk1*i8)*t15 1511 & - (w2*dk2+2.0d0*w4)*t17 1512 & - (w2*dk3-2.0d0*dk2*i8)*t18 1513 w1 = di9 - di4 1514 w2 = di9 - di7 1515 w3 = di6 * k3 1516 w4 = di8 * k3 1517 dmi(j,i) = dmi(j,i) - k0*(w1*t5+w2*t8) 1518 & + 2.0d0*k0*(di5*t6+di6*t7+di8*t9) 1519 & - (w1*k1+2.0d0*w3)*t11 1520 & - (w1*k2+2.0d0*(w4-k1*di5))*t12 1521 & - (w1*k3-2.0d0*k1*di6)*t13 1522 & - (w2*k1+2.0d0*(w3-k2*di5))*t14 1523 & + 2.0d0*(k3*di5+k2*di6+k1*di8)*t15 1524 & - (w2*k2+2.0d0*w4)*t17 1525 & - (w2*k3-2.0d0*k2*di8)*t18 1526c w1 = di9 - di4 1527c w2 = di9 - di7 1528 w3 = di6 * v3 1529 w4 = di8 * v3 1530 dpi(j,i) = dpi(j,i) - (w1*v1+2.0d0*w3)*t11 1531 & - (w1*v2+2.0d0*(w4-v1*di5))*t12 1532 & - (w1*v3-2.0d0*v1*di6)*t13 1533 & - (w2*v1+2.0d0*(w3-v2*di5))*t14 1534 & + 2.0d0*(v3*di5+v2*di6+v1*di8)*t15 1535 & - (w2*v2+2.0d0*w4)*t17 1536 & - (w2*v3-2.0d0*v2*di8)*t18 1537 end if 1538c 1539c get the (dM2/dx)*T*M1 terms for quadrupole at second site 1540c 1541 if (kquad) then 1542 dk4 = dpole(5,j,i,kk) 1543 dk5 = dpole(6,j,i,kk) 1544 dk6 = dpole(7,j,i,kk) 1545 dk7 = dpole(9,j,i,kk) 1546 dk8 = dpole(10,j,i,kk) 1547 dk9 = dpole(13,j,i,kk) 1548 w1 = k9 - k4 1549 w2 = k9 - k7 1550 w3 = k6 * di3 1551 w4 = k8 * di3 1552 dmi(j,i) = dmi(j,i) 1553 & + (w1*di1+2.0d0*w3)*t11 1554 & + (w1*di2+2.0d0*(w4-k5*di1))*t12 1555 & + (w1*di3-2.0d0*k6*di1)*t13 1556 & + (w2*di1+2.0d0*(w3-k5*di2))*t14 1557 & - 2.0d0*(k5*di3+k6*di2+k8*di1)*t15 1558 & + (w2*di2+2.0d0*w4)*t17 1559 & + (w2*di3-2.0d0*k8*di2)*t18 1560 w1 = dk9 - dk4 1561 w2 = dk9 - dk7 1562 w3 = dk6 * i3 1563 w4 = dk8 * i3 1564 dmk(j,i) = dmk(j,i) - i0*(w1*t5+w2*t8) 1565 & + 2.0d0*i0*(dk5*t6+dk6*t7+dk8*t9) 1566 & + (w1*i1+2.0d0*w3)*t11 1567 & + (w1*i2+2.0d0*(w4-dk5*i1))*t12 1568 & + (w1*i3-2.0d0*dk6*i1)*t13 1569 & + (w2*i1+2.0d0*(w3-dk5*i2))*t14 1570 & - 2.0d0*(dk8*i1+dk5*i3+dk6*i2)*t15 1571 & + (w2*i2+2.0d0*w4)*t17 1572 & + (w2*i3-2.0d0*dk8*i2)*t18 1573c w1 = dk9 - dk4 1574c w2 = dk9 - dk7 1575 w3 = dk6 * u3 1576 w4 = dk8 * u3 1577 dpk(j,i) = dpk(j,i) + (w1*u1+2.0d0*w3)*t11 1578 & + (w1*u2+2.0d0*(w4-u1*dk5))*t12 1579 & + (w1*u3-2.0d0*u1*dk6)*t13 1580 & + (w2*u1+2.0d0*(w3-u2*dk5))*t14 1581 & - 2.0d0*(u1*dk8+u2*dk6+u3*dk5)*t15 1582 & + (w2*u2+2.0d0*w4)*t17 1583 & + (w2*u3-2.0d0*u2*dk8)*t18 1584 end if 1585c 1586c get the (dM2/dx)*T*M1 terms for quadrupoles at both sites 1587c 1588 if (iquad .and. kquad) then 1589 w1 = di9 - di4 1590 w2 = di9 - di7 1591 w3 = di5*k9 + 2.0d0*(di6*k8+di8*k6) 1592 w4 = di6 * k6 1593 w5 = di6 * k9 1594 w6 = di8 * k8 1595 w7 = di8 * k9 1596 dmi(j,i) = dmi(j,i) + (w1*(k9-k4)-4.0d0*w4)*t21 1597 & + 2.0d0*(di5*k4-w1*k5-w3)*t22 1598 & + 2.0d0*(di6*k4-w1*k6-w5)*t23 1599 & + (4.0d0*(di5*k5-w4-w6)-w1*k7-w2*k4 1600 & +(2.0d0*di9-di4-di7)*k9)*t24 1601 & + 2.0d0*(di8*k4-w1*k8-w7 1602 & +2.0d0*(di6*k5+di5*k6))*t25 1603 & + 2.0d0*(di5*k7-w2*k5-w3)*t27 1604 & + 2.0d0*(di6*k7-w2*k6-w5 1605 & +2.0d0*(di8*k5+di5*k8))*t28 1606 & + (w2*(k9-k7)-4.0d0*w6)*t31 1607 & + 2.0d0*(di8*k7-w2*k8-w7)*t32 1608 w1 = dk9 - dk4 1609 w2 = dk9 - dk7 1610 w3 = dk5*i9 + 2.0d0*(dk6*i8+dk8*i6) 1611 w4 = dk6 * i6 1612 w5 = dk6 * i9 1613 w6 = dk8 * i8 1614 w7 = dk8 * i9 1615 dmk(j,i) = dmk(j,i) + (w1*(i9-i4)-4.0d0*w4)*t21 1616 & + 2.0d0*(dk5*i4-w1*i5-w3)*t22 1617 & + 2.0d0*(dk6*i4-w1*i6-w5)*t23 1618 & + (4.0d0*(dk5*i5-w4-w6)-w1*i7-w2*i4 1619 & +(2.0d0*dk9-dk4-dk7)*i9)*t24 1620 & + 2.0d0*(dk8*i4-w1*i8-w7 1621 & +2.0d0*(dk6*i5+dk5*i6))*t25 1622 & + 2.0d0*(dk5*i7-w2*i5-w3)*t27 1623 & + 2.0d0*(dk6*i7-w2*i6-w5 1624 & +2.0d0*(dk8*i5+dk5*i8))*t28 1625 & + (w2*(i9-i7)-4.0d0*w6)*t31 1626 & + 2.0d0*(dk8*i7-w2*i8-w7)*t32 1627 end if 1628 end do 1629 end do 1630c 1631c get indD-M, indD-D and indD-Q polarization energy and gradients 1632c 1633 k3u3 = k3 * u3 1634 w1 = k0 * u1 1635 w2 = k0 * u2 1636 w3 = k0 * u3 1637 w4 = k3u3 - k1*u1 1638 w5 = k1*u2 + k2*u1 1639 w6 = k3*u1 + k1*u3 1640 w7 = k3u3 - k2*u2 1641 w8 = k2*u3 + k3*u2 1642 ei = -w1*t2 - w2*t3 - w3*t4 + w4*t5 - w5*t6 1643 & - w6*t7 + w7*t8 - w8*t9 1644 dp(1) = -w1*t5 - w2*t6 - w3*t7 + w4*t11 - w5*t12 1645 & - w6*t13 + w7*t14 - w8*t15 1646 dp(2) = -w1*t6 - w2*t8 - w3*t9 + w4*t12 - w5*t14 1647 & - w6*t15 + w7*t17 - w8*t18 1648 dp(3) = -w1*t7 - w2*t9 - w3*t10 + w4*t13 - w5*t15 1649 & - w6*t16 + w7*t18 - w8*t19 1650 if (kquad) then 1651 k6u3 = k6 * u3 1652 k8u3 = k8 * u3 1653 w1 = -k4k9*u1 + 2.0d0*k6u3 1654 w2 = -k4k9*u2 + 2.0d0*(k8u3-k5*u1) 1655 w3 = -k4k9*u3 - 2.0d0*k6*u1 1656 w4 = -k7k9*u1 + 2.0d0*(k6u3-k5*u2) 1657 w5 = 2.0d0 * (k5*u3 + k6*u2 + k8*u1) 1658 w6 = -k7k9*u2 + 2.0d0*k8u3 1659 w7 = -k7k9*u3 - 2.0d0*k8*u2 1660 ei = ei + w1*t11 + w2*t12 + w3*t13 + w4*t14 1661 & - w5*t15 + w6*t17 + w7*t18 1662 dp(1) = dp(1) + w1*t21 + w2*t22 + w3*t23 1663 & + w4*t24 - w5*t25 + w6*t27 + w7*t28 1664 dp(2) = dp(2) + w1*t22 + w2*t24 + w3*t25 1665 & + w4*t27 - w5*t28 + w6*t31 + w7*t32 1666 dp(3) = dp(3) + w1*t23 + w2*t25 + w3*t26 1667 & + w4*t28 - w5*t29 + w6*t32 + w7*t33 1668 end if 1669c 1670c get M-indD, D-indD and Q-indD polarization energy and gradients 1671c 1672 i3v3 = i3 * v3 1673 w1 = i0 * v1 1674 w2 = i0 * v2 1675 w3 = i0 * v3 1676 w4 = i3v3 - i1*v1 1677 w5 = i1*v2 + i2*v1 1678 w6 = i3*v1 + i1*v3 1679 w7 = i3v3 - i2*v2 1680 w8 = i2*v3 + i3*v2 1681 ek = w1*t2 + w2*t3 + w3*t4 + w4*t5 - w5*t6 1682 & - w6*t7 + w7*t8 - w8*t9 1683 dp(1) = dp(1) + w1*t5 + w2*t6 + w3*t7 + w4*t11 - w5*t12 1684 & - w6*t13 + w7*t14 - w8*t15 1685 dp(2) = dp(2) + w1*t6 + w2*t8 + w3*t9 + w4*t12 - w5*t14 1686 & - w6*t15 + w7*t17 - w8*t18 1687 dp(3) = dp(3) + w1*t7 + w2*t9 + w3*t10 + w4*t13 - w5*t15 1688 & - w6*t16 + w7*t18 - w8*t19 1689 if (iquad) then 1690 i6v3 = i6 * v3 1691 i8v3 = i8 * v3 1692 w1 = i4i9*v1 - 2.0d0*i6v3 1693 w2 = i4i9*v2 + 2.0d0*(i5*v1-i8v3) 1694 w3 = i4i9*v3 + 2.0d0*i6*v1 1695 w4 = i7i9*v1 + 2.0d0*(i5*v2-i6v3) 1696 w5 = 2.0d0 * (i5*v3 + i6*v2 + i8*v1) 1697 w6 = i7i9*v2 - 2.0d0*i8v3 1698 w7 = i7i9*v3 + 2.0d0*i8*v2 1699 ek = ek + w1*t11 + w2*t12 + w3*t13 + w4*t14 1700 & + w5*t15 + w6*t17 + w7*t18 1701 dp(1) = dp(1) + w1*t21 + w2*t22 + w3*t23 1702 & + w4*t24 + w5*t25 + w6*t27 + w7*t28 1703 dp(2) = dp(2) + w1*t22 + w2*t24 + w3*t25 1704 & + w4*t27 + w5*t28 + w6*t31 + w7*t32 1705 dp(3) = dp(3) + w1*t23 + w2*t25 + w3*t26 1706 & + w4*t28 + w5*t29 + w6*t32 + w7*t33 1707 end if 1708c 1709c compute the mutual polarization induced dipole gradient terms 1710c 1711 if (poltyp .eq. 'MUTUAL') then 1712 w1 = -u1*v1 + u3*v3 1713 w2 = -u2*v1 - u1*v2 1714 w3 = -u3*v1 - u1*v3 1715 w4 = -u2*v2 + u3*v3 1716 w5 = -u3*v2 - u2*v3 1717 utu = w1*t5 + w2*t6 + w3*t7 + w4*t8 + w5*t9 1718 dp(1) = dp(1) + w1*t11 + w2*t12 + w3*t13 + w4*t14 + w5*t15 1719 dp(2) = dp(2) + w1*t12 + w2*t14 + w3*t15 + w4*t17 + w5*t18 1720 dp(3) = dp(3) + w1*t13 + w2*t15 + w3*t16 + w4*t18 + w5*t19 1721 end if 1722c 1723c apply a damping factor to polarization energy and derivatives 1724c 1725 damp = pdamp(ii) * pdamp(kk) 1726 if (damp .ne. 0.0d0) then 1727 term = -pgamma * (r/damp)**3 1728 if (term .gt. -50.0d0) then 1729 term = exp(term) 1730 ddamp = (3.0d0*pgamma*r2/damp**3) * term 1731 damp = 1.0d0 - term 1732 de = (ei+ek+utu) * ddamp 1733 ei = ei * damp 1734 ek = ek * damp 1735 dp(1) = dp(1)*damp + de*(xr/r) 1736 dp(2) = dp(2)*damp + de*(yr/r) 1737 dp(3) = dp(3)*damp + de*(zr/r) 1738 do i = 1, 3 1739 do j = 1, 3 1740 dpi(j,i) = dpi(j,i) * damp 1741 dpk(j,i) = dpk(j,i) * damp 1742 end do 1743 end do 1744 end if 1745 end if 1746c 1747c final polarization energy is half of value computed above 1748c 1749 ei = 0.5d0 * ei 1750 ek = 0.5d0 * ek 1751 return 1752 end 1753c 1754c 1755c ################################################################ 1756c ## ## 1757c ## subroutine empole2 -- multipole & polarization Hessian ## 1758c ## ## 1759c ################################################################ 1760c 1761c 1762c "empole2" calculates second derivatives of the multipole 1763c and dipole polarization energy for a single atom at a time 1764c 1765c 1766 subroutine empole2 (i) 1767 implicit none 1768 include 'sizes.i' 1769 include 'atoms.i' 1770 include 'deriv.i' 1771 include 'hessn.i' 1772 integer i,j,k 1773 real*8 eps,old 1774 real*8 d0(3,maxatm) 1775c 1776c 1777c set the stepsize to be used for numerical derivatives 1778c 1779 eps = 1.0d-7 1780c 1781c get multipole first derivatives for the base structure 1782c 1783 call empole2a (i) 1784 do k = 1, n 1785 do j = 1, 3 1786 d0(j,k) = dem(j,k) + dep(j,k) 1787 end do 1788 end do 1789c 1790c find numerical x-components via perturbed structures 1791c 1792 old = x(i) 1793 x(i) = x(i) + eps 1794 call empole2a (i) 1795 x(i) = old 1796 do k = 1, n 1797 do j = 1, 3 1798 hessx(j,k) = hessx(j,k) + (dem(j,k)+dep(j,k)-d0(j,k))/eps 1799 end do 1800 end do 1801c 1802c find numerical y-components via perturbed structures 1803c 1804 old = y(i) 1805 y(i) = y(i) + eps 1806 call empole2a (i) 1807 y(i) = old 1808 do k = 1, n 1809 do j = 1, 3 1810 hessy(j,k) = hessy(j,k) + (dem(j,k)+dep(j,k)-d0(j,k))/eps 1811 end do 1812 end do 1813c 1814c find numerical z-components via perturbed structures 1815c 1816 old = z(i) 1817 z(i) = z(i) + eps 1818 call empole2a (i) 1819 z(i) = old 1820 do k = 1, n 1821 do j = 1, 3 1822 hessz(j,k) = hessz(j,k) + (dem(j,k)+dep(j,k)-d0(j,k))/eps 1823 end do 1824 end do 1825 return 1826 end 1827c 1828c 1829c ################################################################# 1830c ## ## 1831c ## subroutine empole2a -- mpole & polar Hessian; numerical ## 1832c ## ## 1833c ################################################################# 1834c 1835c 1836c "empole2a" computes multipole and dipole polarization first 1837c derivatives for a single atom with respect to Cartesian 1838c coordinates; used to get finite difference second derivatives 1839c 1840c note that since polarization effects are many body, it is not 1841c really correct to neglect interactions where "iatom" is not 1842c directly involved as a multipole site or local coordinate axis; 1843c however, other sites are neglected in this version via the 1844c "ipart" and "kpart" checks to quickly get approximate values 1845c 1846c also, in the current version, the induced dipoles are not 1847c recomputed every time an atom is moved during computation of 1848c the numerical Hessian resulting in a further approximation 1849c 1850c 1851 subroutine empole2a (iatom) 1852 implicit none 1853 include 'sizes.i' 1854 include 'atoms.i' 1855 include 'bound.i' 1856 include 'cell.i' 1857 include 'chgpot.i' 1858 include 'couple.i' 1859 include 'deriv.i' 1860 include 'group.i' 1861 include 'mplpot.i' 1862 include 'mpole.i' 1863 include 'polar.i' 1864 include 'polgrp.i' 1865 include 'polpot.i' 1866 include 'shunt.i' 1867 include 'units.i' 1868 include 'usage.i' 1869 integer i,j,k,m 1870 integer jj,iatom 1871 integer ii,iz,ix 1872 integer kk,kz,kx 1873 real*8 eik,ei,ek,de 1874 real*8 f,fikm,fikp,fgrp 1875 real*8 shift,taper,dtaper 1876 real*8 trans,dtrans 1877 real*8 xi,yi,zi,xr,yr,zr 1878 real*8 r,r2,r3,r4,r5,r6,r7 1879 real*8 a(3,3),d(3,3,3,3) 1880 real*8 rpi(13),rpk(13) 1881 real*8 indi(3),indk(3) 1882 real*8 dm(3),dp(3),utu 1883 real*8 dmi(3,3),dmk(3,3) 1884 real*8 dpi(3,3),dpk(3,3) 1885 real*8 mscale(maxatm) 1886 real*8 pscale(maxatm) 1887 logical proceed,iuse,kuse 1888 logical ipart,kpart 1889c 1890c 1891c zero out the multipole and polarization first derivatives 1892c 1893 do i = 1, n 1894 do j = 1, 3 1895 dem(j,i) = 0.0d0 1896 dep(j,i) = 0.0d0 1897 end do 1898 end do 1899c 1900c set conversion factor and switching function coefficients 1901c 1902 f = electric / dielec 1903 call switch ('MPOLE') 1904c 1905c rotate multipole components and get rotation derivatives 1906c 1907 do i = 1, npole 1908 call rotmat (i,a) 1909 call rotsite (i,a) 1910 call drotmat (i,d) 1911 do j = 1, 3 1912 do k = 1, 3 1913 call drotpole (i,a,d,j,k) 1914 end do 1915 end do 1916 end do 1917c 1918c compute the induced dipoles at each polarizable atom; 1919c currently, we assume this was done in a calling routine 1920c 1921c call induce 1922c 1923c compute and partition multipole interaction energy 1924c 1925 do ii = 1, npole-1 1926 i = ipole(ii) 1927 xi = x(i) 1928 yi = y(i) 1929 zi = z(i) 1930 iz = zaxis(ii) 1931 ix = xaxis(ii) 1932 do j = 1, maxpole 1933 rpi(j) = rpole(j,ii) 1934 end do 1935 do j = 1, 3 1936 indi(j) = uind(j,ii) 1937 end do 1938 ipart = (i.eq.iatom .or. iz.eq.iatom .or. ix.eq.iatom) 1939 iuse = (use(i) .or. use(iz) .or. use(ix)) 1940 do j = ii+1, npole 1941 mscale(ipole(j)) = 1.0d0 1942 pscale(ipole(j)) = 1.0d0 1943 end do 1944 do j = 1, n12(i) 1945 mscale(i12(j,i)) = m2scale 1946 end do 1947 do j = 1, n13(i) 1948 mscale(i13(j,i)) = m3scale 1949 end do 1950 do j = 1, n14(i) 1951 mscale(i14(j,i)) = m4scale 1952 end do 1953 do j = 1, n15(i) 1954 mscale(i15(j,i)) = m5scale 1955 end do 1956 do j = 1, np11(i) 1957 pscale(ip11(j,i)) = p1scale 1958 end do 1959 do j = 1, np12(i) 1960 pscale(ip12(j,i)) = p2scale 1961 end do 1962 do j = 1, np13(i) 1963 pscale(ip13(j,i)) = p3scale 1964 end do 1965 do j = 1, np14(i) 1966 pscale(ip14(j,i)) = p4scale 1967 end do 1968c 1969c decide whether to compute the current interaction 1970c 1971 do kk = ii+1, npole 1972 k = ipole(kk) 1973 kz = zaxis(kk) 1974 kx = xaxis(kk) 1975 kpart = (k.eq.iatom .or. kz.eq.iatom .or. kx.eq.iatom) 1976 kuse = (use(k) .or. use(kz) .or. use(kx)) 1977 if (use_group) call groups (proceed,fgrp,2,i,k,0,0,0) 1978 proceed = .true. 1979 proceed = (ipart .or. kpart) 1980 if (proceed) proceed = (iuse .or. kuse) 1981c 1982c compute the energy contribution for this interaction 1983c 1984 if (proceed) then 1985 xr = x(k) - xi 1986 yr = y(k) - yi 1987 zr = z(k) - zi 1988 if (use_image) call image (xr,yr,zr,0) 1989 r2 = xr*xr + yr*yr + zr*zr 1990 if (r2 .le. off2) then 1991 do j = 1, maxpole 1992 rpk(j) = rpole(j,kk) 1993 end do 1994 do j = 1, 3 1995 indk(j) = uind(j,kk) 1996 end do 1997 r = sqrt(r2) 1998 call empik1 (ii,kk,xr,yr,zr,r,r2,rpi, 1999 & rpk,indi,indk,eik,ei,ek, 2000 & dm,dmi,dmk,dp,dpi,dpk,utu) 2001 fikm = f * mscale(k) 2002 fikp = f * pscale(k) 2003 eik = fikm * eik 2004 ei = fikp * ei 2005 ek = fikp * ek 2006 utu = fikp * utu 2007 do j = 1, 3 2008 dm(j) = fikm * dm(j) 2009 dp(j) = fikp * dp(j) 2010 do jj = 1, 3 2011 dmi(jj,j) = fikm * dmi(jj,j) 2012 dmk(jj,j) = fikm * dmk(jj,j) 2013 dpi(jj,j) = fikp * dpi(jj,j) 2014 dpk(jj,j) = fikp * dpk(jj,j) 2015 end do 2016 end do 2017c 2018c use shifted energy switching if near the cutoff distance 2019c 2020 if (r2 .gt. cut2) then 2021 fikm = fikm * rpi(1) * rpk(1) 2022 shift = fikm / (0.5d0*(off+cut)) 2023 eik = eik - shift 2024 r3 = r2 * r 2025 r4 = r2 * r2 2026 r5 = r2 * r3 2027 r6 = r3 * r3 2028 r7 = r3 * r4 2029 taper = c5*r5 + c4*r4 + c3*r3 2030 & + c2*r2 + c1*r + c0 2031 dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 2032 & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 2033 trans = fikm * (f7*r7 + f6*r6 + f5*r5 + f4*r4 2034 & + f3*r3 + f2*r2 + f1*r + f0) 2035 dtrans = fikm * (7.0d0*f7*r6 + 6.0d0*f6*r5 2036 & + 5.0d0*f5*r4 + 4.0d0*f4*r3 2037 & + 3.0d0*f3*r2 + 2.0d0*f2*r + f1) 2038 de = eik*dtaper + dtrans 2039 dm(1) = dm(1)*taper + de*(xr/r) 2040 dm(2) = dm(2)*taper + de*(yr/r) 2041 dm(3) = dm(3)*taper + de*(zr/r) 2042 de = (2.0d0*(ei+ek)+utu) * dtaper 2043 dp(1) = dp(1)*taper + de*(xr/r) 2044 dp(2) = dp(2)*taper + de*(yr/r) 2045 dp(3) = dp(3)*taper + de*(zr/r) 2046 do j = 1, 3 2047 do jj = 1, 3 2048 dmi(jj,j) = dmi(jj,j) * taper 2049 dmk(jj,j) = dmk(jj,j) * taper 2050 dpi(jj,j) = dpi(jj,j) * taper 2051 dpk(jj,j) = dpk(jj,j) * taper 2052 end do 2053 end do 2054 end if 2055c 2056c scale the interaction based on its group membership; 2057c polarization cannot be group scaled as it is not pairwise 2058c 2059 if (use_group) then 2060 do j = 1, 3 2061 dm(j) = dm(j) * fgrp 2062c dp(j) = dp(j) * fgrp 2063 do jj = 1, 3 2064 dmi(jj,j) = dmi(jj,j) * fgrp 2065 dmk(jj,j) = dmk(jj,j) * fgrp 2066c dpi(jj,j) = dpi(jj,j) * fgrp 2067c dpk(jj,j) = dpk(jj,j) * fgrp 2068 end do 2069 end do 2070 end if 2071c 2072c increment the multipole first derivative expressions 2073c 2074 dem(1,i) = dem(1,i) - dm(1) + dmi(1,1) 2075 dem(2,i) = dem(2,i) - dm(2) + dmi(1,2) 2076 dem(3,i) = dem(3,i) - dm(3) + dmi(1,3) 2077 dem(1,iz) = dem(1,iz) + dmi(2,1) 2078 dem(2,iz) = dem(2,iz) + dmi(2,2) 2079 dem(3,iz) = dem(3,iz) + dmi(2,3) 2080 dem(1,ix) = dem(1,ix) + dmi(3,1) 2081 dem(2,ix) = dem(2,ix) + dmi(3,2) 2082 dem(3,ix) = dem(3,ix) + dmi(3,3) 2083 dem(1,k) = dem(1,k) + dm(1) + dmk(1,1) 2084 dem(2,k) = dem(2,k) + dm(2) + dmk(1,2) 2085 dem(3,k) = dem(3,k) + dm(3) + dmk(1,3) 2086 dem(1,kz) = dem(1,kz) + dmk(2,1) 2087 dem(2,kz) = dem(2,kz) + dmk(2,2) 2088 dem(3,kz) = dem(3,kz) + dmk(2,3) 2089 dem(1,kx) = dem(1,kx) + dmk(3,1) 2090 dem(2,kx) = dem(2,kx) + dmk(3,2) 2091 dem(3,kx) = dem(3,kx) + dmk(3,3) 2092c 2093c increment the polarization first derivative expressions 2094c 2095 dep(1,i) = dep(1,i) - dp(1) + dpi(1,1) 2096 dep(2,i) = dep(2,i) - dp(2) + dpi(1,2) 2097 dep(3,i) = dep(3,i) - dp(3) + dpi(1,3) 2098 dep(1,iz) = dep(1,iz) + dpi(2,1) 2099 dep(2,iz) = dep(2,iz) + dpi(2,2) 2100 dep(3,iz) = dep(3,iz) + dpi(2,3) 2101 dep(1,ix) = dep(1,ix) + dpi(3,1) 2102 dep(2,ix) = dep(2,ix) + dpi(3,2) 2103 dep(3,ix) = dep(3,ix) + dpi(3,3) 2104 dep(1,k) = dep(1,k) + dp(1) + dpk(1,1) 2105 dep(2,k) = dep(2,k) + dp(2) + dpk(1,2) 2106 dep(3,k) = dep(3,k) + dp(3) + dpk(1,3) 2107 dep(1,kz) = dep(1,kz) + dpk(2,1) 2108 dep(2,kz) = dep(2,kz) + dpk(2,2) 2109 dep(3,kz) = dep(3,kz) + dpk(2,3) 2110 dep(1,kx) = dep(1,kx) + dpk(3,1) 2111 dep(2,kx) = dep(2,kx) + dpk(3,2) 2112 dep(3,kx) = dep(3,kx) + dpk(3,3) 2113 end if 2114 end if 2115 end do 2116 end do 2117c 2118c for periodic boundary conditions with large cutoffs 2119c neighbors must be found by the replicates method 2120c 2121 if (.not. use_replica) return 2122c 2123c calculate interaction energy with other unit cells 2124c 2125 do ii = 1, npole 2126 i = ipole(ii) 2127 xi = x(i) 2128 yi = y(i) 2129 zi = z(i) 2130 iz = zaxis(ii) 2131 ix = xaxis(ii) 2132 do j = 1, maxpole 2133 rpi(j) = rpole(j,ii) 2134 end do 2135 do j = 1, 3 2136 indi(j) = uind(j,ii) 2137 end do 2138 ipart = (i.eq.iatom .or. iz.eq.iatom .or. ix.eq.iatom) 2139 iuse = (use(i) .or. use(iz) .or. use(ix)) 2140 do j = ii, npole 2141 mscale(ipole(j)) = 1.0d0 2142 pscale(ipole(j)) = 1.0d0 2143 end do 2144 do j = 1, n12(i) 2145 mscale(i12(j,i)) = m2scale 2146 end do 2147 do j = 1, n13(i) 2148 mscale(i13(j,i)) = m3scale 2149 end do 2150 do j = 1, n14(i) 2151 mscale(i14(j,i)) = m4scale 2152 end do 2153 do j = 1, n15(i) 2154 mscale(i15(j,i)) = m5scale 2155 end do 2156 do j = 1, np11(i) 2157 pscale(ip11(j,i)) = p1scale 2158 end do 2159 do j = 1, np12(i) 2160 pscale(ip12(j,i)) = p2scale 2161 end do 2162 do j = 1, np13(i) 2163 pscale(ip13(j,i)) = p3scale 2164 end do 2165 do j = 1, np14(i) 2166 pscale(ip14(j,i)) = p4scale 2167 end do 2168c 2169c decide whether to compute the current interaction 2170c 2171 do kk = ii, npole 2172 k = ipole(kk) 2173 kz = zaxis(kk) 2174 kx = xaxis(kk) 2175 kpart = (k.eq.iatom .or. kz.eq.iatom .or. kx.eq.iatom) 2176 kuse = (use(k) .or. use(kz) .or. use(kx)) 2177 if (use_group) call groups (proceed,fgrp,2,i,k,0,0,0) 2178 proceed = .true. 2179 proceed = (ipart .or. kpart) 2180 if (proceed) proceed = (iuse .or. kuse) 2181c 2182c compute the energy contribution for this interaction 2183c 2184 if (proceed) then 2185 do m = 2, ncell 2186 xr = x(k) - xi 2187 yr = y(k) - yi 2188 zr = z(k) - zi 2189 call image (xr,yr,zr,m) 2190 r2 = xr*xr + yr*yr + zr*zr 2191 if (r2 .le. off2) then 2192 do j = 1, maxpole 2193 rpk(j) = rpole(j,kk) 2194 end do 2195 do j = 1, 3 2196 indk(j) = uind(j,kk) 2197 end do 2198 r = sqrt(r2) 2199 call empik1 (ii,kk,xr,yr,zr,r,r2,rpi, 2200 & rpk,indi,indk,eik,ei,ek, 2201 & dm,dmi,dmk,dp,dpi,dpk,utu) 2202 fikm = f 2203 fikp = f 2204 if (use_polymer) then 2205 if (r2 .le. polycut2) then 2206 fikm = fikm * mscale(kk) 2207 fikp = fikp * pscale(kk) 2208 end if 2209 end if 2210 eik = fikm * eik 2211 ei = fikp * ei 2212 ek = fikp * ek 2213 utu = fikp * utu 2214 do j = 1, 3 2215 dm(j) = fikm * dm(j) 2216 dp(j) = fikp * dp(j) 2217 do jj = 1, 3 2218 dmi(jj,j) = fikm * dmi(jj,j) 2219 dmk(jj,j) = fikm * dmk(jj,j) 2220 dpi(jj,j) = fikp * dpi(jj,j) 2221 dpk(jj,j) = fikp * dpk(jj,j) 2222 end do 2223 end do 2224c 2225c use shifted energy switching if near the cutoff distance 2226c 2227 if (r2 .gt. cut2) then 2228 fikm = fikm * rpi(1) * rpk(1) 2229 shift = fikm / (0.5d0*(off+cut)) 2230 eik = eik - shift 2231 r3 = r2 * r 2232 r4 = r2 * r2 2233 r5 = r2 * r3 2234 r6 = r3 * r3 2235 r7 = r3 * r4 2236 taper = c5*r5 + c4*r4 + c3*r3 2237 & + c2*r2 + c1*r + c0 2238 dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 2239 & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 2240 trans = fikm * (f7*r7 + f6*r6 + f5*r5 + f4*r4 2241 & + f3*r3 + f2*r2 + f1*r + f0) 2242 dtrans = fikm * (7.0d0*f7*r6 + 6.0d0*f6*r5 2243 & + 5.0d0*f5*r4 + 4.0d0*f4*r3 2244 & + 3.0d0*f3*r2 + 2.0d0*f2*r + f1) 2245 de = eik*dtaper + dtrans 2246 dm(1) = dm(1)*taper + de*(xr/r) 2247 dm(2) = dm(2)*taper + de*(yr/r) 2248 dm(3) = dm(3)*taper + de*(zr/r) 2249 de = (2.0d0*(ei+ek)+utu) * dtaper 2250 dp(1) = dp(1)*taper + de*(xr/r) 2251 dp(2) = dp(2)*taper + de*(yr/r) 2252 dp(3) = dp(3)*taper + de*(zr/r) 2253 do j = 1, 3 2254 do jj = 1, 3 2255 dmi(jj,j) = dmi(jj,j) * taper 2256 dmk(jj,j) = dmk(jj,j) * taper 2257 dpi(jj,j) = dpi(jj,j) * taper 2258 dpk(jj,j) = dpk(jj,j) * taper 2259 end do 2260 end do 2261 end if 2262c 2263c scale the interaction based on its group membership; 2264c polarization cannot be group scaled as it is not pairwise 2265c 2266 if (use_group) then 2267 do j = 1, 3 2268 dm(j) = dm(j) * fgrp 2269c dp(j) = dp(j) * fgrp 2270 do jj = 1, 3 2271 dmi(jj,j) = dmi(jj,j) * fgrp 2272 dmk(jj,j) = dmk(jj,j) * fgrp 2273c dpi(jj,j) = dpi(jj,j) * fgrp 2274c dpk(jj,j) = dpk(jj,j) * fgrp 2275 end do 2276 end do 2277 end if 2278c 2279c increment the multipole first derivative expressions 2280c 2281 dem(1,i) = dem(1,i) - dm(1) + dmi(1,1) 2282 dem(2,i) = dem(2,i) - dm(2) + dmi(1,2) 2283 dem(3,i) = dem(3,i) - dm(3) + dmi(1,3) 2284 dem(1,iz) = dem(1,iz) + dmi(2,1) 2285 dem(2,iz) = dem(2,iz) + dmi(2,2) 2286 dem(3,iz) = dem(3,iz) + dmi(2,3) 2287 dem(1,ix) = dem(1,ix) + dmi(3,1) 2288 dem(2,ix) = dem(2,ix) + dmi(3,2) 2289 dem(3,ix) = dem(3,ix) + dmi(3,3) 2290 if (i .ne. k) then 2291 dem(1,k) = dem(1,k) + dm(1) + dmk(1,1) 2292 dem(2,k) = dem(2,k) + dm(2) + dmk(1,2) 2293 dem(3,k) = dem(3,k) + dm(3) + dmk(1,3) 2294 dem(1,kz) = dem(1,kz) + dmk(2,1) 2295 dem(2,kz) = dem(2,kz) + dmk(2,2) 2296 dem(3,kz) = dem(3,kz) + dmk(2,3) 2297 dem(1,kx) = dem(1,kx) + dmk(3,1) 2298 dem(2,kx) = dem(2,kx) + dmk(3,2) 2299 dem(3,kx) = dem(3,kx) + dmk(3,3) 2300 end if 2301c 2302c increment the polarization first derivative expressions 2303c 2304 dep(1,i) = dep(1,i) - dp(1) + dpi(1,1) 2305 dep(2,i) = dep(2,i) - dp(2) + dpi(1,2) 2306 dep(3,i) = dep(3,i) - dp(3) + dpi(1,3) 2307 dep(1,iz) = dep(1,iz) + dpi(2,1) 2308 dep(2,iz) = dep(2,iz) + dpi(2,2) 2309 dep(3,iz) = dep(3,iz) + dpi(2,3) 2310 dep(1,ix) = dep(1,ix) + dpi(3,1) 2311 dep(2,ix) = dep(2,ix) + dpi(3,2) 2312 dep(3,ix) = dep(3,ix) + dpi(3,3) 2313 if (i .ne. k) then 2314 dep(1,k) = dep(1,k) + dp(1) + dpk(1,1) 2315 dep(2,k) = dep(2,k) + dp(2) + dpk(1,2) 2316 dep(3,k) = dep(3,k) + dp(3) + dpk(1,3) 2317 dep(1,kz) = dep(1,kz) + dpk(2,1) 2318 dep(2,kz) = dep(2,kz) + dpk(2,2) 2319 dep(3,kz) = dep(3,kz) + dpk(2,3) 2320 dep(1,kx) = dep(1,kx) + dpk(3,1) 2321 dep(2,kx) = dep(2,kx) + dpk(3,2) 2322 dep(3,kx) = dep(3,kx) + dpk(3,3) 2323 end if 2324 end if 2325 end do 2326 end if 2327 end do 2328 end do 2329 return 2330 end 2331c 2332c 2333c ############################################################### 2334c ## ## 2335c ## subroutine empole3a -- double loop multipole analysis ## 2336c ## ## 2337c ############################################################### 2338c 2339c 2340c "empole3a" calculates the electrostatic energy due to 2341c atomic multipole interactions and dipole polarizability, and 2342c partitions the energy among the atoms using a double loop 2343c 2344c 2345 subroutine empole3a 2346 implicit none 2347 include 'sizes.i' 2348 include 'action.i' 2349 include 'analyz.i' 2350 include 'atmtyp.i' 2351 include 'atoms.i' 2352 include 'bound.i' 2353 include 'cell.i' 2354 include 'chgpot.i' 2355 include 'couple.i' 2356 include 'energi.i' 2357 include 'group.i' 2358 include 'inform.i' 2359 include 'iounit.i' 2360 include 'moment.i' 2361 include 'mpole.i' 2362 include 'polar.i' 2363 include 'shunt.i' 2364 include 'units.i' 2365 include 'usage.i' 2366 integer i,j,k,m 2367 integer ii,iz,ix 2368 integer kk,kz,kx 2369 integer skip(maxatm) 2370 real*8 eik,ei,ek,a(3,3) 2371 real*8 shift,taper,trans 2372 real*8 f,fik,fgrp 2373 real*8 xr,yr,zr,r 2374 real*8 r2,r3,r4,r5,r6,r7 2375 real*8 rpi(13),rpk(13) 2376 real*8 indi(3),indk(3) 2377 real*8 weight,xcenter,ycenter,zcenter 2378 real*8 totchg,xsum,ysum,zsum 2379 logical header,huge,proceed,iuse,kuse 2380c 2381c 2382c zero out multipole and polarization energy and partitioning 2383c 2384 nem = 0 2385 nep = 0 2386 em = 0.0d0 2387 ep = 0.0d0 2388 do i = 1, n 2389 aem(i) = 0.0d0 2390 aep(i) = 0.0d0 2391 end do 2392 header = .true. 2393c 2394c zero out the list of atoms to be skipped 2395c 2396 do i = 1, n 2397 skip(i) = 0 2398 end do 2399c 2400c set conversion factor and switching function coefficients 2401c 2402 f = electric / dielec 2403 call switch ('MPOLE') 2404c 2405c rotate the multipole components into the global frame 2406c 2407 call rotpole 2408c 2409c compute the induced dipoles at each atom 2410c 2411 call induce 2412c 2413c find the center of mass of the set of active atoms 2414c 2415 weight = 0.0d0 2416 xcenter = 0.0d0 2417 ycenter = 0.0d0 2418 zcenter = 0.0d0 2419 do i = 1, n 2420 if (use(i)) then 2421 weight = weight + mass(i) 2422 xcenter = xcenter + x(i)*mass(i) 2423 ycenter = ycenter + y(i)*mass(i) 2424 zcenter = zcenter + z(i)*mass(i) 2425 end if 2426 end do 2427 if (weight .ne. 0.0d0) then 2428 xcenter = xcenter / weight 2429 ycenter = ycenter / weight 2430 zcenter = zcenter / weight 2431 end if 2432c 2433c get net charge and dipole components over active atoms 2434c 2435 totchg = 0.0d0 2436 xsum = 0.0d0 2437 ysum = 0.0d0 2438 zsum = 0.0d0 2439 do i = 1, npole 2440 k = ipole(i) 2441 if (use(k)) then 2442 totchg = totchg + rpole(1,i) 2443 xsum = xsum + (x(k)-xcenter)*rpole(1,i) 2444 ysum = ysum + (y(k)-ycenter)*rpole(1,i) 2445 zsum = zsum + (z(k)-zcenter)*rpole(1,i) 2446 xsum = xsum + rpole(2,i) + uind(1,i) 2447 ysum = ysum + rpole(3,i) + uind(2,i) 2448 zsum = zsum + rpole(4,i) + uind(3,i) 2449 end if 2450 end do 2451 netchg = netchg + totchg 2452 xdipole = xdipole + debye*xsum 2453 ydipole = ydipole + debye*ysum 2454 zdipole = zdipole + debye*zsum 2455c 2456c calculate the multipole interaction energy term 2457c 2458 do ii = 1, npole-1 2459 i = ipole(ii) 2460 iz = zaxis(ii) 2461 ix = xaxis(ii) 2462 iuse = (use(i) .or. use(iz) .or. use(ix)) 2463 do j = 1, n12(i) 2464 skip(i12(j,i)) = i * chg12use 2465 end do 2466 do j = 1, n13(i) 2467 skip(i13(j,i)) = i * chg13use 2468 end do 2469 do j = 1, n14(i) 2470 skip(i14(j,i)) = i * chg14use 2471 end do 2472 do j = 1, maxpole 2473 rpi(j) = rpole(j,ii) 2474 end do 2475 do j = 1, 3 2476 indi(j) = uind(j,ii) 2477 end do 2478 do kk = ii+1, npole 2479 k = ipole(kk) 2480 kz = zaxis(kk) 2481 kx = xaxis(kk) 2482 kuse = (use(k) .or. use(kz) .or. use(kx)) 2483c 2484c decide whether to compute the current interaction 2485c 2486 proceed = .true. 2487 if (use_group) call groups (proceed,fgrp,2,i,k,0,0,0) 2488 if (proceed) proceed = (iuse .or. kuse) 2489 if (proceed) proceed = (skip(k) .ne. i) 2490c 2491c compute the energy contribution for this interaction 2492c 2493 if (proceed) then 2494 xr = x(k) - x(i) 2495 yr = y(k) - y(i) 2496 zr = z(k) - z(i) 2497 if (use_image) call image (xr,yr,zr,0) 2498 r2 = xr*xr + yr*yr + zr*zr 2499 if (r2 .le. off2) then 2500 do j = 1, maxpole 2501 rpk(j) = rpole(j,kk) 2502 end do 2503 do j = 1, 3 2504 indk(j) = uind(j,kk) 2505 end do 2506 r = sqrt(r2) 2507 call empik (ii,kk,xr,yr,zr,r,rpi,rpk, 2508 & indi,indk,eik,ei,ek) 2509 fik = f 2510 if (skip(k) .eq. -i) fik = fik / chgscale 2511 eik = fik * eik 2512 ei = fik * ei 2513 ek = fik * ek 2514c 2515c use shifted energy switching if near the cutoff distance 2516c 2517 fik = fik * rpi(1) * rpk(1) 2518 shift = fik / (0.5d0*(off+cut)) 2519 eik = eik - shift 2520 if (r2 .gt. cut2) then 2521 r3 = r2 * r 2522 r4 = r2 * r2 2523 r5 = r2 * r3 2524 r6 = r3 * r3 2525 r7 = r3 * r4 2526 taper = c5*r5 + c4*r4 + c3*r3 2527 & + c2*r2 + c1*r + c0 2528 trans = fik * (f7*r7 + f6*r6 + f5*r5 + f4*r4 2529 & + f3*r3 + f2*r2 + f1*r + f0) 2530 eik = eik * taper + trans 2531 ei = ei * taper 2532 ek = ek * taper 2533 end if 2534c 2535c scale the interaction based on its group membership 2536c 2537 if (use_group) then 2538 eik = eik * fgrp 2539 ei = ei * fgrp 2540 ek = ek * fgrp 2541 end if 2542c 2543c increment the overall multipole and polarization energies 2544c 2545 nem = nem + 1 2546 em = em + eik 2547 aem(i) = aem(i) + 0.5d0*eik 2548 aem(k) = aem(k) + 0.5d0*eik 2549 nep = nep + 2 2550 ep = ep + ei + ek 2551 aep(i) = aep(i) + 0.5d0*(ei+ek) 2552 aep(k) = aep(k) + 0.5d0*(ei+ek) 2553c 2554c print a warning if the energy of this interaction is large 2555c 2556 huge = (max(abs(eik),abs(ei),abs(ek)) .gt. 100.0d0) 2557 if (debug .or. (verbose.and.huge)) then 2558 if (header) then 2559 header = .false. 2560 write (iout,10) 2561 10 format (/,' Individual Multipole and', 2562 & ' Polarization Interactions :', 2563 & //,' Type',13x,'Atom Names', 2564 & 9x,'Distance',6x,'Energies', 2565 & ' (MPole, Pol1, Pol2)',/) 2566 end if 2567 write (iout,20) i,name(i),k,name(k),r,eik,ei,ek 2568 20 format (' M-Pole',5x,i5,'-',a3,1x,i5,'-',a3, 2569 & 5x,f8.4,3f12.4) 2570 end if 2571 end if 2572 end if 2573 end do 2574 end do 2575c 2576c for periodic boundary conditions with large cutoffs 2577c neighbors must be found by the replicates method 2578c 2579 if (.not. use_replica) return 2580c 2581c calculate interaction energy with other unit cells 2582c 2583 do ii = 1, npole 2584 i = ipole(ii) 2585 iz = zaxis(ii) 2586 ix = xaxis(ii) 2587 iuse = (use(i) .or. use(iz) .or. use(ix)) 2588 do j = 1, maxpole 2589 rpi(j) = rpole(j,ii) 2590 end do 2591 do j = 1, 3 2592 indi(j) = uind(j,ii) 2593 end do 2594 do kk = ii, npole 2595 k = ipole(kk) 2596 kz = zaxis(kk) 2597 kx = xaxis(kk) 2598 kuse = (use(k) .or. use(kz) .or. use(kx)) 2599c 2600c decide whether to compute the current interaction 2601c 2602 proceed = .true. 2603 if (use_group) call groups (proceed,fgrp,2,i,k,0,0,0) 2604 if (proceed) proceed = (iuse .or. kuse) 2605c 2606c compute the energy contribution for this interaction 2607c 2608 if (proceed) then 2609 do m = 2, ncell 2610 xr = x(k) - x(i) 2611 yr = y(k) - y(i) 2612 zr = z(k) - z(i) 2613 call image (xr,yr,zr,m) 2614 r2 = xr*xr + yr*yr + zr*zr 2615 if (r2 .le. off2) then 2616 do j = 1, maxpole 2617 rpk(j) = rpole(j,kk) 2618 end do 2619 do j = 1, 3 2620 indk(j) = uind(j,kk) 2621 end do 2622 r = sqrt(r2) 2623 call empik (ii,kk,xr,yr,zr,r,rpi,rpk, 2624 & indi,indk,eik,ei,ek) 2625 fik = f 2626 eik = fik * eik 2627 ei = fik * ei 2628 ek = fik * ek 2629c 2630c use shifted energy switching if near the cutoff distance 2631c 2632 fik = fik * rpi(1) * rpk(1) 2633 shift = fik / (0.5d0*(off+cut)) 2634 eik = eik - shift 2635 if (r2 .gt. cut2) then 2636 r3 = r2 * r 2637 r4 = r2 * r2 2638 r5 = r2 * r3 2639 r6 = r3 * r3 2640 r7 = r3 * r4 2641 taper = c5*r5 + c4*r4 + c3*r3 2642 & + c2*r2 + c1*r + c0 2643 trans = fik * (f7*r7 + f6*r6 + f5*r5 + f4*r4 2644 & + f3*r3 + f2*r2 + f1*r + f0) 2645 eik = eik * taper + trans 2646 ei = ei * taper 2647 ek = ek * taper 2648 end if 2649c 2650c scale the interaction based on its group membership 2651c 2652 if (use_group) then 2653 eik = eik * fgrp 2654 ei = ei * fgrp 2655 ek = ek * fgrp 2656 end if 2657c 2658c increment the overall multipole and polarization energies 2659c 2660 nem = nem + 1 2661 nep = nep + 2 2662 if (i .eq. k) then 2663 em = em + 0.5d0*eik 2664 aem(i) = aem(i) + 0.5d0*eik 2665 ep = ep + ei 2666 aep(i) = aep(i) + ei 2667 else 2668 em = em + eik 2669 aem(i) = aem(i) + 0.5d0*eik 2670 aem(k) = aem(k) + 0.5d0*eik 2671 ep = ep + ei + ek 2672 aep(i) = aep(i) + 0.5d0*(ei+ek) 2673 aep(k) = aep(k) + 0.5d0*(ei+ek) 2674 end if 2675c 2676c print a warning if the energy of this interaction is large 2677c 2678 huge = (max(abs(eik),abs(ei),abs(ek)) .gt. 100.0d0) 2679 if (debug .or. (verbose.and.huge)) then 2680 if (header) then 2681 header = .false. 2682 write (iout,30) 2683 30 format (/,' Individual Multipole and', 2684 & ' Polarization Interactions :', 2685 & //,' Type',13x,'Atom Names', 2686 & 9x,'Distance',6x,'Energies', 2687 & ' (MPole, Pol1, Pol2)',/) 2688 end if 2689 write (iout,40) i,name(i),k,name(k),r,eik,ei,ek 2690 40 format (' M-Pole',5x,i5,'-',a3,1x,i5,'-',a3, 2691 & 1x,'(X)',1x,f8.4,3f12.4) 2692 end if 2693 end if 2694 end do 2695 end if 2696 end do 2697 end do 2698 return 2699 end 2700