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=8),dimension(:,:),allocatable :: t_restart1 !(5,maxprocs)
22 integer,dimension(:),allocatable :: iset_restart1 !(maxprocs)
24 real(kind=8),dimension(:),allocatable :: totT_cache,EK_cache,&
25 potE_cache,t_bath_cache,Uconst_cache !(max_cache_traj)
26 real(kind=8),dimension(:,:),allocatable :: qfrag_cache !(50,max_cache_traj)
27 real(kind=8),dimension(:,:),allocatable :: qpair_cache !(100,max_cache_traj)
28 real(kind=8),dimension(:,:),allocatable :: ugamma_cache,&
29 utheta_cache,uscdiff_cache !(maxfrag_back,max_cache_traj)
30 real(kind=8),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
115 write(iout,*) "jestesmy na poczatku MREMD"
119 if(me.eq.king.or..not.out1file) then
120 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
121 write (iout,*) "NREP=",nrep
124 write(iout,*) "jestesmy na poczatku MREMD"
126 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
127 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
129 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
131 !d print *,'MREMD',nodes
132 !d print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
133 !de write (iout,*) "Start MREMD: me",me," t_bath",t_bath
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,nres2*nodes))
198 if(.not.allocated(d_restart2)) allocate(d_restart2(3,nres2*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
352 stdfp=dsqrt(2*Rb*t_bath/d_time)
354 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
357 ! print *,'irep',me,t_bath
359 if (me.eq.king .or. .not. out1file) &
360 write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
361 call rescale_weights(t_bath)
365 !------copy MD--------------
366 ! The driver for molecular dynamics subroutines
367 !------------------------------------------------
373 if(me.eq.king.or..not.out1file) &
374 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
380 ! Determine the inverse of the inertia matrix.
381 call setup_MD_matrices
385 if (me.eq.king .or. .not. out1file) &
386 write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
387 stdfp=dsqrt(2*Rb*t_bath/d_time)
389 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
391 call rescale_weights(t_bath)
395 t_MDsetup = MPI_Wtime()-tt0
397 t_MDsetup = tcpu()-tt0
400 ! Entering the MD loop
406 if (lang.eq.2 .or. lang.eq.3) then
410 call sd_verlet_p_setup
412 call sd_verlet_ciccotti_setup
416 pfric0_mat(i,j,0)=pfric_mat(i,j)
417 afric0_mat(i,j,0)=afric_mat(i,j)
418 vfric0_mat(i,j,0)=vfric_mat(i,j)
419 prand0_mat(i,j,0)=prand_mat(i,j)
420 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
421 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
426 flag_stoch(i)=.false.
430 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
432 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
436 else if (lang.eq.1 .or. lang.eq.4) then
440 if (me.eq.king .or. .not. out1file) &
441 write(iout,*) 'Setup time',time00-walltime
444 t_langsetup=MPI_Wtime()-tt0
447 t_langsetup=tcpu()-tt0
453 do while(.not.end_of_run)
455 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
456 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
458 if (lang.gt.0 .and. surfarea .and. &
459 mod(itime,reset_fricmat).eq.0) then
460 if (lang.eq.2 .or. lang.eq.3) then
464 call sd_verlet_p_setup
466 call sd_verlet_ciccotti_setup
470 pfric0_mat(i,j,0)=pfric_mat(i,j)
471 afric0_mat(i,j,0)=afric_mat(i,j)
472 vfric0_mat(i,j,0)=vfric_mat(i,j)
473 prand0_mat(i,j,0)=prand_mat(i,j)
474 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
475 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
480 flag_stoch(i)=.false.
483 else if (lang.eq.1 .or. lang.eq.4) then
486 write (iout,'(a,i10)') &
487 "Friction matrix reset based on surface area, itime",itime
489 if (reset_vel .and. tbf .and. lang.eq.0 &
490 .and. mod(itime,count_reset_vel).eq.0) then
492 if (me.eq.king .or. .not. out1file) &
493 write(iout,'(a,f20.2)') &
494 "Velocities reset to random values, time",totT
497 d_t_old(j,i)=d_t(j,i)
501 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
505 d_t(j,0)=d_t(j,0)-vcm(j)
508 kinetic_T=2.0d0/(dimen3*Rb)*EK
509 scalfac=dsqrt(T_bath/kinetic_T)
510 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
513 d_t_old(j,i)=scalfac*d_t(j,i)
519 ! Time-reversible RESPA algorithm
520 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
521 call RESPA_step(itime)
523 ! Variable time step algorithm.
524 call velverlet_step(itime)
528 call brown_step(itime)
530 print *,"Brown dynamics not here!"
532 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
538 if (mod(itime,ntwe).eq.0) call statout(itime)
540 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
541 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
543 call hairpin(.true.,nharp,iharp)
544 call secondary2(.true.)
545 call pdbout(potE,tytul,ipdb)
550 if (mod(itime,ntwx).eq.0.and.traj1file) then
551 if(ntwx_cache.lt.max_cache_traj_use) then
552 ntwx_cache=ntwx_cache+1
554 if (max_cache_traj_use.ne.1) &
555 print *,itime,"processor ",me," over cache ",ntwx_cache
558 totT_cache(i)=totT_cache(i+1)
559 EK_cache(i)=EK_cache(i+1)
560 potE_cache(i)=potE_cache(i+1)
561 t_bath_cache(i)=t_bath_cache(i+1)
562 Uconst_cache(i)=Uconst_cache(i+1)
563 iset_cache(i)=iset_cache(i+1)
566 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
569 qpair_cache(ii,i)=qpair_cache(ii,i+1)
572 utheta_cache(ii,i)=utheta_cache(ii,i+1)
573 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
574 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
580 c_cache(j,ii,i)=c_cache(j,ii,i+1)
586 totT_cache(ntwx_cache)=totT
587 EK_cache(ntwx_cache)=EK
588 potE_cache(ntwx_cache)=potE
589 t_bath_cache(ntwx_cache)=t_bath
590 Uconst_cache(ntwx_cache)=Uconst
591 iset_cache(ntwx_cache)=iset
594 qfrag_cache(i,ntwx_cache)=qfrag(i)
597 qpair_cache(i,ntwx_cache)=qpair(i)
600 utheta_cache(i,ntwx_cache)=utheta(i)
601 ugamma_cache(i,ntwx_cache)=ugamma(i)
602 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
607 c_cache(j,i,ntwx_cache)=c(j,i)
612 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
613 .and..not.restart1file) then
616 open(irest1,file=mremd_rst_name,status='unknown')
617 write (irest1,*) "i2rep"
618 write (irest1,*) (i2rep(i),i=0,nodes-1)
619 write (irest1,*) "ifirst"
620 write (irest1,*) (ifirst(i),i=1,remd_m(1))
622 write (irest1,*) "nupa",il
623 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
624 write (irest1,*) "ndowna",il
625 write (irest1,*) ndowna(0,il),&
626 (ndowna(i,il),i=1,ndowna(0,il))
629 write (irest1,*) "nset"
630 write (irest1,*) nset
631 write (irest1,*) "mset"
632 write (irest1,*) (mset(i),i=1,nset)
633 write (irest1,*) "i2set"
634 write (irest1,*) (i2set(i),i=0,nodes-1)
635 write (irest1,*) "i_index"
639 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
647 open(irest2,file=rest2name,status='unknown')
648 write(irest2,*) totT,EK,potE,totE,t_bath
650 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
653 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
656 write (irest2,*) iset
663 ! forced synchronization
664 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
665 .and. .not. mremdsync) then
667 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
669 call mpi_recv(itime_master, 1, MPI_INTEGER,&
670 0,101,CG_COMM, status, ierr)
671 call mpi_barrier(CG_COMM, ierr)
672 !deb if (out1file.or.traj1file) then
673 !deb call mpi_gather(itime,1,mpi_integer,
674 !deb & icache_all,1,mpi_integer,king,
677 call mpi_gather(ntwx_cache,1,mpi_integer,&
678 icache_all,1,mpi_integer,king,&
681 write(iout,*) 'REMD synchro at',itime_master,itime
682 if (itime_master.ge.n_timestep .or. ovrtim()) &
684 !time call flush(iout)
689 if ((mod(itime,nstex).eq.0.and.me.eq.king &
690 .or.end_of_run.and.me.eq.king ) &
691 .and. .not. mremdsync ) then
694 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
695 CG_COMM, ireqi(i), ierr)
696 !d write(iout,*) 'REMD synchro with',i
699 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
700 call mpi_barrier(CG_COMM, ierr)
702 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
703 if (out1file.or.traj1file) then
704 !deb call mpi_gather(itime,1,mpi_integer,
705 !deb & itime_all,1,mpi_integer,king,
707 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
708 !deb & (itime_all(i),i=1,nodes)
710 !deb imin_itime=itime_all(1)
712 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
714 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
715 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
716 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
717 call mpi_gather(ntwx_cache,1,mpi_integer,&
718 icache_all,1,mpi_integer,king,&
720 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
721 ! & (icache_all(i),i=1,nodes)
722 ii_write=icache_all(1)
724 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
726 ! write(iout,*) "MIN ii_write=",ii_write
729 !time call flush(iout)
731 if(mremdsync .and. mod(itime,nstex).eq.0) then
733 if (me.eq.king .or. .not. out1file) &
734 write(iout,*) 'REMD synchro at',itime
737 call mpi_gather(ntwx_cache,1,mpi_integer,&
738 icache_all,1,mpi_integer,king,&
741 write(iout,'(a19,8000i8)') ' ntwx_cache',&
742 (icache_all(i),i=1,nodes)
743 ii_write=icache_all(1)
745 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
747 write(iout,*) "MIN ii_write=",ii_write
753 ! Update the time safety limiy
754 if (time001-time00.gt.safety) then
755 safety=time001-time00+600
756 write (iout,*) "****** SAFETY increased to",safety," s"
758 if (ovrtim()) end_of_run=.true.
760 if(synflag.and..not.end_of_run) then
764 write(iout,*) 'REMD before',me,t_bath
766 ! call mpi_gather(t_bath,1,mpi_double_precision,
767 ! & remd_t_bath,1,mpi_double_precision,king,
769 potEcomp(n_ene+1)=t_bath
771 potEcomp(n_ene+2)=iset
772 if (iset.lt.nset) then
776 potEcomp(n_ene+3)=Uconst
783 potEcomp(n_ene+4)=Uconst
787 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
788 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
791 call mpi_gather(elow,1,mpi_double_precision,&
792 elowi,1,mpi_double_precision,king,&
794 call mpi_gather(ehigh,1,mpi_double_precision,&
795 ehighi,1,mpi_double_precision,king,&
800 if (me.eq.king .or. .not. out1file) then
801 write(iout,*) 'REMD gather times=',time03-time01 &
805 if (restart1file) call write1rst(i_index)
808 if (me.eq.king .or. .not. out1file) then
809 write(iout,*) 'REMD writing rst time=',time04-time03
812 if (traj1file) call write1traj
814 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
815 !deb & icache_all,1,mpi_integer,king,
817 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
818 !deb & (icache_all(i),i=1,nodes)
823 if (me.eq.king .or. .not. out1file) then
824 write(iout,*) 'REMD writing traj time=',time05-time04
831 remd_t_bath(i)=remd_ene(n_ene+1,i)
832 iremd_iset(i)=remd_ene(n_ene+2,i)
835 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
837 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
841 write(iout,*) 'REMD exchange temp,ene'
843 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
844 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)))
860 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
861 if(i.gt.0.and.nupa(0,i).gt.0) then
863 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
865 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
867 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
869 ! do while (iex.eq.i)
870 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
871 iex=nupa(iran_num(1,int(nupa(0,i))),i)
873 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
875 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
877 ! Swap temperatures between conformations i and iex with recalculating the free energies
878 ! following temperature changes.
879 ene_iex_iex=remd_ene(0,iex)
880 ene_i_i=remd_ene(0,i)
881 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
882 ! & " iex",iex," ene_iex_iex",ene_iex_iex
883 ! write (iout,*) "rescaling weights with temperature",
886 call rescale_weights(remd_t_bath(i))
888 ! write (iout,*) "0,iex",remd_t_bath(i)
889 ! call enerprint(remd_ene(0,iex))
891 call sum_energy(remd_ene(0,iex),.false.)
892 ene_iex_i=remd_ene(0,iex)
893 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
895 ! write (iout,*) "0,i",remd_t_bath(i)
896 ! call enerprint(remd_ene(0,i))
898 call sum_energy(remd_ene(0,i),.false.)
899 ! write (iout,*) "ene_i_i",remd_ene(0,i)
901 ! write (iout,*) "rescaling weights with temperature",
903 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
904 write (iout,*) "ERROR: inconsistent energies:",i,&
905 ene_i_i,remd_ene(0,i)
907 call rescale_weights(remd_t_bath(iex))
909 ! write (iout,*) "0,i",remd_t_bath(iex)
910 ! call enerprint(remd_ene(0,i))
912 call sum_energy(remd_ene(0,i),.false.)
913 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
915 ene_i_iex=remd_ene(0,i)
917 ! write (iout,*) "0,iex",remd_t_bath(iex)
918 ! call enerprint(remd_ene(0,iex))
920 call sum_energy(remd_ene(0,iex),.false.)
921 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
922 write (iout,*) "ERROR: inconsistent energies:",iex,&
923 ene_iex_iex,remd_ene(0,iex)
925 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
926 ! write (iout,*) "i",i," iex",iex
927 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
928 ! & " ene_i_iex",ene_i_iex,
929 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
931 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
932 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
934 ! write(iout,*) 'delta',delta
935 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
936 ! & (remd_ene(i)-remd_ene(iex))/Rb/
937 ! & (remd_t_bath(i)*remd_t_bath(iex))
939 if (delta .gt. 50.0d0) then
945 else if (delta.lt.-50.0d0) then
954 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
955 xxx=ran_number(0.0d0,1.0d0)
956 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
958 if (delta .gt. xxx) then
960 remd_t_bath(i)=remd_t_bath(iex)
962 remd_ene(0,i)=ene_i_iex
963 remd_ene(0,iex)=ene_iex_i
969 ehighi(i)=ehighi(iex)
976 nupa(k,i)=nupa(k,iex)
979 ndowna(k,i)=ndowna(k,iex)
983 if (ifirst(il).eq.i) ifirst(il)=iex
985 if (nupa(k,il).eq.i) then
987 elseif (nupa(k,il).eq.iex) then
992 if (ndowna(k,il).eq.i) then
994 elseif (ndowna(k,il).eq.iex) then
1000 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1002 i2rep(i-1)=i2rep(iex-1)
1005 ! write(iout,*) 'exchange',i,iex
1006 ! write (iout,'(a8,100i4)') "@ ifirst",
1007 ! & (ifirst(k),k=1,remd_m(1))
1009 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1010 ! & (nupa(k,il),k=1,nupa(0,il))
1011 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1012 ! & (ndowna(k,il),k=1,ndowna(0,il))
1017 remd_ene(0,iex)=ene_iex_iex
1018 remd_ene(0,i)=ene_i_i
1024 !d write (iout,*) "exchange completed"
1028 !d write(iout,*) "########",ii
1030 i_temp=iran_num(1,nrep)
1031 i_mult=iran_num(1,remd_m(i_temp))
1032 i_iset=iran_num(1,nset)
1033 i_mset=iran_num(1,mset(i_iset))
1034 i=i_index(i_temp,i_mult,i_iset,i_mset)
1036 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1039 !d write(iout,*) "i_dir=",i_dir
1041 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1044 i_mult1=iran_num(1,remd_m(i_temp1))
1046 i_mset1=iran_num(1,mset(i_iset1))
1047 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1049 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1052 i_mult1=iran_num(1,remd_m(i_temp1))
1054 i_mset1=iran_num(1,mset(i_iset1))
1055 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1056 econstr_temp_i=remd_ene(20,i)
1057 econstr_temp_iex=remd_ene(20,iex)
1058 remd_ene(20,i)=remd_ene(n_ene+3,i)
1059 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1061 elseif(remd_m(i_temp+1).gt.0.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)
1077 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1080 ! Swap temperatures between conformations i and iex with recalculating the free energies
1081 ! following temperature changes.
1082 ene_iex_iex=remd_ene(0,iex)
1083 ene_i_i=remd_ene(0,i)
1084 !o write (iout,*) "rescaling weights with temperature",
1086 call rescale_weights(remd_t_bath(i))
1088 call sum_energy(remd_ene(0,iex),.false.)
1089 ene_iex_i=remd_ene(0,iex)
1090 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1091 ! call sum_energy(remd_ene(0,i),.false.)
1092 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1093 ! write (iout,*) "rescaling weights with temperature",
1094 ! & remd_t_bath(iex)
1095 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1096 ! write (iout,*) "ERROR: inconsistent energies:",i,
1097 ! & ene_i_i,remd_ene(0,i)
1099 call rescale_weights(remd_t_bath(iex))
1100 call sum_energy(remd_ene(0,i),.false.)
1101 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1102 ene_i_iex=remd_ene(0,i)
1103 ! call sum_energy(remd_ene(0,iex),.false.)
1104 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1105 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1106 ! & ene_iex_iex,remd_ene(0,iex)
1108 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1109 ! write (iout,*) "i",i," iex",iex
1110 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1111 !d & " ene_i_iex",ene_i_iex,
1112 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1113 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1114 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1116 !d write(iout,*) 'delta',delta
1117 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1118 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1119 ! & (remd_t_bath(i)*remd_t_bath(iex))
1120 if (delta .gt. 50.0d0) then
1125 if (i_dir.eq.1.or.i_dir.eq.3) &
1126 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1127 if (i_dir.eq.2.or.i_dir.eq.3) &
1128 iremd_tot_usa(int(i2set(i-1)))= &
1129 iremd_tot_usa(int(i2set(i-1)))+1
1130 xxx=ran_number(0.0d0,1.0d0)
1131 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1132 if (delta .gt. xxx) then
1134 remd_t_bath(i)=remd_t_bath(iex)
1135 remd_t_bath(iex)=tmp
1138 iremd_iset(i)=iremd_iset(iex)
1139 iremd_iset(iex)=itmp
1141 remd_ene(0,i)=ene_i_iex
1142 remd_ene(0,iex)=ene_iex_i
1144 if (i_dir.eq.1.or.i_dir.eq.3) &
1145 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1148 i2rep(i-1)=i2rep(iex-1)
1151 if (i_dir.eq.2.or.i_dir.eq.3) &
1152 iremd_acc_usa(int(i2set(i-1)))= &
1153 iremd_acc_usa(int(i2set(i-1)))+1
1156 i2set(i-1)=i2set(iex-1)
1159 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1160 i_index(i_temp,i_mult,i_iset,i_mset)= &
1161 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1162 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1165 remd_ene(0,iex)=ene_iex_iex
1166 remd_ene(0,i)=ene_i_i
1167 remd_ene(20,iex)=econstr_temp_iex
1168 remd_ene(20,i)=econstr_temp_i
1172 !d do il1=1,mset(il)
1175 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1188 !-------------------------------------
1189 write (iout,*) "NREP",nrep
1191 if(iremd_tot(i).ne.0) &
1192 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1193 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1198 if(iremd_tot_usa(i).ne.0) &
1199 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1200 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1206 !d write (iout,'(a6,100i4)') "ifirst",
1207 !d & (ifirst(i),i=1,remd_m(1))
1209 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1210 !d & (nupa(i,il),i=1,nupa(0,il))
1211 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1212 !d & (ndowna(i,il),i=1,ndowna(0,il))
1217 !d write (iout,*) "Before scatter"
1219 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1220 t_bath,1,mpi_double_precision,king,&
1222 !d write (iout,*) "After scatter"
1225 call mpi_scatter(iremd_iset,1,mpi_integer,&
1226 iset,1,mpi_integer,king,&
1230 if (me.eq.king .or. .not. out1file) then
1231 write(iout,*) 'REMD scatter time=',time07-time06
1235 call mpi_scatter(elowi,1,mpi_double_precision,&
1236 elow,1,mpi_double_precision,king,&
1238 call mpi_scatter(ehighi,1,mpi_double_precision,&
1239 ehigh,1,mpi_double_precision,king,&
1242 call rescale_weights(t_bath)
1243 !o write (iout,*) "Processor",me,
1244 !o & " rescaling weights with temperature",t_bath
1246 stdfp=dsqrt(2*Rb*t_bath/d_time)
1248 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1251 !de write(iout,*) 'REMD after',me,t_bath
1253 if (me.eq.king .or. .not. out1file) then
1254 write(iout,*) 'REMD exchange time=',time08-time00
1260 if (restart1file) then
1261 if (me.eq.king .or. .not. out1file) &
1262 write(iout,*) 'writing restart at the end of run'
1263 call write1rst(i_index)
1266 if (traj1file) call write1traj
1268 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1269 !deb & icache_all,1,mpi_integer,king,
1270 !deb & CG_COMM,ierr)
1271 !deb write(iout,'(a40,8000i8)')
1272 !deb & ' ntwx_cache after traj1file at the end',
1273 !deb & (icache_all(i),i=1,nodes)
1278 t_MD=MPI_Wtime()-tt0
1282 if (me.eq.king .or. .not. out1file) then
1283 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1285 'MD calculations setup:',t_MDsetup,&
1286 'Energy & gradient evaluation:',t_enegrad,&
1287 'Stochastic MD setup:',t_langsetup,&
1288 'Stochastic MD step setup:',t_sdsetup,&
1290 write (iout,'(/28(1h=),a25,27(1h=))') &
1291 ' End of MD calculation '
1293 !el common /przechowalnia/
1294 deallocate(d_restart1)
1295 deallocate(d_restart2)
1299 end subroutine MREMD
1300 !-----------------------------------------------------------------------------
1301 subroutine write1rst(i_index)
1304 ! implicit real*8 (a-h,o-z)
1305 ! include 'DIMENSIONS'
1307 ! include 'COMMON.MD'
1308 ! include 'COMMON.IOUNITS'
1309 ! include 'COMMON.REMD'
1310 ! include 'COMMON.SETUP'
1311 ! include 'COMMON.CHAIN'
1312 ! include 'COMMON.SBRIDGE'
1313 ! include 'COMMON.INTERACT'
1315 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1316 !el d_restart2(3,2*nres*maxprocs)
1317 real(kind=4) :: r_d(3,2*nres)
1318 real(kind=4) :: t5_restart1(5)
1319 integer :: iret,itmp
1320 integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1321 !el common /przechowalnia/ d_restart1,d_restart2
1322 integer :: i,j,il,il1,ierr,ixdrf
1327 t5_restart1(4)=t_bath
1328 t5_restart1(5)=Uconst
1330 call mpi_gather(t5_restart1,5,mpi_real,&
1331 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1339 call mpi_gather(r_d,3*2*nres,mpi_real,&
1340 d_restart1,3*2*nres,mpi_real,king,&
1349 call mpi_gather(r_d,3*2*nres,mpi_real,&
1350 d_restart2,3*2*nres,mpi_real,king,&
1355 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1357 call xdrfint_(ixdrf, i2rep(i), iret)
1360 call xdrfint_(ixdrf, ifirst(i), iret)
1364 call xdrfint_(ixdrf, nupa(i,il), iret)
1368 call xdrfint_(ixdrf, ndowna(i,il), iret)
1374 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1381 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1388 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1394 call xdrfint_(ixdrf, nset, iret)
1396 call xdrfint_(ixdrf,mset(i), iret)
1399 call xdrfint_(ixdrf,i2set(i), iret)
1405 itmp=i_index(i,j,il,il1)
1406 call xdrfint_(ixdrf,itmp, iret)
1413 call xdrfclose_(ixdrf, iret)
1415 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1417 call xdrfint(ixdrf, i2rep(i), iret)
1420 call xdrfint(ixdrf, ifirst(i), iret)
1424 call xdrfint(ixdrf, nupa(i,il), iret)
1428 call xdrfint(ixdrf, ndowna(i,il), iret)
1434 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1441 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1448 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1455 call xdrfint(ixdrf, nset, iret)
1457 call xdrfint(ixdrf,mset(i), iret)
1460 call xdrfint(ixdrf,i2set(i), iret)
1466 itmp=i_index(i,j,il,il1)
1467 call xdrfint(ixdrf,itmp, iret)
1474 call xdrfclose(ixdrf, iret)
1478 end subroutine write1rst
1479 !-----------------------------------------------------------------------------
1480 subroutine write1traj
1482 ! implicit real*8 (a-h,o-z)
1483 ! include 'DIMENSIONS'
1485 ! include 'COMMON.MD'
1486 ! include 'COMMON.IOUNITS'
1487 ! include 'COMMON.REMD'
1488 ! include 'COMMON.SETUP'
1489 ! include 'COMMON.CHAIN'
1490 ! include 'COMMON.SBRIDGE'
1491 ! include 'COMMON.INTERACT'
1493 real(kind=4) :: t5_restart1(5)
1494 integer :: iret,itmp
1495 real(kind=4) :: xcoord(3,2*nres+2),prec
1496 real(kind=4) :: r_qfrag(50),r_qpair(100)
1497 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1498 real(kind=4) :: p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1499 real(kind=4) :: p_utheta(50*maxprocs),p_ugamma(100*maxprocs),&
1500 p_uscdiff(100*maxprocs)
1501 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1502 real(kind=4) :: r_c(3,2*nres+2)
1503 !el common /przechowalnia/ p_c
1505 integer :: i,j,il,ierr,ii,ixdrf
1507 call mpi_bcast(ii_write,1,mpi_integer,&
1511 print *,'traj1file',me,ii_write,ntwx_cache
1515 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1517 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1520 t5_restart1(1)=totT_cache(ii)
1521 t5_restart1(2)=EK_cache(ii)
1522 t5_restart1(3)=potE_cache(ii)
1523 t5_restart1(4)=t_bath_cache(ii)
1524 t5_restart1(5)=Uconst_cache(ii)
1525 call mpi_gather(t5_restart1,5,mpi_real,&
1526 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1528 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1529 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1532 r_qfrag(i)=qfrag_cache(i,ii)
1535 r_qpair(i)=qpair_cache(i,ii)
1538 r_utheta(i)=utheta_cache(i,ii)
1539 r_ugamma(i)=ugamma_cache(i,ii)
1540 r_uscdiff(i)=uscdiff_cache(i,ii)
1543 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1544 p_qfrag,nfrag,mpi_real,king,&
1546 call mpi_gather(r_qpair,npair,mpi_real,&
1547 p_qpair,npair,mpi_real,king,&
1549 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1550 p_utheta,nfrag_back,mpi_real,king,&
1552 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1553 p_ugamma,nfrag_back,mpi_real,king,&
1555 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1556 p_uscdiff,nfrag_back,mpi_real,king,&
1560 write (iout,*) "p_qfrag"
1562 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1564 write (iout,*) "p_qpair"
1566 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1572 r_c(j,i)=c_cache(j,i,ii)
1576 call mpi_gather(r_c,3*2*nres,mpi_real,&
1577 p_c,3*2*nres,mpi_real,king,&
1583 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1584 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1585 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1586 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1587 call xdrfint_(ixdrf, nss, iret)
1590 call xdrfint(ixdrf, idssb(j)+nres, iret)
1591 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1593 call xdrfint_(ixdrf, ihpb(j), iret)
1594 call xdrfint_(ixdrf, jhpb(j), iret)
1597 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1598 call xdrfint_(ixdrf, iset_restart1(il), iret)
1600 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1603 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1606 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1607 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1608 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1613 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1618 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1622 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1626 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1627 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1628 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1629 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1630 call xdrfint(ixdrf, nss, iret)
1633 call xdrfint(ixdrf, idssb(j)+nres, iret)
1634 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1636 call xdrfint(ixdrf, ihpb(j), iret)
1637 call xdrfint(ixdrf, jhpb(j), iret)
1640 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1641 call xdrfint(ixdrf, iset_restart1(il), iret)
1643 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1646 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1649 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1650 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1651 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1656 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1661 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1665 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1671 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1673 if(me.eq.king) call xdrfclose(ixdrf, iret)
1675 do i=1,ntwx_cache-ii_write
1677 totT_cache(i)=totT_cache(ii_write+i)
1678 EK_cache(i)=EK_cache(ii_write+i)
1679 potE_cache(i)=potE_cache(ii_write+i)
1680 t_bath_cache(i)=t_bath_cache(ii_write+i)
1681 Uconst_cache(i)=Uconst_cache(ii_write+i)
1682 iset_cache(i)=iset_cache(ii_write+i)
1685 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1688 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1691 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1692 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1693 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1698 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1702 ntwx_cache=ntwx_cache-ii_write
1704 end subroutine write1traj
1705 !-----------------------------------------------------------------------------
1706 subroutine read1restart(i_index)
1708 ! implicit real*8 (a-h,o-z)
1709 ! include 'DIMENSIONS'
1711 ! include 'COMMON.MD'
1712 ! include 'COMMON.IOUNITS'
1713 ! include 'COMMON.REMD'
1714 ! include 'COMMON.SETUP'
1715 ! include 'COMMON.CHAIN'
1716 ! include 'COMMON.SBRIDGE'
1717 ! include 'COMMON.INTERACT'
1718 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1719 real(kind=4) :: r_d(3,2*nres),t5_restart1(5)
1720 integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1721 !el common /przechowalnia/ d_restart1
1722 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1724 write (*,*) "Processor",me," called read1restart"
1727 open(irest2,file=mremd_rst_name,status='unknown')
1728 read(irest2,*,err=334) i
1729 write(iout,*) "Reading old rst in ASCI format"
1731 call read1restart_old
1735 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1738 call xdrfint_(ixdrf, i2rep(i), iret)
1741 call xdrfint_(ixdrf, ifirst(i), iret)
1744 call xdrfint_(ixdrf, nupa(0,il), iret)
1746 call xdrfint_(ixdrf, nupa(i,il), iret)
1749 call xdrfint_(ixdrf, ndowna(0,il), iret)
1751 call xdrfint_(ixdrf, ndowna(i,il), iret)
1756 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1760 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1763 call xdrfint(ixdrf, i2rep(i), iret)
1766 call xdrfint(ixdrf, ifirst(i), iret)
1769 call xdrfint(ixdrf, nupa(0,il), iret)
1771 call xdrfint(ixdrf, nupa(i,il), iret)
1774 call xdrfint(ixdrf, ndowna(0,il), iret)
1776 call xdrfint(ixdrf, ndowna(i,il), iret)
1781 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1786 call mpi_scatter(t_restart1,5,mpi_real,&
1787 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1791 t_bath=t5_restart1(4)
1796 ! read(irest2,'(3e15.5)')
1797 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1800 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1802 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1808 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1809 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1819 ! read(irest2,'(3e15.5)')
1820 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1823 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1825 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1831 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1832 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1843 call xdrfint_(ixdrf, nset, iret)
1845 call xdrfint_(ixdrf,mset(i), iret)
1848 call xdrfint_(ixdrf,i2set(i), iret)
1854 call xdrfint_(ixdrf,itmp, iret)
1855 i_index(i,j,il,il1)=itmp
1863 call xdrfint(ixdrf, nset, iret)
1865 call xdrfint(ixdrf,mset(i), iret)
1868 call xdrfint(ixdrf,i2set(i), iret)
1874 call xdrfint(ixdrf,itmp, iret)
1875 i_index(i,j,il,il1)=itmp
1882 call mpi_scatter(i2set,1,mpi_integer,&
1883 iset,1,mpi_integer,king,&
1888 if(me.eq.king) close(irest2)
1890 end subroutine read1restart
1891 !-----------------------------------------------------------------------------
1892 subroutine read1restart_old
1894 ! implicit real*8 (a-h,o-z)
1895 ! include 'DIMENSIONS'
1897 ! include 'COMMON.MD'
1898 ! include 'COMMON.IOUNITS'
1899 ! include 'COMMON.REMD'
1900 ! include 'COMMON.SETUP'
1901 ! include 'COMMON.CHAIN'
1902 ! include 'COMMON.SBRIDGE'
1903 ! include 'COMMON.INTERACT'
1904 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1905 real(kind=4) :: r_d(3,2*nres),t5_restart1(5)
1906 !el common /przechowalnia/ d_restart1
1908 integer :: i,j,il,ierr
1911 open(irest2,file=mremd_rst_name,status='unknown')
1912 read (irest2,*) (i2rep(i),i=0,nodes-1)
1913 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1915 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1916 read (irest2,*) ndowna(0,il),&
1917 (ndowna(i,il),i=1,ndowna(0,il))
1920 read(irest2,*) (t_restart1(j,il),j=1,4)
1923 call mpi_scatter(t_restart1,5,mpi_real,&
1924 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1928 t_bath=t5_restart1(4)
1933 read(irest2,'(3e15.5)') &
1934 (d_restart1(j,i+2*nres*il),j=1,3)
1938 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1939 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1949 read(irest2,'(3e15.5)') &
1950 (d_restart1(j,i+2*nres*il),j=1,3)
1954 call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
1955 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1961 if(me.eq.king) close(irest2)
1963 end subroutine read1restart_old
1964 !----------------------------------------------------------------
1965 subroutine alloc_MREMD_arrays
1967 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
1968 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1)) !(ntyp1))
1970 ! common /remdcommon/ in io: read_REMDpar
1971 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
1972 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
1973 ! common /remdrestart/
1974 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
1976 allocate(i2set(0:2*nodes)) !(0:maxprocs)
1977 allocate(ifirst(0:nodes)) !(maxprocs)
1978 allocate(nupa(0:nodes,0:2*nodes))
1979 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
1980 allocate(t_restart1(5,nodes)) !(5,maxprocs)
1981 allocate(iset_restart1(nodes)) !(maxprocs)
1982 ! common /traj1cache/
1983 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
1984 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
1985 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
1986 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
1987 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
1988 allocate(ugamma_cache(nfrag_back,max_cache_traj))
1989 allocate(utheta_cache(nfrag_back,max_cache_traj))
1990 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
1991 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
1992 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
1995 end subroutine alloc_MREMD_arrays
1996 !-----------------------------------------------------------------------------
1997 !-----------------------------------------------------------------------------