subroutine denature(nres,numberofatoms,kindofatom,           &
                    resnum,coord,CAatomnumber,               &
                    pdbCA,ndenaturing,ndone,                 &
                    intervalnum,den_angle_scale,rmsdistance, &
                    ndensweep,denlogfile,mydenoutfile)
use defs; use eyes; use Ran_Lux_Mod; implicit none
integer(i4b), intent(in) :: nres,numberofatoms
character(len=4), dimension(:,:), intent(in) :: kindofatom
integer(i4b),dimension(:), intent(in) :: resnum
real(dp), dimension(:,:), intent(inout) :: coord 
integer(i4b), dimension(:), intent(in) :: CAatomnumber
real(dp), dimension(:,:), intent(in) :: pdbCA
integer(i4b), intent(in) :: ndenaturing,ndone,intervalnum
! We get intervalnum alternating intervals starting
! with an inverse MC followed by a random interval.
! So intervalnum = 1 gives only one streching interval,
! while intervalnum = 5 gives IMC,R,IMC,R,IMC. 
real(dp), intent(in) :: den_angle_scale
real(dp), intent(out) :: rmsdistance 
integer(i4b), intent(out) :: ndensweep 
character(len=*), intent(in) :: denlogfile,mydenoutfile
integer(i4b) :: sweep,length,residue,imar,     &
                ndecade,sweepdecade
real(sp), dimension(:), allocatable :: RVec
real(dp), dimension(3,numberofatoms) :: coord0
real(dp), dimension(3,nres) :: CA,CA0
real(dp) :: dphi,dpsi,bestdistance
length = 2*(nres-1); allocate(RVec(length))
coord0 = coord; call CAtrace(nres,coord,CAatomnumber,CA)
CA0 = CA
call distance(nres,CA,pdbCA,rmsdistance)
bestdistance = rmsdistance
ndecade = ndenaturing/intervalnum
do sweep = 1, ndenaturing
 sweepdecade = (sweep-1)/ndecade
 if ( mod(sweepdecade,2) == 1 ) go to 3 ! go to random loop
! otherwise start reverse Monte Carlo loop
 call RanLux(RVec); imar = 0! get random numbers
 imar = imar + 1; dpsi = den_angle_scale*(RVec(imar) - 0.5_dp)
 residue = 1
 call psi_rot(residue,coord,dpsi,CAatomnumber,numberofatoms)
 call CAtrace(nres,coord,CAatomnumber,CA)
 call distance(nres,CA,pdbCA,rmsdistance) 
 if ( rmsdistance <= bestdistance ) then
  coord = coord0; CA = CA0   ! go back to previous state
 else
  bestdistance = rmsdistance ! keep new state
  coord0 = coord; CA0 = CA
 end if
 do residue = 2, nres - 1
  if ( kindofatom(2,CAatomnumber(residue)) == ' PRO' ) go to 1
  imar = imar + 1; dphi = den_angle_scale*(RVec(imar) - 0.5_dp)
  call phi_rot(residue,coord,dphi,CAatomnumber,numberofatoms)
  call CAtrace(nres,coord,CAatomnumber,CA)
  call distance(nres,CA,pdbCA,rmsdistance)
  if ( rmsdistance <= bestdistance ) then
   coord = coord0; CA = CA0   ! go back to previous state
  else
   bestdistance = rmsdistance ! keep new state
   coord0 = coord; CA0 = CA
  end if
1 imar = imar + 1; dpsi = den_angle_scale*(RVec(imar) - 0.5_dp)
  call psi_rot(residue,coord,dpsi,CAatomnumber,numberofatoms)
  call CAtrace(nres,coord,CAatomnumber,CA)
  call distance(nres,CA,pdbCA,rmsdistance)
  if ( rmsdistance <= bestdistance ) then
   coord = coord0; CA = CA0   ! go back to previous state
  else
   bestdistance = rmsdistance ! keep new state
   coord0 = coord; CA0 = CA
  end if
 end do  ! end of protein
 if ( kindofatom(2,CAatomnumber(residue)) == ' PRO' ) go to 2
 imar = imar + 1; dphi =den_angle_scale*(RVec(imar) - 0.5_dp)
 call phi_rot(nres,coord,dphi,CAatomnumber,numberofatoms)
 call CAtrace(nres,coord,CAatomnumber,CA)
 call distance(nres,CA,pdbCA,rmsdistance)
 ndensweep = sweep
 if ( rmsdistance <= bestdistance ) then
  coord = coord0; CA = CA0   ! go back to previous state
 else
  bestdistance = rmsdistance ! keep new state
  coord0 = coord; CA0 = CA
 end if 
2 continue
! end of reverse Monte Carlo loop
go to 9 ! skipping really random loop
! start of really random loop
3 call RanLux(RVec); imar = 0! get random numbers
 imar = imar + 1; dpsi = den_angle_scale*(RVec(imar) - 0.5_dp)
 residue = 1
 call psi_rot(residue,coord,dpsi,CAatomnumber,numberofatoms)
! call CAtrace(nres,coord,CAatomnumber,CA)
! call distance(nres,CA,pdbCA,rmsdistance) 
 do residue = 2, nres - 1
  if ( kindofatom(2,CAatomnumber(residue)) == ' PRO' ) go to 5
  imar = imar + 1; dphi = den_angle_scale*(RVec(imar) - 0.5_dp)
  call phi_rot(residue,coord,dphi,CAatomnumber,numberofatoms)
!  call CAtrace(nres,coord,CAatomnumber,CA)
!  call distance(nres,CA,pdbCA,rmsdistance)
5 imar = imar + 1; dpsi = den_angle_scale*(RVec(imar) - 0.5_dp)
  call psi_rot(residue,coord,dpsi,CAatomnumber,numberofatoms)
  call CAtrace(nres,coord,CAatomnumber,CA)
  call distance(nres,CA,pdbCA,rmsdistance)
 end do  ! end of protein
 if ( kindofatom(2,CAatomnumber(residue)) == ' PRO' ) go to 6
 imar = imar + 1; dphi =den_angle_scale*(RVec(imar) - 0.5_dp)
 call phi_rot(nres,coord,dphi,CAatomnumber,numberofatoms)
 call CAtrace(nres,coord,CAatomnumber,CA)
 call distance(nres,CA,pdbCA,rmsdistance)
 bestdistance = rmsdistance
 ndensweep = sweep
6 continue
9 continue
 if ( mod(sweep,50) == 0 ) then
 open(8,file=denlogfile,position='append')
 write(8,'(i7,1x,f10.6)') sweep+ndone, bestdistance
 close(8)
 call writepdb(mydenoutfile,numberofatoms,kindofatom, &
                    resnum,coord)
 end if
end do  ! end of denaturing 
call CAtrace(nres,coord,CAatomnumber,CA)
call distance(nres,CA,pdbCA,rmsdistance)
bestdistance = rmsdistance
deallocate(RVec)
call writepdb(mydenoutfile,numberofatoms,kindofatom, &
                   resnum,coord)
end subroutine denature