4 implicit real*8 (a-h,o-z)
7 include 'COMMON.CONTROL'
11 include 'COMMON.LANGEVIN'
13 include 'COMMON.LANGEVIN.lang0'
15 include 'COMMON.CHAIN'
16 include 'COMMON.DERIV'
18 include 'COMMON.LOCAL'
19 include 'COMMON.INTERACT'
20 include 'COMMON.IOUNITS'
21 include 'COMMON.NAMES'
22 include 'COMMON.TIME1'
24 include 'COMMON.SETUP'
27 double precision cm(3),L(3),vcm(3)
28 double precision energia(0:n_ene)
29 double precision remd_t_bath(maxprocs)
30 integer iremd_iset(maxprocs)
32 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
33 double precision remd_ene(0:n_ene+4,maxprocs),t_bath_old,e_tmp
34 integer iremd_acc(maxprocs),iremd_tot(maxprocs)
35 integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
41 cold integer nup(0:maxprocs),ndown(0:maxprocs)
42 integer rep2i(0:maxprocs),ireqi(maxprocs)
43 integer icache_all(maxprocs)
44 integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
45 logical synflag,end_of_run,file_exist /.false./,ovrtim
52 if(me.eq.king.or..not.out1file) then
53 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
54 write (iout,*) "NREP=",nrep
58 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
59 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
61 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
63 cd print *,'MREMD',nodes,homol_nset
64 cd print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
65 cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath
74 if(homol_nset.gt.1) then
82 if(usampl) i_econstr=20
87 do il1=1,max0(mset(il),1)
103 if(me.eq.king.or..not.out1file) then
104 write(iout,*) "i2rep",(i2rep(i),i=0,nodes-1)
105 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
106 write(iout,*) "i,j,il,il1,i_index(i,j,il,il1)"
111 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
118 c print *,'i2rep',me,i2rep(me)
119 c print *,'rep2i',(rep2i(i),i=0,nrep)
121 cold if (i2rep(me).eq.nrep) then
124 cold nup(0)=remd_m(i2rep(me)+1)
125 cold k=rep2i(int(i2rep(me)))+1
132 cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
134 cold if (i2rep(me).eq.1) then
137 cold ndown(0)=remd_m(i2rep(me)-1)
138 cold k=rep2i(i2rep(me)-2)+1
145 cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
148 write (*,*) "Processor",me," rest",rest,"
149 & restart1fie",restart1file
150 if(rest.and.restart1file) then
152 & inquire(file=mremd_rst_name,exist=file_exist)
153 cd write (*,*) me," Before broadcast: file_exist",file_exist
154 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
156 cd write (*,*) me," After broadcast: file_exist",file_exist
158 if(me.eq.king.or..not.out1file)
159 & write (iout,*) 'Master is reading restart1file'
160 call read1restart(i_index)
162 if(me.eq.king.or..not.out1file)
163 & write (iout,*) 'WARNING : no restart1file'
166 if(me.eq.king.or..not.out1file) then
167 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
168 write(iout,*) "i_index"
173 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
182 if (rest.and..not.restart1file)
183 & inquire(file=mremd_rst_name,exist=file_exist)
184 if(.not.file_exist.and.rest.and..not.restart1file)
185 & write(iout,*) 'WARNING : no restart file',mremd_rst_name
186 IF (rest.and.file_exist.and..not.restart1file) THEN
187 write (iout,*) 'Master is reading restart file',
189 open(irest2,file=mremd_rst_name,status='unknown')
191 read (irest2,*) (i2rep(i),i=0,nodes-1)
193 read (irest2,*) (ifirst(i),i=1,remd_m(1))
196 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
198 read (irest2,*) ndowna(0,il),
199 & (ndowna(i,il),i=1,ndowna(0,il))
201 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
205 read (irest2,*) (mset(i),i=1,nset)
207 read (irest2,*) (i2set(i),i=0,nodes-1)
212 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
217 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
218 write(iout,*) "i_index"
223 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
232 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
233 write (iout,'(a6,1000i5)') "ifirst",
234 & (ifirst(i),i=1,remd_m(1))
236 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
237 & (nupa(i,il),i=1,nupa(0,il))
238 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
239 & (ndowna(i,il),i=1,ndowna(0,il))
241 ELSE IF (.not.(rest.and.file_exist)) THEN
247 if (i2rep(il-1).eq.nrep) then
250 nupa(0,il)=remd_m(i2rep(il-1)+1)
251 k=rep2i(int(i2rep(il-1)))+1
257 if (i2rep(il-1).eq.1) then
260 ndowna(0,il)=remd_m(i2rep(il-1)-1)
261 k=rep2i(i2rep(il-1)-2)+1
269 write (iout,'(a6,100i4)') "ifirst",
270 & (ifirst(i),i=1,remd_m(1))
272 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
273 & (nupa(i,il),i=1,nupa(0,il))
274 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
275 & (ndowna(i,il),i=1,ndowna(0,il))
281 c t_bath=retmin+(retmax-retmin)*me/(nodes-1)
282 if(.not.(rest.and.file_exist.and.restart1file)) then
283 if (me .eq. king) then
286 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
288 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
289 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
292 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
294 c broadcast iset to slaves
295 if (nfgtasks.gt.1) then
296 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
297 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
299 if (hremd.gt.0) call set_hweights(iset)
300 if(me.eq.king.or..not.out1file)
301 & write(iout,*) me,"iset=",iset,"t_bath=",t_bath
304 stdfp=dsqrt(2*Rb*t_bath/d_time)
306 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
309 c print *,'irep',me,t_bath
311 if (me.eq.king .or. .not. out1file)
312 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
313 call rescale_weights(t_bath)
317 c------copy MD--------------
318 c The driver for molecular dynamics subroutines
319 c------------------------------------------------
325 if(me.eq.king.or..not.out1file)
326 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
332 c Determine the inverse of the inertia matrix.
333 call setup_MD_matrices
337 if (me.eq.king .or. .not. out1file)
338 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
339 stdfp=dsqrt(2*Rb*t_bath/d_time)
341 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
343 if (lang.gt.0 .and. .not.surfarea) then
345 stdforcp(i)=stdfp*dsqrt(gamp)
348 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
350 elseif (lang.gt.0 .and. surfarea ) then
353 call rescale_weights(t_bath)
357 t_MDsetup = MPI_Wtime()-tt0
359 t_MDsetup = tcpu()-tt0
362 c Entering the MD loop
368 if (lang.eq.2 .or. lang.eq.3) then
372 call sd_verlet_p_setup
374 call sd_verlet_ciccotti_setup
378 pfric0_mat(i,j,0)=pfric_mat(i,j)
379 afric0_mat(i,j,0)=afric_mat(i,j)
380 vfric0_mat(i,j,0)=vfric_mat(i,j)
381 prand0_mat(i,j,0)=prand_mat(i,j)
382 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
383 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
388 flag_stoch(i)=.false.
392 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
394 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
398 else if (lang.eq.1 .or. lang.eq.4) then
402 if (me.eq.king .or. .not. out1file)
403 & write(iout,*) 'Setup time',time00-walltime
404 ctime call flush(iout)
406 t_langsetup=MPI_Wtime()-tt0
409 t_langsetup=tcpu()-tt0
414 do while(.not.end_of_run)
416 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
417 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
419 if (lang.gt.0 .and. surfarea .and.
420 & mod(itime,reset_fricmat).eq.0) then
421 if (lang.eq.2 .or. lang.eq.3) then
425 call sd_verlet_p_setup
427 call sd_verlet_ciccotti_setup
431 pfric0_mat(i,j,0)=pfric_mat(i,j)
432 afric0_mat(i,j,0)=afric_mat(i,j)
433 vfric0_mat(i,j,0)=vfric_mat(i,j)
434 prand0_mat(i,j,0)=prand_mat(i,j)
435 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
436 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
441 flag_stoch(i)=.false.
444 else if (lang.eq.1 .or. lang.eq.4) then
447 write (iout,'(a,i10)')
448 & "Friction matrix reset based on surface area, itime",itime
450 if (reset_vel .and. tbf .and. lang.eq.0
451 & .and. mod(itime,count_reset_vel).eq.0) then
453 if (me.eq.king .or. .not. out1file)
454 & write(iout,'(a,f20.2)')
455 & "Velocities reset to random values, time",totT
458 d_t_old(j,i)=d_t(j,i)
462 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
466 d_t(j,0)=d_t(j,0)-vcm(j)
469 kinetic_T=2.0d0/(dimen3*Rb)*EK
470 scalfac=dsqrt(T_bath/kinetic_T)
471 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
474 d_t_old(j,i)=scalfac*d_t(j,i)
480 c Time-reversible RESPA algorithm
481 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
482 call RESPA_step(itime)
484 c Variable time step algorithm.
485 call velverlet_step(itime)
489 call brown_step(itime)
491 print *,"Brown dynamics not here!"
493 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
498 if(hmc.gt.0 .and. mod(itime,hmc).eq.0) then
503 if (mod(itime,ntwe).eq.0) call statout(itime)
505 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
506 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
508 call pdbout(potE,tytul,ipdb)
513 if (mod(itime,ntwx).eq.0.and.traj1file) then
514 if(ntwx_cache.lt.max_cache_traj_use) then
515 ntwx_cache=ntwx_cache+1
517 if (max_cache_traj_use.ne.1)
518 & print *,itime,"processor ",me," over cache ",ntwx_cache
521 totT_cache(i)=totT_cache(i+1)
522 EK_cache(i)=EK_cache(i+1)
523 potE_cache(i)=potE_cache(i+1)
524 t_bath_cache(i)=t_bath_cache(i+1)
525 Uconst_cache(i)=Uconst_cache(i+1)
526 iset_cache(i)=iset_cache(i+1)
529 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
532 qpair_cache(ii,i)=qpair_cache(ii,i+1)
535 utheta_cache(ii,i)=utheta_cache(ii,i+1)
536 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
537 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
543 c_cache(j,ii,i)=c_cache(j,ii,i+1)
549 totT_cache(ntwx_cache)=totT
550 EK_cache(ntwx_cache)=EK
551 potE_cache(ntwx_cache)=potE
552 t_bath_cache(ntwx_cache)=t_bath
553 Uconst_cache(ntwx_cache)=Uconst
554 iset_cache(ntwx_cache)=iset
557 qfrag_cache(i,ntwx_cache)=qfrag(i)
560 qpair_cache(i,ntwx_cache)=qpair(i)
563 utheta_cache(i,ntwx_cache)=utheta(i)
564 ugamma_cache(i,ntwx_cache)=ugamma(i)
565 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
570 c_cache(j,i,ntwx_cache)=c(j,i)
575 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
576 & .and..not.restart1file) then
579 open(irest1,file=mremd_rst_name,status='unknown')
580 write (irest1,*) "i2rep"
581 write (irest1,*) (i2rep(i),i=0,nodes-1)
582 write (irest1,*) "ifirst"
583 write (irest1,*) (ifirst(i),i=1,remd_m(1))
585 write (irest1,*) "nupa",il
586 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
587 write (irest1,*) "ndowna",il
588 write (irest1,*) ndowna(0,il),
589 & (ndowna(i,il),i=1,ndowna(0,il))
591 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
592 write (irest1,*) "nset"
593 write (irest1,*) nset
594 write (irest1,*) "mset"
595 write (irest1,*) (mset(i),i=1,nset)
596 write (irest1,*) "i2set"
597 write (irest1,*) (i2set(i),i=0,nodes-1)
598 write (irest1,*) "i_index"
602 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
610 open(irest2,file=rest2name,status='unknown')
611 write(irest2,*) totT,EK,potE,totE,t_bath
613 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
616 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
618 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
619 write (irest2,*) iset
626 c forced synchronization
627 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
628 & .and. .not. mremdsync) then
630 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
632 call mpi_recv(itime_master, 1, MPI_INTEGER,
633 & 0,101,CG_COMM, status, ierr)
634 call mpi_barrier(CG_COMM, ierr)
635 cdeb if (out1file.or.traj1file) then
636 cdeb call mpi_gather(itime,1,mpi_integer,
637 cdeb & icache_all,1,mpi_integer,king,
640 & call mpi_gather(ntwx_cache,1,mpi_integer,
641 & icache_all,1,mpi_integer,king,
644 & write(iout,*) 'REMD synchro at',itime_master,itime
645 if (itime_master.ge.n_timestep .or. ovrtim())
647 ctime call flush(iout)
652 if ((mod(itime,nstex).eq.0.and.me.eq.king
653 & .or.end_of_run.and.me.eq.king )
654 & .and. .not. mremdsync ) then
658 call mpi_isend(itime,1,MPI_INTEGER,i,101,
659 & CG_COMM, ireqi(i), ierr)
660 cd write(iout,*) 'REMD synchro with',i
663 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
664 call mpi_barrier(CG_COMM, ierr)
666 write(iout,*) 'REMD synchro at',itime,'time=',time01-time01_
667 if (out1file.or.traj1file) then
668 cdeb call mpi_gather(itime,1,mpi_integer,
669 cdeb & itime_all,1,mpi_integer,king,
671 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
672 cdeb & (itime_all(i),i=1,nodes)
674 cdeb imin_itime=itime_all(1)
676 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
678 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
679 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
680 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
681 call mpi_gather(ntwx_cache,1,mpi_integer,
682 & icache_all,1,mpi_integer,king,
684 c write(iout,'(a19,8000i8)') ' ntwx_cache',
685 c & (icache_all(i),i=1,nodes)
686 ii_write=icache_all(1)
688 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
690 c write(iout,*) "MIN ii_write=",ii_write
693 ctime call flush(iout)
695 if(mremdsync .and. mod(itime,nstex).eq.0) then
697 if (me.eq.king .or. .not. out1file)
698 & write(iout,*) 'REMD synchro at',itime
701 call mpi_gather(ntwx_cache,1,mpi_integer,
702 & icache_all,1,mpi_integer,king,
705 write(iout,'(a19,8000i8)') ' ntwx_cache',
706 & (icache_all(i),i=1,nodes)
707 ii_write=icache_all(1)
709 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
711 write(iout,*) "MIN ii_write=",ii_write
714 ctest call flush(iout)
717 c Update the time safety limiy
718 if (time001-time00.gt.safety) then
719 safety=time001-time00+600
720 if (me.eq.king .or. .not. out1file)
721 & write (iout,*) "****** SAFETY increased to",safety," s"
723 if (ovrtim() .and. me.eq.king) end_of_run=.true.
724 call MPI_Bcast(end_of_run,1,MPI_LOGICAL,king,CG_COMM,IERR)
726 if(synflag.and..not.end_of_run) then
730 c write(iout,*) 'REMD before',me,t_bath
732 c call mpi_gather(t_bath,1,mpi_double_precision,
733 c & remd_t_bath,1,mpi_double_precision,king,
735 potEcomp(n_ene+1)=t_bath
737 if (usampl.or.homol_nset.gt.1) then
738 potEcomp(n_ene+2)=iset
739 c write(iout,*) potEcomp(n_ene+1),potEcomp(n_ene+2),iset,nset
740 if (iset.lt.nset) then
743 if (homol_nset.gt.1) then
744 c broadcast iset to slaves and reduce energy
745 if (nfgtasks.gt.1) then
746 call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
747 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
748 call e_modeller(e_tmp)
749 c write(iout,*) "iset+1 before reduce",e_tmp
750 call MPI_Barrier(FG_COMM,IERR)
751 call MPI_Reduce(e_tmp,potEcomp(n_ene+3),1,
752 & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
754 call e_modeller(potEcomp(n_ene+3))
756 c write(iout,*) "iset+1",potEcomp(n_ene+3)
759 potEcomp(n_ene+3)=Uconst
762 c broadcast iset to slaves
763 if (nfgtasks.gt.1) then
764 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
765 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
768 potEcomp(n_ene+3)=0.0
773 if (homol_nset.gt.1) then
774 c broadcast iset to slaves and reduce energy
775 if (nfgtasks.gt.1) then
776 call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
777 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
778 call e_modeller(e_tmp)
779 c write(iout,*) "iset-1 before reduce",e_tmp
780 call MPI_Barrier(FG_COMM,IERR)
781 call MPI_Reduce(e_tmp,potEcomp(n_ene+4),1,
782 & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
784 call e_modeller(potEcomp(n_ene+4))
786 c write(iout,*) "iset-1",potEcomp(n_ene+4)
789 potEcomp(n_ene+4)=Uconst
792 c broadcast iset to slaves
793 if (nfgtasks.gt.1) then
794 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
795 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
798 potEcomp(n_ene+4)=0.0
801 if(hremd.gt.0) potEcomp(n_ene+2)=iset
802 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
803 & remd_ene(0,1),n_ene+5,mpi_double_precision,king,
806 call mpi_gather(elow,1,mpi_double_precision,
807 & elowi,1,mpi_double_precision,king,
809 call mpi_gather(ehigh,1,mpi_double_precision,
810 & ehighi,1,mpi_double_precision,king,
815 if (me.eq.king .or. .not. out1file) then
816 write(iout,*) 'REMD gather times=',time03-time01
820 if (restart1file) call write1rst(i_index)
823 if (me.eq.king .or. .not. out1file) then
824 write(iout,*) 'REMD writing rst time=',time04-time03
827 if (traj1file) call write1traj
829 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
830 cdeb & icache_all,1,mpi_integer,king,
832 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
833 cdeb & (icache_all(i),i=1,nodes)
838 if (me.eq.king .or. .not. out1file) then
839 write(iout,*) 'REMD writing traj time=',time05-time04
840 ctime call flush(iout)
845 if(homol_nset.gt.1) write(iout,*)
846 & 'energy_c temperature iset energy_c(iset+1) energy_c(iset-1)'
848 remd_t_bath(i)=remd_ene(n_ene+1,i)
849 iremd_iset(i)=remd_ene(n_ene+2,i)
851 & write(iout,'(i4,f10.3,f6.0,i3,2f10.3)')
852 & i,remd_ene(i_econstr,i),
853 & remd_ene(n_ene+1,i),iremd_iset(i),
854 & remd_ene(n_ene+3,i),remd_ene(n_ene+4,i)
858 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
860 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
864 write(iout,*) 'REMD exchange temp,ene'
866 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
867 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
871 c-------------------------------------
872 IF(.not.usampl.and.hremd.eq.0.and.homol_nset.le.1) THEN
874 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
876 ctime call flush(iout)
877 write (iout,*) "remd_m(1)",remd_m(1)
880 i=ifirst(iran_num(1,remd_m(1)))
884 ctime call flush(iout)
889 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
891 if(i.gt.0.and.nupa(0,i).gt.0) then
893 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
895 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
897 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
899 c do while (iex.eq.i)
900 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
901 iex=nupa(iran_num(1,int(nupa(0,i))),i)
903 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
905 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
907 c Swap temperatures between conformations i and iex with recalculating the free energies
908 c following temperature changes.
909 ene_iex_iex=remd_ene(0,iex)
910 ene_i_i=remd_ene(0,i)
911 c write (iout,*) "i",i," ene_i_i",ene_i_i,
912 c & " iex",iex," ene_iex_iex",ene_iex_iex
913 c write (iout,*) "rescaling weights with temperature",
916 call rescale_weights(remd_t_bath(i))
918 c write (iout,*) "0,iex",remd_t_bath(i)
919 c call enerprint(remd_ene(0,iex))
921 call sum_energy(remd_ene(0,iex),.false.)
922 ene_iex_i=remd_ene(0,iex)
923 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
925 c write (iout,*) "0,i",remd_t_bath(i)
926 c call enerprint(remd_ene(0,i))
928 call sum_energy(remd_ene(0,i),.false.)
929 c write (iout,*) "ene_i_i",remd_ene(0,i)
931 c write (iout,*) "rescaling weights with temperature",
933 if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
934 write (iout,*) "ERROR: inconsistent energies:",i,
935 & ene_i_i,remd_ene(0,i)
937 call rescale_weights(remd_t_bath(iex))
939 c write (iout,*) "0,i",remd_t_bath(iex)
940 c call enerprint(remd_ene(0,i))
942 call sum_energy(remd_ene(0,i),.false.)
943 c write (iout,*) "ene_i_iex",remd_ene(0,i)
945 ene_i_iex=remd_ene(0,i)
947 c write (iout,*) "0,iex",remd_t_bath(iex)
948 c call enerprint(remd_ene(0,iex))
950 call sum_energy(remd_ene(0,iex),.false.)
951 if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
952 write (iout,*) "ERROR: inconsistent energies:",iex,
953 & ene_iex_iex,remd_ene(0,iex)
955 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
956 c write (iout,*) "i",i," iex",iex
957 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
958 c & " ene_i_iex",ene_i_iex,
959 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
961 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
962 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
964 c write(iout,*) 'delta',delta
965 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
966 c & (remd_ene(i)-remd_ene(iex))/Rb/
967 c & (remd_t_bath(i)*remd_t_bath(iex))
969 if (delta .gt. 50.0d0) then
975 else if (delta.lt.-50.0d0) then
984 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
985 xxx=ran_number(0.0d0,1.0d0)
986 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
988 if (delta .gt. xxx) then
990 remd_t_bath(i)=remd_t_bath(iex)
992 remd_ene(0,i)=ene_i_iex
993 remd_ene(0,iex)=ene_iex_i
999 ehighi(i)=ehighi(iex)
1006 nupa(k,i)=nupa(k,iex)
1009 ndowna(k,i)=ndowna(k,iex)
1013 if (ifirst(il).eq.i) ifirst(il)=iex
1015 if (nupa(k,il).eq.i) then
1017 elseif (nupa(k,il).eq.iex) then
1022 if (ndowna(k,il).eq.i) then
1024 elseif (ndowna(k,il).eq.iex) then
1030 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1032 i2rep(i-1)=i2rep(iex-1)
1035 c write(iout,*) 'exchange',i,iex
1036 c write (iout,'(a8,100i4)') "@ ifirst",
1037 c & (ifirst(k),k=1,remd_m(1))
1039 c write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1040 c & (nupa(k,il),k=1,nupa(0,il))
1041 c write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1042 c & (ndowna(k,il),k=1,ndowna(0,il))
1047 remd_ene(0,iex)=ene_iex_iex
1048 remd_ene(0,i)=ene_i_i
1054 cd write (iout,*) "exchange completed"
1056 ELSEIF (usampl.or.homol_nset.gt.1) THEN
1058 c write(iout,*) "########",ii
1060 i_temp=iran_num(1,nrep)
1061 i_mult=iran_num(1,remd_m(i_temp))
1062 i_iset=iran_num(1,nset)
1063 i_mset=iran_num(1,mset(i_iset))
1064 i=i_index(i_temp,i_mult,i_iset,i_mset)
1066 c write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1068 if(t_exchange_only)then
1073 c write(iout,*) "i_dir=",i_dir
1075 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1078 i_mult1=iran_num(1,remd_m(i_temp1))
1080 i_mset1=iran_num(1,mset(i_iset1))
1081 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1083 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1086 i_mult1=iran_num(1,remd_m(i_temp1))
1088 i_mset1=iran_num(1,mset(i_iset1))
1089 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1091 econstr_temp_i=remd_ene(i_econstr,i)
1092 econstr_temp_iex=remd_ene(i_econstr,iex)
1093 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1094 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1096 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1099 i_mult1=iran_num(1,remd_m(i_temp1))
1101 i_mset1=iran_num(1,mset(i_iset1))
1102 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1103 econstr_temp_i=remd_ene(i_econstr,i)
1104 econstr_temp_iex=remd_ene(i_econstr,iex)
1105 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1106 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1112 c write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1113 ctime call flush(iout)
1115 c Swap temperatures between conformations i and iex with recalculating the free energies
1116 c following temperature changes.
1117 ene_iex_iex=remd_ene(0,iex)
1118 ene_i_i=remd_ene(0,i)
1119 co write (iout,*) "rescaling weights with temperature",
1121 call rescale_weights(remd_t_bath(i))
1123 call sum_energy(remd_ene(0,iex),.false.)
1124 ene_iex_i=remd_ene(0,iex)
1126 c ERROR only makes sense for dir =1
1127 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
1128 c call sum_energy(remd_ene(0,i),.false.)
1129 c write (iout,*) "ene_i_i",remd_ene(0,i)
1130 c write (iout,*) "rescaling weights with temperature",
1131 c & remd_t_bath(iex)
1132 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1133 c write (iout,*) "ERROR: inconsistent energies i:",i,
1134 c & ene_i_i,remd_ene(0,i)
1137 call rescale_weights(remd_t_bath(iex))
1138 call sum_energy(remd_ene(0,i),.false.)
1139 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1140 ene_i_iex=remd_ene(0,i)
1142 c ERROR only makes sense for dir =1
1143 c call sum_energy(remd_ene(0,iex),.false.)
1144 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1145 c write (iout,*) "ERROR: inconsistent energies iex:",iex,
1146 c & ene_iex_iex,remd_ene(0,iex)
1148 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1149 c write (iout,*) "i",i," iex",iex
1150 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1151 c & " ene_i_iex",ene_i_iex,
1152 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1154 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1155 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1157 c write(iout,*) 'delta',delta
1158 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1159 c & (remd_ene(i)-remd_ene(iex))/Rb/
1160 c & (remd_t_bath(i)*remd_t_bath(iex))
1161 if (delta .gt. 50.0d0) then
1166 if (i_dir.eq.1.or.i_dir.eq.3)
1167 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1168 if (i_dir.eq.2.or.i_dir.eq.3)
1169 & iremd_tot_usa(int(i2set(i-1)))=
1170 & iremd_tot_usa(int(i2set(i-1)))+1
1171 xxx=ran_number(0.0d0,1.0d0)
1172 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1173 if (delta .gt. xxx) then
1175 remd_t_bath(i)=remd_t_bath(iex)
1176 remd_t_bath(iex)=tmp
1179 iremd_iset(i)=iremd_iset(iex)
1180 iremd_iset(iex)=itmp
1182 remd_ene(0,i)=ene_i_iex
1183 remd_ene(0,iex)=ene_iex_i
1185 if (i_dir.eq.1.or.i_dir.eq.3)
1186 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1189 i2rep(i-1)=i2rep(iex-1)
1192 if (i_dir.eq.2.or.i_dir.eq.3)
1193 & iremd_acc_usa(int(i2set(i-1)))=
1194 & iremd_acc_usa(int(i2set(i-1)))+1
1197 i2set(i-1)=i2set(iex-1)
1200 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1201 i_index(i_temp,i_mult,i_iset,i_mset)=
1202 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1203 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1206 remd_ene(0,iex)=ene_iex_iex
1207 remd_ene(0,i)=ene_i_i
1208 remd_ene(i_econstr,iex)=econstr_temp_iex
1209 remd_ene(i_econstr,i)=econstr_temp_i
1213 cd do il1=1,mset(il)
1216 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1226 ELSEIF (hremd.gt.0) THEN
1228 cd write(iout,*) "########",ii
1230 i_temp=iran_num(1,nrep)
1231 i_mult=iran_num(1,remd_m(i_temp))
1232 i_iset=iran_num(1,nset)
1234 i=i_index(i_temp,i_mult,i_iset,i_mset)
1236 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1238 if(t_exchange_only)then
1244 cd write(iout,*) "i_dir=",i_dir
1246 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1249 i_mult1=iran_num(1,remd_m(i_temp1))
1252 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1254 elseif(i_dir.eq.2)then
1257 i_mult1=iran_num(1,remd_m(i_temp1))
1258 i_iset1=iran_num(1,hremd)
1259 do while(i_iset1.eq.i_iset)
1260 i_iset1=iran_num(1,hremd)
1263 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1265 elseif(remd_m(i_temp+1).gt.0)then
1268 i_mult1=iran_num(1,remd_m(i_temp1))
1269 i_iset1=iran_num(1,hremd)
1270 do while(i_iset1.eq.i_iset)
1271 i_iset1=iran_num(1,hremd)
1274 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1280 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1281 ctime call flush(iout)
1283 c Swap temperatures between conformations i and iex with recalculating the free energies
1284 c following temperature changes.
1285 ene_iex_iex=remd_ene(0,iex)
1286 ene_i_i=remd_ene(0,i)
1288 call set_hweights(i_iset)
1289 call rescale_weights(remd_t_bath(i))
1290 call sum_energy(remd_ene(0,iex),.false.)
1291 ene_iex_i=remd_ene(0,iex)
1293 call set_hweights(i_iset1)
1294 call rescale_weights(remd_t_bath(iex))
1295 call sum_energy(remd_ene(0,i),.false.)
1296 ene_i_iex=remd_ene(0,i)
1298 cd write(iout,*) ene_iex_iex,ene_i_i,ene_iex_i,ene_i_iex
1300 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1301 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1304 if (delta .gt. 50.0d0) then
1310 if (i_dir.eq.1.or.i_dir.eq.3)
1311 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1312 if (i_dir.eq.2.or.i_dir.eq.3)
1313 & iremd_tot_usa(int(i2set(i-1)))=
1314 & iremd_tot_usa(int(i2set(i-1)))+1
1315 xxx=ran_number(0.0d0,1.0d0)
1316 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1317 if (delta .gt. xxx) then
1319 cd write (iout,*) "exchange"
1321 remd_t_bath(i)=remd_t_bath(iex)
1322 remd_t_bath(iex)=tmp
1325 iremd_iset(i)=iremd_iset(iex)
1326 iremd_iset(iex)=itmp
1328 remd_ene(0,i)=ene_i_iex
1329 remd_ene(0,iex)=ene_iex_i
1331 if (i_dir.eq.1.or.i_dir.eq.3)
1332 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1335 i2rep(i-1)=i2rep(iex-1)
1338 if (i_dir.eq.2.or.i_dir.eq.3)
1339 & iremd_acc_usa(int(i2set(i-1)))=
1340 & iremd_acc_usa(int(i2set(i-1)))+1
1343 i2set(i-1)=i2set(iex-1)
1346 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1347 i_index(i_temp,i_mult,i_iset,i_mset)=
1348 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1349 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1352 cd do il1=1,mset(il)
1355 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1362 remd_ene(0,iex)=ene_iex_iex
1363 remd_ene(0,i)=ene_i_i
1374 c-------------------------------------
1375 write (iout,*) "NREP",nrep
1377 if(iremd_tot(i).ne.0)
1378 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1379 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1382 if(usampl.or.homol_nset.gt.1) then
1384 if(iremd_tot_usa(i).ne.0)
1385 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1386 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1392 if(iremd_tot_usa(i).ne.0)
1393 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_hremd',i,
1394 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1399 ctime call flush(iout)
1401 cd write (iout,'(a6,100i4)') "ifirst",
1402 cd & (ifirst(i),i=1,remd_m(1))
1404 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1405 cd & (nupa(i,il),i=1,nupa(0,il))
1406 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1407 cd & (ndowna(i,il),i=1,ndowna(0,il))
1412 cd write (iout,*) "Before scatter"
1414 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1415 & t_bath,1,mpi_double_precision,king,
1417 cd write (iout,*) "After scatter"
1419 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
1420 call mpi_scatter(iremd_iset,1,mpi_integer,
1421 & iset,1,mpi_integer,king,
1423 c 8/31/2015 Correction by AL: send new iset to slaves
1424 if (nfgtasks.gt.1) then
1425 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
1426 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
1432 if (me.eq.king .or. .not. out1file) then
1433 write(iout,*) 'REMD scatter time=',time07-time06
1437 call mpi_scatter(elowi,1,mpi_double_precision,
1438 & elow,1,mpi_double_precision,king,
1440 call mpi_scatter(ehighi,1,mpi_double_precision,
1441 & ehigh,1,mpi_double_precision,king,
1445 if(hremd.gt.0) call set_hweights(iset)
1446 call rescale_weights(t_bath)
1447 co write (iout,*) "Processor",me,
1448 co & " rescaling weights with temperature",t_bath
1450 stdfp=dsqrt(2*Rb*t_bath/d_time)
1452 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1456 stdforcp(i)=stdforcp(i)*sqrt(t_bath/t_bath_old)
1459 stdforcsc(i)=stdforcsc(i)*sqrt(t_bath/t_bath_old)
1462 cde write(iout,*) 'REMD after',me,t_bath
1464 if (me.eq.king .or. .not. out1file) then
1465 write(iout,*) 'REMD exchange time=',time08-time02
1466 ctime call flush(iout)
1471 if (restart1file) then
1472 if (me.eq.king .or. .not. out1file)
1473 & write(iout,*) 'writing restart at the end of run'
1474 call write1rst(i_index)
1477 if (traj1file) call write1traj
1479 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1480 cdeb & icache_all,1,mpi_integer,king,
1481 cdeb & CG_COMM,ierr)
1482 cdeb write(iout,'(a40,8000i8)')
1483 cdeb & ' ntwx_cache after traj1file at the end',
1484 cdeb & (icache_all(i),i=1,nodes)
1489 t_MD=MPI_Wtime()-tt0
1493 if (me.eq.king .or. .not. out1file) then
1494 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1496 & 'MD calculations setup:',t_MDsetup,
1497 & 'Energy & gradient evaluation:',t_enegrad,
1498 & 'Stochastic MD setup:',t_langsetup,
1499 & 'Stochastic MD step setup:',t_sdsetup,
1501 write (iout,'(/28(1h=),a25,27(1h=))')
1502 & ' End of MD calculation '
1503 if(hmc.gt.0) write (iout,*) 'HMC acceptance ratio',
1504 & n_timestep*1.0d0/hmc/hmc_acc
1509 c-----------------------------------------------------------------------
1510 subroutine write1rst(i_index)
1511 implicit real*8 (a-h,o-z)
1512 include 'DIMENSIONS'
1515 include 'COMMON.IOUNITS'
1516 include 'COMMON.REMD'
1517 include 'COMMON.SETUP'
1518 include 'COMMON.CHAIN'
1519 include 'COMMON.SBRIDGE'
1520 include 'COMMON.INTERACT'
1521 include 'COMMON.CONTROL'
1523 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1524 & d_restart2(3,2*maxres*maxprocs)
1528 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1529 common /przechowalnia/ d_restart1,d_restart2
1534 t5_restart1(4)=t_bath
1535 t5_restart1(5)=Uconst
1537 call mpi_gather(t5_restart1,5,mpi_real,
1538 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1546 call mpi_gather(r_d,3*2*nres,mpi_real,
1547 & d_restart1,3*2*nres,mpi_real,king,
1556 call mpi_gather(r_d,3*2*nres,mpi_real,
1557 & d_restart2,3*2*nres,mpi_real,king,
1562 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1564 call xdrfint_(ixdrf, i2rep(i), iret)
1567 call xdrfint_(ixdrf, ifirst(i), iret)
1571 call xdrfint_(ixdrf, nupa(i,il), iret)
1575 call xdrfint_(ixdrf, ndowna(i,il), iret)
1581 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1588 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1595 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1600 if(usampl.or.homol_nset.gt.1) then
1601 call xdrfint_(ixdrf, nset, iret)
1603 call xdrfint_(ixdrf,mset(i), iret)
1606 call xdrfint_(ixdrf,i2set(i), iret)
1612 itmp=i_index(i,j,il,il1)
1613 call xdrfint_(ixdrf,itmp, iret)
1620 call xdrfclose_(ixdrf, iret)
1622 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1624 call xdrfint(ixdrf, i2rep(i), iret)
1627 call xdrfint(ixdrf, ifirst(i), iret)
1631 call xdrfint(ixdrf, nupa(i,il), iret)
1635 call xdrfint(ixdrf, ndowna(i,il), iret)
1641 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1648 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1655 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1661 if(usampl.or.homol_nset.gt.1) then
1662 call xdrfint(ixdrf, nset, iret)
1664 call xdrfint(ixdrf,mset(i), iret)
1667 call xdrfint(ixdrf,i2set(i), iret)
1673 itmp=i_index(i,j,il,il1)
1674 call xdrfint(ixdrf,itmp, iret)
1681 call xdrfclose(ixdrf, iret)
1688 subroutine write1traj
1689 implicit real*8 (a-h,o-z)
1690 include 'DIMENSIONS'
1693 include 'COMMON.IOUNITS'
1694 include 'COMMON.REMD'
1695 include 'COMMON.SETUP'
1696 include 'COMMON.CHAIN'
1697 include 'COMMON.SBRIDGE'
1698 include 'COMMON.INTERACT'
1702 real xcoord(3,maxres2+2),prec
1703 real r_qfrag(50),r_qpair(100)
1704 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1705 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1706 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1707 & p_uscdiff(100*maxprocs)
1708 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1709 common /przechowalnia/ p_c
1711 call mpi_bcast(ii_write,1,mpi_integer,
1712 & king,CG_COMM,ierr)
1715 print *,'traj1file',me,ii_write,ntwx_cache
1719 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1721 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1724 t5_restart1(1)=totT_cache(ii)
1725 t5_restart1(2)=EK_cache(ii)
1726 t5_restart1(3)=potE_cache(ii)
1727 t5_restart1(4)=t_bath_cache(ii)
1728 t5_restart1(5)=Uconst_cache(ii)
1729 call mpi_gather(t5_restart1,5,mpi_real,
1730 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1732 call mpi_gather(iset_cache(ii),1,mpi_integer,
1733 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1736 r_qfrag(i)=qfrag_cache(i,ii)
1739 r_qpair(i)=qpair_cache(i,ii)
1742 r_utheta(i)=utheta_cache(i,ii)
1743 r_ugamma(i)=ugamma_cache(i,ii)
1744 r_uscdiff(i)=uscdiff_cache(i,ii)
1747 call mpi_gather(r_qfrag,nfrag,mpi_real,
1748 & p_qfrag,nfrag,mpi_real,king,
1750 call mpi_gather(r_qpair,npair,mpi_real,
1751 & p_qpair,npair,mpi_real,king,
1753 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1754 & p_utheta,nfrag_back,mpi_real,king,
1756 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1757 & p_ugamma,nfrag_back,mpi_real,king,
1759 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1760 & p_uscdiff,nfrag_back,mpi_real,king,
1764 write (iout,*) "p_qfrag"
1766 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1768 write (iout,*) "p_qpair"
1770 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1772 ctime call flush(iout)
1776 r_c(j,i)=c_cache(j,i,ii)
1780 call mpi_gather(r_c,3*2*nres,mpi_real,
1781 & p_c,3*2*nres,mpi_real,king,
1787 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1788 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1789 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1790 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1791 call xdrfint_(ixdrf, nss, iret)
1794 call xdrfint_(ixdrf, idssb(j)+nres, iret)
1795 call xdrfint_(ixdrf, jdssb(j)+nres, iret)
1797 call xdrfint_(ixdrf, ihpb(j), iret)
1798 call xdrfint_(ixdrf, jhpb(j), iret)
1801 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1802 call xdrfint_(ixdrf, iset_restart1(il), iret)
1804 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1807 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1810 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1811 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1812 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1817 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1822 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1826 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1830 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1831 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1832 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1833 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1834 call xdrfint(ixdrf, nss, iret)
1837 call xdrfint(ixdrf, idssb(j)+nres, iret)
1838 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1840 call xdrfint(ixdrf, ihpb(j), iret)
1841 call xdrfint(ixdrf, jhpb(j), iret)
1844 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1845 call xdrfint(ixdrf, iset_restart1(il), iret)
1847 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1850 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1853 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1854 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1855 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1860 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1865 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1869 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1875 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1877 if(me.eq.king) call xdrfclose(ixdrf, iret)
1879 do i=1,ntwx_cache-ii_write
1881 totT_cache(i)=totT_cache(ii_write+i)
1882 EK_cache(i)=EK_cache(ii_write+i)
1883 potE_cache(i)=potE_cache(ii_write+i)
1884 t_bath_cache(i)=t_bath_cache(ii_write+i)
1885 Uconst_cache(i)=Uconst_cache(ii_write+i)
1886 iset_cache(i)=iset_cache(ii_write+i)
1889 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1892 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1895 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1896 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1897 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1902 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1906 ntwx_cache=ntwx_cache-ii_write
1911 subroutine read1restart(i_index)
1912 implicit real*8 (a-h,o-z)
1913 include 'DIMENSIONS'
1916 include 'COMMON.IOUNITS'
1917 include 'COMMON.REMD'
1918 include 'COMMON.SETUP'
1919 include 'COMMON.CHAIN'
1920 include 'COMMON.SBRIDGE'
1921 include 'COMMON.INTERACT'
1922 include 'COMMON.CONTROL'
1923 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1926 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1927 common /przechowalnia/ d_restart1
1928 write (*,*) "Processor",me," called read1restart"
1931 open(irest2,file=mremd_rst_name,status='unknown')
1932 read(irest2,*,err=334) i
1933 write(iout,*) "Reading old rst in ASCI format"
1935 call read1restart_old
1939 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1942 call xdrfint_(ixdrf, i2rep(i), iret)
1945 call xdrfint_(ixdrf, ifirst(i), iret)
1948 call xdrfint_(ixdrf, nupa(0,il), iret)
1950 call xdrfint_(ixdrf, nupa(i,il), iret)
1953 call xdrfint_(ixdrf, ndowna(0,il), iret)
1955 call xdrfint_(ixdrf, ndowna(i,il), iret)
1960 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1964 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1967 call xdrfint(ixdrf, i2rep(i), iret)
1970 call xdrfint(ixdrf, ifirst(i), iret)
1973 call xdrfint(ixdrf, nupa(0,il), iret)
1975 call xdrfint(ixdrf, nupa(i,il), iret)
1978 call xdrfint(ixdrf, ndowna(0,il), iret)
1980 call xdrfint(ixdrf, ndowna(i,il), iret)
1985 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1990 call mpi_scatter(t_restart1,5,mpi_real,
1991 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1995 t_bath=t5_restart1(4)
2000 c read(irest2,'(3e15.5)')
2001 c & (d_restart1(j,i+2*nres*il),j=1,3)
2004 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
2006 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
2011 write (iout,*) "Conformation read",il
2013 write (iout,'(i5,3f10.5,5x,3f10.5)')
2014 & i,(d_restart1(j,i+2*nres*il),j=1,3),
2015 & (d_restart1(j,nres+i+2*nres*il),j=1,3)
2020 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2021 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2031 c read(irest2,'(3e15.5)')
2032 c & (d_restart1(j,i+2*nres*il),j=1,3)
2035 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
2037 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
2043 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2044 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2052 if(usampl.or.homol_nset.gt.1) then
2055 call xdrfint_(ixdrf, nset, iret)
2057 call xdrfint_(ixdrf,mset(i), iret)
2060 call xdrfint_(ixdrf,i2set(i), iret)
2066 call xdrfint_(ixdrf,itmp, iret)
2067 i_index(i,j,il,il1)=itmp
2075 call xdrfint(ixdrf, nset, iret)
2077 call xdrfint(ixdrf,mset(i), iret)
2080 call xdrfint(ixdrf,i2set(i), iret)
2086 call xdrfint(ixdrf,itmp, iret)
2087 i_index(i,j,il,il1)=itmp
2094 c Corrected AL 8/19/2014: each processor needs whole iset array not only its
2096 c call mpi_scatter(i2set,1,mpi_integer,
2097 c & iset,1,mpi_integer,king,
2099 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
2102 c broadcast iset to slaves
2103 if (nfgtasks.gt.1) then
2104 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
2105 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
2110 if(me.eq.king) close(irest2)
2114 subroutine read1restart_old
2115 implicit real*8 (a-h,o-z)
2116 include 'DIMENSIONS'
2119 include 'COMMON.IOUNITS'
2120 include 'COMMON.REMD'
2121 include 'COMMON.SETUP'
2122 include 'COMMON.CHAIN'
2123 include 'COMMON.SBRIDGE'
2124 include 'COMMON.INTERACT'
2125 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
2127 common /przechowalnia/ d_restart1
2129 open(irest2,file=mremd_rst_name,status='unknown')
2130 read (irest2,*) (i2rep(i),i=0,nodes-1)
2131 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2133 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2134 read (irest2,*) ndowna(0,il),
2135 & (ndowna(i,il),i=1,ndowna(0,il))
2138 read(irest2,*) (t_restart1(j,il),j=1,4)
2141 call mpi_scatter(t_restart1,5,mpi_real,
2142 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2146 t_bath=t5_restart1(4)
2151 read(irest2,'(3e15.5)')
2152 & (d_restart1(j,i+2*nres*il),j=1,3)
2156 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2157 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2167 read(irest2,'(3e15.5)')
2168 & (d_restart1(j,i+2*nres*il),j=1,3)
2172 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2173 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2179 if(me.eq.king) close(irest2)
2182 c-------------------------------------------------------------------
2183 subroutine set_hweights(iiset)
2184 implicit real*8 (a-h,o-z)
2186 include 'DIMENSIONS'
2187 include 'COMMON.FFIELD'
2188 include 'COMMON.REMD'
2191 weights(i)=hweights(iiset,i)