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