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(maxprocs)
86 integer :: iremd_iset(maxprocs)
87 integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
88 real(kind=8) :: remd_ene(0:n_ene+4,maxprocs)
89 integer :: iremd_acc(maxprocs),iremd_tot(maxprocs)
90 integer :: iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
91 integer :: rstcount !el ilen,
93 character(len=50) :: tytul
96 !old integer nup(0:maxprocs),ndown(0:maxprocs)
97 integer :: rep2i(0:maxprocs),ireqi(maxprocs)
98 integer :: icache_all(maxprocs)
99 integer :: status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
100 logical :: synflag, end_of_run, file_exist = .false.!, ovrtim
102 real(kind=8) :: delta,time00,time01,time001,time02,time03,time04,&
103 time05,time06,time07,time08,tt0,scalfac,ene_iex_iex,&
104 ene_i_i,ene_iex_i,ene_i_iex,xxx,tmp,econstr_temp_i,&
106 integer :: k,il,il1,i,j,nharp,ii,ierr,itime_master,irr,iex,&
107 i_set_temp,itmp,i_temp,i_mult,i_iset,i_mset,i_dir,i_temp1,&
108 i_mult1,i_iset1,i_mset1,ierror
109 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
110 !deb imin_itime_old=0
118 if(me.eq.king.or..not.out1file) then
119 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
120 write (iout,*) "NREP=",nrep
124 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
125 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
127 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
129 !d print *,'MREMD',nodes
130 !d print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
131 !de write (iout,*) "Start MREMD: me",me," t_bath",t_bath
135 do il1=1,max0(mset(il),1)
145 i_index(i,j,il,il1)=k
151 if(me.eq.king.or..not.out1file) then
152 write(iout,*) (i2rep(i),i=0,nodes-1)
153 write(iout,*) (i2set(i),i=0,nodes-1)
158 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
165 ! print *,'i2rep',me,i2rep(me)
166 ! print *,'rep2i',(rep2i(i),i=0,nrep)
168 !old if (i2rep(me).eq.nrep) then
171 !old nup(0)=remd_m(i2rep(me)+1)
172 !old k=rep2i(int(i2rep(me)))+1
179 !d print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
181 !old if (i2rep(me).eq.1) then
184 !old ndown(0)=remd_m(i2rep(me)-1)
185 !old k=rep2i(i2rep(me)-2)+1
192 !d print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
194 !el common /przechowalnia/
195 if(.not.allocated(d_restart1)) allocate(d_restart1(3,nres2*nodes))
196 if(.not.allocated(d_restart2)) allocate(d_restart2(3,nres2*nodes))
197 if(.not.allocated(p_c)) allocate(p_c(3,(nres2+2)*nodes))
200 write (*,*) "Processor",me," rest",rest,&
201 "restart1fie",restart1file
202 if(rest.and.restart1file) then
204 inquire(file=mremd_rst_name,exist=file_exist)
205 !d write (*,*) me," Before broadcast: file_exist",file_exist
206 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
208 !d write (*,*) me," After broadcast: file_exist",file_exist
210 if(me.eq.king.or..not.out1file) &
211 write (iout,*) 'Master is reading restart1file'
212 call read1restart(i_index)
214 if(me.eq.king.or..not.out1file) &
215 write (iout,*) 'WARNING : no restart1file'
218 if(me.eq.king.or..not.out1file) then
219 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
220 write(iout,*) "i_index"
225 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
234 if (rest.and..not.restart1file) &
235 inquire(file=mremd_rst_name,exist=file_exist)
236 if(.not.file_exist.and.rest.and..not.restart1file) &
237 write(iout,*) 'WARNING : no restart file',mremd_rst_name
238 IF (rest.and.file_exist.and..not.restart1file) THEN
239 write (iout,*) 'Master is reading restart file',&
241 open(irest2,file=mremd_rst_name,status='unknown')
243 read (irest2,*) (i2rep(i),i=0,nodes-1)
245 read (irest2,*) (ifirst(i),i=1,remd_m(1))
248 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
250 read (irest2,*) ndowna(0,il),&
251 (ndowna(i,il),i=1,ndowna(0,il))
257 read (irest2,*) (mset(i),i=1,nset)
259 read (irest2,*) (i2set(i),i=0,nodes-1)
264 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
269 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
270 write(iout,*) "i_index"
275 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
284 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
285 write (iout,'(a6,1000i5)') "ifirst",&
286 (ifirst(i),i=1,remd_m(1))
288 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
289 (nupa(i,il),i=1,nupa(0,il))
290 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
291 (ndowna(i,il),i=1,ndowna(0,il))
293 ELSE IF (.not.(rest.and.file_exist)) THEN
299 if (i2rep(il-1).eq.nrep) then
302 nupa(0,il)=remd_m(i2rep(il-1)+1)
303 k=rep2i(int(i2rep(il-1)))+1
309 if (i2rep(il-1).eq.1) then
312 ndowna(0,il)=remd_m(i2rep(il-1)-1)
313 k=rep2i(i2rep(il-1)-2)+1
321 write (iout,'(a6,100i4)') "ifirst",&
322 (ifirst(i),i=1,remd_m(1))
324 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
325 (nupa(i,il),i=1,nupa(0,il))
326 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
327 (ndowna(i,il),i=1,ndowna(0,il))
333 ! t_bath=retmin+(retmax-retmin)*me/(nodes-1)
334 if(.not.(rest.and.file_exist.and.restart1file)) then
335 if (me .eq. king) then
338 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
340 !d print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
341 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
346 if(me.eq.king.or..not.out1file) &
347 write(iout,*) me,"iset=",iset,"t_bath=",t_bath
350 stdfp=dsqrt(2*Rb*t_bath/d_time)
352 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
355 ! print *,'irep',me,t_bath
357 if (me.eq.king .or. .not. out1file) &
358 write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
359 call rescale_weights(t_bath)
363 !------copy MD--------------
364 ! The driver for molecular dynamics subroutines
365 !------------------------------------------------
371 if(me.eq.king.or..not.out1file) &
372 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
378 ! Determine the inverse of the inertia matrix.
379 call setup_MD_matrices
383 if (me.eq.king .or. .not. out1file) &
384 write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
385 stdfp=dsqrt(2*Rb*t_bath/d_time)
387 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
389 call rescale_weights(t_bath)
393 t_MDsetup = MPI_Wtime()-tt0
395 t_MDsetup = tcpu()-tt0
398 ! Entering the MD loop
404 if (lang.eq.2 .or. lang.eq.3) then
408 call sd_verlet_p_setup
410 call sd_verlet_ciccotti_setup
414 pfric0_mat(i,j,0)=pfric_mat(i,j)
415 afric0_mat(i,j,0)=afric_mat(i,j)
416 vfric0_mat(i,j,0)=vfric_mat(i,j)
417 prand0_mat(i,j,0)=prand_mat(i,j)
418 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
419 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
424 flag_stoch(i)=.false.
428 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
430 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
434 else if (lang.eq.1 .or. lang.eq.4) then
438 if (me.eq.king .or. .not. out1file) &
439 write(iout,*) 'Setup time',time00-walltime
442 t_langsetup=MPI_Wtime()-tt0
445 t_langsetup=tcpu()-tt0
451 do while(.not.end_of_run)
453 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
454 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
456 if (lang.gt.0 .and. surfarea .and. &
457 mod(itime,reset_fricmat).eq.0) then
458 if (lang.eq.2 .or. lang.eq.3) then
462 call sd_verlet_p_setup
464 call sd_verlet_ciccotti_setup
468 pfric0_mat(i,j,0)=pfric_mat(i,j)
469 afric0_mat(i,j,0)=afric_mat(i,j)
470 vfric0_mat(i,j,0)=vfric_mat(i,j)
471 prand0_mat(i,j,0)=prand_mat(i,j)
472 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
473 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
478 flag_stoch(i)=.false.
481 else if (lang.eq.1 .or. lang.eq.4) then
484 write (iout,'(a,i10)') &
485 "Friction matrix reset based on surface area, itime",itime
487 if (reset_vel .and. tbf .and. lang.eq.0 &
488 .and. mod(itime,count_reset_vel).eq.0) then
490 if (me.eq.king .or. .not. out1file) &
491 write(iout,'(a,f20.2)') &
492 "Velocities reset to random values, time",totT
495 d_t_old(j,i)=d_t(j,i)
499 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
503 d_t(j,0)=d_t(j,0)-vcm(j)
506 kinetic_T=2.0d0/(dimen3*Rb)*EK
507 scalfac=dsqrt(T_bath/kinetic_T)
508 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
511 d_t_old(j,i)=scalfac*d_t(j,i)
517 ! Time-reversible RESPA algorithm
518 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
519 call RESPA_step(itime)
521 ! Variable time step algorithm.
522 call velverlet_step(itime)
526 call brown_step(itime)
528 print *,"Brown dynamics not here!"
530 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
536 if (mod(itime,ntwe).eq.0) call statout(itime)
538 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
539 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
541 call hairpin(.true.,nharp,iharp)
542 call secondary2(.true.)
543 call pdbout(potE,tytul,ipdb)
548 if (mod(itime,ntwx).eq.0.and.traj1file) then
549 if(ntwx_cache.lt.max_cache_traj_use) then
550 ntwx_cache=ntwx_cache+1
552 if (max_cache_traj_use.ne.1) &
553 print *,itime,"processor ",me," over cache ",ntwx_cache
556 totT_cache(i)=totT_cache(i+1)
557 EK_cache(i)=EK_cache(i+1)
558 potE_cache(i)=potE_cache(i+1)
559 t_bath_cache(i)=t_bath_cache(i+1)
560 Uconst_cache(i)=Uconst_cache(i+1)
561 iset_cache(i)=iset_cache(i+1)
564 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
567 qpair_cache(ii,i)=qpair_cache(ii,i+1)
570 utheta_cache(ii,i)=utheta_cache(ii,i+1)
571 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
572 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
578 c_cache(j,ii,i)=c_cache(j,ii,i+1)
584 totT_cache(ntwx_cache)=totT
585 EK_cache(ntwx_cache)=EK
586 potE_cache(ntwx_cache)=potE
587 t_bath_cache(ntwx_cache)=t_bath
588 Uconst_cache(ntwx_cache)=Uconst
589 iset_cache(ntwx_cache)=iset
592 qfrag_cache(i,ntwx_cache)=qfrag(i)
595 qpair_cache(i,ntwx_cache)=qpair(i)
598 utheta_cache(i,ntwx_cache)=utheta(i)
599 ugamma_cache(i,ntwx_cache)=ugamma(i)
600 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
605 c_cache(j,i,ntwx_cache)=c(j,i)
610 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
611 .and..not.restart1file) then
614 open(irest1,file=mremd_rst_name,status='unknown')
615 write (irest1,*) "i2rep"
616 write (irest1,*) (i2rep(i),i=0,nodes-1)
617 write (irest1,*) "ifirst"
618 write (irest1,*) (ifirst(i),i=1,remd_m(1))
620 write (irest1,*) "nupa",il
621 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
622 write (irest1,*) "ndowna",il
623 write (irest1,*) ndowna(0,il),&
624 (ndowna(i,il),i=1,ndowna(0,il))
627 write (irest1,*) "nset"
628 write (irest1,*) nset
629 write (irest1,*) "mset"
630 write (irest1,*) (mset(i),i=1,nset)
631 write (irest1,*) "i2set"
632 write (irest1,*) (i2set(i),i=0,nodes-1)
633 write (irest1,*) "i_index"
637 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
645 open(irest2,file=rest2name,status='unknown')
646 write(irest2,*) totT,EK,potE,totE,t_bath
648 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
651 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
654 write (irest2,*) iset
661 ! forced synchronization
662 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
663 .and. .not. mremdsync) then
665 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
667 call mpi_recv(itime_master, 1, MPI_INTEGER,&
668 0,101,CG_COMM, status, ierr)
669 call mpi_barrier(CG_COMM, ierr)
670 !deb if (out1file.or.traj1file) then
671 !deb call mpi_gather(itime,1,mpi_integer,
672 !deb & icache_all,1,mpi_integer,king,
675 call mpi_gather(ntwx_cache,1,mpi_integer,&
676 icache_all,1,mpi_integer,king,&
679 write(iout,*) 'REMD synchro at',itime_master,itime
680 if (itime_master.ge.n_timestep .or. ovrtim()) &
682 !time call flush(iout)
687 if ((mod(itime,nstex).eq.0.and.me.eq.king &
688 .or.end_of_run.and.me.eq.king ) &
689 .and. .not. mremdsync ) then
692 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
693 CG_COMM, ireqi(i), ierr)
694 !d write(iout,*) 'REMD synchro with',i
697 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
698 call mpi_barrier(CG_COMM, ierr)
700 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
701 if (out1file.or.traj1file) then
702 !deb call mpi_gather(itime,1,mpi_integer,
703 !deb & itime_all,1,mpi_integer,king,
705 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
706 !deb & (itime_all(i),i=1,nodes)
708 !deb imin_itime=itime_all(1)
710 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
712 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
713 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
714 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
715 call mpi_gather(ntwx_cache,1,mpi_integer,&
716 icache_all,1,mpi_integer,king,&
718 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
719 ! & (icache_all(i),i=1,nodes)
720 ii_write=icache_all(1)
722 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
724 ! write(iout,*) "MIN ii_write=",ii_write
727 !time call flush(iout)
729 if(mremdsync .and. mod(itime,nstex).eq.0) then
731 if (me.eq.king .or. .not. out1file) &
732 write(iout,*) 'REMD synchro at',itime
735 call mpi_gather(ntwx_cache,1,mpi_integer,&
736 icache_all,1,mpi_integer,king,&
739 write(iout,'(a19,8000i8)') ' ntwx_cache',&
740 (icache_all(i),i=1,nodes)
741 ii_write=icache_all(1)
743 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
745 write(iout,*) "MIN ii_write=",ii_write
751 ! Update the time safety limiy
752 if (time001-time00.gt.safety) then
753 safety=time001-time00+600
754 write (iout,*) "****** SAFETY increased to",safety," s"
756 if (ovrtim()) end_of_run=.true.
758 if(synflag.and..not.end_of_run) then
762 write(iout,*) 'REMD before',me,t_bath
764 ! call mpi_gather(t_bath,1,mpi_double_precision,
765 ! & remd_t_bath,1,mpi_double_precision,king,
767 potEcomp(n_ene+1)=t_bath
769 potEcomp(n_ene+2)=iset
770 if (iset.lt.nset) then
774 potEcomp(n_ene+3)=Uconst
781 potEcomp(n_ene+4)=Uconst
785 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
786 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
789 call mpi_gather(elow,1,mpi_double_precision,&
790 elowi,1,mpi_double_precision,king,&
792 call mpi_gather(ehigh,1,mpi_double_precision,&
793 ehighi,1,mpi_double_precision,king,&
798 if (me.eq.king .or. .not. out1file) then
799 write(iout,*) 'REMD gather times=',time03-time01 &
803 if (restart1file) call write1rst(i_index)
806 if (me.eq.king .or. .not. out1file) then
807 write(iout,*) 'REMD writing rst time=',time04-time03
810 if (traj1file) call write1traj
812 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
813 !deb & icache_all,1,mpi_integer,king,
815 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
816 !deb & (icache_all(i),i=1,nodes)
821 if (me.eq.king .or. .not. out1file) then
822 write(iout,*) 'REMD writing traj time=',time05-time04
829 remd_t_bath(i)=remd_ene(n_ene+1,i)
830 iremd_iset(i)=remd_ene(n_ene+2,i)
834 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
836 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
840 write(iout,*) 'REMD exchange temp,ene'
842 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
843 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
847 !-------------------------------------
849 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
852 write (iout,*) "remd_m(1)",remd_m(1)
854 i=ifirst(iran_num(1,remd_m(1)))
861 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
863 if(i.gt.0.and.nupa(0,i).gt.0) then
865 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
867 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
869 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
871 ! do while (iex.eq.i)
872 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
873 iex=nupa(iran_num(1,int(nupa(0,i))),i)
875 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
877 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
879 ! Swap temperatures between conformations i and iex with recalculating the free energies
880 ! following temperature changes.
881 ene_iex_iex=remd_ene(0,iex)
882 ene_i_i=remd_ene(0,i)
883 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
884 ! & " iex",iex," ene_iex_iex",ene_iex_iex
885 ! write (iout,*) "rescaling weights with temperature",
888 call rescale_weights(remd_t_bath(i))
890 ! write (iout,*) "0,iex",remd_t_bath(i)
891 ! call enerprint(remd_ene(0,iex))
893 call sum_energy(remd_ene(0,iex),.false.)
894 ene_iex_i=remd_ene(0,iex)
895 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
897 ! write (iout,*) "0,i",remd_t_bath(i)
898 ! call enerprint(remd_ene(0,i))
900 call sum_energy(remd_ene(0,i),.false.)
901 ! write (iout,*) "ene_i_i",remd_ene(0,i)
903 ! write (iout,*) "rescaling weights with temperature",
905 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
906 write (iout,*) "ERROR: inconsistent energies:",i,&
907 ene_i_i,remd_ene(0,i)
909 call rescale_weights(remd_t_bath(iex))
911 ! write (iout,*) "0,i",remd_t_bath(iex)
912 ! call enerprint(remd_ene(0,i))
914 call sum_energy(remd_ene(0,i),.false.)
915 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
917 ene_i_iex=remd_ene(0,i)
919 ! write (iout,*) "0,iex",remd_t_bath(iex)
920 ! call enerprint(remd_ene(0,iex))
922 call sum_energy(remd_ene(0,iex),.false.)
923 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
924 write (iout,*) "ERROR: inconsistent energies:",iex,&
925 ene_iex_iex,remd_ene(0,iex)
927 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
928 ! write (iout,*) "i",i," iex",iex
929 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
930 ! & " ene_i_iex",ene_i_iex,
931 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
933 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
934 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
936 ! write(iout,*) 'delta',delta
937 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
938 ! & (remd_ene(i)-remd_ene(iex))/Rb/
939 ! & (remd_t_bath(i)*remd_t_bath(iex))
941 if (delta .gt. 50.0d0) then
947 else if (delta.lt.-50.0d0) then
956 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
957 xxx=ran_number(0.0d0,1.0d0)
958 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
960 if (delta .gt. xxx) then
962 remd_t_bath(i)=remd_t_bath(iex)
964 remd_ene(0,i)=ene_i_iex
965 remd_ene(0,iex)=ene_iex_i
971 ehighi(i)=ehighi(iex)
978 nupa(k,i)=nupa(k,iex)
981 ndowna(k,i)=ndowna(k,iex)
985 if (ifirst(il).eq.i) ifirst(il)=iex
987 if (nupa(k,il).eq.i) then
989 elseif (nupa(k,il).eq.iex) then
994 if (ndowna(k,il).eq.i) then
996 elseif (ndowna(k,il).eq.iex) then
1002 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1004 i2rep(i-1)=i2rep(iex-1)
1007 ! write(iout,*) 'exchange',i,iex
1008 ! write (iout,'(a8,100i4)') "@ ifirst",
1009 ! & (ifirst(k),k=1,remd_m(1))
1011 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1012 ! & (nupa(k,il),k=1,nupa(0,il))
1013 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1014 ! & (ndowna(k,il),k=1,ndowna(0,il))
1019 remd_ene(0,iex)=ene_iex_iex
1020 remd_ene(0,i)=ene_i_i
1026 !d write (iout,*) "exchange completed"
1030 !d write(iout,*) "########",ii
1032 i_temp=iran_num(1,nrep)
1033 i_mult=iran_num(1,remd_m(i_temp))
1034 i_iset=iran_num(1,nset)
1035 i_mset=iran_num(1,mset(i_iset))
1036 i=i_index(i_temp,i_mult,i_iset,i_mset)
1038 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1041 !d write(iout,*) "i_dir=",i_dir
1043 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1046 i_mult1=iran_num(1,remd_m(i_temp1))
1048 i_mset1=iran_num(1,mset(i_iset1))
1049 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1051 elseif(i_dir.eq.2 .and. mset(i_iset+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)
1058 econstr_temp_i=remd_ene(20,i)
1059 econstr_temp_iex=remd_ene(20,iex)
1060 remd_ene(20,i)=remd_ene(n_ene+3,i)
1061 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1063 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1066 i_mult1=iran_num(1,remd_m(i_temp1))
1068 i_mset1=iran_num(1,mset(i_iset1))
1069 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1070 econstr_temp_i=remd_ene(20,i)
1071 econstr_temp_iex=remd_ene(20,iex)
1072 remd_ene(20,i)=remd_ene(n_ene+3,i)
1073 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1079 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1082 ! Swap temperatures between conformations i and iex with recalculating the free energies
1083 ! following temperature changes.
1084 ene_iex_iex=remd_ene(0,iex)
1085 ene_i_i=remd_ene(0,i)
1086 !o write (iout,*) "rescaling weights with temperature",
1088 call rescale_weights(remd_t_bath(i))
1090 call sum_energy(remd_ene(0,iex),.false.)
1091 ene_iex_i=remd_ene(0,iex)
1092 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1093 ! call sum_energy(remd_ene(0,i),.false.)
1094 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1095 ! write (iout,*) "rescaling weights with temperature",
1096 ! & remd_t_bath(iex)
1097 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1098 ! write (iout,*) "ERROR: inconsistent energies:",i,
1099 ! & ene_i_i,remd_ene(0,i)
1101 call rescale_weights(remd_t_bath(iex))
1102 call sum_energy(remd_ene(0,i),.false.)
1103 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1104 ene_i_iex=remd_ene(0,i)
1105 ! call sum_energy(remd_ene(0,iex),.false.)
1106 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1107 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1108 ! & ene_iex_iex,remd_ene(0,iex)
1110 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1111 ! write (iout,*) "i",i," iex",iex
1112 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1113 !d & " ene_i_iex",ene_i_iex,
1114 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1115 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1116 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1118 !d write(iout,*) 'delta',delta
1119 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1120 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1121 ! & (remd_t_bath(i)*remd_t_bath(iex))
1122 if (delta .gt. 50.0d0) then
1127 if (i_dir.eq.1.or.i_dir.eq.3) &
1128 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1129 if (i_dir.eq.2.or.i_dir.eq.3) &
1130 iremd_tot_usa(int(i2set(i-1)))= &
1131 iremd_tot_usa(int(i2set(i-1)))+1
1132 xxx=ran_number(0.0d0,1.0d0)
1133 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1134 if (delta .gt. xxx) then
1136 remd_t_bath(i)=remd_t_bath(iex)
1137 remd_t_bath(iex)=tmp
1140 iremd_iset(i)=iremd_iset(iex)
1141 iremd_iset(iex)=itmp
1143 remd_ene(0,i)=ene_i_iex
1144 remd_ene(0,iex)=ene_iex_i
1146 if (i_dir.eq.1.or.i_dir.eq.3) &
1147 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1150 i2rep(i-1)=i2rep(iex-1)
1153 if (i_dir.eq.2.or.i_dir.eq.3) &
1154 iremd_acc_usa(int(i2set(i-1)))= &
1155 iremd_acc_usa(int(i2set(i-1)))+1
1158 i2set(i-1)=i2set(iex-1)
1161 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1162 i_index(i_temp,i_mult,i_iset,i_mset)= &
1163 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1164 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1167 remd_ene(0,iex)=ene_iex_iex
1168 remd_ene(0,i)=ene_i_i
1169 remd_ene(20,iex)=econstr_temp_iex
1170 remd_ene(20,i)=econstr_temp_i
1174 !d do il1=1,mset(il)
1177 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1190 !-------------------------------------
1191 write (iout,*) "NREP",nrep
1193 if(iremd_tot(i).ne.0) &
1194 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1195 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1200 if(iremd_tot_usa(i).ne.0) &
1201 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1202 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1208 !d write (iout,'(a6,100i4)') "ifirst",
1209 !d & (ifirst(i),i=1,remd_m(1))
1211 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1212 !d & (nupa(i,il),i=1,nupa(0,il))
1213 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1214 !d & (ndowna(i,il),i=1,ndowna(0,il))
1219 !d write (iout,*) "Before scatter"
1222 if (me.eq.king) then
1223 write (iout,*) "t_bath before scatter",remd_t_bath
1227 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1228 t_bath,1,mpi_double_precision,king,&
1230 !d write (iout,*) "After scatter"
1233 call mpi_scatter(iremd_iset,1,mpi_integer,&
1234 iset,1,mpi_integer,king,&
1238 if (me.eq.king .or. .not. out1file) then
1239 write(iout,*) 'REMD scatter time=',time07-time06
1243 call mpi_scatter(elowi,1,mpi_double_precision,&
1244 elow,1,mpi_double_precision,king,&
1246 call mpi_scatter(ehighi,1,mpi_double_precision,&
1247 ehigh,1,mpi_double_precision,king,&
1250 call rescale_weights(t_bath)
1251 !o write (iout,*) "Processor",me,
1252 !o & " rescaling weights with temperature",t_bath
1254 stdfp=dsqrt(2*Rb*t_bath/d_time)
1256 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1259 !de write(iout,*) 'REMD after',me,t_bath
1261 if (me.eq.king .or. .not. out1file) then
1262 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1263 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1264 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1265 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1266 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1267 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1268 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1269 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1270 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1276 if (restart1file) then
1277 if (me.eq.king .or. .not. out1file) &
1278 write(iout,*) 'writing restart at the end of run'
1279 call write1rst(i_index)
1282 if (traj1file) call write1traj
1284 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1285 !deb & icache_all,1,mpi_integer,king,
1286 !deb & CG_COMM,ierr)
1287 !deb write(iout,'(a40,8000i8)')
1288 !deb & ' ntwx_cache after traj1file at the end',
1289 !deb & (icache_all(i),i=1,nodes)
1294 t_MD=MPI_Wtime()-tt0
1298 if (me.eq.king .or. .not. out1file) then
1299 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1301 'MD calculations setup:',t_MDsetup,&
1302 'Energy & gradient evaluation:',t_enegrad,&
1303 'Stochastic MD setup:',t_langsetup,&
1304 'Stochastic MD step setup:',t_sdsetup,&
1306 write (iout,'(/28(1h=),a25,27(1h=))') &
1307 ' End of MD calculation '
1309 !el common /przechowalnia/
1310 ! deallocate(d_restart1)
1311 ! deallocate(d_restart2)
1315 end subroutine MREMD
1316 !-----------------------------------------------------------------------------
1317 subroutine write1rst(i_index)
1320 ! implicit real*8 (a-h,o-z)
1321 ! include 'DIMENSIONS'
1323 ! include 'COMMON.MD'
1324 ! include 'COMMON.IOUNITS'
1325 ! include 'COMMON.REMD'
1326 ! include 'COMMON.SETUP'
1327 ! include 'COMMON.CHAIN'
1328 ! include 'COMMON.SBRIDGE'
1329 ! include 'COMMON.INTERACT'
1331 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1332 !el d_restart2(3,2*nres*maxprocs)
1333 real(kind=4) :: r_d(3,2*nres)
1334 real(kind=4) :: t5_restart1(5)
1335 integer :: iret,itmp
1336 integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1337 !el common /przechowalnia/ d_restart1,d_restart2
1338 integer :: i,j,il,il1,ierr,ixdrf
1343 t5_restart1(4)=t_bath
1344 t5_restart1(5)=Uconst
1346 call mpi_gather(t5_restart1,5,mpi_real,&
1347 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1355 call mpi_gather(r_d,3*2*nres,mpi_real,&
1356 d_restart1,3*2*nres,mpi_real,king,&
1365 call mpi_gather(r_d,3*2*nres,mpi_real,&
1366 d_restart2,3*2*nres,mpi_real,king,&
1371 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1373 call xdrfint_(ixdrf, i2rep(i), iret)
1376 call xdrfint_(ixdrf, ifirst(i), iret)
1380 call xdrfint_(ixdrf, nupa(i,il), iret)
1384 call xdrfint_(ixdrf, ndowna(i,il), iret)
1390 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1397 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1404 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1410 call xdrfint_(ixdrf, nset, iret)
1412 call xdrfint_(ixdrf,mset(i), iret)
1415 call xdrfint_(ixdrf,i2set(i), iret)
1421 itmp=i_index(i,j,il,il1)
1422 call xdrfint_(ixdrf,itmp, iret)
1429 call xdrfclose_(ixdrf, iret)
1431 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1433 call xdrfint(ixdrf, i2rep(i), iret)
1436 call xdrfint(ixdrf, ifirst(i), iret)
1440 call xdrfint(ixdrf, nupa(i,il), iret)
1444 call xdrfint(ixdrf, ndowna(i,il), iret)
1450 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1457 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1464 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1471 call xdrfint(ixdrf, nset, iret)
1473 call xdrfint(ixdrf,mset(i), iret)
1476 call xdrfint(ixdrf,i2set(i), iret)
1482 itmp=i_index(i,j,il,il1)
1483 call xdrfint(ixdrf,itmp, iret)
1490 call xdrfclose(ixdrf, iret)
1494 end subroutine write1rst
1495 !-----------------------------------------------------------------------------
1496 subroutine write1traj
1498 ! implicit real*8 (a-h,o-z)
1499 ! include 'DIMENSIONS'
1501 ! include 'COMMON.MD'
1502 ! include 'COMMON.IOUNITS'
1503 ! include 'COMMON.REMD'
1504 ! include 'COMMON.SETUP'
1505 ! include 'COMMON.CHAIN'
1506 ! include 'COMMON.SBRIDGE'
1507 ! include 'COMMON.INTERACT'
1509 real(kind=4) :: t5_restart1(5)
1510 integer :: iret,itmp
1511 real(kind=4) :: xcoord(3,2*nres+2),prec
1512 real(kind=4) :: r_qfrag(50),r_qpair(100)
1513 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1514 real(kind=4) :: p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1515 real(kind=4) :: p_utheta(50*maxprocs),p_ugamma(100*maxprocs),&
1516 p_uscdiff(100*maxprocs)
1517 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1518 real(kind=4) :: r_c(3,2*nres+2)
1519 !el common /przechowalnia/ p_c
1521 integer :: i,j,il,ierr,ii,ixdrf
1523 call mpi_bcast(ii_write,1,mpi_integer,&
1527 print *,'traj1file',me,ii_write,ntwx_cache
1531 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1533 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1536 ! write (iout,*) "before gather write1traj: from node",ii
1538 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1540 t5_restart1(1)=totT_cache(ii)
1541 t5_restart1(2)=EK_cache(ii)
1542 t5_restart1(3)=potE_cache(ii)
1543 t5_restart1(4)=t_bath_cache(ii)
1544 t5_restart1(5)=Uconst_cache(ii)
1545 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1547 call mpi_gather(t5_restart1,5,mpi_real,&
1548 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1550 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1553 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1554 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1557 r_qfrag(i)=qfrag_cache(i,ii)
1560 r_qpair(i)=qpair_cache(i,ii)
1563 r_utheta(i)=utheta_cache(i,ii)
1564 r_ugamma(i)=ugamma_cache(i,ii)
1565 r_uscdiff(i)=uscdiff_cache(i,ii)
1568 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1569 p_qfrag,nfrag,mpi_real,king,&
1571 call mpi_gather(r_qpair,npair,mpi_real,&
1572 p_qpair,npair,mpi_real,king,&
1574 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1575 p_utheta,nfrag_back,mpi_real,king,&
1577 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1578 p_ugamma,nfrag_back,mpi_real,king,&
1580 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1581 p_uscdiff,nfrag_back,mpi_real,king,&
1585 write (iout,*) "p_qfrag"
1587 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1589 write (iout,*) "p_qpair"
1591 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1597 r_c(j,i)=c_cache(j,i,ii)
1601 call mpi_gather(r_c,3*2*nres,mpi_real,&
1602 p_c,3*2*nres,mpi_real,king,&
1608 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1609 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1610 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1611 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1612 call xdrfint_(ixdrf, nss, iret)
1615 call xdrfint(ixdrf, idssb(j)+nres, iret)
1616 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1618 call xdrfint_(ixdrf, ihpb(j), iret)
1619 call xdrfint_(ixdrf, jhpb(j), iret)
1622 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1623 call xdrfint_(ixdrf, iset_restart1(il), iret)
1625 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1628 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1631 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1632 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1633 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1638 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1643 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1647 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1651 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1652 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1653 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1654 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1655 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1656 call xdrfint(ixdrf, nss, iret)
1659 call xdrfint(ixdrf, idssb(j)+nres, iret)
1660 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1662 call xdrfint(ixdrf, ihpb(j), iret)
1663 call xdrfint(ixdrf, jhpb(j), iret)
1666 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1667 call xdrfint(ixdrf, iset_restart1(il), iret)
1669 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1672 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1675 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1676 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1677 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1682 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1687 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1691 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1697 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1699 if(me.eq.king) call xdrfclose(ixdrf, iret)
1701 do i=1,ntwx_cache-ii_write
1703 totT_cache(i)=totT_cache(ii_write+i)
1704 EK_cache(i)=EK_cache(ii_write+i)
1705 potE_cache(i)=potE_cache(ii_write+i)
1706 t_bath_cache(i)=t_bath_cache(ii_write+i)
1707 Uconst_cache(i)=Uconst_cache(ii_write+i)
1708 iset_cache(i)=iset_cache(ii_write+i)
1711 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1714 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1717 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1718 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1719 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1724 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1728 ntwx_cache=ntwx_cache-ii_write
1730 end subroutine write1traj
1731 !-----------------------------------------------------------------------------
1732 subroutine read1restart(i_index)
1734 ! implicit real*8 (a-h,o-z)
1735 ! include 'DIMENSIONS'
1737 ! include 'COMMON.MD'
1738 ! include 'COMMON.IOUNITS'
1739 ! include 'COMMON.REMD'
1740 ! include 'COMMON.SETUP'
1741 ! include 'COMMON.CHAIN'
1742 ! include 'COMMON.SBRIDGE'
1743 ! include 'COMMON.INTERACT'
1744 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1745 real(kind=4) :: r_d(3,2*nres),t5_restart1(5)
1746 integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1747 !el common /przechowalnia/ d_restart1
1748 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1750 write (*,*) "Processor",me," called read1restart"
1753 open(irest2,file=mremd_rst_name,status='unknown')
1754 read(irest2,*,err=334) i
1755 write(iout,*) "Reading old rst in ASCI format"
1757 call read1restart_old
1761 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1764 call xdrfint_(ixdrf, i2rep(i), iret)
1767 call xdrfint_(ixdrf, ifirst(i), iret)
1770 call xdrfint_(ixdrf, nupa(0,il), iret)
1772 call xdrfint_(ixdrf, nupa(i,il), iret)
1775 call xdrfint_(ixdrf, ndowna(0,il), iret)
1777 call xdrfint_(ixdrf, ndowna(i,il), iret)
1782 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1786 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1789 call xdrfint(ixdrf, i2rep(i), iret)
1792 call xdrfint(ixdrf, ifirst(i), iret)
1795 call xdrfint(ixdrf, nupa(0,il), iret)
1797 call xdrfint(ixdrf, nupa(i,il), iret)
1800 call xdrfint(ixdrf, ndowna(0,il), iret)
1802 call xdrfint(ixdrf, ndowna(i,il), iret)
1807 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1812 call mpi_scatter(t_restart1,5,mpi_real,&
1813 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1817 t_bath=t5_restart1(4)
1822 ! read(irest2,'(3e15.5)')
1823 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1826 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1828 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1834 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1835 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1845 ! read(irest2,'(3e15.5)')
1846 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1849 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1851 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1857 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1858 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1869 call xdrfint_(ixdrf, nset, iret)
1871 call xdrfint_(ixdrf,mset(i), iret)
1874 call xdrfint_(ixdrf,i2set(i), iret)
1880 call xdrfint_(ixdrf,itmp, iret)
1881 i_index(i,j,il,il1)=itmp
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
1908 call mpi_scatter(i2set,1,mpi_integer,&
1909 iset,1,mpi_integer,king,&
1914 if(me.eq.king) close(irest2)
1916 end subroutine read1restart
1917 !-----------------------------------------------------------------------------
1918 subroutine read1restart_old
1920 ! implicit real*8 (a-h,o-z)
1921 ! include 'DIMENSIONS'
1923 ! include 'COMMON.MD'
1924 ! include 'COMMON.IOUNITS'
1925 ! include 'COMMON.REMD'
1926 ! include 'COMMON.SETUP'
1927 ! include 'COMMON.CHAIN'
1928 ! include 'COMMON.SBRIDGE'
1929 ! include 'COMMON.INTERACT'
1930 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1931 real(kind=4) :: r_d(3,2*nres),t5_restart1(5)
1932 !el common /przechowalnia/ d_restart1
1934 integer :: i,j,il,ierr
1937 open(irest2,file=mremd_rst_name,status='unknown')
1938 read (irest2,*) (i2rep(i),i=0,nodes-1)
1939 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1941 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1942 read (irest2,*) ndowna(0,il),&
1943 (ndowna(i,il),i=1,ndowna(0,il))
1946 read(irest2,*) (t_restart1(j,il),j=1,4)
1949 call mpi_scatter(t_restart1,5,mpi_real,&
1950 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1954 t_bath=t5_restart1(4)
1959 read(irest2,'(3e15.5)') &
1960 (d_restart1(j,i+2*nres*il),j=1,3)
1964 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1965 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1975 read(irest2,'(3e15.5)') &
1976 (d_restart1(j,i+2*nres*il),j=1,3)
1980 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1981 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1987 if(me.eq.king) close(irest2)
1989 end subroutine read1restart_old
1990 !----------------------------------------------------------------
1991 subroutine alloc_MREMD_arrays
1993 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
1994 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1)) !(ntyp1))
1996 ! common /remdcommon/ in io: read_REMDpar
1997 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
1998 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
1999 ! common /remdrestart/
2000 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2002 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2003 allocate(ifirst(0:nodes)) !(maxprocs)
2004 allocate(nupa(0:nodes,0:2*nodes))
2005 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2006 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2007 allocate(iset_restart1(nodes)) !(maxprocs)
2008 ! common /traj1cache/
2009 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2010 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2011 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2012 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2013 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2014 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2015 allocate(utheta_cache(nfrag_back,max_cache_traj))
2016 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2017 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2018 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2021 end subroutine alloc_MREMD_arrays
2022 !-----------------------------------------------------------------------------
2023 !-----------------------------------------------------------------------------