2 !-----------------------------------------------------------------------------
10 ! use control_data, only:maxprocs
14 !-----------------------------------------------------------------------------
16 ! common /remdrestart/
17 integer(kind=2),dimension(:),allocatable :: i2set !(0:maxprocs)
18 integer(kind=2),dimension(:),allocatable :: ifirst !(maxprocs)
19 integer(kind=2),dimension(:,:),allocatable :: nupa,&
20 ndowna !(0:maxprocs/4,0:maxprocs)
21 real(kind=4),dimension(:,:),allocatable :: t_restart1 !(5,maxprocs)
22 integer,dimension(:),allocatable :: iset_restart1 !(maxprocs)
24 real(kind=4),dimension(:),allocatable :: totT_cache,EK_cache,&
25 potE_cache,t_bath_cache,Uconst_cache !(max_cache_traj)
26 real(kind=4),dimension(:,:),allocatable :: qfrag_cache !(50,max_cache_traj)
27 real(kind=4),dimension(:,:),allocatable :: qpair_cache !(100,max_cache_traj)
28 real(kind=4),dimension(:,:),allocatable :: ugamma_cache,&
29 utheta_cache,uscdiff_cache !(maxfrag_back,max_cache_traj)
30 real(kind=4),dimension(:,:,:),allocatable :: c_cache !(3,maxres2+2,max_cache_traj)
31 integer :: ntwx_cache,ii_write !,max_cache_traj_use
32 integer,dimension(:),allocatable :: iset_cache !(max_cache_traj)
33 !-----------------------------------------------------------------------------
34 ! common /przechowalnia/
35 real(kind=4),dimension(:,:),allocatable :: d_restart1 !(3,2*nres*maxprocs)
36 real(kind=4),dimension(:,:),allocatable :: d_restart2 !(3,2*nres*maxprocs)
37 real(kind=4),dimension(:,:),allocatable :: p_c !(3,(nres2+2)*maxprocs)
38 !-----------------------------------------------------------------------------
41 !-----------------------------------------------------------------------------
43 !-----------------------------------------------------------------------------
45 !-----------------------------------------------------------------------------
47 !-----------------------------------------------------------------------------
52 use control, only:tcpu,ovrtim
53 use io_base, only:ilen
56 use random, only: iran_num,ran_number
57 use compare, only:hairpin,secondary2
58 use io, only:cartout,statout
59 ! implicit real*8 (a-h,o-z)
60 ! include 'DIMENSIONS'
62 ! include 'COMMON.CONTROL'
63 ! include 'COMMON.VAR'
66 ! include 'COMMON.LANGEVIN'
68 ! include 'COMMON.LANGEVIN.lang0'
70 ! include 'COMMON.CHAIN'
71 ! include 'COMMON.DERIV'
72 ! include 'COMMON.GEO'
73 ! include 'COMMON.LOCAL'
74 ! include 'COMMON.INTERACT'
75 ! include 'COMMON.IOUNITS'
76 ! include 'COMMON.NAMES'
77 ! include 'COMMON.TIME1'
78 ! include 'COMMON.REMD'
79 ! include 'COMMON.SETUP'
80 ! include 'COMMON.MUCA'
81 ! include 'COMMON.HAIRPIN'
83 real(kind=8),dimension(3) :: L,vcm
84 real(kind=8) :: energia(0:n_ene)
86 real(kind=8) :: remd_t_bath(Nprocs) !(maxprocs)
87 integer :: iremd_iset(Nprocs) !(maxprocs)
88 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
89 ! (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
90 real(kind=8) :: remd_ene(0:n_ene+4,Nprocs) !(0:n_ene+4,maxprocs)
91 integer :: iremd_acc(Nprocs),iremd_tot(Nprocs) !(maxprocs)
92 integer :: iremd_acc_usa(Nprocs),iremd_tot_usa(Nprocs) !(maxprocs)
93 integer :: rstcount !el ilen,
95 character(len=50) :: tytul
98 !old integer nup(0:maxprocs),ndown(0:maxprocs)
99 integer :: rep2i(0:Nprocs),ireqi(Nprocs) !(maxprocs)
100 integer :: icache_all(Nprocs) !(maxprocs)
101 integer :: status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,Nprocs)! (MPI_STATUS_SIZE,maxprocs)
102 logical :: synflag, end_of_run, file_exist = .false.!, ovrtim
104 real(kind=8) :: delta,time00,time01,time001,time02,time03,time04,&
105 time05,time06,time07,time08,tt0,scalfac,ene_iex_iex,&
106 ene_i_i,ene_iex_i,ene_i_iex,xxx,tmp,econstr_temp_i,&
108 integer :: k,il,il1,i,j,nharp,ii,ierr,itime_master,irr,iex,&
109 i_set_temp,itmp,i_temp,i_mult,i_iset,i_mset,i_dir,i_temp1,&
110 i_mult1,i_iset1,i_mset1,ierror
111 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
112 !deb imin_itime_old=0
114 WRITE(iout,*) "JUST AFTER CALL"
115 ! if (.not.allocated(remd_ene)) allocate(remd_ene(0:n_ene+4,Nprocs))
122 if(me.eq.king.or..not.out1file) then
123 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
124 write (iout,*) "NREP=",nrep
128 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
129 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
131 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
133 print *,'MREMD',nodes
134 print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
135 write (iout,*) "Start MREMD: me",me," t_bath",t_bath
136 print *,"NSET=",nset, "MSET=", mset
140 do il1=1,max0(mset(il),1)
150 i_index(i,j,il,il1)=k
156 if(me.eq.king.or..not.out1file) then
157 write(iout,*) (i2rep(i),i=0,nodes-1)
158 write(iout,*) (i2set(i),i=0,nodes-1)
163 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
170 ! print *,'i2rep',me,i2rep(me)
171 ! print *,'rep2i',(rep2i(i),i=0,nrep)
173 !old if (i2rep(me).eq.nrep) then
176 !old nup(0)=remd_m(i2rep(me)+1)
177 !old k=rep2i(int(i2rep(me)))+1
184 !d print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
186 !old if (i2rep(me).eq.1) then
189 !old ndown(0)=remd_m(i2rep(me)-1)
190 !old k=rep2i(i2rep(me)-2)+1
197 !d print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
199 !el common /przechowalnia/
200 if(.not.allocated(d_restart1)) allocate(d_restart1(3,0:(nres2+1)*nodes))
201 if(.not.allocated(d_restart2)) allocate(d_restart2(3,0:(nres2+1)*nodes))
202 if(.not.allocated(p_c)) allocate(p_c(3,(nres2+2)*nodes))
205 write (*,*) "Processor",me," rest",rest,&
206 "restart1fie",restart1file
207 if(rest.and.restart1file) then
209 inquire(file=mremd_rst_name,exist=file_exist)
210 !d write (*,*) me," Before broadcast: file_exist",file_exist
211 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
213 !d write (*,*) me," After broadcast: file_exist",file_exist
215 if(me.eq.king.or..not.out1file) &
216 write (iout,*) 'Master is reading restart1file'
217 call read1restart(i_index)
219 if(me.eq.king.or..not.out1file) &
220 write (iout,*) 'WARNING : no restart1file'
223 if(me.eq.king.or..not.out1file) then
224 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
225 write(iout,*) "i_index"
230 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
239 if (rest.and..not.restart1file) &
240 inquire(file=mremd_rst_name,exist=file_exist)
241 if(.not.file_exist.and.rest.and..not.restart1file) &
242 write(iout,*) 'WARNING : no restart file',mremd_rst_name
243 IF (rest.and.file_exist.and..not.restart1file) THEN
244 write (iout,*) 'Master is reading restart file',&
246 open(irest2,file=mremd_rst_name,status='unknown')
248 read (irest2,*) (i2rep(i),i=0,nodes-1)
250 read (irest2,*) (ifirst(i),i=1,remd_m(1))
253 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
255 read (irest2,*) ndowna(0,il),&
256 (ndowna(i,il),i=1,ndowna(0,il))
262 read (irest2,*) (mset(i),i=1,nset)
264 read (irest2,*) (i2set(i),i=0,nodes-1)
269 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
274 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
275 write(iout,*) "i_index"
280 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
289 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
290 write (iout,'(a6,1000i5)') "ifirst",&
291 (ifirst(i),i=1,remd_m(1))
293 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
294 (nupa(i,il),i=1,nupa(0,il))
295 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
296 (ndowna(i,il),i=1,ndowna(0,il))
298 ELSE IF (.not.(rest.and.file_exist)) THEN
304 if (i2rep(il-1).eq.nrep) then
307 nupa(0,il)=remd_m(i2rep(il-1)+1)
308 k=rep2i(int(i2rep(il-1)))+1
314 if (i2rep(il-1).eq.1) then
317 ndowna(0,il)=remd_m(i2rep(il-1)-1)
318 k=rep2i(i2rep(il-1)-2)+1
326 write (iout,'(a6,100i4)') "ifirst",&
327 (ifirst(i),i=1,remd_m(1))
329 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
330 (nupa(i,il),i=1,nupa(0,il))
331 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
332 (ndowna(i,il),i=1,ndowna(0,il))
338 ! t_bath=retmin+(retmax-retmin)*me/(nodes-1)
339 if(.not.(rest.and.file_exist.and.restart1file)) then
340 if (me .eq. king) then
343 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
345 !d print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
346 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
351 if(me.eq.king.or..not.out1file) &
352 write(iout,*) me,"iset=",iset,"t_bath=",t_bath
356 stdfp(i)=dsqrt(2*Rb*t_bath/d_time)
360 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
364 ! print *,'irep',me,t_bath
366 if (me.eq.king .or. .not. out1file) &
367 write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
368 call rescale_weights(t_bath)
372 !------copy MD--------------
373 ! The driver for molecular dynamics subroutines
374 !------------------------------------------------
380 if(me.eq.king.or..not.out1file) &
381 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
387 ! Determine the inverse of the inertia matrix.
388 call setup_MD_matrices
392 if (me.eq.king .or. .not. out1file) &
393 write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
395 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
397 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
400 call rescale_weights(t_bath)
404 t_MDsetup = MPI_Wtime()-tt0
406 t_MDsetup = tcpu()-tt0
409 ! Entering the MD loop
415 if (lang.eq.2 .or. lang.eq.3) then
419 call sd_verlet_p_setup
421 call sd_verlet_ciccotti_setup
425 pfric0_mat(i,j,0)=pfric_mat(i,j)
426 afric0_mat(i,j,0)=afric_mat(i,j)
427 vfric0_mat(i,j,0)=vfric_mat(i,j)
428 prand0_mat(i,j,0)=prand_mat(i,j)
429 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
430 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
435 flag_stoch(i)=.false.
439 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
441 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
445 else if (lang.eq.1 .or. lang.eq.4) then
449 if (me.eq.king .or. .not. out1file) &
450 write(iout,*) 'Setup time',time00-walltime
453 t_langsetup=MPI_Wtime()-tt0
456 t_langsetup=tcpu()-tt0
462 do while(.not.end_of_run)
464 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
465 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
467 if (lang.gt.0 .and. surfarea .and. &
468 mod(itime,reset_fricmat).eq.0) then
469 if (lang.eq.2 .or. lang.eq.3) then
473 call sd_verlet_p_setup
475 call sd_verlet_ciccotti_setup
479 pfric0_mat(i,j,0)=pfric_mat(i,j)
480 afric0_mat(i,j,0)=afric_mat(i,j)
481 vfric0_mat(i,j,0)=vfric_mat(i,j)
482 prand0_mat(i,j,0)=prand_mat(i,j)
483 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
484 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
489 flag_stoch(i)=.false.
492 else if (lang.eq.1 .or. lang.eq.4) then
495 write (iout,'(a,i10)') &
496 "Friction matrix reset based on surface area, itime",itime
498 if (reset_vel .and. tbf .and. lang.eq.0 &
499 .and. mod(itime,count_reset_vel).eq.0) then
501 if (me.eq.king .or. .not. out1file) &
502 write(iout,'(a,f20.2)') &
503 "Velocities reset to random values, time",totT
506 d_t_old(j,i)=d_t(j,i)
510 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
514 d_t(j,0)=d_t(j,0)-vcm(j)
517 kinetic_T=2.0d0/(dimen3*Rb)*EK
518 scalfac=dsqrt(T_bath/kinetic_T)
519 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
522 d_t_old(j,i)=scalfac*d_t(j,i)
528 ! Time-reversible RESPA algorithm
529 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
530 call RESPA_step(itime)
532 ! Variable time step algorithm.
533 call velverlet_step(itime)
537 call brown_step(itime)
539 print *,"Brown dynamics not here!"
541 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
547 if (mod(itime,ntwe).eq.0) then
549 call enerprint(potEcomp)
552 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
553 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
555 call hairpin(.true.,nharp,iharp)
556 call secondary2(.true.)
557 call pdbout(potE,tytul,ipdb)
562 if (mod(itime,ntwx).eq.0.and.traj1file) then
563 if(ntwx_cache.lt.max_cache_traj_use) then
564 ntwx_cache=ntwx_cache+1
566 if (max_cache_traj_use.ne.1) &
567 print *,itime,"processor ",me," over cache ",ntwx_cache
570 totT_cache(i)=totT_cache(i+1)
571 EK_cache(i)=EK_cache(i+1)
572 potE_cache(i)=potE_cache(i+1)
573 t_bath_cache(i)=t_bath_cache(i+1)
574 Uconst_cache(i)=Uconst_cache(i+1)
575 iset_cache(i)=iset_cache(i+1)
578 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
581 qpair_cache(ii,i)=qpair_cache(ii,i+1)
584 utheta_cache(ii,i)=utheta_cache(ii,i+1)
585 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
586 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
592 c_cache(j,ii,i)=c_cache(j,ii,i+1)
598 totT_cache(ntwx_cache)=totT
599 EK_cache(ntwx_cache)=EK
600 potE_cache(ntwx_cache)=potE
601 t_bath_cache(ntwx_cache)=t_bath
602 Uconst_cache(ntwx_cache)=Uconst
603 iset_cache(ntwx_cache)=iset
606 qfrag_cache(i,ntwx_cache)=qfrag(i)
609 qpair_cache(i,ntwx_cache)=qpair(i)
612 utheta_cache(i,ntwx_cache)=utheta(i)
613 ugamma_cache(i,ntwx_cache)=ugamma(i)
614 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
619 c_cache(j,i,ntwx_cache)=c(j,i)
624 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
625 .and..not.restart1file) then
628 open(irest1,file=mremd_rst_name,status='unknown')
629 write (irest1,*) "i2rep"
630 write (irest1,*) (i2rep(i),i=0,nodes-1)
631 write (irest1,*) "ifirst"
632 write (irest1,*) (ifirst(i),i=1,remd_m(1))
634 write (irest1,*) "nupa",il
635 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
636 write (irest1,*) "ndowna",il
637 write (irest1,*) ndowna(0,il),&
638 (ndowna(i,il),i=1,ndowna(0,il))
641 write (irest1,*) "nset"
642 write (irest1,*) nset
643 write (irest1,*) "mset"
644 write (irest1,*) (mset(i),i=1,nset)
645 write (irest1,*) "i2set"
646 write (irest1,*) (i2set(i),i=0,nodes-1)
647 write (irest1,*) "i_index"
651 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
659 open(irest2,file=rest2name,status='unknown')
660 write(irest2,*) totT,EK,potE,totE,t_bath
662 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
665 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
668 write (irest2,*) iset
675 ! forced synchronization
676 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
677 .and. .not. mremdsync) then
679 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
681 call mpi_recv(itime_master, 1, MPI_INTEGER,&
682 0,101,CG_COMM, status, ierr)
683 call mpi_barrier(CG_COMM, ierr)
684 !deb if (out1file.or.traj1file) then
685 !deb call mpi_gather(itime,1,mpi_integer,
686 !deb & icache_all,1,mpi_integer,king,
689 call mpi_gather(ntwx_cache,1,mpi_integer,&
690 icache_all,1,mpi_integer,king,&
693 write(iout,*) 'REMD synchro at3',itime_master,itime
694 if (itime_master.ge.n_timestep .or. ovrtim()) &
696 !time call flush(iout)
701 if ((mod(itime,nstex).eq.0.and.me.eq.king &
702 .or.end_of_run.and.me.eq.king ) &
703 .and. .not. mremdsync ) then
706 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
707 CG_COMM, ireqi(i), ierr)
708 !d write(iout,*) 'REMD synchro with',i
711 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
712 call mpi_barrier(CG_COMM, ierr)
714 write(iout,*) 'REMD synchro at2',itime,'time=',time01-time00
715 if (out1file.or.traj1file) then
716 !deb call mpi_gather(itime,1,mpi_integer,
717 !deb & itime_all,1,mpi_integer,king,
719 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
720 !deb & (itime_all(i),i=1,nodes)
722 !deb imin_itime=itime_all(1)
724 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
726 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
727 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
728 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
729 call mpi_gather(ntwx_cache,1,mpi_integer,&
730 icache_all,1,mpi_integer,king,&
732 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
733 ! & (icache_all(i),i=1,nodes)
734 ii_write=icache_all(1)
736 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
738 ! write(iout,*) "MIN ii_write=",ii_write
741 !time call flush(iout)
743 if(mremdsync .and. mod(itime,nstex).eq.0) then
745 if (me.eq.king .or. .not. out1file) &
746 write(iout,*) 'REMD synchro at1',itime,ntwx_cache,Nprocs,nodes
747 !!!!!! TRIAL OF MINIM SYNC
750 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
751 CG_COMM, ireqi(i), ierr)
752 !d write(iout,*) 'REMD synchro with',i
755 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
756 call mpi_barrier(CG_COMM, ierr)
760 write(iout,*) icache_all
762 write(iout,*) "before mpi_gather ntwx_cache"
763 call mpi_gather(ntwx_cache,1,mpi_integer,&
764 icache_all(1),1,mpi_integer,king,& ! CONSULT WITH ADAM
766 write(iout,*) "after mpi_gather ntwx_cache"
769 write(iout,'(a19,8000i8)') ' ntwx_cache',&
770 (icache_all(i),i=1,nodes)
771 ii_write=icache_all(1)
773 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
775 write(iout,*) "MIN ii_write=",ii_write
781 ! Update the time safety limiy
782 if (time001-time00.gt.safety) then
783 safety=time001-time00+600
784 write (iout,*) "****** SAFETY increased to",safety," s"
786 if (ovrtim()) end_of_run=.true.
788 if(synflag.and..not.end_of_run) then
792 write(iout,*) 'REMD before',me,t_bath
794 ! call mpi_gather(t_bath,1,mpi_double_precision,
795 ! & remd_t_bath,1,mpi_double_precision,king,
797 potEcomp(n_ene+1)=t_bath
799 potEcomp(n_ene+2)=iset
800 if (iset.lt.nset) then
804 potEcomp(n_ene+3)=Uconst
811 potEcomp(n_ene+4)=Uconst
815 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
816 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
819 call mpi_gather(elow,1,mpi_double_precision,&
820 elowi,1,mpi_double_precision,king,&
822 call mpi_gather(ehigh,1,mpi_double_precision,&
823 ehighi,1,mpi_double_precision,king,&
828 if (me.eq.king .or. .not. out1file) then
829 write(iout,*) 'REMD gather times=',time03-time01 &
833 if (restart1file) call write1rst(i_index)
836 if (me.eq.king .or. .not. out1file) then
837 write(iout,*) 'REMD writing rst time=',time04-time03
840 if (traj1file) call write1traj
842 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
843 !deb & icache_all,1,mpi_integer,king,
845 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
846 !deb & (icache_all(i),i=1,nodes)
851 if (me.eq.king .or. .not. out1file) then
852 write(iout,*) 'REMD writing traj time=',time05-time04
859 remd_t_bath(i)=remd_ene(n_ene+1,i)
860 iremd_iset(i)=remd_ene(n_ene+2,i)
864 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
866 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
870 write(iout,*) 'REMD exchange temp,ene'
872 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
873 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
877 !-------------------------------------
879 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
882 write (iout,*) "remd_m(1)",remd_m(1)
884 i=ifirst(iran_num(1,remd_m(1)))
891 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
893 if(i.gt.0.and.nupa(0,i).gt.0) then
895 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
897 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
899 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
901 ! do while (iex.eq.i)
902 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
903 iex=nupa(iran_num(1,int(nupa(0,i))),i)
905 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
907 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
909 ! Swap temperatures between conformations i and iex with recalculating the free energies
910 ! following temperature changes.
911 ene_iex_iex=remd_ene(0,iex)
912 ene_i_i=remd_ene(0,i)
913 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
914 ! & " iex",iex," ene_iex_iex",ene_iex_iex
915 ! write (iout,*) "rescaling weights with temperature",
918 call rescale_weights(remd_t_bath(i))
920 ! write (iout,*) "0,iex",remd_t_bath(i)
921 ! call enerprint(remd_ene(0,iex))
923 call sum_energy(remd_ene(0,iex),.false.)
924 ene_iex_i=remd_ene(0,iex)
925 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
927 ! write (iout,*) "0,i",remd_t_bath(i)
928 ! call enerprint(remd_ene(0,i))
930 call sum_energy(remd_ene(0,i),.false.)
931 ! write (iout,*) "ene_i_i",remd_ene(0,i)
933 ! write (iout,*) "rescaling weights with temperature",
935 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
936 write (iout,*) "ERROR: inconsistent energies:",i,&
937 ene_i_i,remd_ene(0,i)
939 call rescale_weights(remd_t_bath(iex))
941 ! write (iout,*) "0,i",remd_t_bath(iex)
942 ! call enerprint(remd_ene(0,i))
944 call sum_energy(remd_ene(0,i),.false.)
945 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
947 ene_i_iex=remd_ene(0,i)
949 ! write (iout,*) "0,iex",remd_t_bath(iex)
950 call enerprint(remd_ene(0,iex))
952 call sum_energy(remd_ene(0,iex),.false.)
953 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
954 write (iout,*) "ERROR: inconsistent energies:",iex,&
955 ene_iex_iex,remd_ene(0,iex)
957 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
958 ! write (iout,*) "i",i," iex",iex
959 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
960 ! & " ene_i_iex",ene_i_iex,
961 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
963 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
964 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
966 ! write(iout,*) 'delta',delta
967 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
968 ! & (remd_ene(i)-remd_ene(iex))/Rb/
969 ! & (remd_t_bath(i)*remd_t_bath(iex))
971 if (delta .gt. 50.0d0) then
977 else if (delta.lt.-50.0d0) then
986 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
987 xxx=ran_number(0.0d0,1.0d0)
988 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
990 if (delta .gt. xxx) then
992 remd_t_bath(i)=remd_t_bath(iex)
994 remd_ene(0,i)=ene_i_iex
995 remd_ene(0,iex)=ene_iex_i
1001 ehighi(i)=ehighi(iex)
1008 nupa(k,i)=nupa(k,iex)
1011 ndowna(k,i)=ndowna(k,iex)
1015 if (ifirst(il).eq.i) ifirst(il)=iex
1017 if (nupa(k,il).eq.i) then
1019 elseif (nupa(k,il).eq.iex) then
1024 if (ndowna(k,il).eq.i) then
1026 elseif (ndowna(k,il).eq.iex) then
1032 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1034 i2rep(i-1)=i2rep(iex-1)
1037 ! write(iout,*) 'exchange',i,iex
1038 ! write (iout,'(a8,100i4)') "@ ifirst",
1039 ! & (ifirst(k),k=1,remd_m(1))
1041 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1042 ! & (nupa(k,il),k=1,nupa(0,il))
1043 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1044 ! & (ndowna(k,il),k=1,ndowna(0,il))
1049 remd_ene(0,iex)=ene_iex_iex
1050 remd_ene(0,i)=ene_i_i
1056 !d write (iout,*) "exchange completed"
1060 !d write(iout,*) "########",ii
1062 i_temp=iran_num(1,nrep)
1063 i_mult=iran_num(1,remd_m(i_temp))
1064 i_iset=iran_num(1,nset)
1065 i_mset=iran_num(1,mset(i_iset))
1066 i=i_index(i_temp,i_mult,i_iset,i_mset)
1068 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1071 !d write(iout,*) "i_dir=",i_dir
1073 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1076 i_mult1=iran_num(1,remd_m(i_temp1))
1078 i_mset1=iran_num(1,mset(i_iset1))
1079 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1081 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1084 i_mult1=iran_num(1,remd_m(i_temp1))
1086 i_mset1=iran_num(1,mset(i_iset1))
1087 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1088 econstr_temp_i=remd_ene(20,i)
1089 econstr_temp_iex=remd_ene(20,iex)
1090 remd_ene(20,i)=remd_ene(n_ene+3,i)
1091 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1093 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1096 i_mult1=iran_num(1,remd_m(i_temp1))
1098 i_mset1=iran_num(1,mset(i_iset1))
1099 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1100 econstr_temp_i=remd_ene(20,i)
1101 econstr_temp_iex=remd_ene(20,iex)
1102 remd_ene(20,i)=remd_ene(n_ene+3,i)
1103 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1109 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1112 ! Swap temperatures between conformations i and iex with recalculating the free energies
1113 ! following temperature changes.
1114 ene_iex_iex=remd_ene(0,iex)
1115 ene_i_i=remd_ene(0,i)
1116 !o write (iout,*) "rescaling weights with temperature",
1118 call rescale_weights(remd_t_bath(i))
1120 call sum_energy(remd_ene(0,iex),.false.)
1121 ene_iex_i=remd_ene(0,iex)
1122 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1123 ! call sum_energy(remd_ene(0,i),.false.)
1124 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1125 ! write (iout,*) "rescaling weights with temperature",
1126 ! & remd_t_bath(iex)
1127 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1128 ! write (iout,*) "ERROR: inconsistent energies:",i,
1129 ! & ene_i_i,remd_ene(0,i)
1131 call rescale_weights(remd_t_bath(iex))
1132 call sum_energy(remd_ene(0,i),.false.)
1133 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1134 ene_i_iex=remd_ene(0,i)
1135 ! call sum_energy(remd_ene(0,iex),.false.)
1136 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1137 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1138 ! & ene_iex_iex,remd_ene(0,iex)
1140 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1141 ! write (iout,*) "i",i," iex",iex
1142 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1143 !d & " ene_i_iex",ene_i_iex,
1144 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1145 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1146 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1148 !d write(iout,*) 'delta',delta
1149 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1150 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1151 ! & (remd_t_bath(i)*remd_t_bath(iex))
1152 if (delta .gt. 50.0d0) then
1157 if (i_dir.eq.1.or.i_dir.eq.3) &
1158 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1159 if (i_dir.eq.2.or.i_dir.eq.3) &
1160 iremd_tot_usa(int(i2set(i-1)))= &
1161 iremd_tot_usa(int(i2set(i-1)))+1
1162 xxx=ran_number(0.0d0,1.0d0)
1163 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1164 if (delta .gt. xxx) then
1166 remd_t_bath(i)=remd_t_bath(iex)
1167 remd_t_bath(iex)=tmp
1170 iremd_iset(i)=iremd_iset(iex)
1171 iremd_iset(iex)=itmp
1173 remd_ene(0,i)=ene_i_iex
1174 remd_ene(0,iex)=ene_iex_i
1176 if (i_dir.eq.1.or.i_dir.eq.3) &
1177 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1180 i2rep(i-1)=i2rep(iex-1)
1183 if (i_dir.eq.2.or.i_dir.eq.3) &
1184 iremd_acc_usa(int(i2set(i-1)))= &
1185 iremd_acc_usa(int(i2set(i-1)))+1
1188 i2set(i-1)=i2set(iex-1)
1191 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1192 i_index(i_temp,i_mult,i_iset,i_mset)= &
1193 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1194 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1197 remd_ene(0,iex)=ene_iex_iex
1198 remd_ene(0,i)=ene_i_i
1199 remd_ene(20,iex)=econstr_temp_iex
1200 remd_ene(20,i)=econstr_temp_i
1204 !d do il1=1,mset(il)
1207 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1220 !-------------------------------------
1221 write (iout,*) "NREP",nrep
1223 if(iremd_tot(i).ne.0) &
1224 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1225 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1230 if(iremd_tot_usa(i).ne.0) &
1231 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1232 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1238 !d write (iout,'(a6,100i4)') "ifirst",
1239 !d & (ifirst(i),i=1,remd_m(1))
1241 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1242 !d & (nupa(i,il),i=1,nupa(0,il))
1243 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1244 !d & (ndowna(i,il),i=1,ndowna(0,il))
1249 !d write (iout,*) "Before scatter"
1252 if (me.eq.king) then
1253 write (iout,*) "t_bath before scatter",remd_t_bath
1257 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1258 t_bath,1,mpi_double_precision,king,&
1260 !d write (iout,*) "After scatter"
1263 call mpi_scatter(iremd_iset,1,mpi_integer,&
1264 iset,1,mpi_integer,king,&
1268 if (me.eq.king .or. .not. out1file) then
1269 write(iout,*) 'REMD scatter time=',time07-time06
1273 call mpi_scatter(elowi,1,mpi_double_precision,&
1274 elow,1,mpi_double_precision,king,&
1276 call mpi_scatter(ehighi,1,mpi_double_precision,&
1277 ehigh,1,mpi_double_precision,king,&
1280 call rescale_weights(t_bath)
1281 !o write (iout,*) "Processor",me,
1282 !o & " rescaling weights with temperature",t_bath
1284 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
1286 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
1290 !de write(iout,*) 'REMD after',me,t_bath
1292 if (me.eq.king .or. .not. out1file) then
1293 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1294 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1295 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1296 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1297 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1298 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1299 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1300 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1301 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1307 if (restart1file) then
1308 if (me.eq.king .or. .not. out1file) &
1309 write(iout,*) 'writing restart at the end of run'
1310 call write1rst(i_index)
1313 if (traj1file) call write1traj
1315 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1316 !deb & icache_all,1,mpi_integer,king,
1317 !deb & CG_COMM,ierr)
1318 !deb write(iout,'(a40,8000i8)')
1319 !deb & ' ntwx_cache after traj1file at the end',
1320 !deb & (icache_all(i),i=1,nodes)
1325 t_MD=MPI_Wtime()-tt0
1329 if (me.eq.king .or. .not. out1file) then
1330 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1332 'MD calculations setup:',t_MDsetup,&
1333 'Energy & gradient evaluation:',t_enegrad,&
1334 'Stochastic MD setup:',t_langsetup,&
1335 'Stochastic MD step setup:',t_sdsetup,&
1337 write (iout,'(/28(1h=),a25,27(1h=))') &
1338 ' End of MD calculation '
1340 !el common /przechowalnia/
1341 ! deallocate(d_restart1)
1342 ! deallocate(d_restart2)
1346 end subroutine MREMD
1347 !-----------------------------------------------------------------------------
1348 subroutine write1rst(i_index)
1351 ! implicit real*8 (a-h,o-z)
1352 ! include 'DIMENSIONS'
1354 ! include 'COMMON.MD'
1355 ! include 'COMMON.IOUNITS'
1356 ! include 'COMMON.REMD'
1357 ! include 'COMMON.SETUP'
1358 ! include 'COMMON.CHAIN'
1359 ! include 'COMMON.SBRIDGE'
1360 ! include 'COMMON.INTERACT'
1362 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1363 !el d_restart2(3,2*nres*maxprocs)
1364 real(kind=4) :: r_d(3,0:2*nres)
1365 real(kind=4) :: t5_restart1(5)
1366 integer :: iret,itmp
1367 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1368 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1370 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1371 !el common /przechowalnia/ d_restart1,d_restart2
1372 integer :: i,j,il,il1,ierr,ixdrf
1377 t5_restart1(4)=t_bath
1378 t5_restart1(5)=Uconst
1380 call mpi_gather(t5_restart1,5,mpi_real,&
1381 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1389 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1390 d_restart1,3*2*nres+3,mpi_real,king,&
1402 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1403 d_restart2,3*2*nres+3,mpi_real,king,&
1408 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1410 call xdrfint_(ixdrf, i2rep(i), iret)
1413 call xdrfint_(ixdrf, ifirst(i), iret)
1417 call xdrfint_(ixdrf, nupa(i,il), iret)
1421 call xdrfint_(ixdrf, ndowna(i,il), iret)
1427 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1434 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1441 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1447 call xdrfint_(ixdrf, nset, iret)
1449 call xdrfint_(ixdrf,mset(i), iret)
1452 call xdrfint_(ixdrf,i2set(i), iret)
1458 itmp=i_index(i,j,il,il1)
1459 call xdrfint_(ixdrf,itmp, iret)
1466 call xdrfclose_(ixdrf, iret)
1468 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1470 call xdrfint(ixdrf, i2rep(i), iret)
1473 call xdrfint(ixdrf, ifirst(i), iret)
1477 call xdrfint(ixdrf, nupa(i,il), iret)
1481 call xdrfint(ixdrf, ndowna(i,il), iret)
1487 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1494 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1501 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1508 call xdrfint(ixdrf, nset, iret)
1510 call xdrfint(ixdrf,mset(i), iret)
1513 call xdrfint(ixdrf,i2set(i), iret)
1519 itmp=i_index(i,j,il,il1)
1520 call xdrfint(ixdrf,itmp, iret)
1527 call xdrfclose(ixdrf, iret)
1531 end subroutine write1rst
1532 !-----------------------------------------------------------------------------
1533 subroutine write1traj
1535 ! implicit real*8 (a-h,o-z)
1536 ! include 'DIMENSIONS'
1538 ! include 'COMMON.MD'
1539 ! include 'COMMON.IOUNITS'
1540 ! include 'COMMON.REMD'
1541 ! include 'COMMON.SETUP'
1542 ! include 'COMMON.CHAIN'
1543 ! include 'COMMON.SBRIDGE'
1544 ! include 'COMMON.INTERACT'
1546 real(kind=4) :: t5_restart1(5)
1547 integer :: iret,itmp
1548 real(kind=4) :: xcoord(3,2*nres+2),prec
1549 real(kind=4) :: r_qfrag(50),r_qpair(100)
1550 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1551 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1552 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1553 p_uscdiff(100*Nprocs) !(100*maxprocs)
1554 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1555 real(kind=4) :: r_c(3,2*nres+2)
1556 !el common /przechowalnia/ p_c
1558 integer :: i,j,il,ierr,ii,ixdrf
1560 call mpi_bcast(ii_write,1,mpi_integer,&
1564 ! print *,'traj1file',me,ii_write,ntwx_cache
1568 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1570 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1573 ! write (iout,*) "before gather write1traj: from node",ii
1575 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1577 t5_restart1(1)=totT_cache(ii)
1578 t5_restart1(2)=EK_cache(ii)
1579 t5_restart1(3)=potE_cache(ii)
1580 t5_restart1(4)=t_bath_cache(ii)
1581 t5_restart1(5)=Uconst_cache(ii)
1582 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1584 call mpi_gather(t5_restart1,5,mpi_real,&
1585 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1587 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1590 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1591 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1594 r_qfrag(i)=qfrag_cache(i,ii)
1597 r_qpair(i)=qpair_cache(i,ii)
1600 r_utheta(i)=utheta_cache(i,ii)
1601 r_ugamma(i)=ugamma_cache(i,ii)
1602 r_uscdiff(i)=uscdiff_cache(i,ii)
1605 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1606 p_qfrag,nfrag,mpi_real,king,&
1608 call mpi_gather(r_qpair,npair,mpi_real,&
1609 p_qpair,npair,mpi_real,king,&
1611 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1612 p_utheta,nfrag_back,mpi_real,king,&
1614 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1615 p_ugamma,nfrag_back,mpi_real,king,&
1617 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1618 p_uscdiff,nfrag_back,mpi_real,king,&
1622 write (iout,*) "p_qfrag"
1624 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1626 write (iout,*) "p_qpair"
1628 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1634 r_c(j,i)=c_cache(j,i,ii)
1638 call mpi_gather(r_c,3*2*nres,mpi_real,&
1639 p_c,3*2*nres,mpi_real,king,&
1645 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1646 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1647 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1648 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1649 call xdrfint_(ixdrf, nss, iret)
1652 call xdrfint(ixdrf, idssb(j)+nres, iret)
1653 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1655 call xdrfint_(ixdrf, ihpb(j), iret)
1656 call xdrfint_(ixdrf, jhpb(j), iret)
1659 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1660 call xdrfint_(ixdrf, iset_restart1(il), iret)
1662 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1665 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1668 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1669 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1670 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1675 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1680 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1684 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1688 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1689 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1690 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1691 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1692 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1693 call xdrfint(ixdrf, nss, iret)
1696 call xdrfint(ixdrf, idssb(j)+nres, iret)
1697 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1699 call xdrfint(ixdrf, ihpb(j), iret)
1700 call xdrfint(ixdrf, jhpb(j), iret)
1703 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1704 call xdrfint(ixdrf, iset_restart1(il), iret)
1706 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1709 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1712 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1713 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1714 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1719 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1724 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1728 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1734 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1736 if(me.eq.king) call xdrfclose(ixdrf, iret)
1738 do i=1,ntwx_cache-ii_write
1740 totT_cache(i)=totT_cache(ii_write+i)
1741 EK_cache(i)=EK_cache(ii_write+i)
1742 potE_cache(i)=potE_cache(ii_write+i)
1743 t_bath_cache(i)=t_bath_cache(ii_write+i)
1744 Uconst_cache(i)=Uconst_cache(ii_write+i)
1745 iset_cache(i)=iset_cache(ii_write+i)
1748 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1751 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1754 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1755 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1756 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1761 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1765 ntwx_cache=ntwx_cache-ii_write
1767 end subroutine write1traj
1768 !-----------------------------------------------------------------------------
1769 subroutine read1restart(i_index)
1771 ! implicit real*8 (a-h,o-z)
1772 ! include 'DIMENSIONS'
1774 ! include 'COMMON.MD'
1775 ! include 'COMMON.IOUNITS'
1776 ! include 'COMMON.REMD'
1777 ! include 'COMMON.SETUP'
1778 ! include 'COMMON.CHAIN'
1779 ! include 'COMMON.SBRIDGE'
1780 ! include 'COMMON.INTERACT'
1781 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1782 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1783 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1784 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1786 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1787 !el common /przechowalnia/ d_restart1
1788 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1790 write (*,*) "Processor",me," called read1restart"
1793 open(irest2,file=mremd_rst_name,status='unknown')
1794 read(irest2,*,err=334) i
1795 write(iout,*) "Reading old rst in ASCI format"
1797 call read1restart_old
1801 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1804 call xdrfint_(ixdrf, i2rep(i), iret)
1807 call xdrfint_(ixdrf, ifirst(i), iret)
1810 call xdrfint_(ixdrf, nupa(0,il), iret)
1812 call xdrfint_(ixdrf, nupa(i,il), iret)
1815 call xdrfint_(ixdrf, ndowna(0,il), iret)
1817 call xdrfint_(ixdrf, ndowna(i,il), iret)
1822 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1826 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1829 call xdrfint(ixdrf, i2rep(i), iret)
1832 call xdrfint(ixdrf, ifirst(i), iret)
1835 call xdrfint(ixdrf, nupa(0,il), iret)
1837 call xdrfint(ixdrf, nupa(i,il), iret)
1840 call xdrfint(ixdrf, ndowna(0,il), iret)
1842 call xdrfint(ixdrf, ndowna(i,il), iret)
1847 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1852 call mpi_scatter(t_restart1,5,mpi_real,&
1853 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1857 t_bath=t5_restart1(4)
1862 ! read(irest2,'(3e15.5)')
1863 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1866 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1868 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1874 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1875 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1885 ! read(irest2,'(3e15.5)')
1886 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1889 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1891 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1897 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1898 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1909 call xdrfint_(ixdrf, nset, iret)
1911 call xdrfint_(ixdrf,mset(i), iret)
1914 call xdrfint_(ixdrf,i2set(i), iret)
1920 call xdrfint_(ixdrf,itmp, iret)
1921 i_index(i,j,il,il1)=itmp
1929 call xdrfint(ixdrf, nset, iret)
1931 call xdrfint(ixdrf,mset(i), iret)
1934 call xdrfint(ixdrf,i2set(i), iret)
1940 call xdrfint(ixdrf,itmp, iret)
1941 i_index(i,j,il,il1)=itmp
1948 call mpi_scatter(i2set,1,mpi_integer,&
1949 iset,1,mpi_integer,king,&
1954 if(me.eq.king) close(irest2)
1956 end subroutine read1restart
1957 !-----------------------------------------------------------------------------
1958 subroutine read1restart_old
1960 ! implicit real*8 (a-h,o-z)
1961 ! include 'DIMENSIONS'
1963 ! include 'COMMON.MD'
1964 ! include 'COMMON.IOUNITS'
1965 ! include 'COMMON.REMD'
1966 ! include 'COMMON.SETUP'
1967 ! include 'COMMON.CHAIN'
1968 ! include 'COMMON.SBRIDGE'
1969 ! include 'COMMON.INTERACT'
1970 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1971 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1972 !el common /przechowalnia/ d_restart1
1974 integer :: i,j,il,ierr
1977 open(irest2,file=mremd_rst_name,status='unknown')
1978 read (irest2,*) (i2rep(i),i=0,nodes-1)
1979 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1981 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1982 read (irest2,*) ndowna(0,il),&
1983 (ndowna(i,il),i=1,ndowna(0,il))
1986 read(irest2,*) (t_restart1(j,il),j=1,4)
1989 call mpi_scatter(t_restart1,5,mpi_real,&
1990 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1994 t_bath=t5_restart1(4)
1999 read(irest2,'(3e15.5)') &
2000 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2004 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2005 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2015 read(irest2,'(3e15.5)') &
2016 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2020 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2021 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2027 if(me.eq.king) close(irest2)
2029 end subroutine read1restart_old
2030 !----------------------------------------------------------------
2031 subroutine alloc_MREMD_arrays
2033 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2034 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
2036 ! common /remdcommon/ in io: read_REMDpar
2037 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2038 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2039 ! common /remdrestart/
2040 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2042 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2043 allocate(ifirst(0:nodes)) !(maxprocs)
2044 allocate(nupa(0:nodes,0:2*nodes))
2045 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2046 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2047 allocate(iset_restart1(nodes)) !(maxprocs)
2048 ! common /traj1cache/
2049 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2050 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2051 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2052 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2053 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2054 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2055 allocate(utheta_cache(nfrag_back,max_cache_traj))
2056 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2057 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2058 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2061 end subroutine alloc_MREMD_arrays
2062 !-----------------------------------------------------------------------------
2063 !-----------------------------------------------------------------------------