source/unres/src_MD-M/MREMD.F

   1 c#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() .and. me.eq.king) end_of_run=.true.
 704           call MPI_Bcast(end_of_run,1,MPI_LOGICAL,king,CG_COMM,IERR)
 705         endif
 706         if(synflag.and..not.end_of_run) then
 707            time02=MPI_WTIME()
 708            synflag=.false.
 709
 710 c           write(iout,*) 'REMD before',me,t_bath
 711
 712 c           call mpi_gather(t_bath,1,mpi_double_precision,
 713 c     &             remd_t_bath,1,mpi_double_precision,king,
 714 c     &             CG_COMM,ierr)
 715            potEcomp(n_ene+1)=t_bath
 716            t_bath_old=t_bath
 717            if (usampl.or.homol_nset.gt.1) then
 718              potEcomp(n_ene+2)=iset
 719              if (iset.lt.nset) then
 720                i_set_temp=iset
 721                iset=iset+1
 722                if (homol_nset.gt.1) then
 723 c broadcast iset to slaves and reduce energy
 724                 if (nfgtasks.gt.1) then
 725                  call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
 726                  call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
 727                  call e_modeller(e_tmp)
 728 c                 write(iout,*) "iset+1 before reduce",e_tmp
 729                  call MPI_Barrier(FG_COMM,IERR)
 730                  call MPI_Reduce(e_tmp,potEcomp(n_ene+3),1,
 731      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
 732                 else
 733                  call e_modeller(potEcomp(n_ene+3))
 734                 endif
 735 c                write(iout,*) "iset+1",potEcomp(n_ene+3)
 736                else
 737                 call EconstrQ
 738                 potEcomp(n_ene+3)=Uconst
 739                endif
 740                iset=i_set_temp
 741 c broadcast iset to slaves
 742                if (nfgtasks.gt.1) then
 743                  call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
 744                  call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
 745                endif
 746              else
 747               potEcomp(n_ene+3)=0.0
 748              endif
 749              if (iset.gt.1) then
 750                i_set_temp=iset
 751                iset=iset-1
 752                if (homol_nset.gt.1) then
 753 c broadcast iset to slaves and reduce energy
 754                 if (nfgtasks.gt.1) then
 755                  call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
 756                  call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
 757                  call e_modeller(e_tmp)
 758 c                 write(iout,*) "iset-1 before reduce",e_tmp
 759                  call MPI_Barrier(FG_COMM,IERR)
 760                  call MPI_Reduce(e_tmp,potEcomp(n_ene+4),1,
 761      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
 762                 else
 763                  call e_modeller(potEcomp(n_ene+4))
 764                 endif
 765 c                write(iout,*) "iset-1",potEcomp(n_ene+4)
 766                else
 767                 call EconstrQ
 768                 potEcomp(n_ene+4)=Uconst
 769                endif
 770                iset=i_set_temp
 771 c broadcast iset to slaves
 772                if (nfgtasks.gt.1) then
 773                  call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
 774                  call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
 775                endif
 776              else
 777                potEcomp(n_ene+4)=0.0
 778              endif
 779            endif
 780            call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
 781      &             remd_ene(0,1),n_ene+5,mpi_double_precision,king,
 782      &             CG_COMM,ierr)
 783            if(lmuca) then
 784             call mpi_gather(elow,1,mpi_double_precision,
 785      &             elowi,1,mpi_double_precision,king,
 786      &             CG_COMM,ierr)
 787             call mpi_gather(ehigh,1,mpi_double_precision,
 788      &             ehighi,1,mpi_double_precision,king,
 789      &             CG_COMM,ierr)
 790            endif
 791
 792           time03=MPI_WTIME()
 793           if (me.eq.king .or. .not. out1file) then
 794             write(iout,*) 'REMD gather times=',time03-time01
 795      &                                        ,time03-time02
 796           endif
 797
 798           if (restart1file) call write1rst(i_index)
 799
 800           time04=MPI_WTIME()
 801           if (me.eq.king .or. .not. out1file) then
 802             write(iout,*) 'REMD writing rst time=',time04-time03
 803           endif
 804
 805           if (traj1file) call write1traj
 806 cd debugging
 807 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
 808 cdeb     &             icache_all,1,mpi_integer,king,
 809 cdeb     &             CG_COMM,ierr)
 810 cdeb            write(iout,'(a19,8000i8)') '  ntwx_cache after traj1file',
 811 cdeb     &                    (icache_all(i),i=1,nodes)
 812 cd end
 813
 814
 815           time05=MPI_WTIME()
 816           if (me.eq.king .or. .not. out1file) then
 817             write(iout,*) 'REMD writing traj time=',time05-time04
 818             call flush(iout)
 819           endif
 820
 821
 822           if (me.eq.king) then
 823            if(homol_nset.gt.1) write(iout,*)
 824      &     'energy_c temperature iset energy_c(iset+1) energy_c(iset-1)'
 825             do i=1,nodes
 826                remd_t_bath(i)=remd_ene(n_ene+1,i)
 827                iremd_iset(i)=remd_ene(n_ene+2,i)
 828                if(homol_nset.gt.1)
 829      &                write(iout,'(i4,f10.3,f6.0,i3,2f10.3)')
 830      &                i,remd_ene(i_econstr,i),
 831      &                remd_ene(n_ene+1,i),iremd_iset(i),
 832      &                remd_ene(n_ene+3,i),remd_ene(n_ene+4,i)
 833             enddo
 834 #ifdef DEBUG
 835             if(lmuca) then
 836 co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
 837              do i=1,nodes
 838                write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
 839      &            elowi(i),ehighi(i)
 840              enddo
 841             else
 842               write(iout,*) 'REMD exchange temp,ene'
 843               do i=1,nodes
 844                 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
 845                 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
 846               enddo
 847             endif
 848 #endif
 849 c-------------------------------------
 850            IF(.not.usampl.and.hremd.eq.0.and.homol_nset.le.1) THEN
 851 #ifdef DEBUG
 852             write (iout,*) "Enter exchnge, remd_m",remd_m(1),
 853      &        " nodes",nodes
 854 ctime            call flush(iout)
 855             write (iout,*) "remd_m(1)",remd_m(1)
 856 #endif
 857             do irr=1,remd_m(1)
 858                i=ifirst(iran_num(1,remd_m(1)))
 859 #ifdef DEBUG
 860              write (iout,*) "i",i
 861 #endif
 862 ctime             call flush(iout)
 863
 864              do ii=1,nodes-1
 865
 866 #ifdef DEBUG
 867               write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
 868 #endif
 869              if(i.gt.0.and.nupa(0,i).gt.0) then
 870               iex=i
 871 c              if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
 872 c                write (iout,*)
 873 c     &  "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
 874 c                call flush(iout)
 875 c                call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
 876 c              endif
 877 c              do while (iex.eq.i)
 878 c                write (iout,*) "upper",nupa(int(nupa(0,i)),i)
 879                 iex=nupa(iran_num(1,int(nupa(0,i))),i)
 880 c              enddo
 881 c              write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
 882               if (lmuca) then
 883                call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
 884               else
 885 c Swap temperatures between conformations i and iex with recalculating the free energies
 886 c following temperature changes.
 887                ene_iex_iex=remd_ene(0,iex)
 888                ene_i_i=remd_ene(0,i)
 889 c               write (iout,*) "i",i," ene_i_i",ene_i_i,
 890 c     &          " iex",iex," ene_iex_iex",ene_iex_iex
 891 c               write (iout,*) "rescaling weights with temperature",
 892 c     &          remd_t_bath(i)
 893 c               call flush(iout)
 894                call rescale_weights(remd_t_bath(i))
 895
 896 c               write (iout,*) "0,iex",remd_t_bath(i)
 897 c               call enerprint(remd_ene(0,iex))
 898
 899                call sum_energy(remd_ene(0,iex),.false.)
 900                ene_iex_i=remd_ene(0,iex)
 901 c               write (iout,*) "ene_iex_i",remd_ene(0,iex)
 902
 903 c               write (iout,*) "0,i",remd_t_bath(i)
 904 c               call enerprint(remd_ene(0,i))
 905
 906                call sum_energy(remd_ene(0,i),.false.)
 907 c               write (iout,*) "ene_i_i",remd_ene(0,i)
 908 c               call flush(iout)
 909 c               write (iout,*) "rescaling weights with temperature",
 910 c     &          remd_t_bath(iex)
 911                if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
 912                 write (iout,*) "ERROR: inconsistent energies:",i,
 913      &            ene_i_i,remd_ene(0,i)
 914                endif
 915                call rescale_weights(remd_t_bath(iex))
 916
 917 c               write (iout,*) "0,i",remd_t_bath(iex)
 918 c               call enerprint(remd_ene(0,i))
 919
 920                call sum_energy(remd_ene(0,i),.false.)
 921 c               write (iout,*) "ene_i_iex",remd_ene(0,i)
 922 c               call flush(iout)
 923                ene_i_iex=remd_ene(0,i)
 924
 925 c               write (iout,*) "0,iex",remd_t_bath(iex)
 926 c               call enerprint(remd_ene(0,iex))
 927
 928                call sum_energy(remd_ene(0,iex),.false.)
 929                if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
 930                 write (iout,*) "ERROR: inconsistent energies:",iex,
 931      &            ene_iex_iex,remd_ene(0,iex)
 932                endif
 933 c               write (iout,*) "ene_iex_iex",remd_ene(0,iex)
 934 c               write (iout,*) "i",i," iex",iex
 935 c               write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
 936 c     &           " ene_i_iex",ene_i_iex,
 937 c     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
 938 c               call flush(iout)
 939                delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
 940      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
 941                delta=-delta
 942 c               write(iout,*) 'delta',delta
 943 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
 944 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
 945 c     &              (remd_t_bath(i)*remd_t_bath(iex))
 946               endif
 947               if (delta .gt. 50.0d0) then
 948                 delta=0.0d0
 949               else
 950 #ifdef OSF
 951                 if(isnan(delta))then
 952                   delta=0.0d0
 953                 else if (delta.lt.-50.0d0) then
 954                   delta=dexp(50.0d0)
 955                 else
 956                   delta=dexp(-delta)
 957                 endif
 958 #else
 959                 delta=dexp(-delta)
 960 #endif
 961               endif
 962               iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
 963               xxx=ran_number(0.0d0,1.0d0)
 964 c              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
 965 c              call flush(iout)
 966               if (delta .gt. xxx) then
 967                 tmp=remd_t_bath(i)
 968                 remd_t_bath(i)=remd_t_bath(iex)
 969                 remd_t_bath(iex)=tmp
 970                 remd_ene(0,i)=ene_i_iex
 971                 remd_ene(0,iex)=ene_iex_i
 972                 if(lmuca) then
 973                   tmp=elowi(i)
 974                   elowi(i)=elowi(iex)
 975                   elowi(iex)=tmp
 976                   tmp=ehighi(i)
 977                   ehighi(i)=ehighi(iex)
 978                   ehighi(iex)=tmp
 979                 endif
 980
 981
 982                 do k=0,nodes
 983                   itmp=nupa(k,i)
 984                   nupa(k,i)=nupa(k,iex)
 985                   nupa(k,iex)=itmp
 986                   itmp=ndowna(k,i)
 987                   ndowna(k,i)=ndowna(k,iex)
 988                   ndowna(k,iex)=itmp
 989                 enddo
 990                 do il=1,nodes
 991                  if (ifirst(il).eq.i) ifirst(il)=iex
 992                  do k=1,nupa(0,il)
 993                   if (nupa(k,il).eq.i) then
 994                      nupa(k,il)=iex
 995                   elseif (nupa(k,il).eq.iex) then
 996                      nupa(k,il)=i
 997                   endif
 998                  enddo
 999                  do k=1,ndowna(0,il)
1000                   if (ndowna(k,il).eq.i) then
1001                      ndowna(k,il)=iex
1002                   elseif (ndowna(k,il).eq.iex) then
1003                      ndowna(k,il)=i
1004                   endif
1005                  enddo
1006                 enddo
1007
1008                 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1009                 itmp=i2rep(i-1)
1010                 i2rep(i-1)=i2rep(iex-1)
1011                 i2rep(iex-1)=itmp
1012
1013 c                write(iout,*) 'exchange',i,iex
1014 c                write (iout,'(a8,100i4)') "@ ifirst",
1015 c     &                    (ifirst(k),k=1,remd_m(1))
1016 c                do il=1,nodes
1017 c                 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1018 c     &                    (nupa(k,il),k=1,nupa(0,il))
1019 c                 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1020 c     &                    (ndowna(k,il),k=1,ndowna(0,il))
1021 c                enddo
1022 c                call flush(iout)
1023
1024               else
1025                remd_ene(0,iex)=ene_iex_iex
1026                remd_ene(0,i)=ene_i_i
1027                i=iex
1028               endif
1029             endif
1030            enddo
1031            enddo
1032 cd           write (iout,*) "exchange completed"
1033 cd           call flush(iout)
1034         ELSEIF (usampl.or.homol_nset.gt.1) THEN
1035           do ii=1,nodes
1036 c            write(iout,*) "########",ii
1037
1038             i_temp=iran_num(1,nrep)
1039             i_mult=iran_num(1,remd_m(i_temp))
1040             i_iset=iran_num(1,nset)
1041             i_mset=iran_num(1,mset(i_iset))
1042             i=i_index(i_temp,i_mult,i_iset,i_mset)
1043
1044 c            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1045
1046             if(t_exchange_only)then
1047              i_dir=1
1048             else
1049              i_dir=iran_num(1,3)
1050             endif
1051 c            write(iout,*) "i_dir=",i_dir
1052
1053             if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1054
1055                i_temp1=i_temp+1
1056                i_mult1=iran_num(1,remd_m(i_temp1))
1057                i_iset1=i_iset
1058                i_mset1=iran_num(1,mset(i_iset1))
1059                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1060
1061             elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1062
1063                i_temp1=i_temp
1064                i_mult1=iran_num(1,remd_m(i_temp1))
1065                i_iset1=i_iset+1
1066                i_mset1=iran_num(1,mset(i_iset1))
1067                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1068
1069                econstr_temp_i=remd_ene(i_econstr,i)
1070                econstr_temp_iex=remd_ene(i_econstr,iex)
1071                remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1072                remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1073
1074             elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1075
1076                i_temp1=i_temp+1
1077                i_mult1=iran_num(1,remd_m(i_temp1))
1078                i_iset1=i_iset+1
1079                i_mset1=iran_num(1,mset(i_iset1))
1080                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1081                econstr_temp_i=remd_ene(i_econstr,i)
1082                econstr_temp_iex=remd_ene(i_econstr,iex)
1083                remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1084                remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1085
1086             else
1087                goto 444
1088             endif
1089
1090 c            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1091 ctime            call flush(iout)
1092
1093 c Swap temperatures between conformations i and iex with recalculating the free energies
1094 c following temperature changes.
1095               ene_iex_iex=remd_ene(0,iex)
1096               ene_i_i=remd_ene(0,i)
1097 co              write (iout,*) "rescaling weights with temperature",
1098 co     &          remd_t_bath(i)
1099               call rescale_weights(remd_t_bath(i))
1100
1101               call sum_energy(remd_ene(0,iex),.false.)
1102               ene_iex_i=remd_ene(0,iex)
1103 cdebug
1104 c ERROR only makes sense for dir =1
1105 c              write (iout,*) "ene_iex_i",remd_ene(0,iex)
1106 c              call sum_energy(remd_ene(0,i),.false.)
1107 c              write (iout,*) "ene_i_i",remd_ene(0,i)
1108 c              write (iout,*) "rescaling weights with temperature",
1109 c     &          remd_t_bath(iex)
1110 c              if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1111 c                write (iout,*) "ERROR: inconsistent energies i:",i,
1112 c     &            ene_i_i,remd_ene(0,i)
1113 c              endif
1114 cdebug_end
1115               call rescale_weights(remd_t_bath(iex))
1116               call sum_energy(remd_ene(0,i),.false.)
1117 cd              write (iout,*) "ene_i_iex",remd_ene(0,i)
1118               ene_i_iex=remd_ene(0,i)
1119 cdebug
1120 c ERROR only makes sense for dir =1
1121 c              call sum_energy(remd_ene(0,iex),.false.)
1122 c              if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1123 c                write (iout,*) "ERROR: inconsistent energies iex:",iex,
1124 c     &            ene_iex_iex,remd_ene(0,iex)
1125 c              endif
1126 c              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1127 c              write (iout,*) "i",i," iex",iex
1128 c              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1129 c     &           " ene_i_iex",ene_i_iex,
1130 c     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1131 cdebug_end
1132               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1133      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1134               delta=-delta
1135 c              write(iout,*) 'delta',delta
1136 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
1137 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
1138 c     &              (remd_t_bath(i)*remd_t_bath(iex))
1139               if (delta .gt. 50.0d0) then
1140                 delta=0.0d0
1141               else
1142                 delta=dexp(-delta)
1143               endif
1144               if (i_dir.eq.1.or.i_dir.eq.3)
1145      &         iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1146               if (i_dir.eq.2.or.i_dir.eq.3)
1147      &          iremd_tot_usa(int(i2set(i-1)))=
1148      &                 iremd_tot_usa(int(i2set(i-1)))+1
1149               xxx=ran_number(0.0d0,1.0d0)
1150 c              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1151               if (delta .gt. xxx) then
1152                 tmp=remd_t_bath(i)
1153                 remd_t_bath(i)=remd_t_bath(iex)
1154                 remd_t_bath(iex)=tmp
1155
1156                 itmp=iremd_iset(i)
1157                 iremd_iset(i)=iremd_iset(iex)
1158                 iremd_iset(iex)=itmp
1159
1160                 remd_ene(0,i)=ene_i_iex
1161                 remd_ene(0,iex)=ene_iex_i
1162
1163                 if (i_dir.eq.1.or.i_dir.eq.3)
1164      &           iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1165
1166                 itmp=i2rep(i-1)
1167                 i2rep(i-1)=i2rep(iex-1)
1168                 i2rep(iex-1)=itmp
1169
1170                 if (i_dir.eq.2.or.i_dir.eq.3)
1171      &           iremd_acc_usa(int(i2set(i-1)))=
1172      &                 iremd_acc_usa(int(i2set(i-1)))+1
1173
1174                 itmp=i2set(i-1)
1175                 i2set(i-1)=i2set(iex-1)
1176                 i2set(iex-1)=itmp
1177
1178                 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1179                 i_index(i_temp,i_mult,i_iset,i_mset)=
1180      &                i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1181                 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1182
1183               else
1184                remd_ene(0,iex)=ene_iex_iex
1185                remd_ene(0,i)=ene_i_i
1186                remd_ene(i_econstr,iex)=econstr_temp_iex
1187                remd_ene(i_econstr,i)=econstr_temp_i
1188               endif
1189
1190 cd      do il=1,nset
1191 cd       do il1=1,mset(il)
1192 cd        do i=1,nrep
1193 cd         do j=1,remd_m(i)
1194 cd          write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1195 cd         enddo
1196 cd        enddo
1197 cd       enddo
1198 cd      enddo
1199
1200  444      continue
1201
1202           enddo
1203
1204
1205         ENDIF
1206
1207 c-------------------------------------
1208              write (iout,*) "NREP",nrep
1209              do i=1,nrep
1210               if(iremd_tot(i).ne.0)
1211      &          write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1212      &           ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1213              enddo
1214
1215              if(usampl.or.homol_nset.gt.1) then
1216               do i=1,nset
1217                if(iremd_tot_usa(i).ne.0)
1218      &           write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1219      &         iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1220               enddo
1221              endif
1222
1223              call flush(iout)
1224
1225 cd              write (iout,'(a6,100i4)') "ifirst",
1226 cd     &                    (ifirst(i),i=1,remd_m(1))
1227 cd              do il=1,nodes
1228 cd               write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1229 cd     &                    (nupa(i,il),i=1,nupa(0,il))
1230 cd               write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1231 cd     &                    (ndowna(i,il),i=1,ndowna(0,il))
1232 cd              enddo
1233             endif
1234
1235          time06=MPI_WTIME()
1236 cd         write (iout,*) "Before scatter"
1237 cd         call flush(iout)
1238          call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1239      &           t_bath,1,mpi_double_precision,king,
1240      &           CG_COMM,ierr)
1241 cd         write (iout,*) "After scatter"
1242 cd         call flush(iout)
1243          if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
1244           call mpi_scatter(iremd_iset,1,mpi_integer,
1245      &           iset,1,mpi_integer,king,
1246      &           CG_COMM,ierr)
1247 c 8/31/2015 Correction by AL: send new iset to slaves
1248           if (nfgtasks.gt.1) then
1249            call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
1250            call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
1251           endif
1252
1253          endif
1254
1255          time07=MPI_WTIME()
1256           if (me.eq.king .or. .not. out1file) then
1257             write(iout,*) 'REMD scatter time=',time07-time06
1258           endif
1259
1260          if(lmuca) then
1261            call mpi_scatter(elowi,1,mpi_double_precision,
1262      &           elow,1,mpi_double_precision,king,
1263      &           CG_COMM,ierr)
1264            call mpi_scatter(ehighi,1,mpi_double_precision,
1265      &           ehigh,1,mpi_double_precision,king,
1266      &           CG_COMM,ierr)
1267          endif
1268          call rescale_weights(t_bath)
1269 co         write (iout,*) "Processor",me,
1270 co     &    " rescaling weights with temperature",t_bath
1271
1272          stdfp=dsqrt(2*Rb*t_bath/d_time)
1273          do i=1,ntyp
1274            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1275          enddo
1276
1277 cde         write(iout,*) 'REMD after',me,t_bath
1278            time08=MPI_WTIME()
1279            if (me.eq.king .or. .not. out1file) then
1280             write(iout,*) 'REMD exchange time=',time08-time02
1281 ctime            call flush(iout)
1282            endif
1283         endif
1284       enddo
1285
1286       if (restart1file) then
1287           if (me.eq.king .or. .not. out1file)
1288      &      write(iout,*) 'writing restart at the end of run'
1289            call write1rst(i_index)
1290       endif
1291
1292       if (traj1file) call write1traj
1293 cd debugging
1294 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
1295 cdeb     &             icache_all,1,mpi_integer,king,
1296 cdeb     &             CG_COMM,ierr)
1297 cdeb            write(iout,'(a40,8000i8)')
1298 cdeb     &             '  ntwx_cache after traj1file at the end',
1299 cdeb     &             (icache_all(i),i=1,nodes)
1300 cd end
1301
1302
1303 #ifdef MPI
1304       t_MD=MPI_Wtime()-tt0
1305 #else
1306       t_MD=tcpu()-tt0
1307 #endif
1308       if (me.eq.king .or. .not. out1file) then
1309        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1310      &  '  Timing  ',
1311      & 'MD calculations setup:',t_MDsetup,
1312      & 'Energy & gradient evaluation:',t_enegrad,
1313      & 'Stochastic MD setup:',t_langsetup,
1314      & 'Stochastic MD step setup:',t_sdsetup,
1315      & 'MD steps:',t_MD
1316        write (iout,'(/28(1h=),a25,27(1h=))')
1317      & '  End of MD calculation  '
1318       endif
1319       return
1320       end
1321
1322 c-----------------------------------------------------------------------
1323       subroutine write1rst(i_index)
1324       implicit real*8 (a-h,o-z)
1325       include 'DIMENSIONS'
1326       include 'mpif.h'
1327       include 'COMMON.MD'
1328       include 'COMMON.IOUNITS'
1329       include 'COMMON.REMD'
1330       include 'COMMON.SETUP'
1331       include 'COMMON.CHAIN'
1332       include 'COMMON.SBRIDGE'
1333       include 'COMMON.INTERACT'
1334       include 'COMMON.CONTROL'
1335
1336       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1337      &     d_restart2(3,2*maxres*maxprocs)
1338       real t5_restart1(5)
1339       integer iret,itmp
1340       integer*2 i_index
1341      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1342        common /przechowalnia/ d_restart1,d_restart2
1343
1344        t5_restart1(1)=totT
1345        t5_restart1(2)=EK
1346        t5_restart1(3)=potE
1347        t5_restart1(4)=t_bath
1348        t5_restart1(5)=Uconst
1349
1350        call mpi_gather(t5_restart1,5,mpi_real,
1351      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1352
1353
1354        do i=1,2*nres
1355          do j=1,3
1356            r_d(j,i)=d_t(j,i)
1357          enddo
1358        enddo
1359        call mpi_gather(r_d,3*2*nres,mpi_real,
1360      &           d_restart1,3*2*nres,mpi_real,king,
1361      &           CG_COMM,ierr)
1362
1363
1364        do i=1,2*nres
1365          do j=1,3
1366            r_d(j,i)=dc(j,i)
1367          enddo
1368        enddo
1369        call mpi_gather(r_d,3*2*nres,mpi_real,
1370      &           d_restart2,3*2*nres,mpi_real,king,
1371      &           CG_COMM,ierr)
1372
1373        if(me.eq.king) then
1374 #ifdef AIX
1375          call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1376          do i=0,nodes-1
1377           call xdrfint_(ixdrf, i2rep(i), iret)
1378          enddo
1379          do i=1,remd_m(1)
1380           call xdrfint_(ixdrf, ifirst(i), iret)
1381          enddo
1382          do il=1,nodes
1383               do i=0,nupa(0,il)
1384                call xdrfint_(ixdrf, nupa(i,il), iret)
1385               enddo
1386
1387               do i=0,ndowna(0,il)
1388                call xdrfint_(ixdrf, ndowna(i,il), iret)
1389               enddo
1390          enddo
1391
1392          do il=1,nodes
1393            do j=1,4
1394             call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1395            enddo
1396          enddo
1397
1398          do il=0,nodes-1
1399            do i=1,2*nres
1400             do j=1,3
1401              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1402             enddo
1403            enddo
1404          enddo
1405          do il=0,nodes-1
1406            do i=1,2*nres
1407             do j=1,3
1408              call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1409             enddo
1410            enddo
1411          enddo
1412
1413          if(usampl.or.homol_nset.gt.1) then
1414            call xdrfint_(ixdrf, nset, iret)
1415            do i=1,nset
1416              call xdrfint_(ixdrf,mset(i), iret)
1417            enddo
1418            do i=0,nodes-1
1419              call xdrfint_(ixdrf,i2set(i), iret)
1420            enddo
1421            do il=1,nset
1422              do il1=1,mset(il)
1423                do i=1,nrep
1424                  do j=1,remd_m(i)
1425                    itmp=i_index(i,j,il,il1)
1426                    call xdrfint_(ixdrf,itmp, iret)
1427                  enddo
1428                enddo
1429              enddo
1430            enddo
1431
1432          endif
1433          call xdrfclose_(ixdrf, iret)
1434 #else
1435          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1436          do i=0,nodes-1
1437           call xdrfint(ixdrf, i2rep(i), iret)
1438          enddo
1439          do i=1,remd_m(1)
1440           call xdrfint(ixdrf, ifirst(i), iret)
1441          enddo
1442          do il=1,nodes
1443               do i=0,nupa(0,il)
1444                call xdrfint(ixdrf, nupa(i,il), iret)
1445               enddo
1446
1447               do i=0,ndowna(0,il)
1448                call xdrfint(ixdrf, ndowna(i,il), iret)
1449               enddo
1450          enddo
1451
1452          do il=1,nodes
1453            do j=1,4
1454             call xdrffloat(ixdrf, t_restart1(j,il), iret)
1455            enddo
1456          enddo
1457
1458          do il=0,nodes-1
1459            do i=1,2*nres
1460             do j=1,3
1461              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1462             enddo
1463            enddo
1464          enddo
1465          do il=0,nodes-1
1466            do i=1,2*nres
1467             do j=1,3
1468              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1469             enddo
1470            enddo
1471          enddo
1472
1473
1474              if(usampl.or.homol_nset.gt.1) then
1475               call xdrfint(ixdrf, nset, iret)
1476               do i=1,nset
1477                 call xdrfint(ixdrf,mset(i), iret)
1478               enddo
1479               do i=0,nodes-1
1480                 call xdrfint(ixdrf,i2set(i), iret)
1481               enddo
1482               do il=1,nset
1483                do il1=1,mset(il)
1484                 do i=1,nrep
1485                  do j=1,remd_m(i)
1486                    itmp=i_index(i,j,il,il1)
1487                    call xdrfint(ixdrf,itmp, iret)
1488                  enddo
1489                 enddo
1490                enddo
1491               enddo
1492
1493              endif
1494          call xdrfclose(ixdrf, iret)
1495 #endif
1496        endif
1497       return
1498       end
1499
1500
1501       subroutine write1traj
1502       implicit real*8 (a-h,o-z)
1503       include 'DIMENSIONS'
1504       include 'mpif.h'
1505       include 'COMMON.MD'
1506       include 'COMMON.IOUNITS'
1507       include 'COMMON.REMD'
1508       include 'COMMON.SETUP'
1509       include 'COMMON.CHAIN'
1510       include 'COMMON.SBRIDGE'
1511       include 'COMMON.INTERACT'
1512
1513       real t5_restart1(5)
1514       integer iret,itmp
1515       real xcoord(3,maxres2+2),prec
1516       real r_qfrag(50),r_qpair(100)
1517       real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1518       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1519       real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1520      &     p_uscdiff(100*maxprocs)
1521       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1522       common /przechowalnia/ p_c
1523
1524       call mpi_bcast(ii_write,1,mpi_integer,
1525      &           king,CG_COMM,ierr)
1526
1527 c debugging
1528       print *,'traj1file',me,ii_write,ntwx_cache
1529 c end debugging
1530
1531 #ifdef AIX
1532       if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1533 #else
1534       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1535 #endif
1536       do ii=1,ii_write
1537        t5_restart1(1)=totT_cache(ii)
1538        t5_restart1(2)=EK_cache(ii)
1539        t5_restart1(3)=potE_cache(ii)
1540        t5_restart1(4)=t_bath_cache(ii)
1541        t5_restart1(5)=Uconst_cache(ii)
1542        call mpi_gather(t5_restart1,5,mpi_real,
1543      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1544
1545        call mpi_gather(iset_cache(ii),1,mpi_integer,
1546      &      iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1547
1548           do i=1,nfrag
1549            r_qfrag(i)=qfrag_cache(i,ii)
1550           enddo
1551           do i=1,npair
1552            r_qpair(i)=qpair_cache(i,ii)
1553           enddo
1554           do i=1,nfrag_back
1555            r_utheta(i)=utheta_cache(i,ii)
1556            r_ugamma(i)=ugamma_cache(i,ii)
1557            r_uscdiff(i)=uscdiff_cache(i,ii)
1558           enddo
1559
1560         call mpi_gather(r_qfrag,nfrag,mpi_real,
1561      &           p_qfrag,nfrag,mpi_real,king,
1562      &           CG_COMM,ierr)
1563         call mpi_gather(r_qpair,npair,mpi_real,
1564      &           p_qpair,npair,mpi_real,king,
1565      &           CG_COMM,ierr)
1566         call mpi_gather(r_utheta,nfrag_back,mpi_real,
1567      &           p_utheta,nfrag_back,mpi_real,king,
1568      &           CG_COMM,ierr)
1569         call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1570      &           p_ugamma,nfrag_back,mpi_real,king,
1571      &           CG_COMM,ierr)
1572         call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1573      &           p_uscdiff,nfrag_back,mpi_real,king,
1574      &           CG_COMM,ierr)
1575
1576 #ifdef DEBUG
1577         write (iout,*) "p_qfrag"
1578         do i=1,nodes
1579           write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1580         enddo
1581         write (iout,*) "p_qpair"
1582         do i=1,nodes
1583           write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1584         enddo
1585         call flush(iout)
1586 #endif
1587         do i=1,nres*2
1588          do j=1,3
1589           r_c(j,i)=c_cache(j,i,ii)
1590          enddo
1591         enddo
1592
1593         call mpi_gather(r_c,3*2*nres,mpi_real,
1594      &           p_c,3*2*nres,mpi_real,king,
1595      &           CG_COMM,ierr)
1596
1597        if(me.eq.king) then
1598 #ifdef AIX
1599          do il=1,nodes
1600           call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1601           call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1602           call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1603           call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1604           call xdrfint_(ixdrf, nss, iret)
1605           do j=1,nss
1606            if (dyn_ss) then
1607             call xdrfint_(ixdrf, idssb(j)+nres, iret)
1608             call xdrfint_(ixdrf, jdssb(j)+nres, iret)
1609            else
1610             call xdrfint_(ixdrf, ihpb(j), iret)
1611             call xdrfint_(ixdrf, jhpb(j), iret)
1612            endif
1613           enddo
1614           call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1615           call xdrfint_(ixdrf, iset_restart1(il), iret)
1616           do i=1,nfrag
1617            call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1618           enddo
1619           do i=1,npair
1620            call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1621           enddo
1622           do i=1,nfrag_back
1623            call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1624            call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1625            call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1626           enddo
1627           prec=10000.0
1628           do i=1,nres
1629            do j=1,3
1630             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1631            enddo
1632           enddo
1633           do i=nnt,nct
1634            do j=1,3
1635             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1636            enddo
1637           enddo
1638           itmp=nres+nct-nnt+1
1639           call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1640          enddo
1641 #else
1642          do il=1,nodes
1643           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1644           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1645           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1646           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1647           call xdrfint(ixdrf, nss, iret)
1648           do j=1,nss
1649            if (dyn_ss) then
1650             call xdrfint(ixdrf, idssb(j)+nres, iret)
1651             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1652            else
1653             call xdrfint(ixdrf, ihpb(j), iret)
1654             call xdrfint(ixdrf, jhpb(j), iret)
1655            endif
1656           enddo
1657           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1658           call xdrfint(ixdrf, iset_restart1(il), iret)
1659           do i=1,nfrag
1660            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1661           enddo
1662           do i=1,npair
1663            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1664           enddo
1665           do i=1,nfrag_back
1666            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1667            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1668            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1669           enddo
1670           prec=10000.0
1671           do i=1,nres
1672            do j=1,3
1673             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1674            enddo
1675           enddo
1676           do i=nnt,nct
1677            do j=1,3
1678             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1679            enddo
1680           enddo
1681           itmp=nres+nct-nnt+1
1682           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1683          enddo
1684 #endif
1685        endif
1686       enddo
1687 #ifdef AIX
1688       if(me.eq.king) call xdrfclose_(ixdrf, iret)
1689 #else
1690       if(me.eq.king) call xdrfclose(ixdrf, iret)
1691 #endif
1692       do i=1,ntwx_cache-ii_write
1693
1694             totT_cache(i)=totT_cache(ii_write+i)
1695             EK_cache(i)=EK_cache(ii_write+i)
1696             potE_cache(i)=potE_cache(ii_write+i)
1697             t_bath_cache(i)=t_bath_cache(ii_write+i)
1698             Uconst_cache(i)=Uconst_cache(ii_write+i)
1699             iset_cache(i)=iset_cache(ii_write+i)
1700
1701             do ii=1,nfrag
1702              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1703             enddo
1704             do ii=1,npair
1705              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1706             enddo
1707             do ii=1,nfrag_back
1708               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1709               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1710               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1711             enddo
1712
1713             do ii=1,nres*2
1714              do j=1,3
1715               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1716              enddo
1717             enddo
1718       enddo
1719       ntwx_cache=ntwx_cache-ii_write
1720       return
1721       end
1722
1723
1724       subroutine read1restart(i_index)
1725       implicit real*8 (a-h,o-z)
1726       include 'DIMENSIONS'
1727       include 'mpif.h'
1728       include 'COMMON.MD'
1729       include 'COMMON.IOUNITS'
1730       include 'COMMON.REMD'
1731       include 'COMMON.SETUP'
1732       include 'COMMON.CHAIN'
1733       include 'COMMON.SBRIDGE'
1734       include 'COMMON.INTERACT'
1735       include 'COMMON.CONTROL'
1736       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1737      &                 t5_restart1(5)
1738       integer*2 i_index
1739      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1740       common /przechowalnia/ d_restart1
1741       write (*,*) "Processor",me," called read1restart"
1742
1743          if(me.eq.king)then
1744               open(irest2,file=mremd_rst_name,status='unknown')
1745               read(irest2,*,err=334) i
1746               write(iout,*) "Reading old rst in ASCI format"
1747               close(irest2)
1748                call read1restart_old
1749                return
1750  334          continue
1751 #ifdef AIX
1752               call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1753
1754               do i=0,nodes-1
1755                call xdrfint_(ixdrf, i2rep(i), iret)
1756               enddo
1757               do i=1,remd_m(1)
1758                call xdrfint_(ixdrf, ifirst(i), iret)
1759               enddo
1760              do il=1,nodes
1761               call xdrfint_(ixdrf, nupa(0,il), iret)
1762               do i=1,nupa(0,il)
1763                call xdrfint_(ixdrf, nupa(i,il), iret)
1764               enddo
1765
1766               call xdrfint_(ixdrf, ndowna(0,il), iret)
1767               do i=1,ndowna(0,il)
1768                call xdrfint_(ixdrf, ndowna(i,il), iret)
1769               enddo
1770              enddo
1771              do il=1,nodes
1772                do j=1,4
1773                 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1774                enddo
1775              enddo
1776 #else
1777               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1778
1779               do i=0,nodes-1
1780                call xdrfint(ixdrf, i2rep(i), iret)
1781               enddo
1782               do i=1,remd_m(1)
1783                call xdrfint(ixdrf, ifirst(i), iret)
1784               enddo
1785              do il=1,nodes
1786               call xdrfint(ixdrf, nupa(0,il), iret)
1787               do i=1,nupa(0,il)
1788                call xdrfint(ixdrf, nupa(i,il), iret)
1789               enddo
1790
1791               call xdrfint(ixdrf, ndowna(0,il), iret)
1792               do i=1,ndowna(0,il)
1793                call xdrfint(ixdrf, ndowna(i,il), iret)
1794               enddo
1795              enddo
1796              do il=1,nodes
1797                do j=1,4
1798                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1799                enddo
1800              enddo
1801 #endif
1802          endif
1803          call mpi_scatter(t_restart1,5,mpi_real,
1804      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1805          totT=t5_restart1(1)
1806          EK=t5_restart1(2)
1807          potE=t5_restart1(3)
1808          t_bath=t5_restart1(4)
1809
1810          if(me.eq.king)then
1811               do il=0,nodes-1
1812                do i=1,2*nres
1813 c                read(irest2,'(3e15.5)')
1814 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1815             do j=1,3
1816 #ifdef AIX
1817              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1818 #else
1819              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1820 #endif
1821             enddo
1822                enddo
1823 #ifdef DEBUG
1824             write (iout,*) "Conformation read",il
1825             do i=1,nres
1826               write (iout,'(i5,3f10.5,5x,3f10.5)')
1827      &          i,(d_restart1(j,i+2*nres*il),j=1,3),
1828      &            (d_restart1(j,nres+i+2*nres*il),j=1,3)
1829             enddo
1830 #endif
1831               enddo
1832          endif
1833          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1834      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1835
1836          do i=1,2*nres
1837            do j=1,3
1838             d_t(j,i)=r_d(j,i)
1839            enddo
1840          enddo
1841          if(me.eq.king)then
1842               do il=0,nodes-1
1843                do i=1,2*nres
1844 c                read(irest2,'(3e15.5)')
1845 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1846             do j=1,3
1847 #ifdef AIX
1848              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1849 #else
1850              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1851 #endif
1852             enddo
1853                enddo
1854               enddo
1855          endif
1856          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1857      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1858          do i=1,2*nres
1859            do j=1,3
1860             dc(j,i)=r_d(j,i)
1861            enddo
1862          enddo
1863
1864
1865            if(usampl.or.homol_nset.gt.1) then
1866 #ifdef AIX
1867              if(me.eq.king)then
1868               call xdrfint_(ixdrf, nset, iret)
1869               do i=1,nset
1870                 call xdrfint_(ixdrf,mset(i), iret)
1871               enddo
1872               do i=0,nodes-1
1873                 call xdrfint_(ixdrf,i2set(i), iret)
1874               enddo
1875               do il=1,nset
1876                do il1=1,mset(il)
1877                 do i=1,nrep
1878                  do j=1,remd_m(i)
1879                    call xdrfint_(ixdrf,itmp, iret)
1880                    i_index(i,j,il,il1)=itmp
1881                  enddo
1882                 enddo
1883                enddo
1884               enddo
1885              endif
1886 #else
1887              if(me.eq.king)then
1888               call xdrfint(ixdrf, nset, iret)
1889               do i=1,nset
1890                 call xdrfint(ixdrf,mset(i), iret)
1891               enddo
1892               do i=0,nodes-1
1893                 call xdrfint(ixdrf,i2set(i), iret)
1894               enddo
1895               do il=1,nset
1896                do il1=1,mset(il)
1897                 do i=1,nrep
1898                  do j=1,remd_m(i)
1899                    call xdrfint(ixdrf,itmp, iret)
1900                    i_index(i,j,il,il1)=itmp
1901                  enddo
1902                 enddo
1903                enddo
1904               enddo
1905              endif
1906 #endif
1907 c Corrected AL 8/19/2014: each processor needs whole iset array not only its
1908 c own element
1909 c              call mpi_scatter(i2set,1,mpi_integer,
1910 c     &           iset,1,mpi_integer,king,
1911 c     &           CG_COMM,ierr)
1912               call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1913      &         CG_COMM,ierr)
1914               iset=i2set(me)
1915 c broadcast iset to slaves
1916               if (nfgtasks.gt.1) then
1917                call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
1918                call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
1919               endif
1920            endif
1921
1922
1923         if(me.eq.king) close(irest2)
1924         return
1925         end
1926
1927       subroutine read1restart_old
1928       implicit real*8 (a-h,o-z)
1929       include 'DIMENSIONS'
1930       include 'mpif.h'
1931       include 'COMMON.MD'
1932       include 'COMMON.IOUNITS'
1933       include 'COMMON.REMD'
1934       include 'COMMON.SETUP'
1935       include 'COMMON.CHAIN'
1936       include 'COMMON.SBRIDGE'
1937       include 'COMMON.INTERACT'
1938       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1939      &                 t5_restart1(5)
1940       common /przechowalnia/ d_restart1
1941          if(me.eq.king)then
1942              open(irest2,file=mremd_rst_name,status='unknown')
1943              read (irest2,*) (i2rep(i),i=0,nodes-1)
1944              read (irest2,*) (ifirst(i),i=1,remd_m(1))
1945              do il=1,nodes
1946               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1947               read (irest2,*) ndowna(0,il),
1948      &                    (ndowna(i,il),i=1,ndowna(0,il))
1949              enddo
1950              do il=1,nodes
1951                read(irest2,*) (t_restart1(j,il),j=1,4)
1952              enddo
1953          endif
1954          call mpi_scatter(t_restart1,5,mpi_real,
1955      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1956          totT=t5_restart1(1)
1957          EK=t5_restart1(2)
1958          potE=t5_restart1(3)
1959          t_bath=t5_restart1(4)
1960
1961          if(me.eq.king)then
1962               do il=0,nodes-1
1963                do i=1,2*nres
1964                 read(irest2,'(3e15.5)')
1965      &                (d_restart1(j,i+2*nres*il),j=1,3)
1966                enddo
1967               enddo
1968          endif
1969          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1970      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1971
1972          do i=1,2*nres
1973            do j=1,3
1974             d_t(j,i)=r_d(j,i)
1975            enddo
1976          enddo
1977          if(me.eq.king)then
1978               do il=0,nodes-1
1979                do i=1,2*nres
1980                 read(irest2,'(3e15.5)')
1981      &                (d_restart1(j,i+2*nres*il),j=1,3)
1982                enddo
1983               enddo
1984          endif
1985          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1986      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1987          do i=1,2*nres
1988            do j=1,3
1989             dc(j,i)=r_d(j,i)
1990            enddo
1991          enddo
1992         if(me.eq.king) close(irest2)
1993         return
1994         end
1995 c------------------------------------------