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