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)   !-----------------------------------------------------------------------   ! This subroutine calculates the least square rmsd of two coordinate   ! sets coord1(3,n) and coord2(3,n) using a method based on quaternion.   ! It then calculate the rotation matrix U and the centers of coord, and uses   ! them to align the two molecules.   ! This new version does not use all atoms to calculate the RMSD and rotation matrix   ! but only those which are 'true' in ListAt   !-----------------------------------------------------------------------   ! Input parameters:   ! Na : number of atoms   ! x,y,z : coordinates of the first geometry   ! x2,y2,z2 : coordinates of the second geometry   ! ListAt : Logical, T if this atom should be considered in the RMSD calc and   ! alignment.   !   ! output parameters:   ! U : Real, (3x3) matrix for rotation   ! x2,y2,z2 : coordinates of the second geometry ALIGNED to the first one   !-----------------------------------------------------------------------   Use VarTypes   IMPLICIT NONE   integer(KINT) :: na, NaT   Real(KREAL) :: x(na),y(Na),z(Na)   Real(KREAL) :: x2(Na),y2(Na),z2(Na)   Real(KREAL) :: U(3,3), rmsd   LOGICAL Debug, ListAt(Na)   REAL(KREAL) :: Coord1(3,Na), Coord2(3,Na)   Real(KREAL) :: xc1,yc1,zc1, xc2,yc2,zc2   integer(KINT) :: i, j, ia   Real(KREAL) :: x_norm, y_norm, lambda   Real(KREAL) :: Rmatrix(3,3)   Real(KREAL) :: S(4, 4)   Real(KREAL) :: EigVec(4, 4), EigVal(4)   INTERFACE   function valid(string) result (isValid)   CHARACTER(*), intent(in) :: string   logical :: isValid   END function VALID   END INTERFACE   debug=valid('CalcRmsd').OR.valid('align').OR.valid('alignpartial')   ! calculate the barycenters, centroidal coordinates, and the norms   x_norm = 0.0d0   y_norm = 0.0d0   xc1=0.   yc1=0.   zc1=0.   xc2=0.   yc2=0.   zc2=0.   if (debug) THEN   WRITE(*,*) "DBG AlignPartial X2,Y2,Z2 before align"   DO Ia=1,na   WRITE(*,*) x2(ia),y2(ia),z2(ia)   END DO   END IF   NAt=0   do ia=1,na   IF (listat(ia)) THEN   Nat=NAt+1   xc1=xc1+x(ia)   xc2=xc2+x2(ia)   yc1=yc1+y(ia)   yc2=yc2+y2(ia)   zc1=zc1+z(ia)   zc2=zc2+z2(ia)   ! if (debug) WRITE(*,'(A,I5,4(1X,F10.4))') 'ia...',ia,x(ia),   ! & x2(ia),xc1,xc2   END IF   END DO   if (Nat.EQ.0) THEN   WRITE(*,*) "ERROR, AlignPartial: ListAt is empty :-)"   END IF   xc1=xc1/dble(nat)   yc1=yc1/dble(nat)   zc1=zc1/dble(nat)   xc2=xc2/dble(nat)   yc2=yc2/dble(nat)   zc2=zc2/dble(nat)   IF (Debug) WRITE(*,*) 'Found :', Nat, ' atoms'   IF (debug) WRITE(*,*) (ListAt(I),I=1,Na)   IF (debug) WRITE(*,'(1X,A,3(1X,F10.4))') 'Center1',xc1,yc1,zc1   IF (debug) WRITE(*,'(1X,A,3(1X,F10.4))') 'Center2',xc2,yc2,zc2   do i=1,na   Coord1(1,i)=x(i)-xc1   Coord1(2,i)=y(i)-yc1   Coord1(3,i)=z(i)-zc1   Coord2(1,i)=x2(i)-xc2   Coord2(2,i)=y2(i)-yc2   Coord2(3,i)=z2(i)-zc2   IF (ListAt(I)) THEN   x_norm=x_norm+Coord1(1,i)**2+Coord1(2,i)**2+Coord1(3,i)**2   y_norm=y_norm+Coord2(1,i)**2+Coord2(2,i)**2+Coord2(3,i)**2   END IF   end do   ! calculate the R matrix   do i = 1, 3   do j = 1, 3   Rmatrix(i,j)=0.   do ia=1,na   If (ListAt(ia)) Rmatrix(i,j) = Rmatrix(i,j)+Coord1(i,ia)*Coord2(j,ia)   END DO   end do   end do   IF (debug) THEN   WRITE(*,*) "R matrix"   DO I=1,3   WRITE(*,*) (RMatrix(I,j),j=1,3)   END DO   END IF   ! S matrix   S(1, 1) = Rmatrix(1, 1) + Rmatrix(2, 2) + Rmatrix(3, 3)   S(2, 1) = Rmatrix(2, 3) - Rmatrix(3, 2)   S(3, 1) = Rmatrix(3, 1) - Rmatrix(1, 3)   S(4, 1) = Rmatrix(1, 2) - Rmatrix(2, 1)   S(1, 2) = S(2, 1)   S(2, 2) = Rmatrix(1, 1) - Rmatrix(2, 2) - Rmatrix(3, 3)   S(3, 2) = Rmatrix(1, 2) + Rmatrix(2, 1)   S(4, 2) = Rmatrix(1, 3) + Rmatrix(3, 1)   S(1, 3) = S(3, 1)   S(2, 3) = S(3, 2)   S(3, 3) =-Rmatrix(1, 1) + Rmatrix(2, 2) - Rmatrix(3, 3)   S(4, 3) = Rmatrix(2, 3) + Rmatrix(3, 2)   S(1, 4) = S(4, 1)   S(2, 4) = S(4, 2)   S(3, 4) = S(4, 3)   S(4, 4) =-Rmatrix(1, 1) - Rmatrix(2, 2) + Rmatrix(3, 3)   ! PFL : I use my usual Jacobi diagonalisation   ! Calculate eigenvalues and eigenvectors, and   ! take the maximum eigenvalue lambda and the corresponding eigenvector q.   IF (debug) THEN   WRITE(*,*) "S matrix"   DO I=1,4   WRITE(*,*) (S(I,j),j=1,4)   END DO   END IF   Call Jacobi(S,4,EigVal,EigVec,4)   Call Trie(4,EigVal,EigVec,4)   Lambda=EigVal(4)   ! RMS Deviation   rmsd=sqrt(max(0.0d0,((x_norm+y_norm)-2.0d0*lambda))/dble(nat))   if (debug) WRITE(*,*) "AlignPartial rmsd:",rmsd   Call rotation_matrix(EigVec(1,4),U)   IF (debug) THEN   WRITE(*,*) "U matrix"   DO I=1,3   WRITE(*,*) (U(I,j),j=1,3)   END DO   END IF   DO I=1,na   x2(i)=Coord2(1,i)*U(1,1)+Coord2(2,i)*U(2,1) + Coord2(3,i)*U(3,1) +xc1   y2(i)=Coord2(1,i)*U(1,2)+Coord2(2,i)*U(2,2) + Coord2(3,i)*U(3,2) +yc1   z2(i)=Coord2(1,i)*U(1,3)+Coord2(2,i)*U(2,3) + Coord2(3,i)*U(3,3) +zc1   END DO   if (debug) THEN   WRITE(*,*) "DBG AlignPartial X2,Y2,Z2 after align"   DO Ia=1,na   WRITE(*,*) x2(ia),y2(ia),z2(ia)   END DO   END IF  END subroutine AlignPartial