source/unres/src_MD-M/MREMD.F

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