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)
86 real(kind=8) :: remd_t_bath(Nprocs) !(maxprocs)
87 integer :: iremd_iset(Nprocs) !(maxprocs)
88 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
89 ! (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
90 real(kind=8) :: remd_ene(0:n_ene+4,Nprocs) !(0:n_ene+4,maxprocs)
91 integer :: iremd_acc(Nprocs),iremd_tot(Nprocs) !(maxprocs)
92 integer :: iremd_acc_usa(Nprocs),iremd_tot_usa(Nprocs) !(maxprocs)
93 integer :: rstcount !el ilen,
95 character(len=50) :: tytul
98 !old integer nup(0:maxprocs),ndown(0:maxprocs)
99 integer :: rep2i(0:Nprocs),ireqi(Nprocs) !(maxprocs)
100 integer :: icache_all(Nprocs) !(maxprocs)
101 integer :: status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,Nprocs)! (MPI_STATUS_SIZE,maxprocs)
102 logical :: synflag, end_of_run, file_exist = .false.!, ovrtim
104 real(kind=8) :: delta,time00,time01,time001,time02,time03,time04,&
105 time05,time06,time07,time08,tt0,scalfac,ene_iex_iex,&
106 ene_i_i,ene_iex_i,ene_i_iex,xxx,tmp,econstr_temp_i,&
108 integer :: k,il,il1,i,j,nharp,ii,ierr,itime_master,irr,iex,&
109 i_set_temp,itmp,i_temp,i_mult,i_iset,i_mset,i_dir,i_temp1,&
110 i_mult1,i_iset1,i_mset1,ierror
111 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
112 !deb imin_itime_old=0
114 WRITE(iout,*) "JUST AFTER CALL"
115 ! if (.not.allocated(remd_ene)) allocate(remd_ene(0:n_ene+4,Nprocs))
122 if(me.eq.king.or..not.out1file) then
123 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
124 write (iout,*) "NREP=",nrep
128 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
129 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
131 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
133 print *,'MREMD',nodes
134 print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
135 write (iout,*) "Start MREMD: me",me," t_bath",t_bath
136 print *,"NSET=",nset, "MSET=", mset
140 do il1=1,max0(mset(il),1)
150 i_index(i,j,il,il1)=k
156 if(me.eq.king.or..not.out1file) then
157 write(iout,*) (i2rep(i),i=0,nodes-1)
158 write(iout,*) (i2set(i),i=0,nodes-1)
163 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
170 ! print *,'i2rep',me,i2rep(me)
171 ! print *,'rep2i',(rep2i(i),i=0,nrep)
173 !old if (i2rep(me).eq.nrep) then
176 !old nup(0)=remd_m(i2rep(me)+1)
177 !old k=rep2i(int(i2rep(me)))+1
184 !d print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
186 !old if (i2rep(me).eq.1) then
189 !old ndown(0)=remd_m(i2rep(me)-1)
190 !old k=rep2i(i2rep(me)-2)+1
197 !d print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
199 !el common /przechowalnia/
200 if(.not.allocated(d_restart1)) allocate(d_restart1(3,0:(nres2+1)*nodes))
201 if(.not.allocated(d_restart2)) allocate(d_restart2(3,0:(nres2+1)*nodes))
202 if(.not.allocated(p_c)) allocate(p_c(3,(nres2+2)*nodes))
205 write (*,*) "Processor",me," rest",rest,&
206 "restart1fie",restart1file
207 if(rest.and.restart1file) then
209 inquire(file=mremd_rst_name,exist=file_exist)
210 !d write (*,*) me," Before broadcast: file_exist",file_exist
211 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
213 !d write (*,*) me," After broadcast: file_exist",file_exist
215 if(me.eq.king.or..not.out1file) &
216 write (iout,*) 'Master is reading restart1file'
217 call read1restart(i_index)
219 if(me.eq.king.or..not.out1file) &
220 write (iout,*) 'WARNING : no restart1file'
223 if(me.eq.king.or..not.out1file) then
224 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
225 write(iout,*) "i_index"
230 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
239 if (rest.and..not.restart1file) &
240 inquire(file=mremd_rst_name,exist=file_exist)
241 if(.not.file_exist.and.rest.and..not.restart1file) &
242 write(iout,*) 'WARNING : no restart file',mremd_rst_name
243 IF (rest.and.file_exist.and..not.restart1file) THEN
244 write (iout,*) 'Master is reading restart file',&
246 open(irest2,file=mremd_rst_name,status='unknown')
248 read (irest2,*) (i2rep(i),i=0,nodes-1)
250 read (irest2,*) (ifirst(i),i=1,remd_m(1))
253 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
255 read (irest2,*) ndowna(0,il),&
256 (ndowna(i,il),i=1,ndowna(0,il))
262 read (irest2,*) (mset(i),i=1,nset)
264 read (irest2,*) (i2set(i),i=0,nodes-1)
269 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
274 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
275 write(iout,*) "i_index"
280 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
289 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
290 write (iout,'(a6,1000i5)') "ifirst",&
291 (ifirst(i),i=1,remd_m(1))
293 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
294 (nupa(i,il),i=1,nupa(0,il))
295 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
296 (ndowna(i,il),i=1,ndowna(0,il))
298 ELSE IF (.not.(rest.and.file_exist)) THEN
304 if (i2rep(il-1).eq.nrep) then
307 nupa(0,il)=remd_m(i2rep(il-1)+1)
308 k=rep2i(int(i2rep(il-1)))+1
314 if (i2rep(il-1).eq.1) then
317 ndowna(0,il)=remd_m(i2rep(il-1)-1)
318 k=rep2i(i2rep(il-1)-2)+1
326 write (iout,'(a6,100i4)') "ifirst",&
327 (ifirst(i),i=1,remd_m(1))
329 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
330 (nupa(i,il),i=1,nupa(0,il))
331 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
332 (ndowna(i,il),i=1,ndowna(0,il))
338 ! t_bath=retmin+(retmax-retmin)*me/(nodes-1)
339 if(.not.(rest.and.file_exist.and.restart1file)) then
340 if (me .eq. king) then
343 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
345 !d print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
346 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
351 if(me.eq.king.or..not.out1file) &
352 write(iout,*) me,"iset=",iset,"t_bath=",t_bath
356 stdfp(i)=dsqrt(2*Rb*t_bath/d_time)
360 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
364 ! print *,'irep',me,t_bath
366 if (me.eq.king .or. .not. out1file) &
367 write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
368 call rescale_weights(t_bath)
372 !------copy MD--------------
373 ! The driver for molecular dynamics subroutines
374 !------------------------------------------------
380 if(me.eq.king.or..not.out1file) &
381 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
387 ! Determine the inverse of the inertia matrix.
388 call setup_MD_matrices
392 if (me.eq.king .or. .not. out1file) &
393 write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
395 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
397 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
400 call rescale_weights(t_bath)
404 t_MDsetup = MPI_Wtime()-tt0
406 t_MDsetup = tcpu()-tt0
409 ! Entering the MD loop
415 if (lang.eq.2 .or. lang.eq.3) then
419 call sd_verlet_p_setup
421 call sd_verlet_ciccotti_setup
425 pfric0_mat(i,j,0)=pfric_mat(i,j)
426 afric0_mat(i,j,0)=afric_mat(i,j)
427 vfric0_mat(i,j,0)=vfric_mat(i,j)
428 prand0_mat(i,j,0)=prand_mat(i,j)
429 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
430 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
435 flag_stoch(i)=.false.
439 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
441 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
445 else if (lang.eq.1 .or. lang.eq.4) then
449 if (me.eq.king .or. .not. out1file) &
450 write(iout,*) 'Setup time',time00-walltime
453 t_langsetup=MPI_Wtime()-tt0
456 t_langsetup=tcpu()-tt0
462 do while(.not.end_of_run)
464 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
465 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
467 if (lang.gt.0 .and. surfarea .and. &
468 mod(itime,reset_fricmat).eq.0) then
469 if (lang.eq.2 .or. lang.eq.3) then
473 call sd_verlet_p_setup
475 call sd_verlet_ciccotti_setup
479 pfric0_mat(i,j,0)=pfric_mat(i,j)
480 afric0_mat(i,j,0)=afric_mat(i,j)
481 vfric0_mat(i,j,0)=vfric_mat(i,j)
482 prand0_mat(i,j,0)=prand_mat(i,j)
483 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
484 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
489 flag_stoch(i)=.false.
492 else if (lang.eq.1 .or. lang.eq.4) then
495 write (iout,'(a,i10)') &
496 "Friction matrix reset based on surface area, itime",itime
498 if (reset_vel .and. tbf .and. lang.eq.0 &
499 .and. mod(itime,count_reset_vel).eq.0) then
501 if (me.eq.king .or. .not. out1file) &
502 write(iout,'(a,f20.2)') &
503 "Velocities reset to random values, time",totT
506 d_t_old(j,i)=d_t(j,i)
510 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
514 d_t(j,0)=d_t(j,0)-vcm(j)
517 kinetic_T=2.0d0/(dimen3*Rb)*EK
518 scalfac=dsqrt(T_bath/kinetic_T)
519 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
522 d_t_old(j,i)=scalfac*d_t(j,i)
528 ! Time-reversible RESPA algorithm
529 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
530 call RESPA_step(itime)
532 ! Variable time step algorithm.
533 call velverlet_step(itime)
537 call brown_step(itime)
539 print *,"Brown dynamics not here!"
541 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
547 if (mod(itime,ntwe).eq.0) then
549 call enerprint(potEcomp)
552 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
553 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
555 call hairpin(.true.,nharp,iharp)
556 call secondary2(.true.)
557 call pdbout(potE,tytul,ipdb)
562 if (mod(itime,ntwx).eq.0.and.traj1file) then
563 if(ntwx_cache.lt.max_cache_traj_use) then
564 ntwx_cache=ntwx_cache+1
566 if (max_cache_traj_use.ne.1) &
567 print *,itime,"processor ",me," over cache ",ntwx_cache
570 totT_cache(i)=totT_cache(i+1)
571 EK_cache(i)=EK_cache(i+1)
572 potE_cache(i)=potE_cache(i+1)
573 t_bath_cache(i)=t_bath_cache(i+1)
574 Uconst_cache(i)=Uconst_cache(i+1)
575 iset_cache(i)=iset_cache(i+1)
578 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
581 qpair_cache(ii,i)=qpair_cache(ii,i+1)
584 utheta_cache(ii,i)=utheta_cache(ii,i+1)
585 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
586 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
592 c_cache(j,ii,i)=c_cache(j,ii,i+1)
598 totT_cache(ntwx_cache)=totT
599 EK_cache(ntwx_cache)=EK
600 potE_cache(ntwx_cache)=potE
601 t_bath_cache(ntwx_cache)=t_bath
602 Uconst_cache(ntwx_cache)=Uconst
603 iset_cache(ntwx_cache)=iset
606 qfrag_cache(i,ntwx_cache)=qfrag(i)
609 qpair_cache(i,ntwx_cache)=qpair(i)
612 utheta_cache(i,ntwx_cache)=utheta(i)
613 ugamma_cache(i,ntwx_cache)=ugamma(i)
614 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
619 c_cache(j,i,ntwx_cache)=c(j,i)
624 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
625 .and..not.restart1file) then
628 open(irest1,file=mremd_rst_name,status='unknown')
629 write (irest1,*) "i2rep"
630 write (irest1,*) (i2rep(i),i=0,nodes-1)
631 write (irest1,*) "ifirst"
632 write (irest1,*) (ifirst(i),i=1,remd_m(1))
634 write (irest1,*) "nupa",il
635 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
636 write (irest1,*) "ndowna",il
637 write (irest1,*) ndowna(0,il),&
638 (ndowna(i,il),i=1,ndowna(0,il))
641 write (irest1,*) "nset"
642 write (irest1,*) nset
643 write (irest1,*) "mset"
644 write (irest1,*) (mset(i),i=1,nset)
645 write (irest1,*) "i2set"
646 write (irest1,*) (i2set(i),i=0,nodes-1)
647 write (irest1,*) "i_index"
651 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
659 open(irest2,file=rest2name,status='unknown')
660 write(irest2,*) totT,EK,potE,totE,t_bath
662 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
665 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
668 write (irest2,*) iset
675 ! forced synchronization
676 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
677 .and. .not. mremdsync) then
679 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
681 call mpi_recv(itime_master, 1, MPI_INTEGER,&
682 0,101,CG_COMM, status, ierr)
683 call mpi_barrier(CG_COMM, ierr)
684 !deb if (out1file.or.traj1file) then
685 !deb call mpi_gather(itime,1,mpi_integer,
686 !deb & icache_all,1,mpi_integer,king,
689 call mpi_gather(ntwx_cache,1,mpi_integer,&
690 icache_all,1,mpi_integer,king,&
693 write(iout,*) 'REMD synchro at3',itime_master,itime
694 if (itime_master.ge.n_timestep .or. ovrtim()) &
696 !time call flush(iout)
701 if ((mod(itime,nstex).eq.0.and.me.eq.king &
702 .or.end_of_run.and.me.eq.king ) &
703 .and. .not. mremdsync ) then
706 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
707 CG_COMM, ireqi(i), ierr)
708 !d write(iout,*) 'REMD synchro with',i
711 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
712 call mpi_barrier(CG_COMM, ierr)
714 write(iout,*) 'REMD synchro at2',itime,'time=',time01-time00
715 if (out1file.or.traj1file) then
716 !deb call mpi_gather(itime,1,mpi_integer,
717 !deb & itime_all,1,mpi_integer,king,
719 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
720 !deb & (itime_all(i),i=1,nodes)
722 !deb imin_itime=itime_all(1)
724 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
726 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
727 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
728 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
729 call mpi_gather(ntwx_cache,1,mpi_integer,&
730 icache_all,1,mpi_integer,king,&
732 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
733 ! & (icache_all(i),i=1,nodes)
734 ii_write=icache_all(1)
736 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
738 ! write(iout,*) "MIN ii_write=",ii_write
741 !time call flush(iout)
743 if(mremdsync .and. mod(itime,nstex).eq.0) then
745 if (me.eq.king .or. .not. out1file) &
746 write(iout,*) 'REMD synchro at1',itime,ntwx_cache,Nprocs,nodes
747 write(iout,*) icache_all
749 write(iout,*) "before mpi_gather ntwx_cache"
750 call mpi_gather(ntwx_cache,1,mpi_integer,&
751 icache_all(1),1,mpi_integer,king,& ! CONSULT WITH ADAM
753 write(iout,*) "after mpi_gather ntwx_cache"
756 write(iout,'(a19,8000i8)') ' ntwx_cache',&
757 (icache_all(i),i=1,nodes)
758 ii_write=icache_all(1)
760 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
762 write(iout,*) "MIN ii_write=",ii_write
768 ! Update the time safety limiy
769 if (time001-time00.gt.safety) then
770 safety=time001-time00+600
771 write (iout,*) "****** SAFETY increased to",safety," s"
773 if (ovrtim()) end_of_run=.true.
775 if(synflag.and..not.end_of_run) then
779 write(iout,*) 'REMD before',me,t_bath
781 ! call mpi_gather(t_bath,1,mpi_double_precision,
782 ! & remd_t_bath,1,mpi_double_precision,king,
784 potEcomp(n_ene+1)=t_bath
786 potEcomp(n_ene+2)=iset
787 if (iset.lt.nset) then
791 potEcomp(n_ene+3)=Uconst
798 potEcomp(n_ene+4)=Uconst
802 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
803 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
806 call mpi_gather(elow,1,mpi_double_precision,&
807 elowi,1,mpi_double_precision,king,&
809 call mpi_gather(ehigh,1,mpi_double_precision,&
810 ehighi,1,mpi_double_precision,king,&
815 if (me.eq.king .or. .not. out1file) then
816 write(iout,*) 'REMD gather times=',time03-time01 &
820 if (restart1file) call write1rst(i_index)
823 if (me.eq.king .or. .not. out1file) then
824 write(iout,*) 'REMD writing rst time=',time04-time03
827 if (traj1file) call write1traj
829 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
830 !deb & icache_all,1,mpi_integer,king,
832 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
833 !deb & (icache_all(i),i=1,nodes)
838 if (me.eq.king .or. .not. out1file) then
839 write(iout,*) 'REMD writing traj time=',time05-time04
846 remd_t_bath(i)=remd_ene(n_ene+1,i)
847 iremd_iset(i)=remd_ene(n_ene+2,i)
851 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
853 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
857 write(iout,*) 'REMD exchange temp,ene'
859 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
860 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
864 !-------------------------------------
866 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
869 write (iout,*) "remd_m(1)",remd_m(1)
871 i=ifirst(iran_num(1,remd_m(1)))
878 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
880 if(i.gt.0.and.nupa(0,i).gt.0) then
882 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
884 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
886 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
888 ! do while (iex.eq.i)
889 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
890 iex=nupa(iran_num(1,int(nupa(0,i))),i)
892 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
894 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
896 ! Swap temperatures between conformations i and iex with recalculating the free energies
897 ! following temperature changes.
898 ene_iex_iex=remd_ene(0,iex)
899 ene_i_i=remd_ene(0,i)
900 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
901 ! & " iex",iex," ene_iex_iex",ene_iex_iex
902 ! write (iout,*) "rescaling weights with temperature",
905 call rescale_weights(remd_t_bath(i))
907 ! write (iout,*) "0,iex",remd_t_bath(i)
908 ! call enerprint(remd_ene(0,iex))
910 call sum_energy(remd_ene(0,iex),.false.)
911 ene_iex_i=remd_ene(0,iex)
912 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
914 ! write (iout,*) "0,i",remd_t_bath(i)
915 ! call enerprint(remd_ene(0,i))
917 call sum_energy(remd_ene(0,i),.false.)
918 ! write (iout,*) "ene_i_i",remd_ene(0,i)
920 ! write (iout,*) "rescaling weights with temperature",
922 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
923 write (iout,*) "ERROR: inconsistent energies:",i,&
924 ene_i_i,remd_ene(0,i)
926 call rescale_weights(remd_t_bath(iex))
928 ! write (iout,*) "0,i",remd_t_bath(iex)
929 ! call enerprint(remd_ene(0,i))
931 call sum_energy(remd_ene(0,i),.false.)
932 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
934 ene_i_iex=remd_ene(0,i)
936 ! write (iout,*) "0,iex",remd_t_bath(iex)
937 ! call enerprint(remd_ene(0,iex))
939 call sum_energy(remd_ene(0,iex),.false.)
940 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
941 write (iout,*) "ERROR: inconsistent energies:",iex,&
942 ene_iex_iex,remd_ene(0,iex)
944 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
945 ! write (iout,*) "i",i," iex",iex
946 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
947 ! & " ene_i_iex",ene_i_iex,
948 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
950 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
951 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
953 ! write(iout,*) 'delta',delta
954 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
955 ! & (remd_ene(i)-remd_ene(iex))/Rb/
956 ! & (remd_t_bath(i)*remd_t_bath(iex))
958 if (delta .gt. 50.0d0) then
964 else if (delta.lt.-50.0d0) then
973 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
974 xxx=ran_number(0.0d0,1.0d0)
975 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
977 if (delta .gt. xxx) then
979 remd_t_bath(i)=remd_t_bath(iex)
981 remd_ene(0,i)=ene_i_iex
982 remd_ene(0,iex)=ene_iex_i
988 ehighi(i)=ehighi(iex)
995 nupa(k,i)=nupa(k,iex)
998 ndowna(k,i)=ndowna(k,iex)
1002 if (ifirst(il).eq.i) ifirst(il)=iex
1004 if (nupa(k,il).eq.i) then
1006 elseif (nupa(k,il).eq.iex) then
1011 if (ndowna(k,il).eq.i) then
1013 elseif (ndowna(k,il).eq.iex) then
1019 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1021 i2rep(i-1)=i2rep(iex-1)
1024 ! write(iout,*) 'exchange',i,iex
1025 ! write (iout,'(a8,100i4)') "@ ifirst",
1026 ! & (ifirst(k),k=1,remd_m(1))
1028 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1029 ! & (nupa(k,il),k=1,nupa(0,il))
1030 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1031 ! & (ndowna(k,il),k=1,ndowna(0,il))
1036 remd_ene(0,iex)=ene_iex_iex
1037 remd_ene(0,i)=ene_i_i
1043 !d write (iout,*) "exchange completed"
1047 !d write(iout,*) "########",ii
1049 i_temp=iran_num(1,nrep)
1050 i_mult=iran_num(1,remd_m(i_temp))
1051 i_iset=iran_num(1,nset)
1052 i_mset=iran_num(1,mset(i_iset))
1053 i=i_index(i_temp,i_mult,i_iset,i_mset)
1055 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1058 !d write(iout,*) "i_dir=",i_dir
1060 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1063 i_mult1=iran_num(1,remd_m(i_temp1))
1065 i_mset1=iran_num(1,mset(i_iset1))
1066 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1068 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1071 i_mult1=iran_num(1,remd_m(i_temp1))
1073 i_mset1=iran_num(1,mset(i_iset1))
1074 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1075 econstr_temp_i=remd_ene(20,i)
1076 econstr_temp_iex=remd_ene(20,iex)
1077 remd_ene(20,i)=remd_ene(n_ene+3,i)
1078 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1080 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1083 i_mult1=iran_num(1,remd_m(i_temp1))
1085 i_mset1=iran_num(1,mset(i_iset1))
1086 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1087 econstr_temp_i=remd_ene(20,i)
1088 econstr_temp_iex=remd_ene(20,iex)
1089 remd_ene(20,i)=remd_ene(n_ene+3,i)
1090 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1096 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1099 ! Swap temperatures between conformations i and iex with recalculating the free energies
1100 ! following temperature changes.
1101 ene_iex_iex=remd_ene(0,iex)
1102 ene_i_i=remd_ene(0,i)
1103 !o write (iout,*) "rescaling weights with temperature",
1105 call rescale_weights(remd_t_bath(i))
1107 call sum_energy(remd_ene(0,iex),.false.)
1108 ene_iex_i=remd_ene(0,iex)
1109 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1110 ! call sum_energy(remd_ene(0,i),.false.)
1111 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1112 ! write (iout,*) "rescaling weights with temperature",
1113 ! & remd_t_bath(iex)
1114 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1115 ! write (iout,*) "ERROR: inconsistent energies:",i,
1116 ! & ene_i_i,remd_ene(0,i)
1118 call rescale_weights(remd_t_bath(iex))
1119 call sum_energy(remd_ene(0,i),.false.)
1120 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1121 ene_i_iex=remd_ene(0,i)
1122 ! call sum_energy(remd_ene(0,iex),.false.)
1123 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1124 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1125 ! & ene_iex_iex,remd_ene(0,iex)
1127 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1128 ! write (iout,*) "i",i," iex",iex
1129 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1130 !d & " ene_i_iex",ene_i_iex,
1131 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1132 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1133 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1135 !d write(iout,*) 'delta',delta
1136 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1137 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1138 ! & (remd_t_bath(i)*remd_t_bath(iex))
1139 if (delta .gt. 50.0d0) then
1144 if (i_dir.eq.1.or.i_dir.eq.3) &
1145 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1146 if (i_dir.eq.2.or.i_dir.eq.3) &
1147 iremd_tot_usa(int(i2set(i-1)))= &
1148 iremd_tot_usa(int(i2set(i-1)))+1
1149 xxx=ran_number(0.0d0,1.0d0)
1150 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1151 if (delta .gt. xxx) then
1153 remd_t_bath(i)=remd_t_bath(iex)
1154 remd_t_bath(iex)=tmp
1157 iremd_iset(i)=iremd_iset(iex)
1158 iremd_iset(iex)=itmp
1160 remd_ene(0,i)=ene_i_iex
1161 remd_ene(0,iex)=ene_iex_i
1163 if (i_dir.eq.1.or.i_dir.eq.3) &
1164 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1167 i2rep(i-1)=i2rep(iex-1)
1170 if (i_dir.eq.2.or.i_dir.eq.3) &
1171 iremd_acc_usa(int(i2set(i-1)))= &
1172 iremd_acc_usa(int(i2set(i-1)))+1
1175 i2set(i-1)=i2set(iex-1)
1178 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1179 i_index(i_temp,i_mult,i_iset,i_mset)= &
1180 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1181 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1184 remd_ene(0,iex)=ene_iex_iex
1185 remd_ene(0,i)=ene_i_i
1186 remd_ene(20,iex)=econstr_temp_iex
1187 remd_ene(20,i)=econstr_temp_i
1191 !d do il1=1,mset(il)
1194 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1207 !-------------------------------------
1208 write (iout,*) "NREP",nrep
1210 if(iremd_tot(i).ne.0) &
1211 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1212 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1217 if(iremd_tot_usa(i).ne.0) &
1218 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1219 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1225 !d write (iout,'(a6,100i4)') "ifirst",
1226 !d & (ifirst(i),i=1,remd_m(1))
1228 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1229 !d & (nupa(i,il),i=1,nupa(0,il))
1230 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1231 !d & (ndowna(i,il),i=1,ndowna(0,il))
1236 !d write (iout,*) "Before scatter"
1239 if (me.eq.king) then
1240 write (iout,*) "t_bath before scatter",remd_t_bath
1244 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1245 t_bath,1,mpi_double_precision,king,&
1247 !d write (iout,*) "After scatter"
1250 call mpi_scatter(iremd_iset,1,mpi_integer,&
1251 iset,1,mpi_integer,king,&
1255 if (me.eq.king .or. .not. out1file) then
1256 write(iout,*) 'REMD scatter time=',time07-time06
1260 call mpi_scatter(elowi,1,mpi_double_precision,&
1261 elow,1,mpi_double_precision,king,&
1263 call mpi_scatter(ehighi,1,mpi_double_precision,&
1264 ehigh,1,mpi_double_precision,king,&
1267 call rescale_weights(t_bath)
1268 !o write (iout,*) "Processor",me,
1269 !o & " rescaling weights with temperature",t_bath
1271 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
1273 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
1277 !de write(iout,*) 'REMD after',me,t_bath
1279 if (me.eq.king .or. .not. out1file) then
1280 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1281 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1282 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1283 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1284 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1285 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1286 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1287 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1288 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1294 if (restart1file) then
1295 if (me.eq.king .or. .not. out1file) &
1296 write(iout,*) 'writing restart at the end of run'
1297 call write1rst(i_index)
1300 if (traj1file) call write1traj
1302 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1303 !deb & icache_all,1,mpi_integer,king,
1304 !deb & CG_COMM,ierr)
1305 !deb write(iout,'(a40,8000i8)')
1306 !deb & ' ntwx_cache after traj1file at the end',
1307 !deb & (icache_all(i),i=1,nodes)
1312 t_MD=MPI_Wtime()-tt0
1316 if (me.eq.king .or. .not. out1file) then
1317 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1319 'MD calculations setup:',t_MDsetup,&
1320 'Energy & gradient evaluation:',t_enegrad,&
1321 'Stochastic MD setup:',t_langsetup,&
1322 'Stochastic MD step setup:',t_sdsetup,&
1324 write (iout,'(/28(1h=),a25,27(1h=))') &
1325 ' End of MD calculation '
1327 !el common /przechowalnia/
1328 ! deallocate(d_restart1)
1329 ! deallocate(d_restart2)
1333 end subroutine MREMD
1334 !-----------------------------------------------------------------------------
1335 subroutine write1rst(i_index)
1338 ! implicit real*8 (a-h,o-z)
1339 ! include 'DIMENSIONS'
1341 ! include 'COMMON.MD'
1342 ! include 'COMMON.IOUNITS'
1343 ! include 'COMMON.REMD'
1344 ! include 'COMMON.SETUP'
1345 ! include 'COMMON.CHAIN'
1346 ! include 'COMMON.SBRIDGE'
1347 ! include 'COMMON.INTERACT'
1349 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1350 !el d_restart2(3,2*nres*maxprocs)
1351 real(kind=4) :: r_d(3,0:2*nres)
1352 real(kind=4) :: t5_restart1(5)
1353 integer :: iret,itmp
1354 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1355 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1357 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1358 !el common /przechowalnia/ d_restart1,d_restart2
1359 integer :: i,j,il,il1,ierr,ixdrf
1364 t5_restart1(4)=t_bath
1365 t5_restart1(5)=Uconst
1367 call mpi_gather(t5_restart1,5,mpi_real,&
1368 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1376 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1377 d_restart1,3*2*nres+3,mpi_real,king,&
1389 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1390 d_restart2,3*2*nres+3,mpi_real,king,&
1395 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1397 call xdrfint_(ixdrf, i2rep(i), iret)
1400 call xdrfint_(ixdrf, ifirst(i), iret)
1404 call xdrfint_(ixdrf, nupa(i,il), iret)
1408 call xdrfint_(ixdrf, ndowna(i,il), iret)
1414 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1421 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1428 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1434 call xdrfint_(ixdrf, nset, iret)
1436 call xdrfint_(ixdrf,mset(i), iret)
1439 call xdrfint_(ixdrf,i2set(i), iret)
1445 itmp=i_index(i,j,il,il1)
1446 call xdrfint_(ixdrf,itmp, iret)
1453 call xdrfclose_(ixdrf, iret)
1455 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1457 call xdrfint(ixdrf, i2rep(i), iret)
1460 call xdrfint(ixdrf, ifirst(i), iret)
1464 call xdrfint(ixdrf, nupa(i,il), iret)
1468 call xdrfint(ixdrf, ndowna(i,il), iret)
1474 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1481 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1488 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1495 call xdrfint(ixdrf, nset, iret)
1497 call xdrfint(ixdrf,mset(i), iret)
1500 call xdrfint(ixdrf,i2set(i), iret)
1506 itmp=i_index(i,j,il,il1)
1507 call xdrfint(ixdrf,itmp, iret)
1514 call xdrfclose(ixdrf, iret)
1518 end subroutine write1rst
1519 !-----------------------------------------------------------------------------
1520 subroutine write1traj
1522 ! implicit real*8 (a-h,o-z)
1523 ! include 'DIMENSIONS'
1525 ! include 'COMMON.MD'
1526 ! include 'COMMON.IOUNITS'
1527 ! include 'COMMON.REMD'
1528 ! include 'COMMON.SETUP'
1529 ! include 'COMMON.CHAIN'
1530 ! include 'COMMON.SBRIDGE'
1531 ! include 'COMMON.INTERACT'
1533 real(kind=4) :: t5_restart1(5)
1534 integer :: iret,itmp
1535 real(kind=4) :: xcoord(3,2*nres+2),prec
1536 real(kind=4) :: r_qfrag(50),r_qpair(100)
1537 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1538 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1539 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1540 p_uscdiff(100*Nprocs) !(100*maxprocs)
1541 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1542 real(kind=4) :: r_c(3,2*nres+2)
1543 !el common /przechowalnia/ p_c
1545 integer :: i,j,il,ierr,ii,ixdrf
1547 call mpi_bcast(ii_write,1,mpi_integer,&
1551 print *,'traj1file',me,ii_write,ntwx_cache
1555 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1557 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1560 ! write (iout,*) "before gather write1traj: from node",ii
1562 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1564 t5_restart1(1)=totT_cache(ii)
1565 t5_restart1(2)=EK_cache(ii)
1566 t5_restart1(3)=potE_cache(ii)
1567 t5_restart1(4)=t_bath_cache(ii)
1568 t5_restart1(5)=Uconst_cache(ii)
1569 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1571 call mpi_gather(t5_restart1,5,mpi_real,&
1572 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1574 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1577 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1578 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1581 r_qfrag(i)=qfrag_cache(i,ii)
1584 r_qpair(i)=qpair_cache(i,ii)
1587 r_utheta(i)=utheta_cache(i,ii)
1588 r_ugamma(i)=ugamma_cache(i,ii)
1589 r_uscdiff(i)=uscdiff_cache(i,ii)
1592 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1593 p_qfrag,nfrag,mpi_real,king,&
1595 call mpi_gather(r_qpair,npair,mpi_real,&
1596 p_qpair,npair,mpi_real,king,&
1598 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1599 p_utheta,nfrag_back,mpi_real,king,&
1601 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1602 p_ugamma,nfrag_back,mpi_real,king,&
1604 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1605 p_uscdiff,nfrag_back,mpi_real,king,&
1609 write (iout,*) "p_qfrag"
1611 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1613 write (iout,*) "p_qpair"
1615 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1621 r_c(j,i)=c_cache(j,i,ii)
1625 call mpi_gather(r_c,3*2*nres,mpi_real,&
1626 p_c,3*2*nres,mpi_real,king,&
1632 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1633 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1634 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1635 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1636 call xdrfint_(ixdrf, nss, iret)
1639 call xdrfint(ixdrf, idssb(j)+nres, iret)
1640 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1642 call xdrfint_(ixdrf, ihpb(j), iret)
1643 call xdrfint_(ixdrf, jhpb(j), iret)
1646 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1647 call xdrfint_(ixdrf, iset_restart1(il), iret)
1649 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1652 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1655 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1656 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1657 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1662 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1667 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1671 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1675 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1676 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1677 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1678 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1679 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1680 call xdrfint(ixdrf, nss, iret)
1683 call xdrfint(ixdrf, idssb(j)+nres, iret)
1684 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1686 call xdrfint(ixdrf, ihpb(j), iret)
1687 call xdrfint(ixdrf, jhpb(j), iret)
1690 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1691 call xdrfint(ixdrf, iset_restart1(il), iret)
1693 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1696 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1699 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1700 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1701 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1706 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1711 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1715 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1721 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1723 if(me.eq.king) call xdrfclose(ixdrf, iret)
1725 do i=1,ntwx_cache-ii_write
1727 totT_cache(i)=totT_cache(ii_write+i)
1728 EK_cache(i)=EK_cache(ii_write+i)
1729 potE_cache(i)=potE_cache(ii_write+i)
1730 t_bath_cache(i)=t_bath_cache(ii_write+i)
1731 Uconst_cache(i)=Uconst_cache(ii_write+i)
1732 iset_cache(i)=iset_cache(ii_write+i)
1735 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1738 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1741 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1742 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1743 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1748 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1752 ntwx_cache=ntwx_cache-ii_write
1754 end subroutine write1traj
1755 !-----------------------------------------------------------------------------
1756 subroutine read1restart(i_index)
1758 ! implicit real*8 (a-h,o-z)
1759 ! include 'DIMENSIONS'
1761 ! include 'COMMON.MD'
1762 ! include 'COMMON.IOUNITS'
1763 ! include 'COMMON.REMD'
1764 ! include 'COMMON.SETUP'
1765 ! include 'COMMON.CHAIN'
1766 ! include 'COMMON.SBRIDGE'
1767 ! include 'COMMON.INTERACT'
1768 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1769 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1770 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1771 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1773 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1774 !el common /przechowalnia/ d_restart1
1775 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1777 write (*,*) "Processor",me," called read1restart"
1780 open(irest2,file=mremd_rst_name,status='unknown')
1781 read(irest2,*,err=334) i
1782 write(iout,*) "Reading old rst in ASCI format"
1784 call read1restart_old
1788 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1791 call xdrfint_(ixdrf, i2rep(i), iret)
1794 call xdrfint_(ixdrf, ifirst(i), iret)
1797 call xdrfint_(ixdrf, nupa(0,il), iret)
1799 call xdrfint_(ixdrf, nupa(i,il), iret)
1802 call xdrfint_(ixdrf, ndowna(0,il), iret)
1804 call xdrfint_(ixdrf, ndowna(i,il), iret)
1809 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1813 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1816 call xdrfint(ixdrf, i2rep(i), iret)
1819 call xdrfint(ixdrf, ifirst(i), iret)
1822 call xdrfint(ixdrf, nupa(0,il), iret)
1824 call xdrfint(ixdrf, nupa(i,il), iret)
1827 call xdrfint(ixdrf, ndowna(0,il), iret)
1829 call xdrfint(ixdrf, ndowna(i,il), iret)
1834 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1839 call mpi_scatter(t_restart1,5,mpi_real,&
1840 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1844 t_bath=t5_restart1(4)
1849 ! read(irest2,'(3e15.5)')
1850 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1853 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1855 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1861 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1862 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1872 ! read(irest2,'(3e15.5)')
1873 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1876 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1878 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1884 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1885 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1896 call xdrfint_(ixdrf, nset, iret)
1898 call xdrfint_(ixdrf,mset(i), iret)
1901 call xdrfint_(ixdrf,i2set(i), iret)
1907 call xdrfint_(ixdrf,itmp, iret)
1908 i_index(i,j,il,il1)=itmp
1916 call xdrfint(ixdrf, nset, iret)
1918 call xdrfint(ixdrf,mset(i), iret)
1921 call xdrfint(ixdrf,i2set(i), iret)
1927 call xdrfint(ixdrf,itmp, iret)
1928 i_index(i,j,il,il1)=itmp
1935 call mpi_scatter(i2set,1,mpi_integer,&
1936 iset,1,mpi_integer,king,&
1941 if(me.eq.king) close(irest2)
1943 end subroutine read1restart
1944 !-----------------------------------------------------------------------------
1945 subroutine read1restart_old
1947 ! implicit real*8 (a-h,o-z)
1948 ! include 'DIMENSIONS'
1950 ! include 'COMMON.MD'
1951 ! include 'COMMON.IOUNITS'
1952 ! include 'COMMON.REMD'
1953 ! include 'COMMON.SETUP'
1954 ! include 'COMMON.CHAIN'
1955 ! include 'COMMON.SBRIDGE'
1956 ! include 'COMMON.INTERACT'
1957 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1958 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1959 !el common /przechowalnia/ d_restart1
1961 integer :: i,j,il,ierr
1964 open(irest2,file=mremd_rst_name,status='unknown')
1965 read (irest2,*) (i2rep(i),i=0,nodes-1)
1966 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1968 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1969 read (irest2,*) ndowna(0,il),&
1970 (ndowna(i,il),i=1,ndowna(0,il))
1973 read(irest2,*) (t_restart1(j,il),j=1,4)
1976 call mpi_scatter(t_restart1,5,mpi_real,&
1977 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1981 t_bath=t5_restart1(4)
1986 read(irest2,'(3e15.5)') &
1987 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
1991 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1992 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2002 read(irest2,'(3e15.5)') &
2003 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2007 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2008 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2014 if(me.eq.king) close(irest2)
2016 end subroutine read1restart_old
2017 !----------------------------------------------------------------
2018 subroutine alloc_MREMD_arrays
2020 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2021 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
2023 ! common /remdcommon/ in io: read_REMDpar
2024 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2025 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2026 ! common /remdrestart/
2027 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2029 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2030 allocate(ifirst(0:nodes)) !(maxprocs)
2031 allocate(nupa(0:nodes,0:2*nodes))
2032 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2033 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2034 allocate(iset_restart1(nodes)) !(maxprocs)
2035 ! common /traj1cache/
2036 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2037 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2038 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2039 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2040 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2041 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2042 allocate(utheta_cache(nfrag_back,max_cache_traj))
2043 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2044 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2045 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2048 end subroutine alloc_MREMD_arrays
2049 !-----------------------------------------------------------------------------
2050 !-----------------------------------------------------------------------------