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