MREMD.F.safe

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