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