source/unres/src_MD/MREMD.F

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