source/unres/src_MD-M/MREMD.F.drabinka

   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 cold      integer nup(0:maxprocs),ndown(0:maxprocs)
  35       integer rep2i(0:maxprocs),ireqi(maxprocs)
  36       integer itime_all(maxprocs)
  37       integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
  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 cold       if (i2rep(me).eq.nrep) then
  68 cold        nup(0)=0
  69 cold       else
  70 cold        nup(0)=remd_m(i2rep(me)+1)
  71 cold        k=rep2i(int(i2rep(me)))+1
  72 cold        do i=1,nup(0)
  73 cold         nup(i)=k
  74 cold         k=k+1
  75 cold        enddo
  76 cold       endif
  77
  78 cd       print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
  79
  80 cold       if (i2rep(me).eq.1) then
  81 cold        ndown(0)=0
  82 cold       else
  83 cold        ndown(0)=remd_m(i2rep(me)-1)
  84 cold        k=rep2i(i2rep(me)-2)+1
  85 cold        do i=1,ndown(0)
  86 cold         ndown(i)=k
  87 cold         k=k+1
  88 cold        enddo
  89 cold       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 IF (.not.(rest.and.file_exist)) THEN
 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(int(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(int(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 (mod(itime,ntwx).eq.0.and.traj1file) then
 344           if(ntwx_cache.lt.max_cache_traj) then
 345             ntwx_cache=ntwx_cache+1
 346
 347             totT_cache(ntwx_cache)=totT
 348             EK_cache(ntwx_cache)=EK
 349             potE_cache(ntwx_cache)=potE
 350             t_bath_cache(ntwx_cache)=t_bath
 351             Uconst_cache(ntwx_cache)=Uconst
 352
 353             do i=1,nfrag
 354              qfrag_cache(i,ntwx_cache)=qfrag(i)
 355             enddo
 356             do i=1,npair
 357              qpair_cache(i,ntwx_cache)=qpair(i)
 358             enddo
 359
 360             do i=1,nres*2
 361              do j=1,3
 362               c_cache(j,i,ntwx_cache)=c(j,i)
 363              enddo
 364             enddo
 365           else
 366            print *,itime,"processor ",me," over cache ",ntwx_cache
 367           endif
 368         endif
 369         if ((rstcount.eq.1000.or.itime.eq.n_timestep)
 370      &                         .and..not.restart1file) then
 371
 372            if(me.eq.king) then
 373              open(irest1,file=mremd_rst_name,status='unknown')
 374              write (irest1,*) "i2rep"
 375              write (irest1,*) (i2rep(i),i=0,nodes-1)
 376              write (irest1,*) "ifirst"
 377              write (irest1,*) (ifirst(i),i=1,remd_m(1))
 378              do il=1,nodes
 379               write (irest1,*) "nupa",il
 380               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 381               write (irest1,*) "ndowna",il
 382               write (irest1,*) ndowna(0,il),
 383      &                   (ndowna(i,il),i=1,ndowna(0,il))
 384              enddo
 385              close(irest1)
 386            endif
 387            open(irest2,file=rest2name,status='unknown')
 388            write(irest2,*) totT,EK,potE,totE,t_bath
 389            do i=1,2*nres
 390             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 391            enddo
 392            do i=1,2*nres
 393             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
 394            enddo
 395           close(irest2)
 396           rstcount=0
 397         endif
 398
 399 c REMD - exchange
 400 c forced synchronization
 401         if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
 402      &                                .and. .not. mremdsync) then
 403             synflag=.false.
 404             call mpi_iprobe(0,101,mpi_comm_world,synflag,status,ierr)
 405             if (synflag) then
 406                call mpi_recv(itime_master, 1, MPI_INTEGER,
 407      &                             0,101,mpi_comm_world, status, ierr)
 408                call mpi_barrier(mpi_comm_world, ierr)
 409                if (out1file.or.traj1file) then
 410                 call mpi_gather(itime,1,mpi_integer,
 411      &             itime_all,1,mpi_integer,king,
 412      &             mpi_comm_world,ierr)
 413                endif
 414                if (.not.out1file)
 415      &               write(iout,*) 'REMD synchro at',itime_master,itime
 416                if(itime_master.ge.n_timestep) end_of_run=.true.
 417 ctime               call flush(iout)
 418             endif
 419         endif
 420
 421 c REMD - exchange
 422         if ((mod(itime,nstex).eq.0.and.me.eq.king
 423      &                  .or.end_of_run.and.me.eq.king )
 424      &       .and. .not. mremdsync ) then
 425            synflag=.true.
 426            time00=MPI_WTIME()
 427            do i=1,nodes-1
 428               call mpi_isend(itime,1,MPI_INTEGER,i,101,
 429      &                                mpi_comm_world, ireqi(i), ierr)
 430 cd            write(iout,*) 'REMD synchro with',i
 431 cd            call flush(iout)
 432            enddo
 433            call mpi_waitall(nodes-1,ireqi,statusi,ierr)
 434            call mpi_barrier(mpi_comm_world, ierr)
 435            time01=MPI_WTIME()
 436            write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
 437            if (out1file.or.traj1file) then
 438             call mpi_gather(itime,1,mpi_integer,
 439      &             itime_all,1,mpi_integer,king,
 440      &             mpi_comm_world,ierr)
 441 ctime            write(iout,'(a19,8000i8)') ' REMD synchro itime',
 442 ctime     &                    (itime_all(i),i=1,nodes)
 443             if(traj1file) then
 444              imin_itime=itime_all(1)
 445              do i=2,nodes
 446                if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
 447              enddo
 448              ii_write=(imin_itime-imin_itime_old)/ntwx
 449              imin_itime_old=int(imin_itime/ntwx)*ntwx
 450              write(iout,*) imin_itime,imin_itime_old,ii_write
 451             endif
 452            endif
 453 ctime           call flush(iout)
 454         endif
 455         if(mremdsync .and. mod(itime,nstex).eq.0) then
 456            synflag=.true.
 457            if (me.eq.king .or. .not. out1file)
 458      &      write(iout,*) 'REMD synchro at',itime
 459            call flush(iout)
 460         endif
 461         if(synflag.and..not.end_of_run) then
 462            time02=MPI_WTIME()
 463            synflag=.false.
 464
 465 cd           write(iout,*) 'REMD before',me,t_bath
 466
 467 c           call mpi_gather(t_bath,1,mpi_double_precision,
 468 c     &             remd_t_bath,1,mpi_double_precision,king,
 469 c     &             mpi_comm_world,ierr)
 470            potEcomp(n_ene+1)=t_bath
 471            potEcomp(n_ene+2)=iset
 472            call mpi_gather(potEcomp(0),n_ene+3,mpi_double_precision,
 473      &             remd_ene(0,1),n_ene+3,mpi_double_precision,king,
 474      &             mpi_comm_world,ierr)
 475            if(lmuca) then
 476             call mpi_gather(elow,1,mpi_double_precision,
 477      &             elowi,1,mpi_double_precision,king,
 478      &             mpi_comm_world,ierr)
 479             call mpi_gather(ehigh,1,mpi_double_precision,
 480      &             ehighi,1,mpi_double_precision,king,
 481      &             mpi_comm_world,ierr)
 482            endif
 483
 484           time03=MPI_WTIME()
 485           if (me.eq.king .or. .not. out1file) then
 486             write(iout,*) 'REMD gather times=',time03-time01
 487      &                                        ,time03-time02
 488           endif
 489
 490           if (restart1file) call write1rst
 491
 492           time04=MPI_WTIME()
 493           if (me.eq.king .or. .not. out1file) then
 494             write(iout,*) 'REMD writing rst time=',time04-time03
 495           endif
 496
 497           if (traj1file) call write1traj
 498
 499           time05=MPI_WTIME()
 500           if (me.eq.king .or. .not. out1file) then
 501             write(iout,*) 'REMD writing traj time=',time05-time04
 502           endif
 503
 504
 505           if (me.eq.king) then
 506             do i=1,nodes
 507                remd_t_bath(i)=remd_ene(n_ene+1,i)
 508             enddo
 509             if(lmuca) then
 510 co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
 511              do i=1,nodes
 512                write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
 513      &            elowi(i),ehighi(i)
 514              enddo
 515             else
 516 cd              write(iout,*) 'REMD exchange temp,ene'
 517 c             do i=1,nodes
 518 co              write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
 519 c              write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
 520 c             enddo
 521             endif
 522 c-------------------------------------
 523            do irr=1,remd_m(1)
 524            i=ifirst(iran_num(1,remd_m(1)))
 525            do ii=1,nodes-1
 526
 527             if(i.gt.0.and.nupa(0,i).gt.0) then
 528               iex=nupa(iran_num(1,int(nupa(0,i))),i)
 529              if (lmuca) then
 530                call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
 531              else
 532 c Swap temperatures between conformations i and iex with recalculating the free energies
 533 c following temperature changes.
 534               ene_iex_iex=remd_ene(0,iex)
 535               ene_i_i=remd_ene(0,i)
 536 co              write (iout,*) "rescaling weights with temperature",
 537 co     &          remd_t_bath(i)
 538               call rescale_weights(remd_t_bath(i))
 539               call sum_energy(remd_ene(0,iex))
 540               ene_iex_i=remd_ene(0,iex)
 541 cd              write (iout,*) "ene_iex_i",remd_ene(0,iex)
 542               call sum_energy(remd_ene(0,i))
 543 cd              write (iout,*) "ene_i_i",remd_ene(0,i)
 544 c              write (iout,*) "rescaling weights with temperature",
 545 c     &          remd_t_bath(iex)
 546               if (real(ene_i_i).ne.real(remd_ene(0,i))) then
 547                 write (iout,*) "ERROR: inconsistent energies:",i,
 548      &            ene_i_i,remd_ene(0,i)
 549               endif
 550               call rescale_weights(remd_t_bath(iex))
 551               call sum_energy(remd_ene(0,i))
 552 c              write (iout,*) "ene_i_iex",remd_ene(0,i)
 553               ene_i_iex=remd_ene(0,i)
 554               call sum_energy(remd_ene(0,iex))
 555               if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
 556                 write (iout,*) "ERROR: inconsistent energies:",iex,
 557      &            ene_iex_iex,remd_ene(0,iex)
 558               endif
 559 c              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
 560 c              write (iout,*) "i",i," iex",iex
 561 co              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
 562 co     &           " ene_i_iex",ene_i_iex,
 563 co     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
 564               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
 565      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
 566               delta=-delta
 567 c              write(iout,*) 'delta',delta
 568 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
 569 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
 570 c     &              (remd_t_bath(i)*remd_t_bath(iex))
 571              endif
 572               if (delta .gt. 50.0d0) then
 573                 delta=0.0d0
 574               else
 575 #ifdef OSF
 576                 if(isnan(delta))then
 577                   delta=0.0d0
 578                 else if (delta.lt.-50.0d0) then
 579                   delta=dexp(50.0d0)
 580                 else
 581                 delta=dexp(-delta)
 582               endif
 583 #else
 584                 delta=dexp(-delta)
 585 #endif
 586               endif
 587               iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
 588               xxx=ran_number(0.0d0,1.0d0)
 589 co              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
 590               if (delta .gt. xxx) then
 591                 tmp=remd_t_bath(i)
 592                 remd_t_bath(i)=remd_t_bath(iex)
 593                 remd_t_bath(iex)=tmp
 594                 remd_ene(0,i)=ene_i_iex
 595                 remd_ene(0,iex)=ene_iex_i
 596                 if(lmuca) then
 597                   tmp=elowi(i)
 598                   elowi(i)=elowi(iex)
 599                   elowi(iex)=tmp
 600                   tmp=ehighi(i)
 601                   ehighi(i)=ehighi(iex)
 602                   ehighi(iex)=tmp
 603                 endif
 604
 605
 606                 do k=0,nodes
 607                   itmp=nupa(k,i)
 608                   nupa(k,i)=nupa(k,iex)
 609                   nupa(k,iex)=itmp
 610                   itmp=ndowna(k,i)
 611                   ndowna(k,i)=ndowna(k,iex)
 612                   ndowna(k,iex)=itmp
 613                 enddo
 614                 do il=1,nodes
 615                  if (ifirst(il).eq.i) ifirst(il)=iex
 616                  do k=1,nupa(0,il)
 617                   if (nupa(k,il).eq.i) then
 618                      nupa(k,il)=iex
 619                   elseif (nupa(k,il).eq.iex) then
 620                      nupa(k,il)=i
 621                   endif
 622                  enddo
 623                  do k=1,ndowna(0,il)
 624                   if (ndowna(k,il).eq.i) then
 625                      ndowna(k,il)=iex
 626                   elseif (ndowna(k,il).eq.iex) then
 627                      ndowna(k,il)=i
 628                   endif
 629                  enddo
 630                 enddo
 631
 632                 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
 633                 itmp=i2rep(i-1)
 634                 i2rep(i-1)=i2rep(iex-1)
 635                 i2rep(iex-1)=itmp
 636
 637 cd                write(iout,*) 'exchange',i,iex
 638 cd              write (iout,'(a8,100i4)') "@ ifirst",
 639 cd     &                    (ifirst(k),k=1,remd_m(1))
 640 cd              do il=1,nodes
 641 cd               write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
 642 cd     &                    (nupa(k,il),k=1,nupa(0,il))
 643 cd               write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
 644 cd     &                    (ndowna(k,il),k=1,ndowna(0,il))
 645 cd              enddo
 646
 647               else
 648                remd_ene(0,iex)=ene_iex_iex
 649                remd_ene(0,i)=ene_i_i
 650                i=iex
 651               endif
 652             endif
 653            enddo
 654            enddo
 655 c-------------------------------------
 656              do i=1,nrep
 657               if(iremd_tot(i).ne.0)
 658      &          write(iout,'(a3,i4,2f12.5)') 'ACC',i,remd_t(i)
 659      &           ,iremd_acc(i)/(1.0*iremd_tot(i))
 660              enddo
 661
 662 cd              write (iout,'(a6,100i4)') "ifirst",
 663 cd     &                    (ifirst(i),i=1,remd_m(1))
 664 cd              do il=1,nodes
 665 cd               write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
 666 cd     &                    (nupa(i,il),i=1,nupa(0,il))
 667 cd               write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
 668 cd     &                    (ndowna(i,il),i=1,ndowna(0,il))
 669 cd              enddo
 670           endif
 671
 672          time06=MPI_WTIME()
 673          call mpi_scatter(remd_t_bath,1,mpi_double_precision,
 674      &           t_bath,1,mpi_double_precision,king,
 675      &           mpi_comm_world,ierr)
 676          time07=MPI_WTIME()
 677           if (me.eq.king .or. .not. out1file) then
 678             write(iout,*) 'REMD scatter time=',time07-time06
 679           endif
 680
 681          if(lmuca) then
 682            call mpi_scatter(elowi,1,mpi_double_precision,
 683      &           elow,1,mpi_double_precision,king,
 684      &           mpi_comm_world,ierr)
 685            call mpi_scatter(ehighi,1,mpi_double_precision,
 686      &           ehigh,1,mpi_double_precision,king,
 687      &           mpi_comm_world,ierr)
 688          endif
 689          call rescale_weights(t_bath)
 690 co         write (iout,*) "Processor",me,
 691 co     &    " rescaling weights with temperature",t_bath
 692
 693          stdfp=dsqrt(2*Rb*t_bath/d_time)
 694          do i=1,ntyp
 695            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 696          enddo
 697
 698 cde         write(iout,*) 'REMD after',me,t_bath
 699            time08=MPI_WTIME()
 700            if (me.eq.king .or. .not. out1file) then
 701             write(iout,*) 'REMD exchange time=',time08-time00
 702             call flush(iout)
 703            endif
 704         endif
 705       enddo
 706
 707       if (restart1file) then
 708           if (me.eq.king .or. .not. out1file)
 709      &      write(iout,*) 'writing restart at the end of run'
 710            call write1rst
 711       endif
 712
 713       if (traj1file) call write1traj
 714
 715
 716       t_MD=tcpu()-tt0
 717       if (me.eq.king .or. .not. out1file) then
 718        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
 719      &  '  Timing  ',
 720      & 'MD calculations setup:',t_MDsetup,
 721      & 'Energy & gradient evaluation:',t_enegrad,
 722      & 'Stochastic MD setup:',t_langsetup,
 723      & 'Stochastic MD step setup:',t_sdsetup,
 724      & 'MD steps:',t_MD
 725        write (iout,'(/28(1h=),a25,27(1h=))')
 726      & '  End of MD calculation  '
 727       endif
 728       return
 729       end
 730
 731       subroutine write1rst_
 732       implicit real*8 (a-h,o-z)
 733       include 'DIMENSIONS'
 734       include 'mpif.h'
 735       include 'COMMON.MD'
 736       include 'COMMON.IOUNITS'
 737       include 'COMMON.REMD'
 738       include 'COMMON.SETUP'
 739       include 'COMMON.CHAIN'
 740       include 'COMMON.SBRIDGE'
 741       include 'COMMON.INTERACT'
 742
 743       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
 744      &     d_restart2(3,2*maxres*maxprocs)
 745       real t5_restart1(5)
 746       integer iret,itmp
 747
 748
 749        t5_restart1(1)=totT
 750        t5_restart1(2)=EK
 751        t5_restart1(3)=potE
 752        t5_restart1(4)=t_bath
 753        t5_restart1(5)=Uconst
 754
 755        call mpi_gather(t5_restart1,5,mpi_real,
 756      &      t_restart1,5,mpi_real,king,mpi_comm_world,ierr)
 757
 758
 759        do i=1,2*nres
 760          do j=1,3
 761            r_d(j,i)=d_t(j,i)
 762          enddo
 763        enddo
 764        call mpi_gather(r_d,3*2*nres,mpi_real,
 765      &           d_restart1,3*2*nres,mpi_real,king,
 766      &           mpi_comm_world,ierr)
 767
 768
 769        do i=1,2*nres
 770          do j=1,3
 771            r_d(j,i)=dc(j,i)
 772          enddo
 773        enddo
 774        call mpi_gather(r_d,3*2*nres,mpi_real,
 775      &           d_restart2,3*2*nres,mpi_real,king,
 776      &           mpi_comm_world,ierr)
 777
 778        if(me.eq.king) then
 779          open(irest1,file=mremd_rst_name,status='unknown')
 780          write (irest1,*) (i2rep(i),i=0,nodes-1)
 781          write (irest1,*) (ifirst(i),i=1,remd_m(1))
 782          do il=1,nodes
 783               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 784               write (irest1,*) ndowna(0,il),
 785      &                   (ndowna(i,il),i=1,ndowna(0,il))
 786          enddo
 787
 788          do il=1,nodes
 789            write(irest1,*) (t_restart1(j,il),j=1,4)
 790          enddo
 791
 792          do il=0,nodes-1
 793            do i=1,2*nres
 794              write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
 795            enddo
 796          enddo
 797          do il=0,nodes-1
 798            do i=1,2*nres
 799              write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
 800            enddo
 801          enddo
 802          close(irest1)
 803        endif
 804
 805       return
 806       end
 807
 808       subroutine write1rst
 809       implicit real*8 (a-h,o-z)
 810       include 'DIMENSIONS'
 811       include 'mpif.h'
 812       include 'COMMON.MD'
 813       include 'COMMON.IOUNITS'
 814       include 'COMMON.REMD'
 815       include 'COMMON.SETUP'
 816       include 'COMMON.CHAIN'
 817       include 'COMMON.SBRIDGE'
 818       include 'COMMON.INTERACT'
 819
 820       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
 821      &     d_restart2(3,2*maxres*maxprocs)
 822       real t5_restart1(5)
 823       integer iret,itmp
 824
 825
 826        t5_restart1(1)=totT
 827        t5_restart1(2)=EK
 828        t5_restart1(3)=potE
 829        t5_restart1(4)=t_bath
 830        t5_restart1(5)=Uconst
 831
 832        call mpi_gather(t5_restart1,5,mpi_real,
 833      &      t_restart1,5,mpi_real,king,mpi_comm_world,ierr)
 834
 835
 836        do i=1,2*nres
 837          do j=1,3
 838            r_d(j,i)=d_t(j,i)
 839          enddo
 840        enddo
 841        call mpi_gather(r_d,3*2*nres,mpi_real,
 842      &           d_restart1,3*2*nres,mpi_real,king,
 843      &           mpi_comm_world,ierr)
 844
 845
 846        do i=1,2*nres
 847          do j=1,3
 848            r_d(j,i)=dc(j,i)
 849          enddo
 850        enddo
 851        call mpi_gather(r_d,3*2*nres,mpi_real,
 852      &           d_restart2,3*2*nres,mpi_real,king,
 853      &           mpi_comm_world,ierr)
 854
 855        if(me.eq.king) then
 856 c         open(irest1,file=mremd_rst_name,status='unknown')
 857          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
 858 c         write (irest1,*) (i2rep(i),i=0,nodes-1)
 859          do i=0,nodes-1
 860           call xdrfint(ixdrf, i2rep(i), iret)
 861          enddo
 862 c         write (irest1,*) (ifirst(i),i=1,remd_m(1))
 863          do i=1,remd_m(1)
 864           call xdrfint(ixdrf, ifirst(i), iret)
 865          enddo
 866          do il=1,nodes
 867 c              write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
 868               do i=0,nupa(0,il)
 869                call xdrfint(ixdrf, nupa(i,il), iret)
 870               enddo
 871
 872 c              write (irest1,*) ndowna(0,il),
 873 c     &                   (ndowna(i,il),i=1,ndowna(0,il))
 874               do i=0,ndowna(0,il)
 875                call xdrfint(ixdrf, ndowna(i,il), iret)
 876               enddo
 877          enddo
 878
 879          do il=1,nodes
 880 c           write(irest1,*) (t_restart1(j,il),j=1,4)
 881            do j=1,4
 882             call xdrffloat(ixdrf, t_restart1(j,il), iret)
 883            enddo
 884          enddo
 885
 886          do il=0,nodes-1
 887            do i=1,2*nres
 888 c             write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
 889             do j=1,3
 890              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
 891             enddo
 892            enddo
 893          enddo
 894          do il=0,nodes-1
 895            do i=1,2*nres
 896 c             write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
 897             do j=1,3
 898              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
 899             enddo
 900            enddo
 901          enddo
 902 c         close(irest1)
 903          call xdrfclose(ixdrf, iret)
 904        endif
 905
 906       return
 907       end
 908
 909
 910       subroutine write1traj
 911       implicit real*8 (a-h,o-z)
 912       include 'DIMENSIONS'
 913       include 'mpif.h'
 914       include 'COMMON.MD'
 915       include 'COMMON.IOUNITS'
 916       include 'COMMON.REMD'
 917       include 'COMMON.SETUP'
 918       include 'COMMON.CHAIN'
 919       include 'COMMON.SBRIDGE'
 920       include 'COMMON.INTERACT'
 921
 922       real t5_restart1(5)
 923       integer iret,itmp
 924       real xcoord(3,maxres2+2),prec
 925       real r_qfrag(50),r_qpair(100)
 926       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
 927       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
 928
 929       call mpi_bcast(ii_write,1,mpi_integer,
 930      &           king,mpi_comm_world,ierr)
 931
 932       print *,'traj1file',me,ii_write,ntwx_cache
 933
 934       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
 935       do ii=1,ii_write
 936        t5_restart1(1)=totT_cache(ii)
 937        t5_restart1(2)=EK_cache(ii)
 938        t5_restart1(3)=potE_cache(ii)
 939        t5_restart1(4)=t_bath_cache(ii)
 940        t5_restart1(5)=Uconst_cache(ii)
 941        call mpi_gather(t5_restart1,5,mpi_real,
 942      &      t_restart1,5,mpi_real,king,mpi_comm_world,ierr)
 943
 944           do i=1,nfrag
 945            r_qfrag(i)=qfrag_cache(i,ii)
 946           enddo
 947           do i=1,npair
 948            r_qpair(i)=qpair_cache(i,ii)
 949           enddo
 950
 951         call mpi_gather(r_qfrag,nfrag,mpi_real,
 952      &           p_qfrag,nfrag,mpi_real,king,
 953      &           mpi_comm_world,ierr)
 954         call mpi_gather(r_qpair,npair,mpi_real,
 955      &           p_qpair,npair,mpi_real,king,
 956      &           mpi_comm_world,ierr)
 957
 958         do i=1,nres*2
 959          do j=1,3
 960           r_c(j,i)=c_cache(j,i,ii)
 961          enddo
 962         enddo
 963
 964         call mpi_gather(r_c,3*2*nres,mpi_real,
 965      &           p_c,3*2*nres,mpi_real,king,
 966      &           mpi_comm_world,ierr)
 967
 968        if(me.eq.king) then
 969
 970          do il=1,nodes
 971           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
 972           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
 973           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
 974           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
 975           call xdrfint(ixdrf, nss, iret)
 976           do j=1,nss
 977            call xdrfint(ixdrf, ihpb(j), iret)
 978            call xdrfint(ixdrf, jhpb(j), iret)
 979           enddo
 980           call xdrfint(ixdrf, nfrag+npair, iret)
 981           do i=1,nfrag
 982            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
 983           enddo
 984           do i=1,npair
 985            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
 986           enddo
 987           prec=10000.0
 988           do i=1,nres
 989            do j=1,3
 990             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
 991            enddo
 992           enddo
 993           do i=nnt,nct
 994            do j=1,3
 995             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
 996            enddo
 997           enddo
 998           itmp=nres+nct-nnt+1
 999           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1000          enddo
1001        endif
1002       enddo
1003       if(me.eq.king) call xdrfclose(ixdrf, iret)
1004       do i=1,ntwx_cache-ii_write
1005
1006             totT_cache(i)=totT_cache(ii_write+i)
1007             EK_cache(i)=EK_cache(ii_write+i)
1008             potE_cache(i)=potE_cache(ii_write+i)
1009             t_bath_cache(i)=t_bath_cache(ii_write+i)
1010             Uconst_cache(i)=Uconst_cache(ii_write+i)
1011
1012             do ii=1,nfrag
1013              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1014             enddo
1015             do ii=1,npair
1016              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1017             enddo
1018
1019             do ii=1,nres*2
1020              do j=1,3
1021               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1022              enddo
1023             enddo
1024       enddo
1025       ntwx_cache=ntwx_cache-ii_write
1026       return
1027       end
1028
1029
1030       subroutine read1restart
1031       implicit real*8 (a-h,o-z)
1032       include 'DIMENSIONS'
1033       include 'mpif.h'
1034       include 'COMMON.MD'
1035       include 'COMMON.IOUNITS'
1036       include 'COMMON.REMD'
1037       include 'COMMON.SETUP'
1038       include 'COMMON.CHAIN'
1039       include 'COMMON.SBRIDGE'
1040       include 'COMMON.INTERACT'
1041       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1042      &                 t5_restart1(5)
1043
1044          if(me.eq.king)then
1045               open(irest2,file=mremd_rst_name,status='unknown')
1046               read(irest2,*,err=334) i
1047               write(iout,*) "Reading old rst in ASCI format"
1048               close(irest2)
1049                call read1restart_old
1050                return
1051  334          continue
1052               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1053
1054 c             read (irest2,*) (i2rep(i),i=0,nodes-1)
1055               do i=0,nodes-1
1056                call xdrfint(ixdrf, i2rep(i), iret)
1057               enddo
1058 c             read (irest2,*) (ifirst(i),i=1,remd_m(1))
1059               do i=1,remd_m(1)
1060                call xdrfint(ixdrf, ifirst(i), iret)
1061               enddo
1062              do il=1,nodes
1063 c              read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1064               call xdrfint(ixdrf, nupa(0,il), iret)
1065               do i=1,nupa(0,il)
1066                call xdrfint(ixdrf, nupa(i,il), iret)
1067               enddo
1068
1069 c              read (irest2,*) ndowna(0,il),
1070 c     &                    (ndowna(i,il),i=1,ndowna(0,il))
1071               call xdrfint(ixdrf, ndowna(0,il), iret)
1072               do i=1,ndowna(0,il)
1073                call xdrfint(ixdrf, ndowna(i,il), iret)
1074               enddo
1075              enddo
1076              do il=1,nodes
1077 c               read(irest2,*) (t_restart1(j,il),j=1,4)
1078                do j=1,4
1079                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1080                enddo
1081              enddo
1082          endif
1083          call mpi_scatter(t_restart1,5,mpi_real,
1084      &           t5_restart1,5,mpi_real,king,mpi_comm_world,ierr)
1085          totT=t5_restart1(1)
1086          EK=t5_restart1(2)
1087          potE=t5_restart1(3)
1088          t_bath=t5_restart1(4)
1089
1090          if(me.eq.king)then
1091               do il=0,nodes-1
1092                do i=1,2*nres
1093 c                read(irest2,'(3e15.5)')
1094 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1095             do j=1,3
1096              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1097             enddo
1098                enddo
1099               enddo
1100          endif
1101          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1102      &           r_d,3*2*nres,mpi_real,king,mpi_comm_world,ierr)
1103
1104          do i=1,2*nres
1105            do j=1,3
1106             d_t(j,i)=r_d(j,i)
1107            enddo
1108          enddo
1109          if(me.eq.king)then
1110               do il=0,nodes-1
1111                do i=1,2*nres
1112 c                read(irest2,'(3e15.5)')
1113 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1114             do j=1,3
1115              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1116             enddo
1117                enddo
1118               enddo
1119          endif
1120          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1121      &           r_d,3*2*nres,mpi_real,king,mpi_comm_world,ierr)
1122          do i=1,2*nres
1123            do j=1,3
1124             dc(j,i)=r_d(j,i)
1125            enddo
1126          enddo
1127         if(me.eq.king) close(irest2)
1128         return
1129         end
1130
1131       subroutine read1restart_old
1132       implicit real*8 (a-h,o-z)
1133       include 'DIMENSIONS'
1134       include 'mpif.h'
1135       include 'COMMON.MD'
1136       include 'COMMON.IOUNITS'
1137       include 'COMMON.REMD'
1138       include 'COMMON.SETUP'
1139       include 'COMMON.CHAIN'
1140       include 'COMMON.SBRIDGE'
1141       include 'COMMON.INTERACT'
1142       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1143      &                 t5_restart1(5)
1144
1145          if(me.eq.king)then
1146              open(irest2,file=mremd_rst_name,status='unknown')
1147              read (irest2,*) (i2rep(i),i=0,nodes-1)
1148              read (irest2,*) (ifirst(i),i=1,remd_m(1))
1149              do il=1,nodes
1150               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1151               read (irest2,*) ndowna(0,il),
1152      &                    (ndowna(i,il),i=1,ndowna(0,il))
1153              enddo
1154              do il=1,nodes
1155                read(irest2,*) (t_restart1(j,il),j=1,4)
1156              enddo
1157          endif
1158          call mpi_scatter(t_restart1,5,mpi_real,
1159      &           t5_restart1,5,mpi_real,king,mpi_comm_world,ierr)
1160          totT=t5_restart1(1)
1161          EK=t5_restart1(2)
1162          potE=t5_restart1(3)
1163          t_bath=t5_restart1(4)
1164
1165          if(me.eq.king)then
1166               do il=0,nodes-1
1167                do i=1,2*nres
1168                 read(irest2,'(3e15.5)')
1169      &                (d_restart1(j,i+2*nres*il),j=1,3)
1170                enddo
1171               enddo
1172          endif
1173          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1174      &           r_d,3*2*nres,mpi_real,king,mpi_comm_world,ierr)
1175
1176          do i=1,2*nres
1177            do j=1,3
1178             d_t(j,i)=r_d(j,i)
1179            enddo
1180          enddo
1181          if(me.eq.king)then
1182               do il=0,nodes-1
1183                do i=1,2*nres
1184                 read(irest2,'(3e15.5)')
1185      &                (d_restart1(j,i+2*nres*il),j=1,3)
1186                enddo
1187               enddo
1188          endif
1189          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1190      &           r_d,3*2*nres,mpi_real,king,mpi_comm_world,ierr)
1191          do i=1,2*nres
1192            do j=1,3
1193             dc(j,i)=r_d(j,i)
1194            enddo
1195          enddo
1196         if(me.eq.king) close(irest2)
1197         return
1198         end
1199