subroutine distance(nres,CA,pdbCA,rmsdistance)
! This code computes the mean rms distance between
! the alpha carbons CA of the folded protein
! and those pdbCA of the pdb-file.  After some algebra
! the code calls DGESVD (not DSYEV) which finds the
! singular values of the matrix A. 
 use defs; implicit none 
 integer(i4b), intent(in) :: nres
 real(dp), dimension(:,:), intent(in) :: CA,pdbCA
 real(dp), intent(out) :: rmsdistance
 
 real(dp), dimension(3) :: P, V
 real(dp), dimension(:,:), allocatable :: r,q
 real(dp) :: nd2_cnst,y
 integer(i4b) :: i,j,k,LWORKnew
 character(len=1) :: JOBU, JOBVT
 integer(i4b) :: M,N,LDA,LDU,LDVT,INFO
 real(dp), dimension(3,3) :: A
 real(dp), dimension(3) :: S
 real(dp), dimension(3,3) :: U,VT
 integer(i4b), parameter :: LWORK = 393
 real(dp), dimension(LWORK) :: WORK 

 allocate(r(3,nres),q(3,nres))
 A = 0.0_dp

! The idea is to minimize the quantity
! nres*rmsdistance**2 = sum( r*r + q*q -2q*R*r ) 
!                   = sum(r*r + q*q) -2Tr(sqrt(A*A^T))
! where A(i,j) = sum_k r_i^k q_j^k
! The mean rms distance d is then
! d = sqrt{ sum(r*r + q*q) -2Tr[sqrt(A*A^T)]/nres }.

! First the center V of the folding protein
! and that P of the pdb protein:

 V = 0.0_dp; P = 0.0_dp
 do i = 1, nres
 V = V + CA(:,i)
 P = P + pdbCA(:,i)
 end do
 V = V/real(nres,dp); P = P/real(nres,dp)

! Now make the relative coordinates:
 do i = 1, nres
 r(:,i) = CA(:,i) - V; q(:,i) = pdbCA(:,i) - P
 end do

! This is the constant part of nres*rmsdistance**2
 nd2_cnst = 0.0_dp
 do i = 1, nres
 nd2_cnst = nd2_cnst + dot_product(r(:,i),r(:,i)) &
            + dot_product(q(:,i),q(:,i))
 end do

! Now make the matrix A(i,j) = sum_k r_k(i) * q_k(j)
 do j = 1, 3; do i = 1, 3; do k = 1, nres
 A(i,j) = A(i,j) + r(i,k)*q(j,k)
 enddo; enddo; enddo

! Stuff for DGESVD  A <=> A 
 JOBU = 'N'; JOBVT = 'N'; M = 3; N = 3; LDA = 3
 LDU = 3; LDVT = 3
 call DGESVD ( JOBU, JOBVT, M, N, A, LDA, S, U, &
               LDU, VT, LDVT, WORK, LWORK, INFO )

 LWORKnew = WORK(1)
! The singular values are the elements of S
! Now make y = trace(S) = Tr[sqrt(A*A^T)]
 y = S(1) + S(2) + S(3)

 rmsdistance = nd2_cnst - 2.0_dp*y
 rmsdistance = max(rmsdistance/real(nres,dp),0.0_dp)
 rmsdistance = sqrt(rmsdistance)

 deallocate(r,q)
 if ( INFO /= 0 ) then
  write(6,*) 'INFO =', INFO
  stop 
 else if ( LWORKnew /= LWORK ) then
  write(6,*) 'LWORK should be =', LWORKnew,', not',LWORK
  stop   
 end if
 end subroutine distance