subroutine pro1psi(residue,numberofatoms,dtheta,CAatomnumber,coord)
! When acid_{residue+1} = PRO, with no other nearby PROs,
! finds psi1,phi2,psi2,psi3 with RRRR ~ 1
use defs; use eyes; implicit none
integer(i4b), intent(in) :: residue,numberofatoms
real(dp), intent(in) :: dtheta
integer(i4b), dimension(:), intent(in) :: CAatomnumber
real(dp), dimension(:,:), intent(inout) :: coord  
! internal variables:
integer(i4b) :: amino
real(dp),dimension(3) :: CA1,C1,CA2,C2,N3,CA3,C3,axis1,axis2, &
                         axis3,axis4
real(dp),dimension(4) :: angles
real(dp),dimension(3,3) :: Axes
real(dp) :: norm,maxfactor,dphi,dpsi

! STUFF FOR DGESV
integer(i4b), parameter :: N = 3,NHRS=1,LDA=3,LDB=3
integer(i4b),dimension(3) :: IPIV
! real(dp),dimension(3,3) :: A
real(dp),dimension(3,1) :: B
integer(i4b) :: INFO 

CA1 = coord(:,CAatomnumber(residue))
C1 = coord(:,CAatomnumber(residue)+1) 
CA2 = coord(:,CAatomnumber(residue+1))
C2 = coord(:,CAatomnumber(residue+1)+1)
N3 = coord(:,CAatomnumber(residue+2)-1)
CA3 = coord(:,CAatomnumber(residue+2))
C3 = coord(:,CAatomnumber(residue+2)+1)

! First find the values of the four small angles
axis1 = C1 - CA1; axis1 = axis1/sqrt(dot_product(axis1,axis1))!psi1
axis2 = C2 - CA2; axis2 = axis2/sqrt(dot_product(axis2,axis2))!psi2
axis3 = CA3 - N3; axis3 = axis3/sqrt(dot_product(axis3,axis3))!phi3
axis4 = C3 - CA3; axis4 = axis4/sqrt(dot_product(axis4,axis4))!psi3  
Axes(:,1) = axis1; Axes(:,2) = axis2;  Axes(:,3) = axis3
! A = Axes
B(:,1) = axis4
call DGESV ( N, NHRS, Axes, LDA, IPIV, B, LDB, INFO )
!write(6,*) 'INFO =',INFO
!write(6,*) 'X =', B
!zero = B(1,1)*axis1 + B(2,1)*axis2 + B(3,1)*axis3 - axis4
!write(6,*) 'zero =', zero
! At finite angles the linear rule gives 

angles(1) = B(1,1); angles(2) = B(2,1)
angles(3) = B(3,1); angles(4) = -1.0_dp
norm = dot_product(angles,angles)
maxfactor = dtheta/sqrt(norm)
angles = maxfactor*angles
dpsi = angles(1); amino = residue      ! psi1
call psi_rot(amino,coord,dpsi,Caatomnumber,numberofatoms)
dpsi = angles(2); amino = residue + 1  ! psi2
call psi_rot(amino,coord,dpsi,Caatomnumber,numberofatoms)
dphi = angles(3); amino = residue + 2  ! phi3
call phi_rot(amino,coord,dphi,Caatomnumber,numberofatoms)
dpsi = angles(4); amino = residue + 2  ! psi3  
call psi_rot(amino,coord,dpsi,Caatomnumber,numberofatoms)

end subroutine pro1psi