Statistics
| Revision:

root / src / AlignPartial.f90 @ 10

History | View | Annotate | Download (5.2 kB)

1
subroutine AlignPartial(na, x,y,z,x2,y2,z2,ListAt,U)
2
  !-----------------------------------------------------------------------
3
  !  This subroutine calculates the least square rmsd of two coordinate
4
  !  sets coord1(3,n) and coord2(3,n) using a method based on quaternion.
5
  !  It then calculate the  rotation matrix U and the centers of coord, and uses
6
  ! them to align the two molecules.
7
  ! This new version does not use all atoms to calculate the RMSD and rotation matrix
8
  ! but only those which are 'true' in ListAt
9
  !-----------------------------------------------------------------------
10
  ! Input parameters:
11
  ! Na        : number of atoms
12
  ! x,y,z     : coordinates of the first geometry
13
  ! x2,y2,z2  : coordinates of the second geometry
14
  ! ListAt    : Logical, T if this atom should be considered in the RMSD calc and 
15
  !             alignment.
16
  !
17
  ! output parameters:
18
  ! U         : Real, (3x3) matrix for rotation
19
  ! x2,y2,z2  : coordinates of the second geometry ALIGNED to the first one
20
  !-----------------------------------------------------------------------
21

    
22
  Use VarTypes
23
  IMPLICIT NONE
24

    
25

    
26
  integer(KINT) :: na, NaT
27
  Real(KREAL) ::  x(na),y(Na),z(Na)
28
  Real(KREAL) ::  x2(Na),y2(Na),z2(Na)
29
  Real(KREAL) ::  U(3,3), rmsd
30
  LOGICAL  Debug, ListAt(Na)
31

    
32
  REAL(KREAL) :: Coord1(3,Na), Coord2(3,Na)
33
  Real(KREAL) ::  xc1,yc1,zc1, xc2,yc2,zc2
34

    
35

    
36
  integer(KINT) :: i, j, ia
37
  Real(KREAL) :: x_norm, y_norm, lambda
38
  Real(KREAL) :: Rmatrix(3,3)
39
  Real(KREAL) :: S(4, 4)
40
  Real(KREAL) :: EigVec(4, 4), EigVal(4)
41

    
42
  INTERFACE
43
     function valid(string) result (isValid)
44
       CHARACTER(*), intent(in) :: string
45
       logical                  :: isValid
46
     END function VALID
47
  END INTERFACE
48

    
49

    
50
  debug=valid('CalcRmsd').OR.valid('align').OR.valid('alignpartial')
51

    
52
  ! calculate the barycenters, centroidal coordinates, and the norms
53
  x_norm = 0.0d0
54
  y_norm = 0.0d0
55
  xc1=0.
56
  yc1=0.
57
  zc1=0.
58
  xc2=0.
59
  yc2=0.
60
  zc2=0.
61

    
62
  if (debug) THEN
63
     WRITE(*,*) "DBG AlignPartial X2,Y2,Z2 before align"
64
     DO Ia=1,na
65
        WRITE(*,*) x2(ia),y2(ia),z2(ia)
66
     END DO
67
  END IF
68

    
69
  NAt=0
70
  do ia=1,na
71
     IF (listat(ia)) THEN
72
        Nat=NAt+1
73
        xc1=xc1+x(ia)
74
        xc2=xc2+x2(ia)
75
        yc1=yc1+y(ia)
76
        yc2=yc2+y2(ia)
77
        zc1=zc1+z(ia)
78
        zc2=zc2+z2(ia)
79
        !     if (debug) WRITE(*,'(A,I5,4(1X,F10.4))') 'ia...',ia,x(ia),
80
        !     &                        x2(ia),xc1,xc2
81
     END IF
82
  END DO
83

    
84
  if (Nat.EQ.0) THEN
85
     WRITE(*,*) "ERROR, AlignPartial: ListAt is empty :-)"
86
  END IF
87
  xc1=xc1/dble(nat)
88
  yc1=yc1/dble(nat)
89
  zc1=zc1/dble(nat)
90
  xc2=xc2/dble(nat)
91
  yc2=yc2/dble(nat)
92
  zc2=zc2/dble(nat)
93

    
94
  IF (Debug) WRITE(*,*) 'Found :', Nat, ' atoms'
95
  IF (debug) WRITE(*,*) (ListAt(I),I=1,Na)
96
  IF (debug) WRITE(*,'(1X,A,3(1X,F10.4))') 'Center1',xc1,yc1,zc1
97
  IF (debug) WRITE(*,'(1X,A,3(1X,F10.4))') 'Center2',xc2,yc2,zc2
98
  do i=1,na
99
     Coord1(1,i)=x(i)-xc1
100
     Coord1(2,i)=y(i)-yc1
101
     Coord1(3,i)=z(i)-zc1
102
     Coord2(1,i)=x2(i)-xc2
103
     Coord2(2,i)=y2(i)-yc2
104
     Coord2(3,i)=z2(i)-zc2
105
     IF (ListAt(I)) THEN
106
        x_norm=x_norm+Coord1(1,i)**2+Coord1(2,i)**2+Coord1(3,i)**2
107
        y_norm=y_norm+Coord2(1,i)**2+Coord2(2,i)**2+Coord2(3,i)**2
108
     END IF
109
  end do
110

    
111
  ! calculate the R matrix
112
  do i = 1, 3
113
     do j = 1, 3
114
        Rmatrix(i,j)=0.
115
        do ia=1,na
116
           If (ListAt(ia))  Rmatrix(i,j) = Rmatrix(i,j)+Coord1(i,ia)*Coord2(j,ia)
117
        END DO
118
     end do
119
  end do
120

    
121
  IF (debug) THEN
122
     WRITE(*,*) "R matrix"
123
     DO I=1,3
124
        WRITE(*,*) (RMatrix(I,j),j=1,3)
125
     END DO
126
  END IF
127

    
128

    
129
  ! S matrix
130
  S(1, 1) = Rmatrix(1, 1) + Rmatrix(2, 2) + Rmatrix(3, 3)
131
  S(2, 1) = Rmatrix(2, 3) - Rmatrix(3, 2)
132
  S(3, 1) = Rmatrix(3, 1) - Rmatrix(1, 3)
133
  S(4, 1) = Rmatrix(1, 2) - Rmatrix(2, 1)
134

    
135
  S(1, 2) = S(2, 1)
136
  S(2, 2) = Rmatrix(1, 1) - Rmatrix(2, 2) - Rmatrix(3, 3)
137
  S(3, 2) = Rmatrix(1, 2) + Rmatrix(2, 1)
138
  S(4, 2) = Rmatrix(1, 3) + Rmatrix(3, 1)
139

    
140
  S(1, 3) = S(3, 1)
141
  S(2, 3) = S(3, 2)
142
  S(3, 3) =-Rmatrix(1, 1) + Rmatrix(2, 2) - Rmatrix(3, 3)
143
  S(4, 3) = Rmatrix(2, 3) + Rmatrix(3, 2)
144

    
145
  S(1, 4) = S(4, 1)
146
  S(2, 4) = S(4, 2)
147
  S(3, 4) = S(4, 3)
148
  S(4, 4) =-Rmatrix(1, 1) - Rmatrix(2, 2) + Rmatrix(3, 3) 
149

    
150

    
151
  ! PFL : I use my usual Jacobi diagonalisation
152
  ! Calculate eigenvalues and eigenvectors, and 
153
  ! take the maximum eigenvalue lambda and the corresponding eigenvector q.
154

    
155
  IF (debug) THEN
156
     WRITE(*,*) "S matrix"
157
     DO I=1,4
158
        WRITE(*,*) (S(I,j),j=1,4)
159
     END DO
160
  END IF
161

    
162
  Call Jacobi(S,4,EigVal,EigVec,4)
163

    
164
  Call Trie(4,EigVal,EigVec,4)
165

    
166
  Lambda=EigVal(4)
167

    
168
  ! RMS Deviation
169
  rmsd=sqrt(max(0.0d0,((x_norm+y_norm)-2.0d0*lambda))/dble(nat))
170
  if (debug) WRITE(*,*) "AlignPartial rmsd:",rmsd
171
  Call rotation_matrix(EigVec(1,4),U)
172

    
173

    
174
  IF (debug) THEN
175
     WRITE(*,*) "U matrix"
176
     DO I=1,3
177
        WRITE(*,*) (U(I,j),j=1,3)
178
     END DO
179
  END IF
180

    
181
  DO I=1,na
182
     x2(i)=Coord2(1,i)*U(1,1)+Coord2(2,i)*U(2,1) + Coord2(3,i)*U(3,1) +xc1
183
     y2(i)=Coord2(1,i)*U(1,2)+Coord2(2,i)*U(2,2) + Coord2(3,i)*U(3,2) +yc1
184
     z2(i)=Coord2(1,i)*U(1,3)+Coord2(2,i)*U(2,3) + Coord2(3,i)*U(3,3) +zc1
185
  END DO
186

    
187
  if (debug) THEN
188
     WRITE(*,*) "DBG AlignPartial X2,Y2,Z2 after align"
189
     DO Ia=1,na
190
        WRITE(*,*) x2(ia),y2(ia),z2(ia)
191
     END DO
192
  END IF
193

    
194

    
195
END subroutine AlignPartial