source/unres/src_MD-M/MREMD_nosy1traj.F

   1       subroutine MREMD
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'mpif.h'
   5       include 'COMMON.CONTROL'
   6       include 'COMMON.VAR'
   7       include 'COMMON.MD'
   8 #ifndef LANG0
   9       include 'COMMON.LANGEVIN'
  10 #else
  11       include 'COMMON.LANGEVIN.lang0'
  12 #endif
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.DERIV'
  15       include 'COMMON.GEO'
  16       include 'COMMON.LOCAL'
  17       include 'COMMON.INTERACT'
  18       include 'COMMON.IOUNITS'
  19       include 'COMMON.NAMES'
  20       include 'COMMON.TIME1'
  21       include 'COMMON.REMD'
  22       include 'COMMON.SETUP'
  23       include 'COMMON.MUCA'
  24       include 'COMMON.HAIRPIN'
  25       double precision cm(3),L(3),vcm(3)
  26       double precision energia(0:n_ene)
  27       double precision remd_t_bath(maxprocs)
  28       double precision remd_ene(0:n_ene+1,maxprocs)
  29       integer iremd_acc(maxprocs),iremd_tot(maxprocs)
  30       integer ilen,rstcount
  31       external ilen
  32       character*50 tytul
  33       common /gucio/ cm
  34       integer itime
  35       integer nup(0:maxprocs),ndown(0:maxprocs)
  36       integer rep2i(0:maxprocs)
  37       integer itime_all(maxprocs)
  38       integer status(MPI_STATUS_SIZE)
  39       logical synflag,end_of_run,file_exist
  40
  41       time00=MPI_WTIME()
  42       if(me.eq.king.or..not.out1file)
  43      & write  (iout,*) 'MREMD',nodes,'time before',time00-walltime
  44
  45       synflag=.false.
  46       mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
  47
  48 cd      print *,'MREMD',nodes
  49
  50 cd      print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
  51
  52 cde      write (iout,*) "Start MREMD: me",me," t_bath",t_bath
  53       k=0
  54       rep2i(k)=-1
  55       do i=1,nrep
  56         iremd_acc(i)=0
  57         iremd_tot(i)=0
  58         do j=1,remd_m(i)
  59           i2rep(k)=i
  60           rep2i(i)=k
  61           k=k+1
  62         enddo
  63       enddo
  64
  65 c      print *,'i2rep',me,i2rep(me)
  66 c      print *,'rep2i',(rep2i(i),i=0,nrep)
  67
  68        if (i2rep(me).eq.nrep) then
  69         nup(0)=0
  70        else
  71         nup(0)=remd_m(i2rep(me)+1)
  72         k=rep2i(i2rep(me))+1
  73         do i=1,nup(0)
  74          nup(i)=k
  75          k=k+1
  76         enddo
  77        endif
  78
  79 cd       print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
  80
  81        if (i2rep(me).eq.1) then
  82         ndown(0)=0
  83        else
  84         ndown(0)=remd_m(i2rep(me)-1)
  85         k=rep2i(i2rep(me)-2)+1
  86         do i=1,ndown(0)
  87          ndown(i)=k
  88          k=k+1
  89         enddo
  90        endif
  91
  92 cd       print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
  93
  94
  95        if(rest.and.restart1file) then
  96            inquire(file=mremd_rst_name,exist=file_exist)
  97            if(file_exist) call read1restart
  98        endif
  99
 100        if(me.eq.king) then
 101         if (rest.and..not.restart1file)
 102      &          inquire(file=mremd_rst_name,exist=file_exist)
 103         IF (rest.and.file_exist.and..not.restart1file) THEN
 104              open(irest2,file=mremd_rst_name,status='unknown')
 105              read (irest2,*)
 106              read (irest2,*) (i2rep(i),i=0,nodes-1)
 107              read (irest2,*)
 108              read (irest2,*) (ifirst(i),i=1,remd_m(1))
 109              do il=1,nodes
 110               read (irest2,*)
 111               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 112               read (irest2,*)
 113               read (irest2,*) ndowna(0,il),
 114      &                    (ndowna(i,il),i=1,ndowna(0,il))
 115              enddo
 116              close(irest2)
 117
 118              write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
 119              write (iout,'(a6,1000i5)') "ifirst",
 120      &                    (ifirst(i),i=1,remd_m(1))
 121              do il=1,nodes
 122               write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
 123      &                    (nupa(i,il),i=1,nupa(0,il))
 124               write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
 125      &                    (ndowna(i,il),i=1,ndowna(0,il))
 126              enddo
 127         ELSE
 128          do il=1,remd_m(1)
 129           ifirst(il)=il
 130          enddo
 131
 132          do il=1,nodes
 133
 134           if (i2rep(il-1).eq.nrep) then
 135            nupa(0,il)=0
 136           else
 137            nupa(0,il)=remd_m(i2rep(il-1)+1)
 138            k=rep2i(i2rep(il-1))+1
 139            do i=1,nupa(0,il)
 140             nupa(i,il)=k+1
 141             k=k+1
 142            enddo
 143           endif
 144
 145           if (i2rep(il-1).eq.1) then
 146            ndowna(0,il)=0
 147           else
 148            ndowna(0,il)=remd_m(i2rep(il-1)-1)
 149            k=rep2i(i2rep(il-1)-2)+1
 150            do i=1,ndowna(0,il)
 151             ndowna(i,il)=k+1
 152             k=k+1
 153            enddo
 154           endif
 155
 156          enddo
 157
 158 c        write (iout,'(a6,100i4)') "ifirst",
 159 c     &                    (ifirst(i),i=1,remd_m(1))
 160 c        do il=1,nodes
 161 c         write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
 162 c     &                    (nupa(i,il),i=1,nupa(0,il))
 163 c         write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
 164 c     &                    (ndowna(i,il),i=1,ndowna(0,il))
 165 c        enddo
 166
 167         ENDIF
 168        endif
 169 c
 170 c      t_bath=retmin+(retmax-retmin)*me/(nodes-1)
 171        if(.not.(rest.and.file_exist.and.restart1file)) then
 172          if (me .eq. king) then
 173             t_bath=retmin
 174          else
 175             t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
 176          endif
 177 cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
 178          if (remd_tlist) t_bath=remd_t(i2rep(me))
 179        endif
 180 c
 181         stdfp=dsqrt(2*Rb*t_bath/d_time)
 182         do i=1,ntyp
 183           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 184         enddo
 185
 186 c      print *,'irep',me,t_bath
 187        if (.not.rest) then
 188         if (me.eq.king .or. .not. out1file)
 189      &   write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
 190         call rescale_weights(t_bath)
 191        endif
 192
 193
 194 c------copy MD--------------
 195 c  The driver for molecular dynamics subroutines
 196 c------------------------------------------------
 197       t_MDsetup=0.0d0
 198       t_langsetup=0.0d0
 199       t_MD=0.0d0
 200       t_enegrad=0.0d0
 201       t_sdsetup=0.0d0
 202       if(me.eq.king.or..not.out1file)
 203      & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
 204       tt0 = tcpu()
 205 c Determine the inverse of the inertia matrix.
 206       call setup_MD_matrices
 207 c Initialize MD
 208       call init_MD
 209       if (rest) then
 210        if (me.eq.king .or. .not. out1file)
 211      &  write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
 212        stdfp=dsqrt(2*Rb*t_bath/d_time)
 213        do i=1,ntyp
 214           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 215        enddo
 216        call rescale_weights(t_bath)
 217       endif
 218
 219       t_MDsetup = tcpu()-tt0
 220       rstcount=0
 221 c   Entering the MD loop
 222       tt0 = tcpu()
 223       if (lang.eq.2 .or. lang.eq.3) then
 224         call setup_fricmat
 225         if (lang.eq.2) then
 226           call sd_verlet_p_setup
 227         else
 228           call sd_verlet_ciccotti_setup
 229         endif
 230 #ifndef   LANG0
 231         do i=1,dimen
 232           do j=1,dimen
 233             pfric0_mat(i,j,0)=pfric_mat(i,j)
 234             afric0_mat(i,j,0)=afric_mat(i,j)
 235             vfric0_mat(i,j,0)=vfric_mat(i,j)
 236             prand0_mat(i,j,0)=prand_mat(i,j)
 237             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 238             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 239           enddo
 240         enddo
 241 #endif
 242         flag_stoch(0)=.true.
 243         do i=1,maxflag_stoch
 244           flag_stoch(i)=.false.
 245         enddo
 246       else if (lang.eq.1 .or. lang.eq.4) then
 247         call setup_fricmat
 248       endif
 249       time00=MPI_WTIME()
 250       if (me.eq.king .or. .not. out1file)
 251      & write(iout,*) 'Setup time',time00-walltime
 252       call flush(iout)
 253       t_langsetup=tcpu()-tt0
 254       tt0=tcpu()
 255       itime=0
 256       end_of_run=.false.
 257       do while(.not.end_of_run)
 258         itime=itime+1
 259         if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
 260         if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
 261         rstcount=rstcount+1
 262         if (lang.gt.0 .and. surfarea .and.
 263      &      mod(itime,reset_fricmat).eq.0) then
 264           if (lang.eq.2 .or. lang.eq.3) then
 265             call setup_fricmat
 266             if (lang.eq.2) then
 267               call sd_verlet_p_setup
 268             else
 269               call sd_verlet_ciccotti_setup
 270             endif
 271 #ifndef   LANG0
 272             do i=1,dimen
 273               do j=1,dimen
 274                 pfric0_mat(i,j,0)=pfric_mat(i,j)
 275                 afric0_mat(i,j,0)=afric_mat(i,j)
 276                 vfric0_mat(i,j,0)=vfric_mat(i,j)
 277                 prand0_mat(i,j,0)=prand_mat(i,j)
 278                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 279                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 280               enddo
 281             enddo
 282 #endif
 283             flag_stoch(0)=.true.
 284             do i=1,maxflag_stoch
 285               flag_stoch(i)=.false.
 286             enddo
 287           else if (lang.eq.1 .or. lang.eq.4) then
 288             call setup_fricmat
 289           endif
 290           write (iout,'(a,i10)')
 291      &      "Friction matrix reset based on surface area, itime",itime
 292         endif
 293         if (reset_vel .and. tbf .and. lang.eq.0
 294      &      .and. mod(itime,count_reset_vel).eq.0) then
 295           call random_vel
 296           if (me.eq.king .or. .not. out1file)
 297      &     write(iout,'(a,f20.2)')
 298      &     "Velocities reset to random values, time",totT
 299           do i=0,2*nres
 300             do j=1,3
 301               d_t_old(j,i)=d_t(j,i)
 302             enddo
 303           enddo
 304         endif
 305         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 306           call inertia_tensor
 307           call vcm_vel(vcm)
 308           do j=1,3
 309              d_t(j,0)=d_t(j,0)-vcm(j)
 310           enddo
 311           call kinetic(EK)
 312           kinetic_T=2.0d0/(dimen*Rb)*EK
 313           scalfac=dsqrt(T_bath/kinetic_T)
 314 cd          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
 315           do i=0,2*nres
 316             do j=1,3
 317               d_t_old(j,i)=scalfac*d_t(j,i)
 318             enddo
 319           enddo
 320         endif
 321         if (lang.ne.4) then
 322           if (RESPA) then
 323 c Time-reversible RESPA algorithm
 324 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 325             call RESPA_step(itime)
 326           else
 327 c Variable time step algorithm.
 328             call velverlet_step(itime)
 329           endif
 330         else
 331           call brown_step(itime)
 332         endif
 333         if(ntwe.ne.0) then
 334           if (mod(itime,ntwe).eq.0) call statout(itime)
 335         endif
 336         if (mod(itime,ntwx).eq.0.and..not.traj1file) then
 337           write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
 338           if(mdpdb) then
 339              call hairpin(.true.,nharp,iharp)
 340              call secondary2(.true.)
 341              call pdbout(potE,tytul,ipdb)
 342           else
 343              call cartout(totT)
 344           endif
 345         endif
 346         if ((rstcount.eq.1000.or.itime.eq.n_timestep)
 347      &                         .and..not.restart1file) then
 348
 349            if(me.eq.king) then
 350              open(irest1,file=mremd_rst_name,status='unknown')
 351              write (irest1,*) "i2rep"
 352              write (irest1,*) (i2rep(i),i=0,nodes-1)
 353              write (irest1,*) "ifirst"
 354              write (irest1,*) (ifirst(i),i=1,remd_m(1))
 355              do il=1,nodes
 356               write (irest1,*) "nupa",il
 357               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 358               write (irest1,*) "ndowna",il
 359               write (irest1,*) ndowna(0,il),
 360      &                   (ndowna(i,il),i=1,ndowna(0,il))
 361              enddo
 362              close(irest1)
 363            endif
 364            open(irest2,file=rest2name,status='unknown')
 365            write(irest2,*) totT,EK,potE,totE,t_bath
 366            do i=1,2*nres
 367             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 368            enddo
 369            do i=1,2*nres
 370             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
 371            enddo
 372           close(irest2)
 373           rstcount=0
 374         endif
 375
 376 c REMD - exchange
 377 c forced synchronization
 378         if (mod(itime,100).eq.0 .and. me.ne.king
 379      &                                .and. .not. mremdsync) then
 380             synflag=.false.
 381             call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
 382             if (synflag) then
 383                call mpi_recv(itime_master, 1, MPI_INTEGER,
 384      &                             0,101,CG_COMM, status, ierr)
 385                if (.not. out1file) then
 386                 write(iout,*) 'REMD synchro at',itime_master,itime
 387                else
 388                 call mpi_gather(itime,1,mpi_integer,
 389      &             itime_all,1,mpi_integer,king,
 390      &             CG_COMM,ierr)
 391                endif
 392                if(itime_master.ge.n_timestep) end_of_run=.true.
 393                call flush(iout)
 394             endif
 395         endif
 396
 397 c REMD - exchange
 398         if ((mod(itime,nstex).eq.0.and.me.eq.king
 399      &                  .or.end_of_run.and.me.eq.king )
 400      &       .and. .not. mremdsync ) then
 401            synflag=.true.
 402            time00=MPI_WTIME()
 403            do i=1,nodes-1
 404             call mpi_isend(itime,1,MPI_INTEGER,i,101,
 405      &                                   CG_COMM, ireq, ierr)
 406 cd            write(iout,*) 'REMD synchro with',i
 407 cd            call flush(iout)
 408            enddo
 409            time02=MPI_WTIME()
 410            write(iout,*) 'REMD synchro at',itime,'time=',time02-time00
 411            if (out1file) then
 412             call mpi_gather(itime,1,mpi_integer,
 413      &             itime_all,1,mpi_integer,king,
 414      &             CG_COMM,ierr)
 415             write(iout,'(a18,8000i8)') 'REMD synchro itime',
 416      &                    (itime_all(i),i=1,nodes)
 417            endif
 418            call flush(iout)
 419         endif
 420         if(mremdsync .and. mod(itime,nstex).eq.0) then
 421            synflag=.true.
 422            if (me.eq.king .or. .not. out1file)
 423      &      write(iout,*) 'REMD synchro at',itime
 424            call flush(iout)
 425         endif
 426         if(synflag.and..not.end_of_run) then
 427            time00=MPI_WTIME()
 428            synflag=.false.
 429
 430 cd           write(iout,*) 'REMD before',me,t_bath
 431
 432 c           call mpi_gather(t_bath,1,mpi_double_precision,
 433 c     &             remd_t_bath,1,mpi_double_precision,king,
 434 c     &             CG_COMM,ierr)
 435            potEcomp(n_ene+1)=t_bath
 436            call mpi_gather(potEcomp(0),n_ene+2,mpi_double_precision,
 437      &             remd_ene(0,1),n_ene+2,mpi_double_precision,king,
 438      &             CG_COMM,ierr)
 439            if(lmuca) then
 440             call mpi_gather(elow,1,mpi_double_precision,
 441      &             elowi,1,mpi_double_precision,king,
 442      &             CG_COMM,ierr)
 443             call mpi_gather(ehigh,1,mpi_double_precision,
 444      &             ehighi,1,mpi_double_precision,king,
 445      &             CG_COMM,ierr)
 446            endif
 447
 448           if (restart1file) call write1rst
 449           if (traj1file) call write1traj
 450
 451           if (me.eq.king) then
 452             do i=1,nodes
 453                remd_t_bath(i)=remd_ene(n_ene+1,i)
 454             enddo
 455             if(lmuca) then
 456 co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
 457              do i=1,nodes
 458                write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
 459      &            elowi(i),ehighi(i)
 460              enddo
 461             else
 462 cd              write(iout,*) 'REMD exchange temp,ene'
 463 c             do i=1,nodes
 464 co              write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
 465 c              write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
 466 c             enddo
 467             endif
 468 c-------------------------------------
 469            do irr=1,remd_m(1)
 470            i=ifirst(iran_num(1,remd_m(1)))
 471            do ii=1,nodes-1
 472
 473             if(i.gt.0.and.nupa(0,i).gt.0) then
 474               iex=nupa(iran_num(1,int(nupa(0,i))),i)
 475              if (lmuca) then
 476                call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
 477              else
 478 c Swap temperatures between conformations i and iex with recalculating the free energies
 479 c following temperature changes.
 480               ene_iex_iex=remd_ene(0,iex)
 481               ene_i_i=remd_ene(0,i)
 482 co              write (iout,*) "rescaling weights with temperature",
 483 co     &          remd_t_bath(i)
 484               call rescale_weights(remd_t_bath(i))
 485               call sum_energy(remd_ene(0,iex),.false.)
 486               ene_iex_i=remd_ene(0,iex)
 487 cd              write (iout,*) "ene_iex_i",remd_ene(0,iex)
 488               call sum_energy(remd_ene(0,i),.false.)
 489 cd              write (iout,*) "ene_i_i",remd_ene(0,i)
 490 c              write (iout,*) "rescaling weights with temperature",
 491 c     &          remd_t_bath(iex)
 492               if (real(ene_i_i).ne.real(remd_ene(0,i))) then
 493                 write (iout,*) "ERROR: inconsistent energies:",i,
 494      &            ene_i_i,remd_ene(0,i)
 495               endif
 496               call rescale_weights(remd_t_bath(iex))
 497               call sum_energy(remd_ene(0,i),.false.)
 498 c              write (iout,*) "ene_i_iex",remd_ene(0,i)
 499               ene_i_iex=remd_ene(0,i)
 500               call sum_energy(remd_ene(0,iex),.false.)
 501               if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
 502                 write (iout,*) "ERROR: inconsistent energies:",iex,
 503      &            ene_iex_iex,remd_ene(0,iex)
 504               endif
 505 c              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
 506 c              write (iout,*) "i",i," iex",iex
 507 co              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
 508 co     &           " ene_i_iex",ene_i_iex,
 509 co     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
 510               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
 511      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
 512               delta=-delta
 513 c              write(iout,*) 'delta',delta
 514 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
 515 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
 516 c     &              (remd_t_bath(i)*remd_t_bath(iex))
 517              endif
 518               if (delta .gt. 50.0d0) then
 519                 delta=0.0d0
 520               else
 521 #ifdef OSF
 522                 if(isnan(delta))then
 523                   delta=0.0d0
 524                 else if (delta.lt.-50.0d0) then
 525                   delta=dexp(50.0d0)
 526                 else
 527                 delta=dexp(-delta)
 528               endif
 529 #else
 530                 delta=dexp(-delta)
 531 #endif
 532               endif
 533               iremd_tot(i2rep(i-1))=iremd_tot(i2rep(i-1))+1
 534               xxx=ran_number(0.0d0,1.0d0)
 535 co              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
 536               if (delta .gt. xxx) then
 537                 tmp=remd_t_bath(i)
 538                 remd_t_bath(i)=remd_t_bath(iex)
 539                 remd_t_bath(iex)=tmp
 540                 remd_ene(0,i)=ene_i_iex
 541                 remd_ene(0,iex)=ene_iex_i
 542                 if(lmuca) then
 543                   tmp=elowi(i)
 544                   elowi(i)=elowi(iex)
 545                   elowi(iex)=tmp
 546                   tmp=ehighi(i)
 547                   ehighi(i)=ehighi(iex)
 548                   ehighi(iex)=tmp
 549                 endif
 550
 551
 552                 do k=0,nodes
 553                   itmp=nupa(k,i)
 554                   nupa(k,i)=nupa(k,iex)
 555                   nupa(k,iex)=itmp
 556                   itmp=ndowna(k,i)
 557                   ndowna(k,i)=ndowna(k,iex)
 558                   ndowna(k,iex)=itmp
 559                 enddo
 560                 do il=1,nodes
 561                  if (ifirst(il).eq.i) ifirst(il)=iex
 562                  do k=1,nupa(0,il)
 563                   if (nupa(k,il).eq.i) then
 564                      nupa(k,il)=iex
 565                   elseif (nupa(k,il).eq.iex) then
 566                      nupa(k,il)=i
 567                   endif
 568                  enddo
 569                  do k=1,ndowna(0,il)
 570                   if (ndowna(k,il).eq.i) then
 571                      ndowna(k,il)=iex
 572                   elseif (ndowna(k,il).eq.iex) then
 573                      ndowna(k,il)=i
 574                   endif
 575                  enddo
 576                 enddo
 577
 578                 iremd_acc(i2rep(i-1))=iremd_acc(i2rep(i-1))+1
 579                 itmp=i2rep(i-1)
 580                 i2rep(i-1)=i2rep(iex-1)
 581                 i2rep(iex-1)=itmp
 582
 583 cd                write(iout,*) 'exchange',i,iex
 584 cd              write (iout,'(a8,100i4)') "@ ifirst",
 585 cd     &                    (ifirst(k),k=1,remd_m(1))
 586 cd              do il=1,nodes
 587 cd               write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
 588 cd     &                    (nupa(k,il),k=1,nupa(0,il))
 589 cd               write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
 590 cd     &                    (ndowna(k,il),k=1,ndowna(0,il))
 591 cd              enddo
 592
 593               else
 594                remd_ene(0,iex)=ene_iex_iex
 595                remd_ene(0,i)=ene_i_i
 596                i=iex
 597               endif
 598             endif
 599            enddo
 600            enddo
 601 c-------------------------------------
 602              do i=1,nrep
 603               if(iremd_tot(i).ne.0)
 604      &          write(iout,'(a3,i4,2f12.5)') 'ACC',i,remd_t(i)
 605      &           ,iremd_acc(i)/(1.0*iremd_tot(i))
 606              enddo
 607
 608 cd              write (iout,'(a6,100i4)') "ifirst",
 609 cd     &                    (ifirst(i),i=1,remd_m(1))
 610 cd              do il=1,nodes
 611 cd               write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
 612 cd     &                    (nupa(i,il),i=1,nupa(0,il))
 613 cd               write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
 614 cd     &                    (ndowna(i,il),i=1,ndowna(0,il))
 615 cd              enddo
 616           endif
 617
 618          call mpi_scatter(remd_t_bath,1,mpi_double_precision,
 619      &           t_bath,1,mpi_double_precision,king,
 620      &           CG_COMM,ierr)
 621          if(lmuca) then
 622            call mpi_scatter(elowi,1,mpi_double_precision,
 623      &           elow,1,mpi_double_precision,king,
 624      &           CG_COMM,ierr)
 625            call mpi_scatter(ehighi,1,mpi_double_precision,
 626      &           ehigh,1,mpi_double_precision,king,
 627      &           CG_COMM,ierr)
 628          endif
 629          call rescale_weights(t_bath)
 630 co         write (iout,*) "Processor",me,
 631 co     &    " rescaling weights with temperature",t_bath
 632
 633          stdfp=dsqrt(2*Rb*t_bath/d_time)
 634          do i=1,ntyp
 635            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 636          enddo
 637
 638 cde         write(iout,*) 'REMD after',me,t_bath
 639            time02=MPI_WTIME()
 640            if (me.eq.king .or. .not. out1file) then
 641             write(iout,*) 'REMD exchange time=',time02-time00
 642             call flush(iout)
 643            endif
 644         endif
 645       enddo
 646
 647       if (restart1file) then
 648           if (me.eq.king .or. .not. out1file)
 649      &      write(iout,*) 'writing restart at the end of run'
 650            call write1rst
 651       endif
 652
 653       if (traj1file) call write1traj
 654
 655
 656       t_MD=tcpu()-tt0
 657       if (me.eq.king .or. .not. out1file) then
 658        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
 659      &  '  Timing  ',
 660      & 'MD calculations setup:',t_MDsetup,
 661      & 'Energy & gradient evaluation:',t_enegrad,
 662      & 'Stochastic MD setup:',t_langsetup,
 663      & 'Stochastic MD step setup:',t_sdsetup,
 664      & 'MD steps:',t_MD
 665        write (iout,'(/28(1h=),a25,27(1h=))')
 666      & '  End of MD calculation  '
 667       endif
 668       return
 669       end
 670
 671       subroutine write1rst
 672       implicit real*8 (a-h,o-z)
 673       include 'DIMENSIONS'
 674       include 'mpif.h'
 675       include 'COMMON.MD'
 676       include 'COMMON.IOUNITS'
 677       include 'COMMON.REMD'
 678       include 'COMMON.SETUP'
 679       include 'COMMON.CHAIN'
 680       include 'COMMON.SBRIDGE'
 681       include 'COMMON.INTERACT'
 682
 683       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
 684      &     d_restart2(3,2*maxres*maxprocs)
 685       real t5_restart1(5)
 686       integer iret,itmp
 687
 688
 689        t5_restart1(1)=totT
 690        t5_restart1(2)=EK
 691        t5_restart1(3)=potE
 692        t5_restart1(4)=t_bath
 693        t5_restart1(5)=Uconst
 694
 695        call mpi_gather(t5_restart1,5,mpi_real,
 696      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
 697
 698
 699        do i=1,2*nres
 700          do j=1,3
 701            r_d(j,i)=d_t(j,i)
 702          enddo
 703        enddo
 704        call mpi_gather(r_d,3*2*nres,mpi_real,
 705      &           d_restart1,3*2*nres,mpi_real,king,
 706      &           CG_COMM,ierr)
 707
 708
 709        do i=1,2*nres
 710          do j=1,3
 711            r_d(j,i)=dc(j,i)
 712          enddo
 713        enddo
 714        call mpi_gather(r_d,3*2*nres,mpi_real,
 715      &           d_restart2,3*2*nres,mpi_real,king,
 716      &           CG_COMM,ierr)
 717
 718        if(me.eq.king) then
 719          open(irest1,file=mremd_rst_name,status='unknown')
 720          write (irest1,*) (i2rep(i),i=0,nodes-1)
 721          write (irest1,*) (ifirst(i),i=1,remd_m(1))
 722          do il=1,nodes
 723               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 724               write (irest1,*) ndowna(0,il),
 725      &                   (ndowna(i,il),i=1,ndowna(0,il))
 726          enddo
 727
 728          do il=1,nodes
 729            write(irest1,*) (t_restart1(j,il),j=1,4)
 730          enddo
 731
 732          do il=0,nodes-1
 733            do i=1,2*nres
 734              write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
 735            enddo
 736          enddo
 737          do il=0,nodes-1
 738            do i=1,2*nres
 739              write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
 740            enddo
 741          enddo
 742          close(irest1)
 743        endif
 744
 745
 746       return
 747       end
 748
 749       subroutine write1traj
 750       implicit real*8 (a-h,o-z)
 751       include 'DIMENSIONS'
 752       include 'mpif.h'
 753       include 'COMMON.MD'
 754       include 'COMMON.IOUNITS'
 755       include 'COMMON.REMD'
 756       include 'COMMON.SETUP'
 757       include 'COMMON.CHAIN'
 758       include 'COMMON.SBRIDGE'
 759       include 'COMMON.INTERACT'
 760
 761       real t5_restart1(5)
 762       integer iret,itmp
 763       real xcoord(3,maxres2+2),prec
 764       real r_qfrag(50),r_qpair(100)
 765       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
 766       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
 767
 768       if(.not.restart1file) then
 769        t5_restart1(1)=totT
 770        t5_restart1(2)=EK
 771        t5_restart1(3)=potE
 772        t5_restart1(4)=t_bath
 773        t5_restart1(5)=Uconst
 774        call mpi_gather(t5_restart1,5,mpi_real,
 775      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
 776       endif
 777
 778           do i=1,nfrag
 779            r_qfrag(i)=qfrag(i)
 780           enddo
 781           do i=1,npair
 782            r_qpair(i)=qpair(i)
 783           enddo
 784
 785         call mpi_gather(r_qfrag,nfrag,mpi_real,
 786      &           p_qfrag,nfrag,mpi_real,king,
 787      &           CG_COMM,ierr)
 788         call mpi_gather(qpair,nfrag,mpi_real,
 789      &           p_qpair,nfrag,mpi_real,king,
 790      &           CG_COMM,ierr)
 791
 792         do i=1,nres*2
 793          do j=1,3
 794           r_c(j,i)=c(j,i)
 795          enddo
 796         enddo
 797
 798         call mpi_gather(r_c,3*2*nres,mpi_real,
 799      &           p_c,3*2*nres,mpi_real,king,
 800      &           CG_COMM,ierr)
 801
 802        if(me.eq.king) then
 803          call xdrfopen(ixdrf,cartname, "a", iret)
 804          do il=1,nodes
 805           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
 806           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
 807           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
 808           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
 809           call xdrfint(ixdrf, nss, iret)
 810           do j=1,nss
 811            call xdrfint(ixdrf, ihpb(j), iret)
 812            call xdrfint(ixdrf, jhpb(j), iret)
 813           enddo
 814           call xdrfint(ixdrf, nfrag+npair, iret)
 815           do i=1,nfrag
 816            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
 817           enddo
 818           do i=1,npair
 819            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
 820           enddo
 821           prec=10000.0
 822           do i=1,nres
 823            do j=1,3
 824             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
 825            enddo
 826           enddo
 827           do i=nnt,nct
 828            do j=1,3
 829             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
 830            enddo
 831           enddo
 832           itmp=nres+nct-nnt+1
 833           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
 834          enddo
 835          call xdrfclose(ixdrf, iret)
 836        endif
 837
 838       return
 839       end
 840
 841
 842       subroutine read1restart
 843       implicit real*8 (a-h,o-z)
 844       include 'DIMENSIONS'
 845       include 'mpif.h'
 846       include 'COMMON.MD'
 847       include 'COMMON.IOUNITS'
 848       include 'COMMON.REMD'
 849       include 'COMMON.SETUP'
 850       include 'COMMON.CHAIN'
 851       include 'COMMON.SBRIDGE'
 852       include 'COMMON.INTERACT'
 853       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
 854      &                 t5_restart1(5)
 855
 856          if(me.eq.king)then
 857              open(irest2,file=mremd_rst_name,status='unknown')
 858              read (irest2,*) (i2rep(i),i=0,nodes-1)
 859              read (irest2,*) (ifirst(i),i=1,remd_m(1))
 860              do il=1,nodes
 861               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 862               read (irest2,*) ndowna(0,il),
 863      &                    (ndowna(i,il),i=1,ndowna(0,il))
 864              enddo
 865              do il=1,nodes
 866                read(irest2,*) (t_restart1(j,il),j=1,4)
 867              enddo
 868          endif
 869          call mpi_scatter(t_restart1,5,mpi_real,
 870      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
 871          totT=t5_restart1(1)
 872          EK=t5_restart1(2)
 873          potE=t5_restart1(3)
 874          t_bath=t5_restart1(4)
 875
 876          if(me.eq.king)then
 877               do il=0,nodes-1
 878                do i=1,2*nres
 879                 read(irest2,'(3e15.5)')
 880      &                (d_restart1(j,i+2*nres*il),j=1,3)
 881                enddo
 882               enddo
 883          endif
 884          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
 885      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
 886
 887          do i=1,2*nres
 888            do j=1,3
 889             d_t(j,i)=r_d(j,i)
 890            enddo
 891          enddo
 892          if(me.eq.king)then
 893               do il=0,nodes-1
 894                do i=1,2*nres
 895                 read(irest2,'(3e15.5)')
 896      &                (d_restart1(j,i+2*nres*il),j=1,3)
 897                enddo
 898               enddo
 899          endif
 900          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
 901      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
 902          do i=1,2*nres
 903            do j=1,3
 904             dc(j,i)=r_d(j,i)
 905            enddo
 906          enddo
 907         if(me.eq.king) close(irest2)
 908         return
 909         end
 910