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) :: 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 print *,'MREMD',nodes
131 print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
132 write (iout,*) "Start MREMD: me",me," t_bath",t_bath
133 print *,"NSET=",nset, "MSET=", mset
137 do il1=1,max0(mset(il),1)
147 i_index(i,j,il,il1)=k
153 if(me.eq.king.or..not.out1file) then
154 write(iout,*) (i2rep(i),i=0,nodes-1)
155 write(iout,*) (i2set(i),i=0,nodes-1)
160 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
167 ! print *,'i2rep',me,i2rep(me)
168 ! print *,'rep2i',(rep2i(i),i=0,nrep)
170 !old if (i2rep(me).eq.nrep) then
173 !old nup(0)=remd_m(i2rep(me)+1)
174 !old k=rep2i(int(i2rep(me)))+1
181 !d print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
183 !old if (i2rep(me).eq.1) then
186 !old ndown(0)=remd_m(i2rep(me)-1)
187 !old k=rep2i(i2rep(me)-2)+1
194 !d print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
196 !el common /przechowalnia/
197 if(.not.allocated(d_restart1)) allocate(d_restart1(3,0:(nres2+1)*nodes))
198 if(.not.allocated(d_restart2)) allocate(d_restart2(3,0:(nres2+1)*nodes))
199 if(.not.allocated(p_c)) allocate(p_c(3,(nres2+2)*nodes))
202 write (*,*) "Processor",me," rest",rest,&
203 "restart1fie",restart1file
204 if(rest.and.restart1file) then
206 inquire(file=mremd_rst_name,exist=file_exist)
207 !d write (*,*) me," Before broadcast: file_exist",file_exist
208 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
210 !d write (*,*) me," After broadcast: file_exist",file_exist
212 if(me.eq.king.or..not.out1file) &
213 write (iout,*) 'Master is reading restart1file'
214 call read1restart(i_index)
216 if(me.eq.king.or..not.out1file) &
217 write (iout,*) 'WARNING : no restart1file'
220 if(me.eq.king.or..not.out1file) then
221 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
222 write(iout,*) "i_index"
227 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
236 if (rest.and..not.restart1file) &
237 inquire(file=mremd_rst_name,exist=file_exist)
238 if(.not.file_exist.and.rest.and..not.restart1file) &
239 write(iout,*) 'WARNING : no restart file',mremd_rst_name
240 IF (rest.and.file_exist.and..not.restart1file) THEN
241 write (iout,*) 'Master is reading restart file',&
243 open(irest2,file=mremd_rst_name,status='unknown')
245 read (irest2,*) (i2rep(i),i=0,nodes-1)
247 read (irest2,*) (ifirst(i),i=1,remd_m(1))
250 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
252 read (irest2,*) ndowna(0,il),&
253 (ndowna(i,il),i=1,ndowna(0,il))
259 read (irest2,*) (mset(i),i=1,nset)
261 read (irest2,*) (i2set(i),i=0,nodes-1)
266 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
271 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
272 write(iout,*) "i_index"
277 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
286 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
287 write (iout,'(a6,1000i5)') "ifirst",&
288 (ifirst(i),i=1,remd_m(1))
290 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
291 (nupa(i,il),i=1,nupa(0,il))
292 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
293 (ndowna(i,il),i=1,ndowna(0,il))
295 ELSE IF (.not.(rest.and.file_exist)) THEN
301 if (i2rep(il-1).eq.nrep) then
304 nupa(0,il)=remd_m(i2rep(il-1)+1)
305 k=rep2i(int(i2rep(il-1)))+1
311 if (i2rep(il-1).eq.1) then
314 ndowna(0,il)=remd_m(i2rep(il-1)-1)
315 k=rep2i(i2rep(il-1)-2)+1
323 write (iout,'(a6,100i4)') "ifirst",&
324 (ifirst(i),i=1,remd_m(1))
326 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
327 (nupa(i,il),i=1,nupa(0,il))
328 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
329 (ndowna(i,il),i=1,ndowna(0,il))
335 ! t_bath=retmin+(retmax-retmin)*me/(nodes-1)
336 if(.not.(rest.and.file_exist.and.restart1file)) then
337 if (me .eq. king) then
340 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
342 !d print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
343 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
348 if(me.eq.king.or..not.out1file) &
349 write(iout,*) me,"iset=",iset,"t_bath=",t_bath
353 stdfp(i)=dsqrt(2*Rb*t_bath/d_time)
357 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
361 ! print *,'irep',me,t_bath
363 if (me.eq.king .or. .not. out1file) &
364 write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
365 call rescale_weights(t_bath)
369 !------copy MD--------------
370 ! The driver for molecular dynamics subroutines
371 !------------------------------------------------
377 if(me.eq.king.or..not.out1file) &
378 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
384 ! Determine the inverse of the inertia matrix.
385 call setup_MD_matrices
389 if (me.eq.king .or. .not. out1file) &
390 write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
392 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
394 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
397 call rescale_weights(t_bath)
401 t_MDsetup = MPI_Wtime()-tt0
403 t_MDsetup = tcpu()-tt0
406 ! Entering the MD loop
412 if (lang.eq.2 .or. lang.eq.3) then
416 call sd_verlet_p_setup
418 call sd_verlet_ciccotti_setup
422 pfric0_mat(i,j,0)=pfric_mat(i,j)
423 afric0_mat(i,j,0)=afric_mat(i,j)
424 vfric0_mat(i,j,0)=vfric_mat(i,j)
425 prand0_mat(i,j,0)=prand_mat(i,j)
426 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
427 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
432 flag_stoch(i)=.false.
436 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
438 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
442 else if (lang.eq.1 .or. lang.eq.4) then
446 if (me.eq.king .or. .not. out1file) &
447 write(iout,*) 'Setup time',time00-walltime
450 t_langsetup=MPI_Wtime()-tt0
453 t_langsetup=tcpu()-tt0
459 do while(.not.end_of_run)
461 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
462 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
464 if (lang.gt.0 .and. surfarea .and. &
465 mod(itime,reset_fricmat).eq.0) then
466 if (lang.eq.2 .or. lang.eq.3) then
470 call sd_verlet_p_setup
472 call sd_verlet_ciccotti_setup
476 pfric0_mat(i,j,0)=pfric_mat(i,j)
477 afric0_mat(i,j,0)=afric_mat(i,j)
478 vfric0_mat(i,j,0)=vfric_mat(i,j)
479 prand0_mat(i,j,0)=prand_mat(i,j)
480 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
481 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
486 flag_stoch(i)=.false.
489 else if (lang.eq.1 .or. lang.eq.4) then
492 write (iout,'(a,i10)') &
493 "Friction matrix reset based on surface area, itime",itime
495 if (reset_vel .and. tbf .and. lang.eq.0 &
496 .and. mod(itime,count_reset_vel).eq.0) then
498 if (me.eq.king .or. .not. out1file) &
499 write(iout,'(a,f20.2)') &
500 "Velocities reset to random values, time",totT
503 d_t_old(j,i)=d_t(j,i)
507 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
511 d_t(j,0)=d_t(j,0)-vcm(j)
514 kinetic_T=2.0d0/(dimen3*Rb)*EK
515 scalfac=dsqrt(T_bath/kinetic_T)
516 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
519 d_t_old(j,i)=scalfac*d_t(j,i)
525 ! Time-reversible RESPA algorithm
526 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
527 call RESPA_step(itime)
529 ! Variable time step algorithm.
530 call velverlet_step(itime)
534 call brown_step(itime)
536 print *,"Brown dynamics not here!"
538 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
544 if (mod(itime,ntwe).eq.0) call statout(itime)
546 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
547 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
549 call hairpin(.true.,nharp,iharp)
550 call secondary2(.true.)
551 call pdbout(potE,tytul,ipdb)
556 if (mod(itime,ntwx).eq.0.and.traj1file) then
557 if(ntwx_cache.lt.max_cache_traj_use) then
558 ntwx_cache=ntwx_cache+1
560 if (max_cache_traj_use.ne.1) &
561 print *,itime,"processor ",me," over cache ",ntwx_cache
564 totT_cache(i)=totT_cache(i+1)
565 EK_cache(i)=EK_cache(i+1)
566 potE_cache(i)=potE_cache(i+1)
567 t_bath_cache(i)=t_bath_cache(i+1)
568 Uconst_cache(i)=Uconst_cache(i+1)
569 iset_cache(i)=iset_cache(i+1)
572 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
575 qpair_cache(ii,i)=qpair_cache(ii,i+1)
578 utheta_cache(ii,i)=utheta_cache(ii,i+1)
579 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
580 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
586 c_cache(j,ii,i)=c_cache(j,ii,i+1)
592 totT_cache(ntwx_cache)=totT
593 EK_cache(ntwx_cache)=EK
594 potE_cache(ntwx_cache)=potE
595 t_bath_cache(ntwx_cache)=t_bath
596 Uconst_cache(ntwx_cache)=Uconst
597 iset_cache(ntwx_cache)=iset
600 qfrag_cache(i,ntwx_cache)=qfrag(i)
603 qpair_cache(i,ntwx_cache)=qpair(i)
606 utheta_cache(i,ntwx_cache)=utheta(i)
607 ugamma_cache(i,ntwx_cache)=ugamma(i)
608 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
613 c_cache(j,i,ntwx_cache)=c(j,i)
618 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
619 .and..not.restart1file) then
622 open(irest1,file=mremd_rst_name,status='unknown')
623 write (irest1,*) "i2rep"
624 write (irest1,*) (i2rep(i),i=0,nodes-1)
625 write (irest1,*) "ifirst"
626 write (irest1,*) (ifirst(i),i=1,remd_m(1))
628 write (irest1,*) "nupa",il
629 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
630 write (irest1,*) "ndowna",il
631 write (irest1,*) ndowna(0,il),&
632 (ndowna(i,il),i=1,ndowna(0,il))
635 write (irest1,*) "nset"
636 write (irest1,*) nset
637 write (irest1,*) "mset"
638 write (irest1,*) (mset(i),i=1,nset)
639 write (irest1,*) "i2set"
640 write (irest1,*) (i2set(i),i=0,nodes-1)
641 write (irest1,*) "i_index"
645 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
653 open(irest2,file=rest2name,status='unknown')
654 write(irest2,*) totT,EK,potE,totE,t_bath
656 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
659 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
662 write (irest2,*) iset
669 ! forced synchronization
670 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
671 .and. .not. mremdsync) then
673 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
675 call mpi_recv(itime_master, 1, MPI_INTEGER,&
676 0,101,CG_COMM, status, ierr)
677 call mpi_barrier(CG_COMM, ierr)
678 !deb if (out1file.or.traj1file) then
679 !deb call mpi_gather(itime,1,mpi_integer,
680 !deb & icache_all,1,mpi_integer,king,
683 call mpi_gather(ntwx_cache,1,mpi_integer,&
684 icache_all,1,mpi_integer,king,&
687 write(iout,*) 'REMD synchro at',itime_master,itime
688 if (itime_master.ge.n_timestep .or. ovrtim()) &
690 !time call flush(iout)
695 if ((mod(itime,nstex).eq.0.and.me.eq.king &
696 .or.end_of_run.and.me.eq.king ) &
697 .and. .not. mremdsync ) then
700 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
701 CG_COMM, ireqi(i), ierr)
702 !d write(iout,*) 'REMD synchro with',i
705 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
706 call mpi_barrier(CG_COMM, ierr)
708 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
709 if (out1file.or.traj1file) then
710 !deb call mpi_gather(itime,1,mpi_integer,
711 !deb & itime_all,1,mpi_integer,king,
713 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
714 !deb & (itime_all(i),i=1,nodes)
716 !deb imin_itime=itime_all(1)
718 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
720 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
721 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
722 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
723 call mpi_gather(ntwx_cache,1,mpi_integer,&
724 icache_all,1,mpi_integer,king,&
726 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
727 ! & (icache_all(i),i=1,nodes)
728 ii_write=icache_all(1)
730 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
732 ! write(iout,*) "MIN ii_write=",ii_write
735 !time call flush(iout)
737 if(mremdsync .and. mod(itime,nstex).eq.0) then
739 if (me.eq.king .or. .not. out1file) &
740 write(iout,*) 'REMD synchro at',itime
743 call mpi_gather(ntwx_cache,1,mpi_integer,&
744 icache_all,1,mpi_integer,king,&
747 write(iout,'(a19,8000i8)') ' ntwx_cache',&
748 (icache_all(i),i=1,nodes)
749 ii_write=icache_all(1)
751 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
753 write(iout,*) "MIN ii_write=",ii_write
759 ! Update the time safety limiy
760 if (time001-time00.gt.safety) then
761 safety=time001-time00+600
762 write (iout,*) "****** SAFETY increased to",safety," s"
764 if (ovrtim()) end_of_run=.true.
766 if(synflag.and..not.end_of_run) then
770 write(iout,*) 'REMD before',me,t_bath
772 ! call mpi_gather(t_bath,1,mpi_double_precision,
773 ! & remd_t_bath,1,mpi_double_precision,king,
775 potEcomp(n_ene+1)=t_bath
777 potEcomp(n_ene+2)=iset
778 if (iset.lt.nset) then
782 potEcomp(n_ene+3)=Uconst
789 potEcomp(n_ene+4)=Uconst
793 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
794 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
797 call mpi_gather(elow,1,mpi_double_precision,&
798 elowi,1,mpi_double_precision,king,&
800 call mpi_gather(ehigh,1,mpi_double_precision,&
801 ehighi,1,mpi_double_precision,king,&
806 if (me.eq.king .or. .not. out1file) then
807 write(iout,*) 'REMD gather times=',time03-time01 &
811 if (restart1file) call write1rst(i_index)
814 if (me.eq.king .or. .not. out1file) then
815 write(iout,*) 'REMD writing rst time=',time04-time03
818 if (traj1file) call write1traj
820 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
821 !deb & icache_all,1,mpi_integer,king,
823 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
824 !deb & (icache_all(i),i=1,nodes)
829 if (me.eq.king .or. .not. out1file) then
830 write(iout,*) 'REMD writing traj time=',time05-time04
837 remd_t_bath(i)=remd_ene(n_ene+1,i)
838 iremd_iset(i)=remd_ene(n_ene+2,i)
842 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
844 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
848 write(iout,*) 'REMD exchange temp,ene'
850 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
851 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
855 !-------------------------------------
857 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
860 write (iout,*) "remd_m(1)",remd_m(1)
862 i=ifirst(iran_num(1,remd_m(1)))
869 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
871 if(i.gt.0.and.nupa(0,i).gt.0) then
873 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
875 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
877 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
879 ! do while (iex.eq.i)
880 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
881 iex=nupa(iran_num(1,int(nupa(0,i))),i)
883 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
885 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
887 ! Swap temperatures between conformations i and iex with recalculating the free energies
888 ! following temperature changes.
889 ene_iex_iex=remd_ene(0,iex)
890 ene_i_i=remd_ene(0,i)
891 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
892 ! & " iex",iex," ene_iex_iex",ene_iex_iex
893 ! write (iout,*) "rescaling weights with temperature",
896 call rescale_weights(remd_t_bath(i))
898 ! write (iout,*) "0,iex",remd_t_bath(i)
899 ! call enerprint(remd_ene(0,iex))
901 call sum_energy(remd_ene(0,iex),.false.)
902 ene_iex_i=remd_ene(0,iex)
903 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
905 ! write (iout,*) "0,i",remd_t_bath(i)
906 ! call enerprint(remd_ene(0,i))
908 call sum_energy(remd_ene(0,i),.false.)
909 ! write (iout,*) "ene_i_i",remd_ene(0,i)
911 ! write (iout,*) "rescaling weights with temperature",
913 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
914 write (iout,*) "ERROR: inconsistent energies:",i,&
915 ene_i_i,remd_ene(0,i)
917 call rescale_weights(remd_t_bath(iex))
919 ! write (iout,*) "0,i",remd_t_bath(iex)
920 ! call enerprint(remd_ene(0,i))
922 call sum_energy(remd_ene(0,i),.false.)
923 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
925 ene_i_iex=remd_ene(0,i)
927 ! write (iout,*) "0,iex",remd_t_bath(iex)
928 ! call enerprint(remd_ene(0,iex))
930 call sum_energy(remd_ene(0,iex),.false.)
931 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
932 write (iout,*) "ERROR: inconsistent energies:",iex,&
933 ene_iex_iex,remd_ene(0,iex)
935 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
936 ! write (iout,*) "i",i," iex",iex
937 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
938 ! & " ene_i_iex",ene_i_iex,
939 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
941 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
942 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
944 ! write(iout,*) 'delta',delta
945 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
946 ! & (remd_ene(i)-remd_ene(iex))/Rb/
947 ! & (remd_t_bath(i)*remd_t_bath(iex))
949 if (delta .gt. 50.0d0) then
955 else if (delta.lt.-50.0d0) then
964 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
965 xxx=ran_number(0.0d0,1.0d0)
966 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
968 if (delta .gt. xxx) then
970 remd_t_bath(i)=remd_t_bath(iex)
972 remd_ene(0,i)=ene_i_iex
973 remd_ene(0,iex)=ene_iex_i
979 ehighi(i)=ehighi(iex)
986 nupa(k,i)=nupa(k,iex)
989 ndowna(k,i)=ndowna(k,iex)
993 if (ifirst(il).eq.i) ifirst(il)=iex
995 if (nupa(k,il).eq.i) then
997 elseif (nupa(k,il).eq.iex) then
1002 if (ndowna(k,il).eq.i) then
1004 elseif (ndowna(k,il).eq.iex) then
1010 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1012 i2rep(i-1)=i2rep(iex-1)
1015 ! write(iout,*) 'exchange',i,iex
1016 ! write (iout,'(a8,100i4)') "@ ifirst",
1017 ! & (ifirst(k),k=1,remd_m(1))
1019 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1020 ! & (nupa(k,il),k=1,nupa(0,il))
1021 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1022 ! & (ndowna(k,il),k=1,ndowna(0,il))
1027 remd_ene(0,iex)=ene_iex_iex
1028 remd_ene(0,i)=ene_i_i
1034 !d write (iout,*) "exchange completed"
1038 !d write(iout,*) "########",ii
1040 i_temp=iran_num(1,nrep)
1041 i_mult=iran_num(1,remd_m(i_temp))
1042 i_iset=iran_num(1,nset)
1043 i_mset=iran_num(1,mset(i_iset))
1044 i=i_index(i_temp,i_mult,i_iset,i_mset)
1046 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1049 !d write(iout,*) "i_dir=",i_dir
1051 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1054 i_mult1=iran_num(1,remd_m(i_temp1))
1056 i_mset1=iran_num(1,mset(i_iset1))
1057 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1059 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1062 i_mult1=iran_num(1,remd_m(i_temp1))
1064 i_mset1=iran_num(1,mset(i_iset1))
1065 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1066 econstr_temp_i=remd_ene(20,i)
1067 econstr_temp_iex=remd_ene(20,iex)
1068 remd_ene(20,i)=remd_ene(n_ene+3,i)
1069 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1071 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1074 i_mult1=iran_num(1,remd_m(i_temp1))
1076 i_mset1=iran_num(1,mset(i_iset1))
1077 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1078 econstr_temp_i=remd_ene(20,i)
1079 econstr_temp_iex=remd_ene(20,iex)
1080 remd_ene(20,i)=remd_ene(n_ene+3,i)
1081 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1087 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1090 ! Swap temperatures between conformations i and iex with recalculating the free energies
1091 ! following temperature changes.
1092 ene_iex_iex=remd_ene(0,iex)
1093 ene_i_i=remd_ene(0,i)
1094 !o write (iout,*) "rescaling weights with temperature",
1096 call rescale_weights(remd_t_bath(i))
1098 call sum_energy(remd_ene(0,iex),.false.)
1099 ene_iex_i=remd_ene(0,iex)
1100 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1101 ! call sum_energy(remd_ene(0,i),.false.)
1102 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1103 ! write (iout,*) "rescaling weights with temperature",
1104 ! & remd_t_bath(iex)
1105 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1106 ! write (iout,*) "ERROR: inconsistent energies:",i,
1107 ! & ene_i_i,remd_ene(0,i)
1109 call rescale_weights(remd_t_bath(iex))
1110 call sum_energy(remd_ene(0,i),.false.)
1111 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1112 ene_i_iex=remd_ene(0,i)
1113 ! call sum_energy(remd_ene(0,iex),.false.)
1114 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1115 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1116 ! & ene_iex_iex,remd_ene(0,iex)
1118 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1119 ! write (iout,*) "i",i," iex",iex
1120 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1121 !d & " ene_i_iex",ene_i_iex,
1122 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1123 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1124 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1126 !d write(iout,*) 'delta',delta
1127 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1128 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1129 ! & (remd_t_bath(i)*remd_t_bath(iex))
1130 if (delta .gt. 50.0d0) then
1135 if (i_dir.eq.1.or.i_dir.eq.3) &
1136 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1137 if (i_dir.eq.2.or.i_dir.eq.3) &
1138 iremd_tot_usa(int(i2set(i-1)))= &
1139 iremd_tot_usa(int(i2set(i-1)))+1
1140 xxx=ran_number(0.0d0,1.0d0)
1141 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1142 if (delta .gt. xxx) then
1144 remd_t_bath(i)=remd_t_bath(iex)
1145 remd_t_bath(iex)=tmp
1148 iremd_iset(i)=iremd_iset(iex)
1149 iremd_iset(iex)=itmp
1151 remd_ene(0,i)=ene_i_iex
1152 remd_ene(0,iex)=ene_iex_i
1154 if (i_dir.eq.1.or.i_dir.eq.3) &
1155 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1158 i2rep(i-1)=i2rep(iex-1)
1161 if (i_dir.eq.2.or.i_dir.eq.3) &
1162 iremd_acc_usa(int(i2set(i-1)))= &
1163 iremd_acc_usa(int(i2set(i-1)))+1
1166 i2set(i-1)=i2set(iex-1)
1169 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1170 i_index(i_temp,i_mult,i_iset,i_mset)= &
1171 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1172 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1175 remd_ene(0,iex)=ene_iex_iex
1176 remd_ene(0,i)=ene_i_i
1177 remd_ene(20,iex)=econstr_temp_iex
1178 remd_ene(20,i)=econstr_temp_i
1182 !d do il1=1,mset(il)
1185 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1198 !-------------------------------------
1199 write (iout,*) "NREP",nrep
1201 if(iremd_tot(i).ne.0) &
1202 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1203 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1208 if(iremd_tot_usa(i).ne.0) &
1209 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1210 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1216 !d write (iout,'(a6,100i4)') "ifirst",
1217 !d & (ifirst(i),i=1,remd_m(1))
1219 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1220 !d & (nupa(i,il),i=1,nupa(0,il))
1221 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1222 !d & (ndowna(i,il),i=1,ndowna(0,il))
1227 !d write (iout,*) "Before scatter"
1230 if (me.eq.king) then
1231 write (iout,*) "t_bath before scatter",remd_t_bath
1235 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1236 t_bath,1,mpi_double_precision,king,&
1238 !d write (iout,*) "After scatter"
1241 call mpi_scatter(iremd_iset,1,mpi_integer,&
1242 iset,1,mpi_integer,king,&
1246 if (me.eq.king .or. .not. out1file) then
1247 write(iout,*) 'REMD scatter time=',time07-time06
1251 call mpi_scatter(elowi,1,mpi_double_precision,&
1252 elow,1,mpi_double_precision,king,&
1254 call mpi_scatter(ehighi,1,mpi_double_precision,&
1255 ehigh,1,mpi_double_precision,king,&
1258 call rescale_weights(t_bath)
1259 !o write (iout,*) "Processor",me,
1260 !o & " rescaling weights with temperature",t_bath
1262 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
1264 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
1268 !de write(iout,*) 'REMD after',me,t_bath
1270 if (me.eq.king .or. .not. out1file) then
1271 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1272 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1273 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1274 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1275 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1276 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1277 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1278 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1279 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1285 if (restart1file) then
1286 if (me.eq.king .or. .not. out1file) &
1287 write(iout,*) 'writing restart at the end of run'
1288 call write1rst(i_index)
1291 if (traj1file) call write1traj
1293 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1294 !deb & icache_all,1,mpi_integer,king,
1295 !deb & CG_COMM,ierr)
1296 !deb write(iout,'(a40,8000i8)')
1297 !deb & ' ntwx_cache after traj1file at the end',
1298 !deb & (icache_all(i),i=1,nodes)
1303 t_MD=MPI_Wtime()-tt0
1307 if (me.eq.king .or. .not. out1file) then
1308 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1310 'MD calculations setup:',t_MDsetup,&
1311 'Energy & gradient evaluation:',t_enegrad,&
1312 'Stochastic MD setup:',t_langsetup,&
1313 'Stochastic MD step setup:',t_sdsetup,&
1315 write (iout,'(/28(1h=),a25,27(1h=))') &
1316 ' End of MD calculation '
1318 !el common /przechowalnia/
1319 ! deallocate(d_restart1)
1320 ! deallocate(d_restart2)
1324 end subroutine MREMD
1325 !-----------------------------------------------------------------------------
1326 subroutine write1rst(i_index)
1329 ! implicit real*8 (a-h,o-z)
1330 ! include 'DIMENSIONS'
1332 ! include 'COMMON.MD'
1333 ! include 'COMMON.IOUNITS'
1334 ! include 'COMMON.REMD'
1335 ! include 'COMMON.SETUP'
1336 ! include 'COMMON.CHAIN'
1337 ! include 'COMMON.SBRIDGE'
1338 ! include 'COMMON.INTERACT'
1340 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1341 !el d_restart2(3,2*nres*maxprocs)
1342 real(kind=4) :: r_d(3,0:2*nres)
1343 real(kind=4) :: t5_restart1(5)
1344 integer :: iret,itmp
1345 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1346 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1348 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1349 !el common /przechowalnia/ d_restart1,d_restart2
1350 integer :: i,j,il,il1,ierr,ixdrf
1355 t5_restart1(4)=t_bath
1356 t5_restart1(5)=Uconst
1358 call mpi_gather(t5_restart1,5,mpi_real,&
1359 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1367 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1368 d_restart1,3*2*nres+3,mpi_real,king,&
1380 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1381 d_restart2,3*2*nres+3,mpi_real,king,&
1386 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1388 call xdrfint_(ixdrf, i2rep(i), iret)
1391 call xdrfint_(ixdrf, ifirst(i), iret)
1395 call xdrfint_(ixdrf, nupa(i,il), iret)
1399 call xdrfint_(ixdrf, ndowna(i,il), iret)
1405 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1412 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1419 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1425 call xdrfint_(ixdrf, nset, iret)
1427 call xdrfint_(ixdrf,mset(i), iret)
1430 call xdrfint_(ixdrf,i2set(i), iret)
1436 itmp=i_index(i,j,il,il1)
1437 call xdrfint_(ixdrf,itmp, iret)
1444 call xdrfclose_(ixdrf, iret)
1446 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1448 call xdrfint(ixdrf, i2rep(i), iret)
1451 call xdrfint(ixdrf, ifirst(i), iret)
1455 call xdrfint(ixdrf, nupa(i,il), iret)
1459 call xdrfint(ixdrf, ndowna(i,il), iret)
1465 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1472 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1479 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1486 call xdrfint(ixdrf, nset, iret)
1488 call xdrfint(ixdrf,mset(i), iret)
1491 call xdrfint(ixdrf,i2set(i), iret)
1497 itmp=i_index(i,j,il,il1)
1498 call xdrfint(ixdrf,itmp, iret)
1505 call xdrfclose(ixdrf, iret)
1509 end subroutine write1rst
1510 !-----------------------------------------------------------------------------
1511 subroutine write1traj
1513 ! implicit real*8 (a-h,o-z)
1514 ! include 'DIMENSIONS'
1516 ! include 'COMMON.MD'
1517 ! include 'COMMON.IOUNITS'
1518 ! include 'COMMON.REMD'
1519 ! include 'COMMON.SETUP'
1520 ! include 'COMMON.CHAIN'
1521 ! include 'COMMON.SBRIDGE'
1522 ! include 'COMMON.INTERACT'
1524 real(kind=4) :: t5_restart1(5)
1525 integer :: iret,itmp
1526 real(kind=4) :: xcoord(3,2*nres+2),prec
1527 real(kind=4) :: r_qfrag(50),r_qpair(100)
1528 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1529 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1530 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1531 p_uscdiff(100*Nprocs) !(100*maxprocs)
1532 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1533 real(kind=4) :: r_c(3,2*nres+2)
1534 !el common /przechowalnia/ p_c
1536 integer :: i,j,il,ierr,ii,ixdrf
1538 call mpi_bcast(ii_write,1,mpi_integer,&
1542 print *,'traj1file',me,ii_write,ntwx_cache
1546 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1548 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1551 ! write (iout,*) "before gather write1traj: from node",ii
1553 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1555 t5_restart1(1)=totT_cache(ii)
1556 t5_restart1(2)=EK_cache(ii)
1557 t5_restart1(3)=potE_cache(ii)
1558 t5_restart1(4)=t_bath_cache(ii)
1559 t5_restart1(5)=Uconst_cache(ii)
1560 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1562 call mpi_gather(t5_restart1,5,mpi_real,&
1563 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1565 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1568 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1569 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1572 r_qfrag(i)=qfrag_cache(i,ii)
1575 r_qpair(i)=qpair_cache(i,ii)
1578 r_utheta(i)=utheta_cache(i,ii)
1579 r_ugamma(i)=ugamma_cache(i,ii)
1580 r_uscdiff(i)=uscdiff_cache(i,ii)
1583 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1584 p_qfrag,nfrag,mpi_real,king,&
1586 call mpi_gather(r_qpair,npair,mpi_real,&
1587 p_qpair,npair,mpi_real,king,&
1589 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1590 p_utheta,nfrag_back,mpi_real,king,&
1592 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1593 p_ugamma,nfrag_back,mpi_real,king,&
1595 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1596 p_uscdiff,nfrag_back,mpi_real,king,&
1600 write (iout,*) "p_qfrag"
1602 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1604 write (iout,*) "p_qpair"
1606 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1612 r_c(j,i)=c_cache(j,i,ii)
1616 call mpi_gather(r_c,3*2*nres,mpi_real,&
1617 p_c,3*2*nres,mpi_real,king,&
1623 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1624 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1625 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1626 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1627 call xdrfint_(ixdrf, nss, iret)
1630 call xdrfint(ixdrf, idssb(j)+nres, iret)
1631 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1633 call xdrfint_(ixdrf, ihpb(j), iret)
1634 call xdrfint_(ixdrf, jhpb(j), iret)
1637 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1638 call xdrfint_(ixdrf, iset_restart1(il), iret)
1640 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1643 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1646 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1647 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1648 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1653 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1658 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1662 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1666 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1667 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1668 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1669 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1670 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1671 call xdrfint(ixdrf, nss, iret)
1674 call xdrfint(ixdrf, idssb(j)+nres, iret)
1675 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1677 call xdrfint(ixdrf, ihpb(j), iret)
1678 call xdrfint(ixdrf, jhpb(j), iret)
1681 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1682 call xdrfint(ixdrf, iset_restart1(il), iret)
1684 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1687 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1690 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1691 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1692 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1697 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1702 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1706 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1712 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1714 if(me.eq.king) call xdrfclose(ixdrf, iret)
1716 do i=1,ntwx_cache-ii_write
1718 totT_cache(i)=totT_cache(ii_write+i)
1719 EK_cache(i)=EK_cache(ii_write+i)
1720 potE_cache(i)=potE_cache(ii_write+i)
1721 t_bath_cache(i)=t_bath_cache(ii_write+i)
1722 Uconst_cache(i)=Uconst_cache(ii_write+i)
1723 iset_cache(i)=iset_cache(ii_write+i)
1726 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1729 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1732 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1733 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1734 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1739 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1743 ntwx_cache=ntwx_cache-ii_write
1745 end subroutine write1traj
1746 !-----------------------------------------------------------------------------
1747 subroutine read1restart(i_index)
1749 ! implicit real*8 (a-h,o-z)
1750 ! include 'DIMENSIONS'
1752 ! include 'COMMON.MD'
1753 ! include 'COMMON.IOUNITS'
1754 ! include 'COMMON.REMD'
1755 ! include 'COMMON.SETUP'
1756 ! include 'COMMON.CHAIN'
1757 ! include 'COMMON.SBRIDGE'
1758 ! include 'COMMON.INTERACT'
1759 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1760 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1761 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1762 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1764 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1765 !el common /przechowalnia/ d_restart1
1766 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1768 write (*,*) "Processor",me," called read1restart"
1771 open(irest2,file=mremd_rst_name,status='unknown')
1772 read(irest2,*,err=334) i
1773 write(iout,*) "Reading old rst in ASCI format"
1775 call read1restart_old
1779 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1782 call xdrfint_(ixdrf, i2rep(i), iret)
1785 call xdrfint_(ixdrf, ifirst(i), iret)
1788 call xdrfint_(ixdrf, nupa(0,il), iret)
1790 call xdrfint_(ixdrf, nupa(i,il), iret)
1793 call xdrfint_(ixdrf, ndowna(0,il), iret)
1795 call xdrfint_(ixdrf, ndowna(i,il), iret)
1800 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1804 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1807 call xdrfint(ixdrf, i2rep(i), iret)
1810 call xdrfint(ixdrf, ifirst(i), iret)
1813 call xdrfint(ixdrf, nupa(0,il), iret)
1815 call xdrfint(ixdrf, nupa(i,il), iret)
1818 call xdrfint(ixdrf, ndowna(0,il), iret)
1820 call xdrfint(ixdrf, ndowna(i,il), iret)
1825 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1830 call mpi_scatter(t_restart1,5,mpi_real,&
1831 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1835 t_bath=t5_restart1(4)
1840 ! read(irest2,'(3e15.5)')
1841 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1844 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1846 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1852 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1853 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1863 ! read(irest2,'(3e15.5)')
1864 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1867 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1869 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1875 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1876 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1887 call xdrfint_(ixdrf, nset, iret)
1889 call xdrfint_(ixdrf,mset(i), iret)
1892 call xdrfint_(ixdrf,i2set(i), iret)
1898 call xdrfint_(ixdrf,itmp, iret)
1899 i_index(i,j,il,il1)=itmp
1907 call xdrfint(ixdrf, nset, iret)
1909 call xdrfint(ixdrf,mset(i), iret)
1912 call xdrfint(ixdrf,i2set(i), iret)
1918 call xdrfint(ixdrf,itmp, iret)
1919 i_index(i,j,il,il1)=itmp
1926 call mpi_scatter(i2set,1,mpi_integer,&
1927 iset,1,mpi_integer,king,&
1932 if(me.eq.king) close(irest2)
1934 end subroutine read1restart
1935 !-----------------------------------------------------------------------------
1936 subroutine read1restart_old
1938 ! implicit real*8 (a-h,o-z)
1939 ! include 'DIMENSIONS'
1941 ! include 'COMMON.MD'
1942 ! include 'COMMON.IOUNITS'
1943 ! include 'COMMON.REMD'
1944 ! include 'COMMON.SETUP'
1945 ! include 'COMMON.CHAIN'
1946 ! include 'COMMON.SBRIDGE'
1947 ! include 'COMMON.INTERACT'
1948 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1949 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1950 !el common /przechowalnia/ d_restart1
1952 integer :: i,j,il,ierr
1955 open(irest2,file=mremd_rst_name,status='unknown')
1956 read (irest2,*) (i2rep(i),i=0,nodes-1)
1957 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1959 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1960 read (irest2,*) ndowna(0,il),&
1961 (ndowna(i,il),i=1,ndowna(0,il))
1964 read(irest2,*) (t_restart1(j,il),j=1,4)
1967 call mpi_scatter(t_restart1,5,mpi_real,&
1968 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1972 t_bath=t5_restart1(4)
1977 read(irest2,'(3e15.5)') &
1978 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
1982 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1983 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1993 read(irest2,'(3e15.5)') &
1994 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
1998 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1999 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2005 if(me.eq.king) close(irest2)
2007 end subroutine read1restart_old
2008 !----------------------------------------------------------------
2009 subroutine alloc_MREMD_arrays
2011 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2012 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
2014 ! common /remdcommon/ in io: read_REMDpar
2015 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2016 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2017 ! common /remdrestart/
2018 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2020 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2021 allocate(ifirst(0:nodes)) !(maxprocs)
2022 allocate(nupa(0:nodes,0:2*nodes))
2023 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2024 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2025 allocate(iset_restart1(nodes)) !(maxprocs)
2026 ! common /traj1cache/
2027 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2028 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2029 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2030 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2031 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2032 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2033 allocate(utheta_cache(nfrag_back,max_cache_traj))
2034 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2035 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2036 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2039 end subroutine alloc_MREMD_arrays
2040 !-----------------------------------------------------------------------------
2041 !-----------------------------------------------------------------------------