[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [ccp4bb]: problems with
On Thu, 28 Mar 2002, OnLineHelpForm wrote:
> A rather naive question, I'm afraid. I'm trying to determine the inverse
> of a large number of rotation/translation operations for some
> coordinates, each operation represented as a 3x3 matrix and translation vector. I
> understand that if A -> A' by a particular operation, then the rotation
> matrix for A' -> A is the transpose matrix. But what how is the tranlation vector
> of the inverse operation calculated ? It is not simply the reverse
> vector. Also, are there any CCP4 programs which automate this procedure,
> rather than do them all by hand (or specifying INVERT in dm !)
i have written such a utility, i presume others have as well. the secret
is that in determining the inverse, you apply the translation BEFORE doing
the rotation.
cheers,
=======================================================================
"Supernaturalism is not, on the whole, friendly to science and technology,
and there is no use in pretending that it is." - L. Sprague de Camp
=======================================================================
David J. Schuller
modern man in a post-modern world
MacCHESS, Cornell University
djs63@cornell.edu
program sqmatinv
c djms july 1993
c read in SQUASH type 4x4 matrix
c matrix is R11 R12 R13 T1
c R21 R22 R23 T2
c R31 R32 R33 T3
c 0 0 0 1
c invert transformation
c if y = Ax + t
c x = By + u = A^-1(y) + A^-1(-t)
c
real mat(4,4), invmat(4,4)
write(6,*) ' Input 4x4 SQUASH type transformation:'
do i = 1, 4
read(5,*) (mat(i,j),j=1,4)
end do
c
write(6,*)
write(6,*) ' Starting transformation is:'
do i = 1, 4
write(6,'(1x,4F12.6)') (mat(i,j), j=1,4)
end do
c
call trans3inv( mat, invmat)
c
write(6,*)
write(6,*) ' Inverse transformation is:'
do i = 1, 4
write(6,'(1x,4F12.6)') (invmat(i,j), j=1,4)
end do
end
c==========================================================================
subroutine mat3mul(a,b,c)
C**************************************************************************
c 3x3 matrix multiply, c = a * b
c put default values in 4th row and column
c matrix is R11 R12 R13 T1
c R21 R22 R23 T2
c R31 R32 R33 T3
c 0 0 0 1
real a(4,4), b(4,4), c(4,4)
integer i, j, k
c
do k = 1, 3
do i = 1, 3
c(i, k) = 0.0
do j = 1, 3
c(i,k) = c(i,k) + a(i,j) * b(j,k)
end do
end do
c(4,k) = 0.0
c(k,4) = 0.0
end do
c(4,4) = 1.0
return
end
c==========================================================================
subroutine mat3tpose(a, b)
C**************************************************************************
c find transpose of 3x3 matrix
c b is transpose of a
c put default values in 4th row and column
c matrix is R11 R12 R13 T1
c R21 R22 R23 T2
c R31 R32 R33 T3
c 0 0 0 1
real a(4,4), b(4,4)
integer i, j
c
c find 3x3 transpose
do j = 1, 3
do i = 1, 3
b(i, j) = a(j, i)
end do
b(4,j) = 0.0
b(j,4) = 0.0
end do
b(4,4) = 1.0
return
end
c==========================================================================
subroutine trans3inv(a, b)
C**************************************************************************
c find inverse of transformation
c a is 3x3 rotation matrix, at is translation
c b is 3x3 inverse matrix, bt is translation
c y = Ax + at
c x = 1/Ay + (-at)/A
c matrix is R11 R12 R13 T1
c R21 R22 R23 T2
c R31 R32 R33 T3
c 0 0 0 1
c
real a(4,4), b(4,4)
real ct(3)
integer i, j
c
c find inverse b of matrix a
call mat3inv(a, b)
c
c find translation bt = b(-at)
do i = 1, 3
ct(i) = -1.0 * a(i,4)
end do
do i=1,3
b(i,4) = b(i,1) * ct(1) + b(i,2) * ct(2) + b(i,3) * ct(3)
b(4,i) = 0.0
end do
b(4,4) = 1.0
return
end
c==========================================================================
subroutine trans3mul(b, a, c)
C**************************************************************************
c given transformation y = Ax + at and
c transformation z = By + bt
c find z(x) = BAx + Bt + u = Cx + ct
c put defaults in 4th row
c matrix is R11 R12 R13 T1
c R21 R22 R23 T2
c R31 R32 R33 T3
c 0 0 0 1
c
real a(4,4), b(4,4), c(4,4)
real ct(3)
integer i
c
c find BA
call mat3mul(b, a, c)
c
c find ct = B(at) + bt
do i = 1, 3
c(i,4) = b(i,1) * a(1,4) + b(i,2) * a(2,4) +
& b(i,3) * a(3,4) + b(i,4)
c(4,1) = 0.0
end do
c(4,4) = 1.0
c
return
end
c==========================================================================
subroutine cross( xyz1, xyz2, xyz3 )
C**************************************************************************
c find cross product of xyz1, xyz2; return as xyz3
real xyz1(3), xyz2(3), xyz3(3)
xyz3(1) = ( xyz1(2) * xyz2(3) ) - ( xyz1(3) * xyz2(2) )
xyz3(2) = ( xyz1(3) * xyz2(1) ) - ( xyz1(1) * xyz2(3) )
xyz3(3) = ( xyz1(1) * xyz2(2) ) - ( xyz1(2) * xyz2(1) )
return
end
c==========================================================================
real function dot (xyz1, xyz2)
C**************************************************************************
c find dot product of vectors xyz1, xyz2
real xyz1(3), xyz2(3)
dot = (xyz1(1)*xyz2(1)) + (xyz1(2)*xyz2(2)) + (xyz1(3)*xyz2(3))
return
end
c==========================================================================
real function mat3det(a)
C**************************************************************************
c find determinant of 3x3 matrix A
c matrix is R11 R12 R13 T1
c R21 R22 R23 T2
c R31 R32 R33 T3
c 0 0 0 1
real a(4,4)
real row1(3), row2(3), row3(3), temp(3)
integer i
real dot
c find determinant of a = row1 dot row2 X row3
do i = 1, 3
row1(i) = a(1, i)
row2(i) = a(2, i)
row3(i) = a(3, i)
end do
call cross( row2, row3, temp)
mat3det = dot( row1, temp)
return
end
c==========================================================================
subroutine mat3adjoint(a, b)
C**************************************************************************
c find adjoint of 3x3 matrix
c used in inverting matrices
c b = adj(a)
c put defaults in 4th row or column
c matrix is R11 R12 R13 T1
c R21 R22 R23 T2
c R31 R32 R33 T3
c 0 0 0 1
real a(4,4), b(4,4)
real c(3,3)
real mat3minor
integer i, j
do i = 1, 3
do j = 1, 3
c assign 3x3 rotational matrix from input 4x4
c(i,j) = a(i,j)
end do
end do
do i = 1, 3
do j = 1, 3
c calculate each element of adjunct
b(i,j) = mat3minor(c, j, i)
end do
b(4,i) = 0.0
b(i,4) = 0.0
end do
b(4,4) = 1.0
return
end
c==========================================================================
real function mat3minor(a, i, j)
C**************************************************************************
c find minorij(a); is the determinant of the 2x2 with
c row i and col j missing
c rows and columns expanded cyclically, not collapsed;
c no negative signs needed to balance
c put defaults in 4th row or column
c matrix a is R11 R12 R13
c R21 R22 R23
c R31 R32 R33
real a(3,3)
real minor(2,2)
integer i, j, k, l, row, col
do k = 1, 2
do l = 1, 2
row = i + k
if (row .GT. 3) row = row - 3
col = j + l
if (col .GT. 3) col = col - 3
minor(k, l) = a(row, col)
end do
end do
mat3minor = minor(1,1) * minor(2,2) - minor(1,2) * minor(2,1)
return
end
c==========================================================================
subroutine mat3inv(a, b)
C**************************************************************************
c find inverse b of matrix a
c general case: inverse is adjoint(a)/determinant(a)
c special cases:
c if columns of matrix are ORTHOGONAL(dot(c1,c2)= 0, etc)
c and length of each column is 1.0,
c inverse(A) = transpose(A)
c if matrix is upper diagonal(a21, a31, a32 all = 0)
c inverse is upper diagonal
c
c it will probably be quicker to calculate the general case than
c to test for special cases.
c put defaults in 4th row or column
c matrix is R11 R12 R13 T1
c R21 R22 R23 T2
c R31 R32 R33 T3
c 0 0 0 1
c
real a(4,4), b(4,4)
real deta, mat3det
integer i, j
c find determinant of a
deta = mat3det(a)
if ( abs(deta) .LT. 0.00001) then
c if determinant is 0.0, matrix is not invertible
write(6,*) ' Determinant of matrix is 0.0, within error limits'
write(6,*) ' Cannot find inverse'
STOP ' Matrix determinant is zero'
endif
c inverse of a is 1/det(a) * adjoint(a)
call mat3adjoint(a, b)
do i = 1, 3
do j = 1, 3
b(i, j) = b(i, j) / deta
end do
b(i,4) = 0.0
b(4,i) = 0.0
end do
b(4,4) = 1.0
return
end
c==========================================================================
subroutine mat3brk(cx, q, qt)
C**************************************************************************
c convert unit cell dimensions cx( a, b, c, alpha, beta, gamma)
c to orthogonalization matrix q
c v is volume of unit cell
c Brookhaven convention, "Protein Data Bank Atomic Coordinate and
c Bibliographic entry Format Description", appendix A, p. 23
real cx(6)
real a, b, c, alpha, beta, gamma
real v, degtor
integer i, j
real q(4,4), qt(4,4)
c convert angles in degrees to radians
degtor = acos( -1.0) / 180.0
a = cx(1)
b = cx(2)
c = cx(3)
alpha = degtor * cx(4)
beta = degtor * cx(5)
gamma = degtor * cx(6)
c v is unit cell volume
v = a * b * c * sqrt(1 + 2 * cos(alpha) * cos(beta) * cos(gamma) -
& cos(alpha)*cos(alpha) - cos(beta)*cos(beta) -
& cos(gamma) * cos(gamma))
c
q(1, 1) = a
q(2, 1) = 0.0
q(3, 1) = 0.0
q(1, 2) = b * cos(gamma)
q(2, 2) = b * sin(gamma)
q(3, 2) = 0.0
q(1, 3) = c * cos(beta)
q(2,3) = c * (cos(alpha) - cos(beta) * cos(gamma)) /sin(gamma)
q(3, 3) = v / ( a * b * sin(gamma) )
c
do i = 1, 3
q(i,4) = 0.0
q(4,i) = 0.0
end do
q(4,4) = 1.0
c
c get inverse matrix
call mat3inv(q, qt)
c
return
end