Statistics
| Revision:

root / src / CalcRmsd.f90 @ 11

History | View | Annotate | Download (7.2 kB)

1

    
2
! This subroutine calculates RMSD (Root mean square deviation) using quaternions.
3
! It is baNatsed on the F90 routine bu E. Coutsias
4
! http://www.math.unm.edu/~vageli/homepage.html
5
! I (PFL) have just translated it, and I have changed the diagonalization
6
! subroutine.
7
! I also made some changes to make it suitable for Cart package.
8
!----------------------------------------------------------------------
9
!----------------------------------------------------------------------
10
! Copyright (C) 2004, 2005 Chaok Seok, Evangelos Coutsias and Ken Dill
11
!      UCSF, Univeristy of New Mexico, Seoul National University
12
! Witten by Chaok Seok and Evangelos Coutsias 2004.
13
                                                                                
14
! This library is free software; you can redistribute it and/or
15
! modify it under the terms of the GNU Lesser General Public
16
! License as published by the Free Software Foundation; either
17
! version 2.1 of the License, or (at your option) any later version.
18
!
19
                                                                                
20
! This library is distributed in the hope that it will be useful,
21
! but WITHOUT ANY WARRANTY; without even the implied warranty of
22
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
23
! Lesser General Public License for more details.
24
!
25
                                                                                
26
! You should have received a copy of the GNU Lesser General Public
27
! License along with this library; if not, write to the Free Software
28
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
29
!----------------------------------------------------------------------------
30

    
31
      subroutine CalcRmsd(Nat,x0,y0,z0, x2,y2,z2,U,rmsd,FRot,FAlign)
32
!-----------------------------------------------------------------------
33
!  This subroutine calculates the least square rmsd of two coordinate
34
!  sets coord1(3,n) and coord2(3,n) using a method based on quaternion.
35
!  It then calculate the rotation matrix U and the centers of coord, and uses
36
! them to align the two molecules.
37
! x0(Nat), y0(Nat), z0(Nat) are the coordinates of the reference geometry.
38
! x2(Nat), y2(Nat), z2(Nat) are the coordinates of the geometry which is to
39
! be aligned with the reference geometry. 
40
! The rotation matrix U has INTENT(OUT) in subroutine rotation_matrix(...),
41
! which is called in this CalcRmsd(...).
42
! rmsd: the root mean square deviation and calculated in this subroutine
43
! CalcRmsd(...).
44
! FRot: if .TRUE. the rotation matrix is calculated. This is the matrix that rotates
45
! the first molecule onto the second one.
46
! FAlgin: if .TRUE. then the second molecule is aligned on the first one.
47
!-----------------------------------------------------------------------
48

    
49
        use VarTypes
50

    
51
      IMPLICIT NONE
52

    
53
      INTEGER(KINT) :: Nat
54
      real(KREAL) ::  x0(Nat),y0(Nat),z0(Nat)
55
      real(KREAL) ::  x2(Nat),y2(Nat),z2(Nat)
56
      real(KREAL) ::  U(3,3), rmsd
57
      LOGICAL  FRot,FAlign,Debug
58

    
59
      REAL(KREAL),ALLOCATABLE :: Coord1(:,:), Coord2(:,:) ! 3,Nat
60
      real(KREAL) ::  x0c1,y0c1,z0c1, xc2,yc2,zc2
61
      
62

    
63
      INTEGER(KINT) :: i, j, ia
64
      real(KREAL) :: x_norm, y_norm, lambda
65
      real(KREAL) :: Rmatrix(3,3)
66
      real(KREAL) :: S(4, 4)
67
      real(KREAL) :: EigVec(4, 4), EigVal(4)
68

    
69
      INTERFACE
70
       function valid(string) result (isValid)
71
         CHARACTER(*), intent(in) :: string
72
         logical                  :: isValid
73
         END function VALID
74
      END INTERFACE
75

    
76

    
77
      debug=valid('CalcRmsd').OR.valid('align').OR.valid('alignpartial')
78

    
79
      if (debug) Call Header("Entering CalcRmsd")
80

    
81
      if (debug) THEN
82
         WRITE(*,*) "NAt=",NAt
83
         WRITE(*,*) x0
84
         WRITE(*,*) y0
85
         WRITE(*,*) z0
86
      END IF
87

    
88
      ALLOCATE(Coord1(3,Nat),Coord2(3,Nat))
89
  ! calculate the barycenters, centroidal coordinates, and the norms
90
      x_norm = 0.0d0
91
      y_norm = 0.0d0
92
      x0c1=0.
93
      y0c1=0.
94
      z0c1=0.
95
      xc2=0.
96
      yc2=0.
97
      zc2=0.
98
      do ia=1,Nat
99
         if (debug) WRITE(*,'(A,I5,4(1X,F10.4))') 'ia...',ia,x0(ia), &
100
                             x2(ia),x0c1,xc2
101
         x0c1=x0c1+x0(ia)
102
         xc2=xc2+x2(ia)
103
         y0c1=y0c1+y0(ia)
104
         yc2=yc2+y2(ia)
105
         z0c1=z0c1+z0(ia)
106
         zc2=zc2+z2(ia)
107
         if (debug) WRITE(*,'(A,I5,4(1X,F10.4))') 'ia...',ia,x0(ia), &
108
                             x2(ia),x0c1,xc2
109
      END DO
110
      x0c1=x0c1/dble(Nat)
111
      y0c1=y0c1/dble(Nat)
112
      z0c1=z0c1/dble(Nat)
113
      xc2=xc2/dble(Nat)
114
      yc2=yc2/dble(Nat)
115
      zc2=zc2/dble(Nat)
116

    
117
      IF (debug) WRITE(*,'(1X,A,3(1X,F10.4))') 'Center1',x0c1,y0c1,z0c1
118
      IF (debug) WRITE(*,'(1X,A,3(1X,F10.4))') 'Center2',xc2,yc2,zc2
119
      do i=1,Nat
120
         Coord1(1,i)=x0(i)-x0c1
121
         Coord1(2,i)=y0(i)-y0c1
122
         Coord1(3,i)=z0(i)-z0c1
123
         Coord2(1,i)=x2(i)-xc2
124
         Coord2(2,i)=y2(i)-yc2
125
         Coord2(3,i)=z2(i)-zc2
126
         x_norm=x_norm+Coord1(1,i)**2+Coord1(2,i)**2+Coord1(3,i)**2
127
         y_norm=y_norm+Coord2(1,i)**2+Coord2(2,i)**2+Coord2(3,i)**2
128
      end do
129

    
130
      IF (debug) THEN
131
         WRITE(*,*) "R matrix"
132
         DO I=1,3
133
            WRITE(*,*) (RMatrix(I,j),j=1,3)
134
         END DO
135
      END IF
136

    
137
  ! calculate the R matrix
138
      do i = 1, 3
139
         do j = 1, 3
140
            Rmatrix(i,j)=0.
141
            do ia=1,Nat
142
               Rmatrix(i,j) = Rmatrix(i,j)+Coord1(i,ia)*Coord2(j,ia)
143
            END DO
144
         end do
145
      end do
146

    
147
      IF (debug) THEN
148
         WRITE(*,*) "R matrix"
149
         DO I=1,3
150
            WRITE(*,*) (RMatrix(I,j),j=1,3)
151
         END DO
152
      END IF
153

    
154

    
155
  ! S matrix
156
      S(1, 1) = Rmatrix(1, 1) + Rmatrix(2, 2) + Rmatrix(3, 3)
157
      S(2, 1) = Rmatrix(2, 3) - Rmatrix(3, 2)
158
      S(3, 1) = Rmatrix(3, 1) - Rmatrix(1, 3)
159
      S(4, 1) = Rmatrix(1, 2) - Rmatrix(2, 1)
160

    
161
      S(1, 2) = S(2, 1)
162
      S(2, 2) = Rmatrix(1, 1) - Rmatrix(2, 2) - Rmatrix(3, 3)
163
      S(3, 2) = Rmatrix(1, 2) + Rmatrix(2, 1)
164
      S(4, 2) = Rmatrix(1, 3) + Rmatrix(3, 1)
165

    
166
      S(1, 3) = S(3, 1)
167
      S(2, 3) = S(3, 2)
168
      S(3, 3) =-Rmatrix(1, 1) + Rmatrix(2, 2) - Rmatrix(3, 3)
169
      S(4, 3) = Rmatrix(2, 3) + Rmatrix(3, 2)
170

    
171
      S(1, 4) = S(4, 1)
172
      S(2, 4) = S(4, 2)
173
      S(3, 4) = S(4, 3)
174
      S(4, 4) =-Rmatrix(1, 1) - Rmatrix(2, 2) + Rmatrix(3, 3) 
175

    
176

    
177
! PFL : I use my usual Jacobi diagonalisation
178
  ! Calculate eigenvalues and eigenvectors, and 
179
  ! take the maximum eigenvalue lambda and the corresponding eigenvector q.
180

    
181
      IF (debug) THEN
182
         WRITE(*,*) "S matrix"
183
         DO I=1,4
184
            WRITE(*,*) (S(I,j),j=1,4)
185
         END DO
186
      END IF
187

    
188
      Call Jacobi(S,4,EigVal,EigVec,4)
189

    
190
      Call Trie(4,EigVal,EigVec,4)
191

    
192
      Lambda=EigVal(4)
193

    
194
! RMS Deviation
195
       rmsd=sqrt(max(0.0d0,((x_norm+y_norm)-2.0d0*lambda))/dble(Nat))
196
! The subroutine rotation_matrix(...) constructs rotation matrix U as output
197
! from the input quaternion EigVec(1,4).
198
      if (FRot.OR.FAlign) Call rotation_matrix(EigVec(1,4),U)
199
      IF (FAlign) THEN
200
         DO I=1,Nat
201
            x2(i)=Coord2(1,i)*U(1,1)+Coord2(2,i)*U(2,1) &
202
                +Coord2(3,i)*U(3,1) +x0c1
203
            y2(i)=Coord2(1,i)*U(1,2)+Coord2(2,i)*U(2,2) &
204
                +Coord2(3,i)*U(3,2) +y0c1
205
            z2(i)=Coord2(1,i)*U(1,3)+Coord2(2,i)*U(2,3) &
206
                +Coord2(3,i)*U(3,3) +z0c1
207
         END DO
208
      END IF
209

    
210
      DEALLOCATE(Coord1,Coord2)
211

    
212
    END subroutine CalcRmsd