Merge branch 'devel' into AFM
[unres.git] / source / unres / src_MD-NEWSC / unres.F
diff --git a/source/unres/src_MD-NEWSC/unres.F b/source/unres/src_MD-NEWSC/unres.F
new file mode 100644 (file)
index 0000000..e4f54cb
--- /dev/null
@@ -0,0 +1,799 @@
+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'
+      include 'COMMON.REMD'
+      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  03/28/12 23:29 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
+      if (modecalc.eq.-2) then
+        call test
+        stop
+      else if (modecalc.eq.-1) then
+        write(iout,*) "call check_sc_map next"
+        call check_bond
+        stop
+      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
+      else if (modecalc.eq.1) then
+        call exec_regularize
+      else if (modecalc.eq.2) then
+        call exec_thread
+      else if (modecalc.eq.3 .or. modecalc .eq.6) then
+        call exec_MC
+      else if (modecalc.eq.4) then
+        call exec_mult_eeval_or_minim
+      else if (modecalc.eq.5) then
+        call exec_checkgrad
+      else if (ModeCalc.eq.7) then
+        call exec_map
+      else if (ModeCalc.eq.8) then
+        call exec_CSA
+      else if (modecalc.eq.11) then
+        call exec_softreg
+      else if (modecalc.eq.12) then
+        call exec_MD
+      else if (modecalc.eq.14) then
+        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_MD
+      include 'DIMENSIONS'
+#ifdef MPI
+      include "mpif.h"
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.CONTROL'
+      include 'COMMON.IOUNITS'
+      if (me.eq.king .or. .not. out1file)
+     &   write (iout,*) "Calling chainbuild"
+      call chainbuild
+      call MD
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine exec_MREMD
+      include 'DIMENSIONS'
+#ifdef MPI
+      include "mpif.h"
+      include 'COMMON.SETUP'
+      include 'COMMON.CONTROL'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.REMD'
+      if (me.eq.king .or. .not. out1file)
+     &   write (iout,*) "Calling chainbuild"
+      call chainbuild
+      if (me.eq.king .or. .not. out1file)
+     &   write (iout,*) "Calling REMD"
+      if (remd_mlist) then 
+        call MREMD
+      else
+        do i=1,nrep
+          remd_m(i)=1
+        enddo
+        call MREMD
+      endif
+#else
+      write (iout,*) "MREMD works on parallel machines only"
+#endif
+      return
+      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.REMD'
+      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
+#ifdef MPI
+      time00=MPI_Wtime()
+#else
+      time00=tcpu()
+#endif
+      call chainbuild_cart
+      if (split_ene) then
+       print *,"Processor",myrank," after chainbuild"
+       icall=1
+       call etotal_long(energy_long(0))
+       write (iout,*) "Printing long range energy"
+       call enerprint(energy_long(0))
+       call etotal_short(energy_short(0))
+       write (iout,*) "Printing short range energy"
+       call enerprint(energy_short(0))
+       do i=0,n_ene
+         energy(i)=energy_long(i)+energy_short(i)
+         write (iout,*) i,energy_long(i),energy_short(i),energy(i)
+       enddo
+       write (iout,*) "Printing long+short range energy"
+       call enerprint(energy(0))
+      endif
+      call etotal(energy(0))
+#ifdef MPI
+      time_ene=MPI_Wtime()-time00
+#else
+      time_ene=tcpu()-time00
+#endif
+      write (iout,*) "Time for energy evaluation",time_ene
+      print *,"after etotal"
+      etota = energy(0)
+      etot =etota
+      call enerprint(energy(0))
+      call hairpin(.true.,nharp,iharp)
+      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
+        call hairpin(.true.,nharp,iharp)
+        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
+c---------------------------------------------------------------------------
+      subroutine exec_regularize
+      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.REMD'
+      include 'COMMON.MD'
+      include 'COMMON.SBRIDGE'
+      double precision energy(0:n_ene)
+
+      call gen_dist_constr
+      call sc_conf
+      call intout
+      call regularize(nct-nnt+1,etot,rms,cref(1,nnt),iretcode)
+      call etotal(energy(0))
+      energy(0)=energy(0)-energy(14)
+      etot=energy(0)
+      call enerprint(energy(0))
+      call intout
+      call briefout(0,etot)
+      if (outpdb) call pdbout(etot,titel(:32),ipdb)
+      if (outmol2) call mol2out(etot,titel(:32))
+      if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+      write (iout,'(a,i3)') 'SUMSL return code:',iretcode
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine exec_thread
+      include 'DIMENSIONS'
+#ifdef MP
+      include "mpif.h"
+#endif
+      include "COMMON.SETUP"
+      call thread_seq
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine exec_MC
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      character*10 nodeinfo
+      double precision varia(maxvar)
+#ifdef MPI
+      include "mpif.h"
+#endif
+      include "COMMON.SETUP"
+      include 'COMMON.CONTROL'
+      call mcm_setup
+      if (minim) then
+#ifdef MPI
+        if (modecalc.eq.3) then
+          call do_mcm(ipar)
+        else
+          call entmcm
+        endif
+#else
+        if (modecalc.eq.3) then
+          call do_mcm(ipar)
+        else
+          call entmcm
+        endif
+#endif
+      else
+        call monte_carlo
+      endif
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine exec_mult_eeval_or_minim
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+      dimension muster(mpi_status_size)
+#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.REMD'
+      include 'COMMON.MD'
+      include 'COMMON.SBRIDGE'
+      double precision varia(maxvar)
+      dimension ind(6)
+      double precision energy(0:n_ene)
+      logical eof
+      eof=.false.
+#ifdef MPI
+      if(me.ne.king) then
+        call minim_mcmf
+        return
+      endif
+
+      close (intin)
+      open(intin,file=intinname,status='old')
+      write (istat,'(a5,30a12)')"#    ",
+     &  (wname(print_order(i)),i=1,nprint_ene)
+      if (refstr) then
+        write (istat,'(a5,30a12)')"#    ",
+     &   (ename(print_order(i)),i=1,nprint_ene),
+     &   "ETOT total","RMSD","nat.contact","nnt.contact","cont.order"
+      else
+        write (istat,'(a5,30a12)')"#    ",
+     &    (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)
+#ifdef MPI
+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)
+#endif
+            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.)
+            write (istat,'(i5,30(f12.3))') iconf,
+     &      (energy(print_order(i)),i=1,nprint_ene),etot,
+     &       rms,frac,frac_nn,co
+cjlee end
+          else
+            write (istat,'(i5,30(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
+c      do n=1,nzsc
+      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)
+#ifdef MPI
+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)
+#endif
+            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(it,etot)
+c          if (minim) call briefout(it,etot)
+          if (refstr) then 
+            call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+            write (istat,'(i5,30(f12.3))') iconf,
+     &     (energy(print_order(i)),i=1,nprint_ene),etot,
+     &     rms,frac,frac_nn,co
+          else
+            write (istat,'(i5,30(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)
+#ifdef MPI
+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)
+#endif
+            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(it,etot)
+        if (refstr) then 
+          call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+          write (istat,'(i5,30(f12.3))') iconf,
+     &   (energy(print_order(i)),i=1,nprint_ene),etot,
+     &   rms,frac,frac_nn,co
+        else
+          write (istat,'(i5,30(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
+#else
+      close (intin)
+      open(intin,file=intinname,status='old')
+      write (istat,'(a5,20a12)')"#    ",
+     &   (wname(print_order(i)),i=1,nprint_ene)
+      write (istat,'("#    ",20(1pe12.4))')
+     &   (weights(print_order(i)),i=1,nprint_ene)
+      if (refstr) then
+        write (istat,'(a5,20a12)')"#    ",
+     &   (ename(print_order(i)),i=1,nprint_ene),
+     &   "ETOT total","RMSD","nat.contact","nnt.contact"
+      else
+        write (istat,'(a5,14a12)')"#    ",
+     &   (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
+      endif
+      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)
+#ifdef MPI
+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)
+#endif
+            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
+        if (minim) call minimize(etot,varia,iretcode,nfun)
+        call etotal(energy(0))
+
+        etot=energy(0)
+        call enerprint(energy(0))
+        if (minim) call briefout(it,etot) 
+        if (refstr) then 
+          call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+          write (istat,'(i5,18(f12.3))') iconf,
+     &   (energy(print_order(i)),i=1,nprint_ene),
+     &   etot,rms,frac,frac_nn,co
+cjlee end
+        else
+          write (istat,'(i5,14(f12.3))') iconf,
+     &   (energy(print_order(i)),i=1,nprint_ene),etot
+        endif
+      enddo
+   11 continue
+ 1100 continue
+#endif
+      return
+      end
+c---------------------------------------------------------------------------
+      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'
+      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
+      call read_fragments
+      read(inp,*) t_bath
+c!      t_bath = 300
+      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_map
+C Energy maps
+      call map_read
+      call map
+      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
+csa      call together
+      write (iout,*) "CSA is not supported in this version"
+#else
+csa      write (iout,*) "CSA works on parallel machines only"
+      write (iout,*) "CSA is not supported in this version"
+#endif
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine exec_softreg
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CONTROL'
+      double precision energy(0:max_ene)
+      logical debug /.false./
+      call chainbuild
+      call etotal(energy(0))
+      call enerprint(energy(0))
+      if (.not.lsecondary) then
+        write(iout,*) 'Calling secondary structure recognition'
+        call secondary2(debug)
+      else
+        write(iout,*) 'Using secondary structure supplied in pdb'
+      endif
+
+      call softreg
+
+      call etotal(energy(0))
+      etot=energy(0)
+      call enerprint(energy(0))
+      call intout
+      call briefout(0,etot)
+      call secondary2(.true.)
+      if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+      return
+      end