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