Merge branch 'devel' into feature-ga
[unres.git] / source / unres / src_MD / src / MREMD.F
diff --git a/source/unres/src_MD/src/MREMD.F b/source/unres/src_MD/src/MREMD.F
deleted file mode 100644 (file)
index 2d184b6..0000000
+++ /dev/null
@@ -1,2088 +0,0 @@
-#ifdef MPI
-      subroutine MREMD
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'mpif.h'
-      include 'COMMON.CONTROL'
-      include 'COMMON.VAR'
-      include 'COMMON.MD'
-#ifndef LANG0
-      include 'COMMON.LANGEVIN'
-#else
-      include 'COMMON.LANGEVIN.lang0'
-#endif
-      include 'COMMON.CHAIN'
-      include 'COMMON.DERIV'
-      include 'COMMON.GEO'
-      include 'COMMON.LOCAL'
-      include 'COMMON.INTERACT'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.NAMES'
-      include 'COMMON.TIME1'
-      include 'COMMON.REMD'
-      include 'COMMON.SETUP'
-      include 'COMMON.MUCA'
-      integer ERRCODE
-      double precision cm(3),L(3),vcm(3)
-      double precision energia(0:n_ene)
-      double precision remd_t_bath(maxprocs)
-      integer iremd_iset(maxprocs)
-      integer*2 i_index
-     &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
-      double precision remd_ene(0:n_ene+4,maxprocs),t_bath_old
-      integer iremd_acc(maxprocs),iremd_tot(maxprocs)
-      integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
-      integer ilen,rstcount
-      external ilen
-      character*50 tytul
-      common /gucio/ cm
-      integer itime
-cold      integer nup(0:maxprocs),ndown(0:maxprocs)
-      integer rep2i(0:maxprocs),ireqi(maxprocs)
-      integer icache_all(maxprocs)
-      integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
-      logical synflag,end_of_run,file_exist /.false./,ovrtim
-      real ene_tol /1.0e-5/
-
-cdeb      imin_itime_old=0
-      ntwx_cache=0
-      time00=MPI_WTIME()
-      time01=time00
-      if(me.eq.king.or..not.out1file) then
-       write  (iout,*) 'MREMD',nodes,'time before',time00-walltime
-       write (iout,*) "NREP=",nrep
-      endif
-
-      synflag=.false.
-      if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
-        call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
-      endif
-      mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
-
-cd      print *,'MREMD',nodes
-cd      print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
-cde      write (iout,*) "Start MREMD: me",me," t_bath",t_bath
-
-      if(hremd.gt.0) then
-         nset=hremd
-         do i=1,nset
-          mset(i)=1
-         enddo
-      endif
-
-      k=0
-      rep2i(k)=-1
-      do il=1,max0(nset,1)
-       do il1=1,max0(mset(il),1)
-        do i=1,nrep
-         iremd_acc(i)=0
-         iremd_acc_usa(i)=0
-         iremd_tot(i)=0
-         do j=1,remd_m(i)
-          i2rep(k)=i
-          i2set(k)=il
-          rep2i(i)=k
-          k=k+1
-          i_index(i,j,il,il1)=k
-         enddo
-        enddo
-       enddo
-      enddo
-
-      if(me.eq.king.or..not.out1file) then
-       write(iout,*) "i2rep",(i2rep(i),i=0,nodes-1)
-       write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
-       write(iout,*) "i,j,il,il1,i_index(i,j,il,il1)"
-       do il=1,nset
-        do il1=1,mset(il)
-         do i=1,nrep
-          do j=1,remd_m(i)
-           write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
-          enddo
-         enddo
-        enddo
-       enddo
-      endif
-
-c      print *,'i2rep',me,i2rep(me)
-c      print *,'rep2i',(rep2i(i),i=0,nrep)
-
-cold       if (i2rep(me).eq.nrep) then
-cold        nup(0)=0
-cold       else
-cold        nup(0)=remd_m(i2rep(me)+1)
-cold        k=rep2i(int(i2rep(me)))+1
-cold        do i=1,nup(0)
-cold         nup(i)=k
-cold         k=k+1
-cold        enddo
-cold       endif
-
-cd       print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
-
-cold       if (i2rep(me).eq.1) then
-cold        ndown(0)=0
-cold       else
-cold        ndown(0)=remd_m(i2rep(me)-1)
-cold        k=rep2i(i2rep(me)-2)+1
-cold        do i=1,ndown(0)
-cold         ndown(i)=k
-cold         k=k+1
-cold        enddo
-cold       endif
-
-cd       print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
-
-       
-       write (*,*) "Processor",me," rest",rest,"
-     &   restart1fie",restart1file
-       if(rest.and.restart1file) then 
-           if (me.eq.king)
-     &     inquire(file=mremd_rst_name,exist=file_exist)
-cd           write (*,*) me," Before broadcast: file_exist",file_exist
-           call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
-     &          IERR)
-cd           write (*,*) me," After broadcast: file_exist",file_exist
-           if(file_exist) then 
-             if(me.eq.king.or..not.out1file)
-     &            write  (iout,*) 'Master is reading restart1file'
-             call read1restart(i_index)
-           else
-             if(me.eq.king.or..not.out1file)
-     &            write  (iout,*) 'WARNING : no restart1file'
-           endif
-
-           if(me.eq.king.or..not.out1file) then
-              write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
-              write(iout,*) "i_index"
-              do il=1,nset
-               do il1=1,mset(il)
-                do i=1,nrep
-                 do j=1,remd_m(i)
-                  write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
-                 enddo
-                enddo
-               enddo
-              enddo
-           endif 
-       endif
-
-       if(me.eq.king) then
-        if (rest.and..not.restart1file) 
-     &          inquire(file=mremd_rst_name,exist=file_exist)
-        if(.not.file_exist.and.rest.and..not.restart1file) 
-     &       write(iout,*) 'WARNING : no restart file',mremd_rst_name
-        IF (rest.and.file_exist.and..not.restart1file) THEN
-             write  (iout,*) 'Master is reading restart file',
-     &                        mremd_rst_name
-             open(irest2,file=mremd_rst_name,status='unknown')
-             read (irest2,*) 
-             read (irest2,*) (i2rep(i),i=0,nodes-1)
-             read (irest2,*) 
-             read (irest2,*) (ifirst(i),i=1,remd_m(1))
-             do il=1,nodes
-              read (irest2,*) 
-              read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
-              read (irest2,*) 
-              read (irest2,*) ndowna(0,il),
-     &                    (ndowna(i,il),i=1,ndowna(0,il))
-             enddo
-             if(usampl.or.hremd.gt.0) then
-              read (irest2,*)
-              read (irest2,*) nset
-              read (irest2,*) 
-              read (irest2,*) (mset(i),i=1,nset)
-              read (irest2,*) 
-              read (irest2,*) (i2set(i),i=0,nodes-1)
-              read (irest2,*) 
-              do il=1,nset
-               do il1=1,mset(il)
-                do i=1,nrep
-                  read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
-                enddo
-               enddo
-              enddo
-
-              write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
-              write(iout,*) "i_index"
-              do il=1,nset
-               do il1=1,mset(il)
-                do i=1,nrep
-                 do j=1,remd_m(i)
-                  write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
-                 enddo
-                enddo
-               enddo
-              enddo
-             endif
-
-             close(irest2)
-
-             write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
-             write (iout,'(a6,1000i5)') "ifirst",
-     &                    (ifirst(i),i=1,remd_m(1))
-             do il=1,nodes
-              write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
-     &                    (nupa(i,il),i=1,nupa(0,il))
-              write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
-     &                    (ndowna(i,il),i=1,ndowna(0,il))
-             enddo
-        ELSE IF (.not.(rest.and.file_exist)) THEN
-         do il=1,remd_m(1)
-          ifirst(il)=il
-         enddo
-
-         do il=1,nodes
-          if (i2rep(il-1).eq.nrep) then
-           nupa(0,il)=0
-          else
-           nupa(0,il)=remd_m(i2rep(il-1)+1)
-           k=rep2i(int(i2rep(il-1)))+1
-           do i=1,nupa(0,il)
-            nupa(i,il)=k+1
-            k=k+1
-           enddo
-          endif
-          if (i2rep(il-1).eq.1) then
-           ndowna(0,il)=0
-          else
-           ndowna(0,il)=remd_m(i2rep(il-1)-1)
-           k=rep2i(i2rep(il-1)-2)+1
-           do i=1,ndowna(0,il)
-            ndowna(i,il)=k+1
-            k=k+1
-           enddo
-          endif
-         enddo
-        
-        write (iout,'(a6,100i4)') "ifirst",
-     &                    (ifirst(i),i=1,remd_m(1))
-        do il=1,nodes
-         write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
-     &                    (nupa(i,il),i=1,nupa(0,il))
-         write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
-     &                    (ndowna(i,il),i=1,ndowna(0,il))
-        enddo
-        
-        ENDIF
-       endif
-c
-c      t_bath=retmin+(retmax-retmin)*me/(nodes-1)
-       if(.not.(rest.and.file_exist.and.restart1file)) then
-         if (me .eq. king) then 
-           t_bath=retmin
-         else 
-            t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
-         endif
-cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
-         if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
-
-       endif
-       if(usampl.or.hremd.gt.0) then
-          iset=i2set(me)
-          if (hremd.gt.0) call set_hweights(iset)
-          if(me.eq.king.or..not.out1file) 
-     &     write(iout,*) me,"iset=",iset,"t_bath=",t_bath
-       endif        
-c
-       stdfp=dsqrt(2*Rb*t_bath/d_time)
-       do i=1,ntyp
-          stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
-       enddo 
-
-c      print *,'irep',me,t_bath
-       if (.not.rest) then  
-        if (me.eq.king .or. .not. out1file)
-     &   write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
-        call rescale_weights(t_bath)
-       endif
-
-
-c------copy MD--------------
-c  The driver for molecular dynamics subroutines
-c------------------------------------------------
-      t_MDsetup=0.0d0
-      t_langsetup=0.0d0
-      t_MD=0.0d0
-      t_enegrad=0.0d0
-      t_sdsetup=0.0d0
-      if(me.eq.king.or..not.out1file)
-     & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
-#ifdef MPI
-      tt0 = MPI_Wtime()
-#else
-      tt0 = tcpu()
-#endif
-c Determine the inverse of the inertia matrix.
-      call setup_MD_matrices
-c Initialize MD
-      call init_MD
-      if (rest) then  
-       if (me.eq.king .or. .not. out1file)
-     &  write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
-       stdfp=dsqrt(2*Rb*t_bath/d_time)
-       do i=1,ntyp
-          stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
-       enddo 
-       if (lang.gt.0 .and. .not.surfarea) then
-         do i=nnt,nct-1
-          stdforcp(i)=stdfp*dsqrt(gamp)
-         enddo
-         do i=nnt,nct
-          stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
-         enddo
-       elseif (lang.gt.0 .and. surfarea ) then
-          call setup_fricmat
-       endif
-       call rescale_weights(t_bath)
-      endif
-
-#ifdef MPI
-      t_MDsetup = MPI_Wtime()-tt0
-#else
-      t_MDsetup = tcpu()-tt0
-#endif
-      rstcount=0 
-c   Entering the MD loop       
-#ifdef MPI
-      tt0 = MPI_Wtime()
-#else
-      tt0 = tcpu()
-#endif
-      if (lang.eq.2 .or. lang.eq.3) then
-#ifndef   LANG0
-        call setup_fricmat
-        if (lang.eq.2) then
-          call sd_verlet_p_setup
-        else
-          call sd_verlet_ciccotti_setup
-        endif
-        do i=1,dimen
-          do j=1,dimen
-            pfric0_mat(i,j,0)=pfric_mat(i,j)
-            afric0_mat(i,j,0)=afric_mat(i,j)
-            vfric0_mat(i,j,0)=vfric_mat(i,j)
-            prand0_mat(i,j,0)=prand_mat(i,j)
-            vrand0_mat1(i,j,0)=vrand_mat1(i,j)
-            vrand0_mat2(i,j,0)=vrand_mat2(i,j)
-          enddo
-        enddo
-        flag_stoch(0)=.true.
-        do i=1,maxflag_stoch
-          flag_stoch(i)=.false.
-        enddo
-#else
-        write (iout,*)
-     &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
-#ifdef MPI
-        call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
-#endif
-        stop
-#endif
-      else if (lang.eq.1 .or. lang.eq.4) then
-        call setup_fricmat
-      endif
-      time00=MPI_WTIME()
-      if (me.eq.king .or. .not. out1file)
-     & write(iout,*) 'Setup time',time00-walltime
-ctime      call flush(iout)
-#ifdef MPI
-      t_langsetup=MPI_Wtime()-tt0
-      tt0=MPI_Wtime()
-#else
-      t_langsetup=tcpu()-tt0
-      tt0=tcpu()
-#endif
-      itime=0
-      end_of_run=.false.
-      do while(.not.end_of_run)
-        itime=itime+1
-        if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
-        if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
-        rstcount=rstcount+1
-        if (lang.gt.0 .and. surfarea .and. 
-     &      mod(itime,reset_fricmat).eq.0) then
-          if (lang.eq.2 .or. lang.eq.3) then
-#ifndef   LANG0
-            call setup_fricmat
-            if (lang.eq.2) then
-              call sd_verlet_p_setup
-            else
-              call sd_verlet_ciccotti_setup
-            endif
-            do i=1,dimen
-              do j=1,dimen
-                pfric0_mat(i,j,0)=pfric_mat(i,j)
-                afric0_mat(i,j,0)=afric_mat(i,j)
-                vfric0_mat(i,j,0)=vfric_mat(i,j)
-                prand0_mat(i,j,0)=prand_mat(i,j)
-                vrand0_mat1(i,j,0)=vrand_mat1(i,j)
-                vrand0_mat2(i,j,0)=vrand_mat2(i,j)
-              enddo
-            enddo
-            flag_stoch(0)=.true.
-            do i=1,maxflag_stoch
-              flag_stoch(i)=.false.
-            enddo   
-#endif
-          else if (lang.eq.1 .or. lang.eq.4) then
-            call setup_fricmat
-          endif
-          write (iout,'(a,i10)') 
-     &      "Friction matrix reset based on surface area, itime",itime
-        endif
-        if (reset_vel .and. tbf .and. lang.eq.0 
-     &      .and. mod(itime,count_reset_vel).eq.0) then
-          call random_vel
-          if (me.eq.king .or. .not. out1file)
-     &     write(iout,'(a,f20.2)') 
-     &     "Velocities reset to random values, time",totT      
-          do i=0,2*nres
-            do j=1,3
-              d_t_old(j,i)=d_t(j,i)
-            enddo
-          enddo
-        endif
-               if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
-          call inertia_tensor  
-          call vcm_vel(vcm)
-          do j=1,3
-             d_t(j,0)=d_t(j,0)-vcm(j)
-          enddo
-          call kinetic(EK)
-          kinetic_T=2.0d0/(dimen3*Rb)*EK
-          scalfac=dsqrt(T_bath/kinetic_T)
-cd          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT    
-          do i=0,2*nres
-            do j=1,3
-              d_t_old(j,i)=scalfac*d_t(j,i)
-            enddo
-          enddo
-        endif  
-        if (lang.ne.4) then
-          if (RESPA) then
-c Time-reversible RESPA algorithm 
-c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
-            call RESPA_step(itime)
-          else
-c Variable time step algorithm.
-            call velverlet_step(itime)
-          endif
-        else
-#ifdef BROWN
-          call brown_step(itime)
-#else
-          print *,"Brown dynamics not here!"
-#ifdef MPI
-          call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
-#endif
-          stop
-#endif
-        endif
-        if(ntwe.ne.0) then
-          if (mod(itime,ntwe).eq.0) call statout(itime)
-        endif
-        if (mod(itime,ntwx).eq.0.and..not.traj1file) then
-          write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
-          if(mdpdb) then
-             call pdbout(potE,tytul,ipdb)
-          else 
-             call cartout(totT)
-          endif
-        endif
-        if (mod(itime,ntwx).eq.0.and.traj1file) then
-          if(ntwx_cache.lt.max_cache_traj_use) then
-            ntwx_cache=ntwx_cache+1
-          else
-           if (max_cache_traj_use.ne.1)
-     &      print *,itime,"processor ",me," over cache ",ntwx_cache
-           do i=1,ntwx_cache-1
-
-            totT_cache(i)=totT_cache(i+1)
-            EK_cache(i)=EK_cache(i+1)
-            potE_cache(i)=potE_cache(i+1)
-            t_bath_cache(i)=t_bath_cache(i+1)
-            Uconst_cache(i)=Uconst_cache(i+1)
-            iset_cache(i)=iset_cache(i+1)
-
-            do ii=1,nfrag
-             qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
-            enddo
-            do ii=1,npair
-             qpair_cache(ii,i)=qpair_cache(ii,i+1)
-            enddo
-            do ii=1,nfrag_back
-              utheta_cache(ii,i)=utheta_cache(ii,i+1)
-              ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
-              uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
-            enddo
-
-
-            do ii=1,nres*2
-             do j=1,3
-              c_cache(j,ii,i)=c_cache(j,ii,i+1)
-             enddo
-            enddo
-           enddo
-          endif
-
-            totT_cache(ntwx_cache)=totT
-            EK_cache(ntwx_cache)=EK
-            potE_cache(ntwx_cache)=potE
-            t_bath_cache(ntwx_cache)=t_bath
-            Uconst_cache(ntwx_cache)=Uconst
-            iset_cache(ntwx_cache)=iset
-
-            do i=1,nfrag
-             qfrag_cache(i,ntwx_cache)=qfrag(i)
-            enddo
-            do i=1,npair
-             qpair_cache(i,ntwx_cache)=qpair(i)
-            enddo
-            do i=1,nfrag_back
-              utheta_cache(i,ntwx_cache)=utheta(i)
-              ugamma_cache(i,ntwx_cache)=ugamma(i)
-              uscdiff_cache(i,ntwx_cache)=uscdiff(i)
-            enddo
-
-            do i=1,nres*2
-             do j=1,3
-              c_cache(j,i,ntwx_cache)=c(j,i)
-             enddo
-            enddo
-
-        endif
-        if ((rstcount.eq.1000.or.itime.eq.n_timestep)
-     &                         .and..not.restart1file) then
-
-           if(me.eq.king) then
-             open(irest1,file=mremd_rst_name,status='unknown')
-             write (irest1,*) "i2rep"
-             write (irest1,*) (i2rep(i),i=0,nodes-1)
-             write (irest1,*) "ifirst"
-             write (irest1,*) (ifirst(i),i=1,remd_m(1))
-             do il=1,nodes
-              write (irest1,*) "nupa",il
-              write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
-              write (irest1,*) "ndowna",il
-              write (irest1,*) ndowna(0,il),
-     &                   (ndowna(i,il),i=1,ndowna(0,il))
-             enddo
-             if(usampl.or.hremd.gt.0) then
-              write (irest1,*) "nset"
-              write (irest1,*) nset
-              write (irest1,*) "mset"
-              write (irest1,*) (mset(i),i=1,nset)
-              write (irest1,*) "i2set"
-              write (irest1,*) (i2set(i),i=0,nodes-1)
-              write (irest1,*) "i_index"
-              do il=1,nset
-               do il1=1,mset(il)
-                do i=1,nrep
-                  write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
-                enddo
-               enddo
-              enddo
-
-             endif
-             close(irest1)
-           endif
-           open(irest2,file=rest2name,status='unknown')
-           write(irest2,*) totT,EK,potE,totE,t_bath
-           do i=1,2*nres
-            write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
-           enddo
-           do i=1,2*nres
-            write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
-           enddo
-           if(usampl.or.hremd.gt.0) then
-             write (irest2,*) iset
-           endif
-          close(irest2)
-          rstcount=0
-        endif 
-
-c REMD - exchange
-c forced synchronization
-        if (mod(itime,i_sync_step).eq.0 .and. me.ne.king 
-     &                                .and. .not. mremdsync) then 
-            synflag=.false.
-            call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
-            if (synflag) then 
-               call mpi_recv(itime_master, 1, MPI_INTEGER,
-     &                             0,101,CG_COMM, status, ierr)
-               call mpi_barrier(CG_COMM, ierr)
-cdeb               if (out1file.or.traj1file) then
-cdeb                call mpi_gather(itime,1,mpi_integer,
-cdeb     &             icache_all,1,mpi_integer,king,
-cdeb     &             CG_COMM,ierr)                 
-               if(traj1file)
-     &          call mpi_gather(ntwx_cache,1,mpi_integer,
-     &             icache_all,1,mpi_integer,king,
-     &             CG_COMM,ierr)
-               if (.not.out1file)
-     &               write(iout,*) 'REMD synchro at',itime_master,itime
-               if (itime_master.ge.n_timestep .or. ovrtim()) 
-     &            end_of_run=.true.
-ctime               call flush(iout)
-            endif
-        endif
-
-c REMD - exchange
-        if ((mod(itime,nstex).eq.0.and.me.eq.king
-     &                  .or.end_of_run.and.me.eq.king )
-     &       .and. .not. mremdsync ) then
-           synflag=.true.
-           time01_=MPI_WTIME()
-           do i=1,nodes-1
-              call mpi_isend(itime,1,MPI_INTEGER,i,101,
-     &                                CG_COMM, ireqi(i), ierr)
-cd            write(iout,*) 'REMD synchro with',i
-cd            call flush(iout)
-           enddo
-           call mpi_waitall(nodes-1,ireqi,statusi,ierr)
-           call mpi_barrier(CG_COMM, ierr)
-           time01=MPI_WTIME()
-           write(iout,*) 'REMD synchro at',itime,'time=',time01-time01_
-           if (out1file.or.traj1file) then
-cdeb            call mpi_gather(itime,1,mpi_integer,
-cdeb     &             itime_all,1,mpi_integer,king,
-cdeb     &             CG_COMM,ierr)
-cdeb            write(iout,'(a19,8000i8)') ' REMD synchro itime',
-cdeb     &                    (itime_all(i),i=1,nodes)
-            if(traj1file) then
-cdeb             imin_itime=itime_all(1)
-cdeb             do i=2,nodes
-cdeb               if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
-cdeb             enddo
-cdeb             ii_write=(imin_itime-imin_itime_old)/ntwx
-cdeb             imin_itime_old=int(imin_itime/ntwx)*ntwx
-cdeb             write(iout,*) imin_itime,imin_itime_old,ii_write
-             call mpi_gather(ntwx_cache,1,mpi_integer,
-     &             icache_all,1,mpi_integer,king,
-     &             CG_COMM,ierr)
-c             write(iout,'(a19,8000i8)') '     ntwx_cache',
-c     &                    (icache_all(i),i=1,nodes)
-             ii_write=icache_all(1)
-             do i=2,nodes
-               if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
-             enddo
-c             write(iout,*) "MIN ii_write=",ii_write
-            endif
-           endif
-ctime           call flush(iout)
-        endif
-        if(mremdsync .and. mod(itime,nstex).eq.0) then
-           synflag=.true.
-           if (me.eq.king .or. .not. out1file)
-     &      write(iout,*) 'REMD synchro at',itime
-
-            if(traj1file) then
-             call mpi_gather(ntwx_cache,1,mpi_integer,
-     &             icache_all,1,mpi_integer,king,
-     &             CG_COMM,ierr)
-             if (me.eq.king) then
-               write(iout,'(a19,8000i8)') '     ntwx_cache',
-     &                    (icache_all(i),i=1,nodes)
-               ii_write=icache_all(1)
-               do i=2,nodes
-                 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
-               enddo
-               write(iout,*) "MIN ii_write=",ii_write
-             endif
-            endif
-ctest           call flush(iout)
-        endif
-        if (synflag) then
-c Update the time safety limiy
-          if (time001-time00.gt.safety) then
-            safety=time001-time00+600
-             if (me.eq.king .or. .not. out1file)
-     &       write (iout,*) "****** SAFETY increased to",safety," s"
-          endif
-          if (ovrtim()) end_of_run=.true.
-        endif
-        if(synflag.and..not.end_of_run) then
-           time02=MPI_WTIME()
-           synflag=.false.
-
-c           write(iout,*) 'REMD before',me,t_bath
-
-c           call mpi_gather(t_bath,1,mpi_double_precision,
-c     &             remd_t_bath,1,mpi_double_precision,king,
-c     &             CG_COMM,ierr)
-           potEcomp(n_ene+1)=t_bath
-           t_bath_old=t_bath
-           if (usampl) then
-             potEcomp(n_ene+2)=iset
-             if (iset.lt.nset) then
-               i_set_temp=iset
-               iset=iset+1
-               call EconstrQ
-               potEcomp(n_ene+3)=Uconst
-               iset=i_set_temp
-             endif
-             if (iset.gt.1) then
-               i_set_temp=iset
-               iset=iset-1
-               call EconstrQ
-               potEcomp(n_ene+4)=Uconst 
-               iset=i_set_temp
-             endif
-           endif
-           if(hremd.gt.0) potEcomp(n_ene+2)=iset   
-           call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
-     &             remd_ene(0,1),n_ene+5,mpi_double_precision,king,
-     &             CG_COMM,ierr)
-           if(lmuca) then 
-            call mpi_gather(elow,1,mpi_double_precision,
-     &             elowi,1,mpi_double_precision,king,
-     &             CG_COMM,ierr)
-            call mpi_gather(ehigh,1,mpi_double_precision,
-     &             ehighi,1,mpi_double_precision,king,
-     &             CG_COMM,ierr)
-           endif
-
-          time03=MPI_WTIME()
-          if (me.eq.king .or. .not. out1file) then
-            write(iout,*) 'REMD gather times=',time03-time01
-     &                                        ,time03-time02
-          endif
-
-          if (restart1file) call write1rst(i_index)
-
-          time04=MPI_WTIME()
-          if (me.eq.king .or. .not. out1file) then
-            write(iout,*) 'REMD writing rst time=',time04-time03
-          endif
-
-          if (traj1file) call write1traj
-cd debugging
-cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
-cdeb     &             icache_all,1,mpi_integer,king,
-cdeb     &             CG_COMM,ierr)
-cdeb            write(iout,'(a19,8000i8)') '  ntwx_cache after traj1file',
-cdeb     &                    (icache_all(i),i=1,nodes)
-cd end
-
-
-          time05=MPI_WTIME()
-          if (me.eq.king .or. .not. out1file) then
-            write(iout,*) 'REMD writing traj time=',time05-time04
-ctime            call flush(iout)
-          endif
-
-
-          if (me.eq.king) then
-            do i=1,nodes
-               remd_t_bath(i)=remd_ene(n_ene+1,i)
-               iremd_iset(i)=remd_ene(n_ene+2,i)
-            enddo
-            if(lmuca) then
-co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
-             do i=1,nodes
-               write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
-     &            elowi(i),ehighi(i)       
-             enddo
-            else
-              write(iout,*) 'REMD exchange temp,ene'
-              do i=1,nodes
-                write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
-                write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
-              enddo
-            endif
-c-------------------------------------           
-           IF(.not.usampl.and.hremd.eq.0) THEN
-            write (iout,*) "Enter exchnge, remd_m",remd_m(1),
-     &        " nodes",nodes
-ctime            call flush(iout)
-            write (iout,*) "remd_m(1)",remd_m(1)
-            do irr=1,remd_m(1)
-               i=ifirst(iran_num(1,remd_m(1)))
-             write (iout,*) "i",i
-ctime             call flush(iout)
-
-             do ii=1,nodes-1
-
-              write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
-             if(i.gt.0.and.nupa(0,i).gt.0) then
-              iex=i
-c              if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
-c                write (iout,*) 
-c     &  "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
-c                call flush(iout)
-c                call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
-c              endif
-c              do while (iex.eq.i)
-c                write (iout,*) "upper",nupa(int(nupa(0,i)),i)
-                iex=nupa(iran_num(1,int(nupa(0,i))),i)
-c              enddo
-c              write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
-              if (lmuca) then
-               call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
-              else
-c Swap temperatures between conformations i and iex with recalculating the free energies
-c following temperature changes.
-               ene_iex_iex=remd_ene(0,iex)
-               ene_i_i=remd_ene(0,i)
-c               write (iout,*) "i",i," ene_i_i",ene_i_i,
-c     &          " iex",iex," ene_iex_iex",ene_iex_iex
-c               write (iout,*) "rescaling weights with temperature",
-c     &          remd_t_bath(i)
-c               call flush(iout)
-               call rescale_weights(remd_t_bath(i))
-
-c               write (iout,*) "0,iex",remd_t_bath(i)
-c               call enerprint(remd_ene(0,iex))
-
-               call sum_energy(remd_ene(0,iex),.false.)
-               ene_iex_i=remd_ene(0,iex)
-c               write (iout,*) "ene_iex_i",remd_ene(0,iex)
-
-c               write (iout,*) "0,i",remd_t_bath(i)
-c               call enerprint(remd_ene(0,i))
-
-               call sum_energy(remd_ene(0,i),.false.)
-c               write (iout,*) "ene_i_i",remd_ene(0,i)
-c               call flush(iout)
-c               write (iout,*) "rescaling weights with temperature",
-c     &          remd_t_bath(iex)
-               if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
-                write (iout,*) "ERROR: inconsistent energies:",i,
-     &            ene_i_i,remd_ene(0,i)
-               endif
-               call rescale_weights(remd_t_bath(iex))
-
-c               write (iout,*) "0,i",remd_t_bath(iex)
-c               call enerprint(remd_ene(0,i))
-
-               call sum_energy(remd_ene(0,i),.false.)
-c               write (iout,*) "ene_i_iex",remd_ene(0,i)
-c               call flush(iout)
-               ene_i_iex=remd_ene(0,i)
-
-c               write (iout,*) "0,iex",remd_t_bath(iex)
-c               call enerprint(remd_ene(0,iex))
-
-               call sum_energy(remd_ene(0,iex),.false.)
-               if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
-                write (iout,*) "ERROR: inconsistent energies:",iex,
-     &            ene_iex_iex,remd_ene(0,iex)
-               endif
-c               write (iout,*) "ene_iex_iex",remd_ene(0,iex)
-c               write (iout,*) "i",i," iex",iex
-c               write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
-c     &           " ene_i_iex",ene_i_iex,
-c     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
-c               call flush(iout)
-               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
-     &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
-               delta=-delta
-c               write(iout,*) 'delta',delta
-c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
-c     &              (remd_ene(i)-remd_ene(iex))/Rb/
-c     &              (remd_t_bath(i)*remd_t_bath(iex))
-              endif
-              if (delta .gt. 50.0d0) then
-                delta=0.0d0
-              else
-#ifdef OSF 
-                if(isnan(delta))then
-                  delta=0.0d0
-                else if (delta.lt.-50.0d0) then
-                  delta=dexp(50.0d0)
-                else
-                  delta=dexp(-delta)
-                endif
-#else
-                delta=dexp(-delta)
-#endif
-              endif
-              iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
-              xxx=ran_number(0.0d0,1.0d0)
-c              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
-c              call flush(iout)
-              if (delta .gt. xxx) then
-                tmp=remd_t_bath(i)       
-                remd_t_bath(i)=remd_t_bath(iex)
-                remd_t_bath(iex)=tmp
-                remd_ene(0,i)=ene_i_iex
-                remd_ene(0,iex)=ene_iex_i
-                if(lmuca) then
-                  tmp=elowi(i)
-                  elowi(i)=elowi(iex)
-                  elowi(iex)=tmp  
-                  tmp=ehighi(i)
-                  ehighi(i)=ehighi(iex)
-                  ehighi(iex)=tmp  
-                endif
-
-
-                do k=0,nodes
-                  itmp=nupa(k,i)
-                  nupa(k,i)=nupa(k,iex)
-                  nupa(k,iex)=itmp
-                  itmp=ndowna(k,i)
-                  ndowna(k,i)=ndowna(k,iex)
-                  ndowna(k,iex)=itmp
-                enddo
-                do il=1,nodes
-                 if (ifirst(il).eq.i) ifirst(il)=iex
-                 do k=1,nupa(0,il)
-                  if (nupa(k,il).eq.i) then 
-                     nupa(k,il)=iex
-                  elseif (nupa(k,il).eq.iex) then 
-                     nupa(k,il)=i
-                  endif
-                 enddo
-                 do k=1,ndowna(0,il)
-                  if (ndowna(k,il).eq.i) then 
-                     ndowna(k,il)=iex
-                  elseif (ndowna(k,il).eq.iex) then 
-                     ndowna(k,il)=i
-                  endif
-                 enddo
-                enddo
-
-                iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
-                itmp=i2rep(i-1)
-                i2rep(i-1)=i2rep(iex-1)
-                i2rep(iex-1)=itmp
-
-c                write(iout,*) 'exchange',i,iex
-c                write (iout,'(a8,100i4)') "@ ifirst",
-c     &                    (ifirst(k),k=1,remd_m(1))
-c                do il=1,nodes
-c                 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
-c     &                    (nupa(k,il),k=1,nupa(0,il))
-c                 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
-c     &                    (ndowna(k,il),k=1,ndowna(0,il))
-c                enddo
-c                call flush(iout) 
-
-              else
-               remd_ene(0,iex)=ene_iex_iex
-               remd_ene(0,i)=ene_i_i
-               i=iex
-              endif 
-            endif
-           enddo
-           enddo
-cd           write (iout,*) "exchange completed"
-cd           call flush(iout) 
-        ELSEIF (usampl) THEN
-          do ii=1,nodes  
-cd            write(iout,*) "########",ii
-
-            i_temp=iran_num(1,nrep)
-            i_mult=iran_num(1,remd_m(i_temp))
-            i_iset=iran_num(1,nset)
-            i_mset=iran_num(1,mset(i_iset))
-            i=i_index(i_temp,i_mult,i_iset,i_mset)
-
-cd            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
-
-            if(t_exchange_only)then
-             i_dir=1
-            else
-             i_dir=iran_num(1,3)
-            endif
-cd            write(iout,*) "i_dir=",i_dir
-
-            if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then            
-               
-               i_temp1=i_temp+1
-               i_mult1=iran_num(1,remd_m(i_temp1))
-               i_iset1=i_iset
-               i_mset1=iran_num(1,mset(i_iset1))
-               iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-
-            elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
-
-               i_temp1=i_temp
-               i_mult1=iran_num(1,remd_m(i_temp1))
-               i_iset1=i_iset+1
-               i_mset1=iran_num(1,mset(i_iset1))
-               iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-               econstr_temp_i=remd_ene(20,i)
-               econstr_temp_iex=remd_ene(20,iex)
-               remd_ene(20,i)=remd_ene(n_ene+3,i)
-               remd_ene(20,iex)=remd_ene(n_ene+4,iex)
-
-            elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
-
-               i_temp1=i_temp+1
-               i_mult1=iran_num(1,remd_m(i_temp1))
-               i_iset1=i_iset+1
-               i_mset1=iran_num(1,mset(i_iset1))
-               iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-               econstr_temp_i=remd_ene(20,i)
-               econstr_temp_iex=remd_ene(20,iex)
-               remd_ene(20,i)=remd_ene(n_ene+3,i)
-               remd_ene(20,iex)=remd_ene(n_ene+4,iex)
-
-            else
-               goto 444 
-            endif
-cd            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
-ctime            call flush(iout)
-
-c Swap temperatures between conformations i and iex with recalculating the free energies
-c following temperature changes.
-              ene_iex_iex=remd_ene(0,iex)
-              ene_i_i=remd_ene(0,i)
-co              write (iout,*) "rescaling weights with temperature",
-co     &          remd_t_bath(i)
-              call rescale_weights(remd_t_bath(i))
-              
-              call sum_energy(remd_ene(0,iex),.false.)
-              ene_iex_i=remd_ene(0,iex)
-cd              write (iout,*) "ene_iex_i",remd_ene(0,iex)
-c              call sum_energy(remd_ene(0,i),.false.)
-cd              write (iout,*) "ene_i_i",remd_ene(0,i)
-c              write (iout,*) "rescaling weights with temperature",
-c     &          remd_t_bath(iex)
-c              if (real(ene_i_i).ne.real(remd_ene(0,i))) then
-c                write (iout,*) "ERROR: inconsistent energies:",i,
-c     &            ene_i_i,remd_ene(0,i)
-c              endif
-              call rescale_weights(remd_t_bath(iex))
-              call sum_energy(remd_ene(0,i),.false.)
-cd              write (iout,*) "ene_i_iex",remd_ene(0,i)
-              ene_i_iex=remd_ene(0,i)
-c              call sum_energy(remd_ene(0,iex),.false.)
-c              if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
-c                write (iout,*) "ERROR: inconsistent energies:",iex,
-c     &            ene_iex_iex,remd_ene(0,iex)
-c              endif
-cd              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
-c              write (iout,*) "i",i," iex",iex
-cd              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
-cd     &           " ene_i_iex",ene_i_iex,
-cd     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
-              delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
-     &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
-              delta=-delta
-cd              write(iout,*) 'delta',delta
-c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
-c     &              (remd_ene(i)-remd_ene(iex))/Rb/
-c     &              (remd_t_bath(i)*remd_t_bath(iex))
-              if (delta .gt. 50.0d0) then
-                delta=0.0d0
-              else
-                delta=dexp(-delta)
-              endif
-              if (i_dir.eq.1.or.i_dir.eq.3)
-     &         iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
-              if (i_dir.eq.2.or.i_dir.eq.3)
-     &          iremd_tot_usa(int(i2set(i-1)))=
-     &                 iremd_tot_usa(int(i2set(i-1)))+1
-              xxx=ran_number(0.0d0,1.0d0)
-cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
-              if (delta .gt. xxx) then
-                tmp=remd_t_bath(i)       
-                remd_t_bath(i)=remd_t_bath(iex)
-                remd_t_bath(iex)=tmp
-
-                itmp=iremd_iset(i)       
-                iremd_iset(i)=iremd_iset(iex)
-                iremd_iset(iex)=itmp
-
-                remd_ene(0,i)=ene_i_iex
-                remd_ene(0,iex)=ene_iex_i
-
-                if (i_dir.eq.1.or.i_dir.eq.3) 
-     &           iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
-
-                itmp=i2rep(i-1)
-                i2rep(i-1)=i2rep(iex-1)
-                i2rep(iex-1)=itmp
-
-                if (i_dir.eq.2.or.i_dir.eq.3) 
-     &           iremd_acc_usa(int(i2set(i-1)))=
-     &                 iremd_acc_usa(int(i2set(i-1)))+1
-
-                itmp=i2set(i-1)
-                i2set(i-1)=i2set(iex-1)
-                i2set(iex-1)=itmp
-        
-                itmp=i_index(i_temp,i_mult,i_iset,i_mset)
-                i_index(i_temp,i_mult,i_iset,i_mset)=
-     &                i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-                i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
-
-              else
-               remd_ene(0,iex)=ene_iex_iex
-               remd_ene(0,i)=ene_i_i
-               remd_ene(20,iex)=econstr_temp_iex
-               remd_ene(20,i)=econstr_temp_i
-              endif
-
-cd      do il=1,nset
-cd       do il1=1,mset(il)
-cd        do i=1,nrep
-cd         do j=1,remd_m(i)
-cd          write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
-cd         enddo
-cd        enddo
-cd       enddo
-cd      enddo
-
- 444      continue           
-
-          enddo
-
-        ELSEIF (hremd.gt.0) THEN
-          do ii=1,nodes  
-cd            write(iout,*) "########",ii
-
-            i_temp=iran_num(1,nrep)
-            i_mult=iran_num(1,remd_m(i_temp))
-            i_iset=iran_num(1,nset)
-            i_mset=1
-            i=i_index(i_temp,i_mult,i_iset,i_mset)
-
-cd            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
-
-            if(t_exchange_only)then
-             i_dir=1
-            else
-             i_dir=iran_num(1,3)
-            endif
-
-cd            write(iout,*) "i_dir=",i_dir
-
-            if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then            
-               
-               i_temp1=i_temp+1
-               i_mult1=iran_num(1,remd_m(i_temp1))
-               i_iset1=i_iset
-               i_mset1=1
-               iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-
-            elseif(i_dir.eq.2)then
-
-               i_temp1=i_temp
-               i_mult1=iran_num(1,remd_m(i_temp1))
-               i_iset1=iran_num(1,hremd)
-               do while(i_iset1.eq.i_iset)
-                 i_iset1=iran_num(1,hremd)
-               enddo
-               i_mset1=1
-               iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-
-            elseif(remd_m(i_temp+1).gt.0)then
-
-               i_temp1=i_temp+1
-               i_mult1=iran_num(1,remd_m(i_temp1))
-               i_iset1=iran_num(1,hremd)
-               do while(i_iset1.eq.i_iset)
-                 i_iset1=iran_num(1,hremd)
-               enddo
-               i_mset1=1
-               iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-
-            else
-               goto 445 
-            endif
-cd            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
-ctime            call flush(iout)
-
-c Swap temperatures between conformations i and iex with recalculating the free energies
-c following temperature changes.
-              ene_iex_iex=remd_ene(0,iex)
-              ene_i_i=remd_ene(0,i)
-
-              call set_hweights(i_iset)
-              call rescale_weights(remd_t_bath(i))
-              call sum_energy(remd_ene(0,iex),.false.)
-              ene_iex_i=remd_ene(0,iex)
-
-              call set_hweights(i_iset1)
-              call rescale_weights(remd_t_bath(iex))
-              call sum_energy(remd_ene(0,i),.false.)
-              ene_i_iex=remd_ene(0,i)
-
-cd              write(iout,*)  ene_iex_iex,ene_i_i,ene_iex_i,ene_i_iex
-
-              delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
-     &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
-              delta=-delta
-
-              if (delta .gt. 50.0d0) then
-                delta=0.0d0
-              else
-                delta=dexp(-delta)
-              endif
-
-              if (i_dir.eq.1.or.i_dir.eq.3)
-     &         iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
-              if (i_dir.eq.2.or.i_dir.eq.3)
-     &          iremd_tot_usa(int(i2set(i-1)))=
-     &                 iremd_tot_usa(int(i2set(i-1)))+1
-              xxx=ran_number(0.0d0,1.0d0)
-cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
-              if (delta .gt. xxx) then
-
-cd                write (iout,*) "exchange"
-                tmp=remd_t_bath(i)       
-                remd_t_bath(i)=remd_t_bath(iex)
-                remd_t_bath(iex)=tmp
-
-                itmp=iremd_iset(i)       
-                iremd_iset(i)=iremd_iset(iex)
-                iremd_iset(iex)=itmp
-
-                remd_ene(0,i)=ene_i_iex
-                remd_ene(0,iex)=ene_iex_i
-
-                if (i_dir.eq.1.or.i_dir.eq.3) 
-     &           iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
-
-                itmp=i2rep(i-1)
-                i2rep(i-1)=i2rep(iex-1)
-                i2rep(iex-1)=itmp
-
-                if (i_dir.eq.2.or.i_dir.eq.3) 
-     &           iremd_acc_usa(int(i2set(i-1)))=
-     &                 iremd_acc_usa(int(i2set(i-1)))+1
-
-                itmp=i2set(i-1)
-                i2set(i-1)=i2set(iex-1)
-                i2set(iex-1)=itmp
-        
-                itmp=i_index(i_temp,i_mult,i_iset,i_mset)
-                i_index(i_temp,i_mult,i_iset,i_mset)=
-     &                i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-                i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
-
-cd       do il=1,nset
-cd        do il1=1,mset(il)
-cd         do i=1,nrep
-cd          do j=1,remd_m(i)
-cd            write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
-cd          enddo
-cd         enddo
-cd        enddo
-cd       enddo
-
-              else
-               remd_ene(0,iex)=ene_iex_iex
-               remd_ene(0,i)=ene_i_i
-              endif
-
-
-
- 445      continue           
-
-          enddo
-
-        ENDIF
-
-c-------------------------------------
-             write (iout,*) "NREP",nrep
-             do i=1,nrep
-              if(iremd_tot(i).ne.0)
-     &          write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
-     &           ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
-             enddo
-
-             if(usampl) then
-              do i=1,nset
-               if(iremd_tot_usa(i).ne.0)
-     &           write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
-     &         iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
-              enddo
-             endif
-
-             if(hremd.gt.0) then
-              do i=1,nset
-               if(iremd_tot_usa(i).ne.0)
-     &           write(iout,'(a10,i4,f12.5,i8)') 'ACC_hremd',i,
-     &         iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
-              enddo
-             endif
-
-
-ctime             call flush(iout)
-
-cd              write (iout,'(a6,100i4)') "ifirst",
-cd     &                    (ifirst(i),i=1,remd_m(1))
-cd              do il=1,nodes
-cd               write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
-cd     &                    (nupa(i,il),i=1,nupa(0,il))
-cd               write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
-cd     &                    (ndowna(i,il),i=1,ndowna(0,il))
-cd              enddo
-            endif
-
-         time06=MPI_WTIME()
-cd         write (iout,*) "Before scatter"
-cd         call flush(iout)
-         call mpi_scatter(remd_t_bath,1,mpi_double_precision,
-     &           t_bath,1,mpi_double_precision,king,
-     &           CG_COMM,ierr) 
-cd         write (iout,*) "After scatter"
-cd         call flush(iout)
-         if(usampl.or.hremd.gt.0)
-     &    call mpi_scatter(iremd_iset,1,mpi_integer,
-     &           iset,1,mpi_integer,king,
-     &           CG_COMM,ierr) 
-
-         time07=MPI_WTIME()
-          if (me.eq.king .or. .not. out1file) then
-            write(iout,*) 'REMD scatter time=',time07-time06
-          endif
-
-         if(lmuca) then
-           call mpi_scatter(elowi,1,mpi_double_precision,
-     &           elow,1,mpi_double_precision,king,
-     &           CG_COMM,ierr) 
-           call mpi_scatter(ehighi,1,mpi_double_precision,
-     &           ehigh,1,mpi_double_precision,king,
-     &           CG_COMM,ierr) 
-         endif
-
-         if(hremd.gt.0) call set_hweights(iset)
-         call rescale_weights(t_bath)
-co         write (iout,*) "Processor",me,
-co     &    " rescaling weights with temperature",t_bath
-
-         stdfp=dsqrt(2*Rb*t_bath/d_time)
-         do i=1,ntyp
-           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
-         enddo
-         if (lang.gt.0) then
-           do i=nnt,nct-1
-            stdforcp(i)=stdforcp(i)*sqrt(t_bath/t_bath_old)
-           enddo
-           do i=nnt,nct
-            stdforcsc(i)=stdforcsc(i)*sqrt(t_bath/t_bath_old)
-           enddo
-         endif
-cde         write(iout,*) 'REMD after',me,t_bath
-           time08=MPI_WTIME()
-           if (me.eq.king .or. .not. out1file) then
-            write(iout,*) 'REMD exchange time=',time08-time02
-ctime            call flush(iout)
-           endif
-        endif
-      enddo
-
-      if (restart1file) then 
-          if (me.eq.king .or. .not. out1file)
-     &      write(iout,*) 'writing restart at the end of run'
-           call write1rst(i_index)
-      endif
-
-      if (traj1file) call write1traj
-cd debugging
-cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
-cdeb     &             icache_all,1,mpi_integer,king,
-cdeb     &             CG_COMM,ierr)
-cdeb            write(iout,'(a40,8000i8)') 
-cdeb     &             '  ntwx_cache after traj1file at the end',
-cdeb     &             (icache_all(i),i=1,nodes)
-cd end
-
-
-#ifdef MPI
-      t_MD=MPI_Wtime()-tt0
-#else
-      t_MD=tcpu()-tt0
-#endif
-      if (me.eq.king .or. .not. out1file) then
-       write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') 
-     &  '  Timing  ',
-     & 'MD calculations setup:',t_MDsetup,
-     & 'Energy & gradient evaluation:',t_enegrad,
-     & 'Stochastic MD setup:',t_langsetup,
-     & 'Stochastic MD step setup:',t_sdsetup,
-     & 'MD steps:',t_MD
-       write (iout,'(/28(1h=),a25,27(1h=))') 
-     & '  End of MD calculation  '
-      endif
-      return
-      end
-
-c-----------------------------------------------------------------------
-      subroutine write1rst(i_index)
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'mpif.h'
-      include 'COMMON.MD'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.REMD'
-      include 'COMMON.SETUP'
-      include 'COMMON.CHAIN'
-      include 'COMMON.SBRIDGE'
-      include 'COMMON.INTERACT'
-               
-      real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
-     &     d_restart2(3,2*maxres*maxprocs)
-      real t5_restart1(5)
-      integer iret,itmp
-      integer*2 i_index
-     &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
-       common /przechowalnia/ d_restart1,d_restart2
-
-       t5_restart1(1)=totT
-       t5_restart1(2)=EK
-       t5_restart1(3)=potE
-       t5_restart1(4)=t_bath
-       t5_restart1(5)=Uconst
-       
-       call mpi_gather(t5_restart1,5,mpi_real,
-     &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
-
-
-       do i=1,2*nres
-         do j=1,3
-           r_d(j,i)=d_t(j,i)
-         enddo
-       enddo
-       call mpi_gather(r_d,3*2*nres,mpi_real,
-     &           d_restart1,3*2*nres,mpi_real,king,
-     &           CG_COMM,ierr)
-
-
-       do i=1,2*nres
-         do j=1,3
-           r_d(j,i)=dc(j,i)
-         enddo
-       enddo
-       call mpi_gather(r_d,3*2*nres,mpi_real,
-     &           d_restart2,3*2*nres,mpi_real,king,
-     &           CG_COMM,ierr)
-
-       if(me.eq.king) then
-#ifdef AIX
-         call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
-         do i=0,nodes-1
-          call xdrfint_(ixdrf, i2rep(i), iret)
-         enddo
-         do i=1,remd_m(1)
-          call xdrfint_(ixdrf, ifirst(i), iret)
-         enddo
-         do il=1,nodes
-              do i=0,nupa(0,il)
-               call xdrfint_(ixdrf, nupa(i,il), iret)
-              enddo
-
-              do i=0,ndowna(0,il)
-               call xdrfint_(ixdrf, ndowna(i,il), iret)
-              enddo
-         enddo
-
-         do il=1,nodes
-           do j=1,4
-            call xdrffloat_(ixdrf, t_restart1(j,il), iret)
-           enddo
-         enddo
-
-         do il=0,nodes-1
-           do i=1,2*nres
-            do j=1,3
-             call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
-            enddo
-           enddo
-         enddo
-         do il=0,nodes-1
-           do i=1,2*nres
-            do j=1,3
-             call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
-            enddo
-           enddo
-         enddo
-
-         if(usampl) then
-           call xdrfint_(ixdrf, nset, iret)
-           do i=1,nset
-             call xdrfint_(ixdrf,mset(i), iret)
-           enddo
-           do i=0,nodes-1
-             call xdrfint_(ixdrf,i2set(i), iret)
-           enddo
-           do il=1,nset
-             do il1=1,mset(il)
-               do i=1,nrep
-                 do j=1,remd_m(i)
-                   itmp=i_index(i,j,il,il1)
-                   call xdrfint_(ixdrf,itmp, iret)
-                 enddo
-               enddo
-             enddo
-           enddo
-           
-         endif
-         call xdrfclose_(ixdrf, iret)
-#else
-         call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
-         do i=0,nodes-1
-          call xdrfint(ixdrf, i2rep(i), iret)
-         enddo
-         do i=1,remd_m(1)
-          call xdrfint(ixdrf, ifirst(i), iret)
-         enddo
-         do il=1,nodes
-              do i=0,nupa(0,il)
-               call xdrfint(ixdrf, nupa(i,il), iret)
-              enddo
-
-              do i=0,ndowna(0,il)
-               call xdrfint(ixdrf, ndowna(i,il), iret)
-              enddo
-         enddo
-
-         do il=1,nodes
-           do j=1,4
-            call xdrffloat(ixdrf, t_restart1(j,il), iret)
-           enddo
-         enddo
-
-         do il=0,nodes-1
-           do i=1,2*nres
-            do j=1,3
-             call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
-            enddo
-           enddo
-         enddo
-         do il=0,nodes-1
-           do i=1,2*nres
-            do j=1,3
-             call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
-            enddo
-           enddo
-         enddo
-
-
-             if(usampl) then
-              call xdrfint(ixdrf, nset, iret)
-              do i=1,nset
-                call xdrfint(ixdrf,mset(i), iret)
-              enddo
-              do i=0,nodes-1
-                call xdrfint(ixdrf,i2set(i), iret)
-              enddo
-              do il=1,nset
-               do il1=1,mset(il)
-                do i=1,nrep
-                 do j=1,remd_m(i)
-                   itmp=i_index(i,j,il,il1)
-                   call xdrfint(ixdrf,itmp, iret)
-                 enddo
-                enddo
-               enddo
-              enddo
-           
-             endif
-         call xdrfclose(ixdrf, iret)
-#endif
-       endif
-      return
-      end
-
-
-      subroutine write1traj
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'mpif.h'
-      include 'COMMON.MD'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.REMD'
-      include 'COMMON.SETUP'
-      include 'COMMON.CHAIN'
-      include 'COMMON.SBRIDGE'
-      include 'COMMON.INTERACT'
-               
-      real t5_restart1(5)
-      integer iret,itmp
-      real xcoord(3,maxres2+2),prec
-      real r_qfrag(50),r_qpair(100)
-      real r_utheta(50),r_ugamma(100),r_uscdiff(100)
-      real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
-      real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
-     &     p_uscdiff(100*maxprocs)
-      real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
-      common /przechowalnia/ p_c
-
-      call mpi_bcast(ii_write,1,mpi_integer,
-     &           king,CG_COMM,ierr)
-
-c debugging
-      print *,'traj1file',me,ii_write,ntwx_cache
-c end debugging
-
-#ifdef AIX
-      if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
-#else
-      if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
-#endif
-      do ii=1,ii_write
-       t5_restart1(1)=totT_cache(ii)
-       t5_restart1(2)=EK_cache(ii)
-       t5_restart1(3)=potE_cache(ii)
-       t5_restart1(4)=t_bath_cache(ii)
-       t5_restart1(5)=Uconst_cache(ii)
-       call mpi_gather(t5_restart1,5,mpi_real,
-     &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
-
-       call mpi_gather(iset_cache(ii),1,mpi_integer,
-     &      iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
-
-          do i=1,nfrag
-           r_qfrag(i)=qfrag_cache(i,ii)
-          enddo
-          do i=1,npair
-           r_qpair(i)=qpair_cache(i,ii)
-          enddo
-          do i=1,nfrag_back
-           r_utheta(i)=utheta_cache(i,ii)
-           r_ugamma(i)=ugamma_cache(i,ii)
-           r_uscdiff(i)=uscdiff_cache(i,ii)
-          enddo
-
-        call mpi_gather(r_qfrag,nfrag,mpi_real,
-     &           p_qfrag,nfrag,mpi_real,king,
-     &           CG_COMM,ierr)
-        call mpi_gather(r_qpair,npair,mpi_real,
-     &           p_qpair,npair,mpi_real,king,
-     &           CG_COMM,ierr)
-        call mpi_gather(r_utheta,nfrag_back,mpi_real,
-     &           p_utheta,nfrag_back,mpi_real,king,
-     &           CG_COMM,ierr)
-        call mpi_gather(r_ugamma,nfrag_back,mpi_real,
-     &           p_ugamma,nfrag_back,mpi_real,king,
-     &           CG_COMM,ierr)
-        call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
-     &           p_uscdiff,nfrag_back,mpi_real,king,
-     &           CG_COMM,ierr)
-
-#ifdef DEBUG
-        write (iout,*) "p_qfrag"
-        do i=1,nodes
-          write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
-        enddo
-        write (iout,*) "p_qpair"
-        do i=1,nodes
-          write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
-        enddo
-ctime        call flush(iout)
-#endif
-        do i=1,nres*2
-         do j=1,3
-          r_c(j,i)=c_cache(j,i,ii)
-         enddo
-        enddo
-
-        call mpi_gather(r_c,3*2*nres,mpi_real,
-     &           p_c,3*2*nres,mpi_real,king,
-     &           CG_COMM,ierr)
-
-       if(me.eq.king) then
-#ifdef AIX
-         do il=1,nodes
-          call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
-          call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
-          call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
-          call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
-          call xdrfint_(ixdrf, nss, iret) 
-          do j=1,nss
-           call xdrfint_(ixdrf, ihpb(j), iret)
-           call xdrfint_(ixdrf, jhpb(j), iret)
-          enddo
-          call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
-          call xdrfint_(ixdrf, iset_restart1(il), iret)
-          do i=1,nfrag
-           call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
-          enddo
-          do i=1,npair
-           call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
-          enddo
-          do i=1,nfrag_back
-           call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
-           call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
-           call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
-          enddo
-          prec=10000.0
-          do i=1,nres
-           do j=1,3
-            xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
-           enddo
-          enddo
-          do i=nnt,nct
-           do j=1,3
-            xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
-           enddo
-          enddo
-          itmp=nres+nct-nnt+1
-          call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
-         enddo
-#else
-         do il=1,nodes
-          call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
-          call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
-          call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
-          call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
-          call xdrfint(ixdrf, nss, iret) 
-          do j=1,nss
-           call xdrfint(ixdrf, ihpb(j), iret)
-           call xdrfint(ixdrf, jhpb(j), iret)
-          enddo
-          call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
-          call xdrfint(ixdrf, iset_restart1(il), iret)
-          do i=1,nfrag
-           call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
-          enddo
-          do i=1,npair
-           call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
-          enddo
-          do i=1,nfrag_back
-           call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
-           call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
-           call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
-          enddo
-          prec=10000.0
-          do i=1,nres
-           do j=1,3
-            xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
-           enddo
-          enddo
-          do i=nnt,nct
-           do j=1,3
-            xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
-           enddo
-          enddo
-          itmp=nres+nct-nnt+1
-          call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
-         enddo
-#endif
-       endif
-      enddo
-#ifdef AIX
-      if(me.eq.king) call xdrfclose_(ixdrf, iret)
-#else
-      if(me.eq.king) call xdrfclose(ixdrf, iret)
-#endif
-      do i=1,ntwx_cache-ii_write
-
-            totT_cache(i)=totT_cache(ii_write+i)
-            EK_cache(i)=EK_cache(ii_write+i)
-            potE_cache(i)=potE_cache(ii_write+i)
-            t_bath_cache(i)=t_bath_cache(ii_write+i)
-            Uconst_cache(i)=Uconst_cache(ii_write+i)
-            iset_cache(i)=iset_cache(ii_write+i)
-
-            do ii=1,nfrag
-             qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
-            enddo
-            do ii=1,npair
-             qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
-            enddo
-            do ii=1,nfrag_back
-              utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
-              ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
-              uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
-            enddo
-
-            do ii=1,nres*2
-             do j=1,3
-              c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
-             enddo
-            enddo
-      enddo
-      ntwx_cache=ntwx_cache-ii_write
-      return
-      end
-
-
-      subroutine read1restart(i_index)
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'mpif.h'
-      include 'COMMON.MD'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.REMD'
-      include 'COMMON.SETUP'
-      include 'COMMON.CHAIN'
-      include 'COMMON.SBRIDGE'
-      include 'COMMON.INTERACT'
-      real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
-     &                 t5_restart1(5)
-      integer*2 i_index
-     &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
-      common /przechowalnia/ d_restart1
-      write (*,*) "Processor",me," called read1restart"
-
-         if(me.eq.king)then
-              open(irest2,file=mremd_rst_name,status='unknown')
-              read(irest2,*,err=334) i
-              write(iout,*) "Reading old rst in ASCI format"
-              close(irest2)
-               call read1restart_old
-               return
- 334          continue
-#ifdef AIX
-              call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
-
-              do i=0,nodes-1
-               call xdrfint_(ixdrf, i2rep(i), iret)
-              enddo
-              do i=1,remd_m(1)
-               call xdrfint_(ixdrf, ifirst(i), iret)
-              enddo
-             do il=1,nodes
-              call xdrfint_(ixdrf, nupa(0,il), iret)
-              do i=1,nupa(0,il)
-               call xdrfint_(ixdrf, nupa(i,il), iret)
-              enddo
-
-              call xdrfint_(ixdrf, ndowna(0,il), iret)
-              do i=1,ndowna(0,il)
-               call xdrfint_(ixdrf, ndowna(i,il), iret)
-              enddo
-             enddo
-             do il=1,nodes
-               do j=1,4
-                call xdrffloat_(ixdrf, t_restart1(j,il), iret)
-               enddo
-             enddo
-#else
-              call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
-
-              do i=0,nodes-1
-               call xdrfint(ixdrf, i2rep(i), iret)
-              enddo
-              do i=1,remd_m(1)
-               call xdrfint(ixdrf, ifirst(i), iret)
-              enddo
-             do il=1,nodes
-              call xdrfint(ixdrf, nupa(0,il), iret)
-              do i=1,nupa(0,il)
-               call xdrfint(ixdrf, nupa(i,il), iret)
-              enddo
-
-              call xdrfint(ixdrf, ndowna(0,il), iret)
-              do i=1,ndowna(0,il)
-               call xdrfint(ixdrf, ndowna(i,il), iret)
-              enddo
-             enddo
-             do il=1,nodes
-               do j=1,4
-                call xdrffloat(ixdrf, t_restart1(j,il), iret)
-               enddo
-             enddo
-#endif
-         endif
-         call mpi_scatter(t_restart1,5,mpi_real,
-     &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
-         totT=t5_restart1(1)              
-         EK=t5_restart1(2)
-         potE=t5_restart1(3)
-         t_bath=t5_restart1(4)
-
-         if(me.eq.king)then
-              do il=0,nodes-1
-               do i=1,2*nres
-c                read(irest2,'(3e15.5)') 
-c     &                (d_restart1(j,i+2*nres*il),j=1,3)
-            do j=1,3
-#ifdef AIX
-             call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
-#else
-             call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
-#endif
-            enddo
-               enddo
-              enddo
-         endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,
-     &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
-
-         do i=1,2*nres
-           do j=1,3
-            d_t(j,i)=r_d(j,i)
-           enddo
-         enddo
-         if(me.eq.king)then 
-              do il=0,nodes-1
-               do i=1,2*nres
-c                read(irest2,'(3e15.5)') 
-c     &                (d_restart1(j,i+2*nres*il),j=1,3)
-            do j=1,3
-#ifdef AIX
-             call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
-#else
-             call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
-#endif
-            enddo
-               enddo
-              enddo
-         endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,
-     &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
-         do i=1,2*nres
-           do j=1,3
-            dc(j,i)=r_d(j,i)
-           enddo
-         enddo
-       
-
-           if(usampl) then
-#ifdef AIX
-             if(me.eq.king)then
-              call xdrfint_(ixdrf, nset, iret)
-              do i=1,nset
-                call xdrfint_(ixdrf,mset(i), iret)
-              enddo
-              do i=0,nodes-1
-                call xdrfint_(ixdrf,i2set(i), iret)
-              enddo
-              do il=1,nset
-               do il1=1,mset(il)
-                do i=1,nrep
-                 do j=1,remd_m(i)
-                   call xdrfint_(ixdrf,itmp, iret)
-                   i_index(i,j,il,il1)=itmp
-                 enddo
-                enddo
-               enddo
-              enddo
-             endif
-#else
-             if(me.eq.king)then
-              call xdrfint(ixdrf, nset, iret)
-              do i=1,nset
-                call xdrfint(ixdrf,mset(i), iret)
-              enddo
-              do i=0,nodes-1
-                call xdrfint(ixdrf,i2set(i), iret)
-              enddo
-              do il=1,nset
-               do il1=1,mset(il)
-                do i=1,nrep
-                 do j=1,remd_m(i)
-                   call xdrfint(ixdrf,itmp, iret)
-                   i_index(i,j,il,il1)=itmp
-                 enddo
-                enddo
-               enddo
-              enddo
-             endif
-#endif
-              call mpi_scatter(i2set,1,mpi_integer,
-     &           iset,1,mpi_integer,king,
-     &           CG_COMM,ierr) 
-
-           endif
-
-
-        if(me.eq.king) close(irest2)
-        return
-        end
-
-      subroutine read1restart_old
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'mpif.h'
-      include 'COMMON.MD'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.REMD'
-      include 'COMMON.SETUP'
-      include 'COMMON.CHAIN'
-      include 'COMMON.SBRIDGE'
-      include 'COMMON.INTERACT'
-      real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
-     &                 t5_restart1(5)
-      common /przechowalnia/ d_restart1
-         if(me.eq.king)then
-             open(irest2,file=mremd_rst_name,status='unknown')
-             read (irest2,*) (i2rep(i),i=0,nodes-1)
-             read (irest2,*) (ifirst(i),i=1,remd_m(1))
-             do il=1,nodes
-              read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
-              read (irest2,*) ndowna(0,il),
-     &                    (ndowna(i,il),i=1,ndowna(0,il))
-             enddo
-             do il=1,nodes
-               read(irest2,*) (t_restart1(j,il),j=1,4)
-             enddo
-         endif
-         call mpi_scatter(t_restart1,5,mpi_real,
-     &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
-         totT=t5_restart1(1)              
-         EK=t5_restart1(2)
-         potE=t5_restart1(3)
-         t_bath=t5_restart1(4)
-
-         if(me.eq.king)then
-              do il=0,nodes-1
-               do i=1,2*nres
-                read(irest2,'(3e15.5)') 
-     &                (d_restart1(j,i+2*nres*il),j=1,3)
-               enddo
-              enddo
-         endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,
-     &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
-
-         do i=1,2*nres
-           do j=1,3
-            d_t(j,i)=r_d(j,i)
-           enddo
-         enddo
-         if(me.eq.king)then 
-              do il=0,nodes-1
-               do i=1,2*nres
-                read(irest2,'(3e15.5)') 
-     &                (d_restart1(j,i+2*nres*il),j=1,3)
-               enddo
-              enddo
-         endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,
-     &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
-         do i=1,2*nres
-           do j=1,3
-            dc(j,i)=r_d(j,i)
-           enddo
-         enddo
-        if(me.eq.king) close(irest2)
-        return
-        end
-c-------------------------------------------------------------------
-        subroutine set_hweights(iiset)          
-        implicit real*8 (a-h,o-z)
-        integer i  
-        include 'DIMENSIONS'    
-        include 'COMMON.FFIELD'
-        include 'COMMON.REMD'    
-
-         do i=1,n_ene
-          weights(i)=hweights(iiset,i)
-         enddo
-
-         wsc    =weights(1) 
-         wscp   =weights(2) 
-         welec  =weights(3) 
-         wcorr  =weights(4) 
-         wcorr5 =weights(5) 
-         wcorr6 =weights(6) 
-         wel_loc=weights(7) 
-         wturn3 =weights(8) 
-         wturn4 =weights(9) 
-         wturn6 =weights(10)
-         wang   =weights(11)
-         wscloc =weights(12)
-         wtor   =weights(13)
-         wtor_d =weights(14)
-         wstrain=weights(15)
-         wvdwpp =weights(16)
-         wbond  =weights(17)
-         scal14 =weights(18)
-         wsccor =weights(21)
-
-        return
-        end
-#endif