2 implicit real*8 (a-h,o-z)
5 include 'COMMON.CONTROL'
9 include 'COMMON.LANGEVIN'
11 include 'COMMON.LANGEVIN.lang0'
13 include 'COMMON.CHAIN'
14 include 'COMMON.DERIV'
16 include 'COMMON.LOCAL'
17 include 'COMMON.INTERACT'
18 include 'COMMON.IOUNITS'
19 include 'COMMON.NAMES'
20 include 'COMMON.TIME1'
22 include 'COMMON.SETUP'
24 include 'COMMON.HAIRPIN'
26 double precision cm(3),L(3),vcm(3)
27 double precision energia(0:n_ene)
28 double precision remd_t_bath(maxprocs)
29 integer iremd_iset(maxprocs)
31 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
32 double precision remd_ene(0:n_ene+8,maxprocs)
33 integer iremd_acc(maxprocs),iremd_tot(maxprocs)
34 integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
40 cold integer nup(0:maxprocs),ndown(0:maxprocs)
41 integer rep2i(0:maxprocs),ireqi(maxprocs)
42 integer icache_all(maxprocs)
43 integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
44 logical synflag,end_of_run,file_exist /.false./,ovrtim
50 if(me.eq.king.or..not.out1file) then
51 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
52 write (iout,*) "NREP=",nrep
56 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
57 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
59 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
61 cd print *,'MREMD',nodes
62 cd print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
63 cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath
67 do il1=1,max0(mset(il),1)
83 if(me.eq.king.or..not.out1file) then
84 write(iout,*) (i2rep(i),i=0,nodes-1)
85 write(iout,*) (i2set(i),i=0,nodes-1)
90 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
97 c print *,'i2rep',me,i2rep(me)
98 c print *,'rep2i',(rep2i(i),i=0,nrep)
100 cold if (i2rep(me).eq.nrep) then
103 cold nup(0)=remd_m(i2rep(me)+1)
104 cold k=rep2i(int(i2rep(me)))+1
111 cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
113 cold if (i2rep(me).eq.1) then
116 cold ndown(0)=remd_m(i2rep(me)-1)
117 cold k=rep2i(i2rep(me)-2)+1
124 cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
127 write (*,*) "Processor",me," rest",rest,"
128 & restart1fie",restart1file
129 if(rest.and.restart1file) then
131 & inquire(file=mremd_rst_name,exist=file_exist)
132 cd write (*,*) me," Before broadcast: file_exist",file_exist
133 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
135 cd write (*,*) me," After broadcast: file_exist",file_exist
137 if(me.eq.king.or..not.out1file)
138 & write (iout,*) 'Master is reading restart1file'
139 call read1restart(i_index)
141 if(me.eq.king.or..not.out1file)
142 & write (iout,*) 'WARNING : no restart1file'
145 if(me.eq.king.or..not.out1file) then
146 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
147 write(iout,*) "i_index"
152 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
161 if (rest.and..not.restart1file)
162 & inquire(file=mremd_rst_name,exist=file_exist)
163 if(.not.file_exist.and.rest.and..not.restart1file)
164 & write(iout,*) 'WARNING : no restart file',mremd_rst_name
165 IF (rest.and.file_exist.and..not.restart1file) THEN
166 write (iout,*) 'Master is reading restart file',
168 open(irest2,file=mremd_rst_name,status='unknown')
170 read (irest2,*) (i2rep(i),i=0,nodes-1)
172 read (irest2,*) (ifirst(i),i=1,remd_m(1))
175 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
177 read (irest2,*) ndowna(0,il),
178 & (ndowna(i,il),i=1,ndowna(0,il))
184 read (irest2,*) (mset(i),i=1,nset)
186 read (irest2,*) (i2set(i),i=0,nodes-1)
191 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
196 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
197 write(iout,*) "i_index"
202 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
211 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
212 write (iout,'(a6,1000i5)') "ifirst",
213 & (ifirst(i),i=1,remd_m(1))
215 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
216 & (nupa(i,il),i=1,nupa(0,il))
217 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
218 & (ndowna(i,il),i=1,ndowna(0,il))
220 ELSE IF (.not.(rest.and.file_exist)) THEN
226 if (i2rep(il-1).eq.nrep) then
229 nupa(0,il)=remd_m(i2rep(il-1)+1)
230 k=rep2i(int(i2rep(il-1)))+1
236 if (i2rep(il-1).eq.1) then
239 ndowna(0,il)=remd_m(i2rep(il-1)-1)
240 k=rep2i(i2rep(il-1)-2)+1
248 write (iout,'(a6,100i4)') "ifirst",
249 & (ifirst(i),i=1,remd_m(1))
251 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
252 & (nupa(i,il),i=1,nupa(0,il))
253 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
254 & (ndowna(i,il),i=1,ndowna(0,il))
260 c t_bath=retmin+(retmax-retmin)*me/(nodes-1)
261 if(.not.(rest.and.file_exist.and.restart1file)) then
262 if (me .eq. king) then
265 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
267 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
268 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
273 if(me.eq.king.or..not.out1file)
274 & write(iout,*) me,"iset=",iset,"t_bath=",t_bath
278 stdfp=dsqrt(2*Rb*t_bath/d_time)
280 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
283 c print *,'irep',me,t_bath
285 if (me.eq.king .or. .not. out1file)
286 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
287 call rescale_weights(t_bath)
291 c------copy MD--------------
292 c The driver for molecular dynamics subroutines
293 c------------------------------------------------
299 if(me.eq.king.or..not.out1file)
300 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
306 c Determine the inverse of the inertia matrix.
307 call setup_MD_matrices
311 if (me.eq.king .or. .not. out1file)
312 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
313 stdfp=dsqrt(2*Rb*t_bath/d_time)
315 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
317 call rescale_weights(t_bath)
321 t_MDsetup = MPI_Wtime()-tt0
323 t_MDsetup = tcpu()-tt0
326 c Entering the MD loop
332 if (lang.eq.2 .or. lang.eq.3) then
336 call sd_verlet_p_setup
338 call sd_verlet_ciccotti_setup
342 pfric0_mat(i,j,0)=pfric_mat(i,j)
343 afric0_mat(i,j,0)=afric_mat(i,j)
344 vfric0_mat(i,j,0)=vfric_mat(i,j)
345 prand0_mat(i,j,0)=prand_mat(i,j)
346 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
347 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
352 flag_stoch(i)=.false.
356 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
358 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
362 else if (lang.eq.1 .or. lang.eq.4) then
366 if (me.eq.king .or. .not. out1file)
367 & write(iout,*) 'Setup time',time00-walltime
370 t_langsetup=MPI_Wtime()-tt0
373 t_langsetup=tcpu()-tt0
378 do while(.not.end_of_run)
380 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
381 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
383 if (lang.gt.0 .and. surfarea .and.
384 & mod(itime,reset_fricmat).eq.0) then
385 if (lang.eq.2 .or. lang.eq.3) then
389 call sd_verlet_p_setup
391 call sd_verlet_ciccotti_setup
395 pfric0_mat(i,j,0)=pfric_mat(i,j)
396 afric0_mat(i,j,0)=afric_mat(i,j)
397 vfric0_mat(i,j,0)=vfric_mat(i,j)
398 prand0_mat(i,j,0)=prand_mat(i,j)
399 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
400 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
405 flag_stoch(i)=.false.
408 else if (lang.eq.1 .or. lang.eq.4) then
411 write (iout,'(a,i10)')
412 & "Friction matrix reset based on surface area, itime",itime
414 if (reset_vel .and. tbf .and. lang.eq.0
415 & .and. mod(itime,count_reset_vel).eq.0) then
417 if (me.eq.king .or. .not. out1file)
418 & write(iout,'(a,f20.2)')
419 & "Velocities reset to random values, time",totT
422 d_t_old(j,i)=d_t(j,i)
426 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
430 d_t(j,0)=d_t(j,0)-vcm(j)
433 kinetic_T=2.0d0/(dimen3*Rb)*EK
434 scalfac=dsqrt(T_bath/kinetic_T)
435 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
438 d_t_old(j,i)=scalfac*d_t(j,i)
444 c Time-reversible RESPA algorithm
445 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
446 call RESPA_step(itime)
448 c Variable time step algorithm.
449 call velverlet_step(itime)
453 call brown_step(itime)
455 print *,"Brown dynamics not here!"
457 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
463 if (mod(itime,ntwe).eq.0) call statout(itime)
465 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
466 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
468 call hairpin(.true.,nharp,iharp)
469 call secondary2(.true.)
470 call pdbout(potE,tytul,ipdb)
475 if (mod(itime,ntwx).eq.0.and.traj1file) then
476 if(ntwx_cache.lt.max_cache_traj_use) then
477 ntwx_cache=ntwx_cache+1
479 if (max_cache_traj_use.ne.1)
480 & print *,itime,"processor ",me," over cache ",ntwx_cache
483 totT_cache(i)=totT_cache(i+1)
484 EK_cache(i)=EK_cache(i+1)
485 potE_cache(i)=potE_cache(i+1)
486 t_bath_cache(i)=t_bath_cache(i+1)
487 Uconst_cache(i)=Uconst_cache(i+1)
488 iset_cache(i)=iset_cache(i+1)
491 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
494 qpair_cache(ii,i)=qpair_cache(ii,i+1)
497 utheta_cache(ii,i)=utheta_cache(ii,i+1)
498 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
499 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
505 c_cache(j,ii,i)=c_cache(j,ii,i+1)
511 totT_cache(ntwx_cache)=totT
512 EK_cache(ntwx_cache)=EK
513 potE_cache(ntwx_cache)=potE
514 t_bath_cache(ntwx_cache)=t_bath
515 Uconst_cache(ntwx_cache)=Uconst
516 iset_cache(ntwx_cache)=iset
519 qfrag_cache(i,ntwx_cache)=qfrag(i)
522 qpair_cache(i,ntwx_cache)=qpair(i)
525 utheta_cache(i,ntwx_cache)=utheta(i)
526 ugamma_cache(i,ntwx_cache)=ugamma(i)
527 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
529 C print *,'przed returnbox'
531 C call enerprint(remd_ene(0,i))
534 c_cache(j,i,ntwx_cache)=c(j,i)
539 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
540 & .and..not.restart1file) then
543 open(irest1,file=mremd_rst_name,status='unknown')
544 write (irest1,*) "i2rep"
545 write (irest1,*) (i2rep(i),i=0,nodes-1)
546 write (irest1,*) "ifirst"
547 write (irest1,*) (ifirst(i),i=1,remd_m(1))
549 write (irest1,*) "nupa",il
550 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
551 write (irest1,*) "ndowna",il
552 write (irest1,*) ndowna(0,il),
553 & (ndowna(i,il),i=1,ndowna(0,il))
556 write (irest1,*) "nset"
557 write (irest1,*) nset
558 write (irest1,*) "mset"
559 write (irest1,*) (mset(i),i=1,nset)
560 write (irest1,*) "i2set"
561 write (irest1,*) (i2set(i),i=0,nodes-1)
562 write (irest1,*) "i_index"
566 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
574 open(irest2,file=rest2name,status='unknown')
575 write(irest2,*) totT,EK,potE,totE,t_bath
577 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
580 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
583 write (irest2,*) iset
590 c forced synchronization
591 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
592 & .and. .not. mremdsync) then
594 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
596 call mpi_recv(itime_master, 1, MPI_INTEGER,
597 & 0,101,CG_COMM, status, ierr)
598 call mpi_barrier(CG_COMM, ierr)
599 cdeb if (out1file.or.traj1file) then
600 cdeb call mpi_gather(itime,1,mpi_integer,
601 cdeb & icache_all,1,mpi_integer,king,
604 & call mpi_gather(ntwx_cache,1,mpi_integer,
605 & icache_all,1,mpi_integer,king,
608 & write(iout,*) 'REMD synchro at',itime_master,itime
609 if (itime_master.ge.n_timestep .or. ovrtim())
611 ctime call flush(iout)
616 if ((mod(itime,nstex).eq.0.and.me.eq.king
617 & .or.end_of_run.and.me.eq.king )
618 & .and. .not. mremdsync ) then
621 call mpi_isend(itime,1,MPI_INTEGER,i,101,
622 & CG_COMM, ireqi(i), ierr)
623 cd write(iout,*) 'REMD synchro with',i
626 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
627 call mpi_barrier(CG_COMM, ierr)
629 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
630 if (out1file.or.traj1file) then
631 cdeb call mpi_gather(itime,1,mpi_integer,
632 cdeb & itime_all,1,mpi_integer,king,
634 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
635 cdeb & (itime_all(i),i=1,nodes)
637 cdeb imin_itime=itime_all(1)
639 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
641 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
642 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
643 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
644 call mpi_gather(ntwx_cache,1,mpi_integer,
645 & icache_all,1,mpi_integer,king,
647 c write(iout,'(a19,8000i8)') ' ntwx_cache',
648 c & (icache_all(i),i=1,nodes)
649 ii_write=icache_all(1)
651 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
653 c write(iout,*) "MIN ii_write=",ii_write
656 ctime call flush(iout)
658 if(mremdsync .and. mod(itime,nstex).eq.0) then
660 if (me.eq.king .or. .not. out1file)
661 & write(iout,*) 'REMD synchro at',itime
664 call mpi_gather(ntwx_cache,1,mpi_integer,
665 & icache_all,1,mpi_integer,king,
668 write(iout,'(a19,8000i8)') ' ntwx_cache',
669 & (icache_all(i),i=1,nodes)
670 ii_write=icache_all(1)
672 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
674 write(iout,*) "MIN ii_write=",ii_write
680 c Update the time safety limiy
681 if (time001-time00.gt.safety) then
682 safety=time001-time00+600
683 write (iout,*) "****** SAFETY increased to",safety," s"
685 if (ovrtim()) end_of_run=.true.
687 if(synflag.and..not.end_of_run) then
691 write(iout,*) 'REMD before',me,t_bath
693 c call mpi_gather(t_bath,1,mpi_double_precision,
694 c & remd_t_bath,1,mpi_double_precision,king,
696 potEcomp(n_ene+1)=t_bath
698 potEcomp(n_ene+2)=iset
699 if (iset.lt.nset) then
704 call Econstr_back_qlike
708 potEcomp(n_ene+3)=Uconst+Uconst_back
716 call Econstr_back_qlike
720 potEcomp(n_ene+4)=Uconst+Uconst_back
724 c 9/11/17 Adaptive US
727 write (iout,*) "me ",me," itt",itt
729 if (itt.lt.nrep) then
733 potEcomp(n_ene+5)=Uconst
735 write (iout,*) "t_bath",t_bath_temp,t_bath,
738 if (iset.lt.nset) then
742 potEcomp(n_ene+7)=Uconst
744 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
754 potEcomp(n_ene+6)=Uconst
756 write (iout,*) "t_bath",t_bath_temp,t_bath,
763 potEcomp(n_ene+8)=Uconst
765 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
773 call mpi_gather(potEcomp(0),n_ene+9,mpi_double_precision,
774 & remd_ene(0,1),n_ene+9,mpi_double_precision,king,
777 call mpi_gather(elow,1,mpi_double_precision,
778 & elowi,1,mpi_double_precision,king,
780 call mpi_gather(ehigh,1,mpi_double_precision,
781 & ehighi,1,mpi_double_precision,king,
786 if (me.eq.king .or. .not. out1file) then
787 write(iout,*) 'REMD gather times=',time03-time01
791 if (restart1file) call write1rst(i_index)
794 if (me.eq.king .or. .not. out1file) then
795 write(iout,*) 'REMD writing rst time=',time04-time03
798 if (traj1file) call write1traj
800 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
801 cdeb & icache_all,1,mpi_integer,king,
803 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
804 cdeb & (icache_all(i),i=1,nodes)
809 if (me.eq.king .or. .not. out1file) then
810 write(iout,*) 'REMD writing traj time=',time05-time04
817 remd_t_bath(i)=remd_ene(n_ene+1,i)
818 iremd_iset(i)=remd_ene(n_ene+2,i)
822 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
824 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
828 write(iout,*) 'REMD exchange temp,ene'
830 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
831 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
835 c-------------------------------------
838 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
840 write (iout,*) "remd_m(1)",remd_m(1)
843 i=ifirst(iran_num(1,remd_m(1)))
848 & write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
849 if(i.gt.0.and.nupa(0,i).gt.0) then
851 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
853 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
855 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
857 c do while (iex.eq.i)
858 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
859 iex=nupa(iran_num(1,int(nupa(0,i))),i)
861 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
864 write (iout,*) "i",i,"iex",iex," temperatures",
865 & remd_t_bath(i),remd_t_bath(iex)
868 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
870 c Swap temperatures between conformations i and iex with recalculating the free energies
871 c following temperature changes.
872 ene_iex_iex=remd_ene(0,iex)
873 ene_i_i=remd_ene(0,i)
874 c write (iout,*) "i",i," ene_i_i",ene_i_i,
875 c & " iex",iex," ene_iex_iex",ene_iex_iex
876 c write (iout,*) "rescaling weights with temperature",
879 call rescale_weights(remd_t_bath(i))
881 c write (iout,*) "0,iex",remd_t_bath(i)
882 c call enerprint(remd_ene(0,iex))
884 call sum_energy(remd_ene(0,iex),.false.)
885 ene_iex_i=remd_ene(0,iex)
886 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
888 c write (iout,*) "0,i",remd_t_bath(i)
889 c call enerprint(remd_ene(0,i))
891 call sum_energy(remd_ene(0,i),.false.)
892 c write (iout,*) "ene_i_i",remd_ene(0,i)
894 c write (iout,*) "rescaling weights with temperature",
896 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
897 write (iout,*) "ERROR: inconsistent energies:",i,
898 & ene_i_i,remd_ene(0,i)
900 call rescale_weights(remd_t_bath(iex))
902 c write (iout,*) "0,i",remd_t_bath(iex)
903 c call enerprint(remd_ene(0,i))
905 call sum_energy(remd_ene(0,i),.false.)
906 c write (iout,*) "ene_i_iex",remd_ene(0,i)
908 ene_i_iex=remd_ene(0,i)
910 c write (iout,*) "0,iex",remd_t_bath(iex)
911 c call enerprint(remd_ene(0,iex))
913 call sum_energy(remd_ene(0,iex),.false.)
914 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
915 write (iout,*) "ERROR: inconsistent energies:",iex,
916 & ene_iex_iex,remd_ene(0,iex)
919 write (iout,*) "ene_iex_iex",remd_ene(0,iex)
920 write (iout,*) "i",i," iex",iex
921 write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
922 & " ene_i_iex",ene_i_iex,
923 & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
926 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
927 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
929 c write(iout,*) 'delta',delta
930 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
931 c & (remd_ene(i)-remd_ene(iex))/Rb/
932 c & (remd_t_bath(i)*remd_t_bath(iex))
934 if (delta .gt. 50.0d0) then
940 else if (delta.lt.-50.0d0) then
949 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
950 xxx=ran_number(0.0d0,1.0d0)
951 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
953 if (delta .gt. xxx) then
955 remd_t_bath(i)=remd_t_bath(iex)
957 remd_ene(0,i)=ene_i_iex
958 remd_ene(0,iex)=ene_iex_i
964 ehighi(i)=ehighi(iex)
971 nupa(k,i)=nupa(k,iex)
974 ndowna(k,i)=ndowna(k,iex)
978 if (ifirst(il).eq.i) ifirst(il)=iex
980 if (nupa(k,il).eq.i) then
982 elseif (nupa(k,il).eq.iex) then
987 if (ndowna(k,il).eq.i) then
989 elseif (ndowna(k,il).eq.iex) then
995 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
997 i2rep(i-1)=i2rep(iex-1)
1000 write(iout,*) 'exchange',i,iex
1001 write (iout,'(a8,100i4)') "@ ifirst",
1002 & (ifirst(k),k=1,remd_m(1))
1004 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1005 & (nupa(k,il),k=1,nupa(0,il))
1006 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1007 & (ndowna(k,il),k=1,ndowna(0,il))
1012 remd_ene(0,iex)=ene_iex_iex
1013 remd_ene(0,i)=ene_i_i
1019 cd write (iout,*) "exchange completed"
1023 cd write(iout,*) "########",ii
1025 i_temp=iran_num(1,nrep)
1026 i_mult=iran_num(1,remd_m(i_temp))
1027 i_iset=iran_num(1,nset)
1028 i_mset=iran_num(1,mset(i_iset))
1029 i=i_index(i_temp,i_mult,i_iset,i_mset)
1031 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1034 cd write(iout,*) "i_dir=",i_dir
1036 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1039 i_mult1=iran_num(1,remd_m(i_temp1))
1041 i_mset1=iran_num(1,mset(i_iset1))
1042 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1043 c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned
1044 c on failed replica exchange attempt
1045 econstr_temp_i=remd_ene(20,i)
1046 econstr_temp_iex=remd_ene(20,iex)
1047 c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
1049 remd_ene(20,i)=remd_ene(n_ene+5,i)
1050 remd_ene(20,iex)=remd_ene(n_ene+6,iex)
1052 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1055 i_mult1=iran_num(1,remd_m(i_temp1))
1057 i_mset1=iran_num(1,mset(i_iset1))
1058 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1059 econstr_temp_i=remd_ene(20,i)
1060 econstr_temp_iex=remd_ene(20,iex)
1061 remd_ene(20,i)=remd_ene(n_ene+3,i)
1062 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1064 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1067 i_mult1=iran_num(1,remd_m(i_temp1))
1069 i_mset1=iran_num(1,mset(i_iset1))
1070 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1071 econstr_temp_i=remd_ene(20,i)
1072 econstr_temp_iex=remd_ene(20,iex)
1074 remd_ene(20,i)=remd_ene(n_ene+7,i)
1075 remd_ene(20,iex)=remd_ene(n_ene+8,iex)
1077 remd_ene(20,i)=remd_ene(n_ene+3,i)
1078 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1084 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1087 c Swap temperatures between conformations i and iex with recalculating the free energies
1088 c following temperature changes.
1090 write (iout,*) "i_dir",i_dir," i",i," iex",iex,
1091 & " t_bath",remd_t_bath(i),remd_t_bath(iex)
1093 ene_iex_iex=remd_ene(0,iex)
1094 ene_i_i=remd_ene(0,i)
1095 co write (iout,*) "rescaling weights with temperature",
1097 call rescale_weights(remd_t_bath(i))
1099 call sum_energy(remd_ene(0,iex),.false.)
1100 ene_iex_i=remd_ene(0,iex)
1101 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
1102 c call sum_energy(remd_ene(0,i),.false.)
1103 cd write (iout,*) "ene_i_i",remd_ene(0,i)
1104 c write (iout,*) "rescaling weights with temperature",
1105 c & remd_t_bath(iex)
1106 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1107 c write (iout,*) "ERROR: inconsistent energies:",i,
1108 c & ene_i_i,remd_ene(0,i)
1110 call rescale_weights(remd_t_bath(iex))
1111 call sum_energy(remd_ene(0,i),.false.)
1112 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1113 ene_i_iex=remd_ene(0,i)
1114 c call sum_energy(remd_ene(0,iex),.false.)
1115 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1116 c write (iout,*) "ERROR: inconsistent energies:",iex,
1117 c & ene_iex_iex,remd_ene(0,iex)
1119 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1120 c write (iout,*) "i",i," iex",iex
1121 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1122 cd & " ene_i_iex",ene_i_iex,
1123 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1124 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1125 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1127 cd write(iout,*) 'delta',delta
1128 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1129 c & (remd_ene(i)-remd_ene(iex))/Rb/
1130 c & (remd_t_bath(i)*remd_t_bath(iex))
1131 if (delta .gt. 50.0d0) then
1136 if (i_dir.eq.1.or.i_dir.eq.3)
1137 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1138 if (i_dir.eq.2.or.i_dir.eq.3)
1139 & iremd_tot_usa(int(i2set(i-1)))=
1140 & iremd_tot_usa(int(i2set(i-1)))+1
1141 xxx=ran_number(0.0d0,1.0d0)
1142 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1143 if (delta .gt. xxx) then
1145 remd_t_bath(i)=remd_t_bath(iex)
1146 remd_t_bath(iex)=tmp
1149 iremd_iset(i)=iremd_iset(iex)
1150 iremd_iset(iex)=itmp
1152 remd_ene(0,i)=ene_i_iex
1153 remd_ene(0,iex)=ene_iex_i
1155 if (i_dir.eq.1.or.i_dir.eq.3)
1156 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1159 i2rep(i-1)=i2rep(iex-1)
1162 if (i_dir.eq.2.or.i_dir.eq.3)
1163 & iremd_acc_usa(int(i2set(i-1)))=
1164 & iremd_acc_usa(int(i2set(i-1)))+1
1167 i2set(i-1)=i2set(iex-1)
1170 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1171 i_index(i_temp,i_mult,i_iset,i_mset)=
1172 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1173 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1176 remd_ene(0,iex)=ene_iex_iex
1177 remd_ene(0,i)=ene_i_i
1178 remd_ene(20,iex)=econstr_temp_iex
1179 remd_ene(20,i)=econstr_temp_i
1183 cd do il1=1,mset(il)
1186 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1199 c-------------------------------------
1200 write (iout,*) "NREP",nrep
1202 if(iremd_tot(i).ne.0)
1203 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1204 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1209 if(iremd_tot_usa(i).ne.0)
1210 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1211 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1217 cd write (iout,'(a6,100i4)') "ifirst",
1218 cd & (ifirst(i),i=1,remd_m(1))
1220 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1221 cd & (nupa(i,il),i=1,nupa(0,il))
1222 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1223 cd & (ndowna(i,il),i=1,ndowna(0,il))
1228 cd write (iout,*) "Before scatter"
1230 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1231 & t_bath,1,mpi_double_precision,king,
1233 cd write (iout,*) "After scatter"
1236 & call mpi_scatter(iremd_iset,1,mpi_integer,
1237 & iset,1,mpi_integer,king,
1241 if (me.eq.king .or. .not. out1file) then
1242 write(iout,*) 'REMD scatter time=',time07-time06
1246 call mpi_scatter(elowi,1,mpi_double_precision,
1247 & elow,1,mpi_double_precision,king,
1249 call mpi_scatter(ehighi,1,mpi_double_precision,
1250 & ehigh,1,mpi_double_precision,king,
1253 call rescale_weights(t_bath)
1254 co write (iout,*) "Processor",me,
1255 co & " rescaling weights with temperature",t_bath
1257 stdfp=dsqrt(2*Rb*t_bath/d_time)
1259 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1262 cde write(iout,*) 'REMD after',me,t_bath
1264 if (me.eq.king .or. .not. out1file) then
1265 write(iout,*) 'REMD exchange time=',time08-time00
1271 if (restart1file) then
1272 if (me.eq.king .or. .not. out1file)
1273 & write(iout,*) 'writing restart at the end of run'
1274 call write1rst(i_index)
1277 if (traj1file) call write1traj
1279 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1280 cdeb & icache_all,1,mpi_integer,king,
1281 cdeb & CG_COMM,ierr)
1282 cdeb write(iout,'(a40,8000i8)')
1283 cdeb & ' ntwx_cache after traj1file at the end',
1284 cdeb & (icache_all(i),i=1,nodes)
1289 t_MD=MPI_Wtime()-tt0
1293 if (me.eq.king .or. .not. out1file) then
1294 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1296 & 'MD calculations setup:',t_MDsetup,
1297 & 'Energy & gradient evaluation:',t_enegrad,
1298 & 'Stochastic MD setup:',t_langsetup,
1299 & 'Stochastic MD step setup:',t_sdsetup,
1301 write (iout,'(/28(1h=),a25,27(1h=))')
1302 & ' End of MD calculation '
1307 c-----------------------------------------------------------------------
1308 subroutine write1rst(i_index)
1309 implicit real*8 (a-h,o-z)
1310 include 'DIMENSIONS'
1313 include 'COMMON.IOUNITS'
1314 include 'COMMON.REMD'
1315 include 'COMMON.SETUP'
1316 include 'COMMON.CHAIN'
1317 include 'COMMON.SBRIDGE'
1318 include 'COMMON.INTERACT'
1320 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1321 & d_restart2(3,2*maxres*maxprocs)
1325 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1326 common /przechowalnia/ d_restart1,d_restart2
1331 t5_restart1(4)=t_bath
1332 t5_restart1(5)=Uconst
1334 call mpi_gather(t5_restart1,5,mpi_real,
1335 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1343 call mpi_gather(r_d,3*2*nres,mpi_real,
1344 & d_restart1,3*2*nres,mpi_real,king,
1353 call mpi_gather(r_d,3*2*nres,mpi_real,
1354 & d_restart2,3*2*nres,mpi_real,king,
1359 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1361 call xdrfint_(ixdrf, i2rep(i), iret)
1364 call xdrfint_(ixdrf, ifirst(i), iret)
1368 call xdrfint_(ixdrf, nupa(i,il), iret)
1372 call xdrfint_(ixdrf, ndowna(i,il), iret)
1378 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1385 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1392 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1398 call xdrfint_(ixdrf, nset, iret)
1400 call xdrfint_(ixdrf,mset(i), iret)
1403 call xdrfint_(ixdrf,i2set(i), iret)
1409 itmp=i_index(i,j,il,il1)
1410 call xdrfint_(ixdrf,itmp, iret)
1417 call xdrfclose_(ixdrf, iret)
1419 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1421 call xdrfint(ixdrf, i2rep(i), iret)
1424 call xdrfint(ixdrf, ifirst(i), iret)
1428 call xdrfint(ixdrf, nupa(i,il), iret)
1432 call xdrfint(ixdrf, ndowna(i,il), iret)
1438 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1445 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1452 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1459 call xdrfint(ixdrf, nset, iret)
1461 call xdrfint(ixdrf,mset(i), iret)
1464 call xdrfint(ixdrf,i2set(i), iret)
1470 itmp=i_index(i,j,il,il1)
1471 call xdrfint(ixdrf,itmp, iret)
1478 call xdrfclose(ixdrf, iret)
1485 subroutine write1traj
1486 implicit real*8 (a-h,o-z)
1487 include 'DIMENSIONS'
1490 include 'COMMON.IOUNITS'
1491 include 'COMMON.REMD'
1492 include 'COMMON.SETUP'
1493 include 'COMMON.CHAIN'
1494 include 'COMMON.SBRIDGE'
1495 include 'COMMON.INTERACT'
1499 real xcoord(3,maxres2+2),prec
1500 real r_qfrag(50),r_qpair(100)
1501 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1502 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1503 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1504 & p_uscdiff(100*maxprocs)
1505 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1506 common /przechowalnia/ p_c
1508 call mpi_bcast(ii_write,1,mpi_integer,
1509 & king,CG_COMM,ierr)
1512 c print *,'traj1file',me,ii_write,ntwx_cache
1516 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1518 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1521 t5_restart1(1)=totT_cache(ii)
1522 t5_restart1(2)=EK_cache(ii)
1523 t5_restart1(3)=potE_cache(ii)
1524 t5_restart1(4)=t_bath_cache(ii)
1525 t5_restart1(5)=Uconst_cache(ii)
1526 call mpi_gather(t5_restart1,5,mpi_real,
1527 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1529 call mpi_gather(iset_cache(ii),1,mpi_integer,
1530 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1533 r_qfrag(i)=qfrag_cache(i,ii)
1536 r_qpair(i)=qpair_cache(i,ii)
1539 r_utheta(i)=utheta_cache(i,ii)
1540 r_ugamma(i)=ugamma_cache(i,ii)
1541 r_uscdiff(i)=uscdiff_cache(i,ii)
1544 call mpi_gather(r_qfrag,nfrag,mpi_real,
1545 & p_qfrag,nfrag,mpi_real,king,
1547 call mpi_gather(r_qpair,npair,mpi_real,
1548 & p_qpair,npair,mpi_real,king,
1550 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1551 & p_utheta,nfrag_back,mpi_real,king,
1553 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1554 & p_ugamma,nfrag_back,mpi_real,king,
1556 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1557 & p_uscdiff,nfrag_back,mpi_real,king,
1561 write (iout,*) "p_qfrag"
1563 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1565 write (iout,*) "p_qpair"
1567 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1573 r_c(j,i)=c_cache(j,i,ii)
1577 call mpi_gather(r_c,3*2*nres,mpi_real,
1578 & p_c,3*2*nres,mpi_real,king,
1584 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1585 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1586 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1587 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1588 call xdrfint_(ixdrf, nss, iret)
1591 call xdrfint(ixdrf, idssb(j)+nres, iret)
1592 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1594 call xdrfint_(ixdrf, ihpb(j), iret)
1595 call xdrfint_(ixdrf, jhpb(j), iret)
1598 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1599 call xdrfint_(ixdrf, iset_restart1(il), iret)
1601 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1604 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1607 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1608 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1609 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1614 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1619 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1623 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1627 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1628 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1629 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1630 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1631 call xdrfint(ixdrf, nss, iret)
1634 call xdrfint(ixdrf, idssb(j)+nres, iret)
1635 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1637 call xdrfint(ixdrf, ihpb(j), iret)
1638 call xdrfint(ixdrf, jhpb(j), iret)
1641 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1642 call xdrfint(ixdrf, iset_restart1(il), iret)
1644 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1647 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1650 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1651 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1652 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1657 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1662 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1666 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1672 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1674 if(me.eq.king) call xdrfclose(ixdrf, iret)
1676 do i=1,ntwx_cache-ii_write
1678 totT_cache(i)=totT_cache(ii_write+i)
1679 EK_cache(i)=EK_cache(ii_write+i)
1680 potE_cache(i)=potE_cache(ii_write+i)
1681 t_bath_cache(i)=t_bath_cache(ii_write+i)
1682 Uconst_cache(i)=Uconst_cache(ii_write+i)
1683 iset_cache(i)=iset_cache(ii_write+i)
1686 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1689 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1692 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1693 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1694 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1699 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1703 ntwx_cache=ntwx_cache-ii_write
1708 subroutine read1restart(i_index)
1709 implicit real*8 (a-h,o-z)
1710 include 'DIMENSIONS'
1713 include 'COMMON.IOUNITS'
1714 include 'COMMON.REMD'
1715 include 'COMMON.SETUP'
1716 include 'COMMON.CHAIN'
1717 include 'COMMON.SBRIDGE'
1718 include 'COMMON.INTERACT'
1719 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1722 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1723 common /przechowalnia/ d_restart1
1724 write (*,*) "Processor",me," called read1restart"
1727 open(irest2,file=mremd_rst_name,status='unknown')
1728 read(irest2,*,err=334) i
1729 write(iout,*) "Reading old rst in ASCI format"
1731 call read1restart_old
1735 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1738 call xdrfint_(ixdrf, i2rep(i), iret)
1741 call xdrfint_(ixdrf, ifirst(i), iret)
1744 call xdrfint_(ixdrf, nupa(0,il), iret)
1746 call xdrfint_(ixdrf, nupa(i,il), iret)
1749 call xdrfint_(ixdrf, ndowna(0,il), iret)
1751 call xdrfint_(ixdrf, ndowna(i,il), iret)
1756 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1760 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1763 call xdrfint(ixdrf, i2rep(i), iret)
1766 call xdrfint(ixdrf, ifirst(i), iret)
1769 call xdrfint(ixdrf, nupa(0,il), iret)
1771 call xdrfint(ixdrf, nupa(i,il), iret)
1774 call xdrfint(ixdrf, ndowna(0,il), iret)
1776 call xdrfint(ixdrf, ndowna(i,il), iret)
1781 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1786 call mpi_scatter(t_restart1,5,mpi_real,
1787 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1791 t_bath=t5_restart1(4)
1796 c read(irest2,'(3e15.5)')
1797 c & (d_restart1(j,i+2*nres*il),j=1,3)
1800 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1802 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1808 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1809 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1819 c read(irest2,'(3e15.5)')
1820 c & (d_restart1(j,i+2*nres*il),j=1,3)
1823 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1825 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1831 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1832 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1843 call xdrfint_(ixdrf, nset, iret)
1845 call xdrfint_(ixdrf,mset(i), iret)
1848 call xdrfint_(ixdrf,i2set(i), iret)
1854 call xdrfint_(ixdrf,itmp, iret)
1855 i_index(i,j,il,il1)=itmp
1863 call xdrfint(ixdrf, nset, iret)
1865 call xdrfint(ixdrf,mset(i), iret)
1868 call xdrfint(ixdrf,i2set(i), iret)
1874 call xdrfint(ixdrf,itmp, iret)
1875 i_index(i,j,il,il1)=itmp
1882 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1884 c call mpi_scatter(i2set,1,mpi_integer,
1885 c & iset,1,mpi_integer,king,
1887 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1894 if(me.eq.king) close(irest2)
1898 subroutine read1restart_old
1899 implicit real*8 (a-h,o-z)
1900 include 'DIMENSIONS'
1903 include 'COMMON.IOUNITS'
1904 include 'COMMON.REMD'
1905 include 'COMMON.SETUP'
1906 include 'COMMON.CHAIN'
1907 include 'COMMON.SBRIDGE'
1908 include 'COMMON.INTERACT'
1909 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1911 common /przechowalnia/ d_restart1
1913 open(irest2,file=mremd_rst_name,status='unknown')
1914 read (irest2,*) (i2rep(i),i=0,nodes-1)
1915 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1917 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1918 read (irest2,*) ndowna(0,il),
1919 & (ndowna(i,il),i=1,ndowna(0,il))
1922 read(irest2,*) (t_restart1(j,il),j=1,4)
1925 call mpi_scatter(t_restart1,5,mpi_real,
1926 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1930 t_bath=t5_restart1(4)
1935 read(irest2,'(3e15.5)')
1936 & (d_restart1(j,i+2*nres*il),j=1,3)
1940 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1941 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1951 read(irest2,'(3e15.5)')
1952 & (d_restart1(j,i+2*nres*il),j=1,3)
1956 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1957 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1963 if(me.eq.king) close(irest2)
1966 c------------------------------------------