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