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