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=1,2*nres
 576             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 577            enddo
 578            do i=1,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*maxprocs),r_d(3,2*maxres),
1236      &     d_restart2(3,2*maxres*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=1,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,mpi_real,
1259      &           d_restart1,3*2*nres,mpi_real,king,
1260      &           CG_COMM,ierr)
1261
1262
1263        do i=1,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,mpi_real,
1269      &           d_restart2,3*2*nres,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=1,2*nres
1299             do j=1,3
1300              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1301             enddo
1302            enddo
1303          enddo
1304          do il=0,nodes-1
1305            do i=1,2*nres
1306             do j=1,3
1307              call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1308             enddo
1309            enddo
1310          enddo
1311
1312          if(usampl) then
1313            call xdrfint_(ixdrf, nset, iret)
1314            do i=1,nset
1315              call xdrfint_(ixdrf,mset(i), iret)
1316            enddo
1317            do i=0,nodes-1
1318              call xdrfint_(ixdrf,i2set(i), iret)
1319            enddo
1320            do il=1,nset
1321              do il1=1,mset(il)
1322                do i=1,nrep
1323                  do j=1,remd_m(i)
1324                    itmp=i_index(i,j,il,il1)
1325                    call xdrfint_(ixdrf,itmp, iret)
1326                  enddo
1327                enddo
1328              enddo
1329            enddo
1330
1331          endif
1332          call xdrfclose_(ixdrf, iret)
1333 #else
1334          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1335          do i=0,nodes-1
1336           call xdrfint(ixdrf, i2rep(i), iret)
1337          enddo
1338          do i=1,remd_m(1)
1339           call xdrfint(ixdrf, ifirst(i), iret)
1340          enddo
1341          do il=1,nodes
1342               do i=0,nupa(0,il)
1343                call xdrfint(ixdrf, nupa(i,il), iret)
1344               enddo
1345
1346               do i=0,ndowna(0,il)
1347                call xdrfint(ixdrf, ndowna(i,il), iret)
1348               enddo
1349          enddo
1350
1351          do il=1,nodes
1352            do j=1,4
1353             call xdrffloat(ixdrf, t_restart1(j,il), iret)
1354            enddo
1355          enddo
1356
1357          do il=0,nodes-1
1358            do i=1,2*nres
1359             do j=1,3
1360              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1361             enddo
1362            enddo
1363          enddo
1364          do il=0,nodes-1
1365            do i=1,2*nres
1366             do j=1,3
1367              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1368             enddo
1369            enddo
1370          enddo
1371
1372
1373              if(usampl) then
1374               call xdrfint(ixdrf, nset, iret)
1375               do i=1,nset
1376                 call xdrfint(ixdrf,mset(i), iret)
1377               enddo
1378               do i=0,nodes-1
1379                 call xdrfint(ixdrf,i2set(i), iret)
1380               enddo
1381               do il=1,nset
1382                do il1=1,mset(il)
1383                 do i=1,nrep
1384                  do j=1,remd_m(i)
1385                    itmp=i_index(i,j,il,il1)
1386                    call xdrfint(ixdrf,itmp, iret)
1387                  enddo
1388                 enddo
1389                enddo
1390               enddo
1391
1392              endif
1393          call xdrfclose(ixdrf, iret)
1394 #endif
1395        endif
1396       return
1397       end
1398
1399
1400       subroutine write1traj
1401       implicit real*8 (a-h,o-z)
1402       include 'DIMENSIONS'
1403       include 'mpif.h'
1404       include 'COMMON.MD'
1405       include 'COMMON.IOUNITS'
1406       include 'COMMON.REMD'
1407       include 'COMMON.SETUP'
1408       include 'COMMON.CHAIN'
1409       include 'COMMON.SBRIDGE'
1410       include 'COMMON.INTERACT'
1411
1412       real t5_restart1(5)
1413       integer iret,itmp
1414       real xcoord(3,maxres2+2),prec
1415       real r_qfrag(50),r_qpair(100)
1416       real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1417       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1418       real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1419      &     p_uscdiff(100*maxprocs)
1420       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1421       common /przechowalnia/ p_c
1422
1423       call mpi_bcast(ii_write,1,mpi_integer,
1424      &           king,CG_COMM,ierr)
1425
1426 c debugging
1427       print *,'traj1file',me,ii_write,ntwx_cache
1428 c end debugging
1429
1430 #ifdef AIX
1431       if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1432 #else
1433       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1434 #endif
1435       do ii=1,ii_write
1436        t5_restart1(1)=totT_cache(ii)
1437        t5_restart1(2)=EK_cache(ii)
1438        t5_restart1(3)=potE_cache(ii)
1439        t5_restart1(4)=t_bath_cache(ii)
1440        t5_restart1(5)=Uconst_cache(ii)
1441        call mpi_gather(t5_restart1,5,mpi_real,
1442      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1443
1444        call mpi_gather(iset_cache(ii),1,mpi_integer,
1445      &      iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1446
1447           do i=1,nfrag
1448            r_qfrag(i)=qfrag_cache(i,ii)
1449           enddo
1450           do i=1,npair
1451            r_qpair(i)=qpair_cache(i,ii)
1452           enddo
1453           do i=1,nfrag_back
1454            r_utheta(i)=utheta_cache(i,ii)
1455            r_ugamma(i)=ugamma_cache(i,ii)
1456            r_uscdiff(i)=uscdiff_cache(i,ii)
1457           enddo
1458
1459         call mpi_gather(r_qfrag,nfrag,mpi_real,
1460      &           p_qfrag,nfrag,mpi_real,king,
1461      &           CG_COMM,ierr)
1462         call mpi_gather(r_qpair,npair,mpi_real,
1463      &           p_qpair,npair,mpi_real,king,
1464      &           CG_COMM,ierr)
1465         call mpi_gather(r_utheta,nfrag_back,mpi_real,
1466      &           p_utheta,nfrag_back,mpi_real,king,
1467      &           CG_COMM,ierr)
1468         call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1469      &           p_ugamma,nfrag_back,mpi_real,king,
1470      &           CG_COMM,ierr)
1471         call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1472      &           p_uscdiff,nfrag_back,mpi_real,king,
1473      &           CG_COMM,ierr)
1474
1475 #ifdef DEBUG
1476         write (iout,*) "p_qfrag"
1477         do i=1,nodes
1478           write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1479         enddo
1480         write (iout,*) "p_qpair"
1481         do i=1,nodes
1482           write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1483         enddo
1484         call flush(iout)
1485 #endif
1486         do i=1,nres*2
1487          do j=1,3
1488           r_c(j,i)=c_cache(j,i,ii)
1489          enddo
1490         enddo
1491
1492         call mpi_gather(r_c,3*2*nres,mpi_real,
1493      &           p_c,3*2*nres,mpi_real,king,
1494      &           CG_COMM,ierr)
1495
1496        if(me.eq.king) then
1497 #ifdef AIX
1498          do il=1,nodes
1499           call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1500           call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1501           call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1502           call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1503           call xdrfint_(ixdrf, nss, iret)
1504           do j=1,nss
1505            if (dyn_ss) then
1506             call xdrfint(ixdrf, idssb(j)+nres, iret)
1507             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1508            else
1509             call xdrfint_(ixdrf, ihpb(j), iret)
1510             call xdrfint_(ixdrf, jhpb(j), iret)
1511            endif
1512           enddo
1513           call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1514           call xdrfint_(ixdrf, iset_restart1(il), iret)
1515           do i=1,nfrag
1516            call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1517           enddo
1518           do i=1,npair
1519            call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1520           enddo
1521           do i=1,nfrag_back
1522            call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1523            call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1524            call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1525           enddo
1526           prec=10000.0
1527           do i=1,nres
1528            do j=1,3
1529             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1530            enddo
1531           enddo
1532           do i=nnt,nct
1533            do j=1,3
1534             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1535            enddo
1536           enddo
1537           itmp=nres+nct-nnt+1
1538           call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1539          enddo
1540 #else
1541          do il=1,nodes
1542           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1543           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1544           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1545           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1546           call xdrfint(ixdrf, nss, iret)
1547           do j=1,nss
1548            if (dyn_ss) then
1549             call xdrfint(ixdrf, idssb(j)+nres, iret)
1550             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1551            else
1552             call xdrfint(ixdrf, ihpb(j), iret)
1553             call xdrfint(ixdrf, jhpb(j), iret)
1554            endif
1555           enddo
1556           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1557           call xdrfint(ixdrf, iset_restart1(il), iret)
1558           do i=1,nfrag
1559            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1560           enddo
1561           do i=1,npair
1562            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1563           enddo
1564           do i=1,nfrag_back
1565            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1566            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1567            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1568           enddo
1569           prec=10000.0
1570           do i=1,nres
1571            do j=1,3
1572             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1573            enddo
1574           enddo
1575           do i=nnt,nct
1576            do j=1,3
1577             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1578            enddo
1579           enddo
1580           itmp=nres+nct-nnt+1
1581           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1582          enddo
1583 #endif
1584        endif
1585       enddo
1586 #ifdef AIX
1587       if(me.eq.king) call xdrfclose_(ixdrf, iret)
1588 #else
1589       if(me.eq.king) call xdrfclose(ixdrf, iret)
1590 #endif
1591       do i=1,ntwx_cache-ii_write
1592
1593             totT_cache(i)=totT_cache(ii_write+i)
1594             EK_cache(i)=EK_cache(ii_write+i)
1595             potE_cache(i)=potE_cache(ii_write+i)
1596             t_bath_cache(i)=t_bath_cache(ii_write+i)
1597             Uconst_cache(i)=Uconst_cache(ii_write+i)
1598             iset_cache(i)=iset_cache(ii_write+i)
1599
1600             do ii=1,nfrag
1601              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1602             enddo
1603             do ii=1,npair
1604              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1605             enddo
1606             do ii=1,nfrag_back
1607               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1608               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1609               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1610             enddo
1611
1612             do ii=1,nres*2
1613              do j=1,3
1614               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1615              enddo
1616             enddo
1617       enddo
1618       ntwx_cache=ntwx_cache-ii_write
1619       return
1620       end
1621
1622
1623       subroutine read1restart(i_index)
1624       implicit real*8 (a-h,o-z)
1625       include 'DIMENSIONS'
1626       include 'mpif.h'
1627       include 'COMMON.MD'
1628       include 'COMMON.IOUNITS'
1629       include 'COMMON.REMD'
1630       include 'COMMON.SETUP'
1631       include 'COMMON.CHAIN'
1632       include 'COMMON.SBRIDGE'
1633       include 'COMMON.INTERACT'
1634       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1635      &                 t5_restart1(5)
1636       integer*2 i_index
1637      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1638       common /przechowalnia/ d_restart1
1639       write (*,*) "Processor",me," called read1restart"
1640
1641          if(me.eq.king)then
1642               open(irest2,file=mremd_rst_name,status='unknown')
1643               read(irest2,*,err=334) i
1644               write(iout,*) "Reading old rst in ASCI format"
1645               close(irest2)
1646                call read1restart_old
1647                return
1648  334          continue
1649 #ifdef AIX
1650               call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1651
1652               do i=0,nodes-1
1653                call xdrfint_(ixdrf, i2rep(i), iret)
1654               enddo
1655               do i=1,remd_m(1)
1656                call xdrfint_(ixdrf, ifirst(i), iret)
1657               enddo
1658              do il=1,nodes
1659               call xdrfint_(ixdrf, nupa(0,il), iret)
1660               do i=1,nupa(0,il)
1661                call xdrfint_(ixdrf, nupa(i,il), iret)
1662               enddo
1663
1664               call xdrfint_(ixdrf, ndowna(0,il), iret)
1665               do i=1,ndowna(0,il)
1666                call xdrfint_(ixdrf, ndowna(i,il), iret)
1667               enddo
1668              enddo
1669              do il=1,nodes
1670                do j=1,4
1671                 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1672                enddo
1673              enddo
1674 #else
1675               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1676
1677               do i=0,nodes-1
1678                call xdrfint(ixdrf, i2rep(i), iret)
1679               enddo
1680               do i=1,remd_m(1)
1681                call xdrfint(ixdrf, ifirst(i), iret)
1682               enddo
1683              do il=1,nodes
1684               call xdrfint(ixdrf, nupa(0,il), iret)
1685               do i=1,nupa(0,il)
1686                call xdrfint(ixdrf, nupa(i,il), iret)
1687               enddo
1688
1689               call xdrfint(ixdrf, ndowna(0,il), iret)
1690               do i=1,ndowna(0,il)
1691                call xdrfint(ixdrf, ndowna(i,il), iret)
1692               enddo
1693              enddo
1694              do il=1,nodes
1695                do j=1,4
1696                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1697                enddo
1698              enddo
1699 #endif
1700          endif
1701          call mpi_scatter(t_restart1,5,mpi_real,
1702      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1703          totT=t5_restart1(1)
1704          EK=t5_restart1(2)
1705          potE=t5_restart1(3)
1706          t_bath=t5_restart1(4)
1707
1708          if(me.eq.king)then
1709               do il=0,nodes-1
1710                do i=1,2*nres
1711 c                read(irest2,'(3e15.5)')
1712 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1713             do j=1,3
1714 #ifdef AIX
1715              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1716 #else
1717              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1718 #endif
1719             enddo
1720                enddo
1721               enddo
1722          endif
1723          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1724      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1725
1726          do i=1,2*nres
1727            do j=1,3
1728             d_t(j,i)=r_d(j,i)
1729            enddo
1730          enddo
1731          if(me.eq.king)then
1732               do il=0,nodes-1
1733                do i=1,2*nres
1734 c                read(irest2,'(3e15.5)')
1735 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1736             do j=1,3
1737 #ifdef AIX
1738              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1739 #else
1740              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1741 #endif
1742             enddo
1743                enddo
1744               enddo
1745          endif
1746          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1747      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1748          do i=1,2*nres
1749            do j=1,3
1750             dc(j,i)=r_d(j,i)
1751            enddo
1752          enddo
1753
1754
1755            if(usampl) then
1756 #ifdef AIX
1757              if(me.eq.king)then
1758               call xdrfint_(ixdrf, nset, iret)
1759               do i=1,nset
1760                 call xdrfint_(ixdrf,mset(i), iret)
1761               enddo
1762               do i=0,nodes-1
1763                 call xdrfint_(ixdrf,i2set(i), iret)
1764               enddo
1765               do il=1,nset
1766                do il1=1,mset(il)
1767                 do i=1,nrep
1768                  do j=1,remd_m(i)
1769                    call xdrfint_(ixdrf,itmp, iret)
1770                    i_index(i,j,il,il1)=itmp
1771                  enddo
1772                 enddo
1773                enddo
1774               enddo
1775              endif
1776 #else
1777              if(me.eq.king)then
1778               call xdrfint(ixdrf, nset, iret)
1779               do i=1,nset
1780                 call xdrfint(ixdrf,mset(i), iret)
1781               enddo
1782               do i=0,nodes-1
1783                 call xdrfint(ixdrf,i2set(i), iret)
1784               enddo
1785               do il=1,nset
1786                do il1=1,mset(il)
1787                 do i=1,nrep
1788                  do j=1,remd_m(i)
1789                    call xdrfint(ixdrf,itmp, iret)
1790                    i_index(i,j,il,il1)=itmp
1791                  enddo
1792                 enddo
1793                enddo
1794               enddo
1795              endif
1796 #endif
1797               call mpi_scatter(i2set,1,mpi_integer,
1798      &           iset,1,mpi_integer,king,
1799      &           CG_COMM,ierr)
1800
1801            endif
1802
1803
1804         if(me.eq.king) close(irest2)
1805         return
1806         end
1807
1808       subroutine read1restart_old
1809       implicit real*8 (a-h,o-z)
1810       include 'DIMENSIONS'
1811       include 'mpif.h'
1812       include 'COMMON.MD'
1813       include 'COMMON.IOUNITS'
1814       include 'COMMON.REMD'
1815       include 'COMMON.SETUP'
1816       include 'COMMON.CHAIN'
1817       include 'COMMON.SBRIDGE'
1818       include 'COMMON.INTERACT'
1819       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1820      &                 t5_restart1(5)
1821       common /przechowalnia/ d_restart1
1822          if(me.eq.king)then
1823              open(irest2,file=mremd_rst_name,status='unknown')
1824              read (irest2,*) (i2rep(i),i=0,nodes-1)
1825              read (irest2,*) (ifirst(i),i=1,remd_m(1))
1826              do il=1,nodes
1827               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1828               read (irest2,*) ndowna(0,il),
1829      &                    (ndowna(i,il),i=1,ndowna(0,il))
1830              enddo
1831              do il=1,nodes
1832                read(irest2,*) (t_restart1(j,il),j=1,4)
1833              enddo
1834          endif
1835          call mpi_scatter(t_restart1,5,mpi_real,
1836      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1837          totT=t5_restart1(1)
1838          EK=t5_restart1(2)
1839          potE=t5_restart1(3)
1840          t_bath=t5_restart1(4)
1841
1842          if(me.eq.king)then
1843               do il=0,nodes-1
1844                do i=1,2*nres
1845                 read(irest2,'(3e15.5)')
1846      &                (d_restart1(j,i+2*nres*il),j=1,3)
1847                enddo
1848               enddo
1849          endif
1850          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1851      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1852
1853          do i=1,2*nres
1854            do j=1,3
1855             d_t(j,i)=r_d(j,i)
1856            enddo
1857          enddo
1858          if(me.eq.king)then
1859               do il=0,nodes-1
1860                do i=1,2*nres
1861                 read(irest2,'(3e15.5)')
1862      &                (d_restart1(j,i+2*nres*il),j=1,3)
1863                enddo
1864               enddo
1865          endif
1866          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1867      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1868          do i=1,2*nres
1869            do j=1,3
1870             dc(j,i)=r_d(j,i)
1871            enddo
1872          enddo
1873         if(me.eq.king) close(irest2)
1874         return
1875         end
1876 c------------------------------------------
1877       subroutine returnbox
1878       include 'DIMENSIONS'
1879       include 'mpif.h'
1880       include 'COMMON.CONTROL'
1881       include 'COMMON.VAR'
1882       include 'COMMON.MD'
1883 #ifndef LANG0
1884       include 'COMMON.LANGEVIN'
1885 #else
1886       include 'COMMON.LANGEVIN.lang0'
1887 #endif
1888       include 'COMMON.CHAIN'
1889       include 'COMMON.DERIV'
1890       include 'COMMON.GEO'
1891       include 'COMMON.LOCAL'
1892       include 'COMMON.INTERACT'
1893       include 'COMMON.IOUNITS'
1894       include 'COMMON.NAMES'
1895       include 'COMMON.TIME1'
1896       include 'COMMON.REMD'
1897       include 'COMMON.SETUP'
1898       include 'COMMON.MUCA'
1899       include 'COMMON.HAIRPIN'
1900         j=1
1901         chain_beg=1
1902         allareout=1
1903 C        do i=1,nres
1904 C       write(*,*) 'initial', i,j,c(j,i)
1905 C        enddo
1906         do i=1,nres-1
1907          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
1908           chain_end=i
1909           if (allareout.eq.1) then
1910             ireturnval=int(c(j,i)/boxxsize)
1911             if (c(j,i).le.0) ireturnval=ireturnval-1
1912             do k=chain_beg,chain_end
1913               c(j,k)=c(j,k)-ireturnval*boxxsize
1914               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
1915             enddo
1916            endif
1917            chain_beg=i+1
1918            allareout=1
1919          else
1920           if (int(c(j,i)/boxxsize).eq.0) allareout=0
1921          endif
1922         enddo
1923          if (allareout.eq.1) then
1924             ireturnval=int(c(j,i)/boxxsize)
1925             if (c(j,i).le.0) ireturnval=ireturnval-1
1926             do k=chain_beg,nres
1927               c(j,k)=c(j,k)-ireturnval*boxxsize
1928               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
1929             enddo
1930           endif
1931 C NO JUMP
1932 C        do i=1,nres
1933 C        write(*,*) 'befor no jump', i,j,c(j,i)
1934 C        enddo
1935         nojumpval=0
1936         do i=2,nres
1937            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1938              difference=abs(c(j,i-1)-c(j,i))
1939 C             print *,'diff', difference
1940              if (difference.gt.boxxsize/2.0) then
1941                 if (c(j,i-1).gt.c(j,i)) then
1942                   nojumpval=1
1943                  else
1944                    nojumpval=-1
1945                  endif
1946               else
1947               nojumpval=0
1948               endif
1949               endif
1950               c(j,i)=c(j,i)+nojumpval*boxxsize
1951               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
1952          enddo
1953        nojumpval=0
1954         do i=2,nres
1955            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1956              difference=abs(c(j,i-1)-c(j,i))
1957              if (difference.gt.boxxsize/2.0) then
1958                 if (c(j,i-1).gt.c(j,i)) then
1959                   nojumpval=1
1960                  else
1961                    nojumpval=-1
1962                  endif
1963               else
1964               nojumpval=0
1965               endif
1966              endif
1967               c(j,i)=c(j,i)+nojumpval*boxxsize
1968               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
1969          enddo
1970
1971 C        do i=1,nres
1972 C        write(*,*) 'after no jump', i,j,c(j,i)
1973 C        enddo
1974
1975 C NOW Y dimension
1976         j=2
1977         chain_beg=1
1978         do i=1,nres-1
1979          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
1980           chain_end=i
1981           if (allareout.eq.1) then
1982             ireturnval=int(c(j,i)/boxysize)
1983             if (c(j,i).le.0) ireturnval=ireturnval-1
1984             do k=chain_beg,chain_end
1985               c(j,k)=c(j,k)-ireturnval*boxysize
1986               c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
1987             enddo
1988            endif
1989            chain_beg=i+1
1990            allareout=1
1991          else
1992           if (int(c(j,i)/boxysize).eq.0) allareout=0
1993          endif
1994         enddo
1995          if (allareout.eq.1) then
1996             ireturnval=int(c(j,i)/boxysize)
1997             if (c(j,i).le.0) ireturnval=ireturnval-1
1998             do k=chain_beg,nres
1999               c(j,k)=c(j,k)-ireturnval*boxysize
2000               c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
2001             enddo
2002           endif
2003         nojumpval=0
2004         do i=2,nres
2005            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2006              difference=abs(c(j,i-1)-c(j,i))
2007              if (difference.gt.boxysize/2.0) then
2008                 if (c(j,i-1).gt.c(j,i)) then
2009                   nojumpval=1
2010                  else
2011                    nojumpval=-1
2012                  endif
2013               else
2014               nojumpval=0
2015               endif
2016            endif
2017               c(j,i)=c(j,i)+nojumpval*boxysize
2018               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
2019          enddo
2020       nojumpval=0
2021         do i=2,nres
2022            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2023              difference=abs(c(j,i-1)-c(j,i))
2024              if (difference.gt.boxysize/2.0) then
2025                 if (c(j,i-1).gt.c(j,i)) then
2026                   nojumpval=1
2027                  else
2028                    nojumpval=-1
2029                  endif
2030               else
2031               nojumpval=0
2032               endif
2033             endif
2034               c(j,i)=c(j,i)+nojumpval*boxysize
2035               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
2036          enddo
2037
2038         j=3
2039         chain_beg=1
2040         do i=1,nres-1
2041          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
2042           chain_end=i
2043           if (allareout.eq.1) then
2044             ireturnval=int(c(j,i)/boxysize)
2045             if (c(j,i).le.0) ireturnval=ireturnval-1
2046             do k=chain_beg,chain_end
2047               c(j,k)=c(j,k)-ireturnval*boxzsize
2048               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
2049             enddo
2050            endif
2051            chain_beg=i+1
2052            allareout=1
2053          else
2054           if (int(c(j,i)/boxzsize).eq.0) allareout=0
2055          endif
2056         enddo
2057          if (allareout.eq.1) then
2058             ireturnval=int(c(j,i)/boxzsize)
2059             if (c(j,i).le.0) ireturnval=ireturnval-1
2060             do k=chain_beg,nres
2061               c(j,k)=c(j,k)-ireturnval*boxzsize
2062               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
2063             enddo
2064           endif
2065         nojumpval=0
2066         do i=2,nres
2067            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2068              difference=abs(c(j,i-1)-c(j,i))
2069              if (difference.gt.(boxzsize/2.0)) then
2070                 if (c(j,i-1).gt.c(j,i)) then
2071                   nojumpval=1
2072                  else
2073                    nojumpval=-1
2074                  endif
2075               else
2076               nojumpval=0
2077               endif
2078             endif
2079               c(j,i)=c(j,i)+nojumpval*boxzsize
2080               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
2081          enddo
2082        nojumpval=0
2083         do i=2,nres
2084            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2085              difference=abs(c(j,i-1)-c(j,i))
2086              if (difference.gt.boxzsize/2.0) then
2087                 if (c(j,i-1).gt.c(j,i)) then
2088                   nojumpval=1
2089                  else
2090                    nojumpval=-1
2091                  endif
2092               else
2093               nojumpval=0
2094               endif
2095             endif
2096               c(j,i)=c(j,i)+nojumpval*boxzsize
2097               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
2098          enddo
2099
2100         return
2101         end