subroutine phi_rot(amino,coord,dphi,CAatomnumber,  &
                         numberofatoms)
!     subroutine phi_rot(amino,ca,cb,cp,dphi,n,o)
! Performs a right-handed rotation by angle dphi about
! the n-Ca axis of the phi angle of amino acid amino.
! This should rotate all atoms with atomnumber greater
! than atomnumber(Ca(amino)). 
      use defs; implicit none
      integer(i4b), intent(in) :: amino
      real(dp), dimension(:,:), intent(inout) :: coord
      real(dp), intent(in) :: dphi
      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,CA,N
      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.
! In general, 
!  Rot = A*sin(theta)/|axis| 
!      + (1 - cos(theta))*C/|axis|^2 + Id*cos(theta)

!     axis = ca(:,amino) - n(:,amino)
      CA = coord(:,CAatomnumber(amino))
      N = coord(:,CAatomnumber(amino)-1)
      axis = CA - N
      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(dphi) + (Id - C)*cos(dphi) + C
! Compensate for the displacement of the origin:     
! A' - N = Rot( A - N )  so A' = Rot(A) + N - Rot(N)
!     trans = n(:,amino) - matmul(Rot,n(:,amino))
      trans = N - matmul(Rot,N)
 do i = CAatomnumber(amino)+1, numberofatoms
   coord(:,i) = matmul(Rot,coord(:,i)) + trans
 end do

end subroutine phi_rot