3 implicit real*8 (a-h,o-z)
6 include 'COMMON.CONTROL'
10 include 'COMMON.LANGEVIN'
12 include 'COMMON.LANGEVIN.lang0'
14 include 'COMMON.CHAIN'
15 include 'COMMON.DERIV'
17 include 'COMMON.LOCAL'
18 include 'COMMON.INTERACT'
19 include 'COMMON.IOUNITS'
20 include 'COMMON.NAMES'
21 include 'COMMON.TIME1'
23 include 'COMMON.SETUP'
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)
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)
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
51 if(me.eq.king.or..not.out1file) then
52 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
53 write (iout,*) "NREP=",nrep
57 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
58 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
60 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
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
76 do il1=1,max0(mset(il),1)
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)"
100 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
107 c print *,'i2rep',me,i2rep(me)
108 c print *,'rep2i',(rep2i(i),i=0,nrep)
110 cold if (i2rep(me).eq.nrep) then
113 cold nup(0)=remd_m(i2rep(me)+1)
114 cold k=rep2i(int(i2rep(me)))+1
121 cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
123 cold if (i2rep(me).eq.1) then
126 cold ndown(0)=remd_m(i2rep(me)-1)
127 cold k=rep2i(i2rep(me)-2)+1
134 cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
137 write (*,*) "Processor",me," rest",rest,"
138 & restart1fie",restart1file
139 if(rest.and.restart1file) then
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,
145 cd write (*,*) me," After broadcast: file_exist",file_exist
147 if(me.eq.king.or..not.out1file)
148 & write (iout,*) 'Master is reading restart1file'
149 call read1restart(i_index)
151 if(me.eq.king.or..not.out1file)
152 & write (iout,*) 'WARNING : no restart1file'
155 if(me.eq.king.or..not.out1file) then
156 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
157 write(iout,*) "i_index"
162 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
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',
178 open(irest2,file=mremd_rst_name,status='unknown')
180 read (irest2,*) (i2rep(i),i=0,nodes-1)
182 read (irest2,*) (ifirst(i),i=1,remd_m(1))
185 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
187 read (irest2,*) ndowna(0,il),
188 & (ndowna(i,il),i=1,ndowna(0,il))
190 if(usampl.or.hremd.gt.0) then
194 read (irest2,*) (mset(i),i=1,nset)
196 read (irest2,*) (i2set(i),i=0,nodes-1)
201 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
206 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
207 write(iout,*) "i_index"
212 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
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))
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))
230 ELSE IF (.not.(rest.and.file_exist)) THEN
236 if (i2rep(il-1).eq.nrep) then
239 nupa(0,il)=remd_m(i2rep(il-1)+1)
240 k=rep2i(int(i2rep(il-1)))+1
246 if (i2rep(il-1).eq.1) then
249 ndowna(0,il)=remd_m(i2rep(il-1)-1)
250 k=rep2i(i2rep(il-1)-2)+1
258 write (iout,'(a6,100i4)') "ifirst",
259 & (ifirst(i),i=1,remd_m(1))
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))
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
275 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
277 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
278 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
281 if(usampl.or.hremd.gt.0) then
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
288 stdfp=dsqrt(2*Rb*t_bath/d_time)
290 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
293 c print *,'irep',me,t_bath
295 if (me.eq.king .or. .not. out1file)
296 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
297 call rescale_weights(t_bath)
301 c------copy MD--------------
302 c The driver for molecular dynamics subroutines
303 c------------------------------------------------
309 if(me.eq.king.or..not.out1file)
310 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
316 c Determine the inverse of the inertia matrix.
317 call setup_MD_matrices
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)
325 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
327 if (lang.gt.0 .and. .not.surfarea) then
329 stdforcp(i)=stdfp*dsqrt(gamp)
332 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
334 elseif (lang.gt.0 .and. surfarea ) then
337 call rescale_weights(t_bath)
341 t_MDsetup = MPI_Wtime()-tt0
343 t_MDsetup = tcpu()-tt0
346 c Entering the MD loop
352 if (lang.eq.2 .or. lang.eq.3) then
356 call sd_verlet_p_setup
358 call sd_verlet_ciccotti_setup
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)
372 flag_stoch(i)=.false.
376 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
378 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
382 else if (lang.eq.1 .or. lang.eq.4) then
386 if (me.eq.king .or. .not. out1file)
387 & write(iout,*) 'Setup time',time00-walltime
388 ctime call flush(iout)
390 t_langsetup=MPI_Wtime()-tt0
393 t_langsetup=tcpu()-tt0
398 do while(.not.end_of_run)
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.
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
409 call sd_verlet_p_setup
411 call sd_verlet_ciccotti_setup
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)
425 flag_stoch(i)=.false.
428 else if (lang.eq.1 .or. lang.eq.4) then
431 write (iout,'(a,i10)')
432 & "Friction matrix reset based on surface area, itime",itime
434 if (reset_vel .and. tbf .and. lang.eq.0
435 & .and. mod(itime,count_reset_vel).eq.0) then
437 if (me.eq.king .or. .not. out1file)
438 & write(iout,'(a,f20.2)')
439 & "Velocities reset to random values, time",totT
442 d_t_old(j,i)=d_t(j,i)
446 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
450 d_t(j,0)=d_t(j,0)-vcm(j)
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
458 d_t_old(j,i)=scalfac*d_t(j,i)
464 c Time-reversible RESPA algorithm
465 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
466 call RESPA_step(itime)
468 c Variable time step algorithm.
469 call velverlet_step(itime)
473 call brown_step(itime)
475 print *,"Brown dynamics not here!"
477 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
483 if (mod(itime,ntwe).eq.0) call statout(itime)
485 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
486 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
488 call pdbout(potE,tytul,ipdb)
493 if (mod(itime,ntwx).eq.0.and.traj1file) then
494 if(ntwx_cache.lt.max_cache_traj_use) then
495 ntwx_cache=ntwx_cache+1
497 if (max_cache_traj_use.ne.1)
498 & print *,itime,"processor ",me," over cache ",ntwx_cache
501 totT_cache(i)=totT_cache(i+1)
502 EK_cache(i)=EK_cache(i+1)
503 potE_cache(i)=potE_cache(i+1)
504 t_bath_cache(i)=t_bath_cache(i+1)
505 Uconst_cache(i)=Uconst_cache(i+1)
506 iset_cache(i)=iset_cache(i+1)
509 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
512 qpair_cache(ii,i)=qpair_cache(ii,i+1)
515 utheta_cache(ii,i)=utheta_cache(ii,i+1)
516 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
517 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
523 c_cache(j,ii,i)=c_cache(j,ii,i+1)
529 totT_cache(ntwx_cache)=totT
530 EK_cache(ntwx_cache)=EK
531 potE_cache(ntwx_cache)=potE
532 t_bath_cache(ntwx_cache)=t_bath
533 Uconst_cache(ntwx_cache)=Uconst
534 iset_cache(ntwx_cache)=iset
537 qfrag_cache(i,ntwx_cache)=qfrag(i)
540 qpair_cache(i,ntwx_cache)=qpair(i)
543 utheta_cache(i,ntwx_cache)=utheta(i)
544 ugamma_cache(i,ntwx_cache)=ugamma(i)
545 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
550 c_cache(j,i,ntwx_cache)=c(j,i)
555 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
556 & .and..not.restart1file) then
559 open(irest1,file=mremd_rst_name,status='unknown')
560 write (irest1,*) "i2rep"
561 write (irest1,*) (i2rep(i),i=0,nodes-1)
562 write (irest1,*) "ifirst"
563 write (irest1,*) (ifirst(i),i=1,remd_m(1))
565 write (irest1,*) "nupa",il
566 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
567 write (irest1,*) "ndowna",il
568 write (irest1,*) ndowna(0,il),
569 & (ndowna(i,il),i=1,ndowna(0,il))
571 if(usampl.or.hremd.gt.0) then
572 write (irest1,*) "nset"
573 write (irest1,*) nset
574 write (irest1,*) "mset"
575 write (irest1,*) (mset(i),i=1,nset)
576 write (irest1,*) "i2set"
577 write (irest1,*) (i2set(i),i=0,nodes-1)
578 write (irest1,*) "i_index"
582 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
590 open(irest2,file=rest2name,status='unknown')
591 write(irest2,*) totT,EK,potE,totE,t_bath
593 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
596 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
598 if(usampl.or.hremd.gt.0) then
599 write (irest2,*) iset
606 c forced synchronization
607 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
608 & .and. .not. mremdsync) then
610 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
612 call mpi_recv(itime_master, 1, MPI_INTEGER,
613 & 0,101,CG_COMM, status, ierr)
614 call mpi_barrier(CG_COMM, ierr)
615 cdeb if (out1file.or.traj1file) then
616 cdeb call mpi_gather(itime,1,mpi_integer,
617 cdeb & icache_all,1,mpi_integer,king,
620 & call mpi_gather(ntwx_cache,1,mpi_integer,
621 & icache_all,1,mpi_integer,king,
624 & write(iout,*) 'REMD synchro at',itime_master,itime
625 if (itime_master.ge.n_timestep .or. ovrtim())
627 ctime call flush(iout)
632 if ((mod(itime,nstex).eq.0.and.me.eq.king
633 & .or.end_of_run.and.me.eq.king )
634 & .and. .not. mremdsync ) then
638 call mpi_isend(itime,1,MPI_INTEGER,i,101,
639 & CG_COMM, ireqi(i), ierr)
640 cd write(iout,*) 'REMD synchro with',i
643 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
644 call mpi_barrier(CG_COMM, ierr)
646 write(iout,*) 'REMD synchro at',itime,'time=',time01-time01_
647 if (out1file.or.traj1file) then
648 cdeb call mpi_gather(itime,1,mpi_integer,
649 cdeb & itime_all,1,mpi_integer,king,
651 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
652 cdeb & (itime_all(i),i=1,nodes)
654 cdeb imin_itime=itime_all(1)
656 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
658 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
659 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
660 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
661 call mpi_gather(ntwx_cache,1,mpi_integer,
662 & icache_all,1,mpi_integer,king,
664 c write(iout,'(a19,8000i8)') ' ntwx_cache',
665 c & (icache_all(i),i=1,nodes)
666 ii_write=icache_all(1)
668 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
670 c write(iout,*) "MIN ii_write=",ii_write
673 ctime call flush(iout)
675 if(mremdsync .and. mod(itime,nstex).eq.0) then
677 if (me.eq.king .or. .not. out1file)
678 & write(iout,*) 'REMD synchro at',itime
681 call mpi_gather(ntwx_cache,1,mpi_integer,
682 & icache_all,1,mpi_integer,king,
685 write(iout,'(a19,8000i8)') ' ntwx_cache',
686 & (icache_all(i),i=1,nodes)
687 ii_write=icache_all(1)
689 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
691 write(iout,*) "MIN ii_write=",ii_write
694 ctest call flush(iout)
697 c Update the time safety limiy
698 if (time001-time00.gt.safety) then
699 safety=time001-time00+600
700 if (me.eq.king .or. .not. out1file)
701 & write (iout,*) "****** SAFETY increased to",safety," s"
703 if (ovrtim()) end_of_run=.true.
705 if(synflag.and..not.end_of_run) then
709 c write(iout,*) 'REMD before',me,t_bath
711 c call mpi_gather(t_bath,1,mpi_double_precision,
712 c & remd_t_bath,1,mpi_double_precision,king,
714 potEcomp(n_ene+1)=t_bath
717 potEcomp(n_ene+2)=iset
718 if (iset.lt.nset) then
722 potEcomp(n_ene+3)=Uconst
729 potEcomp(n_ene+4)=Uconst
733 if(hremd.gt.0) potEcomp(n_ene+2)=iset
734 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
735 & remd_ene(0,1),n_ene+5,mpi_double_precision,king,
738 call mpi_gather(elow,1,mpi_double_precision,
739 & elowi,1,mpi_double_precision,king,
741 call mpi_gather(ehigh,1,mpi_double_precision,
742 & ehighi,1,mpi_double_precision,king,
747 if (me.eq.king .or. .not. out1file) then
748 write(iout,*) 'REMD gather times=',time03-time01
752 if (restart1file) call write1rst(i_index)
755 if (me.eq.king .or. .not. out1file) then
756 write(iout,*) 'REMD writing rst time=',time04-time03
759 if (traj1file) call write1traj
761 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
762 cdeb & icache_all,1,mpi_integer,king,
764 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
765 cdeb & (icache_all(i),i=1,nodes)
770 if (me.eq.king .or. .not. out1file) then
771 write(iout,*) 'REMD writing traj time=',time05-time04
772 ctime call flush(iout)
778 remd_t_bath(i)=remd_ene(n_ene+1,i)
779 iremd_iset(i)=remd_ene(n_ene+2,i)
782 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
784 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
788 write(iout,*) 'REMD exchange temp,ene'
790 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
791 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
794 c-------------------------------------
795 IF(.not.usampl.and.hremd.eq.0) THEN
796 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
798 ctime call flush(iout)
799 write (iout,*) "remd_m(1)",remd_m(1)
801 i=ifirst(iran_num(1,remd_m(1)))
803 ctime call flush(iout)
807 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
808 if(i.gt.0.and.nupa(0,i).gt.0) then
810 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
812 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
814 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
816 c do while (iex.eq.i)
817 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
818 iex=nupa(iran_num(1,int(nupa(0,i))),i)
820 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
822 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
824 c Swap temperatures between conformations i and iex with recalculating the free energies
825 c following temperature changes.
826 ene_iex_iex=remd_ene(0,iex)
827 ene_i_i=remd_ene(0,i)
828 c write (iout,*) "i",i," ene_i_i",ene_i_i,
829 c & " iex",iex," ene_iex_iex",ene_iex_iex
830 c write (iout,*) "rescaling weights with temperature",
833 call rescale_weights(remd_t_bath(i))
835 c write (iout,*) "0,iex",remd_t_bath(i)
836 c call enerprint(remd_ene(0,iex))
838 call sum_energy(remd_ene(0,iex),.false.)
839 ene_iex_i=remd_ene(0,iex)
840 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
842 c write (iout,*) "0,i",remd_t_bath(i)
843 c call enerprint(remd_ene(0,i))
845 call sum_energy(remd_ene(0,i),.false.)
846 c write (iout,*) "ene_i_i",remd_ene(0,i)
848 c write (iout,*) "rescaling weights with temperature",
850 if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
851 write (iout,*) "ERROR: inconsistent energies:",i,
852 & ene_i_i,remd_ene(0,i)
854 call rescale_weights(remd_t_bath(iex))
856 c write (iout,*) "0,i",remd_t_bath(iex)
857 c call enerprint(remd_ene(0,i))
859 call sum_energy(remd_ene(0,i),.false.)
860 c write (iout,*) "ene_i_iex",remd_ene(0,i)
862 ene_i_iex=remd_ene(0,i)
864 c write (iout,*) "0,iex",remd_t_bath(iex)
865 c call enerprint(remd_ene(0,iex))
867 call sum_energy(remd_ene(0,iex),.false.)
868 if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
869 write (iout,*) "ERROR: inconsistent energies:",iex,
870 & ene_iex_iex,remd_ene(0,iex)
872 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
873 c write (iout,*) "i",i," iex",iex
874 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
875 c & " ene_i_iex",ene_i_iex,
876 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
878 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
879 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
881 c write(iout,*) 'delta',delta
882 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
883 c & (remd_ene(i)-remd_ene(iex))/Rb/
884 c & (remd_t_bath(i)*remd_t_bath(iex))
886 if (delta .gt. 50.0d0) then
892 else if (delta.lt.-50.0d0) then
901 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
902 xxx=ran_number(0.0d0,1.0d0)
903 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
905 if (delta .gt. xxx) then
907 remd_t_bath(i)=remd_t_bath(iex)
909 remd_ene(0,i)=ene_i_iex
910 remd_ene(0,iex)=ene_iex_i
916 ehighi(i)=ehighi(iex)
923 nupa(k,i)=nupa(k,iex)
926 ndowna(k,i)=ndowna(k,iex)
930 if (ifirst(il).eq.i) ifirst(il)=iex
932 if (nupa(k,il).eq.i) then
934 elseif (nupa(k,il).eq.iex) then
939 if (ndowna(k,il).eq.i) then
941 elseif (ndowna(k,il).eq.iex) then
947 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
949 i2rep(i-1)=i2rep(iex-1)
952 c write(iout,*) 'exchange',i,iex
953 c write (iout,'(a8,100i4)') "@ ifirst",
954 c & (ifirst(k),k=1,remd_m(1))
956 c write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
957 c & (nupa(k,il),k=1,nupa(0,il))
958 c write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
959 c & (ndowna(k,il),k=1,ndowna(0,il))
964 remd_ene(0,iex)=ene_iex_iex
965 remd_ene(0,i)=ene_i_i
971 cd write (iout,*) "exchange completed"
975 cd write(iout,*) "########",ii
977 i_temp=iran_num(1,nrep)
978 i_mult=iran_num(1,remd_m(i_temp))
979 i_iset=iran_num(1,nset)
980 i_mset=iran_num(1,mset(i_iset))
981 i=i_index(i_temp,i_mult,i_iset,i_mset)
983 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
985 if(t_exchange_only)then
990 cd write(iout,*) "i_dir=",i_dir
992 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
995 i_mult1=iran_num(1,remd_m(i_temp1))
997 i_mset1=iran_num(1,mset(i_iset1))
998 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1000 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1003 i_mult1=iran_num(1,remd_m(i_temp1))
1005 i_mset1=iran_num(1,mset(i_iset1))
1006 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1007 econstr_temp_i=remd_ene(20,i)
1008 econstr_temp_iex=remd_ene(20,iex)
1009 remd_ene(20,i)=remd_ene(n_ene+3,i)
1010 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1012 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1015 i_mult1=iran_num(1,remd_m(i_temp1))
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)
1028 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1029 ctime call flush(iout)
1031 c Swap temperatures between conformations i and iex with recalculating the free energies
1032 c following temperature changes.
1033 ene_iex_iex=remd_ene(0,iex)
1034 ene_i_i=remd_ene(0,i)
1035 co write (iout,*) "rescaling weights with temperature",
1037 call rescale_weights(remd_t_bath(i))
1039 call sum_energy(remd_ene(0,iex),.false.)
1040 ene_iex_i=remd_ene(0,iex)
1041 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
1042 c call sum_energy(remd_ene(0,i),.false.)
1043 cd write (iout,*) "ene_i_i",remd_ene(0,i)
1044 c write (iout,*) "rescaling weights with temperature",
1045 c & remd_t_bath(iex)
1046 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1047 c write (iout,*) "ERROR: inconsistent energies:",i,
1048 c & ene_i_i,remd_ene(0,i)
1050 call rescale_weights(remd_t_bath(iex))
1051 call sum_energy(remd_ene(0,i),.false.)
1052 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1053 ene_i_iex=remd_ene(0,i)
1054 c call sum_energy(remd_ene(0,iex),.false.)
1055 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1056 c write (iout,*) "ERROR: inconsistent energies:",iex,
1057 c & ene_iex_iex,remd_ene(0,iex)
1059 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1060 c write (iout,*) "i",i," iex",iex
1061 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1062 cd & " ene_i_iex",ene_i_iex,
1063 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1064 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1065 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1067 cd write(iout,*) 'delta',delta
1068 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1069 c & (remd_ene(i)-remd_ene(iex))/Rb/
1070 c & (remd_t_bath(i)*remd_t_bath(iex))
1071 if (delta .gt. 50.0d0) then
1076 if (i_dir.eq.1.or.i_dir.eq.3)
1077 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1078 if (i_dir.eq.2.or.i_dir.eq.3)
1079 & iremd_tot_usa(int(i2set(i-1)))=
1080 & iremd_tot_usa(int(i2set(i-1)))+1
1081 xxx=ran_number(0.0d0,1.0d0)
1082 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1083 if (delta .gt. xxx) then
1085 remd_t_bath(i)=remd_t_bath(iex)
1086 remd_t_bath(iex)=tmp
1089 iremd_iset(i)=iremd_iset(iex)
1090 iremd_iset(iex)=itmp
1092 remd_ene(0,i)=ene_i_iex
1093 remd_ene(0,iex)=ene_iex_i
1095 if (i_dir.eq.1.or.i_dir.eq.3)
1096 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1099 i2rep(i-1)=i2rep(iex-1)
1102 if (i_dir.eq.2.or.i_dir.eq.3)
1103 & iremd_acc_usa(int(i2set(i-1)))=
1104 & iremd_acc_usa(int(i2set(i-1)))+1
1107 i2set(i-1)=i2set(iex-1)
1110 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1111 i_index(i_temp,i_mult,i_iset,i_mset)=
1112 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1113 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1116 remd_ene(0,iex)=ene_iex_iex
1117 remd_ene(0,i)=ene_i_i
1118 remd_ene(20,iex)=econstr_temp_iex
1119 remd_ene(20,i)=econstr_temp_i
1123 cd do il1=1,mset(il)
1126 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1136 ELSEIF (hremd.gt.0) THEN
1138 cd write(iout,*) "########",ii
1140 i_temp=iran_num(1,nrep)
1141 i_mult=iran_num(1,remd_m(i_temp))
1142 i_iset=iran_num(1,nset)
1144 i=i_index(i_temp,i_mult,i_iset,i_mset)
1146 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1148 if(t_exchange_only)then
1154 cd write(iout,*) "i_dir=",i_dir
1156 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1159 i_mult1=iran_num(1,remd_m(i_temp1))
1162 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1164 elseif(i_dir.eq.2)then
1167 i_mult1=iran_num(1,remd_m(i_temp1))
1168 i_iset1=iran_num(1,hremd)
1169 do while(i_iset1.eq.i_iset)
1170 i_iset1=iran_num(1,hremd)
1173 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1175 elseif(remd_m(i_temp+1).gt.0)then
1178 i_mult1=iran_num(1,remd_m(i_temp1))
1179 i_iset1=iran_num(1,hremd)
1180 do while(i_iset1.eq.i_iset)
1181 i_iset1=iran_num(1,hremd)
1184 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1190 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1191 ctime call flush(iout)
1193 c Swap temperatures between conformations i and iex with recalculating the free energies
1194 c following temperature changes.
1195 ene_iex_iex=remd_ene(0,iex)
1196 ene_i_i=remd_ene(0,i)
1198 call set_hweights(i_iset)
1199 call rescale_weights(remd_t_bath(i))
1200 call sum_energy(remd_ene(0,iex),.false.)
1201 ene_iex_i=remd_ene(0,iex)
1203 call set_hweights(i_iset1)
1204 call rescale_weights(remd_t_bath(iex))
1205 call sum_energy(remd_ene(0,i),.false.)
1206 ene_i_iex=remd_ene(0,i)
1208 cd write(iout,*) ene_iex_iex,ene_i_i,ene_iex_i,ene_i_iex
1210 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1211 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1214 if (delta .gt. 50.0d0) then
1220 if (i_dir.eq.1.or.i_dir.eq.3)
1221 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1222 if (i_dir.eq.2.or.i_dir.eq.3)
1223 & iremd_tot_usa(int(i2set(i-1)))=
1224 & iremd_tot_usa(int(i2set(i-1)))+1
1225 xxx=ran_number(0.0d0,1.0d0)
1226 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1227 if (delta .gt. xxx) then
1229 cd write (iout,*) "exchange"
1231 remd_t_bath(i)=remd_t_bath(iex)
1232 remd_t_bath(iex)=tmp
1235 iremd_iset(i)=iremd_iset(iex)
1236 iremd_iset(iex)=itmp
1238 remd_ene(0,i)=ene_i_iex
1239 remd_ene(0,iex)=ene_iex_i
1241 if (i_dir.eq.1.or.i_dir.eq.3)
1242 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1245 i2rep(i-1)=i2rep(iex-1)
1248 if (i_dir.eq.2.or.i_dir.eq.3)
1249 & iremd_acc_usa(int(i2set(i-1)))=
1250 & iremd_acc_usa(int(i2set(i-1)))+1
1253 i2set(i-1)=i2set(iex-1)
1256 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1257 i_index(i_temp,i_mult,i_iset,i_mset)=
1258 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1259 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1262 cd do il1=1,mset(il)
1265 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1272 remd_ene(0,iex)=ene_iex_iex
1273 remd_ene(0,i)=ene_i_i
1284 c-------------------------------------
1285 write (iout,*) "NREP",nrep
1287 if(iremd_tot(i).ne.0)
1288 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1289 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1294 if(iremd_tot_usa(i).ne.0)
1295 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1296 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1302 if(iremd_tot_usa(i).ne.0)
1303 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_hremd',i,
1304 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1309 ctime call flush(iout)
1311 cd write (iout,'(a6,100i4)') "ifirst",
1312 cd & (ifirst(i),i=1,remd_m(1))
1314 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1315 cd & (nupa(i,il),i=1,nupa(0,il))
1316 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1317 cd & (ndowna(i,il),i=1,ndowna(0,il))
1322 cd write (iout,*) "Before scatter"
1324 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1325 & t_bath,1,mpi_double_precision,king,
1327 cd write (iout,*) "After scatter"
1329 if(usampl.or.hremd.gt.0)
1330 & call mpi_scatter(iremd_iset,1,mpi_integer,
1331 & iset,1,mpi_integer,king,
1335 if (me.eq.king .or. .not. out1file) then
1336 write(iout,*) 'REMD scatter time=',time07-time06
1340 call mpi_scatter(elowi,1,mpi_double_precision,
1341 & elow,1,mpi_double_precision,king,
1343 call mpi_scatter(ehighi,1,mpi_double_precision,
1344 & ehigh,1,mpi_double_precision,king,
1348 if(hremd.gt.0) call set_hweights(iset)
1349 call rescale_weights(t_bath)
1350 co write (iout,*) "Processor",me,
1351 co & " rescaling weights with temperature",t_bath
1353 stdfp=dsqrt(2*Rb*t_bath/d_time)
1355 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1359 stdforcp(i)=stdforcp(i)*sqrt(t_bath/t_bath_old)
1362 stdforcsc(i)=stdforcsc(i)*sqrt(t_bath/t_bath_old)
1365 cde write(iout,*) 'REMD after',me,t_bath
1367 if (me.eq.king .or. .not. out1file) then
1368 write(iout,*) 'REMD exchange time=',time08-time02
1369 ctime call flush(iout)
1374 if (restart1file) then
1375 if (me.eq.king .or. .not. out1file)
1376 & write(iout,*) 'writing restart at the end of run'
1377 call write1rst(i_index)
1380 if (traj1file) call write1traj
1382 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1383 cdeb & icache_all,1,mpi_integer,king,
1384 cdeb & CG_COMM,ierr)
1385 cdeb write(iout,'(a40,8000i8)')
1386 cdeb & ' ntwx_cache after traj1file at the end',
1387 cdeb & (icache_all(i),i=1,nodes)
1392 t_MD=MPI_Wtime()-tt0
1396 if (me.eq.king .or. .not. out1file) then
1397 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1399 & 'MD calculations setup:',t_MDsetup,
1400 & 'Energy & gradient evaluation:',t_enegrad,
1401 & 'Stochastic MD setup:',t_langsetup,
1402 & 'Stochastic MD step setup:',t_sdsetup,
1404 write (iout,'(/28(1h=),a25,27(1h=))')
1405 & ' End of MD calculation '
1410 c-----------------------------------------------------------------------
1411 subroutine write1rst(i_index)
1412 implicit real*8 (a-h,o-z)
1413 include 'DIMENSIONS'
1416 include 'COMMON.IOUNITS'
1417 include 'COMMON.REMD'
1418 include 'COMMON.SETUP'
1419 include 'COMMON.CHAIN'
1420 include 'COMMON.SBRIDGE'
1421 include 'COMMON.INTERACT'
1423 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1424 & d_restart2(3,2*maxres*maxprocs)
1428 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1429 common /przechowalnia/ d_restart1,d_restart2
1434 t5_restart1(4)=t_bath
1435 t5_restart1(5)=Uconst
1437 call mpi_gather(t5_restart1,5,mpi_real,
1438 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1446 call mpi_gather(r_d,3*2*nres,mpi_real,
1447 & d_restart1,3*2*nres,mpi_real,king,
1456 call mpi_gather(r_d,3*2*nres,mpi_real,
1457 & d_restart2,3*2*nres,mpi_real,king,
1462 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1464 call xdrfint_(ixdrf, i2rep(i), iret)
1467 call xdrfint_(ixdrf, ifirst(i), iret)
1471 call xdrfint_(ixdrf, nupa(i,il), iret)
1475 call xdrfint_(ixdrf, ndowna(i,il), iret)
1481 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1488 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1495 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1501 call xdrfint_(ixdrf, nset, iret)
1503 call xdrfint_(ixdrf,mset(i), iret)
1506 call xdrfint_(ixdrf,i2set(i), iret)
1512 itmp=i_index(i,j,il,il1)
1513 call xdrfint_(ixdrf,itmp, iret)
1520 call xdrfclose_(ixdrf, iret)
1522 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1524 call xdrfint(ixdrf, i2rep(i), iret)
1527 call xdrfint(ixdrf, ifirst(i), iret)
1531 call xdrfint(ixdrf, nupa(i,il), iret)
1535 call xdrfint(ixdrf, ndowna(i,il), iret)
1541 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1548 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1555 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1562 call xdrfint(ixdrf, nset, iret)
1564 call xdrfint(ixdrf,mset(i), iret)
1567 call xdrfint(ixdrf,i2set(i), iret)
1573 itmp=i_index(i,j,il,il1)
1574 call xdrfint(ixdrf,itmp, iret)
1581 call xdrfclose(ixdrf, iret)
1588 subroutine write1traj
1589 implicit real*8 (a-h,o-z)
1590 include 'DIMENSIONS'
1593 include 'COMMON.IOUNITS'
1594 include 'COMMON.REMD'
1595 include 'COMMON.SETUP'
1596 include 'COMMON.CHAIN'
1597 include 'COMMON.SBRIDGE'
1598 include 'COMMON.INTERACT'
1602 real xcoord(3,maxres2+2),prec
1603 real r_qfrag(50),r_qpair(100)
1604 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1605 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1606 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1607 & p_uscdiff(100*maxprocs)
1608 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1609 common /przechowalnia/ p_c
1611 call mpi_bcast(ii_write,1,mpi_integer,
1612 & king,CG_COMM,ierr)
1615 print *,'traj1file',me,ii_write,ntwx_cache
1619 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1621 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1624 t5_restart1(1)=totT_cache(ii)
1625 t5_restart1(2)=EK_cache(ii)
1626 t5_restart1(3)=potE_cache(ii)
1627 t5_restart1(4)=t_bath_cache(ii)
1628 t5_restart1(5)=Uconst_cache(ii)
1629 call mpi_gather(t5_restart1,5,mpi_real,
1630 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1632 call mpi_gather(iset_cache(ii),1,mpi_integer,
1633 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1636 r_qfrag(i)=qfrag_cache(i,ii)
1639 r_qpair(i)=qpair_cache(i,ii)
1642 r_utheta(i)=utheta_cache(i,ii)
1643 r_ugamma(i)=ugamma_cache(i,ii)
1644 r_uscdiff(i)=uscdiff_cache(i,ii)
1647 call mpi_gather(r_qfrag,nfrag,mpi_real,
1648 & p_qfrag,nfrag,mpi_real,king,
1650 call mpi_gather(r_qpair,npair,mpi_real,
1651 & p_qpair,npair,mpi_real,king,
1653 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1654 & p_utheta,nfrag_back,mpi_real,king,
1656 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1657 & p_ugamma,nfrag_back,mpi_real,king,
1659 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1660 & p_uscdiff,nfrag_back,mpi_real,king,
1664 write (iout,*) "p_qfrag"
1666 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1668 write (iout,*) "p_qpair"
1670 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1672 ctime call flush(iout)
1676 r_c(j,i)=c_cache(j,i,ii)
1680 call mpi_gather(r_c,3*2*nres,mpi_real,
1681 & p_c,3*2*nres,mpi_real,king,
1687 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1688 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1689 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1690 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1691 call xdrfint_(ixdrf, nss, iret)
1693 call xdrfint_(ixdrf, ihpb(j), iret)
1694 call xdrfint_(ixdrf, jhpb(j), iret)
1696 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1697 call xdrfint_(ixdrf, iset_restart1(il), iret)
1699 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1702 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1705 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1706 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1707 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1712 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1717 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1721 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1725 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1726 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1727 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1728 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1729 call xdrfint(ixdrf, nss, iret)
1731 call xdrfint(ixdrf, ihpb(j), iret)
1732 call xdrfint(ixdrf, jhpb(j), iret)
1734 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1735 call xdrfint(ixdrf, iset_restart1(il), iret)
1737 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1740 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1743 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1744 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1745 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1750 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1755 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1759 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1765 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1767 if(me.eq.king) call xdrfclose(ixdrf, iret)
1769 do i=1,ntwx_cache-ii_write
1771 totT_cache(i)=totT_cache(ii_write+i)
1772 EK_cache(i)=EK_cache(ii_write+i)
1773 potE_cache(i)=potE_cache(ii_write+i)
1774 t_bath_cache(i)=t_bath_cache(ii_write+i)
1775 Uconst_cache(i)=Uconst_cache(ii_write+i)
1776 iset_cache(i)=iset_cache(ii_write+i)
1779 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1782 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1785 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1786 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1787 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1792 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1796 ntwx_cache=ntwx_cache-ii_write
1801 subroutine read1restart(i_index)
1802 implicit real*8 (a-h,o-z)
1803 include 'DIMENSIONS'
1806 include 'COMMON.IOUNITS'
1807 include 'COMMON.REMD'
1808 include 'COMMON.SETUP'
1809 include 'COMMON.CHAIN'
1810 include 'COMMON.SBRIDGE'
1811 include 'COMMON.INTERACT'
1812 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1815 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1816 common /przechowalnia/ d_restart1
1817 write (*,*) "Processor",me," called read1restart"
1820 open(irest2,file=mremd_rst_name,status='unknown')
1821 read(irest2,*,err=334) i
1822 write(iout,*) "Reading old rst in ASCI format"
1824 call read1restart_old
1828 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1831 call xdrfint_(ixdrf, i2rep(i), iret)
1834 call xdrfint_(ixdrf, ifirst(i), iret)
1837 call xdrfint_(ixdrf, nupa(0,il), iret)
1839 call xdrfint_(ixdrf, nupa(i,il), iret)
1842 call xdrfint_(ixdrf, ndowna(0,il), iret)
1844 call xdrfint_(ixdrf, ndowna(i,il), iret)
1849 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1853 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1856 call xdrfint(ixdrf, i2rep(i), iret)
1859 call xdrfint(ixdrf, ifirst(i), iret)
1862 call xdrfint(ixdrf, nupa(0,il), iret)
1864 call xdrfint(ixdrf, nupa(i,il), iret)
1867 call xdrfint(ixdrf, ndowna(0,il), iret)
1869 call xdrfint(ixdrf, ndowna(i,il), iret)
1874 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1879 call mpi_scatter(t_restart1,5,mpi_real,
1880 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1884 t_bath=t5_restart1(4)
1889 c read(irest2,'(3e15.5)')
1890 c & (d_restart1(j,i+2*nres*il),j=1,3)
1893 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1895 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1901 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1902 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1912 c read(irest2,'(3e15.5)')
1913 c & (d_restart1(j,i+2*nres*il),j=1,3)
1916 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1918 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1924 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1925 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1936 call xdrfint_(ixdrf, nset, iret)
1938 call xdrfint_(ixdrf,mset(i), iret)
1941 call xdrfint_(ixdrf,i2set(i), iret)
1947 call xdrfint_(ixdrf,itmp, iret)
1948 i_index(i,j,il,il1)=itmp
1956 call xdrfint(ixdrf, nset, iret)
1958 call xdrfint(ixdrf,mset(i), iret)
1961 call xdrfint(ixdrf,i2set(i), iret)
1967 call xdrfint(ixdrf,itmp, iret)
1968 i_index(i,j,il,il1)=itmp
1975 call mpi_scatter(i2set,1,mpi_integer,
1976 & iset,1,mpi_integer,king,
1982 if(me.eq.king) close(irest2)
1986 subroutine read1restart_old
1987 implicit real*8 (a-h,o-z)
1988 include 'DIMENSIONS'
1991 include 'COMMON.IOUNITS'
1992 include 'COMMON.REMD'
1993 include 'COMMON.SETUP'
1994 include 'COMMON.CHAIN'
1995 include 'COMMON.SBRIDGE'
1996 include 'COMMON.INTERACT'
1997 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1999 common /przechowalnia/ d_restart1
2001 open(irest2,file=mremd_rst_name,status='unknown')
2002 read (irest2,*) (i2rep(i),i=0,nodes-1)
2003 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2005 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2006 read (irest2,*) ndowna(0,il),
2007 & (ndowna(i,il),i=1,ndowna(0,il))
2010 read(irest2,*) (t_restart1(j,il),j=1,4)
2013 call mpi_scatter(t_restart1,5,mpi_real,
2014 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2018 t_bath=t5_restart1(4)
2023 read(irest2,'(3e15.5)')
2024 & (d_restart1(j,i+2*nres*il),j=1,3)
2028 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2029 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2039 read(irest2,'(3e15.5)')
2040 & (d_restart1(j,i+2*nres*il),j=1,3)
2044 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2045 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2051 if(me.eq.king) close(irest2)
2054 c-------------------------------------------------------------------
2055 subroutine set_hweights(iiset)
2056 implicit real*8 (a-h,o-z)
2058 include 'DIMENSIONS'
2059 include 'COMMON.FFIELD'
2060 include 'COMMON.REMD'
2063 weights(i)=hweights(iiset,i)