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
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)
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)
807 write(iout,'(i4,f10.3,f6.0,i3,2f10.3)')
808 & i,remd_ene(i_econstr,i),
809 & remd_ene(n_ene+1,i),iremd_iset(i),
810 & remd_ene(n_ene+3,i),remd_ene(n_ene+4,i)
814 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
816 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
820 write(iout,*) 'REMD exchange temp,ene'
822 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
823 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
827 c-------------------------------------
828 IF(.not.usampl.and.hremd.eq.0.and.homol_nset.le.1) THEN
830 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
832 ctime call flush(iout)
833 write (iout,*) "remd_m(1)",remd_m(1)
836 i=ifirst(iran_num(1,remd_m(1)))
840 ctime call flush(iout)
845 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
847 if(i.gt.0.and.nupa(0,i).gt.0) then
849 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
851 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
853 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
855 c do while (iex.eq.i)
856 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
857 iex=nupa(iran_num(1,int(nupa(0,i))),i)
859 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
861 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
863 c Swap temperatures between conformations i and iex with recalculating the free energies
864 c following temperature changes.
865 ene_iex_iex=remd_ene(0,iex)
866 ene_i_i=remd_ene(0,i)
867 c write (iout,*) "i",i," ene_i_i",ene_i_i,
868 c & " iex",iex," ene_iex_iex",ene_iex_iex
869 c write (iout,*) "rescaling weights with temperature",
872 call rescale_weights(remd_t_bath(i))
874 c write (iout,*) "0,iex",remd_t_bath(i)
875 c call enerprint(remd_ene(0,iex))
877 call sum_energy(remd_ene(0,iex),.false.)
878 ene_iex_i=remd_ene(0,iex)
879 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
881 c write (iout,*) "0,i",remd_t_bath(i)
882 c call enerprint(remd_ene(0,i))
884 call sum_energy(remd_ene(0,i),.false.)
885 c write (iout,*) "ene_i_i",remd_ene(0,i)
887 c write (iout,*) "rescaling weights with temperature",
889 if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
890 write (iout,*) "ERROR: inconsistent energies:",i,
891 & ene_i_i,remd_ene(0,i)
893 call rescale_weights(remd_t_bath(iex))
895 c write (iout,*) "0,i",remd_t_bath(iex)
896 c call enerprint(remd_ene(0,i))
898 call sum_energy(remd_ene(0,i),.false.)
899 c write (iout,*) "ene_i_iex",remd_ene(0,i)
901 ene_i_iex=remd_ene(0,i)
903 c write (iout,*) "0,iex",remd_t_bath(iex)
904 c call enerprint(remd_ene(0,iex))
906 call sum_energy(remd_ene(0,iex),.false.)
907 if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
908 write (iout,*) "ERROR: inconsistent energies:",iex,
909 & ene_iex_iex,remd_ene(0,iex)
911 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
912 c write (iout,*) "i",i," iex",iex
913 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
914 c & " ene_i_iex",ene_i_iex,
915 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
917 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
918 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
920 c write(iout,*) 'delta',delta
921 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
922 c & (remd_ene(i)-remd_ene(iex))/Rb/
923 c & (remd_t_bath(i)*remd_t_bath(iex))
925 if (delta .gt. 50.0d0) then
931 else if (delta.lt.-50.0d0) then
940 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
941 xxx=ran_number(0.0d0,1.0d0)
942 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
944 if (delta .gt. xxx) then
946 remd_t_bath(i)=remd_t_bath(iex)
948 remd_ene(0,i)=ene_i_iex
949 remd_ene(0,iex)=ene_iex_i
955 ehighi(i)=ehighi(iex)
962 nupa(k,i)=nupa(k,iex)
965 ndowna(k,i)=ndowna(k,iex)
969 if (ifirst(il).eq.i) ifirst(il)=iex
971 if (nupa(k,il).eq.i) then
973 elseif (nupa(k,il).eq.iex) then
978 if (ndowna(k,il).eq.i) then
980 elseif (ndowna(k,il).eq.iex) then
986 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
988 i2rep(i-1)=i2rep(iex-1)
991 c write(iout,*) 'exchange',i,iex
992 c write (iout,'(a8,100i4)') "@ ifirst",
993 c & (ifirst(k),k=1,remd_m(1))
995 c write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
996 c & (nupa(k,il),k=1,nupa(0,il))
997 c write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
998 c & (ndowna(k,il),k=1,ndowna(0,il))
1003 remd_ene(0,iex)=ene_iex_iex
1004 remd_ene(0,i)=ene_i_i
1010 cd write (iout,*) "exchange completed"
1012 ELSEIF (usampl.or.homol_nset.gt.1) THEN
1014 c write(iout,*) "########",ii
1016 i_temp=iran_num(1,nrep)
1017 i_mult=iran_num(1,remd_m(i_temp))
1018 i_iset=iran_num(1,nset)
1019 i_mset=iran_num(1,mset(i_iset))
1020 i=i_index(i_temp,i_mult,i_iset,i_mset)
1022 c write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1024 if(t_exchange_only)then
1029 c write(iout,*) "i_dir=",i_dir
1031 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1034 i_mult1=iran_num(1,remd_m(i_temp1))
1036 i_mset1=iran_num(1,mset(i_iset1))
1037 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1039 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1042 i_mult1=iran_num(1,remd_m(i_temp1))
1044 i_mset1=iran_num(1,mset(i_iset1))
1045 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1047 econstr_temp_i=remd_ene(i_econstr,i)
1048 econstr_temp_iex=remd_ene(i_econstr,iex)
1049 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1050 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1052 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1055 i_mult1=iran_num(1,remd_m(i_temp1))
1057 i_mset1=iran_num(1,mset(i_iset1))
1058 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1059 econstr_temp_i=remd_ene(i_econstr,i)
1060 econstr_temp_iex=remd_ene(i_econstr,iex)
1061 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1062 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1068 c write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1069 ctime call flush(iout)
1071 c Swap temperatures between conformations i and iex with recalculating the free energies
1072 c following temperature changes.
1073 ene_iex_iex=remd_ene(0,iex)
1074 ene_i_i=remd_ene(0,i)
1075 co write (iout,*) "rescaling weights with temperature",
1077 call rescale_weights(remd_t_bath(i))
1079 call sum_energy(remd_ene(0,iex),.false.)
1080 ene_iex_i=remd_ene(0,iex)
1082 c ERROR only makes sense for dir =1
1083 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
1084 c call sum_energy(remd_ene(0,i),.false.)
1085 c write (iout,*) "ene_i_i",remd_ene(0,i)
1086 c write (iout,*) "rescaling weights with temperature",
1087 c & remd_t_bath(iex)
1088 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1089 c write (iout,*) "ERROR: inconsistent energies i:",i,
1090 c & ene_i_i,remd_ene(0,i)
1093 call rescale_weights(remd_t_bath(iex))
1094 call sum_energy(remd_ene(0,i),.false.)
1095 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1096 ene_i_iex=remd_ene(0,i)
1098 c ERROR only makes sense for dir =1
1099 c call sum_energy(remd_ene(0,iex),.false.)
1100 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1101 c write (iout,*) "ERROR: inconsistent energies iex:",iex,
1102 c & ene_iex_iex,remd_ene(0,iex)
1104 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1105 c write (iout,*) "i",i," iex",iex
1106 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1107 c & " ene_i_iex",ene_i_iex,
1108 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1110 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1111 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1113 c write(iout,*) 'delta',delta
1114 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1115 c & (remd_ene(i)-remd_ene(iex))/Rb/
1116 c & (remd_t_bath(i)*remd_t_bath(iex))
1117 if (delta .gt. 50.0d0) then
1122 if (i_dir.eq.1.or.i_dir.eq.3)
1123 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1124 if (i_dir.eq.2.or.i_dir.eq.3)
1125 & iremd_tot_usa(int(i2set(i-1)))=
1126 & iremd_tot_usa(int(i2set(i-1)))+1
1127 xxx=ran_number(0.0d0,1.0d0)
1128 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1129 if (delta .gt. xxx) then
1131 remd_t_bath(i)=remd_t_bath(iex)
1132 remd_t_bath(iex)=tmp
1135 iremd_iset(i)=iremd_iset(iex)
1136 iremd_iset(iex)=itmp
1138 remd_ene(0,i)=ene_i_iex
1139 remd_ene(0,iex)=ene_iex_i
1141 if (i_dir.eq.1.or.i_dir.eq.3)
1142 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1145 i2rep(i-1)=i2rep(iex-1)
1148 if (i_dir.eq.2.or.i_dir.eq.3)
1149 & iremd_acc_usa(int(i2set(i-1)))=
1150 & iremd_acc_usa(int(i2set(i-1)))+1
1153 i2set(i-1)=i2set(iex-1)
1156 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1157 i_index(i_temp,i_mult,i_iset,i_mset)=
1158 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1159 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1162 remd_ene(0,iex)=ene_iex_iex
1163 remd_ene(0,i)=ene_i_i
1164 remd_ene(i_econstr,iex)=econstr_temp_iex
1165 remd_ene(i_econstr,i)=econstr_temp_i
1169 cd do il1=1,mset(il)
1172 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1182 ELSEIF (hremd.gt.0) THEN
1184 cd write(iout,*) "########",ii
1186 i_temp=iran_num(1,nrep)
1187 i_mult=iran_num(1,remd_m(i_temp))
1188 i_iset=iran_num(1,nset)
1190 i=i_index(i_temp,i_mult,i_iset,i_mset)
1192 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1194 if(t_exchange_only)then
1200 cd write(iout,*) "i_dir=",i_dir
1202 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1205 i_mult1=iran_num(1,remd_m(i_temp1))
1208 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1210 elseif(i_dir.eq.2)then
1213 i_mult1=iran_num(1,remd_m(i_temp1))
1214 i_iset1=iran_num(1,hremd)
1215 do while(i_iset1.eq.i_iset)
1216 i_iset1=iran_num(1,hremd)
1219 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1221 elseif(remd_m(i_temp+1).gt.0)then
1224 i_mult1=iran_num(1,remd_m(i_temp1))
1225 i_iset1=iran_num(1,hremd)
1226 do while(i_iset1.eq.i_iset)
1227 i_iset1=iran_num(1,hremd)
1230 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1236 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1237 ctime call flush(iout)
1239 c Swap temperatures between conformations i and iex with recalculating the free energies
1240 c following temperature changes.
1241 ene_iex_iex=remd_ene(0,iex)
1242 ene_i_i=remd_ene(0,i)
1244 call set_hweights(i_iset)
1245 call rescale_weights(remd_t_bath(i))
1246 call sum_energy(remd_ene(0,iex),.false.)
1247 ene_iex_i=remd_ene(0,iex)
1249 call set_hweights(i_iset1)
1250 call rescale_weights(remd_t_bath(iex))
1251 call sum_energy(remd_ene(0,i),.false.)
1252 ene_i_iex=remd_ene(0,i)
1254 cd write(iout,*) ene_iex_iex,ene_i_i,ene_iex_i,ene_i_iex
1256 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1257 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1260 if (delta .gt. 50.0d0) then
1266 if (i_dir.eq.1.or.i_dir.eq.3)
1267 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1268 if (i_dir.eq.2.or.i_dir.eq.3)
1269 & iremd_tot_usa(int(i2set(i-1)))=
1270 & iremd_tot_usa(int(i2set(i-1)))+1
1271 xxx=ran_number(0.0d0,1.0d0)
1272 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1273 if (delta .gt. xxx) then
1275 cd write (iout,*) "exchange"
1277 remd_t_bath(i)=remd_t_bath(iex)
1278 remd_t_bath(iex)=tmp
1281 iremd_iset(i)=iremd_iset(iex)
1282 iremd_iset(iex)=itmp
1284 remd_ene(0,i)=ene_i_iex
1285 remd_ene(0,iex)=ene_iex_i
1287 if (i_dir.eq.1.or.i_dir.eq.3)
1288 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1291 i2rep(i-1)=i2rep(iex-1)
1294 if (i_dir.eq.2.or.i_dir.eq.3)
1295 & iremd_acc_usa(int(i2set(i-1)))=
1296 & iremd_acc_usa(int(i2set(i-1)))+1
1299 i2set(i-1)=i2set(iex-1)
1302 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1303 i_index(i_temp,i_mult,i_iset,i_mset)=
1304 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1305 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1308 cd do il1=1,mset(il)
1311 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1318 remd_ene(0,iex)=ene_iex_iex
1319 remd_ene(0,i)=ene_i_i
1330 c-------------------------------------
1331 write (iout,*) "NREP",nrep
1333 if(iremd_tot(i).ne.0)
1334 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1335 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1338 if(usampl.or.homol_nset.gt.1) then
1340 if(iremd_tot_usa(i).ne.0)
1341 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1342 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1348 if(iremd_tot_usa(i).ne.0)
1349 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_hremd',i,
1350 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1355 ctime call flush(iout)
1357 cd write (iout,'(a6,100i4)') "ifirst",
1358 cd & (ifirst(i),i=1,remd_m(1))
1360 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1361 cd & (nupa(i,il),i=1,nupa(0,il))
1362 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1363 cd & (ndowna(i,il),i=1,ndowna(0,il))
1368 cd write (iout,*) "Before scatter"
1370 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1371 & t_bath,1,mpi_double_precision,king,
1373 cd write (iout,*) "After scatter"
1375 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1)
1376 & call mpi_scatter(iremd_iset,1,mpi_integer,
1377 & iset,1,mpi_integer,king,
1381 if (me.eq.king .or. .not. out1file) then
1382 write(iout,*) 'REMD scatter time=',time07-time06
1386 call mpi_scatter(elowi,1,mpi_double_precision,
1387 & elow,1,mpi_double_precision,king,
1389 call mpi_scatter(ehighi,1,mpi_double_precision,
1390 & ehigh,1,mpi_double_precision,king,
1394 if(hremd.gt.0) call set_hweights(iset)
1395 call rescale_weights(t_bath)
1396 co write (iout,*) "Processor",me,
1397 co & " rescaling weights with temperature",t_bath
1399 stdfp=dsqrt(2*Rb*t_bath/d_time)
1401 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1405 stdforcp(i)=stdforcp(i)*sqrt(t_bath/t_bath_old)
1408 stdforcsc(i)=stdforcsc(i)*sqrt(t_bath/t_bath_old)
1411 cde write(iout,*) 'REMD after',me,t_bath
1413 if (me.eq.king .or. .not. out1file) then
1414 write(iout,*) 'REMD exchange time=',time08-time02
1415 ctime call flush(iout)
1420 if (restart1file) then
1421 if (me.eq.king .or. .not. out1file)
1422 & write(iout,*) 'writing restart at the end of run'
1423 call write1rst(i_index)
1426 if (traj1file) call write1traj
1428 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1429 cdeb & icache_all,1,mpi_integer,king,
1430 cdeb & CG_COMM,ierr)
1431 cdeb write(iout,'(a40,8000i8)')
1432 cdeb & ' ntwx_cache after traj1file at the end',
1433 cdeb & (icache_all(i),i=1,nodes)
1438 t_MD=MPI_Wtime()-tt0
1442 if (me.eq.king .or. .not. out1file) then
1443 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1445 & 'MD calculations setup:',t_MDsetup,
1446 & 'Energy & gradient evaluation:',t_enegrad,
1447 & 'Stochastic MD setup:',t_langsetup,
1448 & 'Stochastic MD step setup:',t_sdsetup,
1450 write (iout,'(/28(1h=),a25,27(1h=))')
1451 & ' End of MD calculation '
1452 if(hmc.gt.0) write (iout,*) 'HMC acceptance ratio',
1453 & n_timestep*1.0d0/hmc/hmc_acc
1458 c-----------------------------------------------------------------------
1459 subroutine write1rst(i_index)
1460 implicit real*8 (a-h,o-z)
1461 include 'DIMENSIONS'
1464 include 'COMMON.IOUNITS'
1465 include 'COMMON.REMD'
1466 include 'COMMON.SETUP'
1467 include 'COMMON.CHAIN'
1468 include 'COMMON.SBRIDGE'
1469 include 'COMMON.INTERACT'
1471 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1472 & d_restart2(3,2*maxres*maxprocs)
1476 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1477 common /przechowalnia/ d_restart1,d_restart2
1482 t5_restart1(4)=t_bath
1483 t5_restart1(5)=Uconst
1485 call mpi_gather(t5_restart1,5,mpi_real,
1486 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1494 call mpi_gather(r_d,3*2*nres,mpi_real,
1495 & d_restart1,3*2*nres,mpi_real,king,
1504 call mpi_gather(r_d,3*2*nres,mpi_real,
1505 & d_restart2,3*2*nres,mpi_real,king,
1510 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1512 call xdrfint_(ixdrf, i2rep(i), iret)
1515 call xdrfint_(ixdrf, ifirst(i), iret)
1519 call xdrfint_(ixdrf, nupa(i,il), iret)
1523 call xdrfint_(ixdrf, ndowna(i,il), iret)
1529 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1536 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1543 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1548 if(usampl.or.homol_nset.gt.1) then
1549 call xdrfint_(ixdrf, nset, iret)
1551 call xdrfint_(ixdrf,mset(i), iret)
1554 call xdrfint_(ixdrf,i2set(i), iret)
1560 itmp=i_index(i,j,il,il1)
1561 call xdrfint_(ixdrf,itmp, iret)
1568 call xdrfclose_(ixdrf, iret)
1570 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1572 call xdrfint(ixdrf, i2rep(i), iret)
1575 call xdrfint(ixdrf, ifirst(i), iret)
1579 call xdrfint(ixdrf, nupa(i,il), iret)
1583 call xdrfint(ixdrf, ndowna(i,il), iret)
1589 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1596 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1603 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1609 if(usampl.or.homol_nset.gt.1) then
1610 call xdrfint(ixdrf, nset, iret)
1612 call xdrfint(ixdrf,mset(i), iret)
1615 call xdrfint(ixdrf,i2set(i), iret)
1621 itmp=i_index(i,j,il,il1)
1622 call xdrfint(ixdrf,itmp, iret)
1629 call xdrfclose(ixdrf, iret)
1636 subroutine write1traj
1637 implicit real*8 (a-h,o-z)
1638 include 'DIMENSIONS'
1641 include 'COMMON.IOUNITS'
1642 include 'COMMON.REMD'
1643 include 'COMMON.SETUP'
1644 include 'COMMON.CHAIN'
1645 include 'COMMON.SBRIDGE'
1646 include 'COMMON.INTERACT'
1650 real xcoord(3,maxres2+2),prec
1651 real r_qfrag(50),r_qpair(100)
1652 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1653 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1654 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1655 & p_uscdiff(100*maxprocs)
1656 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1657 common /przechowalnia/ p_c
1659 call mpi_bcast(ii_write,1,mpi_integer,
1660 & king,CG_COMM,ierr)
1663 print *,'traj1file',me,ii_write,ntwx_cache
1667 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1669 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1672 t5_restart1(1)=totT_cache(ii)
1673 t5_restart1(2)=EK_cache(ii)
1674 t5_restart1(3)=potE_cache(ii)
1675 t5_restart1(4)=t_bath_cache(ii)
1676 t5_restart1(5)=Uconst_cache(ii)
1677 call mpi_gather(t5_restart1,5,mpi_real,
1678 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1680 call mpi_gather(iset_cache(ii),1,mpi_integer,
1681 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1684 r_qfrag(i)=qfrag_cache(i,ii)
1687 r_qpair(i)=qpair_cache(i,ii)
1690 r_utheta(i)=utheta_cache(i,ii)
1691 r_ugamma(i)=ugamma_cache(i,ii)
1692 r_uscdiff(i)=uscdiff_cache(i,ii)
1695 call mpi_gather(r_qfrag,nfrag,mpi_real,
1696 & p_qfrag,nfrag,mpi_real,king,
1698 call mpi_gather(r_qpair,npair,mpi_real,
1699 & p_qpair,npair,mpi_real,king,
1701 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1702 & p_utheta,nfrag_back,mpi_real,king,
1704 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1705 & p_ugamma,nfrag_back,mpi_real,king,
1707 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1708 & p_uscdiff,nfrag_back,mpi_real,king,
1712 write (iout,*) "p_qfrag"
1714 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1716 write (iout,*) "p_qpair"
1718 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1720 ctime call flush(iout)
1724 r_c(j,i)=c_cache(j,i,ii)
1728 call mpi_gather(r_c,3*2*nres,mpi_real,
1729 & p_c,3*2*nres,mpi_real,king,
1735 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1736 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1737 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1738 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1739 call xdrfint_(ixdrf, nss, iret)
1742 call xdrfint_(ixdrf, idssb(j)+nres, iret)
1743 call xdrfint_(ixdrf, jdssb(j)+nres, iret)
1745 call xdrfint_(ixdrf, ihpb(j), iret)
1746 call xdrfint_(ixdrf, jhpb(j), iret)
1749 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1750 call xdrfint_(ixdrf, iset_restart1(il), iret)
1752 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1755 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1758 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1759 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1760 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1765 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1770 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1774 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1778 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1779 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1780 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1781 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1782 call xdrfint(ixdrf, nss, iret)
1785 call xdrfint(ixdrf, idssb(j)+nres, iret)
1786 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1788 call xdrfint(ixdrf, ihpb(j), iret)
1789 call xdrfint(ixdrf, jhpb(j), iret)
1792 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1793 call xdrfint(ixdrf, iset_restart1(il), iret)
1795 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1798 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1801 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1802 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1803 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1808 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1813 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1817 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1823 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1825 if(me.eq.king) call xdrfclose(ixdrf, iret)
1827 do i=1,ntwx_cache-ii_write
1829 totT_cache(i)=totT_cache(ii_write+i)
1830 EK_cache(i)=EK_cache(ii_write+i)
1831 potE_cache(i)=potE_cache(ii_write+i)
1832 t_bath_cache(i)=t_bath_cache(ii_write+i)
1833 Uconst_cache(i)=Uconst_cache(ii_write+i)
1834 iset_cache(i)=iset_cache(ii_write+i)
1837 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1840 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1843 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1844 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1845 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1850 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1854 ntwx_cache=ntwx_cache-ii_write
1859 subroutine read1restart(i_index)
1860 implicit real*8 (a-h,o-z)
1861 include 'DIMENSIONS'
1864 include 'COMMON.IOUNITS'
1865 include 'COMMON.REMD'
1866 include 'COMMON.SETUP'
1867 include 'COMMON.CHAIN'
1868 include 'COMMON.SBRIDGE'
1869 include 'COMMON.INTERACT'
1870 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1873 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1874 common /przechowalnia/ d_restart1
1875 write (*,*) "Processor",me," called read1restart"
1878 open(irest2,file=mremd_rst_name,status='unknown')
1879 read(irest2,*,err=334) i
1880 write(iout,*) "Reading old rst in ASCI format"
1882 call read1restart_old
1886 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1889 call xdrfint_(ixdrf, i2rep(i), iret)
1892 call xdrfint_(ixdrf, ifirst(i), iret)
1895 call xdrfint_(ixdrf, nupa(0,il), iret)
1897 call xdrfint_(ixdrf, nupa(i,il), iret)
1900 call xdrfint_(ixdrf, ndowna(0,il), iret)
1902 call xdrfint_(ixdrf, ndowna(i,il), iret)
1907 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1911 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1914 call xdrfint(ixdrf, i2rep(i), iret)
1917 call xdrfint(ixdrf, ifirst(i), iret)
1920 call xdrfint(ixdrf, nupa(0,il), iret)
1922 call xdrfint(ixdrf, nupa(i,il), iret)
1925 call xdrfint(ixdrf, ndowna(0,il), iret)
1927 call xdrfint(ixdrf, ndowna(i,il), iret)
1932 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1937 call mpi_scatter(t_restart1,5,mpi_real,
1938 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1942 t_bath=t5_restart1(4)
1947 c read(irest2,'(3e15.5)')
1948 c & (d_restart1(j,i+2*nres*il),j=1,3)
1951 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1953 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1959 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1960 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1970 c read(irest2,'(3e15.5)')
1971 c & (d_restart1(j,i+2*nres*il),j=1,3)
1974 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1976 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1982 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1983 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1991 if(usampl.or.homol_nset.gt.1) then
1994 call xdrfint_(ixdrf, nset, iret)
1996 call xdrfint_(ixdrf,mset(i), iret)
1999 call xdrfint_(ixdrf,i2set(i), iret)
2005 call xdrfint_(ixdrf,itmp, iret)
2006 i_index(i,j,il,il1)=itmp
2014 call xdrfint(ixdrf, nset, iret)
2016 call xdrfint(ixdrf,mset(i), iret)
2019 call xdrfint(ixdrf,i2set(i), iret)
2025 call xdrfint(ixdrf,itmp, iret)
2026 i_index(i,j,il,il1)=itmp
2033 c Corrected AL 8/19/2014: each processor needs whole iset array not only its
2035 c call mpi_scatter(i2set,1,mpi_integer,
2036 c & iset,1,mpi_integer,king,
2038 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
2045 if(me.eq.king) close(irest2)
2049 subroutine read1restart_old
2050 implicit real*8 (a-h,o-z)
2051 include 'DIMENSIONS'
2054 include 'COMMON.IOUNITS'
2055 include 'COMMON.REMD'
2056 include 'COMMON.SETUP'
2057 include 'COMMON.CHAIN'
2058 include 'COMMON.SBRIDGE'
2059 include 'COMMON.INTERACT'
2060 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
2062 common /przechowalnia/ d_restart1
2064 open(irest2,file=mremd_rst_name,status='unknown')
2065 read (irest2,*) (i2rep(i),i=0,nodes-1)
2066 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2068 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2069 read (irest2,*) ndowna(0,il),
2070 & (ndowna(i,il),i=1,ndowna(0,il))
2073 read(irest2,*) (t_restart1(j,il),j=1,4)
2076 call mpi_scatter(t_restart1,5,mpi_real,
2077 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2081 t_bath=t5_restart1(4)
2086 read(irest2,'(3e15.5)')
2087 & (d_restart1(j,i+2*nres*il),j=1,3)
2091 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2092 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2102 read(irest2,'(3e15.5)')
2103 & (d_restart1(j,i+2*nres*il),j=1,3)
2107 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2108 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2114 if(me.eq.king) close(irest2)
2117 c-------------------------------------------------------------------
2118 subroutine set_hweights(iiset)
2119 implicit real*8 (a-h,o-z)
2121 include 'DIMENSIONS'
2122 include 'COMMON.FFIELD'
2123 include 'COMMON.REMD'
2126 weights(i)=hweights(iiset,i)