source/unres/src_MD/MREMD.F

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