source/unres/src_MD-M/MREMD.F

   1 c#define DEBUG
   2       subroutine MREMD
   3       implicit real*8 (a-h,o-z)
   4       include 'DIMENSIONS'
   5       include 'mpif.h'
   6       include 'COMMON.CONTROL'
   7       include 'COMMON.VAR'
   8       include 'COMMON.MD'
   9 #ifndef LANG0
  10       include 'COMMON.LANGEVIN'
  11 #else
  12       include 'COMMON.LANGEVIN.lang0'
  13 #endif
  14       include 'COMMON.CHAIN'
  15       include 'COMMON.DERIV'
  16       include 'COMMON.GEO'
  17       include 'COMMON.LOCAL'
  18       include 'COMMON.INTERACT'
  19       include 'COMMON.IOUNITS'
  20       include 'COMMON.NAMES'
  21       include 'COMMON.TIME1'
  22       include 'COMMON.REMD'
  23       include 'COMMON.SETUP'
  24       include 'COMMON.MUCA'
  25       include 'COMMON.HAIRPIN'
  26       integer ERRCODE
  27       double precision cm(3),L(3),vcm(3)
  28       double precision energia(0:n_ene)
  29       double precision remd_t_bath(maxprocs)
  30       integer iremd_iset(maxprocs)
  31       integer*2 i_index
  32      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
  33       double precision remd_ene(0:n_ene+4,maxprocs),t_bath_old,e_tmp
  34       integer iremd_acc(maxprocs),iremd_tot(maxprocs)
  35       integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
  36       integer ilen,rstcount
  37       external ilen
  38       character*50 tytul
  39       common /gucio/ cm
  40       integer itime
  41 cold      integer nup(0:maxprocs),ndown(0:maxprocs)
  42       integer rep2i(0:maxprocs),ireqi(maxprocs)
  43       integer icache_all(maxprocs)
  44       integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
  45       logical synflag,end_of_run,file_exist /.false./,ovrtim
  46       real ene_tol /1.0e-5/
  47
  48 cdeb      imin_itime_old=0
  49       ntwx_cache=0
  50       time00=MPI_WTIME()
  51       time01=time00
  52       if(me.eq.king.or..not.out1file) then
  53        write  (iout,*) 'MREMD',nodes,'time before',time00-walltime
  54        write (iout,*) "NREP=",nrep
  55       endif
  56
  57       synflag=.false.
  58       if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
  59         call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
  60       endif
  61       mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
  62
  63 cd      print *,'MREMD',nodes,homol_nset
  64 cd      print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
  65 cde      write (iout,*) "Start MREMD: me",me," t_bath",t_bath
  66       if(homol_nset.gt.1) then
  67          i_econstr=24
  68          nset=homol_nset
  69          do i=1,nset
  70           mset(i)=1
  71          enddo
  72       endif
  73
  74       if(usampl) i_econstr=20
  75
  76       k=0
  77       rep2i(k)=-1
  78       do il=1,max0(nset,1)
  79        do il1=1,max0(mset(il),1)
  80         do i=1,nrep
  81          iremd_acc(i)=0
  82          iremd_acc_usa(i)=0
  83          iremd_tot(i)=0
  84          do j=1,remd_m(i)
  85           i2rep(k)=i
  86           i2set(k)=il
  87           rep2i(i)=k
  88           k=k+1
  89           i_index(i,j,il,il1)=k
  90          enddo
  91         enddo
  92        enddo
  93       enddo
  94
  95       if(me.eq.king.or..not.out1file) then
  96        write(iout,*) "i2rep",(i2rep(i),i=0,nodes-1)
  97        write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
  98        write(iout,*) "i,j,il,il1,i_index(i,j,il,il1)"
  99        do il=1,nset
 100         do il1=1,mset(il)
 101          do i=1,nrep
 102           do j=1,remd_m(i)
 103            write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
 104           enddo
 105          enddo
 106         enddo
 107        enddo
 108       endif
 109
 110 c      print *,'i2rep',me,i2rep(me)
 111 c      print *,'rep2i',(rep2i(i),i=0,nrep)
 112
 113 cold       if (i2rep(me).eq.nrep) then
 114 cold        nup(0)=0
 115 cold       else
 116 cold        nup(0)=remd_m(i2rep(me)+1)
 117 cold        k=rep2i(int(i2rep(me)))+1
 118 cold        do i=1,nup(0)
 119 cold         nup(i)=k
 120 cold         k=k+1
 121 cold        enddo
 122 cold       endif
 123
 124 cd       print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
 125
 126 cold       if (i2rep(me).eq.1) then
 127 cold        ndown(0)=0
 128 cold       else
 129 cold        ndown(0)=remd_m(i2rep(me)-1)
 130 cold        k=rep2i(i2rep(me)-2)+1
 131 cold        do i=1,ndown(0)
 132 cold         ndown(i)=k
 133 cold         k=k+1
 134 cold        enddo
 135 cold       endif
 136
 137 cd       print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
 138
 139
 140        write (*,*) "Processor",me," rest",rest,"
 141      &   restart1fie",restart1file
 142        if(rest.and.restart1file) then
 143            if (me.eq.king)
 144      &     inquire(file=mremd_rst_name,exist=file_exist)
 145 cd           write (*,*) me," Before broadcast: file_exist",file_exist
 146            call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
 147      &          IERR)
 148 cd           write (*,*) me," After broadcast: file_exist",file_exist
 149            if(file_exist) then
 150              if(me.eq.king.or..not.out1file)
 151      &            write  (iout,*) 'Master is reading restart1file'
 152              call read1restart(i_index)
 153            else
 154              if(me.eq.king.or..not.out1file)
 155      &            write  (iout,*) 'WARNING : no restart1file'
 156            endif
 157
 158            if(me.eq.king.or..not.out1file) then
 159               write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
 160               write(iout,*) "i_index"
 161               do il=1,nset
 162                do il1=1,mset(il)
 163                 do i=1,nrep
 164                  do j=1,remd_m(i)
 165                   write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
 166                  enddo
 167                 enddo
 168                enddo
 169               enddo
 170            endif
 171        endif
 172
 173        if(me.eq.king) then
 174         if (rest.and..not.restart1file)
 175      &          inquire(file=mremd_rst_name,exist=file_exist)
 176         if(.not.file_exist.and.rest.and..not.restart1file)
 177      &       write(iout,*) 'WARNING : no restart file',mremd_rst_name
 178         IF (rest.and.file_exist.and..not.restart1file) THEN
 179              write  (iout,*) 'Master is reading restart file',
 180      &                        mremd_rst_name
 181              open(irest2,file=mremd_rst_name,status='unknown')
 182              read (irest2,*)
 183              read (irest2,*) (i2rep(i),i=0,nodes-1)
 184              read (irest2,*)
 185              read (irest2,*) (ifirst(i),i=1,remd_m(1))
 186              do il=1,nodes
 187               read (irest2,*)
 188               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 189               read (irest2,*)
 190               read (irest2,*) ndowna(0,il),
 191      &                    (ndowna(i,il),i=1,ndowna(0,il))
 192              enddo
 193              if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
 194               read (irest2,*)
 195               read (irest2,*) nset
 196               read (irest2,*)
 197               read (irest2,*) (mset(i),i=1,nset)
 198               read (irest2,*)
 199               read (irest2,*) (i2set(i),i=0,nodes-1)
 200               read (irest2,*)
 201               do il=1,nset
 202                do il1=1,mset(il)
 203                 do i=1,nrep
 204                   read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
 205                 enddo
 206                enddo
 207               enddo
 208
 209               write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
 210               write(iout,*) "i_index"
 211               do il=1,nset
 212                do il1=1,mset(il)
 213                 do i=1,nrep
 214                  do j=1,remd_m(i)
 215                   write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
 216                  enddo
 217                 enddo
 218                enddo
 219               enddo
 220              endif
 221
 222              close(irest2)
 223
 224              write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
 225              write (iout,'(a6,1000i5)') "ifirst",
 226      &                    (ifirst(i),i=1,remd_m(1))
 227              do il=1,nodes
 228               write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
 229      &                    (nupa(i,il),i=1,nupa(0,il))
 230               write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
 231      &                    (ndowna(i,il),i=1,ndowna(0,il))
 232              enddo
 233         ELSE IF (.not.(rest.and.file_exist)) THEN
 234          do il=1,remd_m(1)
 235           ifirst(il)=il
 236          enddo
 237
 238          do il=1,nodes
 239           if (i2rep(il-1).eq.nrep) then
 240            nupa(0,il)=0
 241           else
 242            nupa(0,il)=remd_m(i2rep(il-1)+1)
 243            k=rep2i(int(i2rep(il-1)))+1
 244            do i=1,nupa(0,il)
 245             nupa(i,il)=k+1
 246             k=k+1
 247            enddo
 248           endif
 249           if (i2rep(il-1).eq.1) then
 250            ndowna(0,il)=0
 251           else
 252            ndowna(0,il)=remd_m(i2rep(il-1)-1)
 253            k=rep2i(i2rep(il-1)-2)+1
 254            do i=1,ndowna(0,il)
 255             ndowna(i,il)=k+1
 256             k=k+1
 257            enddo
 258           endif
 259          enddo
 260
 261         write (iout,'(a6,100i4)') "ifirst",
 262      &                    (ifirst(i),i=1,remd_m(1))
 263         do il=1,nodes
 264          write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
 265      &                    (nupa(i,il),i=1,nupa(0,il))
 266          write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
 267      &                    (ndowna(i,il),i=1,ndowna(0,il))
 268         enddo
 269
 270         ENDIF
 271        endif
 272 c
 273 c      t_bath=retmin+(retmax-retmin)*me/(nodes-1)
 274        if(.not.(rest.and.file_exist.and.restart1file)) then
 275          if (me .eq. king) then
 276             t_bath=retmin
 277          else
 278             t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
 279          endif
 280 cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
 281          if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
 282
 283        endif
 284        if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
 285           iset=i2set(me)
 286 c broadcast iset to slaves
 287           if (nfgtasks.gt.1) then
 288            call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
 289            call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
 290           endif
 291           if(me.eq.king.or..not.out1file)
 292      &     write(iout,*) me,"iset=",iset,"t_bath=",t_bath
 293        endif
 294 c
 295        stdfp=dsqrt(2*Rb*t_bath/d_time)
 296        do i=1,ntyp
 297           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 298        enddo
 299
 300 c      print *,'irep',me,t_bath
 301        if (.not.rest) then
 302         if (me.eq.king .or. .not. out1file)
 303      &   write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
 304         call rescale_weights(t_bath)
 305        endif
 306
 307
 308 c------copy MD--------------
 309 c  The driver for molecular dynamics subroutines
 310 c------------------------------------------------
 311       t_MDsetup=0.0d0
 312       t_langsetup=0.0d0
 313       t_MD=0.0d0
 314       t_enegrad=0.0d0
 315       t_sdsetup=0.0d0
 316       if(me.eq.king.or..not.out1file)
 317      & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
 318 #ifdef MPI
 319       tt0 = MPI_Wtime()
 320 #else
 321       tt0 = tcpu()
 322 #endif
 323 c Determine the inverse of the inertia matrix.
 324       call setup_MD_matrices
 325 c Initialize MD
 326       call init_MD
 327       if (rest) then
 328        if (me.eq.king .or. .not. out1file)
 329      &  write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
 330        stdfp=dsqrt(2*Rb*t_bath/d_time)
 331        do i=1,ntyp
 332           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 333        enddo
 334        call rescale_weights(t_bath)
 335       endif
 336
 337 #ifdef MPI
 338       t_MDsetup = MPI_Wtime()-tt0
 339 #else
 340       t_MDsetup = tcpu()-tt0
 341 #endif
 342       rstcount=0
 343 c   Entering the MD loop
 344 #ifdef MPI
 345       tt0 = MPI_Wtime()
 346 #else
 347       tt0 = tcpu()
 348 #endif
 349       if (lang.eq.2 .or. lang.eq.3) then
 350 #ifndef   LANG0
 351         call setup_fricmat
 352         if (lang.eq.2) then
 353           call sd_verlet_p_setup
 354         else
 355           call sd_verlet_ciccotti_setup
 356         endif
 357         do i=1,dimen
 358           do j=1,dimen
 359             pfric0_mat(i,j,0)=pfric_mat(i,j)
 360             afric0_mat(i,j,0)=afric_mat(i,j)
 361             vfric0_mat(i,j,0)=vfric_mat(i,j)
 362             prand0_mat(i,j,0)=prand_mat(i,j)
 363             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 364             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 365           enddo
 366         enddo
 367         flag_stoch(0)=.true.
 368         do i=1,maxflag_stoch
 369           flag_stoch(i)=.false.
 370         enddo
 371 #else
 372         write (iout,*)
 373      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 374 #ifdef MPI
 375         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 376 #endif
 377         stop
 378 #endif
 379       else if (lang.eq.1 .or. lang.eq.4) then
 380         call setup_fricmat
 381       endif
 382       time00=MPI_WTIME()
 383       if (me.eq.king .or. .not. out1file)
 384      & write(iout,*) 'Setup time',time00-walltime
 385       call flush(iout)
 386 #ifdef MPI
 387       t_langsetup=MPI_Wtime()-tt0
 388       tt0=MPI_Wtime()
 389 #else
 390       t_langsetup=tcpu()-tt0
 391       tt0=tcpu()
 392 #endif
 393       itime=0
 394       end_of_run=.false.
 395       do while(.not.end_of_run)
 396         itime=itime+1
 397         if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
 398         if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
 399         rstcount=rstcount+1
 400         if (lang.gt.0 .and. surfarea .and.
 401      &      mod(itime,reset_fricmat).eq.0) then
 402           if (lang.eq.2 .or. lang.eq.3) then
 403 #ifndef   LANG0
 404             call setup_fricmat
 405             if (lang.eq.2) then
 406               call sd_verlet_p_setup
 407             else
 408               call sd_verlet_ciccotti_setup
 409             endif
 410             do i=1,dimen
 411               do j=1,dimen
 412                 pfric0_mat(i,j,0)=pfric_mat(i,j)
 413                 afric0_mat(i,j,0)=afric_mat(i,j)
 414                 vfric0_mat(i,j,0)=vfric_mat(i,j)
 415                 prand0_mat(i,j,0)=prand_mat(i,j)
 416                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 417                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 418               enddo
 419             enddo
 420             flag_stoch(0)=.true.
 421             do i=1,maxflag_stoch
 422               flag_stoch(i)=.false.
 423             enddo
 424 #endif
 425           else if (lang.eq.1 .or. lang.eq.4) then
 426             call setup_fricmat
 427           endif
 428           write (iout,'(a,i10)')
 429      &      "Friction matrix reset based on surface area, itime",itime
 430         endif
 431         if (reset_vel .and. tbf .and. lang.eq.0
 432      &      .and. mod(itime,count_reset_vel).eq.0) then
 433           call random_vel
 434           if (me.eq.king .or. .not. out1file)
 435      &     write(iout,'(a,f20.2)')
 436      &     "Velocities reset to random values, time",totT
 437           do i=0,2*nres
 438             do j=1,3
 439               d_t_old(j,i)=d_t(j,i)
 440             enddo
 441           enddo
 442         endif
 443         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 444           call inertia_tensor
 445           call vcm_vel(vcm)
 446           do j=1,3
 447              d_t(j,0)=d_t(j,0)-vcm(j)
 448           enddo
 449           call kinetic(EK)
 450           kinetic_T=2.0d0/(dimen3*Rb)*EK
 451           scalfac=dsqrt(T_bath/kinetic_T)
 452 cd          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
 453           do i=0,2*nres
 454             do j=1,3
 455               d_t_old(j,i)=scalfac*d_t(j,i)
 456             enddo
 457           enddo
 458         endif
 459         if (lang.ne.4) then
 460           if (RESPA) then
 461 c Time-reversible RESPA algorithm
 462 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 463             call RESPA_step(itime)
 464           else
 465 c Variable time step algorithm.
 466             call velverlet_step(itime)
 467           endif
 468         else
 469 #ifdef BROWN
 470           call brown_step(itime)
 471 #else
 472           print *,"Brown dynamics not here!"
 473 #ifdef MPI
 474           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 475 #endif
 476           stop
 477 #endif
 478         endif
 479         if(ntwe.ne.0) then
 480           if (mod(itime,ntwe).eq.0) call statout(itime)
 481         endif
 482         if (mod(itime,ntwx).eq.0.and..not.traj1file) then
 483           write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
 484           if(mdpdb) then
 485              call hairpin(.true.,nharp,iharp)
 486              call secondary2(.true.)
 487              call pdbout(potE,tytul,ipdb)
 488           else
 489              call cartout(totT)
 490           endif
 491         endif
 492         if (mod(itime,ntwx).eq.0.and.traj1file) then
 493           if(ntwx_cache.lt.max_cache_traj_use) then
 494             ntwx_cache=ntwx_cache+1
 495           else
 496            if (max_cache_traj_use.ne.1)
 497      &      print *,itime,"processor ",me," over cache ",ntwx_cache
 498            do i=1,ntwx_cache-1
 499
 500             totT_cache(i)=totT_cache(i+1)
 501             EK_cache(i)=EK_cache(i+1)
 502             potE_cache(i)=potE_cache(i+1)
 503             t_bath_cache(i)=t_bath_cache(i+1)
 504             Uconst_cache(i)=Uconst_cache(i+1)
 505             iset_cache(i)=iset_cache(i+1)
 506
 507             do ii=1,nfrag
 508              qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
 509             enddo
 510             do ii=1,npair
 511              qpair_cache(ii,i)=qpair_cache(ii,i+1)
 512             enddo
 513             do ii=1,nfrag_back
 514               utheta_cache(ii,i)=utheta_cache(ii,i+1)
 515               ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
 516               uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
 517             enddo
 518
 519
 520             do ii=1,nres*2
 521              do j=1,3
 522               c_cache(j,ii,i)=c_cache(j,ii,i+1)
 523              enddo
 524             enddo
 525            enddo
 526           endif
 527
 528             totT_cache(ntwx_cache)=totT
 529             EK_cache(ntwx_cache)=EK
 530             potE_cache(ntwx_cache)=potE
 531             t_bath_cache(ntwx_cache)=t_bath
 532             Uconst_cache(ntwx_cache)=Uconst
 533             iset_cache(ntwx_cache)=iset
 534
 535             do i=1,nfrag
 536              qfrag_cache(i,ntwx_cache)=qfrag(i)
 537             enddo
 538             do i=1,npair
 539              qpair_cache(i,ntwx_cache)=qpair(i)
 540             enddo
 541             do i=1,nfrag_back
 542               utheta_cache(i,ntwx_cache)=utheta(i)
 543               ugamma_cache(i,ntwx_cache)=ugamma(i)
 544               uscdiff_cache(i,ntwx_cache)=uscdiff(i)
 545             enddo
 546 C            print *,'przed returnbox'
 547             call returnbox
 548 C            call enerprint(remd_ene(0,i))
 549             do i=1,nres*2
 550              do j=1,3
 551               c_cache(j,i,ntwx_cache)=c(j,i)
 552              enddo
 553             enddo
 554
 555         endif
 556         if ((rstcount.eq.1000.or.itime.eq.n_timestep)
 557      &                         .and..not.restart1file) then
 558
 559            if(me.eq.king) then
 560              open(irest1,file=mremd_rst_name,status='unknown')
 561              write (irest1,*) "i2rep"
 562              write (irest1,*) (i2rep(i),i=0,nodes-1)
 563              write (irest1,*) "ifirst"
 564              write (irest1,*) (ifirst(i),i=1,remd_m(1))
 565              do il=1,nodes
 566               write (irest1,*) "nupa",il
 567               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 568               write (irest1,*) "ndowna",il
 569               write (irest1,*) ndowna(0,il),
 570      &                   (ndowna(i,il),i=1,ndowna(0,il))
 571              enddo
 572              if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
 573               write (irest1,*) "nset"
 574               write (irest1,*) nset
 575               write (irest1,*) "mset"
 576               write (irest1,*) (mset(i),i=1,nset)
 577               write (irest1,*) "i2set"
 578               write (irest1,*) (i2set(i),i=0,nodes-1)
 579               write (irest1,*) "i_index"
 580               do il=1,nset
 581                do il1=1,mset(il)
 582                 do i=1,nrep
 583                   write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
 584                 enddo
 585                enddo
 586               enddo
 587
 588              endif
 589              close(irest1)
 590            endif
 591            open(irest2,file=rest2name,status='unknown')
 592            write(irest2,*) totT,EK,potE,totE,t_bath
 593            do i=1,2*nres
 594             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 595            enddo
 596            do i=1,2*nres
 597             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
 598            enddo
 599            if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
 600              write (irest2,*) iset
 601            endif
 602           close(irest2)
 603           rstcount=0
 604         endif
 605
 606 c REMD - exchange
 607 c forced synchronization
 608         if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
 609      &                                .and. .not. mremdsync) then
 610             synflag=.false.
 611             call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
 612             if (synflag) then
 613                call mpi_recv(itime_master, 1, MPI_INTEGER,
 614      &                             0,101,CG_COMM, status, ierr)
 615                call mpi_barrier(CG_COMM, ierr)
 616 cdeb               if (out1file.or.traj1file) then
 617 cdeb                call mpi_gather(itime,1,mpi_integer,
 618 cdeb     &             icache_all,1,mpi_integer,king,
 619 cdeb     &             CG_COMM,ierr)
 620                if(traj1file)
 621      &          call mpi_gather(ntwx_cache,1,mpi_integer,
 622      &             icache_all,1,mpi_integer,king,
 623      &             CG_COMM,ierr)
 624                if (.not.out1file)
 625      &               write(iout,*) 'REMD synchro at',itime_master,itime
 626                if (itime_master.ge.n_timestep .or. ovrtim())
 627      &            end_of_run=.true.
 628 ctime               call flush(iout)
 629             endif
 630         endif
 631
 632 c REMD - exchange
 633         if ((mod(itime,nstex).eq.0.and.me.eq.king
 634      &                  .or.end_of_run.and.me.eq.king )
 635      &       .and. .not. mremdsync ) then
 636            synflag=.true.
 637            time01_=MPI_WTIME()
 638            do i=1,nodes-1
 639               call mpi_isend(itime,1,MPI_INTEGER,i,101,
 640      &                                CG_COMM, ireqi(i), ierr)
 641 cd            write(iout,*) 'REMD synchro with',i
 642 cd            call flush(iout)
 643            enddo
 644            call mpi_waitall(nodes-1,ireqi,statusi,ierr)
 645            call mpi_barrier(CG_COMM, ierr)
 646            time01=MPI_WTIME()
 647            write(iout,*) 'REMD synchro at',itime,'time=',time01-time01_
 648            if (out1file.or.traj1file) then
 649 cdeb            call mpi_gather(itime,1,mpi_integer,
 650 cdeb     &             itime_all,1,mpi_integer,king,
 651 cdeb     &             CG_COMM,ierr)
 652 cdeb            write(iout,'(a19,8000i8)') ' REMD synchro itime',
 653 cdeb     &                    (itime_all(i),i=1,nodes)
 654             if(traj1file) then
 655 cdeb             imin_itime=itime_all(1)
 656 cdeb             do i=2,nodes
 657 cdeb               if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
 658 cdeb             enddo
 659 cdeb             ii_write=(imin_itime-imin_itime_old)/ntwx
 660 cdeb             imin_itime_old=int(imin_itime/ntwx)*ntwx
 661 cdeb             write(iout,*) imin_itime,imin_itime_old,ii_write
 662              call mpi_gather(ntwx_cache,1,mpi_integer,
 663      &             icache_all,1,mpi_integer,king,
 664      &             CG_COMM,ierr)
 665 c             write(iout,'(a19,8000i8)') '     ntwx_cache',
 666 c     &                    (icache_all(i),i=1,nodes)
 667              ii_write=icache_all(1)
 668              do i=2,nodes
 669                if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
 670              enddo
 671 c             write(iout,*) "MIN ii_write=",ii_write
 672             endif
 673            endif
 674 ctime           call flush(iout)
 675         endif
 676         if(mremdsync .and. mod(itime,nstex).eq.0) then
 677            synflag=.true.
 678            if (me.eq.king .or. .not. out1file)
 679      &      write(iout,*) 'REMD synchro at',itime
 680
 681             if(traj1file) then
 682              call mpi_gather(ntwx_cache,1,mpi_integer,
 683      &             icache_all,1,mpi_integer,king,
 684      &             CG_COMM,ierr)
 685              if (me.eq.king) then
 686                write(iout,'(a19,8000i8)') '     ntwx_cache',
 687      &                    (icache_all(i),i=1,nodes)
 688                ii_write=icache_all(1)
 689                do i=2,nodes
 690                  if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
 691                enddo
 692                write(iout,*) "MIN ii_write=",ii_write
 693              endif
 694             endif
 695            call flush(iout)
 696         endif
 697         if (synflag) then
 698 c Update the time safety limiy
 699           if (time001-time00.gt.safety) then
 700             safety=time001-time00+600
 701             write (iout,*) "****** SAFETY increased to",safety," s"
 702           endif
 703           if (ovrtim() .and. me.eq.king) end_of_run=.true.
 704           call MPI_Bcast(end_of_run,1,MPI_LOGICAL,king,CG_COMM,IERR)
 705         endif
 706         if(synflag.and..not.end_of_run) then
 707            time02=MPI_WTIME()
 708            synflag=.false.
 709
 710 c           write(iout,*) 'REMD before',me,t_bath
 711
 712 c           call mpi_gather(t_bath,1,mpi_double_precision,
 713 c     &             remd_t_bath,1,mpi_double_precision,king,
 714 c     &             CG_COMM,ierr)
 715            potEcomp(n_ene+1)=t_bath
 716            t_bath_old=t_bath
 717            if (usampl.or.homol_nset.gt.1) then
 718              potEcomp(n_ene+2)=iset
 719              if (iset.lt.nset) then
 720                i_set_temp=iset
 721                iset=iset+1
 722                if (homol_nset.gt.1) then
 723 c broadcast iset to slaves and reduce energy
 724                 if (nfgtasks.gt.1) then
 725                  call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
 726                  call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
 727                  call e_modeller(e_tmp)
 728 c                 write(iout,*) "iset+1 before reduce",e_tmp
 729                  call MPI_Barrier(FG_COMM,IERR)
 730                  call MPI_Reduce(e_tmp,potEcomp(n_ene+3),1,
 731      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
 732                 else
 733                  call e_modeller(potEcomp(n_ene+3))
 734                 endif
 735 c                write(iout,*) "iset+1",potEcomp(n_ene+3)
 736                else
 737                 call EconstrQ
 738                 potEcomp(n_ene+3)=Uconst
 739                endif
 740                iset=i_set_temp
 741 c broadcast iset to slaves
 742                if (nfgtasks.gt.1) then
 743                  call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
 744                  call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
 745                endif
 746              else
 747               potEcomp(n_ene+3)=0.0
 748              endif
 749              if (iset.gt.1) then
 750                i_set_temp=iset
 751                iset=iset-1
 752                if (homol_nset.gt.1) then
 753 c broadcast iset to slaves and reduce energy
 754                 if (nfgtasks.gt.1) then
 755                  call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
 756                  call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
 757                  call e_modeller(e_tmp)
 758 c                 write(iout,*) "iset-1 before reduce",e_tmp
 759                  call MPI_Barrier(FG_COMM,IERR)
 760                  call MPI_Reduce(e_tmp,potEcomp(n_ene+4),1,
 761      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
 762                 else
 763                  call e_modeller(potEcomp(n_ene+4))
 764                 endif
 765 c                write(iout,*) "iset-1",potEcomp(n_ene+4)
 766                else
 767                 call EconstrQ
 768                 potEcomp(n_ene+4)=Uconst
 769                endif
 770                iset=i_set_temp
 771 c broadcast iset to slaves
 772                if (nfgtasks.gt.1) then
 773                  call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
 774                  call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
 775                endif
 776              else
 777                potEcomp(n_ene+4)=0.0
 778              endif
 779            endif
 780            call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
 781      &             remd_ene(0,1),n_ene+5,mpi_double_precision,king,
 782      &             CG_COMM,ierr)
 783            if(lmuca) then
 784             call mpi_gather(elow,1,mpi_double_precision,
 785      &             elowi,1,mpi_double_precision,king,
 786      &             CG_COMM,ierr)
 787             call mpi_gather(ehigh,1,mpi_double_precision,
 788      &             ehighi,1,mpi_double_precision,king,
 789      &             CG_COMM,ierr)
 790            endif
 791
 792           time03=MPI_WTIME()
 793           if (me.eq.king .or. .not. out1file) then
 794             write(iout,*) 'REMD gather times=',time03-time01
 795      &                                        ,time03-time02
 796           endif
 797
 798           if (restart1file) call write1rst(i_index)
 799
 800           time04=MPI_WTIME()
 801           if (me.eq.king .or. .not. out1file) then
 802             write(iout,*) 'REMD writing rst time=',time04-time03
 803           endif
 804
 805           if (traj1file) call write1traj
 806 cd debugging
 807 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
 808 cdeb     &             icache_all,1,mpi_integer,king,
 809 cdeb     &             CG_COMM,ierr)
 810 cdeb            write(iout,'(a19,8000i8)') '  ntwx_cache after traj1file',
 811 cdeb     &                    (icache_all(i),i=1,nodes)
 812 cd end
 813
 814
 815           time05=MPI_WTIME()
 816           if (me.eq.king .or. .not. out1file) then
 817             write(iout,*) 'REMD writing traj time=',time05-time04
 818             call flush(iout)
 819           endif
 820
 821
 822           if (me.eq.king) then
 823            if(homol_nset.gt.1) write(iout,*)
 824      &     'energy_c temperature iset energy_c(iset+1) energy_c(iset-1)'
 825             do i=1,nodes
 826                remd_t_bath(i)=remd_ene(n_ene+1,i)
 827                iremd_iset(i)=remd_ene(n_ene+2,i)
 828                if(homol_nset.gt.1)
 829      &                write(iout,'(i4,f10.3,f6.0,i3,2f10.3)')
 830      &                i,remd_ene(i_econstr,i),
 831      &                remd_ene(n_ene+1,i),iremd_iset(i),
 832      &                remd_ene(n_ene+3,i),remd_ene(n_ene+4,i)
 833             enddo
 834 #ifdef DEBUG
 835             if(lmuca) then
 836 co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
 837              do i=1,nodes
 838                write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
 839      &            elowi(i),ehighi(i)
 840              enddo
 841             else
 842               write(iout,*) 'REMD exchange temp,ene'
 843               do i=1,nodes
 844                 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
 845                 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
 846               enddo
 847             endif
 848 #endif
 849 c-------------------------------------
 850            IF(.not.usampl.and.hremd.eq.0.and.homol_nset.le.1) THEN
 851 #ifdef DEBUG
 852             write (iout,*) "Enter exchnge, remd_m",remd_m(1),
 853      &        " nodes",nodes
 854 ctime            call flush(iout)
 855             write (iout,*) "remd_m(1)",remd_m(1)
 856 #endif
 857             do irr=1,remd_m(1)
 858                i=ifirst(iran_num(1,remd_m(1)))
 859 #ifdef DEBUG
 860              write (iout,*) "i",i
 861 #endif
 862 ctime             call flush(iout)
 863
 864              do ii=1,nodes-1
 865
 866 #ifdef DEBUG
 867               write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
 868 #endif
 869              if(i.gt.0.and.nupa(0,i).gt.0) then
 870               iex=i
 871 c              if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
 872 c                write (iout,*)
 873 c     &  "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
 874 c                call flush(iout)
 875 c                call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
 876 c              endif
 877 c              do while (iex.eq.i)
 878 c                write (iout,*) "upper",nupa(int(nupa(0,i)),i)
 879                 iex=nupa(iran_num(1,int(nupa(0,i))),i)
 880 c              enddo
 881 c              write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
 882               if (lmuca) then
 883                call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
 884               else
 885 c Swap temperatures between conformations i and iex with recalculating the free energies
 886 c following temperature changes.
 887                ene_iex_iex=remd_ene(0,iex)
 888                ene_i_i=remd_ene(0,i)
 889 c               write (iout,*) "i",i," ene_i_i",ene_i_i,
 890 c     &          " iex",iex," ene_iex_iex",ene_iex_iex
 891 c               write (iout,*) "rescaling weights with temperature",
 892 c     &          remd_t_bath(i)
 893 c               call flush(iout)
 894                call rescale_weights(remd_t_bath(i))
 895
 896 c               write (iout,*) "0,iex",remd_t_bath(i)
 897 c               call enerprint(remd_ene(0,iex))
 898
 899                call sum_energy(remd_ene(0,iex),.false.)
 900                ene_iex_i=remd_ene(0,iex)
 901 c               write (iout,*) "ene_iex_i",remd_ene(0,iex)
 902
 903 c               write (iout,*) "0,i",remd_t_bath(i)
 904 c               call enerprint(remd_ene(0,i))
 905
 906                call sum_energy(remd_ene(0,i),.false.)
 907 c               write (iout,*) "ene_i_i",remd_ene(0,i)
 908 c               call flush(iout)
 909 c               write (iout,*) "rescaling weights with temperature",
 910 c     &          remd_t_bath(iex)
 911                if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
 912                 write (iout,*) "ERROR: inconsistent energies:",i,
 913      &            ene_i_i,remd_ene(0,i)
 914                endif
 915                call rescale_weights(remd_t_bath(iex))
 916
 917 c               write (iout,*) "0,i",remd_t_bath(iex)
 918 c               call enerprint(remd_ene(0,i))
 919
 920                call sum_energy(remd_ene(0,i),.false.)
 921 c               write (iout,*) "ene_i_iex",remd_ene(0,i)
 922 c               call flush(iout)
 923                ene_i_iex=remd_ene(0,i)
 924
 925 c               write (iout,*) "0,iex",remd_t_bath(iex)
 926 c               call enerprint(remd_ene(0,iex))
 927
 928                call sum_energy(remd_ene(0,iex),.false.)
 929                if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
 930                 write (iout,*) "ERROR: inconsistent energies:",iex,
 931      &            ene_iex_iex,remd_ene(0,iex)
 932                endif
 933 c               write (iout,*) "ene_iex_iex",remd_ene(0,iex)
 934 c               write (iout,*) "i",i," iex",iex
 935 c               write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
 936 c     &           " ene_i_iex",ene_i_iex,
 937 c     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
 938 c               call flush(iout)
 939                delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
 940      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
 941                delta=-delta
 942 c               write(iout,*) 'delta',delta
 943 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
 944 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
 945 c     &              (remd_t_bath(i)*remd_t_bath(iex))
 946               endif
 947               if (delta .gt. 50.0d0) then
 948                 delta=0.0d0
 949               else
 950 #ifdef OSF
 951                 if(isnan(delta))then
 952                   delta=0.0d0
 953                 else if (delta.lt.-50.0d0) then
 954                   delta=dexp(50.0d0)
 955                 else
 956                   delta=dexp(-delta)
 957                 endif
 958 #else
 959                 delta=dexp(-delta)
 960 #endif
 961               endif
 962               iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
 963               xxx=ran_number(0.0d0,1.0d0)
 964 c              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
 965 c              call flush(iout)
 966               if (delta .gt. xxx) then
 967                 tmp=remd_t_bath(i)
 968                 remd_t_bath(i)=remd_t_bath(iex)
 969                 remd_t_bath(iex)=tmp
 970                 remd_ene(0,i)=ene_i_iex
 971                 remd_ene(0,iex)=ene_iex_i
 972                 if(lmuca) then
 973                   tmp=elowi(i)
 974                   elowi(i)=elowi(iex)
 975                   elowi(iex)=tmp
 976                   tmp=ehighi(i)
 977                   ehighi(i)=ehighi(iex)
 978                   ehighi(iex)=tmp
 979                 endif
 980
 981
 982                 do k=0,nodes
 983                   itmp=nupa(k,i)
 984                   nupa(k,i)=nupa(k,iex)
 985                   nupa(k,iex)=itmp
 986                   itmp=ndowna(k,i)
 987                   ndowna(k,i)=ndowna(k,iex)
 988                   ndowna(k,iex)=itmp
 989                 enddo
 990                 do il=1,nodes
 991                  if (ifirst(il).eq.i) ifirst(il)=iex
 992                  do k=1,nupa(0,il)
 993                   if (nupa(k,il).eq.i) then
 994                      nupa(k,il)=iex
 995                   elseif (nupa(k,il).eq.iex) then
 996                      nupa(k,il)=i
 997                   endif
 998                  enddo
 999                  do k=1,ndowna(0,il)
1000                   if (ndowna(k,il).eq.i) then
1001                      ndowna(k,il)=iex
1002                   elseif (ndowna(k,il).eq.iex) then
1003                      ndowna(k,il)=i
1004                   endif
1005                  enddo
1006                 enddo
1007
1008                 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1009                 itmp=i2rep(i-1)
1010                 i2rep(i-1)=i2rep(iex-1)
1011                 i2rep(iex-1)=itmp
1012
1013 c                write(iout,*) 'exchange',i,iex
1014 c                write (iout,'(a8,100i4)') "@ ifirst",
1015 c     &                    (ifirst(k),k=1,remd_m(1))
1016 c                do il=1,nodes
1017 c                 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1018 c     &                    (nupa(k,il),k=1,nupa(0,il))
1019 c                 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1020 c     &                    (ndowna(k,il),k=1,ndowna(0,il))
1021 c                enddo
1022 c                call flush(iout)
1023
1024               else
1025                remd_ene(0,iex)=ene_iex_iex
1026                remd_ene(0,i)=ene_i_i
1027                i=iex
1028               endif
1029             endif
1030            enddo
1031            enddo
1032 cd           write (iout,*) "exchange completed"
1033 cd           call flush(iout)
1034         ELSEIF (usampl.or.homol_nset.gt.1) THEN
1035           do ii=1,nodes
1036 c            write(iout,*) "########",ii
1037
1038             i_temp=iran_num(1,nrep)
1039             i_mult=iran_num(1,remd_m(i_temp))
1040             i_iset=iran_num(1,nset)
1041             i_mset=iran_num(1,mset(i_iset))
1042             i=i_index(i_temp,i_mult,i_iset,i_mset)
1043
1044 c            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1045
1046             if(t_exchange_only)then
1047              i_dir=1
1048             else
1049              i_dir=iran_num(1,3)
1050             endif
1051 c            write(iout,*) "i_dir=",i_dir
1052
1053             if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1054
1055                i_temp1=i_temp+1
1056                i_mult1=iran_num(1,remd_m(i_temp1))
1057                i_iset1=i_iset
1058                i_mset1=iran_num(1,mset(i_iset1))
1059                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1060
1061             elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1062
1063                i_temp1=i_temp
1064                i_mult1=iran_num(1,remd_m(i_temp1))
1065                i_iset1=i_iset+1
1066                i_mset1=iran_num(1,mset(i_iset1))
1067                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1068
1069                econstr_temp_i=remd_ene(i_econstr,i)
1070                econstr_temp_iex=remd_ene(i_econstr,iex)
1071                remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1072                remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1073
1074             elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1075
1076                i_temp1=i_temp+1
1077                i_mult1=iran_num(1,remd_m(i_temp1))
1078                i_iset1=i_iset+1
1079                i_mset1=iran_num(1,mset(i_iset1))
1080                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1081                econstr_temp_i=remd_ene(i_econstr,i)
1082                econstr_temp_iex=remd_ene(i_econstr,iex)
1083                remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1084                remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1085
1086             else
1087                goto 444
1088             endif
1089
1090 c            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1091 ctime            call flush(iout)
1092
1093 c Swap temperatures between conformations i and iex with recalculating the free energies
1094 c following temperature changes.
1095               ene_iex_iex=remd_ene(0,iex)
1096               ene_i_i=remd_ene(0,i)
1097 co              write (iout,*) "rescaling weights with temperature",
1098 co     &          remd_t_bath(i)
1099               call rescale_weights(remd_t_bath(i))
1100
1101               call sum_energy(remd_ene(0,iex),.false.)
1102               ene_iex_i=remd_ene(0,iex)
1103 cdebug
1104 c ERROR only makes sense for dir =1
1105 c              write (iout,*) "ene_iex_i",remd_ene(0,iex)
1106 c              call sum_energy(remd_ene(0,i),.false.)
1107 c              write (iout,*) "ene_i_i",remd_ene(0,i)
1108 c              write (iout,*) "rescaling weights with temperature",
1109 c     &          remd_t_bath(iex)
1110 c              if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1111 c                write (iout,*) "ERROR: inconsistent energies i:",i,
1112 c     &            ene_i_i,remd_ene(0,i)
1113 c              endif
1114 cdebug_end
1115               call rescale_weights(remd_t_bath(iex))
1116               call sum_energy(remd_ene(0,i),.false.)
1117 cd              write (iout,*) "ene_i_iex",remd_ene(0,i)
1118               ene_i_iex=remd_ene(0,i)
1119 cdebug
1120 c ERROR only makes sense for dir =1
1121 c              call sum_energy(remd_ene(0,iex),.false.)
1122 c              if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1123 c                write (iout,*) "ERROR: inconsistent energies iex:",iex,
1124 c     &            ene_iex_iex,remd_ene(0,iex)
1125 c              endif
1126 c              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1127 c              write (iout,*) "i",i," iex",iex
1128 c              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1129 c     &           " ene_i_iex",ene_i_iex,
1130 c     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1131 cdebug_end
1132               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1133      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1134               delta=-delta
1135 c              write(iout,*) 'delta',delta
1136 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
1137 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
1138 c     &              (remd_t_bath(i)*remd_t_bath(iex))
1139               if (delta .gt. 50.0d0) then
1140                 delta=0.0d0
1141               else
1142                 delta=dexp(-delta)
1143               endif
1144               if (i_dir.eq.1.or.i_dir.eq.3)
1145      &         iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1146               if (i_dir.eq.2.or.i_dir.eq.3)
1147      &          iremd_tot_usa(int(i2set(i-1)))=
1148      &                 iremd_tot_usa(int(i2set(i-1)))+1
1149               xxx=ran_number(0.0d0,1.0d0)
1150 c              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1151               if (delta .gt. xxx) then
1152                 tmp=remd_t_bath(i)
1153                 remd_t_bath(i)=remd_t_bath(iex)
1154                 remd_t_bath(iex)=tmp
1155
1156                 itmp=iremd_iset(i)
1157                 iremd_iset(i)=iremd_iset(iex)
1158                 iremd_iset(iex)=itmp
1159
1160                 remd_ene(0,i)=ene_i_iex
1161                 remd_ene(0,iex)=ene_iex_i
1162
1163                 if (i_dir.eq.1.or.i_dir.eq.3)
1164      &           iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1165
1166                 itmp=i2rep(i-1)
1167                 i2rep(i-1)=i2rep(iex-1)
1168                 i2rep(iex-1)=itmp
1169
1170                 if (i_dir.eq.2.or.i_dir.eq.3)
1171      &           iremd_acc_usa(int(i2set(i-1)))=
1172      &                 iremd_acc_usa(int(i2set(i-1)))+1
1173
1174                 itmp=i2set(i-1)
1175                 i2set(i-1)=i2set(iex-1)
1176                 i2set(iex-1)=itmp
1177
1178                 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1179                 i_index(i_temp,i_mult,i_iset,i_mset)=
1180      &                i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1181                 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1182
1183               else
1184                remd_ene(0,iex)=ene_iex_iex
1185                remd_ene(0,i)=ene_i_i
1186                remd_ene(i_econstr,iex)=econstr_temp_iex
1187                remd_ene(i_econstr,i)=econstr_temp_i
1188               endif
1189
1190 cd      do il=1,nset
1191 cd       do il1=1,mset(il)
1192 cd        do i=1,nrep
1193 cd         do j=1,remd_m(i)
1194 cd          write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1195 cd         enddo
1196 cd        enddo
1197 cd       enddo
1198 cd      enddo
1199
1200  444      continue
1201
1202           enddo
1203
1204
1205         ENDIF
1206
1207 c-------------------------------------
1208              write (iout,*) "NREP",nrep
1209              do i=1,nrep
1210               if(iremd_tot(i).ne.0)
1211      &          write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1212      &           ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1213              enddo
1214
1215              if(usampl.or.homol_nset.gt.1) then
1216               do i=1,nset
1217                if(iremd_tot_usa(i).ne.0)
1218      &           write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1219      &         iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1220               enddo
1221              endif
1222
1223              call flush(iout)
1224
1225 cd              write (iout,'(a6,100i4)') "ifirst",
1226 cd     &                    (ifirst(i),i=1,remd_m(1))
1227 cd              do il=1,nodes
1228 cd               write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1229 cd     &                    (nupa(i,il),i=1,nupa(0,il))
1230 cd               write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1231 cd     &                    (ndowna(i,il),i=1,ndowna(0,il))
1232 cd              enddo
1233             endif
1234
1235          time06=MPI_WTIME()
1236 cd         write (iout,*) "Before scatter"
1237 cd         call flush(iout)
1238          call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1239      &           t_bath,1,mpi_double_precision,king,
1240      &           CG_COMM,ierr)
1241 cd         write (iout,*) "After scatter"
1242 cd         call flush(iout)
1243          if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
1244           call mpi_scatter(iremd_iset,1,mpi_integer,
1245      &           iset,1,mpi_integer,king,
1246      &           CG_COMM,ierr)
1247 c 8/31/2015 Correction by AL: send new iset to slaves
1248           if (nfgtasks.gt.1) then
1249            call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
1250            call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
1251           endif
1252
1253          endif
1254
1255          time07=MPI_WTIME()
1256           if (me.eq.king .or. .not. out1file) then
1257             write(iout,*) 'REMD scatter time=',time07-time06
1258           endif
1259
1260          if(lmuca) then
1261            call mpi_scatter(elowi,1,mpi_double_precision,
1262      &           elow,1,mpi_double_precision,king,
1263      &           CG_COMM,ierr)
1264            call mpi_scatter(ehighi,1,mpi_double_precision,
1265      &           ehigh,1,mpi_double_precision,king,
1266      &           CG_COMM,ierr)
1267          endif
1268          call rescale_weights(t_bath)
1269 co         write (iout,*) "Processor",me,
1270 co     &    " rescaling weights with temperature",t_bath
1271
1272          stdfp=dsqrt(2*Rb*t_bath/d_time)
1273          do i=1,ntyp
1274            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1275          enddo
1276
1277 c Compute the standard deviations of stochastic forces for Langevin dynamics
1278 c if the friction coefficients do not depend on surface area
1279       if (lang.gt.0 .and. .not.surfarea) then
1280         do i=nnt,nct-1
1281           stdforcp(i)=stdfp*dsqrt(gamp)
1282         enddo
1283         do i=nnt,nct
1284           if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1285      &                *dsqrt(gamsc(iabs(itype(i))))
1286         enddo
1287       endif
1288
1289 cde         write(iout,*) 'REMD after',me,t_bath
1290            time08=MPI_WTIME()
1291            if (me.eq.king .or. .not. out1file) then
1292             write(iout,*) 'REMD exchange time=',time08-time02
1293 ctime            call flush(iout)
1294            endif
1295         endif
1296       enddo
1297
1298       if (restart1file) then
1299           if (me.eq.king .or. .not. out1file)
1300      &      write(iout,*) 'writing restart at the end of run'
1301            call write1rst(i_index)
1302       endif
1303
1304       if (traj1file) call write1traj
1305 cd debugging
1306 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
1307 cdeb     &             icache_all,1,mpi_integer,king,
1308 cdeb     &             CG_COMM,ierr)
1309 cdeb            write(iout,'(a40,8000i8)')
1310 cdeb     &             '  ntwx_cache after traj1file at the end',
1311 cdeb     &             (icache_all(i),i=1,nodes)
1312 cd end
1313
1314
1315 #ifdef MPI
1316       t_MD=MPI_Wtime()-tt0
1317 #else
1318       t_MD=tcpu()-tt0
1319 #endif
1320       if (me.eq.king .or. .not. out1file) then
1321        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1322      &  '  Timing  ',
1323      & 'MD calculations setup:',t_MDsetup,
1324      & 'Energy & gradient evaluation:',t_enegrad,
1325      & 'Stochastic MD setup:',t_langsetup,
1326      & 'Stochastic MD step setup:',t_sdsetup,
1327      & 'MD steps:',t_MD
1328        write (iout,'(/28(1h=),a25,27(1h=))')
1329      & '  End of MD calculation  '
1330       endif
1331       return
1332       end
1333
1334 c-----------------------------------------------------------------------
1335       subroutine write1rst(i_index)
1336       implicit real*8 (a-h,o-z)
1337       include 'DIMENSIONS'
1338       include 'mpif.h'
1339       include 'COMMON.MD'
1340       include 'COMMON.IOUNITS'
1341       include 'COMMON.REMD'
1342       include 'COMMON.SETUP'
1343       include 'COMMON.CHAIN'
1344       include 'COMMON.SBRIDGE'
1345       include 'COMMON.INTERACT'
1346       include 'COMMON.CONTROL'
1347
1348       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1349      &     d_restart2(3,2*maxres*maxprocs)
1350       real t5_restart1(5)
1351       integer iret,itmp
1352       integer*2 i_index
1353      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1354        common /przechowalnia/ d_restart1,d_restart2
1355
1356        t5_restart1(1)=totT
1357        t5_restart1(2)=EK
1358        t5_restart1(3)=potE
1359        t5_restart1(4)=t_bath
1360        t5_restart1(5)=Uconst
1361
1362        call mpi_gather(t5_restart1,5,mpi_real,
1363      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1364
1365
1366        do i=1,2*nres
1367          do j=1,3
1368            r_d(j,i)=d_t(j,i)
1369          enddo
1370        enddo
1371        call mpi_gather(r_d,3*2*nres,mpi_real,
1372      &           d_restart1,3*2*nres,mpi_real,king,
1373      &           CG_COMM,ierr)
1374
1375
1376        do i=1,2*nres
1377          do j=1,3
1378            r_d(j,i)=dc(j,i)
1379          enddo
1380        enddo
1381        call mpi_gather(r_d,3*2*nres,mpi_real,
1382      &           d_restart2,3*2*nres,mpi_real,king,
1383      &           CG_COMM,ierr)
1384
1385        if(me.eq.king) then
1386 #ifdef AIX
1387          call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1388          do i=0,nodes-1
1389           call xdrfint_(ixdrf, i2rep(i), iret)
1390          enddo
1391          do i=1,remd_m(1)
1392           call xdrfint_(ixdrf, ifirst(i), iret)
1393          enddo
1394          do il=1,nodes
1395               do i=0,nupa(0,il)
1396                call xdrfint_(ixdrf, nupa(i,il), iret)
1397               enddo
1398
1399               do i=0,ndowna(0,il)
1400                call xdrfint_(ixdrf, ndowna(i,il), iret)
1401               enddo
1402          enddo
1403
1404          do il=1,nodes
1405            do j=1,4
1406             call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1407            enddo
1408          enddo
1409
1410          do il=0,nodes-1
1411            do i=1,2*nres
1412             do j=1,3
1413              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1414             enddo
1415            enddo
1416          enddo
1417          do il=0,nodes-1
1418            do i=1,2*nres
1419             do j=1,3
1420              call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1421             enddo
1422            enddo
1423          enddo
1424
1425          if(usampl.or.homol_nset.gt.1) then
1426            call xdrfint_(ixdrf, nset, iret)
1427            do i=1,nset
1428              call xdrfint_(ixdrf,mset(i), iret)
1429            enddo
1430            do i=0,nodes-1
1431              call xdrfint_(ixdrf,i2set(i), iret)
1432            enddo
1433            do il=1,nset
1434              do il1=1,mset(il)
1435                do i=1,nrep
1436                  do j=1,remd_m(i)
1437                    itmp=i_index(i,j,il,il1)
1438                    call xdrfint_(ixdrf,itmp, iret)
1439                  enddo
1440                enddo
1441              enddo
1442            enddo
1443
1444          endif
1445          call xdrfclose_(ixdrf, iret)
1446 #else
1447          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1448          do i=0,nodes-1
1449           call xdrfint(ixdrf, i2rep(i), iret)
1450          enddo
1451          do i=1,remd_m(1)
1452           call xdrfint(ixdrf, ifirst(i), iret)
1453          enddo
1454          do il=1,nodes
1455               do i=0,nupa(0,il)
1456                call xdrfint(ixdrf, nupa(i,il), iret)
1457               enddo
1458
1459               do i=0,ndowna(0,il)
1460                call xdrfint(ixdrf, ndowna(i,il), iret)
1461               enddo
1462          enddo
1463
1464          do il=1,nodes
1465            do j=1,4
1466             call xdrffloat(ixdrf, t_restart1(j,il), iret)
1467            enddo
1468          enddo
1469
1470          do il=0,nodes-1
1471            do i=1,2*nres
1472             do j=1,3
1473              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1474             enddo
1475            enddo
1476          enddo
1477          do il=0,nodes-1
1478            do i=1,2*nres
1479             do j=1,3
1480              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1481             enddo
1482            enddo
1483          enddo
1484
1485
1486              if(usampl.or.homol_nset.gt.1) then
1487               call xdrfint(ixdrf, nset, iret)
1488               do i=1,nset
1489                 call xdrfint(ixdrf,mset(i), iret)
1490               enddo
1491               do i=0,nodes-1
1492                 call xdrfint(ixdrf,i2set(i), iret)
1493               enddo
1494               do il=1,nset
1495                do il1=1,mset(il)
1496                 do i=1,nrep
1497                  do j=1,remd_m(i)
1498                    itmp=i_index(i,j,il,il1)
1499                    call xdrfint(ixdrf,itmp, iret)
1500                  enddo
1501                 enddo
1502                enddo
1503               enddo
1504
1505              endif
1506          call xdrfclose(ixdrf, iret)
1507 #endif
1508        endif
1509       return
1510       end
1511
1512
1513       subroutine write1traj
1514       implicit real*8 (a-h,o-z)
1515       include 'DIMENSIONS'
1516       include 'mpif.h'
1517       include 'COMMON.MD'
1518       include 'COMMON.IOUNITS'
1519       include 'COMMON.REMD'
1520       include 'COMMON.SETUP'
1521       include 'COMMON.CHAIN'
1522       include 'COMMON.SBRIDGE'
1523       include 'COMMON.INTERACT'
1524
1525       real t5_restart1(5)
1526       integer iret,itmp
1527       real xcoord(3,maxres2+2),prec
1528       real r_qfrag(50),r_qpair(100)
1529       real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1530       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1531       real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1532      &     p_uscdiff(100*maxprocs)
1533       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1534       common /przechowalnia/ p_c
1535
1536       call mpi_bcast(ii_write,1,mpi_integer,
1537      &           king,CG_COMM,ierr)
1538
1539 c debugging
1540       print *,'traj1file',me,ii_write,ntwx_cache
1541 c end debugging
1542
1543 #ifdef AIX
1544       if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1545 #else
1546       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1547 #endif
1548       do ii=1,ii_write
1549        t5_restart1(1)=totT_cache(ii)
1550        t5_restart1(2)=EK_cache(ii)
1551        t5_restart1(3)=potE_cache(ii)
1552        t5_restart1(4)=t_bath_cache(ii)
1553        t5_restart1(5)=Uconst_cache(ii)
1554        call mpi_gather(t5_restart1,5,mpi_real,
1555      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1556
1557        call mpi_gather(iset_cache(ii),1,mpi_integer,
1558      &      iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1559
1560           do i=1,nfrag
1561            r_qfrag(i)=qfrag_cache(i,ii)
1562           enddo
1563           do i=1,npair
1564            r_qpair(i)=qpair_cache(i,ii)
1565           enddo
1566           do i=1,nfrag_back
1567            r_utheta(i)=utheta_cache(i,ii)
1568            r_ugamma(i)=ugamma_cache(i,ii)
1569            r_uscdiff(i)=uscdiff_cache(i,ii)
1570           enddo
1571
1572         call mpi_gather(r_qfrag,nfrag,mpi_real,
1573      &           p_qfrag,nfrag,mpi_real,king,
1574      &           CG_COMM,ierr)
1575         call mpi_gather(r_qpair,npair,mpi_real,
1576      &           p_qpair,npair,mpi_real,king,
1577      &           CG_COMM,ierr)
1578         call mpi_gather(r_utheta,nfrag_back,mpi_real,
1579      &           p_utheta,nfrag_back,mpi_real,king,
1580      &           CG_COMM,ierr)
1581         call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1582      &           p_ugamma,nfrag_back,mpi_real,king,
1583      &           CG_COMM,ierr)
1584         call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1585      &           p_uscdiff,nfrag_back,mpi_real,king,
1586      &           CG_COMM,ierr)
1587
1588 #ifdef DEBUG
1589         write (iout,*) "p_qfrag"
1590         do i=1,nodes
1591           write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1592         enddo
1593         write (iout,*) "p_qpair"
1594         do i=1,nodes
1595           write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1596         enddo
1597         call flush(iout)
1598 #endif
1599         do i=1,nres*2
1600          do j=1,3
1601           r_c(j,i)=c_cache(j,i,ii)
1602          enddo
1603         enddo
1604
1605         call mpi_gather(r_c,3*2*nres,mpi_real,
1606      &           p_c,3*2*nres,mpi_real,king,
1607      &           CG_COMM,ierr)
1608
1609        if(me.eq.king) then
1610 #ifdef AIX
1611          do il=1,nodes
1612           call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1613           call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1614           call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1615           call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1616           call xdrfint_(ixdrf, nss, iret)
1617           do j=1,nss
1618            if (dyn_ss) then
1619             call xdrfint_(ixdrf, idssb(j)+nres, iret)
1620             call xdrfint_(ixdrf, jdssb(j)+nres, iret)
1621            else
1622             call xdrfint_(ixdrf, ihpb(j), iret)
1623             call xdrfint_(ixdrf, jhpb(j), iret)
1624            endif
1625           enddo
1626           call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1627           call xdrfint_(ixdrf, iset_restart1(il), iret)
1628           do i=1,nfrag
1629            call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1630           enddo
1631           do i=1,npair
1632            call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1633           enddo
1634           do i=1,nfrag_back
1635            call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1636            call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1637            call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1638           enddo
1639           prec=10000.0
1640           do i=1,nres
1641            do j=1,3
1642             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1643            enddo
1644           enddo
1645           do i=nnt,nct
1646            do j=1,3
1647             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1648            enddo
1649           enddo
1650           itmp=nres+nct-nnt+1
1651           call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1652          enddo
1653 #else
1654          do il=1,nodes
1655           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1656           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1657           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1658           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1659           call xdrfint(ixdrf, nss, iret)
1660           do j=1,nss
1661            if (dyn_ss) then
1662             call xdrfint(ixdrf, idssb(j)+nres, iret)
1663             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1664            else
1665             call xdrfint(ixdrf, ihpb(j), iret)
1666             call xdrfint(ixdrf, jhpb(j), iret)
1667            endif
1668           enddo
1669           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1670           call xdrfint(ixdrf, iset_restart1(il), iret)
1671           do i=1,nfrag
1672            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1673           enddo
1674           do i=1,npair
1675            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1676           enddo
1677           do i=1,nfrag_back
1678            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1679            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1680            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1681           enddo
1682           prec=10000.0
1683           do i=1,nres
1684            do j=1,3
1685             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1686            enddo
1687           enddo
1688           do i=nnt,nct
1689            do j=1,3
1690             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1691            enddo
1692           enddo
1693           itmp=nres+nct-nnt+1
1694           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1695          enddo
1696 #endif
1697        endif
1698       enddo
1699 #ifdef AIX
1700       if(me.eq.king) call xdrfclose_(ixdrf, iret)
1701 #else
1702       if(me.eq.king) call xdrfclose(ixdrf, iret)
1703 #endif
1704       do i=1,ntwx_cache-ii_write
1705
1706             totT_cache(i)=totT_cache(ii_write+i)
1707             EK_cache(i)=EK_cache(ii_write+i)
1708             potE_cache(i)=potE_cache(ii_write+i)
1709             t_bath_cache(i)=t_bath_cache(ii_write+i)
1710             Uconst_cache(i)=Uconst_cache(ii_write+i)
1711             iset_cache(i)=iset_cache(ii_write+i)
1712
1713             do ii=1,nfrag
1714              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1715             enddo
1716             do ii=1,npair
1717              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1718             enddo
1719             do ii=1,nfrag_back
1720               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1721               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1722               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1723             enddo
1724
1725             do ii=1,nres*2
1726              do j=1,3
1727               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1728              enddo
1729             enddo
1730       enddo
1731       ntwx_cache=ntwx_cache-ii_write
1732       return
1733       end
1734
1735
1736       subroutine read1restart(i_index)
1737       implicit real*8 (a-h,o-z)
1738       include 'DIMENSIONS'
1739       include 'mpif.h'
1740       include 'COMMON.MD'
1741       include 'COMMON.IOUNITS'
1742       include 'COMMON.REMD'
1743       include 'COMMON.SETUP'
1744       include 'COMMON.CHAIN'
1745       include 'COMMON.SBRIDGE'
1746       include 'COMMON.INTERACT'
1747       include 'COMMON.CONTROL'
1748       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1749      &                 t5_restart1(5)
1750       integer*2 i_index
1751      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1752       common /przechowalnia/ d_restart1
1753       write (*,*) "Processor",me," called read1restart"
1754
1755          if(me.eq.king)then
1756               open(irest2,file=mremd_rst_name,status='unknown')
1757               read(irest2,*,err=334) i
1758               write(iout,*) "Reading old rst in ASCI format"
1759               close(irest2)
1760                call read1restart_old
1761                return
1762  334          continue
1763 #ifdef AIX
1764               call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1765
1766               do i=0,nodes-1
1767                call xdrfint_(ixdrf, i2rep(i), iret)
1768               enddo
1769               do i=1,remd_m(1)
1770                call xdrfint_(ixdrf, ifirst(i), iret)
1771               enddo
1772              do il=1,nodes
1773               call xdrfint_(ixdrf, nupa(0,il), iret)
1774               do i=1,nupa(0,il)
1775                call xdrfint_(ixdrf, nupa(i,il), iret)
1776               enddo
1777
1778               call xdrfint_(ixdrf, ndowna(0,il), iret)
1779               do i=1,ndowna(0,il)
1780                call xdrfint_(ixdrf, ndowna(i,il), iret)
1781               enddo
1782              enddo
1783              do il=1,nodes
1784                do j=1,4
1785                 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1786                enddo
1787              enddo
1788 #else
1789               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1790
1791               do i=0,nodes-1
1792                call xdrfint(ixdrf, i2rep(i), iret)
1793               enddo
1794               do i=1,remd_m(1)
1795                call xdrfint(ixdrf, ifirst(i), iret)
1796               enddo
1797              do il=1,nodes
1798               call xdrfint(ixdrf, nupa(0,il), iret)
1799               do i=1,nupa(0,il)
1800                call xdrfint(ixdrf, nupa(i,il), iret)
1801               enddo
1802
1803               call xdrfint(ixdrf, ndowna(0,il), iret)
1804               do i=1,ndowna(0,il)
1805                call xdrfint(ixdrf, ndowna(i,il), iret)
1806               enddo
1807              enddo
1808              do il=1,nodes
1809                do j=1,4
1810                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1811                enddo
1812              enddo
1813 #endif
1814          endif
1815          call mpi_scatter(t_restart1,5,mpi_real,
1816      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1817          totT=t5_restart1(1)
1818          EK=t5_restart1(2)
1819          potE=t5_restart1(3)
1820          t_bath=t5_restart1(4)
1821
1822          if(me.eq.king)then
1823               do il=0,nodes-1
1824                do i=1,2*nres
1825 c                read(irest2,'(3e15.5)')
1826 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1827             do j=1,3
1828 #ifdef AIX
1829              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1830 #else
1831              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1832 #endif
1833             enddo
1834                enddo
1835 #ifdef DEBUG
1836             write (iout,*) "Conformation read",il
1837             do i=1,nres
1838               write (iout,'(i5,3f10.5,5x,3f10.5)')
1839      &          i,(d_restart1(j,i+2*nres*il),j=1,3),
1840      &            (d_restart1(j,nres+i+2*nres*il),j=1,3)
1841             enddo
1842 #endif
1843               enddo
1844          endif
1845          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1846      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1847
1848          do i=1,2*nres
1849            do j=1,3
1850             d_t(j,i)=r_d(j,i)
1851            enddo
1852          enddo
1853          if(me.eq.king)then
1854               do il=0,nodes-1
1855                do i=1,2*nres
1856 c                read(irest2,'(3e15.5)')
1857 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1858             do j=1,3
1859 #ifdef AIX
1860              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1861 #else
1862              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1863 #endif
1864             enddo
1865                enddo
1866               enddo
1867          endif
1868          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1869      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1870          do i=1,2*nres
1871            do j=1,3
1872             dc(j,i)=r_d(j,i)
1873            enddo
1874          enddo
1875
1876
1877            if(usampl.or.homol_nset.gt.1) then
1878 #ifdef AIX
1879              if(me.eq.king)then
1880               call xdrfint_(ixdrf, nset, iret)
1881               do i=1,nset
1882                 call xdrfint_(ixdrf,mset(i), iret)
1883               enddo
1884               do i=0,nodes-1
1885                 call xdrfint_(ixdrf,i2set(i), iret)
1886               enddo
1887               do il=1,nset
1888                do il1=1,mset(il)
1889                 do i=1,nrep
1890                  do j=1,remd_m(i)
1891                    call xdrfint_(ixdrf,itmp, iret)
1892                    i_index(i,j,il,il1)=itmp
1893                  enddo
1894                 enddo
1895                enddo
1896               enddo
1897              endif
1898 #else
1899              if(me.eq.king)then
1900               call xdrfint(ixdrf, nset, iret)
1901               do i=1,nset
1902                 call xdrfint(ixdrf,mset(i), iret)
1903               enddo
1904               do i=0,nodes-1
1905                 call xdrfint(ixdrf,i2set(i), iret)
1906               enddo
1907               do il=1,nset
1908                do il1=1,mset(il)
1909                 do i=1,nrep
1910                  do j=1,remd_m(i)
1911                    call xdrfint(ixdrf,itmp, iret)
1912                    i_index(i,j,il,il1)=itmp
1913                  enddo
1914                 enddo
1915                enddo
1916               enddo
1917              endif
1918 #endif
1919 c Corrected AL 8/19/2014: each processor needs whole iset array not only its
1920 c own element
1921 c              call mpi_scatter(i2set,1,mpi_integer,
1922 c     &           iset,1,mpi_integer,king,
1923 c     &           CG_COMM,ierr)
1924               call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1925      &         CG_COMM,ierr)
1926               iset=i2set(me)
1927 c broadcast iset to slaves
1928               if (nfgtasks.gt.1) then
1929                call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
1930                call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
1931               endif
1932            endif
1933
1934
1935         if(me.eq.king) close(irest2)
1936         return
1937         end
1938
1939       subroutine read1restart_old
1940       implicit real*8 (a-h,o-z)
1941       include 'DIMENSIONS'
1942       include 'mpif.h'
1943       include 'COMMON.MD'
1944       include 'COMMON.IOUNITS'
1945       include 'COMMON.REMD'
1946       include 'COMMON.SETUP'
1947       include 'COMMON.CHAIN'
1948       include 'COMMON.SBRIDGE'
1949       include 'COMMON.INTERACT'
1950       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1951      &                 t5_restart1(5)
1952       common /przechowalnia/ d_restart1
1953          if(me.eq.king)then
1954              open(irest2,file=mremd_rst_name,status='unknown')
1955              read (irest2,*) (i2rep(i),i=0,nodes-1)
1956              read (irest2,*) (ifirst(i),i=1,remd_m(1))
1957              do il=1,nodes
1958               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1959               read (irest2,*) ndowna(0,il),
1960      &                    (ndowna(i,il),i=1,ndowna(0,il))
1961              enddo
1962              do il=1,nodes
1963                read(irest2,*) (t_restart1(j,il),j=1,4)
1964              enddo
1965          endif
1966          call mpi_scatter(t_restart1,5,mpi_real,
1967      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1968          totT=t5_restart1(1)
1969          EK=t5_restart1(2)
1970          potE=t5_restart1(3)
1971          t_bath=t5_restart1(4)
1972
1973          if(me.eq.king)then
1974               do il=0,nodes-1
1975                do i=1,2*nres
1976                 read(irest2,'(3e15.5)')
1977      &                (d_restart1(j,i+2*nres*il),j=1,3)
1978                enddo
1979               enddo
1980          endif
1981          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1982      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1983
1984          do i=1,2*nres
1985            do j=1,3
1986             d_t(j,i)=r_d(j,i)
1987            enddo
1988          enddo
1989          if(me.eq.king)then
1990               do il=0,nodes-1
1991                do i=1,2*nres
1992                 read(irest2,'(3e15.5)')
1993      &                (d_restart1(j,i+2*nres*il),j=1,3)
1994                enddo
1995               enddo
1996          endif
1997          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1998      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1999          do i=1,2*nres
2000            do j=1,3
2001             dc(j,i)=r_d(j,i)
2002            enddo
2003          enddo
2004         if(me.eq.king) close(irest2)
2005         return
2006         end
2007 c------------------------------------------