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