Chimie Théorique » Opt'n Path

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

  !-----------------------------------------------------------------------

!----------------------------------------------------------------------

!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon,

!  Centre National de la Recherche Scientifique,

!  Université Claude Bernard Lyon 1. All rights reserved.

!

!  This work is registered with the Agency for the Protection of Programs

!  as IDDN.FR.001.100009.000.S.P.2014.000.30625

!

!  Authors: P. Fleurat-Lessard, P. Dayal

!  Contact: optnpath@gmail.com

!

! This file is part of "Opt'n Path".

!

!  "Opt'n Path" is free software: you can redistribute it and/or modify

!  it under the terms of the GNU Affero General Public License as

!  published by the Free Software Foundation, either version 3 of the License,

!  or (at your option) any later version.

!

!  "Opt'n Path" is distributed in the hope that it will be useful,

!  but WITHOUT ANY WARRANTY; without even the implied warranty of

!

!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

!  GNU Affero General Public License for more details.

!

!  You should have received a copy of the GNU Affero General Public License

!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.

!

! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,

! for commercial licensing opportunities.

!----------------------------------------------------------------------

  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