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'
532 c_cache(j,i,ntwx_cache)=c(j,i)
537 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
538 & .and..not.restart1file) then
541 open(irest1,file=mremd_rst_name,status='unknown')
542 write (irest1,*) "i2rep"
543 write (irest1,*) (i2rep(i),i=0,nodes-1)
544 write (irest1,*) "ifirst"
545 write (irest1,*) (ifirst(i),i=1,remd_m(1))
547 write (irest1,*) "nupa",il
548 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
549 write (irest1,*) "ndowna",il
550 write (irest1,*) ndowna(0,il),
551 & (ndowna(i,il),i=1,ndowna(0,il))
554 write (irest1,*) "nset"
555 write (irest1,*) nset
556 write (irest1,*) "mset"
557 write (irest1,*) (mset(i),i=1,nset)
558 write (irest1,*) "i2set"
559 write (irest1,*) (i2set(i),i=0,nodes-1)
560 write (irest1,*) "i_index"
564 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
572 open(irest2,file=rest2name,status='unknown')
573 write(irest2,*) totT,EK,potE,totE,t_bath
575 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
578 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
581 write (irest2,*) iset
588 c forced synchronization
589 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
590 & .and. .not. mremdsync) then
592 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
594 call mpi_recv(itime_master, 1, MPI_INTEGER,
595 & 0,101,CG_COMM, status, ierr)
596 call mpi_barrier(CG_COMM, ierr)
597 cdeb if (out1file.or.traj1file) then
598 cdeb call mpi_gather(itime,1,mpi_integer,
599 cdeb & icache_all,1,mpi_integer,king,
602 & call mpi_gather(ntwx_cache,1,mpi_integer,
603 & icache_all,1,mpi_integer,king,
606 & write(iout,*) 'REMD synchro at',itime_master,itime
607 if (itime_master.ge.n_timestep .or. ovrtim())
609 ctime call flush(iout)
614 if ((mod(itime,nstex).eq.0.and.me.eq.king
615 & .or.end_of_run.and.me.eq.king )
616 & .and. .not. mremdsync ) then
619 call mpi_isend(itime,1,MPI_INTEGER,i,101,
620 & CG_COMM, ireqi(i), ierr)
621 cd write(iout,*) 'REMD synchro with',i
624 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
625 call mpi_barrier(CG_COMM, ierr)
627 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
628 if (out1file.or.traj1file) then
629 cdeb call mpi_gather(itime,1,mpi_integer,
630 cdeb & itime_all,1,mpi_integer,king,
632 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
633 cdeb & (itime_all(i),i=1,nodes)
635 cdeb imin_itime=itime_all(1)
637 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
639 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
640 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
641 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
642 call mpi_gather(ntwx_cache,1,mpi_integer,
643 & icache_all,1,mpi_integer,king,
645 c write(iout,'(a19,8000i8)') ' ntwx_cache',
646 c & (icache_all(i),i=1,nodes)
647 ii_write=icache_all(1)
649 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
651 c write(iout,*) "MIN ii_write=",ii_write
654 ctime call flush(iout)
656 if(mremdsync .and. mod(itime,nstex).eq.0) then
658 if (me.eq.king .or. .not. out1file)
659 & write(iout,*) 'REMD synchro at',itime
662 call mpi_gather(ntwx_cache,1,mpi_integer,
663 & icache_all,1,mpi_integer,king,
666 write(iout,'(a19,8000i8)') ' ntwx_cache',
667 & (icache_all(i),i=1,nodes)
668 ii_write=icache_all(1)
670 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
672 write(iout,*) "MIN ii_write=",ii_write
678 c Update the time safety limiy
679 if (time001-time00.gt.safety) then
680 safety=time001-time00+600
681 write (iout,*) "****** SAFETY increased to",safety," s"
683 if (ovrtim()) end_of_run=.true.
685 if(synflag.and..not.end_of_run) then
689 write(iout,*) 'REMD before',me,t_bath
691 c call mpi_gather(t_bath,1,mpi_double_precision,
692 c & remd_t_bath,1,mpi_double_precision,king,
694 potEcomp(n_ene+1)=t_bath
696 potEcomp(n_ene+2)=iset
697 if (iset.lt.nset) then
701 potEcomp(n_ene+3)=Uconst
708 potEcomp(n_ene+4)=Uconst
712 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
713 & remd_ene(0,1),n_ene+5,mpi_double_precision,king,
716 call mpi_gather(elow,1,mpi_double_precision,
717 & elowi,1,mpi_double_precision,king,
719 call mpi_gather(ehigh,1,mpi_double_precision,
720 & ehighi,1,mpi_double_precision,king,
725 if (me.eq.king .or. .not. out1file) then
726 write(iout,*) 'REMD gather times=',time03-time01
730 if (restart1file) call write1rst(i_index)
733 if (me.eq.king .or. .not. out1file) then
734 write(iout,*) 'REMD writing rst time=',time04-time03
737 if (traj1file) call write1traj
739 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
740 cdeb & icache_all,1,mpi_integer,king,
742 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
743 cdeb & (icache_all(i),i=1,nodes)
748 if (me.eq.king .or. .not. out1file) then
749 write(iout,*) 'REMD writing traj time=',time05-time04
756 remd_t_bath(i)=remd_ene(n_ene+1,i)
757 iremd_iset(i)=remd_ene(n_ene+2,i)
760 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
762 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
766 write(iout,*) 'REMD exchange temp,ene'
768 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
769 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
772 c-------------------------------------
774 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
777 write (iout,*) "remd_m(1)",remd_m(1)
779 i=ifirst(iran_num(1,remd_m(1)))
785 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
786 if(i.gt.0.and.nupa(0,i).gt.0) then
788 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
790 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
792 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
794 c do while (iex.eq.i)
795 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
796 iex=nupa(iran_num(1,int(nupa(0,i))),i)
798 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
800 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
802 c Swap temperatures between conformations i and iex with recalculating the free energies
803 c following temperature changes.
804 ene_iex_iex=remd_ene(0,iex)
805 ene_i_i=remd_ene(0,i)
806 c write (iout,*) "i",i," ene_i_i",ene_i_i,
807 c & " iex",iex," ene_iex_iex",ene_iex_iex
808 c write (iout,*) "rescaling weights with temperature",
811 call rescale_weights(remd_t_bath(i))
813 c write (iout,*) "0,iex",remd_t_bath(i)
814 c call enerprint(remd_ene(0,iex))
816 call sum_energy(remd_ene(0,iex),.false.)
817 ene_iex_i=remd_ene(0,iex)
818 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
820 c write (iout,*) "0,i",remd_t_bath(i)
821 c call enerprint(remd_ene(0,i))
823 call sum_energy(remd_ene(0,i),.false.)
824 c write (iout,*) "ene_i_i",remd_ene(0,i)
826 c write (iout,*) "rescaling weights with temperature",
828 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
829 write (iout,*) "ERROR: inconsistent energies:",i,
830 & ene_i_i,remd_ene(0,i)
832 call rescale_weights(remd_t_bath(iex))
834 c write (iout,*) "0,i",remd_t_bath(iex)
835 c call enerprint(remd_ene(0,i))
837 call sum_energy(remd_ene(0,i),.false.)
838 c write (iout,*) "ene_i_iex",remd_ene(0,i)
840 ene_i_iex=remd_ene(0,i)
842 c write (iout,*) "0,iex",remd_t_bath(iex)
843 c call enerprint(remd_ene(0,iex))
845 call sum_energy(remd_ene(0,iex),.false.)
846 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
847 write (iout,*) "ERROR: inconsistent energies:",iex,
848 & ene_iex_iex,remd_ene(0,iex)
850 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
851 c write (iout,*) "i",i," iex",iex
852 c write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
853 c & " ene_i_iex",ene_i_iex,
854 c & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
856 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
857 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
859 c write(iout,*) 'delta',delta
860 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
861 c & (remd_ene(i)-remd_ene(iex))/Rb/
862 c & (remd_t_bath(i)*remd_t_bath(iex))
864 if (delta .gt. 50.0d0) then
870 else if (delta.lt.-50.0d0) then
879 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
880 xxx=ran_number(0.0d0,1.0d0)
881 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
883 if (delta .gt. xxx) then
885 remd_t_bath(i)=remd_t_bath(iex)
887 remd_ene(0,i)=ene_i_iex
888 remd_ene(0,iex)=ene_iex_i
894 ehighi(i)=ehighi(iex)
901 nupa(k,i)=nupa(k,iex)
904 ndowna(k,i)=ndowna(k,iex)
908 if (ifirst(il).eq.i) ifirst(il)=iex
910 if (nupa(k,il).eq.i) then
912 elseif (nupa(k,il).eq.iex) then
917 if (ndowna(k,il).eq.i) then
919 elseif (ndowna(k,il).eq.iex) then
925 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
927 i2rep(i-1)=i2rep(iex-1)
930 c write(iout,*) 'exchange',i,iex
931 c write (iout,'(a8,100i4)') "@ ifirst",
932 c & (ifirst(k),k=1,remd_m(1))
934 c write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
935 c & (nupa(k,il),k=1,nupa(0,il))
936 c write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
937 c & (ndowna(k,il),k=1,ndowna(0,il))
942 remd_ene(0,iex)=ene_iex_iex
943 remd_ene(0,i)=ene_i_i
949 cd write (iout,*) "exchange completed"
953 cd write(iout,*) "########",ii
955 i_temp=iran_num(1,nrep)
956 i_mult=iran_num(1,remd_m(i_temp))
957 i_iset=iran_num(1,nset)
958 i_mset=iran_num(1,mset(i_iset))
959 i=i_index(i_temp,i_mult,i_iset,i_mset)
961 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
964 cd write(iout,*) "i_dir=",i_dir
966 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
969 i_mult1=iran_num(1,remd_m(i_temp1))
971 i_mset1=iran_num(1,mset(i_iset1))
972 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
974 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
977 i_mult1=iran_num(1,remd_m(i_temp1))
979 i_mset1=iran_num(1,mset(i_iset1))
980 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
981 econstr_temp_i=remd_ene(20,i)
982 econstr_temp_iex=remd_ene(20,iex)
983 remd_ene(20,i)=remd_ene(n_ene+3,i)
984 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
986 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
989 i_mult1=iran_num(1,remd_m(i_temp1))
991 i_mset1=iran_num(1,mset(i_iset1))
992 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
993 econstr_temp_i=remd_ene(20,i)
994 econstr_temp_iex=remd_ene(20,iex)
995 remd_ene(20,i)=remd_ene(n_ene+3,i)
996 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1002 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1005 c Swap temperatures between conformations i and iex with recalculating the free energies
1006 c following temperature changes.
1007 ene_iex_iex=remd_ene(0,iex)
1008 ene_i_i=remd_ene(0,i)
1009 co write (iout,*) "rescaling weights with temperature",
1011 call rescale_weights(remd_t_bath(i))
1013 call sum_energy(remd_ene(0,iex),.false.)
1014 ene_iex_i=remd_ene(0,iex)
1015 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
1016 c call sum_energy(remd_ene(0,i),.false.)
1017 cd write (iout,*) "ene_i_i",remd_ene(0,i)
1018 c write (iout,*) "rescaling weights with temperature",
1019 c & remd_t_bath(iex)
1020 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1021 c write (iout,*) "ERROR: inconsistent energies:",i,
1022 c & ene_i_i,remd_ene(0,i)
1024 call rescale_weights(remd_t_bath(iex))
1025 call sum_energy(remd_ene(0,i),.false.)
1026 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1027 ene_i_iex=remd_ene(0,i)
1028 c call sum_energy(remd_ene(0,iex),.false.)
1029 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1030 c write (iout,*) "ERROR: inconsistent energies:",iex,
1031 c & ene_iex_iex,remd_ene(0,iex)
1033 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1034 c write (iout,*) "i",i," iex",iex
1035 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1036 cd & " ene_i_iex",ene_i_iex,
1037 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1038 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1039 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1041 cd write(iout,*) 'delta',delta
1042 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1043 c & (remd_ene(i)-remd_ene(iex))/Rb/
1044 c & (remd_t_bath(i)*remd_t_bath(iex))
1045 if (delta .gt. 50.0d0) then
1050 if (i_dir.eq.1.or.i_dir.eq.3)
1051 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1052 if (i_dir.eq.2.or.i_dir.eq.3)
1053 & iremd_tot_usa(int(i2set(i-1)))=
1054 & iremd_tot_usa(int(i2set(i-1)))+1
1055 xxx=ran_number(0.0d0,1.0d0)
1056 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1057 if (delta .gt. xxx) then
1059 remd_t_bath(i)=remd_t_bath(iex)
1060 remd_t_bath(iex)=tmp
1063 iremd_iset(i)=iremd_iset(iex)
1064 iremd_iset(iex)=itmp
1066 remd_ene(0,i)=ene_i_iex
1067 remd_ene(0,iex)=ene_iex_i
1069 if (i_dir.eq.1.or.i_dir.eq.3)
1070 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1073 i2rep(i-1)=i2rep(iex-1)
1076 if (i_dir.eq.2.or.i_dir.eq.3)
1077 & iremd_acc_usa(int(i2set(i-1)))=
1078 & iremd_acc_usa(int(i2set(i-1)))+1
1081 i2set(i-1)=i2set(iex-1)
1084 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1085 i_index(i_temp,i_mult,i_iset,i_mset)=
1086 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1087 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1090 remd_ene(0,iex)=ene_iex_iex
1091 remd_ene(0,i)=ene_i_i
1092 remd_ene(20,iex)=econstr_temp_iex
1093 remd_ene(20,i)=econstr_temp_i
1097 cd do il1=1,mset(il)
1100 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1113 c-------------------------------------
1114 write (iout,*) "NREP",nrep
1116 if(iremd_tot(i).ne.0)
1117 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1118 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1123 if(iremd_tot_usa(i).ne.0)
1124 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1125 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1131 cd write (iout,'(a6,100i4)') "ifirst",
1132 cd & (ifirst(i),i=1,remd_m(1))
1134 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1135 cd & (nupa(i,il),i=1,nupa(0,il))
1136 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1137 cd & (ndowna(i,il),i=1,ndowna(0,il))
1142 cd write (iout,*) "Before scatter"
1144 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1145 & t_bath,1,mpi_double_precision,king,
1147 cd write (iout,*) "After scatter"
1150 & call mpi_scatter(iremd_iset,1,mpi_integer,
1151 & iset,1,mpi_integer,king,
1155 if (me.eq.king .or. .not. out1file) then
1156 write(iout,*) 'REMD scatter time=',time07-time06
1160 call mpi_scatter(elowi,1,mpi_double_precision,
1161 & elow,1,mpi_double_precision,king,
1163 call mpi_scatter(ehighi,1,mpi_double_precision,
1164 & ehigh,1,mpi_double_precision,king,
1167 call rescale_weights(t_bath)
1168 co write (iout,*) "Processor",me,
1169 co & " rescaling weights with temperature",t_bath
1171 stdfp=dsqrt(2*Rb*t_bath/d_time)
1173 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1176 cde write(iout,*) 'REMD after',me,t_bath
1178 if (me.eq.king .or. .not. out1file) then
1179 write(iout,*) 'REMD exchange time=',time08-time00
1185 if (restart1file) then
1186 if (me.eq.king .or. .not. out1file)
1187 & write(iout,*) 'writing restart at the end of run'
1188 call write1rst(i_index)
1191 if (traj1file) call write1traj
1193 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1194 cdeb & icache_all,1,mpi_integer,king,
1195 cdeb & CG_COMM,ierr)
1196 cdeb write(iout,'(a40,8000i8)')
1197 cdeb & ' ntwx_cache after traj1file at the end',
1198 cdeb & (icache_all(i),i=1,nodes)
1203 t_MD=MPI_Wtime()-tt0
1207 if (me.eq.king .or. .not. out1file) then
1208 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1210 & 'MD calculations setup:',t_MDsetup,
1211 & 'Energy & gradient evaluation:',t_enegrad,
1212 & 'Stochastic MD setup:',t_langsetup,
1213 & 'Stochastic MD step setup:',t_sdsetup,
1215 write (iout,'(/28(1h=),a25,27(1h=))')
1216 & ' End of MD calculation '
1221 c-----------------------------------------------------------------------
1222 subroutine write1rst(i_index)
1223 implicit real*8 (a-h,o-z)
1224 include 'DIMENSIONS'
1227 include 'COMMON.IOUNITS'
1228 include 'COMMON.REMD'
1229 include 'COMMON.SETUP'
1230 include 'COMMON.CHAIN'
1231 include 'COMMON.SBRIDGE'
1232 include 'COMMON.INTERACT'
1234 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1235 & d_restart2(3,2*maxres*maxprocs)
1239 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1240 common /przechowalnia/ d_restart1,d_restart2
1245 t5_restart1(4)=t_bath
1246 t5_restart1(5)=Uconst
1248 call mpi_gather(t5_restart1,5,mpi_real,
1249 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1257 call mpi_gather(r_d,3*2*nres,mpi_real,
1258 & d_restart1,3*2*nres,mpi_real,king,
1267 call mpi_gather(r_d,3*2*nres,mpi_real,
1268 & d_restart2,3*2*nres,mpi_real,king,
1273 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1275 call xdrfint_(ixdrf, i2rep(i), iret)
1278 call xdrfint_(ixdrf, ifirst(i), iret)
1282 call xdrfint_(ixdrf, nupa(i,il), iret)
1286 call xdrfint_(ixdrf, ndowna(i,il), iret)
1292 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1299 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1306 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1312 call xdrfint_(ixdrf, nset, iret)
1314 call xdrfint_(ixdrf,mset(i), iret)
1317 call xdrfint_(ixdrf,i2set(i), iret)
1323 itmp=i_index(i,j,il,il1)
1324 call xdrfint_(ixdrf,itmp, iret)
1331 call xdrfclose_(ixdrf, iret)
1333 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1335 call xdrfint(ixdrf, i2rep(i), iret)
1338 call xdrfint(ixdrf, ifirst(i), iret)
1342 call xdrfint(ixdrf, nupa(i,il), iret)
1346 call xdrfint(ixdrf, ndowna(i,il), iret)
1352 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1359 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1366 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1373 call xdrfint(ixdrf, nset, iret)
1375 call xdrfint(ixdrf,mset(i), iret)
1378 call xdrfint(ixdrf,i2set(i), iret)
1384 itmp=i_index(i,j,il,il1)
1385 call xdrfint(ixdrf,itmp, iret)
1392 call xdrfclose(ixdrf, iret)
1399 subroutine write1traj
1400 implicit real*8 (a-h,o-z)
1401 include 'DIMENSIONS'
1404 include 'COMMON.IOUNITS'
1405 include 'COMMON.REMD'
1406 include 'COMMON.SETUP'
1407 include 'COMMON.CHAIN'
1408 include 'COMMON.SBRIDGE'
1409 include 'COMMON.INTERACT'
1413 real xcoord(3,maxres2+2),prec
1414 real r_qfrag(50),r_qpair(100)
1415 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1416 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1417 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1418 & p_uscdiff(100*maxprocs)
1419 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1420 common /przechowalnia/ p_c
1422 call mpi_bcast(ii_write,1,mpi_integer,
1423 & king,CG_COMM,ierr)
1426 print *,'traj1file',me,ii_write,ntwx_cache
1430 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1432 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1435 t5_restart1(1)=totT_cache(ii)
1436 t5_restart1(2)=EK_cache(ii)
1437 t5_restart1(3)=potE_cache(ii)
1438 t5_restart1(4)=t_bath_cache(ii)
1439 t5_restart1(5)=Uconst_cache(ii)
1440 call mpi_gather(t5_restart1,5,mpi_real,
1441 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1443 call mpi_gather(iset_cache(ii),1,mpi_integer,
1444 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1447 r_qfrag(i)=qfrag_cache(i,ii)
1450 r_qpair(i)=qpair_cache(i,ii)
1453 r_utheta(i)=utheta_cache(i,ii)
1454 r_ugamma(i)=ugamma_cache(i,ii)
1455 r_uscdiff(i)=uscdiff_cache(i,ii)
1458 call mpi_gather(r_qfrag,nfrag,mpi_real,
1459 & p_qfrag,nfrag,mpi_real,king,
1461 call mpi_gather(r_qpair,npair,mpi_real,
1462 & p_qpair,npair,mpi_real,king,
1464 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1465 & p_utheta,nfrag_back,mpi_real,king,
1467 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1468 & p_ugamma,nfrag_back,mpi_real,king,
1470 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1471 & p_uscdiff,nfrag_back,mpi_real,king,
1475 write (iout,*) "p_qfrag"
1477 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1479 write (iout,*) "p_qpair"
1481 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1487 r_c(j,i)=c_cache(j,i,ii)
1491 call mpi_gather(r_c,3*2*nres,mpi_real,
1492 & p_c,3*2*nres,mpi_real,king,
1498 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1499 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1500 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1501 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1502 call xdrfint_(ixdrf, nss, iret)
1504 call xdrfint_(ixdrf, ihpb(j), iret)
1505 call xdrfint_(ixdrf, jhpb(j), iret)
1507 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1508 call xdrfint_(ixdrf, iset_restart1(il), iret)
1510 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1513 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1516 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1517 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1518 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1523 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1528 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1532 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1536 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1537 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1538 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1539 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1540 call xdrfint(ixdrf, nss, iret)
1542 call xdrfint(ixdrf, ihpb(j), iret)
1543 call xdrfint(ixdrf, jhpb(j), iret)
1545 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1546 call xdrfint(ixdrf, iset_restart1(il), iret)
1548 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1551 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1554 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1555 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1556 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1561 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1566 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1570 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1576 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1578 if(me.eq.king) call xdrfclose(ixdrf, iret)
1580 do i=1,ntwx_cache-ii_write
1582 totT_cache(i)=totT_cache(ii_write+i)
1583 EK_cache(i)=EK_cache(ii_write+i)
1584 potE_cache(i)=potE_cache(ii_write+i)
1585 t_bath_cache(i)=t_bath_cache(ii_write+i)
1586 Uconst_cache(i)=Uconst_cache(ii_write+i)
1587 iset_cache(i)=iset_cache(ii_write+i)
1590 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1593 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1596 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1597 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1598 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1603 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1607 ntwx_cache=ntwx_cache-ii_write
1612 subroutine read1restart(i_index)
1613 implicit real*8 (a-h,o-z)
1614 include 'DIMENSIONS'
1617 include 'COMMON.IOUNITS'
1618 include 'COMMON.REMD'
1619 include 'COMMON.SETUP'
1620 include 'COMMON.CHAIN'
1621 include 'COMMON.SBRIDGE'
1622 include 'COMMON.INTERACT'
1623 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1626 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1627 common /przechowalnia/ d_restart1
1628 write (*,*) "Processor",me," called read1restart"
1631 open(irest2,file=mremd_rst_name,status='unknown')
1632 read(irest2,*,err=334) i
1633 write(iout,*) "Reading old rst in ASCI format"
1635 call read1restart_old
1639 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1642 call xdrfint_(ixdrf, i2rep(i), iret)
1645 call xdrfint_(ixdrf, ifirst(i), iret)
1648 call xdrfint_(ixdrf, nupa(0,il), iret)
1650 call xdrfint_(ixdrf, nupa(i,il), iret)
1653 call xdrfint_(ixdrf, ndowna(0,il), iret)
1655 call xdrfint_(ixdrf, ndowna(i,il), iret)
1660 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
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)
1690 call mpi_scatter(t_restart1,5,mpi_real,
1691 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1695 t_bath=t5_restart1(4)
1700 c read(irest2,'(3e15.5)')
1701 c & (d_restart1(j,i+2*nres*il),j=1,3)
1704 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1706 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1712 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1713 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1723 c read(irest2,'(3e15.5)')
1724 c & (d_restart1(j,i+2*nres*il),j=1,3)
1727 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1729 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1735 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1736 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1747 call xdrfint_(ixdrf, nset, iret)
1749 call xdrfint_(ixdrf,mset(i), iret)
1752 call xdrfint_(ixdrf,i2set(i), iret)
1758 call xdrfint_(ixdrf,itmp, iret)
1759 i_index(i,j,il,il1)=itmp
1767 call xdrfint(ixdrf, nset, iret)
1769 call xdrfint(ixdrf,mset(i), iret)
1772 call xdrfint(ixdrf,i2set(i), iret)
1778 call xdrfint(ixdrf,itmp, iret)
1779 i_index(i,j,il,il1)=itmp
1786 call mpi_scatter(i2set,1,mpi_integer,
1787 & iset,1,mpi_integer,king,
1793 if(me.eq.king) close(irest2)
1797 subroutine read1restart_old
1798 implicit real*8 (a-h,o-z)
1799 include 'DIMENSIONS'
1802 include 'COMMON.IOUNITS'
1803 include 'COMMON.REMD'
1804 include 'COMMON.SETUP'
1805 include 'COMMON.CHAIN'
1806 include 'COMMON.SBRIDGE'
1807 include 'COMMON.INTERACT'
1808 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1810 common /przechowalnia/ d_restart1
1812 open(irest2,file=mremd_rst_name,status='unknown')
1813 read (irest2,*) (i2rep(i),i=0,nodes-1)
1814 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1816 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1817 read (irest2,*) ndowna(0,il),
1818 & (ndowna(i,il),i=1,ndowna(0,il))
1821 read(irest2,*) (t_restart1(j,il),j=1,4)
1824 call mpi_scatter(t_restart1,5,mpi_real,
1825 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1829 t_bath=t5_restart1(4)
1834 read(irest2,'(3e15.5)')
1835 & (d_restart1(j,i+2*nres*il),j=1,3)
1839 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1840 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1850 read(irest2,'(3e15.5)')
1851 & (d_restart1(j,i+2*nres*il),j=1,3)
1855 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1856 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1862 if(me.eq.king) close(irest2)
1865 c------------------------------------------
1866 subroutine returnbox
1867 include 'DIMENSIONS'
1869 include 'COMMON.CONTROL'
1870 include 'COMMON.VAR'
1873 include 'COMMON.LANGEVIN'
1875 include 'COMMON.LANGEVIN.lang0'
1877 include 'COMMON.CHAIN'
1878 include 'COMMON.DERIV'
1879 include 'COMMON.GEO'
1880 include 'COMMON.LOCAL'
1881 include 'COMMON.INTERACT'
1882 include 'COMMON.IOUNITS'
1883 include 'COMMON.NAMES'
1884 include 'COMMON.TIME1'
1885 include 'COMMON.REMD'
1886 include 'COMMON.SETUP'
1887 include 'COMMON.MUCA'
1888 include 'COMMON.HAIRPIN'
1892 C write(*,*) 'initial', i,j,c(j,i)
1895 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
1897 if (allareout.eq.1) then
1898 ireturnval=int(c(j,i)/boxxsize)
1899 if (c(j,i).le.0) ireturnval=ireturnval-1
1900 do k=chain_beg,chain_end
1901 c(j,k)=c(j,k)-ireturnval*boxxsize
1902 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
1908 if (int(c(j,i)/boxxsize).eq.0) allareout=0
1911 if (allareout.eq.1) then
1912 ireturnval=int(c(j,i)/boxxsize)
1913 if (c(j,i).le.0) ireturnval=ireturnval-1
1915 c(j,k)=c(j,k)-ireturnval*boxxsize
1916 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
1921 C write(*,*) 'befor no jump', i,j,c(j,i)
1925 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1926 difference=abs(c(j,i-1)-c(j,i))
1927 C print *,'diff', difference
1928 if (difference.gt.boxxsize/2.0) then
1929 if (c(j,i-1).gt.c(j,i)) then
1938 c(j,k)=c(j,k)+nojumpval*boxxsize
1939 c(j,k+nres)=c(j,k+nres)+nojumpval*boxxsize
1943 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1944 difference=abs(c(j,i-1)-c(j,i))
1945 if (difference.gt.boxxsize/2.0) then
1946 if (c(j,i-1).gt.c(j,i)) then
1955 c(j,k)=c(j,k)+nojumpval*boxxsize
1956 c(j,k+nres)=c(j,k+nres)+nojumpval*boxxsize
1960 write(*,*) 'after no jump', i,j,c(j,i)
1967 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
1969 if (allareout.eq.1) then
1970 ireturnval=int(c(j,i)/boxysize)
1971 if (c(j,i).le.0) ireturnval=ireturnval-1
1972 do k=chain_beg,chain_end
1973 c(j,k)=c(j,k)-ireturnval*boxysize
1974 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
1980 if (int(c(j,i)/boxysize).eq.0) allareout=0
1983 if (allareout.eq.1) then
1984 ireturnval=int(c(j,i)/boxysize)
1985 if (c(j,i).le.0) ireturnval=ireturnval-1
1987 c(j,k)=c(j,k)-ireturnval*boxysize
1988 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
1993 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
1994 difference=abs(c(j,i-1)-c(j,i))
1995 if (difference.gt.boxysize/2.0) then
1996 if (c(j,i-1).gt.c(j,i)) then
2005 c(j,k)=c(j,k)+nojumpval*boxysize
2006 c(j,k+nres)=c(j,k+nres)+nojumpval*boxysize
2010 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2011 difference=abs(c(j,i-1)-c(j,i))
2012 if (difference.gt.boxysize/2.0) then
2013 if (c(j,i-1).gt.c(j,i)) then
2022 c(j,k)=c(j,k)+nojumpval*boxysize
2023 c(j,k+nres)=c(j,k+nres)+nojumpval*boxysize
2029 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
2031 if (allareout.eq.1) then
2032 ireturnval=int(c(j,i)/boxysize)
2033 if (c(j,i).le.0) ireturnval=ireturnval-1
2034 do k=chain_beg,chain_end
2035 c(j,k)=c(j,k)-ireturnval*boxzsize
2036 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
2042 if (int(c(j,i)/boxzsize).eq.0) allareout=0
2045 if (allareout.eq.1) then
2046 ireturnval=int(c(j,i)/boxzsize)
2047 if (c(j,i).le.0) ireturnval=ireturnval-1
2049 c(j,k)=c(j,k)-ireturnval*boxzsize
2050 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
2055 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2056 difference=abs(c(j,i-1)-c(j,i))
2057 if (difference.gt.(boxzsize/2.0)) then
2058 if (c(j,i-1).gt.c(j,i)) then
2067 c(j,k)=c(j,k)+nojumpval*boxzsize
2068 c(j,k+nres)=c(j,k+nres)+nojumpval*boxzsize
2072 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
2073 difference=abs(c(j,i-1)-c(j,i))
2074 if (difference.gt.boxzsize/2.0) then
2075 if (c(j,i-1).gt.c(j,i)) then
2084 c(j,k)=c(j,k)+nojumpval*boxzsize
2085 c(j,k+nres)=c(j,k+nres)+nojumpval*boxzsize