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*nodes))
197 if(.not.allocated(d_restart2)) allocate(d_restart2(3,nres2*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,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,mpi_real,&
1358 d_restart1,3*2*nres,mpi_real,king,&
1367 call mpi_gather(r_d,3*2*nres,mpi_real,&
1368 d_restart2,3*2*nres,mpi_real,king,&
1373 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1375 call xdrfint_(ixdrf, i2rep(i), iret)
1378 call xdrfint_(ixdrf, ifirst(i), iret)
1382 call xdrfint_(ixdrf, nupa(i,il), iret)
1386 call xdrfint_(ixdrf, ndowna(i,il), iret)
1392 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1399 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1406 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1412 call xdrfint_(ixdrf, nset, iret)
1414 call xdrfint_(ixdrf,mset(i), iret)
1417 call xdrfint_(ixdrf,i2set(i), iret)
1423 itmp=i_index(i,j,il,il1)
1424 call xdrfint_(ixdrf,itmp, iret)
1431 call xdrfclose_(ixdrf, iret)
1433 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1435 call xdrfint(ixdrf, i2rep(i), iret)
1438 call xdrfint(ixdrf, ifirst(i), iret)
1442 call xdrfint(ixdrf, nupa(i,il), iret)
1446 call xdrfint(ixdrf, ndowna(i,il), iret)
1452 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1459 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1466 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1473 call xdrfint(ixdrf, nset, iret)
1475 call xdrfint(ixdrf,mset(i), iret)
1478 call xdrfint(ixdrf,i2set(i), iret)
1484 itmp=i_index(i,j,il,il1)
1485 call xdrfint(ixdrf,itmp, iret)
1492 call xdrfclose(ixdrf, iret)
1496 end subroutine write1rst
1497 !-----------------------------------------------------------------------------
1498 subroutine write1traj
1500 ! implicit real*8 (a-h,o-z)
1501 ! include 'DIMENSIONS'
1503 ! include 'COMMON.MD'
1504 ! include 'COMMON.IOUNITS'
1505 ! include 'COMMON.REMD'
1506 ! include 'COMMON.SETUP'
1507 ! include 'COMMON.CHAIN'
1508 ! include 'COMMON.SBRIDGE'
1509 ! include 'COMMON.INTERACT'
1511 real(kind=4) :: t5_restart1(5)
1512 integer :: iret,itmp
1513 real(kind=4) :: xcoord(3,2*nres+2),prec
1514 real(kind=4) :: r_qfrag(50),r_qpair(100)
1515 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1516 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1517 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1518 p_uscdiff(100*Nprocs) !(100*maxprocs)
1519 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1520 real(kind=4) :: r_c(3,2*nres+2)
1521 !el common /przechowalnia/ p_c
1523 integer :: i,j,il,ierr,ii,ixdrf
1525 call mpi_bcast(ii_write,1,mpi_integer,&
1529 print *,'traj1file',me,ii_write,ntwx_cache
1533 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1535 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1538 ! write (iout,*) "before gather write1traj: from node",ii
1540 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1542 t5_restart1(1)=totT_cache(ii)
1543 t5_restart1(2)=EK_cache(ii)
1544 t5_restart1(3)=potE_cache(ii)
1545 t5_restart1(4)=t_bath_cache(ii)
1546 t5_restart1(5)=Uconst_cache(ii)
1547 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1549 call mpi_gather(t5_restart1,5,mpi_real,&
1550 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1552 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1555 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1556 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1559 r_qfrag(i)=qfrag_cache(i,ii)
1562 r_qpair(i)=qpair_cache(i,ii)
1565 r_utheta(i)=utheta_cache(i,ii)
1566 r_ugamma(i)=ugamma_cache(i,ii)
1567 r_uscdiff(i)=uscdiff_cache(i,ii)
1570 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1571 p_qfrag,nfrag,mpi_real,king,&
1573 call mpi_gather(r_qpair,npair,mpi_real,&
1574 p_qpair,npair,mpi_real,king,&
1576 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1577 p_utheta,nfrag_back,mpi_real,king,&
1579 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1580 p_ugamma,nfrag_back,mpi_real,king,&
1582 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1583 p_uscdiff,nfrag_back,mpi_real,king,&
1587 write (iout,*) "p_qfrag"
1589 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1591 write (iout,*) "p_qpair"
1593 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1599 r_c(j,i)=c_cache(j,i,ii)
1603 call mpi_gather(r_c,3*2*nres,mpi_real,&
1604 p_c,3*2*nres,mpi_real,king,&
1610 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1611 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1612 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1613 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1614 call xdrfint_(ixdrf, nss, iret)
1617 call xdrfint(ixdrf, idssb(j)+nres, iret)
1618 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1620 call xdrfint_(ixdrf, ihpb(j), iret)
1621 call xdrfint_(ixdrf, jhpb(j), iret)
1624 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1625 call xdrfint_(ixdrf, iset_restart1(il), iret)
1627 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1630 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1633 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1634 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1635 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1640 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1645 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1649 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1653 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1654 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1655 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1656 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1657 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1658 call xdrfint(ixdrf, nss, iret)
1661 call xdrfint(ixdrf, idssb(j)+nres, iret)
1662 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1664 call xdrfint(ixdrf, ihpb(j), iret)
1665 call xdrfint(ixdrf, jhpb(j), iret)
1668 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1669 call xdrfint(ixdrf, iset_restart1(il), iret)
1671 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1674 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1677 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1678 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1679 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1684 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1689 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1693 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1699 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1701 if(me.eq.king) call xdrfclose(ixdrf, iret)
1703 do i=1,ntwx_cache-ii_write
1705 totT_cache(i)=totT_cache(ii_write+i)
1706 EK_cache(i)=EK_cache(ii_write+i)
1707 potE_cache(i)=potE_cache(ii_write+i)
1708 t_bath_cache(i)=t_bath_cache(ii_write+i)
1709 Uconst_cache(i)=Uconst_cache(ii_write+i)
1710 iset_cache(i)=iset_cache(ii_write+i)
1713 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1716 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1719 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1720 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1721 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1726 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1730 ntwx_cache=ntwx_cache-ii_write
1732 end subroutine write1traj
1733 !-----------------------------------------------------------------------------
1734 subroutine read1restart(i_index)
1736 ! implicit real*8 (a-h,o-z)
1737 ! include 'DIMENSIONS'
1739 ! include 'COMMON.MD'
1740 ! include 'COMMON.IOUNITS'
1741 ! include 'COMMON.REMD'
1742 ! include 'COMMON.SETUP'
1743 ! include 'COMMON.CHAIN'
1744 ! include 'COMMON.SBRIDGE'
1745 ! include 'COMMON.INTERACT'
1746 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1747 real(kind=4) :: r_d(3,2*nres),t5_restart1(5)
1748 integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1749 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1750 !el common /przechowalnia/ d_restart1
1751 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1753 write (*,*) "Processor",me," called read1restart"
1756 open(irest2,file=mremd_rst_name,status='unknown')
1757 read(irest2,*,err=334) i
1758 write(iout,*) "Reading old rst in ASCI format"
1760 call read1restart_old
1764 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1767 call xdrfint_(ixdrf, i2rep(i), iret)
1770 call xdrfint_(ixdrf, ifirst(i), iret)
1773 call xdrfint_(ixdrf, nupa(0,il), iret)
1775 call xdrfint_(ixdrf, nupa(i,il), iret)
1778 call xdrfint_(ixdrf, ndowna(0,il), iret)
1780 call xdrfint_(ixdrf, ndowna(i,il), iret)
1785 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1789 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1792 call xdrfint(ixdrf, i2rep(i), iret)
1795 call xdrfint(ixdrf, ifirst(i), iret)
1798 call xdrfint(ixdrf, nupa(0,il), iret)
1800 call xdrfint(ixdrf, nupa(i,il), iret)
1803 call xdrfint(ixdrf, ndowna(0,il), iret)
1805 call xdrfint(ixdrf, ndowna(i,il), iret)
1810 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1815 call mpi_scatter(t_restart1,5,mpi_real,&
1816 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1820 t_bath=t5_restart1(4)
1825 ! read(irest2,'(3e15.5)')
1826 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1829 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1831 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1837 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1838 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1848 ! read(irest2,'(3e15.5)')
1849 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1852 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1854 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1860 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1861 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1872 call xdrfint_(ixdrf, nset, iret)
1874 call xdrfint_(ixdrf,mset(i), iret)
1877 call xdrfint_(ixdrf,i2set(i), iret)
1883 call xdrfint_(ixdrf,itmp, iret)
1884 i_index(i,j,il,il1)=itmp
1892 call xdrfint(ixdrf, nset, iret)
1894 call xdrfint(ixdrf,mset(i), iret)
1897 call xdrfint(ixdrf,i2set(i), iret)
1903 call xdrfint(ixdrf,itmp, iret)
1904 i_index(i,j,il,il1)=itmp
1911 call mpi_scatter(i2set,1,mpi_integer,&
1912 iset,1,mpi_integer,king,&
1917 if(me.eq.king) close(irest2)
1919 end subroutine read1restart
1920 !-----------------------------------------------------------------------------
1921 subroutine read1restart_old
1923 ! implicit real*8 (a-h,o-z)
1924 ! include 'DIMENSIONS'
1926 ! include 'COMMON.MD'
1927 ! include 'COMMON.IOUNITS'
1928 ! include 'COMMON.REMD'
1929 ! include 'COMMON.SETUP'
1930 ! include 'COMMON.CHAIN'
1931 ! include 'COMMON.SBRIDGE'
1932 ! include 'COMMON.INTERACT'
1933 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1934 real(kind=4) :: r_d(3,2*nres),t5_restart1(5)
1935 !el common /przechowalnia/ d_restart1
1937 integer :: i,j,il,ierr
1940 open(irest2,file=mremd_rst_name,status='unknown')
1941 read (irest2,*) (i2rep(i),i=0,nodes-1)
1942 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1944 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1945 read (irest2,*) ndowna(0,il),&
1946 (ndowna(i,il),i=1,ndowna(0,il))
1949 read(irest2,*) (t_restart1(j,il),j=1,4)
1952 call mpi_scatter(t_restart1,5,mpi_real,&
1953 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1957 t_bath=t5_restart1(4)
1962 read(irest2,'(3e15.5)') &
1963 (d_restart1(j,i+2*nres*il),j=1,3)
1967 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1968 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1978 read(irest2,'(3e15.5)') &
1979 (d_restart1(j,i+2*nres*il),j=1,3)
1983 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1984 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1990 if(me.eq.king) close(irest2)
1992 end subroutine read1restart_old
1993 !----------------------------------------------------------------
1994 subroutine alloc_MREMD_arrays
1996 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
1997 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1)) !(ntyp1))
1999 ! common /remdcommon/ in io: read_REMDpar
2000 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2001 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2002 ! common /remdrestart/
2003 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2005 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2006 allocate(ifirst(0:nodes)) !(maxprocs)
2007 allocate(nupa(0:nodes,0:2*nodes))
2008 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2009 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2010 allocate(iset_restart1(nodes)) !(maxprocs)
2011 ! common /traj1cache/
2012 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2013 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2014 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2015 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2016 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2017 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2018 allocate(utheta_cache(nfrag_back,max_cache_traj))
2019 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2020 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2021 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2024 end subroutine alloc_MREMD_arrays
2025 !-----------------------------------------------------------------------------
2026 !-----------------------------------------------------------------------------