source/unres/src_MD/MREMD.F

   1 #define DEBUG
   2 #ifdef MPI
   3       subroutine MREMD
   4       implicit real*8 (a-h,o-z)
   5       include 'DIMENSIONS'
   6       include 'mpif.h'
   7       include 'COMMON.CONTROL'
   8       include 'COMMON.VAR'
   9       include 'COMMON.MD'
  10 #ifndef LANG0
  11       include 'COMMON.LANGEVIN'
  12 #else
  13       include 'COMMON.LANGEVIN.lang0'
  14 #endif
  15       include 'COMMON.CHAIN'
  16       include 'COMMON.DERIV'
  17       include 'COMMON.GEO'
  18       include 'COMMON.LOCAL'
  19       include 'COMMON.INTERACT'
  20       include 'COMMON.IOUNITS'
  21       include 'COMMON.NAMES'
  22       include 'COMMON.TIME1'
  23       include 'COMMON.REMD'
  24       include 'COMMON.SETUP'
  25       include 'COMMON.MUCA'
  26       integer ERRCODE
  27       double precision cm(3),L(3),vcm(3)
  28       double precision energia(0:n_ene)
  29       double precision remd_t_bath(maxprocs)
  30       integer iremd_iset(maxprocs)
  31       integer*2 i_index
  32      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
  33       double precision remd_ene(0:n_ene+4,maxprocs),t_bath_old,e_tmp
  34       integer iremd_acc(maxprocs),iremd_tot(maxprocs)
  35       integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
  36       integer ilen,rstcount
  37       external ilen
  38       character*50 tytul
  39       common /gucio/ cm
  40       integer itime
  41 cold      integer nup(0:maxprocs),ndown(0:maxprocs)
  42       integer rep2i(0:maxprocs),ireqi(maxprocs)
  43       integer icache_all(maxprocs)
  44       integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
  45       logical synflag,end_of_run,file_exist /.false./,ovrtim
  46       real ene_tol /1.0e-5/
  47
  48 cdeb      imin_itime_old=0
  49       ntwx_cache=0
  50       time00=MPI_WTIME()
  51       time01=time00
  52       if(me.eq.king.or..not.out1file) then
  53        write  (iout,*) 'MREMD',nodes,'time before',time00-walltime
  54        write (iout,*) "NREP=",nrep
  55       endif
  56
  57       synflag=.false.
  58       if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
  59         call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
  60       endif
  61       mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
  62
  63 cd      print *,'MREMD',nodes,homol_nset
  64 cd      print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
  65 cde      write (iout,*) "Start MREMD: me",me," t_bath",t_bath
  66
  67       if(hremd.gt.0) then
  68          nset=hremd
  69          do i=1,nset
  70           mset(i)=1
  71          enddo
  72       endif
  73
  74       if(homol_nset.gt.1) then
  75          i_econstr=24
  76          nset=homol_nset
  77          do i=1,nset
  78           mset(i)=1
  79          enddo
  80       endif
  81
  82       if(usampl) i_econstr=20
  83
  84       k=0
  85       rep2i(k)=-1
  86       do il=1,max0(nset,1)
  87        do il1=1,max0(mset(il),1)
  88         do i=1,nrep
  89          iremd_acc(i)=0
  90          iremd_acc_usa(i)=0
  91          iremd_tot(i)=0
  92          do j=1,remd_m(i)
  93           i2rep(k)=i
  94           i2set(k)=il
  95           rep2i(i)=k
  96           k=k+1
  97           i_index(i,j,il,il1)=k
  98          enddo
  99         enddo
 100        enddo
 101       enddo
 102
 103       if(me.eq.king.or..not.out1file) then
 104        write(iout,*) "i2rep",(i2rep(i),i=0,nodes-1)
 105        write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
 106        write(iout,*) "i,j,il,il1,i_index(i,j,il,il1)"
 107        do il=1,nset
 108         do il1=1,mset(il)
 109          do i=1,nrep
 110           do j=1,remd_m(i)
 111            write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
 112           enddo
 113          enddo
 114         enddo
 115        enddo
 116       endif
 117
 118 c      print *,'i2rep',me,i2rep(me)
 119 c      print *,'rep2i',(rep2i(i),i=0,nrep)
 120
 121 cold       if (i2rep(me).eq.nrep) then
 122 cold        nup(0)=0
 123 cold       else
 124 cold        nup(0)=remd_m(i2rep(me)+1)
 125 cold        k=rep2i(int(i2rep(me)))+1
 126 cold        do i=1,nup(0)
 127 cold         nup(i)=k
 128 cold         k=k+1
 129 cold        enddo
 130 cold       endif
 131
 132 cd       print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
 133
 134 cold       if (i2rep(me).eq.1) then
 135 cold        ndown(0)=0
 136 cold       else
 137 cold        ndown(0)=remd_m(i2rep(me)-1)
 138 cold        k=rep2i(i2rep(me)-2)+1
 139 cold        do i=1,ndown(0)
 140 cold         ndown(i)=k
 141 cold         k=k+1
 142 cold        enddo
 143 cold       endif
 144
 145 cd       print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
 146
 147
 148        write (*,*) "Processor",me," rest",rest,"
 149      &   restart1fie",restart1file
 150        if(rest.and.restart1file) then
 151            if (me.eq.king)
 152      &     inquire(file=mremd_rst_name,exist=file_exist)
 153 cd           write (*,*) me," Before broadcast: file_exist",file_exist
 154            call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
 155      &          IERR)
 156 cd           write (*,*) me," After broadcast: file_exist",file_exist
 157            if(file_exist) then
 158              if(me.eq.king.or..not.out1file)
 159      &            write  (iout,*) 'Master is reading restart1file'
 160              call read1restart(i_index)
 161            else
 162              if(me.eq.king.or..not.out1file)
 163      &            write  (iout,*) 'WARNING : no restart1file'
 164            endif
 165
 166            if(me.eq.king.or..not.out1file) then
 167               write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
 168               write(iout,*) "i_index"
 169               do il=1,nset
 170                do il1=1,mset(il)
 171                 do i=1,nrep
 172                  do j=1,remd_m(i)
 173                   write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
 174                  enddo
 175                 enddo
 176                enddo
 177               enddo
 178            endif
 179        endif
 180
 181        if(me.eq.king) then
 182         if (rest.and..not.restart1file)
 183      &          inquire(file=mremd_rst_name,exist=file_exist)
 184         if(.not.file_exist.and.rest.and..not.restart1file)
 185      &       write(iout,*) 'WARNING : no restart file',mremd_rst_name
 186         IF (rest.and.file_exist.and..not.restart1file) THEN
 187              write  (iout,*) 'Master is reading restart file',
 188      &                        mremd_rst_name
 189              open(irest2,file=mremd_rst_name,status='unknown')
 190              read (irest2,*)
 191              read (irest2,*) (i2rep(i),i=0,nodes-1)
 192              read (irest2,*)
 193              read (irest2,*) (ifirst(i),i=1,remd_m(1))
 194              do il=1,nodes
 195               read (irest2,*)
 196               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 197               read (irest2,*)
 198               read (irest2,*) ndowna(0,il),
 199      &                    (ndowna(i,il),i=1,ndowna(0,il))
 200              enddo
 201              if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
 202               read (irest2,*)
 203               read (irest2,*) nset
 204               read (irest2,*)
 205               read (irest2,*) (mset(i),i=1,nset)
 206               read (irest2,*)
 207               read (irest2,*) (i2set(i),i=0,nodes-1)
 208               read (irest2,*)
 209               do il=1,nset
 210                do il1=1,mset(il)
 211                 do i=1,nrep
 212                   read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
 213                 enddo
 214                enddo
 215               enddo
 216
 217               write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
 218               write(iout,*) "i_index"
 219               do il=1,nset
 220                do il1=1,mset(il)
 221                 do i=1,nrep
 222                  do j=1,remd_m(i)
 223                   write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
 224                  enddo
 225                 enddo
 226                enddo
 227               enddo
 228              endif
 229
 230              close(irest2)
 231
 232              write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
 233              write (iout,'(a6,1000i5)') "ifirst",
 234      &                    (ifirst(i),i=1,remd_m(1))
 235              do il=1,nodes
 236               write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
 237      &                    (nupa(i,il),i=1,nupa(0,il))
 238               write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
 239      &                    (ndowna(i,il),i=1,ndowna(0,il))
 240              enddo
 241         ELSE IF (.not.(rest.and.file_exist)) THEN
 242          do il=1,remd_m(1)
 243           ifirst(il)=il
 244          enddo
 245
 246          do il=1,nodes
 247           if (i2rep(il-1).eq.nrep) then
 248            nupa(0,il)=0
 249           else
 250            nupa(0,il)=remd_m(i2rep(il-1)+1)
 251            k=rep2i(int(i2rep(il-1)))+1
 252            do i=1,nupa(0,il)
 253             nupa(i,il)=k+1
 254             k=k+1
 255            enddo
 256           endif
 257           if (i2rep(il-1).eq.1) then
 258            ndowna(0,il)=0
 259           else
 260            ndowna(0,il)=remd_m(i2rep(il-1)-1)
 261            k=rep2i(i2rep(il-1)-2)+1
 262            do i=1,ndowna(0,il)
 263             ndowna(i,il)=k+1
 264             k=k+1
 265            enddo
 266           endif
 267          enddo
 268
 269         write (iout,'(a6,100i4)') "ifirst",
 270      &                    (ifirst(i),i=1,remd_m(1))
 271         do il=1,nodes
 272          write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
 273      &                    (nupa(i,il),i=1,nupa(0,il))
 274          write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
 275      &                    (ndowna(i,il),i=1,ndowna(0,il))
 276         enddo
 277
 278         ENDIF
 279        endif
 280 c
 281 c      t_bath=retmin+(retmax-retmin)*me/(nodes-1)
 282        if(.not.(rest.and.file_exist.and.restart1file)) then
 283          if (me .eq. king) then
 284             t_bath=retmin
 285          else
 286             t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
 287          endif
 288 cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
 289          if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
 290
 291        endif
 292        if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
 293           iset=i2set(me)
 294 c broadcast iset to slaves
 295           if (nfgtasks.gt.1) then
 296            call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
 297            call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
 298           endif
 299           if (hremd.gt.0) call set_hweights(iset)
 300           if(me.eq.king.or..not.out1file)
 301      &     write(iout,*) me,"iset=",iset,"t_bath=",t_bath
 302        endif
 303 c
 304        stdfp=dsqrt(2*Rb*t_bath/d_time)
 305        do i=1,ntyp
 306           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 307        enddo
 308
 309 c      print *,'irep',me,t_bath
 310        if (.not.rest) then
 311         if (me.eq.king .or. .not. out1file)
 312      &   write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
 313         call rescale_weights(t_bath)
 314        endif
 315
 316
 317 c------copy MD--------------
 318 c  The driver for molecular dynamics subroutines
 319 c------------------------------------------------
 320       t_MDsetup=0.0d0
 321       t_langsetup=0.0d0
 322       t_MD=0.0d0
 323       t_enegrad=0.0d0
 324       t_sdsetup=0.0d0
 325       if(me.eq.king.or..not.out1file)
 326      & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
 327 #ifdef MPI
 328       tt0 = MPI_Wtime()
 329 #else
 330       tt0 = tcpu()
 331 #endif
 332 c Determine the inverse of the inertia matrix.
 333       call setup_MD_matrices
 334 c Initialize MD
 335       call init_MD
 336       if (rest) then
 337        if (me.eq.king .or. .not. out1file)
 338      &  write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
 339        stdfp=dsqrt(2*Rb*t_bath/d_time)
 340        do i=1,ntyp
 341           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 342        enddo
 343        if (lang.gt.0 .and. .not.surfarea) then
 344          do i=nnt,nct-1
 345           stdforcp(i)=stdfp*dsqrt(gamp)
 346          enddo
 347          do i=nnt,nct
 348           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
 349          enddo
 350        elseif (lang.gt.0 .and. surfarea ) then
 351           call setup_fricmat
 352        endif
 353        call rescale_weights(t_bath)
 354       endif
 355
 356 #ifdef MPI
 357       t_MDsetup = MPI_Wtime()-tt0
 358 #else
 359       t_MDsetup = tcpu()-tt0
 360 #endif
 361       rstcount=0
 362 c   Entering the MD loop
 363 #ifdef MPI
 364       tt0 = MPI_Wtime()
 365 #else
 366       tt0 = tcpu()
 367 #endif
 368       if (lang.eq.2 .or. lang.eq.3) then
 369 #ifndef   LANG0
 370         call setup_fricmat
 371         if (lang.eq.2) then
 372           call sd_verlet_p_setup
 373         else
 374           call sd_verlet_ciccotti_setup
 375         endif
 376         do i=1,dimen
 377           do j=1,dimen
 378             pfric0_mat(i,j,0)=pfric_mat(i,j)
 379             afric0_mat(i,j,0)=afric_mat(i,j)
 380             vfric0_mat(i,j,0)=vfric_mat(i,j)
 381             prand0_mat(i,j,0)=prand_mat(i,j)
 382             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 383             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 384           enddo
 385         enddo
 386         flag_stoch(0)=.true.
 387         do i=1,maxflag_stoch
 388           flag_stoch(i)=.false.
 389         enddo
 390 #else
 391         write (iout,*)
 392      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 393 #ifdef MPI
 394         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 395 #endif
 396         stop
 397 #endif
 398       else if (lang.eq.1 .or. lang.eq.4) then
 399         call setup_fricmat
 400       endif
 401       time00=MPI_WTIME()
 402       if (me.eq.king .or. .not. out1file)
 403      & write(iout,*) 'Setup time',time00-walltime
 404 ctime      call flush(iout)
 405 #ifdef MPI
 406       t_langsetup=MPI_Wtime()-tt0
 407       tt0=MPI_Wtime()
 408 #else
 409       t_langsetup=tcpu()-tt0
 410       tt0=tcpu()
 411 #endif
 412       itime=0
 413       end_of_run=.false.
 414       do while(.not.end_of_run)
 415         itime=itime+1
 416         if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
 417         if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
 418         rstcount=rstcount+1
 419         if (lang.gt.0 .and. surfarea .and.
 420      &      mod(itime,reset_fricmat).eq.0) then
 421           if (lang.eq.2 .or. lang.eq.3) then
 422 #ifndef   LANG0
 423             call setup_fricmat
 424             if (lang.eq.2) then
 425               call sd_verlet_p_setup
 426             else
 427               call sd_verlet_ciccotti_setup
 428             endif
 429             do i=1,dimen
 430               do j=1,dimen
 431                 pfric0_mat(i,j,0)=pfric_mat(i,j)
 432                 afric0_mat(i,j,0)=afric_mat(i,j)
 433                 vfric0_mat(i,j,0)=vfric_mat(i,j)
 434                 prand0_mat(i,j,0)=prand_mat(i,j)
 435                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 436                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 437               enddo
 438             enddo
 439             flag_stoch(0)=.true.
 440             do i=1,maxflag_stoch
 441               flag_stoch(i)=.false.
 442             enddo
 443 #endif
 444           else if (lang.eq.1 .or. lang.eq.4) then
 445             call setup_fricmat
 446           endif
 447           write (iout,'(a,i10)')
 448      &      "Friction matrix reset based on surface area, itime",itime
 449         endif
 450         if (reset_vel .and. tbf .and. lang.eq.0
 451      &      .and. mod(itime,count_reset_vel).eq.0) then
 452           call random_vel
 453           if (me.eq.king .or. .not. out1file)
 454      &     write(iout,'(a,f20.2)')
 455      &     "Velocities reset to random values, time",totT
 456           do i=0,2*nres
 457             do j=1,3
 458               d_t_old(j,i)=d_t(j,i)
 459             enddo
 460           enddo
 461         endif
 462         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 463           call inertia_tensor
 464           call vcm_vel(vcm)
 465           do j=1,3
 466              d_t(j,0)=d_t(j,0)-vcm(j)
 467           enddo
 468           call kinetic(EK)
 469           kinetic_T=2.0d0/(dimen3*Rb)*EK
 470           scalfac=dsqrt(T_bath/kinetic_T)
 471 cd          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
 472           do i=0,2*nres
 473             do j=1,3
 474               d_t_old(j,i)=scalfac*d_t(j,i)
 475             enddo
 476           enddo
 477         endif
 478         if (lang.ne.4) then
 479           if (RESPA) then
 480 c Time-reversible RESPA algorithm
 481 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 482             call RESPA_step(itime)
 483           else
 484 c Variable time step algorithm.
 485             call velverlet_step(itime)
 486           endif
 487         else
 488 #ifdef BROWN
 489           call brown_step(itime)
 490 #else
 491           print *,"Brown dynamics not here!"
 492 #ifdef MPI
 493           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 494 #endif
 495           stop
 496 #endif
 497         endif
 498         if(hmc.gt.0 .and. mod(itime,hmc).eq.0) then
 499           call statout(itime)
 500           call hmc_test(itime)
 501         endif
 502         if(ntwe.ne.0) then
 503           if (mod(itime,ntwe).eq.0) call statout(itime)
 504         endif
 505         if (mod(itime,ntwx).eq.0.and..not.traj1file) then
 506           write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
 507           if(mdpdb) then
 508              call pdbout(potE,tytul,ipdb)
 509           else
 510              call cartout(totT)
 511           endif
 512         endif
 513         if (mod(itime,ntwx).eq.0.and.traj1file) then
 514           if(ntwx_cache.lt.max_cache_traj_use) then
 515             ntwx_cache=ntwx_cache+1
 516           else
 517            if (max_cache_traj_use.ne.1)
 518      &      print *,itime,"processor ",me," over cache ",ntwx_cache
 519            do i=1,ntwx_cache-1
 520
 521             totT_cache(i)=totT_cache(i+1)
 522             EK_cache(i)=EK_cache(i+1)
 523             potE_cache(i)=potE_cache(i+1)
 524             t_bath_cache(i)=t_bath_cache(i+1)
 525             Uconst_cache(i)=Uconst_cache(i+1)
 526             iset_cache(i)=iset_cache(i+1)
 527
 528             do ii=1,nfrag
 529              qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
 530             enddo
 531             do ii=1,npair
 532              qpair_cache(ii,i)=qpair_cache(ii,i+1)
 533             enddo
 534             do ii=1,nfrag_back
 535               utheta_cache(ii,i)=utheta_cache(ii,i+1)
 536               ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
 537               uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
 538             enddo
 539
 540
 541             do ii=1,nres*2
 542              do j=1,3
 543               c_cache(j,ii,i)=c_cache(j,ii,i+1)
 544              enddo
 545             enddo
 546            enddo
 547           endif
 548
 549             totT_cache(ntwx_cache)=totT
 550             EK_cache(ntwx_cache)=EK
 551             potE_cache(ntwx_cache)=potE
 552             t_bath_cache(ntwx_cache)=t_bath
 553             Uconst_cache(ntwx_cache)=Uconst
 554             iset_cache(ntwx_cache)=iset
 555
 556             do i=1,nfrag
 557              qfrag_cache(i,ntwx_cache)=qfrag(i)
 558             enddo
 559             do i=1,npair
 560              qpair_cache(i,ntwx_cache)=qpair(i)
 561             enddo
 562             do i=1,nfrag_back
 563               utheta_cache(i,ntwx_cache)=utheta(i)
 564               ugamma_cache(i,ntwx_cache)=ugamma(i)
 565               uscdiff_cache(i,ntwx_cache)=uscdiff(i)
 566             enddo
 567
 568             do i=1,nres*2
 569              do j=1,3
 570               c_cache(j,i,ntwx_cache)=c(j,i)
 571              enddo
 572             enddo
 573
 574         endif
 575         if ((rstcount.eq.1000.or.itime.eq.n_timestep)
 576      &                         .and..not.restart1file) then
 577
 578            if(me.eq.king) then
 579              open(irest1,file=mremd_rst_name,status='unknown')
 580              write (irest1,*) "i2rep"
 581              write (irest1,*) (i2rep(i),i=0,nodes-1)
 582              write (irest1,*) "ifirst"
 583              write (irest1,*) (ifirst(i),i=1,remd_m(1))
 584              do il=1,nodes
 585               write (irest1,*) "nupa",il
 586               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 587               write (irest1,*) "ndowna",il
 588               write (irest1,*) ndowna(0,il),
 589      &                   (ndowna(i,il),i=1,ndowna(0,il))
 590              enddo
 591              if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
 592               write (irest1,*) "nset"
 593               write (irest1,*) nset
 594               write (irest1,*) "mset"
 595               write (irest1,*) (mset(i),i=1,nset)
 596               write (irest1,*) "i2set"
 597               write (irest1,*) (i2set(i),i=0,nodes-1)
 598               write (irest1,*) "i_index"
 599               do il=1,nset
 600                do il1=1,mset(il)
 601                 do i=1,nrep
 602                   write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
 603                 enddo
 604                enddo
 605               enddo
 606
 607              endif
 608              close(irest1)
 609            endif
 610            open(irest2,file=rest2name,status='unknown')
 611            write(irest2,*) totT,EK,potE,totE,t_bath
 612            do i=1,2*nres
 613             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 614            enddo
 615            do i=1,2*nres
 616             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
 617            enddo
 618            if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
 619              write (irest2,*) iset
 620            endif
 621           close(irest2)
 622           rstcount=0
 623         endif
 624
 625 c REMD - exchange
 626 c forced synchronization
 627         if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
 628      &                                .and. .not. mremdsync) then
 629             synflag=.false.
 630             call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
 631             if (synflag) then
 632                call mpi_recv(itime_master, 1, MPI_INTEGER,
 633      &                             0,101,CG_COMM, status, ierr)
 634                call mpi_barrier(CG_COMM, ierr)
 635 cdeb               if (out1file.or.traj1file) then
 636 cdeb                call mpi_gather(itime,1,mpi_integer,
 637 cdeb     &             icache_all,1,mpi_integer,king,
 638 cdeb     &             CG_COMM,ierr)
 639                if(traj1file)
 640      &          call mpi_gather(ntwx_cache,1,mpi_integer,
 641      &             icache_all,1,mpi_integer,king,
 642      &             CG_COMM,ierr)
 643                if (.not.out1file)
 644      &               write(iout,*) 'REMD synchro at',itime_master,itime
 645                if (itime_master.ge.n_timestep .or. ovrtim())
 646      &            end_of_run=.true.
 647 ctime               call flush(iout)
 648             endif
 649         endif
 650
 651 c REMD - exchange
 652         if ((mod(itime,nstex).eq.0.and.me.eq.king
 653      &                  .or.end_of_run.and.me.eq.king )
 654      &       .and. .not. mremdsync ) then
 655            synflag=.true.
 656            time01_=MPI_WTIME()
 657            do i=1,nodes-1
 658               call mpi_isend(itime,1,MPI_INTEGER,i,101,
 659      &                                CG_COMM, ireqi(i), ierr)
 660 cd            write(iout,*) 'REMD synchro with',i
 661 cd            call flush(iout)
 662            enddo
 663            call mpi_waitall(nodes-1,ireqi,statusi,ierr)
 664            call mpi_barrier(CG_COMM, ierr)
 665            time01=MPI_WTIME()
 666            write(iout,*) 'REMD synchro at',itime,'time=',time01-time01_
 667            if (out1file.or.traj1file) then
 668 cdeb            call mpi_gather(itime,1,mpi_integer,
 669 cdeb     &             itime_all,1,mpi_integer,king,
 670 cdeb     &             CG_COMM,ierr)
 671 cdeb            write(iout,'(a19,8000i8)') ' REMD synchro itime',
 672 cdeb     &                    (itime_all(i),i=1,nodes)
 673             if(traj1file) then
 674 cdeb             imin_itime=itime_all(1)
 675 cdeb             do i=2,nodes
 676 cdeb               if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
 677 cdeb             enddo
 678 cdeb             ii_write=(imin_itime-imin_itime_old)/ntwx
 679 cdeb             imin_itime_old=int(imin_itime/ntwx)*ntwx
 680 cdeb             write(iout,*) imin_itime,imin_itime_old,ii_write
 681              call mpi_gather(ntwx_cache,1,mpi_integer,
 682      &             icache_all,1,mpi_integer,king,
 683      &             CG_COMM,ierr)
 684 c             write(iout,'(a19,8000i8)') '     ntwx_cache',
 685 c     &                    (icache_all(i),i=1,nodes)
 686              ii_write=icache_all(1)
 687              do i=2,nodes
 688                if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
 689              enddo
 690 c             write(iout,*) "MIN ii_write=",ii_write
 691             endif
 692            endif
 693 ctime           call flush(iout)
 694         endif
 695         if(mremdsync .and. mod(itime,nstex).eq.0) then
 696            synflag=.true.
 697            if (me.eq.king .or. .not. out1file)
 698      &      write(iout,*) 'REMD synchro at',itime
 699
 700             if(traj1file) then
 701              call mpi_gather(ntwx_cache,1,mpi_integer,
 702      &             icache_all,1,mpi_integer,king,
 703      &             CG_COMM,ierr)
 704              if (me.eq.king) then
 705                write(iout,'(a19,8000i8)') '     ntwx_cache',
 706      &                    (icache_all(i),i=1,nodes)
 707                ii_write=icache_all(1)
 708                do i=2,nodes
 709                  if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
 710                enddo
 711                write(iout,*) "MIN ii_write=",ii_write
 712              endif
 713             endif
 714 ctest           call flush(iout)
 715         endif
 716         if (synflag) then
 717 c Update the time safety limiy
 718           if (time001-time00.gt.safety) then
 719             safety=time001-time00+600
 720              if (me.eq.king .or. .not. out1file)
 721      &       write (iout,*) "****** SAFETY increased to",safety," s"
 722           endif
 723           if (ovrtim()) end_of_run=.true.
 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)
1419      &    call mpi_scatter(iremd_iset,1,mpi_integer,
1420      &           iset,1,mpi_integer,king,
1421      &           CG_COMM,ierr)
1422
1423          time07=MPI_WTIME()
1424           if (me.eq.king .or. .not. out1file) then
1425             write(iout,*) 'REMD scatter time=',time07-time06
1426           endif
1427
1428          if(lmuca) then
1429            call mpi_scatter(elowi,1,mpi_double_precision,
1430      &           elow,1,mpi_double_precision,king,
1431      &           CG_COMM,ierr)
1432            call mpi_scatter(ehighi,1,mpi_double_precision,
1433      &           ehigh,1,mpi_double_precision,king,
1434      &           CG_COMM,ierr)
1435          endif
1436
1437          if(hremd.gt.0) call set_hweights(iset)
1438          call rescale_weights(t_bath)
1439 co         write (iout,*) "Processor",me,
1440 co     &    " rescaling weights with temperature",t_bath
1441
1442          stdfp=dsqrt(2*Rb*t_bath/d_time)
1443          do i=1,ntyp
1444            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1445          enddo
1446          if (lang.gt.0) then
1447            do i=nnt,nct-1
1448             stdforcp(i)=stdforcp(i)*sqrt(t_bath/t_bath_old)
1449            enddo
1450            do i=nnt,nct
1451             stdforcsc(i)=stdforcsc(i)*sqrt(t_bath/t_bath_old)
1452            enddo
1453          endif
1454 cde         write(iout,*) 'REMD after',me,t_bath
1455            time08=MPI_WTIME()
1456            if (me.eq.king .or. .not. out1file) then
1457             write(iout,*) 'REMD exchange time=',time08-time02
1458 ctime            call flush(iout)
1459            endif
1460         endif
1461       enddo
1462
1463       if (restart1file) then
1464           if (me.eq.king .or. .not. out1file)
1465      &      write(iout,*) 'writing restart at the end of run'
1466            call write1rst(i_index)
1467       endif
1468
1469       if (traj1file) call write1traj
1470 cd debugging
1471 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
1472 cdeb     &             icache_all,1,mpi_integer,king,
1473 cdeb     &             CG_COMM,ierr)
1474 cdeb            write(iout,'(a40,8000i8)')
1475 cdeb     &             '  ntwx_cache after traj1file at the end',
1476 cdeb     &             (icache_all(i),i=1,nodes)
1477 cd end
1478
1479
1480 #ifdef MPI
1481       t_MD=MPI_Wtime()-tt0
1482 #else
1483       t_MD=tcpu()-tt0
1484 #endif
1485       if (me.eq.king .or. .not. out1file) then
1486        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1487      &  '  Timing  ',
1488      & 'MD calculations setup:',t_MDsetup,
1489      & 'Energy & gradient evaluation:',t_enegrad,
1490      & 'Stochastic MD setup:',t_langsetup,
1491      & 'Stochastic MD step setup:',t_sdsetup,
1492      & 'MD steps:',t_MD
1493        write (iout,'(/28(1h=),a25,27(1h=))')
1494      & '  End of MD calculation  '
1495        if(hmc.gt.0) write (iout,*) 'HMC acceptance ratio',
1496      &         n_timestep*1.0d0/hmc/hmc_acc
1497       endif
1498       return
1499       end
1500
1501 c-----------------------------------------------------------------------
1502       subroutine write1rst(i_index)
1503       implicit real*8 (a-h,o-z)
1504       include 'DIMENSIONS'
1505       include 'mpif.h'
1506       include 'COMMON.MD'
1507       include 'COMMON.IOUNITS'
1508       include 'COMMON.REMD'
1509       include 'COMMON.SETUP'
1510       include 'COMMON.CHAIN'
1511       include 'COMMON.SBRIDGE'
1512       include 'COMMON.INTERACT'
1513       include 'COMMON.CONTROL'
1514
1515       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1516      &     d_restart2(3,2*maxres*maxprocs)
1517       real t5_restart1(5)
1518       integer iret,itmp
1519       integer*2 i_index
1520      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1521        common /przechowalnia/ d_restart1,d_restart2
1522
1523        t5_restart1(1)=totT
1524        t5_restart1(2)=EK
1525        t5_restart1(3)=potE
1526        t5_restart1(4)=t_bath
1527        t5_restart1(5)=Uconst
1528
1529        call mpi_gather(t5_restart1,5,mpi_real,
1530      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1531
1532
1533        do i=1,2*nres
1534          do j=1,3
1535            r_d(j,i)=d_t(j,i)
1536          enddo
1537        enddo
1538        call mpi_gather(r_d,3*2*nres,mpi_real,
1539      &           d_restart1,3*2*nres,mpi_real,king,
1540      &           CG_COMM,ierr)
1541
1542
1543        do i=1,2*nres
1544          do j=1,3
1545            r_d(j,i)=dc(j,i)
1546          enddo
1547        enddo
1548        call mpi_gather(r_d,3*2*nres,mpi_real,
1549      &           d_restart2,3*2*nres,mpi_real,king,
1550      &           CG_COMM,ierr)
1551
1552        if(me.eq.king) then
1553 #ifdef AIX
1554          call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1555          do i=0,nodes-1
1556           call xdrfint_(ixdrf, i2rep(i), iret)
1557          enddo
1558          do i=1,remd_m(1)
1559           call xdrfint_(ixdrf, ifirst(i), iret)
1560          enddo
1561          do il=1,nodes
1562               do i=0,nupa(0,il)
1563                call xdrfint_(ixdrf, nupa(i,il), iret)
1564               enddo
1565
1566               do i=0,ndowna(0,il)
1567                call xdrfint_(ixdrf, ndowna(i,il), iret)
1568               enddo
1569          enddo
1570
1571          do il=1,nodes
1572            do j=1,4
1573             call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1574            enddo
1575          enddo
1576
1577          do il=0,nodes-1
1578            do i=1,2*nres
1579             do j=1,3
1580              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1581             enddo
1582            enddo
1583          enddo
1584          do il=0,nodes-1
1585            do i=1,2*nres
1586             do j=1,3
1587              call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1588             enddo
1589            enddo
1590          enddo
1591
1592          if(usampl.or.homol_nset.gt.1) then
1593            call xdrfint_(ixdrf, nset, iret)
1594            do i=1,nset
1595              call xdrfint_(ixdrf,mset(i), iret)
1596            enddo
1597            do i=0,nodes-1
1598              call xdrfint_(ixdrf,i2set(i), iret)
1599            enddo
1600            do il=1,nset
1601              do il1=1,mset(il)
1602                do i=1,nrep
1603                  do j=1,remd_m(i)
1604                    itmp=i_index(i,j,il,il1)
1605                    call xdrfint_(ixdrf,itmp, iret)
1606                  enddo
1607                enddo
1608              enddo
1609            enddo
1610
1611          endif
1612          call xdrfclose_(ixdrf, iret)
1613 #else
1614          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1615          do i=0,nodes-1
1616           call xdrfint(ixdrf, i2rep(i), iret)
1617          enddo
1618          do i=1,remd_m(1)
1619           call xdrfint(ixdrf, ifirst(i), iret)
1620          enddo
1621          do il=1,nodes
1622               do i=0,nupa(0,il)
1623                call xdrfint(ixdrf, nupa(i,il), iret)
1624               enddo
1625
1626               do i=0,ndowna(0,il)
1627                call xdrfint(ixdrf, ndowna(i,il), iret)
1628               enddo
1629          enddo
1630
1631          do il=1,nodes
1632            do j=1,4
1633             call xdrffloat(ixdrf, t_restart1(j,il), iret)
1634            enddo
1635          enddo
1636
1637          do il=0,nodes-1
1638            do i=1,2*nres
1639             do j=1,3
1640              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1641             enddo
1642            enddo
1643          enddo
1644          do il=0,nodes-1
1645            do i=1,2*nres
1646             do j=1,3
1647              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1648             enddo
1649            enddo
1650          enddo
1651
1652
1653              if(usampl.or.homol_nset.gt.1) then
1654               call xdrfint(ixdrf, nset, iret)
1655               do i=1,nset
1656                 call xdrfint(ixdrf,mset(i), iret)
1657               enddo
1658               do i=0,nodes-1
1659                 call xdrfint(ixdrf,i2set(i), iret)
1660               enddo
1661               do il=1,nset
1662                do il1=1,mset(il)
1663                 do i=1,nrep
1664                  do j=1,remd_m(i)
1665                    itmp=i_index(i,j,il,il1)
1666                    call xdrfint(ixdrf,itmp, iret)
1667                  enddo
1668                 enddo
1669                enddo
1670               enddo
1671
1672              endif
1673          call xdrfclose(ixdrf, iret)
1674 #endif
1675        endif
1676       return
1677       end
1678
1679
1680       subroutine write1traj
1681       implicit real*8 (a-h,o-z)
1682       include 'DIMENSIONS'
1683       include 'mpif.h'
1684       include 'COMMON.MD'
1685       include 'COMMON.IOUNITS'
1686       include 'COMMON.REMD'
1687       include 'COMMON.SETUP'
1688       include 'COMMON.CHAIN'
1689       include 'COMMON.SBRIDGE'
1690       include 'COMMON.INTERACT'
1691
1692       real t5_restart1(5)
1693       integer iret,itmp
1694       real xcoord(3,maxres2+2),prec
1695       real r_qfrag(50),r_qpair(100)
1696       real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1697       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1698       real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1699      &     p_uscdiff(100*maxprocs)
1700       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1701       common /przechowalnia/ p_c
1702
1703       call mpi_bcast(ii_write,1,mpi_integer,
1704      &           king,CG_COMM,ierr)
1705
1706 c debugging
1707       print *,'traj1file',me,ii_write,ntwx_cache
1708 c end debugging
1709
1710 #ifdef AIX
1711       if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1712 #else
1713       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1714 #endif
1715       do ii=1,ii_write
1716        t5_restart1(1)=totT_cache(ii)
1717        t5_restart1(2)=EK_cache(ii)
1718        t5_restart1(3)=potE_cache(ii)
1719        t5_restart1(4)=t_bath_cache(ii)
1720        t5_restart1(5)=Uconst_cache(ii)
1721        call mpi_gather(t5_restart1,5,mpi_real,
1722      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1723
1724        call mpi_gather(iset_cache(ii),1,mpi_integer,
1725      &      iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1726
1727           do i=1,nfrag
1728            r_qfrag(i)=qfrag_cache(i,ii)
1729           enddo
1730           do i=1,npair
1731            r_qpair(i)=qpair_cache(i,ii)
1732           enddo
1733           do i=1,nfrag_back
1734            r_utheta(i)=utheta_cache(i,ii)
1735            r_ugamma(i)=ugamma_cache(i,ii)
1736            r_uscdiff(i)=uscdiff_cache(i,ii)
1737           enddo
1738
1739         call mpi_gather(r_qfrag,nfrag,mpi_real,
1740      &           p_qfrag,nfrag,mpi_real,king,
1741      &           CG_COMM,ierr)
1742         call mpi_gather(r_qpair,npair,mpi_real,
1743      &           p_qpair,npair,mpi_real,king,
1744      &           CG_COMM,ierr)
1745         call mpi_gather(r_utheta,nfrag_back,mpi_real,
1746      &           p_utheta,nfrag_back,mpi_real,king,
1747      &           CG_COMM,ierr)
1748         call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1749      &           p_ugamma,nfrag_back,mpi_real,king,
1750      &           CG_COMM,ierr)
1751         call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1752      &           p_uscdiff,nfrag_back,mpi_real,king,
1753      &           CG_COMM,ierr)
1754
1755 #ifdef DEBUG
1756         write (iout,*) "p_qfrag"
1757         do i=1,nodes
1758           write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1759         enddo
1760         write (iout,*) "p_qpair"
1761         do i=1,nodes
1762           write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1763         enddo
1764 ctime        call flush(iout)
1765 #endif
1766         do i=1,nres*2
1767          do j=1,3
1768           r_c(j,i)=c_cache(j,i,ii)
1769          enddo
1770         enddo
1771
1772         call mpi_gather(r_c,3*2*nres,mpi_real,
1773      &           p_c,3*2*nres,mpi_real,king,
1774      &           CG_COMM,ierr)
1775
1776        if(me.eq.king) then
1777 #ifdef AIX
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 #else
1821          do il=1,nodes
1822           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1823           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1824           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1825           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1826           call xdrfint(ixdrf, nss, iret)
1827           do j=1,nss
1828            if (dyn_ss) then
1829             call xdrfint(ixdrf, idssb(j)+nres, iret)
1830             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1831            else
1832             call xdrfint(ixdrf, ihpb(j), iret)
1833             call xdrfint(ixdrf, jhpb(j), iret)
1834            endif
1835           enddo
1836           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1837           call xdrfint(ixdrf, iset_restart1(il), iret)
1838           do i=1,nfrag
1839            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1840           enddo
1841           do i=1,npair
1842            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1843           enddo
1844           do i=1,nfrag_back
1845            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1846            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1847            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1848           enddo
1849           prec=10000.0
1850           do i=1,nres
1851            do j=1,3
1852             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1853            enddo
1854           enddo
1855           do i=nnt,nct
1856            do j=1,3
1857             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1858            enddo
1859           enddo
1860           itmp=nres+nct-nnt+1
1861           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1862          enddo
1863 #endif
1864        endif
1865       enddo
1866 #ifdef AIX
1867       if(me.eq.king) call xdrfclose_(ixdrf, iret)
1868 #else
1869       if(me.eq.king) call xdrfclose(ixdrf, iret)
1870 #endif
1871       do i=1,ntwx_cache-ii_write
1872
1873             totT_cache(i)=totT_cache(ii_write+i)
1874             EK_cache(i)=EK_cache(ii_write+i)
1875             potE_cache(i)=potE_cache(ii_write+i)
1876             t_bath_cache(i)=t_bath_cache(ii_write+i)
1877             Uconst_cache(i)=Uconst_cache(ii_write+i)
1878             iset_cache(i)=iset_cache(ii_write+i)
1879
1880             do ii=1,nfrag
1881              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1882             enddo
1883             do ii=1,npair
1884              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1885             enddo
1886             do ii=1,nfrag_back
1887               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1888               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1889               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1890             enddo
1891
1892             do ii=1,nres*2
1893              do j=1,3
1894               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1895              enddo
1896             enddo
1897       enddo
1898       ntwx_cache=ntwx_cache-ii_write
1899       return
1900       end
1901
1902
1903       subroutine read1restart(i_index)
1904       implicit real*8 (a-h,o-z)
1905       include 'DIMENSIONS'
1906       include 'mpif.h'
1907       include 'COMMON.MD'
1908       include 'COMMON.IOUNITS'
1909       include 'COMMON.REMD'
1910       include 'COMMON.SETUP'
1911       include 'COMMON.CHAIN'
1912       include 'COMMON.SBRIDGE'
1913       include 'COMMON.INTERACT'
1914       include 'COMMON.CONTROL'
1915       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1916      &                 t5_restart1(5)
1917       integer*2 i_index
1918      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1919       common /przechowalnia/ d_restart1
1920       write (*,*) "Processor",me," called read1restart"
1921
1922          if(me.eq.king)then
1923               open(irest2,file=mremd_rst_name,status='unknown')
1924               read(irest2,*,err=334) i
1925               write(iout,*) "Reading old rst in ASCI format"
1926               close(irest2)
1927                call read1restart_old
1928                return
1929  334          continue
1930 #ifdef AIX
1931               call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1932
1933               do i=0,nodes-1
1934                call xdrfint_(ixdrf, i2rep(i), iret)
1935               enddo
1936               do i=1,remd_m(1)
1937                call xdrfint_(ixdrf, ifirst(i), iret)
1938               enddo
1939              do il=1,nodes
1940               call xdrfint_(ixdrf, nupa(0,il), iret)
1941               do i=1,nupa(0,il)
1942                call xdrfint_(ixdrf, nupa(i,il), iret)
1943               enddo
1944
1945               call xdrfint_(ixdrf, ndowna(0,il), iret)
1946               do i=1,ndowna(0,il)
1947                call xdrfint_(ixdrf, ndowna(i,il), iret)
1948               enddo
1949              enddo
1950              do il=1,nodes
1951                do j=1,4
1952                 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1953                enddo
1954              enddo
1955 #else
1956               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1957
1958               do i=0,nodes-1
1959                call xdrfint(ixdrf, i2rep(i), iret)
1960               enddo
1961               do i=1,remd_m(1)
1962                call xdrfint(ixdrf, ifirst(i), iret)
1963               enddo
1964              do il=1,nodes
1965               call xdrfint(ixdrf, nupa(0,il), iret)
1966               do i=1,nupa(0,il)
1967                call xdrfint(ixdrf, nupa(i,il), iret)
1968               enddo
1969
1970               call xdrfint(ixdrf, ndowna(0,il), iret)
1971               do i=1,ndowna(0,il)
1972                call xdrfint(ixdrf, ndowna(i,il), iret)
1973               enddo
1974              enddo
1975              do il=1,nodes
1976                do j=1,4
1977                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1978                enddo
1979              enddo
1980 #endif
1981          endif
1982          call mpi_scatter(t_restart1,5,mpi_real,
1983      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1984          totT=t5_restart1(1)
1985          EK=t5_restart1(2)
1986          potE=t5_restart1(3)
1987          t_bath=t5_restart1(4)
1988
1989          if(me.eq.king)then
1990               do il=0,nodes-1
1991                do i=1,2*nres
1992 c                read(irest2,'(3e15.5)')
1993 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1994             do j=1,3
1995 #ifdef AIX
1996              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1997 #else
1998              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1999 #endif
2000             enddo
2001                enddo
2002               enddo
2003          endif
2004          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2005      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2006
2007          do i=1,2*nres
2008            do j=1,3
2009             d_t(j,i)=r_d(j,i)
2010            enddo
2011          enddo
2012          if(me.eq.king)then
2013               do il=0,nodes-1
2014                do i=1,2*nres
2015 c                read(irest2,'(3e15.5)')
2016 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
2017             do j=1,3
2018 #ifdef AIX
2019              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
2020 #else
2021              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
2022 #endif
2023             enddo
2024                enddo
2025               enddo
2026          endif
2027          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2028      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2029          do i=1,2*nres
2030            do j=1,3
2031             dc(j,i)=r_d(j,i)
2032            enddo
2033          enddo
2034
2035
2036            if(usampl.or.homol_nset.gt.1) then
2037 #ifdef AIX
2038              if(me.eq.king)then
2039               call xdrfint_(ixdrf, nset, iret)
2040               do i=1,nset
2041                 call xdrfint_(ixdrf,mset(i), iret)
2042               enddo
2043               do i=0,nodes-1
2044                 call xdrfint_(ixdrf,i2set(i), iret)
2045               enddo
2046               do il=1,nset
2047                do il1=1,mset(il)
2048                 do i=1,nrep
2049                  do j=1,remd_m(i)
2050                    call xdrfint_(ixdrf,itmp, iret)
2051                    i_index(i,j,il,il1)=itmp
2052                  enddo
2053                 enddo
2054                enddo
2055               enddo
2056              endif
2057 #else
2058              if(me.eq.king)then
2059               call xdrfint(ixdrf, nset, iret)
2060               do i=1,nset
2061                 call xdrfint(ixdrf,mset(i), iret)
2062               enddo
2063               do i=0,nodes-1
2064                 call xdrfint(ixdrf,i2set(i), iret)
2065               enddo
2066               do il=1,nset
2067                do il1=1,mset(il)
2068                 do i=1,nrep
2069                  do j=1,remd_m(i)
2070                    call xdrfint(ixdrf,itmp, iret)
2071                    i_index(i,j,il,il1)=itmp
2072                  enddo
2073                 enddo
2074                enddo
2075               enddo
2076              endif
2077 #endif
2078 c Corrected AL 8/19/2014: each processor needs whole iset array not only its
2079 c own element
2080 c              call mpi_scatter(i2set,1,mpi_integer,
2081 c     &           iset,1,mpi_integer,king,
2082 c     &           CG_COMM,ierr)
2083               call mpi_bcast(i2set(0),nodes,mpi_integer,king,
2084      &         CG_COMM,ierr)
2085               iset=i2set(me)
2086 c broadcast iset to slaves
2087               if (nfgtasks.gt.1) then
2088                call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
2089                call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
2090               endif
2091            endif
2092
2093
2094         if(me.eq.king) close(irest2)
2095         return
2096         end
2097
2098       subroutine read1restart_old
2099       implicit real*8 (a-h,o-z)
2100       include 'DIMENSIONS'
2101       include 'mpif.h'
2102       include 'COMMON.MD'
2103       include 'COMMON.IOUNITS'
2104       include 'COMMON.REMD'
2105       include 'COMMON.SETUP'
2106       include 'COMMON.CHAIN'
2107       include 'COMMON.SBRIDGE'
2108       include 'COMMON.INTERACT'
2109       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
2110      &                 t5_restart1(5)
2111       common /przechowalnia/ d_restart1
2112          if(me.eq.king)then
2113              open(irest2,file=mremd_rst_name,status='unknown')
2114              read (irest2,*) (i2rep(i),i=0,nodes-1)
2115              read (irest2,*) (ifirst(i),i=1,remd_m(1))
2116              do il=1,nodes
2117               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2118               read (irest2,*) ndowna(0,il),
2119      &                    (ndowna(i,il),i=1,ndowna(0,il))
2120              enddo
2121              do il=1,nodes
2122                read(irest2,*) (t_restart1(j,il),j=1,4)
2123              enddo
2124          endif
2125          call mpi_scatter(t_restart1,5,mpi_real,
2126      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2127          totT=t5_restart1(1)
2128          EK=t5_restart1(2)
2129          potE=t5_restart1(3)
2130          t_bath=t5_restart1(4)
2131
2132          if(me.eq.king)then
2133               do il=0,nodes-1
2134                do i=1,2*nres
2135                 read(irest2,'(3e15.5)')
2136      &                (d_restart1(j,i+2*nres*il),j=1,3)
2137                enddo
2138               enddo
2139          endif
2140          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2141      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2142
2143          do i=1,2*nres
2144            do j=1,3
2145             d_t(j,i)=r_d(j,i)
2146            enddo
2147          enddo
2148          if(me.eq.king)then
2149               do il=0,nodes-1
2150                do i=1,2*nres
2151                 read(irest2,'(3e15.5)')
2152      &                (d_restart1(j,i+2*nres*il),j=1,3)
2153                enddo
2154               enddo
2155          endif
2156          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2157      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2158          do i=1,2*nres
2159            do j=1,3
2160             dc(j,i)=r_d(j,i)
2161            enddo
2162          enddo
2163         if(me.eq.king) close(irest2)
2164         return
2165         end
2166 c-------------------------------------------------------------------
2167         subroutine set_hweights(iiset)
2168         implicit real*8 (a-h,o-z)
2169         integer i
2170         include 'DIMENSIONS'
2171         include 'COMMON.FFIELD'
2172         include 'COMMON.REMD'
2173
2174          do i=1,n_ene
2175           weights(i)=hweights(iiset,i)
2176          enddo
2177
2178          wsc    =weights(1)
2179          wscp   =weights(2)
2180          welec  =weights(3)
2181          wcorr  =weights(4)
2182          wcorr5 =weights(5)
2183          wcorr6 =weights(6)
2184          wel_loc=weights(7)
2185          wturn3 =weights(8)
2186          wturn4 =weights(9)
2187          wturn6 =weights(10)
2188          wang   =weights(11)
2189          wscloc =weights(12)
2190          wtor   =weights(13)
2191          wtor_d =weights(14)
2192          wstrain=weights(15)
2193          wvdwpp =weights(16)
2194          wbond  =weights(17)
2195          scal14 =weights(18)
2196          wsccor =weights(21)
2197
2198         return
2199         end
2200 #endif