1c 2c 3c ############################################################# 4c ## COPYRIGHT (C) 1999 by Pengyu Ren & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ############################################################# 7c 8c ################################################################# 9c ## ## 10c ## subroutine empole2 -- atomic multipole Hessian elements ## 11c ## ## 12c ################################################################# 13c 14c 15c "empole2" calculates second derivatives of the multipole energy 16c for a single atom at a time 17c 18c 19 subroutine empole2 (i) 20 use atoms 21 use deriv 22 use energi 23 use hessn 24 use mpole 25 use potent 26 implicit none 27 integer i,j,k 28 integer nlist 29 integer, allocatable :: list(:) 30 real*8 eps,old 31 real*8, allocatable :: d0(:,:) 32 logical prior 33 logical twosided 34c 35c 36c set the default stepsize and accuracy control flags 37c 38 eps = 1.0d-5 39 twosided = .false. 40 if (n .le. 300) twosided = .true. 41c 42c perform dynamic allocation of some local arrays 43c 44 allocate (list(npole)) 45 allocate (d0(3,n)) 46c 47c perform dynamic allocation of some global arrays 48c 49 prior = .false. 50 if (allocated(dem)) then 51 prior = .true. 52 if (size(dem) .lt. 3*n) deallocate (dem) 53 end if 54 if (.not. allocated(dem)) allocate (dem(3,n)) 55c 56c find the multipole definitions involving the current atom 57c 58 nlist = 0 59 do k = 1, npole 60 if (ipole(k).eq.i .or. zaxis(k).eq.i 61 & .or. xaxis(k).eq.i .or. abs(yaxis(k)).eq.i) then 62 nlist = nlist + 1 63 list(nlist) = k 64 end if 65 end do 66c 67c get multipole first derivatives for the base structure 68c 69 if (.not. twosided) then 70 call empole2a (nlist,list) 71 do k = 1, n 72 do j = 1, 3 73 d0(j,k) = dem(j,k) 74 end do 75 end do 76 end if 77c 78c find numerical x-components via perturbed structures 79c 80 old = x(i) 81 if (twosided) then 82 x(i) = x(i) - 0.5d0*eps 83 call empole2a (nlist,list) 84 do k = 1, n 85 do j = 1, 3 86 d0(j,k) = dem(j,k) 87 end do 88 end do 89 end if 90 x(i) = x(i) + eps 91 call empole2a (nlist,list) 92 x(i) = old 93 do k = 1, n 94 do j = 1, 3 95 hessx(j,k) = hessx(j,k) + (dem(j,k)-d0(j,k))/eps 96 end do 97 end do 98c 99c find numerical y-components via perturbed structures 100c 101 old = y(i) 102 if (twosided) then 103 y(i) = y(i) - 0.5d0*eps 104 call empole2a (nlist,list) 105 do k = 1, n 106 do j = 1, 3 107 d0(j,k) = dem(j,k) 108 end do 109 end do 110 end if 111 y(i) = y(i) + eps 112 call empole2a (nlist,list) 113 y(i) = old 114 do k = 1, n 115 do j = 1, 3 116 hessy(j,k) = hessy(j,k) + (dem(j,k)-d0(j,k))/eps 117 end do 118 end do 119c 120c find numerical z-components via perturbed structures 121c 122 old = z(i) 123 if (twosided) then 124 z(i) = z(i) - 0.5d0*eps 125 call empole2a (nlist,list) 126 do k = 1, n 127 do j = 1, 3 128 d0(j,k) = dem(j,k) 129 end do 130 end do 131 end if 132 z(i) = z(i) + eps 133 call empole2a (nlist,list) 134 z(i) = old 135 do k = 1, n 136 do j = 1, 3 137 hessz(j,k) = hessz(j,k) + (dem(j,k)-d0(j,k))/eps 138 end do 139 end do 140c 141c perform deallocation of some global arrays 142c 143 if (.not. prior) deallocate (dem) 144c 145c perform deallocation of some local arrays 146c 147 deallocate (list) 148 deallocate (d0) 149 return 150 end 151c 152c 153c ############################################################## 154c ## ## 155c ## subroutine empole2a -- multipole derivatives utility ## 156c ## ## 157c ############################################################## 158c 159c 160c "empole2a" computes multipole first derivatives for a single 161c atom; used to get finite difference second derivatives 162c 163c 164 subroutine empole2a (nlist,list) 165 use atoms 166 use bound 167 use boxes 168 use cell 169 use chgpen 170 use chgpot 171 use couple 172 use deriv 173 use group 174 use limits 175 use molcul 176 use mplpot 177 use mpole 178 use potent 179 use shunt 180 use usage 181 implicit none 182 integer i,j,k 183 integer ii,kk,iii 184 integer nlist,jcell 185 integer ix,iy,iz 186 integer kx,ky,kz 187 integer list(*) 188 real*8 de,f,fgrp 189 real*8 xi,yi,zi 190 real*8 xr,yr,zr 191 real*8 r,r2,rr1,rr3 192 real*8 rr5,rr7,rr9,rr11 193 real*8 rr3i,rr5i,rr7i 194 real*8 rr3k,rr5k,rr7k 195 real*8 rr3ik,rr5ik,rr7ik 196 real*8 rr9ik,rr11ik 197 real*8 ci,dix,diy,diz 198 real*8 qixx,qixy,qixz 199 real*8 qiyy,qiyz,qizz 200 real*8 ck,dkx,dky,dkz 201 real*8 qkxx,qkxy,qkxz 202 real*8 qkyy,qkyz,qkzz 203 real*8 dir,dkr,dik,qik 204 real*8 qix,qiy,qiz,qir 205 real*8 qkx,qky,qkz,qkr 206 real*8 diqk,dkqi,qiqk 207 real*8 dirx,diry,dirz 208 real*8 dkrx,dkry,dkrz 209 real*8 dikx,diky,dikz 210 real*8 qirx,qiry,qirz 211 real*8 qkrx,qkry,qkrz 212 real*8 qikx,qiky,qikz 213 real*8 qixk,qiyk,qizk 214 real*8 qkxi,qkyi,qkzi 215 real*8 qikrx,qikry,qikrz 216 real*8 qkirx,qkiry,qkirz 217 real*8 diqkx,diqky,diqkz 218 real*8 dkqix,dkqiy,dkqiz 219 real*8 diqkrx,diqkry,diqkrz 220 real*8 dkqirx,dkqiry,dkqirz 221 real*8 dqikx,dqiky,dqikz 222 real*8 corei,corek 223 real*8 vali,valk 224 real*8 alphai,alphak 225 real*8 term1,term2,term3 226 real*8 term4,term5,term6 227 real*8 term1i,term2i,term3i 228 real*8 term1k,term2k,term3k 229 real*8 term1ik,term2ik,term3ik 230 real*8 term4ik,term5ik 231 real*8 poti,potk 232 real*8 frcx,frcy,frcz 233 real*8 ttmi(3),ttmk(3) 234 real*8 fix(3),fiy(3),fiz(3) 235 real*8 dmpi(9),dmpk(9) 236 real*8 dmpik(11) 237 real*8, allocatable :: mscale(:) 238 real*8, allocatable :: tem(:,:) 239 real*8, allocatable :: pot(:) 240 real*8, allocatable :: decfx(:) 241 real*8, allocatable :: decfy(:) 242 real*8, allocatable :: decfz(:) 243 logical proceed 244 logical usei,usek 245 character*6 mode 246c 247c 248c zero out the multipole first derivative components 249c 250 do i = 1, n 251 do j = 1, 3 252 dem(j,i) = 0.0d0 253 end do 254 end do 255 if (nlist .eq. 0) return 256c 257c perform dynamic allocation of some local arrays 258c 259 allocate (mscale(n)) 260 allocate (tem(3,n)) 261 allocate (pot(n)) 262 allocate (decfx(n)) 263 allocate (decfy(n)) 264 allocate (decfz(n)) 265c 266c initialize scaling, torque and potential arrays 267c 268 do i = 1, n 269 mscale(i) = 1.0d0 270 do j = 1, 3 271 tem(j,i) = 0.0d0 272 end do 273 pot(i) = 0.0d0 274 end do 275c 276c set conversion factor, cutoff and switching coefficients 277c 278 f = electric / dielec 279 mode = 'MPOLE' 280 call switch (mode) 281c 282c alter partial charges and multipoles for charge flux 283c 284 if (use_chgflx) call alterchg 285c 286c check the sign of multipole components at chiral sites 287c 288 call chkpole 289c 290c rotate the multipole components into the global frame 291c 292 call rotpole 293c 294c compute components of the multipole interaction gradient 295c 296 do iii = 1, nlist 297 ii = list(iii) 298 i = ipole(ii) 299 iz = zaxis(ii) 300 ix = xaxis(ii) 301 iy = abs(yaxis(ii)) 302 xi = x(i) 303 yi = y(i) 304 zi = z(i) 305 ci = rpole(1,ii) 306 dix = rpole(2,ii) 307 diy = rpole(3,ii) 308 diz = rpole(4,ii) 309 qixx = rpole(5,ii) 310 qixy = rpole(6,ii) 311 qixz = rpole(7,ii) 312 qiyy = rpole(9,ii) 313 qiyz = rpole(10,ii) 314 qizz = rpole(13,ii) 315 if (use_chgpen) then 316 corei = pcore(ii) 317 vali = pval(ii) 318 alphai = palpha(ii) 319 end if 320 usei = (use(i) .or. use(iz) .or. use(ix) .or. use(iy)) 321c 322c set exclusion coefficients for connected atoms 323c 324 do j = 1, n12(i) 325 mscale(i12(j,i)) = m2scale 326 end do 327 do j = 1, n13(i) 328 mscale(i13(j,i)) = m3scale 329 end do 330 do j = 1, n14(i) 331 mscale(i14(j,i)) = m4scale 332 end do 333 do j = 1, n15(i) 334 mscale(i15(j,i)) = m5scale 335 end do 336c 337c evaluate all sites within the cutoff distance 338c 339 do kk = 1, npole 340 if (kk .eq. ii) goto 10 341 k = ipole(kk) 342 kz = zaxis(kk) 343 kx = xaxis(kk) 344 ky = abs(yaxis(kk)) 345 usek = (use(k) .or. use(kz) .or. use(kx) .or. use(ky)) 346 proceed = .true. 347 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 348 if (.not. use_intra) proceed = .true. 349 if (proceed) proceed = (usei .or. usek) 350 if (.not. proceed) goto 10 351 xr = x(k) - xi 352 yr = y(k) - yi 353 zr = z(k) - zi 354 if (use_bounds) call image (xr,yr,zr) 355 r2 = xr*xr + yr*yr + zr*zr 356 if (r2 .le. off2) then 357 r = sqrt(r2) 358 ck = rpole(1,kk) 359 dkx = rpole(2,kk) 360 dky = rpole(3,kk) 361 dkz = rpole(4,kk) 362 qkxx = rpole(5,kk) 363 qkxy = rpole(6,kk) 364 qkxz = rpole(7,kk) 365 qkyy = rpole(9,kk) 366 qkyz = rpole(10,kk) 367 qkzz = rpole(13,kk) 368c 369c intermediates involving moments and separation distance 370c 371 dir = dix*xr + diy*yr + diz*zr 372 qix = qixx*xr + qixy*yr + qixz*zr 373 qiy = qixy*xr + qiyy*yr + qiyz*zr 374 qiz = qixz*xr + qiyz*yr + qizz*zr 375 qir = qix*xr + qiy*yr + qiz*zr 376 dkr = dkx*xr + dky*yr + dkz*zr 377 qkx = qkxx*xr + qkxy*yr + qkxz*zr 378 qky = qkxy*xr + qkyy*yr + qkyz*zr 379 qkz = qkxz*xr + qkyz*yr + qkzz*zr 380 qkr = qkx*xr + qky*yr + qkz*zr 381 dik = dix*dkx + diy*dky + diz*dkz 382 qik = qix*qkx + qiy*qky + qiz*qkz 383 diqk = dix*qkx + diy*qky + diz*qkz 384 dkqi = dkx*qix + dky*qiy + dkz*qiz 385 qiqk = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 386 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 387c 388c additional intermediates involving moments and distance 389c 390 dirx = diy*zr - diz*yr 391 diry = diz*xr - dix*zr 392 dirz = dix*yr - diy*xr 393 dkrx = dky*zr - dkz*yr 394 dkry = dkz*xr - dkx*zr 395 dkrz = dkx*yr - dky*xr 396 dikx = diy*dkz - diz*dky 397 diky = diz*dkx - dix*dkz 398 dikz = dix*dky - diy*dkx 399 qirx = qiz*yr - qiy*zr 400 qiry = qix*zr - qiz*xr 401 qirz = qiy*xr - qix*yr 402 qkrx = qkz*yr - qky*zr 403 qkry = qkx*zr - qkz*xr 404 qkrz = qky*xr - qkx*yr 405 qikx = qky*qiz - qkz*qiy 406 qiky = qkz*qix - qkx*qiz 407 qikz = qkx*qiy - qky*qix 408 qixk = qixx*qkx + qixy*qky + qixz*qkz 409 qiyk = qixy*qkx + qiyy*qky + qiyz*qkz 410 qizk = qixz*qkx + qiyz*qky + qizz*qkz 411 qkxi = qkxx*qix + qkxy*qiy + qkxz*qiz 412 qkyi = qkxy*qix + qkyy*qiy + qkyz*qiz 413 qkzi = qkxz*qix + qkyz*qiy + qkzz*qiz 414 qikrx = qizk*yr - qiyk*zr 415 qikry = qixk*zr - qizk*xr 416 qikrz = qiyk*xr - qixk*yr 417 qkirx = qkzi*yr - qkyi*zr 418 qkiry = qkxi*zr - qkzi*xr 419 qkirz = qkyi*xr - qkxi*yr 420 diqkx = dix*qkxx + diy*qkxy + diz*qkxz 421 diqky = dix*qkxy + diy*qkyy + diz*qkyz 422 diqkz = dix*qkxz + diy*qkyz + diz*qkzz 423 dkqix = dkx*qixx + dky*qixy + dkz*qixz 424 dkqiy = dkx*qixy + dky*qiyy + dkz*qiyz 425 dkqiz = dkx*qixz + dky*qiyz + dkz*qizz 426 diqkrx = diqkz*yr - diqky*zr 427 diqkry = diqkx*zr - diqkz*xr 428 diqkrz = diqky*xr - diqkx*yr 429 dkqirx = dkqiz*yr - dkqiy*zr 430 dkqiry = dkqix*zr - dkqiz*xr 431 dkqirz = dkqiy*xr - dkqix*yr 432 dqikx = diy*qkz - diz*qky + dky*qiz - dkz*qiy 433 & - 2.0d0*(qixy*qkxz+qiyy*qkyz+qiyz*qkzz 434 & -qixz*qkxy-qiyz*qkyy-qizz*qkyz) 435 dqiky = diz*qkx - dix*qkz + dkz*qix - dkx*qiz 436 & - 2.0d0*(qixz*qkxx+qiyz*qkxy+qizz*qkxz 437 & -qixx*qkxz-qixy*qkyz-qixz*qkzz) 438 dqikz = dix*qky - diy*qkx + dkx*qiy - dky*qix 439 & - 2.0d0*(qixx*qkxy+qixy*qkyy+qixz*qkyz 440 & -qixy*qkxx-qiyy*qkxy-qiyz*qkxz) 441c 442c get reciprocal distance terms for this interaction 443c 444 rr1 = f * mscale(k) / r 445 rr3 = rr1 / r2 446 rr5 = 3.0d0 * rr3 / r2 447 rr7 = 5.0d0 * rr5 / r2 448 rr9 = 7.0d0 * rr7 / r2 449 rr11 = 9.0d0 * rr9 / r2 450c 451c find damped multipole intermediates for force and torque 452c 453 if (use_chgpen) then 454 corek = pcore(kk) 455 valk = pval(kk) 456 alphak = palpha(kk) 457 term1 = corei*corek 458 term1i = corek*vali 459 term2i = corek*dir 460 term3i = corek*qir 461 term1k = corei*valk 462 term2k = -corei*dkr 463 term3k = corei*qkr 464 term1ik = vali*valk 465 term2ik = valk*dir - vali*dkr + dik 466 term3ik = vali*qkr + valk*qir - dir*dkr 467 & + 2.0d0*(dkqi-diqk+qiqk) 468 term4ik = dir*qkr - dkr*qir - 4.0d0*qik 469 term5ik = qir*qkr 470 call damppole (r,11,alphai,alphak, 471 & dmpi,dmpk,dmpik) 472 rr3i = dmpi(3)*rr3 473 rr5i = dmpi(5)*rr5 474 rr7i = dmpi(7)*rr7 475 rr3k = dmpk(3)*rr3 476 rr5k = dmpk(5)*rr5 477 rr7k = dmpk(7)*rr7 478 rr3ik = dmpik(3)*rr3 479 rr5ik = dmpik(5)*rr5 480 rr7ik = dmpik(7)*rr7 481 rr9ik = dmpik(9)*rr9 482 rr11ik = dmpik(11)*rr11 483 de = term1*rr3 + term4ik*rr9ik + term5ik*rr11ik 484 & + term1i*rr3i + term1k*rr3k + term1ik*rr3ik 485 & + term2i*rr5i + term2k*rr5k + term2ik*rr5ik 486 & + term3i*rr7i + term3k*rr7k + term3ik*rr7ik 487 term1 = -corek*rr3i - valk*rr3ik 488 & + dkr*rr5ik - qkr*rr7ik 489 term2 = corei*rr3k + vali*rr3ik 490 & + dir*rr5ik + qir*rr7ik 491 term3 = 2.0d0 * rr5ik 492 term4 = -2.0d0 * (corek*rr5i+valk*rr5ik 493 & -dkr*rr7ik+qkr*rr9ik) 494 term5 = -2.0d0 * (corei*rr5k+vali*rr5ik 495 & +dir*rr7ik+qir*rr9ik) 496 term6 = 4.0d0 * rr7ik 497c 498c find standard multipole intermediates for force and torque 499c 500 else 501 term1 = ci*ck 502 term2 = ck*dir - ci*dkr + dik 503 term3 = ci*qkr + ck*qir - dir*dkr 504 & + 2.0d0*(dkqi-diqk+qiqk) 505 term4 = dir*qkr - dkr*qir - 4.0d0*qik 506 term5 = qir*qkr 507 de = term1*rr3 + term2*rr5 + term3*rr7 508 & + term4*rr9 + term5*rr11 509 term1 = -ck*rr3 + dkr*rr5 - qkr*rr7 510 term2 = ci*rr3 + dir*rr5 + qir*rr7 511 term3 = 2.0d0 * rr5 512 term4 = 2.0d0 * (-ck*rr5+dkr*rr7-qkr*rr9) 513 term5 = 2.0d0 * (-ci*rr5-dir*rr7-qir*rr9) 514 term6 = 4.0d0 * rr7 515 end if 516c 517c store the potential at each site for use in charge flux 518c 519 if (use_chgflx) then 520 if (use_chgpen) then 521 term1i = corek*dmpi(1) + valk*dmpik(1) 522 term1k = corei*dmpk(1) + vali*dmpik(1) 523 term2i = -dkr * dmpik(3) 524 term2k = dir * dmpik(3) 525 term3i = qkr * dmpik(5) 526 term3k = qir * dmpik(5) 527 poti = term1i*rr1 + term2i*rr3 + term3i*rr5 528 potk = term1k*rr1 + term2k*rr3 + term3k*rr5 529 else 530 poti = ck*rr1 - dkr*rr3 + qkr*rr5 531 potk = ci*rr1 + dir*rr3 + qir*rr5 532 end if 533 pot(i) = pot(i) + poti 534 pot(k) = pot(k) + potk 535 end if 536c 537c compute the force components for this interaction 538c 539 frcx = de*xr + term1*dix + term2*dkx 540 & + term3*(diqkx-dkqix) + term4*qix 541 & + term5*qkx + term6*(qixk+qkxi) 542 frcy = de*yr + term1*diy + term2*dky 543 & + term3*(diqky-dkqiy) + term4*qiy 544 & + term5*qky + term6*(qiyk+qkyi) 545 frcz = de*zr + term1*diz + term2*dkz 546 & + term3*(diqkz-dkqiz) + term4*qiz 547 & + term5*qkz + term6*(qizk+qkzi) 548c 549c compute the torque components for this interaction 550c 551 if (use_chgpen) rr3 = rr3ik 552 ttmi(1) = -rr3*dikx + term1*dirx 553 & + term3*(dqikx+dkqirx) 554 & - term4*qirx - term6*(qikrx+qikx) 555 ttmi(2) = -rr3*diky + term1*diry 556 & + term3*(dqiky+dkqiry) 557 & - term4*qiry - term6*(qikry+qiky) 558 ttmi(3) = -rr3*dikz + term1*dirz 559 & + term3*(dqikz+dkqirz) 560 & - term4*qirz - term6*(qikrz+qikz) 561 ttmk(1) = rr3*dikx + term2*dkrx 562 & - term3*(dqikx+diqkrx) 563 & - term5*qkrx - term6*(qkirx-qikx) 564 ttmk(2) = rr3*diky + term2*dkry 565 & - term3*(dqiky+diqkry) 566 & - term5*qkry - term6*(qkiry-qiky) 567 ttmk(3) = rr3*dikz + term2*dkrz 568 & - term3*(dqikz+diqkrz) 569 & - term5*qkrz - term6*(qkirz-qikz) 570c 571c force and torque components scaled by group membership 572c 573 if (use_group) then 574 frcx = fgrp * frcx 575 frcy = fgrp * frcy 576 frcz = fgrp * frcz 577 do j = 1, 3 578 ttmi(j) = fgrp * ttmi(j) 579 ttmk(j) = fgrp * ttmk(j) 580 end do 581 end if 582c 583c increment force-based gradient and torque on first site 584c 585 dem(1,i) = dem(1,i) + frcx 586 dem(2,i) = dem(2,i) + frcy 587 dem(3,i) = dem(3,i) + frcz 588 tem(1,i) = tem(1,i) + ttmi(1) 589 tem(2,i) = tem(2,i) + ttmi(2) 590 tem(3,i) = tem(3,i) + ttmi(3) 591c 592c increment force-based gradient and torque on second site 593c 594 dem(1,k) = dem(1,k) - frcx 595 dem(2,k) = dem(2,k) - frcy 596 dem(3,k) = dem(3,k) - frcz 597 tem(1,k) = tem(1,k) + ttmk(1) 598 tem(2,k) = tem(2,k) + ttmk(2) 599 tem(3,k) = tem(3,k) + ttmk(3) 600 end if 601 10 continue 602 end do 603c 604c reset exclusion coefficients for connected atoms 605c 606 do j = 1, n12(i) 607 mscale(i12(j,i)) = 1.0d0 608 end do 609 do j = 1, n13(i) 610 mscale(i13(j,i)) = 1.0d0 611 end do 612 do j = 1, n14(i) 613 mscale(i14(j,i)) = 1.0d0 614 end do 615 do j = 1, n15(i) 616 mscale(i15(j,i)) = 1.0d0 617 end do 618 end do 619c 620c for periodic boundary conditions with large cutoffs 621c neighbors must be found by the replicates method 622c 623 if (use_replica) then 624c 625c calculate interaction with other unit cells 626c 627 do iii = 1, nlist 628 ii = list(iii) 629 i = ipole(ii) 630 iz = zaxis(ii) 631 ix = xaxis(ii) 632 iy = abs(yaxis(ii)) 633 xi = x(i) 634 yi = y(i) 635 zi = z(i) 636 ci = rpole(1,ii) 637 dix = rpole(2,ii) 638 diy = rpole(3,ii) 639 diz = rpole(4,ii) 640 qixx = rpole(5,ii) 641 qixy = rpole(6,ii) 642 qixz = rpole(7,ii) 643 qiyy = rpole(9,ii) 644 qiyz = rpole(10,ii) 645 qizz = rpole(13,ii) 646 if (use_chgpen) then 647 corei = pcore(ii) 648 vali = pval(ii) 649 alphai = palpha(ii) 650 end if 651 usei = (use(i) .or. use(iz) .or. use(ix) .or. use(iy)) 652c 653c set exclusion coefficients for connected atoms 654c 655 do j = 1, n12(i) 656 mscale(i12(j,i)) = m2scale 657 end do 658 do j = 1, n13(i) 659 mscale(i13(j,i)) = m3scale 660 end do 661 do j = 1, n14(i) 662 mscale(i14(j,i)) = m4scale 663 end do 664 do j = 1, n15(i) 665 mscale(i15(j,i)) = m5scale 666 end do 667c 668c evaluate all sites within the cutoff distance 669c 670 do kk = 1, npole 671 k = ipole(kk) 672 kz = zaxis(kk) 673 kx = xaxis(kk) 674 ky = abs(yaxis(kk)) 675 usek = (use(k) .or. use(kz) .or. use(kx) .or. use(ky)) 676 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 677 proceed = .true. 678 if (proceed) proceed = (usei .or. usek) 679 if (.not. proceed) goto 20 680 do jcell = 2, ncell 681 xr = x(k) - xi 682 yr = y(k) - yi 683 zr = z(k) - zi 684 call imager (xr,yr,zr,jcell) 685 r2 = xr*xr + yr*yr + zr*zr 686 if (.not. (use_polymer .and. r2.le.polycut2)) then 687 mscale(k) = 1.0d0 688 end if 689 if (r2 .le. off2) then 690 r = sqrt(r2) 691 ck = rpole(1,kk) 692 dkx = rpole(2,kk) 693 dky = rpole(3,kk) 694 dkz = rpole(4,kk) 695 qkxx = rpole(5,kk) 696 qkxy = rpole(6,kk) 697 qkxz = rpole(7,kk) 698 qkyy = rpole(9,kk) 699 qkyz = rpole(10,kk) 700 qkzz = rpole(13,kk) 701c 702c intermediates involving moments and separation distance 703c 704 dir = dix*xr + diy*yr + diz*zr 705 qix = qixx*xr + qixy*yr + qixz*zr 706 qiy = qixy*xr + qiyy*yr + qiyz*zr 707 qiz = qixz*xr + qiyz*yr + qizz*zr 708 qir = qix*xr + qiy*yr + qiz*zr 709 dkr = dkx*xr + dky*yr + dkz*zr 710 qkx = qkxx*xr + qkxy*yr + qkxz*zr 711 qky = qkxy*xr + qkyy*yr + qkyz*zr 712 qkz = qkxz*xr + qkyz*yr + qkzz*zr 713 qkr = qkx*xr + qky*yr + qkz*zr 714 dik = dix*dkx + diy*dky + diz*dkz 715 qik = qix*qkx + qiy*qky + qiz*qkz 716 diqk = dix*qkx + diy*qky + diz*qkz 717 dkqi = dkx*qix + dky*qiy + dkz*qiz 718 qiqk = 2.0d0*(qixy*qkxy+qixz*qkxz+qiyz*qkyz) 719 & + qixx*qkxx + qiyy*qkyy + qizz*qkzz 720c 721c additional intermediates involving moments and distance 722c 723 dirx = diy*zr - diz*yr 724 diry = diz*xr - dix*zr 725 dirz = dix*yr - diy*xr 726 dkrx = dky*zr - dkz*yr 727 dkry = dkz*xr - dkx*zr 728 dkrz = dkx*yr - dky*xr 729 dikx = diy*dkz - diz*dky 730 diky = diz*dkx - dix*dkz 731 dikz = dix*dky - diy*dkx 732 qirx = qiz*yr - qiy*zr 733 qiry = qix*zr - qiz*xr 734 qirz = qiy*xr - qix*yr 735 qkrx = qkz*yr - qky*zr 736 qkry = qkx*zr - qkz*xr 737 qkrz = qky*xr - qkx*yr 738 qikx = qky*qiz - qkz*qiy 739 qiky = qkz*qix - qkx*qiz 740 qikz = qkx*qiy - qky*qix 741 qixk = qixx*qkx + qixy*qky + qixz*qkz 742 qiyk = qixy*qkx + qiyy*qky + qiyz*qkz 743 qizk = qixz*qkx + qiyz*qky + qizz*qkz 744 qkxi = qkxx*qix + qkxy*qiy + qkxz*qiz 745 qkyi = qkxy*qix + qkyy*qiy + qkyz*qiz 746 qkzi = qkxz*qix + qkyz*qiy + qkzz*qiz 747 qikrx = qizk*yr - qiyk*zr 748 qikry = qixk*zr - qizk*xr 749 qikrz = qiyk*xr - qixk*yr 750 qkirx = qkzi*yr - qkyi*zr 751 qkiry = qkxi*zr - qkzi*xr 752 qkirz = qkyi*xr - qkxi*yr 753 diqkx = dix*qkxx + diy*qkxy + diz*qkxz 754 diqky = dix*qkxy + diy*qkyy + diz*qkyz 755 diqkz = dix*qkxz + diy*qkyz + diz*qkzz 756 dkqix = dkx*qixx + dky*qixy + dkz*qixz 757 dkqiy = dkx*qixy + dky*qiyy + dkz*qiyz 758 dkqiz = dkx*qixz + dky*qiyz + dkz*qizz 759 diqkrx = diqkz*yr - diqky*zr 760 diqkry = diqkx*zr - diqkz*xr 761 diqkrz = diqky*xr - diqkx*yr 762 dkqirx = dkqiz*yr - dkqiy*zr 763 dkqiry = dkqix*zr - dkqiz*xr 764 dkqirz = dkqiy*xr - dkqix*yr 765 dqikx = diy*qkz - diz*qky + dky*qiz - dkz*qiy 766 & - 2.0d0*(qixy*qkxz+qiyy*qkyz+qiyz*qkzz 767 & -qixz*qkxy-qiyz*qkyy-qizz*qkyz) 768 dqiky = diz*qkx - dix*qkz + dkz*qix - dkx*qiz 769 & - 2.0d0*(qixz*qkxx+qiyz*qkxy+qizz*qkxz 770 & -qixx*qkxz-qixy*qkyz-qixz*qkzz) 771 dqikz = dix*qky - diy*qkx + dkx*qiy - dky*qix 772 & - 2.0d0*(qixx*qkxy+qixy*qkyy+qixz*qkyz 773 & -qixy*qkxx-qiyy*qkxy-qiyz*qkxz) 774c 775c get reciprocal distance terms for this interaction 776c 777 rr1 = f * mscale(k) / r 778 rr3 = rr1 / r2 779 rr5 = 3.0d0 * rr3 / r2 780 rr7 = 5.0d0 * rr5 / r2 781 rr9 = 7.0d0 * rr7 / r2 782 rr11 = 9.0d0 * rr9 / r2 783c 784c find damped multipole intermediates for force and torque 785c 786 if (use_chgpen) then 787 corek = pcore(kk) 788 valk = pval(kk) 789 alphak = palpha(kk) 790 term1 = corei*corek 791 term1i = corek*vali 792 term2i = corek*dir 793 term3i = corek*qir 794 term1k = corei*valk 795 term2k = -corei*dkr 796 term3k = corei*qkr 797 term1ik = vali*valk 798 term2ik = valk*dir - vali*dkr + dik 799 term3ik = vali*qkr + valk*qir - dir*dkr 800 & + 2.0d0*(dkqi-diqk+qiqk) 801 term4ik = dir*qkr - dkr*qir - 4.0d0*qik 802 term5ik = qir*qkr 803 call damppole (r,11,alphai,alphak, 804 & dmpi,dmpk,dmpik) 805 rr3i = dmpi(3)*rr3 806 rr5i = dmpi(5)*rr5 807 rr7i = dmpi(7)*rr7 808 rr3k = dmpk(3)*rr3 809 rr5k = dmpk(5)*rr5 810 rr7k = dmpk(7)*rr7 811 rr3ik = dmpik(3)*rr3 812 rr5ik = dmpik(5)*rr5 813 rr7ik = dmpik(7)*rr7 814 rr9ik = dmpik(9)*rr9 815 rr11ik = dmpik(11)*rr11 816 de = term1*rr3 + term4ik*rr9ik + term5ik*rr11ik 817 & + term1i*rr3i + term1k*rr3k + term1ik*rr3ik 818 & + term2i*rr5i + term2k*rr5k + term2ik*rr5ik 819 & + term3i*rr7i + term3k*rr7k + term3ik*rr7ik 820 term1 = -corek*rr3i - valk*rr3ik 821 & + dkr*rr5ik - qkr*rr7ik 822 term2 = corei*rr3k + vali*rr3ik 823 & + dir*rr5ik + qir*rr7ik 824 term3 = 2.0d0 * rr5ik 825 term4 = -2.0d0 * (corek*rr5i+valk*rr5ik 826 & -dkr*rr7ik+qkr*rr9ik) 827 term5 = -2.0d0 * (corei*rr5k+vali*rr5ik 828 & +dir*rr7ik+qir*rr9ik) 829 term6 = 4.0d0 * rr7ik 830 rr3 = rr3ik 831c 832c find standard multipole intermediates for force and torque 833c 834 else 835 term1 = ci*ck 836 term2 = ck*dir - ci*dkr + dik 837 term3 = ci*qkr + ck*qir - dir*dkr 838 & + 2.0d0*(dkqi-diqk+qiqk) 839 term4 = dir*qkr - dkr*qir - 4.0d0*qik 840 term5 = qir*qkr 841 de = term1*rr3 + term2*rr5 + term3*rr7 842 & + term4*rr9 + term5*rr11 843 term1 = -ck*rr3 + dkr*rr5 - qkr*rr7 844 term2 = ci*rr3 + dir*rr5 + qir*rr7 845 term3 = 2.0d0 * rr5 846 term4 = 2.0d0 * (-ck*rr5+dkr*rr7-qkr*rr9) 847 term5 = 2.0d0 * (-ci*rr5-dir*rr7-qir*rr9) 848 term6 = 4.0d0 * rr7 849 end if 850c 851c store the potential at each site for use in charge flux 852c 853 if (use_chgflx) then 854 if (use_chgpen) then 855 term1i = corek*dmpi(1) + valk*dmpik(1) 856 term1k = corei*dmpk(1) + vali*dmpik(1) 857 term2i = -dkr * dmpik(3) 858 term2k = dir * dmpik(3) 859 term3i = qkr * dmpik(5) 860 term3k = qir * dmpik(5) 861 poti = term1i*rr1 + term2i*rr3 + term3i*rr5 862 potk = term1k*rr1 + term2k*rr3 + term3k*rr5 863 else 864 poti = ck*rr1 - dkr*rr3 + qkr*rr5 865 potk = ci*rr1 + dir*rr3 + qir*rr5 866 end if 867 pot(i) = pot(i) + poti 868 pot(k) = pot(k) + potk 869 end if 870c 871c compute the force components for this interaction 872c 873 frcx = de*xr + term1*dix + term2*dkx 874 & + term3*(diqkx-dkqix) + term4*qix 875 & + term5*qkx + term6*(qixk+qkxi) 876 frcy = de*yr + term1*diy + term2*dky 877 & + term3*(diqky-dkqiy) + term4*qiy 878 & + term5*qky + term6*(qiyk+qkyi) 879 frcz = de*zr + term1*diz + term2*dkz 880 & + term3*(diqkz-dkqiz) + term4*qiz 881 & + term5*qkz + term6*(qizk+qkzi) 882c 883c compute the torque components for this interaction 884c 885 ttmi(1) = -rr3*dikx + term1*dirx 886 & + term3*(dqikx+dkqirx) 887 & - term4*qirx - term6*(qikrx+qikx) 888 ttmi(2) = -rr3*diky + term1*diry 889 & + term3*(dqiky+dkqiry) 890 & - term4*qiry - term6*(qikry+qiky) 891 ttmi(3) = -rr3*dikz + term1*dirz 892 & + term3*(dqikz+dkqirz) 893 & - term4*qirz - term6*(qikrz+qikz) 894 ttmk(1) = rr3*dikx + term2*dkrx 895 & - term3*(dqikx+diqkrx) 896 & - term5*qkrx - term6*(qkirx-qikx) 897 ttmk(2) = rr3*diky + term2*dkry 898 & - term3*(dqiky+diqkry) 899 & - term5*qkry - term6*(qkiry-qiky) 900 ttmk(3) = rr3*dikz + term2*dkrz 901 & - term3*(dqikz+diqkrz) 902 & - term5*qkrz - term6*(qkirz-qikz) 903c 904c force and torque scaled for self-interactions and groups 905c 906 if (i .eq. k) then 907 frcx = 0.5d0 * frcx 908 frcy = 0.5d0 * frcy 909 frcz = 0.5d0 * frcz 910 do j = 1, 3 911 ttmi(j) = 0.5d0 * ttmi(j) 912 ttmk(j) = 0.5d0 * ttmk(j) 913 end do 914 end if 915 if (use_group) then 916 frcx = fgrp * frcx 917 frcy = fgrp * frcy 918 frcz = fgrp * frcz 919 do j = 1, 3 920 ttmi(j) = fgrp * ttmi(j) 921 ttmk(j) = fgrp * ttmk(j) 922 end do 923 end if 924c 925c increment force-based gradient and torque on first site 926c 927 dem(1,i) = dem(1,i) + frcx 928 dem(2,i) = dem(2,i) + frcy 929 dem(3,i) = dem(3,i) + frcz 930 tem(1,i) = tem(1,i) + ttmi(1) 931 tem(2,i) = tem(2,i) + ttmi(2) 932 tem(3,i) = tem(3,i) + ttmi(3) 933c 934c increment force-based gradient and torque on second site 935c 936 dem(1,k) = dem(1,k) - frcx 937 dem(2,k) = dem(2,k) - frcy 938 dem(3,k) = dem(3,k) - frcz 939 tem(1,k) = tem(1,k) + ttmk(1) 940 tem(2,k) = tem(2,k) + ttmk(2) 941 tem(3,k) = tem(3,k) + ttmk(3) 942 end if 943 end do 944 20 continue 945 end do 946c 947c reset exclusion coefficients for connected atoms 948c 949 do j = 1, n12(i) 950 mscale(i12(j,i)) = 1.0d0 951 end do 952 do j = 1, n13(i) 953 mscale(i13(j,i)) = 1.0d0 954 end do 955 do j = 1, n14(i) 956 mscale(i14(j,i)) = 1.0d0 957 end do 958 do j = 1, n15(i) 959 mscale(i15(j,i)) = 1.0d0 960 end do 961 end do 962 end if 963c 964c resolve site torques then increment multipole forces 965c 966 do ii = 1, npole 967 i = ipole(ii) 968 call torque (ii,tem(1,i),fix,fiy,fiz,dem) 969 end do 970c 971c modify the gradient and virial for charge flux 972c 973 if (use_chgflx) then 974 call dcflux (pot,decfx,decfy,decfz) 975 do ii = 1, npole 976 i = ipole(ii) 977 frcx = decfx(i) 978 frcy = decfy(i) 979 frcz = decfz(i) 980 dem(1,i) = dem(1,i) + frcx 981 dem(2,i) = dem(2,i) + frcy 982 dem(3,i) = dem(3,i) + frcz 983 end do 984 end if 985c 986c perform deallocation of some local arrays 987c 988 deallocate (mscale) 989 deallocate (tem) 990 deallocate (pot) 991 deallocate (decfx) 992 deallocate (decfy) 993 deallocate (decfz) 994 return 995 end 996