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