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,homol_nset
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
73 if(homol_nset.gt.1) then
81 if(usampl) i_econstr=20
86 do il1=1,max0(mset(il),1)
102 if(me.eq.king.or..not.out1file) then
103 write(iout,*) "i2rep",(i2rep(i),i=0,nodes-1)
104 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
105 write(iout,*) "i,j,il,il1,i_index(i,j,il,il1)"
110 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
117 c print *,'i2rep',me,i2rep(me)
118 c print *,'rep2i',(rep2i(i),i=0,nrep)
120 cold if (i2rep(me).eq.nrep) then
123 cold nup(0)=remd_m(i2rep(me)+1)
124 cold k=rep2i(int(i2rep(me)))+1
131 cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
133 cold if (i2rep(me).eq.1) then
136 cold ndown(0)=remd_m(i2rep(me)-1)
137 cold k=rep2i(i2rep(me)-2)+1
144 cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
147 write (*,*) "Processor",me," rest",rest,"
148 & restart1fie",restart1file
149 if(rest.and.restart1file) then
151 & inquire(file=mremd_rst_name,exist=file_exist)
152 cd write (*,*) me," Before broadcast: file_exist",file_exist
153 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
155 cd write (*,*) me," After broadcast: file_exist",file_exist
157 if(me.eq.king.or..not.out1file)
158 & write (iout,*) 'Master is reading restart1file'
159 call read1restart(i_index)
161 if(me.eq.king.or..not.out1file)
162 & write (iout,*) 'WARNING : no restart1file'
165 if(me.eq.king.or..not.out1file) then
166 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
167 write(iout,*) "i_index"
172 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
181 if (rest.and..not.restart1file)
182 & inquire(file=mremd_rst_name,exist=file_exist)
183 if(.not.file_exist.and.rest.and..not.restart1file)
184 & write(iout,*) 'WARNING : no restart file',mremd_rst_name
185 IF (rest.and.file_exist.and..not.restart1file) THEN
186 write (iout,*) 'Master is reading restart file',
188 open(irest2,file=mremd_rst_name,status='unknown')
190 read (irest2,*) (i2rep(i),i=0,nodes-1)
192 read (irest2,*) (ifirst(i),i=1,remd_m(1))
195 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
197 read (irest2,*) ndowna(0,il),
198 & (ndowna(i,il),i=1,ndowna(0,il))
200 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
204 read (irest2,*) (mset(i),i=1,nset)
206 read (irest2,*) (i2set(i),i=0,nodes-1)
211 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
216 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
217 write(iout,*) "i_index"
222 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
231 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
232 write (iout,'(a6,1000i5)') "ifirst",
233 & (ifirst(i),i=1,remd_m(1))
235 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
236 & (nupa(i,il),i=1,nupa(0,il))
237 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
238 & (ndowna(i,il),i=1,ndowna(0,il))
240 ELSE IF (.not.(rest.and.file_exist)) THEN
246 if (i2rep(il-1).eq.nrep) then
249 nupa(0,il)=remd_m(i2rep(il-1)+1)
250 k=rep2i(int(i2rep(il-1)))+1
256 if (i2rep(il-1).eq.1) then
259 ndowna(0,il)=remd_m(i2rep(il-1)-1)
260 k=rep2i(i2rep(il-1)-2)+1
268 write (iout,'(a6,100i4)') "ifirst",
269 & (ifirst(i),i=1,remd_m(1))
271 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
272 & (nupa(i,il),i=1,nupa(0,il))
273 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
274 & (ndowna(i,il),i=1,ndowna(0,il))
280 c t_bath=retmin+(retmax-retmin)*me/(nodes-1)
281 if(.not.(rest.and.file_exist.and.restart1file)) then
282 if (me .eq. king) then
285 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
287 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
288 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
291 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
293 if (hremd.gt.0) call set_hweights(iset)
294 if(me.eq.king.or..not.out1file)
295 & write(iout,*) me,"iset=",iset,"t_bath=",t_bath
298 stdfp=dsqrt(2*Rb*t_bath/d_time)
300 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
303 c print *,'irep',me,t_bath
305 if (me.eq.king .or. .not. out1file)
306 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
307 call rescale_weights(t_bath)
311 c------copy MD--------------
312 c The driver for molecular dynamics subroutines
313 c------------------------------------------------
319 if(me.eq.king.or..not.out1file)
320 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
326 c Determine the inverse of the inertia matrix.
327 call setup_MD_matrices
331 if (me.eq.king .or. .not. out1file)
332 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
333 stdfp=dsqrt(2*Rb*t_bath/d_time)
335 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
337 if (lang.gt.0 .and. .not.surfarea) then
339 stdforcp(i)=stdfp*dsqrt(gamp)
342 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
344 elseif (lang.gt.0 .and. surfarea ) then
347 call rescale_weights(t_bath)
351 t_MDsetup = MPI_Wtime()-tt0
353 t_MDsetup = tcpu()-tt0
356 c Entering the MD loop
362 if (lang.eq.2 .or. lang.eq.3) then
366 call sd_verlet_p_setup
368 call sd_verlet_ciccotti_setup
372 pfric0_mat(i,j,0)=pfric_mat(i,j)
373 afric0_mat(i,j,0)=afric_mat(i,j)
374 vfric0_mat(i,j,0)=vfric_mat(i,j)
375 prand0_mat(i,j,0)=prand_mat(i,j)
376 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
377 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
382 flag_stoch(i)=.false.
386 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
388 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
392 else if (lang.eq.1 .or. lang.eq.4) then
396 if (me.eq.king .or. .not. out1file)
397 & write(iout,*) 'Setup time',time00-walltime
398 ctime call flush(iout)
400 t_langsetup=MPI_Wtime()-tt0
403 t_langsetup=tcpu()-tt0
408 do while(.not.end_of_run)
410 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
411 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
413 if (lang.gt.0 .and. surfarea .and.
414 & mod(itime,reset_fricmat).eq.0) then
415 if (lang.eq.2 .or. lang.eq.3) then
419 call sd_verlet_p_setup
421 call sd_verlet_ciccotti_setup
425 pfric0_mat(i,j,0)=pfric_mat(i,j)
426 afric0_mat(i,j,0)=afric_mat(i,j)
427 vfric0_mat(i,j,0)=vfric_mat(i,j)
428 prand0_mat(i,j,0)=prand_mat(i,j)
429 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
430 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
435 flag_stoch(i)=.false.
438 else if (lang.eq.1 .or. lang.eq.4) then
441 write (iout,'(a,i10)')
442 & "Friction matrix reset based on surface area, itime",itime
444 if (reset_vel .and. tbf .and. lang.eq.0
445 & .and. mod(itime,count_reset_vel).eq.0) then
447 if (me.eq.king .or. .not. out1file)
448 & write(iout,'(a,f20.2)')
449 & "Velocities reset to random values, time",totT
452 d_t_old(j,i)=d_t(j,i)
456 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
460 d_t(j,0)=d_t(j,0)-vcm(j)
463 kinetic_T=2.0d0/(dimen3*Rb)*EK
464 scalfac=dsqrt(T_bath/kinetic_T)
465 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
468 d_t_old(j,i)=scalfac*d_t(j,i)
474 c Time-reversible RESPA algorithm
475 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
476 call RESPA_step(itime)
478 c Variable time step algorithm.
479 call velverlet_step(itime)
483 call brown_step(itime)
485 print *,"Brown dynamics not here!"
487 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
492 if(hmc.gt.0 .and. mod(itime,hmc).eq.0) then
497 if (mod(itime,ntwe).eq.0) call statout(itime)
499 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
500 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
502 call pdbout(potE,tytul,ipdb)
507 if (mod(itime,ntwx).eq.0.and.traj1file) then
508 if(ntwx_cache.lt.max_cache_traj_use) then
509 ntwx_cache=ntwx_cache+1
511 if (max_cache_traj_use.ne.1)
512 & print *,itime,"processor ",me," over cache ",ntwx_cache
515 totT_cache(i)=totT_cache(i+1)
516 EK_cache(i)=EK_cache(i+1)
517 potE_cache(i)=potE_cache(i+1)
518 t_bath_cache(i)=t_bath_cache(i+1)
519 Uconst_cache(i)=Uconst_cache(i+1)
520 iset_cache(i)=iset_cache(i+1)
523 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
526 qpair_cache(ii,i)=qpair_cache(ii,i+1)
529 utheta_cache(ii,i)=utheta_cache(ii,i+1)
530 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
531 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
537 c_cache(j,ii,i)=c_cache(j,ii,i+1)
543 totT_cache(ntwx_cache)=totT
544 EK_cache(ntwx_cache)=EK
545 potE_cache(ntwx_cache)=potE
546 t_bath_cache(ntwx_cache)=t_bath
547 Uconst_cache(ntwx_cache)=Uconst
548 iset_cache(ntwx_cache)=iset
551 qfrag_cache(i,ntwx_cache)=qfrag(i)
554 qpair_cache(i,ntwx_cache)=qpair(i)
557 utheta_cache(i,ntwx_cache)=utheta(i)
558 ugamma_cache(i,ntwx_cache)=ugamma(i)
559 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
564 c_cache(j,i,ntwx_cache)=c(j,i)
569 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
570 & .and..not.restart1file) then
573 open(irest1,file=mremd_rst_name,status='unknown')
574 write (irest1,*) "i2rep"
575 write (irest1,*) (i2rep(i),i=0,nodes-1)
576 write (irest1,*) "ifirst"
577 write (irest1,*) (ifirst(i),i=1,remd_m(1))
579 write (irest1,*) "nupa",il
580 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
581 write (irest1,*) "ndowna",il
582 write (irest1,*) ndowna(0,il),
583 & (ndowna(i,il),i=1,ndowna(0,il))
585 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
586 write (irest1,*) "nset"
587 write (irest1,*) nset
588 write (irest1,*) "mset"
589 write (irest1,*) (mset(i),i=1,nset)
590 write (irest1,*) "i2set"
591 write (irest1,*) (i2set(i),i=0,nodes-1)
592 write (irest1,*) "i_index"
596 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
604 open(irest2,file=rest2name,status='unknown')
605 write(irest2,*) totT,EK,potE,totE,t_bath
607 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
610 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
612 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
613 write (irest2,*) iset
620 c forced synchronization
621 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
622 & .and. .not. mremdsync) then
624 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
626 call mpi_recv(itime_master, 1, MPI_INTEGER,
627 & 0,101,CG_COMM, status, ierr)
628 call mpi_barrier(CG_COMM, ierr)
629 cdeb if (out1file.or.traj1file) then
630 cdeb call mpi_gather(itime,1,mpi_integer,
631 cdeb & icache_all,1,mpi_integer,king,
634 & call mpi_gather(ntwx_cache,1,mpi_integer,
635 & icache_all,1,mpi_integer,king,
638 & write(iout,*) 'REMD synchro at',itime_master,itime
639 if (itime_master.ge.n_timestep .or. ovrtim())
641 ctime call flush(iout)
646 if ((mod(itime,nstex).eq.0.and.me.eq.king
647 & .or.end_of_run.and.me.eq.king )
648 & .and. .not. mremdsync ) then
652 call mpi_isend(itime,1,MPI_INTEGER,i,101,
653 & CG_COMM, ireqi(i), ierr)
654 cd write(iout,*) 'REMD synchro with',i
657 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
658 call mpi_barrier(CG_COMM, ierr)
660 write(iout,*) 'REMD synchro at',itime,'time=',time01-time01_
661 if (out1file.or.traj1file) then
662 cdeb call mpi_gather(itime,1,mpi_integer,
663 cdeb & itime_all,1,mpi_integer,king,
665 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
666 cdeb & (itime_all(i),i=1,nodes)
668 cdeb imin_itime=itime_all(1)
670 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
672 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
673 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
674 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
675 call mpi_gather(ntwx_cache,1,mpi_integer,
676 & icache_all,1,mpi_integer,king,
678 c write(iout,'(a19,8000i8)') ' ntwx_cache',
679 c & (icache_all(i),i=1,nodes)
680 ii_write=icache_all(1)
682 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
684 c write(iout,*) "MIN ii_write=",ii_write
687 ctime call flush(iout)
689 if(mremdsync .and. mod(itime,nstex).eq.0) then
691 if (me.eq.king .or. .not. out1file)
692 & write(iout,*) 'REMD synchro at',itime
695 call mpi_gather(ntwx_cache,1,mpi_integer,
696 & icache_all,1,mpi_integer,king,
699 write(iout,'(a19,8000i8)') ' ntwx_cache',
700 & (icache_all(i),i=1,nodes)
701 ii_write=icache_all(1)
703 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
705 write(iout,*) "MIN ii_write=",ii_write
708 ctest call flush(iout)
711 c Update the time safety limiy
712 if (time001-time00.gt.safety) then
713 safety=time001-time00+600
714 if (me.eq.king .or. .not. out1file)
715 & write (iout,*) "****** SAFETY increased to",safety," s"
717 if (ovrtim()) end_of_run=.true.
719 if(synflag.and..not.end_of_run) then
723 c write(iout,*) 'REMD before',me,t_bath
725 c call mpi_gather(t_bath,1,mpi_double_precision,
726 c & remd_t_bath,1,mpi_double_precision,king,
728 potEcomp(n_ene+1)=t_bath
730 if (usampl.or.homol_nset.gt.1) then
731 potEcomp(n_ene+2)=iset
732 c write(iout,*) potEcomp(n_ene+1),potEcomp(n_ene+2),iset,nset
733 if (iset.lt.nset) then
736 if (homol_nset.gt.1) then
737 call e_modeller(potEcomp(n_ene+3))
738 c write(iout,*) "iset+1",potEcomp(n_ene+3)
741 potEcomp(n_ene+3)=Uconst
748 if (homol_nset.gt.1) then
749 call e_modeller(potEcomp(n_ene+4))
750 c write(iout,*) "iset-1",potEcomp(n_ene+4)
753 potEcomp(n_ene+4)=Uconst
758 if(hremd.gt.0) potEcomp(n_ene+2)=iset
759 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
760 & remd_ene(0,1),n_ene+5,mpi_double_precision,king,
763 call mpi_gather(elow,1,mpi_double_precision,
764 & elowi,1,mpi_double_precision,king,
766 call mpi_gather(ehigh,1,mpi_double_precision,
767 & ehighi,1,mpi_double_precision,king,
772 if (me.eq.king .or. .not. out1file) then
773 write(iout,*) 'REMD gather times=',time03-time01
777 if (restart1file) call write1rst(i_index)
780 if (me.eq.king .or. .not. out1file) then
781 write(iout,*) 'REMD writing rst time=',time04-time03
784 if (traj1file) call write1traj
786 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
787 cdeb & icache_all,1,mpi_integer,king,
789 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
790 cdeb & (icache_all(i),i=1,nodes)
795 if (me.eq.king .or. .not. out1file) then
796 write(iout,*) 'REMD writing traj time=',time05-time04
797 ctime call flush(iout)
802 if(homol_nset.gt.1) write(iout,*)
803 & 'energy_c temperature iset energy_c(iset+1) energy_c(iset-1)'
805 remd_t_bath(i)=remd_ene(n_ene+1,i)
806 iremd_iset(i)=remd_ene(n_ene+2,i)
808 & write(iout,'(i4,f10.3,f6.0,i3,2f10.3)')
809 & i,remd_ene(i_econstr,i),
810 & remd_ene(n_ene+1,i),iremd_iset(i),
811 & remd_ene(n_ene+3,i),remd_ene(n_ene+4,i)
815 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
817 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
821 write(iout,*) 'REMD exchange temp,ene'
823 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
824 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
828 c-------------------------------------
829 IF(.not.usampl.and.hremd.eq.0.and.homol_nset.le.1) THEN
831 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
833 ctime call flush(iout)
834 write (iout,*) "remd_m(1)",remd_m(1)
837 i=ifirst(iran_num(1,remd_m(1)))
841 ctime call flush(iout)
846 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
848 if(i.gt.0.and.nupa(0,i).gt.0) then
850 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
852 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
854 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
856 c do while (iex.eq.i)
857 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
858 iex=nupa(iran_num(1,int(nupa(0,i))),i)
860 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
862 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
864 c Swap temperatures between conformations i and iex with recalculating the free energies
865 c following temperature changes.
866 ene_iex_iex=remd_ene(0,iex)
867 ene_i_i=remd_ene(0,i)
868 c write (iout,*) "i",i," ene_i_i",ene_i_i,
869 c & " iex",iex," ene_iex_iex",ene_iex_iex
870 c write (iout,*) "rescaling weights with temperature",
873 call rescale_weights(remd_t_bath(i))
875 c write (iout,*) "0,iex",remd_t_bath(i)
876 c call enerprint(remd_ene(0,iex))
878 call sum_energy(remd_ene(0,iex),.false.)
879 ene_iex_i=remd_ene(0,iex)
880 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
882 c write (iout,*) "0,i",remd_t_bath(i)
883 c call enerprint(remd_ene(0,i))
885 call sum_energy(remd_ene(0,i),.false.)
886 c write (iout,*) "ene_i_i",remd_ene(0,i)
888 c write (iout,*) "rescaling weights with temperature",
890 if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
891 write (iout,*) "ERROR: inconsistent energies:",i,
892 & ene_i_i,remd_ene(0,i)
894 call rescale_weights(remd_t_bath(iex))
896 c write (iout,*) "0,i",remd_t_bath(iex)
897 c call enerprint(remd_ene(0,i))
899 call sum_energy(remd_ene(0,i),.false.)
900 c write (iout,*) "ene_i_iex",remd_ene(0,i)
902 ene_i_iex=remd_ene(0,i)
904 c write (iout,*) "0,iex",remd_t_bath(iex)
905 c call enerprint(remd_ene(0,iex))
907 call sum_energy(remd_ene(0,iex),.false.)
908 if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
909 write (iout,*) "ERROR: inconsistent energies:",iex,
910 & ene_iex_iex,remd_ene(0,iex)
912 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
913 c write (iout,*) "i",i," iex",iex
914 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
915 c & " ene_i_iex",ene_i_iex,
916 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
918 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
919 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
921 c write(iout,*) 'delta',delta
922 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
923 c & (remd_ene(i)-remd_ene(iex))/Rb/
924 c & (remd_t_bath(i)*remd_t_bath(iex))
926 if (delta .gt. 50.0d0) then
932 else if (delta.lt.-50.0d0) then
941 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
942 xxx=ran_number(0.0d0,1.0d0)
943 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
945 if (delta .gt. xxx) then
947 remd_t_bath(i)=remd_t_bath(iex)
949 remd_ene(0,i)=ene_i_iex
950 remd_ene(0,iex)=ene_iex_i
956 ehighi(i)=ehighi(iex)
963 nupa(k,i)=nupa(k,iex)
966 ndowna(k,i)=ndowna(k,iex)
970 if (ifirst(il).eq.i) ifirst(il)=iex
972 if (nupa(k,il).eq.i) then
974 elseif (nupa(k,il).eq.iex) then
979 if (ndowna(k,il).eq.i) then
981 elseif (ndowna(k,il).eq.iex) then
987 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
989 i2rep(i-1)=i2rep(iex-1)
992 c write(iout,*) 'exchange',i,iex
993 c write (iout,'(a8,100i4)') "@ ifirst",
994 c & (ifirst(k),k=1,remd_m(1))
996 c write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
997 c & (nupa(k,il),k=1,nupa(0,il))
998 c write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
999 c & (ndowna(k,il),k=1,ndowna(0,il))
1004 remd_ene(0,iex)=ene_iex_iex
1005 remd_ene(0,i)=ene_i_i
1011 cd write (iout,*) "exchange completed"
1013 ELSEIF (usampl.or.homol_nset.gt.1) THEN
1015 c write(iout,*) "########",ii
1017 i_temp=iran_num(1,nrep)
1018 i_mult=iran_num(1,remd_m(i_temp))
1019 i_iset=iran_num(1,nset)
1020 i_mset=iran_num(1,mset(i_iset))
1021 i=i_index(i_temp,i_mult,i_iset,i_mset)
1023 c write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1025 if(t_exchange_only)then
1030 c write(iout,*) "i_dir=",i_dir
1032 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1035 i_mult1=iran_num(1,remd_m(i_temp1))
1037 i_mset1=iran_num(1,mset(i_iset1))
1038 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1040 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1043 i_mult1=iran_num(1,remd_m(i_temp1))
1045 i_mset1=iran_num(1,mset(i_iset1))
1046 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1048 econstr_temp_i=remd_ene(i_econstr,i)
1049 econstr_temp_iex=remd_ene(i_econstr,iex)
1050 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1051 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1053 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1056 i_mult1=iran_num(1,remd_m(i_temp1))
1058 i_mset1=iran_num(1,mset(i_iset1))
1059 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1060 econstr_temp_i=remd_ene(i_econstr,i)
1061 econstr_temp_iex=remd_ene(i_econstr,iex)
1062 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1063 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1069 c write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1070 ctime call flush(iout)
1072 c Swap temperatures between conformations i and iex with recalculating the free energies
1073 c following temperature changes.
1074 ene_iex_iex=remd_ene(0,iex)
1075 ene_i_i=remd_ene(0,i)
1076 co write (iout,*) "rescaling weights with temperature",
1078 call rescale_weights(remd_t_bath(i))
1080 call sum_energy(remd_ene(0,iex),.false.)
1081 ene_iex_i=remd_ene(0,iex)
1083 c ERROR only makes sense for dir =1
1084 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
1085 c call sum_energy(remd_ene(0,i),.false.)
1086 c write (iout,*) "ene_i_i",remd_ene(0,i)
1087 c write (iout,*) "rescaling weights with temperature",
1088 c & remd_t_bath(iex)
1089 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1090 c write (iout,*) "ERROR: inconsistent energies i:",i,
1091 c & ene_i_i,remd_ene(0,i)
1094 call rescale_weights(remd_t_bath(iex))
1095 call sum_energy(remd_ene(0,i),.false.)
1096 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1097 ene_i_iex=remd_ene(0,i)
1099 c ERROR only makes sense for dir =1
1100 c call sum_energy(remd_ene(0,iex),.false.)
1101 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1102 c write (iout,*) "ERROR: inconsistent energies iex:",iex,
1103 c & ene_iex_iex,remd_ene(0,iex)
1105 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1106 c write (iout,*) "i",i," iex",iex
1107 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1108 c & " ene_i_iex",ene_i_iex,
1109 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1111 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1112 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1114 c write(iout,*) 'delta',delta
1115 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1116 c & (remd_ene(i)-remd_ene(iex))/Rb/
1117 c & (remd_t_bath(i)*remd_t_bath(iex))
1118 if (delta .gt. 50.0d0) then
1123 if (i_dir.eq.1.or.i_dir.eq.3)
1124 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1125 if (i_dir.eq.2.or.i_dir.eq.3)
1126 & iremd_tot_usa(int(i2set(i-1)))=
1127 & iremd_tot_usa(int(i2set(i-1)))+1
1128 xxx=ran_number(0.0d0,1.0d0)
1129 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1130 if (delta .gt. xxx) then
1132 remd_t_bath(i)=remd_t_bath(iex)
1133 remd_t_bath(iex)=tmp
1136 iremd_iset(i)=iremd_iset(iex)
1137 iremd_iset(iex)=itmp
1139 remd_ene(0,i)=ene_i_iex
1140 remd_ene(0,iex)=ene_iex_i
1142 if (i_dir.eq.1.or.i_dir.eq.3)
1143 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1146 i2rep(i-1)=i2rep(iex-1)
1149 if (i_dir.eq.2.or.i_dir.eq.3)
1150 & iremd_acc_usa(int(i2set(i-1)))=
1151 & iremd_acc_usa(int(i2set(i-1)))+1
1154 i2set(i-1)=i2set(iex-1)
1157 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1158 i_index(i_temp,i_mult,i_iset,i_mset)=
1159 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1160 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1163 remd_ene(0,iex)=ene_iex_iex
1164 remd_ene(0,i)=ene_i_i
1165 remd_ene(i_econstr,iex)=econstr_temp_iex
1166 remd_ene(i_econstr,i)=econstr_temp_i
1170 cd do il1=1,mset(il)
1173 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1183 ELSEIF (hremd.gt.0) THEN
1185 cd write(iout,*) "########",ii
1187 i_temp=iran_num(1,nrep)
1188 i_mult=iran_num(1,remd_m(i_temp))
1189 i_iset=iran_num(1,nset)
1191 i=i_index(i_temp,i_mult,i_iset,i_mset)
1193 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1195 if(t_exchange_only)then
1201 cd write(iout,*) "i_dir=",i_dir
1203 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1206 i_mult1=iran_num(1,remd_m(i_temp1))
1209 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1211 elseif(i_dir.eq.2)then
1214 i_mult1=iran_num(1,remd_m(i_temp1))
1215 i_iset1=iran_num(1,hremd)
1216 do while(i_iset1.eq.i_iset)
1217 i_iset1=iran_num(1,hremd)
1220 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1222 elseif(remd_m(i_temp+1).gt.0)then
1225 i_mult1=iran_num(1,remd_m(i_temp1))
1226 i_iset1=iran_num(1,hremd)
1227 do while(i_iset1.eq.i_iset)
1228 i_iset1=iran_num(1,hremd)
1231 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1237 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1238 ctime call flush(iout)
1240 c Swap temperatures between conformations i and iex with recalculating the free energies
1241 c following temperature changes.
1242 ene_iex_iex=remd_ene(0,iex)
1243 ene_i_i=remd_ene(0,i)
1245 call set_hweights(i_iset)
1246 call rescale_weights(remd_t_bath(i))
1247 call sum_energy(remd_ene(0,iex),.false.)
1248 ene_iex_i=remd_ene(0,iex)
1250 call set_hweights(i_iset1)
1251 call rescale_weights(remd_t_bath(iex))
1252 call sum_energy(remd_ene(0,i),.false.)
1253 ene_i_iex=remd_ene(0,i)
1255 cd write(iout,*) ene_iex_iex,ene_i_i,ene_iex_i,ene_i_iex
1257 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1258 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1261 if (delta .gt. 50.0d0) then
1267 if (i_dir.eq.1.or.i_dir.eq.3)
1268 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1269 if (i_dir.eq.2.or.i_dir.eq.3)
1270 & iremd_tot_usa(int(i2set(i-1)))=
1271 & iremd_tot_usa(int(i2set(i-1)))+1
1272 xxx=ran_number(0.0d0,1.0d0)
1273 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1274 if (delta .gt. xxx) then
1276 cd write (iout,*) "exchange"
1278 remd_t_bath(i)=remd_t_bath(iex)
1279 remd_t_bath(iex)=tmp
1282 iremd_iset(i)=iremd_iset(iex)
1283 iremd_iset(iex)=itmp
1285 remd_ene(0,i)=ene_i_iex
1286 remd_ene(0,iex)=ene_iex_i
1288 if (i_dir.eq.1.or.i_dir.eq.3)
1289 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1292 i2rep(i-1)=i2rep(iex-1)
1295 if (i_dir.eq.2.or.i_dir.eq.3)
1296 & iremd_acc_usa(int(i2set(i-1)))=
1297 & iremd_acc_usa(int(i2set(i-1)))+1
1300 i2set(i-1)=i2set(iex-1)
1303 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1304 i_index(i_temp,i_mult,i_iset,i_mset)=
1305 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1306 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1309 cd do il1=1,mset(il)
1312 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1319 remd_ene(0,iex)=ene_iex_iex
1320 remd_ene(0,i)=ene_i_i
1331 c-------------------------------------
1332 write (iout,*) "NREP",nrep
1334 if(iremd_tot(i).ne.0)
1335 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1336 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1339 if(usampl.or.homol_nset.gt.1) then
1341 if(iremd_tot_usa(i).ne.0)
1342 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1343 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1349 if(iremd_tot_usa(i).ne.0)
1350 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_hremd',i,
1351 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1356 ctime call flush(iout)
1358 cd write (iout,'(a6,100i4)') "ifirst",
1359 cd & (ifirst(i),i=1,remd_m(1))
1361 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1362 cd & (nupa(i,il),i=1,nupa(0,il))
1363 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1364 cd & (ndowna(i,il),i=1,ndowna(0,il))
1369 cd write (iout,*) "Before scatter"
1371 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1372 & t_bath,1,mpi_double_precision,king,
1374 cd write (iout,*) "After scatter"
1376 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1)
1377 & call mpi_scatter(iremd_iset,1,mpi_integer,
1378 & iset,1,mpi_integer,king,
1382 if (me.eq.king .or. .not. out1file) then
1383 write(iout,*) 'REMD scatter time=',time07-time06
1387 call mpi_scatter(elowi,1,mpi_double_precision,
1388 & elow,1,mpi_double_precision,king,
1390 call mpi_scatter(ehighi,1,mpi_double_precision,
1391 & ehigh,1,mpi_double_precision,king,
1395 if(hremd.gt.0) call set_hweights(iset)
1396 call rescale_weights(t_bath)
1397 co write (iout,*) "Processor",me,
1398 co & " rescaling weights with temperature",t_bath
1400 stdfp=dsqrt(2*Rb*t_bath/d_time)
1402 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1406 stdforcp(i)=stdforcp(i)*sqrt(t_bath/t_bath_old)
1409 stdforcsc(i)=stdforcsc(i)*sqrt(t_bath/t_bath_old)
1412 cde write(iout,*) 'REMD after',me,t_bath
1414 if (me.eq.king .or. .not. out1file) then
1415 write(iout,*) 'REMD exchange time=',time08-time02
1416 ctime call flush(iout)
1421 if (restart1file) then
1422 if (me.eq.king .or. .not. out1file)
1423 & write(iout,*) 'writing restart at the end of run'
1424 call write1rst(i_index)
1427 if (traj1file) call write1traj
1429 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1430 cdeb & icache_all,1,mpi_integer,king,
1431 cdeb & CG_COMM,ierr)
1432 cdeb write(iout,'(a40,8000i8)')
1433 cdeb & ' ntwx_cache after traj1file at the end',
1434 cdeb & (icache_all(i),i=1,nodes)
1439 t_MD=MPI_Wtime()-tt0
1443 if (me.eq.king .or. .not. out1file) then
1444 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1446 & 'MD calculations setup:',t_MDsetup,
1447 & 'Energy & gradient evaluation:',t_enegrad,
1448 & 'Stochastic MD setup:',t_langsetup,
1449 & 'Stochastic MD step setup:',t_sdsetup,
1451 write (iout,'(/28(1h=),a25,27(1h=))')
1452 & ' End of MD calculation '
1453 if(hmc.gt.0) write (iout,*) 'HMC acceptance ratio',
1454 & n_timestep*1.0d0/hmc/hmc_acc
1459 c-----------------------------------------------------------------------
1460 subroutine write1rst(i_index)
1461 implicit real*8 (a-h,o-z)
1462 include 'DIMENSIONS'
1465 include 'COMMON.IOUNITS'
1466 include 'COMMON.REMD'
1467 include 'COMMON.SETUP'
1468 include 'COMMON.CHAIN'
1469 include 'COMMON.SBRIDGE'
1470 include 'COMMON.INTERACT'
1471 include 'COMMON.CONTROL'
1473 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1474 & d_restart2(3,2*maxres*maxprocs)
1478 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1479 common /przechowalnia/ d_restart1,d_restart2
1484 t5_restart1(4)=t_bath
1485 t5_restart1(5)=Uconst
1487 call mpi_gather(t5_restart1,5,mpi_real,
1488 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1496 call mpi_gather(r_d,3*2*nres,mpi_real,
1497 & d_restart1,3*2*nres,mpi_real,king,
1506 call mpi_gather(r_d,3*2*nres,mpi_real,
1507 & d_restart2,3*2*nres,mpi_real,king,
1512 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1514 call xdrfint_(ixdrf, i2rep(i), iret)
1517 call xdrfint_(ixdrf, ifirst(i), iret)
1521 call xdrfint_(ixdrf, nupa(i,il), iret)
1525 call xdrfint_(ixdrf, ndowna(i,il), iret)
1531 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1538 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1545 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1550 if(usampl.or.homol_nset.gt.1) then
1551 call xdrfint_(ixdrf, nset, iret)
1553 call xdrfint_(ixdrf,mset(i), iret)
1556 call xdrfint_(ixdrf,i2set(i), iret)
1562 itmp=i_index(i,j,il,il1)
1563 call xdrfint_(ixdrf,itmp, iret)
1570 call xdrfclose_(ixdrf, iret)
1572 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1574 call xdrfint(ixdrf, i2rep(i), iret)
1577 call xdrfint(ixdrf, ifirst(i), iret)
1581 call xdrfint(ixdrf, nupa(i,il), iret)
1585 call xdrfint(ixdrf, ndowna(i,il), iret)
1591 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1598 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1605 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1611 if(usampl.or.homol_nset.gt.1) then
1612 call xdrfint(ixdrf, nset, iret)
1614 call xdrfint(ixdrf,mset(i), iret)
1617 call xdrfint(ixdrf,i2set(i), iret)
1623 itmp=i_index(i,j,il,il1)
1624 call xdrfint(ixdrf,itmp, iret)
1631 call xdrfclose(ixdrf, iret)
1638 subroutine write1traj
1639 implicit real*8 (a-h,o-z)
1640 include 'DIMENSIONS'
1643 include 'COMMON.IOUNITS'
1644 include 'COMMON.REMD'
1645 include 'COMMON.SETUP'
1646 include 'COMMON.CHAIN'
1647 include 'COMMON.SBRIDGE'
1648 include 'COMMON.INTERACT'
1652 real xcoord(3,maxres2+2),prec
1653 real r_qfrag(50),r_qpair(100)
1654 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1655 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1656 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1657 & p_uscdiff(100*maxprocs)
1658 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1659 common /przechowalnia/ p_c
1661 call mpi_bcast(ii_write,1,mpi_integer,
1662 & king,CG_COMM,ierr)
1665 print *,'traj1file',me,ii_write,ntwx_cache
1669 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1671 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1674 t5_restart1(1)=totT_cache(ii)
1675 t5_restart1(2)=EK_cache(ii)
1676 t5_restart1(3)=potE_cache(ii)
1677 t5_restart1(4)=t_bath_cache(ii)
1678 t5_restart1(5)=Uconst_cache(ii)
1679 call mpi_gather(t5_restart1,5,mpi_real,
1680 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1682 call mpi_gather(iset_cache(ii),1,mpi_integer,
1683 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1686 r_qfrag(i)=qfrag_cache(i,ii)
1689 r_qpair(i)=qpair_cache(i,ii)
1692 r_utheta(i)=utheta_cache(i,ii)
1693 r_ugamma(i)=ugamma_cache(i,ii)
1694 r_uscdiff(i)=uscdiff_cache(i,ii)
1697 call mpi_gather(r_qfrag,nfrag,mpi_real,
1698 & p_qfrag,nfrag,mpi_real,king,
1700 call mpi_gather(r_qpair,npair,mpi_real,
1701 & p_qpair,npair,mpi_real,king,
1703 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1704 & p_utheta,nfrag_back,mpi_real,king,
1706 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1707 & p_ugamma,nfrag_back,mpi_real,king,
1709 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1710 & p_uscdiff,nfrag_back,mpi_real,king,
1714 write (iout,*) "p_qfrag"
1716 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1718 write (iout,*) "p_qpair"
1720 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1722 ctime call flush(iout)
1726 r_c(j,i)=c_cache(j,i,ii)
1730 call mpi_gather(r_c,3*2*nres,mpi_real,
1731 & p_c,3*2*nres,mpi_real,king,
1737 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1738 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1739 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1740 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1741 call xdrfint_(ixdrf, nss, iret)
1744 call xdrfint_(ixdrf, idssb(j)+nres, iret)
1745 call xdrfint_(ixdrf, jdssb(j)+nres, iret)
1747 call xdrfint_(ixdrf, ihpb(j), iret)
1748 call xdrfint_(ixdrf, jhpb(j), iret)
1751 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1752 call xdrfint_(ixdrf, iset_restart1(il), iret)
1754 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1757 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1760 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1761 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1762 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1767 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1772 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1776 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1780 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1781 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1782 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1783 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1784 call xdrfint(ixdrf, nss, iret)
1787 call xdrfint(ixdrf, idssb(j)+nres, iret)
1788 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1790 call xdrfint(ixdrf, ihpb(j), iret)
1791 call xdrfint(ixdrf, jhpb(j), iret)
1794 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1795 call xdrfint(ixdrf, iset_restart1(il), iret)
1797 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1800 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1803 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1804 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1805 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1810 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1815 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1819 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1825 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1827 if(me.eq.king) call xdrfclose(ixdrf, iret)
1829 do i=1,ntwx_cache-ii_write
1831 totT_cache(i)=totT_cache(ii_write+i)
1832 EK_cache(i)=EK_cache(ii_write+i)
1833 potE_cache(i)=potE_cache(ii_write+i)
1834 t_bath_cache(i)=t_bath_cache(ii_write+i)
1835 Uconst_cache(i)=Uconst_cache(ii_write+i)
1836 iset_cache(i)=iset_cache(ii_write+i)
1839 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1842 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1845 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1846 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1847 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1852 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1856 ntwx_cache=ntwx_cache-ii_write
1861 subroutine read1restart(i_index)
1862 implicit real*8 (a-h,o-z)
1863 include 'DIMENSIONS'
1866 include 'COMMON.IOUNITS'
1867 include 'COMMON.REMD'
1868 include 'COMMON.SETUP'
1869 include 'COMMON.CHAIN'
1870 include 'COMMON.SBRIDGE'
1871 include 'COMMON.INTERACT'
1872 include 'COMMON.CONTROL'
1873 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1876 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1877 common /przechowalnia/ d_restart1
1878 write (*,*) "Processor",me," called read1restart"
1881 open(irest2,file=mremd_rst_name,status='unknown')
1882 read(irest2,*,err=334) i
1883 write(iout,*) "Reading old rst in ASCI format"
1885 call read1restart_old
1889 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1892 call xdrfint_(ixdrf, i2rep(i), iret)
1895 call xdrfint_(ixdrf, ifirst(i), iret)
1898 call xdrfint_(ixdrf, nupa(0,il), iret)
1900 call xdrfint_(ixdrf, nupa(i,il), iret)
1903 call xdrfint_(ixdrf, ndowna(0,il), iret)
1905 call xdrfint_(ixdrf, ndowna(i,il), iret)
1910 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1914 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1917 call xdrfint(ixdrf, i2rep(i), iret)
1920 call xdrfint(ixdrf, ifirst(i), iret)
1923 call xdrfint(ixdrf, nupa(0,il), iret)
1925 call xdrfint(ixdrf, nupa(i,il), iret)
1928 call xdrfint(ixdrf, ndowna(0,il), iret)
1930 call xdrfint(ixdrf, ndowna(i,il), iret)
1935 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1940 call mpi_scatter(t_restart1,5,mpi_real,
1941 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1945 t_bath=t5_restart1(4)
1950 c read(irest2,'(3e15.5)')
1951 c & (d_restart1(j,i+2*nres*il),j=1,3)
1954 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1956 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1962 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1963 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1973 c read(irest2,'(3e15.5)')
1974 c & (d_restart1(j,i+2*nres*il),j=1,3)
1977 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1979 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1985 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1986 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1994 if(usampl.or.homol_nset.gt.1) then
1997 call xdrfint_(ixdrf, nset, iret)
1999 call xdrfint_(ixdrf,mset(i), iret)
2002 call xdrfint_(ixdrf,i2set(i), iret)
2008 call xdrfint_(ixdrf,itmp, iret)
2009 i_index(i,j,il,il1)=itmp
2017 call xdrfint(ixdrf, nset, iret)
2019 call xdrfint(ixdrf,mset(i), iret)
2022 call xdrfint(ixdrf,i2set(i), iret)
2028 call xdrfint(ixdrf,itmp, iret)
2029 i_index(i,j,il,il1)=itmp
2037 c Corrected AL 8/19/2014: each processor needs whole iset array not only its
2039 Corrected AL 8/19/2014: each processor needs whole iset array not only its
2040 c>>>>>>> 66693b0684c228404e7aadcffe6d2a8c9f489063
2042 c call mpi_scatter(i2set,1,mpi_integer,
2043 c & iset,1,mpi_integer,king,
2045 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
2052 if(me.eq.king) close(irest2)
2056 subroutine read1restart_old
2057 implicit real*8 (a-h,o-z)
2058 include 'DIMENSIONS'
2061 include 'COMMON.IOUNITS'
2062 include 'COMMON.REMD'
2063 include 'COMMON.SETUP'
2064 include 'COMMON.CHAIN'
2065 include 'COMMON.SBRIDGE'
2066 include 'COMMON.INTERACT'
2067 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
2069 common /przechowalnia/ d_restart1
2071 open(irest2,file=mremd_rst_name,status='unknown')
2072 read (irest2,*) (i2rep(i),i=0,nodes-1)
2073 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2075 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2076 read (irest2,*) ndowna(0,il),
2077 & (ndowna(i,il),i=1,ndowna(0,il))
2080 read(irest2,*) (t_restart1(j,il),j=1,4)
2083 call mpi_scatter(t_restart1,5,mpi_real,
2084 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2088 t_bath=t5_restart1(4)
2093 read(irest2,'(3e15.5)')
2094 & (d_restart1(j,i+2*nres*il),j=1,3)
2098 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2099 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2109 read(irest2,'(3e15.5)')
2110 & (d_restart1(j,i+2*nres*il),j=1,3)
2114 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2115 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2121 if(me.eq.king) close(irest2)
2124 c-------------------------------------------------------------------
2125 subroutine set_hweights(iiset)
2126 implicit real*8 (a-h,o-z)
2128 include 'DIMENSIONS'
2129 include 'COMMON.FFIELD'
2130 include 'COMMON.REMD'
2133 weights(i)=hweights(iiset,i)