source/unres/src-HCD-5D/MREMD.F

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