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