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