subroutine psi_rot(amino,coord,dpsi,CAatomnumber,numberofatoms)
!      subroutine psi_rot(amino,ca,cb,cp,dpsi,n,o)
! Performs a right-handed rotation by angle dpsi about
! the Ca-C' axis of the psi angle of the amino-th amino acid.
! This should rotate all atoms in subsequent residues.
! Thus this rotation does not apply to the last residue.
      use defs; implicit none
      integer(i4b), intent(in) :: amino
      real(dp), dimension(:,:), intent(inout) :: coord
      real(dp), intent(in) :: dpsi
      integer(i4b), dimension(:), intent(in) :: CAatomnumber
      integer(i4b), intent(in) :: numberofatoms
      real(dp), dimension(3,3) :: A,C,Id,Rot
      real(dp), dimension(3) :: axis,trans,Cp,CA
      integer(i4b) :: i,j
! The rule is:      
!  Rot = A*sin(theta) + (Id-C)*cos(theta) + C
! where A(i,j) = -axis * epsilon(k,i,j)
! and C(i,j) = axis(i)*axis(j) 
! in which axis is a unit vector.

      Cp = coord(:,CAatomnumber(amino)+1)
      CA = coord(:,CAatomnumber(amino))
      axis = Cp - CA
!     axis = cp(:,amino) - ca(:,amino)
      axis = axis/sqrt(dot_product(axis,axis))
      
      Id = 0.0_dp
      Id(1,1) = 1.0_dp; Id(2,2) = 1.0_dp; Id(3,3) = 1.0_dp
      
      A = 0.0_dp
      A(1,2) = - axis(3); A(1,3) = axis(2)
      A(2,1) = axis(3); A(2,3) = - axis(1)
      A(3,1) = - axis(2); A(3,2) = axis(1)

      do j = 1, 3
      do i = 1, 3
        C(i,j) = axis(i)*axis(j)
      end do
      end do
      
      Rot = A*sin(dpsi) + (Id - C)*cos(dpsi) + C
! Compensate for the displacement of the origin:     
! A' - N = Rot( A - N )  so A' = Rot(A) + N - Rot(N)  
      trans = CA - matmul(Rot,CA)    
!     trans = ca(:,amino) - matmul(Rot,ca(:,amino))
 do i = CAatomnumber(amino)+2, numberofatoms
   coord(:,i) = matmul(Rot,coord(:,i)) + trans
 end do
end subroutine psi_rot