source/unres/src_MD/MREMD.F

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