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            call xdrfint_(ixdrf, ihpb(j), iret)
1506            call xdrfint_(ixdrf, jhpb(j), iret)
1507           enddo
1508           call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1509           call xdrfint_(ixdrf, iset_restart1(il), iret)
1510           do i=1,nfrag
1511            call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1512           enddo
1513           do i=1,npair
1514            call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1515           enddo
1516           do i=1,nfrag_back
1517            call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1518            call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1519            call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1520           enddo
1521           prec=10000.0
1522           do i=1,nres
1523            do j=1,3
1524             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1525            enddo
1526           enddo
1527           do i=nnt,nct
1528            do j=1,3
1529             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1530            enddo
1531           enddo
1532           itmp=nres+nct-nnt+1
1533           call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1534          enddo
1535 #else
1536          do il=1,nodes
1537           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1538           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1539           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1540           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1541           call xdrfint(ixdrf, nss, iret)
1542           do j=1,nss
1543            call xdrfint(ixdrf, ihpb(j), iret)
1544            call xdrfint(ixdrf, jhpb(j), iret)
1545           enddo
1546           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1547           call xdrfint(ixdrf, iset_restart1(il), iret)
1548           do i=1,nfrag
1549            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1550           enddo
1551           do i=1,npair
1552            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1553           enddo
1554           do i=1,nfrag_back
1555            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1556            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1557            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1558           enddo
1559           prec=10000.0
1560           do i=1,nres
1561            do j=1,3
1562             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1563            enddo
1564           enddo
1565           do i=nnt,nct
1566            do j=1,3
1567             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1568            enddo
1569           enddo
1570           itmp=nres+nct-nnt+1
1571           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1572          enddo
1573 #endif
1574        endif
1575       enddo
1576 #ifdef AIX
1577       if(me.eq.king) call xdrfclose_(ixdrf, iret)
1578 #else
1579       if(me.eq.king) call xdrfclose(ixdrf, iret)
1580 #endif
1581       do i=1,ntwx_cache-ii_write
1582
1583             totT_cache(i)=totT_cache(ii_write+i)
1584             EK_cache(i)=EK_cache(ii_write+i)
1585             potE_cache(i)=potE_cache(ii_write+i)
1586             t_bath_cache(i)=t_bath_cache(ii_write+i)
1587             Uconst_cache(i)=Uconst_cache(ii_write+i)
1588             iset_cache(i)=iset_cache(ii_write+i)
1589
1590             do ii=1,nfrag
1591              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1592             enddo
1593             do ii=1,npair
1594              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1595             enddo
1596             do ii=1,nfrag_back
1597               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1598               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1599               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1600             enddo
1601
1602             do ii=1,nres*2
1603              do j=1,3
1604               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1605              enddo
1606             enddo
1607       enddo
1608       ntwx_cache=ntwx_cache-ii_write
1609       return
1610       end
1611
1612
1613       subroutine read1restart(i_index)
1614       implicit real*8 (a-h,o-z)
1615       include 'DIMENSIONS'
1616       include 'mpif.h'
1617       include 'COMMON.MD'
1618       include 'COMMON.IOUNITS'
1619       include 'COMMON.REMD'
1620       include 'COMMON.SETUP'
1621       include 'COMMON.CHAIN'
1622       include 'COMMON.SBRIDGE'
1623       include 'COMMON.INTERACT'
1624       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1625      &                 t5_restart1(5)
1626       integer*2 i_index
1627      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1628       common /przechowalnia/ d_restart1
1629       write (*,*) "Processor",me," called read1restart"
1630
1631          if(me.eq.king)then
1632               open(irest2,file=mremd_rst_name,status='unknown')
1633               read(irest2,*,err=334) i
1634               write(iout,*) "Reading old rst in ASCI format"
1635               close(irest2)
1636                call read1restart_old
1637                return
1638  334          continue
1639 #ifdef AIX
1640               call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1641
1642               do i=0,nodes-1
1643                call xdrfint_(ixdrf, i2rep(i), iret)
1644               enddo
1645               do i=1,remd_m(1)
1646                call xdrfint_(ixdrf, ifirst(i), iret)
1647               enddo
1648              do il=1,nodes
1649               call xdrfint_(ixdrf, nupa(0,il), iret)
1650               do i=1,nupa(0,il)
1651                call xdrfint_(ixdrf, nupa(i,il), iret)
1652               enddo
1653
1654               call xdrfint_(ixdrf, ndowna(0,il), iret)
1655               do i=1,ndowna(0,il)
1656                call xdrfint_(ixdrf, ndowna(i,il), iret)
1657               enddo
1658              enddo
1659              do il=1,nodes
1660                do j=1,4
1661                 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1662                enddo
1663              enddo
1664 #else
1665               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1666
1667               do i=0,nodes-1
1668                call xdrfint(ixdrf, i2rep(i), iret)
1669               enddo
1670               do i=1,remd_m(1)
1671                call xdrfint(ixdrf, ifirst(i), iret)
1672               enddo
1673              do il=1,nodes
1674               call xdrfint(ixdrf, nupa(0,il), iret)
1675               do i=1,nupa(0,il)
1676                call xdrfint(ixdrf, nupa(i,il), iret)
1677               enddo
1678
1679               call xdrfint(ixdrf, ndowna(0,il), iret)
1680               do i=1,ndowna(0,il)
1681                call xdrfint(ixdrf, ndowna(i,il), iret)
1682               enddo
1683              enddo
1684              do il=1,nodes
1685                do j=1,4
1686                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1687                enddo
1688              enddo
1689 #endif
1690          endif
1691          call mpi_scatter(t_restart1,5,mpi_real,
1692      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1693          totT=t5_restart1(1)
1694          EK=t5_restart1(2)
1695          potE=t5_restart1(3)
1696          t_bath=t5_restart1(4)
1697
1698          if(me.eq.king)then
1699               do il=0,nodes-1
1700                do i=1,2*nres
1701 c                read(irest2,'(3e15.5)')
1702 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1703             do j=1,3
1704 #ifdef AIX
1705              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1706 #else
1707              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1708 #endif
1709             enddo
1710                enddo
1711               enddo
1712          endif
1713          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1714      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1715
1716          do i=1,2*nres
1717            do j=1,3
1718             d_t(j,i)=r_d(j,i)
1719            enddo
1720          enddo
1721          if(me.eq.king)then
1722               do il=0,nodes-1
1723                do i=1,2*nres
1724 c                read(irest2,'(3e15.5)')
1725 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1726             do j=1,3
1727 #ifdef AIX
1728              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1729 #else
1730              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1731 #endif
1732             enddo
1733                enddo
1734               enddo
1735          endif
1736          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1737      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1738          do i=1,2*nres
1739            do j=1,3
1740             dc(j,i)=r_d(j,i)
1741            enddo
1742          enddo
1743
1744
1745            if(usampl) then
1746 #ifdef AIX
1747              if(me.eq.king)then
1748               call xdrfint_(ixdrf, nset, iret)
1749               do i=1,nset
1750                 call xdrfint_(ixdrf,mset(i), iret)
1751               enddo
1752               do i=0,nodes-1
1753                 call xdrfint_(ixdrf,i2set(i), iret)
1754               enddo
1755               do il=1,nset
1756                do il1=1,mset(il)
1757                 do i=1,nrep
1758                  do j=1,remd_m(i)
1759                    call xdrfint_(ixdrf,itmp, iret)
1760                    i_index(i,j,il,il1)=itmp
1761                  enddo
1762                 enddo
1763                enddo
1764               enddo
1765              endif
1766 #else
1767              if(me.eq.king)then
1768               call xdrfint(ixdrf, nset, iret)
1769               do i=1,nset
1770                 call xdrfint(ixdrf,mset(i), iret)
1771               enddo
1772               do i=0,nodes-1
1773                 call xdrfint(ixdrf,i2set(i), iret)
1774               enddo
1775               do il=1,nset
1776                do il1=1,mset(il)
1777                 do i=1,nrep
1778                  do j=1,remd_m(i)
1779                    call xdrfint(ixdrf,itmp, iret)
1780                    i_index(i,j,il,il1)=itmp
1781                  enddo
1782                 enddo
1783                enddo
1784               enddo
1785              endif
1786 #endif
1787               call mpi_scatter(i2set,1,mpi_integer,
1788      &           iset,1,mpi_integer,king,
1789      &           CG_COMM,ierr)
1790
1791            endif
1792
1793
1794         if(me.eq.king) close(irest2)
1795         return
1796         end
1797
1798       subroutine read1restart_old
1799       implicit real*8 (a-h,o-z)
1800       include 'DIMENSIONS'
1801       include 'mpif.h'
1802       include 'COMMON.MD'
1803       include 'COMMON.IOUNITS'
1804       include 'COMMON.REMD'
1805       include 'COMMON.SETUP'
1806       include 'COMMON.CHAIN'
1807       include 'COMMON.SBRIDGE'
1808       include 'COMMON.INTERACT'
1809       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1810      &                 t5_restart1(5)
1811       common /przechowalnia/ d_restart1
1812          if(me.eq.king)then
1813              open(irest2,file=mremd_rst_name,status='unknown')
1814              read (irest2,*) (i2rep(i),i=0,nodes-1)
1815              read (irest2,*) (ifirst(i),i=1,remd_m(1))
1816              do il=1,nodes
1817               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1818               read (irest2,*) ndowna(0,il),
1819      &                    (ndowna(i,il),i=1,ndowna(0,il))
1820              enddo
1821              do il=1,nodes
1822                read(irest2,*) (t_restart1(j,il),j=1,4)
1823              enddo
1824          endif
1825          call mpi_scatter(t_restart1,5,mpi_real,
1826      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1827          totT=t5_restart1(1)
1828          EK=t5_restart1(2)
1829          potE=t5_restart1(3)
1830          t_bath=t5_restart1(4)
1831
1832          if(me.eq.king)then
1833               do il=0,nodes-1
1834                do i=1,2*nres
1835                 read(irest2,'(3e15.5)')
1836      &                (d_restart1(j,i+2*nres*il),j=1,3)
1837                enddo
1838               enddo
1839          endif
1840          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1841      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1842
1843          do i=1,2*nres
1844            do j=1,3
1845             d_t(j,i)=r_d(j,i)
1846            enddo
1847          enddo
1848          if(me.eq.king)then
1849               do il=0,nodes-1
1850                do i=1,2*nres
1851                 read(irest2,'(3e15.5)')
1852      &                (d_restart1(j,i+2*nres*il),j=1,3)
1853                enddo
1854               enddo
1855          endif
1856          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1857      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1858          do i=1,2*nres
1859            do j=1,3
1860             dc(j,i)=r_d(j,i)
1861            enddo
1862          enddo
1863         if(me.eq.king) close(irest2)
1864         return
1865         end
1866 c------------------------------------------
1867       subroutine returnbox
1868       include 'DIMENSIONS'
1869       include 'mpif.h'
1870       include 'COMMON.CONTROL'
1871       include 'COMMON.VAR'
1872       include 'COMMON.MD'
1873 #ifndef LANG0
1874       include 'COMMON.LANGEVIN'
1875 #else
1876       include 'COMMON.LANGEVIN.lang0'
1877 #endif
1878       include 'COMMON.CHAIN'
1879       include 'COMMON.DERIV'
1880       include 'COMMON.GEO'
1881       include 'COMMON.LOCAL'
1882       include 'COMMON.INTERACT'
1883       include 'COMMON.IOUNITS'
1884       include 'COMMON.NAMES'
1885       include 'COMMON.TIME1'
1886       include 'COMMON.REMD'
1887       include 'COMMON.SETUP'
1888       include 'COMMON.MUCA'
1889       include 'COMMON.HAIRPIN'
1890         j=1
1891         chain_beg=1
1892 C        do i=1,nres
1893 C       write(*,*) 'initial', i,j,c(j,i)
1894 C        enddo
1895         do i=1,nres-1
1896          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
1897           chain_end=i
1898           if (allareout.eq.1) then
1899             ireturnval=int(c(j,i)/boxxsize)
1900             if (c(j,i).le.0) ireturnval=ireturnval-1
1901             do k=chain_beg,chain_end
1902               c(j,k)=c(j,k)-ireturnval*boxxsize
1903               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
1904             enddo
1905            endif
1906            chain_beg=i+1
1907            allareout=1
1908          else
1909           if (int(c(j,i)/boxxsize).eq.0) allareout=0
1910          endif
1911         enddo
1912          if (allareout.eq.1) then
1913             ireturnval=int(c(j,i)/boxxsize)
1914             if (c(j,i).le.0) ireturnval=ireturnval-1
1915             do k=chain_beg,nres
1916               c(j,k)=c(j,k)-ireturnval*boxxsize
1917               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
1918             enddo
1919           endif
1920 C NO JUMP
1921 C        do i=1,nres
1922 C        write(*,*) 'befor no jump', i,j,c(j,i)
1923 C        enddo
1924         nojumpval=0
1925         do i=2,nres
1926            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1927              difference=abs(c(j,i-1)-c(j,i))
1928 C             print *,'diff', difference
1929              if (difference.gt.boxxsize/2.0) then
1930                 if (c(j,i-1).gt.c(j,i)) then
1931                   nojumpval=1
1932                  else
1933                    nojumpval=-1
1934                  endif
1935               else
1936               nojumpval=0
1937               endif
1938               endif
1939               c(j,i)=c(j,i)+nojumpval*boxxsize
1940               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
1941          enddo
1942        nojumpval=0
1943         do i=2,nres
1944            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1945              difference=abs(c(j,i-1)-c(j,i))
1946              if (difference.gt.boxxsize/2.0) then
1947                 if (c(j,i-1).gt.c(j,i)) then
1948                   nojumpval=1
1949                  else
1950                    nojumpval=-1
1951                  endif
1952               else
1953               nojumpval=0
1954               endif
1955              endif
1956               c(j,i)=c(j,i)+nojumpval*boxxsize
1957               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
1958          enddo
1959
1960 C        do i=1,nres
1961 C        write(*,*) 'after no jump', i,j,c(j,i)
1962 C        enddo
1963
1964 C NOW Y dimension
1965         j=2
1966         chain_beg=1
1967         do i=1,nres-1
1968          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
1969           chain_end=i
1970           if (allareout.eq.1) then
1971             ireturnval=int(c(j,i)/boxysize)
1972             if (c(j,i).le.0) ireturnval=ireturnval-1
1973             do k=chain_beg,chain_end
1974               c(j,k)=c(j,k)-ireturnval*boxysize
1975               c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
1976             enddo
1977            endif
1978            chain_beg=i+1
1979            allareout=1
1980          else
1981           if (int(c(j,i)/boxysize).eq.0) allareout=0
1982          endif
1983         enddo
1984          if (allareout.eq.1) then
1985             ireturnval=int(c(j,i)/boxysize)
1986             if (c(j,i).le.0) ireturnval=ireturnval-1
1987             do k=chain_beg,nres
1988               c(j,k)=c(j,k)-ireturnval*boxysize
1989               c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
1990             enddo
1991           endif
1992         nojumpval=0
1993         do i=2,nres
1994            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1995              difference=abs(c(j,i-1)-c(j,i))
1996              if (difference.gt.boxysize/2.0) then
1997                 if (c(j,i-1).gt.c(j,i)) then
1998                   nojumpval=1
1999                  else
2000                    nojumpval=-1
2001                  endif
2002               else
2003               nojumpval=0
2004               endif
2005            endif
2006               c(j,i)=c(j,i)+nojumpval*boxysize
2007               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
2008          enddo
2009       nojumpval=0
2010         do i=2,nres
2011            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2012              difference=abs(c(j,i-1)-c(j,i))
2013              if (difference.gt.boxysize/2.0) then
2014                 if (c(j,i-1).gt.c(j,i)) then
2015                   nojumpval=1
2016                  else
2017                    nojumpval=-1
2018                  endif
2019               else
2020               nojumpval=0
2021               endif
2022             endif
2023               c(j,i)=c(j,i)+nojumpval*boxysize
2024               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
2025          enddo
2026
2027         j=3
2028         chain_beg=1
2029         do i=1,nres-1
2030          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
2031           chain_end=i
2032           if (allareout.eq.1) then
2033             ireturnval=int(c(j,i)/boxysize)
2034             if (c(j,i).le.0) ireturnval=ireturnval-1
2035             do k=chain_beg,chain_end
2036               c(j,k)=c(j,k)-ireturnval*boxzsize
2037               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
2038             enddo
2039            endif
2040            chain_beg=i+1
2041            allareout=1
2042          else
2043           if (int(c(j,i)/boxzsize).eq.0) allareout=0
2044          endif
2045         enddo
2046          if (allareout.eq.1) then
2047             ireturnval=int(c(j,i)/boxzsize)
2048             if (c(j,i).le.0) ireturnval=ireturnval-1
2049             do k=chain_beg,nres
2050               c(j,k)=c(j,k)-ireturnval*boxzsize
2051               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
2052             enddo
2053           endif
2054         nojumpval=0
2055         do i=2,nres
2056            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2057              difference=abs(c(j,i-1)-c(j,i))
2058              if (difference.gt.(boxzsize/2.0)) then
2059                 if (c(j,i-1).gt.c(j,i)) then
2060                   nojumpval=1
2061                  else
2062                    nojumpval=-1
2063                  endif
2064               else
2065               nojumpval=0
2066               endif
2067             endif
2068               c(j,i)=c(j,i)+nojumpval*boxzsize
2069               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
2070          enddo
2071        nojumpval=0
2072         do i=2,nres
2073            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2074              difference=abs(c(j,i-1)-c(j,i))
2075              if (difference.gt.boxzsize/2.0) then
2076                 if (c(j,i-1).gt.c(j,i)) then
2077                   nojumpval=1
2078                  else
2079                    nojumpval=-1
2080                  endif
2081               else
2082               nojumpval=0
2083               endif
2084             endif
2085               c(j,i)=c(j,i)+nojumpval*boxzsize
2086               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
2087          enddo
2088
2089         return
2090         end