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'
25 include 'COMMON.HAIRPIN'
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
66 if(homol_nset.gt.1) then
74 if(usampl) i_econstr=20
79 do il1=1,max0(mset(il),1)
95 if(me.eq.king.or..not.out1file) then
96 write(iout,*) "i2rep",(i2rep(i),i=0,nodes-1)
97 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
98 write(iout,*) "i,j,il,il1,i_index(i,j,il,il1)"
103 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
110 c print *,'i2rep',me,i2rep(me)
111 c print *,'rep2i',(rep2i(i),i=0,nrep)
113 cold if (i2rep(me).eq.nrep) then
116 cold nup(0)=remd_m(i2rep(me)+1)
117 cold k=rep2i(int(i2rep(me)))+1
124 cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
126 cold if (i2rep(me).eq.1) then
129 cold ndown(0)=remd_m(i2rep(me)-1)
130 cold k=rep2i(i2rep(me)-2)+1
137 cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
140 write (*,*) "Processor",me," rest",rest,"
141 & restart1fie",restart1file
142 if(rest.and.restart1file) then
144 & inquire(file=mremd_rst_name,exist=file_exist)
145 cd write (*,*) me," Before broadcast: file_exist",file_exist
146 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
148 cd write (*,*) me," After broadcast: file_exist",file_exist
150 if(me.eq.king.or..not.out1file)
151 & write (iout,*) 'Master is reading restart1file'
152 call read1restart(i_index)
154 if(me.eq.king.or..not.out1file)
155 & write (iout,*) 'WARNING : no restart1file'
158 if(me.eq.king.or..not.out1file) then
159 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
160 write(iout,*) "i_index"
165 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
174 if (rest.and..not.restart1file)
175 & inquire(file=mremd_rst_name,exist=file_exist)
176 if(.not.file_exist.and.rest.and..not.restart1file)
177 & write(iout,*) 'WARNING : no restart file',mremd_rst_name
178 IF (rest.and.file_exist.and..not.restart1file) THEN
179 write (iout,*) 'Master is reading restart file',
181 open(irest2,file=mremd_rst_name,status='unknown')
183 read (irest2,*) (i2rep(i),i=0,nodes-1)
185 read (irest2,*) (ifirst(i),i=1,remd_m(1))
188 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
190 read (irest2,*) ndowna(0,il),
191 & (ndowna(i,il),i=1,ndowna(0,il))
193 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
197 read (irest2,*) (mset(i),i=1,nset)
199 read (irest2,*) (i2set(i),i=0,nodes-1)
204 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
209 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
210 write(iout,*) "i_index"
215 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
224 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
225 write (iout,'(a6,1000i5)') "ifirst",
226 & (ifirst(i),i=1,remd_m(1))
228 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
229 & (nupa(i,il),i=1,nupa(0,il))
230 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
231 & (ndowna(i,il),i=1,ndowna(0,il))
233 ELSE IF (.not.(rest.and.file_exist)) THEN
239 if (i2rep(il-1).eq.nrep) then
242 nupa(0,il)=remd_m(i2rep(il-1)+1)
243 k=rep2i(int(i2rep(il-1)))+1
249 if (i2rep(il-1).eq.1) then
252 ndowna(0,il)=remd_m(i2rep(il-1)-1)
253 k=rep2i(i2rep(il-1)-2)+1
261 write (iout,'(a6,100i4)') "ifirst",
262 & (ifirst(i),i=1,remd_m(1))
264 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
265 & (nupa(i,il),i=1,nupa(0,il))
266 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
267 & (ndowna(i,il),i=1,ndowna(0,il))
273 c t_bath=retmin+(retmax-retmin)*me/(nodes-1)
274 if(.not.(rest.and.file_exist.and.restart1file)) then
275 if (me .eq. king) then
278 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
280 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
281 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
284 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
286 c broadcast iset to slaves
287 if (nfgtasks.gt.1) then
288 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
289 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
291 if(me.eq.king.or..not.out1file)
292 & write(iout,*) me,"iset=",iset,"t_bath=",t_bath
295 stdfp=dsqrt(2*Rb*t_bath/d_time)
297 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
300 c print *,'irep',me,t_bath
302 if (me.eq.king .or. .not. out1file)
303 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
304 call rescale_weights(t_bath)
308 c------copy MD--------------
309 c The driver for molecular dynamics subroutines
310 c------------------------------------------------
316 if(me.eq.king.or..not.out1file)
317 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
323 c Determine the inverse of the inertia matrix.
324 call setup_MD_matrices
328 if (me.eq.king .or. .not. out1file)
329 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
330 stdfp=dsqrt(2*Rb*t_bath/d_time)
332 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
334 call rescale_weights(t_bath)
338 t_MDsetup = MPI_Wtime()-tt0
340 t_MDsetup = tcpu()-tt0
343 c Entering the MD loop
349 if (lang.eq.2 .or. lang.eq.3) then
353 call sd_verlet_p_setup
355 call sd_verlet_ciccotti_setup
359 pfric0_mat(i,j,0)=pfric_mat(i,j)
360 afric0_mat(i,j,0)=afric_mat(i,j)
361 vfric0_mat(i,j,0)=vfric_mat(i,j)
362 prand0_mat(i,j,0)=prand_mat(i,j)
363 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
364 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
369 flag_stoch(i)=.false.
373 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
375 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
379 else if (lang.eq.1 .or. lang.eq.4) then
383 if (me.eq.king .or. .not. out1file)
384 & write(iout,*) 'Setup time',time00-walltime
387 t_langsetup=MPI_Wtime()-tt0
390 t_langsetup=tcpu()-tt0
395 do while(.not.end_of_run)
397 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
398 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
400 if (lang.gt.0 .and. surfarea .and.
401 & mod(itime,reset_fricmat).eq.0) then
402 if (lang.eq.2 .or. lang.eq.3) then
406 call sd_verlet_p_setup
408 call sd_verlet_ciccotti_setup
412 pfric0_mat(i,j,0)=pfric_mat(i,j)
413 afric0_mat(i,j,0)=afric_mat(i,j)
414 vfric0_mat(i,j,0)=vfric_mat(i,j)
415 prand0_mat(i,j,0)=prand_mat(i,j)
416 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
417 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
422 flag_stoch(i)=.false.
425 else if (lang.eq.1 .or. lang.eq.4) then
428 write (iout,'(a,i10)')
429 & "Friction matrix reset based on surface area, itime",itime
431 if (reset_vel .and. tbf .and. lang.eq.0
432 & .and. mod(itime,count_reset_vel).eq.0) then
434 if (me.eq.king .or. .not. out1file)
435 & write(iout,'(a,f20.2)')
436 & "Velocities reset to random values, time",totT
439 d_t_old(j,i)=d_t(j,i)
443 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
447 d_t(j,0)=d_t(j,0)-vcm(j)
450 kinetic_T=2.0d0/(dimen3*Rb)*EK
451 scalfac=dsqrt(T_bath/kinetic_T)
452 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
455 d_t_old(j,i)=scalfac*d_t(j,i)
461 c Time-reversible RESPA algorithm
462 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
463 call RESPA_step(itime)
465 c Variable time step algorithm.
466 call velverlet_step(itime)
470 call brown_step(itime)
472 print *,"Brown dynamics not here!"
474 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
480 if (mod(itime,ntwe).eq.0) call statout(itime)
482 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
483 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
485 call hairpin(.true.,nharp,iharp)
486 call secondary2(.true.)
487 call pdbout(potE,tytul,ipdb)
492 if (mod(itime,ntwx).eq.0.and.traj1file) then
493 if(ntwx_cache.lt.max_cache_traj_use) then
494 ntwx_cache=ntwx_cache+1
496 if (max_cache_traj_use.ne.1)
497 & print *,itime,"processor ",me," over cache ",ntwx_cache
500 totT_cache(i)=totT_cache(i+1)
501 EK_cache(i)=EK_cache(i+1)
502 potE_cache(i)=potE_cache(i+1)
503 t_bath_cache(i)=t_bath_cache(i+1)
504 Uconst_cache(i)=Uconst_cache(i+1)
505 iset_cache(i)=iset_cache(i+1)
508 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
511 qpair_cache(ii,i)=qpair_cache(ii,i+1)
514 utheta_cache(ii,i)=utheta_cache(ii,i+1)
515 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
516 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
522 c_cache(j,ii,i)=c_cache(j,ii,i+1)
528 totT_cache(ntwx_cache)=totT
529 EK_cache(ntwx_cache)=EK
530 potE_cache(ntwx_cache)=potE
531 t_bath_cache(ntwx_cache)=t_bath
532 Uconst_cache(ntwx_cache)=Uconst
533 iset_cache(ntwx_cache)=iset
536 qfrag_cache(i,ntwx_cache)=qfrag(i)
539 qpair_cache(i,ntwx_cache)=qpair(i)
542 utheta_cache(i,ntwx_cache)=utheta(i)
543 ugamma_cache(i,ntwx_cache)=ugamma(i)
544 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
546 C print *,'przed returnbox'
548 C call enerprint(remd_ene(0,i))
551 c_cache(j,i,ntwx_cache)=c(j,i)
556 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
557 & .and..not.restart1file) then
560 open(irest1,file=mremd_rst_name,status='unknown')
561 write (irest1,*) "i2rep"
562 write (irest1,*) (i2rep(i),i=0,nodes-1)
563 write (irest1,*) "ifirst"
564 write (irest1,*) (ifirst(i),i=1,remd_m(1))
566 write (irest1,*) "nupa",il
567 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
568 write (irest1,*) "ndowna",il
569 write (irest1,*) ndowna(0,il),
570 & (ndowna(i,il),i=1,ndowna(0,il))
572 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
573 write (irest1,*) "nset"
574 write (irest1,*) nset
575 write (irest1,*) "mset"
576 write (irest1,*) (mset(i),i=1,nset)
577 write (irest1,*) "i2set"
578 write (irest1,*) (i2set(i),i=0,nodes-1)
579 write (irest1,*) "i_index"
583 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
591 open(irest2,file=rest2name,status='unknown')
592 write(irest2,*) totT,EK,potE,totE,t_bath
594 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
597 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
599 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
600 write (irest2,*) iset
607 c forced synchronization
608 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
609 & .and. .not. mremdsync) then
611 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
613 call mpi_recv(itime_master, 1, MPI_INTEGER,
614 & 0,101,CG_COMM, status, ierr)
615 call mpi_barrier(CG_COMM, ierr)
616 cdeb if (out1file.or.traj1file) then
617 cdeb call mpi_gather(itime,1,mpi_integer,
618 cdeb & icache_all,1,mpi_integer,king,
621 & call mpi_gather(ntwx_cache,1,mpi_integer,
622 & icache_all,1,mpi_integer,king,
625 & write(iout,*) 'REMD synchro at',itime_master,itime
626 if (itime_master.ge.n_timestep .or. ovrtim())
628 ctime call flush(iout)
633 if ((mod(itime,nstex).eq.0.and.me.eq.king
634 & .or.end_of_run.and.me.eq.king )
635 & .and. .not. mremdsync ) then
639 call mpi_isend(itime,1,MPI_INTEGER,i,101,
640 & CG_COMM, ireqi(i), ierr)
641 cd write(iout,*) 'REMD synchro with',i
644 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
645 call mpi_barrier(CG_COMM, ierr)
647 write(iout,*) 'REMD synchro at',itime,'time=',time01-time01_
648 if (out1file.or.traj1file) then
649 cdeb call mpi_gather(itime,1,mpi_integer,
650 cdeb & itime_all,1,mpi_integer,king,
652 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
653 cdeb & (itime_all(i),i=1,nodes)
655 cdeb imin_itime=itime_all(1)
657 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
659 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
660 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
661 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
662 call mpi_gather(ntwx_cache,1,mpi_integer,
663 & icache_all,1,mpi_integer,king,
665 c write(iout,'(a19,8000i8)') ' ntwx_cache',
666 c & (icache_all(i),i=1,nodes)
667 ii_write=icache_all(1)
669 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
671 c write(iout,*) "MIN ii_write=",ii_write
674 ctime call flush(iout)
676 if(mremdsync .and. mod(itime,nstex).eq.0) then
678 if (me.eq.king .or. .not. out1file)
679 & write(iout,*) 'REMD synchro at',itime
682 call mpi_gather(ntwx_cache,1,mpi_integer,
683 & icache_all,1,mpi_integer,king,
686 write(iout,'(a19,8000i8)') ' ntwx_cache',
687 & (icache_all(i),i=1,nodes)
688 ii_write=icache_all(1)
690 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
692 write(iout,*) "MIN ii_write=",ii_write
698 c Update the time safety limiy
699 if (time001-time00.gt.safety) then
700 safety=time001-time00+600
701 write (iout,*) "****** SAFETY increased to",safety," s"
703 if (ovrtim()) end_of_run=.true.
705 if(synflag.and..not.end_of_run) then
709 write(iout,*) 'REMD before',me,t_bath
711 c call mpi_gather(t_bath,1,mpi_double_precision,
712 c & remd_t_bath,1,mpi_double_precision,king,
714 potEcomp(n_ene+1)=t_bath
716 if (usampl.or.homol_nset.gt.1) then
717 potEcomp(n_ene+2)=iset
718 if (iset.lt.nset) then
721 if (homol_nset.gt.1) then
722 c broadcast iset to slaves and reduce energy
723 if (nfgtasks.gt.1) then
724 call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
725 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
726 call e_modeller(e_tmp)
727 c write(iout,*) "iset+1 before reduce",e_tmp
728 call MPI_Barrier(FG_COMM,IERR)
729 call MPI_Reduce(e_tmp,potEcomp(n_ene+3),1,
730 & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
732 call e_modeller(potEcomp(n_ene+3))
734 c write(iout,*) "iset+1",potEcomp(n_ene+3)
737 potEcomp(n_ene+3)=Uconst
740 c broadcast iset to slaves
741 if (nfgtasks.gt.1) then
742 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
743 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
746 potEcomp(n_ene+3)=0.0
751 if (homol_nset.gt.1) then
752 c broadcast iset to slaves and reduce energy
753 if (nfgtasks.gt.1) then
754 call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR)
755 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
756 call e_modeller(e_tmp)
757 c write(iout,*) "iset-1 before reduce",e_tmp
758 call MPI_Barrier(FG_COMM,IERR)
759 call MPI_Reduce(e_tmp,potEcomp(n_ene+4),1,
760 & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
762 call e_modeller(potEcomp(n_ene+4))
764 c write(iout,*) "iset-1",potEcomp(n_ene+4)
767 potEcomp(n_ene+4)=Uconst
770 c broadcast iset to slaves
771 if (nfgtasks.gt.1) then
772 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
773 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
776 potEcomp(n_ene+4)=0.0
779 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
780 & remd_ene(0,1),n_ene+5,mpi_double_precision,king,
783 call mpi_gather(elow,1,mpi_double_precision,
784 & elowi,1,mpi_double_precision,king,
786 call mpi_gather(ehigh,1,mpi_double_precision,
787 & ehighi,1,mpi_double_precision,king,
792 if (me.eq.king .or. .not. out1file) then
793 write(iout,*) 'REMD gather times=',time03-time01
797 if (restart1file) call write1rst(i_index)
800 if (me.eq.king .or. .not. out1file) then
801 write(iout,*) 'REMD writing rst time=',time04-time03
804 if (traj1file) call write1traj
806 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
807 cdeb & icache_all,1,mpi_integer,king,
809 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
810 cdeb & (icache_all(i),i=1,nodes)
815 if (me.eq.king .or. .not. out1file) then
816 write(iout,*) 'REMD writing traj time=',time05-time04
822 if(homol_nset.gt.1) write(iout,*)
823 & 'energy_c temperature iset energy_c(iset+1) energy_c(iset-1)'
825 remd_t_bath(i)=remd_ene(n_ene+1,i)
826 iremd_iset(i)=remd_ene(n_ene+2,i)
828 & write(iout,'(i4,f10.3,f6.0,i3,2f10.3)')
829 & i,remd_ene(i_econstr,i),
830 & remd_ene(n_ene+1,i),iremd_iset(i),
831 & remd_ene(n_ene+3,i),remd_ene(n_ene+4,i)
835 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
837 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
841 write(iout,*) 'REMD exchange temp,ene'
843 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
844 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
848 c-------------------------------------
849 IF(.not.usampl.and.hremd.eq.0.and.homol_nset.le.1) THEN
851 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
853 ctime call flush(iout)
854 write (iout,*) "remd_m(1)",remd_m(1)
857 i=ifirst(iran_num(1,remd_m(1)))
861 ctime call flush(iout)
866 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
868 if(i.gt.0.and.nupa(0,i).gt.0) then
870 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
872 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
874 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
876 c do while (iex.eq.i)
877 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
878 iex=nupa(iran_num(1,int(nupa(0,i))),i)
880 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
882 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
884 c Swap temperatures between conformations i and iex with recalculating the free energies
885 c following temperature changes.
886 ene_iex_iex=remd_ene(0,iex)
887 ene_i_i=remd_ene(0,i)
888 c write (iout,*) "i",i," ene_i_i",ene_i_i,
889 c & " iex",iex," ene_iex_iex",ene_iex_iex
890 c write (iout,*) "rescaling weights with temperature",
893 call rescale_weights(remd_t_bath(i))
895 c write (iout,*) "0,iex",remd_t_bath(i)
896 c call enerprint(remd_ene(0,iex))
898 call sum_energy(remd_ene(0,iex),.false.)
899 ene_iex_i=remd_ene(0,iex)
900 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
902 c write (iout,*) "0,i",remd_t_bath(i)
903 c call enerprint(remd_ene(0,i))
905 call sum_energy(remd_ene(0,i),.false.)
906 c write (iout,*) "ene_i_i",remd_ene(0,i)
908 c write (iout,*) "rescaling weights with temperature",
910 if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then
911 write (iout,*) "ERROR: inconsistent energies:",i,
912 & ene_i_i,remd_ene(0,i)
914 call rescale_weights(remd_t_bath(iex))
916 c write (iout,*) "0,i",remd_t_bath(iex)
917 c call enerprint(remd_ene(0,i))
919 call sum_energy(remd_ene(0,i),.false.)
920 c write (iout,*) "ene_i_iex",remd_ene(0,i)
922 ene_i_iex=remd_ene(0,i)
924 c write (iout,*) "0,iex",remd_t_bath(iex)
925 c call enerprint(remd_ene(0,iex))
927 call sum_energy(remd_ene(0,iex),.false.)
928 if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then
929 write (iout,*) "ERROR: inconsistent energies:",iex,
930 & ene_iex_iex,remd_ene(0,iex)
932 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
933 c write (iout,*) "i",i," iex",iex
934 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
935 c & " ene_i_iex",ene_i_iex,
936 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
938 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
939 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
941 c write(iout,*) 'delta',delta
942 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
943 c & (remd_ene(i)-remd_ene(iex))/Rb/
944 c & (remd_t_bath(i)*remd_t_bath(iex))
946 if (delta .gt. 50.0d0) then
952 else if (delta.lt.-50.0d0) then
961 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
962 xxx=ran_number(0.0d0,1.0d0)
963 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
965 if (delta .gt. xxx) then
967 remd_t_bath(i)=remd_t_bath(iex)
969 remd_ene(0,i)=ene_i_iex
970 remd_ene(0,iex)=ene_iex_i
976 ehighi(i)=ehighi(iex)
983 nupa(k,i)=nupa(k,iex)
986 ndowna(k,i)=ndowna(k,iex)
990 if (ifirst(il).eq.i) ifirst(il)=iex
992 if (nupa(k,il).eq.i) then
994 elseif (nupa(k,il).eq.iex) then
999 if (ndowna(k,il).eq.i) then
1001 elseif (ndowna(k,il).eq.iex) then
1007 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1009 i2rep(i-1)=i2rep(iex-1)
1012 c write(iout,*) 'exchange',i,iex
1013 c write (iout,'(a8,100i4)') "@ ifirst",
1014 c & (ifirst(k),k=1,remd_m(1))
1016 c write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1017 c & (nupa(k,il),k=1,nupa(0,il))
1018 c write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1019 c & (ndowna(k,il),k=1,ndowna(0,il))
1024 remd_ene(0,iex)=ene_iex_iex
1025 remd_ene(0,i)=ene_i_i
1031 cd write (iout,*) "exchange completed"
1033 ELSEIF (usampl.or.homol_nset.gt.1) THEN
1035 c write(iout,*) "########",ii
1037 i_temp=iran_num(1,nrep)
1038 i_mult=iran_num(1,remd_m(i_temp))
1039 i_iset=iran_num(1,nset)
1040 i_mset=iran_num(1,mset(i_iset))
1041 i=i_index(i_temp,i_mult,i_iset,i_mset)
1043 c write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1045 if(t_exchange_only)then
1050 c write(iout,*) "i_dir=",i_dir
1052 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1055 i_mult1=iran_num(1,remd_m(i_temp1))
1057 i_mset1=iran_num(1,mset(i_iset1))
1058 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1060 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1063 i_mult1=iran_num(1,remd_m(i_temp1))
1065 i_mset1=iran_num(1,mset(i_iset1))
1066 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1068 econstr_temp_i=remd_ene(i_econstr,i)
1069 econstr_temp_iex=remd_ene(i_econstr,iex)
1070 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1071 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1073 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1076 i_mult1=iran_num(1,remd_m(i_temp1))
1078 i_mset1=iran_num(1,mset(i_iset1))
1079 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1080 econstr_temp_i=remd_ene(i_econstr,i)
1081 econstr_temp_iex=remd_ene(i_econstr,iex)
1082 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1083 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1089 c write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1090 ctime call flush(iout)
1092 c Swap temperatures between conformations i and iex with recalculating the free energies
1093 c following temperature changes.
1094 ene_iex_iex=remd_ene(0,iex)
1095 ene_i_i=remd_ene(0,i)
1096 co write (iout,*) "rescaling weights with temperature",
1098 call rescale_weights(remd_t_bath(i))
1100 call sum_energy(remd_ene(0,iex),.false.)
1101 ene_iex_i=remd_ene(0,iex)
1103 c ERROR only makes sense for dir =1
1104 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
1105 c call sum_energy(remd_ene(0,i),.false.)
1106 c write (iout,*) "ene_i_i",remd_ene(0,i)
1107 c write (iout,*) "rescaling weights with temperature",
1108 c & remd_t_bath(iex)
1109 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1110 c write (iout,*) "ERROR: inconsistent energies i:",i,
1111 c & ene_i_i,remd_ene(0,i)
1114 call rescale_weights(remd_t_bath(iex))
1115 call sum_energy(remd_ene(0,i),.false.)
1116 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1117 ene_i_iex=remd_ene(0,i)
1119 c ERROR only makes sense for dir =1
1120 c call sum_energy(remd_ene(0,iex),.false.)
1121 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1122 c write (iout,*) "ERROR: inconsistent energies iex:",iex,
1123 c & ene_iex_iex,remd_ene(0,iex)
1125 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1126 c write (iout,*) "i",i," iex",iex
1127 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1128 c & " ene_i_iex",ene_i_iex,
1129 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1131 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1132 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1134 c write(iout,*) 'delta',delta
1135 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1136 c & (remd_ene(i)-remd_ene(iex))/Rb/
1137 c & (remd_t_bath(i)*remd_t_bath(iex))
1138 if (delta .gt. 50.0d0) then
1143 if (i_dir.eq.1.or.i_dir.eq.3)
1144 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1145 if (i_dir.eq.2.or.i_dir.eq.3)
1146 & iremd_tot_usa(int(i2set(i-1)))=
1147 & iremd_tot_usa(int(i2set(i-1)))+1
1148 xxx=ran_number(0.0d0,1.0d0)
1149 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1150 if (delta .gt. xxx) then
1152 remd_t_bath(i)=remd_t_bath(iex)
1153 remd_t_bath(iex)=tmp
1156 iremd_iset(i)=iremd_iset(iex)
1157 iremd_iset(iex)=itmp
1159 remd_ene(0,i)=ene_i_iex
1160 remd_ene(0,iex)=ene_iex_i
1162 if (i_dir.eq.1.or.i_dir.eq.3)
1163 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1166 i2rep(i-1)=i2rep(iex-1)
1169 if (i_dir.eq.2.or.i_dir.eq.3)
1170 & iremd_acc_usa(int(i2set(i-1)))=
1171 & iremd_acc_usa(int(i2set(i-1)))+1
1174 i2set(i-1)=i2set(iex-1)
1177 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1178 i_index(i_temp,i_mult,i_iset,i_mset)=
1179 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1180 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1183 remd_ene(0,iex)=ene_iex_iex
1184 remd_ene(0,i)=ene_i_i
1185 remd_ene(i_econstr,iex)=econstr_temp_iex
1186 remd_ene(i_econstr,i)=econstr_temp_i
1190 cd do il1=1,mset(il)
1193 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1206 c-------------------------------------
1207 write (iout,*) "NREP",nrep
1209 if(iremd_tot(i).ne.0)
1210 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1211 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1214 if(usampl.or.homol_nset.gt.1) then
1216 if(iremd_tot_usa(i).ne.0)
1217 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1218 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1224 cd write (iout,'(a6,100i4)') "ifirst",
1225 cd & (ifirst(i),i=1,remd_m(1))
1227 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1228 cd & (nupa(i,il),i=1,nupa(0,il))
1229 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1230 cd & (ndowna(i,il),i=1,ndowna(0,il))
1235 cd write (iout,*) "Before scatter"
1237 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1238 & t_bath,1,mpi_double_precision,king,
1240 cd write (iout,*) "After scatter"
1242 if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
1243 call mpi_scatter(iremd_iset,1,mpi_integer,
1244 & iset,1,mpi_integer,king,
1246 c 8/31/2015 Correction by AL: send new iset to slaves
1247 if (nfgtasks.gt.1) then
1248 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
1249 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
1255 if (me.eq.king .or. .not. out1file) then
1256 write(iout,*) 'REMD scatter time=',time07-time06
1260 call mpi_scatter(elowi,1,mpi_double_precision,
1261 & elow,1,mpi_double_precision,king,
1263 call mpi_scatter(ehighi,1,mpi_double_precision,
1264 & ehigh,1,mpi_double_precision,king,
1267 call rescale_weights(t_bath)
1268 co write (iout,*) "Processor",me,
1269 co & " rescaling weights with temperature",t_bath
1271 stdfp=dsqrt(2*Rb*t_bath/d_time)
1273 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1276 cde write(iout,*) 'REMD after',me,t_bath
1278 if (me.eq.king .or. .not. out1file) then
1279 write(iout,*) 'REMD exchange time=',time08-time02
1280 ctime call flush(iout)
1285 if (restart1file) then
1286 if (me.eq.king .or. .not. out1file)
1287 & write(iout,*) 'writing restart at the end of run'
1288 call write1rst(i_index)
1291 if (traj1file) call write1traj
1293 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1294 cdeb & icache_all,1,mpi_integer,king,
1295 cdeb & CG_COMM,ierr)
1296 cdeb write(iout,'(a40,8000i8)')
1297 cdeb & ' ntwx_cache after traj1file at the end',
1298 cdeb & (icache_all(i),i=1,nodes)
1303 t_MD=MPI_Wtime()-tt0
1307 if (me.eq.king .or. .not. out1file) then
1308 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1310 & 'MD calculations setup:',t_MDsetup,
1311 & 'Energy & gradient evaluation:',t_enegrad,
1312 & 'Stochastic MD setup:',t_langsetup,
1313 & 'Stochastic MD step setup:',t_sdsetup,
1315 write (iout,'(/28(1h=),a25,27(1h=))')
1316 & ' End of MD calculation '
1321 c-----------------------------------------------------------------------
1322 subroutine write1rst(i_index)
1323 implicit real*8 (a-h,o-z)
1324 include 'DIMENSIONS'
1327 include 'COMMON.IOUNITS'
1328 include 'COMMON.REMD'
1329 include 'COMMON.SETUP'
1330 include 'COMMON.CHAIN'
1331 include 'COMMON.SBRIDGE'
1332 include 'COMMON.INTERACT'
1333 include 'COMMON.CONTROL'
1335 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1336 & d_restart2(3,2*maxres*maxprocs)
1340 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1341 common /przechowalnia/ d_restart1,d_restart2
1346 t5_restart1(4)=t_bath
1347 t5_restart1(5)=Uconst
1349 call mpi_gather(t5_restart1,5,mpi_real,
1350 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1358 call mpi_gather(r_d,3*2*nres,mpi_real,
1359 & d_restart1,3*2*nres,mpi_real,king,
1368 call mpi_gather(r_d,3*2*nres,mpi_real,
1369 & d_restart2,3*2*nres,mpi_real,king,
1374 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1376 call xdrfint_(ixdrf, i2rep(i), iret)
1379 call xdrfint_(ixdrf, ifirst(i), iret)
1383 call xdrfint_(ixdrf, nupa(i,il), iret)
1387 call xdrfint_(ixdrf, ndowna(i,il), iret)
1393 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1400 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1407 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1412 if(usampl.or.homol_nset.gt.1) then
1413 call xdrfint_(ixdrf, nset, iret)
1415 call xdrfint_(ixdrf,mset(i), iret)
1418 call xdrfint_(ixdrf,i2set(i), iret)
1424 itmp=i_index(i,j,il,il1)
1425 call xdrfint_(ixdrf,itmp, iret)
1432 call xdrfclose_(ixdrf, iret)
1434 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1436 call xdrfint(ixdrf, i2rep(i), iret)
1439 call xdrfint(ixdrf, ifirst(i), iret)
1443 call xdrfint(ixdrf, nupa(i,il), iret)
1447 call xdrfint(ixdrf, ndowna(i,il), iret)
1453 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1460 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1467 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1473 if(usampl.or.homol_nset.gt.1) then
1474 call xdrfint(ixdrf, nset, iret)
1476 call xdrfint(ixdrf,mset(i), iret)
1479 call xdrfint(ixdrf,i2set(i), iret)
1485 itmp=i_index(i,j,il,il1)
1486 call xdrfint(ixdrf,itmp, iret)
1493 call xdrfclose(ixdrf, iret)
1500 subroutine write1traj
1501 implicit real*8 (a-h,o-z)
1502 include 'DIMENSIONS'
1505 include 'COMMON.IOUNITS'
1506 include 'COMMON.REMD'
1507 include 'COMMON.SETUP'
1508 include 'COMMON.CHAIN'
1509 include 'COMMON.SBRIDGE'
1510 include 'COMMON.INTERACT'
1514 real xcoord(3,maxres2+2),prec
1515 real r_qfrag(50),r_qpair(100)
1516 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1517 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1518 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1519 & p_uscdiff(100*maxprocs)
1520 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1521 common /przechowalnia/ p_c
1523 call mpi_bcast(ii_write,1,mpi_integer,
1524 & king,CG_COMM,ierr)
1527 print *,'traj1file',me,ii_write,ntwx_cache
1531 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1533 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1536 t5_restart1(1)=totT_cache(ii)
1537 t5_restart1(2)=EK_cache(ii)
1538 t5_restart1(3)=potE_cache(ii)
1539 t5_restart1(4)=t_bath_cache(ii)
1540 t5_restart1(5)=Uconst_cache(ii)
1541 call mpi_gather(t5_restart1,5,mpi_real,
1542 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1544 call mpi_gather(iset_cache(ii),1,mpi_integer,
1545 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1548 r_qfrag(i)=qfrag_cache(i,ii)
1551 r_qpair(i)=qpair_cache(i,ii)
1554 r_utheta(i)=utheta_cache(i,ii)
1555 r_ugamma(i)=ugamma_cache(i,ii)
1556 r_uscdiff(i)=uscdiff_cache(i,ii)
1559 call mpi_gather(r_qfrag,nfrag,mpi_real,
1560 & p_qfrag,nfrag,mpi_real,king,
1562 call mpi_gather(r_qpair,npair,mpi_real,
1563 & p_qpair,npair,mpi_real,king,
1565 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1566 & p_utheta,nfrag_back,mpi_real,king,
1568 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1569 & p_ugamma,nfrag_back,mpi_real,king,
1571 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1572 & p_uscdiff,nfrag_back,mpi_real,king,
1576 write (iout,*) "p_qfrag"
1578 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1580 write (iout,*) "p_qpair"
1582 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1588 r_c(j,i)=c_cache(j,i,ii)
1592 call mpi_gather(r_c,3*2*nres,mpi_real,
1593 & p_c,3*2*nres,mpi_real,king,
1599 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1600 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1601 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1602 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1603 call xdrfint_(ixdrf, nss, iret)
1606 call xdrfint_(ixdrf, idssb(j)+nres, iret)
1607 call xdrfint_(ixdrf, jdssb(j)+nres, iret)
1609 call xdrfint_(ixdrf, ihpb(j), iret)
1610 call xdrfint_(ixdrf, jhpb(j), iret)
1613 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1614 call xdrfint_(ixdrf, iset_restart1(il), iret)
1616 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1619 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1622 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1623 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1624 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1629 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1634 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1638 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1642 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1643 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1644 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1645 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1646 call xdrfint(ixdrf, nss, iret)
1649 call xdrfint(ixdrf, idssb(j)+nres, iret)
1650 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1652 call xdrfint(ixdrf, ihpb(j), iret)
1653 call xdrfint(ixdrf, jhpb(j), iret)
1656 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1657 call xdrfint(ixdrf, iset_restart1(il), iret)
1659 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1662 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1665 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1666 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1667 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1672 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1677 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1681 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1687 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1689 if(me.eq.king) call xdrfclose(ixdrf, iret)
1691 do i=1,ntwx_cache-ii_write
1693 totT_cache(i)=totT_cache(ii_write+i)
1694 EK_cache(i)=EK_cache(ii_write+i)
1695 potE_cache(i)=potE_cache(ii_write+i)
1696 t_bath_cache(i)=t_bath_cache(ii_write+i)
1697 Uconst_cache(i)=Uconst_cache(ii_write+i)
1698 iset_cache(i)=iset_cache(ii_write+i)
1701 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1704 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1707 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1708 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1709 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1714 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1718 ntwx_cache=ntwx_cache-ii_write
1723 subroutine read1restart(i_index)
1724 implicit real*8 (a-h,o-z)
1725 include 'DIMENSIONS'
1728 include 'COMMON.IOUNITS'
1729 include 'COMMON.REMD'
1730 include 'COMMON.SETUP'
1731 include 'COMMON.CHAIN'
1732 include 'COMMON.SBRIDGE'
1733 include 'COMMON.INTERACT'
1734 include 'COMMON.CONTROL'
1735 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1738 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1739 common /przechowalnia/ d_restart1
1740 write (*,*) "Processor",me," called read1restart"
1743 open(irest2,file=mremd_rst_name,status='unknown')
1744 read(irest2,*,err=334) i
1745 write(iout,*) "Reading old rst in ASCI format"
1747 call read1restart_old
1751 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1754 call xdrfint_(ixdrf, i2rep(i), iret)
1757 call xdrfint_(ixdrf, ifirst(i), iret)
1760 call xdrfint_(ixdrf, nupa(0,il), iret)
1762 call xdrfint_(ixdrf, nupa(i,il), iret)
1765 call xdrfint_(ixdrf, ndowna(0,il), iret)
1767 call xdrfint_(ixdrf, ndowna(i,il), iret)
1772 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1776 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1779 call xdrfint(ixdrf, i2rep(i), iret)
1782 call xdrfint(ixdrf, ifirst(i), iret)
1785 call xdrfint(ixdrf, nupa(0,il), iret)
1787 call xdrfint(ixdrf, nupa(i,il), iret)
1790 call xdrfint(ixdrf, ndowna(0,il), iret)
1792 call xdrfint(ixdrf, ndowna(i,il), iret)
1797 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1802 call mpi_scatter(t_restart1,5,mpi_real,
1803 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1807 t_bath=t5_restart1(4)
1812 c read(irest2,'(3e15.5)')
1813 c & (d_restart1(j,i+2*nres*il),j=1,3)
1816 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1818 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1823 write (iout,*) "Conformation read",il
1825 write (iout,'(i5,3f10.5,5x,3f10.5)')
1826 & i,(d_restart1(j,i+2*nres*il),j=1,3),
1827 & (d_restart1(j,nres+i+2*nres*il),j=1,3)
1832 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1833 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1843 c read(irest2,'(3e15.5)')
1844 c & (d_restart1(j,i+2*nres*il),j=1,3)
1847 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1849 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1855 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1856 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1864 if(usampl.or.homol_nset.gt.1) then
1867 call xdrfint_(ixdrf, nset, iret)
1869 call xdrfint_(ixdrf,mset(i), iret)
1872 call xdrfint_(ixdrf,i2set(i), iret)
1878 call xdrfint_(ixdrf,itmp, iret)
1879 i_index(i,j,il,il1)=itmp
1887 call xdrfint(ixdrf, nset, iret)
1889 call xdrfint(ixdrf,mset(i), iret)
1892 call xdrfint(ixdrf,i2set(i), iret)
1898 call xdrfint(ixdrf,itmp, iret)
1899 i_index(i,j,il,il1)=itmp
1906 c Corrected AL 8/19/2014: each processor needs whole iset array not only its
1908 c call mpi_scatter(i2set,1,mpi_integer,
1909 c & iset,1,mpi_integer,king,
1911 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1914 c broadcast iset to slaves
1915 if (nfgtasks.gt.1) then
1916 call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
1917 call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR)
1922 if(me.eq.king) close(irest2)
1926 subroutine read1restart_old
1927 implicit real*8 (a-h,o-z)
1928 include 'DIMENSIONS'
1931 include 'COMMON.IOUNITS'
1932 include 'COMMON.REMD'
1933 include 'COMMON.SETUP'
1934 include 'COMMON.CHAIN'
1935 include 'COMMON.SBRIDGE'
1936 include 'COMMON.INTERACT'
1937 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1939 common /przechowalnia/ d_restart1
1941 open(irest2,file=mremd_rst_name,status='unknown')
1942 read (irest2,*) (i2rep(i),i=0,nodes-1)
1943 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1945 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1946 read (irest2,*) ndowna(0,il),
1947 & (ndowna(i,il),i=1,ndowna(0,il))
1950 read(irest2,*) (t_restart1(j,il),j=1,4)
1953 call mpi_scatter(t_restart1,5,mpi_real,
1954 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1958 t_bath=t5_restart1(4)
1963 read(irest2,'(3e15.5)')
1964 & (d_restart1(j,i+2*nres*il),j=1,3)
1968 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1969 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1979 read(irest2,'(3e15.5)')
1980 & (d_restart1(j,i+2*nres*il),j=1,3)
1984 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1985 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1991 if(me.eq.king) close(irest2)
1994 c------------------------------------------