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
1472       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1473      &     d_restart2(3,2*maxres*maxprocs)
1474       real t5_restart1(5)
1475       integer iret,itmp
1476       integer*2 i_index
1477      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1478        common /przechowalnia/ d_restart1,d_restart2
1479
1480        t5_restart1(1)=totT
1481        t5_restart1(2)=EK
1482        t5_restart1(3)=potE
1483        t5_restart1(4)=t_bath
1484        t5_restart1(5)=Uconst
1485
1486        call mpi_gather(t5_restart1,5,mpi_real,
1487      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1488
1489
1490        do i=1,2*nres
1491          do j=1,3
1492            r_d(j,i)=d_t(j,i)
1493          enddo
1494        enddo
1495        call mpi_gather(r_d,3*2*nres,mpi_real,
1496      &           d_restart1,3*2*nres,mpi_real,king,
1497      &           CG_COMM,ierr)
1498
1499
1500        do i=1,2*nres
1501          do j=1,3
1502            r_d(j,i)=dc(j,i)
1503          enddo
1504        enddo
1505        call mpi_gather(r_d,3*2*nres,mpi_real,
1506      &           d_restart2,3*2*nres,mpi_real,king,
1507      &           CG_COMM,ierr)
1508
1509        if(me.eq.king) then
1510 #ifdef AIX
1511          call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1512          do i=0,nodes-1
1513           call xdrfint_(ixdrf, i2rep(i), iret)
1514          enddo
1515          do i=1,remd_m(1)
1516           call xdrfint_(ixdrf, ifirst(i), iret)
1517          enddo
1518          do il=1,nodes
1519               do i=0,nupa(0,il)
1520                call xdrfint_(ixdrf, nupa(i,il), iret)
1521               enddo
1522
1523               do i=0,ndowna(0,il)
1524                call xdrfint_(ixdrf, ndowna(i,il), iret)
1525               enddo
1526          enddo
1527
1528          do il=1,nodes
1529            do j=1,4
1530             call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1531            enddo
1532          enddo
1533
1534          do il=0,nodes-1
1535            do i=1,2*nres
1536             do j=1,3
1537              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1538             enddo
1539            enddo
1540          enddo
1541          do il=0,nodes-1
1542            do i=1,2*nres
1543             do j=1,3
1544              call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1545             enddo
1546            enddo
1547          enddo
1548
1549          if(usampl.or.homol_nset.gt.1) then
1550            call xdrfint_(ixdrf, nset, iret)
1551            do i=1,nset
1552              call xdrfint_(ixdrf,mset(i), iret)
1553            enddo
1554            do i=0,nodes-1
1555              call xdrfint_(ixdrf,i2set(i), iret)
1556            enddo
1557            do il=1,nset
1558              do il1=1,mset(il)
1559                do i=1,nrep
1560                  do j=1,remd_m(i)
1561                    itmp=i_index(i,j,il,il1)
1562                    call xdrfint_(ixdrf,itmp, iret)
1563                  enddo
1564                enddo
1565              enddo
1566            enddo
1567
1568          endif
1569          call xdrfclose_(ixdrf, iret)
1570 #else
1571          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1572          do i=0,nodes-1
1573           call xdrfint(ixdrf, i2rep(i), iret)
1574          enddo
1575          do i=1,remd_m(1)
1576           call xdrfint(ixdrf, ifirst(i), iret)
1577          enddo
1578          do il=1,nodes
1579               do i=0,nupa(0,il)
1580                call xdrfint(ixdrf, nupa(i,il), iret)
1581               enddo
1582
1583               do i=0,ndowna(0,il)
1584                call xdrfint(ixdrf, ndowna(i,il), iret)
1585               enddo
1586          enddo
1587
1588          do il=1,nodes
1589            do j=1,4
1590             call xdrffloat(ixdrf, t_restart1(j,il), iret)
1591            enddo
1592          enddo
1593
1594          do il=0,nodes-1
1595            do i=1,2*nres
1596             do j=1,3
1597              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1598             enddo
1599            enddo
1600          enddo
1601          do il=0,nodes-1
1602            do i=1,2*nres
1603             do j=1,3
1604              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1605             enddo
1606            enddo
1607          enddo
1608
1609
1610              if(usampl.or.homol_nset.gt.1) then
1611               call xdrfint(ixdrf, nset, iret)
1612               do i=1,nset
1613                 call xdrfint(ixdrf,mset(i), iret)
1614               enddo
1615               do i=0,nodes-1
1616                 call xdrfint(ixdrf,i2set(i), iret)
1617               enddo
1618               do il=1,nset
1619                do il1=1,mset(il)
1620                 do i=1,nrep
1621                  do j=1,remd_m(i)
1622                    itmp=i_index(i,j,il,il1)
1623                    call xdrfint(ixdrf,itmp, iret)
1624                  enddo
1625                 enddo
1626                enddo
1627               enddo
1628
1629              endif
1630          call xdrfclose(ixdrf, iret)
1631 #endif
1632        endif
1633       return
1634       end
1635
1636
1637       subroutine write1traj
1638       implicit real*8 (a-h,o-z)
1639       include 'DIMENSIONS'
1640       include 'mpif.h'
1641       include 'COMMON.MD'
1642       include 'COMMON.IOUNITS'
1643       include 'COMMON.REMD'
1644       include 'COMMON.SETUP'
1645       include 'COMMON.CHAIN'
1646       include 'COMMON.SBRIDGE'
1647       include 'COMMON.INTERACT'
1648
1649       real t5_restart1(5)
1650       integer iret,itmp
1651       real xcoord(3,maxres2+2),prec
1652       real r_qfrag(50),r_qpair(100)
1653       real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1654       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1655       real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1656      &     p_uscdiff(100*maxprocs)
1657       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1658       common /przechowalnia/ p_c
1659
1660       call mpi_bcast(ii_write,1,mpi_integer,
1661      &           king,CG_COMM,ierr)
1662
1663 c debugging
1664       print *,'traj1file',me,ii_write,ntwx_cache
1665 c end debugging
1666
1667 #ifdef AIX
1668       if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1669 #else
1670       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1671 #endif
1672       do ii=1,ii_write
1673        t5_restart1(1)=totT_cache(ii)
1674        t5_restart1(2)=EK_cache(ii)
1675        t5_restart1(3)=potE_cache(ii)
1676        t5_restart1(4)=t_bath_cache(ii)
1677        t5_restart1(5)=Uconst_cache(ii)
1678        call mpi_gather(t5_restart1,5,mpi_real,
1679      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1680
1681        call mpi_gather(iset_cache(ii),1,mpi_integer,
1682      &      iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1683
1684           do i=1,nfrag
1685            r_qfrag(i)=qfrag_cache(i,ii)
1686           enddo
1687           do i=1,npair
1688            r_qpair(i)=qpair_cache(i,ii)
1689           enddo
1690           do i=1,nfrag_back
1691            r_utheta(i)=utheta_cache(i,ii)
1692            r_ugamma(i)=ugamma_cache(i,ii)
1693            r_uscdiff(i)=uscdiff_cache(i,ii)
1694           enddo
1695
1696         call mpi_gather(r_qfrag,nfrag,mpi_real,
1697      &           p_qfrag,nfrag,mpi_real,king,
1698      &           CG_COMM,ierr)
1699         call mpi_gather(r_qpair,npair,mpi_real,
1700      &           p_qpair,npair,mpi_real,king,
1701      &           CG_COMM,ierr)
1702         call mpi_gather(r_utheta,nfrag_back,mpi_real,
1703      &           p_utheta,nfrag_back,mpi_real,king,
1704      &           CG_COMM,ierr)
1705         call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1706      &           p_ugamma,nfrag_back,mpi_real,king,
1707      &           CG_COMM,ierr)
1708         call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1709      &           p_uscdiff,nfrag_back,mpi_real,king,
1710      &           CG_COMM,ierr)
1711
1712 #ifdef DEBUG
1713         write (iout,*) "p_qfrag"
1714         do i=1,nodes
1715           write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1716         enddo
1717         write (iout,*) "p_qpair"
1718         do i=1,nodes
1719           write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1720         enddo
1721 ctime        call flush(iout)
1722 #endif
1723         do i=1,nres*2
1724          do j=1,3
1725           r_c(j,i)=c_cache(j,i,ii)
1726          enddo
1727         enddo
1728
1729         call mpi_gather(r_c,3*2*nres,mpi_real,
1730      &           p_c,3*2*nres,mpi_real,king,
1731      &           CG_COMM,ierr)
1732
1733        if(me.eq.king) then
1734 #ifdef AIX
1735          do il=1,nodes
1736           call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1737           call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1738           call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1739           call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1740           call xdrfint_(ixdrf, nss, iret)
1741           do j=1,nss
1742            if (dyn_ss) then
1743             call xdrfint_(ixdrf, idssb(j)+nres, iret)
1744             call xdrfint_(ixdrf, jdssb(j)+nres, iret)
1745            else
1746             call xdrfint_(ixdrf, ihpb(j), iret)
1747             call xdrfint_(ixdrf, jhpb(j), iret)
1748            endif
1749           enddo
1750           call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1751           call xdrfint_(ixdrf, iset_restart1(il), iret)
1752           do i=1,nfrag
1753            call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1754           enddo
1755           do i=1,npair
1756            call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1757           enddo
1758           do i=1,nfrag_back
1759            call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1760            call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1761            call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1762           enddo
1763           prec=10000.0
1764           do i=1,nres
1765            do j=1,3
1766             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1767            enddo
1768           enddo
1769           do i=nnt,nct
1770            do j=1,3
1771             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1772            enddo
1773           enddo
1774           itmp=nres+nct-nnt+1
1775           call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1776          enddo
1777 #else
1778          do il=1,nodes
1779           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1780           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1781           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1782           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1783           call xdrfint(ixdrf, nss, iret)
1784           do j=1,nss
1785            if (dyn_ss) then
1786             call xdrfint(ixdrf, idssb(j)+nres, iret)
1787             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1788            else
1789             call xdrfint(ixdrf, ihpb(j), iret)
1790             call xdrfint(ixdrf, jhpb(j), iret)
1791            endif
1792           enddo
1793           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1794           call xdrfint(ixdrf, iset_restart1(il), iret)
1795           do i=1,nfrag
1796            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1797           enddo
1798           do i=1,npair
1799            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1800           enddo
1801           do i=1,nfrag_back
1802            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1803            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1804            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1805           enddo
1806           prec=10000.0
1807           do i=1,nres
1808            do j=1,3
1809             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1810            enddo
1811           enddo
1812           do i=nnt,nct
1813            do j=1,3
1814             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1815            enddo
1816           enddo
1817           itmp=nres+nct-nnt+1
1818           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1819          enddo
1820 #endif
1821        endif
1822       enddo
1823 #ifdef AIX
1824       if(me.eq.king) call xdrfclose_(ixdrf, iret)
1825 #else
1826       if(me.eq.king) call xdrfclose(ixdrf, iret)
1827 #endif
1828       do i=1,ntwx_cache-ii_write
1829
1830             totT_cache(i)=totT_cache(ii_write+i)
1831             EK_cache(i)=EK_cache(ii_write+i)
1832             potE_cache(i)=potE_cache(ii_write+i)
1833             t_bath_cache(i)=t_bath_cache(ii_write+i)
1834             Uconst_cache(i)=Uconst_cache(ii_write+i)
1835             iset_cache(i)=iset_cache(ii_write+i)
1836
1837             do ii=1,nfrag
1838              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1839             enddo
1840             do ii=1,npair
1841              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1842             enddo
1843             do ii=1,nfrag_back
1844               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1845               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1846               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1847             enddo
1848
1849             do ii=1,nres*2
1850              do j=1,3
1851               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1852              enddo
1853             enddo
1854       enddo
1855       ntwx_cache=ntwx_cache-ii_write
1856       return
1857       end
1858
1859
1860       subroutine read1restart(i_index)
1861       implicit real*8 (a-h,o-z)
1862       include 'DIMENSIONS'
1863       include 'mpif.h'
1864       include 'COMMON.MD'
1865       include 'COMMON.IOUNITS'
1866       include 'COMMON.REMD'
1867       include 'COMMON.SETUP'
1868       include 'COMMON.CHAIN'
1869       include 'COMMON.SBRIDGE'
1870       include 'COMMON.INTERACT'
1871       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1872      &                 t5_restart1(5)
1873       integer*2 i_index
1874      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1875       common /przechowalnia/ d_restart1
1876       write (*,*) "Processor",me," called read1restart"
1877
1878          if(me.eq.king)then
1879               open(irest2,file=mremd_rst_name,status='unknown')
1880               read(irest2,*,err=334) i
1881               write(iout,*) "Reading old rst in ASCI format"
1882               close(irest2)
1883                call read1restart_old
1884                return
1885  334          continue
1886 #ifdef AIX
1887               call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1888
1889               do i=0,nodes-1
1890                call xdrfint_(ixdrf, i2rep(i), iret)
1891               enddo
1892               do i=1,remd_m(1)
1893                call xdrfint_(ixdrf, ifirst(i), iret)
1894               enddo
1895              do il=1,nodes
1896               call xdrfint_(ixdrf, nupa(0,il), iret)
1897               do i=1,nupa(0,il)
1898                call xdrfint_(ixdrf, nupa(i,il), iret)
1899               enddo
1900
1901               call xdrfint_(ixdrf, ndowna(0,il), iret)
1902               do i=1,ndowna(0,il)
1903                call xdrfint_(ixdrf, ndowna(i,il), iret)
1904               enddo
1905              enddo
1906              do il=1,nodes
1907                do j=1,4
1908                 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1909                enddo
1910              enddo
1911 #else
1912               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1913
1914               do i=0,nodes-1
1915                call xdrfint(ixdrf, i2rep(i), iret)
1916               enddo
1917               do i=1,remd_m(1)
1918                call xdrfint(ixdrf, ifirst(i), iret)
1919               enddo
1920              do il=1,nodes
1921               call xdrfint(ixdrf, nupa(0,il), iret)
1922               do i=1,nupa(0,il)
1923                call xdrfint(ixdrf, nupa(i,il), iret)
1924               enddo
1925
1926               call xdrfint(ixdrf, ndowna(0,il), iret)
1927               do i=1,ndowna(0,il)
1928                call xdrfint(ixdrf, ndowna(i,il), iret)
1929               enddo
1930              enddo
1931              do il=1,nodes
1932                do j=1,4
1933                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1934                enddo
1935              enddo
1936 #endif
1937          endif
1938          call mpi_scatter(t_restart1,5,mpi_real,
1939      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1940          totT=t5_restart1(1)
1941          EK=t5_restart1(2)
1942          potE=t5_restart1(3)
1943          t_bath=t5_restart1(4)
1944
1945          if(me.eq.king)then
1946               do il=0,nodes-1
1947                do i=1,2*nres
1948 c                read(irest2,'(3e15.5)')
1949 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1950             do j=1,3
1951 #ifdef AIX
1952              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1953 #else
1954              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1955 #endif
1956             enddo
1957                enddo
1958               enddo
1959          endif
1960          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1961      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1962
1963          do i=1,2*nres
1964            do j=1,3
1965             d_t(j,i)=r_d(j,i)
1966            enddo
1967          enddo
1968          if(me.eq.king)then
1969               do il=0,nodes-1
1970                do i=1,2*nres
1971 c                read(irest2,'(3e15.5)')
1972 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1973             do j=1,3
1974 #ifdef AIX
1975              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1976 #else
1977              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1978 #endif
1979             enddo
1980                enddo
1981               enddo
1982          endif
1983          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1984      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1985          do i=1,2*nres
1986            do j=1,3
1987             dc(j,i)=r_d(j,i)
1988            enddo
1989          enddo
1990
1991
1992            if(usampl.or.homol_nset.gt.1) then
1993 #ifdef AIX
1994              if(me.eq.king)then
1995               call xdrfint_(ixdrf, nset, iret)
1996               do i=1,nset
1997                 call xdrfint_(ixdrf,mset(i), iret)
1998               enddo
1999               do i=0,nodes-1
2000                 call xdrfint_(ixdrf,i2set(i), iret)
2001               enddo
2002               do il=1,nset
2003                do il1=1,mset(il)
2004                 do i=1,nrep
2005                  do j=1,remd_m(i)
2006                    call xdrfint_(ixdrf,itmp, iret)
2007                    i_index(i,j,il,il1)=itmp
2008                  enddo
2009                 enddo
2010                enddo
2011               enddo
2012              endif
2013 #else
2014              if(me.eq.king)then
2015               call xdrfint(ixdrf, nset, iret)
2016               do i=1,nset
2017                 call xdrfint(ixdrf,mset(i), iret)
2018               enddo
2019               do i=0,nodes-1
2020                 call xdrfint(ixdrf,i2set(i), iret)
2021               enddo
2022               do il=1,nset
2023                do il1=1,mset(il)
2024                 do i=1,nrep
2025                  do j=1,remd_m(i)
2026                    call xdrfint(ixdrf,itmp, iret)
2027                    i_index(i,j,il,il1)=itmp
2028                  enddo
2029                 enddo
2030                enddo
2031               enddo
2032              endif
2033 #endif
2034 c Corrected AL 8/19/2014: each processor needs whole iset array not only its
2035 c own element
2036 c              call mpi_scatter(i2set,1,mpi_integer,
2037 c     &           iset,1,mpi_integer,king,
2038 c     &           CG_COMM,ierr)
2039               call mpi_bcast(i2set(0),nodes,mpi_integer,king,
2040      &         CG_COMM,ierr)
2041               iset=i2set(me)
2042
2043            endif
2044
2045
2046         if(me.eq.king) close(irest2)
2047         return
2048         end
2049
2050       subroutine read1restart_old
2051       implicit real*8 (a-h,o-z)
2052       include 'DIMENSIONS'
2053       include 'mpif.h'
2054       include 'COMMON.MD'
2055       include 'COMMON.IOUNITS'
2056       include 'COMMON.REMD'
2057       include 'COMMON.SETUP'
2058       include 'COMMON.CHAIN'
2059       include 'COMMON.SBRIDGE'
2060       include 'COMMON.INTERACT'
2061       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
2062      &                 t5_restart1(5)
2063       common /przechowalnia/ d_restart1
2064          if(me.eq.king)then
2065              open(irest2,file=mremd_rst_name,status='unknown')
2066              read (irest2,*) (i2rep(i),i=0,nodes-1)
2067              read (irest2,*) (ifirst(i),i=1,remd_m(1))
2068              do il=1,nodes
2069               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2070               read (irest2,*) ndowna(0,il),
2071      &                    (ndowna(i,il),i=1,ndowna(0,il))
2072              enddo
2073              do il=1,nodes
2074                read(irest2,*) (t_restart1(j,il),j=1,4)
2075              enddo
2076          endif
2077          call mpi_scatter(t_restart1,5,mpi_real,
2078      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2079          totT=t5_restart1(1)
2080          EK=t5_restart1(2)
2081          potE=t5_restart1(3)
2082          t_bath=t5_restart1(4)
2083
2084          if(me.eq.king)then
2085               do il=0,nodes-1
2086                do i=1,2*nres
2087                 read(irest2,'(3e15.5)')
2088      &                (d_restart1(j,i+2*nres*il),j=1,3)
2089                enddo
2090               enddo
2091          endif
2092          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2093      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2094
2095          do i=1,2*nres
2096            do j=1,3
2097             d_t(j,i)=r_d(j,i)
2098            enddo
2099          enddo
2100          if(me.eq.king)then
2101               do il=0,nodes-1
2102                do i=1,2*nres
2103                 read(irest2,'(3e15.5)')
2104      &                (d_restart1(j,i+2*nres*il),j=1,3)
2105                enddo
2106               enddo
2107          endif
2108          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2109      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2110          do i=1,2*nres
2111            do j=1,3
2112             dc(j,i)=r_d(j,i)
2113            enddo
2114          enddo
2115         if(me.eq.king) close(irest2)
2116         return
2117         end
2118 c-------------------------------------------------------------------
2119         subroutine set_hweights(iiset)
2120         implicit real*8 (a-h,o-z)
2121         integer i
2122         include 'DIMENSIONS'
2123         include 'COMMON.FFIELD'
2124         include 'COMMON.REMD'
2125
2126          do i=1,n_ene
2127           weights(i)=hweights(iiset,i)
2128          enddo
2129
2130          wsc    =weights(1)
2131          wscp   =weights(2)
2132          welec  =weights(3)
2133          wcorr  =weights(4)
2134          wcorr5 =weights(5)
2135          wcorr6 =weights(6)
2136          wel_loc=weights(7)
2137          wturn3 =weights(8)
2138          wturn4 =weights(9)
2139          wturn6 =weights(10)
2140          wang   =weights(11)
2141          wscloc =weights(12)
2142          wtor   =weights(13)
2143          wtor_d =weights(14)
2144          wstrain=weights(15)
2145          wvdwpp =weights(16)
2146          wbond  =weights(17)
2147          scal14 =weights(18)
2148          wsccor =weights(21)
2149
2150         return
2151         end
2152 #endif