1c 2c 3c ############################################################# 4c ## COPYRIGHT (C) 1999 by Pengyu Ren & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ############################################################# 7c 8c ################################################################ 9c ## ## 10c ## subroutine empole1 -- mpole/polar energy & derivatives ## 11c ## ## 12c ################################################################ 13c 14c 15c "empole1" calculates the multipole and dipole polarization 16c energy and derivatives with respect to Cartesian coordinates 17c 18c 19 subroutine empole1 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 empole1d 29 else 30 call empole1c 31 end if 32 else 33 if (use_mlist) then 34 call empole1b 35 else 36 call empole1a 37 end if 38 end if 39 return 40 end 41c 42c 43c ################################################################## 44c ## ## 45c ## subroutine empole1a -- double loop multipole derivatives ## 46c ## ## 47c ################################################################## 48c 49c 50c "empole1a" calculates the multipole energy and derivatives with 51c respect to Cartesian coordinates using a pairwise double loop 52c 53c 54 subroutine empole1a 55 use sizes 56 use atoms 57 use bound 58 use cell 59 use chgpot 60 use couple 61 use deriv 62 use energi 63 use group 64 use inter 65 use molcul 66 use mplpot 67 use mpole 68 use shunt 69 use usage 70 use virial 71 implicit none 72 integer i,j,k 73 integer ii,kk,jcell 74 integer ix,iy,iz 75 integer kx,ky,kz 76 integer iax,iay,iaz 77 real*8 e,de,f,fgrp 78 real*8 xi,yi,zi 79 real*8 xr,yr,zr 80 real*8 xix,yix,zix 81 real*8 xiy,yiy,ziy 82 real*8 xiz,yiz,ziz 83 real*8 r,r2,rr1,rr3 84 real*8 rr5,rr7,rr9,rr11 85 real*8 ci,dix,diy,diz 86 real*8 qixx,qixy,qixz 87 real*8 qiyy,qiyz,qizz 88 real*8 ck,dkx,dky,dkz 89 real*8 qkxx,qkxy,qkxz 90 real*8 qkyy,qkyz,qkzz 91 real*8 dikx,diky,dikz 92 real*8 dirx,diry,dirz 93 real*8 dkrx,dkry,dkrz 94 real*8 qrix,qriy,qriz 95 real*8 qrkx,qrky,qrkz 96 real*8 qrixr,qriyr,qrizr 97 real*8 qrkxr,qrkyr,qrkzr 98 real*8 qrrx,qrry,qrrz 99 real*8 qikrx,qikry,qikrz 100 real*8 qkirx,qkiry,qkirz 101 real*8 qikrxr,qikryr,qikrzr 102 real*8 qkirxr,qkiryr,qkirzr 103 real*8 diqkx,diqky,diqkz 104 real*8 dkqix,dkqiy,dkqiz 105 real*8 diqkxr,diqkyr,diqkzr 106 real*8 dkqixr,dkqiyr,dkqizr 107 real*8 dqiqkx,dqiqky,dqiqkz 108 real*8 dri,drk,qrri,qrrk 109 real*8 diqrk,dkqri 110 real*8 dik,qik,qrrik 111 real*8 term1,term2,term3 112 real*8 term4,term5,term6 113 real*8 frcx,frcy,frcz 114 real*8 vxx,vyy,vzz 115 real*8 vxy,vxz,vyz 116 real*8 ttmi(3),ttmk(3) 117 real*8 fix(3),fiy(3),fiz(3) 118 real*8, allocatable :: mscale(:) 119 real*8, allocatable :: tem(:,:) 120 logical proceed,usei,usek 121 character*6 mode 122c 123c 124c zero out the atomic multipole energy and derivatives 125c 126 em = 0.0d0 127 do i = 1, n 128 do j = 1, 3 129 dem(j,i) = 0.0d0 130 end do 131 end do 132 if (npole .eq. 0) return 133c 134c check the sign of multipole components at chiral sites 135c 136 call chkpole 137c 138c rotate the multipole components into the global frame 139c 140 call rotpole 141c 142c perform dynamic allocation of some local arrays 143c 144 allocate (mscale(n)) 145 allocate (tem(3,n)) 146c 147c initialize connected atom scaling and torque arrays 148c 149 do i = 1, n 150 mscale(i) = 1.0d0 151 do j = 1, 3 152 tem(j,i) = 0.0d0 153 end do 154 end do 155c 156c set conversion factor, cutoff and switching coefficients 157c 158 f = electric / dielec 159 mode = 'MPOLE' 160 call switch (mode) 161c 162c compute the multipole interaction energy and gradient 163c 164 do i = 1, npole-1 165 ii = ipole(i) 166 iz = zaxis(i) 167 ix = xaxis(i) 168 iy = yaxis(i) 169 xi = x(ii) 170 yi = y(ii) 171 zi = z(ii) 172 ci = rpole(1,i) 173 dix = rpole(2,i) 174 diy = rpole(3,i) 175 diz = rpole(4,i) 176 qixx = rpole(5,i) 177 qixy = rpole(6,i) 178 qixz = rpole(7,i) 179 qiyy = rpole(9,i) 180 qiyz = rpole(10,i) 181 qizz = rpole(13,i) 182 usei = (use(ii) .or. use(iz) .or. use(ix) .or. use(iy)) 183 do j = 1, n12(ii) 184 mscale(i12(j,ii)) = m2scale 185 end do 186 do j = 1, n13(ii) 187 mscale(i13(j,ii)) = m3scale 188 end do 189 do j = 1, n14(ii) 190 mscale(i14(j,ii)) = m4scale 191 end do 192 do j = 1, n15(ii) 193 mscale(i15(j,ii)) = m5scale 194 end do 195c 196c evaluate all sites within the cutoff distance 197c 198 do k = i+1, npole 199 kk = ipole(k) 200 kz = zaxis(k) 201 kx = xaxis(k) 202 ky = yaxis(k) 203 usek = (use(kk) .or. use(kz) .or. use(kx) .or. use(ky)) 204 proceed = .true. 205 if (use_group) call groups (proceed,fgrp,ii,kk,0,0,0,0) 206 if (.not. use_intra) proceed = .true. 207 if (proceed) proceed = (usei .or. usek) 208 if (.not. proceed) goto 10 209 xr = x(kk) - xi 210 yr = y(kk) - yi 211 zr = z(kk) - zi 212 if (use_bounds) call image (xr,yr,zr) 213 r2 = xr*xr + yr*yr + zr*zr 214 if (r2 .le. off2) then 215 r = sqrt(r2) 216 ck = rpole(1,k) 217 dkx = rpole(2,k) 218 dky = rpole(3,k) 219 dkz = rpole(4,k) 220 qkxx = rpole(5,k) 221 qkxy = rpole(6,k) 222 qkxz = rpole(7,k) 223 qkyy = rpole(9,k) 224 qkyz = rpole(10,k) 225 qkzz = rpole(13,k) 226c 227c get reciprocal distance terms for this interaction 228c 229 rr1 = f * mscale(kk) / r 230 rr3 = rr1 / r2 231 rr5 = 3.0d0 * rr3 / r2 232 rr7 = 5.0d0 * rr5 / r2 233 rr9 = 7.0d0 * rr7 / r2 234 rr11 = 9.0d0 * rr9 / r2 235c 236c intermediates involving moments and distance separation 237c 238 dikx = diy*dkz - diz*dky 239 diky = diz*dkx - dix*dkz 240 dikz = dix*dky - diy*dkx 241 dirx = diy*zr - diz*yr 242 diry = diz*xr - dix*zr 243 dirz = dix*yr - diy*xr 244 dkrx = dky*zr - dkz*yr 245 dkry = dkz*xr - dkx*zr 246 dkrz = dkx*yr - dky*xr 247 dri = dix*xr + diy*yr + diz*zr 248 drk = dkx*xr + dky*yr + dkz*zr 249 dik = dix*dkx + diy*dky + diz*dkz 250 qrix = qixx*xr + qixy*yr + qixz*zr 251 qriy = qixy*xr + qiyy*yr + qiyz*zr 252 qriz = qixz*xr + qiyz*yr + qizz*zr 253 qrkx = qkxx*xr + qkxy*yr + qkxz*zr 254 qrky = qkxy*xr + qkyy*yr + qkyz*zr 255 qrkz = qkxz*xr + qkyz*yr + qkzz*zr 256 qrri = qrix*xr + qriy*yr + qriz*zr 257 qrrk = qrkx*xr + qrky*yr + qrkz*zr 258 qrrik = qrix*qrkx + qriy*qrky + qriz*qrkz 259 qik = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 260 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 261 qrixr = qriz*yr - qriy*zr 262 qriyr = qrix*zr - qriz*xr 263 qrizr = qriy*xr - qrix*yr 264 qrkxr = qrkz*yr - qrky*zr 265 qrkyr = qrkx*zr - qrkz*xr 266 qrkzr = qrky*xr - qrkx*yr 267 qrrx = qrky*qriz - qrkz*qriy 268 qrry = qrkz*qrix - qrkx*qriz 269 qrrz = qrkx*qriy - qrky*qrix 270 qikrx = qixx*qrkx + qixy*qrky + qixz*qrkz 271 qikry = qixy*qrkx + qiyy*qrky + qiyz*qrkz 272 qikrz = qixz*qrkx + qiyz*qrky + qizz*qrkz 273 qkirx = qkxx*qrix + qkxy*qriy + qkxz*qriz 274 qkiry = qkxy*qrix + qkyy*qriy + qkyz*qriz 275 qkirz = qkxz*qrix + qkyz*qriy + qkzz*qriz 276 qikrxr = qikrz*yr - qikry*zr 277 qikryr = qikrx*zr - qikrz*xr 278 qikrzr = qikry*xr - qikrx*yr 279 qkirxr = qkirz*yr - qkiry*zr 280 qkiryr = qkirx*zr - qkirz*xr 281 qkirzr = qkiry*xr - qkirx*yr 282 diqkx = dix*qkxx + diy*qkxy + diz*qkxz 283 diqky = dix*qkxy + diy*qkyy + diz*qkyz 284 diqkz = dix*qkxz + diy*qkyz + diz*qkzz 285 dkqix = dkx*qixx + dky*qixy + dkz*qixz 286 dkqiy = dkx*qixy + dky*qiyy + dkz*qiyz 287 dkqiz = dkx*qixz + dky*qiyz + dkz*qizz 288 diqrk = dix*qrkx + diy*qrky + diz*qrkz 289 dkqri = dkx*qrix + dky*qriy + dkz*qriz 290 diqkxr = diqkz*yr - diqky*zr 291 diqkyr = diqkx*zr - diqkz*xr 292 diqkzr = diqky*xr - diqkx*yr 293 dkqixr = dkqiz*yr - dkqiy*zr 294 dkqiyr = dkqix*zr - dkqiz*xr 295 dkqizr = dkqiy*xr - dkqix*yr 296 dqiqkx = diy*qrkz - diz*qrky + dky*qriz - dkz*qriy 297 & - 2.0d0*(qixy*qkxz+qiyy*qkyz+qiyz*qkzz 298 & -qixz*qkxy-qiyz*qkyy-qizz*qkyz) 299 dqiqky = diz*qrkx - dix*qrkz + dkz*qrix - dkx*qriz 300 & - 2.0d0*(qixz*qkxx+qiyz*qkxy+qizz*qkxz 301 & -qixx*qkxz-qixy*qkyz-qixz*qkzz) 302 dqiqkz = dix*qrky - diy*qrkx + dkx*qriy - dky*qrix 303 & - 2.0d0*(qixx*qkxy+qixy*qkyy+qixz*qkyz 304 & -qixy*qkxx-qiyy*qkxy-qiyz*qkxz) 305c 306c calculate intermediate terms for multipole energy 307c 308 term1 = ci*ck 309 term2 = ck*dri - ci*drk + dik 310 term3 = ci*qrrk + ck*qrri - dri*drk 311 & + 2.0d0*(dkqri-diqrk+qik) 312 term4 = dri*qrrk - drk*qrri - 4.0d0*qrrik 313 term5 = qrri*qrrk 314c 315c compute the energy contribution for this interaction 316c 317 e = term1*rr1 + term2*rr3 + term3*rr5 318 & + term4*rr7 + term5*rr9 319 em = em + e 320 if (molcule(ii) .ne. molcule(kk)) 321 & einter = einter + e 322c 323c calculate intermediate terms for force and torque 324c 325 de = term1*rr3 + term2*rr5 + term3*rr7 326 & + term4*rr9 + term5*rr11 327 term1 = -ck*rr3 + drk*rr5 - qrrk*rr7 328 term2 = ci*rr3 + dri*rr5 + qrri*rr7 329 term3 = 2.0d0 * rr5 330 term4 = 2.0d0 * (-ck*rr5+drk*rr7-qrrk*rr9) 331 term5 = 2.0d0 * (-ci*rr5-dri*rr7-qrri*rr9) 332 term6 = 4.0d0 * rr7 333c 334c compute the force components for this interaction 335c 336 frcx = de*xr + term1*dix + term2*dkx 337 & + term3*(diqkx-dkqix) + term4*qrix 338 & + term5*qrkx + term6*(qikrx+qkirx) 339 frcy = de*yr + term1*diy + term2*dky 340 & + term3*(diqky-dkqiy) + term4*qriy 341 & + term5*qrky + term6*(qikry+qkiry) 342 frcz = de*zr + term1*diz + term2*dkz 343 & + term3*(diqkz-dkqiz) + term4*qriz 344 & + term5*qrkz + term6*(qikrz+qkirz) 345c 346c compute the torque components for this interaction 347c 348 ttmi(1) = -rr3*dikx + term1*dirx + term3*(dqiqkx+dkqixr) 349 & - term4*qrixr - term6*(qikrxr+qrrx) 350 ttmi(2) = -rr3*diky + term1*diry + term3*(dqiqky+dkqiyr) 351 & - term4*qriyr - term6*(qikryr+qrry) 352 ttmi(3) = -rr3*dikz + term1*dirz + term3*(dqiqkz+dkqizr) 353 & - term4*qrizr - term6*(qikrzr+qrrz) 354 ttmk(1) = rr3*dikx + term2*dkrx - term3*(dqiqkx+diqkxr) 355 & - term5*qrkxr - term6*(qkirxr-qrrx) 356 ttmk(2) = rr3*diky + term2*dkry - term3*(dqiqky+diqkyr) 357 & - term5*qrkyr - term6*(qkiryr-qrry) 358 ttmk(3) = rr3*dikz + term2*dkrz - term3*(dqiqkz+diqkzr) 359 & - term5*qrkzr - term6*(qkirzr-qrrz) 360c 361c force and torque components scaled by group membership 362c 363 if (use_group) then 364 frcx = fgrp * frcx 365 frcy = fgrp * frcy 366 frcz = fgrp * frcz 367 do j = 1, 3 368 ttmi(j) = fgrp * ttmi(j) 369 ttmk(j) = fgrp * ttmk(j) 370 end do 371 end if 372c 373c increment force-based gradient and torque on first site 374c 375 dem(1,ii) = dem(1,ii) + frcx 376 dem(2,ii) = dem(2,ii) + frcy 377 dem(3,ii) = dem(3,ii) + frcz 378 tem(1,i) = tem(1,i) + ttmi(1) 379 tem(2,i) = tem(2,i) + ttmi(2) 380 tem(3,i) = tem(3,i) + ttmi(3) 381c 382c increment force-based gradient and torque on second site 383c 384 dem(1,kk) = dem(1,kk) - frcx 385 dem(2,kk) = dem(2,kk) - frcy 386 dem(3,kk) = dem(3,kk) - frcz 387 tem(1,k) = tem(1,k) + ttmk(1) 388 tem(2,k) = tem(2,k) + ttmk(2) 389 tem(3,k) = tem(3,k) + ttmk(3) 390c 391c increment the virial due to pairwise Cartesian forces 392c 393 vxx = -xr * frcx 394 vxy = -yr * frcx 395 vxz = -zr * frcx 396 vyy = -yr * frcy 397 vyz = -zr * frcy 398 vzz = -zr * frcz 399 vir(1,1) = vir(1,1) + vxx 400 vir(2,1) = vir(2,1) + vxy 401 vir(3,1) = vir(3,1) + vxz 402 vir(1,2) = vir(1,2) + vxy 403 vir(2,2) = vir(2,2) + vyy 404 vir(3,2) = vir(3,2) + vyz 405 vir(1,3) = vir(1,3) + vxz 406 vir(2,3) = vir(2,3) + vyz 407 vir(3,3) = vir(3,3) + vzz 408 end if 409 10 continue 410 end do 411c 412c reset exclusion coefficients for connected atoms 413c 414 do j = 1, n12(ii) 415 mscale(i12(j,ii)) = 1.0d0 416 end do 417 do j = 1, n13(ii) 418 mscale(i13(j,ii)) = 1.0d0 419 end do 420 do j = 1, n14(ii) 421 mscale(i14(j,ii)) = 1.0d0 422 end do 423 do j = 1, n15(ii) 424 mscale(i15(j,ii)) = 1.0d0 425 end do 426 end do 427c 428c for periodic boundary conditions with large cutoffs 429c neighbors must be found by the replicates method 430c 431 if (use_replica) then 432c 433c calculate interaction with other unit cells 434c 435 do i = 1, npole 436 ii = ipole(i) 437 iz = zaxis(i) 438 ix = xaxis(i) 439 iy = yaxis(i) 440 xi = x(ii) 441 yi = y(ii) 442 zi = z(ii) 443 ci = rpole(1,i) 444 dix = rpole(2,i) 445 diy = rpole(3,i) 446 diz = rpole(4,i) 447 qixx = rpole(5,i) 448 qixy = rpole(6,i) 449 qixz = rpole(7,i) 450 qiyy = rpole(9,i) 451 qiyz = rpole(10,i) 452 qizz = rpole(13,i) 453 usei = (use(ii) .or. use(iz) .or. use(ix) .or. use(iy)) 454 do j = 1, n12(ii) 455 mscale(i12(j,ii)) = m2scale 456 end do 457 do j = 1, n13(ii) 458 mscale(i13(j,ii)) = m3scale 459 end do 460 do j = 1, n14(ii) 461 mscale(i14(j,ii)) = m4scale 462 end do 463 do j = 1, n15(ii) 464 mscale(i15(j,ii)) = m5scale 465 end do 466c 467c evaluate all sites within the cutoff distance 468c 469 do k = i, npole 470 kk = ipole(k) 471 kz = zaxis(k) 472 kx = xaxis(k) 473 ky = yaxis(k) 474 usek = (use(kk) .or. use(kz) .or. use(kx) .or. use(ky)) 475 if (use_group) call groups (proceed,fgrp,ii,kk,0,0,0,0) 476 proceed = .true. 477 if (proceed) proceed = (usei .or. usek) 478 if (.not. proceed) goto 20 479 do jcell = 2, ncell 480 xr = x(kk) - xi 481 yr = y(kk) - yi 482 zr = z(kk) - zi 483 call imager (xr,yr,zr,jcell) 484 r2 = xr*xr + yr*yr + zr*zr 485 if (.not. (use_polymer .and. r2.le.polycut2)) then 486 mscale(kk) = 1.0d0 487 end if 488 if (r2 .le. off2) then 489 r = sqrt(r2) 490 ck = rpole(1,k) 491 dkx = rpole(2,k) 492 dky = rpole(3,k) 493 dkz = rpole(4,k) 494 qkxx = rpole(5,k) 495 qkxy = rpole(6,k) 496 qkxz = rpole(7,k) 497 qkyy = rpole(9,k) 498 qkyz = rpole(10,k) 499 qkzz = rpole(13,k) 500c 501c get reciprocal distance terms for this interaction 502c 503 rr1 = f * mscale(kk) / r 504 rr3 = rr1 / r2 505 rr5 = 3.0d0 * rr3 / r2 506 rr7 = 5.0d0 * rr5 / r2 507 rr9 = 7.0d0 * rr7 / r2 508 rr11 = 9.0d0 * rr9 / r2 509c 510c intermediates involving moments and distance separation 511c 512 dikx = diy*dkz - diz*dky 513 diky = diz*dkx - dix*dkz 514 dikz = dix*dky - diy*dkx 515 dirx = diy*zr - diz*yr 516 diry = diz*xr - dix*zr 517 dirz = dix*yr - diy*xr 518 dkrx = dky*zr - dkz*yr 519 dkry = dkz*xr - dkx*zr 520 dkrz = dkx*yr - dky*xr 521 dri = dix*xr + diy*yr + diz*zr 522 drk = dkx*xr + dky*yr + dkz*zr 523 dik = dix*dkx + diy*dky + diz*dkz 524 qrix = qixx*xr + qixy*yr + qixz*zr 525 qriy = qixy*xr + qiyy*yr + qiyz*zr 526 qriz = qixz*xr + qiyz*yr + qizz*zr 527 qrkx = qkxx*xr + qkxy*yr + qkxz*zr 528 qrky = qkxy*xr + qkyy*yr + qkyz*zr 529 qrkz = qkxz*xr + qkyz*yr + qkzz*zr 530 qrri = qrix*xr + qriy*yr + qriz*zr 531 qrrk = qrkx*xr + qrky*yr + qrkz*zr 532 qrrik = qrix*qrkx + qriy*qrky + qriz*qrkz 533 qik = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 534 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 535 qrixr = qriz*yr - qriy*zr 536 qriyr = qrix*zr - qriz*xr 537 qrizr = qriy*xr - qrix*yr 538 qrkxr = qrkz*yr - qrky*zr 539 qrkyr = qrkx*zr - qrkz*xr 540 qrkzr = qrky*xr - qrkx*yr 541 qrrx = qrky*qriz - qrkz*qriy 542 qrry = qrkz*qrix - qrkx*qriz 543 qrrz = qrkx*qriy - qrky*qrix 544 qikrx = qixx*qrkx + qixy*qrky + qixz*qrkz 545 qikry = qixy*qrkx + qiyy*qrky + qiyz*qrkz 546 qikrz = qixz*qrkx + qiyz*qrky + qizz*qrkz 547 qkirx = qkxx*qrix + qkxy*qriy + qkxz*qriz 548 qkiry = qkxy*qrix + qkyy*qriy + qkyz*qriz 549 qkirz = qkxz*qrix + qkyz*qriy + qkzz*qriz 550 qikrxr = qikrz*yr - qikry*zr 551 qikryr = qikrx*zr - qikrz*xr 552 qikrzr = qikry*xr - qikrx*yr 553 qkirxr = qkirz*yr - qkiry*zr 554 qkiryr = qkirx*zr - qkirz*xr 555 qkirzr = qkiry*xr - qkirx*yr 556 diqkx = dix*qkxx + diy*qkxy + diz*qkxz 557 diqky = dix*qkxy + diy*qkyy + diz*qkyz 558 diqkz = dix*qkxz + diy*qkyz + diz*qkzz 559 dkqix = dkx*qixx + dky*qixy + dkz*qixz 560 dkqiy = dkx*qixy + dky*qiyy + dkz*qiyz 561 dkqiz = dkx*qixz + dky*qiyz + dkz*qizz 562 diqrk = dix*qrkx + diy*qrky + diz*qrkz 563 dkqri = dkx*qrix + dky*qriy + dkz*qriz 564 diqkxr = diqkz*yr - diqky*zr 565 diqkyr = diqkx*zr - diqkz*xr 566 diqkzr = diqky*xr - diqkx*yr 567 dkqixr = dkqiz*yr - dkqiy*zr 568 dkqiyr = dkqix*zr - dkqiz*xr 569 dkqizr = dkqiy*xr - dkqix*yr 570 dqiqkx = diy*qrkz - diz*qrky + dky*qriz - dkz*qriy 571 & - 2.0d0*(qixy*qkxz+qiyy*qkyz+qiyz*qkzz 572 & -qixz*qkxy-qiyz*qkyy-qizz*qkyz) 573 dqiqky = diz*qrkx - dix*qrkz + dkz*qrix - dkx*qriz 574 & - 2.0d0*(qixz*qkxx+qiyz*qkxy+qizz*qkxz 575 & -qixx*qkxz-qixy*qkyz-qixz*qkzz) 576 dqiqkz = dix*qrky - diy*qrkx + dkx*qriy - dky*qrix 577 & - 2.0d0*(qixx*qkxy+qixy*qkyy+qixz*qkyz 578 & -qixy*qkxx-qiyy*qkxy-qiyz*qkxz) 579c 580c calculate intermediate terms for multipole energy 581c 582 term1 = ci*ck 583 term2 = ck*dri - ci*drk + dik 584 term3 = ci*qrrk + ck*qrri - dri*drk 585 & + 2.0d0*(dkqri-diqrk+qik) 586 term4 = dri*qrrk - drk*qrri - 4.0d0*qrrik 587 term5 = qrri*qrrk 588c 589c compute the energy contribution for this interaction 590c 591 e = term1*rr1 + term2*rr3 + term3*rr5 592 & + term4*rr7 + term5*rr9 593 if (ii .eq. kk) e = 0.5d0 * e 594 em = em + e 595 einter = einter + e 596c 597c calculate intermediate terms for force and torque 598c 599 de = term1*rr3 + term2*rr5 + term3*rr7 600 & + term4*rr9 + term5*rr11 601 term1 = -ck*rr3 + drk*rr5 - qrrk*rr7 602 term2 = ci*rr3 + dri*rr5 + qrri*rr7 603 term3 = 2.0d0 * rr5 604 term4 = 2.0d0 * (-ck*rr5+drk*rr7-qrrk*rr9) 605 term5 = 2.0d0 * (-ci*rr5-dri*rr7-qrri*rr9) 606 term6 = 4.0d0 * rr7 607c 608c compute the force components for this interaction 609c 610 frcx = de*xr + term1*dix + term2*dkx 611 & + term3*(diqkx-dkqix) + term4*qrix 612 & + term5*qrkx + term6*(qikrx+qkirx) 613 frcy = de*yr + term1*diy + term2*dky 614 & + term3*(diqky-dkqiy) + term4*qriy 615 & + term5*qrky + term6*(qikry+qkiry) 616 frcz = de*zr + term1*diz + term2*dkz 617 & + term3*(diqkz-dkqiz) + term4*qriz 618 & + term5*qrkz + term6*(qikrz+qkirz) 619c 620c compute the torque components for this interaction 621c 622 ttmi(1) = -rr3*dikx + term1*dirx + term3*(dqiqkx+dkqixr) 623 & - term4*qrixr - term6*(qikrxr+qrrx) 624 ttmi(2) = -rr3*diky + term1*diry + term3*(dqiqky+dkqiyr) 625 & - term4*qriyr - term6*(qikryr+qrry) 626 ttmi(3) = -rr3*dikz + term1*dirz + term3*(dqiqkz+dkqizr) 627 & - term4*qrizr - term6*(qikrzr+qrrz) 628 ttmk(1) = rr3*dikx + term2*dkrx - term3*(dqiqkx+diqkxr) 629 & - term5*qrkxr - term6*(qkirxr-qrrx) 630 ttmk(2) = rr3*diky + term2*dkry - term3*(dqiqky+diqkyr) 631 & - term5*qrkyr - term6*(qkiryr-qrry) 632 ttmk(3) = rr3*dikz + term2*dkrz - term3*(dqiqkz+diqkzr) 633 & - term5*qrkzr - term6*(qkirzr-qrrz) 634c 635c force and torque scaled for self-interactions and groups 636c 637 if (ii .eq. kk) then 638 frcx = 0.5d0 * frcx 639 frcy = 0.5d0 * frcy 640 frcz = 0.5d0 * frcz 641 do j = 1, 3 642 ttmi(j) = 0.5d0 * ttmi(j) 643 ttmk(j) = 0.5d0 * ttmk(j) 644 end do 645 end if 646 if (use_group) then 647 frcx = fgrp * frcx 648 frcy = fgrp * frcy 649 frcz = fgrp * frcz 650 do j = 1, 3 651 ttmi(j) = fgrp * ttmi(j) 652 ttmk(j) = fgrp * ttmk(j) 653 end do 654 end if 655c 656c increment force-based gradient and torque on first site 657c 658 dem(1,ii) = dem(1,ii) + frcx 659 dem(2,ii) = dem(2,ii) + frcy 660 dem(3,ii) = dem(3,ii) + frcz 661 tem(1,i) = tem(1,i) + ttmi(1) 662 tem(2,i) = tem(2,i) + ttmi(2) 663 tem(3,i) = tem(3,i) + ttmi(3) 664c 665c increment force-based gradient and torque on second site 666c 667 dem(1,kk) = dem(1,kk) - frcx 668 dem(2,kk) = dem(2,kk) - frcy 669 dem(3,kk) = dem(3,kk) - frcz 670 tem(1,k) = tem(1,k) + ttmk(1) 671 tem(2,k) = tem(2,k) + ttmk(2) 672 tem(3,k) = tem(3,k) + ttmk(3) 673c 674c increment the virial due to pairwise Cartesian forces 675c 676 vxx = -xr * frcx 677 vxy = -yr * frcx 678 vxz = -zr * frcx 679 vyy = -yr * frcy 680 vyz = -zr * frcy 681 vzz = -zr * frcz 682 vir(1,1) = vir(1,1) + vxx 683 vir(2,1) = vir(2,1) + vxy 684 vir(3,1) = vir(3,1) + vxz 685 vir(1,2) = vir(1,2) + vxy 686 vir(2,2) = vir(2,2) + vyy 687 vir(3,2) = vir(3,2) + vyz 688 vir(1,3) = vir(1,3) + vxz 689 vir(2,3) = vir(2,3) + vyz 690 vir(3,3) = vir(3,3) + vzz 691 end if 692 end do 693 20 continue 694 end do 695c 696c reset exclusion coefficients for connected atoms 697c 698 do j = 1, n12(ii) 699 mscale(i12(j,ii)) = 1.0d0 700 end do 701 do j = 1, n13(ii) 702 mscale(i13(j,ii)) = 1.0d0 703 end do 704 do j = 1, n14(ii) 705 mscale(i14(j,ii)) = 1.0d0 706 end do 707 do j = 1, n15(ii) 708 mscale(i15(j,ii)) = 1.0d0 709 end do 710 end do 711 end if 712c 713c resolve site torques then increment forces and virial 714c 715 do i = 1, npole 716 call torque (i,tem(1,i),fix,fiy,fiz,dem) 717 ii = ipole(i) 718 iaz = zaxis(i) 719 iax = xaxis(i) 720 iay = yaxis(i) 721 if (iaz .eq. 0) iaz = ii 722 if (iax .eq. 0) iax = ii 723 if (iay .eq. 0) iay = ii 724 xiz = x(iaz) - x(ii) 725 yiz = y(iaz) - y(ii) 726 ziz = z(iaz) - z(ii) 727 xix = x(iax) - x(ii) 728 yix = y(iax) - y(ii) 729 zix = z(iax) - z(ii) 730 xiy = x(iay) - x(ii) 731 yiy = y(iay) - y(ii) 732 ziy = z(iay) - z(ii) 733 vxx = xix*fix(1) + xiy*fiy(1) + xiz*fiz(1) 734 vxy = 0.5d0 * (yix*fix(1) + yiy*fiy(1) + yiz*fiz(1) 735 & + xix*fix(2) + xiy*fiy(2) + xiz*fiz(2)) 736 vxz = 0.5d0 * (zix*fix(1) + ziy*fiy(1) + ziz*fiz(1) 737 & + xix*fix(3) + xiy*fiy(3) + xiz*fiz(3)) 738 vyy = yix*fix(2) + yiy*fiy(2) + yiz*fiz(2) 739 vyz = 0.5d0 * (zix*fix(2) + ziy*fiy(2) + ziz*fiz(2) 740 & + yix*fix(3) + yiy*fiy(3) + yiz*fiz(3)) 741 vzz = zix*fix(3) + ziy*fiy(3) + ziz*fiz(3) 742 vir(1,1) = vir(1,1) + vxx 743 vir(2,1) = vir(2,1) + vxy 744 vir(3,1) = vir(3,1) + vxz 745 vir(1,2) = vir(1,2) + vxy 746 vir(2,2) = vir(2,2) + vyy 747 vir(3,2) = vir(3,2) + vyz 748 vir(1,3) = vir(1,3) + vxz 749 vir(2,3) = vir(2,3) + vyz 750 vir(3,3) = vir(3,3) + vzz 751 end do 752c 753c perform deallocation of some local arrays 754c 755 deallocate (mscale) 756 deallocate (tem) 757 return 758 end 759c 760c 761c ############################################################### 762c ## ## 763c ## subroutine empole1b -- neighbor list multipole derivs ## 764c ## ## 765c ############################################################### 766c 767c 768c "empole1b" calculates the multipole energy and derivatives 769c with respect to Cartesian coordinates using a neighbor list 770c 771c 772 subroutine empole1b 773 use sizes 774 use atoms 775 use bound 776 use chgpot 777 use couple 778 use deriv 779 use energi 780 use group 781 use inter 782 use molcul 783 use mplpot 784 use mpole 785 use neigh 786 use shunt 787 use usage 788 use virial 789 implicit none 790 integer i,j,k 791 integer ii,kk,kkk 792 integer ix,iy,iz 793 integer kx,ky,kz 794 integer iax,iay,iaz 795 real*8 e,de,f,fgrp 796 real*8 xi,yi,zi 797 real*8 xr,yr,zr 798 real*8 xix,yix,zix 799 real*8 xiy,yiy,ziy 800 real*8 xiz,yiz,ziz 801 real*8 r,r2,rr1,rr3 802 real*8 rr5,rr7,rr9,rr11 803 real*8 ci,dix,diy,diz 804 real*8 qixx,qixy,qixz 805 real*8 qiyy,qiyz,qizz 806 real*8 ck,dkx,dky,dkz 807 real*8 qkxx,qkxy,qkxz 808 real*8 qkyy,qkyz,qkzz 809 real*8 dikx,diky,dikz 810 real*8 dirx,diry,dirz 811 real*8 dkrx,dkry,dkrz 812 real*8 qrix,qriy,qriz 813 real*8 qrkx,qrky,qrkz 814 real*8 qrixr,qriyr,qrizr 815 real*8 qrkxr,qrkyr,qrkzr 816 real*8 qrrx,qrry,qrrz 817 real*8 qikrx,qikry,qikrz 818 real*8 qkirx,qkiry,qkirz 819 real*8 qikrxr,qikryr,qikrzr 820 real*8 qkirxr,qkiryr,qkirzr 821 real*8 diqkx,diqky,diqkz 822 real*8 dkqix,dkqiy,dkqiz 823 real*8 diqkxr,diqkyr,diqkzr 824 real*8 dkqixr,dkqiyr,dkqizr 825 real*8 dqiqkx,dqiqky,dqiqkz 826 real*8 dri,drk,qrri,qrrk 827 real*8 diqrk,dkqri 828 real*8 dik,qik,qrrik 829 real*8 term1,term2,term3 830 real*8 term4,term5,term6 831 real*8 frcx,frcy,frcz 832 real*8 vxx,vyy,vzz 833 real*8 vxy,vxz,vyz 834 real*8 ttmi(3),ttmk(3) 835 real*8 fix(3),fiy(3),fiz(3) 836 real*8, allocatable :: mscale(:) 837 real*8, allocatable :: tem(:,:) 838 logical proceed,usei,usek 839 character*6 mode 840c 841c 842c zero out the atomic multipole energy and derivatives 843c 844 em = 0.0d0 845 do i = 1, n 846 do j = 1, 3 847 dem(j,i) = 0.0d0 848 end do 849 end do 850 if (npole .eq. 0) return 851c 852c check the sign of multipole components at chiral sites 853c 854 call chkpole 855c 856c rotate the multipole components into the global frame 857c 858 call rotpole 859c 860c perform dynamic allocation of some local arrays 861c 862 allocate (mscale(n)) 863 allocate (tem(3,n)) 864c 865c initialize connected atom scaling and torque arrays 866c 867 do i = 1, n 868 mscale(i) = 1.0d0 869 do j = 1, 3 870 tem(j,i) = 0.0d0 871 end do 872 end do 873c 874c set conversion factor, cutoff and scaling coefficients 875c 876 f = electric / dielec 877 mode = 'MPOLE' 878 call switch (mode) 879c 880c OpenMP directives for the major loop structure 881c 882!$OMP PARALLEL default(private) 883!$OMP& shared(npole,ipole,x,y,z,xaxis,yaxis,zaxis,rpole,use,n12, 884!$OMP& i12,n13,i13,n14,i14,n15,i15,m2scale,m3scale,m4scale,m5scale, 885!$OMP& nelst,elst,use_group,use_intra,use_bounds,off2,f,molcule) 886!$OMP& firstprivate(mscale) shared (em,einter,dem,tem,vir) 887!$OMP DO reduction(+:em,einter,dem,tem,vir) schedule(guided) 888c 889c compute the multipole interaction energy and gradient 890c 891 do i = 1, npole 892 ii = ipole(i) 893 iz = zaxis(i) 894 ix = xaxis(i) 895 iy = yaxis(i) 896 xi = x(ii) 897 yi = y(ii) 898 zi = z(ii) 899 ci = rpole(1,i) 900 dix = rpole(2,i) 901 diy = rpole(3,i) 902 diz = rpole(4,i) 903 qixx = rpole(5,i) 904 qixy = rpole(6,i) 905 qixz = rpole(7,i) 906 qiyy = rpole(9,i) 907 qiyz = rpole(10,i) 908 qizz = rpole(13,i) 909 usei = (use(ii) .or. use(iz) .or. use(ix) .or. use(iy)) 910 do j = 1, n12(ii) 911 mscale(i12(j,ii)) = m2scale 912 end do 913 do j = 1, n13(ii) 914 mscale(i13(j,ii)) = m3scale 915 end do 916 do j = 1, n14(ii) 917 mscale(i14(j,ii)) = m4scale 918 end do 919 do j = 1, n15(ii) 920 mscale(i15(j,ii)) = m5scale 921 end do 922c 923c evaluate all sites within the cutoff distance 924c 925 do kkk = 1, nelst(i) 926 k = elst(kkk,i) 927 kk = ipole(k) 928 kz = zaxis(k) 929 kx = xaxis(k) 930 ky = yaxis(k) 931 usek = (use(kk) .or. use(kz) .or. use(kx) .or. use(ky)) 932 proceed = .true. 933 if (use_group) call groups (proceed,fgrp,ii,kk,0,0,0,0) 934 if (.not. use_intra) proceed = .true. 935 if (proceed) proceed = (usei .or. usek) 936 if (.not. proceed) goto 10 937 xr = x(kk) - xi 938 yr = y(kk) - yi 939 zr = z(kk) - zi 940 if (use_bounds) call image (xr,yr,zr) 941 r2 = xr*xr + yr*yr + zr*zr 942 if (r2 .le. off2) then 943 r = sqrt(r2) 944 ck = rpole(1,k) 945 dkx = rpole(2,k) 946 dky = rpole(3,k) 947 dkz = rpole(4,k) 948 qkxx = rpole(5,k) 949 qkxy = rpole(6,k) 950 qkxz = rpole(7,k) 951 qkyy = rpole(9,k) 952 qkyz = rpole(10,k) 953 qkzz = rpole(13,k) 954c 955c get reciprocal distance terms for this interaction 956c 957 rr1 = f * mscale(kk) / r 958 rr3 = rr1 / r2 959 rr5 = 3.0d0 * rr3 / r2 960 rr7 = 5.0d0 * rr5 / r2 961 rr9 = 7.0d0 * rr7 / r2 962 rr11 = 9.0d0 * rr9 / r2 963c 964c intermediates involving moments and distance separation 965c 966 dikx = diy*dkz - diz*dky 967 diky = diz*dkx - dix*dkz 968 dikz = dix*dky - diy*dkx 969 dirx = diy*zr - diz*yr 970 diry = diz*xr - dix*zr 971 dirz = dix*yr - diy*xr 972 dkrx = dky*zr - dkz*yr 973 dkry = dkz*xr - dkx*zr 974 dkrz = dkx*yr - dky*xr 975 dri = dix*xr + diy*yr + diz*zr 976 drk = dkx*xr + dky*yr + dkz*zr 977 dik = dix*dkx + diy*dky + diz*dkz 978 qrix = qixx*xr + qixy*yr + qixz*zr 979 qriy = qixy*xr + qiyy*yr + qiyz*zr 980 qriz = qixz*xr + qiyz*yr + qizz*zr 981 qrkx = qkxx*xr + qkxy*yr + qkxz*zr 982 qrky = qkxy*xr + qkyy*yr + qkyz*zr 983 qrkz = qkxz*xr + qkyz*yr + qkzz*zr 984 qrri = qrix*xr + qriy*yr + qriz*zr 985 qrrk = qrkx*xr + qrky*yr + qrkz*zr 986 qrrik = qrix*qrkx + qriy*qrky + qriz*qrkz 987 qik = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 988 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 989 qrixr = qriz*yr - qriy*zr 990 qriyr = qrix*zr - qriz*xr 991 qrizr = qriy*xr - qrix*yr 992 qrkxr = qrkz*yr - qrky*zr 993 qrkyr = qrkx*zr - qrkz*xr 994 qrkzr = qrky*xr - qrkx*yr 995 qrrx = qrky*qriz - qrkz*qriy 996 qrry = qrkz*qrix - qrkx*qriz 997 qrrz = qrkx*qriy - qrky*qrix 998 qikrx = qixx*qrkx + qixy*qrky + qixz*qrkz 999 qikry = qixy*qrkx + qiyy*qrky + qiyz*qrkz 1000 qikrz = qixz*qrkx + qiyz*qrky + qizz*qrkz 1001 qkirx = qkxx*qrix + qkxy*qriy + qkxz*qriz 1002 qkiry = qkxy*qrix + qkyy*qriy + qkyz*qriz 1003 qkirz = qkxz*qrix + qkyz*qriy + qkzz*qriz 1004 qikrxr = qikrz*yr - qikry*zr 1005 qikryr = qikrx*zr - qikrz*xr 1006 qikrzr = qikry*xr - qikrx*yr 1007 qkirxr = qkirz*yr - qkiry*zr 1008 qkiryr = qkirx*zr - qkirz*xr 1009 qkirzr = qkiry*xr - qkirx*yr 1010 diqkx = dix*qkxx + diy*qkxy + diz*qkxz 1011 diqky = dix*qkxy + diy*qkyy + diz*qkyz 1012 diqkz = dix*qkxz + diy*qkyz + diz*qkzz 1013 dkqix = dkx*qixx + dky*qixy + dkz*qixz 1014 dkqiy = dkx*qixy + dky*qiyy + dkz*qiyz 1015 dkqiz = dkx*qixz + dky*qiyz + dkz*qizz 1016 diqrk = dix*qrkx + diy*qrky + diz*qrkz 1017 dkqri = dkx*qrix + dky*qriy + dkz*qriz 1018 diqkxr = diqkz*yr - diqky*zr 1019 diqkyr = diqkx*zr - diqkz*xr 1020 diqkzr = diqky*xr - diqkx*yr 1021 dkqixr = dkqiz*yr - dkqiy*zr 1022 dkqiyr = dkqix*zr - dkqiz*xr 1023 dkqizr = dkqiy*xr - dkqix*yr 1024 dqiqkx = diy*qrkz - diz*qrky + dky*qriz - dkz*qriy 1025 & - 2.0d0*(qixy*qkxz+qiyy*qkyz+qiyz*qkzz 1026 & -qixz*qkxy-qiyz*qkyy-qizz*qkyz) 1027 dqiqky = diz*qrkx - dix*qrkz + dkz*qrix - dkx*qriz 1028 & - 2.0d0*(qixz*qkxx+qiyz*qkxy+qizz*qkxz 1029 & -qixx*qkxz-qixy*qkyz-qixz*qkzz) 1030 dqiqkz = dix*qrky - diy*qrkx + dkx*qriy - dky*qrix 1031 & - 2.0d0*(qixx*qkxy+qixy*qkyy+qixz*qkyz 1032 & -qixy*qkxx-qiyy*qkxy-qiyz*qkxz) 1033c 1034c calculate intermediate terms for multipole energy 1035c 1036 term1 = ci*ck 1037 term2 = ck*dri - ci*drk + dik 1038 term3 = ci*qrrk + ck*qrri - dri*drk 1039 & + 2.0d0*(dkqri-diqrk+qik) 1040 term4 = dri*qrrk - drk*qrri - 4.0d0*qrrik 1041 term5 = qrri*qrrk 1042c 1043c compute the energy contribution for this interaction 1044c 1045 e = term1*rr1 + term2*rr3 + term3*rr5 1046 & + term4*rr7 + term5*rr9 1047 if (use_group) e = e * fgrp 1048 em = em + e 1049 if (molcule(ii) .ne. molcule(kk)) 1050 & einter = einter + e 1051c 1052c calculate intermediate terms for force and torque 1053c 1054 de = term1*rr3 + term2*rr5 + term3*rr7 1055 & + term4*rr9 + term5*rr11 1056 term1 = -ck*rr3 + drk*rr5 - qrrk*rr7 1057 term2 = ci*rr3 + dri*rr5 + qrri*rr7 1058 term3 = 2.0d0 * rr5 1059 term4 = 2.0d0 * (-ck*rr5+drk*rr7-qrrk*rr9) 1060 term5 = 2.0d0 * (-ci*rr5-dri*rr7-qrri*rr9) 1061 term6 = 4.0d0 * rr7 1062c 1063c compute the force components for this interaction 1064c 1065 frcx = de*xr + term1*dix + term2*dkx 1066 & + term3*(diqkx-dkqix) + term4*qrix 1067 & + term5*qrkx + term6*(qikrx+qkirx) 1068 frcy = de*yr + term1*diy + term2*dky 1069 & + term3*(diqky-dkqiy) + term4*qriy 1070 & + term5*qrky + term6*(qikry+qkiry) 1071 frcz = de*zr + term1*diz + term2*dkz 1072 & + term3*(diqkz-dkqiz) + term4*qriz 1073 & + term5*qrkz + term6*(qikrz+qkirz) 1074c 1075c compute the torque components for this interaction 1076c 1077 ttmi(1) = -rr3*dikx + term1*dirx + term3*(dqiqkx+dkqixr) 1078 & - term4*qrixr - term6*(qikrxr+qrrx) 1079 ttmi(2) = -rr3*diky + term1*diry + term3*(dqiqky+dkqiyr) 1080 & - term4*qriyr - term6*(qikryr+qrry) 1081 ttmi(3) = -rr3*dikz + term1*dirz + term3*(dqiqkz+dkqizr) 1082 & - term4*qrizr - term6*(qikrzr+qrrz) 1083 ttmk(1) = rr3*dikx + term2*dkrx - term3*(dqiqkx+diqkxr) 1084 & - term5*qrkxr - term6*(qkirxr-qrrx) 1085 ttmk(2) = rr3*diky + term2*dkry - term3*(dqiqky+diqkyr) 1086 & - term5*qrkyr - term6*(qkiryr-qrry) 1087 ttmk(3) = rr3*dikz + term2*dkrz - term3*(dqiqkz+diqkzr) 1088 & - term5*qrkzr - term6*(qkirzr-qrrz) 1089c 1090c increment force-based gradient and torque on first site 1091c 1092 dem(1,ii) = dem(1,ii) + frcx 1093 dem(2,ii) = dem(2,ii) + frcy 1094 dem(3,ii) = dem(3,ii) + frcz 1095 tem(1,i) = tem(1,i) + ttmi(1) 1096 tem(2,i) = tem(2,i) + ttmi(2) 1097 tem(3,i) = tem(3,i) + ttmi(3) 1098c 1099c increment force-based gradient and torque on second site 1100c 1101 dem(1,kk) = dem(1,kk) - frcx 1102 dem(2,kk) = dem(2,kk) - frcy 1103 dem(3,kk) = dem(3,kk) - frcz 1104 tem(1,k) = tem(1,k) + ttmk(1) 1105 tem(2,k) = tem(2,k) + ttmk(2) 1106 tem(3,k) = tem(3,k) + ttmk(3) 1107c 1108c increment the virial due to pairwise Cartesian forces 1109c 1110 vxx = -xr * frcx 1111 vxy = -yr * frcx 1112 vxz = -zr * frcx 1113 vyy = -yr * frcy 1114 vyz = -zr * frcy 1115 vzz = -zr * frcz 1116 vir(1,1) = vir(1,1) + vxx 1117 vir(2,1) = vir(2,1) + vxy 1118 vir(3,1) = vir(3,1) + vxz 1119 vir(1,2) = vir(1,2) + vxy 1120 vir(2,2) = vir(2,2) + vyy 1121 vir(3,2) = vir(3,2) + vyz 1122 vir(1,3) = vir(1,3) + vxz 1123 vir(2,3) = vir(2,3) + vyz 1124 vir(3,3) = vir(3,3) + vzz 1125 end if 1126 10 continue 1127 end do 1128c 1129c reset exclusion coefficients for connected atoms 1130c 1131 do j = 1, n12(ii) 1132 mscale(i12(j,ii)) = 1.0d0 1133 end do 1134 do j = 1, n13(ii) 1135 mscale(i13(j,ii)) = 1.0d0 1136 end do 1137 do j = 1, n14(ii) 1138 mscale(i14(j,ii)) = 1.0d0 1139 end do 1140 do j = 1, n15(ii) 1141 mscale(i15(j,ii)) = 1.0d0 1142 end do 1143 end do 1144c 1145c OpenMP directives for the major loop structure 1146c 1147!$OMP END DO 1148!$OMP DO reduction(+:dem,vir) schedule(guided) 1149c 1150c resolve site torques then increment forces and virial 1151c 1152 do i = 1, npole 1153 call torque (i,tem(1,i),fix,fiy,fiz,dem) 1154 ii = ipole(i) 1155 iaz = zaxis(i) 1156 iax = xaxis(i) 1157 iay = yaxis(i) 1158 if (iaz .eq. 0) iaz = ii 1159 if (iax .eq. 0) iax = ii 1160 if (iay .eq. 0) iay = ii 1161 xiz = x(iaz) - x(ii) 1162 yiz = y(iaz) - y(ii) 1163 ziz = z(iaz) - z(ii) 1164 xix = x(iax) - x(ii) 1165 yix = y(iax) - y(ii) 1166 zix = z(iax) - z(ii) 1167 xiy = x(iay) - x(ii) 1168 yiy = y(iay) - y(ii) 1169 ziy = z(iay) - z(ii) 1170 vxx = xix*fix(1) + xiy*fiy(1) + xiz*fiz(1) 1171 vxy = 0.5d0 * (yix*fix(1) + yiy*fiy(1) + yiz*fiz(1) 1172 & + xix*fix(2) + xiy*fiy(2) + xiz*fiz(2)) 1173 vxz = 0.5d0 * (zix*fix(1) + ziy*fiy(1) + ziz*fiz(1) 1174 & + xix*fix(3) + xiy*fiy(3) + xiz*fiz(3)) 1175 vyy = yix*fix(2) + yiy*fiy(2) + yiz*fiz(2) 1176 vyz = 0.5d0 * (zix*fix(2) + ziy*fiy(2) + ziz*fiz(2) 1177 & + yix*fix(3) + yiy*fiy(3) + yiz*fiz(3)) 1178 vzz = zix*fix(3) + ziy*fiy(3) + ziz*fiz(3) 1179 vir(1,1) = vir(1,1) + vxx 1180 vir(2,1) = vir(2,1) + vxy 1181 vir(3,1) = vir(3,1) + vxz 1182 vir(1,2) = vir(1,2) + vxy 1183 vir(2,2) = vir(2,2) + vyy 1184 vir(3,2) = vir(3,2) + vyz 1185 vir(1,3) = vir(1,3) + vxz 1186 vir(2,3) = vir(2,3) + vyz 1187 vir(3,3) = vir(3,3) + vzz 1188 end do 1189c 1190c OpenMP directives for the major loop structure 1191c 1192!$OMP END DO 1193!$OMP END PARALLEL 1194c 1195c perform deallocation of some local arrays 1196c 1197 deallocate (mscale) 1198 deallocate (tem) 1199 return 1200 end 1201c 1202c 1203c ################################################################ 1204c ## ## 1205c ## subroutine empole1c -- Ewald multipole derivs via loop ## 1206c ## ## 1207c ################################################################ 1208c 1209c 1210c "empole1c" calculates the multipole energy and derivatives 1211c with respect to Cartesian coordinates using particle mesh 1212c Ewald summation and a double loop 1213c 1214c 1215 subroutine empole1c 1216 use sizes 1217 use atoms 1218 use boxes 1219 use chgpot 1220 use deriv 1221 use energi 1222 use ewald 1223 use math 1224 use mpole 1225 use virial 1226 implicit none 1227 integer i,j,ii 1228 real*8 e,f 1229 real*8 term,fterm 1230 real*8 cii,dii,qii 1231 real*8 xd,yd,zd 1232 real*8 xq,yq,zq 1233 real*8 xv,yv,zv,vterm 1234 real*8 ci,dix,diy,diz 1235 real*8 qixx,qixy,qixz 1236 real*8 qiyy,qiyz,qizz 1237 real*8 xdfield,ydfield 1238 real*8 zdfield 1239 real*8 trq(3),frcx(3) 1240 real*8 frcy(3),frcz(3) 1241c 1242c 1243c zero out the atomic multipole energy and derivatives 1244c 1245 em = 0.0d0 1246 do i = 1, n 1247 do j = 1, 3 1248 dem(j,i) = 0.0d0 1249 end do 1250 end do 1251 if (npole .eq. 0) return 1252c 1253c set the energy unit conversion factor 1254c 1255 f = electric / dielec 1256c 1257c check the sign of multipole components at chiral sites 1258c 1259 call chkpole 1260c 1261c rotate the multipole components into the global frame 1262c 1263 call rotpole 1264c 1265c compute the real space part of the Ewald summation 1266c 1267 call emreal1c 1268c 1269c compute the reciprocal space part of the Ewald summation 1270c 1271 call emrecip1 1272c 1273c compute the Ewald self-energy term over all the atoms 1274c 1275 term = 2.0d0 * aewald * aewald 1276 fterm = -f * aewald / sqrtpi 1277 do i = 1, npole 1278 ci = rpole(1,i) 1279 dix = rpole(2,i) 1280 diy = rpole(3,i) 1281 diz = rpole(4,i) 1282 qixx = rpole(5,i) 1283 qixy = rpole(6,i) 1284 qixz = rpole(7,i) 1285 qiyy = rpole(9,i) 1286 qiyz = rpole(10,i) 1287 qizz = rpole(13,i) 1288 cii = ci*ci 1289 dii = dix*dix + diy*diy + diz*diz 1290 qii = 2.0d0*(qixy*qixy+qixz*qixz+qiyz*qiyz) 1291 & + qixx*qixx + qiyy*qiyy + qizz*qizz 1292 e = fterm * (cii + term*(dii/3.0d0+2.0d0*term*qii/5.0d0)) 1293 em = em + e 1294 end do 1295c 1296c compute the cell dipole boundary correction term 1297c 1298 if (boundary .eq. 'VACUUM') then 1299 xd = 0.0d0 1300 yd = 0.0d0 1301 zd = 0.0d0 1302 do i = 1, npole 1303 ii = ipole(i) 1304 xd = xd + rpole(2,i) + rpole(1,i)*x(ii) 1305 yd = yd + rpole(3,i) + rpole(1,i)*y(ii) 1306 zd = zd + rpole(4,i) + rpole(1,i)*z(ii) 1307 end do 1308 term = (2.0d0/3.0d0) * f * (pi/volbox) 1309 em = em + term*(xd*xd+yd*yd+zd*zd) 1310 do i = 1, npole 1311 ii = ipole(i) 1312 dem(1,ii) = dem(1,ii) + 2.0d0*term*rpole(1,i)*xd 1313 dem(2,ii) = dem(2,ii) + 2.0d0*term*rpole(1,i)*yd 1314 dem(3,ii) = dem(3,ii) + 2.0d0*term*rpole(1,i)*zd 1315 end do 1316 xdfield = -2.0d0 * term * xd 1317 ydfield = -2.0d0 * term * yd 1318 zdfield = -2.0d0 * term * zd 1319 do i = 1, npole 1320 trq(1) = rpole(3,i)*zdfield - rpole(4,i)*ydfield 1321 trq(2) = rpole(4,i)*xdfield - rpole(2,i)*zdfield 1322 trq(3) = rpole(2,i)*ydfield - rpole(3,i)*xdfield 1323 call torque (i,trq,frcx,frcy,frcz,dem) 1324 end do 1325c 1326c boundary correction to virial due to overall cell dipole 1327c 1328 xd = 0.0d0 1329 yd = 0.0d0 1330 zd = 0.0d0 1331 xq = 0.0d0 1332 yq = 0.0d0 1333 zq = 0.0d0 1334 do i = 1, npole 1335 ii = ipole(i) 1336 xd = xd + rpole(2,i) 1337 yd = yd + rpole(3,i) 1338 zd = zd + rpole(4,i) 1339 xq = xq + rpole(1,i)*x(ii) 1340 yq = yq + rpole(1,i)*y(ii) 1341 zq = zq + rpole(1,i)*z(ii) 1342 end do 1343 xv = xd * xq 1344 yv = yd * yq 1345 zv = zd * zq 1346 vterm = term * (xd*xd + yd*yd + zd*zd + 2.0d0*(xv+yv+zv) 1347 & + xq*xq + yq*yq + zq*zq) 1348 vir(1,1) = vir(1,1) + 2.0d0*term*(xq*xq+xv) + vterm 1349 vir(2,1) = vir(2,1) + 2.0d0*term*(xq*yq+xv) 1350 vir(3,1) = vir(3,1) + 2.0d0*term*(xq*zq+xv) 1351 vir(1,2) = vir(1,2) + 2.0d0*term*(yq*xq+yv) 1352 vir(2,2) = vir(2,2) + 2.0d0*term*(yq*yq+yv) + vterm 1353 vir(3,2) = vir(3,2) + 2.0d0*term*(yq*zq+yv) 1354 vir(1,3) = vir(1,3) + 2.0d0*term*(zq*xq+zv) 1355 vir(2,3) = vir(2,3) + 2.0d0*term*(zq*yq+zv) 1356 vir(3,3) = vir(3,3) + 2.0d0*term*(zq*zq+zv) + vterm 1357 end if 1358 return 1359 end 1360c 1361c 1362c ################################################################# 1363c ## ## 1364c ## subroutine emreal1c -- Ewald real space derivs via loop ## 1365c ## ## 1366c ################################################################# 1367c 1368c 1369c "emreal1c" evaluates the real space portion of the Ewald 1370c summation energy and gradient due to multipole interactions 1371c via a double loop 1372c 1373c 1374 subroutine emreal1c 1375 use sizes 1376 use atoms 1377 use bound 1378 use cell 1379 use chgpot 1380 use couple 1381 use deriv 1382 use energi 1383 use ewald 1384 use inter 1385 use math 1386 use molcul 1387 use mplpot 1388 use mpole 1389 use shunt 1390 use virial 1391 implicit none 1392 integer i,j,k 1393 integer ii,kk,jcell 1394 integer iax,iay,iaz 1395 real*8 e,efull,de,f 1396 real*8 bfac,erfc 1397 real*8 alsq2,alsq2n 1398 real*8 exp2a,ralpha 1399 real*8 scalekk 1400 real*8 xi,yi,zi 1401 real*8 xr,yr,zr 1402 real*8 xix,yix,zix 1403 real*8 xiy,yiy,ziy 1404 real*8 xiz,yiz,ziz 1405 real*8 r,r2,rr1,rr3 1406 real*8 rr5,rr7,rr9,rr11 1407 real*8 ci,dix,diy,diz 1408 real*8 qixx,qixy,qixz 1409 real*8 qiyy,qiyz,qizz 1410 real*8 ck,dkx,dky,dkz 1411 real*8 qkxx,qkxy,qkxz 1412 real*8 qkyy,qkyz,qkzz 1413 real*8 dikx,diky,dikz 1414 real*8 dirx,diry,dirz 1415 real*8 dkrx,dkry,dkrz 1416 real*8 qrix,qriy,qriz 1417 real*8 qrkx,qrky,qrkz 1418 real*8 qrixr,qriyr,qrizr 1419 real*8 qrkxr,qrkyr,qrkzr 1420 real*8 qrrx,qrry,qrrz 1421 real*8 qikrx,qikry,qikrz 1422 real*8 qkirx,qkiry,qkirz 1423 real*8 qikrxr,qikryr,qikrzr 1424 real*8 qkirxr,qkiryr,qkirzr 1425 real*8 diqkx,diqky,diqkz 1426 real*8 dkqix,dkqiy,dkqiz 1427 real*8 diqkxr,diqkyr,diqkzr 1428 real*8 dkqixr,dkqiyr,dkqizr 1429 real*8 dqiqkx,dqiqky,dqiqkz 1430 real*8 dri,drk,qrri,qrrk 1431 real*8 diqrk,dkqri 1432 real*8 dik,qik,qrrik 1433 real*8 term1,term2,term3 1434 real*8 term4,term5,term6 1435 real*8 frcx,frcy,frcz 1436 real*8 vxx,vyy,vzz 1437 real*8 vxy,vxz,vyz 1438 real*8 ttmi(3),ttmk(3) 1439 real*8 fix(3),fiy(3),fiz(3) 1440 real*8 bn(0:5) 1441 real*8, allocatable :: mscale(:) 1442 real*8, allocatable :: tem(:,:) 1443 character*6 mode 1444 external erfc 1445c 1446c 1447c perform dynamic allocation of some local arrays 1448c 1449 allocate (mscale(n)) 1450 allocate (tem(3,n)) 1451c 1452c initialize connected atom scaling and torque arrays 1453c 1454 do i = 1, n 1455 mscale(i) = 1.0d0 1456 do j = 1, 3 1457 tem(j,i) = 0.0d0 1458 end do 1459 end do 1460c 1461c set conversion factor, cutoff and switching coefficients 1462c 1463 f = electric / dielec 1464 mode = 'EWALD' 1465 call switch (mode) 1466c 1467c compute the real space portion of the Ewald summation 1468c 1469 do i = 1, npole-1 1470 ii = ipole(i) 1471 xi = x(ii) 1472 yi = y(ii) 1473 zi = z(ii) 1474 ci = rpole(1,i) 1475 dix = rpole(2,i) 1476 diy = rpole(3,i) 1477 diz = rpole(4,i) 1478 qixx = rpole(5,i) 1479 qixy = rpole(6,i) 1480 qixz = rpole(7,i) 1481 qiyy = rpole(9,i) 1482 qiyz = rpole(10,i) 1483 qizz = rpole(13,i) 1484 do j = 1, n12(ii) 1485 mscale(i12(j,ii)) = m2scale 1486 end do 1487 do j = 1, n13(ii) 1488 mscale(i13(j,ii)) = m3scale 1489 end do 1490 do j = 1, n14(ii) 1491 mscale(i14(j,ii)) = m4scale 1492 end do 1493 do j = 1, n15(ii) 1494 mscale(i15(j,ii)) = m5scale 1495 end do 1496c 1497c evaluate all sites within the cutoff distance 1498c 1499 do k = i+1, npole 1500 kk = ipole(k) 1501 xr = x(kk) - xi 1502 yr = y(kk) - yi 1503 zr = z(kk) - zi 1504 if (use_bounds) call image (xr,yr,zr) 1505 r2 = xr*xr + yr*yr + zr*zr 1506 if (r2 .le. off2) then 1507 r = sqrt(r2) 1508 ck = rpole(1,k) 1509 dkx = rpole(2,k) 1510 dky = rpole(3,k) 1511 dkz = rpole(4,k) 1512 qkxx = rpole(5,k) 1513 qkxy = rpole(6,k) 1514 qkxz = rpole(7,k) 1515 qkyy = rpole(9,k) 1516 qkyz = rpole(10,k) 1517 qkzz = rpole(13,k) 1518c 1519c get reciprocal distance terms for this interaction 1520c 1521 rr1 = f / r 1522 rr3 = rr1 / r2 1523 rr5 = 3.0d0 * rr3 / r2 1524 rr7 = 5.0d0 * rr5 / r2 1525 rr9 = 7.0d0 * rr7 / r2 1526 rr11 = 9.0d0 * rr9 / r2 1527c 1528c calculate the real space Ewald error function terms 1529c 1530 ralpha = aewald * r 1531 bn(0) = erfc(ralpha) / r 1532 alsq2 = 2.0d0 * aewald**2 1533 alsq2n = 0.0d0 1534 if (aewald .gt. 0.0d0) alsq2n = 1.0d0 / (sqrtpi*aewald) 1535 exp2a = exp(-ralpha**2) 1536 do j = 1, 5 1537 bfac = dble(j+j-1) 1538 alsq2n = alsq2 * alsq2n 1539 bn(j) = (bfac*bn(j-1)+alsq2n*exp2a) / r2 1540 end do 1541 do j = 0, 5 1542 bn(j) = f * bn(j) 1543 end do 1544c 1545c intermediates involving moments and distance separation 1546c 1547 dikx = diy*dkz - diz*dky 1548 diky = diz*dkx - dix*dkz 1549 dikz = dix*dky - diy*dkx 1550 dirx = diy*zr - diz*yr 1551 diry = diz*xr - dix*zr 1552 dirz = dix*yr - diy*xr 1553 dkrx = dky*zr - dkz*yr 1554 dkry = dkz*xr - dkx*zr 1555 dkrz = dkx*yr - dky*xr 1556 dri = dix*xr + diy*yr + diz*zr 1557 drk = dkx*xr + dky*yr + dkz*zr 1558 dik = dix*dkx + diy*dky + diz*dkz 1559 qrix = qixx*xr + qixy*yr + qixz*zr 1560 qriy = qixy*xr + qiyy*yr + qiyz*zr 1561 qriz = qixz*xr + qiyz*yr + qizz*zr 1562 qrkx = qkxx*xr + qkxy*yr + qkxz*zr 1563 qrky = qkxy*xr + qkyy*yr + qkyz*zr 1564 qrkz = qkxz*xr + qkyz*yr + qkzz*zr 1565 qrri = qrix*xr + qriy*yr + qriz*zr 1566 qrrk = qrkx*xr + qrky*yr + qrkz*zr 1567 qrrik = qrix*qrkx + qriy*qrky + qriz*qrkz 1568 qik = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 1569 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 1570 qrixr = qriz*yr - qriy*zr 1571 qriyr = qrix*zr - qriz*xr 1572 qrizr = qriy*xr - qrix*yr 1573 qrkxr = qrkz*yr - qrky*zr 1574 qrkyr = qrkx*zr - qrkz*xr 1575 qrkzr = qrky*xr - qrkx*yr 1576 qrrx = qrky*qriz - qrkz*qriy 1577 qrry = qrkz*qrix - qrkx*qriz 1578 qrrz = qrkx*qriy - qrky*qrix 1579 qikrx = qixx*qrkx + qixy*qrky + qixz*qrkz 1580 qikry = qixy*qrkx + qiyy*qrky + qiyz*qrkz 1581 qikrz = qixz*qrkx + qiyz*qrky + qizz*qrkz 1582 qkirx = qkxx*qrix + qkxy*qriy + qkxz*qriz 1583 qkiry = qkxy*qrix + qkyy*qriy + qkyz*qriz 1584 qkirz = qkxz*qrix + qkyz*qriy + qkzz*qriz 1585 qikrxr = qikrz*yr - qikry*zr 1586 qikryr = qikrx*zr - qikrz*xr 1587 qikrzr = qikry*xr - qikrx*yr 1588 qkirxr = qkirz*yr - qkiry*zr 1589 qkiryr = qkirx*zr - qkirz*xr 1590 qkirzr = qkiry*xr - qkirx*yr 1591 diqkx = dix*qkxx + diy*qkxy + diz*qkxz 1592 diqky = dix*qkxy + diy*qkyy + diz*qkyz 1593 diqkz = dix*qkxz + diy*qkyz + diz*qkzz 1594 dkqix = dkx*qixx + dky*qixy + dkz*qixz 1595 dkqiy = dkx*qixy + dky*qiyy + dkz*qiyz 1596 dkqiz = dkx*qixz + dky*qiyz + dkz*qizz 1597 diqrk = dix*qrkx + diy*qrky + diz*qrkz 1598 dkqri = dkx*qrix + dky*qriy + dkz*qriz 1599 diqkxr = diqkz*yr - diqky*zr 1600 diqkyr = diqkx*zr - diqkz*xr 1601 diqkzr = diqky*xr - diqkx*yr 1602 dkqixr = dkqiz*yr - dkqiy*zr 1603 dkqiyr = dkqix*zr - dkqiz*xr 1604 dkqizr = dkqiy*xr - dkqix*yr 1605 dqiqkx = diy*qrkz - diz*qrky + dky*qriz - dkz*qriy 1606 & - 2.0d0*(qixy*qkxz+qiyy*qkyz+qiyz*qkzz 1607 & -qixz*qkxy-qiyz*qkyy-qizz*qkyz) 1608 dqiqky = diz*qrkx - dix*qrkz + dkz*qrix - dkx*qriz 1609 & - 2.0d0*(qixz*qkxx+qiyz*qkxy+qizz*qkxz 1610 & -qixx*qkxz-qixy*qkyz-qixz*qkzz) 1611 dqiqkz = dix*qrky - diy*qrkx + dkx*qriy - dky*qrix 1612 & - 2.0d0*(qixx*qkxy+qixy*qkyy+qixz*qkyz 1613 & -qixy*qkxx-qiyy*qkxy-qiyz*qkxz) 1614c 1615c calculate intermediate terms for multipole energy 1616c 1617 term1 = ci*ck 1618 term2 = ck*dri - ci*drk + dik 1619 term3 = ci*qrrk + ck*qrri - dri*drk 1620 & + 2.0d0*(dkqri-diqrk+qik) 1621 term4 = dri*qrrk - drk*qrri - 4.0d0*qrrik 1622 term5 = qrri*qrrk 1623c 1624c compute the full energy without any Ewald scaling 1625c 1626 efull = term1*rr1 + term2*rr3 + term3*rr5 1627 & + term4*rr7 + term5*rr9 1628 efull = efull * mscale(kk) 1629 if (molcule(ii) .ne. molcule(kk)) 1630 & einter = einter + efull 1631c 1632c modify distances to account for Ewald and exclusions 1633c 1634 scalekk = 1.0d0 - mscale(kk) 1635 rr1 = bn(0) - scalekk*rr1 1636 rr3 = bn(1) - scalekk*rr3 1637 rr5 = bn(2) - scalekk*rr5 1638 rr7 = bn(3) - scalekk*rr7 1639 rr9 = bn(4) - scalekk*rr9 1640 rr11 = bn(5) - scalekk*rr11 1641c 1642c compute the energy contribution for this interaction 1643c 1644 e = term1*rr1 + term2*rr3 + term3*rr5 1645 & + term4*rr7 + term5*rr9 1646 em = em + e 1647c 1648c calculate intermediate terms for force and torque 1649c 1650 de = term1*rr3 + term2*rr5 + term3*rr7 1651 & + term4*rr9 + term5*rr11 1652 term1 = -ck*rr3 + drk*rr5 - qrrk*rr7 1653 term2 = ci*rr3 + dri*rr5 + qrri*rr7 1654 term3 = 2.0d0 * rr5 1655 term4 = 2.0d0 * (-ck*rr5+drk*rr7-qrrk*rr9) 1656 term5 = 2.0d0 * (-ci*rr5-dri*rr7-qrri*rr9) 1657 term6 = 4.0d0 * rr7 1658c 1659c compute the force components for this interaction 1660c 1661 frcx = de*xr + term1*dix + term2*dkx 1662 & + term3*(diqkx-dkqix) + term4*qrix 1663 & + term5*qrkx + term6*(qikrx+qkirx) 1664 frcy = de*yr + term1*diy + term2*dky 1665 & + term3*(diqky-dkqiy) + term4*qriy 1666 & + term5*qrky + term6*(qikry+qkiry) 1667 frcz = de*zr + term1*diz + term2*dkz 1668 & + term3*(diqkz-dkqiz) + term4*qriz 1669 & + term5*qrkz + term6*(qikrz+qkirz) 1670c 1671c compute the torque components for this interaction 1672c 1673 ttmi(1) = -rr3*dikx + term1*dirx + term3*(dqiqkx+dkqixr) 1674 & - term4*qrixr - term6*(qikrxr+qrrx) 1675 ttmi(2) = -rr3*diky + term1*diry + term3*(dqiqky+dkqiyr) 1676 & - term4*qriyr - term6*(qikryr+qrry) 1677 ttmi(3) = -rr3*dikz + term1*dirz + term3*(dqiqkz+dkqizr) 1678 & - term4*qrizr - term6*(qikrzr+qrrz) 1679 ttmk(1) = rr3*dikx + term2*dkrx - term3*(dqiqkx+diqkxr) 1680 & - term5*qrkxr - term6*(qkirxr-qrrx) 1681 ttmk(2) = rr3*diky + term2*dkry - term3*(dqiqky+diqkyr) 1682 & - term5*qrkyr - term6*(qkiryr-qrry) 1683 ttmk(3) = rr3*dikz + term2*dkrz - term3*(dqiqkz+diqkzr) 1684 & - term5*qrkzr - term6*(qkirzr-qrrz) 1685c 1686c increment force-based gradient and torque on first site 1687c 1688 dem(1,ii) = dem(1,ii) + frcx 1689 dem(2,ii) = dem(2,ii) + frcy 1690 dem(3,ii) = dem(3,ii) + frcz 1691 tem(1,i) = tem(1,i) + ttmi(1) 1692 tem(2,i) = tem(2,i) + ttmi(2) 1693 tem(3,i) = tem(3,i) + ttmi(3) 1694c 1695c increment force-based gradient and torque on second site 1696c 1697 dem(1,kk) = dem(1,kk) - frcx 1698 dem(2,kk) = dem(2,kk) - frcy 1699 dem(3,kk) = dem(3,kk) - frcz 1700 tem(1,k) = tem(1,k) + ttmk(1) 1701 tem(2,k) = tem(2,k) + ttmk(2) 1702 tem(3,k) = tem(3,k) + ttmk(3) 1703c 1704c increment the virial due to pairwise Cartesian forces 1705c 1706 vxx = -xr * frcx 1707 vxy = -yr * frcx 1708 vxz = -zr * frcx 1709 vyy = -yr * frcy 1710 vyz = -zr * frcy 1711 vzz = -zr * frcz 1712 vir(1,1) = vir(1,1) + vxx 1713 vir(2,1) = vir(2,1) + vxy 1714 vir(3,1) = vir(3,1) + vxz 1715 vir(1,2) = vir(1,2) + vxy 1716 vir(2,2) = vir(2,2) + vyy 1717 vir(3,2) = vir(3,2) + vyz 1718 vir(1,3) = vir(1,3) + vxz 1719 vir(2,3) = vir(2,3) + vyz 1720 vir(3,3) = vir(3,3) + vzz 1721 end if 1722 end do 1723c 1724c reset exclusion coefficients for connected atoms 1725c 1726 do j = 1, n12(ii) 1727 mscale(i12(j,ii)) = 1.0d0 1728 end do 1729 do j = 1, n13(ii) 1730 mscale(i13(j,ii)) = 1.0d0 1731 end do 1732 do j = 1, n14(ii) 1733 mscale(i14(j,ii)) = 1.0d0 1734 end do 1735 do j = 1, n15(ii) 1736 mscale(i15(j,ii)) = 1.0d0 1737 end do 1738 end do 1739c 1740c for periodic boundary conditions with large cutoffs 1741c neighbors must be found by the replicates method 1742c 1743 if (use_replica) then 1744c 1745c calculate interaction with other unit cells 1746c 1747 do i = 1, npole 1748 ii = ipole(i) 1749 xi = x(ii) 1750 yi = y(ii) 1751 zi = z(ii) 1752 ci = rpole(1,i) 1753 dix = rpole(2,i) 1754 diy = rpole(3,i) 1755 diz = rpole(4,i) 1756 qixx = rpole(5,i) 1757 qixy = rpole(6,i) 1758 qixz = rpole(7,i) 1759 qiyy = rpole(9,i) 1760 qiyz = rpole(10,i) 1761 qizz = rpole(13,i) 1762 do j = 1, n12(ii) 1763 mscale(i12(j,ii)) = m2scale 1764 end do 1765 do j = 1, n13(ii) 1766 mscale(i13(j,ii)) = m3scale 1767 end do 1768 do j = 1, n14(ii) 1769 mscale(i14(j,ii)) = m4scale 1770 end do 1771 do j = 1, n15(ii) 1772 mscale(i15(j,ii)) = m5scale 1773 end do 1774c 1775c evaluate all sites within the cutoff distance 1776c 1777 do k = i, npole 1778 kk = ipole(k) 1779 do jcell = 2, ncell 1780 xr = x(kk) - xi 1781 yr = y(kk) - yi 1782 zr = z(kk) - zi 1783 call imager (xr,yr,zr,jcell) 1784 r2 = xr*xr + yr*yr + zr*zr 1785 if (.not. (use_polymer .and. r2.le.polycut2)) then 1786 mscale(kk) = 1.0d0 1787 end if 1788 if (r2 .le. off2) then 1789 r = sqrt(r2) 1790 ck = rpole(1,k) 1791 dkx = rpole(2,k) 1792 dky = rpole(3,k) 1793 dkz = rpole(4,k) 1794 qkxx = rpole(5,k) 1795 qkxy = rpole(6,k) 1796 qkxz = rpole(7,k) 1797 qkyy = rpole(9,k) 1798 qkyz = rpole(10,k) 1799 qkzz = rpole(13,k) 1800c 1801c get reciprocal distance terms for this interaction 1802c 1803 rr1 = f / r 1804 rr3 = rr1 / r2 1805 rr5 = 3.0d0 * rr3 / r2 1806 rr7 = 5.0d0 * rr5 / r2 1807 rr9 = 7.0d0 * rr7 / r2 1808 rr11 = 9.0d0 * rr9 / r2 1809c 1810c calculate the real space Ewald error function terms 1811c 1812 ralpha = aewald * r 1813 bn(0) = erfc(ralpha) / r 1814 alsq2 = 2.0d0 * aewald**2 1815 alsq2n = 0.0d0 1816 if (aewald .gt. 0.0d0) alsq2n = 1.0d0 / (sqrtpi*aewald) 1817 exp2a = exp(-ralpha**2) 1818 do j = 1, 5 1819 bfac = dble(j+j-1) 1820 alsq2n = alsq2 * alsq2n 1821 bn(j) = (bfac*bn(j-1)+alsq2n*exp2a) / r2 1822 end do 1823 do j = 0, 5 1824 bn(j) = f * bn(j) 1825 end do 1826c 1827c intermediates involving moments and distance separation 1828c 1829 dikx = diy*dkz - diz*dky 1830 diky = diz*dkx - dix*dkz 1831 dikz = dix*dky - diy*dkx 1832 dirx = diy*zr - diz*yr 1833 diry = diz*xr - dix*zr 1834 dirz = dix*yr - diy*xr 1835 dkrx = dky*zr - dkz*yr 1836 dkry = dkz*xr - dkx*zr 1837 dkrz = dkx*yr - dky*xr 1838 dri = dix*xr + diy*yr + diz*zr 1839 drk = dkx*xr + dky*yr + dkz*zr 1840 dik = dix*dkx + diy*dky + diz*dkz 1841 qrix = qixx*xr + qixy*yr + qixz*zr 1842 qriy = qixy*xr + qiyy*yr + qiyz*zr 1843 qriz = qixz*xr + qiyz*yr + qizz*zr 1844 qrkx = qkxx*xr + qkxy*yr + qkxz*zr 1845 qrky = qkxy*xr + qkyy*yr + qkyz*zr 1846 qrkz = qkxz*xr + qkyz*yr + qkzz*zr 1847 qrri = qrix*xr + qriy*yr + qriz*zr 1848 qrrk = qrkx*xr + qrky*yr + qrkz*zr 1849 qrrik = qrix*qrkx + qriy*qrky + qriz*qrkz 1850 qik = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 1851 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 1852 qrixr = qriz*yr - qriy*zr 1853 qriyr = qrix*zr - qriz*xr 1854 qrizr = qriy*xr - qrix*yr 1855 qrkxr = qrkz*yr - qrky*zr 1856 qrkyr = qrkx*zr - qrkz*xr 1857 qrkzr = qrky*xr - qrkx*yr 1858 qrrx = qrky*qriz - qrkz*qriy 1859 qrry = qrkz*qrix - qrkx*qriz 1860 qrrz = qrkx*qriy - qrky*qrix 1861 qikrx = qixx*qrkx + qixy*qrky + qixz*qrkz 1862 qikry = qixy*qrkx + qiyy*qrky + qiyz*qrkz 1863 qikrz = qixz*qrkx + qiyz*qrky + qizz*qrkz 1864 qkirx = qkxx*qrix + qkxy*qriy + qkxz*qriz 1865 qkiry = qkxy*qrix + qkyy*qriy + qkyz*qriz 1866 qkirz = qkxz*qrix + qkyz*qriy + qkzz*qriz 1867 qikrxr = qikrz*yr - qikry*zr 1868 qikryr = qikrx*zr - qikrz*xr 1869 qikrzr = qikry*xr - qikrx*yr 1870 qkirxr = qkirz*yr - qkiry*zr 1871 qkiryr = qkirx*zr - qkirz*xr 1872 qkirzr = qkiry*xr - qkirx*yr 1873 diqkx = dix*qkxx + diy*qkxy + diz*qkxz 1874 diqky = dix*qkxy + diy*qkyy + diz*qkyz 1875 diqkz = dix*qkxz + diy*qkyz + diz*qkzz 1876 dkqix = dkx*qixx + dky*qixy + dkz*qixz 1877 dkqiy = dkx*qixy + dky*qiyy + dkz*qiyz 1878 dkqiz = dkx*qixz + dky*qiyz + dkz*qizz 1879 diqrk = dix*qrkx + diy*qrky + diz*qrkz 1880 dkqri = dkx*qrix + dky*qriy + dkz*qriz 1881 diqkxr = diqkz*yr - diqky*zr 1882 diqkyr = diqkx*zr - diqkz*xr 1883 diqkzr = diqky*xr - diqkx*yr 1884 dkqixr = dkqiz*yr - dkqiy*zr 1885 dkqiyr = dkqix*zr - dkqiz*xr 1886 dkqizr = dkqiy*xr - dkqix*yr 1887 dqiqkx = diy*qrkz - diz*qrky + dky*qriz - dkz*qriy 1888 & - 2.0d0*(qixy*qkxz+qiyy*qkyz+qiyz*qkzz 1889 & -qixz*qkxy-qiyz*qkyy-qizz*qkyz) 1890 dqiqky = diz*qrkx - dix*qrkz + dkz*qrix - dkx*qriz 1891 & - 2.0d0*(qixz*qkxx+qiyz*qkxy+qizz*qkxz 1892 & -qixx*qkxz-qixy*qkyz-qixz*qkzz) 1893 dqiqkz = dix*qrky - diy*qrkx + dkx*qriy - dky*qrix 1894 & - 2.0d0*(qixx*qkxy+qixy*qkyy+qixz*qkyz 1895 & -qixy*qkxx-qiyy*qkxy-qiyz*qkxz) 1896c 1897c calculate intermediate terms for multipole energy 1898c 1899 term1 = ci*ck 1900 term2 = ck*dri - ci*drk + dik 1901 term3 = ci*qrrk + ck*qrri - dri*drk 1902 & + 2.0d0*(dkqri-diqrk+qik) 1903 term4 = dri*qrrk - drk*qrri - 4.0d0*qrrik 1904 term5 = qrri*qrrk 1905c 1906c compute the full energy without any Ewald scaling 1907c 1908 efull = term1*rr1 + term2*rr3 + term3*rr5 1909 & + term4*rr7 + term5*rr9 1910 efull = efull * mscale(kk) 1911 if (ii .eq. kk) efull = 0.5d0 * e 1912 einter = einter + efull 1913c 1914c modify distances to account for Ewald and exclusions 1915c 1916 scalekk = 1.0d0 - mscale(kk) 1917 rr1 = bn(0) - scalekk*rr1 1918 rr3 = bn(1) - scalekk*rr3 1919 rr5 = bn(2) - scalekk*rr5 1920 rr7 = bn(3) - scalekk*rr7 1921 rr9 = bn(4) - scalekk*rr9 1922 rr11 = bn(5) - scalekk*rr11 1923c 1924c compute the energy contribution for this interaction 1925c 1926 e = term1*rr1 + term2*rr3 + term3*rr5 1927 & + term4*rr7 + term5*rr9 1928 if (ii .eq. kk) e = 0.5d0 * e 1929 em = em + e 1930c 1931c calculate intermediate terms for force and torque 1932c 1933 de = term1*rr3 + term2*rr5 + term3*rr7 1934 & + term4*rr9 + term5*rr11 1935 term1 = -ck*rr3 + drk*rr5 - qrrk*rr7 1936 term2 = ci*rr3 + dri*rr5 + qrri*rr7 1937 term3 = 2.0d0 * rr5 1938 term4 = 2.0d0 * (-ck*rr5+drk*rr7-qrrk*rr9) 1939 term5 = 2.0d0 * (-ci*rr5-dri*rr7-qrri*rr9) 1940 term6 = 4.0d0 * rr7 1941c 1942c compute the force components for this interaction 1943c 1944 frcx = de*xr + term1*dix + term2*dkx 1945 & + term3*(diqkx-dkqix) + term4*qrix 1946 & + term5*qrkx + term6*(qikrx+qkirx) 1947 frcy = de*yr + term1*diy + term2*dky 1948 & + term3*(diqky-dkqiy) + term4*qriy 1949 & + term5*qrky + term6*(qikry+qkiry) 1950 frcz = de*zr + term1*diz + term2*dkz 1951 & + term3*(diqkz-dkqiz) + term4*qriz 1952 & + term5*qrkz + term6*(qikrz+qkirz) 1953c 1954c compute the torque components for this interaction 1955c 1956 ttmi(1) = -rr3*dikx + term1*dirx + term3*(dqiqkx+dkqixr) 1957 & - term4*qrixr - term6*(qikrxr+qrrx) 1958 ttmi(2) = -rr3*diky + term1*diry + term3*(dqiqky+dkqiyr) 1959 & - term4*qriyr - term6*(qikryr+qrry) 1960 ttmi(3) = -rr3*dikz + term1*dirz + term3*(dqiqkz+dkqizr) 1961 & - term4*qrizr - term6*(qikrzr+qrrz) 1962 ttmk(1) = rr3*dikx + term2*dkrx - term3*(dqiqkx+diqkxr) 1963 & - term5*qrkxr - term6*(qkirxr-qrrx) 1964 ttmk(2) = rr3*diky + term2*dkry - term3*(dqiqky+diqkyr) 1965 & - term5*qrkyr - term6*(qkiryr-qrry) 1966 ttmk(3) = rr3*dikz + term2*dkrz - term3*(dqiqkz+diqkzr) 1967 & - term5*qrkzr - term6*(qkirzr-qrrz) 1968c 1969c force and torque components scaled for self-interactions 1970c 1971 if (ii .eq. kk) then 1972 frcx = 0.5d0 * frcx 1973 frcy = 0.5d0 * frcy 1974 frcz = 0.5d0 * frcz 1975 do j = 1, 3 1976 ttmi(j) = 0.5d0 * ttmi(j) 1977 ttmk(j) = 0.5d0 * ttmk(j) 1978 end do 1979 end if 1980c 1981c increment force-based gradient and torque on first site 1982c 1983 dem(1,ii) = dem(1,ii) + frcx 1984 dem(2,ii) = dem(2,ii) + frcy 1985 dem(3,ii) = dem(3,ii) + frcz 1986 tem(1,i) = tem(1,i) + ttmi(1) 1987 tem(2,i) = tem(2,i) + ttmi(2) 1988 tem(3,i) = tem(3,i) + ttmi(3) 1989c 1990c increment force-based gradient and torque on second site 1991c 1992 dem(1,kk) = dem(1,kk) - frcx 1993 dem(2,kk) = dem(2,kk) - frcy 1994 dem(3,kk) = dem(3,kk) - frcz 1995 tem(1,k) = tem(1,k) + ttmk(1) 1996 tem(2,k) = tem(2,k) + ttmk(2) 1997 tem(3,k) = tem(3,k) + ttmk(3) 1998c 1999c increment the virial due to pairwise Cartesian forces 2000c 2001 vxx = -xr * frcx 2002 vxy = -yr * frcx 2003 vxz = -zr * frcx 2004 vyy = -yr * frcy 2005 vyz = -zr * frcy 2006 vzz = -zr * frcz 2007 vir(1,1) = vir(1,1) + vxx 2008 vir(2,1) = vir(2,1) + vxy 2009 vir(3,1) = vir(3,1) + vxz 2010 vir(1,2) = vir(1,2) + vxy 2011 vir(2,2) = vir(2,2) + vyy 2012 vir(3,2) = vir(3,2) + vyz 2013 vir(1,3) = vir(1,3) + vxz 2014 vir(2,3) = vir(2,3) + vyz 2015 vir(3,3) = vir(3,3) + vzz 2016 end if 2017 end do 2018 end do 2019c 2020c reset exclusion coefficients for connected atoms 2021c 2022 do j = 1, n12(ii) 2023 mscale(i12(j,ii)) = 1.0d0 2024 end do 2025 do j = 1, n13(ii) 2026 mscale(i13(j,ii)) = 1.0d0 2027 end do 2028 do j = 1, n14(ii) 2029 mscale(i14(j,ii)) = 1.0d0 2030 end do 2031 do j = 1, n15(ii) 2032 mscale(i15(j,ii)) = 1.0d0 2033 end do 2034 end do 2035 end if 2036c 2037c resolve site torques then increment forces and virial 2038c 2039 do i = 1, npole 2040 call torque (i,tem(1,i),fix,fiy,fiz,dem) 2041 ii = ipole(i) 2042 iaz = zaxis(i) 2043 iax = xaxis(i) 2044 iay = yaxis(i) 2045 if (iaz .eq. 0) iaz = ii 2046 if (iax .eq. 0) iax = ii 2047 if (iay .eq. 0) iay = ii 2048 xiz = x(iaz) - x(ii) 2049 yiz = y(iaz) - y(ii) 2050 ziz = z(iaz) - z(ii) 2051 xix = x(iax) - x(ii) 2052 yix = y(iax) - y(ii) 2053 zix = z(iax) - z(ii) 2054 xiy = x(iay) - x(ii) 2055 yiy = y(iay) - y(ii) 2056 ziy = z(iay) - z(ii) 2057 2058c vxx = xix*fix(1) + xiy*fiy(1) + xiz*fiz(1) 2059c vxy = 0.5d0 * (yix*fix(1) + yiy*fiy(1) + yiz*fiz(1) 2060c & + xix*fix(2) + xiy*fiy(2) + xiz*fiz(2)) 2061c vxz = 0.5d0 * (zix*fix(1) + ziy*fiy(1) + ziz*fiz(1) 2062c & + xix*fix(3) + xiy*fiy(3) + xiz*fiz(3)) 2063c vyy = yix*fix(2) + yiy*fiy(2) + yiz*fiz(2) 2064c vyz = 0.5d0 * (zix*fix(2) + ziy*fiy(2) + ziz*fiz(2) 2065c & + yix*fix(3) + yiy*fiy(3) + yiz*fiz(3)) 2066c vzz = zix*fix(3) + ziy*fiy(3) + ziz*fiz(3) 2067 2068 vxx = xix*fix(1) + xiy*fiy(1) + xiz*fiz(1) 2069c vxy = 0.5d0*(yix*fix(1) + yiy*fiy(1) + yiz*fiz(1) 2070c & + xix*fix(2) + xiy*fiy(2) + xiz*fiz(2)) 2071 vxy = 0.5d0*(xiy*fix(1) + yiy*fiy(1) + ziy*fiz(1) 2072 & + xix*fix(2) + yix*fiy(2) + zix*fiz(2)) 2073c vxz = 0.5d0*(zix*fix(1) + ziy*fiy(1) + ziz*fiz(1) 2074c & + xix*fix(3) + xiy*fiy(3) + xiz*fiz(3)) 2075 vxz = 0.5d0*(xiz*fix(1) + yiz*fiy(1) + ziz*fiz(1) 2076 & + xix*fix(3) + yix*fiy(3) + zix*fiz(3)) 2077 vyy = yix*fix(2) + yiy*fiy(2) + yiz*fiz(2) 2078c vyz = 0.5d0*(zix*fix(2) + ziy*fiy(2) + ziz*fiz(2) 2079c & + yix*fix(3) + yiy*fiy(3) + yiz*fiz(3)) 2080 vyz = 0.5d0*(xiz*fix(2) + yiz*fiy(2) + ziz*fiz(2) 2081 & + xiy*fix(3) + yiy*fiy(3) + ziy*fiz(3)) 2082 vzz = zix*fix(3) + ziy*fiy(3) + ziz*fiz(3) 2083 2084 vir(1,1) = vir(1,1) + vxx 2085 vir(2,1) = vir(2,1) + vxy 2086 vir(3,1) = vir(3,1) + vxz 2087 vir(1,2) = vir(1,2) + vxy 2088 vir(2,2) = vir(2,2) + vyy 2089 vir(3,2) = vir(3,2) + vyz 2090 vir(1,3) = vir(1,3) + vxz 2091 vir(2,3) = vir(2,3) + vyz 2092 vir(3,3) = vir(3,3) + vzz 2093 end do 2094c 2095c perform deallocation of some local arrays 2096c 2097 deallocate (mscale) 2098 deallocate (tem) 2099 return 2100 end 2101c 2102c 2103c ################################################################ 2104c ## ## 2105c ## subroutine empole1d -- Ewald multipole derivs via list ## 2106c ## ## 2107c ################################################################ 2108c 2109c 2110c "empole1d" calculates the multipole energy and derivatives 2111c with respect to Cartesian coordinates using particle mesh Ewald 2112c summation and a neighbor list 2113c 2114c 2115 subroutine empole1d 2116 use sizes 2117 use atoms 2118 use boxes 2119 use chgpot 2120 use deriv 2121 use energi 2122 use ewald 2123 use math 2124 use mpole 2125 use virial 2126 implicit none 2127 integer i,j,ii 2128 real*8 e,f 2129 real*8 term,fterm 2130 real*8 cii,dii,qii 2131 real*8 xd,yd,zd 2132 real*8 xq,yq,zq 2133 real*8 xv,yv,zv,vterm 2134 real*8 ci,dix,diy,diz 2135 real*8 qixx,qixy,qixz 2136 real*8 qiyy,qiyz,qizz 2137 real*8 xdfield,ydfield 2138 real*8 zdfield 2139 real*8 trq(3),frcx(3) 2140 real*8 frcy(3),frcz(3) 2141c 2142c 2143c zero out the atomic multipole energy and derivatives 2144c 2145 em = 0.0d0 2146 do i = 1, n 2147 do j = 1, 3 2148 dem(j,i) = 0.0d0 2149 end do 2150 end do 2151 if (npole .eq. 0) return 2152c 2153c set the energy unit conversion factor 2154c 2155 f = electric / dielec 2156c 2157c check the sign of multipole components at chiral sites 2158c 2159 call chkpole 2160c 2161c rotate the multipole components into the global frame 2162c 2163 call rotpole 2164c 2165c compute the real space part of the Ewald summation 2166c 2167 call emreal1d 2168c 2169c compute the reciprocal space part of the Ewald summation 2170c 2171 call emrecip1 2172c 2173c compute the Ewald self-energy term over all the atoms 2174c 2175 term = 2.0d0 * aewald * aewald 2176 fterm = -f * aewald / sqrtpi 2177 do i = 1, npole 2178 ci = rpole(1,i) 2179 dix = rpole(2,i) 2180 diy = rpole(3,i) 2181 diz = rpole(4,i) 2182 qixx = rpole(5,i) 2183 qixy = rpole(6,i) 2184 qixz = rpole(7,i) 2185 qiyy = rpole(9,i) 2186 qiyz = rpole(10,i) 2187 qizz = rpole(13,i) 2188 cii = ci*ci 2189 dii = dix*dix + diy*diy + diz*diz 2190 qii = 2.0d0*(qixy*qixy+qixz*qixz+qiyz*qiyz) 2191 & + qixx*qixx + qiyy*qiyy + qizz*qizz 2192 e = fterm * (cii + term*(dii/3.0d0+2.0d0*term*qii/5.0d0)) 2193 em = em + e 2194 end do 2195c 2196c compute the cell dipole boundary correction term 2197c 2198 if (boundary .eq. 'VACUUM') then 2199 xd = 0.0d0 2200 yd = 0.0d0 2201 zd = 0.0d0 2202 do i = 1, npole 2203 ii = ipole(i) 2204 xd = xd + rpole(2,i) + rpole(1,i)*x(ii) 2205 yd = yd + rpole(3,i) + rpole(1,i)*y(ii) 2206 zd = zd + rpole(4,i) + rpole(1,i)*z(ii) 2207 end do 2208 term = (2.0d0/3.0d0) * f * (pi/volbox) 2209 em = em + term*(xd*xd+yd*yd+zd*zd) 2210 do i = 1, npole 2211 ii = ipole(i) 2212 dem(1,ii) = dem(1,ii) + 2.0d0*term*rpole(1,i)*xd 2213 dem(2,ii) = dem(2,ii) + 2.0d0*term*rpole(1,i)*yd 2214 dem(3,ii) = dem(3,ii) + 2.0d0*term*rpole(1,i)*zd 2215 end do 2216 xdfield = -2.0d0 * term * xd 2217 ydfield = -2.0d0 * term * yd 2218 zdfield = -2.0d0 * term * zd 2219 do i = 1, npole 2220 trq(1) = rpole(3,i)*zdfield - rpole(4,i)*ydfield 2221 trq(2) = rpole(4,i)*xdfield - rpole(2,i)*zdfield 2222 trq(3) = rpole(2,i)*ydfield - rpole(3,i)*xdfield 2223 call torque (i,trq,frcx,frcy,frcz,dem) 2224 end do 2225c 2226c boundary correction to virial due to overall cell dipole 2227c 2228 xd = 0.0d0 2229 yd = 0.0d0 2230 zd = 0.0d0 2231 xq = 0.0d0 2232 yq = 0.0d0 2233 zq = 0.0d0 2234 do i = 1, npole 2235 ii = ipole(i) 2236 xd = xd + rpole(2,i) 2237 yd = yd + rpole(3,i) 2238 zd = zd + rpole(4,i) 2239 xq = xq + rpole(1,i)*x(ii) 2240 yq = yq + rpole(1,i)*y(ii) 2241 zq = zq + rpole(1,i)*z(ii) 2242 end do 2243 xv = xd * xq 2244 yv = yd * yq 2245 zv = zd * zq 2246 vterm = term * (xd*xd + yd*yd + zd*zd + 2.0d0*(xv+yv+zv) 2247 & + xq*xq + yq*yq + zq*zq) 2248 vir(1,1) = vir(1,1) + 2.0d0*term*(xq*xq+xv) + vterm 2249 vir(2,1) = vir(2,1) + 2.0d0*term*(xq*yq+xv) 2250 vir(3,1) = vir(3,1) + 2.0d0*term*(xq*zq+xv) 2251 vir(1,2) = vir(1,2) + 2.0d0*term*(yq*xq+yv) 2252 vir(2,2) = vir(2,2) + 2.0d0*term*(yq*yq+yv) + vterm 2253 vir(3,2) = vir(3,2) + 2.0d0*term*(yq*zq+yv) 2254 vir(1,3) = vir(1,3) + 2.0d0*term*(zq*xq+zv) 2255 vir(2,3) = vir(2,3) + 2.0d0*term*(zq*yq+zv) 2256 vir(3,3) = vir(3,3) + 2.0d0*term*(zq*zq+zv) + vterm 2257 end if 2258 return 2259 end 2260c 2261c 2262c ################################################################# 2263c ## ## 2264c ## subroutine emreal1d -- Ewald real space derivs via list ## 2265c ## ## 2266c ################################################################# 2267c 2268c 2269c "emreal1d" evaluates the real space portion of the Ewald 2270c summation energy and gradient due to multipole interactions 2271c via a neighbor list 2272c 2273c 2274 subroutine emreal1d 2275 use sizes 2276 use atoms 2277 use bound 2278 use chgpot 2279 use couple 2280 use deriv 2281 use energi 2282 use ewald 2283 use inter 2284 use math 2285 use molcul 2286 use mplpot 2287 use mpole 2288 use neigh 2289 use shunt 2290 use virial 2291 implicit none 2292 integer i,j,k 2293 integer ii,kk,kkk 2294 integer iax,iay,iaz 2295 real*8 e,efull,de,f 2296 real*8 bfac,erfc 2297 real*8 alsq2,alsq2n 2298 real*8 exp2a,ralpha 2299 real*8 scalekk 2300 real*8 xi,yi,zi 2301 real*8 xr,yr,zr 2302 real*8 xix,yix,zix 2303 real*8 xiy,yiy,ziy 2304 real*8 xiz,yiz,ziz 2305 real*8 r,r2,rr1,rr3 2306 real*8 rr5,rr7,rr9,rr11 2307 real*8 ci,dix,diy,diz 2308 real*8 qixx,qixy,qixz 2309 real*8 qiyy,qiyz,qizz 2310 real*8 ck,dkx,dky,dkz 2311 real*8 qkxx,qkxy,qkxz 2312 real*8 qkyy,qkyz,qkzz 2313 real*8 dikx,diky,dikz 2314 real*8 dirx,diry,dirz 2315 real*8 dkrx,dkry,dkrz 2316 real*8 qrix,qriy,qriz 2317 real*8 qrkx,qrky,qrkz 2318 real*8 qrixr,qriyr,qrizr 2319 real*8 qrkxr,qrkyr,qrkzr 2320 real*8 qrrx,qrry,qrrz 2321 real*8 qikrx,qikry,qikrz 2322 real*8 qkirx,qkiry,qkirz 2323 real*8 qikrxr,qikryr,qikrzr 2324 real*8 qkirxr,qkiryr,qkirzr 2325 real*8 diqkx,diqky,diqkz 2326 real*8 dkqix,dkqiy,dkqiz 2327 real*8 diqkxr,diqkyr,diqkzr 2328 real*8 dkqixr,dkqiyr,dkqizr 2329 real*8 dqiqkx,dqiqky,dqiqkz 2330 real*8 dri,drk,qrri,qrrk 2331 real*8 diqrk,dkqri 2332 real*8 dik,qik,qrrik 2333 real*8 term1,term2,term3 2334 real*8 term4,term5,term6 2335 real*8 frcx,frcy,frcz 2336 real*8 vxx,vyy,vzz 2337 real*8 vxy,vxz,vyz 2338 real*8 ttmi(3),ttmk(3) 2339 real*8 fix(3),fiy(3),fiz(3) 2340 real*8 bn(0:5) 2341 real*8, allocatable :: mscale(:) 2342 real*8, allocatable :: tem(:,:) 2343 character*6 mode 2344 external erfc 2345c 2346c 2347c perform dynamic allocation of some local arrays 2348c 2349 allocate (mscale(n)) 2350 allocate (tem(3,n)) 2351c 2352c initialize connected atom scaling and torque arrays 2353c 2354 do i = 1, n 2355 mscale(i) = 1.0d0 2356 do j = 1, 3 2357 tem(j,i) = 0.0d0 2358 end do 2359 end do 2360c 2361c set conversion factor, cutoff and switching coefficients 2362c 2363 f = electric / dielec 2364 mode = 'EWALD' 2365 call switch (mode) 2366c 2367c OpenMP directives for the major loop structure 2368c 2369!$OMP PARALLEL default(private) 2370!$OMP& shared(npole,ipole,x,y,z,rpole,n12,i12,n13,i13,n14,i14, 2371!$OMP& n15,i15,m2scale,m3scale,m4scale,m5scale,nelst,elst, 2372!$OMP& use_bounds,f,off2,aewald,molcule,xaxis,yaxis,zaxis) 2373!$OMP& firstprivate(mscale) shared (em,einter,dem,tem,vir) 2374!$OMP DO reduction(+:em,einter,dem,tem,vir) schedule(guided) 2375c 2376c compute the real space portion of the Ewald summation 2377c 2378 do i = 1, npole 2379 ii = ipole(i) 2380 xi = x(ii) 2381 yi = y(ii) 2382 zi = z(ii) 2383 ci = rpole(1,i) 2384 dix = rpole(2,i) 2385 diy = rpole(3,i) 2386 diz = rpole(4,i) 2387 qixx = rpole(5,i) 2388 qixy = rpole(6,i) 2389 qixz = rpole(7,i) 2390 qiyy = rpole(9,i) 2391 qiyz = rpole(10,i) 2392 qizz = rpole(13,i) 2393 do j = 1, n12(ii) 2394 mscale(i12(j,ii)) = m2scale 2395 end do 2396 do j = 1, n13(ii) 2397 mscale(i13(j,ii)) = m3scale 2398 end do 2399 do j = 1, n14(ii) 2400 mscale(i14(j,ii)) = m4scale 2401 end do 2402 do j = 1, n15(ii) 2403 mscale(i15(j,ii)) = m5scale 2404 end do 2405c 2406c evaluate all sites within the cutoff distance 2407c 2408 do kkk = 1, nelst(i) 2409 k = elst(kkk,i) 2410 kk = ipole(k) 2411 xr = x(kk) - xi 2412 yr = y(kk) - yi 2413 zr = z(kk) - zi 2414 if (use_bounds) call image (xr,yr,zr) 2415 r2 = xr*xr + yr*yr + zr*zr 2416 if (r2 .le. off2) then 2417 r = sqrt(r2) 2418 ck = rpole(1,k) 2419 dkx = rpole(2,k) 2420 dky = rpole(3,k) 2421 dkz = rpole(4,k) 2422 qkxx = rpole(5,k) 2423 qkxy = rpole(6,k) 2424 qkxz = rpole(7,k) 2425 qkyy = rpole(9,k) 2426 qkyz = rpole(10,k) 2427 qkzz = rpole(13,k) 2428c 2429c get reciprocal distance terms for this interaction 2430c 2431 rr1 = f / r 2432 rr3 = rr1 / r2 2433 rr5 = 3.0d0 * rr3 / r2 2434 rr7 = 5.0d0 * rr5 / r2 2435 rr9 = 7.0d0 * rr7 / r2 2436 rr11 = 9.0d0 * rr9 / r2 2437c 2438c calculate the real space Ewald error function terms 2439c 2440 ralpha = aewald * r 2441 bn(0) = erfc(ralpha) / r 2442 alsq2 = 2.0d0 * aewald**2 2443 alsq2n = 0.0d0 2444 if (aewald .gt. 0.0d0) alsq2n = 1.0d0 / (sqrtpi*aewald) 2445 exp2a = exp(-ralpha**2) 2446 do j = 1, 5 2447 bfac = dble(j+j-1) 2448 alsq2n = alsq2 * alsq2n 2449 bn(j) = (bfac*bn(j-1)+alsq2n*exp2a) / r2 2450 end do 2451 do j = 0, 5 2452 bn(j) = f * bn(j) 2453 end do 2454c 2455c intermediates involving moments and distance separation 2456c 2457 dikx = diy*dkz - diz*dky 2458 diky = diz*dkx - dix*dkz 2459 dikz = dix*dky - diy*dkx 2460 dirx = diy*zr - diz*yr 2461 diry = diz*xr - dix*zr 2462 dirz = dix*yr - diy*xr 2463 dkrx = dky*zr - dkz*yr 2464 dkry = dkz*xr - dkx*zr 2465 dkrz = dkx*yr - dky*xr 2466 dri = dix*xr + diy*yr + diz*zr 2467 drk = dkx*xr + dky*yr + dkz*zr 2468 dik = dix*dkx + diy*dky + diz*dkz 2469 qrix = qixx*xr + qixy*yr + qixz*zr 2470 qriy = qixy*xr + qiyy*yr + qiyz*zr 2471 qriz = qixz*xr + qiyz*yr + qizz*zr 2472 qrkx = qkxx*xr + qkxy*yr + qkxz*zr 2473 qrky = qkxy*xr + qkyy*yr + qkyz*zr 2474 qrkz = qkxz*xr + qkyz*yr + qkzz*zr 2475 qrri = qrix*xr + qriy*yr + qriz*zr 2476 qrrk = qrkx*xr + qrky*yr + qrkz*zr 2477 qrrik = qrix*qrkx + qriy*qrky + qriz*qrkz 2478 qik = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 2479 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 2480 qrixr = qriz*yr - qriy*zr 2481 qriyr = qrix*zr - qriz*xr 2482 qrizr = qriy*xr - qrix*yr 2483 qrkxr = qrkz*yr - qrky*zr 2484 qrkyr = qrkx*zr - qrkz*xr 2485 qrkzr = qrky*xr - qrkx*yr 2486 qrrx = qrky*qriz - qrkz*qriy 2487 qrry = qrkz*qrix - qrkx*qriz 2488 qrrz = qrkx*qriy - qrky*qrix 2489 qikrx = qixx*qrkx + qixy*qrky + qixz*qrkz 2490 qikry = qixy*qrkx + qiyy*qrky + qiyz*qrkz 2491 qikrz = qixz*qrkx + qiyz*qrky + qizz*qrkz 2492 qkirx = qkxx*qrix + qkxy*qriy + qkxz*qriz 2493 qkiry = qkxy*qrix + qkyy*qriy + qkyz*qriz 2494 qkirz = qkxz*qrix + qkyz*qriy + qkzz*qriz 2495 qikrxr = qikrz*yr - qikry*zr 2496 qikryr = qikrx*zr - qikrz*xr 2497 qikrzr = qikry*xr - qikrx*yr 2498 qkirxr = qkirz*yr - qkiry*zr 2499 qkiryr = qkirx*zr - qkirz*xr 2500 qkirzr = qkiry*xr - qkirx*yr 2501 diqkx = dix*qkxx + diy*qkxy + diz*qkxz 2502 diqky = dix*qkxy + diy*qkyy + diz*qkyz 2503 diqkz = dix*qkxz + diy*qkyz + diz*qkzz 2504 dkqix = dkx*qixx + dky*qixy + dkz*qixz 2505 dkqiy = dkx*qixy + dky*qiyy + dkz*qiyz 2506 dkqiz = dkx*qixz + dky*qiyz + dkz*qizz 2507 diqrk = dix*qrkx + diy*qrky + diz*qrkz 2508 dkqri = dkx*qrix + dky*qriy + dkz*qriz 2509 diqkxr = diqkz*yr - diqky*zr 2510 diqkyr = diqkx*zr - diqkz*xr 2511 diqkzr = diqky*xr - diqkx*yr 2512 dkqixr = dkqiz*yr - dkqiy*zr 2513 dkqiyr = dkqix*zr - dkqiz*xr 2514 dkqizr = dkqiy*xr - dkqix*yr 2515 dqiqkx = diy*qrkz - diz*qrky + dky*qriz - dkz*qriy 2516 & - 2.0d0*(qixy*qkxz+qiyy*qkyz+qiyz*qkzz 2517 & -qixz*qkxy-qiyz*qkyy-qizz*qkyz) 2518 dqiqky = diz*qrkx - dix*qrkz + dkz*qrix - dkx*qriz 2519 & - 2.0d0*(qixz*qkxx+qiyz*qkxy+qizz*qkxz 2520 & -qixx*qkxz-qixy*qkyz-qixz*qkzz) 2521 dqiqkz = dix*qrky - diy*qrkx + dkx*qriy - dky*qrix 2522 & - 2.0d0*(qixx*qkxy+qixy*qkyy+qixz*qkyz 2523 & -qixy*qkxx-qiyy*qkxy-qiyz*qkxz) 2524c 2525c calculate intermediate terms for multipole energy 2526c 2527 term1 = ci*ck 2528 term2 = ck*dri - ci*drk + dik 2529 term3 = ci*qrrk + ck*qrri - dri*drk 2530 & + 2.0d0*(dkqri-diqrk+qik) 2531 term4 = dri*qrrk - drk*qrri - 4.0d0*qrrik 2532 term5 = qrri*qrrk 2533c 2534c compute the full energy without any Ewald scaling 2535c 2536 efull = term1*rr1 + term2*rr3 + term3*rr5 2537 & + term4*rr7 + term5*rr9 2538 efull = efull * mscale(kk) 2539 if (molcule(ii) .ne. molcule(kk)) 2540 & einter = einter + efull 2541c 2542c modify distances to account for Ewald and exclusions 2543c 2544 scalekk = 1.0d0 - mscale(kk) 2545 rr1 = bn(0) - scalekk*rr1 2546 rr3 = bn(1) - scalekk*rr3 2547 rr5 = bn(2) - scalekk*rr5 2548 rr7 = bn(3) - scalekk*rr7 2549 rr9 = bn(4) - scalekk*rr9 2550 rr11 = bn(5) - scalekk*rr11 2551c 2552c compute the energy contributions for this interaction 2553c 2554 e = term1*rr1 + term2*rr3 + term3*rr5 2555 & + term4*rr7 + term5*rr9 2556 em = em + e 2557c 2558c calculate intermediate terms for force and torque 2559c 2560 de = term1*rr3 + term2*rr5 + term3*rr7 2561 & + term4*rr9 + term5*rr11 2562 term1 = -ck*rr3 + drk*rr5 - qrrk*rr7 2563 term2 = ci*rr3 + dri*rr5 + qrri*rr7 2564 term3 = 2.0d0 * rr5 2565 term4 = 2.0d0 * (-ck*rr5+drk*rr7-qrrk*rr9) 2566 term5 = 2.0d0 * (-ci*rr5-dri*rr7-qrri*rr9) 2567 term6 = 4.0d0 * rr7 2568c 2569c compute the force components for this interaction 2570c 2571 frcx = de*xr + term1*dix + term2*dkx 2572 & + term3*(diqkx-dkqix) + term4*qrix 2573 & + term5*qrkx + term6*(qikrx+qkirx) 2574 frcy = de*yr + term1*diy + term2*dky 2575 & + term3*(diqky-dkqiy) + term4*qriy 2576 & + term5*qrky + term6*(qikry+qkiry) 2577 frcz = de*zr + term1*diz + term2*dkz 2578 & + term3*(diqkz-dkqiz) + term4*qriz 2579 & + term5*qrkz + term6*(qikrz+qkirz) 2580c 2581c compute the torque components for this interaction 2582c 2583 ttmi(1) = -rr3*dikx + term1*dirx + term3*(dqiqkx+dkqixr) 2584 & - term4*qrixr - term6*(qikrxr+qrrx) 2585 ttmi(2) = -rr3*diky + term1*diry + term3*(dqiqky+dkqiyr) 2586 & - term4*qriyr - term6*(qikryr+qrry) 2587 ttmi(3) = -rr3*dikz + term1*dirz + term3*(dqiqkz+dkqizr) 2588 & - term4*qrizr - term6*(qikrzr+qrrz) 2589 ttmk(1) = rr3*dikx + term2*dkrx - term3*(dqiqkx+diqkxr) 2590 & - term5*qrkxr - term6*(qkirxr-qrrx) 2591 ttmk(2) = rr3*diky + term2*dkry - term3*(dqiqky+diqkyr) 2592 & - term5*qrkyr - term6*(qkiryr-qrry) 2593 ttmk(3) = rr3*dikz + term2*dkrz - term3*(dqiqkz+diqkzr) 2594 & - term5*qrkzr - term6*(qkirzr-qrrz) 2595c 2596c increment force-based gradient and torque on first site 2597c 2598 dem(1,ii) = dem(1,ii) + frcx 2599 dem(2,ii) = dem(2,ii) + frcy 2600 dem(3,ii) = dem(3,ii) + frcz 2601 tem(1,i) = tem(1,i) + ttmi(1) 2602 tem(2,i) = tem(2,i) + ttmi(2) 2603 tem(3,i) = tem(3,i) + ttmi(3) 2604c 2605c increment force-based gradient and torque on second site 2606c 2607 dem(1,kk) = dem(1,kk) - frcx 2608 dem(2,kk) = dem(2,kk) - frcy 2609 dem(3,kk) = dem(3,kk) - frcz 2610 tem(1,k) = tem(1,k) + ttmk(1) 2611 tem(2,k) = tem(2,k) + ttmk(2) 2612 tem(3,k) = tem(3,k) + ttmk(3) 2613c 2614c increment the virial due to pairwise Cartesian forces 2615c 2616 vxx = -xr * frcx 2617 vxy = -0.5d0 * (yr*frcx + xr*frcy) 2618 vxz = -0.5d0 * (zr*frcx + xr*frcz) 2619 vyy = -yr * frcy 2620 vyz = -0.5d0 * (zr*frcy + yr*frcz) 2621 vzz = -zr * frcz 2622 vir(1,1) = vir(1,1) + vxx 2623 vir(2,1) = vir(2,1) + vxy 2624 vir(3,1) = vir(3,1) + vxz 2625 vir(1,2) = vir(1,2) + vxy 2626 vir(2,2) = vir(2,2) + vyy 2627 vir(3,2) = vir(3,2) + vyz 2628 vir(1,3) = vir(1,3) + vxz 2629 vir(2,3) = vir(2,3) + vyz 2630 vir(3,3) = vir(3,3) + vzz 2631 end if 2632 end do 2633c 2634c reset exclusion coefficients for connected atoms 2635c 2636 do j = 1, n12(ii) 2637 mscale(i12(j,ii)) = 1.0d0 2638 end do 2639 do j = 1, n13(ii) 2640 mscale(i13(j,ii)) = 1.0d0 2641 end do 2642 do j = 1, n14(ii) 2643 mscale(i14(j,ii)) = 1.0d0 2644 end do 2645 do j = 1, n15(ii) 2646 mscale(i15(j,ii)) = 1.0d0 2647 end do 2648 end do 2649c 2650c OpenMP directives for the major loop structure 2651c 2652!$OMP END DO 2653!$OMP DO reduction(+:dem,vir) schedule(guided) 2654c 2655c resolve site torques then increment forces and virial 2656c 2657 do i = 1, npole 2658 call torque (i,tem(1,i),fix,fiy,fiz,dem) 2659 ii = ipole(i) 2660 iaz = zaxis(i) 2661 iax = xaxis(i) 2662 iay = yaxis(i) 2663 if (iaz .eq. 0) iaz = ii 2664 if (iax .eq. 0) iax = ii 2665 if (iay .eq. 0) iay = ii 2666 xiz = x(iaz) - x(ii) 2667 yiz = y(iaz) - y(ii) 2668 ziz = z(iaz) - z(ii) 2669 xix = x(iax) - x(ii) 2670 yix = y(iax) - y(ii) 2671 zix = z(iax) - z(ii) 2672 xiy = x(iay) - x(ii) 2673 yiy = y(iay) - y(ii) 2674 ziy = z(iay) - z(ii) 2675 2676c vxx = xix*fix(1) + xiy*fiy(1) + xiz*fiz(1) 2677c vxy = 0.5d0 * (yix*fix(1) + yiy*fiy(1) + yiz*fiz(1) 2678c & + xix*fix(2) + xiy*fiy(2) + xiz*fiz(2)) 2679c vxz = 0.5d0 * (zix*fix(1) + ziy*fiy(1) + ziz*fiz(1) 2680c & + xix*fix(3) + xiy*fiy(3) + xiz*fiz(3)) 2681c vyy = yix*fix(2) + yiy*fiy(2) + yiz*fiz(2) 2682c vyz = 0.5d0 * (zix*fix(2) + ziy*fiy(2) + ziz*fiz(2) 2683c & + yix*fix(3) + yiy*fiy(3) + yiz*fiz(3)) 2684c vzz = zix*fix(3) + ziy*fiy(3) + ziz*fiz(3) 2685 2686 vxx = xix*fix(1) + xiy*fiy(1) + xiz*fiz(1) 2687c vxy = 0.5d0*(yix*fix(1) + yiy*fiy(1) + yiz*fiz(1) 2688c & + xix*fix(2) + xiy*fiy(2) + xiz*fiz(2)) 2689 vxy = 0.5d0*(xiy*fix(1) + yiy*fiy(1) + ziy*fiz(1) 2690 & + xix*fix(2) + yix*fiy(2) + zix*fiz(2)) 2691c vxz = 0.5d0*(zix*fix(1) + ziy*fiy(1) + ziz*fiz(1) 2692c & + xix*fix(3) + xiy*fiy(3) + xiz*fiz(3)) 2693 vxz = 0.5d0*(xiz*fix(1) + yiz*fiy(1) + ziz*fiz(1) 2694 & + xix*fix(3) + yix*fiy(3) + zix*fiz(3)) 2695 vyy = yix*fix(2) + yiy*fiy(2) + yiz*fiz(2) 2696c vyz = 0.5d0*(zix*fix(2) + ziy*fiy(2) + ziz*fiz(2) 2697c & + yix*fix(3) + yiy*fiy(3) + yiz*fiz(3)) 2698 vyz = 0.5d0*(xiz*fix(2) + yiz*fiy(2) + ziz*fiz(2) 2699 & + xiy*fix(3) + yiy*fiy(3) + ziy*fiz(3)) 2700 vzz = zix*fix(3) + ziy*fiy(3) + ziz*fiz(3) 2701 2702 vir(1,1) = vir(1,1) + vxx 2703 vir(2,1) = vir(2,1) + vxy 2704 vir(3,1) = vir(3,1) + vxz 2705 vir(1,2) = vir(1,2) + vxy 2706 vir(2,2) = vir(2,2) + vyy 2707 vir(3,2) = vir(3,2) + vyz 2708 vir(1,3) = vir(1,3) + vxz 2709 vir(2,3) = vir(2,3) + vyz 2710 vir(3,3) = vir(3,3) + vzz 2711 end do 2712c 2713c OpenMP directives for the major loop structure 2714c 2715!$OMP END DO 2716!$OMP END PARALLEL 2717c 2718c perform deallocation of some local arrays 2719c 2720 deallocate (mscale) 2721 deallocate (tem) 2722 return 2723 end 2724c 2725c 2726c #################################################################### 2727c ## ## 2728c ## subroutine emrecip1 -- PME recip multipole energy & derivs ## 2729c ## ## 2730c #################################################################### 2731c 2732c 2733c "emrecip1" evaluates the reciprocal space portion of the particle 2734c mesh Ewald summation energy and gradient due to multipoles 2735c 2736c literature reference: 2737c 2738c C. Sagui, L. G. Pedersen and T. A. Darden, "Towards an Accurate 2739c Representation of Electrostatics in Classical Force Fields: 2740c Efficient Implementation of Multipolar Interactions in 2741c Biomolecular Simulations", Journal of Chemical Physics, 120, 2742c 73-87 (2004) 2743c 2744c modifications for nonperiodic systems suggested by Tom Darden 2745c during May 2007 2746c 2747c 2748 subroutine emrecip1 2749 use sizes 2750 use atoms 2751 use bound 2752 use boxes 2753 use chgpot 2754 use deriv 2755 use energi 2756 use ewald 2757 use math 2758 use mpole 2759 use mrecip 2760 use pme 2761 use virial 2762 implicit none 2763 integer i,j,k,ii 2764 integer k1,k2,k3 2765 integer m1,m2,m3 2766 integer iax,iay,iaz 2767 integer ntot,nff 2768 integer nf1,nf2,nf3 2769 integer deriv1(10) 2770 integer deriv2(10) 2771 integer deriv3(10) 2772 real*8 e,eterm,f 2773 real*8 r1,r2,r3 2774 real*8 h1,h2,h3 2775 real*8 f1,f2,f3 2776 real*8 xix,yix,zix 2777 real*8 xiy,yiy,ziy 2778 real*8 xiz,yiz,ziz 2779 real*8 vxx,vyy,vzz 2780 real*8 vxy,vxz,vyz 2781 real*8 volterm,denom 2782 real*8 hsq,expterm 2783 real*8 term,pterm 2784 real*8 vterm,struc2 2785 real*8 trq(3),fix(3) 2786 real*8 fiy(3),fiz(3) 2787c 2788c indices into the electrostatic field array 2789c 2790 data deriv1 / 2, 5, 8, 9, 11, 16, 18, 14, 15, 20 / 2791 data deriv2 / 3, 8, 6, 10, 14, 12, 19, 16, 20, 17 / 2792 data deriv3 / 4, 9, 10, 7, 15, 17, 13, 20, 18, 19 / 2793c 2794c 2795c return if the Ewald coefficient is zero 2796c 2797 if (aewald .lt. 1.0d-6) return 2798 f = electric / dielec 2799c 2800c perform dynamic allocation of some global arrays 2801c 2802 if (allocated(cmp)) then 2803 if (size(cmp) .lt. 10*npole) then 2804 deallocate (cmp) 2805 deallocate (fmp) 2806 deallocate (cphi) 2807 deallocate (fphi) 2808 end if 2809 end if 2810 if (.not. allocated(cmp)) then 2811 allocate (cmp(10,npole)) 2812 allocate (fmp(10,npole)) 2813 allocate (cphi(10,npole)) 2814 allocate (fphi(20,npole)) 2815 end if 2816c 2817c zero out the temporary virial accumulation variables 2818c 2819 vxx = 0.0d0 2820 vxy = 0.0d0 2821 vxz = 0.0d0 2822 vyy = 0.0d0 2823 vyz = 0.0d0 2824 vzz = 0.0d0 2825c 2826c copy multipole moments and coordinates to local storage 2827c 2828 do i = 1, npole 2829 cmp(1,i) = rpole(1,i) 2830 cmp(2,i) = rpole(2,i) 2831 cmp(3,i) = rpole(3,i) 2832 cmp(4,i) = rpole(4,i) 2833 cmp(5,i) = rpole(5,i) 2834 cmp(6,i) = rpole(9,i) 2835 cmp(7,i) = rpole(13,i) 2836 cmp(8,i) = 2.0d0 * rpole(6,i) 2837 cmp(9,i) = 2.0d0 * rpole(7,i) 2838 cmp(10,i) = 2.0d0 * rpole(10,i) 2839 end do 2840c 2841c compute the arrays of B-spline coefficients 2842c 2843 call bspline_fill 2844 call table_fill 2845c 2846c assign permanent multipoles to PME grid and perform 2847c the 3-D FFT forward transformation 2848c 2849 call cmp_to_fmp (cmp,fmp) 2850 call grid_mpole (fmp) 2851 call fftfront 2852c 2853c initialize variables required for the scalar summation 2854c 2855 ntot = nfft1 * nfft2 * nfft3 2856 pterm = (pi/aewald)**2 2857 volterm = pi * volbox 2858 nff = nfft1 * nfft2 2859 nf1 = (nfft1+1) / 2 2860 nf2 = (nfft2+1) / 2 2861 nf3 = (nfft3+1) / 2 2862c 2863c make the scalar summation over reciprocal lattice 2864c 2865 do i = 1, ntot-1 2866 k3 = i/nff + 1 2867 j = i - (k3-1)*nff 2868 k2 = j/nfft1 + 1 2869 k1 = j - (k2-1)*nfft1 + 1 2870 m1 = k1 - 1 2871 m2 = k2 - 1 2872 m3 = k3 - 1 2873 if (k1 .gt. nf1) m1 = m1 - nfft1 2874 if (k2 .gt. nf2) m2 = m2 - nfft2 2875 if (k3 .gt. nf3) m3 = m3 - nfft3 2876 r1 = dble(m1) 2877 r2 = dble(m2) 2878 r3 = dble(m3) 2879 h1 = recip(1,1)*r1 + recip(1,2)*r2 + recip(1,3)*r3 2880 h2 = recip(2,1)*r1 + recip(2,2)*r2 + recip(2,3)*r3 2881 h3 = recip(3,1)*r1 + recip(3,2)*r2 + recip(3,3)*r3 2882 hsq = h1*h1 + h2*h2 + h3*h3 2883 term = -pterm * hsq 2884 expterm = 0.0d0 2885 if (term .gt. -50.0d0) then 2886 denom = volterm*hsq*bsmod1(k1)*bsmod2(k2)*bsmod3(k3) 2887 expterm = exp(term) / denom 2888 if (.not. use_bounds) then 2889 expterm = expterm * (1.0d0-cos(pi*xbox*sqrt(hsq))) 2890 else if (octahedron) then 2891 if (mod(m1+m2+m3,2) .ne. 0) expterm = 0.0d0 2892 end if 2893 struc2 = qgrid(1,k1,k2,k3)**2 + qgrid(2,k1,k2,k3)**2 2894 eterm = 0.5d0 * electric * expterm * struc2 2895 vterm = (2.0d0/hsq) * (1.0d0-term) * eterm 2896 vxx = vxx + h1*h1*vterm - eterm 2897 vxy = vxy + h1*h2*vterm 2898 vxz = vxz + h1*h3*vterm 2899 vyy = vyy + h2*h2*vterm - eterm 2900 vyz = vyz + h2*h3*vterm 2901 vzz = vzz + h3*h3*vterm - eterm 2902 end if 2903 qfac(k1,k2,k3) = expterm 2904 end do 2905c 2906c save the virial for use in polarization computation 2907c 2908 vmxx = vxx 2909 vmxy = vxy 2910 vmxz = vxz 2911 vmyy = vyy 2912 vmyz = vyz 2913 vmzz = vzz 2914c 2915c account for zeroth grid point for nonperiodic system 2916c 2917 qfac(1,1,1) = 0.0d0 2918 if (.not. use_bounds) then 2919 expterm = 0.5d0 * pi / xbox 2920 struc2 = qgrid(1,1,1,1)**2 + qgrid(2,1,1,1)**2 2921 e = f * expterm * struc2 2922 em = em + e 2923 qfac(1,1,1) = expterm 2924 end if 2925c 2926c complete the transformation of the PME grid 2927c 2928 do k = 1, nfft3 2929 do j = 1, nfft2 2930 do i = 1, nfft1 2931 term = qfac(i,j,k) 2932 qgrid(1,i,j,k) = term * qgrid(1,i,j,k) 2933 qgrid(2,i,j,k) = term * qgrid(2,i,j,k) 2934 end do 2935 end do 2936 end do 2937c 2938c perform 3-D FFT backward transform and get potential 2939c 2940 call fftback 2941 call fphi_mpole (fphi) 2942 do i = 1, npole 2943 do j = 1, 20 2944 fphi(j,i) = f * fphi(j,i) 2945 end do 2946 end do 2947 call fphi_to_cphi (fphi,cphi) 2948c 2949c increment the permanent multipole energy and gradient 2950c 2951 e = 0.0d0 2952 do i = 1, npole 2953 f1 = 0.0d0 2954 f2 = 0.0d0 2955 f3 = 0.0d0 2956 do k = 1, 10 2957 e = e + fmp(k,i)*fphi(k,i) 2958 f1 = f1 + fmp(k,i)*fphi(deriv1(k),i) 2959 f2 = f2 + fmp(k,i)*fphi(deriv2(k),i) 2960 f3 = f3 + fmp(k,i)*fphi(deriv3(k),i) 2961 end do 2962 f1 = dble(nfft1) * f1 2963 f2 = dble(nfft2) * f2 2964 f3 = dble(nfft3) * f3 2965 h1 = recip(1,1)*f1 + recip(1,2)*f2 + recip(1,3)*f3 2966 h2 = recip(2,1)*f1 + recip(2,2)*f2 + recip(2,3)*f3 2967 h3 = recip(3,1)*f1 + recip(3,2)*f2 + recip(3,3)*f3 2968 ii = ipole(i) 2969 dem(1,ii) = dem(1,ii) + h1 2970 dem(2,ii) = dem(2,ii) + h2 2971 dem(3,ii) = dem(3,ii) + h3 2972 end do 2973 e = 0.5d0 * e 2974 em = em + e 2975c 2976c increment the permanent multipole virial contributions 2977c 2978 do i = 1, npole 2979 vxx = vxx - cmp(2,i)*cphi(2,i) - 2.0d0*cmp(5,i)*cphi(5,i) 2980 & - cmp(8,i)*cphi(8,i) - cmp(9,i)*cphi(9,i) 2981 vxy = vxy - 0.5d0*(cmp(3,i)*cphi(2,i)+cmp(2,i)*cphi(3,i)) 2982 & - (cmp(5,i)+cmp(6,i))*cphi(8,i) 2983 & - 0.5d0*cmp(8,i)*(cphi(5,i)+cphi(6,i)) 2984 & - 0.5d0*(cmp(9,i)*cphi(10,i)+cmp(10,i)*cphi(9,i)) 2985 vxz = vxz - 0.5d0*(cmp(4,i)*cphi(2,i)+cmp(2,i)*cphi(4,i)) 2986 & - (cmp(5,i)+cmp(7,i))*cphi(9,i) 2987 & - 0.5d0*cmp(9,i)*(cphi(5,i)+cphi(7,i)) 2988 & - 0.5d0*(cmp(8,i)*cphi(10,i)+cmp(10,i)*cphi(8,i)) 2989 vyy = vyy - cmp(3,i)*cphi(3,i) - 2.0d0*cmp(6,i)*cphi(6,i) 2990 & - cmp(8,i)*cphi(8,i) - cmp(10,i)*cphi(10,i) 2991 vyz = vyz - 0.5d0*(cmp(4,i)*cphi(3,i)+cmp(3,i)*cphi(4,i)) 2992 & - (cmp(6,i)+cmp(7,i))*cphi(10,i) 2993 & - 0.5d0*cmp(10,i)*(cphi(6,i)+cphi(7,i)) 2994 & - 0.5d0*(cmp(8,i)*cphi(9,i)+cmp(9,i)*cphi(8,i)) 2995 vzz = vzz - cmp(4,i)*cphi(4,i) - 2.0d0*cmp(7,i)*cphi(7,i) 2996 & - cmp(9,i)*cphi(9,i) - cmp(10,i)*cphi(10,i) 2997 end do 2998c 2999c resolve site torques then increment forces and virial 3000c 3001 do i = 1, npole 3002 trq(1) = cmp(4,i)*cphi(3,i) - cmp(3,i)*cphi(4,i) 3003 & + 2.0d0*(cmp(7,i)-cmp(6,i))*cphi(10,i) 3004 & + cmp(9,i)*cphi(8,i) + cmp(10,i)*cphi(6,i) 3005 & - cmp(8,i)*cphi(9,i) - cmp(10,i)*cphi(7,i) 3006 trq(2) = cmp(2,i)*cphi(4,i) - cmp(4,i)*cphi(2,i) 3007 & + 2.0d0*(cmp(5,i)-cmp(7,i))*cphi(9,i) 3008 & + cmp(8,i)*cphi(10,i) + cmp(9,i)*cphi(7,i) 3009 & - cmp(9,i)*cphi(5,i) - cmp(10,i)*cphi(8,i) 3010 trq(3) = cmp(3,i)*cphi(2,i) - cmp(2,i)*cphi(3,i) 3011 & + 2.0d0*(cmp(6,i)-cmp(5,i))*cphi(8,i) 3012 & + cmp(8,i)*cphi(5,i) + cmp(10,i)*cphi(9,i) 3013 & - cmp(8,i)*cphi(6,i) - cmp(9,i)*cphi(10,i) 3014 call torque (i,trq,fix,fiy,fiz,dem) 3015 ii = ipole(i) 3016 iaz = zaxis(i) 3017 iax = xaxis(i) 3018 iay = yaxis(i) 3019 if (iaz .eq. 0) iaz = ii 3020 if (iax .eq. 0) iax = ii 3021 if (iay .eq. 0) iay = ii 3022 xiz = x(iaz) - x(ii) 3023 yiz = y(iaz) - y(ii) 3024 ziz = z(iaz) - z(ii) 3025 xix = x(iax) - x(ii) 3026 yix = y(iax) - y(ii) 3027 zix = z(iax) - z(ii) 3028 xiy = x(iay) - x(ii) 3029 yiy = y(iay) - y(ii) 3030 ziy = z(iay) - z(ii) 3031 vxx = vxx + xix*fix(1) + xiy*fiy(1) + xiz*fiz(1) 3032c vxy = vxy + 0.5d0*(yix*fix(1) + yiy*fiy(1) + yiz*fiz(1) 3033c & + xix*fix(2) + xiy*fiy(2) + xiz*fiz(2)) 3034 vxy = vxy + 0.5d0*(xiy*fix(1) + yiy*fiy(1) + ziy*fiz(1) 3035 & + xix*fix(2) + yix*fiy(2) + zix*fiz(2)) 3036c vxz = vxz + 0.5d0*(zix*fix(1) + ziy*fiy(1) + ziz*fiz(1) 3037c & + xix*fix(3) + xiy*fiy(3) + xiz*fiz(3)) 3038 vxz = vxz + 0.5d0*(xiz*fix(1) + yiz*fiy(1) + ziz*fiz(1) 3039 & + xix*fix(3) + yix*fiy(3) + zix*fiz(3)) 3040 vyy = vyy + yix*fix(2) + yiy*fiy(2) + yiz*fiz(2) 3041c vyz = vyz + 0.5d0*(zix*fix(2) + ziy*fiy(2) + ziz*fiz(2) 3042c & + yix*fix(3) + yiy*fiy(3) + yiz*fiz(3)) 3043 vyz = vyz + 0.5d0*(xiz*fix(2) + yiz*fiy(2) + ziz*fiz(2) 3044 & + xiy*fix(3) + yiy*fiy(3) + ziy*fiz(3)) 3045 vzz = vzz + zix*fix(3) + ziy*fiy(3) + ziz*fiz(3) 3046 end do 3047c 3048c increment the total internal virial tensor components 3049c 3050 vir(1,1) = vir(1,1) + vxx 3051 vir(2,1) = vir(2,1) + vxy 3052 vir(3,1) = vir(3,1) + vxz 3053 vir(1,2) = vir(1,2) + vxy 3054 vir(2,2) = vir(2,2) + vyy 3055 vir(3,2) = vir(3,2) + vyz 3056 vir(1,3) = vir(1,3) + vxz 3057 vir(2,3) = vir(2,3) + vyz 3058 vir(3,3) = vir(3,3) + vzz 3059 return 3060 end 3061