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