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
 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            if (dyn_ss) then
1504             call xdrfint(ixdrf, idssb(j)+nres, iret)
1505             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1506            else
1507             call xdrfint_(ixdrf, ihpb(j), iret)
1508             call xdrfint_(ixdrf, jhpb(j), iret)
1509            endif
1510           enddo
1511           call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1512           call xdrfint_(ixdrf, iset_restart1(il), iret)
1513           do i=1,nfrag
1514            call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1515           enddo
1516           do i=1,npair
1517            call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1518           enddo
1519           do i=1,nfrag_back
1520            call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1521            call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1522            call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1523           enddo
1524           prec=10000.0
1525           do i=1,nres
1526            do j=1,3
1527             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1528            enddo
1529           enddo
1530           do i=nnt,nct
1531            do j=1,3
1532             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1533            enddo
1534           enddo
1535           itmp=nres+nct-nnt+1
1536           call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1537          enddo
1538 #else
1539          do il=1,nodes
1540           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1541           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1542           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1543           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1544           call xdrfint(ixdrf, nss, iret)
1545           do j=1,nss
1546            if (dyn_ss) then
1547             call xdrfint(ixdrf, idssb(j)+nres, iret)
1548             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1549            else
1550             call xdrfint(ixdrf, ihpb(j), iret)
1551             call xdrfint(ixdrf, jhpb(j), iret)
1552            endif
1553           enddo
1554           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1555           call xdrfint(ixdrf, iset_restart1(il), iret)
1556           do i=1,nfrag
1557            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1558           enddo
1559           do i=1,npair
1560            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1561           enddo
1562           do i=1,nfrag_back
1563            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1564            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1565            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1566           enddo
1567           prec=10000.0
1568           do i=1,nres
1569            do j=1,3
1570             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1571            enddo
1572           enddo
1573           do i=nnt,nct
1574            do j=1,3
1575             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1576            enddo
1577           enddo
1578           itmp=nres+nct-nnt+1
1579           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1580          enddo
1581 #endif
1582        endif
1583       enddo
1584 #ifdef AIX
1585       if(me.eq.king) call xdrfclose_(ixdrf, iret)
1586 #else
1587       if(me.eq.king) call xdrfclose(ixdrf, iret)
1588 #endif
1589       do i=1,ntwx_cache-ii_write
1590
1591             totT_cache(i)=totT_cache(ii_write+i)
1592             EK_cache(i)=EK_cache(ii_write+i)
1593             potE_cache(i)=potE_cache(ii_write+i)
1594             t_bath_cache(i)=t_bath_cache(ii_write+i)
1595             Uconst_cache(i)=Uconst_cache(ii_write+i)
1596             iset_cache(i)=iset_cache(ii_write+i)
1597
1598             do ii=1,nfrag
1599              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1600             enddo
1601             do ii=1,npair
1602              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1603             enddo
1604             do ii=1,nfrag_back
1605               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1606               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1607               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1608             enddo
1609
1610             do ii=1,nres*2
1611              do j=1,3
1612               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1613              enddo
1614             enddo
1615       enddo
1616       ntwx_cache=ntwx_cache-ii_write
1617       return
1618       end
1619
1620
1621       subroutine read1restart(i_index)
1622       implicit real*8 (a-h,o-z)
1623       include 'DIMENSIONS'
1624       include 'mpif.h'
1625       include 'COMMON.MD'
1626       include 'COMMON.IOUNITS'
1627       include 'COMMON.REMD'
1628       include 'COMMON.SETUP'
1629       include 'COMMON.CHAIN'
1630       include 'COMMON.SBRIDGE'
1631       include 'COMMON.INTERACT'
1632       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1633      &                 t5_restart1(5)
1634       integer*2 i_index
1635      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1636       common /przechowalnia/ d_restart1
1637       write (*,*) "Processor",me," called read1restart"
1638
1639          if(me.eq.king)then
1640               open(irest2,file=mremd_rst_name,status='unknown')
1641               read(irest2,*,err=334) i
1642               write(iout,*) "Reading old rst in ASCI format"
1643               close(irest2)
1644                call read1restart_old
1645                return
1646  334          continue
1647 #ifdef AIX
1648               call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1649
1650               do i=0,nodes-1
1651                call xdrfint_(ixdrf, i2rep(i), iret)
1652               enddo
1653               do i=1,remd_m(1)
1654                call xdrfint_(ixdrf, ifirst(i), iret)
1655               enddo
1656              do il=1,nodes
1657               call xdrfint_(ixdrf, nupa(0,il), iret)
1658               do i=1,nupa(0,il)
1659                call xdrfint_(ixdrf, nupa(i,il), iret)
1660               enddo
1661
1662               call xdrfint_(ixdrf, ndowna(0,il), iret)
1663               do i=1,ndowna(0,il)
1664                call xdrfint_(ixdrf, ndowna(i,il), iret)
1665               enddo
1666              enddo
1667              do il=1,nodes
1668                do j=1,4
1669                 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1670                enddo
1671              enddo
1672 #else
1673               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1674
1675               do i=0,nodes-1
1676                call xdrfint(ixdrf, i2rep(i), iret)
1677               enddo
1678               do i=1,remd_m(1)
1679                call xdrfint(ixdrf, ifirst(i), iret)
1680               enddo
1681              do il=1,nodes
1682               call xdrfint(ixdrf, nupa(0,il), iret)
1683               do i=1,nupa(0,il)
1684                call xdrfint(ixdrf, nupa(i,il), iret)
1685               enddo
1686
1687               call xdrfint(ixdrf, ndowna(0,il), iret)
1688               do i=1,ndowna(0,il)
1689                call xdrfint(ixdrf, ndowna(i,il), iret)
1690               enddo
1691              enddo
1692              do il=1,nodes
1693                do j=1,4
1694                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1695                enddo
1696              enddo
1697 #endif
1698          endif
1699          call mpi_scatter(t_restart1,5,mpi_real,
1700      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1701          totT=t5_restart1(1)
1702          EK=t5_restart1(2)
1703          potE=t5_restart1(3)
1704          t_bath=t5_restart1(4)
1705
1706          if(me.eq.king)then
1707               do il=0,nodes-1
1708                do i=1,2*nres
1709 c                read(irest2,'(3e15.5)')
1710 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1711             do j=1,3
1712 #ifdef AIX
1713              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1714 #else
1715              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1716 #endif
1717             enddo
1718                enddo
1719               enddo
1720          endif
1721          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1722      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1723
1724          do i=1,2*nres
1725            do j=1,3
1726             d_t(j,i)=r_d(j,i)
1727            enddo
1728          enddo
1729          if(me.eq.king)then
1730               do il=0,nodes-1
1731                do i=1,2*nres
1732 c                read(irest2,'(3e15.5)')
1733 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1734             do j=1,3
1735 #ifdef AIX
1736              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1737 #else
1738              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1739 #endif
1740             enddo
1741                enddo
1742               enddo
1743          endif
1744          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1745      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1746          do i=1,2*nres
1747            do j=1,3
1748             dc(j,i)=r_d(j,i)
1749            enddo
1750          enddo
1751
1752
1753            if(usampl) then
1754 #ifdef AIX
1755              if(me.eq.king)then
1756               call xdrfint_(ixdrf, nset, iret)
1757               do i=1,nset
1758                 call xdrfint_(ixdrf,mset(i), iret)
1759               enddo
1760               do i=0,nodes-1
1761                 call xdrfint_(ixdrf,i2set(i), iret)
1762               enddo
1763               do il=1,nset
1764                do il1=1,mset(il)
1765                 do i=1,nrep
1766                  do j=1,remd_m(i)
1767                    call xdrfint_(ixdrf,itmp, iret)
1768                    i_index(i,j,il,il1)=itmp
1769                  enddo
1770                 enddo
1771                enddo
1772               enddo
1773              endif
1774 #else
1775              if(me.eq.king)then
1776               call xdrfint(ixdrf, nset, iret)
1777               do i=1,nset
1778                 call xdrfint(ixdrf,mset(i), iret)
1779               enddo
1780               do i=0,nodes-1
1781                 call xdrfint(ixdrf,i2set(i), iret)
1782               enddo
1783               do il=1,nset
1784                do il1=1,mset(il)
1785                 do i=1,nrep
1786                  do j=1,remd_m(i)
1787                    call xdrfint(ixdrf,itmp, iret)
1788                    i_index(i,j,il,il1)=itmp
1789                  enddo
1790                 enddo
1791                enddo
1792               enddo
1793              endif
1794 #endif
1795 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1796 c own element
1797 c              call mpi_scatter(i2set,1,mpi_integer,
1798 c     &           iset,1,mpi_integer,king,
1799 c     &           CG_COMM,ierr)
1800               call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1801      &         CG_COMM,ierr)
1802               iset=i2set(me)
1803            endif
1804
1805
1806         if(me.eq.king) close(irest2)
1807         return
1808         end
1809
1810       subroutine read1restart_old
1811       implicit real*8 (a-h,o-z)
1812       include 'DIMENSIONS'
1813       include 'mpif.h'
1814       include 'COMMON.MD'
1815       include 'COMMON.IOUNITS'
1816       include 'COMMON.REMD'
1817       include 'COMMON.SETUP'
1818       include 'COMMON.CHAIN'
1819       include 'COMMON.SBRIDGE'
1820       include 'COMMON.INTERACT'
1821       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1822      &                 t5_restart1(5)
1823       common /przechowalnia/ d_restart1
1824          if(me.eq.king)then
1825              open(irest2,file=mremd_rst_name,status='unknown')
1826              read (irest2,*) (i2rep(i),i=0,nodes-1)
1827              read (irest2,*) (ifirst(i),i=1,remd_m(1))
1828              do il=1,nodes
1829               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1830               read (irest2,*) ndowna(0,il),
1831      &                    (ndowna(i,il),i=1,ndowna(0,il))
1832              enddo
1833              do il=1,nodes
1834                read(irest2,*) (t_restart1(j,il),j=1,4)
1835              enddo
1836          endif
1837          call mpi_scatter(t_restart1,5,mpi_real,
1838      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1839          totT=t5_restart1(1)
1840          EK=t5_restart1(2)
1841          potE=t5_restart1(3)
1842          t_bath=t5_restart1(4)
1843
1844          if(me.eq.king)then
1845               do il=0,nodes-1
1846                do i=1,2*nres
1847                 read(irest2,'(3e15.5)')
1848      &                (d_restart1(j,i+2*nres*il),j=1,3)
1849                enddo
1850               enddo
1851          endif
1852          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1853      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1854
1855          do i=1,2*nres
1856            do j=1,3
1857             d_t(j,i)=r_d(j,i)
1858            enddo
1859          enddo
1860          if(me.eq.king)then
1861               do il=0,nodes-1
1862                do i=1,2*nres
1863                 read(irest2,'(3e15.5)')
1864      &                (d_restart1(j,i+2*nres*il),j=1,3)
1865                enddo
1866               enddo
1867          endif
1868          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1869      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1870          do i=1,2*nres
1871            do j=1,3
1872             dc(j,i)=r_d(j,i)
1873            enddo
1874          enddo
1875         if(me.eq.king) close(irest2)
1876         return
1877         end
1878