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