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)
C
c      if (modecalc.eq.-2) then
c        call test
c        stop
c      else if (modecalc.eq.-1) then
c        write(iout,*) "call check_sc_map next"
c        call check_bond
c        stop
c      endif
#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
c      else if (modecalc.eq.1) then
c        call exec_regularize
c      else if (modecalc.eq.2) then
c        call exec_thread
c      else if (modecalc.eq.3 .or. modecalc .eq.6) then
c        call exec_MC
      else if (modecalc.eq.4) then
        call exec_mult_eeval_or_minim
      else if (modecalc.eq.5) then
         call exec_checkgrad
c      else if (ModeCalc.eq.7) then
c        call exec_map
      else if (ModeCalc.eq.8) then
        call exec_CSA
c      else if (modecalc.eq.11) then
c        call exec_softreg
c      else if (modecalc.eq.12) then
c        call exec_MD
c      else if (modecalc.eq.14) then
c        call exec_MREMD
      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'
c      include 'COMMON.REMD'
c      include 'COMMON.MD'
      include 'COMMON.SBRIDGE'
      common /srutu/ icall
      double precision energy(0:n_ene)
      double precision energy_long(0:n_ene),energy_short(0:n_ene)
      double precision varia(maxvar)
      if (indpdb.eq.0) call chainbuild
c      time00=MPI_Wtime()
      call chainbuild_cart
c      if (split_ene) then
c       print *,"Processor",myrank," after chainbuild"
c       icall=1
c       call etotal_long(energy_long(0))
c       write (iout,*) "Printing long range energy"
c       call enerprint(energy_long(0))
c       call etotal_short(energy_short(0))
c       write (iout,*) "Printing short range energy"
c       call enerprint(energy_short(0))
c       do i=0,n_ene
c         energy(i)=energy_long(i)+energy_short(i)
c         write (iout,*) i,energy_long(i),energy_short(i),energy(i)
c       enddo
c       write (iout,*) "Printing long+short range energy"
c       call enerprint(energy(0))
c      endif
      call etotal(energy(0))
c      time_ene=MPI_Wtime()-time00
      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

        if (dccart) then
          print *, 'Calling MINIM_DC'
c          time1=MPI_WTIME()
          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.'
c          time1=MPI_WTIME()
          call minimize(etot,varia,iretcode,nfun)
        endif
        print *,'SUMSL return code is',iretcode,' eval ',nfun
c        evals=nfun/(MPI_WTIME()-time1)
        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)
        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
      else
        print *,'refstr=',refstr
        if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
        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
c      read(inp,*) t_bath
c      call rescale_weights(t_bath)
      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
c---------------------------------------------------------------------------
      subroutine exec_CSA
#ifdef MPI
      include "mpif.h"
#endif
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
C Conformational Space Annealling programmed by Jooyoung Lee.
C This method works only with parallel machines!
#ifdef MPI
      call together
#else
      write (iout,*) "CSA works on parallel machines only"
#endif
      return
      end
c---------------------------------------------------------------------------
#ifdef MPI
      subroutine exec_mult_eeval_or_minim
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'mpif.h'
      integer muster(mpi_status_size)
      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'
      double precision varia(maxvar)
      integer ind(6)
      double precision energy(0:n_ene)
      logical eof
      eof=.false.

      if(me.ne.king) then
        call minim_mcmf
        return
      endif

      close (intin)
      open(intin,file=intinname,status='old')
      write (istat,'(a5,100a12)')"#    ",
     &  (wname(print_order(i)),i=1,nprint_ene)
      if (refstr) then
        write (istat,'(a5,100a12)')"#    ",
     &   (ename(print_order(i)),i=1,nprint_ene),
     &   "ETOT total","RMSD","nat.contact","nnt.contact",
     &   "cont.order","TMscore"
      else
        write (istat,'(a5,100a12)')"#    ",
     &    (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
      endif

      if (.not.minim) then
        do while (.not. eof)
          if (read_cart) then
            read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene
            call read_x(intin,*11)
c Broadcast the order to compute internal coordinates to the slaves.
            if (nfgtasks.gt.1)
     &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
            call int_from_cart1(.false.)
          else
            read (intin,'(i5)',end=1100,err=1100) iconf
            call read_angles(intin,*11)
            call geom_to_var(nvar,varia)
            call chainbuild
          endif
          write (iout,'(a,i7)') 'Conformation #',iconf
          call etotal(energy(0))
          call briefout(iconf,energy(0))
          call enerprint(energy(0))
          etot=energy(0)
          if (refstr) then 
            call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
            call calc_tmscore(tm,.true.)
            write (istat,'(i5,100(f12.3))') iconf,
     &      (energy(print_order(i)),i=1,nprint_ene),etot,
     &       rms,frac,frac_nn,co,tm
          else
            write (istat,'(i5,100(f12.3))') iconf,
     &     (energy(print_order(i)),i=1,nprint_ene),etot
          endif
        enddo
1100    continue
        goto 1101
      endif

      mm=0
      imm=0
      nft=0
      ene0=0.0d0
      n=0
      iconf=0
      do while (.not. eof)
        mm=mm+1
        if (mm.lt.nodes) then
          if (read_cart) then
            read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
            call read_x(intin,*11)
c Broadcast the order to compute internal coordinates to the slaves.
            if (nfgtasks.gt.1) 
     &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
            call int_from_cart1(.false.)
          else
            read (intin,'(i5)',end=11,err=11) iconf
            call read_angles(intin,*11)
            call geom_to_var(nvar,varia)
            call chainbuild
          endif

          n=n+1
          write (iout,*) 'Conformation #',iconf,' read'
         imm=imm+1
         ind(1)=1
         ind(2)=n
         ind(3)=0
         ind(4)=0
         ind(5)=0
         ind(6)=0
         ene0=0.0d0
         call mpi_send(ind,6,mpi_integer,mm,idint,CG_COMM,
     *                  ierr)
         call mpi_send(varia,nvar,mpi_double_precision,mm,
     *                  idreal,CG_COMM,ierr)
         call mpi_send(ene0,1,mpi_double_precision,mm,
     *                  idreal,CG_COMM,ierr)
c         print *,'task ',n,' sent to worker ',mm,nvar
        else
         call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
     *                 CG_COMM,muster,ierr)
         man=muster(mpi_source)
c         print *,'receiving result from worker ',man,' (',iii1,iii,')'
         call mpi_recv(varia,nvar,mpi_double_precision, 
     *               man,idreal,CG_COMM,muster,ierr)
         call mpi_recv(ene,1,
     *               mpi_double_precision,man,idreal,
     *               CG_COMM,muster,ierr)
         call mpi_recv(ene0,1,
     *               mpi_double_precision,man,idreal,
     *               CG_COMM,muster,ierr)
c         print *,'result received from worker ',man,' sending now'

          call var_to_geom(nvar,varia)
          call chainbuild
          call etotal(energy(0))
          iconf=ind(2)
          write (iout,*)
          write (iout,*)
          write (iout,*) 'Conformation #',iconf," sumsl return code ",
     &                      ind(5)

          etot=energy(0)
          call enerprint(energy(0))
          call briefout(iconf,etot)
          if (refstr) then 
            call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
            call calc_tmscore(tm,.true.)
            write (istat,'(i5,100(f12.3))') iconf,
     &     (energy(print_order(i)),i=1,nprint_ene),etot,
     &     rms,frac,frac_nn,co,tm
          else
            write (istat,'(i5,100(f12.3))') iconf,
     &     (energy(print_order(i)),i=1,nprint_ene),etot
          endif

          imm=imm-1
          if (read_cart) then
            read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
            call read_x(intin,*11)
c Broadcast the order to compute internal coordinates to the slaves.
            if (nfgtasks.gt.1)
     &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
            call int_from_cart1(.false.)
          else
            read (intin,'(i5)',end=11,err=11) iconf
            call read_angles(intin,*11)
            call geom_to_var(nvar,varia)
            call chainbuild
          endif
          n=n+1
          write (iout,*) 'Conformation #',iconf,' read'
          imm=imm+1
          ind(1)=1
          ind(2)=n
          ind(3)=0
          ind(4)=0
          ind(5)=0
          ind(6)=0
          call mpi_send(ind,6,mpi_integer,man,idint,CG_COMM,
     *                  ierr)
          call mpi_send(varia,nvar,mpi_double_precision,man, 
     *                  idreal,CG_COMM,ierr)
          call mpi_send(ene0,1,mpi_double_precision,man,
     *                  idreal,CG_COMM,ierr)
          nf_mcmf=nf_mcmf+ind(4)
          nmin=nmin+1
        endif
      enddo
11    continue
      do j=1,imm
        call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
     *               CG_COMM,muster,ierr)
        man=muster(mpi_source)
        call mpi_recv(varia,nvar,mpi_double_precision, 
     *               man,idreal,CG_COMM,muster,ierr)
        call mpi_recv(ene,1,
     *               mpi_double_precision,man,idreal,
     *               CG_COMM,muster,ierr)
        call mpi_recv(ene0,1,
     *               mpi_double_precision,man,idreal,
     *               CG_COMM,muster,ierr)

        call var_to_geom(nvar,varia)
        call chainbuild
        call etotal(energy(0))
        iconf=ind(2)
        write (iout,*)
        write (iout,*)
        write (iout,*) 'Conformation #',iconf," sumsl return code ",
     &                  ind(5)

        etot=energy(0)
        call enerprint(energy(0))
        call briefout(iconf,etot)
        if (refstr) then 
          call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
          call calc_tmscore(tm,.true.)
          write (istat,'(i5,100(f12.3))') iconf,
     &   (energy(print_order(i)),i=1,nprint_ene),etot,
     &   rms,frac,frac_nn,co,tm
        else
          write (istat,'(i5,100(f12.3))') iconf,
     &    (energy(print_order(i)),i=1,nprint_ene),etot
        endif
        nmin=nmin+1
      enddo
1101  continue
      do i=1, nodes-1
         ind(1)=0
         ind(2)=0
         ind(3)=0
         ind(4)=0
         ind(5)=0
         ind(6)=0
         call mpi_send(ind,6,mpi_integer,i,idint,CG_COMM,
     *                  ierr)
      enddo
      return
      end
#else
      subroutine exec_mult_eeval_or_minim
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      write (iout,*) "Unsupported option in serial version"
      return
      end
#endif
c---------------------------------------------------------------------------