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