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,max0(Nprocs/200,1),max0(Nprocs/200,1))
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 integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,max0(Nprocs/200,1),max0(Nprocs/200,1))
1347 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1348 !el common /przechowalnia/ d_restart1,d_restart2
1349 integer :: i,j,il,il1,ierr,ixdrf
1354 t5_restart1(4)=t_bath
1355 t5_restart1(5)=Uconst
1357 call mpi_gather(t5_restart1,5,mpi_real,&
1358 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1366 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1367 d_restart1,3*2*nres+3,mpi_real,king,&
1379 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1380 d_restart2,3*2*nres+3,mpi_real,king,&
1385 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1387 call xdrfint_(ixdrf, i2rep(i), iret)
1390 call xdrfint_(ixdrf, ifirst(i), iret)
1394 call xdrfint_(ixdrf, nupa(i,il), iret)
1398 call xdrfint_(ixdrf, ndowna(i,il), iret)
1404 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1411 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1418 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1424 call xdrfint_(ixdrf, nset, iret)
1426 call xdrfint_(ixdrf,mset(i), iret)
1429 call xdrfint_(ixdrf,i2set(i), iret)
1435 itmp=i_index(i,j,il,il1)
1436 call xdrfint_(ixdrf,itmp, iret)
1443 call xdrfclose_(ixdrf, iret)
1445 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1447 call xdrfint(ixdrf, i2rep(i), iret)
1450 call xdrfint(ixdrf, ifirst(i), iret)
1454 call xdrfint(ixdrf, nupa(i,il), iret)
1458 call xdrfint(ixdrf, ndowna(i,il), iret)
1464 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1471 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1478 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1485 call xdrfint(ixdrf, nset, iret)
1487 call xdrfint(ixdrf,mset(i), iret)
1490 call xdrfint(ixdrf,i2set(i), iret)
1496 itmp=i_index(i,j,il,il1)
1497 call xdrfint(ixdrf,itmp, iret)
1504 call xdrfclose(ixdrf, iret)
1508 end subroutine write1rst
1509 !-----------------------------------------------------------------------------
1510 subroutine write1traj
1512 ! implicit real*8 (a-h,o-z)
1513 ! include 'DIMENSIONS'
1515 ! include 'COMMON.MD'
1516 ! include 'COMMON.IOUNITS'
1517 ! include 'COMMON.REMD'
1518 ! include 'COMMON.SETUP'
1519 ! include 'COMMON.CHAIN'
1520 ! include 'COMMON.SBRIDGE'
1521 ! include 'COMMON.INTERACT'
1523 real(kind=4) :: t5_restart1(5)
1524 integer :: iret,itmp
1525 real(kind=4) :: xcoord(3,2*nres+2),prec
1526 real(kind=4) :: r_qfrag(50),r_qpair(100)
1527 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1528 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1529 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1530 p_uscdiff(100*Nprocs) !(100*maxprocs)
1531 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1532 real(kind=4) :: r_c(3,2*nres+2)
1533 !el common /przechowalnia/ p_c
1535 integer :: i,j,il,ierr,ii,ixdrf
1537 call mpi_bcast(ii_write,1,mpi_integer,&
1541 print *,'traj1file',me,ii_write,ntwx_cache
1545 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1547 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1550 ! write (iout,*) "before gather write1traj: from node",ii
1552 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1554 t5_restart1(1)=totT_cache(ii)
1555 t5_restart1(2)=EK_cache(ii)
1556 t5_restart1(3)=potE_cache(ii)
1557 t5_restart1(4)=t_bath_cache(ii)
1558 t5_restart1(5)=Uconst_cache(ii)
1559 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1561 call mpi_gather(t5_restart1,5,mpi_real,&
1562 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1564 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1567 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1568 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1571 r_qfrag(i)=qfrag_cache(i,ii)
1574 r_qpair(i)=qpair_cache(i,ii)
1577 r_utheta(i)=utheta_cache(i,ii)
1578 r_ugamma(i)=ugamma_cache(i,ii)
1579 r_uscdiff(i)=uscdiff_cache(i,ii)
1582 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1583 p_qfrag,nfrag,mpi_real,king,&
1585 call mpi_gather(r_qpair,npair,mpi_real,&
1586 p_qpair,npair,mpi_real,king,&
1588 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1589 p_utheta,nfrag_back,mpi_real,king,&
1591 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1592 p_ugamma,nfrag_back,mpi_real,king,&
1594 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1595 p_uscdiff,nfrag_back,mpi_real,king,&
1599 write (iout,*) "p_qfrag"
1601 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1603 write (iout,*) "p_qpair"
1605 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1611 r_c(j,i)=c_cache(j,i,ii)
1615 call mpi_gather(r_c,3*2*nres,mpi_real,&
1616 p_c,3*2*nres,mpi_real,king,&
1622 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1623 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1624 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1625 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1626 call xdrfint_(ixdrf, nss, iret)
1629 call xdrfint(ixdrf, idssb(j)+nres, iret)
1630 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1632 call xdrfint_(ixdrf, ihpb(j), iret)
1633 call xdrfint_(ixdrf, jhpb(j), iret)
1636 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1637 call xdrfint_(ixdrf, iset_restart1(il), iret)
1639 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1642 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1645 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1646 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1647 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1652 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1657 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1661 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1665 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1666 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1667 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1668 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1669 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1670 call xdrfint(ixdrf, nss, iret)
1673 call xdrfint(ixdrf, idssb(j)+nres, iret)
1674 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1676 call xdrfint(ixdrf, ihpb(j), iret)
1677 call xdrfint(ixdrf, jhpb(j), iret)
1680 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1681 call xdrfint(ixdrf, iset_restart1(il), iret)
1683 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1686 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1689 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1690 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1691 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1696 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1701 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1705 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1711 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1713 if(me.eq.king) call xdrfclose(ixdrf, iret)
1715 do i=1,ntwx_cache-ii_write
1717 totT_cache(i)=totT_cache(ii_write+i)
1718 EK_cache(i)=EK_cache(ii_write+i)
1719 potE_cache(i)=potE_cache(ii_write+i)
1720 t_bath_cache(i)=t_bath_cache(ii_write+i)
1721 Uconst_cache(i)=Uconst_cache(ii_write+i)
1722 iset_cache(i)=iset_cache(ii_write+i)
1725 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1728 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1731 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1732 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1733 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1738 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1742 ntwx_cache=ntwx_cache-ii_write
1744 end subroutine write1traj
1745 !-----------------------------------------------------------------------------
1746 subroutine read1restart(i_index)
1748 ! implicit real*8 (a-h,o-z)
1749 ! include 'DIMENSIONS'
1751 ! include 'COMMON.MD'
1752 ! include 'COMMON.IOUNITS'
1753 ! include 'COMMON.REMD'
1754 ! include 'COMMON.SETUP'
1755 ! include 'COMMON.CHAIN'
1756 ! include 'COMMON.SBRIDGE'
1757 ! include 'COMMON.INTERACT'
1758 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1759 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1760 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1761 integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,max0(Nprocs/200,1),max0(Nprocs/200,1))
1763 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1764 !el common /przechowalnia/ d_restart1
1765 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1767 write (*,*) "Processor",me," called read1restart"
1770 open(irest2,file=mremd_rst_name,status='unknown')
1771 read(irest2,*,err=334) i
1772 write(iout,*) "Reading old rst in ASCI format"
1774 call read1restart_old
1778 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1781 call xdrfint_(ixdrf, i2rep(i), iret)
1784 call xdrfint_(ixdrf, ifirst(i), iret)
1787 call xdrfint_(ixdrf, nupa(0,il), iret)
1789 call xdrfint_(ixdrf, nupa(i,il), iret)
1792 call xdrfint_(ixdrf, ndowna(0,il), iret)
1794 call xdrfint_(ixdrf, ndowna(i,il), iret)
1799 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1803 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1806 call xdrfint(ixdrf, i2rep(i), iret)
1809 call xdrfint(ixdrf, ifirst(i), iret)
1812 call xdrfint(ixdrf, nupa(0,il), iret)
1814 call xdrfint(ixdrf, nupa(i,il), iret)
1817 call xdrfint(ixdrf, ndowna(0,il), iret)
1819 call xdrfint(ixdrf, ndowna(i,il), iret)
1824 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1829 call mpi_scatter(t_restart1,5,mpi_real,&
1830 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1834 t_bath=t5_restart1(4)
1839 ! read(irest2,'(3e15.5)')
1840 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1843 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1845 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1851 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1852 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1862 ! read(irest2,'(3e15.5)')
1863 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1866 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1868 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1874 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1875 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1886 call xdrfint_(ixdrf, nset, iret)
1888 call xdrfint_(ixdrf,mset(i), iret)
1891 call xdrfint_(ixdrf,i2set(i), iret)
1897 call xdrfint_(ixdrf,itmp, iret)
1898 i_index(i,j,il,il1)=itmp
1906 call xdrfint(ixdrf, nset, iret)
1908 call xdrfint(ixdrf,mset(i), iret)
1911 call xdrfint(ixdrf,i2set(i), iret)
1917 call xdrfint(ixdrf,itmp, iret)
1918 i_index(i,j,il,il1)=itmp
1925 call mpi_scatter(i2set,1,mpi_integer,&
1926 iset,1,mpi_integer,king,&
1931 if(me.eq.king) close(irest2)
1933 end subroutine read1restart
1934 !-----------------------------------------------------------------------------
1935 subroutine read1restart_old
1937 ! implicit real*8 (a-h,o-z)
1938 ! include 'DIMENSIONS'
1940 ! include 'COMMON.MD'
1941 ! include 'COMMON.IOUNITS'
1942 ! include 'COMMON.REMD'
1943 ! include 'COMMON.SETUP'
1944 ! include 'COMMON.CHAIN'
1945 ! include 'COMMON.SBRIDGE'
1946 ! include 'COMMON.INTERACT'
1947 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1948 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1949 !el common /przechowalnia/ d_restart1
1951 integer :: i,j,il,ierr
1954 open(irest2,file=mremd_rst_name,status='unknown')
1955 read (irest2,*) (i2rep(i),i=0,nodes-1)
1956 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1958 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1959 read (irest2,*) ndowna(0,il),&
1960 (ndowna(i,il),i=1,ndowna(0,il))
1963 read(irest2,*) (t_restart1(j,il),j=1,4)
1966 call mpi_scatter(t_restart1,5,mpi_real,&
1967 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1971 t_bath=t5_restart1(4)
1976 read(irest2,'(3e15.5)') &
1977 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
1981 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1982 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1992 read(irest2,'(3e15.5)') &
1993 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
1997 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1998 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2004 if(me.eq.king) close(irest2)
2006 end subroutine read1restart_old
2007 !----------------------------------------------------------------
2008 subroutine alloc_MREMD_arrays
2010 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2011 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
2013 ! common /remdcommon/ in io: read_REMDpar
2014 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2015 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2016 ! common /remdrestart/
2017 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2019 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2020 allocate(ifirst(0:nodes)) !(maxprocs)
2021 allocate(nupa(0:nodes,0:2*nodes))
2022 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2023 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2024 allocate(iset_restart1(nodes)) !(maxprocs)
2025 ! common /traj1cache/
2026 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2027 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2028 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2029 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2030 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2031 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2032 allocate(utheta_cache(nfrag_back,max_cache_traj))
2033 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2034 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2035 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2038 end subroutine alloc_MREMD_arrays
2039 !-----------------------------------------------------------------------------
2040 !-----------------------------------------------------------------------------