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