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
351 stdfp=dsqrt(2*Rb*t_bath/d_time)
353 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
356 ! print *,'irep',me,t_bath
358 if (me.eq.king .or. .not. out1file) &
359 write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
360 call rescale_weights(t_bath)
364 !------copy MD--------------
365 ! The driver for molecular dynamics subroutines
366 !------------------------------------------------
372 if(me.eq.king.or..not.out1file) &
373 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
379 ! Determine the inverse of the inertia matrix.
380 call setup_MD_matrices
384 if (me.eq.king .or. .not. out1file) &
385 write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
386 stdfp=dsqrt(2*Rb*t_bath/d_time)
388 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
390 call rescale_weights(t_bath)
394 t_MDsetup = MPI_Wtime()-tt0
396 t_MDsetup = tcpu()-tt0
399 ! Entering the MD loop
405 if (lang.eq.2 .or. lang.eq.3) then
409 call sd_verlet_p_setup
411 call sd_verlet_ciccotti_setup
415 pfric0_mat(i,j,0)=pfric_mat(i,j)
416 afric0_mat(i,j,0)=afric_mat(i,j)
417 vfric0_mat(i,j,0)=vfric_mat(i,j)
418 prand0_mat(i,j,0)=prand_mat(i,j)
419 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
420 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
425 flag_stoch(i)=.false.
429 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
431 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
435 else if (lang.eq.1 .or. lang.eq.4) then
439 if (me.eq.king .or. .not. out1file) &
440 write(iout,*) 'Setup time',time00-walltime
443 t_langsetup=MPI_Wtime()-tt0
446 t_langsetup=tcpu()-tt0
452 do while(.not.end_of_run)
454 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
455 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
457 if (lang.gt.0 .and. surfarea .and. &
458 mod(itime,reset_fricmat).eq.0) then
459 if (lang.eq.2 .or. lang.eq.3) then
463 call sd_verlet_p_setup
465 call sd_verlet_ciccotti_setup
469 pfric0_mat(i,j,0)=pfric_mat(i,j)
470 afric0_mat(i,j,0)=afric_mat(i,j)
471 vfric0_mat(i,j,0)=vfric_mat(i,j)
472 prand0_mat(i,j,0)=prand_mat(i,j)
473 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
474 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
479 flag_stoch(i)=.false.
482 else if (lang.eq.1 .or. lang.eq.4) then
485 write (iout,'(a,i10)') &
486 "Friction matrix reset based on surface area, itime",itime
488 if (reset_vel .and. tbf .and. lang.eq.0 &
489 .and. mod(itime,count_reset_vel).eq.0) then
491 if (me.eq.king .or. .not. out1file) &
492 write(iout,'(a,f20.2)') &
493 "Velocities reset to random values, time",totT
496 d_t_old(j,i)=d_t(j,i)
500 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
504 d_t(j,0)=d_t(j,0)-vcm(j)
507 kinetic_T=2.0d0/(dimen3*Rb)*EK
508 scalfac=dsqrt(T_bath/kinetic_T)
509 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
512 d_t_old(j,i)=scalfac*d_t(j,i)
518 ! Time-reversible RESPA algorithm
519 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
520 call RESPA_step(itime)
522 ! Variable time step algorithm.
523 call velverlet_step(itime)
527 call brown_step(itime)
529 print *,"Brown dynamics not here!"
531 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
537 if (mod(itime,ntwe).eq.0) call statout(itime)
539 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
540 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
542 call hairpin(.true.,nharp,iharp)
543 call secondary2(.true.)
544 call pdbout(potE,tytul,ipdb)
549 if (mod(itime,ntwx).eq.0.and.traj1file) then
550 if(ntwx_cache.lt.max_cache_traj_use) then
551 ntwx_cache=ntwx_cache+1
553 if (max_cache_traj_use.ne.1) &
554 print *,itime,"processor ",me," over cache ",ntwx_cache
557 totT_cache(i)=totT_cache(i+1)
558 EK_cache(i)=EK_cache(i+1)
559 potE_cache(i)=potE_cache(i+1)
560 t_bath_cache(i)=t_bath_cache(i+1)
561 Uconst_cache(i)=Uconst_cache(i+1)
562 iset_cache(i)=iset_cache(i+1)
565 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
568 qpair_cache(ii,i)=qpair_cache(ii,i+1)
571 utheta_cache(ii,i)=utheta_cache(ii,i+1)
572 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
573 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
579 c_cache(j,ii,i)=c_cache(j,ii,i+1)
585 totT_cache(ntwx_cache)=totT
586 EK_cache(ntwx_cache)=EK
587 potE_cache(ntwx_cache)=potE
588 t_bath_cache(ntwx_cache)=t_bath
589 Uconst_cache(ntwx_cache)=Uconst
590 iset_cache(ntwx_cache)=iset
593 qfrag_cache(i,ntwx_cache)=qfrag(i)
596 qpair_cache(i,ntwx_cache)=qpair(i)
599 utheta_cache(i,ntwx_cache)=utheta(i)
600 ugamma_cache(i,ntwx_cache)=ugamma(i)
601 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
606 c_cache(j,i,ntwx_cache)=c(j,i)
611 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
612 .and..not.restart1file) then
615 open(irest1,file=mremd_rst_name,status='unknown')
616 write (irest1,*) "i2rep"
617 write (irest1,*) (i2rep(i),i=0,nodes-1)
618 write (irest1,*) "ifirst"
619 write (irest1,*) (ifirst(i),i=1,remd_m(1))
621 write (irest1,*) "nupa",il
622 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
623 write (irest1,*) "ndowna",il
624 write (irest1,*) ndowna(0,il),&
625 (ndowna(i,il),i=1,ndowna(0,il))
628 write (irest1,*) "nset"
629 write (irest1,*) nset
630 write (irest1,*) "mset"
631 write (irest1,*) (mset(i),i=1,nset)
632 write (irest1,*) "i2set"
633 write (irest1,*) (i2set(i),i=0,nodes-1)
634 write (irest1,*) "i_index"
638 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
646 open(irest2,file=rest2name,status='unknown')
647 write(irest2,*) totT,EK,potE,totE,t_bath
649 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
652 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
655 write (irest2,*) iset
662 ! forced synchronization
663 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
664 .and. .not. mremdsync) then
666 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
668 call mpi_recv(itime_master, 1, MPI_INTEGER,&
669 0,101,CG_COMM, status, ierr)
670 call mpi_barrier(CG_COMM, ierr)
671 !deb if (out1file.or.traj1file) then
672 !deb call mpi_gather(itime,1,mpi_integer,
673 !deb & icache_all,1,mpi_integer,king,
676 call mpi_gather(ntwx_cache,1,mpi_integer,&
677 icache_all,1,mpi_integer,king,&
680 write(iout,*) 'REMD synchro at',itime_master,itime
681 if (itime_master.ge.n_timestep .or. ovrtim()) &
683 !time call flush(iout)
688 if ((mod(itime,nstex).eq.0.and.me.eq.king &
689 .or.end_of_run.and.me.eq.king ) &
690 .and. .not. mremdsync ) then
693 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
694 CG_COMM, ireqi(i), ierr)
695 !d write(iout,*) 'REMD synchro with',i
698 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
699 call mpi_barrier(CG_COMM, ierr)
701 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
702 if (out1file.or.traj1file) then
703 !deb call mpi_gather(itime,1,mpi_integer,
704 !deb & itime_all,1,mpi_integer,king,
706 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
707 !deb & (itime_all(i),i=1,nodes)
709 !deb imin_itime=itime_all(1)
711 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
713 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
714 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
715 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
716 call mpi_gather(ntwx_cache,1,mpi_integer,&
717 icache_all,1,mpi_integer,king,&
719 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
720 ! & (icache_all(i),i=1,nodes)
721 ii_write=icache_all(1)
723 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
725 ! write(iout,*) "MIN ii_write=",ii_write
728 !time call flush(iout)
730 if(mremdsync .and. mod(itime,nstex).eq.0) then
732 if (me.eq.king .or. .not. out1file) &
733 write(iout,*) 'REMD synchro at',itime
736 call mpi_gather(ntwx_cache,1,mpi_integer,&
737 icache_all,1,mpi_integer,king,&
740 write(iout,'(a19,8000i8)') ' ntwx_cache',&
741 (icache_all(i),i=1,nodes)
742 ii_write=icache_all(1)
744 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
746 write(iout,*) "MIN ii_write=",ii_write
752 ! Update the time safety limiy
753 if (time001-time00.gt.safety) then
754 safety=time001-time00+600
755 write (iout,*) "****** SAFETY increased to",safety," s"
757 if (ovrtim()) end_of_run=.true.
759 if(synflag.and..not.end_of_run) then
763 write(iout,*) 'REMD before',me,t_bath
765 ! call mpi_gather(t_bath,1,mpi_double_precision,
766 ! & remd_t_bath,1,mpi_double_precision,king,
768 potEcomp(n_ene+1)=t_bath
770 potEcomp(n_ene+2)=iset
771 if (iset.lt.nset) then
775 potEcomp(n_ene+3)=Uconst
782 potEcomp(n_ene+4)=Uconst
786 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
787 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
790 call mpi_gather(elow,1,mpi_double_precision,&
791 elowi,1,mpi_double_precision,king,&
793 call mpi_gather(ehigh,1,mpi_double_precision,&
794 ehighi,1,mpi_double_precision,king,&
799 if (me.eq.king .or. .not. out1file) then
800 write(iout,*) 'REMD gather times=',time03-time01 &
804 if (restart1file) call write1rst(i_index)
807 if (me.eq.king .or. .not. out1file) then
808 write(iout,*) 'REMD writing rst time=',time04-time03
811 if (traj1file) call write1traj
813 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
814 !deb & icache_all,1,mpi_integer,king,
816 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
817 !deb & (icache_all(i),i=1,nodes)
822 if (me.eq.king .or. .not. out1file) then
823 write(iout,*) 'REMD writing traj time=',time05-time04
830 remd_t_bath(i)=remd_ene(n_ene+1,i)
831 iremd_iset(i)=remd_ene(n_ene+2,i)
835 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
837 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
841 write(iout,*) 'REMD exchange temp,ene'
843 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
844 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
848 !-------------------------------------
850 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
853 write (iout,*) "remd_m(1)",remd_m(1)
855 i=ifirst(iran_num(1,remd_m(1)))
862 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
864 if(i.gt.0.and.nupa(0,i).gt.0) then
866 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
868 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
870 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
872 ! do while (iex.eq.i)
873 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
874 iex=nupa(iran_num(1,int(nupa(0,i))),i)
876 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
878 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
880 ! Swap temperatures between conformations i and iex with recalculating the free energies
881 ! following temperature changes.
882 ene_iex_iex=remd_ene(0,iex)
883 ene_i_i=remd_ene(0,i)
884 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
885 ! & " iex",iex," ene_iex_iex",ene_iex_iex
886 ! write (iout,*) "rescaling weights with temperature",
889 call rescale_weights(remd_t_bath(i))
891 ! write (iout,*) "0,iex",remd_t_bath(i)
892 ! call enerprint(remd_ene(0,iex))
894 call sum_energy(remd_ene(0,iex),.false.)
895 ene_iex_i=remd_ene(0,iex)
896 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
898 ! write (iout,*) "0,i",remd_t_bath(i)
899 ! call enerprint(remd_ene(0,i))
901 call sum_energy(remd_ene(0,i),.false.)
902 ! write (iout,*) "ene_i_i",remd_ene(0,i)
904 ! write (iout,*) "rescaling weights with temperature",
906 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
907 write (iout,*) "ERROR: inconsistent energies:",i,&
908 ene_i_i,remd_ene(0,i)
910 call rescale_weights(remd_t_bath(iex))
912 ! write (iout,*) "0,i",remd_t_bath(iex)
913 ! call enerprint(remd_ene(0,i))
915 call sum_energy(remd_ene(0,i),.false.)
916 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
918 ene_i_iex=remd_ene(0,i)
920 ! write (iout,*) "0,iex",remd_t_bath(iex)
921 ! call enerprint(remd_ene(0,iex))
923 call sum_energy(remd_ene(0,iex),.false.)
924 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
925 write (iout,*) "ERROR: inconsistent energies:",iex,&
926 ene_iex_iex,remd_ene(0,iex)
928 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
929 ! write (iout,*) "i",i," iex",iex
930 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
931 ! & " ene_i_iex",ene_i_iex,
932 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
934 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
935 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
937 ! write(iout,*) 'delta',delta
938 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
939 ! & (remd_ene(i)-remd_ene(iex))/Rb/
940 ! & (remd_t_bath(i)*remd_t_bath(iex))
942 if (delta .gt. 50.0d0) then
948 else if (delta.lt.-50.0d0) then
957 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
958 xxx=ran_number(0.0d0,1.0d0)
959 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
961 if (delta .gt. xxx) then
963 remd_t_bath(i)=remd_t_bath(iex)
965 remd_ene(0,i)=ene_i_iex
966 remd_ene(0,iex)=ene_iex_i
972 ehighi(i)=ehighi(iex)
979 nupa(k,i)=nupa(k,iex)
982 ndowna(k,i)=ndowna(k,iex)
986 if (ifirst(il).eq.i) ifirst(il)=iex
988 if (nupa(k,il).eq.i) then
990 elseif (nupa(k,il).eq.iex) then
995 if (ndowna(k,il).eq.i) then
997 elseif (ndowna(k,il).eq.iex) then
1003 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1005 i2rep(i-1)=i2rep(iex-1)
1008 ! write(iout,*) 'exchange',i,iex
1009 ! write (iout,'(a8,100i4)') "@ ifirst",
1010 ! & (ifirst(k),k=1,remd_m(1))
1012 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1013 ! & (nupa(k,il),k=1,nupa(0,il))
1014 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1015 ! & (ndowna(k,il),k=1,ndowna(0,il))
1020 remd_ene(0,iex)=ene_iex_iex
1021 remd_ene(0,i)=ene_i_i
1027 !d write (iout,*) "exchange completed"
1031 !d write(iout,*) "########",ii
1033 i_temp=iran_num(1,nrep)
1034 i_mult=iran_num(1,remd_m(i_temp))
1035 i_iset=iran_num(1,nset)
1036 i_mset=iran_num(1,mset(i_iset))
1037 i=i_index(i_temp,i_mult,i_iset,i_mset)
1039 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1042 !d write(iout,*) "i_dir=",i_dir
1044 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1047 i_mult1=iran_num(1,remd_m(i_temp1))
1049 i_mset1=iran_num(1,mset(i_iset1))
1050 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1052 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1055 i_mult1=iran_num(1,remd_m(i_temp1))
1057 i_mset1=iran_num(1,mset(i_iset1))
1058 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1059 econstr_temp_i=remd_ene(20,i)
1060 econstr_temp_iex=remd_ene(20,iex)
1061 remd_ene(20,i)=remd_ene(n_ene+3,i)
1062 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1064 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1067 i_mult1=iran_num(1,remd_m(i_temp1))
1069 i_mset1=iran_num(1,mset(i_iset1))
1070 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1071 econstr_temp_i=remd_ene(20,i)
1072 econstr_temp_iex=remd_ene(20,iex)
1073 remd_ene(20,i)=remd_ene(n_ene+3,i)
1074 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1080 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1083 ! Swap temperatures between conformations i and iex with recalculating the free energies
1084 ! following temperature changes.
1085 ene_iex_iex=remd_ene(0,iex)
1086 ene_i_i=remd_ene(0,i)
1087 !o write (iout,*) "rescaling weights with temperature",
1089 call rescale_weights(remd_t_bath(i))
1091 call sum_energy(remd_ene(0,iex),.false.)
1092 ene_iex_i=remd_ene(0,iex)
1093 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1094 ! call sum_energy(remd_ene(0,i),.false.)
1095 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1096 ! write (iout,*) "rescaling weights with temperature",
1097 ! & remd_t_bath(iex)
1098 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1099 ! write (iout,*) "ERROR: inconsistent energies:",i,
1100 ! & ene_i_i,remd_ene(0,i)
1102 call rescale_weights(remd_t_bath(iex))
1103 call sum_energy(remd_ene(0,i),.false.)
1104 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1105 ene_i_iex=remd_ene(0,i)
1106 ! call sum_energy(remd_ene(0,iex),.false.)
1107 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1108 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1109 ! & ene_iex_iex,remd_ene(0,iex)
1111 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1112 ! write (iout,*) "i",i," iex",iex
1113 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1114 !d & " ene_i_iex",ene_i_iex,
1115 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1116 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1117 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1119 !d write(iout,*) 'delta',delta
1120 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1121 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1122 ! & (remd_t_bath(i)*remd_t_bath(iex))
1123 if (delta .gt. 50.0d0) then
1128 if (i_dir.eq.1.or.i_dir.eq.3) &
1129 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1130 if (i_dir.eq.2.or.i_dir.eq.3) &
1131 iremd_tot_usa(int(i2set(i-1)))= &
1132 iremd_tot_usa(int(i2set(i-1)))+1
1133 xxx=ran_number(0.0d0,1.0d0)
1134 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1135 if (delta .gt. xxx) then
1137 remd_t_bath(i)=remd_t_bath(iex)
1138 remd_t_bath(iex)=tmp
1141 iremd_iset(i)=iremd_iset(iex)
1142 iremd_iset(iex)=itmp
1144 remd_ene(0,i)=ene_i_iex
1145 remd_ene(0,iex)=ene_iex_i
1147 if (i_dir.eq.1.or.i_dir.eq.3) &
1148 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1151 i2rep(i-1)=i2rep(iex-1)
1154 if (i_dir.eq.2.or.i_dir.eq.3) &
1155 iremd_acc_usa(int(i2set(i-1)))= &
1156 iremd_acc_usa(int(i2set(i-1)))+1
1159 i2set(i-1)=i2set(iex-1)
1162 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1163 i_index(i_temp,i_mult,i_iset,i_mset)= &
1164 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1165 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1168 remd_ene(0,iex)=ene_iex_iex
1169 remd_ene(0,i)=ene_i_i
1170 remd_ene(20,iex)=econstr_temp_iex
1171 remd_ene(20,i)=econstr_temp_i
1175 !d do il1=1,mset(il)
1178 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1191 !-------------------------------------
1192 write (iout,*) "NREP",nrep
1194 if(iremd_tot(i).ne.0) &
1195 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1196 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1201 if(iremd_tot_usa(i).ne.0) &
1202 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1203 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1209 !d write (iout,'(a6,100i4)') "ifirst",
1210 !d & (ifirst(i),i=1,remd_m(1))
1212 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1213 !d & (nupa(i,il),i=1,nupa(0,il))
1214 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1215 !d & (ndowna(i,il),i=1,ndowna(0,il))
1220 !d write (iout,*) "Before scatter"
1223 if (me.eq.king) then
1224 write (iout,*) "t_bath before scatter",remd_t_bath
1228 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1229 t_bath,1,mpi_double_precision,king,&
1231 !d write (iout,*) "After scatter"
1234 call mpi_scatter(iremd_iset,1,mpi_integer,&
1235 iset,1,mpi_integer,king,&
1239 if (me.eq.king .or. .not. out1file) then
1240 write(iout,*) 'REMD scatter time=',time07-time06
1244 call mpi_scatter(elowi,1,mpi_double_precision,&
1245 elow,1,mpi_double_precision,king,&
1247 call mpi_scatter(ehighi,1,mpi_double_precision,&
1248 ehigh,1,mpi_double_precision,king,&
1251 call rescale_weights(t_bath)
1252 !o write (iout,*) "Processor",me,
1253 !o & " rescaling weights with temperature",t_bath
1255 stdfp=dsqrt(2*Rb*t_bath/d_time)
1257 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1260 !de write(iout,*) 'REMD after',me,t_bath
1262 if (me.eq.king .or. .not. out1file) then
1263 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1264 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1265 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1266 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1267 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1268 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1269 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1270 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1271 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1277 if (restart1file) then
1278 if (me.eq.king .or. .not. out1file) &
1279 write(iout,*) 'writing restart at the end of run'
1280 call write1rst(i_index)
1283 if (traj1file) call write1traj
1285 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1286 !deb & icache_all,1,mpi_integer,king,
1287 !deb & CG_COMM,ierr)
1288 !deb write(iout,'(a40,8000i8)')
1289 !deb & ' ntwx_cache after traj1file at the end',
1290 !deb & (icache_all(i),i=1,nodes)
1295 t_MD=MPI_Wtime()-tt0
1299 if (me.eq.king .or. .not. out1file) then
1300 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1302 'MD calculations setup:',t_MDsetup,&
1303 'Energy & gradient evaluation:',t_enegrad,&
1304 'Stochastic MD setup:',t_langsetup,&
1305 'Stochastic MD step setup:',t_sdsetup,&
1307 write (iout,'(/28(1h=),a25,27(1h=))') &
1308 ' End of MD calculation '
1310 !el common /przechowalnia/
1311 ! deallocate(d_restart1)
1312 ! deallocate(d_restart2)
1316 end subroutine MREMD
1317 !-----------------------------------------------------------------------------
1318 subroutine write1rst(i_index)
1321 ! implicit real*8 (a-h,o-z)
1322 ! include 'DIMENSIONS'
1324 ! include 'COMMON.MD'
1325 ! include 'COMMON.IOUNITS'
1326 ! include 'COMMON.REMD'
1327 ! include 'COMMON.SETUP'
1328 ! include 'COMMON.CHAIN'
1329 ! include 'COMMON.SBRIDGE'
1330 ! include 'COMMON.INTERACT'
1332 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1333 !el d_restart2(3,2*nres*maxprocs)
1334 real(kind=4) :: r_d(3,0:2*nres)
1335 real(kind=4) :: t5_restart1(5)
1336 integer :: iret,itmp
1337 integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1338 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1339 !el common /przechowalnia/ d_restart1,d_restart2
1340 integer :: i,j,il,il1,ierr,ixdrf
1345 t5_restart1(4)=t_bath
1346 t5_restart1(5)=Uconst
1348 call mpi_gather(t5_restart1,5,mpi_real,&
1349 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1357 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1358 d_restart1,3*2*nres+3,mpi_real,king,&
1370 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1371 d_restart2,3*2*nres+3,mpi_real,king,&
1376 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1378 call xdrfint_(ixdrf, i2rep(i), iret)
1381 call xdrfint_(ixdrf, ifirst(i), iret)
1385 call xdrfint_(ixdrf, nupa(i,il), iret)
1389 call xdrfint_(ixdrf, ndowna(i,il), iret)
1395 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1402 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1409 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1415 call xdrfint_(ixdrf, nset, iret)
1417 call xdrfint_(ixdrf,mset(i), iret)
1420 call xdrfint_(ixdrf,i2set(i), iret)
1426 itmp=i_index(i,j,il,il1)
1427 call xdrfint_(ixdrf,itmp, iret)
1434 call xdrfclose_(ixdrf, iret)
1436 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1438 call xdrfint(ixdrf, i2rep(i), iret)
1441 call xdrfint(ixdrf, ifirst(i), iret)
1445 call xdrfint(ixdrf, nupa(i,il), iret)
1449 call xdrfint(ixdrf, ndowna(i,il), iret)
1455 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1462 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1469 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1476 call xdrfint(ixdrf, nset, iret)
1478 call xdrfint(ixdrf,mset(i), iret)
1481 call xdrfint(ixdrf,i2set(i), iret)
1487 itmp=i_index(i,j,il,il1)
1488 call xdrfint(ixdrf,itmp, iret)
1495 call xdrfclose(ixdrf, iret)
1499 end subroutine write1rst
1500 !-----------------------------------------------------------------------------
1501 subroutine write1traj
1503 ! implicit real*8 (a-h,o-z)
1504 ! include 'DIMENSIONS'
1506 ! include 'COMMON.MD'
1507 ! include 'COMMON.IOUNITS'
1508 ! include 'COMMON.REMD'
1509 ! include 'COMMON.SETUP'
1510 ! include 'COMMON.CHAIN'
1511 ! include 'COMMON.SBRIDGE'
1512 ! include 'COMMON.INTERACT'
1514 real(kind=4) :: t5_restart1(5)
1515 integer :: iret,itmp
1516 real(kind=4) :: xcoord(3,2*nres+2),prec
1517 real(kind=4) :: r_qfrag(50),r_qpair(100)
1518 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1519 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1520 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1521 p_uscdiff(100*Nprocs) !(100*maxprocs)
1522 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1523 real(kind=4) :: r_c(3,2*nres+2)
1524 !el common /przechowalnia/ p_c
1526 integer :: i,j,il,ierr,ii,ixdrf
1528 call mpi_bcast(ii_write,1,mpi_integer,&
1532 print *,'traj1file',me,ii_write,ntwx_cache
1536 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1538 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1541 ! write (iout,*) "before gather write1traj: from node",ii
1543 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1545 t5_restart1(1)=totT_cache(ii)
1546 t5_restart1(2)=EK_cache(ii)
1547 t5_restart1(3)=potE_cache(ii)
1548 t5_restart1(4)=t_bath_cache(ii)
1549 t5_restart1(5)=Uconst_cache(ii)
1550 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1552 call mpi_gather(t5_restart1,5,mpi_real,&
1553 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1555 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1558 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1559 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1562 r_qfrag(i)=qfrag_cache(i,ii)
1565 r_qpair(i)=qpair_cache(i,ii)
1568 r_utheta(i)=utheta_cache(i,ii)
1569 r_ugamma(i)=ugamma_cache(i,ii)
1570 r_uscdiff(i)=uscdiff_cache(i,ii)
1573 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1574 p_qfrag,nfrag,mpi_real,king,&
1576 call mpi_gather(r_qpair,npair,mpi_real,&
1577 p_qpair,npair,mpi_real,king,&
1579 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1580 p_utheta,nfrag_back,mpi_real,king,&
1582 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1583 p_ugamma,nfrag_back,mpi_real,king,&
1585 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1586 p_uscdiff,nfrag_back,mpi_real,king,&
1590 write (iout,*) "p_qfrag"
1592 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1594 write (iout,*) "p_qpair"
1596 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1602 r_c(j,i)=c_cache(j,i,ii)
1606 call mpi_gather(r_c,3*2*nres,mpi_real,&
1607 p_c,3*2*nres,mpi_real,king,&
1613 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1614 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1615 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1616 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1617 call xdrfint_(ixdrf, nss, iret)
1620 call xdrfint(ixdrf, idssb(j)+nres, iret)
1621 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1623 call xdrfint_(ixdrf, ihpb(j), iret)
1624 call xdrfint_(ixdrf, jhpb(j), iret)
1627 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1628 call xdrfint_(ixdrf, iset_restart1(il), iret)
1630 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1633 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1636 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1637 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1638 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1643 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1648 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1652 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1656 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1657 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1658 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1659 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1660 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1661 call xdrfint(ixdrf, nss, iret)
1664 call xdrfint(ixdrf, idssb(j)+nres, iret)
1665 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1667 call xdrfint(ixdrf, ihpb(j), iret)
1668 call xdrfint(ixdrf, jhpb(j), iret)
1671 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1672 call xdrfint(ixdrf, iset_restart1(il), iret)
1674 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1677 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1680 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1681 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1682 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1687 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1692 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1696 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1702 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1704 if(me.eq.king) call xdrfclose(ixdrf, iret)
1706 do i=1,ntwx_cache-ii_write
1708 totT_cache(i)=totT_cache(ii_write+i)
1709 EK_cache(i)=EK_cache(ii_write+i)
1710 potE_cache(i)=potE_cache(ii_write+i)
1711 t_bath_cache(i)=t_bath_cache(ii_write+i)
1712 Uconst_cache(i)=Uconst_cache(ii_write+i)
1713 iset_cache(i)=iset_cache(ii_write+i)
1716 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1719 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1722 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1723 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1724 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1729 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1733 ntwx_cache=ntwx_cache-ii_write
1735 end subroutine write1traj
1736 !-----------------------------------------------------------------------------
1737 subroutine read1restart(i_index)
1739 ! implicit real*8 (a-h,o-z)
1740 ! include 'DIMENSIONS'
1742 ! include 'COMMON.MD'
1743 ! include 'COMMON.IOUNITS'
1744 ! include 'COMMON.REMD'
1745 ! include 'COMMON.SETUP'
1746 ! include 'COMMON.CHAIN'
1747 ! include 'COMMON.SBRIDGE'
1748 ! include 'COMMON.INTERACT'
1749 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1750 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1751 integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1752 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1753 !el common /przechowalnia/ d_restart1
1754 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1756 write (*,*) "Processor",me," called read1restart"
1759 open(irest2,file=mremd_rst_name,status='unknown')
1760 read(irest2,*,err=334) i
1761 write(iout,*) "Reading old rst in ASCI format"
1763 call read1restart_old
1767 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1770 call xdrfint_(ixdrf, i2rep(i), iret)
1773 call xdrfint_(ixdrf, ifirst(i), iret)
1776 call xdrfint_(ixdrf, nupa(0,il), iret)
1778 call xdrfint_(ixdrf, nupa(i,il), iret)
1781 call xdrfint_(ixdrf, ndowna(0,il), iret)
1783 call xdrfint_(ixdrf, ndowna(i,il), iret)
1788 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1792 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1795 call xdrfint(ixdrf, i2rep(i), iret)
1798 call xdrfint(ixdrf, ifirst(i), iret)
1801 call xdrfint(ixdrf, nupa(0,il), iret)
1803 call xdrfint(ixdrf, nupa(i,il), iret)
1806 call xdrfint(ixdrf, ndowna(0,il), iret)
1808 call xdrfint(ixdrf, ndowna(i,il), iret)
1813 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1818 call mpi_scatter(t_restart1,5,mpi_real,&
1819 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1823 t_bath=t5_restart1(4)
1828 ! read(irest2,'(3e15.5)')
1829 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1832 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1834 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1840 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1841 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1851 ! read(irest2,'(3e15.5)')
1852 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1855 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1857 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1863 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1864 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1875 call xdrfint_(ixdrf, nset, iret)
1877 call xdrfint_(ixdrf,mset(i), iret)
1880 call xdrfint_(ixdrf,i2set(i), iret)
1886 call xdrfint_(ixdrf,itmp, iret)
1887 i_index(i,j,il,il1)=itmp
1895 call xdrfint(ixdrf, nset, iret)
1897 call xdrfint(ixdrf,mset(i), iret)
1900 call xdrfint(ixdrf,i2set(i), iret)
1906 call xdrfint(ixdrf,itmp, iret)
1907 i_index(i,j,il,il1)=itmp
1914 call mpi_scatter(i2set,1,mpi_integer,&
1915 iset,1,mpi_integer,king,&
1920 if(me.eq.king) close(irest2)
1922 end subroutine read1restart
1923 !-----------------------------------------------------------------------------
1924 subroutine read1restart_old
1926 ! implicit real*8 (a-h,o-z)
1927 ! include 'DIMENSIONS'
1929 ! include 'COMMON.MD'
1930 ! include 'COMMON.IOUNITS'
1931 ! include 'COMMON.REMD'
1932 ! include 'COMMON.SETUP'
1933 ! include 'COMMON.CHAIN'
1934 ! include 'COMMON.SBRIDGE'
1935 ! include 'COMMON.INTERACT'
1936 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1937 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1938 !el common /przechowalnia/ d_restart1
1940 integer :: i,j,il,ierr
1943 open(irest2,file=mremd_rst_name,status='unknown')
1944 read (irest2,*) (i2rep(i),i=0,nodes-1)
1945 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1947 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1948 read (irest2,*) ndowna(0,il),&
1949 (ndowna(i,il),i=1,ndowna(0,il))
1952 read(irest2,*) (t_restart1(j,il),j=1,4)
1955 call mpi_scatter(t_restart1,5,mpi_real,&
1956 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1960 t_bath=t5_restart1(4)
1965 read(irest2,'(3e15.5)') &
1966 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
1970 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1971 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1981 read(irest2,'(3e15.5)') &
1982 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
1986 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1987 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1993 if(me.eq.king) close(irest2)
1995 end subroutine read1restart_old
1996 !----------------------------------------------------------------
1997 subroutine alloc_MREMD_arrays
1999 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2000 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1)) !(ntyp1))
2002 ! common /remdcommon/ in io: read_REMDpar
2003 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2004 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2005 ! common /remdrestart/
2006 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2008 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2009 allocate(ifirst(0:nodes)) !(maxprocs)
2010 allocate(nupa(0:nodes,0:2*nodes))
2011 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2012 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2013 allocate(iset_restart1(nodes)) !(maxprocs)
2014 ! common /traj1cache/
2015 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2016 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2017 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2018 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2019 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2020 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2021 allocate(utheta_cache(nfrag_back,max_cache_traj))
2022 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2023 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2024 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2027 end subroutine alloc_MREMD_arrays
2028 !-----------------------------------------------------------------------------
2029 !-----------------------------------------------------------------------------