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