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