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+4,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
277 stdfp=dsqrt(2*Rb*t_bath/d_time)
279 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
282 c print *,'irep',me,t_bath
284 if (me.eq.king .or. .not. out1file)
285 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
286 call rescale_weights(t_bath)
290 c------copy MD--------------
291 c The driver for molecular dynamics subroutines
292 c------------------------------------------------
298 if(me.eq.king.or..not.out1file)
299 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
305 c Determine the inverse of the inertia matrix.
306 call setup_MD_matrices
310 if (me.eq.king .or. .not. out1file)
311 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
312 stdfp=dsqrt(2*Rb*t_bath/d_time)
314 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
316 call rescale_weights(t_bath)
320 t_MDsetup = MPI_Wtime()-tt0
322 t_MDsetup = tcpu()-tt0
325 c Entering the MD loop
331 if (lang.eq.2 .or. lang.eq.3) then
335 call sd_verlet_p_setup
337 call sd_verlet_ciccotti_setup
341 pfric0_mat(i,j,0)=pfric_mat(i,j)
342 afric0_mat(i,j,0)=afric_mat(i,j)
343 vfric0_mat(i,j,0)=vfric_mat(i,j)
344 prand0_mat(i,j,0)=prand_mat(i,j)
345 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
346 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
351 flag_stoch(i)=.false.
355 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
357 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
361 else if (lang.eq.1 .or. lang.eq.4) then
365 if (me.eq.king .or. .not. out1file)
366 & write(iout,*) 'Setup time',time00-walltime
369 t_langsetup=MPI_Wtime()-tt0
372 t_langsetup=tcpu()-tt0
377 do while(.not.end_of_run)
379 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
380 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
382 if (lang.gt.0 .and. surfarea .and.
383 & mod(itime,reset_fricmat).eq.0) then
384 if (lang.eq.2 .or. lang.eq.3) then
388 call sd_verlet_p_setup
390 call sd_verlet_ciccotti_setup
394 pfric0_mat(i,j,0)=pfric_mat(i,j)
395 afric0_mat(i,j,0)=afric_mat(i,j)
396 vfric0_mat(i,j,0)=vfric_mat(i,j)
397 prand0_mat(i,j,0)=prand_mat(i,j)
398 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
399 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
404 flag_stoch(i)=.false.
407 else if (lang.eq.1 .or. lang.eq.4) then
410 write (iout,'(a,i10)')
411 & "Friction matrix reset based on surface area, itime",itime
413 if (reset_vel .and. tbf .and. lang.eq.0
414 & .and. mod(itime,count_reset_vel).eq.0) then
416 if (me.eq.king .or. .not. out1file)
417 & write(iout,'(a,f20.2)')
418 & "Velocities reset to random values, time",totT
421 d_t_old(j,i)=d_t(j,i)
425 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
429 d_t(j,0)=d_t(j,0)-vcm(j)
432 kinetic_T=2.0d0/(dimen3*Rb)*EK
433 scalfac=dsqrt(T_bath/kinetic_T)
434 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
437 d_t_old(j,i)=scalfac*d_t(j,i)
443 c Time-reversible RESPA algorithm
444 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
445 call RESPA_step(itime)
447 c Variable time step algorithm.
448 call velverlet_step(itime)
452 call brown_step(itime)
454 print *,"Brown dynamics not here!"
456 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
462 if (mod(itime,ntwe).eq.0) call statout(itime)
464 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
465 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
467 call hairpin(.true.,nharp,iharp)
468 call secondary2(.true.)
469 call pdbout(potE,tytul,ipdb)
474 if (mod(itime,ntwx).eq.0.and.traj1file) then
475 if(ntwx_cache.lt.max_cache_traj_use) then
476 ntwx_cache=ntwx_cache+1
478 if (max_cache_traj_use.ne.1)
479 & print *,itime,"processor ",me," over cache ",ntwx_cache
482 totT_cache(i)=totT_cache(i+1)
483 EK_cache(i)=EK_cache(i+1)
484 potE_cache(i)=potE_cache(i+1)
485 t_bath_cache(i)=t_bath_cache(i+1)
486 Uconst_cache(i)=Uconst_cache(i+1)
487 iset_cache(i)=iset_cache(i+1)
490 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
493 qpair_cache(ii,i)=qpair_cache(ii,i+1)
496 utheta_cache(ii,i)=utheta_cache(ii,i+1)
497 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
498 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
504 c_cache(j,ii,i)=c_cache(j,ii,i+1)
510 totT_cache(ntwx_cache)=totT
511 EK_cache(ntwx_cache)=EK
512 potE_cache(ntwx_cache)=potE
513 t_bath_cache(ntwx_cache)=t_bath
514 Uconst_cache(ntwx_cache)=Uconst
515 iset_cache(ntwx_cache)=iset
518 qfrag_cache(i,ntwx_cache)=qfrag(i)
521 qpair_cache(i,ntwx_cache)=qpair(i)
524 utheta_cache(i,ntwx_cache)=utheta(i)
525 ugamma_cache(i,ntwx_cache)=ugamma(i)
526 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
528 C print *,'przed returnbox'
530 C call enerprint(remd_ene(0,i))
533 c_cache(j,i,ntwx_cache)=c(j,i)
538 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
539 & .and..not.restart1file) then
542 open(irest1,file=mremd_rst_name,status='unknown')
543 write (irest1,*) "i2rep"
544 write (irest1,*) (i2rep(i),i=0,nodes-1)
545 write (irest1,*) "ifirst"
546 write (irest1,*) (ifirst(i),i=1,remd_m(1))
548 write (irest1,*) "nupa",il
549 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
550 write (irest1,*) "ndowna",il
551 write (irest1,*) ndowna(0,il),
552 & (ndowna(i,il),i=1,ndowna(0,il))
555 write (irest1,*) "nset"
556 write (irest1,*) nset
557 write (irest1,*) "mset"
558 write (irest1,*) (mset(i),i=1,nset)
559 write (irest1,*) "i2set"
560 write (irest1,*) (i2set(i),i=0,nodes-1)
561 write (irest1,*) "i_index"
565 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
573 open(irest2,file=rest2name,status='unknown')
574 write(irest2,*) totT,EK,potE,totE,t_bath
576 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
579 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
582 write (irest2,*) iset
589 c forced synchronization
590 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
591 & .and. .not. mremdsync) then
593 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
595 call mpi_recv(itime_master, 1, MPI_INTEGER,
596 & 0,101,CG_COMM, status, ierr)
597 call mpi_barrier(CG_COMM, ierr)
598 cdeb if (out1file.or.traj1file) then
599 cdeb call mpi_gather(itime,1,mpi_integer,
600 cdeb & icache_all,1,mpi_integer,king,
603 & call mpi_gather(ntwx_cache,1,mpi_integer,
604 & icache_all,1,mpi_integer,king,
607 & write(iout,*) 'REMD synchro at',itime_master,itime
608 if (itime_master.ge.n_timestep .or. ovrtim())
610 ctime call flush(iout)
615 if ((mod(itime,nstex).eq.0.and.me.eq.king
616 & .or.end_of_run.and.me.eq.king )
617 & .and. .not. mremdsync ) then
620 call mpi_isend(itime,1,MPI_INTEGER,i,101,
621 & CG_COMM, ireqi(i), ierr)
622 cd write(iout,*) 'REMD synchro with',i
625 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
626 call mpi_barrier(CG_COMM, ierr)
628 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
629 if (out1file.or.traj1file) then
630 cdeb call mpi_gather(itime,1,mpi_integer,
631 cdeb & itime_all,1,mpi_integer,king,
633 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
634 cdeb & (itime_all(i),i=1,nodes)
636 cdeb imin_itime=itime_all(1)
638 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
640 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
641 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
642 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
643 call mpi_gather(ntwx_cache,1,mpi_integer,
644 & icache_all,1,mpi_integer,king,
646 c write(iout,'(a19,8000i8)') ' ntwx_cache',
647 c & (icache_all(i),i=1,nodes)
648 ii_write=icache_all(1)
650 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
652 c write(iout,*) "MIN ii_write=",ii_write
655 ctime call flush(iout)
657 if(mremdsync .and. mod(itime,nstex).eq.0) then
659 if (me.eq.king .or. .not. out1file)
660 & write(iout,*) 'REMD synchro at',itime
663 call mpi_gather(ntwx_cache,1,mpi_integer,
664 & icache_all,1,mpi_integer,king,
667 write(iout,'(a19,8000i8)') ' ntwx_cache',
668 & (icache_all(i),i=1,nodes)
669 ii_write=icache_all(1)
671 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
673 write(iout,*) "MIN ii_write=",ii_write
679 c Update the time safety limiy
680 if (time001-time00.gt.safety) then
681 safety=time001-time00+600
682 write (iout,*) "****** SAFETY increased to",safety," s"
684 if (ovrtim()) end_of_run=.true.
686 if(synflag.and..not.end_of_run) then
690 write(iout,*) 'REMD before',me,t_bath
692 c call mpi_gather(t_bath,1,mpi_double_precision,
693 c & remd_t_bath,1,mpi_double_precision,king,
695 potEcomp(n_ene+1)=t_bath
697 potEcomp(n_ene+2)=iset
698 if (iset.lt.nset) then
702 potEcomp(n_ene+3)=Uconst
709 potEcomp(n_ene+4)=Uconst
713 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
714 & remd_ene(0,1),n_ene+5,mpi_double_precision,king,
717 call mpi_gather(elow,1,mpi_double_precision,
718 & elowi,1,mpi_double_precision,king,
720 call mpi_gather(ehigh,1,mpi_double_precision,
721 & ehighi,1,mpi_double_precision,king,
726 if (me.eq.king .or. .not. out1file) then
727 write(iout,*) 'REMD gather times=',time03-time01
731 if (restart1file) call write1rst(i_index)
734 if (me.eq.king .or. .not. out1file) then
735 write(iout,*) 'REMD writing rst time=',time04-time03
738 if (traj1file) call write1traj
740 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
741 cdeb & icache_all,1,mpi_integer,king,
743 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
744 cdeb & (icache_all(i),i=1,nodes)
749 if (me.eq.king .or. .not. out1file) then
750 write(iout,*) 'REMD writing traj time=',time05-time04
757 remd_t_bath(i)=remd_ene(n_ene+1,i)
758 iremd_iset(i)=remd_ene(n_ene+2,i)
761 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
763 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
767 write(iout,*) 'REMD exchange temp,ene'
769 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
770 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
773 c-------------------------------------
775 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
778 write (iout,*) "remd_m(1)",remd_m(1)
780 i=ifirst(iran_num(1,remd_m(1)))
786 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
787 if(i.gt.0.and.nupa(0,i).gt.0) then
789 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
791 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
793 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
795 c do while (iex.eq.i)
796 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
797 iex=nupa(iran_num(1,int(nupa(0,i))),i)
799 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
801 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
803 c Swap temperatures between conformations i and iex with recalculating the free energies
804 c following temperature changes.
805 ene_iex_iex=remd_ene(0,iex)
806 ene_i_i=remd_ene(0,i)
807 c write (iout,*) "i",i," ene_i_i",ene_i_i,
808 c & " iex",iex," ene_iex_iex",ene_iex_iex
809 c write (iout,*) "rescaling weights with temperature",
812 call rescale_weights(remd_t_bath(i))
814 c write (iout,*) "0,iex",remd_t_bath(i)
815 c call enerprint(remd_ene(0,iex))
817 call sum_energy(remd_ene(0,iex),.false.)
818 ene_iex_i=remd_ene(0,iex)
819 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
821 c write (iout,*) "0,i",remd_t_bath(i)
822 c call enerprint(remd_ene(0,i))
824 call sum_energy(remd_ene(0,i),.false.)
825 c write (iout,*) "ene_i_i",remd_ene(0,i)
827 c write (iout,*) "rescaling weights with temperature",
829 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
830 write (iout,*) "ERROR: inconsistent energies:",i,
831 & ene_i_i,remd_ene(0,i)
833 call rescale_weights(remd_t_bath(iex))
835 c write (iout,*) "0,i",remd_t_bath(iex)
836 call enerprint(remd_ene(0,i))
838 call sum_energy(remd_ene(0,i),.false.)
839 c write (iout,*) "ene_i_iex",remd_ene(0,i)
841 ene_i_iex=remd_ene(0,i)
843 c write (iout,*) "0,iex",remd_t_bath(iex)
844 c call enerprint(remd_ene(0,iex))
846 call sum_energy(remd_ene(0,iex),.false.)
847 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
848 write (iout,*) "ERROR: inconsistent energies:",iex,
849 & ene_iex_iex,remd_ene(0,iex)
851 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
852 c write (iout,*) "i",i," iex",iex
853 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
854 c & " ene_i_iex",ene_i_iex,
855 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
857 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
858 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
860 c write(iout,*) 'delta',delta
861 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
862 c & (remd_ene(i)-remd_ene(iex))/Rb/
863 c & (remd_t_bath(i)*remd_t_bath(iex))
865 if (delta .gt. 50.0d0) then
871 else if (delta.lt.-50.0d0) then
880 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
881 xxx=ran_number(0.0d0,1.0d0)
882 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
884 if (delta .gt. xxx) then
886 remd_t_bath(i)=remd_t_bath(iex)
888 remd_ene(0,i)=ene_i_iex
889 remd_ene(0,iex)=ene_iex_i
895 ehighi(i)=ehighi(iex)
902 nupa(k,i)=nupa(k,iex)
905 ndowna(k,i)=ndowna(k,iex)
909 if (ifirst(il).eq.i) ifirst(il)=iex
911 if (nupa(k,il).eq.i) then
913 elseif (nupa(k,il).eq.iex) then
918 if (ndowna(k,il).eq.i) then
920 elseif (ndowna(k,il).eq.iex) then
926 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
928 i2rep(i-1)=i2rep(iex-1)
931 c write(iout,*) 'exchange',i,iex
932 c write (iout,'(a8,100i4)') "@ ifirst",
933 c & (ifirst(k),k=1,remd_m(1))
935 c write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
936 c & (nupa(k,il),k=1,nupa(0,il))
937 c write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
938 c & (ndowna(k,il),k=1,ndowna(0,il))
943 remd_ene(0,iex)=ene_iex_iex
944 remd_ene(0,i)=ene_i_i
950 cd write (iout,*) "exchange completed"
954 cd write(iout,*) "########",ii
956 i_temp=iran_num(1,nrep)
957 i_mult=iran_num(1,remd_m(i_temp))
958 i_iset=iran_num(1,nset)
959 i_mset=iran_num(1,mset(i_iset))
960 i=i_index(i_temp,i_mult,i_iset,i_mset)
962 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
965 cd write(iout,*) "i_dir=",i_dir
967 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
970 i_mult1=iran_num(1,remd_m(i_temp1))
972 i_mset1=iran_num(1,mset(i_iset1))
973 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
975 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
978 i_mult1=iran_num(1,remd_m(i_temp1))
980 i_mset1=iran_num(1,mset(i_iset1))
981 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
982 econstr_temp_i=remd_ene(20,i)
983 econstr_temp_iex=remd_ene(20,iex)
984 remd_ene(20,i)=remd_ene(n_ene+3,i)
985 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
987 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
990 i_mult1=iran_num(1,remd_m(i_temp1))
992 i_mset1=iran_num(1,mset(i_iset1))
993 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
994 econstr_temp_i=remd_ene(20,i)
995 econstr_temp_iex=remd_ene(20,iex)
996 remd_ene(20,i)=remd_ene(n_ene+3,i)
997 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1003 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1006 c Swap temperatures between conformations i and iex with recalculating the free energies
1007 c following temperature changes.
1008 ene_iex_iex=remd_ene(0,iex)
1009 ene_i_i=remd_ene(0,i)
1010 co write (iout,*) "rescaling weights with temperature",
1012 call rescale_weights(remd_t_bath(i))
1014 call sum_energy(remd_ene(0,iex),.false.)
1015 ene_iex_i=remd_ene(0,iex)
1016 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
1017 c call sum_energy(remd_ene(0,i),.false.)
1018 cd write (iout,*) "ene_i_i",remd_ene(0,i)
1019 c write (iout,*) "rescaling weights with temperature",
1020 c & remd_t_bath(iex)
1021 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1022 c write (iout,*) "ERROR: inconsistent energies:",i,
1023 c & ene_i_i,remd_ene(0,i)
1025 call rescale_weights(remd_t_bath(iex))
1026 call sum_energy(remd_ene(0,i),.false.)
1027 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1028 ene_i_iex=remd_ene(0,i)
1029 c call sum_energy(remd_ene(0,iex),.false.)
1030 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1031 c write (iout,*) "ERROR: inconsistent energies:",iex,
1032 c & ene_iex_iex,remd_ene(0,iex)
1034 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1035 c write (iout,*) "i",i," iex",iex
1036 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1037 cd & " ene_i_iex",ene_i_iex,
1038 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1039 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1040 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1042 cd write(iout,*) 'delta',delta
1043 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1044 c & (remd_ene(i)-remd_ene(iex))/Rb/
1045 c & (remd_t_bath(i)*remd_t_bath(iex))
1046 if (delta .gt. 50.0d0) then
1051 if (i_dir.eq.1.or.i_dir.eq.3)
1052 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1053 if (i_dir.eq.2.or.i_dir.eq.3)
1054 & iremd_tot_usa(int(i2set(i-1)))=
1055 & iremd_tot_usa(int(i2set(i-1)))+1
1056 xxx=ran_number(0.0d0,1.0d0)
1057 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1058 if (delta .gt. xxx) then
1060 remd_t_bath(i)=remd_t_bath(iex)
1061 remd_t_bath(iex)=tmp
1064 iremd_iset(i)=iremd_iset(iex)
1065 iremd_iset(iex)=itmp
1067 remd_ene(0,i)=ene_i_iex
1068 remd_ene(0,iex)=ene_iex_i
1070 if (i_dir.eq.1.or.i_dir.eq.3)
1071 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1074 i2rep(i-1)=i2rep(iex-1)
1077 if (i_dir.eq.2.or.i_dir.eq.3)
1078 & iremd_acc_usa(int(i2set(i-1)))=
1079 & iremd_acc_usa(int(i2set(i-1)))+1
1082 i2set(i-1)=i2set(iex-1)
1085 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1086 i_index(i_temp,i_mult,i_iset,i_mset)=
1087 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1088 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1091 remd_ene(0,iex)=ene_iex_iex
1092 remd_ene(0,i)=ene_i_i
1093 remd_ene(20,iex)=econstr_temp_iex
1094 remd_ene(20,i)=econstr_temp_i
1098 cd do il1=1,mset(il)
1101 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1114 c-------------------------------------
1115 write (iout,*) "NREP",nrep
1117 if(iremd_tot(i).ne.0)
1118 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1119 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1124 if(iremd_tot_usa(i).ne.0)
1125 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1126 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1132 cd write (iout,'(a6,100i4)') "ifirst",
1133 cd & (ifirst(i),i=1,remd_m(1))
1135 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1136 cd & (nupa(i,il),i=1,nupa(0,il))
1137 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1138 cd & (ndowna(i,il),i=1,ndowna(0,il))
1143 cd write (iout,*) "Before scatter"
1145 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1146 & t_bath,1,mpi_double_precision,king,
1148 cd write (iout,*) "After scatter"
1151 & call mpi_scatter(iremd_iset,1,mpi_integer,
1152 & iset,1,mpi_integer,king,
1156 if (me.eq.king .or. .not. out1file) then
1157 write(iout,*) 'REMD scatter time=',time07-time06
1161 call mpi_scatter(elowi,1,mpi_double_precision,
1162 & elow,1,mpi_double_precision,king,
1164 call mpi_scatter(ehighi,1,mpi_double_precision,
1165 & ehigh,1,mpi_double_precision,king,
1168 call rescale_weights(t_bath)
1169 co write (iout,*) "Processor",me,
1170 co & " rescaling weights with temperature",t_bath
1172 stdfp=dsqrt(2*Rb*t_bath/d_time)
1174 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1177 cde write(iout,*) 'REMD after',me,t_bath
1179 if (me.eq.king .or. .not. out1file) then
1180 write(iout,*) 'REMD exchange time=',time08-time00
1186 if (restart1file) then
1187 if (me.eq.king .or. .not. out1file)
1188 & write(iout,*) 'writing restart at the end of run'
1189 call write1rst(i_index)
1192 if (traj1file) call write1traj
1194 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1195 cdeb & icache_all,1,mpi_integer,king,
1196 cdeb & CG_COMM,ierr)
1197 cdeb write(iout,'(a40,8000i8)')
1198 cdeb & ' ntwx_cache after traj1file at the end',
1199 cdeb & (icache_all(i),i=1,nodes)
1204 t_MD=MPI_Wtime()-tt0
1208 if (me.eq.king .or. .not. out1file) then
1209 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1211 & 'MD calculations setup:',t_MDsetup,
1212 & 'Energy & gradient evaluation:',t_enegrad,
1213 & 'Stochastic MD setup:',t_langsetup,
1214 & 'Stochastic MD step setup:',t_sdsetup,
1216 write (iout,'(/28(1h=),a25,27(1h=))')
1217 & ' End of MD calculation '
1222 c-----------------------------------------------------------------------
1223 subroutine write1rst(i_index)
1224 implicit real*8 (a-h,o-z)
1225 include 'DIMENSIONS'
1228 include 'COMMON.IOUNITS'
1229 include 'COMMON.REMD'
1230 include 'COMMON.SETUP'
1231 include 'COMMON.CHAIN'
1232 include 'COMMON.SBRIDGE'
1233 include 'COMMON.INTERACT'
1235 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1236 & d_restart2(3,2*maxres*maxprocs)
1240 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1241 common /przechowalnia/ d_restart1,d_restart2
1246 t5_restart1(4)=t_bath
1247 t5_restart1(5)=Uconst
1249 call mpi_gather(t5_restart1,5,mpi_real,
1250 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1258 call mpi_gather(r_d,3*2*nres,mpi_real,
1259 & d_restart1,3*2*nres,mpi_real,king,
1268 call mpi_gather(r_d,3*2*nres,mpi_real,
1269 & d_restart2,3*2*nres,mpi_real,king,
1274 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1276 call xdrfint_(ixdrf, i2rep(i), iret)
1279 call xdrfint_(ixdrf, ifirst(i), iret)
1283 call xdrfint_(ixdrf, nupa(i,il), iret)
1287 call xdrfint_(ixdrf, ndowna(i,il), iret)
1293 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1300 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1307 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1313 call xdrfint_(ixdrf, nset, iret)
1315 call xdrfint_(ixdrf,mset(i), iret)
1318 call xdrfint_(ixdrf,i2set(i), iret)
1324 itmp=i_index(i,j,il,il1)
1325 call xdrfint_(ixdrf,itmp, iret)
1332 call xdrfclose_(ixdrf, iret)
1334 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1336 call xdrfint(ixdrf, i2rep(i), iret)
1339 call xdrfint(ixdrf, ifirst(i), iret)
1343 call xdrfint(ixdrf, nupa(i,il), iret)
1347 call xdrfint(ixdrf, ndowna(i,il), iret)
1353 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1360 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1367 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1374 call xdrfint(ixdrf, nset, iret)
1376 call xdrfint(ixdrf,mset(i), iret)
1379 call xdrfint(ixdrf,i2set(i), iret)
1385 itmp=i_index(i,j,il,il1)
1386 call xdrfint(ixdrf,itmp, iret)
1393 call xdrfclose(ixdrf, iret)
1400 subroutine write1traj
1401 implicit real*8 (a-h,o-z)
1402 include 'DIMENSIONS'
1405 include 'COMMON.IOUNITS'
1406 include 'COMMON.REMD'
1407 include 'COMMON.SETUP'
1408 include 'COMMON.CHAIN'
1409 include 'COMMON.SBRIDGE'
1410 include 'COMMON.INTERACT'
1414 real xcoord(3,maxres2+2),prec
1415 real r_qfrag(50),r_qpair(100)
1416 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1417 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1418 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1419 & p_uscdiff(100*maxprocs)
1420 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1421 common /przechowalnia/ p_c
1423 call mpi_bcast(ii_write,1,mpi_integer,
1424 & king,CG_COMM,ierr)
1427 print *,'traj1file',me,ii_write,ntwx_cache
1431 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1433 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1436 t5_restart1(1)=totT_cache(ii)
1437 t5_restart1(2)=EK_cache(ii)
1438 t5_restart1(3)=potE_cache(ii)
1439 t5_restart1(4)=t_bath_cache(ii)
1440 t5_restart1(5)=Uconst_cache(ii)
1441 call mpi_gather(t5_restart1,5,mpi_real,
1442 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1444 call mpi_gather(iset_cache(ii),1,mpi_integer,
1445 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1448 r_qfrag(i)=qfrag_cache(i,ii)
1451 r_qpair(i)=qpair_cache(i,ii)
1454 r_utheta(i)=utheta_cache(i,ii)
1455 r_ugamma(i)=ugamma_cache(i,ii)
1456 r_uscdiff(i)=uscdiff_cache(i,ii)
1459 call mpi_gather(r_qfrag,nfrag,mpi_real,
1460 & p_qfrag,nfrag,mpi_real,king,
1462 call mpi_gather(r_qpair,npair,mpi_real,
1463 & p_qpair,npair,mpi_real,king,
1465 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1466 & p_utheta,nfrag_back,mpi_real,king,
1468 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1469 & p_ugamma,nfrag_back,mpi_real,king,
1471 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1472 & p_uscdiff,nfrag_back,mpi_real,king,
1476 write (iout,*) "p_qfrag"
1478 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1480 write (iout,*) "p_qpair"
1482 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1488 r_c(j,i)=c_cache(j,i,ii)
1492 call mpi_gather(r_c,3*2*nres,mpi_real,
1493 & p_c,3*2*nres,mpi_real,king,
1499 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1500 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1501 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1502 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1503 call xdrfint_(ixdrf, nss, iret)
1505 call xdrfint_(ixdrf, ihpb(j), iret)
1506 call xdrfint_(ixdrf, jhpb(j), iret)
1508 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1509 call xdrfint_(ixdrf, iset_restart1(il), iret)
1511 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1514 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1517 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1518 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1519 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1524 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1529 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1533 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1537 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1538 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1539 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1540 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1541 call xdrfint(ixdrf, nss, iret)
1543 call xdrfint(ixdrf, ihpb(j), iret)
1544 call xdrfint(ixdrf, jhpb(j), iret)
1546 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1547 call xdrfint(ixdrf, iset_restart1(il), iret)
1549 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1552 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1555 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1556 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1557 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1562 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1567 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1571 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1577 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1579 if(me.eq.king) call xdrfclose(ixdrf, iret)
1581 do i=1,ntwx_cache-ii_write
1583 totT_cache(i)=totT_cache(ii_write+i)
1584 EK_cache(i)=EK_cache(ii_write+i)
1585 potE_cache(i)=potE_cache(ii_write+i)
1586 t_bath_cache(i)=t_bath_cache(ii_write+i)
1587 Uconst_cache(i)=Uconst_cache(ii_write+i)
1588 iset_cache(i)=iset_cache(ii_write+i)
1591 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1594 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1597 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1598 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1599 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1604 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1608 ntwx_cache=ntwx_cache-ii_write
1613 subroutine read1restart(i_index)
1614 implicit real*8 (a-h,o-z)
1615 include 'DIMENSIONS'
1618 include 'COMMON.IOUNITS'
1619 include 'COMMON.REMD'
1620 include 'COMMON.SETUP'
1621 include 'COMMON.CHAIN'
1622 include 'COMMON.SBRIDGE'
1623 include 'COMMON.INTERACT'
1624 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1627 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1628 common /przechowalnia/ d_restart1
1629 write (*,*) "Processor",me," called read1restart"
1632 open(irest2,file=mremd_rst_name,status='unknown')
1633 read(irest2,*,err=334) i
1634 write(iout,*) "Reading old rst in ASCI format"
1636 call read1restart_old
1640 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1643 call xdrfint_(ixdrf, i2rep(i), iret)
1646 call xdrfint_(ixdrf, ifirst(i), iret)
1649 call xdrfint_(ixdrf, nupa(0,il), iret)
1651 call xdrfint_(ixdrf, nupa(i,il), iret)
1654 call xdrfint_(ixdrf, ndowna(0,il), iret)
1656 call xdrfint_(ixdrf, ndowna(i,il), iret)
1661 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1665 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1668 call xdrfint(ixdrf, i2rep(i), iret)
1671 call xdrfint(ixdrf, ifirst(i), iret)
1674 call xdrfint(ixdrf, nupa(0,il), iret)
1676 call xdrfint(ixdrf, nupa(i,il), iret)
1679 call xdrfint(ixdrf, ndowna(0,il), iret)
1681 call xdrfint(ixdrf, ndowna(i,il), iret)
1686 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1691 call mpi_scatter(t_restart1,5,mpi_real,
1692 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1696 t_bath=t5_restart1(4)
1701 c read(irest2,'(3e15.5)')
1702 c & (d_restart1(j,i+2*nres*il),j=1,3)
1705 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1707 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1713 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1714 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1724 c read(irest2,'(3e15.5)')
1725 c & (d_restart1(j,i+2*nres*il),j=1,3)
1728 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1730 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1736 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1737 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1748 call xdrfint_(ixdrf, nset, iret)
1750 call xdrfint_(ixdrf,mset(i), iret)
1753 call xdrfint_(ixdrf,i2set(i), iret)
1759 call xdrfint_(ixdrf,itmp, iret)
1760 i_index(i,j,il,il1)=itmp
1768 call xdrfint(ixdrf, nset, iret)
1770 call xdrfint(ixdrf,mset(i), iret)
1773 call xdrfint(ixdrf,i2set(i), iret)
1779 call xdrfint(ixdrf,itmp, iret)
1780 i_index(i,j,il,il1)=itmp
1787 call mpi_scatter(i2set,1,mpi_integer,
1788 & iset,1,mpi_integer,king,
1794 if(me.eq.king) close(irest2)
1798 subroutine read1restart_old
1799 implicit real*8 (a-h,o-z)
1800 include 'DIMENSIONS'
1803 include 'COMMON.IOUNITS'
1804 include 'COMMON.REMD'
1805 include 'COMMON.SETUP'
1806 include 'COMMON.CHAIN'
1807 include 'COMMON.SBRIDGE'
1808 include 'COMMON.INTERACT'
1809 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1811 common /przechowalnia/ d_restart1
1813 open(irest2,file=mremd_rst_name,status='unknown')
1814 read (irest2,*) (i2rep(i),i=0,nodes-1)
1815 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1817 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1818 read (irest2,*) ndowna(0,il),
1819 & (ndowna(i,il),i=1,ndowna(0,il))
1822 read(irest2,*) (t_restart1(j,il),j=1,4)
1825 call mpi_scatter(t_restart1,5,mpi_real,
1826 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1830 t_bath=t5_restart1(4)
1835 read(irest2,'(3e15.5)')
1836 & (d_restart1(j,i+2*nres*il),j=1,3)
1840 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1841 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1851 read(irest2,'(3e15.5)')
1852 & (d_restart1(j,i+2*nres*il),j=1,3)
1856 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1857 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1863 if(me.eq.king) close(irest2)
1866 c------------------------------------------
1867 subroutine returnbox
1868 include 'DIMENSIONS'
1870 include 'COMMON.CONTROL'
1871 include 'COMMON.VAR'
1874 include 'COMMON.LANGEVIN'
1876 include 'COMMON.LANGEVIN.lang0'
1878 include 'COMMON.CHAIN'
1879 include 'COMMON.DERIV'
1880 include 'COMMON.GEO'
1881 include 'COMMON.LOCAL'
1882 include 'COMMON.INTERACT'
1883 include 'COMMON.IOUNITS'
1884 include 'COMMON.NAMES'
1885 include 'COMMON.TIME1'
1886 include 'COMMON.REMD'
1887 include 'COMMON.SETUP'
1888 include 'COMMON.MUCA'
1889 include 'COMMON.HAIRPIN'
1893 C write(*,*) 'initial', i,j,c(j,i)
1896 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
1898 if (allareout.eq.1) then
1899 ireturnval=int(c(j,i)/boxxsize)
1900 if (c(j,i).le.0) ireturnval=ireturnval-1
1901 do k=chain_beg,chain_end
1902 c(j,k)=c(j,k)-ireturnval*boxxsize
1903 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
1909 if (int(c(j,i)/boxxsize).eq.0) allareout=0
1912 if (allareout.eq.1) then
1913 ireturnval=int(c(j,i)/boxxsize)
1914 if (c(j,i).le.0) ireturnval=ireturnval-1
1916 c(j,k)=c(j,k)-ireturnval*boxxsize
1917 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
1922 C write(*,*) 'befor no jump', i,j,c(j,i)
1926 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1927 difference=abs(c(j,i-1)-c(j,i))
1928 C print *,'diff', difference
1929 if (difference.gt.boxxsize/2.0) then
1930 if (c(j,i-1).gt.c(j,i)) then
1939 c(j,i)=c(j,i)+nojumpval*boxxsize
1940 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
1944 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1945 difference=abs(c(j,i-1)-c(j,i))
1946 if (difference.gt.boxxsize/2.0) then
1947 if (c(j,i-1).gt.c(j,i)) then
1956 c(j,i)=c(j,i)+nojumpval*boxxsize
1957 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
1961 C write(*,*) 'after no jump', i,j,c(j,i)
1968 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
1970 if (allareout.eq.1) then
1971 ireturnval=int(c(j,i)/boxysize)
1972 if (c(j,i).le.0) ireturnval=ireturnval-1
1973 do k=chain_beg,chain_end
1974 c(j,k)=c(j,k)-ireturnval*boxysize
1975 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
1981 if (int(c(j,i)/boxysize).eq.0) allareout=0
1984 if (allareout.eq.1) then
1985 ireturnval=int(c(j,i)/boxysize)
1986 if (c(j,i).le.0) ireturnval=ireturnval-1
1988 c(j,k)=c(j,k)-ireturnval*boxysize
1989 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
1994 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1995 difference=abs(c(j,i-1)-c(j,i))
1996 if (difference.gt.boxysize/2.0) then
1997 if (c(j,i-1).gt.c(j,i)) then
2006 c(j,i)=c(j,i)+nojumpval*boxysize
2007 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
2011 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2012 difference=abs(c(j,i-1)-c(j,i))
2013 if (difference.gt.boxysize/2.0) then
2014 if (c(j,i-1).gt.c(j,i)) then
2023 c(j,i)=c(j,i)+nojumpval*boxysize
2024 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
2030 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
2032 if (allareout.eq.1) then
2033 ireturnval=int(c(j,i)/boxysize)
2034 if (c(j,i).le.0) ireturnval=ireturnval-1
2035 do k=chain_beg,chain_end
2036 c(j,k)=c(j,k)-ireturnval*boxzsize
2037 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
2043 if (int(c(j,i)/boxzsize).eq.0) allareout=0
2046 if (allareout.eq.1) then
2047 ireturnval=int(c(j,i)/boxzsize)
2048 if (c(j,i).le.0) ireturnval=ireturnval-1
2050 c(j,k)=c(j,k)-ireturnval*boxzsize
2051 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
2056 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2057 difference=abs(c(j,i-1)-c(j,i))
2058 if (difference.gt.(boxzsize/2.0)) then
2059 if (c(j,i-1).gt.c(j,i)) then
2068 c(j,i)=c(j,i)+nojumpval*boxzsize
2069 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
2073 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2074 difference=abs(c(j,i-1)-c(j,i))
2075 if (difference.gt.boxzsize/2.0) then
2076 if (c(j,i-1).gt.c(j,i)) then
2085 c(j,i)=c(j,i)+nojumpval*boxzsize
2086 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize