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