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 #ifdef DEBUG
 786             if(lmuca) then
 787 co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
 788              do i=1,nodes
 789                write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
 790      &            elowi(i),ehighi(i)
 791              enddo
 792             else
 793               write(iout,*) 'REMD exchange temp,ene'
 794               do i=1,nodes
 795                 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
 796                 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
 797               enddo
 798             endif
 799 #endif
 800 c-------------------------------------
 801            IF(.not.usampl.and.hremd.eq.0) THEN
 802 #ifdef DEBUG
 803             write (iout,*) "Enter exchnge, remd_m",remd_m(1),
 804      &        " nodes",nodes
 805 ctime            call flush(iout)
 806             write (iout,*) "remd_m(1)",remd_m(1)
 807 #endif
 808             do irr=1,remd_m(1)
 809                i=ifirst(iran_num(1,remd_m(1)))
 810 #ifdef DEBUG
 811              write (iout,*) "i",i
 812 #endif
 813 ctime             call flush(iout)
 814
 815              do ii=1,nodes-1
 816
 817 #ifdef DEBUG
 818               write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
 819 #endif
 820              if(i.gt.0.and.nupa(0,i).gt.0) then
 821               iex=i
 822 c              if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
 823 c                write (iout,*)
 824 c     &  "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
 825 c                call flush(iout)
 826 c                call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
 827 c              endif
 828 c              do while (iex.eq.i)
 829 c                write (iout,*) "upper",nupa(int(nupa(0,i)),i)
 830                 iex=nupa(iran_num(1,int(nupa(0,i))),i)
 831 c              enddo
 832 c              write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
 833               if (lmuca) then
 834                call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
 835               else
 836 c Swap temperatures between conformations i and iex with recalculating the free energies
 837 c following temperature changes.
 838                ene_iex_iex=remd_ene(0,iex)
 839                ene_i_i=remd_ene(0,i)
 840 c               write (iout,*) "i",i," ene_i_i",ene_i_i,
 841 c     &          " iex",iex," ene_iex_iex",ene_iex_iex
 842 c               write (iout,*) "rescaling weights with temperature",
 843 c     &          remd_t_bath(i)
 844 c               call flush(iout)
 845                call rescale_weights(remd_t_bath(i))
 846
 847 c               write (iout,*) "0,iex",remd_t_bath(i)
 848 c               call enerprint(remd_ene(0,iex))
 849
 850                call sum_energy(remd_ene(0,iex),.false.)
 851                ene_iex_i=remd_ene(0,iex)
 852 c               write (iout,*) "ene_iex_i",remd_ene(0,iex)
 853
 854 c               write (iout,*) "0,i",remd_t_bath(i)
 855 c               call enerprint(remd_ene(0,i))
 856
 857                call sum_energy(remd_ene(0,i),.false.)
 858 c               write (iout,*) "ene_i_i",remd_ene(0,i)
 859 c               call flush(iout)
 860 c               write (iout,*) "rescaling weights with temperature",
 861 c     &          remd_t_bath(iex)
 862                if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
 863                 write (iout,*) "ERROR: inconsistent energies:",i,
 864      &            ene_i_i,remd_ene(0,i)
 865                endif
 866                call rescale_weights(remd_t_bath(iex))
 867
 868 c               write (iout,*) "0,i",remd_t_bath(iex)
 869 c               call enerprint(remd_ene(0,i))
 870
 871                call sum_energy(remd_ene(0,i),.false.)
 872 c               write (iout,*) "ene_i_iex",remd_ene(0,i)
 873 c               call flush(iout)
 874                ene_i_iex=remd_ene(0,i)
 875
 876 c               write (iout,*) "0,iex",remd_t_bath(iex)
 877 c               call enerprint(remd_ene(0,iex))
 878
 879                call sum_energy(remd_ene(0,iex),.false.)
 880                if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
 881                 write (iout,*) "ERROR: inconsistent energies:",iex,
 882      &            ene_iex_iex,remd_ene(0,iex)
 883                endif
 884 c               write (iout,*) "ene_iex_iex",remd_ene(0,iex)
 885 c               write (iout,*) "i",i," iex",iex
 886 c               write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
 887 c     &           " ene_i_iex",ene_i_iex,
 888 c     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
 889 c               call flush(iout)
 890                delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
 891      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
 892                delta=-delta
 893 c               write(iout,*) 'delta',delta
 894 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
 895 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
 896 c     &              (remd_t_bath(i)*remd_t_bath(iex))
 897               endif
 898               if (delta .gt. 50.0d0) then
 899                 delta=0.0d0
 900               else
 901 #ifdef OSF
 902                 if(isnan(delta))then
 903                   delta=0.0d0
 904                 else if (delta.lt.-50.0d0) then
 905                   delta=dexp(50.0d0)
 906                 else
 907                   delta=dexp(-delta)
 908                 endif
 909 #else
 910                 delta=dexp(-delta)
 911 #endif
 912               endif
 913               iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
 914               xxx=ran_number(0.0d0,1.0d0)
 915 c              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
 916 c              call flush(iout)
 917               if (delta .gt. xxx) then
 918                 tmp=remd_t_bath(i)
 919                 remd_t_bath(i)=remd_t_bath(iex)
 920                 remd_t_bath(iex)=tmp
 921                 remd_ene(0,i)=ene_i_iex
 922                 remd_ene(0,iex)=ene_iex_i
 923                 if(lmuca) then
 924                   tmp=elowi(i)
 925                   elowi(i)=elowi(iex)
 926                   elowi(iex)=tmp
 927                   tmp=ehighi(i)
 928                   ehighi(i)=ehighi(iex)
 929                   ehighi(iex)=tmp
 930                 endif
 931
 932
 933                 do k=0,nodes
 934                   itmp=nupa(k,i)
 935                   nupa(k,i)=nupa(k,iex)
 936                   nupa(k,iex)=itmp
 937                   itmp=ndowna(k,i)
 938                   ndowna(k,i)=ndowna(k,iex)
 939                   ndowna(k,iex)=itmp
 940                 enddo
 941                 do il=1,nodes
 942                  if (ifirst(il).eq.i) ifirst(il)=iex
 943                  do k=1,nupa(0,il)
 944                   if (nupa(k,il).eq.i) then
 945                      nupa(k,il)=iex
 946                   elseif (nupa(k,il).eq.iex) then
 947                      nupa(k,il)=i
 948                   endif
 949                  enddo
 950                  do k=1,ndowna(0,il)
 951                   if (ndowna(k,il).eq.i) then
 952                      ndowna(k,il)=iex
 953                   elseif (ndowna(k,il).eq.iex) then
 954                      ndowna(k,il)=i
 955                   endif
 956                  enddo
 957                 enddo
 958
 959                 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
 960                 itmp=i2rep(i-1)
 961                 i2rep(i-1)=i2rep(iex-1)
 962                 i2rep(iex-1)=itmp
 963
 964 c                write(iout,*) 'exchange',i,iex
 965 c                write (iout,'(a8,100i4)') "@ ifirst",
 966 c     &                    (ifirst(k),k=1,remd_m(1))
 967 c                do il=1,nodes
 968 c                 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
 969 c     &                    (nupa(k,il),k=1,nupa(0,il))
 970 c                 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
 971 c     &                    (ndowna(k,il),k=1,ndowna(0,il))
 972 c                enddo
 973 c                call flush(iout)
 974
 975               else
 976                remd_ene(0,iex)=ene_iex_iex
 977                remd_ene(0,i)=ene_i_i
 978                i=iex
 979               endif
 980             endif
 981            enddo
 982            enddo
 983 cd           write (iout,*) "exchange completed"
 984 cd           call flush(iout)
 985         ELSEIF (usampl) THEN
 986           do ii=1,nodes
 987 cd            write(iout,*) "########",ii
 988
 989             i_temp=iran_num(1,nrep)
 990             i_mult=iran_num(1,remd_m(i_temp))
 991             i_iset=iran_num(1,nset)
 992             i_mset=iran_num(1,mset(i_iset))
 993             i=i_index(i_temp,i_mult,i_iset,i_mset)
 994
 995 cd            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
 996
 997             if(t_exchange_only)then
 998              i_dir=1
 999             else
1000              i_dir=iran_num(1,3)
1001             endif
1002 cd            write(iout,*) "i_dir=",i_dir
1003
1004             if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1005
1006                i_temp1=i_temp+1
1007                i_mult1=iran_num(1,remd_m(i_temp1))
1008                i_iset1=i_iset
1009                i_mset1=iran_num(1,mset(i_iset1))
1010                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1011
1012             elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1013
1014                i_temp1=i_temp
1015                i_mult1=iran_num(1,remd_m(i_temp1))
1016                i_iset1=i_iset+1
1017                i_mset1=iran_num(1,mset(i_iset1))
1018                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1019                econstr_temp_i=remd_ene(20,i)
1020                econstr_temp_iex=remd_ene(20,iex)
1021                remd_ene(20,i)=remd_ene(n_ene+3,i)
1022                remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1023
1024             elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1025
1026                i_temp1=i_temp+1
1027                i_mult1=iran_num(1,remd_m(i_temp1))
1028                i_iset1=i_iset+1
1029                i_mset1=iran_num(1,mset(i_iset1))
1030                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1031                econstr_temp_i=remd_ene(20,i)
1032                econstr_temp_iex=remd_ene(20,iex)
1033                remd_ene(20,i)=remd_ene(n_ene+3,i)
1034                remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1035
1036             else
1037                goto 444
1038             endif
1039
1040 cd            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1041 ctime            call flush(iout)
1042
1043 c Swap temperatures between conformations i and iex with recalculating the free energies
1044 c following temperature changes.
1045               ene_iex_iex=remd_ene(0,iex)
1046               ene_i_i=remd_ene(0,i)
1047 co              write (iout,*) "rescaling weights with temperature",
1048 co     &          remd_t_bath(i)
1049               call rescale_weights(remd_t_bath(i))
1050
1051               call sum_energy(remd_ene(0,iex),.false.)
1052               ene_iex_i=remd_ene(0,iex)
1053 cd              write (iout,*) "ene_iex_i",remd_ene(0,iex)
1054 c              call sum_energy(remd_ene(0,i),.false.)
1055 cd              write (iout,*) "ene_i_i",remd_ene(0,i)
1056 c              write (iout,*) "rescaling weights with temperature",
1057 c     &          remd_t_bath(iex)
1058 c              if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1059 c                write (iout,*) "ERROR: inconsistent energies:",i,
1060 c     &            ene_i_i,remd_ene(0,i)
1061 c              endif
1062               call rescale_weights(remd_t_bath(iex))
1063               call sum_energy(remd_ene(0,i),.false.)
1064 cd              write (iout,*) "ene_i_iex",remd_ene(0,i)
1065               ene_i_iex=remd_ene(0,i)
1066 c              call sum_energy(remd_ene(0,iex),.false.)
1067 c              if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1068 c                write (iout,*) "ERROR: inconsistent energies:",iex,
1069 c     &            ene_iex_iex,remd_ene(0,iex)
1070 c              endif
1071 cd              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1072 c              write (iout,*) "i",i," iex",iex
1073 cd              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1074 cd     &           " ene_i_iex",ene_i_iex,
1075 cd     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1076               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1077      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1078               delta=-delta
1079 cd              write(iout,*) 'delta',delta
1080 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
1081 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
1082 c     &              (remd_t_bath(i)*remd_t_bath(iex))
1083               if (delta .gt. 50.0d0) then
1084                 delta=0.0d0
1085               else
1086                 delta=dexp(-delta)
1087               endif
1088               if (i_dir.eq.1.or.i_dir.eq.3)
1089      &         iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1090               if (i_dir.eq.2.or.i_dir.eq.3)
1091      &          iremd_tot_usa(int(i2set(i-1)))=
1092      &                 iremd_tot_usa(int(i2set(i-1)))+1
1093               xxx=ran_number(0.0d0,1.0d0)
1094 cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1095               if (delta .gt. xxx) then
1096                 tmp=remd_t_bath(i)
1097                 remd_t_bath(i)=remd_t_bath(iex)
1098                 remd_t_bath(iex)=tmp
1099
1100                 itmp=iremd_iset(i)
1101                 iremd_iset(i)=iremd_iset(iex)
1102                 iremd_iset(iex)=itmp
1103
1104                 remd_ene(0,i)=ene_i_iex
1105                 remd_ene(0,iex)=ene_iex_i
1106
1107                 if (i_dir.eq.1.or.i_dir.eq.3)
1108      &           iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1109
1110                 itmp=i2rep(i-1)
1111                 i2rep(i-1)=i2rep(iex-1)
1112                 i2rep(iex-1)=itmp
1113
1114                 if (i_dir.eq.2.or.i_dir.eq.3)
1115      &           iremd_acc_usa(int(i2set(i-1)))=
1116      &                 iremd_acc_usa(int(i2set(i-1)))+1
1117
1118                 itmp=i2set(i-1)
1119                 i2set(i-1)=i2set(iex-1)
1120                 i2set(iex-1)=itmp
1121
1122                 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1123                 i_index(i_temp,i_mult,i_iset,i_mset)=
1124      &                i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1125                 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1126
1127               else
1128                remd_ene(0,iex)=ene_iex_iex
1129                remd_ene(0,i)=ene_i_i
1130                remd_ene(20,iex)=econstr_temp_iex
1131                remd_ene(20,i)=econstr_temp_i
1132               endif
1133
1134 cd      do il=1,nset
1135 cd       do il1=1,mset(il)
1136 cd        do i=1,nrep
1137 cd         do j=1,remd_m(i)
1138 cd          write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1139 cd         enddo
1140 cd        enddo
1141 cd       enddo
1142 cd      enddo
1143
1144  444      continue
1145
1146           enddo
1147
1148         ELSEIF (hremd.gt.0) THEN
1149           do ii=1,nodes
1150 cd            write(iout,*) "########",ii
1151
1152             i_temp=iran_num(1,nrep)
1153             i_mult=iran_num(1,remd_m(i_temp))
1154             i_iset=iran_num(1,nset)
1155             i_mset=1
1156             i=i_index(i_temp,i_mult,i_iset,i_mset)
1157
1158 cd            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1159
1160             if(t_exchange_only)then
1161              i_dir=1
1162             else
1163              i_dir=iran_num(1,3)
1164             endif
1165
1166 cd            write(iout,*) "i_dir=",i_dir
1167
1168             if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1169
1170                i_temp1=i_temp+1
1171                i_mult1=iran_num(1,remd_m(i_temp1))
1172                i_iset1=i_iset
1173                i_mset1=1
1174                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1175
1176             elseif(i_dir.eq.2)then
1177
1178                i_temp1=i_temp
1179                i_mult1=iran_num(1,remd_m(i_temp1))
1180                i_iset1=iran_num(1,hremd)
1181                do while(i_iset1.eq.i_iset)
1182                  i_iset1=iran_num(1,hremd)
1183                enddo
1184                i_mset1=1
1185                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1186
1187             elseif(remd_m(i_temp+1).gt.0)then
1188
1189                i_temp1=i_temp+1
1190                i_mult1=iran_num(1,remd_m(i_temp1))
1191                i_iset1=iran_num(1,hremd)
1192                do while(i_iset1.eq.i_iset)
1193                  i_iset1=iran_num(1,hremd)
1194                enddo
1195                i_mset1=1
1196                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1197
1198             else
1199                goto 445
1200             endif
1201
1202 cd            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1203 ctime            call flush(iout)
1204
1205 c Swap temperatures between conformations i and iex with recalculating the free energies
1206 c following temperature changes.
1207               ene_iex_iex=remd_ene(0,iex)
1208               ene_i_i=remd_ene(0,i)
1209
1210               call set_hweights(i_iset)
1211               call rescale_weights(remd_t_bath(i))
1212               call sum_energy(remd_ene(0,iex),.false.)
1213               ene_iex_i=remd_ene(0,iex)
1214
1215               call set_hweights(i_iset1)
1216               call rescale_weights(remd_t_bath(iex))
1217               call sum_energy(remd_ene(0,i),.false.)
1218               ene_i_iex=remd_ene(0,i)
1219
1220 cd              write(iout,*)  ene_iex_iex,ene_i_i,ene_iex_i,ene_i_iex
1221
1222               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1223      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1224               delta=-delta
1225
1226               if (delta .gt. 50.0d0) then
1227                 delta=0.0d0
1228               else
1229                 delta=dexp(-delta)
1230               endif
1231
1232               if (i_dir.eq.1.or.i_dir.eq.3)
1233      &         iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1234               if (i_dir.eq.2.or.i_dir.eq.3)
1235      &          iremd_tot_usa(int(i2set(i-1)))=
1236      &                 iremd_tot_usa(int(i2set(i-1)))+1
1237               xxx=ran_number(0.0d0,1.0d0)
1238 cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1239               if (delta .gt. xxx) then
1240
1241 cd                write (iout,*) "exchange"
1242                 tmp=remd_t_bath(i)
1243                 remd_t_bath(i)=remd_t_bath(iex)
1244                 remd_t_bath(iex)=tmp
1245
1246                 itmp=iremd_iset(i)
1247                 iremd_iset(i)=iremd_iset(iex)
1248                 iremd_iset(iex)=itmp
1249
1250                 remd_ene(0,i)=ene_i_iex
1251                 remd_ene(0,iex)=ene_iex_i
1252
1253                 if (i_dir.eq.1.or.i_dir.eq.3)
1254      &           iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1255
1256                 itmp=i2rep(i-1)
1257                 i2rep(i-1)=i2rep(iex-1)
1258                 i2rep(iex-1)=itmp
1259
1260                 if (i_dir.eq.2.or.i_dir.eq.3)
1261      &           iremd_acc_usa(int(i2set(i-1)))=
1262      &                 iremd_acc_usa(int(i2set(i-1)))+1
1263
1264                 itmp=i2set(i-1)
1265                 i2set(i-1)=i2set(iex-1)
1266                 i2set(iex-1)=itmp
1267
1268                 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1269                 i_index(i_temp,i_mult,i_iset,i_mset)=
1270      &                i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1271                 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1272
1273 cd       do il=1,nset
1274 cd        do il1=1,mset(il)
1275 cd         do i=1,nrep
1276 cd          do j=1,remd_m(i)
1277 cd            write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1278 cd          enddo
1279 cd         enddo
1280 cd        enddo
1281 cd       enddo
1282
1283               else
1284                remd_ene(0,iex)=ene_iex_iex
1285                remd_ene(0,i)=ene_i_i
1286               endif
1287
1288
1289
1290  445      continue
1291
1292           enddo
1293
1294         ENDIF
1295
1296 c-------------------------------------
1297              write (iout,*) "NREP",nrep
1298              do i=1,nrep
1299               if(iremd_tot(i).ne.0)
1300      &          write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1301      &           ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1302              enddo
1303
1304              if(usampl) then
1305               do i=1,nset
1306                if(iremd_tot_usa(i).ne.0)
1307      &           write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1308      &         iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1309               enddo
1310              endif
1311
1312              if(hremd.gt.0) then
1313               do i=1,nset
1314                if(iremd_tot_usa(i).ne.0)
1315      &           write(iout,'(a10,i4,f12.5,i8)') 'ACC_hremd',i,
1316      &         iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1317               enddo
1318              endif
1319
1320
1321 ctime             call flush(iout)
1322
1323 cd              write (iout,'(a6,100i4)') "ifirst",
1324 cd     &                    (ifirst(i),i=1,remd_m(1))
1325 cd              do il=1,nodes
1326 cd               write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1327 cd     &                    (nupa(i,il),i=1,nupa(0,il))
1328 cd               write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1329 cd     &                    (ndowna(i,il),i=1,ndowna(0,il))
1330 cd              enddo
1331             endif
1332
1333          time06=MPI_WTIME()
1334 cd         write (iout,*) "Before scatter"
1335 cd         call flush(iout)
1336          call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1337      &           t_bath,1,mpi_double_precision,king,
1338      &           CG_COMM,ierr)
1339 cd         write (iout,*) "After scatter"
1340 cd         call flush(iout)
1341          if(usampl.or.hremd.gt.0)
1342      &    call mpi_scatter(iremd_iset,1,mpi_integer,
1343      &           iset,1,mpi_integer,king,
1344      &           CG_COMM,ierr)
1345
1346          time07=MPI_WTIME()
1347           if (me.eq.king .or. .not. out1file) then
1348             write(iout,*) 'REMD scatter time=',time07-time06
1349           endif
1350
1351          if(lmuca) then
1352            call mpi_scatter(elowi,1,mpi_double_precision,
1353      &           elow,1,mpi_double_precision,king,
1354      &           CG_COMM,ierr)
1355            call mpi_scatter(ehighi,1,mpi_double_precision,
1356      &           ehigh,1,mpi_double_precision,king,
1357      &           CG_COMM,ierr)
1358          endif
1359
1360          if(hremd.gt.0) call set_hweights(iset)
1361          call rescale_weights(t_bath)
1362 co         write (iout,*) "Processor",me,
1363 co     &    " rescaling weights with temperature",t_bath
1364
1365          stdfp=dsqrt(2*Rb*t_bath/d_time)
1366          do i=1,ntyp
1367            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1368          enddo
1369          if (lang.gt.0) then
1370            do i=nnt,nct-1
1371             stdforcp(i)=stdforcp(i)*sqrt(t_bath/t_bath_old)
1372            enddo
1373            do i=nnt,nct
1374             stdforcsc(i)=stdforcsc(i)*sqrt(t_bath/t_bath_old)
1375            enddo
1376          endif
1377 cde         write(iout,*) 'REMD after',me,t_bath
1378            time08=MPI_WTIME()
1379            if (me.eq.king .or. .not. out1file) then
1380             write(iout,*) 'REMD exchange time=',time08-time02
1381 ctime            call flush(iout)
1382            endif
1383         endif
1384       enddo
1385
1386       if (restart1file) then
1387           if (me.eq.king .or. .not. out1file)
1388      &      write(iout,*) 'writing restart at the end of run'
1389            call write1rst(i_index)
1390       endif
1391
1392       if (traj1file) call write1traj
1393 cd debugging
1394 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
1395 cdeb     &             icache_all,1,mpi_integer,king,
1396 cdeb     &             CG_COMM,ierr)
1397 cdeb            write(iout,'(a40,8000i8)')
1398 cdeb     &             '  ntwx_cache after traj1file at the end',
1399 cdeb     &             (icache_all(i),i=1,nodes)
1400 cd end
1401
1402
1403 #ifdef MPI
1404       t_MD=MPI_Wtime()-tt0
1405 #else
1406       t_MD=tcpu()-tt0
1407 #endif
1408       if (me.eq.king .or. .not. out1file) then
1409        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1410      &  '  Timing  ',
1411      & 'MD calculations setup:',t_MDsetup,
1412      & 'Energy & gradient evaluation:',t_enegrad,
1413      & 'Stochastic MD setup:',t_langsetup,
1414      & 'Stochastic MD step setup:',t_sdsetup,
1415      & 'MD steps:',t_MD
1416        write (iout,'(/28(1h=),a25,27(1h=))')
1417      & '  End of MD calculation  '
1418        if(hmc.gt.0) write (iout,*) 'HMC acceptance ratio',
1419      &         n_timestep*1.0d0/hmc/hmc_acc
1420       endif
1421       return
1422       end
1423
1424 c-----------------------------------------------------------------------
1425       subroutine write1rst(i_index)
1426       implicit real*8 (a-h,o-z)
1427       include 'DIMENSIONS'
1428       include 'mpif.h'
1429       include 'COMMON.MD'
1430       include 'COMMON.IOUNITS'
1431       include 'COMMON.REMD'
1432       include 'COMMON.SETUP'
1433       include 'COMMON.CHAIN'
1434       include 'COMMON.SBRIDGE'
1435       include 'COMMON.INTERACT'
1436
1437       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1438      &     d_restart2(3,2*maxres*maxprocs)
1439       real t5_restart1(5)
1440       integer iret,itmp
1441       integer*2 i_index
1442      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1443        common /przechowalnia/ d_restart1,d_restart2
1444
1445        t5_restart1(1)=totT
1446        t5_restart1(2)=EK
1447        t5_restart1(3)=potE
1448        t5_restart1(4)=t_bath
1449        t5_restart1(5)=Uconst
1450
1451        call mpi_gather(t5_restart1,5,mpi_real,
1452      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1453
1454
1455        do i=1,2*nres
1456          do j=1,3
1457            r_d(j,i)=d_t(j,i)
1458          enddo
1459        enddo
1460        call mpi_gather(r_d,3*2*nres,mpi_real,
1461      &           d_restart1,3*2*nres,mpi_real,king,
1462      &           CG_COMM,ierr)
1463
1464
1465        do i=1,2*nres
1466          do j=1,3
1467            r_d(j,i)=dc(j,i)
1468          enddo
1469        enddo
1470        call mpi_gather(r_d,3*2*nres,mpi_real,
1471      &           d_restart2,3*2*nres,mpi_real,king,
1472      &           CG_COMM,ierr)
1473
1474        if(me.eq.king) then
1475 #ifdef AIX
1476          call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1477          do i=0,nodes-1
1478           call xdrfint_(ixdrf, i2rep(i), iret)
1479          enddo
1480          do i=1,remd_m(1)
1481           call xdrfint_(ixdrf, ifirst(i), iret)
1482          enddo
1483          do il=1,nodes
1484               do i=0,nupa(0,il)
1485                call xdrfint_(ixdrf, nupa(i,il), iret)
1486               enddo
1487
1488               do i=0,ndowna(0,il)
1489                call xdrfint_(ixdrf, ndowna(i,il), iret)
1490               enddo
1491          enddo
1492
1493          do il=1,nodes
1494            do j=1,4
1495             call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1496            enddo
1497          enddo
1498
1499          do il=0,nodes-1
1500            do i=1,2*nres
1501             do j=1,3
1502              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1503             enddo
1504            enddo
1505          enddo
1506          do il=0,nodes-1
1507            do i=1,2*nres
1508             do j=1,3
1509              call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1510             enddo
1511            enddo
1512          enddo
1513
1514          if(usampl) then
1515            call xdrfint_(ixdrf, nset, iret)
1516            do i=1,nset
1517              call xdrfint_(ixdrf,mset(i), iret)
1518            enddo
1519            do i=0,nodes-1
1520              call xdrfint_(ixdrf,i2set(i), iret)
1521            enddo
1522            do il=1,nset
1523              do il1=1,mset(il)
1524                do i=1,nrep
1525                  do j=1,remd_m(i)
1526                    itmp=i_index(i,j,il,il1)
1527                    call xdrfint_(ixdrf,itmp, iret)
1528                  enddo
1529                enddo
1530              enddo
1531            enddo
1532
1533          endif
1534          call xdrfclose_(ixdrf, iret)
1535 #else
1536          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1537          do i=0,nodes-1
1538           call xdrfint(ixdrf, i2rep(i), iret)
1539          enddo
1540          do i=1,remd_m(1)
1541           call xdrfint(ixdrf, ifirst(i), iret)
1542          enddo
1543          do il=1,nodes
1544               do i=0,nupa(0,il)
1545                call xdrfint(ixdrf, nupa(i,il), iret)
1546               enddo
1547
1548               do i=0,ndowna(0,il)
1549                call xdrfint(ixdrf, ndowna(i,il), iret)
1550               enddo
1551          enddo
1552
1553          do il=1,nodes
1554            do j=1,4
1555             call xdrffloat(ixdrf, t_restart1(j,il), iret)
1556            enddo
1557          enddo
1558
1559          do il=0,nodes-1
1560            do i=1,2*nres
1561             do j=1,3
1562              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1563             enddo
1564            enddo
1565          enddo
1566          do il=0,nodes-1
1567            do i=1,2*nres
1568             do j=1,3
1569              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1570             enddo
1571            enddo
1572          enddo
1573
1574
1575              if(usampl) then
1576               call xdrfint(ixdrf, nset, iret)
1577               do i=1,nset
1578                 call xdrfint(ixdrf,mset(i), iret)
1579               enddo
1580               do i=0,nodes-1
1581                 call xdrfint(ixdrf,i2set(i), iret)
1582               enddo
1583               do il=1,nset
1584                do il1=1,mset(il)
1585                 do i=1,nrep
1586                  do j=1,remd_m(i)
1587                    itmp=i_index(i,j,il,il1)
1588                    call xdrfint(ixdrf,itmp, iret)
1589                  enddo
1590                 enddo
1591                enddo
1592               enddo
1593
1594              endif
1595          call xdrfclose(ixdrf, iret)
1596 #endif
1597        endif
1598       return
1599       end
1600
1601
1602       subroutine write1traj
1603       implicit real*8 (a-h,o-z)
1604       include 'DIMENSIONS'
1605       include 'mpif.h'
1606       include 'COMMON.MD'
1607       include 'COMMON.IOUNITS'
1608       include 'COMMON.REMD'
1609       include 'COMMON.SETUP'
1610       include 'COMMON.CHAIN'
1611       include 'COMMON.SBRIDGE'
1612       include 'COMMON.INTERACT'
1613
1614       real t5_restart1(5)
1615       integer iret,itmp
1616       real xcoord(3,maxres2+2),prec
1617       real r_qfrag(50),r_qpair(100)
1618       real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1619       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1620       real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1621      &     p_uscdiff(100*maxprocs)
1622       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1623       common /przechowalnia/ p_c
1624
1625       call mpi_bcast(ii_write,1,mpi_integer,
1626      &           king,CG_COMM,ierr)
1627
1628 c debugging
1629       print *,'traj1file',me,ii_write,ntwx_cache
1630 c end debugging
1631
1632 #ifdef AIX
1633       if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1634 #else
1635       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1636 #endif
1637       do ii=1,ii_write
1638        t5_restart1(1)=totT_cache(ii)
1639        t5_restart1(2)=EK_cache(ii)
1640        t5_restart1(3)=potE_cache(ii)
1641        t5_restart1(4)=t_bath_cache(ii)
1642        t5_restart1(5)=Uconst_cache(ii)
1643        call mpi_gather(t5_restart1,5,mpi_real,
1644      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1645
1646        call mpi_gather(iset_cache(ii),1,mpi_integer,
1647      &      iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1648
1649           do i=1,nfrag
1650            r_qfrag(i)=qfrag_cache(i,ii)
1651           enddo
1652           do i=1,npair
1653            r_qpair(i)=qpair_cache(i,ii)
1654           enddo
1655           do i=1,nfrag_back
1656            r_utheta(i)=utheta_cache(i,ii)
1657            r_ugamma(i)=ugamma_cache(i,ii)
1658            r_uscdiff(i)=uscdiff_cache(i,ii)
1659           enddo
1660
1661         call mpi_gather(r_qfrag,nfrag,mpi_real,
1662      &           p_qfrag,nfrag,mpi_real,king,
1663      &           CG_COMM,ierr)
1664         call mpi_gather(r_qpair,npair,mpi_real,
1665      &           p_qpair,npair,mpi_real,king,
1666      &           CG_COMM,ierr)
1667         call mpi_gather(r_utheta,nfrag_back,mpi_real,
1668      &           p_utheta,nfrag_back,mpi_real,king,
1669      &           CG_COMM,ierr)
1670         call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1671      &           p_ugamma,nfrag_back,mpi_real,king,
1672      &           CG_COMM,ierr)
1673         call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1674      &           p_uscdiff,nfrag_back,mpi_real,king,
1675      &           CG_COMM,ierr)
1676
1677 #ifdef DEBUG
1678         write (iout,*) "p_qfrag"
1679         do i=1,nodes
1680           write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1681         enddo
1682         write (iout,*) "p_qpair"
1683         do i=1,nodes
1684           write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1685         enddo
1686 ctime        call flush(iout)
1687 #endif
1688         do i=1,nres*2
1689          do j=1,3
1690           r_c(j,i)=c_cache(j,i,ii)
1691          enddo
1692         enddo
1693
1694         call mpi_gather(r_c,3*2*nres,mpi_real,
1695      &           p_c,3*2*nres,mpi_real,king,
1696      &           CG_COMM,ierr)
1697
1698        if(me.eq.king) then
1699 #ifdef AIX
1700          do il=1,nodes
1701           call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1702           call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1703           call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1704           call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1705           call xdrfint_(ixdrf, nss, iret)
1706           do j=1,nss
1707           do j=1,nss
1708            if (dyn_ss) then
1709             call xdrfint(ixdrf, idssb(j)+nres, iret)
1710             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1711            else
1712             call xdrfint_(ixdrf, ihpb(j), iret)
1713             call xdrfint_(ixdrf, jhpb(j), iret)
1714            endif
1715           enddo
1716           call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1717           call xdrfint_(ixdrf, iset_restart1(il), iret)
1718           do i=1,nfrag
1719            call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1720           enddo
1721           do i=1,npair
1722            call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1723           enddo
1724           do i=1,nfrag_back
1725            call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1726            call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1727            call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1728           enddo
1729           prec=10000.0
1730           do i=1,nres
1731            do j=1,3
1732             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1733            enddo
1734           enddo
1735           do i=nnt,nct
1736            do j=1,3
1737             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1738            enddo
1739           enddo
1740           itmp=nres+nct-nnt+1
1741           call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1742          enddo
1743 #else
1744          do il=1,nodes
1745           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1746           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1747           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1748           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1749           call xdrfint(ixdrf, nss, iret)
1750           do j=1,nss
1751            if (dyn_ss) then
1752             call xdrfint(ixdrf, idssb(j)+nres, iret)
1753             call xdrfint(ixdrf, jdssb(j)+nres, iret)
1754            else
1755             call xdrfint(ixdrf, ihpb(j), iret)
1756             call xdrfint(ixdrf, jhpb(j), iret)
1757            endif
1758           enddo
1759           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1760           call xdrfint(ixdrf, iset_restart1(il), iret)
1761           do i=1,nfrag
1762            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1763           enddo
1764           do i=1,npair
1765            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1766           enddo
1767           do i=1,nfrag_back
1768            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1769            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1770            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1771           enddo
1772           prec=10000.0
1773           do i=1,nres
1774            do j=1,3
1775             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1776            enddo
1777           enddo
1778           do i=nnt,nct
1779            do j=1,3
1780             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1781            enddo
1782           enddo
1783           itmp=nres+nct-nnt+1
1784           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1785          enddo
1786 #endif
1787        endif
1788       enddo
1789 #ifdef AIX
1790       if(me.eq.king) call xdrfclose_(ixdrf, iret)
1791 #else
1792       if(me.eq.king) call xdrfclose(ixdrf, iret)
1793 #endif
1794       do i=1,ntwx_cache-ii_write
1795
1796             totT_cache(i)=totT_cache(ii_write+i)
1797             EK_cache(i)=EK_cache(ii_write+i)
1798             potE_cache(i)=potE_cache(ii_write+i)
1799             t_bath_cache(i)=t_bath_cache(ii_write+i)
1800             Uconst_cache(i)=Uconst_cache(ii_write+i)
1801             iset_cache(i)=iset_cache(ii_write+i)
1802
1803             do ii=1,nfrag
1804              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1805             enddo
1806             do ii=1,npair
1807              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1808             enddo
1809             do ii=1,nfrag_back
1810               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1811               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1812               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1813             enddo
1814
1815             do ii=1,nres*2
1816              do j=1,3
1817               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1818              enddo
1819             enddo
1820       enddo
1821       ntwx_cache=ntwx_cache-ii_write
1822       return
1823       end
1824
1825
1826       subroutine read1restart(i_index)
1827       implicit real*8 (a-h,o-z)
1828       include 'DIMENSIONS'
1829       include 'mpif.h'
1830       include 'COMMON.MD'
1831       include 'COMMON.IOUNITS'
1832       include 'COMMON.REMD'
1833       include 'COMMON.SETUP'
1834       include 'COMMON.CHAIN'
1835       include 'COMMON.SBRIDGE'
1836       include 'COMMON.INTERACT'
1837       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1838      &                 t5_restart1(5)
1839       integer*2 i_index
1840      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1841       common /przechowalnia/ d_restart1
1842       write (*,*) "Processor",me," called read1restart"
1843
1844          if(me.eq.king)then
1845               open(irest2,file=mremd_rst_name,status='unknown')
1846               read(irest2,*,err=334) i
1847               write(iout,*) "Reading old rst in ASCI format"
1848               close(irest2)
1849                call read1restart_old
1850                return
1851  334          continue
1852 #ifdef AIX
1853               call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1854
1855               do i=0,nodes-1
1856                call xdrfint_(ixdrf, i2rep(i), iret)
1857               enddo
1858               do i=1,remd_m(1)
1859                call xdrfint_(ixdrf, ifirst(i), iret)
1860               enddo
1861              do il=1,nodes
1862               call xdrfint_(ixdrf, nupa(0,il), iret)
1863               do i=1,nupa(0,il)
1864                call xdrfint_(ixdrf, nupa(i,il), iret)
1865               enddo
1866
1867               call xdrfint_(ixdrf, ndowna(0,il), iret)
1868               do i=1,ndowna(0,il)
1869                call xdrfint_(ixdrf, ndowna(i,il), iret)
1870               enddo
1871              enddo
1872              do il=1,nodes
1873                do j=1,4
1874                 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1875                enddo
1876              enddo
1877 #else
1878               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1879
1880               do i=0,nodes-1
1881                call xdrfint(ixdrf, i2rep(i), iret)
1882               enddo
1883               do i=1,remd_m(1)
1884                call xdrfint(ixdrf, ifirst(i), iret)
1885               enddo
1886              do il=1,nodes
1887               call xdrfint(ixdrf, nupa(0,il), iret)
1888               do i=1,nupa(0,il)
1889                call xdrfint(ixdrf, nupa(i,il), iret)
1890               enddo
1891
1892               call xdrfint(ixdrf, ndowna(0,il), iret)
1893               do i=1,ndowna(0,il)
1894                call xdrfint(ixdrf, ndowna(i,il), iret)
1895               enddo
1896              enddo
1897              do il=1,nodes
1898                do j=1,4
1899                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1900                enddo
1901              enddo
1902 #endif
1903          endif
1904          call mpi_scatter(t_restart1,5,mpi_real,
1905      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1906          totT=t5_restart1(1)
1907          EK=t5_restart1(2)
1908          potE=t5_restart1(3)
1909          t_bath=t5_restart1(4)
1910
1911          if(me.eq.king)then
1912               do il=0,nodes-1
1913                do i=1,2*nres
1914 c                read(irest2,'(3e15.5)')
1915 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1916             do j=1,3
1917 #ifdef AIX
1918              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1919 #else
1920              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1921 #endif
1922             enddo
1923                enddo
1924               enddo
1925          endif
1926          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1927      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1928
1929          do i=1,2*nres
1930            do j=1,3
1931             d_t(j,i)=r_d(j,i)
1932            enddo
1933          enddo
1934          if(me.eq.king)then
1935               do il=0,nodes-1
1936                do i=1,2*nres
1937 c                read(irest2,'(3e15.5)')
1938 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1939             do j=1,3
1940 #ifdef AIX
1941              call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1942 #else
1943              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1944 #endif
1945             enddo
1946                enddo
1947               enddo
1948          endif
1949          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1950      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1951          do i=1,2*nres
1952            do j=1,3
1953             dc(j,i)=r_d(j,i)
1954            enddo
1955          enddo
1956
1957
1958            if(usampl) then
1959 #ifdef AIX
1960              if(me.eq.king)then
1961               call xdrfint_(ixdrf, nset, iret)
1962               do i=1,nset
1963                 call xdrfint_(ixdrf,mset(i), iret)
1964               enddo
1965               do i=0,nodes-1
1966                 call xdrfint_(ixdrf,i2set(i), iret)
1967               enddo
1968               do il=1,nset
1969                do il1=1,mset(il)
1970                 do i=1,nrep
1971                  do j=1,remd_m(i)
1972                    call xdrfint_(ixdrf,itmp, iret)
1973                    i_index(i,j,il,il1)=itmp
1974                  enddo
1975                 enddo
1976                enddo
1977               enddo
1978              endif
1979 #else
1980              if(me.eq.king)then
1981               call xdrfint(ixdrf, nset, iret)
1982               do i=1,nset
1983                 call xdrfint(ixdrf,mset(i), iret)
1984               enddo
1985               do i=0,nodes-1
1986                 call xdrfint(ixdrf,i2set(i), iret)
1987               enddo
1988               do il=1,nset
1989                do il1=1,mset(il)
1990                 do i=1,nrep
1991                  do j=1,remd_m(i)
1992                    call xdrfint(ixdrf,itmp, iret)
1993                    i_index(i,j,il,il1)=itmp
1994                  enddo
1995                 enddo
1996                enddo
1997               enddo
1998              endif
1999 #endif
2000               call mpi_scatter(i2set,1,mpi_integer,
2001      &           iset,1,mpi_integer,king,
2002      &           CG_COMM,ierr)
2003
2004            endif
2005
2006
2007         if(me.eq.king) close(irest2)
2008         return
2009         end
2010
2011       subroutine read1restart_old
2012       implicit real*8 (a-h,o-z)
2013       include 'DIMENSIONS'
2014       include 'mpif.h'
2015       include 'COMMON.MD'
2016       include 'COMMON.IOUNITS'
2017       include 'COMMON.REMD'
2018       include 'COMMON.SETUP'
2019       include 'COMMON.CHAIN'
2020       include 'COMMON.SBRIDGE'
2021       include 'COMMON.INTERACT'
2022       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
2023      &                 t5_restart1(5)
2024       common /przechowalnia/ d_restart1
2025          if(me.eq.king)then
2026              open(irest2,file=mremd_rst_name,status='unknown')
2027              read (irest2,*) (i2rep(i),i=0,nodes-1)
2028              read (irest2,*) (ifirst(i),i=1,remd_m(1))
2029              do il=1,nodes
2030               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2031               read (irest2,*) ndowna(0,il),
2032      &                    (ndowna(i,il),i=1,ndowna(0,il))
2033              enddo
2034              do il=1,nodes
2035                read(irest2,*) (t_restart1(j,il),j=1,4)
2036              enddo
2037          endif
2038          call mpi_scatter(t_restart1,5,mpi_real,
2039      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2040          totT=t5_restart1(1)
2041          EK=t5_restart1(2)
2042          potE=t5_restart1(3)
2043          t_bath=t5_restart1(4)
2044
2045          if(me.eq.king)then
2046               do il=0,nodes-1
2047                do i=1,2*nres
2048                 read(irest2,'(3e15.5)')
2049      &                (d_restart1(j,i+2*nres*il),j=1,3)
2050                enddo
2051               enddo
2052          endif
2053          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2054      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2055
2056          do i=1,2*nres
2057            do j=1,3
2058             d_t(j,i)=r_d(j,i)
2059            enddo
2060          enddo
2061          if(me.eq.king)then
2062               do il=0,nodes-1
2063                do i=1,2*nres
2064                 read(irest2,'(3e15.5)')
2065      &                (d_restart1(j,i+2*nres*il),j=1,3)
2066                enddo
2067               enddo
2068          endif
2069          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2070      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2071          do i=1,2*nres
2072            do j=1,3
2073             dc(j,i)=r_d(j,i)
2074            enddo
2075          enddo
2076         if(me.eq.king) close(irest2)
2077         return
2078         end
2079 c-------------------------------------------------------------------
2080         subroutine set_hweights(iiset)
2081         implicit real*8 (a-h,o-z)
2082         integer i
2083         include 'DIMENSIONS'
2084         include 'COMMON.FFIELD'
2085         include 'COMMON.REMD'
2086
2087          do i=1,n_ene
2088           weights(i)=hweights(iiset,i)
2089          enddo
2090
2091          wsc    =weights(1)
2092          wscp   =weights(2)
2093          welec  =weights(3)
2094          wcorr  =weights(4)
2095          wcorr5 =weights(5)
2096          wcorr6 =weights(6)
2097          wel_loc=weights(7)
2098          wturn3 =weights(8)
2099          wturn4 =weights(9)
2100          wturn6 =weights(10)
2101          wang   =weights(11)
2102          wscloc =weights(12)
2103          wtor   =weights(13)
2104          wtor_d =weights(14)
2105          wstrain=weights(15)
2106          wvdwpp =weights(16)
2107          wbond  =weights(17)
2108          scal14 =weights(18)
2109          wsccor =weights(21)
2110
2111         return
2112         end
2113 #endif