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
