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