1c 2c 3c ############################################################## 4c ## COPYRIGHT (C) 1997 by Rohit Pappu & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ############################################################## 7c 8c ############################################################### 9c ## ## 10c ## subroutine orient -- rigid body reference coordinates ## 11c ## ## 12c ############################################################### 13c 14c 15c "orient" computes a set of reference Cartesian coordinates 16c in standard orientation for each rigid body atom group 17c 18c 19 subroutine orient 20 use atoms 21 use group 22 use rigid 23 implicit none 24 integer i,j,k 25 integer init,stop 26 real*8 xcm,ycm,zcm 27 real*8 phi,theta,psi 28 real*8 xterm,yterm,zterm 29 real*8 cphi,ctheta,cpsi 30 real*8 sphi,stheta,spsi 31 real*8 a(3,3) 32c 33c 34c perform dynamic allocation of some global arrays 35c 36 if (.not. allocated(xrb)) allocate (xrb(n)) 37 if (.not. allocated(yrb)) allocate (yrb(n)) 38 if (.not. allocated(zrb)) allocate (zrb(n)) 39 if (.not. allocated(rbc)) allocate (rbc(6,ngrp)) 40c 41c use current coordinates as default reference coordinates 42c 43 do i = 1, n 44 xrb(i) = x(i) 45 yrb(i) = y(i) 46 zrb(i) = z(i) 47 end do 48c 49c compute the rigid body coordinates for each atom group 50c 51 call xyzrigid 52c 53c get the center of mass and Euler angles for each group 54c 55 do i = 1, ngrp 56 xcm = rbc(1,i) 57 ycm = rbc(2,i) 58 zcm = rbc(3,i) 59 phi = rbc(4,i) 60 theta = rbc(5,i) 61 psi = rbc(6,i) 62 cphi = cos(phi) 63 sphi = sin(phi) 64 ctheta = cos(theta) 65 stheta = sin(theta) 66 cpsi = cos(psi) 67 spsi = sin(psi) 68c 69c construct the rotation matrix from Euler angle values 70c 71 a(1,1) = ctheta * cphi 72 a(2,1) = spsi*stheta*cphi - cpsi*sphi 73 a(3,1) = cpsi*stheta*cphi + spsi*sphi 74 a(1,2) = ctheta * sphi 75 a(2,2) = spsi*stheta*sphi + cpsi*cphi 76 a(3,2) = cpsi*stheta*sphi - spsi*cphi 77 a(1,3) = -stheta 78 a(2,3) = ctheta * spsi 79 a(3,3) = ctheta * cpsi 80c 81c translate and rotate each atom group into inertial frame 82c 83 init = igrp(1,i) 84 stop = igrp(2,i) 85 do j = init, stop 86 k = kgrp(j) 87 xterm = x(k) - xcm 88 yterm = y(k) - ycm 89 zterm = z(k) - zcm 90 xrb(k) = a(1,1)*xterm + a(1,2)*yterm + a(1,3)*zterm 91 yrb(k) = a(2,1)*xterm + a(2,2)*yterm + a(2,3)*zterm 92 zrb(k) = a(3,1)*xterm + a(3,2)*yterm + a(3,3)*zterm 93 end do 94 end do 95 return 96 end 97c 98c 99c ################################################################# 100c ## ## 101c ## subroutine xyzrigid -- determine rigid body coordinates ## 102c ## ## 103c ################################################################# 104c 105c 106c "xyzrigid" computes the center of mass and Euler angle rigid 107c body coordinates for each atom group in the system 108c 109c literature reference: 110c 111c Herbert Goldstein, "Classical Mechanics, 2nd Edition", 112c Addison-Wesley, Reading, MA, 1980; see the Euler angle 113c xyz convention in Appendix B 114c 115c 116 subroutine xyzrigid 117 use atoms 118 use atomid 119 use group 120 use rigid 121 implicit none 122 integer i,j,k,m 123 integer init,stop 124 real*8 xcm,ycm,zcm 125 real*8 phi,theta,psi 126 real*8 weigh,total,dot 127 real*8 xx,xy,xz,yy,yz,zz 128 real*8 xterm,yterm,zterm 129 real*8 moment(3) 130 real*8 vec(3,3) 131 real*8 tensor(3,3) 132 real*8 a(3,3) 133c 134c 135c get the first and last atom in the current group 136c 137 do i = 1, ngrp 138 init = igrp(1,i) 139 stop = igrp(2,i) 140c 141c compute the position of the group center of mass 142c 143 total = 0.0d0 144 xcm = 0.0d0 145 ycm = 0.0d0 146 zcm = 0.0d0 147 do j = init, stop 148 k = kgrp(j) 149 weigh = mass(k) 150 total = total + weigh 151 xcm = xcm + x(k)*weigh 152 ycm = ycm + y(k)*weigh 153 zcm = zcm + z(k)*weigh 154 end do 155 if (total .ne. 0.0d0) then 156 xcm = xcm / total 157 ycm = ycm / total 158 zcm = zcm / total 159 end if 160c 161c compute and then diagonalize the inertia tensor 162c 163 xx = 0.0d0 164 xy = 0.0d0 165 xz = 0.0d0 166 yy = 0.0d0 167 yz = 0.0d0 168 zz = 0.0d0 169 do j = init, stop 170 k = kgrp(j) 171 weigh = mass(k) 172 xterm = x(k) - xcm 173 yterm = y(k) - ycm 174 zterm = z(k) - zcm 175 xx = xx + xterm*xterm*weigh 176 xy = xy + xterm*yterm*weigh 177 xz = xz + xterm*zterm*weigh 178 yy = yy + yterm*yterm*weigh 179 yz = yz + yterm*zterm*weigh 180 zz = zz + zterm*zterm*weigh 181 end do 182 tensor(1,1) = yy + zz 183 tensor(2,1) = -xy 184 tensor(3,1) = -xz 185 tensor(1,2) = -xy 186 tensor(2,2) = xx + zz 187 tensor(3,2) = -yz 188 tensor(1,3) = -xz 189 tensor(2,3) = -yz 190 tensor(3,3) = xx + yy 191 call jacobi (3,tensor,moment,vec) 192c 193c select the direction for each principle moment axis 194c 195 do m = 1, 2 196 do j = init, stop 197 k = kgrp(j) 198 xterm = vec(1,m) * (x(k)-xcm) 199 yterm = vec(2,m) * (y(k)-ycm) 200 zterm = vec(3,m) * (z(k)-zcm) 201 dot = xterm + yterm + zterm 202 if (dot .lt. 0.0d0) then 203 vec(1,m) = -vec(1,m) 204 vec(2,m) = -vec(2,m) 205 vec(3,m) = -vec(3,m) 206 end if 207 if (dot .ne. 0.0d0) goto 10 208 end do 209 10 continue 210 end do 211c 212c moment axes must give a right-handed coordinate system 213c 214 xterm = vec(1,1) * (vec(2,2)*vec(3,3)-vec(2,3)*vec(3,2)) 215 yterm = vec(2,1) * (vec(1,3)*vec(3,2)-vec(1,2)*vec(3,3)) 216 zterm = vec(3,1) * (vec(1,2)*vec(2,3)-vec(1,3)*vec(2,2)) 217 dot = xterm + yterm + zterm 218 if (dot .lt. 0.0d0) then 219 do j = 1, 3 220 vec(j,3) = -vec(j,3) 221 end do 222 end if 223c 224c principal moment axes form rows of Euler rotation matrix 225c 226 do k = 1, 3 227 do j = 1, 3 228 a(k,j) = vec(j,k) 229 end do 230 end do 231c 232c compute Euler angles consistent with the rotation matrix 233c 234 call roteuler (a,phi,theta,psi) 235c 236c set the rigid body coordinates for each atom group 237c 238 rbc(1,i) = xcm 239 rbc(2,i) = ycm 240 rbc(3,i) = zcm 241 rbc(4,i) = phi 242 rbc(5,i) = theta 243 rbc(6,i) = psi 244 end do 245 return 246 end 247c 248c 249c ################################################################# 250c ## ## 251c ## subroutine roteuler -- rotation matrix to Euler angles ## 252c ## ## 253c ################################################################# 254c 255c 256c "roteuler" computes a set of Euler angle values consistent 257c with an input rotation matrix 258c 259c 260 subroutine roteuler (a,phi,theta,psi) 261 use math 262 implicit none 263 integer i 264 real*8 phi,theta,psi,eps 265 real*8 cphi,ctheta,cpsi 266 real*8 sphi,stheta,spsi 267 real*8 a(3,3),b(3) 268 logical flip(3) 269c 270c 271c set the tolerance for Euler angles and rotation elements 272c 273 eps = 1.0d-7 274c 275c get a trial value of theta from a single rotation element 276c 277 theta = asin(min(1.0d0,max(-1.0d0,-a(1,3)))) 278 ctheta = cos(theta) 279 stheta = -a(1,3) 280c 281c set the phi/psi difference when theta is either 90 or -90 282c 283 if (abs(ctheta) .le. eps) then 284 phi = 0.0d0 285 if (abs(a(3,1)) .lt. eps) then 286 psi = asin(min(1.0d0,max(-1.0d0,-a(2,1)/a(1,3)))) 287 else if (abs(a(2,1)) .lt. eps) then 288 psi = acos(min(1.0d0,max(-1.0d0,-a(3,1)/a(1,3)))) 289 else 290 psi = atan(a(2,1)/a(3,1)) 291 end if 292c 293c set the phi and psi values for all other theta values 294c 295 else 296 if (abs(a(1,1)) .lt. eps) then 297 phi = asin(min(1.0d0,max(-1.0d0,a(1,2)/ctheta))) 298 else if (abs(a(1,2)) .lt. eps) then 299 phi = acos(min(1.0d0,max(-1.0d0,a(1,1)/ctheta))) 300 else 301 phi = atan(a(1,2)/a(1,1)) 302 end if 303 if (abs(a(3,3)) .lt. eps) then 304 psi = asin(min(1.0d0,max(-1.0d0,a(2,3)/ctheta))) 305 else if (abs(a(2,3)) .lt. eps) then 306 psi = acos(min(1.0d0,max(-1.0d0,a(3,3)/ctheta))) 307 else 308 psi = atan(a(2,3)/a(3,3)) 309 end if 310 end if 311c 312c find sine and cosine of the trial phi and psi values 313c 314 cphi = cos(phi) 315 sphi = sin(phi) 316 cpsi = cos(psi) 317 spsi = sin(psi) 318c 319c reconstruct the diagonal of the rotation matrix 320c 321 b(1) = ctheta * cphi 322 b(2) = spsi*stheta*sphi + cpsi*cphi 323 b(3) = ctheta * cpsi 324c 325c compare the correct matrix diagonal to rebuilt diagonal 326c 327 do i = 1, 3 328 flip(i) = .false. 329 if (abs(a(i,i)-b(i)) .gt. eps) flip(i) = .true. 330 end do 331c 332c alter Euler angles to get correct rotation matrix values 333c 334 if (flip(1) .and. flip(2)) phi = phi - sign(pi,phi) 335 if (flip(1) .and. flip(3)) theta = -theta + sign(pi,theta) 336 if (flip(2) .and. flip(3)) psi = psi - sign(pi,psi) 337c 338c convert maximum negative angles to positive values 339c 340 if (phi .le. -pi) phi = pi 341 if (theta .le. -pi) theta = pi 342 if (psi .le. -pi) psi = pi 343 return 344 end 345c 346c 347c ############################################################### 348c ## ## 349c ## subroutine rigidxyz -- rigid body to Cartesian coords ## 350c ## ## 351c ############################################################### 352c 353c 354c "rigidxyz" computes Cartesian coordinates for a rigid body 355c group via rotation and translation of reference coordinates 356c 357c literature reference: 358c 359c Herbert Goldstein, "Classical Mechanics, 2nd Edition", 360c Addison-Wesley, Reading, MA, 1980; see the Euler angle 361c xyz convention in Appendix B 362c 363c 364 subroutine rigidxyz 365 use atoms 366 use group 367 use rigid 368 implicit none 369 integer i,j,k 370 integer init,stop 371 real*8 xcm,ycm,zcm 372 real*8 phi,theta,psi 373 real*8 xterm,yterm,zterm 374 real*8 cphi,ctheta,cpsi 375 real*8 sphi,stheta,spsi 376 real*8 a(3,3) 377c 378c 379c get the center of mass and Euler angles for each group 380c 381 do i = 1, ngrp 382 xcm = rbc(1,i) 383 ycm = rbc(2,i) 384 zcm = rbc(3,i) 385 phi = rbc(4,i) 386 theta = rbc(5,i) 387 psi = rbc(6,i) 388 cphi = cos(phi) 389 sphi = sin(phi) 390 ctheta = cos(theta) 391 stheta = sin(theta) 392 cpsi = cos(psi) 393 spsi = sin(psi) 394c 395c construct the rotation matrix from Euler angle values 396c 397 a(1,1) = ctheta * cphi 398 a(2,1) = spsi*stheta*cphi - cpsi*sphi 399 a(3,1) = cpsi*stheta*cphi + spsi*sphi 400 a(1,2) = ctheta * sphi 401 a(2,2) = spsi*stheta*sphi + cpsi*cphi 402 a(3,2) = cpsi*stheta*sphi - spsi*cphi 403 a(1,3) = -stheta 404 a(2,3) = ctheta * spsi 405 a(3,3) = ctheta * cpsi 406c 407c rotate and translate reference coordinates into global frame 408c 409 init = igrp(1,i) 410 stop = igrp(2,i) 411 do j = init, stop 412 k = kgrp(j) 413 xterm = xrb(k) 414 yterm = yrb(k) 415 zterm = zrb(k) 416 x(k) = a(1,1)*xterm + a(2,1)*yterm + a(3,1)*zterm + xcm 417 y(k) = a(1,2)*xterm + a(2,2)*yterm + a(3,2)*zterm + ycm 418 z(k) = a(1,3)*xterm + a(2,3)*yterm + a(3,3)*zterm + zcm 419 end do 420 end do 421 return 422 end 423