http://www.ccl.net/cca/software/SOURCES/FORTRAN/fitest/fitest.shtml
CCL fitest.f
```C
C David J. Heisterberg
C The Ohio Supercomputer Center
C 1224 Kinnear Rd.
C Columbus, OH  43212-1163
C (614)292-6036
C djh@ccl.net    djh@ohstpy.bitnet    ohstpy::djh
C
C ANALYZE
C analyze the x to y fit and the fitting matrix
C
SUBROUTINE analyz (n, x, y, w, u)
IMPLICIT CHARACTER (A-Z)
C
INTEGER n
DOUBLEPRECISION x (3, n)
DOUBLEPRECISION y (3, n)
DOUBLEPRECISION w (n)
DOUBLEPRECISION u (3, 3)
C
DOUBLEPRECISION err
DOUBLEPRECISION wnorm
DOUBLEPRECISION urnorm (3)
DOUBLEPRECISION ucnorm (3)
INTEGER i, j
C
C find the mean sum of squares error
C
err = 0.0D0
wnorm = 0.0D0
DO 10000 i = 1, n
err = err + w (i) * ((y (1, i) - x (1, i)) ** 2 +
&                       (y (2, i) - x (2, i)) ** 2 +
&                       (y (3, i) - x (3, i)) ** 2)
wnorm = wnorm + w (i)
10000   CONTINUE
err = SQRT (err / wnorm)
C
C find the row and column norms of u
C
ucnorm (1) = u (1, 1) ** 2 + u (2, 1) ** 2 + u (3, 1) ** 2
ucnorm (2) = u (1, 2) ** 2 + u (2, 2) ** 2 + u (3, 2) ** 2
ucnorm (3) = u (1, 3) ** 2 + u (2, 3) ** 2 + u (3, 3) ** 2
urnorm (1) = u (1, 1) ** 2 + u (1, 2) ** 2 + u (1, 3) ** 2
urnorm (2) = u (2, 1) ** 2 + u (2, 2) ** 2 + u (2, 3) ** 2
urnorm (3) = u (3, 1) ** 2 + u (3, 2) ** 2 + u (3, 3) ** 2
C
C write the error and u norms
C
WRITE (*, *)
DO 11000 i = 1, 3
WRITE (*, 1000) (u (i, j), j = 1, 3), urnorm (i)
11000   CONTINUE
WRITE (*, *)
WRITE (*, 1000) (ucnorm (j), j = 1, 3), err
C
RETURN
C
1000 FORMAT (1X, 3E18.10, 4X, E18.10)
C
END
C
C CENTER
C center a molecule about its weighted centroid or other origin
C
SUBROUTINE center (n, x, w, io, o, t)
IMPLICIT CHARACTER (A-Z)
C
INTEGER n
DOUBLEPRECISION x (3, n)
DOUBLEPRECISION w (n)
LOGICAL io
DOUBLEPRECISION o (3)
DOUBLEPRECISION t (3)
C
DOUBLEPRECISION wnorm
INTEGER i
C
IF (io) THEN
t (1) = o (1)
t (2) = o (2)
t (3) = o (3)
ELSE
t (1) = 0.0D0
t (2) = 0.0D0
t (3) = 0.0D0
wnorm = 0.0D0
DO 10000 i = 1, n
t (1) = t (1) + x (1, i) * SQRT (w (i))
t (2) = t (2) + x (2, i) * SQRT (w (i))
t (3) = t (3) + x (3, i) * SQRT (w (i))
wnorm = wnorm + SQRT (w (i))
10000     CONTINUE
t (1) = t (1) / wnorm
t (2) = t (2) / wnorm
t (3) = t (3) / wnorm
ENDIF
C
DO 10100 i = 1, n
x (1, i) = x (1, i) - t (1)
x (2, i) = x (2, i) - t (2)
x (3, i) = x (3, i) - t (3)
10100   CONTINUE
C
RETURN
END
C
C FITEST
C rigid fit test driver
C
PROGRAM fitest
IMPLICIT CHARACTER (A-Z)
C
INTEGER n
DOUBLEPRECISION angle
DOUBLEPRECISION pert
PARAMETER (n = 401)
PARAMETER (angle = 1.0D0)
PARAMETER (pert = 0.2D0)
C
DOUBLEPRECISION y (3, n)
DOUBLEPRECISION x (3, n)
DOUBLEPRECISION z (3, n)
DOUBLEPRECISION w (n)
DOUBLEPRECISION o (3)
DOUBLEPRECISION q (4)
DOUBLEPRECISION u (3, 3)
INTEGER ierr
C
C seed the random number generator
C call ranini for vax
C
CALL ranset (123571)
*     CALL ranini ()
C
C generate reference and test molecules and weights
C
CALL genxyw (angle, pert, n, x, y, w)
CALL center (n, y, w, .FALSE., o, o)
CALL center (n, x, w, .FALSE., o, o)
CALL identm (3, 3, u)
CALL analyz (n, x, y, w, u)
C
C fit x to y
C
WRITE (*, *)
WRITE (*, *) 'fitting x to y'
CALL qtrfit (n, x, y, w, q, u, ierr)
CALL rotmol (n, x, z, u)
CALL analyz (n, z, y, w, u)
C
C fit y to x
C
WRITE (*, *)
WRITE (*, *) 'fitting y to x'
CALL qtrfit (n, y, x, w, q, u, ierr)
CALL rotmol (n, y, z, u)
CALL analyz (n, z, x, w, u)
C
END
C
C GENROT
C Generate a left rotation matrix from a normalized rotation axis
C and cosine and sine of the rotation angle.
C
C INPUT
C   a      - rotation axis
C   cost   - cosine of rotation angle
C   sint   - sine of rotation angle
C
C OUTPUT
C   u      - the rotation matrix
C
SUBROUTINE genrot (a, cost, sint, u)
IMPLICIT CHARACTER (A-Z)
C
DOUBLEPRECISION a (3)
DOUBLEPRECISION cost, sint
DOUBLEPRECISION u (3, 3)
C
u (1, 1) = a (1) ** 2 + (1.0D0 - a (1) ** 2) * cost
u (2, 1) = a (1) * a (2) * (1.0D0 - cost) + a (3) * sint
u (3, 1) = a (1) * a (3) * (1.0D0 - cost) - a (2) * sint
C
u (1, 2) = a (1) * a (2) * (1.0D0 - cost) - a (3) * sint
u (2, 2) = a (2) ** 2 + (1.0D0 - a (2) ** 2) * cost
u (3, 2) = a (2) * a (3) * (1.0D0 - cost) + a (1) * sint
C
u (1, 3) = a (1) * a (3) * (1.0D0 - cost) + a (2) * sint
u (2, 3) = a (2) * a (3) * (1.0D0 - cost) - a (1) * sint
u (3, 3) = a (3) ** 2 + (1.0D0 - a (3) ** 2) * cost
C
RETURN
END
C
C GENXYW
C generate random reference and test molecules and weights
C
SUBROUTINE genxyw (angle, pert, n, x, y, w)
IMPLICIT CHARACTER (A-Z)
C
DOUBLEPRECISION angle
DOUBLEPRECISION pert
INTEGER n
DOUBLEPRECISION x (3, n)
DOUBLEPRECISION y (3, n)
DOUBLEPRECISION w (n)
C
DOUBLEPRECISION u (3, 3)
INTEGER i
C
DOUBLEPRECISION random
EXTERNAL random
C
CALL ranmol (n, y)
CALL ranrot (angle, u)
CALL rotmol (n, y, x, u)
CALL pertur (pert, n, x)
CALL ranrot (angle, u)
CALL rotmol (n, x, x, u)
C
DO 10000 i = 1, n
w (i) = random ()
10000   CONTINUE
C
RETURN
END
C
C IDENTM
C generate an identity matrix
C
SUBROUTINE identm (n, np, u)
IMPLICIT CHARACTER (A-Z)
C
INTEGER n, np
DOUBLEPRECISION u (np, n)
C
INTEGER i, j
C
DO 10000 j = 1, n
DO 10010 i = 1, n
u (i, j) = 0.0D0
10010     CONTINUE
u (j, j) = 1.0D0
10000   CONTINUE
C
RETURN
END
C
C JACOBI
C Jacobi diagonalizer with sorted output.  Same calling sequence as
C EISPACK routine, but must specify nrot!
C
SUBROUTINE jacobi (a, n, np, d, v, nrot)
IMPLICIT CHARACTER (A-Z)
C
INTEGER n, np, nrot
DOUBLEPRECISION a (np, n)
DOUBLEPRECISION d (n)
DOUBLEPRECISION v (np, n)
C
DOUBLEPRECISION onorm, dnorm
DOUBLEPRECISION b, dma, q, t, c, s
DOUBLEPRECISION atemp, vtemp, dtemp
INTEGER i, j, k, l
C
DO 10000 j = 1, n
DO 10010 i = 1, n
v (i, j) = 0.0D0
10010     CONTINUE
v (j, j) = 1.0D0
d (j) = a (j, j)
10000   CONTINUE
C
DO 20000 l = 1, nrot
dnorm = 0.0D0
onorm = 0.0D0
DO 20100 j = 1, n
dnorm = dnorm + ABS (d (j))
DO 20110 i = 1, j - 1
onorm = onorm + ABS (a (i, j))
20110       CONTINUE
20100     CONTINUE
IF (onorm / dnorm .LE. 0.0D0) GOTO 19999
DO 21000 j = 2, n
DO 21010 i = 1, j - 1
b = a (i, j)
IF (ABS (b) .GT. 0.0D0) THEN
dma = d (j) - d (i)
IF (ABS (dma) + ABS (b) .LE. ABS (dma)) THEN
t = b / dma
ELSE
q = 0.5D0 * dma / b
t = SIGN (1.0D0 / (ABS (q) + SQRT (1.0D0 + q * q)), q)
ENDIF
c = 1.0D0 / SQRT (t * t + 1.0D0)
s = t * c
a (i, j) = 0.0D0
DO 21110 k = 1, i - 1
atemp    = c * a (k, i) - s * a (k, j)
a (k, j) = s * a (k, i) + c * a (k, j)
a (k, i) = atemp
21110           CONTINUE
DO 21120 k = i + 1, j - 1
atemp    = c * a (i, k) - s * a (k, j)
a (k, j) = s * a (i, k) + c * a (k, j)
a (i, k) = atemp
21120           CONTINUE
DO 21130 k = j + 1, n
atemp    = c * a (i, k) - s * a (j, k)
a (j, k) = s * a (i, k) + c * a (j, k)
a (i, k) = atemp
21130           CONTINUE
DO 21140 k = 1, n
vtemp    = c * v (k, i) - s * v (k, j)
v (k, j) = s * v (k, i) + c * v (k, j)
v (k, i) = vtemp
21140           CONTINUE
dtemp = c * c * d (i) + s * s * d (j) -
&                2.0D0 * c * s * b
d (j) = s * s * d (i) + c * c * d (j) +
&                2.0D0 * c * s * b
d (i) = dtemp
ENDIF
21010       CONTINUE
21000     CONTINUE
20000   CONTINUE
19999 CONTINUE
nrot = l
C
DO 30000 j = 1, n - 1
k = j
dtemp = d (k)
DO 30100 i = j + 1, n
IF (d (i) .LT. dtemp) THEN
k = i
dtemp = d (k)
ENDIF
30100     CONTINUE
IF (k .GT. j) THEN
d (k) = d (j)
d (j) = dtemp
DO 30200 i = 1, n
dtemp    = v (i, k)
v (i, k) = v (i, j)
v (i, j) = dtemp
30200       CONTINUE
ENDIF
30000   CONTINUE
C
RETURN
END
C
C PERTUR
C apply a random perturbation
C
SUBROUTINE pertur (pert, n, x)
IMPLICIT CHARACTER (A-Z)
C
DOUBLEPRECISION pert
INTEGER n
DOUBLEPRECISION x (3, n)
C
INTEGER i, j
C
DOUBLEPRECISION random
EXTERNAL random
C
DO 10000 j = 1, 3
DO 10010 i = 1, n
x (j, i) = x (j, i) *
&               (1.0D0 + 2.0D0 * pert * (0.5D0 - random ()))
10010     CONTINUE
10000   CONTINUE
C
RETURN
END
C
C Q2MAT
C Generate a left rotation matrix from a normalized quaternion
C
C INPUT
C   q      - normalized quaternion
C
C OUTPUT
C   u      - the rotation matrix
C
SUBROUTINE q2mat (q, u)
IMPLICIT NONE
C
DOUBLEPRECISION q (0 : 3)
DOUBLEPRECISION u (3, 3)
C
u (1, 1) = q (0) ** 2 + q (1) ** 2 - q (2) ** 2 - q (3) ** 2
u (2, 1) = 2.0D0 * (q (1) * q (2) - q (0) * q (3))
u (3, 1) = 2.0D0 * (q (1) * q (3) + q (0) * q (2))
C
u (1, 2) = 2.0D0 * (q (2) * q (1) + q (0) * q (3))
u (2, 2) = q (0) ** 2 - q (1) ** 2 + q (2) ** 2 - q (3) ** 2
u (3, 2) = 2.0D0 * (q (2) * q (3) - q (0) * q (1))
C
u (1, 3) = 2.0D0 * (q (3) * q (1) - q (0) * q (2))
u (2, 3) = 2.0D0 * (q (3) * q (2) + q (0) * q (1))
u (3, 3) = q (0) ** 2 - q (1) ** 2 - q (2) ** 2 + q (3) ** 2
C
RETURN
END
C
C QTRFIT
C Find the quaternion, q, [and left rotation matrix, u] that minimizes
C
C   |qTXq - Y| ^ 2    [|uX - Y| ^ 2]
C
C This is equivalent to maximizing Re (qTXTqY).
C
C This is equivalent to finding the largest eigenvalue and corresponding
C eigenvector of the matrix
C
C   [A2   AUx  AUy  AUz ]
C   [AUx  Ux2  UxUy UzUx]
C   [AUy  UxUy Uy2  UyUz]
C   [AUz  UzUx UyUz Uz2 ]
C
C where
C
C   A2   = Xx Yx + Xy Yy + Xz Yz
C   Ux2  = Xx Yx - Xy Yy - Xz Yz
C   Uy2  = Xy Yy - Xz Yz - Xx Yx
C   Uz2  = Xz Yz - Xx Yx - Xy Yy
C   AUx  = Xz Yy - Xy Yz
C   AUy  = Xx Yz - Xz Yx
C   AUz  = Xy Yx - Xx Yy
C   UxUy = Xx Yy + Xy Yx
C   UyUz = Xy Yz + Xz Yy
C   UzUx = Xz Yx + Xx Yz
C
C The left rotation matrix, u, is obtained from q by
C
C   u = qT1q
C
C INPUT
C   n      - number of points
C   x      - test vector
C   y      - reference vector
C   w      - weight vector
C
C OUTPUT
C   q      - the best-fit quaternion
C   u      - the best-fit left rotation matrix
C   nr     - number of jacobi sweeps required
C
SUBROUTINE qtrfit (n, x, y, w, q, u, nr)
IMPLICIT CHARACTER (A-Z)
C
INTEGER n
DOUBLEPRECISION x (3, n)
DOUBLEPRECISION y (3, n)
DOUBLEPRECISION w (n)
DOUBLEPRECISION q (0 : 3)
DOUBLEPRECISION u (3, 3)
INTEGER nr
C
DOUBLEPRECISION xxyx, xxyy, xxyz
DOUBLEPRECISION xyyx, xyyy, xyyz
DOUBLEPRECISION xzyx, xzyy, xzyz
DOUBLEPRECISION c (0 : 3, 0 : 3), v (0 : 3, 0 : 3)
DOUBLEPRECISION d (0 : 3)
INTEGER i
C
C generate the upper triangle of the quadratic form matrix
C
xxyx = 0.0D0
xxyy = 0.0D0
xxyz = 0.0D0
xyyx = 0.0D0
xyyy = 0.0D0
xyyz = 0.0D0
xzyx = 0.0D0
xzyy = 0.0D0
xzyz = 0.0D0
DO 11000 i = 1, n
xxyx = xxyx + x (1, i) * y (1, i) * w (i)
xxyy = xxyy + x (1, i) * y (2, i) * w (i)
xxyz = xxyz + x (1, i) * y (3, i) * w (i)
xyyx = xyyx + x (2, i) * y (1, i) * w (i)
xyyy = xyyy + x (2, i) * y (2, i) * w (i)
xyyz = xyyz + x (2, i) * y (3, i) * w (i)
xzyx = xzyx + x (3, i) * y (1, i) * w (i)
xzyy = xzyy + x (3, i) * y (2, i) * w (i)
xzyz = xzyz + x (3, i) * y (3, i) * w (i)
11000   CONTINUE
C
c (0, 0) = xxyx + xyyy + xzyz
C
c (0, 1) = xzyy - xyyz
c (1, 1) = xxyx - xyyy - xzyz
C
c (0, 2) = xxyz - xzyx
c (1, 2) = xxyy + xyyx
c (2, 2) = xyyy - xzyz - xxyx
C
c (0, 3) = xyyx - xxyy
c (1, 3) = xzyx + xxyz
c (2, 3) = xyyz + xzyy
c (3, 3) = xzyz - xxyx - xyyy
C
C diagonalize c
C
nr = 16
CALL jacobi (c, 4, 4, d, v, nr)
C
C extract the desired quaternion
C
q (0) = v (0, 3)
q (1) = v (1, 3)
q (2) = v (2, 3)
q (3) = v (3, 3)
C
C generate the rotation matrix
C
CALL q2mat (q, u)
C
RETURN
END
C
C RANDOM
C random number generator after Knuth
C
DOUBLEPRECISION FUNCTION random ()
IMPLICIT CHARACTER (A-Z)
C
INTEGER ntable
INTEGER im1, im2, im3
INTEGER ia1, ia2, ia3
INTEGER ic1, ic2, ic3
DOUBLEPRECISION rm1, rm2, rm3
PARAMETER (ntable = 97)
PARAMETER (im1 = 714025, im2 = 214326, im3 = 139968)
PARAMETER (ia1 = 1366,   ia2 = 3613,   ia3 = 3877)
PARAMETER (ic1 = 150889, ic2 = 45289,  ic3 = 29573)
PARAMETER (rm1 = 1.0D0 / im1)
PARAMETER (rm2 = 1.0D0 / im2)
PARAMETER (rm3 = 1.0D0 / im3)
C
INTEGER iseed0, iseed1, iseed2, iseed3
DOUBLEPRECISION r (ntable)
INTEGER i
C
COMMON /rtable/ iseed0, iseed1, iseed2, iseed3, r
SAVE /rtable/
C
iseed1 = MOD (ia1 * iseed1 + ic1, im1)
iseed2 = MOD (ia2 * iseed2 + ic2, im2)
iseed3 = MOD (ia3 * iseed3 + ic3, im3)
C
i = 1 + (ntable * iseed3) * rm3
random = r (i)
r (i) = (iseed1 + iseed2 * rm2) * rm1
C
RETURN
END
INTEGER FUNCTION ranget ()
IMPLICIT CHARACTER (A-Z)
C
INTEGER ntable
INTEGER im1, im2, im3
INTEGER ia1, ia2, ia3
INTEGER ic1, ic2, ic3
DOUBLEPRECISION rm1, rm2, rm3
PARAMETER (ntable = 97)
PARAMETER (im1 = 714025, im2 = 214326, im3 = 139968)
PARAMETER (ia1 = 1366,   ia2 = 3613,   ia3 = 3877)
PARAMETER (ic1 = 150889, ic2 = 45289,  ic3 = 29573)
PARAMETER (rm1 = 1.0D0 / im1)
PARAMETER (rm2 = 1.0D0 / im2)
PARAMETER (rm3 = 1.0D0 / im3)
C
INTEGER iseed0, iseed1, iseed2, iseed3
DOUBLEPRECISION r (ntable)
C
COMMON /rtable/ iseed0, iseed1, iseed2, iseed3, r
SAVE /rtable/
C
ranget = iseed0
C
RETURN
END
C
C RANINI
C seed the random number generator
C
*     SUBROUTINE ranini ()
*     IMPLICIT CHARACTER (A-Z)
C
*     INTEGER bintim (2)
C
*     CALL sys\$gettim (bintim)
*     CALL ranset (IAND (bintim (1), '7FFFFFFF'X))
C
*     RETURN
*     END
C
C RANMOL
C generate a random molecule
C
SUBROUTINE ranmol (n, x)
IMPLICIT CHARACTER (A-Z)
C
INTEGER n
DOUBLEPRECISION x (3, n)
C
INTEGER i, j
C
DOUBLEPRECISION random
EXTERNAL random
C
C use nested loops to get same result for scalar/vector
C
DO 10000 j = 1, 3
DO 10010 i = 1, n
x (j, i) = random ()
10010     CONTINUE
10000   CONTINUE
C
RETURN
END
C
C RANROT
C generate a random (almost) rotation matrix
C
SUBROUTINE ranrot (angle, u)
IMPLICIT CHARACTER (A-Z)
C
DOUBLEPRECISION angle
DOUBLEPRECISION u (3, 3)
C
DOUBLEPRECISION a (3)
DOUBLEPRECISION anorm
DOUBLEPRECISION theta
DOUBLEPRECISION pi
C
DOUBLEPRECISION random
EXTERNAL random
C
pi = 4.0D0 * ATAN (1.0D0)
C
a (1) = 0.5D0 - random ()
a (2) = 0.5D0 - random ()
a (3) = 0.5D0 - random ()
C
anorm = SQRT (a (1) ** 2 + a (2) ** 2 + a (3) ** 2)
a (1) = a (1) / anorm
a (2) = a (2) / anorm
a (3) = a (3) / anorM
C
theta = angle * pi * random ()
C
CALL genrot (a, COS (theta), SIN (theta), u)
C
RETURN
END
SUBROUTINE ranset (iseed)
IMPLICIT CHARACTER (A-Z)
C
INTEGER ntable
INTEGER im1, im2, im3
INTEGER ia1, ia2, ia3
INTEGER ic1, ic2, ic3
DOUBLEPRECISION rm1, rm2, rm3
PARAMETER (ntable = 97)
PARAMETER (im1 = 714025, im2 = 214326, im3 = 139968)
PARAMETER (ia1 = 1366,   ia2 = 3613,   ia3 = 3877)
PARAMETER (ic1 = 150889, ic2 = 45289,  ic3 = 29573)
PARAMETER (rm1 = 1.0D0 / im1)
PARAMETER (rm2 = 1.0D0 / im2)
PARAMETER (rm3 = 1.0D0 / im3)
C
INTEGER iseed
C
INTEGER iseed0, iseed1, iseed2, iseed3
DOUBLEPRECISION r (ntable)
INTEGER i
C
COMMON /rtable/ iseed0, iseed1, iseed2, iseed3, r
SAVE /rtable/
C
iseed0 = iseed
iseed1 = MOD (iseed0, im1)
iseed2 = MOD (iseed0, im2)
iseed3 = MOD (iseed0, im3)
C
DO 10000 i = 1, ntable
iseed1 = MOD (ia1 * iseed1 + ic1, im1)
iseed2 = MOD (ia2 * iseed2 + ic2, im2)
iseed3 = MOD (ia3 * iseed3 + ic3, im3)
r (i) = (iseed1 + iseed2 * rm2) * rm1
10000   CONTINUE
C
RETURN
END
C
C ROTMOL
C rotate a molecule
C
SUBROUTINE rotmol (n, x, y, u)
IMPLICIT CHARACTER (A-Z)
C
INTEGER n
DOUBLEPRECISION x (3, n)
DOUBLEPRECISION y (3, n)
DOUBLEPRECISION u (3, 3)
C
DOUBLEPRECISION yx, yy, yz
INTEGER i
C
DO 10000 i = 1, n
yx = u(1, 1) * x(1, i) + u(1, 2) * x(2, i) + u(1, 3) * x(3, i)
yy = u(2, 1) * x(1, i) + u(2, 2) * x(2, i) + u(2, 3) * x(3, i)
yz = u(3, 1) * x(1, i) + u(3, 2) * x(2, i) + u(3, 3) * x(3, i)
C
y (1, i) = yx
y (2, i) = yy
y (3, i) = yz
10000   CONTINUE
C
RETURN
END
```
 Modified: Wed Aug 21 16:00:00 1991 GMT Page accessed 9847 times since Sat Apr 17 21:34:37 1999 GMT