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