2 !-----------------------------------------------------------------------------
10 ! use control_data, only:maxprocs
14 !-----------------------------------------------------------------------------
16 ! common /remdrestart/
17 integer(kind=2),dimension(:),allocatable :: i2set !(0:maxprocs)
18 integer(kind=2),dimension(:),allocatable :: ifirst !(maxprocs)
19 integer(kind=2),dimension(:,:),allocatable :: nupa,&
20 ndowna !(0:maxprocs/4,0:maxprocs)
21 real(kind=4),dimension(:,:),allocatable :: t_restart1 !(5,maxprocs)
22 integer,dimension(:),allocatable :: iset_restart1 !(maxprocs)
24 real(kind=4),dimension(:),allocatable :: totT_cache,EK_cache,&
25 potE_cache,t_bath_cache,Uconst_cache !(max_cache_traj)
26 real(kind=4),dimension(:,:),allocatable :: qfrag_cache !(50,max_cache_traj)
27 real(kind=4),dimension(:,:),allocatable :: qpair_cache !(100,max_cache_traj)
28 real(kind=4),dimension(:,:),allocatable :: ugamma_cache,&
29 utheta_cache,uscdiff_cache !(maxfrag_back,max_cache_traj)
30 real(kind=4),dimension(:,:,:),allocatable :: c_cache !(3,maxres2+2,max_cache_traj)
31 integer :: ntwx_cache,ii_write !,max_cache_traj_use
32 integer,dimension(:),allocatable :: iset_cache !(max_cache_traj)
33 !-----------------------------------------------------------------------------
34 ! common /przechowalnia/
35 real(kind=4),dimension(:,:),allocatable :: d_restart1 !(3,2*nres*maxprocs)
36 real(kind=4),dimension(:,:),allocatable :: d_restart2 !(3,2*nres*maxprocs)
37 real(kind=4),dimension(:,:),allocatable :: p_c !(3,(nres2+2)*maxprocs)
38 !-----------------------------------------------------------------------------
41 !-----------------------------------------------------------------------------
43 !-----------------------------------------------------------------------------
45 !-----------------------------------------------------------------------------
47 !-----------------------------------------------------------------------------
52 use control, only:tcpu,ovrtim
53 use io_base, only:ilen
56 use random, only: iran_num,ran_number
57 use compare, only:hairpin,secondary2
58 use io, only:cartout,statout
59 ! implicit real*8 (a-h,o-z)
60 ! include 'DIMENSIONS'
62 ! include 'COMMON.CONTROL'
63 ! include 'COMMON.VAR'
66 ! include 'COMMON.LANGEVIN'
68 ! include 'COMMON.LANGEVIN.lang0'
70 ! include 'COMMON.CHAIN'
71 ! include 'COMMON.DERIV'
72 ! include 'COMMON.GEO'
73 ! include 'COMMON.LOCAL'
74 ! include 'COMMON.INTERACT'
75 ! include 'COMMON.IOUNITS'
76 ! include 'COMMON.NAMES'
77 ! include 'COMMON.TIME1'
78 ! include 'COMMON.REMD'
79 ! include 'COMMON.SETUP'
80 ! include 'COMMON.MUCA'
81 ! include 'COMMON.HAIRPIN'
83 real(kind=8),dimension(3) :: L,vcm
84 real(kind=8) :: energia(0:n_ene)
85 real(kind=8) :: remd_t_bath(Nprocs) !(maxprocs)
86 integer :: iremd_iset(Nprocs) !(maxprocs)
87 integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
88 ! (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
89 real(kind=8) :: remd_ene(0:n_ene+4,Nprocs) !(0:n_ene+4,maxprocs)
90 integer :: iremd_acc(Nprocs),iremd_tot(Nprocs) !(maxprocs)
91 integer :: iremd_acc_usa(Nprocs),iremd_tot_usa(Nprocs) !(maxprocs)
92 integer :: rstcount !el ilen,
94 character(len=50) :: tytul
97 !old integer nup(0:maxprocs),ndown(0:maxprocs)
98 integer :: rep2i(0:Nprocs),ireqi(Nprocs) !(maxprocs)
99 integer :: icache_all(Nprocs) !(maxprocs)
100 integer :: status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,Nprocs)! (MPI_STATUS_SIZE,maxprocs)
101 logical :: synflag, end_of_run, file_exist = .false.!, ovrtim
103 real(kind=8) :: delta,time00,time01,time001,time02,time03,time04,&
104 time05,time06,time07,time08,tt0,scalfac,ene_iex_iex,&
105 ene_i_i,ene_iex_i,ene_i_iex,xxx,tmp,econstr_temp_i,&
107 integer :: k,il,il1,i,j,nharp,ii,ierr,itime_master,irr,iex,&
108 i_set_temp,itmp,i_temp,i_mult,i_iset,i_mset,i_dir,i_temp1,&
109 i_mult1,i_iset1,i_mset1,ierror
110 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
111 !deb imin_itime_old=0
119 if(me.eq.king.or..not.out1file) then
120 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
121 write (iout,*) "NREP=",nrep
125 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
126 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
128 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
130 !d print *,'MREMD',nodes
131 !d print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
132 !de write (iout,*) "Start MREMD: me",me," t_bath",t_bath
136 do il1=1,max0(mset(il),1)
146 i_index(i,j,il,il1)=k
152 if(me.eq.king.or..not.out1file) then
153 write(iout,*) (i2rep(i),i=0,nodes-1)
154 write(iout,*) (i2set(i),i=0,nodes-1)
159 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
166 ! print *,'i2rep',me,i2rep(me)
167 ! print *,'rep2i',(rep2i(i),i=0,nrep)
169 !old if (i2rep(me).eq.nrep) then
172 !old nup(0)=remd_m(i2rep(me)+1)
173 !old k=rep2i(int(i2rep(me)))+1
180 !d print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
182 !old if (i2rep(me).eq.1) then
185 !old ndown(0)=remd_m(i2rep(me)-1)
186 !old k=rep2i(i2rep(me)-2)+1
193 !d print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
195 !el common /przechowalnia/
196 if(.not.allocated(d_restart1)) allocate(d_restart1(3,(nres2+1)*nodes))
197 if(.not.allocated(d_restart2)) allocate(d_restart2(3,(nres2+1)*nodes))
198 if(.not.allocated(p_c)) allocate(p_c(3,(nres2+2)*nodes))
201 write (*,*) "Processor",me," rest",rest,&
202 "restart1fie",restart1file
203 if(rest.and.restart1file) then
205 inquire(file=mremd_rst_name,exist=file_exist)
206 !d write (*,*) me," Before broadcast: file_exist",file_exist
207 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
209 !d write (*,*) me," After broadcast: file_exist",file_exist
211 if(me.eq.king.or..not.out1file) &
212 write (iout,*) 'Master is reading restart1file'
213 call read1restart(i_index)
215 if(me.eq.king.or..not.out1file) &
216 write (iout,*) 'WARNING : no restart1file'
219 if(me.eq.king.or..not.out1file) then
220 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
221 write(iout,*) "i_index"
226 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
235 if (rest.and..not.restart1file) &
236 inquire(file=mremd_rst_name,exist=file_exist)
237 if(.not.file_exist.and.rest.and..not.restart1file) &
238 write(iout,*) 'WARNING : no restart file',mremd_rst_name
239 IF (rest.and.file_exist.and..not.restart1file) THEN
240 write (iout,*) 'Master is reading restart file',&
242 open(irest2,file=mremd_rst_name,status='unknown')
244 read (irest2,*) (i2rep(i),i=0,nodes-1)
246 read (irest2,*) (ifirst(i),i=1,remd_m(1))
249 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
251 read (irest2,*) ndowna(0,il),&
252 (ndowna(i,il),i=1,ndowna(0,il))
258 read (irest2,*) (mset(i),i=1,nset)
260 read (irest2,*) (i2set(i),i=0,nodes-1)
265 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
270 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
271 write(iout,*) "i_index"
276 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
285 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
286 write (iout,'(a6,1000i5)') "ifirst",&
287 (ifirst(i),i=1,remd_m(1))
289 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
290 (nupa(i,il),i=1,nupa(0,il))
291 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
292 (ndowna(i,il),i=1,ndowna(0,il))
294 ELSE IF (.not.(rest.and.file_exist)) THEN
300 if (i2rep(il-1).eq.nrep) then
303 nupa(0,il)=remd_m(i2rep(il-1)+1)
304 k=rep2i(int(i2rep(il-1)))+1
310 if (i2rep(il-1).eq.1) then
313 ndowna(0,il)=remd_m(i2rep(il-1)-1)
314 k=rep2i(i2rep(il-1)-2)+1
322 write (iout,'(a6,100i4)') "ifirst",&
323 (ifirst(i),i=1,remd_m(1))
325 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
326 (nupa(i,il),i=1,nupa(0,il))
327 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
328 (ndowna(i,il),i=1,ndowna(0,il))
334 ! t_bath=retmin+(retmax-retmin)*me/(nodes-1)
335 if(.not.(rest.and.file_exist.and.restart1file)) then
336 if (me .eq. king) then
339 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
341 !d print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
342 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
347 if(me.eq.king.or..not.out1file) &
348 write(iout,*) me,"iset=",iset,"t_bath=",t_bath
352 stdfp(i)=dsqrt(2*Rb*t_bath/d_time)
356 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
360 ! print *,'irep',me,t_bath
362 if (me.eq.king .or. .not. out1file) &
363 write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
364 call rescale_weights(t_bath)
368 !------copy MD--------------
369 ! The driver for molecular dynamics subroutines
370 !------------------------------------------------
376 if(me.eq.king.or..not.out1file) &
377 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
383 ! Determine the inverse of the inertia matrix.
384 call setup_MD_matrices
388 if (me.eq.king .or. .not. out1file) &
389 write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
391 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
393 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
396 call rescale_weights(t_bath)
400 t_MDsetup = MPI_Wtime()-tt0
402 t_MDsetup = tcpu()-tt0
405 ! Entering the MD loop
411 if (lang.eq.2 .or. lang.eq.3) then
415 call sd_verlet_p_setup
417 call sd_verlet_ciccotti_setup
421 pfric0_mat(i,j,0)=pfric_mat(i,j)
422 afric0_mat(i,j,0)=afric_mat(i,j)
423 vfric0_mat(i,j,0)=vfric_mat(i,j)
424 prand0_mat(i,j,0)=prand_mat(i,j)
425 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
426 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
431 flag_stoch(i)=.false.
435 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
437 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
441 else if (lang.eq.1 .or. lang.eq.4) then
445 if (me.eq.king .or. .not. out1file) &
446 write(iout,*) 'Setup time',time00-walltime
449 t_langsetup=MPI_Wtime()-tt0
452 t_langsetup=tcpu()-tt0
458 do while(.not.end_of_run)
460 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
461 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
463 if (lang.gt.0 .and. surfarea .and. &
464 mod(itime,reset_fricmat).eq.0) then
465 if (lang.eq.2 .or. lang.eq.3) then
469 call sd_verlet_p_setup
471 call sd_verlet_ciccotti_setup
475 pfric0_mat(i,j,0)=pfric_mat(i,j)
476 afric0_mat(i,j,0)=afric_mat(i,j)
477 vfric0_mat(i,j,0)=vfric_mat(i,j)
478 prand0_mat(i,j,0)=prand_mat(i,j)
479 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
480 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
485 flag_stoch(i)=.false.
488 else if (lang.eq.1 .or. lang.eq.4) then
491 write (iout,'(a,i10)') &
492 "Friction matrix reset based on surface area, itime",itime
494 if (reset_vel .and. tbf .and. lang.eq.0 &
495 .and. mod(itime,count_reset_vel).eq.0) then
497 if (me.eq.king .or. .not. out1file) &
498 write(iout,'(a,f20.2)') &
499 "Velocities reset to random values, time",totT
502 d_t_old(j,i)=d_t(j,i)
506 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
510 d_t(j,0)=d_t(j,0)-vcm(j)
513 kinetic_T=2.0d0/(dimen3*Rb)*EK
514 scalfac=dsqrt(T_bath/kinetic_T)
515 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
518 d_t_old(j,i)=scalfac*d_t(j,i)
524 ! Time-reversible RESPA algorithm
525 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
526 call RESPA_step(itime)
528 ! Variable time step algorithm.
529 call velverlet_step(itime)
533 call brown_step(itime)
535 print *,"Brown dynamics not here!"
537 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
543 if (mod(itime,ntwe).eq.0) call statout(itime)
545 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
546 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
548 call hairpin(.true.,nharp,iharp)
549 call secondary2(.true.)
550 call pdbout(potE,tytul,ipdb)
555 if (mod(itime,ntwx).eq.0.and.traj1file) then
556 if(ntwx_cache.lt.max_cache_traj_use) then
557 ntwx_cache=ntwx_cache+1
559 if (max_cache_traj_use.ne.1) &
560 print *,itime,"processor ",me," over cache ",ntwx_cache
563 totT_cache(i)=totT_cache(i+1)
564 EK_cache(i)=EK_cache(i+1)
565 potE_cache(i)=potE_cache(i+1)
566 t_bath_cache(i)=t_bath_cache(i+1)
567 Uconst_cache(i)=Uconst_cache(i+1)
568 iset_cache(i)=iset_cache(i+1)
571 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
574 qpair_cache(ii,i)=qpair_cache(ii,i+1)
577 utheta_cache(ii,i)=utheta_cache(ii,i+1)
578 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
579 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
585 c_cache(j,ii,i)=c_cache(j,ii,i+1)
591 totT_cache(ntwx_cache)=totT
592 EK_cache(ntwx_cache)=EK
593 potE_cache(ntwx_cache)=potE
594 t_bath_cache(ntwx_cache)=t_bath
595 Uconst_cache(ntwx_cache)=Uconst
596 iset_cache(ntwx_cache)=iset
599 qfrag_cache(i,ntwx_cache)=qfrag(i)
602 qpair_cache(i,ntwx_cache)=qpair(i)
605 utheta_cache(i,ntwx_cache)=utheta(i)
606 ugamma_cache(i,ntwx_cache)=ugamma(i)
607 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
612 c_cache(j,i,ntwx_cache)=c(j,i)
617 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
618 .and..not.restart1file) then
621 open(irest1,file=mremd_rst_name,status='unknown')
622 write (irest1,*) "i2rep"
623 write (irest1,*) (i2rep(i),i=0,nodes-1)
624 write (irest1,*) "ifirst"
625 write (irest1,*) (ifirst(i),i=1,remd_m(1))
627 write (irest1,*) "nupa",il
628 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
629 write (irest1,*) "ndowna",il
630 write (irest1,*) ndowna(0,il),&
631 (ndowna(i,il),i=1,ndowna(0,il))
634 write (irest1,*) "nset"
635 write (irest1,*) nset
636 write (irest1,*) "mset"
637 write (irest1,*) (mset(i),i=1,nset)
638 write (irest1,*) "i2set"
639 write (irest1,*) (i2set(i),i=0,nodes-1)
640 write (irest1,*) "i_index"
644 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
652 open(irest2,file=rest2name,status='unknown')
653 write(irest2,*) totT,EK,potE,totE,t_bath
655 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
658 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
661 write (irest2,*) iset
668 ! forced synchronization
669 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
670 .and. .not. mremdsync) then
672 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
674 call mpi_recv(itime_master, 1, MPI_INTEGER,&
675 0,101,CG_COMM, status, ierr)
676 call mpi_barrier(CG_COMM, ierr)
677 !deb if (out1file.or.traj1file) then
678 !deb call mpi_gather(itime,1,mpi_integer,
679 !deb & icache_all,1,mpi_integer,king,
682 call mpi_gather(ntwx_cache,1,mpi_integer,&
683 icache_all,1,mpi_integer,king,&
686 write(iout,*) 'REMD synchro at',itime_master,itime
687 if (itime_master.ge.n_timestep .or. ovrtim()) &
689 !time call flush(iout)
694 if ((mod(itime,nstex).eq.0.and.me.eq.king &
695 .or.end_of_run.and.me.eq.king ) &
696 .and. .not. mremdsync ) then
699 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
700 CG_COMM, ireqi(i), ierr)
701 !d write(iout,*) 'REMD synchro with',i
704 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
705 call mpi_barrier(CG_COMM, ierr)
707 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
708 if (out1file.or.traj1file) then
709 !deb call mpi_gather(itime,1,mpi_integer,
710 !deb & itime_all,1,mpi_integer,king,
712 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
713 !deb & (itime_all(i),i=1,nodes)
715 !deb imin_itime=itime_all(1)
717 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
719 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
720 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
721 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
722 call mpi_gather(ntwx_cache,1,mpi_integer,&
723 icache_all,1,mpi_integer,king,&
725 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
726 ! & (icache_all(i),i=1,nodes)
727 ii_write=icache_all(1)
729 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
731 ! write(iout,*) "MIN ii_write=",ii_write
734 !time call flush(iout)
736 if(mremdsync .and. mod(itime,nstex).eq.0) then
738 if (me.eq.king .or. .not. out1file) &
739 write(iout,*) 'REMD synchro at',itime
742 call mpi_gather(ntwx_cache,1,mpi_integer,&
743 icache_all,1,mpi_integer,king,&
746 write(iout,'(a19,8000i8)') ' ntwx_cache',&
747 (icache_all(i),i=1,nodes)
748 ii_write=icache_all(1)
750 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
752 write(iout,*) "MIN ii_write=",ii_write
758 ! Update the time safety limiy
759 if (time001-time00.gt.safety) then
760 safety=time001-time00+600
761 write (iout,*) "****** SAFETY increased to",safety," s"
763 if (ovrtim()) end_of_run=.true.
765 if(synflag.and..not.end_of_run) then
769 write(iout,*) 'REMD before',me,t_bath
771 ! call mpi_gather(t_bath,1,mpi_double_precision,
772 ! & remd_t_bath,1,mpi_double_precision,king,
774 potEcomp(n_ene+1)=t_bath
776 potEcomp(n_ene+2)=iset
777 if (iset.lt.nset) then
781 potEcomp(n_ene+3)=Uconst
788 potEcomp(n_ene+4)=Uconst
792 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
793 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
796 call mpi_gather(elow,1,mpi_double_precision,&
797 elowi,1,mpi_double_precision,king,&
799 call mpi_gather(ehigh,1,mpi_double_precision,&
800 ehighi,1,mpi_double_precision,king,&
805 if (me.eq.king .or. .not. out1file) then
806 write(iout,*) 'REMD gather times=',time03-time01 &
810 if (restart1file) call write1rst(i_index)
813 if (me.eq.king .or. .not. out1file) then
814 write(iout,*) 'REMD writing rst time=',time04-time03
817 if (traj1file) call write1traj
819 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
820 !deb & icache_all,1,mpi_integer,king,
822 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
823 !deb & (icache_all(i),i=1,nodes)
828 if (me.eq.king .or. .not. out1file) then
829 write(iout,*) 'REMD writing traj time=',time05-time04
836 remd_t_bath(i)=remd_ene(n_ene+1,i)
837 iremd_iset(i)=remd_ene(n_ene+2,i)
841 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
843 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
847 write(iout,*) 'REMD exchange temp,ene'
849 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
850 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
854 !-------------------------------------
856 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
859 write (iout,*) "remd_m(1)",remd_m(1)
861 i=ifirst(iran_num(1,remd_m(1)))
868 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
870 if(i.gt.0.and.nupa(0,i).gt.0) then
872 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
874 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
876 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
878 ! do while (iex.eq.i)
879 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
880 iex=nupa(iran_num(1,int(nupa(0,i))),i)
882 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
884 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
886 ! Swap temperatures between conformations i and iex with recalculating the free energies
887 ! following temperature changes.
888 ene_iex_iex=remd_ene(0,iex)
889 ene_i_i=remd_ene(0,i)
890 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
891 ! & " iex",iex," ene_iex_iex",ene_iex_iex
892 ! write (iout,*) "rescaling weights with temperature",
895 call rescale_weights(remd_t_bath(i))
897 ! write (iout,*) "0,iex",remd_t_bath(i)
898 ! call enerprint(remd_ene(0,iex))
900 call sum_energy(remd_ene(0,iex),.false.)
901 ene_iex_i=remd_ene(0,iex)
902 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
904 ! write (iout,*) "0,i",remd_t_bath(i)
905 ! call enerprint(remd_ene(0,i))
907 call sum_energy(remd_ene(0,i),.false.)
908 ! write (iout,*) "ene_i_i",remd_ene(0,i)
910 ! write (iout,*) "rescaling weights with temperature",
912 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
913 write (iout,*) "ERROR: inconsistent energies:",i,&
914 ene_i_i,remd_ene(0,i)
916 call rescale_weights(remd_t_bath(iex))
918 ! write (iout,*) "0,i",remd_t_bath(iex)
919 ! call enerprint(remd_ene(0,i))
921 call sum_energy(remd_ene(0,i),.false.)
922 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
924 ene_i_iex=remd_ene(0,i)
926 ! write (iout,*) "0,iex",remd_t_bath(iex)
927 ! call enerprint(remd_ene(0,iex))
929 call sum_energy(remd_ene(0,iex),.false.)
930 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
931 write (iout,*) "ERROR: inconsistent energies:",iex,&
932 ene_iex_iex,remd_ene(0,iex)
934 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
935 ! write (iout,*) "i",i," iex",iex
936 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
937 ! & " ene_i_iex",ene_i_iex,
938 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
940 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
941 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
943 ! write(iout,*) 'delta',delta
944 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
945 ! & (remd_ene(i)-remd_ene(iex))/Rb/
946 ! & (remd_t_bath(i)*remd_t_bath(iex))
948 if (delta .gt. 50.0d0) then
954 else if (delta.lt.-50.0d0) then
963 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
964 xxx=ran_number(0.0d0,1.0d0)
965 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
967 if (delta .gt. xxx) then
969 remd_t_bath(i)=remd_t_bath(iex)
971 remd_ene(0,i)=ene_i_iex
972 remd_ene(0,iex)=ene_iex_i
978 ehighi(i)=ehighi(iex)
985 nupa(k,i)=nupa(k,iex)
988 ndowna(k,i)=ndowna(k,iex)
992 if (ifirst(il).eq.i) ifirst(il)=iex
994 if (nupa(k,il).eq.i) then
996 elseif (nupa(k,il).eq.iex) then
1001 if (ndowna(k,il).eq.i) then
1003 elseif (ndowna(k,il).eq.iex) then
1009 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1011 i2rep(i-1)=i2rep(iex-1)
1014 ! write(iout,*) 'exchange',i,iex
1015 ! write (iout,'(a8,100i4)') "@ ifirst",
1016 ! & (ifirst(k),k=1,remd_m(1))
1018 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1019 ! & (nupa(k,il),k=1,nupa(0,il))
1020 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1021 ! & (ndowna(k,il),k=1,ndowna(0,il))
1026 remd_ene(0,iex)=ene_iex_iex
1027 remd_ene(0,i)=ene_i_i
1033 !d write (iout,*) "exchange completed"
1037 !d write(iout,*) "########",ii
1039 i_temp=iran_num(1,nrep)
1040 i_mult=iran_num(1,remd_m(i_temp))
1041 i_iset=iran_num(1,nset)
1042 i_mset=iran_num(1,mset(i_iset))
1043 i=i_index(i_temp,i_mult,i_iset,i_mset)
1045 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1048 !d write(iout,*) "i_dir=",i_dir
1050 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1053 i_mult1=iran_num(1,remd_m(i_temp1))
1055 i_mset1=iran_num(1,mset(i_iset1))
1056 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1058 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1061 i_mult1=iran_num(1,remd_m(i_temp1))
1063 i_mset1=iran_num(1,mset(i_iset1))
1064 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1065 econstr_temp_i=remd_ene(20,i)
1066 econstr_temp_iex=remd_ene(20,iex)
1067 remd_ene(20,i)=remd_ene(n_ene+3,i)
1068 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1070 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1073 i_mult1=iran_num(1,remd_m(i_temp1))
1075 i_mset1=iran_num(1,mset(i_iset1))
1076 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1077 econstr_temp_i=remd_ene(20,i)
1078 econstr_temp_iex=remd_ene(20,iex)
1079 remd_ene(20,i)=remd_ene(n_ene+3,i)
1080 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1086 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1089 ! Swap temperatures between conformations i and iex with recalculating the free energies
1090 ! following temperature changes.
1091 ene_iex_iex=remd_ene(0,iex)
1092 ene_i_i=remd_ene(0,i)
1093 !o write (iout,*) "rescaling weights with temperature",
1095 call rescale_weights(remd_t_bath(i))
1097 call sum_energy(remd_ene(0,iex),.false.)
1098 ene_iex_i=remd_ene(0,iex)
1099 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1100 ! call sum_energy(remd_ene(0,i),.false.)
1101 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1102 ! write (iout,*) "rescaling weights with temperature",
1103 ! & remd_t_bath(iex)
1104 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1105 ! write (iout,*) "ERROR: inconsistent energies:",i,
1106 ! & ene_i_i,remd_ene(0,i)
1108 call rescale_weights(remd_t_bath(iex))
1109 call sum_energy(remd_ene(0,i),.false.)
1110 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1111 ene_i_iex=remd_ene(0,i)
1112 ! call sum_energy(remd_ene(0,iex),.false.)
1113 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1114 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1115 ! & ene_iex_iex,remd_ene(0,iex)
1117 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1118 ! write (iout,*) "i",i," iex",iex
1119 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1120 !d & " ene_i_iex",ene_i_iex,
1121 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1122 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1123 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1125 !d write(iout,*) 'delta',delta
1126 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1127 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1128 ! & (remd_t_bath(i)*remd_t_bath(iex))
1129 if (delta .gt. 50.0d0) then
1134 if (i_dir.eq.1.or.i_dir.eq.3) &
1135 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1136 if (i_dir.eq.2.or.i_dir.eq.3) &
1137 iremd_tot_usa(int(i2set(i-1)))= &
1138 iremd_tot_usa(int(i2set(i-1)))+1
1139 xxx=ran_number(0.0d0,1.0d0)
1140 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1141 if (delta .gt. xxx) then
1143 remd_t_bath(i)=remd_t_bath(iex)
1144 remd_t_bath(iex)=tmp
1147 iremd_iset(i)=iremd_iset(iex)
1148 iremd_iset(iex)=itmp
1150 remd_ene(0,i)=ene_i_iex
1151 remd_ene(0,iex)=ene_iex_i
1153 if (i_dir.eq.1.or.i_dir.eq.3) &
1154 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1157 i2rep(i-1)=i2rep(iex-1)
1160 if (i_dir.eq.2.or.i_dir.eq.3) &
1161 iremd_acc_usa(int(i2set(i-1)))= &
1162 iremd_acc_usa(int(i2set(i-1)))+1
1165 i2set(i-1)=i2set(iex-1)
1168 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1169 i_index(i_temp,i_mult,i_iset,i_mset)= &
1170 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1171 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1174 remd_ene(0,iex)=ene_iex_iex
1175 remd_ene(0,i)=ene_i_i
1176 remd_ene(20,iex)=econstr_temp_iex
1177 remd_ene(20,i)=econstr_temp_i
1181 !d do il1=1,mset(il)
1184 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1197 !-------------------------------------
1198 write (iout,*) "NREP",nrep
1200 if(iremd_tot(i).ne.0) &
1201 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1202 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1207 if(iremd_tot_usa(i).ne.0) &
1208 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1209 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1215 !d write (iout,'(a6,100i4)') "ifirst",
1216 !d & (ifirst(i),i=1,remd_m(1))
1218 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1219 !d & (nupa(i,il),i=1,nupa(0,il))
1220 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1221 !d & (ndowna(i,il),i=1,ndowna(0,il))
1226 !d write (iout,*) "Before scatter"
1229 if (me.eq.king) then
1230 write (iout,*) "t_bath before scatter",remd_t_bath
1234 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1235 t_bath,1,mpi_double_precision,king,&
1237 !d write (iout,*) "After scatter"
1240 call mpi_scatter(iremd_iset,1,mpi_integer,&
1241 iset,1,mpi_integer,king,&
1245 if (me.eq.king .or. .not. out1file) then
1246 write(iout,*) 'REMD scatter time=',time07-time06
1250 call mpi_scatter(elowi,1,mpi_double_precision,&
1251 elow,1,mpi_double_precision,king,&
1253 call mpi_scatter(ehighi,1,mpi_double_precision,&
1254 ehigh,1,mpi_double_precision,king,&
1257 call rescale_weights(t_bath)
1258 !o write (iout,*) "Processor",me,
1259 !o & " rescaling weights with temperature",t_bath
1261 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
1263 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
1267 !de write(iout,*) 'REMD after',me,t_bath
1269 if (me.eq.king .or. .not. out1file) then
1270 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1271 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1272 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1273 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1274 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1275 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1276 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1277 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1278 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1284 if (restart1file) then
1285 if (me.eq.king .or. .not. out1file) &
1286 write(iout,*) 'writing restart at the end of run'
1287 call write1rst(i_index)
1290 if (traj1file) call write1traj
1292 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1293 !deb & icache_all,1,mpi_integer,king,
1294 !deb & CG_COMM,ierr)
1295 !deb write(iout,'(a40,8000i8)')
1296 !deb & ' ntwx_cache after traj1file at the end',
1297 !deb & (icache_all(i),i=1,nodes)
1302 t_MD=MPI_Wtime()-tt0
1306 if (me.eq.king .or. .not. out1file) then
1307 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1309 'MD calculations setup:',t_MDsetup,&
1310 'Energy & gradient evaluation:',t_enegrad,&
1311 'Stochastic MD setup:',t_langsetup,&
1312 'Stochastic MD step setup:',t_sdsetup,&
1314 write (iout,'(/28(1h=),a25,27(1h=))') &
1315 ' End of MD calculation '
1317 !el common /przechowalnia/
1318 ! deallocate(d_restart1)
1319 ! deallocate(d_restart2)
1323 end subroutine MREMD
1324 !-----------------------------------------------------------------------------
1325 subroutine write1rst(i_index)
1328 ! implicit real*8 (a-h,o-z)
1329 ! include 'DIMENSIONS'
1331 ! include 'COMMON.MD'
1332 ! include 'COMMON.IOUNITS'
1333 ! include 'COMMON.REMD'
1334 ! include 'COMMON.SETUP'
1335 ! include 'COMMON.CHAIN'
1336 ! include 'COMMON.SBRIDGE'
1337 ! include 'COMMON.INTERACT'
1339 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1340 !el d_restart2(3,2*nres*maxprocs)
1341 real(kind=4) :: r_d(3,0:2*nres)
1342 real(kind=4) :: t5_restart1(5)
1343 integer :: iret,itmp
1344 integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1345 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1346 !el common /przechowalnia/ d_restart1,d_restart2
1347 integer :: i,j,il,il1,ierr,ixdrf
1352 t5_restart1(4)=t_bath
1353 t5_restart1(5)=Uconst
1355 call mpi_gather(t5_restart1,5,mpi_real,&
1356 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1364 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1365 d_restart1,3*2*nres+3,mpi_real,king,&
1377 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1378 d_restart2,3*2*nres+3,mpi_real,king,&
1383 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1385 call xdrfint_(ixdrf, i2rep(i), iret)
1388 call xdrfint_(ixdrf, ifirst(i), iret)
1392 call xdrfint_(ixdrf, nupa(i,il), iret)
1396 call xdrfint_(ixdrf, ndowna(i,il), iret)
1402 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1409 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1416 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1422 call xdrfint_(ixdrf, nset, iret)
1424 call xdrfint_(ixdrf,mset(i), iret)
1427 call xdrfint_(ixdrf,i2set(i), iret)
1433 itmp=i_index(i,j,il,il1)
1434 call xdrfint_(ixdrf,itmp, iret)
1441 call xdrfclose_(ixdrf, iret)
1443 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1445 call xdrfint(ixdrf, i2rep(i), iret)
1448 call xdrfint(ixdrf, ifirst(i), iret)
1452 call xdrfint(ixdrf, nupa(i,il), iret)
1456 call xdrfint(ixdrf, ndowna(i,il), iret)
1462 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1469 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1476 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1483 call xdrfint(ixdrf, nset, iret)
1485 call xdrfint(ixdrf,mset(i), iret)
1488 call xdrfint(ixdrf,i2set(i), iret)
1494 itmp=i_index(i,j,il,il1)
1495 call xdrfint(ixdrf,itmp, iret)
1502 call xdrfclose(ixdrf, iret)
1506 end subroutine write1rst
1507 !-----------------------------------------------------------------------------
1508 subroutine write1traj
1510 ! implicit real*8 (a-h,o-z)
1511 ! include 'DIMENSIONS'
1513 ! include 'COMMON.MD'
1514 ! include 'COMMON.IOUNITS'
1515 ! include 'COMMON.REMD'
1516 ! include 'COMMON.SETUP'
1517 ! include 'COMMON.CHAIN'
1518 ! include 'COMMON.SBRIDGE'
1519 ! include 'COMMON.INTERACT'
1521 real(kind=4) :: t5_restart1(5)
1522 integer :: iret,itmp
1523 real(kind=4) :: xcoord(3,2*nres+2),prec
1524 real(kind=4) :: r_qfrag(50),r_qpair(100)
1525 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1526 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1527 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1528 p_uscdiff(100*Nprocs) !(100*maxprocs)
1529 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1530 real(kind=4) :: r_c(3,2*nres+2)
1531 !el common /przechowalnia/ p_c
1533 integer :: i,j,il,ierr,ii,ixdrf
1535 call mpi_bcast(ii_write,1,mpi_integer,&
1539 print *,'traj1file',me,ii_write,ntwx_cache
1543 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1545 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1548 ! write (iout,*) "before gather write1traj: from node",ii
1550 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1552 t5_restart1(1)=totT_cache(ii)
1553 t5_restart1(2)=EK_cache(ii)
1554 t5_restart1(3)=potE_cache(ii)
1555 t5_restart1(4)=t_bath_cache(ii)
1556 t5_restart1(5)=Uconst_cache(ii)
1557 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1559 call mpi_gather(t5_restart1,5,mpi_real,&
1560 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1562 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1565 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1566 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1569 r_qfrag(i)=qfrag_cache(i,ii)
1572 r_qpair(i)=qpair_cache(i,ii)
1575 r_utheta(i)=utheta_cache(i,ii)
1576 r_ugamma(i)=ugamma_cache(i,ii)
1577 r_uscdiff(i)=uscdiff_cache(i,ii)
1580 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1581 p_qfrag,nfrag,mpi_real,king,&
1583 call mpi_gather(r_qpair,npair,mpi_real,&
1584 p_qpair,npair,mpi_real,king,&
1586 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1587 p_utheta,nfrag_back,mpi_real,king,&
1589 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1590 p_ugamma,nfrag_back,mpi_real,king,&
1592 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1593 p_uscdiff,nfrag_back,mpi_real,king,&
1597 write (iout,*) "p_qfrag"
1599 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1601 write (iout,*) "p_qpair"
1603 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1609 r_c(j,i)=c_cache(j,i,ii)
1613 call mpi_gather(r_c,3*2*nres,mpi_real,&
1614 p_c,3*2*nres,mpi_real,king,&
1620 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1621 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1622 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1623 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1624 call xdrfint_(ixdrf, nss, iret)
1627 call xdrfint(ixdrf, idssb(j)+nres, iret)
1628 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1630 call xdrfint_(ixdrf, ihpb(j), iret)
1631 call xdrfint_(ixdrf, jhpb(j), iret)
1634 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1635 call xdrfint_(ixdrf, iset_restart1(il), iret)
1637 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1640 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1643 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1644 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1645 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1650 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1655 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1659 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1663 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1664 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1665 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1666 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1667 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1668 call xdrfint(ixdrf, nss, iret)
1671 call xdrfint(ixdrf, idssb(j)+nres, iret)
1672 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1674 call xdrfint(ixdrf, ihpb(j), iret)
1675 call xdrfint(ixdrf, jhpb(j), iret)
1678 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1679 call xdrfint(ixdrf, iset_restart1(il), iret)
1681 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1684 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1687 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1688 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1689 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1694 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1699 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1703 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1709 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1711 if(me.eq.king) call xdrfclose(ixdrf, iret)
1713 do i=1,ntwx_cache-ii_write
1715 totT_cache(i)=totT_cache(ii_write+i)
1716 EK_cache(i)=EK_cache(ii_write+i)
1717 potE_cache(i)=potE_cache(ii_write+i)
1718 t_bath_cache(i)=t_bath_cache(ii_write+i)
1719 Uconst_cache(i)=Uconst_cache(ii_write+i)
1720 iset_cache(i)=iset_cache(ii_write+i)
1723 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1726 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1729 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1730 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1731 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1736 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1740 ntwx_cache=ntwx_cache-ii_write
1742 end subroutine write1traj
1743 !-----------------------------------------------------------------------------
1744 subroutine read1restart(i_index)
1746 ! implicit real*8 (a-h,o-z)
1747 ! include 'DIMENSIONS'
1749 ! include 'COMMON.MD'
1750 ! include 'COMMON.IOUNITS'
1751 ! include 'COMMON.REMD'
1752 ! include 'COMMON.SETUP'
1753 ! include 'COMMON.CHAIN'
1754 ! include 'COMMON.SBRIDGE'
1755 ! include 'COMMON.INTERACT'
1756 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1757 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1758 integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1759 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1760 !el common /przechowalnia/ d_restart1
1761 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1763 write (*,*) "Processor",me," called read1restart"
1766 open(irest2,file=mremd_rst_name,status='unknown')
1767 read(irest2,*,err=334) i
1768 write(iout,*) "Reading old rst in ASCI format"
1770 call read1restart_old
1774 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1777 call xdrfint_(ixdrf, i2rep(i), iret)
1780 call xdrfint_(ixdrf, ifirst(i), iret)
1783 call xdrfint_(ixdrf, nupa(0,il), iret)
1785 call xdrfint_(ixdrf, nupa(i,il), iret)
1788 call xdrfint_(ixdrf, ndowna(0,il), iret)
1790 call xdrfint_(ixdrf, ndowna(i,il), iret)
1795 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1799 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1802 call xdrfint(ixdrf, i2rep(i), iret)
1805 call xdrfint(ixdrf, ifirst(i), iret)
1808 call xdrfint(ixdrf, nupa(0,il), iret)
1810 call xdrfint(ixdrf, nupa(i,il), iret)
1813 call xdrfint(ixdrf, ndowna(0,il), iret)
1815 call xdrfint(ixdrf, ndowna(i,il), iret)
1820 call xdrffloat(ixdrf, t_restart1(j,il), iret)
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)
1839 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1841 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1847 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1848 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1858 ! read(irest2,'(3e15.5)')
1859 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1862 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1864 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1870 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1871 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1882 call xdrfint_(ixdrf, nset, iret)
1884 call xdrfint_(ixdrf,mset(i), iret)
1887 call xdrfint_(ixdrf,i2set(i), iret)
1893 call xdrfint_(ixdrf,itmp, iret)
1894 i_index(i,j,il,il1)=itmp
1902 call xdrfint(ixdrf, nset, iret)
1904 call xdrfint(ixdrf,mset(i), iret)
1907 call xdrfint(ixdrf,i2set(i), iret)
1913 call xdrfint(ixdrf,itmp, iret)
1914 i_index(i,j,il,il1)=itmp
1921 call mpi_scatter(i2set,1,mpi_integer,&
1922 iset,1,mpi_integer,king,&
1927 if(me.eq.king) close(irest2)
1929 end subroutine read1restart
1930 !-----------------------------------------------------------------------------
1931 subroutine read1restart_old
1933 ! implicit real*8 (a-h,o-z)
1934 ! include 'DIMENSIONS'
1936 ! include 'COMMON.MD'
1937 ! include 'COMMON.IOUNITS'
1938 ! include 'COMMON.REMD'
1939 ! include 'COMMON.SETUP'
1940 ! include 'COMMON.CHAIN'
1941 ! include 'COMMON.SBRIDGE'
1942 ! include 'COMMON.INTERACT'
1943 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1944 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1945 !el common /przechowalnia/ d_restart1
1947 integer :: i,j,il,ierr
1950 open(irest2,file=mremd_rst_name,status='unknown')
1951 read (irest2,*) (i2rep(i),i=0,nodes-1)
1952 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1954 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1955 read (irest2,*) ndowna(0,il),&
1956 (ndowna(i,il),i=1,ndowna(0,il))
1959 read(irest2,*) (t_restart1(j,il),j=1,4)
1962 call mpi_scatter(t_restart1,5,mpi_real,&
1963 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1967 t_bath=t5_restart1(4)
1972 read(irest2,'(3e15.5)') &
1973 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
1977 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1978 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1988 read(irest2,'(3e15.5)') &
1989 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
1993 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1994 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2000 if(me.eq.king) close(irest2)
2002 end subroutine read1restart_old
2003 !----------------------------------------------------------------
2004 subroutine alloc_MREMD_arrays
2006 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2007 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
2009 ! common /remdcommon/ in io: read_REMDpar
2010 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2011 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2012 ! common /remdrestart/
2013 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2015 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2016 allocate(ifirst(0:nodes)) !(maxprocs)
2017 allocate(nupa(0:nodes,0:2*nodes))
2018 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2019 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2020 allocate(iset_restart1(nodes)) !(maxprocs)
2021 ! common /traj1cache/
2022 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2023 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2024 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2025 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2026 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2027 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2028 allocate(utheta_cache(nfrag_back,max_cache_traj))
2029 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2030 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2031 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2034 end subroutine alloc_MREMD_arrays
2035 !-----------------------------------------------------------------------------
2036 !-----------------------------------------------------------------------------