source/unres/src-HCD-5D/MREMD.F

   1       subroutine MREMD
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'mpif.h'
   5       include 'COMMON.CONTROL'
   6       include 'COMMON.VAR'
   7       include 'COMMON.MD'
   8 #ifdef FIVEDIAG
   9        include 'COMMON.LAGRANGE.5diag'
  10 #else
  11        include 'COMMON.LAGRANGE'
  12 #endif
  13       include 'COMMON.QRESTR'
  14 #ifndef LANG0
  15       include 'COMMON.LANGEVIN'
  16 #else
  17 #ifdef FIVEDIAG
  18       include 'COMMON.LANGEVIN.lang0.5diag'
  19 #else
  20       include 'COMMON.LANGEVIN.lang0'
  21 #endif
  22 #endif
  23       include 'COMMON.CHAIN'
  24       include 'COMMON.DERIV'
  25       include 'COMMON.GEO'
  26       include 'COMMON.LOCAL'
  27       include 'COMMON.INTERACT'
  28       include 'COMMON.IOUNITS'
  29       include 'COMMON.NAMES'
  30       include 'COMMON.TIME1'
  31       include 'COMMON.REMD'
  32       include 'COMMON.SETUP'
  33       include 'COMMON.MUCA'
  34       include 'COMMON.HAIRPIN'
  35       double precision time00,time01,time02,time03,time04,time05,
  36      & time06,time07,time08,time001,tt0
  37       double precision scalfac
  38       integer i,j,k,il,il1,ii,iex,itmp,i_temp,i_mult,i_iset,i_mset,
  39      & i_dir,i_temp1,i_mult1,i_mset1
  40       integer ERRCODE,ierr,ierror
  41       double precision cm(3),L(3),vcm(3)
  42       double precision energia(0:n_ene)
  43       double precision remd_t_bath(maxprocs)
  44       integer iremd_iset(maxprocs)
  45       integer*2 i_index
  46      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
  47       double precision remd_ene(0:n_ene+8,maxprocs)
  48       integer iremd_acc(maxprocs),iremd_tot(maxprocs)
  49       integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
  50       integer ilen,rstcount
  51       external ilen
  52       character*50 tytul
  53       common /gucio/ cm
  54       integer itime,i_set_temp,itt,itime_master,irr,i_iset1
  55       integer nharp,iharp(4,maxres/3)
  56 cold      integer nup(0:maxprocs),ndown(0:maxprocs)
  57       integer rep2i(0:maxprocs),ireqi(maxprocs)
  58       integer icache_all(maxprocs)
  59       integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
  60       logical synflag,end_of_run,file_exist /.false./,ovrtim,first_pass
  61       double precision t_bath_temp,delta,ene_iex_iex,ene_i_i,ene_iex_i,
  62      & ene_i_iex,xxx,tmp,econstr_temp_iex,econstr_temp_i
  63       integer iran_num
  64       double precision ran_number
  65       integer i_econstr/20/
  66
  67 cdeb      imin_itime_old=0
  68       ntwx_cache=0
  69       time00=MPI_WTIME()
  70       time01=time00
  71       if(me.eq.king.or..not.out1file) then
  72        write  (iout,*) 'MREMD',nodes,'time before',time00-walltime
  73        write (iout,*) "NREP=",nrep
  74       endif
  75
  76       synflag=.false.
  77       if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
  78         call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
  79       endif
  80       mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
  81
  82 cd      print *,'MREMD',nodes
  83 cd      print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
  84 cde      write (iout,*) "Start MREMD: me",me," t_bath",t_bath
  85       k=0
  86       rep2i(k)=-1
  87       do il=1,max0(nset,1)
  88        do il1=1,max0(mset(il),1)
  89         do i=1,nrep
  90          iremd_acc(i)=0
  91          iremd_acc_usa(i)=0
  92          iremd_tot(i)=0
  93          do j=1,remd_m(i)
  94           i2rep(k)=i
  95           i2set(k)=il
  96           rep2i(i)=k
  97           k=k+1
  98           i_index(i,j,il,il1)=k
  99          enddo
 100         enddo
 101        enddo
 102       enddo
 103
 104       if(me.eq.king.or..not.out1file) then
 105        write(iout,*) (i2rep(i),i=0,nodes-1)
 106        write(iout,*) (i2set(i),i=0,nodes-1)
 107        do il=1,nset
 108         do il1=1,mset(il)
 109          do i=1,nrep
 110           do j=1,remd_m(i)
 111            write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
 112           enddo
 113          enddo
 114         enddo
 115        enddo
 116       endif
 117
 118 c      print *,'i2rep',me,i2rep(me)
 119 c      print *,'rep2i',(rep2i(i),i=0,nrep)
 120
 121 cold       if (i2rep(me).eq.nrep) then
 122 cold        nup(0)=0
 123 cold       else
 124 cold        nup(0)=remd_m(i2rep(me)+1)
 125 cold        k=rep2i(int(i2rep(me)))+1
 126 cold        do i=1,nup(0)
 127 cold         nup(i)=k
 128 cold         k=k+1
 129 cold        enddo
 130 cold       endif
 131
 132 cd       print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
 133
 134 cold       if (i2rep(me).eq.1) then
 135 cold        ndown(0)=0
 136 cold       else
 137 cold        ndown(0)=remd_m(i2rep(me)-1)
 138 cold        k=rep2i(i2rep(me)-2)+1
 139 cold        do i=1,ndown(0)
 140 cold         ndown(i)=k
 141 cold         k=k+1
 142 cold        enddo
 143 cold       endif
 144
 145 cd       print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
 146
 147
 148 c       write (*,*) "Processor",me," rest",rest,"
 149 c     &   restart1fie",restart1file
 150        if(rest.and.restart1file) then
 151            if (me.eq.king)
 152      &     inquire(file=mremd_rst_name,exist=file_exist)
 153 cd           write (*,*) me," Before broadcast: file_exist",file_exist
 154            call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
 155      &          IERR)
 156 cd           write (*,*) me," After broadcast: file_exist",file_exist
 157            if(file_exist) then
 158              if(me.eq.king.or..not.out1file)
 159      &            write  (iout,*) 'Master is reading restart1file'
 160              call read1restart(i_index)
 161            else
 162              if(me.eq.king.or..not.out1file)
 163      &            write  (iout,*) 'WARNING : no restart1file'
 164            endif
 165
 166            if(me.eq.king.or..not.out1file) then
 167               write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
 168               write(iout,*) "i_index"
 169               do il=1,nset
 170                do il1=1,mset(il)
 171                 do i=1,nrep
 172                  do j=1,remd_m(i)
 173                   write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
 174                  enddo
 175                 enddo
 176                enddo
 177               enddo
 178            endif
 179            stdfp=dsqrt(2*Rb*t_bath/d_time)
 180            do i=1,ntyp
 181              stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 182            enddo
 183            if (lang.gt.0 .and. .not.surfarea) then
 184              do i=nnt,nct-1
 185                stdforcp(i)=stdfp*dsqrt(gamp)
 186              enddo
 187              do i=nnt,nct
 188               if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
 189      &                   *dsqrt(gamsc(iabs(itype(i))))
 190               enddo
 191            endif
 192
 193        endif
 194
 195        if(me.eq.king) then
 196         if (rest.and..not.restart1file)
 197      &          inquire(file=mremd_rst_name,exist=file_exist)
 198         if(.not.file_exist.and.rest.and..not.restart1file)
 199      &       write(iout,*) 'WARNING : no restart file',mremd_rst_name
 200         IF (rest.and.file_exist.and..not.restart1file) THEN
 201              write  (iout,*) 'Master is reading restart file',
 202      &                        mremd_rst_name
 203              open(irest2,file=mremd_rst_name,status='unknown')
 204              read (irest2,*)
 205              read (irest2,*) (i2rep(i),i=0,nodes-1)
 206              read (irest2,*)
 207              read (irest2,*) (ifirst(i),i=1,remd_m(1))
 208              do il=1,nodes
 209               read (irest2,*)
 210               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 211               read (irest2,*)
 212               read (irest2,*) ndowna(0,il),
 213      &                    (ndowna(i,il),i=1,ndowna(0,il))
 214              enddo
 215              if(usampl) then
 216               read (irest2,*)
 217               read (irest2,*) nset
 218               read (irest2,*)
 219               read (irest2,*) (mset(i),i=1,nset)
 220               read (irest2,*)
 221               read (irest2,*) (i2set(i),i=0,nodes-1)
 222               read (irest2,*)
 223               do il=1,nset
 224                do il1=1,mset(il)
 225                 do i=1,nrep
 226                   read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
 227                 enddo
 228                enddo
 229               enddo
 230
 231               write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
 232               write(iout,*) "i_index"
 233               do il=1,nset
 234                do il1=1,mset(il)
 235                 do i=1,nrep
 236                  do j=1,remd_m(i)
 237                   write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
 238                  enddo
 239                 enddo
 240                enddo
 241               enddo
 242              endif
 243
 244              close(irest2)
 245
 246              write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
 247              write (iout,'(a6,1000i5)') "ifirst",
 248      &                    (ifirst(i),i=1,remd_m(1))
 249              do il=1,nodes
 250               write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
 251      &                    (nupa(i,il),i=1,nupa(0,il))
 252               write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
 253      &                    (ndowna(i,il),i=1,ndowna(0,il))
 254              enddo
 255
 256              stdfp=dsqrt(2*Rb*t_bath/d_time)
 257              do i=1,ntyp
 258                stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 259              enddo
 260              if (lang.gt.0 .and. .not.surfarea) then
 261              do i=nnt,nct-1
 262                stdforcp(i)=stdfp*dsqrt(gamp)
 263              enddo
 264              do i=nnt,nct
 265               if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
 266      &                   *dsqrt(gamsc(iabs(itype(i))))
 267              enddo
 268          endif
 269         ELSE IF (.not.(rest.and.file_exist)) THEN
 270          do il=1,remd_m(1)
 271           ifirst(il)=il
 272          enddo
 273
 274          do il=1,nodes
 275           if (i2rep(il-1).eq.nrep) then
 276            nupa(0,il)=0
 277           else
 278            nupa(0,il)=remd_m(i2rep(il-1)+1)
 279            k=rep2i(int(i2rep(il-1)))+1
 280            do i=1,nupa(0,il)
 281             nupa(i,il)=k+1
 282             k=k+1
 283            enddo
 284           endif
 285           if (i2rep(il-1).eq.1) then
 286            ndowna(0,il)=0
 287           else
 288            ndowna(0,il)=remd_m(i2rep(il-1)-1)
 289            k=rep2i(i2rep(il-1)-2)+1
 290            do i=1,ndowna(0,il)
 291             ndowna(i,il)=k+1
 292             k=k+1
 293            enddo
 294           endif
 295          enddo
 296
 297         write (iout,'(a6,100i4)') "ifirst",
 298      &                    (ifirst(i),i=1,remd_m(1))
 299         do il=1,nodes
 300          write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
 301      &                    (nupa(i,il),i=1,nupa(0,il))
 302          write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
 303      &                    (ndowna(i,il),i=1,ndowna(0,il))
 304         enddo
 305
 306         ENDIF
 307        endif
 308 c
 309 c      t_bath=retmin+(retmax-retmin)*me/(nodes-1)
 310        if(.not.(rest.and.file_exist.and.restart1file)) then
 311          if (me .eq. king) then
 312             t_bath=retmin
 313          else
 314             t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
 315          endif
 316 cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
 317          if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
 318
 319        endif
 320        if(usampl) then
 321           iset=i2set(me)
 322           if(me.eq.king.or..not.out1file)
 323      &     write(iout,*) me,"iset=",iset,"t_bath=",t_bath
 324 c          call flush(iout)
 325        endif
 326 c
 327        stdfp=dsqrt(2*Rb*t_bath/d_time)
 328        do i=1,ntyp
 329           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 330        enddo
 331
 332 c      print *,'irep',me,t_bath
 333        if (.not.rest) then
 334         if (me.eq.king .or. .not. out1file)
 335      &   write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
 336         call rescale_weights(t_bath)
 337        endif
 338
 339
 340 c------copy MD--------------
 341 c  The driver for molecular dynamics subroutines
 342 c------------------------------------------------
 343       t_MDsetup=0.0d0
 344       t_langsetup=0.0d0
 345       t_MD=0.0d0
 346       t_enegrad=0.0d0
 347       t_sdsetup=0.0d0
 348       if(me.eq.king.or..not.out1file)
 349      & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
 350 #ifdef MPI
 351       tt0 = MPI_Wtime()
 352 #else
 353       tt0 = tcpu()
 354 #endif
 355 c Determine the inverse of the inertia matrix.
 356       call setup_MD_matrices
 357 c Initialize MD
 358       call init_MD
 359       if (rest) then
 360        if (me.eq.king .or. .not. out1file)
 361      &  write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
 362        stdfp=dsqrt(2*Rb*t_bath/d_time)
 363        do i=1,ntyp
 364           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 365        enddo
 366        call rescale_weights(t_bath)
 367       endif
 368
 369 #ifdef MPI
 370       t_MDsetup = MPI_Wtime()-tt0
 371 #else
 372       t_MDsetup = tcpu()-tt0
 373 #endif
 374       rstcount=0
 375 c   Entering the MD loop
 376 #ifdef MPI
 377       tt0 = MPI_Wtime()
 378 #else
 379       tt0 = tcpu()
 380 #endif
 381       if (lang.eq.2 .or. lang.eq.3) then
 382 #ifndef   LANG0
 383         call setup_fricmat
 384         if (lang.eq.2) then
 385           call sd_verlet_p_setup
 386         else
 387           call sd_verlet_ciccotti_setup
 388         endif
 389         do i=1,dimen
 390           do j=1,dimen
 391             pfric0_mat(i,j,0)=pfric_mat(i,j)
 392             afric0_mat(i,j,0)=afric_mat(i,j)
 393             vfric0_mat(i,j,0)=vfric_mat(i,j)
 394             prand0_mat(i,j,0)=prand_mat(i,j)
 395             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 396             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 397           enddo
 398         enddo
 399         flag_stoch(0)=.true.
 400         do i=1,maxflag_stoch
 401           flag_stoch(i)=.false.
 402         enddo
 403 #else
 404         write (iout,*)
 405      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 406 #ifdef MPI
 407         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 408 #endif
 409         stop
 410 #endif
 411       else if (lang.eq.1 .or. lang.eq.4) then
 412         call setup_fricmat
 413       endif
 414       time00=MPI_WTIME()
 415       if (me.eq.king .or. .not. out1file)
 416      & write(iout,*) 'Setup time',time00-walltime
 417       call flush(iout)
 418 #ifdef MPI
 419       t_langsetup=MPI_Wtime()-tt0
 420       tt0=MPI_Wtime()
 421 #else
 422       t_langsetup=tcpu()-tt0
 423       tt0=tcpu()
 424 #endif
 425       itime=0
 426       end_of_run=.false.
 427       first_pass=.not.rest
 428 c      write (iout,*) "first_pass",first_pass
 429       do while(.not.end_of_run)
 430         itime=itime+1
 431         if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
 432         if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
 433         rstcount=rstcount+1
 434         if (lang.gt.0 .and. surfarea .and.
 435      &      mod(itime,reset_fricmat).eq.0) then
 436           if (lang.eq.2 .or. lang.eq.3) then
 437 #ifndef   LANG0
 438             call setup_fricmat
 439             if (lang.eq.2) then
 440               call sd_verlet_p_setup
 441             else
 442               call sd_verlet_ciccotti_setup
 443             endif
 444             do i=1,dimen
 445               do j=1,dimen
 446                 pfric0_mat(i,j,0)=pfric_mat(i,j)
 447                 afric0_mat(i,j,0)=afric_mat(i,j)
 448                 vfric0_mat(i,j,0)=vfric_mat(i,j)
 449                 prand0_mat(i,j,0)=prand_mat(i,j)
 450                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 451                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 452               enddo
 453             enddo
 454             flag_stoch(0)=.true.
 455             do i=1,maxflag_stoch
 456               flag_stoch(i)=.false.
 457             enddo
 458 #endif
 459           else if (lang.eq.1 .or. lang.eq.4) then
 460             call setup_fricmat
 461           endif
 462           write (iout,'(a,i10)')
 463      &      "Friction matrix reset based on surface area, itime",itime
 464         endif
 465         if (reset_vel .and. tbf .and. lang.eq.0
 466      &      .and. mod(itime,count_reset_vel).eq.0) then
 467           call random_vel
 468           if (me.eq.king .or. .not. out1file)
 469      &     write(iout,'(a,f20.2)')
 470      &     "Velocities reset to random values, time",totT
 471           do i=0,2*nres
 472             do j=1,3
 473               d_t_old(j,i)=d_t(j,i)
 474             enddo
 475           enddo
 476         endif
 477         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 478           call inertia_tensor
 479           call vcm_vel(vcm)
 480           do j=1,3
 481              d_t(j,0)=d_t(j,0)-vcm(j)
 482           enddo
 483           call kinetic(EK)
 484           kinetic_T=2.0d0/(dimen3*Rb)*EK
 485           scalfac=dsqrt(T_bath/kinetic_T)
 486 cd          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
 487           do i=0,2*nres
 488             do j=1,3
 489               d_t_old(j,i)=scalfac*d_t(j,i)
 490             enddo
 491           enddo
 492         endif
 493         if (lang.ne.4) then
 494           if (RESPA) then
 495 c Time-reversible RESPA algorithm
 496 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 497             call RESPA_step(itime)
 498           else
 499 c Variable time step algorithm.
 500             call velverlet_step(itime)
 501           endif
 502         else
 503 #ifdef BROWN
 504           call brown_step(itime)
 505 #else
 506           print *,"Brown dynamics not here!"
 507 #ifdef MPI
 508           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 509 #endif
 510           stop
 511 #endif
 512         endif
 513         if(ntwe.ne.0) then
 514           if (mod(itime,ntwe).eq.0) call statout(itime)
 515         endif
 516         if (mod(itime,ntwx).eq.0.and..not.traj1file) then
 517           write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
 518           if(mdpdb) then
 519              call hairpin(.true.,nharp,iharp)
 520              call secondary2(.true.)
 521              call pdbout(potE,tytul,ipdb)
 522           else
 523              call cartout(totT)
 524           endif
 525         endif
 526         if (mod(itime,ntwx).eq.0.and.traj1file) then
 527           if(ntwx_cache.lt.max_cache_traj_use) then
 528             ntwx_cache=ntwx_cache+1
 529           else
 530            if (max_cache_traj_use.ne.1)
 531      &      print *,itime,"processor ",me," over cache ",ntwx_cache
 532            do i=1,ntwx_cache-1
 533
 534             totT_cache(i)=totT_cache(i+1)
 535             EK_cache(i)=EK_cache(i+1)
 536             potE_cache(i)=potE_cache(i+1)
 537             t_bath_cache(i)=t_bath_cache(i+1)
 538             Uconst_cache(i)=Uconst_cache(i+1)
 539             iset_cache(i)=iset_cache(i+1)
 540
 541             do ii=1,nfrag
 542              qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
 543             enddo
 544             do ii=1,npair
 545              qpair_cache(ii,i)=qpair_cache(ii,i+1)
 546             enddo
 547             do ii=1,nfrag_back
 548               utheta_cache(ii,i)=utheta_cache(ii,i+1)
 549               ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
 550               uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
 551             enddo
 552
 553
 554             do ii=1,nres*2
 555              do j=1,3
 556               c_cache(j,ii,i)=c_cache(j,ii,i+1)
 557              enddo
 558             enddo
 559            enddo
 560           endif
 561
 562             totT_cache(ntwx_cache)=totT
 563             EK_cache(ntwx_cache)=EK
 564             potE_cache(ntwx_cache)=potE
 565             t_bath_cache(ntwx_cache)=t_bath
 566             Uconst_cache(ntwx_cache)=Uconst
 567             iset_cache(ntwx_cache)=iset
 568
 569             do i=1,nfrag
 570              qfrag_cache(i,ntwx_cache)=qfrag(i)
 571             enddo
 572             do i=1,npair
 573              qpair_cache(i,ntwx_cache)=qpair(i)
 574             enddo
 575             do i=1,nfrag_back
 576               utheta_cache(i,ntwx_cache)=utheta(i)
 577               ugamma_cache(i,ntwx_cache)=ugamma(i)
 578               uscdiff_cache(i,ntwx_cache)=uscdiff(i)
 579             enddo
 580 C            print *,'przed returnbox'
 581             call returnbox
 582 C            call enerprint(remd_ene(0,i))
 583             do i=1,nres*2
 584              do j=1,3
 585               c_cache(j,i,ntwx_cache)=c(j,i)
 586              enddo
 587             enddo
 588
 589         endif
 590         if ((rstcount.eq.1000.or.itime.eq.n_timestep)
 591      &                         .and..not.restart1file) then
 592
 593            if(me.eq.king) then
 594              open(irest1,file=mremd_rst_name,status='unknown')
 595              write (irest1,*) "i2rep"
 596              write (irest1,*) (i2rep(i),i=0,nodes-1)
 597              write (irest1,*) "ifirst"
 598              write (irest1,*) (ifirst(i),i=1,remd_m(1))
 599              do il=1,nodes
 600               write (irest1,*) "nupa",il
 601               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 602               write (irest1,*) "ndowna",il
 603               write (irest1,*) ndowna(0,il),
 604      &                   (ndowna(i,il),i=1,ndowna(0,il))
 605              enddo
 606              if(usampl) then
 607               write (irest1,*) "nset"
 608               write (irest1,*) nset
 609               write (irest1,*) "mset"
 610               write (irest1,*) (mset(i),i=1,nset)
 611               write (irest1,*) "i2set"
 612               write (irest1,*) (i2set(i),i=0,nodes-1)
 613               write (irest1,*) "i_index"
 614               do il=1,nset
 615                do il1=1,mset(il)
 616                 do i=1,nrep
 617                   write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
 618                 enddo
 619                enddo
 620               enddo
 621
 622              endif
 623              close(irest1)
 624            endif
 625            open(irest2,file=rest2name,status='unknown')
 626            write(irest2,*) totT,EK,potE,totE,t_bath
 627            do i=1,2*nres
 628             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 629            enddo
 630            do i=1,2*nres
 631             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
 632            enddo
 633            if(usampl) then
 634              write (irest2,*) iset
 635            endif
 636           close(irest2)
 637           rstcount=0
 638         endif
 639
 640 c REMD - exchange
 641 c forced synchronization
 642         if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
 643      &                                .and. .not. mremdsync) then
 644             synflag=.false.
 645             call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
 646             if (synflag) then
 647                call mpi_recv(itime_master, 1, MPI_INTEGER,
 648      &                             0,101,CG_COMM, status, ierr)
 649                call mpi_barrier(CG_COMM, ierr)
 650 cdeb               if (out1file.or.traj1file) then
 651 cdeb                call mpi_gather(itime,1,mpi_integer,
 652 cdeb     &             icache_all,1,mpi_integer,king,
 653 cdeb     &             CG_COMM,ierr)
 654                if(traj1file)
 655      &          call mpi_gather(ntwx_cache,1,mpi_integer,
 656      &             icache_all,1,mpi_integer,king,
 657      &             CG_COMM,ierr)
 658                if (.not.out1file)
 659      &               write(iout,*) 'REMD synchro at',itime_master,itime
 660                if (itime_master.ge.n_timestep .or. ovrtim())
 661      &            end_of_run=.true.
 662 ctime               call flush(iout)
 663             endif
 664         endif
 665
 666 c REMD - exchange
 667         if ((mod(itime,nstex).eq.0.and.me.eq.king
 668      &                  .or.end_of_run.and.me.eq.king )
 669      &       .and. .not. mremdsync ) then
 670            synflag=.true.
 671            do i=1,nodes-1
 672               call mpi_isend(itime,1,MPI_INTEGER,i,101,
 673      &                                CG_COMM, ireqi(i), ierr)
 674 cd            write(iout,*) 'REMD synchro with',i
 675 cd            call flush(iout)
 676            enddo
 677            call mpi_waitall(nodes-1,ireqi,statusi,ierr)
 678            call mpi_barrier(CG_COMM, ierr)
 679            time01=MPI_WTIME()
 680            write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
 681            if (out1file.or.traj1file) then
 682 cdeb            call mpi_gather(itime,1,mpi_integer,
 683 cdeb     &             itime_all,1,mpi_integer,king,
 684 cdeb     &             CG_COMM,ierr)
 685 cdeb            write(iout,'(a19,8000i8)') ' REMD synchro itime',
 686 cdeb     &                    (itime_all(i),i=1,nodes)
 687             if(traj1file) then
 688 cdeb             imin_itime=itime_all(1)
 689 cdeb             do i=2,nodes
 690 cdeb               if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
 691 cdeb             enddo
 692 cdeb             ii_write=(imin_itime-imin_itime_old)/ntwx
 693 cdeb             imin_itime_old=int(imin_itime/ntwx)*ntwx
 694 cdeb             write(iout,*) imin_itime,imin_itime_old,ii_write
 695              call mpi_gather(ntwx_cache,1,mpi_integer,
 696      &             icache_all,1,mpi_integer,king,
 697      &             CG_COMM,ierr)
 698 c             write(iout,'(a19,8000i8)') '     ntwx_cache',
 699 c     &                    (icache_all(i),i=1,nodes)
 700              ii_write=icache_all(1)
 701              do i=2,nodes
 702                if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
 703              enddo
 704 c             write(iout,*) "MIN ii_write=",ii_write
 705             endif
 706            endif
 707 ctime           call flush(iout)
 708         endif
 709         if(mremdsync .and. mod(itime,nstex).eq.0) then
 710            synflag=.true.
 711            if (me.eq.king .or. .not. out1file)
 712      &      write(iout,*) 'REMD synchro at',itime
 713
 714             if(traj1file) then
 715              call mpi_gather(ntwx_cache,1,mpi_integer,
 716      &             icache_all,1,mpi_integer,king,
 717      &             CG_COMM,ierr)
 718              if (me.eq.king) then
 719                write(iout,'(a19,8000i8)') '     ntwx_cache',
 720      &                    (icache_all(i),i=1,nodes)
 721                ii_write=icache_all(1)
 722                do i=2,nodes
 723                  if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
 724                enddo
 725                write(iout,*) "MIN ii_write=",ii_write
 726              endif
 727             endif
 728 c           call flush(iout)
 729         endif
 730         if (synflag) then
 731 c Update the time safety limiy
 732           if (time001-time00.gt.safety) then
 733             safety=time001-time00+600
 734             write (iout,*) "****** SAFETY increased to",safety," s"
 735           endif
 736           if (ovrtim()) end_of_run=.true.
 737         endif
 738         if(synflag.and..not.end_of_run) then
 739            time02=MPI_WTIME()
 740            synflag=.false.
 741
 742            write(iout,*) 'REMD before',me,t_bath
 743
 744 c           call mpi_gather(t_bath,1,mpi_double_precision,
 745 c     &             remd_t_bath,1,mpi_double_precision,king,
 746 c     &             CG_COMM,ierr)
 747            potEcomp(n_ene+1)=t_bath
 748            if (usampl) then
 749              potEcomp(n_ene+2)=iset
 750              if (iset.lt.nset) then
 751                i_set_temp=iset
 752                iset=iset+1
 753                call EconstrQ
 754                if (loc_qlike) then
 755                  call Econstr_back_qlike
 756                else
 757                  call Econstr_back
 758                endif
 759                potEcomp(n_ene+3)=Uconst+Uconst_back
 760                iset=i_set_temp
 761              endif
 762              if (iset.gt.1) then
 763                i_set_temp=iset
 764                iset=iset-1
 765                call EconstrQ
 766                if (loc_qlike) then
 767                  call Econstr_back_qlike
 768                else
 769                  call Econstr_back
 770                endif
 771                potEcomp(n_ene+4)=Uconst+Uconst_back
 772                iset=i_set_temp
 773              endif
 774              if (adaptive) then
 775 c 9/11/17 Adaptive US
 776                itt=i2rep(me)
 777 #ifdef DEBUG
 778                write (iout,*) "me ",me," itt",itt
 779 #endif
 780                if (itt.lt.nrep) then
 781                  t_bath_temp=t_bath
 782                  t_bath=remd_t(itt+1)
 783                  call EconstrQ
 784                  potEcomp(n_ene+5)=Uconst
 785 #ifdef DEBUG
 786                  write (iout,*) "t_bath",t_bath_temp,t_bath,
 787      &                   " Uconst",Uconst
 788 #endif
 789                  if (iset.lt.nset) then
 790                    i_set_temp=iset
 791                    iset=iset+1
 792                    call EconstrQ
 793                    potEcomp(n_ene+7)=Uconst
 794 #ifdef DEBUG
 795                    write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
 796 #endif
 797                    iset=i_set_temp
 798                  endif
 799                  t_bath=t_bath_temp
 800                endif
 801                if (itt.gt.1) then
 802                  t_bath_temp=t_bath
 803                  t_bath=remd_t(itt-1)
 804                  call EconstrQ
 805                  potEcomp(n_ene+6)=Uconst
 806 #ifdef DEBUG
 807                  write (iout,*) "t_bath",t_bath_temp,t_bath,
 808      &                   " Uconst",Uconst
 809 #endif
 810                  if (iset.gt.1) then
 811                    i_set_temp=iset
 812                    iset=iset-1
 813                    call EconstrQ
 814                    potEcomp(n_ene+8)=Uconst
 815 #ifdef DEBUG
 816                    write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
 817 #endif
 818                    iset=i_set_temp
 819                  endif
 820                  t_bath=t_bath_temp
 821                endif
 822              endif
 823            endif
 824            call mpi_gather(potEcomp(0),n_ene+9,mpi_double_precision,
 825      &             remd_ene(0,1),n_ene+9,mpi_double_precision,king,
 826      &             CG_COMM,ierr)
 827            if(lmuca) then
 828             call mpi_gather(elow,1,mpi_double_precision,
 829      &             elowi,1,mpi_double_precision,king,
 830      &             CG_COMM,ierr)
 831             call mpi_gather(ehigh,1,mpi_double_precision,
 832      &             ehighi,1,mpi_double_precision,king,
 833      &             CG_COMM,ierr)
 834            endif
 835
 836           time03=MPI_WTIME()
 837           if (me.eq.king .or. .not. out1file) then
 838             write(iout,*) 'REMD gather times=',time03-time01
 839      &                                        ,time03-time02
 840           endif
 841
 842           if (restart1file) call write1rst(i_index)
 843
 844           time04=MPI_WTIME()
 845           if (me.eq.king .or. .not. out1file) then
 846             write(iout,*) 'REMD writing rst time=',time04-time03
 847           endif
 848
 849           if (traj1file) call write1traj
 850 cd debugging
 851 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
 852 cdeb     &             icache_all,1,mpi_integer,king,
 853 cdeb     &             CG_COMM,ierr)
 854 cdeb            write(iout,'(a19,8000i8)') '  ntwx_cache after traj1file',
 855 cdeb     &                    (icache_all(i),i=1,nodes)
 856 cd end
 857
 858
 859           time05=MPI_WTIME()
 860           if (me.eq.king .or. .not. out1file) then
 861             write(iout,*) 'REMD writing traj time=',time05-time04
 862 c            call flush(iout)
 863           endif
 864
 865
 866           if (me.eq.king) then
 867             do i=1,nodes
 868                remd_t_bath(i)=remd_ene(n_ene+1,i)
 869                iremd_iset(i)=remd_ene(n_ene+2,i)
 870             enddo
 871             if (mremd_dec) then
 872             if(lmuca) then
 873 co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
 874              do i=1,nodes
 875                write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
 876      &            elowi(i),ehighi(i)
 877              enddo
 878             else
 879               write(iout,*) 'REMD exchange temp,ene'
 880               do i=1,nodes
 881                 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
 882                 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
 883               enddo
 884             endif
 885             endif
 886 c-------------------------------------
 887            IF(.not.usampl) THEN
 888             if (mremd_dec) then
 889             write (iout,*) "Enter exchnge, remd_m",remd_m(1),
 890      &        " nodes",nodes
 891             write (iout,*) "remd_m(1)",remd_m(1)
 892             endif
 893             do irr=1,remd_m(1)
 894                i=ifirst(iran_num(1,remd_m(1)))
 895
 896              do ii=1,nodes-1
 897
 898               if (mremd_dec)
 899      &          write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
 900              if(i.gt.0.and.nupa(0,i).gt.0) then
 901               iex=i
 902 c              if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
 903 c                write (iout,*)
 904 c     &  "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
 905 c                call flush(iout)
 906 c                call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
 907 c              endif
 908 c              do while (iex.eq.i)
 909 c                write (iout,*) "upper",nupa(int(nupa(0,i)),i)
 910                 iex=nupa(iran_num(1,int(nupa(0,i))),i)
 911 c              enddo
 912 c              write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
 913 #ifdef DEBUG
 914 c 8/21/17 AL : debug
 915               write (iout,*) "i",i,"iex",iex," temperatures",
 916      &            remd_t_bath(i),remd_t_bath(iex)
 917 #endif
 918               if (lmuca) then
 919                call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
 920               else
 921 c Swap temperatures between conformations i and iex with recalculating the free energies
 922 c following temperature changes.
 923                ene_iex_iex=remd_ene(0,iex)
 924                ene_i_i=remd_ene(0,i)
 925 c               write (iout,*) "i",i," ene_i_i",ene_i_i,
 926 c     &          " iex",iex," ene_iex_iex",ene_iex_iex
 927 c               write (iout,*) "rescaling weights with temperature",
 928 c     &          remd_t_bath(i)
 929 c               call flush(iout)
 930                call rescale_weights(remd_t_bath(i))
 931
 932 c               write (iout,*) "0,iex",remd_t_bath(i)
 933 c               call enerprint(remd_ene(0,iex))
 934
 935                call sum_energy(remd_ene(0,iex),.false.)
 936                ene_iex_i=remd_ene(0,iex)
 937 c               write (iout,*) "ene_iex_i",remd_ene(0,iex)
 938
 939 c               write (iout,*) "0,i",remd_t_bath(i)
 940 c               call enerprint(remd_ene(0,i))
 941
 942                call sum_energy(remd_ene(0,i),.false.)
 943 c               write (iout,*) "ene_i_i",remd_ene(0,i)
 944 c               call flush(iout)
 945 c               write (iout,*) "rescaling weights with temperature",
 946 c     &          remd_t_bath(iex)
 947 c               write (iout,*) "first_pass",first_pass
 948                if (.not.first_pass.and.
 949      &           real(ene_i_i).ne.real(remd_ene(0,i)))
 950      &         then
 951                 write (iout,*) "ERROR: inconsistent energies:",i,
 952      &            ene_i_i,remd_ene(0,i)
 953                endif
 954                call rescale_weights(remd_t_bath(iex))
 955
 956 c               write (iout,*) "0,i",remd_t_bath(iex)
 957 c               call enerprint(remd_ene(0,i))
 958
 959                call sum_energy(remd_ene(0,i),.false.)
 960 c               write (iout,*) "ene_i_iex",remd_ene(0,i)
 961 c               call flush(iout)
 962                ene_i_iex=remd_ene(0,i)
 963
 964 c               write (iout,*) "0,iex",remd_t_bath(iex)
 965 c               call enerprint(remd_ene(0,iex))
 966
 967                call sum_energy(remd_ene(0,iex),.false.)
 968                if (.not.first_pass.and.
 969      &           real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
 970                 write (iout,*) "ERROR: inconsistent energies:",iex,
 971      &            ene_iex_iex,remd_ene(0,iex)
 972                endif
 973 #ifdef DEBUG
 974                write (iout,*) "ene_iex_iex",remd_ene(0,iex)
 975                write (iout,*) "i",i," iex",iex
 976                write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
 977      &           " ene_i_iex",ene_i_iex,
 978      &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
 979                call flush(iout)
 980 #endif
 981                delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
 982      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
 983                delta=-delta
 984 c               write(iout,*) 'delta',delta
 985 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
 986 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
 987 c     &              (remd_t_bath(i)*remd_t_bath(iex))
 988               endif
 989               if (delta .gt. 50.0d0) then
 990                 delta=0.0d0
 991               else
 992 #ifdef OSF
 993                 if(isnan(delta))then
 994                   delta=0.0d0
 995                 else if (delta.lt.-50.0d0) then
 996                   delta=dexp(50.0d0)
 997                 else
 998                   delta=dexp(-delta)
 999                 endif
1000 #else
1001                 delta=dexp(-delta)
1002 #endif
1003               endif
1004               iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1005               xxx=ran_number(0.0d0,1.0d0)
1006 c              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1007 c              call flush(iout)
1008               if (delta .gt. xxx) then
1009                 tmp=remd_t_bath(i)
1010                 remd_t_bath(i)=remd_t_bath(iex)
1011                 remd_t_bath(iex)=tmp
1012                 remd_ene(0,i)=ene_i_iex
1013                 remd_ene(0,iex)=ene_iex_i
1014                 if(lmuca) then
1015                   tmp=elowi(i)
1016                   elowi(i)=elowi(iex)
1017                   elowi(iex)=tmp
1018                   tmp=ehighi(i)
1019                   ehighi(i)=ehighi(iex)
1020                   ehighi(iex)=tmp
1021                 endif
1022
1023
1024                 do k=0,nodes
1025                   itmp=nupa(k,i)
1026                   nupa(k,i)=nupa(k,iex)
1027                   nupa(k,iex)=itmp
1028                   itmp=ndowna(k,i)
1029                   ndowna(k,i)=ndowna(k,iex)
1030                   ndowna(k,iex)=itmp
1031                 enddo
1032                 do il=1,nodes
1033                  if (ifirst(il).eq.i) ifirst(il)=iex
1034                  do k=1,nupa(0,il)
1035                   if (nupa(k,il).eq.i) then
1036                      nupa(k,il)=iex
1037                   elseif (nupa(k,il).eq.iex) then
1038                      nupa(k,il)=i
1039                   endif
1040                  enddo
1041                  do k=1,ndowna(0,il)
1042                   if (ndowna(k,il).eq.i) then
1043                      ndowna(k,il)=iex
1044                   elseif (ndowna(k,il).eq.iex) then
1045                      ndowna(k,il)=i
1046                   endif
1047                  enddo
1048                 enddo
1049
1050                 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1051                 itmp=i2rep(i-1)
1052                 i2rep(i-1)=i2rep(iex-1)
1053                 i2rep(iex-1)=itmp
1054 #ifdef DEBUG
1055                 write(iout,*) 'exchange',i,iex
1056                 write (iout,'(a8,100i4)') "@ ifirst",
1057      &                    (ifirst(k),k=1,remd_m(1))
1058                 do il=1,nodes
1059                  write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1060      &                    (nupa(k,il),k=1,nupa(0,il))
1061                  write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1062      &                    (ndowna(k,il),k=1,ndowna(0,il))
1063                 enddo
1064                 call flush(iout)
1065 #endif
1066               else
1067                remd_ene(0,iex)=ene_iex_iex
1068                remd_ene(0,i)=ene_i_i
1069                i=iex
1070               endif
1071             endif
1072            enddo
1073            enddo
1074            first_pass=.false.
1075 cd           write (iout,*) "exchange completed"
1076 cd           call flush(iout)
1077         ELSE
1078           do ii=1,nodes
1079 cd            write(iout,*) "########",ii
1080
1081             i_temp=iran_num(1,nrep)
1082             i_mult=iran_num(1,remd_m(i_temp))
1083             i_iset=iran_num(1,nset)
1084             i_mset=iran_num(1,mset(i_iset))
1085             i=i_index(i_temp,i_mult,i_iset,i_mset)
1086
1087 cd            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1088
1089              i_dir=iran_num(1,3)
1090 cd            write(iout,*) "i_dir=",i_dir
1091
1092             if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1093
1094                i_temp1=i_temp+1
1095                i_mult1=iran_num(1,remd_m(i_temp1))
1096                i_iset1=i_iset
1097                i_mset1=iran_num(1,mset(i_iset1))
1098                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1099 c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned
1100 c       on failed replica exchange attempt
1101                econstr_temp_i=remd_ene(i_econstr,i)
1102                econstr_temp_iex=remd_ene(i_econstr,iex)
1103 c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
1104                if (adaptive) then
1105                  remd_ene(i_econstr,i)=remd_ene(n_ene+5,i)
1106                  remd_ene(i_econstr,iex)=remd_ene(n_ene+6,iex)
1107                endif
1108             elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1109
1110                i_temp1=i_temp
1111                i_mult1=iran_num(1,remd_m(i_temp1))
1112                i_iset1=i_iset+1
1113                i_mset1=iran_num(1,mset(i_iset1))
1114                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1115                econstr_temp_i=remd_ene(i_econstr,i)
1116                econstr_temp_iex=remd_ene(i_econstr,iex)
1117                remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1118                remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1119
1120             elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1121
1122                i_temp1=i_temp+1
1123                i_mult1=iran_num(1,remd_m(i_temp1))
1124                i_iset1=i_iset+1
1125                i_mset1=iran_num(1,mset(i_iset1))
1126                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1127                econstr_temp_i=remd_ene(i_econstr,i)
1128                econstr_temp_iex=remd_ene(i_econstr,iex)
1129                if (adaptive) then
1130                  remd_ene(i_econstr,i)=remd_ene(n_ene+7,i)
1131                  remd_ene(i_econstr,iex)=remd_ene(n_ene+8,iex)
1132                else
1133                  remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1134                  remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1135                endif
1136             else
1137                goto 444
1138             endif
1139
1140 cd            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1141 c            call flush(iout)
1142
1143 c Swap temperatures between conformations i and iex with recalculating the free energies
1144 c following temperature changes.
1145 #ifdef DEBUG
1146               write (iout,*) "i_dir",i_dir," i",i," iex",iex,
1147      &           " t_bath",remd_t_bath(i),remd_t_bath(iex)
1148 #endif
1149               ene_iex_iex=remd_ene(0,iex)
1150               ene_i_i=remd_ene(0,i)
1151 co              write (iout,*) "rescaling weights with temperature",
1152 co     &          remd_t_bath(i)
1153               call rescale_weights(remd_t_bath(i))
1154
1155               call sum_energy(remd_ene(0,iex),.false.)
1156               ene_iex_i=remd_ene(0,iex)
1157 cd              write (iout,*) "ene_iex_i",remd_ene(0,iex)
1158 c              call sum_energy(remd_ene(0,i),.false.)
1159 cd              write (iout,*) "ene_i_i",remd_ene(0,i)
1160 c              write (iout,*) "rescaling weights with temperature",
1161 c     &          remd_t_bath(iex)
1162 c              if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1163 c                write (iout,*) "ERROR: inconsistent energies:",i,
1164 c     &            ene_i_i,remd_ene(0,i)
1165 c              endif
1166               call rescale_weights(remd_t_bath(iex))
1167               call sum_energy(remd_ene(0,i),.false.)
1168 cd              write (iout,*) "ene_i_iex",remd_ene(0,i)
1169               ene_i_iex=remd_ene(0,i)
1170 c              call sum_energy(remd_ene(0,iex),.false.)
1171 c              if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1172 c                write (iout,*) "ERROR: inconsistent energies:",iex,
1173 c     &            ene_iex_iex,remd_ene(0,iex)
1174 c              endif
1175 cd              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1176 c              write (iout,*) "i",i," iex",iex
1177 cd              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1178 cd     &           " ene_i_iex",ene_i_iex,
1179 cd     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1180               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1181      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1182               delta=-delta
1183 cd              write(iout,*) 'delta',delta
1184 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
1185 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
1186 c     &              (remd_t_bath(i)*remd_t_bath(iex))
1187               if (delta .gt. 50.0d0) then
1188                 delta=0.0d0
1189               else
1190                 delta=dexp(-delta)
1191               endif
1192               if (i_dir.eq.1.or.i_dir.eq.3)
1193      &         iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1194               if (i_dir.eq.2.or.i_dir.eq.3)
1195      &          iremd_tot_usa(int(i2set(i-1)))=
1196      &                 iremd_tot_usa(int(i2set(i-1)))+1
1197               xxx=ran_number(0.0d0,1.0d0)
1198 cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1199               if (delta .gt. xxx) then
1200                 tmp=remd_t_bath(i)
1201                 remd_t_bath(i)=remd_t_bath(iex)
1202                 remd_t_bath(iex)=tmp
1203
1204                 itmp=iremd_iset(i)
1205                 iremd_iset(i)=iremd_iset(iex)
1206                 iremd_iset(iex)=itmp
1207
1208                 remd_ene(0,i)=ene_i_iex
1209                 remd_ene(0,iex)=ene_iex_i
1210
1211                 if (i_dir.eq.1.or.i_dir.eq.3)
1212      &           iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1213
1214                 itmp=i2rep(i-1)
1215                 i2rep(i-1)=i2rep(iex-1)
1216                 i2rep(iex-1)=itmp
1217
1218                 if (i_dir.eq.2.or.i_dir.eq.3)
1219      &           iremd_acc_usa(int(i2set(i-1)))=
1220      &                 iremd_acc_usa(int(i2set(i-1)))+1
1221
1222                 itmp=i2set(i-1)
1223                 i2set(i-1)=i2set(iex-1)
1224                 i2set(iex-1)=itmp
1225
1226                 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1227                 i_index(i_temp,i_mult,i_iset,i_mset)=
1228      &                i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1229                 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1230
1231               else
1232                remd_ene(0,iex)=ene_iex_iex
1233                remd_ene(0,i)=ene_i_i
1234                remd_ene(i_econstr,iex)=econstr_temp_iex
1235                remd_ene(i_econstr,i)=econstr_temp_i
1236               endif
1237
1238 cd      do il=1,nset
1239 cd       do il1=1,mset(il)
1240 cd        do i=1,nrep
1241 cd         do j=1,remd_m(i)
1242 cd          write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1243 cd         enddo
1244 cd        enddo
1245 cd       enddo
1246 cd      enddo
1247
1248  444      continue
1249
1250           enddo
1251
1252
1253         ENDIF
1254
1255 c-------------------------------------
1256              write (iout,*) "NREP",nrep
1257              do i=1,nrep
1258               if(iremd_tot(i).ne.0)
1259      &          write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1260      &           ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1261              enddo
1262
1263              if(usampl) then
1264               do i=1,nset
1265                if(iremd_tot_usa(i).ne.0)
1266      &           write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1267      &         iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1268               enddo
1269              endif
1270
1271 c             call flush(iout)
1272
1273 cd              write (iout,'(a6,100i4)') "ifirst",
1274 cd     &                    (ifirst(i),i=1,remd_m(1))
1275 cd              do il=1,nodes
1276 cd               write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1277 cd     &                    (nupa(i,il),i=1,nupa(0,il))
1278 cd               write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1279 cd     &                    (ndowna(i,il),i=1,ndowna(0,il))
1280 cd              enddo
1281             endif
1282
1283          time06=MPI_WTIME()
1284 cd         write (iout,*) "Before scatter"
1285 cd         call flush(iout)
1286          call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1287      &           t_bath,1,mpi_double_precision,king,
1288      &           CG_COMM,ierr)
1289 cd         write (iout,*) "After scatter"
1290 cd         call flush(iout)
1291          if(usampl)
1292      &    call mpi_scatter(iremd_iset,1,mpi_integer,
1293      &           iset,1,mpi_integer,king,
1294      &           CG_COMM,ierr)
1295
1296          time07=MPI_WTIME()
1297           if (me.eq.king .or. .not. out1file) then
1298             write(iout,*) 'REMD scatter time=',time07-time06
1299           endif
1300
1301          if(lmuca) then
1302            call mpi_scatter(elowi,1,mpi_double_precision,
1303      &           elow,1,mpi_double_precision,king,
1304      &           CG_COMM,ierr)
1305            call mpi_scatter(ehighi,1,mpi_double_precision,
1306      &           ehigh,1,mpi_double_precision,king,
1307      &           CG_COMM,ierr)
1308          endif
1309          call rescale_weights(t_bath)
1310 co         write (iout,*) "Processor",me,
1311 co     &    " rescaling weights with temperature",t_bath
1312
1313          stdfp=dsqrt(2*Rb*t_bath/d_time)
1314          do i=1,ntyp
1315            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1316          enddo
1317 c Compute the standard deviations of stochastic forces for Langevin dynamics
1318 c if the friction coefficients do not depend on surface area
1319          if (lang.gt.0 .and. .not.surfarea) then
1320            do i=nnt,nct-1
1321              stdforcp(i)=stdfp*dsqrt(gamp)
1322            enddo
1323            do i=nnt,nct
1324              if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1325      &                   *dsqrt(gamsc(iabs(itype(i))))
1326            enddo
1327          endif
1328 cde           write(iout,*) 'REMD after',me,t_bath
1329            time08=MPI_WTIME()
1330            if (me.eq.king .or. .not. out1file) then
1331             write(iout,*) 'REMD exchange time=',time08-time00
1332 c            call flush(iout)
1333            endif
1334         endif
1335       enddo
1336
1337       if (restart1file) then
1338           if (me.eq.king .or. .not. out1file)
1339      &      write(iout,*) 'writing restart at the end of run'
1340            call write1rst(i_index)
1341       endif
1342
1343       if (traj1file) call write1traj
1344 cd debugging
1345 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
1346 cdeb     &             icache_all,1,mpi_integer,king,
1347 cdeb     &             CG_COMM,ierr)
1348 cdeb            write(iout,'(a40,8000i8)')
1349 cdeb     &             '  ntwx_cache after traj1file at the end',
1350 cdeb     &             (icache_all(i),i=1,nodes)
1351 cd end
1352
1353
1354 #ifdef MPI
1355       t_MD=MPI_Wtime()-tt0
1356 #else
1357       t_MD=tcpu()-tt0
1358 #endif
1359       if (me.eq.king .or. .not. out1file) then
1360        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1361      &  '  Timing  ',
1362      & 'MD calculations setup:',t_MDsetup,
1363      & 'Energy & gradient evaluation:',t_enegrad,
1364      & 'Stochastic MD setup:',t_langsetup,
1365      & 'Stochastic MD step setup:',t_sdsetup,
1366      & 'MD steps:',t_MD
1367        write (iout,'(/28(1h=),a25,27(1h=))')
1368      & '  End of MD calculation  '
1369       endif
1370       return
1371       end
1372
1373 c-----------------------------------------------------------------------
1374       subroutine write1rst(i_index)
1375       implicit none
1376       include 'DIMENSIONS'
1377       include 'mpif.h'
1378       include 'COMMON.CONTROL'
1379       include 'COMMON.MD'
1380 #ifdef FIVEDIAG
1381        include 'COMMON.LAGRANGE.5diag'
1382 #else
1383        include 'COMMON.LAGRANGE'
1384 #endif
1385       include 'COMMON.QRESTR'
1386       include 'COMMON.IOUNITS'
1387       include 'COMMON.REMD'
1388       include 'COMMON.SETUP'
1389       include 'COMMON.CHAIN'
1390       include 'COMMON.SBRIDGE'
1391       include 'COMMON.INTERACT'
1392
1393       real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1),
1394      &     d_restart2(3,2*maxres*maxprocs)
1395       real t5_restart1(5)
1396       integer iret,itmp
1397       integer*2 i_index
1398      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1399        common /przechowalnia/ d_restart1,d_restart2
1400       integer i,j,il1,il,ixdrf
1401       integer ierr
1402
1403        t5_restart1(1)=totT
1404        t5_restart1(2)=EK
1405        t5_restart1(3)=potE
1406        t5_restart1(4)=t_bath
1407        t5_restart1(5)=Uconst
1408
1409        call mpi_gather(t5_restart1,5,mpi_real,
1410      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1411
1412
1413        do i=0,2*nres-1
1414          do j=1,3
1415            r_d(j,i)=d_t(j,i)
1416          enddo
1417        enddo
1418        call mpi_gather(r_d,3*2*nres,mpi_real,
1419      &           d_restart1,3*2*nres,mpi_real,king,
1420      &           CG_COMM,ierr)
1421
1422
1423        do i=0,2*nres-1
1424          do j=1,3
1425            r_d(j,i)=dc(j,i)
1426          enddo
1427        enddo
1428        call mpi_gather(r_d,3*2*nres,mpi_real,
1429      &           d_restart2,3*2*nres,mpi_real,king,
1430      &           CG_COMM,ierr)
1431
1432        if(me.eq.king) then
1433 #ifdef AIX
1434          call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1435          do i=0,nodes-1
1436           call xdrfint_(ixdrf, i2rep(i), iret)
1437          enddo
1438          do i=1,remd_m(1)
1439           call xdrfint_(ixdrf, ifirst(i), iret)
1440          enddo
1441          do il=1,nodes
1442               do i=0,nupa(0,il)
1443                call xdrfint_(ixdrf, nupa(i,il), iret)
1444               enddo
1445
1446               do i=0,ndowna(0,il)
1447                call xdrfint_(ixdrf, ndowna(i,il), iret)
1448               enddo
1449          enddo
1450
1451          do il=1,nodes
1452            do j=1,4
1453             call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1454            enddo
1455          enddo
1456
1457          do il=0,nodes-1
1458            do i=1,2*nres
1459             do j=1,3
1460              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1461             enddo
1462            enddo
1463          enddo
1464          do il=0,nodes-1
1465            do i=1,2*nres
1466             do j=1,3
1467              call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1468             enddo
1469            enddo
1470          enddo
1471
1472          if(usampl) then
1473            call xdrfint_(ixdrf, nset, iret)
1474            do i=1,nset
1475              call xdrfint_(ixdrf,mset(i), iret)
1476            enddo
1477            do i=0,nodes-1
1478              call xdrfint_(ixdrf,i2set(i), iret)
1479            enddo
1480            do il=1,nset
1481              do il1=1,mset(il)
1482                do i=1,nrep
1483                  do j=1,remd_m(i)
1484                    itmp=i_index(i,j,il,il1)
1485                    call xdrfint_(ixdrf,itmp, iret)
1486                  enddo
1487                enddo
1488              enddo
1489            enddo
1490
1491          endif
1492          call xdrfclose_(ixdrf, iret)
1493 #else
1494          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1495          do i=0,nodes-1
1496           call xdrfint(ixdrf, i2rep(i), iret)
1497          enddo
1498          do i=1,remd_m(1)
1499           call xdrfint(ixdrf, ifirst(i), iret)
1500          enddo
1501          do il=1,nodes
1502               do i=0,nupa(0,il)
1503                call xdrfint(ixdrf, nupa(i,il), iret)
1504               enddo
1505
1506               do i=0,ndowna(0,il)
1507                call xdrfint(ixdrf, ndowna(i,il), iret)
1508               enddo
1509          enddo
1510
1511          do il=1,nodes
1512            do j=1,4
1513             call xdrffloat(ixdrf, t_restart1(j,il), iret)
1514            enddo
1515          enddo
1516
1517          do il=0,nodes-1
1518            do i=1,2*nres
1519             do j=1,3
1520              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1521             enddo
1522            enddo
1523          enddo
1524          do il=0,nodes-1
1525            do i=1,2*nres
1526             do j=1,3
1527              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1528             enddo
1529            enddo
1530          enddo
1531
1532
1533              if(usampl) then
1534               call xdrfint(ixdrf, nset, iret)
1535               do i=1,nset
1536                 call xdrfint(ixdrf,mset(i), iret)
1537               enddo
1538               do i=0,nodes-1
1539                 call xdrfint(ixdrf,i2set(i), iret)
1540               enddo
1541               do il=1,nset
1542                do il1=1,mset(il)
1543                 do i=1,nrep
1544                  do j=1,remd_m(i)
1545                    itmp=i_index(i,j,il,il1)
1546                    call xdrfint(ixdrf,itmp, iret)
1547                  enddo
1548                 enddo
1549                enddo
1550               enddo
1551
1552              endif
1553          call xdrfclose(ixdrf, iret)
1554 #endif
1555        endif
1556       return
1557       end
1558
1559
1560       subroutine write1traj
1561       implicit none
1562       include 'DIMENSIONS'
1563       include 'mpif.h'
1564       include 'COMMON.MD'
1565       include 'COMMON.QRESTR'
1566       include 'COMMON.IOUNITS'
1567       include 'COMMON.REMD'
1568       include 'COMMON.SETUP'
1569       include 'COMMON.CHAIN'
1570       include 'COMMON.SBRIDGE'
1571       include 'COMMON.INTERACT'
1572
1573       real t5_restart1(5)
1574       integer iret,itmp
1575       real xcoord(3,maxres2+2),prec
1576       real r_qfrag(50),r_qpair(100)
1577       real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1578       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1579       real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1580      &     p_uscdiff(100*maxprocs)
1581       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1582       common /przechowalnia/ p_c
1583       integer ii,i,il,j,ixdrf
1584       integer ierr
1585
1586       call mpi_bcast(ii_write,1,mpi_integer,
1587      &           king,CG_COMM,ierr)
1588
1589 c debugging
1590 c      print *,'traj1file',me,ii_write,ntwx_cache
1591 c end debugging
1592
1593 #ifdef AIX
1594       if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1595 #else
1596       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1597 #endif
1598       do ii=1,ii_write
1599        t5_restart1(1)=totT_cache(ii)
1600        t5_restart1(2)=EK_cache(ii)
1601        t5_restart1(3)=potE_cache(ii)
1602        t5_restart1(4)=t_bath_cache(ii)
1603        t5_restart1(5)=Uconst_cache(ii)
1604        call mpi_gather(t5_restart1,5,mpi_real,
1605      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1606
1607        call mpi_gather(iset_cache(ii),1,mpi_integer,
1608      &      iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1609
1610           do i=1,nfrag
1611            r_qfrag(i)=qfrag_cache(i,ii)
1612           enddo
1613           do i=1,npair
1614            r_qpair(i)=qpair_cache(i,ii)
1615           enddo
1616           do i=1,nfrag_back
1617            r_utheta(i)=utheta_cache(i,ii)
1618            r_ugamma(i)=ugamma_cache(i,ii)
1619            r_uscdiff(i)=uscdiff_cache(i,ii)
1620           enddo
1621
1622         call mpi_gather(r_qfrag,nfrag,mpi_real,
1623      &           p_qfrag,nfrag,mpi_real,king,
1624      &           CG_COMM,ierr)
1625         call mpi_gather(r_qpair,npair,mpi_real,
1626      &           p_qpair,npair,mpi_real,king,
1627      &           CG_COMM,ierr)
1628         call mpi_gather(r_utheta,nfrag_back,mpi_real,
1629      &           p_utheta,nfrag_back,mpi_real,king,
1630      &           CG_COMM,ierr)
1631         call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1632      &           p_ugamma,nfrag_back,mpi_real,king,
1633      &           CG_COMM,ierr)
1634         call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1635      &           p_uscdiff,nfrag_back,mpi_real,king,
1636      &           CG_COMM,ierr)
1637
1638 #ifdef DEBUG
1639         write (iout,*) "p_qfrag"
1640         do i=1,nodes
1641           write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1642         enddo
1643         write (iout,*) "p_qpair"
1644         do i=1,nodes
1645           write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1646         enddo
1647         call flush(iout)
1648 #endif
1649         do i=1,nres*2
1650          do j=1,3
1651           r_c(j,i)=c_cache(j,i,ii)
1652          enddo
1653         enddo
1654
1655         call mpi_gather(r_c,3*2*nres,mpi_real,
1656      &           p_c,3*2*nres,mpi_real,king,
1657      &           CG_COMM,ierr)
1658
1659        if(me.eq.king) then
1660 #ifdef AIX
1661          do il=1,nodes
1662           call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1663           call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1664           call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1665           call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1666           call xdrfint_(ixdrf, nss, iret)
1667           do j=1,nss
1668            if (dyn_ss) then
1669             call xdrfint(ixdrf, idssb(j)+nres, iret)
1670             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1671            else
1672             call xdrfint_(ixdrf, ihpb(j), iret)
1673             call xdrfint_(ixdrf, jhpb(j), iret)
1674            endif
1675           enddo
1676           call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1677           call xdrfint_(ixdrf, iset_restart1(il), iret)
1678           do i=1,nfrag
1679            call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1680           enddo
1681           do i=1,npair
1682            call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1683           enddo
1684           do i=1,nfrag_back
1685            call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1686            call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1687            call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1688           enddo
1689           prec=10000.0
1690           do i=1,nres
1691            do j=1,3
1692             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1693            enddo
1694           enddo
1695           do i=nnt,nct
1696            do j=1,3
1697             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1698            enddo
1699           enddo
1700           itmp=nres+nct-nnt+1
1701           call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1702          enddo
1703 #else
1704          do il=1,nodes
1705           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1706           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1707           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1708           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1709           call xdrfint(ixdrf, nss, iret)
1710           do j=1,nss
1711            if (dyn_ss) then
1712             call xdrfint(ixdrf, idssb(j)+nres, iret)
1713             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1714            else
1715             call xdrfint(ixdrf, ihpb(j), iret)
1716             call xdrfint(ixdrf, jhpb(j), iret)
1717            endif
1718           enddo
1719           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1720           call xdrfint(ixdrf, iset_restart1(il), iret)
1721           do i=1,nfrag
1722            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1723           enddo
1724           do i=1,npair
1725            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1726           enddo
1727           do i=1,nfrag_back
1728            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1729            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1730            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1731           enddo
1732           prec=10000.0
1733           do i=1,nres
1734            do j=1,3
1735             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1736            enddo
1737           enddo
1738           do i=nnt,nct
1739            do j=1,3
1740             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1741            enddo
1742           enddo
1743           itmp=nres+nct-nnt+1
1744           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1745          enddo
1746 #endif
1747        endif
1748       enddo
1749 #ifdef AIX
1750       if(me.eq.king) call xdrfclose_(ixdrf, iret)
1751 #else
1752       if(me.eq.king) call xdrfclose(ixdrf, iret)
1753 #endif
1754       do i=1,ntwx_cache-ii_write
1755
1756             totT_cache(i)=totT_cache(ii_write+i)
1757             EK_cache(i)=EK_cache(ii_write+i)
1758             potE_cache(i)=potE_cache(ii_write+i)
1759             t_bath_cache(i)=t_bath_cache(ii_write+i)
1760             Uconst_cache(i)=Uconst_cache(ii_write+i)
1761             iset_cache(i)=iset_cache(ii_write+i)
1762
1763             do ii=1,nfrag
1764              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1765             enddo
1766             do ii=1,npair
1767              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1768             enddo
1769             do ii=1,nfrag_back
1770               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1771               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1772               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1773             enddo
1774
1775             do ii=1,nres*2
1776              do j=1,3
1777               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1778              enddo
1779             enddo
1780       enddo
1781       ntwx_cache=ntwx_cache-ii_write
1782       return
1783       end
1784
1785
1786       subroutine read1restart(i_index)
1787       implicit none
1788       include 'DIMENSIONS'
1789       include 'mpif.h'
1790       include 'COMMON.CONTROL'
1791       include 'COMMON.MD'
1792 #ifdef FIVEDIAG
1793        include 'COMMON.LAGRANGE.5diag'
1794 #else
1795        include 'COMMON.LAGRANGE'
1796 #endif
1797       include 'COMMON.QRESTR'
1798       include 'COMMON.IOUNITS'
1799       include 'COMMON.REMD'
1800       include 'COMMON.SETUP'
1801       include 'COMMON.CHAIN'
1802       include 'COMMON.SBRIDGE'
1803       include 'COMMON.INTERACT'
1804       real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1),
1805      &                 t5_restart1(5)
1806       integer*2 i_index
1807      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1808       common /przechowalnia/ d_restart1
1809       integer i,j,il,il1,ixdrf,iret,itmp
1810       integer ierr
1811 c      write (*,*) "Processor",me," called read1restart"
1812
1813          if(me.eq.king)then
1814               open(irest2,file=mremd_rst_name,status='unknown')
1815               read(irest2,*,err=334) i
1816               write(iout,*) "Reading old rst in ASCI format"
1817               close(irest2)
1818                call read1restart_old
1819                return
1820  334          continue
1821 #ifdef AIX
1822               call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1823
1824               do i=0,nodes-1
1825                call xdrfint_(ixdrf, i2rep(i), iret)
1826               enddo
1827               do i=1,remd_m(1)
1828                call xdrfint_(ixdrf, ifirst(i), iret)
1829               enddo
1830              do il=1,nodes
1831               call xdrfint_(ixdrf, nupa(0,il), iret)
1832               do i=1,nupa(0,il)
1833                call xdrfint_(ixdrf, nupa(i,il), iret)
1834               enddo
1835
1836               call xdrfint_(ixdrf, ndowna(0,il), iret)
1837               do i=1,ndowna(0,il)
1838                call xdrfint_(ixdrf, ndowna(i,il), iret)
1839               enddo
1840              enddo
1841              do il=1,nodes
1842                do j=1,4
1843                 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1844                enddo
1845              enddo
1846 #else
1847               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1848
1849               do i=0,nodes-1
1850                call xdrfint(ixdrf, i2rep(i), iret)
1851               enddo
1852               do i=1,remd_m(1)
1853                call xdrfint(ixdrf, ifirst(i), iret)
1854               enddo
1855              do il=1,nodes
1856               call xdrfint(ixdrf, nupa(0,il), iret)
1857               do i=1,nupa(0,il)
1858                call xdrfint(ixdrf, nupa(i,il), iret)
1859               enddo
1860
1861               call xdrfint(ixdrf, ndowna(0,il), iret)
1862               do i=1,ndowna(0,il)
1863                call xdrfint(ixdrf, ndowna(i,il), iret)
1864               enddo
1865              enddo
1866              do il=1,nodes
1867                do j=1,4
1868                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1869                enddo
1870              enddo
1871 #endif
1872          endif
1873          call mpi_scatter(t_restart1,5,mpi_real,
1874      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1875          totT=t5_restart1(1)
1876          EK=t5_restart1(2)
1877          potE=t5_restart1(3)
1878          t_bath=t5_restart1(4)
1879
1880          if(me.eq.king)then
1881               do il=0,nodes-1
1882                do i=1,2*nres
1883 c                read(irest2,'(3e15.5)')
1884 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1885             do j=1,3
1886 #ifdef AIX
1887              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1888 #else
1889              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1890 #endif
1891             enddo
1892                enddo
1893               enddo
1894          endif
1895          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1896      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1897
1898          do i=0,2*nres-1
1899            do j=1,3
1900             d_t(j,i)=r_d(j,i)
1901            enddo
1902          enddo
1903          if(me.eq.king)then
1904               do il=0,nodes-1
1905                do i=1,2*nres
1906 c                read(irest2,'(3e15.5)')
1907 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1908             do j=1,3
1909 #ifdef AIX
1910              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1911 #else
1912              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1913 #endif
1914             enddo
1915                enddo
1916               enddo
1917          endif
1918          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1919      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1920          do i=0,2*nres-1
1921            do j=1,3
1922             dc(j,i)=r_d(j,i)
1923            enddo
1924          enddo
1925
1926
1927            if(usampl) then
1928 #ifdef AIX
1929              if(me.eq.king)then
1930               call xdrfint_(ixdrf, nset, iret)
1931               do i=1,nset
1932                 call xdrfint_(ixdrf,mset(i), iret)
1933               enddo
1934               do i=0,nodes-1
1935                 call xdrfint_(ixdrf,i2set(i), iret)
1936               enddo
1937               do il=1,nset
1938                do il1=1,mset(il)
1939                 do i=1,nrep
1940                  do j=1,remd_m(i)
1941                    call xdrfint_(ixdrf,itmp, iret)
1942                    i_index(i,j,il,il1)=itmp
1943                  enddo
1944                 enddo
1945                enddo
1946               enddo
1947              endif
1948 #else
1949              if(me.eq.king)then
1950               call xdrfint(ixdrf, nset, iret)
1951               do i=1,nset
1952                 call xdrfint(ixdrf,mset(i), iret)
1953               enddo
1954               do i=0,nodes-1
1955                 call xdrfint(ixdrf,i2set(i), iret)
1956               enddo
1957               do il=1,nset
1958                do il1=1,mset(il)
1959                 do i=1,nrep
1960                  do j=1,remd_m(i)
1961                    call xdrfint(ixdrf,itmp, iret)
1962                    i_index(i,j,il,il1)=itmp
1963                  enddo
1964                 enddo
1965                enddo
1966               enddo
1967              endif
1968 #endif
1969 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1970 c own element
1971 c              call mpi_scatter(i2set,1,mpi_integer,
1972 c     &           iset,1,mpi_integer,king,
1973 c     &           CG_COMM,ierr)
1974               call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1975      &         CG_COMM,ierr)
1976               iset=i2set(me)
1977
1978            endif
1979
1980
1981         if(me.eq.king) close(irest2)
1982         return
1983         end
1984
1985       subroutine read1restart_old
1986       implicit none
1987       include 'DIMENSIONS'
1988       include 'mpif.h'
1989       include 'COMMON.MD'
1990 #ifdef FIVEDIAG
1991        include 'COMMON.LAGRANGE.5diag'
1992 #else
1993        include 'COMMON.LAGRANGE'
1994 #endif
1995       include 'COMMON.IOUNITS'
1996       include 'COMMON.REMD'
1997       include 'COMMON.SETUP'
1998       include 'COMMON.CHAIN'
1999       include 'COMMON.SBRIDGE'
2000       include 'COMMON.INTERACT'
2001       real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1),
2002      &                 t5_restart1(5)
2003       common /przechowalnia/ d_restart1
2004       integer i,j,il,itmp
2005       integer ierr
2006          if(me.eq.king)then
2007              open(irest2,file=mremd_rst_name,status='unknown')
2008              read (irest2,*) (i2rep(i),i=0,nodes-1)
2009              read (irest2,*) (ifirst(i),i=1,remd_m(1))
2010              do il=1,nodes
2011               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2012               read (irest2,*) ndowna(0,il),
2013      &                    (ndowna(i,il),i=1,ndowna(0,il))
2014              enddo
2015              do il=1,nodes
2016                read(irest2,*) (t_restart1(j,il),j=1,4)
2017              enddo
2018          endif
2019          call mpi_scatter(t_restart1,5,mpi_real,
2020      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2021          totT=t5_restart1(1)
2022          EK=t5_restart1(2)
2023          potE=t5_restart1(3)
2024          t_bath=t5_restart1(4)
2025
2026          if(me.eq.king)then
2027               do il=0,nodes-1
2028                do i=1,2*nres
2029                 read(irest2,'(3e15.5)')
2030      &                (d_restart1(j,i+2*nres*il),j=1,3)
2031                enddo
2032               enddo
2033          endif
2034          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2035      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2036
2037          do i=0,2*nres-1
2038            do j=1,3
2039             d_t(j,i)=r_d(j,i)
2040            enddo
2041          enddo
2042          if(me.eq.king)then
2043               do il=0,nodes-1
2044                do i=1,2*nres
2045                 read(irest2,'(3e15.5)')
2046      &                (d_restart1(j,i+2*nres*il),j=1,3)
2047                enddo
2048               enddo
2049          endif
2050          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2051      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2052          do i=0,2*nres-1
2053            do j=1,3
2054             dc(j,i)=r_d(j,i)
2055            enddo
2056          enddo
2057         if(me.eq.king) close(irest2)
2058         return
2059         end
2060 c------------------------------------------