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