source/unres/src_MD/MREMD.F

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