1c 2c 3c ############################################################ 4c ## COPYRIGHT (C) 1995 by Yong Kong & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ############################################################ 7c 8c ################################################################# 9c ## ## 10c ## subroutine rotpole -- rotate multipoles to global frame ## 11c ## ## 12c ################################################################# 13c 14c 15c "rotpole" constructs the set of atomic multipoles in the global 16c frame by applying the correct rotation matrix for each site 17c 18c 19 subroutine rotpole 20 use mpole 21 implicit none 22 integer i 23 real*8 a(3,3) 24c real*8 d1(3,3) 25c real*8 d2(5,5) 26c logical use_qi 27c 28c 29c rotate the atomic multipoles at each site in turn 30c 31 do i = 1, npole 32 call rotmat (i,a) 33 call rotsite (i,a) 34 end do 35c 36c set QI dipole rotation matrix by permuting the Cartesian 37c rotation matrix to adhere to spherical harmonic ordering 38c 39c do i = 1, npole 40c d1(1,1) = a(3,3) 41c d1(2,1) = a(3,1) 42c d1(3,1) = a(3,2) 43c d1(1,2) = a(1,3) 44c d1(2,2) = a(1,1) 45c d1(3,2) = a(1,2) 46c d1(1,3) = a(2,3) 47c d1(2,3) = a(2,1) 48c d1(3,3) = a(2,2) 49c call shrotmat (d1,d2) 50c call shrotsite (i,d1,d2) 51c end do 52 return 53 end 54c 55c 56c ################################################################ 57c ## ## 58c ## subroutine rotmat -- find global frame rotation matrix ## 59c ## ## 60c ################################################################ 61c 62c 63c "rotmat" finds the rotation matrix that rotates the local 64c coordinate system into the global frame at a multipole site 65c 66c 67 subroutine rotmat (i,a) 68 use atoms 69 use mpole 70 implicit none 71 integer i,ii 72 integer ix,iy,iz 73 real*8 r,dot 74 real*8 xi,yi,zi 75 real*8 dx,dy,dz 76 real*8 dx1,dy1,dz1 77 real*8 dx2,dy2,dz2 78 real*8 dx3,dy3,dz3 79 real*8 a(3,3) 80c 81c 82c get coordinates and frame definition for the multipole site 83c 84 ii = ipole(i) 85 xi = x(ii) 86 yi = y(ii) 87 zi = z(ii) 88 iz = zaxis(i) 89 ix = xaxis(i) 90 iy = abs(yaxis(i)) 91c 92c use the identity matrix as the default rotation matrix 93c 94 a(1,1) = 1.0d0 95 a(2,1) = 0.0d0 96 a(3,1) = 0.0d0 97 a(1,3) = 0.0d0 98 a(2,3) = 0.0d0 99 a(3,3) = 1.0d0 100c 101c z-only frame rotation matrix elements for z-axis only 102c 103 if (polaxe(i) .eq. 'Z-Only') then 104 dx = x(iz) - xi 105 dy = y(iz) - yi 106 dz = z(iz) - zi 107 r = sqrt(dx*dx + dy*dy + dz*dz) 108 a(1,3) = dx / r 109 a(2,3) = dy / r 110 a(3,3) = dz / r 111 dx = 1.0d0 112 dy = 0.0d0 113 dz = 0.0d0 114 dot = a(1,3) 115 if (abs(dot) .gt. 0.866d0) then 116 dx = 0.0d0 117 dy = 1.0d0 118 dot = a(2,3) 119 end if 120 dx = dx - dot*a(1,3) 121 dy = dy - dot*a(2,3) 122 dz = dz - dot*a(3,3) 123 r = sqrt(dx*dx + dy*dy + dz*dz) 124 a(1,1) = dx / r 125 a(2,1) = dy / r 126 a(3,1) = dz / r 127c 128c z-then-x frame rotation matrix elements for z- and x-axes 129c 130 else if (polaxe(i) .eq. 'Z-then-X') then 131 dx = x(iz) - xi 132 dy = y(iz) - yi 133 dz = z(iz) - zi 134 r = sqrt(dx*dx + dy*dy + dz*dz) 135 a(1,3) = dx / r 136 a(2,3) = dy / r 137 a(3,3) = dz / r 138 dx = x(ix) - xi 139 dy = y(ix) - yi 140 dz = z(ix) - zi 141 dot = dx*a(1,3) + dy*a(2,3) + dz*a(3,3) 142 dx = dx - dot*a(1,3) 143 dy = dy - dot*a(2,3) 144 dz = dz - dot*a(3,3) 145 r = sqrt(dx*dx + dy*dy + dz*dz) 146 a(1,1) = dx / r 147 a(2,1) = dy / r 148 a(3,1) = dz / r 149c 150c bisector frame rotation matrix elements for z- and x-axes 151c 152 else if (polaxe(i) .eq. 'Bisector') then 153 dx = x(iz) - xi 154 dy = y(iz) - yi 155 dz = z(iz) - zi 156 r = sqrt(dx*dx + dy*dy + dz*dz) 157 dx1 = dx / r 158 dy1 = dy / r 159 dz1 = dz / r 160 dx = x(ix) - xi 161 dy = y(ix) - yi 162 dz = z(ix) - zi 163 r = sqrt(dx*dx + dy*dy + dz*dz) 164 dx2 = dx / r 165 dy2 = dy / r 166 dz2 = dz / r 167 dx = dx1 + dx2 168 dy = dy1 + dy2 169 dz = dz1 + dz2 170 r = sqrt(dx*dx + dy*dy + dz*dz) 171 a(1,3) = dx / r 172 a(2,3) = dy / r 173 a(3,3) = dz / r 174 dot = dx2*a(1,3) + dy2*a(2,3) + dz2*a(3,3) 175 dx = dx2 - dot*a(1,3) 176 dy = dy2 - dot*a(2,3) 177 dz = dz2 - dot*a(3,3) 178 r = sqrt(dx*dx + dy*dy + dz*dz) 179 a(1,1) = dx / r 180 a(2,1) = dy / r 181 a(3,1) = dz / r 182c 183c z-bisect frame rotation matrix elements for z- and x-axes 184c 185 else if (polaxe(i) .eq. 'Z-Bisect') then 186 dx = x(iz) - xi 187 dy = y(iz) - yi 188 dz = z(iz) - zi 189 r = sqrt(dx*dx + dy*dy + dz*dz) 190 a(1,3) = dx / r 191 a(2,3) = dy / r 192 a(3,3) = dz / r 193 dx = x(ix) - xi 194 dy = y(ix) - yi 195 dz = z(ix) - zi 196 r = sqrt(dx*dx + dy*dy + dz*dz) 197 dx1 = dx / r 198 dy1 = dy / r 199 dz1 = dz / r 200 dx = x(iy) - xi 201 dy = y(iy) - yi 202 dz = z(iy) - zi 203 r = sqrt(dx*dx + dy*dy + dz*dz) 204 dx2 = dx / r 205 dy2 = dy / r 206 dz2 = dz / r 207 dx = dx1 + dx2 208 dy = dy1 + dy2 209 dz = dz1 + dz2 210 r = sqrt(dx*dx + dy*dy + dz*dz) 211 dx = dx / r 212 dy = dy / r 213 dz = dz / r 214 dot = dx*a(1,3) + dy*a(2,3) + dz*a(3,3) 215 dx = dx - dot*a(1,3) 216 dy = dy - dot*a(2,3) 217 dz = dz - dot*a(3,3) 218 r = sqrt(dx*dx + dy*dy + dz*dz) 219 a(1,1) = dx / r 220 a(2,1) = dy / r 221 a(3,1) = dz / r 222c 223c 3-fold frame rotation matrix elements for z- and x-axes 224c 225 else if (polaxe(i) .eq. '3-Fold') then 226 dx = x(iz) - xi 227 dy = y(iz) - yi 228 dz = z(iz) - zi 229 r = sqrt(dx*dx + dy*dy + dz*dz) 230 dx1 = dx / r 231 dy1 = dy / r 232 dz1 = dz / r 233 dx = x(ix) - xi 234 dy = y(ix) - yi 235 dz = z(ix) - zi 236 r = sqrt(dx*dx + dy*dy + dz*dz) 237 dx2 = dx / r 238 dy2 = dy / r 239 dz2 = dz / r 240 dx = x(iy) - xi 241 dy = y(iy) - yi 242 dz = z(iy) - zi 243 r = sqrt(dx*dx + dy*dy + dz*dz) 244 dx3 = dx / r 245 dy3 = dy / r 246 dz3 = dz / r 247 dx = dx1 + dx2 + dx3 248 dy = dy1 + dy2 + dy3 249 dz = dz1 + dz2 + dz3 250 r = sqrt(dx*dx + dy*dy + dz*dz) 251 a(1,3) = dx / r 252 a(2,3) = dy / r 253 a(3,3) = dz / r 254 dot = dx2*a(1,3) + dy2*a(2,3) + dz2*a(3,3) 255 dx = dx2 - dot*a(1,3) 256 dy = dy2 - dot*a(2,3) 257 dz = dz2 - dot*a(3,3) 258 r = sqrt(dx*dx + dy*dy + dz*dz) 259 a(1,1) = dx / r 260 a(2,1) = dy / r 261 a(3,1) = dz / r 262 end if 263c 264c finally, find rotation matrix elements for the y-axis 265c 266 a(1,2) = a(3,1)*a(2,3) - a(2,1)*a(3,3) 267 a(2,2) = a(1,1)*a(3,3) - a(3,1)*a(1,3) 268 a(3,2) = a(2,1)*a(1,3) - a(1,1)*a(2,3) 269 return 270 end 271c 272c 273c ################################################################ 274c ## ## 275c ## subroutine rotsite -- rotate multipoles at single site ## 276c ## ## 277c ################################################################ 278c 279c 280c "rotsite" rotates the local frame atomic multipoles at a 281c specified site into the global coordinate frame by applying 282c a rotation matrix 283c 284c 285 subroutine rotsite (isite,a) 286 use atoms 287 use mpole 288 implicit none 289 integer i,j,k,m 290 integer isite 291 real*8 a(3,3) 292 real*8 mp(3,3) 293 real*8 rp(3,3) 294c 295c 296c monopoles have the same value in any coordinate frame 297c 298 rpole(1,isite) = pole(1,isite) 299c 300c rotate the dipoles to the global coordinate frame 301c 302 do i = 2, 4 303 rpole(i,isite) = 0.0d0 304 do j = 2, 4 305 rpole(i,isite) = rpole(i,isite) + pole(j,isite)*a(i-1,j-1) 306 end do 307 end do 308c 309c rotate the quadrupoles to the global coordinate frame 310c 311 k = 5 312 do i = 1, 3 313 do j = 1, 3 314 mp(i,j) = pole(k,isite) 315 rp(i,j) = 0.0d0 316 k = k + 1 317 end do 318 end do 319 do i = 1, 3 320 do j = 1, 3 321 if (j .lt. i) then 322 rp(i,j) = rp(j,i) 323 else 324 do k = 1, 3 325 do m = 1, 3 326 rp(i,j) = rp(i,j) + a(i,k)*a(j,m)*mp(k,m) 327 end do 328 end do 329 end if 330 end do 331 end do 332 k = 5 333 do i = 1, 3 334 do j = 1, 3 335 rpole(k,isite) = rp(i,j) 336 k = k + 1 337 end do 338 end do 339 return 340 end 341c 342c 343c ############################################################## 344c ## ## 345c ## subroutine shrotmat -- QI quadrupole rotation matrix ## 346c ## ## 347c ############################################################## 348c 349c 350c "shrotmat" finds the rotation matrix that converts spherical 351c harmonic quadrupoles from the local to the global frame given 352c the required dipole rotation matrix 353c 354c 355 subroutine shrotmat (d1,d2) 356 use math 357 implicit none 358 real*8 d1(3,3) 359 real*8 d2(5,5) 360c 361c 362c build the quadrupole rotation matrix 363c 364 d2(1,1) = 0.5d0*(3.0d0*d1(1,1)*d1(1,1) - 1.0d0) 365 d2(2,1) = root3*d1(1,1)*d1(2,1) 366 d2(3,1) = root3*d1(1,1)*d1(3,1) 367 d2(4,1) = 0.5d0*root3*(d1(2,1)*d1(2,1) - d1(3,1)*d1(3,1)) 368 d2(5,1) = root3*d1(2,1)*d1(3,1) 369 d2(1,2) = root3*d1(1,1)*d1(1,2) 370 d2(2,2) = d1(2,1)*d1(1,2) + d1(1,1)*d1(2,2) 371 d2(3,2) = d1(3,1)*d1(1,2) + d1(1,1)*d1(3,2) 372 d2(4,2) = d1(2,1)*d1(2,2) - d1(3,1)*d1(3,2) 373 d2(5,2) = d1(3,1)*d1(2,2) + d1(2,1)*d1(3,2) 374 d2(1,3) = root3*d1(1,1)*d1(1,3) 375 d2(2,3) = d1(2,1)*d1(1,3) + d1(1,1)*d1(2,3) 376 d2(3,3) = d1(3,1)*d1(1,3) + d1(1,1)*d1(3,3) 377 d2(4,3) = d1(2,1)*d1(2,3) - d1(3,1)*d1(3,3) 378 d2(5,3) = d1(3,1)*d1(2,3) + d1(2,1)*d1(3,3) 379 d2(1,4) = 0.5d0*root3*(d1(1,2)*d1(1,2) - d1(1,3)*d1(1,3)) 380 d2(2,4) = d1(1,2)*d1(2,2) - d1(1,3)*d1(2,3) 381 d2(3,4) = d1(1,2)*d1(3,2) - d1(1,3)*d1(3,3) 382 d2(4,4) = 0.5d0*(d1(2,2)*d1(2,2) - d1(3,2)*d1(3,2) 383 & - d1(2,3)*d1(2,3) + d1(3,3)*d1(3,3)) 384 d2(5,4) = d1(2,2)*d1(3,2) - d1(2,3)*d1(3,3) 385 d2(1,5) = root3*d1(1,2)*d1(1,3) 386 d2(2,5) = d1(2,2)*d1(1,3) + d1(1,2)*d1(2,3) 387 d2(3,5) = d1(3,2)*d1(1,3) + d1(1,2)*d1(3,3) 388 d2(4,5) = d1(2,2)*d1(2,3) - d1(3,2)*d1(3,3) 389 d2(5,5) = d1(3,2)*d1(2,3) + d1(2,2)*d1(3,3) 390 return 391 end 392c 393c 394c ################################################################## 395c ## ## 396c ## subroutine shrotsite -- rotate SH mpoles local to global ## 397c ## ## 398c ################################################################## 399c 400c 401c "shrotsite" converts spherical harmonic multipoles from the 402c local to the global frame given required rotation matrices 403c 404c 405 subroutine shrotsite (i,d1,d2) 406 use mpole 407 implicit none 408 integer i 409 real*8 d1(3,3) 410 real*8 d2(5,5) 411c 412c 413c find the global frame spherical harmonic multipoles 414c 415 srpole(1,i) = spole(1,i) 416 srpole(2,i) = spole(2,i)*d1(1,1) + spole(3,i)*d1(2,1) 417 & + spole(4,i)*d1(3,1) 418 srpole(3,i) = spole(2,i)*d1(1,2) + spole(3,i)*d1(2,2) 419 & + spole(4,i)*d1(3,2) 420 srpole(4,i) = spole(2,i)*d1(1,3) + spole(3,i)*d1(2,3) 421 & + spole(4,i)*d1(3,3) 422 srpole(5,i) = spole(5,i)*d2(1,1) + spole(6,i)*d2(2,1) 423 & + spole(7,i)*d2(3,1) + spole(8,i)*d2(4,1) 424 & + spole(9,i)*d2(5,1) 425 srpole(6,i) = spole(5,i)*d2(1,2) + spole(6,i)*d2(2,2) 426 & + spole(7,i)*d2(3,2) + spole(8,i)*d2(4,2) 427 & + spole(9,i)*d2(5,2) 428 srpole(7,i) = spole(5,i)*d2(1,3) + spole(6,i)*d2(2,3) 429 & + spole(7,i)*d2(3,3) + spole(8,i)*d2(4,3) 430 & + spole(9,i)*d2(5,3) 431 srpole(8,i) = spole(5,i)*d2(1,4) + spole(6,i)*d2(2,4) 432 & + spole(7,i)*d2(3,4) + spole(8,i)*d2(4,4) 433 & + spole(9,i)*d2(5,4) 434 srpole(9,i) = spole(5,i)*d2(1,5) + spole(6,i)*d2(2,5) 435 & + spole(7,i)*d2(3,5) + spole(8,i)*d2(4,5) 436 & + spole(9,i)*d2(5,5) 437 return 438 end 439c 440c 441c ############################################################### 442c ## ## 443c ## subroutine qirotmat -- QI interatomic rotation matrix ## 444c ## ## 445c ############################################################### 446c 447c 448c "qirotmat" finds a rotation matrix that describes the 449c interatomic vector 450c 451c 452 subroutine qirotmat (i,k,rinv,d1) 453 use atoms 454 use mpole 455 implicit none 456 integer i,k 457 real*8 rinv,r,dot 458 real*8 dx,dy,dz 459 real*8 a(3,3) 460 real*8 d1(3,3) 461c 462c set the z-axis to be the interatomic vector 463c 464 dx = x(k) - x(i) 465 dy = y(k) - y(i) 466 dz = z(k) - z(i) 467 a(1,3) = dx * rinv 468 a(2,3) = dy * rinv 469 a(3,3) = dz * rinv 470c 471c find an x-axis that is orthogonal to the z-axis 472c 473 if (y(i).ne.y(k) .or. z(i).ne.z(k)) then 474 dx = dx + 1.0d0 475 else 476 dy = dy + 1.0d0 477 end if 478 dot = dx*a(1,3) + dy*a(2,3) + dz*a(3,3) 479 dx = dx - dot*a(1,3) 480 dy = dy - dot*a(2,3) 481 dz = dz - dot*a(3,3) 482 r = 1.d0 / sqrt(dx*dx + dy*dy + dz*dz) 483 a(1,1) = dx * r 484 a(2,1) = dy * r 485 a(3,1) = dz * r 486c 487c get rotation matrix elements for the y-axis 488c 489 a(1,2) = a(3,1)*a(2,3) - a(2,1)*a(3,3) 490 a(2,2) = a(1,1)*a(3,3) - a(3,1)*a(1,3) 491 a(3,2) = a(2,1)*a(1,3) - a(1,1)*a(2,3) 492c 493c reorder to account for spherical harmonics 494c 495 d1(1,1) = a(3,3) 496 d1(2,1) = a(1,3) 497 d1(3,1) = a(2,3) 498 d1(1,2) = a(3,1) 499 d1(2,2) = a(1,1) 500 d1(3,2) = a(2,1) 501 d1(1,3) = a(3,2) 502 d1(2,3) = a(1,2) 503 d1(3,3) = a(2,2) 504 return 505 end 506