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,e_tmp
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 c broadcast iset to slaves
294 if (nfgtasks.gt.1) then
295 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
296 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
298 if (hremd.gt.0) call set_hweights(iset)
299 if(me.eq.king.or..not.out1file)
300 & write(iout,*) me,"iset=",iset,"t_bath=",t_bath
303 stdfp=dsqrt(2*Rb*t_bath/d_time)
305 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
308 c print *,'irep',me,t_bath
310 if (me.eq.king .or. .not. out1file)
311 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
312 call rescale_weights(t_bath)
316 c------copy MD--------------
317 c The driver for molecular dynamics subroutines
318 c------------------------------------------------
324 if(me.eq.king.or..not.out1file)
325 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
331 c Determine the inverse of the inertia matrix.
332 call setup_MD_matrices
336 if (me.eq.king .or. .not. out1file)
337 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
338 stdfp=dsqrt(2*Rb*t_bath/d_time)
340 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
342 if (lang.gt.0 .and. .not.surfarea) then
344 stdforcp(i)=stdfp*dsqrt(gamp)
347 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
349 elseif (lang.gt.0 .and. surfarea ) then
352 call rescale_weights(t_bath)
356 t_MDsetup = MPI_Wtime()-tt0
358 t_MDsetup = tcpu()-tt0
361 c Entering the MD loop
367 if (lang.eq.2 .or. lang.eq.3) then
371 call sd_verlet_p_setup
373 call sd_verlet_ciccotti_setup
377 pfric0_mat(i,j,0)=pfric_mat(i,j)
378 afric0_mat(i,j,0)=afric_mat(i,j)
379 vfric0_mat(i,j,0)=vfric_mat(i,j)
380 prand0_mat(i,j,0)=prand_mat(i,j)
381 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
382 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
387 flag_stoch(i)=.false.
391 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
393 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
397 else if (lang.eq.1 .or. lang.eq.4) then
401 if (me.eq.king .or. .not. out1file)
402 & write(iout,*) 'Setup time',time00-walltime
403 ctime call flush(iout)
405 t_langsetup=MPI_Wtime()-tt0
408 t_langsetup=tcpu()-tt0
413 do while(.not.end_of_run)
415 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
416 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
418 if (lang.gt.0 .and. surfarea .and.
419 & mod(itime,reset_fricmat).eq.0) then
420 if (lang.eq.2 .or. lang.eq.3) then
424 call sd_verlet_p_setup
426 call sd_verlet_ciccotti_setup
430 pfric0_mat(i,j,0)=pfric_mat(i,j)
431 afric0_mat(i,j,0)=afric_mat(i,j)
432 vfric0_mat(i,j,0)=vfric_mat(i,j)
433 prand0_mat(i,j,0)=prand_mat(i,j)
434 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
435 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
440 flag_stoch(i)=.false.
443 else if (lang.eq.1 .or. lang.eq.4) then
446 write (iout,'(a,i10)')
447 & "Friction matrix reset based on surface area, itime",itime
449 if (reset_vel .and. tbf .and. lang.eq.0
450 & .and. mod(itime,count_reset_vel).eq.0) then
452 if (me.eq.king .or. .not. out1file)
453 & write(iout,'(a,f20.2)')
454 & "Velocities reset to random values, time",totT
457 d_t_old(j,i)=d_t(j,i)
461 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
465 d_t(j,0)=d_t(j,0)-vcm(j)
468 kinetic_T=2.0d0/(dimen3*Rb)*EK
469 scalfac=dsqrt(T_bath/kinetic_T)
470 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
473 d_t_old(j,i)=scalfac*d_t(j,i)
479 c Time-reversible RESPA algorithm
480 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
481 call RESPA_step(itime)
483 c Variable time step algorithm.
484 call velverlet_step(itime)
488 call brown_step(itime)
490 print *,"Brown dynamics not here!"
492 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
497 if(hmc.gt.0 .and. mod(itime,hmc).eq.0) then
502 if (mod(itime,ntwe).eq.0) call statout(itime)
504 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
505 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
507 call pdbout(potE,tytul,ipdb)
512 if (mod(itime,ntwx).eq.0.and.traj1file) then
513 if(ntwx_cache.lt.max_cache_traj_use) then
514 ntwx_cache=ntwx_cache+1
516 if (max_cache_traj_use.ne.1)
517 & print *,itime,"processor ",me," over cache ",ntwx_cache
520 totT_cache(i)=totT_cache(i+1)
521 EK_cache(i)=EK_cache(i+1)
522 potE_cache(i)=potE_cache(i+1)
523 t_bath_cache(i)=t_bath_cache(i+1)
524 Uconst_cache(i)=Uconst_cache(i+1)
525 iset_cache(i)=iset_cache(i+1)
528 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
531 qpair_cache(ii,i)=qpair_cache(ii,i+1)
534 utheta_cache(ii,i)=utheta_cache(ii,i+1)
535 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
536 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
542 c_cache(j,ii,i)=c_cache(j,ii,i+1)
548 totT_cache(ntwx_cache)=totT
549 EK_cache(ntwx_cache)=EK
550 potE_cache(ntwx_cache)=potE
551 t_bath_cache(ntwx_cache)=t_bath
552 Uconst_cache(ntwx_cache)=Uconst
553 iset_cache(ntwx_cache)=iset
556 qfrag_cache(i,ntwx_cache)=qfrag(i)
559 qpair_cache(i,ntwx_cache)=qpair(i)
562 utheta_cache(i,ntwx_cache)=utheta(i)
563 ugamma_cache(i,ntwx_cache)=ugamma(i)
564 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
569 c_cache(j,i,ntwx_cache)=c(j,i)
574 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
575 & .and..not.restart1file) then
578 open(irest1,file=mremd_rst_name,status='unknown')
579 write (irest1,*) "i2rep"
580 write (irest1,*) (i2rep(i),i=0,nodes-1)
581 write (irest1,*) "ifirst"
582 write (irest1,*) (ifirst(i),i=1,remd_m(1))
584 write (irest1,*) "nupa",il
585 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
586 write (irest1,*) "ndowna",il
587 write (irest1,*) ndowna(0,il),
588 & (ndowna(i,il),i=1,ndowna(0,il))
590 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
591 write (irest1,*) "nset"
592 write (irest1,*) nset
593 write (irest1,*) "mset"
594 write (irest1,*) (mset(i),i=1,nset)
595 write (irest1,*) "i2set"
596 write (irest1,*) (i2set(i),i=0,nodes-1)
597 write (irest1,*) "i_index"
601 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
609 open(irest2,file=rest2name,status='unknown')
610 write(irest2,*) totT,EK,potE,totE,t_bath
612 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
615 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
617 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
618 write (irest2,*) iset
625 c forced synchronization
626 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
627 & .and. .not. mremdsync) then
629 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
631 call mpi_recv(itime_master, 1, MPI_INTEGER,
632 & 0,101,CG_COMM, status, ierr)
633 call mpi_barrier(CG_COMM, ierr)
634 cdeb if (out1file.or.traj1file) then
635 cdeb call mpi_gather(itime,1,mpi_integer,
636 cdeb & icache_all,1,mpi_integer,king,
639 & call mpi_gather(ntwx_cache,1,mpi_integer,
640 & icache_all,1,mpi_integer,king,
643 & write(iout,*) 'REMD synchro at',itime_master,itime
644 if (itime_master.ge.n_timestep .or. ovrtim())
646 ctime call flush(iout)
651 if ((mod(itime,nstex).eq.0.and.me.eq.king
652 & .or.end_of_run.and.me.eq.king )
653 & .and. .not. mremdsync ) then
657 call mpi_isend(itime,1,MPI_INTEGER,i,101,
658 & CG_COMM, ireqi(i), ierr)
659 cd write(iout,*) 'REMD synchro with',i
662 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
663 call mpi_barrier(CG_COMM, ierr)
665 write(iout,*) 'REMD synchro at',itime,'time=',time01-time01_
666 if (out1file.or.traj1file) then
667 cdeb call mpi_gather(itime,1,mpi_integer,
668 cdeb & itime_all,1,mpi_integer,king,
670 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
671 cdeb & (itime_all(i),i=1,nodes)
673 cdeb imin_itime=itime_all(1)
675 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
677 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
678 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
679 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
680 call mpi_gather(ntwx_cache,1,mpi_integer,
681 & icache_all,1,mpi_integer,king,
683 c write(iout,'(a19,8000i8)') ' ntwx_cache',
684 c & (icache_all(i),i=1,nodes)
685 ii_write=icache_all(1)
687 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
689 c write(iout,*) "MIN ii_write=",ii_write
692 ctime call flush(iout)
694 if(mremdsync .and. mod(itime,nstex).eq.0) then
696 if (me.eq.king .or. .not. out1file)
697 & write(iout,*) 'REMD synchro at',itime
700 call mpi_gather(ntwx_cache,1,mpi_integer,
701 & icache_all,1,mpi_integer,king,
704 write(iout,'(a19,8000i8)') ' ntwx_cache',
705 & (icache_all(i),i=1,nodes)
706 ii_write=icache_all(1)
708 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
710 write(iout,*) "MIN ii_write=",ii_write
713 ctest call flush(iout)
716 c Update the time safety limiy
717 if (time001-time00.gt.safety) then
718 safety=time001-time00+600
719 if (me.eq.king .or. .not. out1file)
720 & write (iout,*) "****** SAFETY increased to",safety," s"
722 if (ovrtim() .and. me.eq.king) end_of_run=.true.
723 call MPI_Bcast(end_of_run,1,MPI_LOGICAL,king,CG_COMM,IERR)
725 if(synflag.and..not.end_of_run) then
729 c write(iout,*) 'REMD before',me,t_bath
731 c call mpi_gather(t_bath,1,mpi_double_precision,
732 c & remd_t_bath,1,mpi_double_precision,king,
734 potEcomp(n_ene+1)=t_bath
736 if (usampl.or.homol_nset.gt.1) then
737 potEcomp(n_ene+2)=iset
738 c write(iout,*) potEcomp(n_ene+1),potEcomp(n_ene+2),iset,nset
739 if (iset.lt.nset) then
742 if (homol_nset.gt.1) then
743 c broadcast iset to slaves and reduce energy
744 if (nfgtasks.gt.1) then
745 call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
746 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
747 call e_modeller(e_tmp)
748 c write(iout,*) "iset+1 before reduce",e_tmp
749 call MPI_Barrier(FG_COMM,IERR)
750 call MPI_Reduce(e_tmp,potEcomp(n_ene+3),1,
751 & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
753 call e_modeller(potEcomp(n_ene+3))
755 c write(iout,*) "iset+1",potEcomp(n_ene+3)
758 potEcomp(n_ene+3)=Uconst
761 c broadcast iset to slaves
762 if (nfgtasks.gt.1) then
763 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
764 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
767 potEcomp(n_ene+3)=0.0
772 if (homol_nset.gt.1) then
773 c broadcast iset to slaves and reduce energy
774 if (nfgtasks.gt.1) then
775 call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
776 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
777 call e_modeller(e_tmp)
778 c write(iout,*) "iset-1 before reduce",e_tmp
779 call MPI_Barrier(FG_COMM,IERR)
780 call MPI_Reduce(e_tmp,potEcomp(n_ene+4),1,
781 & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
783 call e_modeller(potEcomp(n_ene+4))
785 c write(iout,*) "iset-1",potEcomp(n_ene+4)
788 potEcomp(n_ene+4)=Uconst
791 c broadcast iset to slaves
792 if (nfgtasks.gt.1) then
793 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
794 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
797 potEcomp(n_ene+4)=0.0
800 if(hremd.gt.0) potEcomp(n_ene+2)=iset
801 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
802 & remd_ene(0,1),n_ene+5,mpi_double_precision,king,
805 call mpi_gather(elow,1,mpi_double_precision,
806 & elowi,1,mpi_double_precision,king,
808 call mpi_gather(ehigh,1,mpi_double_precision,
809 & ehighi,1,mpi_double_precision,king,
814 if (me.eq.king .or. .not. out1file) then
815 write(iout,*) 'REMD gather times=',time03-time01
819 if (restart1file) call write1rst(i_index)
822 if (me.eq.king .or. .not. out1file) then
823 write(iout,*) 'REMD writing rst time=',time04-time03
826 if (traj1file) call write1traj
828 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
829 cdeb & icache_all,1,mpi_integer,king,
831 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
832 cdeb & (icache_all(i),i=1,nodes)
837 if (me.eq.king .or. .not. out1file) then
838 write(iout,*) 'REMD writing traj time=',time05-time04
839 ctime call flush(iout)
844 if(homol_nset.gt.1) write(iout,*)
845 & 'energy_c temperature iset energy_c(iset+1) energy_c(iset-1)'
847 remd_t_bath(i)=remd_ene(n_ene+1,i)
848 iremd_iset(i)=remd_ene(n_ene+2,i)
850 & write(iout,'(i4,f10.3,f6.0,i3,2f10.3)')
851 & i,remd_ene(i_econstr,i),
852 & remd_ene(n_ene+1,i),iremd_iset(i),
853 & remd_ene(n_ene+3,i),remd_ene(n_ene+4,i)
857 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
859 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
863 write(iout,*) 'REMD exchange temp,ene'
865 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
866 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
870 c-------------------------------------
871 IF(.not.usampl.and.hremd.eq.0.and.homol_nset.le.1) THEN
873 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
875 ctime call flush(iout)
876 write (iout,*) "remd_m(1)",remd_m(1)
879 i=ifirst(iran_num(1,remd_m(1)))
883 ctime call flush(iout)
888 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
890 if(i.gt.0.and.nupa(0,i).gt.0) then
892 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
894 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
896 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
898 c do while (iex.eq.i)
899 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
900 iex=nupa(iran_num(1,int(nupa(0,i))),i)
902 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
904 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
906 c Swap temperatures between conformations i and iex with recalculating the free energies
907 c following temperature changes.
908 ene_iex_iex=remd_ene(0,iex)
909 ene_i_i=remd_ene(0,i)
910 c write (iout,*) "i",i," ene_i_i",ene_i_i,
911 c & " iex",iex," ene_iex_iex",ene_iex_iex
912 c write (iout,*) "rescaling weights with temperature",
915 call rescale_weights(remd_t_bath(i))
917 c write (iout,*) "0,iex",remd_t_bath(i)
918 c call enerprint(remd_ene(0,iex))
920 call sum_energy(remd_ene(0,iex),.false.)
921 ene_iex_i=remd_ene(0,iex)
922 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
924 c write (iout,*) "0,i",remd_t_bath(i)
925 c call enerprint(remd_ene(0,i))
927 call sum_energy(remd_ene(0,i),.false.)
928 c write (iout,*) "ene_i_i",remd_ene(0,i)
930 c write (iout,*) "rescaling weights with temperature",
932 if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
933 write (iout,*) "ERROR: inconsistent energies:",i,
934 & ene_i_i,remd_ene(0,i)
936 call rescale_weights(remd_t_bath(iex))
938 c write (iout,*) "0,i",remd_t_bath(iex)
939 c call enerprint(remd_ene(0,i))
941 call sum_energy(remd_ene(0,i),.false.)
942 c write (iout,*) "ene_i_iex",remd_ene(0,i)
944 ene_i_iex=remd_ene(0,i)
946 c write (iout,*) "0,iex",remd_t_bath(iex)
947 c call enerprint(remd_ene(0,iex))
949 call sum_energy(remd_ene(0,iex),.false.)
950 if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
951 write (iout,*) "ERROR: inconsistent energies:",iex,
952 & ene_iex_iex,remd_ene(0,iex)
954 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
955 c write (iout,*) "i",i," iex",iex
956 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
957 c & " ene_i_iex",ene_i_iex,
958 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
960 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
961 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
963 c write(iout,*) 'delta',delta
964 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
965 c & (remd_ene(i)-remd_ene(iex))/Rb/
966 c & (remd_t_bath(i)*remd_t_bath(iex))
968 if (delta .gt. 50.0d0) then
974 else if (delta.lt.-50.0d0) then
983 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
984 xxx=ran_number(0.0d0,1.0d0)
985 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
987 if (delta .gt. xxx) then
989 remd_t_bath(i)=remd_t_bath(iex)
991 remd_ene(0,i)=ene_i_iex
992 remd_ene(0,iex)=ene_iex_i
998 ehighi(i)=ehighi(iex)
1005 nupa(k,i)=nupa(k,iex)
1008 ndowna(k,i)=ndowna(k,iex)
1012 if (ifirst(il).eq.i) ifirst(il)=iex
1014 if (nupa(k,il).eq.i) then
1016 elseif (nupa(k,il).eq.iex) then
1021 if (ndowna(k,il).eq.i) then
1023 elseif (ndowna(k,il).eq.iex) then
1029 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1031 i2rep(i-1)=i2rep(iex-1)
1034 c write(iout,*) 'exchange',i,iex
1035 c write (iout,'(a8,100i4)') "@ ifirst",
1036 c & (ifirst(k),k=1,remd_m(1))
1038 c write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1039 c & (nupa(k,il),k=1,nupa(0,il))
1040 c write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1041 c & (ndowna(k,il),k=1,ndowna(0,il))
1046 remd_ene(0,iex)=ene_iex_iex
1047 remd_ene(0,i)=ene_i_i
1053 cd write (iout,*) "exchange completed"
1055 ELSEIF (usampl.or.homol_nset.gt.1) THEN
1057 c write(iout,*) "########",ii
1059 i_temp=iran_num(1,nrep)
1060 i_mult=iran_num(1,remd_m(i_temp))
1061 i_iset=iran_num(1,nset)
1062 i_mset=iran_num(1,mset(i_iset))
1063 i=i_index(i_temp,i_mult,i_iset,i_mset)
1065 c write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1067 if(t_exchange_only)then
1072 c write(iout,*) "i_dir=",i_dir
1074 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1077 i_mult1=iran_num(1,remd_m(i_temp1))
1079 i_mset1=iran_num(1,mset(i_iset1))
1080 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1082 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1085 i_mult1=iran_num(1,remd_m(i_temp1))
1087 i_mset1=iran_num(1,mset(i_iset1))
1088 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1090 econstr_temp_i=remd_ene(i_econstr,i)
1091 econstr_temp_iex=remd_ene(i_econstr,iex)
1092 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1093 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1095 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1098 i_mult1=iran_num(1,remd_m(i_temp1))
1100 i_mset1=iran_num(1,mset(i_iset1))
1101 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1102 econstr_temp_i=remd_ene(i_econstr,i)
1103 econstr_temp_iex=remd_ene(i_econstr,iex)
1104 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1105 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1111 c write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1112 ctime call flush(iout)
1114 c Swap temperatures between conformations i and iex with recalculating the free energies
1115 c following temperature changes.
1116 ene_iex_iex=remd_ene(0,iex)
1117 ene_i_i=remd_ene(0,i)
1118 co write (iout,*) "rescaling weights with temperature",
1120 call rescale_weights(remd_t_bath(i))
1122 call sum_energy(remd_ene(0,iex),.false.)
1123 ene_iex_i=remd_ene(0,iex)
1125 c ERROR only makes sense for dir =1
1126 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
1127 c call sum_energy(remd_ene(0,i),.false.)
1128 c write (iout,*) "ene_i_i",remd_ene(0,i)
1129 c write (iout,*) "rescaling weights with temperature",
1130 c & remd_t_bath(iex)
1131 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1132 c write (iout,*) "ERROR: inconsistent energies i:",i,
1133 c & ene_i_i,remd_ene(0,i)
1136 call rescale_weights(remd_t_bath(iex))
1137 call sum_energy(remd_ene(0,i),.false.)
1138 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1139 ene_i_iex=remd_ene(0,i)
1141 c ERROR only makes sense for dir =1
1142 c call sum_energy(remd_ene(0,iex),.false.)
1143 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1144 c write (iout,*) "ERROR: inconsistent energies iex:",iex,
1145 c & ene_iex_iex,remd_ene(0,iex)
1147 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1148 c write (iout,*) "i",i," iex",iex
1149 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1150 c & " ene_i_iex",ene_i_iex,
1151 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1153 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1154 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1156 c write(iout,*) 'delta',delta
1157 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1158 c & (remd_ene(i)-remd_ene(iex))/Rb/
1159 c & (remd_t_bath(i)*remd_t_bath(iex))
1160 if (delta .gt. 50.0d0) then
1165 if (i_dir.eq.1.or.i_dir.eq.3)
1166 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1167 if (i_dir.eq.2.or.i_dir.eq.3)
1168 & iremd_tot_usa(int(i2set(i-1)))=
1169 & iremd_tot_usa(int(i2set(i-1)))+1
1170 xxx=ran_number(0.0d0,1.0d0)
1171 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1172 if (delta .gt. xxx) then
1174 remd_t_bath(i)=remd_t_bath(iex)
1175 remd_t_bath(iex)=tmp
1178 iremd_iset(i)=iremd_iset(iex)
1179 iremd_iset(iex)=itmp
1181 remd_ene(0,i)=ene_i_iex
1182 remd_ene(0,iex)=ene_iex_i
1184 if (i_dir.eq.1.or.i_dir.eq.3)
1185 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1188 i2rep(i-1)=i2rep(iex-1)
1191 if (i_dir.eq.2.or.i_dir.eq.3)
1192 & iremd_acc_usa(int(i2set(i-1)))=
1193 & iremd_acc_usa(int(i2set(i-1)))+1
1196 i2set(i-1)=i2set(iex-1)
1199 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1200 i_index(i_temp,i_mult,i_iset,i_mset)=
1201 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1202 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1205 remd_ene(0,iex)=ene_iex_iex
1206 remd_ene(0,i)=ene_i_i
1207 remd_ene(i_econstr,iex)=econstr_temp_iex
1208 remd_ene(i_econstr,i)=econstr_temp_i
1212 cd do il1=1,mset(il)
1215 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1225 ELSEIF (hremd.gt.0) THEN
1227 cd write(iout,*) "########",ii
1229 i_temp=iran_num(1,nrep)
1230 i_mult=iran_num(1,remd_m(i_temp))
1231 i_iset=iran_num(1,nset)
1233 i=i_index(i_temp,i_mult,i_iset,i_mset)
1235 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1237 if(t_exchange_only)then
1243 cd write(iout,*) "i_dir=",i_dir
1245 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1248 i_mult1=iran_num(1,remd_m(i_temp1))
1251 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1253 elseif(i_dir.eq.2)then
1256 i_mult1=iran_num(1,remd_m(i_temp1))
1257 i_iset1=iran_num(1,hremd)
1258 do while(i_iset1.eq.i_iset)
1259 i_iset1=iran_num(1,hremd)
1262 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1264 elseif(remd_m(i_temp+1).gt.0)then
1267 i_mult1=iran_num(1,remd_m(i_temp1))
1268 i_iset1=iran_num(1,hremd)
1269 do while(i_iset1.eq.i_iset)
1270 i_iset1=iran_num(1,hremd)
1273 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1279 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1280 ctime call flush(iout)
1282 c Swap temperatures between conformations i and iex with recalculating the free energies
1283 c following temperature changes.
1284 ene_iex_iex=remd_ene(0,iex)
1285 ene_i_i=remd_ene(0,i)
1287 call set_hweights(i_iset)
1288 call rescale_weights(remd_t_bath(i))
1289 call sum_energy(remd_ene(0,iex),.false.)
1290 ene_iex_i=remd_ene(0,iex)
1292 call set_hweights(i_iset1)
1293 call rescale_weights(remd_t_bath(iex))
1294 call sum_energy(remd_ene(0,i),.false.)
1295 ene_i_iex=remd_ene(0,i)
1297 cd write(iout,*) ene_iex_iex,ene_i_i,ene_iex_i,ene_i_iex
1299 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1300 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1303 if (delta .gt. 50.0d0) then
1309 if (i_dir.eq.1.or.i_dir.eq.3)
1310 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1311 if (i_dir.eq.2.or.i_dir.eq.3)
1312 & iremd_tot_usa(int(i2set(i-1)))=
1313 & iremd_tot_usa(int(i2set(i-1)))+1
1314 xxx=ran_number(0.0d0,1.0d0)
1315 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1316 if (delta .gt. xxx) then
1318 cd write (iout,*) "exchange"
1320 remd_t_bath(i)=remd_t_bath(iex)
1321 remd_t_bath(iex)=tmp
1324 iremd_iset(i)=iremd_iset(iex)
1325 iremd_iset(iex)=itmp
1327 remd_ene(0,i)=ene_i_iex
1328 remd_ene(0,iex)=ene_iex_i
1330 if (i_dir.eq.1.or.i_dir.eq.3)
1331 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1334 i2rep(i-1)=i2rep(iex-1)
1337 if (i_dir.eq.2.or.i_dir.eq.3)
1338 & iremd_acc_usa(int(i2set(i-1)))=
1339 & iremd_acc_usa(int(i2set(i-1)))+1
1342 i2set(i-1)=i2set(iex-1)
1345 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1346 i_index(i_temp,i_mult,i_iset,i_mset)=
1347 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1348 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1351 cd do il1=1,mset(il)
1354 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1361 remd_ene(0,iex)=ene_iex_iex
1362 remd_ene(0,i)=ene_i_i
1373 c-------------------------------------
1374 write (iout,*) "NREP",nrep
1376 if(iremd_tot(i).ne.0)
1377 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1378 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1381 if(usampl.or.homol_nset.gt.1) then
1383 if(iremd_tot_usa(i).ne.0)
1384 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1385 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1391 if(iremd_tot_usa(i).ne.0)
1392 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_hremd',i,
1393 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1398 ctime call flush(iout)
1400 cd write (iout,'(a6,100i4)') "ifirst",
1401 cd & (ifirst(i),i=1,remd_m(1))
1403 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1404 cd & (nupa(i,il),i=1,nupa(0,il))
1405 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1406 cd & (ndowna(i,il),i=1,ndowna(0,il))
1411 cd write (iout,*) "Before scatter"
1413 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1414 & t_bath,1,mpi_double_precision,king,
1416 cd write (iout,*) "After scatter"
1418 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
1419 call mpi_scatter(iremd_iset,1,mpi_integer,
1420 & iset,1,mpi_integer,king,
1422 c 8/31/2015 Correction by AL: send new iset to slaves
1423 if (nfgtasks.gt.1) then
1424 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
1425 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
1431 if (me.eq.king .or. .not. out1file) then
1432 write(iout,*) 'REMD scatter time=',time07-time06
1436 call mpi_scatter(elowi,1,mpi_double_precision,
1437 & elow,1,mpi_double_precision,king,
1439 call mpi_scatter(ehighi,1,mpi_double_precision,
1440 & ehigh,1,mpi_double_precision,king,
1444 if(hremd.gt.0) call set_hweights(iset)
1445 call rescale_weights(t_bath)
1446 co write (iout,*) "Processor",me,
1447 co & " rescaling weights with temperature",t_bath
1449 stdfp=dsqrt(2*Rb*t_bath/d_time)
1451 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1455 stdforcp(i)=stdforcp(i)*sqrt(t_bath/t_bath_old)
1458 stdforcsc(i)=stdforcsc(i)*sqrt(t_bath/t_bath_old)
1461 cde write(iout,*) 'REMD after',me,t_bath
1463 if (me.eq.king .or. .not. out1file) then
1464 write(iout,*) 'REMD exchange time=',time08-time02
1465 ctime call flush(iout)
1470 if (restart1file) then
1471 if (me.eq.king .or. .not. out1file)
1472 & write(iout,*) 'writing restart at the end of run'
1473 call write1rst(i_index)
1476 if (traj1file) call write1traj
1478 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1479 cdeb & icache_all,1,mpi_integer,king,
1480 cdeb & CG_COMM,ierr)
1481 cdeb write(iout,'(a40,8000i8)')
1482 cdeb & ' ntwx_cache after traj1file at the end',
1483 cdeb & (icache_all(i),i=1,nodes)
1488 t_MD=MPI_Wtime()-tt0
1492 if (me.eq.king .or. .not. out1file) then
1493 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1495 & 'MD calculations setup:',t_MDsetup,
1496 & 'Energy & gradient evaluation:',t_enegrad,
1497 & 'Stochastic MD setup:',t_langsetup,
1498 & 'Stochastic MD step setup:',t_sdsetup,
1500 write (iout,'(/28(1h=),a25,27(1h=))')
1501 & ' End of MD calculation '
1502 if(hmc.gt.0) write (iout,*) 'HMC acceptance ratio',
1503 & n_timestep*1.0d0/hmc/hmc_acc
1508 c-----------------------------------------------------------------------
1509 subroutine write1rst(i_index)
1510 implicit real*8 (a-h,o-z)
1511 include 'DIMENSIONS'
1514 include 'COMMON.IOUNITS'
1515 include 'COMMON.REMD'
1516 include 'COMMON.SETUP'
1517 include 'COMMON.CHAIN'
1518 include 'COMMON.SBRIDGE'
1519 include 'COMMON.INTERACT'
1520 include 'COMMON.CONTROL'
1522 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1523 & d_restart2(3,2*maxres*maxprocs)
1527 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1528 common /przechowalnia/ d_restart1,d_restart2
1533 t5_restart1(4)=t_bath
1534 t5_restart1(5)=Uconst
1536 call mpi_gather(t5_restart1,5,mpi_real,
1537 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1545 call mpi_gather(r_d,3*2*nres,mpi_real,
1546 & d_restart1,3*2*nres,mpi_real,king,
1555 call mpi_gather(r_d,3*2*nres,mpi_real,
1556 & d_restart2,3*2*nres,mpi_real,king,
1561 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1563 call xdrfint_(ixdrf, i2rep(i), iret)
1566 call xdrfint_(ixdrf, ifirst(i), iret)
1570 call xdrfint_(ixdrf, nupa(i,il), iret)
1574 call xdrfint_(ixdrf, ndowna(i,il), iret)
1580 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1587 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1594 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1599 if(usampl.or.homol_nset.gt.1) then
1600 call xdrfint_(ixdrf, nset, iret)
1602 call xdrfint_(ixdrf,mset(i), iret)
1605 call xdrfint_(ixdrf,i2set(i), iret)
1611 itmp=i_index(i,j,il,il1)
1612 call xdrfint_(ixdrf,itmp, iret)
1619 call xdrfclose_(ixdrf, iret)
1621 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1623 call xdrfint(ixdrf, i2rep(i), iret)
1626 call xdrfint(ixdrf, ifirst(i), iret)
1630 call xdrfint(ixdrf, nupa(i,il), iret)
1634 call xdrfint(ixdrf, ndowna(i,il), iret)
1640 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1647 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1654 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1660 if(usampl.or.homol_nset.gt.1) then
1661 call xdrfint(ixdrf, nset, iret)
1663 call xdrfint(ixdrf,mset(i), iret)
1666 call xdrfint(ixdrf,i2set(i), iret)
1672 itmp=i_index(i,j,il,il1)
1673 call xdrfint(ixdrf,itmp, iret)
1680 call xdrfclose(ixdrf, iret)
1687 subroutine write1traj
1688 implicit real*8 (a-h,o-z)
1689 include 'DIMENSIONS'
1692 include 'COMMON.IOUNITS'
1693 include 'COMMON.REMD'
1694 include 'COMMON.SETUP'
1695 include 'COMMON.CHAIN'
1696 include 'COMMON.SBRIDGE'
1697 include 'COMMON.INTERACT'
1701 real xcoord(3,maxres2+2),prec
1702 real r_qfrag(50),r_qpair(100)
1703 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1704 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1705 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1706 & p_uscdiff(100*maxprocs)
1707 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1708 common /przechowalnia/ p_c
1710 call mpi_bcast(ii_write,1,mpi_integer,
1711 & king,CG_COMM,ierr)
1714 print *,'traj1file',me,ii_write,ntwx_cache
1718 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1720 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1723 t5_restart1(1)=totT_cache(ii)
1724 t5_restart1(2)=EK_cache(ii)
1725 t5_restart1(3)=potE_cache(ii)
1726 t5_restart1(4)=t_bath_cache(ii)
1727 t5_restart1(5)=Uconst_cache(ii)
1728 call mpi_gather(t5_restart1,5,mpi_real,
1729 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1731 call mpi_gather(iset_cache(ii),1,mpi_integer,
1732 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1735 r_qfrag(i)=qfrag_cache(i,ii)
1738 r_qpair(i)=qpair_cache(i,ii)
1741 r_utheta(i)=utheta_cache(i,ii)
1742 r_ugamma(i)=ugamma_cache(i,ii)
1743 r_uscdiff(i)=uscdiff_cache(i,ii)
1746 call mpi_gather(r_qfrag,nfrag,mpi_real,
1747 & p_qfrag,nfrag,mpi_real,king,
1749 call mpi_gather(r_qpair,npair,mpi_real,
1750 & p_qpair,npair,mpi_real,king,
1752 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1753 & p_utheta,nfrag_back,mpi_real,king,
1755 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1756 & p_ugamma,nfrag_back,mpi_real,king,
1758 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1759 & p_uscdiff,nfrag_back,mpi_real,king,
1763 write (iout,*) "p_qfrag"
1765 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1767 write (iout,*) "p_qpair"
1769 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1771 ctime call flush(iout)
1775 r_c(j,i)=c_cache(j,i,ii)
1779 call mpi_gather(r_c,3*2*nres,mpi_real,
1780 & p_c,3*2*nres,mpi_real,king,
1786 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1787 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1788 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1789 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1790 call xdrfint_(ixdrf, nss, iret)
1793 call xdrfint_(ixdrf, idssb(j)+nres, iret)
1794 call xdrfint_(ixdrf, jdssb(j)+nres, iret)
1796 call xdrfint_(ixdrf, ihpb(j), iret)
1797 call xdrfint_(ixdrf, jhpb(j), iret)
1800 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1801 call xdrfint_(ixdrf, iset_restart1(il), iret)
1803 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1806 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1809 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1810 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1811 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1816 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1821 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1825 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1829 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1830 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1831 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1832 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1833 call xdrfint(ixdrf, nss, iret)
1836 call xdrfint(ixdrf, idssb(j)+nres, iret)
1837 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1839 call xdrfint(ixdrf, ihpb(j), iret)
1840 call xdrfint(ixdrf, jhpb(j), iret)
1843 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1844 call xdrfint(ixdrf, iset_restart1(il), iret)
1846 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1849 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1852 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1853 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1854 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1859 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1864 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1868 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1874 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1876 if(me.eq.king) call xdrfclose(ixdrf, iret)
1878 do i=1,ntwx_cache-ii_write
1880 totT_cache(i)=totT_cache(ii_write+i)
1881 EK_cache(i)=EK_cache(ii_write+i)
1882 potE_cache(i)=potE_cache(ii_write+i)
1883 t_bath_cache(i)=t_bath_cache(ii_write+i)
1884 Uconst_cache(i)=Uconst_cache(ii_write+i)
1885 iset_cache(i)=iset_cache(ii_write+i)
1888 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1891 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1894 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1895 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1896 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1901 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1905 ntwx_cache=ntwx_cache-ii_write
1910 subroutine read1restart(i_index)
1911 implicit real*8 (a-h,o-z)
1912 include 'DIMENSIONS'
1915 include 'COMMON.IOUNITS'
1916 include 'COMMON.REMD'
1917 include 'COMMON.SETUP'
1918 include 'COMMON.CHAIN'
1919 include 'COMMON.SBRIDGE'
1920 include 'COMMON.INTERACT'
1921 include 'COMMON.CONTROL'
1922 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1925 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1926 common /przechowalnia/ d_restart1
1927 write (*,*) "Processor",me," called read1restart"
1930 open(irest2,file=mremd_rst_name,status='unknown')
1931 read(irest2,*,err=334) i
1932 write(iout,*) "Reading old rst in ASCI format"
1934 call read1restart_old
1938 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1941 call xdrfint_(ixdrf, i2rep(i), iret)
1944 call xdrfint_(ixdrf, ifirst(i), iret)
1947 call xdrfint_(ixdrf, nupa(0,il), iret)
1949 call xdrfint_(ixdrf, nupa(i,il), iret)
1952 call xdrfint_(ixdrf, ndowna(0,il), iret)
1954 call xdrfint_(ixdrf, ndowna(i,il), iret)
1959 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1963 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1966 call xdrfint(ixdrf, i2rep(i), iret)
1969 call xdrfint(ixdrf, ifirst(i), iret)
1972 call xdrfint(ixdrf, nupa(0,il), iret)
1974 call xdrfint(ixdrf, nupa(i,il), iret)
1977 call xdrfint(ixdrf, ndowna(0,il), iret)
1979 call xdrfint(ixdrf, ndowna(i,il), iret)
1984 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1989 call mpi_scatter(t_restart1,5,mpi_real,
1990 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1994 t_bath=t5_restart1(4)
1999 c read(irest2,'(3e15.5)')
2000 c & (d_restart1(j,i+2*nres*il),j=1,3)
2003 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
2005 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
2010 write (iout,*) "Conformation read",il
2012 write (iout,'(i5,3f10.5,5x,3f10.5)')
2013 & i,(d_restart1(j,i+2*nres*il),j=1,3),
2014 & (d_restart1(j,nres+i+2*nres*il),j=1,3)
2019 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2020 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2030 c read(irest2,'(3e15.5)')
2031 c & (d_restart1(j,i+2*nres*il),j=1,3)
2034 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
2036 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
2042 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2043 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2051 if(usampl.or.homol_nset.gt.1) then
2054 call xdrfint_(ixdrf, nset, iret)
2056 call xdrfint_(ixdrf,mset(i), iret)
2059 call xdrfint_(ixdrf,i2set(i), iret)
2065 call xdrfint_(ixdrf,itmp, iret)
2066 i_index(i,j,il,il1)=itmp
2074 call xdrfint(ixdrf, nset, iret)
2076 call xdrfint(ixdrf,mset(i), iret)
2079 call xdrfint(ixdrf,i2set(i), iret)
2085 call xdrfint(ixdrf,itmp, iret)
2086 i_index(i,j,il,il1)=itmp
2093 c Corrected AL 8/19/2014: each processor needs whole iset array not only its
2095 c call mpi_scatter(i2set,1,mpi_integer,
2096 c & iset,1,mpi_integer,king,
2098 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
2101 c broadcast iset to slaves
2102 if (nfgtasks.gt.1) then
2103 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
2104 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
2109 if(me.eq.king) close(irest2)
2113 subroutine read1restart_old
2114 implicit real*8 (a-h,o-z)
2115 include 'DIMENSIONS'
2118 include 'COMMON.IOUNITS'
2119 include 'COMMON.REMD'
2120 include 'COMMON.SETUP'
2121 include 'COMMON.CHAIN'
2122 include 'COMMON.SBRIDGE'
2123 include 'COMMON.INTERACT'
2124 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
2126 common /przechowalnia/ d_restart1
2128 open(irest2,file=mremd_rst_name,status='unknown')
2129 read (irest2,*) (i2rep(i),i=0,nodes-1)
2130 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2132 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2133 read (irest2,*) ndowna(0,il),
2134 & (ndowna(i,il),i=1,ndowna(0,il))
2137 read(irest2,*) (t_restart1(j,il),j=1,4)
2140 call mpi_scatter(t_restart1,5,mpi_real,
2141 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2145 t_bath=t5_restart1(4)
2150 read(irest2,'(3e15.5)')
2151 & (d_restart1(j,i+2*nres*il),j=1,3)
2155 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2156 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2166 read(irest2,'(3e15.5)')
2167 & (d_restart1(j,i+2*nres*il),j=1,3)
2171 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2172 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2178 if(me.eq.king) close(irest2)
2181 c-------------------------------------------------------------------
2182 subroutine set_hweights(iiset)
2183 implicit real*8 (a-h,o-z)
2185 include 'DIMENSIONS'
2186 include 'COMMON.FFIELD'
2187 include 'COMMON.REMD'
2190 weights(i)=hweights(iiset,i)