source/unres/src_MD/MREMD.F

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