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,Nprocs,Nprocs,Nprocs)
88 ! (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
89 real(kind=8),dimension(:,:), allocatable :: 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
113 WRITE(iout,*) "JUST AFTER CALL"
114 if (.not.allocated(remd_ene)) allocate(remd_ene(0:n_ene+30,Nprocs))
121 if(me.eq.king.or..not.out1file) then
122 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
123 write (iout,*) "NREP=",nrep
127 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
128 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
130 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
132 print *,'MREMD',nodes
133 print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
134 write (iout,*) "Start MREMD: me",me," t_bath",t_bath
135 print *,"NSET=",nset, "MSET=", mset
139 do il1=1,max0(mset(il),1)
149 i_index(i,j,il,il1)=k
155 if(me.eq.king.or..not.out1file) then
156 write(iout,*) (i2rep(i),i=0,nodes-1)
157 write(iout,*) (i2set(i),i=0,nodes-1)
162 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
169 ! print *,'i2rep',me,i2rep(me)
170 ! print *,'rep2i',(rep2i(i),i=0,nrep)
172 !old if (i2rep(me).eq.nrep) then
175 !old nup(0)=remd_m(i2rep(me)+1)
176 !old k=rep2i(int(i2rep(me)))+1
183 !d print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
185 !old if (i2rep(me).eq.1) then
188 !old ndown(0)=remd_m(i2rep(me)-1)
189 !old k=rep2i(i2rep(me)-2)+1
196 !d print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
198 !el common /przechowalnia/
199 if(.not.allocated(d_restart1)) allocate(d_restart1(3,0:(nres2+1)*nodes))
200 if(.not.allocated(d_restart2)) allocate(d_restart2(3,0:(nres2+1)*nodes))
201 if(.not.allocated(p_c)) allocate(p_c(3,(nres2+2)*nodes))
204 write (*,*) "Processor",me," rest",rest,&
205 "restart1fie",restart1file
206 if(rest.and.restart1file) then
208 inquire(file=mremd_rst_name,exist=file_exist)
209 !d write (*,*) me," Before broadcast: file_exist",file_exist
210 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
212 !d write (*,*) me," After broadcast: file_exist",file_exist
214 if(me.eq.king.or..not.out1file) &
215 write (iout,*) 'Master is reading restart1file'
216 call read1restart(i_index)
218 if(me.eq.king.or..not.out1file) &
219 write (iout,*) 'WARNING : no restart1file'
222 if(me.eq.king.or..not.out1file) then
223 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
224 write(iout,*) "i_index"
229 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
238 if (rest.and..not.restart1file) &
239 inquire(file=mremd_rst_name,exist=file_exist)
240 if(.not.file_exist.and.rest.and..not.restart1file) &
241 write(iout,*) 'WARNING : no restart file',mremd_rst_name
242 IF (rest.and.file_exist.and..not.restart1file) THEN
243 write (iout,*) 'Master is reading restart file',&
245 open(irest2,file=mremd_rst_name,status='unknown')
247 read (irest2,*) (i2rep(i),i=0,nodes-1)
249 read (irest2,*) (ifirst(i),i=1,remd_m(1))
252 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
254 read (irest2,*) ndowna(0,il),&
255 (ndowna(i,il),i=1,ndowna(0,il))
261 read (irest2,*) (mset(i),i=1,nset)
263 read (irest2,*) (i2set(i),i=0,nodes-1)
268 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
273 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
274 write(iout,*) "i_index"
279 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
288 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
289 write (iout,'(a6,1000i5)') "ifirst",&
290 (ifirst(i),i=1,remd_m(1))
292 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
293 (nupa(i,il),i=1,nupa(0,il))
294 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
295 (ndowna(i,il),i=1,ndowna(0,il))
297 ELSE IF (.not.(rest.and.file_exist)) THEN
303 if (i2rep(il-1).eq.nrep) then
306 nupa(0,il)=remd_m(i2rep(il-1)+1)
307 k=rep2i(int(i2rep(il-1)))+1
313 if (i2rep(il-1).eq.1) then
316 ndowna(0,il)=remd_m(i2rep(il-1)-1)
317 k=rep2i(i2rep(il-1)-2)+1
325 write (iout,'(a6,100i4)') "ifirst",&
326 (ifirst(i),i=1,remd_m(1))
328 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
329 (nupa(i,il),i=1,nupa(0,il))
330 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
331 (ndowna(i,il),i=1,ndowna(0,il))
337 ! t_bath=retmin+(retmax-retmin)*me/(nodes-1)
338 if(.not.(rest.and.file_exist.and.restart1file)) then
339 if (me .eq. king) then
342 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
344 !d print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
345 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
350 if(me.eq.king.or..not.out1file) &
351 write(iout,*) me,"iset=",iset,"t_bath=",t_bath
355 stdfp(i)=dsqrt(2*Rb*t_bath/d_time)
359 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
363 ! print *,'irep',me,t_bath
365 if (me.eq.king .or. .not. out1file) &
366 write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
367 call rescale_weights(t_bath)
371 !------copy MD--------------
372 ! The driver for molecular dynamics subroutines
373 !------------------------------------------------
379 if(me.eq.king.or..not.out1file) &
380 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
386 ! Determine the inverse of the inertia matrix.
387 call setup_MD_matrices
391 if (me.eq.king .or. .not. out1file) &
392 write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
394 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
396 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
399 call rescale_weights(t_bath)
403 t_MDsetup = MPI_Wtime()-tt0
405 t_MDsetup = tcpu()-tt0
408 ! Entering the MD loop
414 if (lang.eq.2 .or. lang.eq.3) then
418 call sd_verlet_p_setup
420 call sd_verlet_ciccotti_setup
424 pfric0_mat(i,j,0)=pfric_mat(i,j)
425 afric0_mat(i,j,0)=afric_mat(i,j)
426 vfric0_mat(i,j,0)=vfric_mat(i,j)
427 prand0_mat(i,j,0)=prand_mat(i,j)
428 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
429 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
434 flag_stoch(i)=.false.
438 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
440 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
444 else if (lang.eq.1 .or. lang.eq.4) then
448 if (me.eq.king .or. .not. out1file) &
449 write(iout,*) 'Setup time',time00-walltime
452 t_langsetup=MPI_Wtime()-tt0
455 t_langsetup=tcpu()-tt0
461 do while(.not.end_of_run)
463 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
464 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
466 if (lang.gt.0 .and. surfarea .and. &
467 mod(itime,reset_fricmat).eq.0) then
468 if (lang.eq.2 .or. lang.eq.3) then
472 call sd_verlet_p_setup
474 call sd_verlet_ciccotti_setup
478 pfric0_mat(i,j,0)=pfric_mat(i,j)
479 afric0_mat(i,j,0)=afric_mat(i,j)
480 vfric0_mat(i,j,0)=vfric_mat(i,j)
481 prand0_mat(i,j,0)=prand_mat(i,j)
482 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
483 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
488 flag_stoch(i)=.false.
491 else if (lang.eq.1 .or. lang.eq.4) then
494 write (iout,'(a,i10)') &
495 "Friction matrix reset based on surface area, itime",itime
497 if (reset_vel .and. tbf .and. lang.eq.0 &
498 .and. mod(itime,count_reset_vel).eq.0) then
500 if (me.eq.king .or. .not. out1file) &
501 write(iout,'(a,f20.2)') &
502 "Velocities reset to random values, time",totT
505 d_t_old(j,i)=d_t(j,i)
509 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
513 d_t(j,0)=d_t(j,0)-vcm(j)
516 kinetic_T=2.0d0/(dimen3*Rb)*EK
517 scalfac=dsqrt(T_bath/kinetic_T)
518 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
521 d_t_old(j,i)=scalfac*d_t(j,i)
527 ! Time-reversible RESPA algorithm
528 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
529 call RESPA_step(itime)
531 ! Variable time step algorithm.
532 call velverlet_step(itime)
536 call brown_step(itime)
538 print *,"Brown dynamics not here!"
540 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
546 if (mod(itime,ntwe).eq.0) call statout(itime)
548 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
549 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
551 call hairpin(.true.,nharp,iharp)
552 call secondary2(.true.)
553 call pdbout(potE,tytul,ipdb)
558 if (mod(itime,ntwx).eq.0.and.traj1file) then
559 if(ntwx_cache.lt.max_cache_traj_use) then
560 ntwx_cache=ntwx_cache+1
562 if (max_cache_traj_use.ne.1) &
563 print *,itime,"processor ",me," over cache ",ntwx_cache
566 totT_cache(i)=totT_cache(i+1)
567 EK_cache(i)=EK_cache(i+1)
568 potE_cache(i)=potE_cache(i+1)
569 t_bath_cache(i)=t_bath_cache(i+1)
570 Uconst_cache(i)=Uconst_cache(i+1)
571 iset_cache(i)=iset_cache(i+1)
574 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
577 qpair_cache(ii,i)=qpair_cache(ii,i+1)
580 utheta_cache(ii,i)=utheta_cache(ii,i+1)
581 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
582 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
588 c_cache(j,ii,i)=c_cache(j,ii,i+1)
594 totT_cache(ntwx_cache)=totT
595 EK_cache(ntwx_cache)=EK
596 potE_cache(ntwx_cache)=potE
597 t_bath_cache(ntwx_cache)=t_bath
598 Uconst_cache(ntwx_cache)=Uconst
599 iset_cache(ntwx_cache)=iset
602 qfrag_cache(i,ntwx_cache)=qfrag(i)
605 qpair_cache(i,ntwx_cache)=qpair(i)
608 utheta_cache(i,ntwx_cache)=utheta(i)
609 ugamma_cache(i,ntwx_cache)=ugamma(i)
610 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
615 c_cache(j,i,ntwx_cache)=c(j,i)
620 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
621 .and..not.restart1file) then
624 open(irest1,file=mremd_rst_name,status='unknown')
625 write (irest1,*) "i2rep"
626 write (irest1,*) (i2rep(i),i=0,nodes-1)
627 write (irest1,*) "ifirst"
628 write (irest1,*) (ifirst(i),i=1,remd_m(1))
630 write (irest1,*) "nupa",il
631 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
632 write (irest1,*) "ndowna",il
633 write (irest1,*) ndowna(0,il),&
634 (ndowna(i,il),i=1,ndowna(0,il))
637 write (irest1,*) "nset"
638 write (irest1,*) nset
639 write (irest1,*) "mset"
640 write (irest1,*) (mset(i),i=1,nset)
641 write (irest1,*) "i2set"
642 write (irest1,*) (i2set(i),i=0,nodes-1)
643 write (irest1,*) "i_index"
647 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
655 open(irest2,file=rest2name,status='unknown')
656 write(irest2,*) totT,EK,potE,totE,t_bath
658 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
661 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
664 write (irest2,*) iset
671 ! forced synchronization
672 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
673 .and. .not. mremdsync) then
675 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
677 call mpi_recv(itime_master, 1, MPI_INTEGER,&
678 0,101,CG_COMM, status, ierr)
679 call mpi_barrier(CG_COMM, ierr)
680 !deb if (out1file.or.traj1file) then
681 !deb call mpi_gather(itime,1,mpi_integer,
682 !deb & icache_all,1,mpi_integer,king,
685 call mpi_gather(ntwx_cache,1,mpi_integer,&
686 icache_all,1,mpi_integer,king,&
689 write(iout,*) 'REMD synchro at',itime_master,itime
690 if (itime_master.ge.n_timestep .or. ovrtim()) &
692 !time call flush(iout)
697 if ((mod(itime,nstex).eq.0.and.me.eq.king &
698 .or.end_of_run.and.me.eq.king ) &
699 .and. .not. mremdsync ) then
702 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
703 CG_COMM, ireqi(i), ierr)
704 !d write(iout,*) 'REMD synchro with',i
707 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
708 call mpi_barrier(CG_COMM, ierr)
710 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
711 if (out1file.or.traj1file) then
712 !deb call mpi_gather(itime,1,mpi_integer,
713 !deb & itime_all,1,mpi_integer,king,
715 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
716 !deb & (itime_all(i),i=1,nodes)
718 !deb imin_itime=itime_all(1)
720 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
722 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
723 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
724 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
725 call mpi_gather(ntwx_cache,1,mpi_integer,&
726 icache_all,1,mpi_integer,king,&
728 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
729 ! & (icache_all(i),i=1,nodes)
730 ii_write=icache_all(1)
732 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
734 ! write(iout,*) "MIN ii_write=",ii_write
737 !time call flush(iout)
739 if(mremdsync .and. mod(itime,nstex).eq.0) then
741 if (me.eq.king .or. .not. out1file) &
742 write(iout,*) 'REMD synchro at',itime
745 call mpi_gather(ntwx_cache,1,mpi_integer,&
746 icache_all,1,mpi_integer,king,&
749 write(iout,'(a19,8000i8)') ' ntwx_cache',&
750 (icache_all(i),i=1,nodes)
751 ii_write=icache_all(1)
753 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
755 write(iout,*) "MIN ii_write=",ii_write
761 ! Update the time safety limiy
762 if (time001-time00.gt.safety) then
763 safety=time001-time00+600
764 write (iout,*) "****** SAFETY increased to",safety," s"
766 if (ovrtim()) end_of_run=.true.
768 if(synflag.and..not.end_of_run) then
772 write(iout,*) 'REMD before',me,t_bath
774 ! call mpi_gather(t_bath,1,mpi_double_precision,
775 ! & remd_t_bath,1,mpi_double_precision,king,
777 potEcomp(n_ene+1)=t_bath
779 potEcomp(n_ene+2)=iset
780 if (iset.lt.nset) then
784 potEcomp(n_ene+3)=Uconst
791 potEcomp(n_ene+4)=Uconst
795 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
796 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
799 call mpi_gather(elow,1,mpi_double_precision,&
800 elowi,1,mpi_double_precision,king,&
802 call mpi_gather(ehigh,1,mpi_double_precision,&
803 ehighi,1,mpi_double_precision,king,&
808 if (me.eq.king .or. .not. out1file) then
809 write(iout,*) 'REMD gather times=',time03-time01 &
813 if (restart1file) call write1rst(i_index)
816 if (me.eq.king .or. .not. out1file) then
817 write(iout,*) 'REMD writing rst time=',time04-time03
820 if (traj1file) call write1traj
822 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
823 !deb & icache_all,1,mpi_integer,king,
825 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
826 !deb & (icache_all(i),i=1,nodes)
831 if (me.eq.king .or. .not. out1file) then
832 write(iout,*) 'REMD writing traj time=',time05-time04
839 remd_t_bath(i)=remd_ene(n_ene+1,i)
840 iremd_iset(i)=remd_ene(n_ene+2,i)
844 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
846 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
850 write(iout,*) 'REMD exchange temp,ene'
852 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
853 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
857 !-------------------------------------
859 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
862 write (iout,*) "remd_m(1)",remd_m(1)
864 i=ifirst(iran_num(1,remd_m(1)))
871 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
873 if(i.gt.0.and.nupa(0,i).gt.0) then
875 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
877 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
879 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
881 ! do while (iex.eq.i)
882 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
883 iex=nupa(iran_num(1,int(nupa(0,i))),i)
885 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
887 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
889 ! Swap temperatures between conformations i and iex with recalculating the free energies
890 ! following temperature changes.
891 ene_iex_iex=remd_ene(0,iex)
892 ene_i_i=remd_ene(0,i)
893 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
894 ! & " iex",iex," ene_iex_iex",ene_iex_iex
895 ! write (iout,*) "rescaling weights with temperature",
898 call rescale_weights(remd_t_bath(i))
900 ! write (iout,*) "0,iex",remd_t_bath(i)
901 ! call enerprint(remd_ene(0,iex))
903 call sum_energy(remd_ene(0,iex),.false.)
904 ene_iex_i=remd_ene(0,iex)
905 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
907 ! write (iout,*) "0,i",remd_t_bath(i)
908 ! call enerprint(remd_ene(0,i))
910 call sum_energy(remd_ene(0,i),.false.)
911 ! write (iout,*) "ene_i_i",remd_ene(0,i)
913 ! write (iout,*) "rescaling weights with temperature",
915 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
916 write (iout,*) "ERROR: inconsistent energies:",i,&
917 ene_i_i,remd_ene(0,i)
919 call rescale_weights(remd_t_bath(iex))
921 ! write (iout,*) "0,i",remd_t_bath(iex)
922 ! call enerprint(remd_ene(0,i))
924 call sum_energy(remd_ene(0,i),.false.)
925 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
927 ene_i_iex=remd_ene(0,i)
929 ! write (iout,*) "0,iex",remd_t_bath(iex)
930 ! call enerprint(remd_ene(0,iex))
932 call sum_energy(remd_ene(0,iex),.false.)
933 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
934 write (iout,*) "ERROR: inconsistent energies:",iex,&
935 ene_iex_iex,remd_ene(0,iex)
937 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
938 ! write (iout,*) "i",i," iex",iex
939 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
940 ! & " ene_i_iex",ene_i_iex,
941 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
943 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
944 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
946 ! write(iout,*) 'delta',delta
947 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
948 ! & (remd_ene(i)-remd_ene(iex))/Rb/
949 ! & (remd_t_bath(i)*remd_t_bath(iex))
951 if (delta .gt. 50.0d0) then
957 else if (delta.lt.-50.0d0) then
966 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
967 xxx=ran_number(0.0d0,1.0d0)
968 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
970 if (delta .gt. xxx) then
972 remd_t_bath(i)=remd_t_bath(iex)
974 remd_ene(0,i)=ene_i_iex
975 remd_ene(0,iex)=ene_iex_i
981 ehighi(i)=ehighi(iex)
988 nupa(k,i)=nupa(k,iex)
991 ndowna(k,i)=ndowna(k,iex)
995 if (ifirst(il).eq.i) ifirst(il)=iex
997 if (nupa(k,il).eq.i) then
999 elseif (nupa(k,il).eq.iex) then
1004 if (ndowna(k,il).eq.i) then
1006 elseif (ndowna(k,il).eq.iex) then
1012 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1014 i2rep(i-1)=i2rep(iex-1)
1017 ! write(iout,*) 'exchange',i,iex
1018 ! write (iout,'(a8,100i4)') "@ ifirst",
1019 ! & (ifirst(k),k=1,remd_m(1))
1021 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1022 ! & (nupa(k,il),k=1,nupa(0,il))
1023 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1024 ! & (ndowna(k,il),k=1,ndowna(0,il))
1029 remd_ene(0,iex)=ene_iex_iex
1030 remd_ene(0,i)=ene_i_i
1036 !d write (iout,*) "exchange completed"
1040 !d write(iout,*) "########",ii
1042 i_temp=iran_num(1,nrep)
1043 i_mult=iran_num(1,remd_m(i_temp))
1044 i_iset=iran_num(1,nset)
1045 i_mset=iran_num(1,mset(i_iset))
1046 i=i_index(i_temp,i_mult,i_iset,i_mset)
1048 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1051 !d write(iout,*) "i_dir=",i_dir
1053 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1056 i_mult1=iran_num(1,remd_m(i_temp1))
1058 i_mset1=iran_num(1,mset(i_iset1))
1059 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1061 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1064 i_mult1=iran_num(1,remd_m(i_temp1))
1066 i_mset1=iran_num(1,mset(i_iset1))
1067 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1068 econstr_temp_i=remd_ene(20,i)
1069 econstr_temp_iex=remd_ene(20,iex)
1070 remd_ene(20,i)=remd_ene(n_ene+3,i)
1071 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1073 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1076 i_mult1=iran_num(1,remd_m(i_temp1))
1078 i_mset1=iran_num(1,mset(i_iset1))
1079 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1080 econstr_temp_i=remd_ene(20,i)
1081 econstr_temp_iex=remd_ene(20,iex)
1082 remd_ene(20,i)=remd_ene(n_ene+3,i)
1083 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1089 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1092 ! Swap temperatures between conformations i and iex with recalculating the free energies
1093 ! following temperature changes.
1094 ene_iex_iex=remd_ene(0,iex)
1095 ene_i_i=remd_ene(0,i)
1096 !o write (iout,*) "rescaling weights with temperature",
1098 call rescale_weights(remd_t_bath(i))
1100 call sum_energy(remd_ene(0,iex),.false.)
1101 ene_iex_i=remd_ene(0,iex)
1102 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1103 ! call sum_energy(remd_ene(0,i),.false.)
1104 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1105 ! write (iout,*) "rescaling weights with temperature",
1106 ! & remd_t_bath(iex)
1107 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1108 ! write (iout,*) "ERROR: inconsistent energies:",i,
1109 ! & ene_i_i,remd_ene(0,i)
1111 call rescale_weights(remd_t_bath(iex))
1112 call sum_energy(remd_ene(0,i),.false.)
1113 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1114 ene_i_iex=remd_ene(0,i)
1115 ! call sum_energy(remd_ene(0,iex),.false.)
1116 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1117 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1118 ! & ene_iex_iex,remd_ene(0,iex)
1120 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1121 ! write (iout,*) "i",i," iex",iex
1122 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1123 !d & " ene_i_iex",ene_i_iex,
1124 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1125 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1126 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1128 !d write(iout,*) 'delta',delta
1129 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1130 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1131 ! & (remd_t_bath(i)*remd_t_bath(iex))
1132 if (delta .gt. 50.0d0) then
1137 if (i_dir.eq.1.or.i_dir.eq.3) &
1138 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1139 if (i_dir.eq.2.or.i_dir.eq.3) &
1140 iremd_tot_usa(int(i2set(i-1)))= &
1141 iremd_tot_usa(int(i2set(i-1)))+1
1142 xxx=ran_number(0.0d0,1.0d0)
1143 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1144 if (delta .gt. xxx) then
1146 remd_t_bath(i)=remd_t_bath(iex)
1147 remd_t_bath(iex)=tmp
1150 iremd_iset(i)=iremd_iset(iex)
1151 iremd_iset(iex)=itmp
1153 remd_ene(0,i)=ene_i_iex
1154 remd_ene(0,iex)=ene_iex_i
1156 if (i_dir.eq.1.or.i_dir.eq.3) &
1157 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1160 i2rep(i-1)=i2rep(iex-1)
1163 if (i_dir.eq.2.or.i_dir.eq.3) &
1164 iremd_acc_usa(int(i2set(i-1)))= &
1165 iremd_acc_usa(int(i2set(i-1)))+1
1168 i2set(i-1)=i2set(iex-1)
1171 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1172 i_index(i_temp,i_mult,i_iset,i_mset)= &
1173 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1174 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1177 remd_ene(0,iex)=ene_iex_iex
1178 remd_ene(0,i)=ene_i_i
1179 remd_ene(20,iex)=econstr_temp_iex
1180 remd_ene(20,i)=econstr_temp_i
1184 !d do il1=1,mset(il)
1187 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1200 !-------------------------------------
1201 write (iout,*) "NREP",nrep
1203 if(iremd_tot(i).ne.0) &
1204 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1205 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1210 if(iremd_tot_usa(i).ne.0) &
1211 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1212 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1218 !d write (iout,'(a6,100i4)') "ifirst",
1219 !d & (ifirst(i),i=1,remd_m(1))
1221 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1222 !d & (nupa(i,il),i=1,nupa(0,il))
1223 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1224 !d & (ndowna(i,il),i=1,ndowna(0,il))
1229 !d write (iout,*) "Before scatter"
1232 if (me.eq.king) then
1233 write (iout,*) "t_bath before scatter",remd_t_bath
1237 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1238 t_bath,1,mpi_double_precision,king,&
1240 !d write (iout,*) "After scatter"
1243 call mpi_scatter(iremd_iset,1,mpi_integer,&
1244 iset,1,mpi_integer,king,&
1248 if (me.eq.king .or. .not. out1file) then
1249 write(iout,*) 'REMD scatter time=',time07-time06
1253 call mpi_scatter(elowi,1,mpi_double_precision,&
1254 elow,1,mpi_double_precision,king,&
1256 call mpi_scatter(ehighi,1,mpi_double_precision,&
1257 ehigh,1,mpi_double_precision,king,&
1260 call rescale_weights(t_bath)
1261 !o write (iout,*) "Processor",me,
1262 !o & " rescaling weights with temperature",t_bath
1264 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
1266 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
1270 !de write(iout,*) 'REMD after',me,t_bath
1272 if (me.eq.king .or. .not. out1file) then
1273 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1274 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1275 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1276 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1277 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1278 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1279 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1280 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1281 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1287 if (restart1file) then
1288 if (me.eq.king .or. .not. out1file) &
1289 write(iout,*) 'writing restart at the end of run'
1290 call write1rst(i_index)
1293 if (traj1file) call write1traj
1295 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1296 !deb & icache_all,1,mpi_integer,king,
1297 !deb & CG_COMM,ierr)
1298 !deb write(iout,'(a40,8000i8)')
1299 !deb & ' ntwx_cache after traj1file at the end',
1300 !deb & (icache_all(i),i=1,nodes)
1305 t_MD=MPI_Wtime()-tt0
1309 if (me.eq.king .or. .not. out1file) then
1310 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1312 'MD calculations setup:',t_MDsetup,&
1313 'Energy & gradient evaluation:',t_enegrad,&
1314 'Stochastic MD setup:',t_langsetup,&
1315 'Stochastic MD step setup:',t_sdsetup,&
1317 write (iout,'(/28(1h=),a25,27(1h=))') &
1318 ' End of MD calculation '
1320 !el common /przechowalnia/
1321 ! deallocate(d_restart1)
1322 ! deallocate(d_restart2)
1326 end subroutine MREMD
1327 !-----------------------------------------------------------------------------
1328 subroutine write1rst(i_index)
1331 ! implicit real*8 (a-h,o-z)
1332 ! include 'DIMENSIONS'
1334 ! include 'COMMON.MD'
1335 ! include 'COMMON.IOUNITS'
1336 ! include 'COMMON.REMD'
1337 ! include 'COMMON.SETUP'
1338 ! include 'COMMON.CHAIN'
1339 ! include 'COMMON.SBRIDGE'
1340 ! include 'COMMON.INTERACT'
1342 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1343 !el d_restart2(3,2*nres*maxprocs)
1344 real(kind=4) :: r_d(3,0:2*nres)
1345 real(kind=4) :: t5_restart1(5)
1346 integer :: iret,itmp
1347 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1348 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1350 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1351 !el common /przechowalnia/ d_restart1,d_restart2
1352 integer :: i,j,il,il1,ierr,ixdrf
1357 t5_restart1(4)=t_bath
1358 t5_restart1(5)=Uconst
1360 call mpi_gather(t5_restart1,5,mpi_real,&
1361 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1369 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1370 d_restart1,3*2*nres+3,mpi_real,king,&
1382 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1383 d_restart2,3*2*nres+3,mpi_real,king,&
1388 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1390 call xdrfint_(ixdrf, i2rep(i), iret)
1393 call xdrfint_(ixdrf, ifirst(i), iret)
1397 call xdrfint_(ixdrf, nupa(i,il), iret)
1401 call xdrfint_(ixdrf, ndowna(i,il), iret)
1407 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1414 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1421 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1427 call xdrfint_(ixdrf, nset, iret)
1429 call xdrfint_(ixdrf,mset(i), iret)
1432 call xdrfint_(ixdrf,i2set(i), iret)
1438 itmp=i_index(i,j,il,il1)
1439 call xdrfint_(ixdrf,itmp, iret)
1446 call xdrfclose_(ixdrf, iret)
1448 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1450 call xdrfint(ixdrf, i2rep(i), iret)
1453 call xdrfint(ixdrf, ifirst(i), iret)
1457 call xdrfint(ixdrf, nupa(i,il), iret)
1461 call xdrfint(ixdrf, ndowna(i,il), iret)
1467 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1474 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1481 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1488 call xdrfint(ixdrf, nset, iret)
1490 call xdrfint(ixdrf,mset(i), iret)
1493 call xdrfint(ixdrf,i2set(i), iret)
1499 itmp=i_index(i,j,il,il1)
1500 call xdrfint(ixdrf,itmp, iret)
1507 call xdrfclose(ixdrf, iret)
1511 end subroutine write1rst
1512 !-----------------------------------------------------------------------------
1513 subroutine write1traj
1515 ! implicit real*8 (a-h,o-z)
1516 ! include 'DIMENSIONS'
1518 ! include 'COMMON.MD'
1519 ! include 'COMMON.IOUNITS'
1520 ! include 'COMMON.REMD'
1521 ! include 'COMMON.SETUP'
1522 ! include 'COMMON.CHAIN'
1523 ! include 'COMMON.SBRIDGE'
1524 ! include 'COMMON.INTERACT'
1526 real(kind=4) :: t5_restart1(5)
1527 integer :: iret,itmp
1528 real(kind=4) :: xcoord(3,2*nres+2),prec
1529 real(kind=4) :: r_qfrag(50),r_qpair(100)
1530 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1531 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1532 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1533 p_uscdiff(100*Nprocs) !(100*maxprocs)
1534 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1535 real(kind=4) :: r_c(3,2*nres+2)
1536 !el common /przechowalnia/ p_c
1538 integer :: i,j,il,ierr,ii,ixdrf
1540 call mpi_bcast(ii_write,1,mpi_integer,&
1544 print *,'traj1file',me,ii_write,ntwx_cache
1548 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1550 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1553 ! write (iout,*) "before gather write1traj: from node",ii
1555 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1557 t5_restart1(1)=totT_cache(ii)
1558 t5_restart1(2)=EK_cache(ii)
1559 t5_restart1(3)=potE_cache(ii)
1560 t5_restart1(4)=t_bath_cache(ii)
1561 t5_restart1(5)=Uconst_cache(ii)
1562 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1564 call mpi_gather(t5_restart1,5,mpi_real,&
1565 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1567 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1570 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1571 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1574 r_qfrag(i)=qfrag_cache(i,ii)
1577 r_qpair(i)=qpair_cache(i,ii)
1580 r_utheta(i)=utheta_cache(i,ii)
1581 r_ugamma(i)=ugamma_cache(i,ii)
1582 r_uscdiff(i)=uscdiff_cache(i,ii)
1585 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1586 p_qfrag,nfrag,mpi_real,king,&
1588 call mpi_gather(r_qpair,npair,mpi_real,&
1589 p_qpair,npair,mpi_real,king,&
1591 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1592 p_utheta,nfrag_back,mpi_real,king,&
1594 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1595 p_ugamma,nfrag_back,mpi_real,king,&
1597 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1598 p_uscdiff,nfrag_back,mpi_real,king,&
1602 write (iout,*) "p_qfrag"
1604 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1606 write (iout,*) "p_qpair"
1608 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1614 r_c(j,i)=c_cache(j,i,ii)
1618 call mpi_gather(r_c,3*2*nres,mpi_real,&
1619 p_c,3*2*nres,mpi_real,king,&
1625 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1626 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1627 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1628 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1629 call xdrfint_(ixdrf, nss, iret)
1632 call xdrfint(ixdrf, idssb(j)+nres, iret)
1633 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1635 call xdrfint_(ixdrf, ihpb(j), iret)
1636 call xdrfint_(ixdrf, jhpb(j), iret)
1639 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1640 call xdrfint_(ixdrf, iset_restart1(il), iret)
1642 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1645 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1648 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1649 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1650 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1655 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1660 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1664 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1668 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1669 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1670 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1671 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1672 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1673 call xdrfint(ixdrf, nss, iret)
1676 call xdrfint(ixdrf, idssb(j)+nres, iret)
1677 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1679 call xdrfint(ixdrf, ihpb(j), iret)
1680 call xdrfint(ixdrf, jhpb(j), iret)
1683 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1684 call xdrfint(ixdrf, iset_restart1(il), iret)
1686 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1689 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1692 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1693 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1694 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1699 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1704 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1708 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1714 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1716 if(me.eq.king) call xdrfclose(ixdrf, iret)
1718 do i=1,ntwx_cache-ii_write
1720 totT_cache(i)=totT_cache(ii_write+i)
1721 EK_cache(i)=EK_cache(ii_write+i)
1722 potE_cache(i)=potE_cache(ii_write+i)
1723 t_bath_cache(i)=t_bath_cache(ii_write+i)
1724 Uconst_cache(i)=Uconst_cache(ii_write+i)
1725 iset_cache(i)=iset_cache(ii_write+i)
1728 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1731 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1734 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1735 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1736 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1741 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1745 ntwx_cache=ntwx_cache-ii_write
1747 end subroutine write1traj
1748 !-----------------------------------------------------------------------------
1749 subroutine read1restart(i_index)
1751 ! implicit real*8 (a-h,o-z)
1752 ! include 'DIMENSIONS'
1754 ! include 'COMMON.MD'
1755 ! include 'COMMON.IOUNITS'
1756 ! include 'COMMON.REMD'
1757 ! include 'COMMON.SETUP'
1758 ! include 'COMMON.CHAIN'
1759 ! include 'COMMON.SBRIDGE'
1760 ! include 'COMMON.INTERACT'
1761 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1762 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1763 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1764 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1766 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1767 !el common /przechowalnia/ d_restart1
1768 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1770 write (*,*) "Processor",me," called read1restart"
1773 open(irest2,file=mremd_rst_name,status='unknown')
1774 read(irest2,*,err=334) i
1775 write(iout,*) "Reading old rst in ASCI format"
1777 call read1restart_old
1781 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1784 call xdrfint_(ixdrf, i2rep(i), iret)
1787 call xdrfint_(ixdrf, ifirst(i), iret)
1790 call xdrfint_(ixdrf, nupa(0,il), iret)
1792 call xdrfint_(ixdrf, nupa(i,il), iret)
1795 call xdrfint_(ixdrf, ndowna(0,il), iret)
1797 call xdrfint_(ixdrf, ndowna(i,il), iret)
1802 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1806 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1809 call xdrfint(ixdrf, i2rep(i), iret)
1812 call xdrfint(ixdrf, ifirst(i), iret)
1815 call xdrfint(ixdrf, nupa(0,il), iret)
1817 call xdrfint(ixdrf, nupa(i,il), iret)
1820 call xdrfint(ixdrf, ndowna(0,il), iret)
1822 call xdrfint(ixdrf, ndowna(i,il), iret)
1827 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1832 call mpi_scatter(t_restart1,5,mpi_real,&
1833 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1837 t_bath=t5_restart1(4)
1842 ! read(irest2,'(3e15.5)')
1843 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1846 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1848 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1854 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1855 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1865 ! read(irest2,'(3e15.5)')
1866 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1869 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1871 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1877 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1878 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1889 call xdrfint_(ixdrf, nset, iret)
1891 call xdrfint_(ixdrf,mset(i), iret)
1894 call xdrfint_(ixdrf,i2set(i), iret)
1900 call xdrfint_(ixdrf,itmp, iret)
1901 i_index(i,j,il,il1)=itmp
1909 call xdrfint(ixdrf, nset, iret)
1911 call xdrfint(ixdrf,mset(i), iret)
1914 call xdrfint(ixdrf,i2set(i), iret)
1920 call xdrfint(ixdrf,itmp, iret)
1921 i_index(i,j,il,il1)=itmp
1928 call mpi_scatter(i2set,1,mpi_integer,&
1929 iset,1,mpi_integer,king,&
1934 if(me.eq.king) close(irest2)
1936 end subroutine read1restart
1937 !-----------------------------------------------------------------------------
1938 subroutine read1restart_old
1940 ! implicit real*8 (a-h,o-z)
1941 ! include 'DIMENSIONS'
1943 ! include 'COMMON.MD'
1944 ! include 'COMMON.IOUNITS'
1945 ! include 'COMMON.REMD'
1946 ! include 'COMMON.SETUP'
1947 ! include 'COMMON.CHAIN'
1948 ! include 'COMMON.SBRIDGE'
1949 ! include 'COMMON.INTERACT'
1950 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1951 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1952 !el common /przechowalnia/ d_restart1
1954 integer :: i,j,il,ierr
1957 open(irest2,file=mremd_rst_name,status='unknown')
1958 read (irest2,*) (i2rep(i),i=0,nodes-1)
1959 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1961 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1962 read (irest2,*) ndowna(0,il),&
1963 (ndowna(i,il),i=1,ndowna(0,il))
1966 read(irest2,*) (t_restart1(j,il),j=1,4)
1969 call mpi_scatter(t_restart1,5,mpi_real,&
1970 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1974 t_bath=t5_restart1(4)
1979 read(irest2,'(3e15.5)') &
1980 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
1984 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1985 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1995 read(irest2,'(3e15.5)') &
1996 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2000 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2001 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2007 if(me.eq.king) close(irest2)
2009 end subroutine read1restart_old
2010 !----------------------------------------------------------------
2011 subroutine alloc_MREMD_arrays
2013 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2014 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
2016 ! common /remdcommon/ in io: read_REMDpar
2017 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2018 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2019 ! common /remdrestart/
2020 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2022 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2023 allocate(ifirst(0:nodes)) !(maxprocs)
2024 allocate(nupa(0:nodes,0:2*nodes))
2025 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2026 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2027 allocate(iset_restart1(nodes)) !(maxprocs)
2028 ! common /traj1cache/
2029 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2030 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2031 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2032 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2033 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2034 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2035 allocate(utheta_cache(nfrag_back,max_cache_traj))
2036 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2037 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2038 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2041 end subroutine alloc_MREMD_arrays
2042 !-----------------------------------------------------------------------------
2043 !-----------------------------------------------------------------------------