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