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 cde         write(iout,*) 'REMD after',me,t_bath
1283            time08=MPI_WTIME()
1284            if (me.eq.king .or. .not. out1file) then
1285             write(iout,*) 'REMD exchange time=',time08-time00
1286             call flush(iout)
1287            endif
1288         endif
1289       enddo
1290
1291       if (restart1file) then
1292           if (me.eq.king .or. .not. out1file)
1293      &      write(iout,*) 'writing restart at the end of run'
1294            call write1rst(i_index)
1295       endif
1296
1297       if (traj1file) call write1traj
1298 cd debugging
1299 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
1300 cdeb     &             icache_all,1,mpi_integer,king,
1301 cdeb     &             CG_COMM,ierr)
1302 cdeb            write(iout,'(a40,8000i8)')
1303 cdeb     &             '  ntwx_cache after traj1file at the end',
1304 cdeb     &             (icache_all(i),i=1,nodes)
1305 cd end
1306
1307
1308 #ifdef MPI
1309       t_MD=MPI_Wtime()-tt0
1310 #else
1311       t_MD=tcpu()-tt0
1312 #endif
1313       if (me.eq.king .or. .not. out1file) then
1314        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1315      &  '  Timing  ',
1316      & 'MD calculations setup:',t_MDsetup,
1317      & 'Energy & gradient evaluation:',t_enegrad,
1318      & 'Stochastic MD setup:',t_langsetup,
1319      & 'Stochastic MD step setup:',t_sdsetup,
1320      & 'MD steps:',t_MD
1321        write (iout,'(/28(1h=),a25,27(1h=))')
1322      & '  End of MD calculation  '
1323       endif
1324       return
1325       end
1326
1327 c-----------------------------------------------------------------------
1328       subroutine write1rst(i_index)
1329       implicit none
1330       include 'DIMENSIONS'
1331       include 'mpif.h'
1332       include 'COMMON.CONTROL'
1333       include 'COMMON.MD'
1334 #ifdef FIVEDIAG
1335        include 'COMMON.LAGRANGE.5diag'
1336 #else
1337        include 'COMMON.LAGRANGE'
1338 #endif
1339       include 'COMMON.QRESTR'
1340       include 'COMMON.IOUNITS'
1341       include 'COMMON.REMD'
1342       include 'COMMON.SETUP'
1343       include 'COMMON.CHAIN'
1344       include 'COMMON.SBRIDGE'
1345       include 'COMMON.INTERACT'
1346
1347       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1348      &     d_restart2(3,2*maxres*maxprocs)
1349       real t5_restart1(5)
1350       integer iret,itmp
1351       integer*2 i_index
1352      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1353        common /przechowalnia/ d_restart1,d_restart2
1354       integer i,j,il1,il,ixdrf
1355       integer ierr
1356
1357        t5_restart1(1)=totT
1358        t5_restart1(2)=EK
1359        t5_restart1(3)=potE
1360        t5_restart1(4)=t_bath
1361        t5_restart1(5)=Uconst
1362
1363        call mpi_gather(t5_restart1,5,mpi_real,
1364      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1365
1366
1367        do i=1,2*nres
1368          do j=1,3
1369            r_d(j,i)=d_t(j,i)
1370          enddo
1371        enddo
1372        call mpi_gather(r_d,3*2*nres,mpi_real,
1373      &           d_restart1,3*2*nres,mpi_real,king,
1374      &           CG_COMM,ierr)
1375
1376
1377        do i=1,2*nres
1378          do j=1,3
1379            r_d(j,i)=dc(j,i)
1380          enddo
1381        enddo
1382        call mpi_gather(r_d,3*2*nres,mpi_real,
1383      &           d_restart2,3*2*nres,mpi_real,king,
1384      &           CG_COMM,ierr)
1385
1386        if(me.eq.king) then
1387 #ifdef AIX
1388          call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1389          do i=0,nodes-1
1390           call xdrfint_(ixdrf, i2rep(i), iret)
1391          enddo
1392          do i=1,remd_m(1)
1393           call xdrfint_(ixdrf, ifirst(i), iret)
1394          enddo
1395          do il=1,nodes
1396               do i=0,nupa(0,il)
1397                call xdrfint_(ixdrf, nupa(i,il), iret)
1398               enddo
1399
1400               do i=0,ndowna(0,il)
1401                call xdrfint_(ixdrf, ndowna(i,il), iret)
1402               enddo
1403          enddo
1404
1405          do il=1,nodes
1406            do j=1,4
1407             call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1408            enddo
1409          enddo
1410
1411          do il=0,nodes-1
1412            do i=1,2*nres
1413             do j=1,3
1414              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1415             enddo
1416            enddo
1417          enddo
1418          do il=0,nodes-1
1419            do i=1,2*nres
1420             do j=1,3
1421              call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1422             enddo
1423            enddo
1424          enddo
1425
1426          if(usampl) then
1427            call xdrfint_(ixdrf, nset, iret)
1428            do i=1,nset
1429              call xdrfint_(ixdrf,mset(i), iret)
1430            enddo
1431            do i=0,nodes-1
1432              call xdrfint_(ixdrf,i2set(i), iret)
1433            enddo
1434            do il=1,nset
1435              do il1=1,mset(il)
1436                do i=1,nrep
1437                  do j=1,remd_m(i)
1438                    itmp=i_index(i,j,il,il1)
1439                    call xdrfint_(ixdrf,itmp, iret)
1440                  enddo
1441                enddo
1442              enddo
1443            enddo
1444
1445          endif
1446          call xdrfclose_(ixdrf, iret)
1447 #else
1448          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1449          do i=0,nodes-1
1450           call xdrfint(ixdrf, i2rep(i), iret)
1451          enddo
1452          do i=1,remd_m(1)
1453           call xdrfint(ixdrf, ifirst(i), iret)
1454          enddo
1455          do il=1,nodes
1456               do i=0,nupa(0,il)
1457                call xdrfint(ixdrf, nupa(i,il), iret)
1458               enddo
1459
1460               do i=0,ndowna(0,il)
1461                call xdrfint(ixdrf, ndowna(i,il), iret)
1462               enddo
1463          enddo
1464
1465          do il=1,nodes
1466            do j=1,4
1467             call xdrffloat(ixdrf, t_restart1(j,il), iret)
1468            enddo
1469          enddo
1470
1471          do il=0,nodes-1
1472            do i=1,2*nres
1473             do j=1,3
1474              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1475             enddo
1476            enddo
1477          enddo
1478          do il=0,nodes-1
1479            do i=1,2*nres
1480             do j=1,3
1481              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1482             enddo
1483            enddo
1484          enddo
1485
1486
1487              if(usampl) then
1488               call xdrfint(ixdrf, nset, iret)
1489               do i=1,nset
1490                 call xdrfint(ixdrf,mset(i), iret)
1491               enddo
1492               do i=0,nodes-1
1493                 call xdrfint(ixdrf,i2set(i), iret)
1494               enddo
1495               do il=1,nset
1496                do il1=1,mset(il)
1497                 do i=1,nrep
1498                  do j=1,remd_m(i)
1499                    itmp=i_index(i,j,il,il1)
1500                    call xdrfint(ixdrf,itmp, iret)
1501                  enddo
1502                 enddo
1503                enddo
1504               enddo
1505
1506              endif
1507          call xdrfclose(ixdrf, iret)
1508 #endif
1509        endif
1510       return
1511       end
1512
1513
1514       subroutine write1traj
1515       implicit none
1516       include 'DIMENSIONS'
1517       include 'mpif.h'
1518       include 'COMMON.MD'
1519       include 'COMMON.QRESTR'
1520       include 'COMMON.IOUNITS'
1521       include 'COMMON.REMD'
1522       include 'COMMON.SETUP'
1523       include 'COMMON.CHAIN'
1524       include 'COMMON.SBRIDGE'
1525       include 'COMMON.INTERACT'
1526
1527       real t5_restart1(5)
1528       integer iret,itmp
1529       real xcoord(3,maxres2+2),prec
1530       real r_qfrag(50),r_qpair(100)
1531       real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1532       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1533       real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1534      &     p_uscdiff(100*maxprocs)
1535       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1536       common /przechowalnia/ p_c
1537       integer ii,i,il,j,ixdrf
1538       integer ierr
1539
1540       call mpi_bcast(ii_write,1,mpi_integer,
1541      &           king,CG_COMM,ierr)
1542
1543 c debugging
1544 c      print *,'traj1file',me,ii_write,ntwx_cache
1545 c end debugging
1546
1547 #ifdef AIX
1548       if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1549 #else
1550       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1551 #endif
1552       do ii=1,ii_write
1553        t5_restart1(1)=totT_cache(ii)
1554        t5_restart1(2)=EK_cache(ii)
1555        t5_restart1(3)=potE_cache(ii)
1556        t5_restart1(4)=t_bath_cache(ii)
1557        t5_restart1(5)=Uconst_cache(ii)
1558        call mpi_gather(t5_restart1,5,mpi_real,
1559      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1560
1561        call mpi_gather(iset_cache(ii),1,mpi_integer,
1562      &      iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1563
1564           do i=1,nfrag
1565            r_qfrag(i)=qfrag_cache(i,ii)
1566           enddo
1567           do i=1,npair
1568            r_qpair(i)=qpair_cache(i,ii)
1569           enddo
1570           do i=1,nfrag_back
1571            r_utheta(i)=utheta_cache(i,ii)
1572            r_ugamma(i)=ugamma_cache(i,ii)
1573            r_uscdiff(i)=uscdiff_cache(i,ii)
1574           enddo
1575
1576         call mpi_gather(r_qfrag,nfrag,mpi_real,
1577      &           p_qfrag,nfrag,mpi_real,king,
1578      &           CG_COMM,ierr)
1579         call mpi_gather(r_qpair,npair,mpi_real,
1580      &           p_qpair,npair,mpi_real,king,
1581      &           CG_COMM,ierr)
1582         call mpi_gather(r_utheta,nfrag_back,mpi_real,
1583      &           p_utheta,nfrag_back,mpi_real,king,
1584      &           CG_COMM,ierr)
1585         call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1586      &           p_ugamma,nfrag_back,mpi_real,king,
1587      &           CG_COMM,ierr)
1588         call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1589      &           p_uscdiff,nfrag_back,mpi_real,king,
1590      &           CG_COMM,ierr)
1591
1592 #ifdef DEBUG
1593         write (iout,*) "p_qfrag"
1594         do i=1,nodes
1595           write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1596         enddo
1597         write (iout,*) "p_qpair"
1598         do i=1,nodes
1599           write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1600         enddo
1601         call flush(iout)
1602 #endif
1603         do i=1,nres*2
1604          do j=1,3
1605           r_c(j,i)=c_cache(j,i,ii)
1606          enddo
1607         enddo
1608
1609         call mpi_gather(r_c,3*2*nres,mpi_real,
1610      &           p_c,3*2*nres,mpi_real,king,
1611      &           CG_COMM,ierr)
1612
1613        if(me.eq.king) then
1614 #ifdef AIX
1615          do il=1,nodes
1616           call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1617           call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1618           call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1619           call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1620           call xdrfint_(ixdrf, nss, iret)
1621           do j=1,nss
1622            if (dyn_ss) then
1623             call xdrfint(ixdrf, idssb(j)+nres, iret)
1624             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1625            else
1626             call xdrfint_(ixdrf, ihpb(j), iret)
1627             call xdrfint_(ixdrf, jhpb(j), iret)
1628            endif
1629           enddo
1630           call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1631           call xdrfint_(ixdrf, iset_restart1(il), iret)
1632           do i=1,nfrag
1633            call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1634           enddo
1635           do i=1,npair
1636            call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1637           enddo
1638           do i=1,nfrag_back
1639            call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1640            call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1641            call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1642           enddo
1643           prec=10000.0
1644           do i=1,nres
1645            do j=1,3
1646             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1647            enddo
1648           enddo
1649           do i=nnt,nct
1650            do j=1,3
1651             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1652            enddo
1653           enddo
1654           itmp=nres+nct-nnt+1
1655           call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1656          enddo
1657 #else
1658          do il=1,nodes
1659           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1660           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1661           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1662           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1663           call xdrfint(ixdrf, nss, iret)
1664           do j=1,nss
1665            if (dyn_ss) then
1666             call xdrfint(ixdrf, idssb(j)+nres, iret)
1667             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1668            else
1669             call xdrfint(ixdrf, ihpb(j), iret)
1670             call xdrfint(ixdrf, jhpb(j), iret)
1671            endif
1672           enddo
1673           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1674           call xdrfint(ixdrf, iset_restart1(il), iret)
1675           do i=1,nfrag
1676            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1677           enddo
1678           do i=1,npair
1679            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1680           enddo
1681           do i=1,nfrag_back
1682            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1683            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1684            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1685           enddo
1686           prec=10000.0
1687           do i=1,nres
1688            do j=1,3
1689             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1690            enddo
1691           enddo
1692           do i=nnt,nct
1693            do j=1,3
1694             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1695            enddo
1696           enddo
1697           itmp=nres+nct-nnt+1
1698           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1699          enddo
1700 #endif
1701        endif
1702       enddo
1703 #ifdef AIX
1704       if(me.eq.king) call xdrfclose_(ixdrf, iret)
1705 #else
1706       if(me.eq.king) call xdrfclose(ixdrf, iret)
1707 #endif
1708       do i=1,ntwx_cache-ii_write
1709
1710             totT_cache(i)=totT_cache(ii_write+i)
1711             EK_cache(i)=EK_cache(ii_write+i)
1712             potE_cache(i)=potE_cache(ii_write+i)
1713             t_bath_cache(i)=t_bath_cache(ii_write+i)
1714             Uconst_cache(i)=Uconst_cache(ii_write+i)
1715             iset_cache(i)=iset_cache(ii_write+i)
1716
1717             do ii=1,nfrag
1718              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1719             enddo
1720             do ii=1,npair
1721              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1722             enddo
1723             do ii=1,nfrag_back
1724               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1725               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1726               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1727             enddo
1728
1729             do ii=1,nres*2
1730              do j=1,3
1731               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1732              enddo
1733             enddo
1734       enddo
1735       ntwx_cache=ntwx_cache-ii_write
1736       return
1737       end
1738
1739
1740       subroutine read1restart(i_index)
1741       implicit none
1742       include 'DIMENSIONS'
1743       include 'mpif.h'
1744       include 'COMMON.CONTROL'
1745       include 'COMMON.MD'
1746 #ifdef FIVEDIAG
1747        include 'COMMON.LAGRANGE.5diag'
1748 #else
1749        include 'COMMON.LAGRANGE'
1750 #endif
1751       include 'COMMON.QRESTR'
1752       include 'COMMON.IOUNITS'
1753       include 'COMMON.REMD'
1754       include 'COMMON.SETUP'
1755       include 'COMMON.CHAIN'
1756       include 'COMMON.SBRIDGE'
1757       include 'COMMON.INTERACT'
1758       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1759      &                 t5_restart1(5)
1760       integer*2 i_index
1761      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1762       common /przechowalnia/ d_restart1
1763       integer i,j,il,il1,ixdrf,iret,itmp
1764       integer ierr
1765       write (*,*) "Processor",me," called read1restart"
1766
1767          if(me.eq.king)then
1768               open(irest2,file=mremd_rst_name,status='unknown')
1769               read(irest2,*,err=334) i
1770               write(iout,*) "Reading old rst in ASCI format"
1771               close(irest2)
1772                call read1restart_old
1773                return
1774  334          continue
1775 #ifdef AIX
1776               call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1777
1778               do i=0,nodes-1
1779                call xdrfint_(ixdrf, i2rep(i), iret)
1780               enddo
1781               do i=1,remd_m(1)
1782                call xdrfint_(ixdrf, ifirst(i), iret)
1783               enddo
1784              do il=1,nodes
1785               call xdrfint_(ixdrf, nupa(0,il), iret)
1786               do i=1,nupa(0,il)
1787                call xdrfint_(ixdrf, nupa(i,il), iret)
1788               enddo
1789
1790               call xdrfint_(ixdrf, ndowna(0,il), iret)
1791               do i=1,ndowna(0,il)
1792                call xdrfint_(ixdrf, ndowna(i,il), iret)
1793               enddo
1794              enddo
1795              do il=1,nodes
1796                do j=1,4
1797                 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1798                enddo
1799              enddo
1800 #else
1801               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1802
1803               do i=0,nodes-1
1804                call xdrfint(ixdrf, i2rep(i), iret)
1805               enddo
1806               do i=1,remd_m(1)
1807                call xdrfint(ixdrf, ifirst(i), iret)
1808               enddo
1809              do il=1,nodes
1810               call xdrfint(ixdrf, nupa(0,il), iret)
1811               do i=1,nupa(0,il)
1812                call xdrfint(ixdrf, nupa(i,il), iret)
1813               enddo
1814
1815               call xdrfint(ixdrf, ndowna(0,il), iret)
1816               do i=1,ndowna(0,il)
1817                call xdrfint(ixdrf, ndowna(i,il), iret)
1818               enddo
1819              enddo
1820              do il=1,nodes
1821                do j=1,4
1822                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1823                enddo
1824              enddo
1825 #endif
1826          endif
1827          call mpi_scatter(t_restart1,5,mpi_real,
1828      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1829          totT=t5_restart1(1)
1830          EK=t5_restart1(2)
1831          potE=t5_restart1(3)
1832          t_bath=t5_restart1(4)
1833
1834          if(me.eq.king)then
1835               do il=0,nodes-1
1836                do i=1,2*nres
1837 c                read(irest2,'(3e15.5)')
1838 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1839             do j=1,3
1840 #ifdef AIX
1841              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1842 #else
1843              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1844 #endif
1845             enddo
1846                enddo
1847               enddo
1848          endif
1849          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1850      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1851
1852          do i=1,2*nres
1853            do j=1,3
1854             d_t(j,i)=r_d(j,i)
1855            enddo
1856          enddo
1857          if(me.eq.king)then
1858               do il=0,nodes-1
1859                do i=1,2*nres
1860 c                read(irest2,'(3e15.5)')
1861 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1862             do j=1,3
1863 #ifdef AIX
1864              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1865 #else
1866              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1867 #endif
1868             enddo
1869                enddo
1870               enddo
1871          endif
1872          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1873      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1874          do i=1,2*nres
1875            do j=1,3
1876             dc(j,i)=r_d(j,i)
1877            enddo
1878          enddo
1879
1880
1881            if(usampl) then
1882 #ifdef AIX
1883              if(me.eq.king)then
1884               call xdrfint_(ixdrf, nset, iret)
1885               do i=1,nset
1886                 call xdrfint_(ixdrf,mset(i), iret)
1887               enddo
1888               do i=0,nodes-1
1889                 call xdrfint_(ixdrf,i2set(i), iret)
1890               enddo
1891               do il=1,nset
1892                do il1=1,mset(il)
1893                 do i=1,nrep
1894                  do j=1,remd_m(i)
1895                    call xdrfint_(ixdrf,itmp, iret)
1896                    i_index(i,j,il,il1)=itmp
1897                  enddo
1898                 enddo
1899                enddo
1900               enddo
1901              endif
1902 #else
1903              if(me.eq.king)then
1904               call xdrfint(ixdrf, nset, iret)
1905               do i=1,nset
1906                 call xdrfint(ixdrf,mset(i), iret)
1907               enddo
1908               do i=0,nodes-1
1909                 call xdrfint(ixdrf,i2set(i), iret)
1910               enddo
1911               do il=1,nset
1912                do il1=1,mset(il)
1913                 do i=1,nrep
1914                  do j=1,remd_m(i)
1915                    call xdrfint(ixdrf,itmp, iret)
1916                    i_index(i,j,il,il1)=itmp
1917                  enddo
1918                 enddo
1919                enddo
1920               enddo
1921              endif
1922 #endif
1923 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1924 c own element
1925 c              call mpi_scatter(i2set,1,mpi_integer,
1926 c     &           iset,1,mpi_integer,king,
1927 c     &           CG_COMM,ierr)
1928               call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1929      &         CG_COMM,ierr)
1930               iset=i2set(me)
1931
1932            endif
1933
1934
1935         if(me.eq.king) close(irest2)
1936         return
1937         end
1938
1939       subroutine read1restart_old
1940       implicit none
1941       include 'DIMENSIONS'
1942       include 'mpif.h'
1943       include 'COMMON.MD'
1944 #ifdef FIVEDIAG
1945        include 'COMMON.LAGRANGE.5diag'
1946 #else
1947        include 'COMMON.LAGRANGE'
1948 #endif
1949       include 'COMMON.IOUNITS'
1950       include 'COMMON.REMD'
1951       include 'COMMON.SETUP'
1952       include 'COMMON.CHAIN'
1953       include 'COMMON.SBRIDGE'
1954       include 'COMMON.INTERACT'
1955       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1956      &                 t5_restart1(5)
1957       common /przechowalnia/ d_restart1
1958       integer i,j,il,itmp
1959       integer ierr
1960          if(me.eq.king)then
1961              open(irest2,file=mremd_rst_name,status='unknown')
1962              read (irest2,*) (i2rep(i),i=0,nodes-1)
1963              read (irest2,*) (ifirst(i),i=1,remd_m(1))
1964              do il=1,nodes
1965               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1966               read (irest2,*) ndowna(0,il),
1967      &                    (ndowna(i,il),i=1,ndowna(0,il))
1968              enddo
1969              do il=1,nodes
1970                read(irest2,*) (t_restart1(j,il),j=1,4)
1971              enddo
1972          endif
1973          call mpi_scatter(t_restart1,5,mpi_real,
1974      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1975          totT=t5_restart1(1)
1976          EK=t5_restart1(2)
1977          potE=t5_restart1(3)
1978          t_bath=t5_restart1(4)
1979
1980          if(me.eq.king)then
1981               do il=0,nodes-1
1982                do i=1,2*nres
1983                 read(irest2,'(3e15.5)')
1984      &                (d_restart1(j,i+2*nres*il),j=1,3)
1985                enddo
1986               enddo
1987          endif
1988          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1989      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1990
1991          do i=1,2*nres
1992            do j=1,3
1993             d_t(j,i)=r_d(j,i)
1994            enddo
1995          enddo
1996          if(me.eq.king)then
1997               do il=0,nodes-1
1998                do i=1,2*nres
1999                 read(irest2,'(3e15.5)')
2000      &                (d_restart1(j,i+2*nres*il),j=1,3)
2001                enddo
2002               enddo
2003          endif
2004          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2005      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2006          do i=1,2*nres
2007            do j=1,3
2008             dc(j,i)=r_d(j,i)
2009            enddo
2010          enddo
2011         if(me.eq.king) close(irest2)
2012         return
2013         end
2014 c------------------------------------------