source/unres/MREMD.f90

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