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 iperm(maxprocs)
30 c 8/20/17 AL: replica permutation table
31 integer iremd_iset(maxprocs)
33 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
34 double precision remd_ene(0:n_ene+4,maxprocs)
35 integer iremd_acc(maxprocs),iremd_tot(maxprocs)
36 integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
42 cold integer nup(0:maxprocs),ndown(0:maxprocs)
43 integer rep2i(0:maxprocs),ireqi(maxprocs)
44 integer icache_all(maxprocs)
45 integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
46 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
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
69 do il1=1,max0(mset(il),1)
85 if(me.eq.king.or..not.out1file) then
86 write(iout,*) (i2rep(i),i=0,nodes-1)
87 write(iout,*) (i2set(i),i=0,nodes-1)
92 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
99 c print *,'i2rep',me,i2rep(me)
100 c print *,'rep2i',(rep2i(i),i=0,nrep)
102 cold if (i2rep(me).eq.nrep) then
105 cold nup(0)=remd_m(i2rep(me)+1)
106 cold k=rep2i(int(i2rep(me)))+1
113 cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
115 cold if (i2rep(me).eq.1) then
118 cold ndown(0)=remd_m(i2rep(me)-1)
119 cold k=rep2i(i2rep(me)-2)+1
126 cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
129 write (*,*) "Processor",me," rest",rest,"
130 & restart1fie",restart1file
131 if(rest.and.restart1file) then
133 & inquire(file=mremd_rst_name,exist=file_exist)
134 cd write (*,*) me," Before broadcast: file_exist",file_exist
135 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
137 cd write (*,*) me," After broadcast: file_exist",file_exist
139 if(me.eq.king.or..not.out1file)
140 & write (iout,*) 'Master is reading restart1file'
141 call read1restart(i_index)
143 if(me.eq.king.or..not.out1file)
144 & write (iout,*) 'WARNING : no restart1file'
147 if(me.eq.king.or..not.out1file) then
148 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
149 write(iout,*) "i_index"
154 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
163 if (rest.and..not.restart1file)
164 & inquire(file=mremd_rst_name,exist=file_exist)
165 if(.not.file_exist.and.rest.and..not.restart1file)
166 & write(iout,*) 'WARNING : no restart file',mremd_rst_name
167 IF (rest.and.file_exist.and..not.restart1file) THEN
168 write (iout,*) 'Master is reading restart file',
170 open(irest2,file=mremd_rst_name,status='unknown')
172 read (irest2,*) (i2rep(i),i=0,nodes-1)
174 read (irest2,*) (ifirst(i),i=1,remd_m(1))
177 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
179 read (irest2,*) ndowna(0,il),
180 & (ndowna(i,il),i=1,ndowna(0,il))
186 read (irest2,*) (mset(i),i=1,nset)
188 read (irest2,*) (i2set(i),i=0,nodes-1)
193 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
198 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
199 write(iout,*) "i_index"
204 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
213 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
214 write (iout,'(a6,1000i5)') "ifirst",
215 & (ifirst(i),i=1,remd_m(1))
217 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
218 & (nupa(i,il),i=1,nupa(0,il))
219 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
220 & (ndowna(i,il),i=1,ndowna(0,il))
222 ELSE IF (.not.(rest.and.file_exist)) THEN
228 if (i2rep(il-1).eq.nrep) then
231 nupa(0,il)=remd_m(i2rep(il-1)+1)
232 k=rep2i(int(i2rep(il-1)))+1
238 if (i2rep(il-1).eq.1) then
241 ndowna(0,il)=remd_m(i2rep(il-1)-1)
242 k=rep2i(i2rep(il-1)-2)+1
250 write (iout,'(a6,100i4)') "ifirst",
251 & (ifirst(i),i=1,remd_m(1))
253 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
254 & (nupa(i,il),i=1,nupa(0,il))
255 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
256 & (ndowna(i,il),i=1,ndowna(0,il))
262 c t_bath=retmin+(retmax-retmin)*me/(nodes-1)
263 if(.not.(rest.and.file_exist.and.restart1file)) then
264 if (me .eq. king) then
267 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
269 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
270 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
275 if(me.eq.king.or..not.out1file)
276 & write(iout,*) me,"iset=",iset,"t_bath=",t_bath
279 stdfp=dsqrt(2*Rb*t_bath/d_time)
281 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
284 c print *,'irep',me,t_bath
286 if (me.eq.king .or. .not. out1file)
287 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
288 call rescale_weights(t_bath)
292 c------copy MD--------------
293 c The driver for molecular dynamics subroutines
294 c------------------------------------------------
300 if(me.eq.king.or..not.out1file)
301 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
307 c Determine the inverse of the inertia matrix.
308 call setup_MD_matrices
312 if (me.eq.king .or. .not. out1file)
313 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
314 stdfp=dsqrt(2*Rb*t_bath/d_time)
316 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
318 call rescale_weights(t_bath)
322 t_MDsetup = MPI_Wtime()-tt0
324 t_MDsetup = tcpu()-tt0
327 c Entering the MD loop
333 if (lang.eq.2 .or. lang.eq.3) then
337 call sd_verlet_p_setup
339 call sd_verlet_ciccotti_setup
343 pfric0_mat(i,j,0)=pfric_mat(i,j)
344 afric0_mat(i,j,0)=afric_mat(i,j)
345 vfric0_mat(i,j,0)=vfric_mat(i,j)
346 prand0_mat(i,j,0)=prand_mat(i,j)
347 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
348 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
353 flag_stoch(i)=.false.
357 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
359 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
363 else if (lang.eq.1 .or. lang.eq.4) then
367 if (me.eq.king .or. .not. out1file)
368 & write(iout,*) 'Setup time',time00-walltime
371 t_langsetup=MPI_Wtime()-tt0
374 t_langsetup=tcpu()-tt0
379 do while(.not.end_of_run)
381 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
382 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
384 if (lang.gt.0 .and. surfarea .and.
385 & mod(itime,reset_fricmat).eq.0) then
386 if (lang.eq.2 .or. lang.eq.3) then
390 call sd_verlet_p_setup
392 call sd_verlet_ciccotti_setup
396 pfric0_mat(i,j,0)=pfric_mat(i,j)
397 afric0_mat(i,j,0)=afric_mat(i,j)
398 vfric0_mat(i,j,0)=vfric_mat(i,j)
399 prand0_mat(i,j,0)=prand_mat(i,j)
400 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
401 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
406 flag_stoch(i)=.false.
409 else if (lang.eq.1 .or. lang.eq.4) then
412 write (iout,'(a,i10)')
413 & "Friction matrix reset based on surface area, itime",itime
415 if (reset_vel .and. tbf .and. lang.eq.0
416 & .and. mod(itime,count_reset_vel).eq.0) then
418 if (me.eq.king .or. .not. out1file)
419 & write(iout,'(a,f20.2)')
420 & "Velocities reset to random values, time",totT
423 d_t_old(j,i)=d_t(j,i)
427 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
431 d_t(j,0)=d_t(j,0)-vcm(j)
434 kinetic_T=2.0d0/(dimen3*Rb)*EK
435 scalfac=dsqrt(T_bath/kinetic_T)
436 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
439 d_t_old(j,i)=scalfac*d_t(j,i)
445 c Time-reversible RESPA algorithm
446 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
447 call RESPA_step(itime)
449 c Variable time step algorithm.
450 call velverlet_step(itime)
454 call brown_step(itime)
456 print *,"Brown dynamics not here!"
458 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
464 if (mod(itime,ntwe).eq.0) call statout(itime)
466 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
467 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
469 call hairpin(.true.,nharp,iharp)
470 call secondary2(.true.)
471 call pdbout(potE,tytul,ipdb)
476 if (mod(itime,ntwx).eq.0.and.traj1file) then
477 if(ntwx_cache.lt.max_cache_traj_use) then
478 ntwx_cache=ntwx_cache+1
480 if (max_cache_traj_use.ne.1)
481 & print *,itime,"processor ",me," over cache ",ntwx_cache
484 totT_cache(i)=totT_cache(i+1)
485 EK_cache(i)=EK_cache(i+1)
486 potE_cache(i)=potE_cache(i+1)
487 t_bath_cache(i)=t_bath_cache(i+1)
488 Uconst_cache(i)=Uconst_cache(i+1)
489 iset_cache(i)=iset_cache(i+1)
492 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
495 qpair_cache(ii,i)=qpair_cache(ii,i+1)
498 utheta_cache(ii,i)=utheta_cache(ii,i+1)
499 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
500 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
506 c_cache(j,ii,i)=c_cache(j,ii,i+1)
512 totT_cache(ntwx_cache)=totT
513 EK_cache(ntwx_cache)=EK
514 potE_cache(ntwx_cache)=potE
515 t_bath_cache(ntwx_cache)=t_bath
516 Uconst_cache(ntwx_cache)=Uconst
517 iset_cache(ntwx_cache)=iset
520 qfrag_cache(i,ntwx_cache)=qfrag(i)
523 qpair_cache(i,ntwx_cache)=qpair(i)
526 utheta_cache(i,ntwx_cache)=utheta(i)
527 ugamma_cache(i,ntwx_cache)=ugamma(i)
528 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
530 C print *,'przed returnbox'
532 C call enerprint(remd_ene(0,i))
535 c_cache(j,i,ntwx_cache)=c(j,i)
540 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
541 & .and..not.restart1file) then
544 open(irest1,file=mremd_rst_name,status='unknown')
545 write (irest1,*) "i2rep"
546 write (irest1,*) (i2rep(i),i=0,nodes-1)
547 write (irest1,*) "ifirst"
548 write (irest1,*) (ifirst(i),i=1,remd_m(1))
550 write (irest1,*) "nupa",il
551 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
552 write (irest1,*) "ndowna",il
553 write (irest1,*) ndowna(0,il),
554 & (ndowna(i,il),i=1,ndowna(0,il))
557 write (irest1,*) "nset"
558 write (irest1,*) nset
559 write (irest1,*) "mset"
560 write (irest1,*) (mset(i),i=1,nset)
561 write (irest1,*) "i2set"
562 write (irest1,*) (i2set(i),i=0,nodes-1)
563 write (irest1,*) "i_index"
567 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
575 open(irest2,file=rest2name,status='unknown')
576 write(irest2,*) totT,EK,potE,totE,t_bath
578 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
581 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
584 write (irest2,*) iset
591 c forced synchronization
592 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
593 & .and. .not. mremdsync) then
595 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
597 call mpi_recv(itime_master, 1, MPI_INTEGER,
598 & 0,101,CG_COMM, status, ierr)
599 call mpi_barrier(CG_COMM, ierr)
600 cdeb if (out1file.or.traj1file) then
601 cdeb call mpi_gather(itime,1,mpi_integer,
602 cdeb & icache_all,1,mpi_integer,king,
605 & call mpi_gather(ntwx_cache,1,mpi_integer,
606 & icache_all,1,mpi_integer,king,
609 & write(iout,*) 'REMD synchro at',itime_master,itime
610 if (itime_master.ge.n_timestep .or. ovrtim())
612 ctime call flush(iout)
617 if ((mod(itime,nstex).eq.0.and.me.eq.king
618 & .or.end_of_run.and.me.eq.king )
619 & .and. .not. mremdsync ) then
622 call mpi_isend(itime,1,MPI_INTEGER,i,101,
623 & CG_COMM, ireqi(i), ierr)
624 cd write(iout,*) 'REMD synchro with',i
627 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
628 call mpi_barrier(CG_COMM, ierr)
630 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
631 if (out1file.or.traj1file) then
632 cdeb call mpi_gather(itime,1,mpi_integer,
633 cdeb & itime_all,1,mpi_integer,king,
635 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
636 cdeb & (itime_all(i),i=1,nodes)
638 cdeb imin_itime=itime_all(1)
640 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
642 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
643 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
644 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
645 call mpi_gather(ntwx_cache,1,mpi_integer,
646 & icache_all,1,mpi_integer,king,
648 c write(iout,'(a19,8000i8)') ' ntwx_cache',
649 c & (icache_all(i),i=1,nodes)
650 ii_write=icache_all(1)
652 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
654 c write(iout,*) "MIN ii_write=",ii_write
657 ctime call flush(iout)
659 if(mremdsync .and. mod(itime,nstex).eq.0) then
661 if (me.eq.king .or. .not. out1file)
662 & write(iout,*) 'REMD synchro at',itime
665 call mpi_gather(ntwx_cache,1,mpi_integer,
666 & icache_all,1,mpi_integer,king,
669 write(iout,'(a19,8000i8)') ' ntwx_cache',
670 & (icache_all(i),i=1,nodes)
671 ii_write=icache_all(1)
673 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
675 write(iout,*) "MIN ii_write=",ii_write
681 c Update the time safety limiy
682 if (time001-time00.gt.safety) then
683 safety=time001-time00+600
684 write (iout,*) "****** SAFETY increased to",safety," s"
686 if (ovrtim()) end_of_run=.true.
688 if(synflag.and..not.end_of_run) then
692 write(iout,*) 'REMD before',me,t_bath
694 c call mpi_gather(t_bath,1,mpi_double_precision,
695 c & remd_t_bath,1,mpi_double_precision,king,
697 potEcomp(n_ene+1)=t_bath
699 potEcomp(n_ene+2)=iset
700 if (iset.lt.nset) then
705 call Econstr_back_qlike
709 potEcomp(n_ene+3)=Uconst+Uconst_back
717 call Econstr_back_qlike
721 potEcomp(n_ene+4)=Uconst+Uconst_back
725 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
726 & remd_ene(0,1),n_ene+5,mpi_double_precision,king,
729 call mpi_gather(elow,1,mpi_double_precision,
730 & elowi,1,mpi_double_precision,king,
732 call mpi_gather(ehigh,1,mpi_double_precision,
733 & ehighi,1,mpi_double_precision,king,
738 if (me.eq.king .or. .not. out1file) then
739 write(iout,*) 'REMD gather times=',time03-time01
743 if (restart1file) call write1rst(i_index)
746 if (me.eq.king .or. .not. out1file) then
747 write(iout,*) 'REMD writing rst time=',time04-time03
750 if (traj1file) call write1traj
752 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
753 cdeb & icache_all,1,mpi_integer,king,
755 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
756 cdeb & (icache_all(i),i=1,nodes)
761 if (me.eq.king .or. .not. out1file) then
762 write(iout,*) 'REMD writing traj time=',time05-time04
769 remd_t_bath(i)=remd_ene(n_ene+1,i)
770 iremd_iset(i)=remd_ene(n_ene+2,i)
774 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
776 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
780 write(iout,*) 'REMD exchange temp,ene'
782 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
783 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
787 c-------------------------------------
790 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
792 write (iout,*) "remd_m(1)",remd_m(1)
795 i=ifirst(iran_num(1,remd_m(1)))
800 & write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
801 if(i.gt.0.and.nupa(0,i).gt.0) then
803 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
805 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
807 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
809 c do while (iex.eq.i)
810 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
811 iex=nupa(iran_num(1,int(nupa(0,i))),i)
813 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
815 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
817 c Swap temperatures between conformations i and iex with recalculating the free energies
818 c following temperature changes.
819 ene_iex_iex=remd_ene(0,iex)
820 ene_i_i=remd_ene(0,i)
821 c write (iout,*) "i",i," ene_i_i",ene_i_i,
822 c & " iex",iex," ene_iex_iex",ene_iex_iex
823 c write (iout,*) "rescaling weights with temperature",
826 call rescale_weights(remd_t_bath(i))
828 c write (iout,*) "0,iex",remd_t_bath(i)
829 c call enerprint(remd_ene(0,iex))
831 call sum_energy(remd_ene(0,iex),.false.)
832 ene_iex_i=remd_ene(0,iex)
833 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
835 c write (iout,*) "0,i",remd_t_bath(i)
836 c call enerprint(remd_ene(0,i))
838 call sum_energy(remd_ene(0,i),.false.)
839 c write (iout,*) "ene_i_i",remd_ene(0,i)
841 c write (iout,*) "rescaling weights with temperature",
843 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
844 write (iout,*) "ERROR: inconsistent energies:",i,
845 & ene_i_i,remd_ene(0,i)
847 call rescale_weights(remd_t_bath(iex))
849 c write (iout,*) "0,i",remd_t_bath(iex)
850 c call enerprint(remd_ene(0,i))
852 call sum_energy(remd_ene(0,i),.false.)
853 c write (iout,*) "ene_i_iex",remd_ene(0,i)
855 ene_i_iex=remd_ene(0,i)
857 c write (iout,*) "0,iex",remd_t_bath(iex)
858 c call enerprint(remd_ene(0,iex))
860 call sum_energy(remd_ene(0,iex),.false.)
861 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
862 write (iout,*) "ERROR: inconsistent energies:",iex,
863 & ene_iex_iex,remd_ene(0,iex)
865 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
866 c write (iout,*) "i",i," iex",iex
867 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
868 c & " ene_i_iex",ene_i_iex,
869 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
871 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
872 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
874 c write(iout,*) 'delta',delta
875 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
876 c & (remd_ene(i)-remd_ene(iex))/Rb/
877 c & (remd_t_bath(i)*remd_t_bath(iex))
879 if (delta .gt. 50.0d0) then
885 else if (delta.lt.-50.0d0) then
894 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
895 xxx=ran_number(0.0d0,1.0d0)
896 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
898 if (delta .gt. xxx) then
900 remd_t_bath(i)=remd_t_bath(iex)
902 remd_ene(0,i)=ene_i_iex
903 remd_ene(0,iex)=ene_iex_i
909 ehighi(i)=ehighi(iex)
916 nupa(k,i)=nupa(k,iex)
919 ndowna(k,i)=ndowna(k,iex)
923 if (ifirst(il).eq.i) ifirst(il)=iex
925 if (nupa(k,il).eq.i) then
927 elseif (nupa(k,il).eq.iex) then
932 if (ndowna(k,il).eq.i) then
934 elseif (ndowna(k,il).eq.iex) then
940 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
942 i2rep(i-1)=i2rep(iex-1)
945 c write(iout,*) 'exchange',i,iex
946 c write (iout,'(a8,100i4)') "@ ifirst",
947 c & (ifirst(k),k=1,remd_m(1))
949 c write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
950 c & (nupa(k,il),k=1,nupa(0,il))
951 c write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
952 c & (ndowna(k,il),k=1,ndowna(0,il))
957 remd_ene(0,iex)=ene_iex_iex
958 remd_ene(0,i)=ene_i_i
964 cd write (iout,*) "exchange completed"
968 cd write(iout,*) "########",ii
970 i_temp=iran_num(1,nrep)
971 i_mult=iran_num(1,remd_m(i_temp))
972 i_iset=iran_num(1,nset)
973 i_mset=iran_num(1,mset(i_iset))
974 i=i_index(i_temp,i_mult,i_iset,i_mset)
976 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
979 cd write(iout,*) "i_dir=",i_dir
981 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
984 i_mult1=iran_num(1,remd_m(i_temp1))
986 i_mset1=iran_num(1,mset(i_iset1))
987 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
989 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
992 i_mult1=iran_num(1,remd_m(i_temp1))
994 i_mset1=iran_num(1,mset(i_iset1))
995 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
996 econstr_temp_i=remd_ene(20,i)
997 econstr_temp_iex=remd_ene(20,iex)
998 remd_ene(20,i)=remd_ene(n_ene+3,i)
999 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1001 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1004 i_mult1=iran_num(1,remd_m(i_temp1))
1006 i_mset1=iran_num(1,mset(i_iset1))
1007 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1008 econstr_temp_i=remd_ene(20,i)
1009 econstr_temp_iex=remd_ene(20,iex)
1010 remd_ene(20,i)=remd_ene(n_ene+3,i)
1011 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1017 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1020 c Swap temperatures between conformations i and iex with recalculating the free energies
1021 c following temperature changes.
1022 ene_iex_iex=remd_ene(0,iex)
1023 ene_i_i=remd_ene(0,i)
1024 co write (iout,*) "rescaling weights with temperature",
1026 call rescale_weights(remd_t_bath(i))
1028 call sum_energy(remd_ene(0,iex),.false.)
1029 ene_iex_i=remd_ene(0,iex)
1030 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
1031 c call sum_energy(remd_ene(0,i),.false.)
1032 cd write (iout,*) "ene_i_i",remd_ene(0,i)
1033 c write (iout,*) "rescaling weights with temperature",
1034 c & remd_t_bath(iex)
1035 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1036 c write (iout,*) "ERROR: inconsistent energies:",i,
1037 c & ene_i_i,remd_ene(0,i)
1039 call rescale_weights(remd_t_bath(iex))
1040 call sum_energy(remd_ene(0,i),.false.)
1041 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1042 ene_i_iex=remd_ene(0,i)
1043 c call sum_energy(remd_ene(0,iex),.false.)
1044 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1045 c write (iout,*) "ERROR: inconsistent energies:",iex,
1046 c & ene_iex_iex,remd_ene(0,iex)
1048 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1049 c write (iout,*) "i",i," iex",iex
1050 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1051 cd & " ene_i_iex",ene_i_iex,
1052 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1053 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1054 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1056 cd write(iout,*) 'delta',delta
1057 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1058 c & (remd_ene(i)-remd_ene(iex))/Rb/
1059 c & (remd_t_bath(i)*remd_t_bath(iex))
1060 if (delta .gt. 50.0d0) then
1065 if (i_dir.eq.1.or.i_dir.eq.3)
1066 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1067 if (i_dir.eq.2.or.i_dir.eq.3)
1068 & iremd_tot_usa(int(i2set(i-1)))=
1069 & iremd_tot_usa(int(i2set(i-1)))+1
1070 xxx=ran_number(0.0d0,1.0d0)
1071 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1072 if (delta .gt. xxx) then
1074 remd_t_bath(i)=remd_t_bath(iex)
1075 remd_t_bath(iex)=tmp
1078 iremd_iset(i)=iremd_iset(iex)
1079 iremd_iset(iex)=itmp
1081 remd_ene(0,i)=ene_i_iex
1082 remd_ene(0,iex)=ene_iex_i
1084 if (i_dir.eq.1.or.i_dir.eq.3)
1085 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1088 i2rep(i-1)=i2rep(iex-1)
1091 if (i_dir.eq.2.or.i_dir.eq.3)
1092 & iremd_acc_usa(int(i2set(i-1)))=
1093 & iremd_acc_usa(int(i2set(i-1)))+1
1096 i2set(i-1)=i2set(iex-1)
1099 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1100 i_index(i_temp,i_mult,i_iset,i_mset)=
1101 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1102 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1105 remd_ene(0,iex)=ene_iex_iex
1106 remd_ene(0,i)=ene_i_i
1107 remd_ene(20,iex)=econstr_temp_iex
1108 remd_ene(20,i)=econstr_temp_i
1112 cd do il1=1,mset(il)
1115 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1128 c-------------------------------------
1129 write (iout,*) "NREP",nrep
1131 if(iremd_tot(i).ne.0)
1132 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1133 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1138 if(iremd_tot_usa(i).ne.0)
1139 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1140 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1146 cd write (iout,'(a6,100i4)') "ifirst",
1147 cd & (ifirst(i),i=1,remd_m(1))
1149 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1150 cd & (nupa(i,il),i=1,nupa(0,il))
1151 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1152 cd & (ndowna(i,il),i=1,ndowna(0,il))
1157 cd write (iout,*) "Before scatter"
1159 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1160 & t_bath,1,mpi_double_precision,king,
1162 cd write (iout,*) "After scatter"
1165 & call mpi_scatter(iremd_iset,1,mpi_integer,
1166 & iset,1,mpi_integer,king,
1170 if (me.eq.king .or. .not. out1file) then
1171 write(iout,*) 'REMD scatter time=',time07-time06
1175 call mpi_scatter(elowi,1,mpi_double_precision,
1176 & elow,1,mpi_double_precision,king,
1178 call mpi_scatter(ehighi,1,mpi_double_precision,
1179 & ehigh,1,mpi_double_precision,king,
1182 call rescale_weights(t_bath)
1183 co write (iout,*) "Processor",me,
1184 co & " rescaling weights with temperature",t_bath
1186 stdfp=dsqrt(2*Rb*t_bath/d_time)
1188 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1191 cde write(iout,*) 'REMD after',me,t_bath
1193 if (me.eq.king .or. .not. out1file) then
1194 write(iout,*) 'REMD exchange time=',time08-time00
1200 if (restart1file) then
1201 if (me.eq.king .or. .not. out1file)
1202 & write(iout,*) 'writing restart at the end of run'
1203 call write1rst(i_index)
1206 if (traj1file) call write1traj
1208 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1209 cdeb & icache_all,1,mpi_integer,king,
1210 cdeb & CG_COMM,ierr)
1211 cdeb write(iout,'(a40,8000i8)')
1212 cdeb & ' ntwx_cache after traj1file at the end',
1213 cdeb & (icache_all(i),i=1,nodes)
1218 t_MD=MPI_Wtime()-tt0
1222 if (me.eq.king .or. .not. out1file) then
1223 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1225 & 'MD calculations setup:',t_MDsetup,
1226 & 'Energy & gradient evaluation:',t_enegrad,
1227 & 'Stochastic MD setup:',t_langsetup,
1228 & 'Stochastic MD step setup:',t_sdsetup,
1230 write (iout,'(/28(1h=),a25,27(1h=))')
1231 & ' End of MD calculation '
1236 c-----------------------------------------------------------------------
1237 subroutine write1rst(i_index)
1238 implicit real*8 (a-h,o-z)
1239 include 'DIMENSIONS'
1242 include 'COMMON.IOUNITS'
1243 include 'COMMON.REMD'
1244 include 'COMMON.SETUP'
1245 include 'COMMON.CHAIN'
1246 include 'COMMON.SBRIDGE'
1247 include 'COMMON.INTERACT'
1249 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1250 & d_restart2(3,2*maxres*maxprocs)
1254 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1255 common /przechowalnia/ d_restart1,d_restart2
1260 t5_restart1(4)=t_bath
1261 t5_restart1(5)=Uconst
1263 call mpi_gather(t5_restart1,5,mpi_real,
1264 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1272 call mpi_gather(r_d,3*2*nres,mpi_real,
1273 & d_restart1,3*2*nres,mpi_real,king,
1282 call mpi_gather(r_d,3*2*nres,mpi_real,
1283 & d_restart2,3*2*nres,mpi_real,king,
1288 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1290 call xdrfint_(ixdrf, i2rep(i), iret)
1293 call xdrfint_(ixdrf, ifirst(i), iret)
1297 call xdrfint_(ixdrf, nupa(i,il), iret)
1301 call xdrfint_(ixdrf, ndowna(i,il), iret)
1307 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1314 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1321 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1327 call xdrfint_(ixdrf, nset, iret)
1329 call xdrfint_(ixdrf,mset(i), iret)
1332 call xdrfint_(ixdrf,i2set(i), iret)
1338 itmp=i_index(i,j,il,il1)
1339 call xdrfint_(ixdrf,itmp, iret)
1346 call xdrfclose_(ixdrf, iret)
1348 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1350 call xdrfint(ixdrf, i2rep(i), iret)
1353 call xdrfint(ixdrf, ifirst(i), iret)
1357 call xdrfint(ixdrf, nupa(i,il), iret)
1361 call xdrfint(ixdrf, ndowna(i,il), iret)
1367 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1374 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1381 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1388 call xdrfint(ixdrf, nset, iret)
1390 call xdrfint(ixdrf,mset(i), iret)
1393 call xdrfint(ixdrf,i2set(i), iret)
1399 itmp=i_index(i,j,il,il1)
1400 call xdrfint(ixdrf,itmp, iret)
1407 call xdrfclose(ixdrf, iret)
1414 subroutine write1traj
1415 implicit real*8 (a-h,o-z)
1416 include 'DIMENSIONS'
1419 include 'COMMON.IOUNITS'
1420 include 'COMMON.REMD'
1421 include 'COMMON.SETUP'
1422 include 'COMMON.CHAIN'
1423 include 'COMMON.SBRIDGE'
1424 include 'COMMON.INTERACT'
1428 real xcoord(3,maxres2+2),prec
1429 real r_qfrag(50),r_qpair(100)
1430 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1431 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1432 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1433 & p_uscdiff(100*maxprocs)
1434 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1435 common /przechowalnia/ p_c
1437 call mpi_bcast(ii_write,1,mpi_integer,
1438 & king,CG_COMM,ierr)
1441 c print *,'traj1file',me,ii_write,ntwx_cache
1445 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1447 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1450 t5_restart1(1)=totT_cache(ii)
1451 t5_restart1(2)=EK_cache(ii)
1452 t5_restart1(3)=potE_cache(ii)
1453 t5_restart1(4)=t_bath_cache(ii)
1454 t5_restart1(5)=Uconst_cache(ii)
1455 call mpi_gather(t5_restart1,5,mpi_real,
1456 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1458 call mpi_gather(iset_cache(ii),1,mpi_integer,
1459 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1462 r_qfrag(i)=qfrag_cache(i,ii)
1465 r_qpair(i)=qpair_cache(i,ii)
1468 r_utheta(i)=utheta_cache(i,ii)
1469 r_ugamma(i)=ugamma_cache(i,ii)
1470 r_uscdiff(i)=uscdiff_cache(i,ii)
1473 call mpi_gather(r_qfrag,nfrag,mpi_real,
1474 & p_qfrag,nfrag,mpi_real,king,
1476 call mpi_gather(r_qpair,npair,mpi_real,
1477 & p_qpair,npair,mpi_real,king,
1479 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1480 & p_utheta,nfrag_back,mpi_real,king,
1482 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1483 & p_ugamma,nfrag_back,mpi_real,king,
1485 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1486 & p_uscdiff,nfrag_back,mpi_real,king,
1490 write (iout,*) "p_qfrag"
1492 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1494 write (iout,*) "p_qpair"
1496 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1502 r_c(j,i)=c_cache(j,i,ii)
1506 call mpi_gather(r_c,3*2*nres,mpi_real,
1507 & p_c,3*2*nres,mpi_real,king,
1513 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1514 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1515 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1516 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1517 call xdrfint_(ixdrf, nss, iret)
1520 call xdrfint(ixdrf, idssb(j)+nres, iret)
1521 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1523 call xdrfint_(ixdrf, ihpb(j), iret)
1524 call xdrfint_(ixdrf, jhpb(j), iret)
1527 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1528 call xdrfint_(ixdrf, iset_restart1(il), iret)
1530 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1533 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1536 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1537 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1538 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1543 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1548 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1552 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1556 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1557 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1558 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1559 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1560 call xdrfint(ixdrf, nss, iret)
1563 call xdrfint(ixdrf, idssb(j)+nres, iret)
1564 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1566 call xdrfint(ixdrf, ihpb(j), iret)
1567 call xdrfint(ixdrf, jhpb(j), iret)
1570 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1571 call xdrfint(ixdrf, iset_restart1(il), iret)
1573 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1576 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1579 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1580 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1581 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1586 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1591 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1595 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1601 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1603 if(me.eq.king) call xdrfclose(ixdrf, iret)
1605 do i=1,ntwx_cache-ii_write
1607 totT_cache(i)=totT_cache(ii_write+i)
1608 EK_cache(i)=EK_cache(ii_write+i)
1609 potE_cache(i)=potE_cache(ii_write+i)
1610 t_bath_cache(i)=t_bath_cache(ii_write+i)
1611 Uconst_cache(i)=Uconst_cache(ii_write+i)
1612 iset_cache(i)=iset_cache(ii_write+i)
1615 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1618 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1621 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1622 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1623 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1628 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1632 ntwx_cache=ntwx_cache-ii_write
1637 subroutine read1restart(i_index)
1638 implicit real*8 (a-h,o-z)
1639 include 'DIMENSIONS'
1642 include 'COMMON.IOUNITS'
1643 include 'COMMON.REMD'
1644 include 'COMMON.SETUP'
1645 include 'COMMON.CHAIN'
1646 include 'COMMON.SBRIDGE'
1647 include 'COMMON.INTERACT'
1648 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1651 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1652 common /przechowalnia/ d_restart1
1653 write (*,*) "Processor",me," called read1restart"
1656 open(irest2,file=mremd_rst_name,status='unknown')
1657 read(irest2,*,err=334) i
1658 write(iout,*) "Reading old rst in ASCI format"
1660 call read1restart_old
1664 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1667 call xdrfint_(ixdrf, i2rep(i), iret)
1670 call xdrfint_(ixdrf, ifirst(i), iret)
1673 call xdrfint_(ixdrf, nupa(0,il), iret)
1675 call xdrfint_(ixdrf, nupa(i,il), iret)
1678 call xdrfint_(ixdrf, ndowna(0,il), iret)
1680 call xdrfint_(ixdrf, ndowna(i,il), iret)
1685 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1689 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1692 call xdrfint(ixdrf, i2rep(i), iret)
1695 call xdrfint(ixdrf, ifirst(i), iret)
1698 call xdrfint(ixdrf, nupa(0,il), iret)
1700 call xdrfint(ixdrf, nupa(i,il), iret)
1703 call xdrfint(ixdrf, ndowna(0,il), iret)
1705 call xdrfint(ixdrf, ndowna(i,il), iret)
1710 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1715 call mpi_scatter(t_restart1,5,mpi_real,
1716 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1720 t_bath=t5_restart1(4)
1725 c read(irest2,'(3e15.5)')
1726 c & (d_restart1(j,i+2*nres*il),j=1,3)
1729 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1731 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1737 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1738 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1748 c read(irest2,'(3e15.5)')
1749 c & (d_restart1(j,i+2*nres*il),j=1,3)
1752 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1754 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1760 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1761 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1772 call xdrfint_(ixdrf, nset, iret)
1774 call xdrfint_(ixdrf,mset(i), iret)
1777 call xdrfint_(ixdrf,i2set(i), iret)
1783 call xdrfint_(ixdrf,itmp, iret)
1784 i_index(i,j,il,il1)=itmp
1792 call xdrfint(ixdrf, nset, iret)
1794 call xdrfint(ixdrf,mset(i), iret)
1797 call xdrfint(ixdrf,i2set(i), iret)
1803 call xdrfint(ixdrf,itmp, iret)
1804 i_index(i,j,il,il1)=itmp
1811 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1813 c call mpi_scatter(i2set,1,mpi_integer,
1814 c & iset,1,mpi_integer,king,
1816 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1823 if(me.eq.king) close(irest2)
1827 subroutine read1restart_old
1828 implicit real*8 (a-h,o-z)
1829 include 'DIMENSIONS'
1832 include 'COMMON.IOUNITS'
1833 include 'COMMON.REMD'
1834 include 'COMMON.SETUP'
1835 include 'COMMON.CHAIN'
1836 include 'COMMON.SBRIDGE'
1837 include 'COMMON.INTERACT'
1838 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1840 common /przechowalnia/ d_restart1
1842 open(irest2,file=mremd_rst_name,status='unknown')
1843 read (irest2,*) (i2rep(i),i=0,nodes-1)
1844 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1846 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1847 read (irest2,*) ndowna(0,il),
1848 & (ndowna(i,il),i=1,ndowna(0,il))
1851 read(irest2,*) (t_restart1(j,il),j=1,4)
1854 call mpi_scatter(t_restart1,5,mpi_real,
1855 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1859 t_bath=t5_restart1(4)
1864 read(irest2,'(3e15.5)')
1865 & (d_restart1(j,i+2*nres*il),j=1,3)
1869 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1870 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1880 read(irest2,'(3e15.5)')
1881 & (d_restart1(j,i+2*nres*il),j=1,3)
1885 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1886 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1892 if(me.eq.king) close(irest2)
1895 c------------------------------------------