source/unres/src-5hdiag-tmp/MREMD.F

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