CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                              C
C                                U N R E S                                     C
C                                                                              C
C Program to carry out conformational search of proteins in an united-residue  C
C approximation.                                                               C
C                                                                              C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'


#ifdef MPI
      include 'mpif.h'
      include 'COMMON.SETUP'
#endif
      include 'COMMON.TIME1'
      include 'COMMON.INTERACT'
      include 'COMMON.NAMES'
      include 'COMMON.GEO'
      include 'COMMON.HEADER'
      include 'COMMON.CONTROL'
      include 'COMMON.CONTACTS'
      include 'COMMON.CHAIN'
      include 'COMMON.VAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.FFIELD'
c      include 'COMMON.REMD'
c      include 'COMMON.MD'
      include 'COMMON.SBRIDGE'
      double precision hrtime,mintime,sectime
      character*64 text_mode_calc(-2:14) /'test',
     & 'SC rotamer distribution',
     & 'Energy evaluation or minimization',
     & 'Regularization of PDB structure',
     & 'Threading of a sequence on PDB structures',
     & 'Monte Carlo (with minimization) ',
     & 'Energy minimization of multiple conformations',
     & 'Checking energy gradient',
     & 'Entropic sampling Monte Carlo (with minimization)',
     & 'Energy map',
     & 'CSA calculations',
     & 'Not used 9',
     & 'Not used 10',
     & 'Soft regularization of PDB structure',
     & 'Mesoscopic molecular dynamics (MD) ',
     & 'Not used 13',
     & 'Replica exchange molecular dynamics (REMD)'/
      external ilen

c      call memmon_print_usage()

      call init_task
      if (me.eq.king)
     & write(iout,*)'### LAST MODIFIED  11/03/09 1:19PM by czarek'  
      if (me.eq.king) call cinfo
C Read force field parameters and job setup data
      call readrtns
      call flush(iout)
C
      if (me.eq.king .or. .not. out1file) then
       write (iout,'(2a/)') 
     & text_mode_calc(modecalc)(:ilen(text_mode_calc(modecalc))),
     & ' calculation.' 
       if (minim) write (iout,'(a)') 
     &  'Conformations will be energy-minimized.'
       write (iout,'(80(1h*)/)') 
      endif
      call flush(iout)
#ifdef MPI
      if (fg_rank.gt.0) then
C Fine-grain slaves just do energy and gradient components.
        call ergastulum ! slave workhouse in Latin
      else
#endif
      if (modecalc.eq.0) then
        call exec_eeval_or_minim
      else if (modecalc.eq.5) then
         call exec_checkgrad
      else
        write (iout,'(a)') 'This calculation type is not supported',
     &   ModeCalc
      endif
#ifdef MPI
      endif
C Finish task.
      if (fg_rank.eq.0) call finish_task
c      call memmon_print_usage()
#ifdef TIMING
       call print_detailed_timing
#endif
      call MPI_Finalize(ierr)
      stop 'Bye Bye...'
#else
      call dajczas(tcpu(),hrtime,mintime,sectime)
      stop '********** Program terminated normally.'
#endif
      end
c---------------------------------------------------------------------------
      subroutine exec_eeval_or_minim
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.SETUP'
      include 'COMMON.TIME1'
      include 'COMMON.INTERACT'
      include 'COMMON.NAMES'
      include 'COMMON.GEO'
      include 'COMMON.HEADER'
      include 'COMMON.CONTROL'
      include 'COMMON.CONTACTS'
      include 'COMMON.CHAIN'
      include 'COMMON.VAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.FFIELD'
      include 'COMMON.SBRIDGE'
      common /srutu/ icall
      double precision energy(0:n_ene),varia(maxvar)
      double precision energy_long(0:n_ene),energy_short(0:n_ene)
      if (indpdb.eq.0) call chainbuild
#ifdef MPI
      time00=MPI_Wtime()
#else
      time00=tcpu()
#endif
      call chainbuild_cart
      call etotal(energy(0))
#ifdef MPI
      time_ene=MPI_Wtime()-time00
#else
      time_ene=tcpu()
#endif
      write (iout,*) "Time for energy evaluation",time_ene
      print *,"after etotal"
      etota = energy(0)
      etot =etota
      call enerprint(energy(0))
c      call hairpin(.true.,nharp,iharp)
c      call secondary2(.true.)
      if (minim) then

crc overlap test
        if (overlapsc) then
          print *, 'Calling OVERLAP_SC'
          call overlap_sc(fail)
        endif

        if (searchsc) then
          call sc_move(2,nres-1,10,1d10,nft_sc,etot)
          print *,'SC_move',nft_sc,etot
          write(iout,*) 'SC_move',nft_sc,etot
        endif

        if (dccart) then
          print *, 'Calling MINIM_DC'
#ifdef MPI
          time1=MPI_WTIME()
#else
          time1=tcpu()
#endif
          call minim_dc(etot,iretcode,nfun)
        else
          if (indpdb.ne.0) then 
            call bond_regular
            call chainbuild
          endif
          call geom_to_var(nvar,varia)
          print *,'Calling MINIMIZE.'
#ifdef MPI
          time1=MPI_WTIME()
#else
          time1=tcpu()
#endif
          call minimize(etot,varia,iretcode,nfun)
        endif
        print *,'SUMSL return code is',iretcode,' eval ',nfun
#ifdef MPI
        evals=nfun/(MPI_WTIME()-time1)
#else
        evals=nfun/(tcpu()-time1)
#endif
        print *,'# eval/s',evals
        print *,'refstr=',refstr
c        call hairpin(.true.,nharp,iharp)
c        call secondary2(.true.)
        call etotal(energy(0))
        etot = energy(0)
        call enerprint(energy(0))

        call intout
        call briefout(0,etot)
c        if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
          write (iout,'(a,i3)') 'SUMSL return code:',iretcode
          write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
          write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
c      else
c        print *,'refstr=',refstr
c        if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
c        call briefout(0,etot)
      endif
      if (outpdb) call pdbout(etot,titel(:32),ipdb)
      if (outmol2) call mol2out(etot,titel(:32))
      return
      end

      subroutine exec_checkgrad
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.SETUP'
      include 'COMMON.TIME1'
      include 'COMMON.INTERACT'
      include 'COMMON.NAMES'
      include 'COMMON.GEO'
      include 'COMMON.HEADER'
      include 'COMMON.CONTROL'
      include 'COMMON.CONTACTS'
      include 'COMMON.CHAIN'
      include 'COMMON.VAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.FFIELD'
c      include 'COMMON.REMD'
      include 'COMMON.MD_'
      include 'COMMON.SBRIDGE'
      common /srutu/ icall
      double precision energy(0:max_ene)
c      do i=2,nres
c        vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0)
c        if (itype(i).ne.10) 
c     &      vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0)
c      enddo
      if (indpdb.eq.0) call chainbuild
c      do i=0,nres
c        do j=1,3
c          dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0)
c        enddo
c      enddo
c      do i=1,nres-1
c        if (itype(i).ne.10) then
c          do j=1,3
c            dc(j,i+nres)=dc(j,i+nres)+ran_number(-0.2d0,0.2d0)
c          enddo
c        endif
c      enddo
c      do j=1,3
c        dc(j,0)=ran_number(-0.2d0,0.2d0)
c      enddo
      usampl=.true.
      totT=1.d0
      eq_time=0.0d0
c      call read_fragments
      call chainbuild_cart
      call cartprint
      call intout
      icall=1
      call etotal(energy(0))
      etot = energy(0)
      call enerprint(energy(0))
      write (iout,*) "Uconst",Uconst," Uconst_back",uconst_back
      print *,'icheckgrad=',icheckgrad
      goto (10,20,30) icheckgrad
  10  call check_ecartint
      return
  20  call check_cartgrad
      return
  30  call check_eint
      return
      end