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