1c 2c 3c ############################################################# 4c ## COPYRIGHT (C) 1999 by Pengyu Ren & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ############################################################# 7c 8c ############################################################# 9c ## ## 10c ## subroutine empole -- atomic multipole moment energy ## 11c ## ## 12c ############################################################# 13c 14c 15c "empole" calculates the electrostatic energy due to atomic 16c multipole interactions 17c 18c 19 subroutine empole 20 use limits 21 implicit none 22c 23c 24c choose the method for summing over multipole interactions 25c 26 if (use_ewald) then 27 if (use_mlist) then 28 call empole0d 29 else 30 call empole0c 31 end if 32 else 33 if (use_mlist) then 34 call empole0b 35 else 36 call empole0a 37 end if 38 end if 39 return 40 end 41c 42c 43c ############################################################# 44c ## ## 45c ## subroutine empole0a -- double loop multipole energy ## 46c ## ## 47c ############################################################# 48c 49c 50c "empole0a" calculates the atomic multipole interaction energy 51c using a double loop 52c 53c 54 subroutine empole0a 55 use atoms 56 use bound 57 use cell 58 use chgpen 59 use chgpot 60 use couple 61 use energi 62 use group 63 use math 64 use mplpot 65 use mpole 66 use polpot 67 use shunt 68 use usage 69 implicit none 70 integer i,j,k 71 integer ii,kk 72 integer ix,iy,iz 73 integer kx,ky,kz 74 real*8 e,f,fgrp 75 real*8 xi,yi,zi 76 real*8 xr,yr,zr 77 real*8 r,r2,rr1,rr3 78 real*8 rr5,rr7,rr9 79 real*8 rr1i,rr3i,rr5i 80 real*8 rr1k,rr3k,rr5k 81 real*8 rr1ik,rr3ik,rr5ik 82 real*8 rr7ik,rr9ik 83 real*8 ci,dix,diy,diz 84 real*8 qixx,qixy,qixz 85 real*8 qiyy,qiyz,qizz 86 real*8 ck,dkx,dky,dkz 87 real*8 qkxx,qkxy,qkxz 88 real*8 qkyy,qkyz,qkzz 89 real*8 dir,dkr,dik,qik 90 real*8 qix,qiy,qiz,qir 91 real*8 qkx,qky,qkz,qkr 92 real*8 diqk,dkqi,qiqk 93 real*8 corei,corek 94 real*8 vali,valk 95 real*8 alphai,alphak 96 real*8 term1,term2,term3 97 real*8 term4,term5 98 real*8 term1i,term2i,term3i 99 real*8 term1k,term2k,term3k 100 real*8 term1ik,term2ik,term3ik 101 real*8 term4ik,term5ik 102 real*8 dmpi(9),dmpk(9) 103 real*8 dmpik(9) 104 real*8, allocatable :: mscale(:) 105 logical proceed,usei,usek 106 character*6 mode 107c 108c 109c zero out the total atomic multipole energy 110c 111 em = 0.0d0 112 if (npole .eq. 0) return 113c 114c check the sign of multipole components at chiral sites 115c 116 call chkpole 117c 118c rotate the multipole components into the global frame 119c 120 call rotpole 121c 122c perform dynamic allocation of some local arrays 123c 124 allocate (mscale(n)) 125c 126c initialize connected atom exclusion coefficients 127c 128 do i = 1, n 129 mscale(i) = 1.0d0 130 end do 131c 132c set conversion factor, cutoff and switching coefficients 133c 134 f = electric / dielec 135 mode = 'MPOLE' 136 call switch (mode) 137c 138c calculate the multipole interaction energy term 139c 140 do ii = 1, npole-1 141 i = ipole(ii) 142 iz = zaxis(ii) 143 ix = xaxis(ii) 144 iy = abs(yaxis(ii)) 145 xi = x(i) 146 yi = y(i) 147 zi = z(i) 148 ci = rpole(1,ii) 149 dix = rpole(2,ii) 150 diy = rpole(3,ii) 151 diz = rpole(4,ii) 152 qixx = rpole(5,ii) 153 qixy = rpole(6,ii) 154 qixz = rpole(7,ii) 155 qiyy = rpole(9,ii) 156 qiyz = rpole(10,ii) 157 qizz = rpole(13,ii) 158 if (use_chgpen) then 159 corei = pcore(ii) 160 vali = pval(ii) 161 alphai = palpha(ii) 162 end if 163 usei = (use(i) .or. use(iz) .or. use(ix) .or. use(iy)) 164c 165c set exclusion coefficients for connected atoms 166c 167 do j = 1, n12(i) 168 mscale(i12(j,i)) = m2scale 169 end do 170 do j = 1, n13(i) 171 mscale(i13(j,i)) = m3scale 172 end do 173 do j = 1, n14(i) 174 mscale(i14(j,i)) = m4scale 175 end do 176 do j = 1, n15(i) 177 mscale(i15(j,i)) = m5scale 178 end do 179c 180c evaluate all sites within the cutoff distance 181c 182 do kk = ii+1, npole 183 k = ipole(kk) 184 kz = zaxis(kk) 185 kx = xaxis(kk) 186 ky = abs(yaxis(kk)) 187 usek = (use(k) .or. use(kz) .or. use(kx) .or. use(ky)) 188 proceed = .true. 189 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 190 if (.not. use_intra) proceed = .true. 191 if (proceed) proceed = (usei .or. usek) 192 if (proceed) then 193 xr = x(k) - xi 194 yr = y(k) - yi 195 zr = z(k) - zi 196 if (use_bounds) call image (xr,yr,zr) 197 r2 = xr*xr + yr* yr + zr*zr 198 if (r2 .le. off2) then 199 r = sqrt(r2) 200 ck = rpole(1,kk) 201 dkx = rpole(2,kk) 202 dky = rpole(3,kk) 203 dkz = rpole(4,kk) 204 qkxx = rpole(5,kk) 205 qkxy = rpole(6,kk) 206 qkxz = rpole(7,kk) 207 qkyy = rpole(9,kk) 208 qkyz = rpole(10,kk) 209 qkzz = rpole(13,kk) 210c 211c intermediates involving moments and separation distance 212c 213 dir = dix*xr + diy*yr + diz*zr 214 qix = qixx*xr + qixy*yr + qixz*zr 215 qiy = qixy*xr + qiyy*yr + qiyz*zr 216 qiz = qixz*xr + qiyz*yr + qizz*zr 217 qir = qix*xr + qiy*yr + qiz*zr 218 dkr = dkx*xr + dky*yr + dkz*zr 219 qkx = qkxx*xr + qkxy*yr + qkxz*zr 220 qky = qkxy*xr + qkyy*yr + qkyz*zr 221 qkz = qkxz*xr + qkyz*yr + qkzz*zr 222 qkr = qkx*xr + qky*yr + qkz*zr 223 dik = dix*dkx + diy*dky + diz*dkz 224 qik = qix*qkx + qiy*qky + qiz*qkz 225 diqk = dix*qkx + diy*qky + diz*qkz 226 dkqi = dkx*qix + dky*qiy + dkz*qiz 227 qiqk = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 228 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 229c 230c get reciprocal distance terms for this interaction 231c 232 rr1 = f * mscale(k) / r 233 rr3 = rr1 / r2 234 rr5 = 3.0d0 * rr3 / r2 235 rr7 = 5.0d0 * rr5 / r2 236 rr9 = 7.0d0 * rr7 / r2 237c 238c find damped multipole intermediates and energy value 239c 240 if (use_chgpen) then 241 corek = pcore(kk) 242 valk = pval(kk) 243 alphak = palpha(kk) 244 term1 = corei*corek 245 term1i = corek*vali 246 term2i = corek*dir 247 term3i = corek*qir 248 term1k = corei*valk 249 term2k = -corei*dkr 250 term3k = corei*qkr 251 term1ik = vali*valk 252 term2ik = valk*dir - vali*dkr + dik 253 term3ik = vali*qkr + valk*qir - dir*dkr 254 & + 2.0d0*(dkqi-diqk+qiqk) 255 term4ik = dir*qkr - dkr*qir - 4.0d0*qik 256 term5ik = qir*qkr 257 call damppole (r,9,alphai,alphak, 258 & dmpi,dmpk,dmpik) 259 rr1i = dmpi(1)*rr1 260 rr3i = dmpi(3)*rr3 261 rr5i = dmpi(5)*rr5 262 rr1k = dmpk(1)*rr1 263 rr3k = dmpk(3)*rr3 264 rr5k = dmpk(5)*rr5 265 rr1ik = dmpik(1)*rr1 266 rr3ik = dmpik(3)*rr3 267 rr5ik = dmpik(5)*rr5 268 rr7ik = dmpik(7)*rr7 269 rr9ik = dmpik(9)*rr9 270 e = term1*rr1 + term1i*rr1i 271 & + term1k*rr1k + term1ik*rr1ik 272 & + term2i*rr3i + term2k*rr3k 273 & + term2ik*rr3ik + term3i*rr5i 274 & + term3k*rr5k + term3ik*rr5ik 275 & + term4ik*rr7ik + term5ik*rr9ik 276c 277c find standard multipole intermediates and energy value 278c 279 else 280 term1 = ci*ck 281 term2 = ck*dir - ci*dkr + dik 282 term3 = ci*qkr + ck*qir - dir*dkr 283 & + 2.0d0*(dkqi-diqk+qiqk) 284 term4 = dir*qkr - dkr*qir - 4.0d0*qik 285 term5 = qir*qkr 286 e = term1*rr1 + term2*rr3 + term3*rr5 287 & + term4*rr7 + term5*rr9 288 end if 289c 290c compute the energy contribution for this interaction 291c 292 if (use_group) e = e * fgrp 293 em = em + e 294 end if 295 end if 296 end do 297c 298c reset exclusion coefficients for connected atoms 299c 300 do j = 1, n12(i) 301 mscale(i12(j,i)) = 1.0d0 302 end do 303 do j = 1, n13(i) 304 mscale(i13(j,i)) = 1.0d0 305 end do 306 do j = 1, n14(i) 307 mscale(i14(j,i)) = 1.0d0 308 end do 309 do j = 1, n15(i) 310 mscale(i15(j,i)) = 1.0d0 311 end do 312 end do 313c 314c for periodic boundary conditions with large cutoffs 315c neighbors must be found by the replicates method 316c 317 if (use_replica) then 318c 319c calculate interaction energy with other unit cells 320c 321 do ii = 1, npole 322 i = ipole(ii) 323 iz = zaxis(ii) 324 ix = xaxis(ii) 325 iy = abs(yaxis(ii)) 326 xi = x(i) 327 yi = y(i) 328 zi = z(i) 329 ci = rpole(1,ii) 330 dix = rpole(2,ii) 331 diy = rpole(3,ii) 332 diz = rpole(4,ii) 333 qixx = rpole(5,ii) 334 qixy = rpole(6,ii) 335 qixz = rpole(7,ii) 336 qiyy = rpole(9,ii) 337 qiyz = rpole(10,ii) 338 qizz = rpole(13,ii) 339 if (use_chgpen) then 340 corei = pcore(ii) 341 vali = pval(ii) 342 alphai = palpha(ii) 343 end if 344 usei = (use(i) .or. use(iz) .or. use(ix) .or. use(iy)) 345c 346c set exclusion coefficients for connected atoms 347c 348 do j = 1, n12(i) 349 mscale(i12(j,i)) = m2scale 350 end do 351 do j = 1, n13(i) 352 mscale(i13(j,i)) = m3scale 353 end do 354 do j = 1, n14(i) 355 mscale(i14(j,i)) = m4scale 356 end do 357 do j = 1, n15(i) 358 mscale(i15(j,i)) = m5scale 359 end do 360c 361c evaluate all sites within the cutoff distance 362c 363 do kk = ii, npole 364 k = ipole(kk) 365 kz = zaxis(kk) 366 kx = xaxis(kk) 367 ky = abs(yaxis(kk)) 368 usek = (use(k) .or. use(kz) .or. use(kx) .or. use(ky)) 369 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 370 proceed = .true. 371 if (proceed) proceed = (usei .or. usek) 372 if (proceed) then 373 do j = 2, ncell 374 xr = x(k) - xi 375 yr = y(k) - yi 376 zr = z(k) - zi 377 call imager (xr,yr,zr,j) 378 r2 = xr*xr + yr* yr + zr*zr 379 if (.not. (use_polymer .and. r2.le.polycut2)) 380 & mscale(k) = 1.0d0 381 if (r2 .le. off2) then 382 r = sqrt(r2) 383 ck = rpole(1,kk) 384 dkx = rpole(2,kk) 385 dky = rpole(3,kk) 386 dkz = rpole(4,kk) 387 qkxx = rpole(5,kk) 388 qkxy = rpole(6,kk) 389 qkxz = rpole(7,kk) 390 qkyy = rpole(9,kk) 391 qkyz = rpole(10,kk) 392 qkzz = rpole(13,kk) 393c 394c intermediates involving moments and separation distance 395c 396 dir = dix*xr + diy*yr + diz*zr 397 qix = qixx*xr + qixy*yr + qixz*zr 398 qiy = qixy*xr + qiyy*yr + qiyz*zr 399 qiz = qixz*xr + qiyz*yr + qizz*zr 400 qir = qix*xr + qiy*yr + qiz*zr 401 dkr = dkx*xr + dky*yr + dkz*zr 402 qkx = qkxx*xr + qkxy*yr + qkxz*zr 403 qky = qkxy*xr + qkyy*yr + qkyz*zr 404 qkz = qkxz*xr + qkyz*yr + qkzz*zr 405 qkr = qkx*xr + qky*yr + qkz*zr 406 dik = dix*dkx + diy*dky + diz*dkz 407 qik = qix*qkx + qiy*qky + qiz*qkz 408 diqk = dix*qkx + diy*qky + diz*qkz 409 dkqi = dkx*qix + dky*qiy + dkz*qiz 410 qiqk = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 411 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 412c 413c get reciprocal distance terms for this interaction 414c 415 rr1 = f * mscale(k) / r 416 rr3 = rr1 / r2 417 rr5 = 3.0d0 * rr3 / r2 418 rr7 = 5.0d0 * rr5 / r2 419 rr9 = 7.0d0 * rr7 / r2 420c 421c find damped multipole intermediates and energy value 422c 423 if (use_chgpen) then 424 corek = pcore(kk) 425 valk = pval(kk) 426 alphak = palpha(kk) 427 term1 = corei*corek 428 term1i = corek*vali 429 term2i = corek*dir 430 term3i = corek*qir 431 term1k = corei*valk 432 term2k = -corei*dkr 433 term3k = corei*qkr 434 term1ik = vali*valk 435 term2ik = valk*dir - vali*dkr + dik 436 term3ik = vali*qkr + valk*qir - dir*dkr 437 & + 2.0d0*(dkqi-diqk+qiqk) 438 term4ik = dir*qkr - dkr*qir - 4.0d0*qik 439 term5ik = qir*qkr 440 call damppole (r,9,alphai,alphak, 441 & dmpi,dmpk,dmpik) 442 rr1i = dmpi(1)*rr1 443 rr3i = dmpi(3)*rr3 444 rr5i = dmpi(5)*rr5 445 rr1k = dmpk(1)*rr1 446 rr3k = dmpk(3)*rr3 447 rr5k = dmpk(5)*rr5 448 rr1ik = dmpik(1)*rr1 449 rr3ik = dmpik(3)*rr3 450 rr5ik = dmpik(5)*rr5 451 rr7ik = dmpik(7)*rr7 452 rr9ik = dmpik(9)*rr9 453 e = term1*rr1 + term1i*rr1i 454 & + term1k*rr1k + term1ik*rr1ik 455 & + term2i*rr3i + term2k*rr3k 456 & + term2ik*rr3ik + term3i*rr5i 457 & + term3k*rr5k + term3ik*rr5ik 458 & + term4ik*rr7ik + term5ik*rr9ik 459c 460c find standard multipole intermediates and energy value 461c 462 else 463 term1 = ci*ck 464 term2 = ck*dir - ci*dkr + dik 465 term3 = ci*qkr + ck*qir - dir*dkr 466 & + 2.0d0*(dkqi-diqk+qiqk) 467 term4 = dir*qkr - dkr*qir - 4.0d0*qik 468 term5 = qir*qkr 469 e = term1*rr1 + term2*rr3 + term3*rr5 470 & + term4*rr7 + term5*rr9 471 end if 472c 473c compute the energy contribution for this interaction 474c 475 if (use_group) e = e * fgrp 476 if (i .eq. k) e = 0.5d0 * e 477 em = em + e 478 end if 479 end do 480 end if 481 end do 482c 483c reset exclusion coefficients for connected atoms 484c 485 do j = 1, n12(i) 486 mscale(i12(j,i)) = 1.0d0 487 end do 488 do j = 1, n13(i) 489 mscale(i13(j,i)) = 1.0d0 490 end do 491 do j = 1, n14(i) 492 mscale(i14(j,i)) = 1.0d0 493 end do 494 do j = 1, n15(i) 495 mscale(i15(j,i)) = 1.0d0 496 end do 497 end do 498 end if 499c 500c perform deallocation of some local arrays 501c 502 deallocate (mscale) 503 return 504 end 505c 506c 507c ############################################################### 508c ## ## 509c ## subroutine empole0b -- neighbor list multipole energy ## 510c ## ## 511c ############################################################### 512c 513c 514c "empole0b" calculates the atomic multipole interaction energy 515c using a neighbor list 516c 517c 518 subroutine empole0b 519 use atoms 520 use bound 521 use chgpen 522 use chgpot 523 use couple 524 use energi 525 use group 526 use math 527 use mplpot 528 use mpole 529 use neigh 530 use shunt 531 use usage 532 implicit none 533 integer i,j,k 534 integer ii,kk,kkk 535 integer ix,iy,iz 536 integer kx,ky,kz 537 real*8 e,f,fgrp 538 real*8 xi,yi,zi 539 real*8 xr,yr,zr 540 real*8 r,r2,rr1,rr3 541 real*8 rr5,rr7,rr9 542 real*8 rr1i,rr3i,rr5i 543 real*8 rr1k,rr3k,rr5k 544 real*8 rr1ik,rr3ik,rr5ik 545 real*8 rr7ik,rr9ik 546 real*8 ci,dix,diy,diz 547 real*8 qixx,qixy,qixz 548 real*8 qiyy,qiyz,qizz 549 real*8 ck,dkx,dky,dkz 550 real*8 qkxx,qkxy,qkxz 551 real*8 qkyy,qkyz,qkzz 552 real*8 dir,dkr,dik,qik 553 real*8 qix,qiy,qiz,qir 554 real*8 qkx,qky,qkz,qkr 555 real*8 diqk,dkqi,qiqk 556 real*8 corei,corek 557 real*8 vali,valk 558 real*8 alphai,alphak 559 real*8 term1,term2,term3 560 real*8 term4,term5 561 real*8 term1i,term2i,term3i 562 real*8 term1k,term2k,term3k 563 real*8 term1ik,term2ik,term3ik 564 real*8 term4ik,term5ik 565 real*8 dmpi(9),dmpk(9) 566 real*8 dmpik(9) 567 real*8, allocatable :: mscale(:) 568 logical proceed,usei,usek 569 character*6 mode 570c 571c 572c zero out the total atomic multipole energy 573c 574 em = 0.0d0 575 if (npole .eq. 0) return 576c 577c check the sign of multipole components at chiral sites 578c 579 call chkpole 580c 581c rotate the multipole components into the global frame 582c 583 call rotpole 584c 585c perform dynamic allocation of some local arrays 586c 587 allocate (mscale(n)) 588c 589c initialize connected atom exclusion coefficients 590c 591 do i = 1, n 592 mscale(i) = 1.0d0 593 end do 594c 595c set conversion factor, cutoff and switching coefficients 596c 597 f = electric / dielec 598 mode = 'MPOLE' 599 call switch (mode) 600c 601c OpenMP directives for the major loop structure 602c 603!$OMP PARALLEL default(private) 604!$OMP& shared(npole,ipole,x,y,z,xaxis,yaxis,zaxis,rpole,pcore,pval, 605!$OMP& palpha,use,n12,i12,n13,i13,n14,i14,n15,i15,m2scale,m3scale, 606!$OMP& m4scale,m5scale,f,nelst,elst,use_chgpen,use_group,use_intra, 607!$OMP& use_bounds,off2) 608!$OMP& firstprivate(mscale) shared (em) 609!$OMP DO reduction(+:em) schedule(guided) 610c 611c compute the real space portion of the Ewald summation 612c 613 do ii = 1, npole 614 i = ipole(ii) 615 iz = zaxis(ii) 616 ix = xaxis(ii) 617 iy = abs(yaxis(ii)) 618 xi = x(i) 619 yi = y(i) 620 zi = z(i) 621 ci = rpole(1,ii) 622 dix = rpole(2,ii) 623 diy = rpole(3,ii) 624 diz = rpole(4,ii) 625 qixx = rpole(5,ii) 626 qixy = rpole(6,ii) 627 qixz = rpole(7,ii) 628 qiyy = rpole(9,ii) 629 qiyz = rpole(10,ii) 630 qizz = rpole(13,ii) 631 if (use_chgpen) then 632 corei = pcore(ii) 633 vali = pval(ii) 634 alphai = palpha(ii) 635 end if 636 usei = (use(i) .or. use(iz) .or. use(ix) .or. use(iy)) 637c 638c set exclusion coefficients for connected atoms 639c 640 do j = 1, n12(i) 641 mscale(i12(j,i)) = m2scale 642 end do 643 do j = 1, n13(i) 644 mscale(i13(j,i)) = m3scale 645 end do 646 do j = 1, n14(i) 647 mscale(i14(j,i)) = m4scale 648 end do 649 do j = 1, n15(i) 650 mscale(i15(j,i)) = m5scale 651 end do 652c 653c evaluate all sites within the cutoff distance 654c 655 do kkk = 1, nelst(ii) 656 kk = elst(kkk,ii) 657 k = ipole(kk) 658 kz = zaxis(kk) 659 kx = xaxis(kk) 660 ky = abs(yaxis(kk)) 661 usek = (use(k) .or. use(kz) .or. use(kx) .or. use(ky)) 662 proceed = .true. 663 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 664 if (.not. use_intra) proceed = .true. 665 if (proceed) proceed = (usei .or. usek) 666 if (proceed) then 667 xr = x(k) - xi 668 yr = y(k) - yi 669 zr = z(k) - zi 670 if (use_bounds) call image (xr,yr,zr) 671 r2 = xr*xr + yr* yr + zr*zr 672 if (r2 .le. off2) then 673 r = sqrt(r2) 674 ck = rpole(1,kk) 675 dkx = rpole(2,kk) 676 dky = rpole(3,kk) 677 dkz = rpole(4,kk) 678 qkxx = rpole(5,kk) 679 qkxy = rpole(6,kk) 680 qkxz = rpole(7,kk) 681 qkyy = rpole(9,kk) 682 qkyz = rpole(10,kk) 683 qkzz = rpole(13,kk) 684c 685c intermediates involving moments and separation distance 686c 687 dir = dix*xr + diy*yr + diz*zr 688 qix = qixx*xr + qixy*yr + qixz*zr 689 qiy = qixy*xr + qiyy*yr + qiyz*zr 690 qiz = qixz*xr + qiyz*yr + qizz*zr 691 qir = qix*xr + qiy*yr + qiz*zr 692 dkr = dkx*xr + dky*yr + dkz*zr 693 qkx = qkxx*xr + qkxy*yr + qkxz*zr 694 qky = qkxy*xr + qkyy*yr + qkyz*zr 695 qkz = qkxz*xr + qkyz*yr + qkzz*zr 696 qkr = qkx*xr + qky*yr + qkz*zr 697 dik = dix*dkx + diy*dky + diz*dkz 698 qik = qix*qkx + qiy*qky + qiz*qkz 699 diqk = dix*qkx + diy*qky + diz*qkz 700 dkqi = dkx*qix + dky*qiy + dkz*qiz 701 qiqk = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 702 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 703c 704c get reciprocal distance terms for this interaction 705c 706 rr1 = f * mscale(k) / r 707 rr3 = rr1 / r2 708 rr5 = 3.0d0 * rr3 / r2 709 rr7 = 5.0d0 * rr5 / r2 710 rr9 = 7.0d0 * rr7 / r2 711c 712c find damped multipole intermediates and energy value 713c 714 if (use_chgpen) then 715 corek = pcore(kk) 716 valk = pval(kk) 717 alphak = palpha(kk) 718 term1 = corei*corek 719 term1i = corek*vali 720 term2i = corek*dir 721 term3i = corek*qir 722 term1k = corei*valk 723 term2k = -corei*dkr 724 term3k = corei*qkr 725 term1ik = vali*valk 726 term2ik = valk*dir - vali*dkr + dik 727 term3ik = vali*qkr + valk*qir - dir*dkr 728 & + 2.0d0*(dkqi-diqk+qiqk) 729 term4ik = dir*qkr - dkr*qir - 4.0d0*qik 730 term5ik = qir*qkr 731 call damppole (r,9,alphai,alphak, 732 & dmpi,dmpk,dmpik) 733 rr1i = dmpi(1)*rr1 734 rr3i = dmpi(3)*rr3 735 rr5i = dmpi(5)*rr5 736 rr1k = dmpk(1)*rr1 737 rr3k = dmpk(3)*rr3 738 rr5k = dmpk(5)*rr5 739 rr1ik = dmpik(1)*rr1 740 rr3ik = dmpik(3)*rr3 741 rr5ik = dmpik(5)*rr5 742 rr7ik = dmpik(7)*rr7 743 rr9ik = dmpik(9)*rr9 744 e = term1*rr1 + term1i*rr1i 745 & + term1k*rr1k + term1ik*rr1ik 746 & + term2i*rr3i + term2k*rr3k 747 & + term2ik*rr3ik + term3i*rr5i 748 & + term3k*rr5k + term3ik*rr5ik 749 & + term4ik*rr7ik + term5ik*rr9ik 750c 751c find standard multipole intermediates and energy value 752c 753 else 754 term1 = ci*ck 755 term2 = ck*dir - ci*dkr + dik 756 term3 = ci*qkr + ck*qir - dir*dkr 757 & + 2.0d0*(dkqi-diqk+qiqk) 758 term4 = dir*qkr - dkr*qir - 4.0d0*qik 759 term5 = qir*qkr 760 e = term1*rr1 + term2*rr3 + term3*rr5 761 & + term4*rr7 + term5*rr9 762 end if 763c 764c compute the energy contribution for this interaction 765c 766 if (use_group) e = e * fgrp 767 em = em + e 768 end if 769 end if 770 end do 771c 772c reset exclusion coefficients for connected atoms 773c 774 do j = 1, n12(i) 775 mscale(i12(j,i)) = 1.0d0 776 end do 777 do j = 1, n13(i) 778 mscale(i13(j,i)) = 1.0d0 779 end do 780 do j = 1, n14(i) 781 mscale(i14(j,i)) = 1.0d0 782 end do 783 do j = 1, n15(i) 784 mscale(i15(j,i)) = 1.0d0 785 end do 786 end do 787c 788c OpenMP directives for the major loop structure 789c 790!$OMP END DO 791!$OMP END PARALLEL 792c 793c perform deallocation of some local arrays 794c 795 deallocate (mscale) 796 return 797 end 798c 799c 800c ################################################################ 801c ## ## 802c ## subroutine empole0c -- Ewald multipole energy via loop ## 803c ## ## 804c ################################################################ 805c 806c 807c "empole0c" calculates the atomic multipole interaction energy 808c using particle mesh Ewald summation and a double loop 809c 810c 811 subroutine empole0c 812 use atoms 813 use boxes 814 use chgpot 815 use energi 816 use ewald 817 use math 818 use mpole 819 use pme 820 implicit none 821 integer i,ii 822 real*8 e,f 823 real*8 term,fterm 824 real*8 cii,dii,qii 825 real*8 xd,yd,zd 826 real*8 ci,dix,diy,diz 827 real*8 qixx,qixy,qixz 828 real*8 qiyy,qiyz,qizz 829c 830c 831c zero out the total atomic multipole energy 832c 833 em = 0.0d0 834 if (npole .eq. 0) return 835c 836c set grid size, spline order and Ewald coefficient 837c 838 nfft1 = nefft1 839 nfft2 = nefft2 840 nfft3 = nefft3 841 bsorder = bseorder 842 aewald = aeewald 843c 844c set the energy unit conversion factor 845c 846 f = electric / dielec 847c 848c check the sign of multipole components at chiral sites 849c 850 call chkpole 851c 852c rotate the multipole components into the global frame 853c 854 call rotpole 855c 856c compute the real space portion of the Ewald summation 857c 858 call emreal0c 859c 860c compute the reciprocal space part of the Ewald summation 861c 862 call emrecip 863c 864c compute the self-energy portion of the Ewald summation 865c 866 term = 2.0d0 * aewald * aewald 867 fterm = -f * aewald / rootpi 868 do ii = 1, npole 869 ci = rpole(1,ii) 870 dix = rpole(2,ii) 871 diy = rpole(3,ii) 872 diz = rpole(4,ii) 873 qixx = rpole(5,ii) 874 qixy = rpole(6,ii) 875 qixz = rpole(7,ii) 876 qiyy = rpole(9,ii) 877 qiyz = rpole(10,ii) 878 qizz = rpole(13,ii) 879 cii = ci*ci 880 dii = dix*dix + diy*diy + diz*diz 881 qii = 2.0d0*(qixy*qixy+qixz*qixz+qiyz*qiyz) 882 & + qixx*qixx + qiyy*qiyy + qizz*qizz 883 e = fterm * (cii + term*(dii/3.0d0+2.0d0*term*qii/5.0d0)) 884 em = em + e 885 end do 886c 887c compute the cell dipole boundary correction term 888c 889 if (boundary .eq. 'VACUUM') then 890 xd = 0.0d0 891 yd = 0.0d0 892 zd = 0.0d0 893 do ii = 1, npole 894 i = ipole(ii) 895 dix = rpole(2,ii) 896 diy = rpole(3,ii) 897 diz = rpole(4,ii) 898 xd = xd + dix + rpole(1,ii)*x(i) 899 yd = yd + diy + rpole(1,ii)*y(i) 900 zd = zd + diz + rpole(1,ii)*z(i) 901 end do 902 term = (2.0d0/3.0d0) * f * (pi/volbox) 903 e = term * (xd*xd+yd*yd+zd*zd) 904 em = em + e 905 end if 906 return 907 end 908c 909c 910c ################################################################# 911c ## ## 912c ## subroutine emreal0c -- real space mpole energy via loop ## 913c ## ## 914c ################################################################# 915c 916c 917c "emreal0c" evaluates the real space portion of the Ewald sum 918c energy due to atomic multipoles using a double loop 919c 920c literature reference: 921c 922c W. Smith, "Point Multipoles in the Ewald Summation (Revisited)", 923c CCP5 Newsletter, 46, 18-30, 1998 [newsletters are available at 924c https://www.ccp5.ac.uk/infoweb/newsletters] 925c 926c 927 subroutine emreal0c 928 use atoms 929 use bound 930 use cell 931 use chgpen 932 use chgpot 933 use couple 934 use energi 935 use math 936 use mplpot 937 use mpole 938 use shunt 939 implicit none 940 integer i,j,k 941 integer ii,kk 942 integer jcell 943 real*8 e,f,scalek 944 real*8 xi,yi,zi 945 real*8 xr,yr,zr 946 real*8 r,r2,rr1,rr3 947 real*8 rr5,rr7,rr9 948 real*8 rr1i,rr3i,rr5i 949 real*8 rr1k,rr3k,rr5k 950 real*8 rr1ik,rr3ik,rr5ik 951 real*8 rr7ik,rr9ik 952 real*8 ci,dix,diy,diz 953 real*8 qixx,qixy,qixz 954 real*8 qiyy,qiyz,qizz 955 real*8 ck,dkx,dky,dkz 956 real*8 qkxx,qkxy,qkxz 957 real*8 qkyy,qkyz,qkzz 958 real*8 dir,dkr,dik,qik 959 real*8 qix,qiy,qiz,qir 960 real*8 qkx,qky,qkz,qkr 961 real*8 diqk,dkqi,qiqk 962 real*8 corei,corek 963 real*8 vali,valk 964 real*8 alphai,alphak 965 real*8 term1,term2,term3 966 real*8 term4,term5 967 real*8 term1i,term2i,term3i 968 real*8 term1k,term2k,term3k 969 real*8 term1ik,term2ik,term3ik 970 real*8 term4ik,term5ik 971 real*8 dmpi(9),dmpk(9) 972 real*8 dmpik(9),dmpe(9) 973 real*8, allocatable :: mscale(:) 974 character*6 mode 975c 976c 977c perform dynamic allocation of some local arrays 978c 979 allocate (mscale(n)) 980c 981c initialize connected atom exclusion coefficients 982c 983 do i = 1, n 984 mscale(i) = 1.0d0 985 end do 986c 987c set conversion factor, cutoff and switching coefficients 988c 989 f = electric / dielec 990 mode = 'EWALD' 991 call switch (mode) 992c 993c compute the real space portion of the Ewald summation 994c 995 do ii = 1, npole-1 996 i = ipole(ii) 997 xi = x(i) 998 yi = y(i) 999 zi = z(i) 1000 ci = rpole(1,ii) 1001 dix = rpole(2,ii) 1002 diy = rpole(3,ii) 1003 diz = rpole(4,ii) 1004 qixx = rpole(5,ii) 1005 qixy = rpole(6,ii) 1006 qixz = rpole(7,ii) 1007 qiyy = rpole(9,ii) 1008 qiyz = rpole(10,ii) 1009 qizz = rpole(13,ii) 1010 if (use_chgpen) then 1011 corei = pcore(ii) 1012 vali = pval(ii) 1013 alphai = palpha(ii) 1014 end if 1015c 1016c set exclusion coefficients for connected atoms 1017c 1018 do j = 1, n12(i) 1019 mscale(i12(j,i)) = m2scale 1020 end do 1021 do j = 1, n13(i) 1022 mscale(i13(j,i)) = m3scale 1023 end do 1024 do j = 1, n14(i) 1025 mscale(i14(j,i)) = m4scale 1026 end do 1027 do j = 1, n15(i) 1028 mscale(i15(j,i)) = m5scale 1029 end do 1030c 1031c evaluate all sites within the cutoff distance 1032c 1033 do kk = ii+1, npole 1034 k = ipole(kk) 1035 xr = x(k) - xi 1036 yr = y(k) - yi 1037 zr = z(k) - zi 1038 if (use_bounds) call image (xr,yr,zr) 1039 r2 = xr*xr + yr* yr + zr*zr 1040 if (r2 .le. off2) then 1041 r = sqrt(r2) 1042 ck = rpole(1,kk) 1043 dkx = rpole(2,kk) 1044 dky = rpole(3,kk) 1045 dkz = rpole(4,kk) 1046 qkxx = rpole(5,kk) 1047 qkxy = rpole(6,kk) 1048 qkxz = rpole(7,kk) 1049 qkyy = rpole(9,kk) 1050 qkyz = rpole(10,kk) 1051 qkzz = rpole(13,kk) 1052c 1053c intermediates involving moments and separation distance 1054c 1055 dir = dix*xr + diy*yr + diz*zr 1056 qix = qixx*xr + qixy*yr + qixz*zr 1057 qiy = qixy*xr + qiyy*yr + qiyz*zr 1058 qiz = qixz*xr + qiyz*yr + qizz*zr 1059 qir = qix*xr + qiy*yr + qiz*zr 1060 dkr = dkx*xr + dky*yr + dkz*zr 1061 qkx = qkxx*xr + qkxy*yr + qkxz*zr 1062 qky = qkxy*xr + qkyy*yr + qkyz*zr 1063 qkz = qkxz*xr + qkyz*yr + qkzz*zr 1064 qkr = qkx*xr + qky*yr + qkz*zr 1065 dik = dix*dkx + diy*dky + diz*dkz 1066 qik = qix*qkx + qiy*qky + qiz*qkz 1067 diqk = dix*qkx + diy*qky + diz*qkz 1068 dkqi = dkx*qix + dky*qiy + dkz*qiz 1069 qiqk = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 1070 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 1071c 1072c get reciprocal distance terms for this interaction 1073c 1074 rr1 = f / r 1075 rr3 = rr1 / r2 1076 rr5 = 3.0d0 * rr3 / r2 1077 rr7 = 5.0d0 * rr5 / r2 1078 rr9 = 7.0d0 * rr7 / r2 1079c 1080c calculate real space Ewald error function damping 1081c 1082 call dampewald (9,r,r2,f,dmpe) 1083c 1084c find damped multipole intermediates and energy value 1085c 1086 if (use_chgpen) then 1087 corek = pcore(kk) 1088 valk = pval(kk) 1089 alphak = palpha(kk) 1090 term1 = corei*corek 1091 term1i = corek*vali 1092 term2i = corek*dir 1093 term3i = corek*qir 1094 term1k = corei*valk 1095 term2k = -corei*dkr 1096 term3k = corei*qkr 1097 term1ik = vali*valk 1098 term2ik = valk*dir - vali*dkr + dik 1099 term3ik = vali*qkr + valk*qir - dir*dkr 1100 & + 2.0d0*(dkqi-diqk+qiqk) 1101 term4ik = dir*qkr - dkr*qir - 4.0d0*qik 1102 term5ik = qir*qkr 1103 call damppole (r,9,alphai,alphak, 1104 & dmpi,dmpk,dmpik) 1105 scalek = mscale(k) 1106 scalek = mscale(k) 1107 rr1i = dmpe(1) - (1.0d0-scalek*dmpi(1))*rr1 1108 rr3i = dmpe(3) - (1.0d0-scalek*dmpi(3))*rr3 1109 rr5i = dmpe(5) - (1.0d0-scalek*dmpi(5))*rr5 1110 rr1k = dmpe(1) - (1.0d0-scalek*dmpk(1))*rr1 1111 rr3k = dmpe(3) - (1.0d0-scalek*dmpk(3))*rr3 1112 rr5k = dmpe(5) - (1.0d0-scalek*dmpk(5))*rr5 1113 rr1ik = dmpe(1) - (1.0d0-scalek*dmpik(1))*rr1 1114 rr3ik = dmpe(3) - (1.0d0-scalek*dmpik(3))*rr3 1115 rr5ik = dmpe(5) - (1.0d0-scalek*dmpik(5))*rr5 1116 rr7ik = dmpe(7) - (1.0d0-scalek*dmpik(7))*rr7 1117 rr9ik = dmpe(9) - (1.0d0-scalek*dmpik(9))*rr9 1118 rr1 = dmpe(1) - (1.0d0-scalek)*rr1 1119 e = term1*rr1 + term4ik*rr7ik + term5ik*rr9ik 1120 & + term1i*rr1i + term1k*rr1k + term1ik*rr1ik 1121 & + term2i*rr3i + term2k*rr3k + term2ik*rr3ik 1122 & + term3i*rr5i + term3k*rr5k + term3ik*rr5ik 1123c 1124c find standard multipole intermediates and energy value 1125c 1126 else 1127 term1 = ci*ck 1128 term2 = ck*dir - ci*dkr + dik 1129 term3 = ci*qkr + ck*qir - dir*dkr 1130 & + 2.0d0*(dkqi-diqk+qiqk) 1131 term4 = dir*qkr - dkr*qir - 4.0d0*qik 1132 term5 = qir*qkr 1133 scalek = 1.0d0 - mscale(k) 1134 rr1 = dmpe(1) - scalek*rr1 1135 rr3 = dmpe(3) - scalek*rr3 1136 rr5 = dmpe(5) - scalek*rr5 1137 rr7 = dmpe(7) - scalek*rr7 1138 rr9 = dmpe(9) - scalek*rr9 1139 e = term1*rr1 + term2*rr3 + term3*rr5 1140 & + term4*rr7 + term5*rr9 1141 end if 1142c 1143c increment the overall multipole energy component 1144c 1145 em = em + e 1146 end if 1147 end do 1148c 1149c reset exclusion coefficients for connected atoms 1150c 1151 do j = 1, n12(i) 1152 mscale(i12(j,i)) = 1.0d0 1153 end do 1154 do j = 1, n13(i) 1155 mscale(i13(j,i)) = 1.0d0 1156 end do 1157 do j = 1, n14(i) 1158 mscale(i14(j,i)) = 1.0d0 1159 end do 1160 do j = 1, n15(i) 1161 mscale(i15(j,i)) = 1.0d0 1162 end do 1163 end do 1164c 1165c for periodic boundary conditions with large cutoffs 1166c neighbors must be found by the replicates method 1167c 1168 if (use_replica) then 1169c 1170c calculate interaction energy with other unit cells 1171c 1172 do ii = 1, npole 1173 i = ipole(ii) 1174 xi = x(i) 1175 yi = y(i) 1176 zi = z(i) 1177 ci = rpole(1,ii) 1178 dix = rpole(2,ii) 1179 diy = rpole(3,ii) 1180 diz = rpole(4,ii) 1181 qixx = rpole(5,ii) 1182 qixy = rpole(6,ii) 1183 qixz = rpole(7,ii) 1184 qiyy = rpole(9,ii) 1185 qiyz = rpole(10,ii) 1186 qizz = rpole(13,ii) 1187 if (use_chgpen) then 1188 corei = pcore(ii) 1189 vali = pval(ii) 1190 alphai = palpha(ii) 1191 end if 1192c 1193c set exclusion coefficients for connected atoms 1194c 1195 do j = 1, n12(i) 1196 mscale(i12(j,i)) = m2scale 1197 end do 1198 do j = 1, n13(i) 1199 mscale(i13(j,i)) = m3scale 1200 end do 1201 do j = 1, n14(i) 1202 mscale(i14(j,i)) = m4scale 1203 end do 1204 do j = 1, n15(i) 1205 mscale(i15(j,i)) = m5scale 1206 end do 1207c 1208c evaluate all sites within the cutoff distance 1209c 1210 do kk = ii, npole 1211 k = ipole(kk) 1212 do jcell = 2, ncell 1213 xr = x(k) - xi 1214 yr = y(k) - yi 1215 zr = z(k) - zi 1216 call imager (xr,yr,zr,jcell) 1217 r2 = xr*xr + yr* yr + zr*zr 1218 if (.not. (use_polymer .and. r2.le.polycut2)) 1219 & mscale(k) = 1.0d0 1220 if (r2 .le. off2) then 1221 r = sqrt(r2) 1222 ck = rpole(1,kk) 1223 dkx = rpole(2,kk) 1224 dky = rpole(3,kk) 1225 dkz = rpole(4,kk) 1226 qkxx = rpole(5,kk) 1227 qkxy = rpole(6,kk) 1228 qkxz = rpole(7,kk) 1229 qkyy = rpole(9,kk) 1230 qkyz = rpole(10,kk) 1231 qkzz = rpole(13,kk) 1232c 1233c intermediates involving moments and separation distance 1234c 1235 dir = dix*xr + diy*yr + diz*zr 1236 qix = qixx*xr + qixy*yr + qixz*zr 1237 qiy = qixy*xr + qiyy*yr + qiyz*zr 1238 qiz = qixz*xr + qiyz*yr + qizz*zr 1239 qir = qix*xr + qiy*yr + qiz*zr 1240 dkr = dkx*xr + dky*yr + dkz*zr 1241 qkx = qkxx*xr + qkxy*yr + qkxz*zr 1242 qky = qkxy*xr + qkyy*yr + qkyz*zr 1243 qkz = qkxz*xr + qkyz*yr + qkzz*zr 1244 qkr = qkx*xr + qky*yr + qkz*zr 1245 dik = dix*dkx + diy*dky + diz*dkz 1246 qik = qix*qkx + qiy*qky + qiz*qkz 1247 diqk = dix*qkx + diy*qky + diz*qkz 1248 dkqi = dkx*qix + dky*qiy + dkz*qiz 1249 qiqk = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 1250 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 1251c 1252c get reciprocal distance terms for this interaction 1253c 1254 rr1 = f / r 1255 rr3 = rr1 / r2 1256 rr5 = 3.0d0 * rr3 / r2 1257 rr7 = 5.0d0 * rr5 / r2 1258 rr9 = 7.0d0 * rr7 / r2 1259c 1260c calculate real space Ewald error function damping 1261c 1262 call dampewald (9,r,r2,f,dmpe) 1263c 1264c find damped multipole intermediates and energy value 1265c 1266 if (use_chgpen) then 1267 corek = pcore(kk) 1268 valk = pval(kk) 1269 alphak = palpha(kk) 1270 term1 = corei*corek 1271 term1i = corek*vali 1272 term2i = corek*dir 1273 term3i = corek*qir 1274 term1k = corei*valk 1275 term2k = -corei*dkr 1276 term3k = corei*qkr 1277 term1ik = vali*valk 1278 term2ik = valk*dir - vali*dkr + dik 1279 term3ik = vali*qkr + valk*qir - dir*dkr 1280 & + 2.0d0*(dkqi-diqk+qiqk) 1281 term4ik = dir*qkr - dkr*qir - 4.0d0*qik 1282 term5ik = qir*qkr 1283 call damppole (r,9,alphai,alphak, 1284 & dmpi,dmpk,dmpik) 1285 scalek = mscale(k) 1286 rr1i = dmpe(1) - (1.0d0-scalek*dmpi(1))*rr1 1287 rr3i = dmpe(3) - (1.0d0-scalek*dmpi(3))*rr3 1288 rr5i = dmpe(5) - (1.0d0-scalek*dmpi(5))*rr5 1289 rr1k = dmpe(1) - (1.0d0-scalek*dmpk(1))*rr1 1290 rr3k = dmpe(3) - (1.0d0-scalek*dmpk(3))*rr3 1291 rr5k = dmpe(5) - (1.0d0-scalek*dmpk(5))*rr5 1292 rr1ik = dmpe(1) - (1.0d0-scalek*dmpik(1))*rr1 1293 rr3ik = dmpe(3) - (1.0d0-scalek*dmpik(3))*rr3 1294 rr5ik = dmpe(5) - (1.0d0-scalek*dmpik(5))*rr5 1295 rr7ik = dmpe(7) - (1.0d0-scalek*dmpik(7))*rr7 1296 rr9ik = dmpe(9) - (1.0d0-scalek*dmpik(9))*rr9 1297 rr1 = dmpe(1) - (1.0d0-scalek)*rr1 1298 e = term1*rr1 + term1i*rr1i 1299 & + term1k*rr1k + term1ik*rr1ik 1300 & + term2i*rr3i + term2k*rr3k 1301 & + term2ik*rr3ik + term3i*rr5i 1302 & + term3k*rr5k + term3ik*rr5ik 1303 & + term4ik*rr7ik + term5ik*rr9ik 1304c 1305c find standard multipole intermediates and energy value 1306c 1307 else 1308 term1 = ci*ck 1309 term2 = ck*dir - ci*dkr + dik 1310 term3 = ci*qkr + ck*qir - dir*dkr 1311 & + 2.0d0*(dkqi-diqk+qiqk) 1312 term4 = dir*qkr - dkr*qir - 4.0d0*qik 1313 term5 = qir*qkr 1314 scalek = 1.0d0 - mscale(k) 1315 rr1 = dmpe(1) - scalek*rr1 1316 rr3 = dmpe(3) - scalek*rr3 1317 rr5 = dmpe(5) - scalek*rr5 1318 rr7 = dmpe(7) - scalek*rr7 1319 rr9 = dmpe(9) - scalek*rr9 1320 e = term1*rr1 + term2*rr3 + term3*rr5 1321 & + term4*rr7 + term5*rr9 1322 end if 1323c 1324c increment the overall multipole energy component 1325c 1326 if (i .eq. k) e = 0.5d0 * e 1327 em = em + e 1328 end if 1329 end do 1330 end do 1331c 1332c reset exclusion coefficients for connected atoms 1333c 1334 do j = 1, n12(i) 1335 mscale(i12(j,i)) = 1.0d0 1336 end do 1337 do j = 1, n13(i) 1338 mscale(i13(j,i)) = 1.0d0 1339 end do 1340 do j = 1, n14(i) 1341 mscale(i14(j,i)) = 1.0d0 1342 end do 1343 do j = 1, n15(i) 1344 mscale(i15(j,i)) = 1.0d0 1345 end do 1346 end do 1347 end if 1348c 1349c perform deallocation of some local arrays 1350c 1351 deallocate (mscale) 1352 return 1353 end 1354c 1355c 1356c ################################################################ 1357c ## ## 1358c ## subroutine empole0d -- Ewald multipole energy via list ## 1359c ## ## 1360c ################################################################ 1361c 1362c 1363c "empole0d" calculates the atomic multipole interaction energy 1364c using particle mesh Ewald summation and a neighbor list 1365c 1366c 1367 subroutine empole0d 1368 use atoms 1369 use boxes 1370 use chgpot 1371 use energi 1372 use ewald 1373 use math 1374 use mpole 1375 use pme 1376 implicit none 1377 integer i,ii 1378 real*8 e,f 1379 real*8 term,fterm 1380 real*8 cii,dii,qii 1381 real*8 xd,yd,zd 1382 real*8 ci,dix,diy,diz 1383 real*8 qixx,qixy,qixz 1384 real*8 qiyy,qiyz,qizz 1385c 1386c 1387c zero out the total atomic multipole energy 1388c 1389 em = 0.0d0 1390 if (npole .eq. 0) return 1391c 1392c set grid size, spline order and Ewald coefficient 1393c 1394 nfft1 = nefft1 1395 nfft2 = nefft2 1396 nfft3 = nefft3 1397 bsorder = bseorder 1398 aewald = aeewald 1399c 1400c set the energy unit conversion factor 1401c 1402 f = electric / dielec 1403c 1404c check the sign of multipole components at chiral sites 1405c 1406 call chkpole 1407c 1408c rotate the multipole components into the global frame 1409c 1410 call rotpole 1411c 1412c compute the real space portion of the Ewald summation 1413c 1414 call emreal0d 1415c 1416c compute the reciprocal space part of the Ewald summation 1417c 1418 call emrecip 1419c 1420c compute the self-energy portion of the Ewald summation 1421c 1422 term = 2.0d0 * aewald * aewald 1423 fterm = -f * aewald / rootpi 1424 do ii = 1, npole 1425 ci = rpole(1,ii) 1426 dix = rpole(2,ii) 1427 diy = rpole(3,ii) 1428 diz = rpole(4,ii) 1429 qixx = rpole(5,ii) 1430 qixy = rpole(6,ii) 1431 qixz = rpole(7,ii) 1432 qiyy = rpole(9,ii) 1433 qiyz = rpole(10,ii) 1434 qizz = rpole(13,ii) 1435 cii = ci*ci 1436 dii = dix*dix + diy*diy + diz*diz 1437 qii = 2.0d0*(qixy*qixy+qixz*qixz+qiyz*qiyz) 1438 & + qixx*qixx + qiyy*qiyy + qizz*qizz 1439 e = fterm * (cii + term*(dii/3.0d0+2.0d0*term*qii/5.0d0)) 1440 em = em + e 1441 end do 1442c 1443c compute the cell dipole boundary correction term 1444c 1445 if (boundary .eq. 'VACUUM') then 1446 xd = 0.0d0 1447 yd = 0.0d0 1448 zd = 0.0d0 1449 do ii = 1, npole 1450 i = ipole(ii) 1451 dix = rpole(2,ii) 1452 diy = rpole(3,ii) 1453 diz = rpole(4,ii) 1454 xd = xd + dix + rpole(1,ii)*x(i) 1455 yd = yd + diy + rpole(1,ii)*y(i) 1456 zd = zd + diz + rpole(1,ii)*z(i) 1457 end do 1458 term = (2.0d0/3.0d0) * f * (pi/volbox) 1459 e = term * (xd*xd+yd*yd+zd*zd) 1460 em = em + e 1461 end if 1462 return 1463 end 1464c 1465c 1466c ################################################################# 1467c ## ## 1468c ## subroutine emreal0d -- real space mpole energy via list ## 1469c ## ## 1470c ################################################################# 1471c 1472c 1473c "emreal0d" evaluates the real space portion of the Ewald sum 1474c energy due to atomic multipoles using a neighbor list 1475c 1476c literature reference: 1477c 1478c W. Smith, "Point Multipoles in the Ewald Summation (Revisited)", 1479c CCP5 Newsletter, 46, 18-30, 1998 [newsletters are available at 1480c https://www.ccp5.ac.uk/infoweb/newsletters] 1481c 1482c 1483 subroutine emreal0d 1484 use atoms 1485 use bound 1486 use chgpen 1487 use chgpot 1488 use couple 1489 use energi 1490 use math 1491 use mplpot 1492 use mpole 1493 use neigh 1494 use shunt 1495 implicit none 1496 integer i,j,k 1497 integer ii,kk,kkk 1498 real*8 e,f,scalek 1499 real*8 xi,yi,zi 1500 real*8 xr,yr,zr 1501 real*8 r,r2,rr1,rr3 1502 real*8 rr5,rr7,rr9 1503 real*8 rr1i,rr3i,rr5i 1504 real*8 rr1k,rr3k,rr5k 1505 real*8 rr1ik,rr3ik,rr5ik 1506 real*8 rr7ik,rr9ik 1507 real*8 ci,dix,diy,diz 1508 real*8 qixx,qixy,qixz 1509 real*8 qiyy,qiyz,qizz 1510 real*8 ck,dkx,dky,dkz 1511 real*8 qkxx,qkxy,qkxz 1512 real*8 qkyy,qkyz,qkzz 1513 real*8 dir,dkr,dik,qik 1514 real*8 qix,qiy,qiz,qir 1515 real*8 qkx,qky,qkz,qkr 1516 real*8 diqk,dkqi,qiqk 1517 real*8 corei,corek 1518 real*8 vali,valk 1519 real*8 alphai,alphak 1520 real*8 term1,term2,term3 1521 real*8 term4,term5 1522 real*8 term1i,term2i,term3i 1523 real*8 term1k,term2k,term3k 1524 real*8 term1ik,term2ik,term3ik 1525 real*8 term4ik,term5ik 1526 real*8 dmpi(9),dmpk(9) 1527 real*8 dmpik(9),dmpe(9) 1528 real*8, allocatable :: mscale(:) 1529 character*6 mode 1530c 1531c 1532c perform dynamic allocation of some local arrays 1533c 1534 allocate (mscale(n)) 1535c 1536c initialize connected atom exclusion coefficients 1537c 1538 do i = 1, n 1539 mscale(i) = 1.0d0 1540 end do 1541c 1542c set conversion factor, cutoff and switching coefficients 1543c 1544 f = electric / dielec 1545 mode = 'EWALD' 1546 call switch (mode) 1547c 1548c OpenMP directives for the major loop structure 1549c 1550!$OMP PARALLEL default(private) 1551!$OMP& shared(npole,ipole,x,y,z,rpole,pcore,pval,palpha,n12,i12, 1552!$OMP& n13,i13,n14,i14,n15,i15,m2scale,m3scale,m4scale,m5scale, 1553!$OMP& f,nelst,elst,use_bounds,use_chgpen,off2) 1554!$OMP& firstprivate(mscale) shared (em) 1555!$OMP DO reduction(+:em) schedule(guided) 1556c 1557c compute the real space portion of the Ewald summation 1558c 1559 do ii = 1, npole 1560 i = ipole(ii) 1561 xi = x(i) 1562 yi = y(i) 1563 zi = z(i) 1564 ci = rpole(1,ii) 1565 dix = rpole(2,ii) 1566 diy = rpole(3,ii) 1567 diz = rpole(4,ii) 1568 qixx = rpole(5,ii) 1569 qixy = rpole(6,ii) 1570 qixz = rpole(7,ii) 1571 qiyy = rpole(9,ii) 1572 qiyz = rpole(10,ii) 1573 qizz = rpole(13,ii) 1574 if (use_chgpen) then 1575 corei = pcore(ii) 1576 vali = pval(ii) 1577 alphai = palpha(ii) 1578 end if 1579c 1580c set exclusion coefficients for connected atoms 1581c 1582 do j = 1, n12(i) 1583 mscale(i12(j,i)) = m2scale 1584 end do 1585 do j = 1, n13(i) 1586 mscale(i13(j,i)) = m3scale 1587 end do 1588 do j = 1, n14(i) 1589 mscale(i14(j,i)) = m4scale 1590 end do 1591 do j = 1, n15(i) 1592 mscale(i15(j,i)) = m5scale 1593 end do 1594c 1595c evaluate all sites within the cutoff distance 1596c 1597 do kkk = 1, nelst(ii) 1598 kk = elst(kkk,ii) 1599 k = ipole(kk) 1600 xr = x(k) - xi 1601 yr = y(k) - yi 1602 zr = z(k) - zi 1603 if (use_bounds) call image (xr,yr,zr) 1604 r2 = xr*xr + yr* yr + zr*zr 1605 if (r2 .le. off2) then 1606 r = sqrt(r2) 1607 ck = rpole(1,kk) 1608 dkx = rpole(2,kk) 1609 dky = rpole(3,kk) 1610 dkz = rpole(4,kk) 1611 qkxx = rpole(5,kk) 1612 qkxy = rpole(6,kk) 1613 qkxz = rpole(7,kk) 1614 qkyy = rpole(9,kk) 1615 qkyz = rpole(10,kk) 1616 qkzz = rpole(13,kk) 1617c 1618c intermediates involving moments and separation distance 1619c 1620 dir = dix*xr + diy*yr + diz*zr 1621 qix = qixx*xr + qixy*yr + qixz*zr 1622 qiy = qixy*xr + qiyy*yr + qiyz*zr 1623 qiz = qixz*xr + qiyz*yr + qizz*zr 1624 qir = qix*xr + qiy*yr + qiz*zr 1625 dkr = dkx*xr + dky*yr + dkz*zr 1626 qkx = qkxx*xr + qkxy*yr + qkxz*zr 1627 qky = qkxy*xr + qkyy*yr + qkyz*zr 1628 qkz = qkxz*xr + qkyz*yr + qkzz*zr 1629 qkr = qkx*xr + qky*yr + qkz*zr 1630 dik = dix*dkx + diy*dky + diz*dkz 1631 qik = qix*qkx + qiy*qky + qiz*qkz 1632 diqk = dix*qkx + diy*qky + diz*qkz 1633 dkqi = dkx*qix + dky*qiy + dkz*qiz 1634 qiqk = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 1635 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 1636c 1637c get reciprocal distance terms for this interaction 1638c 1639 rr1 = f / r 1640 rr3 = rr1 / r2 1641 rr5 = 3.0d0 * rr3 / r2 1642 rr7 = 5.0d0 * rr5 / r2 1643 rr9 = 7.0d0 * rr7 / r2 1644c 1645c calculate real space Ewald error function damping 1646c 1647 call dampewald (9,r,r2,f,dmpe) 1648c 1649c find damped multipole intermediates and energy value 1650c 1651 if (use_chgpen) then 1652 corek = pcore(kk) 1653 valk = pval(kk) 1654 alphak = palpha(kk) 1655 term1 = corei*corek 1656 term1i = corek*vali 1657 term2i = corek*dir 1658 term3i = corek*qir 1659 term1k = corei*valk 1660 term2k = -corei*dkr 1661 term3k = corei*qkr 1662 term1ik = vali*valk 1663 term2ik = valk*dir - vali*dkr + dik 1664 term3ik = vali*qkr + valk*qir - dir*dkr 1665 & + 2.0d0*(dkqi-diqk+qiqk) 1666 term4ik = dir*qkr - dkr*qir - 4.0d0*qik 1667 term5ik = qir*qkr 1668 call damppole (r,9,alphai,alphak, 1669 & dmpi,dmpk,dmpik) 1670 scalek = mscale(k) 1671 rr1i = dmpe(1) - (1.0d0-scalek*dmpi(1))*rr1 1672 rr3i = dmpe(3) - (1.0d0-scalek*dmpi(3))*rr3 1673 rr5i = dmpe(5) - (1.0d0-scalek*dmpi(5))*rr5 1674 rr1k = dmpe(1) - (1.0d0-scalek*dmpk(1))*rr1 1675 rr3k = dmpe(3) - (1.0d0-scalek*dmpk(3))*rr3 1676 rr5k = dmpe(5) - (1.0d0-scalek*dmpk(5))*rr5 1677 rr1ik = dmpe(1) - (1.0d0-scalek*dmpik(1))*rr1 1678 rr3ik = dmpe(3) - (1.0d0-scalek*dmpik(3))*rr3 1679 rr5ik = dmpe(5) - (1.0d0-scalek*dmpik(5))*rr5 1680 rr7ik = dmpe(7) - (1.0d0-scalek*dmpik(7))*rr7 1681 rr9ik = dmpe(9) - (1.0d0-scalek*dmpik(9))*rr9 1682 rr1 = dmpe(1) - (1.0d0-scalek)*rr1 1683 e = term1*rr1 + term4ik*rr7ik + term5ik*rr9ik 1684 & + term1i*rr1i + term1k*rr1k + term1ik*rr1ik 1685 & + term2i*rr3i + term2k*rr3k + term2ik*rr3ik 1686 & + term3i*rr5i + term3k*rr5k + term3ik*rr5ik 1687c 1688c find standard multipole intermediates and energy value 1689c 1690 else 1691 term1 = ci*ck 1692 term2 = ck*dir - ci*dkr + dik 1693 term3 = ci*qkr + ck*qir - dir*dkr 1694 & + 2.0d0*(dkqi-diqk+qiqk) 1695 term4 = dir*qkr - dkr*qir - 4.0d0*qik 1696 term5 = qir*qkr 1697 scalek = 1.0d0 - mscale(k) 1698 rr1 = dmpe(1) - scalek*rr1 1699 rr3 = dmpe(3) - scalek*rr3 1700 rr5 = dmpe(5) - scalek*rr5 1701 rr7 = dmpe(7) - scalek*rr7 1702 rr9 = dmpe(9) - scalek*rr9 1703 e = term1*rr1 + term2*rr3 + term3*rr5 1704 & + term4*rr7 + term5*rr9 1705 end if 1706c 1707c increment the overall multipole energy component 1708c 1709 em = em + e 1710 end if 1711 end do 1712c 1713c reset exclusion coefficients for connected atoms 1714c 1715 do j = 1, n12(i) 1716 mscale(i12(j,i)) = 1.0d0 1717 end do 1718 do j = 1, n13(i) 1719 mscale(i13(j,i)) = 1.0d0 1720 end do 1721 do j = 1, n14(i) 1722 mscale(i14(j,i)) = 1.0d0 1723 end do 1724 do j = 1, n15(i) 1725 mscale(i15(j,i)) = 1.0d0 1726 end do 1727 end do 1728c 1729c OpenMP directives for the major loop structure 1730c 1731!$OMP END DO 1732!$OMP END PARALLEL 1733c 1734c perform deallocation of some local arrays 1735c 1736 deallocate (mscale) 1737 return 1738 end 1739c 1740c 1741c ################################################################ 1742c ## ## 1743c ## subroutine emrecip -- PME recip space multipole energy ## 1744c ## ## 1745c ################################################################ 1746c 1747c 1748c "emrecip" evaluates the reciprocal space portion of the particle 1749c mesh Ewald energy due to atomic multipole interactions 1750c 1751c literature references: 1752c 1753c C. Sagui, L. G. Pedersen and T. A. Darden, "Towards an Accurate 1754c Representation of Electrostatics in Classical Force Fields: 1755c Efficient Implementation of Multipolar Interactions in 1756c Biomolecular Simulations", Journal of Chemical Physics, 120, 1757c 73-87 (2004) 1758c 1759c W. Smith and D. Fincham, "The Ewald Sum in Truncated Octahedral 1760c and Rhombic Dodecahedral Boundary Conditions", Molecular Physics, 1761c 10, 67-71 (1993) 1762c 1763c modifications for nonperiodic systems suggested by Tom Darden 1764c during May 2007 1765c 1766c 1767 subroutine emrecip 1768 use bound 1769 use boxes 1770 use chgpot 1771 use energi 1772 use ewald 1773 use math 1774 use mpole 1775 use mrecip 1776 use pme 1777 use potent 1778 implicit none 1779 integer i,j,k 1780 integer k1,k2,k3 1781 integer m1,m2,m3 1782 integer ntot,nff 1783 integer nf1,nf2,nf3 1784 real*8 e,r1,r2,r3 1785 real*8 f,h1,h2,h3 1786 real*8 volterm,denom 1787 real*8 hsq,expterm 1788 real*8 term,pterm 1789 real*8 struc2 1790c 1791c 1792c return if the Ewald coefficient is zero 1793c 1794 if (aewald .lt. 1.0d-6) return 1795 f = electric / dielec 1796c 1797c perform dynamic allocation of some global arrays 1798c 1799 if (allocated(cmp)) then 1800 if (size(cmp) .lt. 10*npole) deallocate (cmp) 1801 end if 1802 if (allocated(fmp)) then 1803 if (size(fmp) .lt. 10*npole) deallocate (fmp) 1804 end if 1805 if (allocated(fphi)) then 1806 if (size(fphi) .lt. 20*npole) deallocate (fphi) 1807 end if 1808 if (.not. allocated(cmp)) allocate (cmp(10,npole)) 1809 if (.not. allocated(fmp)) allocate (fmp(10,npole)) 1810 if (.not. allocated(fphi)) allocate (fphi(20,npole)) 1811c 1812c perform dynamic allocation of some global arrays 1813c 1814 ntot = nfft1 * nfft2 * nfft3 1815 if (allocated(qgrid)) then 1816 if (size(qgrid) .ne. 2*ntot) call fftclose 1817 end if 1818 if (allocated(qfac)) then 1819 if (size(qfac) .ne. ntot) deallocate (qfac) 1820 end if 1821 if (.not. allocated(qgrid)) call fftsetup 1822 if (.not. allocated(qfac)) allocate (qfac(nfft1,nfft2,nfft3)) 1823c 1824c setup spatial decomposition and B-spline coefficients 1825c 1826 call getchunk 1827 call moduli 1828 call bspline_fill 1829 call table_fill 1830c 1831c copy the multipole moments into local storage areas 1832c 1833 do i = 1, npole 1834 cmp(1,i) = rpole(1,i) 1835 cmp(2,i) = rpole(2,i) 1836 cmp(3,i) = rpole(3,i) 1837 cmp(4,i) = rpole(4,i) 1838 cmp(5,i) = rpole(5,i) 1839 cmp(6,i) = rpole(9,i) 1840 cmp(7,i) = rpole(13,i) 1841 cmp(8,i) = 2.0d0 * rpole(6,i) 1842 cmp(9,i) = 2.0d0 * rpole(7,i) 1843 cmp(10,i) = 2.0d0 * rpole(10,i) 1844 end do 1845c 1846c convert Cartesian multipoles to fractional coordinates 1847c 1848 call cmp_to_fmp (cmp,fmp) 1849c 1850c assign PME grid and perform 3-D FFT forward transform 1851c 1852 call grid_mpole (fmp) 1853 call fftfront 1854c 1855c make the scalar summation over reciprocal lattice 1856c 1857 qfac(1,1,1) = 0.0d0 1858 pterm = (pi/aewald)**2 1859 volterm = pi * volbox 1860 nf1 = (nfft1+1) / 2 1861 nf2 = (nfft2+1) / 2 1862 nf3 = (nfft3+1) / 2 1863 nff = nfft1 * nfft2 1864 ntot = nff * nfft3 1865 do i = 1, ntot-1 1866 k3 = i/nff + 1 1867 j = i - (k3-1)*nff 1868 k2 = j/nfft1 + 1 1869 k1 = j - (k2-1)*nfft1 + 1 1870 m1 = k1 - 1 1871 m2 = k2 - 1 1872 m3 = k3 - 1 1873 if (k1 .gt. nf1) m1 = m1 - nfft1 1874 if (k2 .gt. nf2) m2 = m2 - nfft2 1875 if (k3 .gt. nf3) m3 = m3 - nfft3 1876 r1 = dble(m1) 1877 r2 = dble(m2) 1878 r3 = dble(m3) 1879 h1 = recip(1,1)*r1 + recip(1,2)*r2 + recip(1,3)*r3 1880 h2 = recip(2,1)*r1 + recip(2,2)*r2 + recip(2,3)*r3 1881 h3 = recip(3,1)*r1 + recip(3,2)*r2 + recip(3,3)*r3 1882 hsq = h1*h1 + h2*h2 + h3*h3 1883 term = -pterm * hsq 1884 expterm = 0.0d0 1885 if (term .gt. -50.0d0) then 1886 denom = volterm*hsq*bsmod1(k1)*bsmod2(k2)*bsmod3(k3) 1887 expterm = exp(term) / denom 1888 if (.not. use_bounds) then 1889 expterm = expterm * (1.0d0-cos(pi*xbox*sqrt(hsq))) 1890 else if (nonprism) then 1891 if (mod(m1+m2+m3,2) .ne. 0) expterm = 0.0d0 1892 end if 1893 end if 1894 qfac(k1,k2,k3) = expterm 1895 end do 1896c 1897c account for zeroth grid point for nonperiodic system 1898c 1899 if (.not. use_bounds) then 1900 expterm = 0.5d0 * pi / xbox 1901 qfac(1,1,1) = expterm 1902 struc2 = qgrid(1,1,1,1)**2 + qgrid(2,1,1,1)**2 1903 e = f * expterm * struc2 1904 em = em + e 1905 end if 1906c 1907c complete the transformation of the PME grid 1908c 1909 do k = 1, nfft3 1910 do j = 1, nfft2 1911 do i = 1, nfft1 1912 term = qfac(i,j,k) 1913 qgrid(1,i,j,k) = term * qgrid(1,i,j,k) 1914 qgrid(2,i,j,k) = term * qgrid(2,i,j,k) 1915 end do 1916 end do 1917 end do 1918c 1919c perform 3-D FFT backward transform and get potential 1920c 1921 call fftback 1922 call fphi_mpole (fphi) 1923c 1924c sum over multipoles and increment total multipole energy 1925c 1926 e = 0.0d0 1927 do i = 1, npole 1928 do k = 1, 10 1929 e = e + fmp(k,i)*fphi(k,i) 1930 end do 1931 end do 1932 e = 0.5d0 * f * e 1933 em = em + e 1934 return 1935 end 1936