source/unres/MREMD.f90

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