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()) end_of_run=.true.
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)
1419 & call mpi_scatter(iremd_iset,1,mpi_integer,
1420 & iset,1,mpi_integer,king,
1424 if (me.eq.king .or. .not. out1file) then
1425 write(iout,*) 'REMD scatter time=',time07-time06
1429 call mpi_scatter(elowi,1,mpi_double_precision,
1430 & elow,1,mpi_double_precision,king,
1432 call mpi_scatter(ehighi,1,mpi_double_precision,
1433 & ehigh,1,mpi_double_precision,king,
1437 if(hremd.gt.0) call set_hweights(iset)
1438 call rescale_weights(t_bath)
1439 co write (iout,*) "Processor",me,
1440 co & " rescaling weights with temperature",t_bath
1442 stdfp=dsqrt(2*Rb*t_bath/d_time)
1444 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1448 stdforcp(i)=stdforcp(i)*sqrt(t_bath/t_bath_old)
1451 stdforcsc(i)=stdforcsc(i)*sqrt(t_bath/t_bath_old)
1454 cde write(iout,*) 'REMD after',me,t_bath
1456 if (me.eq.king .or. .not. out1file) then
1457 write(iout,*) 'REMD exchange time=',time08-time02
1458 ctime call flush(iout)
1463 if (restart1file) then
1464 if (me.eq.king .or. .not. out1file)
1465 & write(iout,*) 'writing restart at the end of run'
1466 call write1rst(i_index)
1469 if (traj1file) call write1traj
1471 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1472 cdeb & icache_all,1,mpi_integer,king,
1473 cdeb & CG_COMM,ierr)
1474 cdeb write(iout,'(a40,8000i8)')
1475 cdeb & ' ntwx_cache after traj1file at the end',
1476 cdeb & (icache_all(i),i=1,nodes)
1481 t_MD=MPI_Wtime()-tt0
1485 if (me.eq.king .or. .not. out1file) then
1486 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1488 & 'MD calculations setup:',t_MDsetup,
1489 & 'Energy & gradient evaluation:',t_enegrad,
1490 & 'Stochastic MD setup:',t_langsetup,
1491 & 'Stochastic MD step setup:',t_sdsetup,
1493 write (iout,'(/28(1h=),a25,27(1h=))')
1494 & ' End of MD calculation '
1495 if(hmc.gt.0) write (iout,*) 'HMC acceptance ratio',
1496 & n_timestep*1.0d0/hmc/hmc_acc
1501 c-----------------------------------------------------------------------
1502 subroutine write1rst(i_index)
1503 implicit real*8 (a-h,o-z)
1504 include 'DIMENSIONS'
1507 include 'COMMON.IOUNITS'
1508 include 'COMMON.REMD'
1509 include 'COMMON.SETUP'
1510 include 'COMMON.CHAIN'
1511 include 'COMMON.SBRIDGE'
1512 include 'COMMON.INTERACT'
1513 include 'COMMON.CONTROL'
1515 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1516 & d_restart2(3,2*maxres*maxprocs)
1520 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1521 common /przechowalnia/ d_restart1,d_restart2
1526 t5_restart1(4)=t_bath
1527 t5_restart1(5)=Uconst
1529 call mpi_gather(t5_restart1,5,mpi_real,
1530 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1538 call mpi_gather(r_d,3*2*nres,mpi_real,
1539 & d_restart1,3*2*nres,mpi_real,king,
1548 call mpi_gather(r_d,3*2*nres,mpi_real,
1549 & d_restart2,3*2*nres,mpi_real,king,
1554 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1556 call xdrfint_(ixdrf, i2rep(i), iret)
1559 call xdrfint_(ixdrf, ifirst(i), iret)
1563 call xdrfint_(ixdrf, nupa(i,il), iret)
1567 call xdrfint_(ixdrf, ndowna(i,il), iret)
1573 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1580 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1587 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1592 if(usampl.or.homol_nset.gt.1) then
1593 call xdrfint_(ixdrf, nset, iret)
1595 call xdrfint_(ixdrf,mset(i), iret)
1598 call xdrfint_(ixdrf,i2set(i), iret)
1604 itmp=i_index(i,j,il,il1)
1605 call xdrfint_(ixdrf,itmp, iret)
1612 call xdrfclose_(ixdrf, iret)
1614 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1616 call xdrfint(ixdrf, i2rep(i), iret)
1619 call xdrfint(ixdrf, ifirst(i), iret)
1623 call xdrfint(ixdrf, nupa(i,il), iret)
1627 call xdrfint(ixdrf, ndowna(i,il), iret)
1633 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1640 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1647 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1653 if(usampl.or.homol_nset.gt.1) then
1654 call xdrfint(ixdrf, nset, iret)
1656 call xdrfint(ixdrf,mset(i), iret)
1659 call xdrfint(ixdrf,i2set(i), iret)
1665 itmp=i_index(i,j,il,il1)
1666 call xdrfint(ixdrf,itmp, iret)
1673 call xdrfclose(ixdrf, iret)
1680 subroutine write1traj
1681 implicit real*8 (a-h,o-z)
1682 include 'DIMENSIONS'
1685 include 'COMMON.IOUNITS'
1686 include 'COMMON.REMD'
1687 include 'COMMON.SETUP'
1688 include 'COMMON.CHAIN'
1689 include 'COMMON.SBRIDGE'
1690 include 'COMMON.INTERACT'
1694 real xcoord(3,maxres2+2),prec
1695 real r_qfrag(50),r_qpair(100)
1696 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1697 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1698 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1699 & p_uscdiff(100*maxprocs)
1700 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1701 common /przechowalnia/ p_c
1703 call mpi_bcast(ii_write,1,mpi_integer,
1704 & king,CG_COMM,ierr)
1707 print *,'traj1file',me,ii_write,ntwx_cache
1711 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1713 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1716 t5_restart1(1)=totT_cache(ii)
1717 t5_restart1(2)=EK_cache(ii)
1718 t5_restart1(3)=potE_cache(ii)
1719 t5_restart1(4)=t_bath_cache(ii)
1720 t5_restart1(5)=Uconst_cache(ii)
1721 call mpi_gather(t5_restart1,5,mpi_real,
1722 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1724 call mpi_gather(iset_cache(ii),1,mpi_integer,
1725 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1728 r_qfrag(i)=qfrag_cache(i,ii)
1731 r_qpair(i)=qpair_cache(i,ii)
1734 r_utheta(i)=utheta_cache(i,ii)
1735 r_ugamma(i)=ugamma_cache(i,ii)
1736 r_uscdiff(i)=uscdiff_cache(i,ii)
1739 call mpi_gather(r_qfrag,nfrag,mpi_real,
1740 & p_qfrag,nfrag,mpi_real,king,
1742 call mpi_gather(r_qpair,npair,mpi_real,
1743 & p_qpair,npair,mpi_real,king,
1745 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1746 & p_utheta,nfrag_back,mpi_real,king,
1748 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1749 & p_ugamma,nfrag_back,mpi_real,king,
1751 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1752 & p_uscdiff,nfrag_back,mpi_real,king,
1756 write (iout,*) "p_qfrag"
1758 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1760 write (iout,*) "p_qpair"
1762 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1764 ctime call flush(iout)
1768 r_c(j,i)=c_cache(j,i,ii)
1772 call mpi_gather(r_c,3*2*nres,mpi_real,
1773 & p_c,3*2*nres,mpi_real,king,
1779 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1780 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1781 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1782 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1783 call xdrfint_(ixdrf, nss, iret)
1786 call xdrfint_(ixdrf, idssb(j)+nres, iret)
1787 call xdrfint_(ixdrf, jdssb(j)+nres, iret)
1789 call xdrfint_(ixdrf, ihpb(j), iret)
1790 call xdrfint_(ixdrf, jhpb(j), iret)
1793 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1794 call xdrfint_(ixdrf, iset_restart1(il), iret)
1796 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1799 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1802 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1803 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1804 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1809 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1814 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1818 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1822 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1823 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1824 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1825 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1826 call xdrfint(ixdrf, nss, iret)
1829 call xdrfint(ixdrf, idssb(j)+nres, iret)
1830 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1832 call xdrfint(ixdrf, ihpb(j), iret)
1833 call xdrfint(ixdrf, jhpb(j), iret)
1836 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1837 call xdrfint(ixdrf, iset_restart1(il), iret)
1839 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1842 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1845 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1846 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1847 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1852 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1857 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1861 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1867 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1869 if(me.eq.king) call xdrfclose(ixdrf, iret)
1871 do i=1,ntwx_cache-ii_write
1873 totT_cache(i)=totT_cache(ii_write+i)
1874 EK_cache(i)=EK_cache(ii_write+i)
1875 potE_cache(i)=potE_cache(ii_write+i)
1876 t_bath_cache(i)=t_bath_cache(ii_write+i)
1877 Uconst_cache(i)=Uconst_cache(ii_write+i)
1878 iset_cache(i)=iset_cache(ii_write+i)
1881 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1884 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1887 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1888 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1889 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1894 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1898 ntwx_cache=ntwx_cache-ii_write
1903 subroutine read1restart(i_index)
1904 implicit real*8 (a-h,o-z)
1905 include 'DIMENSIONS'
1908 include 'COMMON.IOUNITS'
1909 include 'COMMON.REMD'
1910 include 'COMMON.SETUP'
1911 include 'COMMON.CHAIN'
1912 include 'COMMON.SBRIDGE'
1913 include 'COMMON.INTERACT'
1914 include 'COMMON.CONTROL'
1915 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1918 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1919 common /przechowalnia/ d_restart1
1920 write (*,*) "Processor",me," called read1restart"
1923 open(irest2,file=mremd_rst_name,status='unknown')
1924 read(irest2,*,err=334) i
1925 write(iout,*) "Reading old rst in ASCI format"
1927 call read1restart_old
1931 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1934 call xdrfint_(ixdrf, i2rep(i), iret)
1937 call xdrfint_(ixdrf, ifirst(i), iret)
1940 call xdrfint_(ixdrf, nupa(0,il), iret)
1942 call xdrfint_(ixdrf, nupa(i,il), iret)
1945 call xdrfint_(ixdrf, ndowna(0,il), iret)
1947 call xdrfint_(ixdrf, ndowna(i,il), iret)
1952 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1956 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1959 call xdrfint(ixdrf, i2rep(i), iret)
1962 call xdrfint(ixdrf, ifirst(i), iret)
1965 call xdrfint(ixdrf, nupa(0,il), iret)
1967 call xdrfint(ixdrf, nupa(i,il), iret)
1970 call xdrfint(ixdrf, ndowna(0,il), iret)
1972 call xdrfint(ixdrf, ndowna(i,il), iret)
1977 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1982 call mpi_scatter(t_restart1,5,mpi_real,
1983 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1987 t_bath=t5_restart1(4)
1992 c read(irest2,'(3e15.5)')
1993 c & (d_restart1(j,i+2*nres*il),j=1,3)
1996 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1998 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
2004 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2005 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2015 c read(irest2,'(3e15.5)')
2016 c & (d_restart1(j,i+2*nres*il),j=1,3)
2019 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
2021 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
2027 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2028 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2036 if(usampl.or.homol_nset.gt.1) then
2039 call xdrfint_(ixdrf, nset, iret)
2041 call xdrfint_(ixdrf,mset(i), iret)
2044 call xdrfint_(ixdrf,i2set(i), iret)
2050 call xdrfint_(ixdrf,itmp, iret)
2051 i_index(i,j,il,il1)=itmp
2059 call xdrfint(ixdrf, nset, iret)
2061 call xdrfint(ixdrf,mset(i), iret)
2064 call xdrfint(ixdrf,i2set(i), iret)
2070 call xdrfint(ixdrf,itmp, iret)
2071 i_index(i,j,il,il1)=itmp
2078 c Corrected AL 8/19/2014: each processor needs whole iset array not only its
2080 c call mpi_scatter(i2set,1,mpi_integer,
2081 c & iset,1,mpi_integer,king,
2083 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
2086 c broadcast iset to slaves
2087 if (nfgtasks.gt.1) then
2088 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
2089 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
2094 if(me.eq.king) close(irest2)
2098 subroutine read1restart_old
2099 implicit real*8 (a-h,o-z)
2100 include 'DIMENSIONS'
2103 include 'COMMON.IOUNITS'
2104 include 'COMMON.REMD'
2105 include 'COMMON.SETUP'
2106 include 'COMMON.CHAIN'
2107 include 'COMMON.SBRIDGE'
2108 include 'COMMON.INTERACT'
2109 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
2111 common /przechowalnia/ d_restart1
2113 open(irest2,file=mremd_rst_name,status='unknown')
2114 read (irest2,*) (i2rep(i),i=0,nodes-1)
2115 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2117 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2118 read (irest2,*) ndowna(0,il),
2119 & (ndowna(i,il),i=1,ndowna(0,il))
2122 read(irest2,*) (t_restart1(j,il),j=1,4)
2125 call mpi_scatter(t_restart1,5,mpi_real,
2126 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2130 t_bath=t5_restart1(4)
2135 read(irest2,'(3e15.5)')
2136 & (d_restart1(j,i+2*nres*il),j=1,3)
2140 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2141 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
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)
2163 if(me.eq.king) close(irest2)
2166 c-------------------------------------------------------------------
2167 subroutine set_hweights(iiset)
2168 implicit real*8 (a-h,o-z)
2170 include 'DIMENSIONS'
2171 include 'COMMON.FFIELD'
2172 include 'COMMON.REMD'
2175 weights(i)=hweights(iiset,i)