+++ /dev/null
- module MREMDyn
-!-----------------------------------------------------------------------------
- use io_units
- use names
- use MPI_data
- use md_data
- use remd_data
- use geometry_data
- use energy_data
- use control_data, only:maxprocs
- use MDyn
-
- implicit none
-!-----------------------------------------------------------------------------
-! commom.remd
-! common /remdrestart/
- integer(kind=2),dimension(:),allocatable :: i2set !(0:maxprocs)
- integer(kind=2),dimension(:),allocatable :: ifirst !(maxprocs)
- integer(kind=2),dimension(:,:),allocatable :: nupa,&
- ndowna !(0:maxprocs/4,0:maxprocs)
- real(kind=4),dimension(:,:),allocatable :: t_restart1 !(5,maxprocs)
- integer,dimension(:),allocatable :: iset_restart1 !(maxprocs)
-! common /traj1cache/
- real(kind=4),dimension(:),allocatable :: totT_cache,EK_cache,&
- potE_cache,t_bath_cache,Uconst_cache !(max_cache_traj)
- real(kind=4),dimension(:,:),allocatable :: qfrag_cache !(50,max_cache_traj)
- real(kind=4),dimension(:,:),allocatable :: qpair_cache !(100,max_cache_traj)
- real(kind=4),dimension(:,:),allocatable :: ugamma_cache,&
- utheta_cache,uscdiff_cache !(maxfrag_back,max_cache_traj)
- real(kind=4),dimension(:,:,:),allocatable :: c_cache !(3,maxres2+2,max_cache_traj)
- integer :: ntwx_cache,ii_write !,max_cache_traj_use
- integer,dimension(:),allocatable :: iset_cache !(max_cache_traj)
-!-----------------------------------------------------------------------------
-! common /przechowalnia/
- real(kind=4),dimension(:,:),allocatable :: d_restart1 !(3,2*nres*maxprocs)
- real(kind=4),dimension(:,:),allocatable :: d_restart2 !(3,2*nres*maxprocs)
- real(kind=4),dimension(:,:),allocatable :: p_c !(3,(nres2+2)*maxprocs)
-!-----------------------------------------------------------------------------
-!
-!
-!-----------------------------------------------------------------------------
- contains
-!-----------------------------------------------------------------------------
-!
-!-----------------------------------------------------------------------------
-! MREMD.F
-!-----------------------------------------------------------------------------
-
- subroutine MREMD
-
- use comm_gucio
- use control, only:tcpu,ovrtim
- use io_base, only:ilen
- use control_data
- use geometry_data
- use random, only: iran_num,ran_number
- use compare, only:hairpin,secondary2
- use io, only:cartout,statout
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- include 'mpif.h'
-! include 'COMMON.CONTROL'
-! include 'COMMON.VAR'
-! include 'COMMON.MD'
-!#ifndef LANG0
-! include 'COMMON.LANGEVIN'
-!#else
-! include 'COMMON.LANGEVIN.lang0'
-!#endif
-! include 'COMMON.CHAIN'
-! include 'COMMON.DERIV'
-! include 'COMMON.GEO'
-! include 'COMMON.LOCAL'
-! include 'COMMON.INTERACT'
-! include 'COMMON.IOUNITS'
-! include 'COMMON.NAMES'
-! include 'COMMON.TIME1'
-! include 'COMMON.REMD'
-! include 'COMMON.SETUP'
-! include 'COMMON.MUCA'
-! include 'COMMON.HAIRPIN'
- integer :: ERRCODE
- real(kind=8),dimension(3) :: L,vcm
- real(kind=8) :: energia(0:n_ene)
- real(kind=8) :: remd_t_bath(maxprocs)
- integer :: iremd_iset(maxprocs)
- integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
- real(kind=8) :: remd_ene(0:n_ene+4,maxprocs)
- integer :: iremd_acc(maxprocs),iremd_tot(maxprocs)
- integer :: iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
- integer :: rstcount !el ilen,
-!el external ilen
- character(len=50) :: tytul
-!el common /gucio/ cm
- integer :: itime
-!old integer nup(0:maxprocs),ndown(0:maxprocs)
- integer :: rep2i(0:maxprocs),ireqi(maxprocs)
- integer :: icache_all(maxprocs)
- integer :: status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
- logical :: synflag, end_of_run, file_exist = .false.!, ovrtim
-
- real(kind=8) :: delta,time00,time01,time001,time02,time03,time04,&
- time05,time06,time07,time08,tt0,scalfac,ene_iex_iex,&
- ene_i_i,ene_iex_i,ene_i_iex,xxx,tmp,econstr_temp_i,&
- econstr_temp_iex
- integer :: k,il,il1,i,j,nharp,ii,ierr,itime_master,irr,iex,&
- i_set_temp,itmp,i_temp,i_mult,i_iset,i_mset,i_dir,i_temp1,&
- i_mult1,i_iset1,i_mset1,ierror
- integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
-!deb imin_itime_old=0
- integer :: nres2 !el
- nres2=2*nres
- time001=0.0d0
-
- ntwx_cache=0
- time00=MPI_WTIME()
- time01=time00
- if(me.eq.king.or..not.out1file) then
- write (iout,*) 'MREMD',nodes,'time before',time00-walltime
- write (iout,*) "NREP=",nrep
- endif
-
- synflag=.false.
- if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
- call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
- endif
- mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
-
-!d print *,'MREMD',nodes
-!d print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
-!de write (iout,*) "Start MREMD: me",me," t_bath",t_bath
- k=0
- rep2i(k)=-1
- do il=1,max0(nset,1)
- do il1=1,max0(mset(il),1)
- do i=1,nrep
- iremd_acc(i)=0
- iremd_acc_usa(i)=0
- iremd_tot(i)=0
- do j=1,remd_m(i)
- i2rep(k)=i
- i2set(k)=il
- rep2i(i)=k
- k=k+1
- i_index(i,j,il,il1)=k
- enddo
- enddo
- enddo
- enddo
-
- if(me.eq.king.or..not.out1file) then
- write(iout,*) (i2rep(i),i=0,nodes-1)
- write(iout,*) (i2set(i),i=0,nodes-1)
- do il=1,nset
- do il1=1,mset(il)
- do i=1,nrep
- do j=1,remd_m(i)
- write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
- enddo
- enddo
- enddo
- enddo
- endif
-
-! print *,'i2rep',me,i2rep(me)
-! print *,'rep2i',(rep2i(i),i=0,nrep)
-
-!old if (i2rep(me).eq.nrep) then
-!old nup(0)=0
-!old else
-!old nup(0)=remd_m(i2rep(me)+1)
-!old k=rep2i(int(i2rep(me)))+1
-!old do i=1,nup(0)
-!old nup(i)=k
-!old k=k+1
-!old enddo
-!old endif
-
-!d print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
-
-!old if (i2rep(me).eq.1) then
-!old ndown(0)=0
-!old else
-!old ndown(0)=remd_m(i2rep(me)-1)
-!old k=rep2i(i2rep(me)-2)+1
-!old do i=1,ndown(0)
-!old ndown(i)=k
-!old k=k+1
-!old enddo
-!old endif
-
-!d print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
-
-!el common /przechowalnia/
- if(.not.allocated(d_restart1)) allocate(d_restart1(3,nres2*nodes))
- if(.not.allocated(d_restart2)) allocate(d_restart2(3,nres2*nodes))
- if(.not.allocated(p_c)) allocate(p_c(3,(nres2+2)*nodes))
-!el-------------
-
- write (*,*) "Processor",me," rest",rest,&
- "restart1fie",restart1file
- if(rest.and.restart1file) then
- if (me.eq.king) &
- inquire(file=mremd_rst_name,exist=file_exist)
-!d write (*,*) me," Before broadcast: file_exist",file_exist
- call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
- IERR)
-!d write (*,*) me," After broadcast: file_exist",file_exist
- if(file_exist) then
- if(me.eq.king.or..not.out1file) &
- write (iout,*) 'Master is reading restart1file'
- call read1restart(i_index)
- else
- if(me.eq.king.or..not.out1file) &
- write (iout,*) 'WARNING : no restart1file'
- endif
-
- if(me.eq.king.or..not.out1file) then
- write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
- write(iout,*) "i_index"
- do il=1,nset
- do il1=1,mset(il)
- do i=1,nrep
- do j=1,remd_m(i)
- write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
- enddo
- enddo
- enddo
- enddo
- endif
- endif
-
- if(me.eq.king) then
- if (rest.and..not.restart1file) &
- inquire(file=mremd_rst_name,exist=file_exist)
- if(.not.file_exist.and.rest.and..not.restart1file) &
- write(iout,*) 'WARNING : no restart file',mremd_rst_name
- IF (rest.and.file_exist.and..not.restart1file) THEN
- write (iout,*) 'Master is reading restart file',&
- mremd_rst_name
- open(irest2,file=mremd_rst_name,status='unknown')
- read (irest2,*)
- read (irest2,*) (i2rep(i),i=0,nodes-1)
- read (irest2,*)
- read (irest2,*) (ifirst(i),i=1,remd_m(1))
- do il=1,nodes
- read (irest2,*)
- read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
- read (irest2,*)
- read (irest2,*) ndowna(0,il),&
- (ndowna(i,il),i=1,ndowna(0,il))
- enddo
- if(usampl) then
- read (irest2,*)
- read (irest2,*) nset
- read (irest2,*)
- read (irest2,*) (mset(i),i=1,nset)
- read (irest2,*)
- read (irest2,*) (i2set(i),i=0,nodes-1)
- read (irest2,*)
- do il=1,nset
- do il1=1,mset(il)
- do i=1,nrep
- read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
- enddo
- enddo
- enddo
-
- write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
- write(iout,*) "i_index"
- do il=1,nset
- do il1=1,mset(il)
- do i=1,nrep
- do j=1,remd_m(i)
- write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
- enddo
- enddo
- enddo
- enddo
- endif
-
- close(irest2)
-
- write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
- write (iout,'(a6,1000i5)') "ifirst",&
- (ifirst(i),i=1,remd_m(1))
- do il=1,nodes
- write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
- (nupa(i,il),i=1,nupa(0,il))
- write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
- (ndowna(i,il),i=1,ndowna(0,il))
- enddo
- ELSE IF (.not.(rest.and.file_exist)) THEN
- do il=1,remd_m(1)
- ifirst(il)=il
- enddo
-
- do il=1,nodes
- if (i2rep(il-1).eq.nrep) then
- nupa(0,il)=0
- else
- nupa(0,il)=remd_m(i2rep(il-1)+1)
- k=rep2i(int(i2rep(il-1)))+1
- do i=1,nupa(0,il)
- nupa(i,il)=k+1
- k=k+1
- enddo
- endif
- if (i2rep(il-1).eq.1) then
- ndowna(0,il)=0
- else
- ndowna(0,il)=remd_m(i2rep(il-1)-1)
- k=rep2i(i2rep(il-1)-2)+1
- do i=1,ndowna(0,il)
- ndowna(i,il)=k+1
- k=k+1
- enddo
- endif
- enddo
-
- write (iout,'(a6,100i4)') "ifirst",&
- (ifirst(i),i=1,remd_m(1))
- do il=1,nodes
- write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
- (nupa(i,il),i=1,nupa(0,il))
- write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
- (ndowna(i,il),i=1,ndowna(0,il))
- enddo
-
- ENDIF
- endif
-!
-! t_bath=retmin+(retmax-retmin)*me/(nodes-1)
- if(.not.(rest.and.file_exist.and.restart1file)) then
- if (me .eq. king) then
- t_bath=retmin
- else
- t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
- endif
-!d print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
- if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
-
- endif
- if(usampl) then
- iset=i2set(me)
- if(me.eq.king.or..not.out1file) &
- write(iout,*) me,"iset=",iset,"t_bath=",t_bath
- endif
-!
- stdfp=dsqrt(2*Rb*t_bath/d_time)
- do i=1,ntyp
- stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
- enddo
-
-! print *,'irep',me,t_bath
- if (.not.rest) then
- if (me.eq.king .or. .not. out1file) &
- write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
- call rescale_weights(t_bath)
- endif
-
-
-!------copy MD--------------
-! The driver for molecular dynamics subroutines
-!------------------------------------------------
- t_MDsetup=0.0d0
- t_langsetup=0.0d0
- t_MD=0.0d0
- t_enegrad=0.0d0
- t_sdsetup=0.0d0
- if(me.eq.king.or..not.out1file) &
- write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
-#ifdef MPI
- tt0 = MPI_Wtime()
-#else
- tt0 = tcpu()
-#endif
-! Determine the inverse of the inertia matrix.
- call setup_MD_matrices
-! Initialize MD
- call init_MD
- if (rest) then
- if (me.eq.king .or. .not. out1file) &
- write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
- stdfp=dsqrt(2*Rb*t_bath/d_time)
- do i=1,ntyp
- stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
- enddo
- call rescale_weights(t_bath)
- endif
-
-#ifdef MPI
- t_MDsetup = MPI_Wtime()-tt0
-#else
- t_MDsetup = tcpu()-tt0
-#endif
- rstcount=0
-! Entering the MD loop
-#ifdef MPI
- tt0 = MPI_Wtime()
-#else
- tt0 = tcpu()
-#endif
- if (lang.eq.2 .or. lang.eq.3) then
-#ifndef LANG0
- call setup_fricmat
- if (lang.eq.2) then
- call sd_verlet_p_setup
- else
- call sd_verlet_ciccotti_setup
- endif
- do i=1,dimen
- do j=1,dimen
- pfric0_mat(i,j,0)=pfric_mat(i,j)
- afric0_mat(i,j,0)=afric_mat(i,j)
- vfric0_mat(i,j,0)=vfric_mat(i,j)
- prand0_mat(i,j,0)=prand_mat(i,j)
- vrand0_mat1(i,j,0)=vrand_mat1(i,j)
- vrand0_mat2(i,j,0)=vrand_mat2(i,j)
- enddo
- enddo
- flag_stoch(0)=.true.
- do i=1,maxflag_stoch
- flag_stoch(i)=.false.
- enddo
-#else
- write (iout,*) &
- "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
-#ifdef MPI
- call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
-#endif
- stop
-#endif
- else if (lang.eq.1 .or. lang.eq.4) then
- call setup_fricmat
- endif
- time00=MPI_WTIME()
- if (me.eq.king .or. .not. out1file) &
- write(iout,*) 'Setup time',time00-walltime
- call flush(iout)
-#ifdef MPI
- t_langsetup=MPI_Wtime()-tt0
- tt0=MPI_Wtime()
-#else
- t_langsetup=tcpu()-tt0
- tt0=tcpu()
-#endif
- itime=0
- end_of_run=.false.
-
- do while(.not.end_of_run)
- itime=itime+1
- if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
- if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
- rstcount=rstcount+1
- if (lang.gt.0 .and. surfarea .and. &
- mod(itime,reset_fricmat).eq.0) then
- if (lang.eq.2 .or. lang.eq.3) then
-#ifndef LANG0
- call setup_fricmat
- if (lang.eq.2) then
- call sd_verlet_p_setup
- else
- call sd_verlet_ciccotti_setup
- endif
- do i=1,dimen
- do j=1,dimen
- pfric0_mat(i,j,0)=pfric_mat(i,j)
- afric0_mat(i,j,0)=afric_mat(i,j)
- vfric0_mat(i,j,0)=vfric_mat(i,j)
- prand0_mat(i,j,0)=prand_mat(i,j)
- vrand0_mat1(i,j,0)=vrand_mat1(i,j)
- vrand0_mat2(i,j,0)=vrand_mat2(i,j)
- enddo
- enddo
- flag_stoch(0)=.true.
- do i=1,maxflag_stoch
- flag_stoch(i)=.false.
- enddo
-#endif
- else if (lang.eq.1 .or. lang.eq.4) then
- call setup_fricmat
- endif
- write (iout,'(a,i10)') &
- "Friction matrix reset based on surface area, itime",itime
- endif
- if (reset_vel .and. tbf .and. lang.eq.0 &
- .and. mod(itime,count_reset_vel).eq.0) then
- call random_vel
- if (me.eq.king .or. .not. out1file) &
- write(iout,'(a,f20.2)') &
- "Velocities reset to random values, time",totT
- do i=0,2*nres
- do j=1,3
- d_t_old(j,i)=d_t(j,i)
- enddo
- enddo
- endif
- if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
- call inertia_tensor
- call vcm_vel(vcm)
- do j=1,3
- d_t(j,0)=d_t(j,0)-vcm(j)
- enddo
- call kinetic(EK)
- kinetic_T=2.0d0/(dimen3*Rb)*EK
- scalfac=dsqrt(T_bath/kinetic_T)
-!d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
- do i=0,2*nres
- do j=1,3
- d_t_old(j,i)=scalfac*d_t(j,i)
- enddo
- enddo
- endif
- if (lang.ne.4) then
- if (RESPA) then
-! Time-reversible RESPA algorithm
-! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
- call RESPA_step(itime)
- else
-! Variable time step algorithm.
- call velverlet_step(itime)
- endif
- else
-#ifdef BROWN
- call brown_step(itime)
-#else
- print *,"Brown dynamics not here!"
-#ifdef MPI
- call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
-#endif
- stop
-#endif
- endif
- if(ntwe.ne.0) then
- if (mod(itime,ntwe).eq.0) call statout(itime)
- endif
- if (mod(itime,ntwx).eq.0.and..not.traj1file) then
- write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
- if(mdpdb) then
- call hairpin(.true.,nharp,iharp)
- call secondary2(.true.)
- call pdbout(potE,tytul,ipdb)
- else
- call cartout(totT)
- endif
- endif
- if (mod(itime,ntwx).eq.0.and.traj1file) then
- if(ntwx_cache.lt.max_cache_traj_use) then
- ntwx_cache=ntwx_cache+1
- else
- if (max_cache_traj_use.ne.1) &
- print *,itime,"processor ",me," over cache ",ntwx_cache
- do i=1,ntwx_cache-1
-
- totT_cache(i)=totT_cache(i+1)
- EK_cache(i)=EK_cache(i+1)
- potE_cache(i)=potE_cache(i+1)
- t_bath_cache(i)=t_bath_cache(i+1)
- Uconst_cache(i)=Uconst_cache(i+1)
- iset_cache(i)=iset_cache(i+1)
-
- do ii=1,nfrag
- qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
- enddo
- do ii=1,npair
- qpair_cache(ii,i)=qpair_cache(ii,i+1)
- enddo
- do ii=1,nfrag_back
- utheta_cache(ii,i)=utheta_cache(ii,i+1)
- ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
- uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
- enddo
-
-
- do ii=1,nres*2
- do j=1,3
- c_cache(j,ii,i)=c_cache(j,ii,i+1)
- enddo
- enddo
- enddo
- endif
-
- totT_cache(ntwx_cache)=totT
- EK_cache(ntwx_cache)=EK
- potE_cache(ntwx_cache)=potE
- t_bath_cache(ntwx_cache)=t_bath
- Uconst_cache(ntwx_cache)=Uconst
- iset_cache(ntwx_cache)=iset
-
- do i=1,nfrag
- qfrag_cache(i,ntwx_cache)=qfrag(i)
- enddo
- do i=1,npair
- qpair_cache(i,ntwx_cache)=qpair(i)
- enddo
- do i=1,nfrag_back
- utheta_cache(i,ntwx_cache)=utheta(i)
- ugamma_cache(i,ntwx_cache)=ugamma(i)
- uscdiff_cache(i,ntwx_cache)=uscdiff(i)
- enddo
-
- do i=1,nres*2
- do j=1,3
- c_cache(j,i,ntwx_cache)=c(j,i)
- enddo
- enddo
-
- endif
- if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
- .and..not.restart1file) then
-
- if(me.eq.king) then
- open(irest1,file=mremd_rst_name,status='unknown')
- write (irest1,*) "i2rep"
- write (irest1,*) (i2rep(i),i=0,nodes-1)
- write (irest1,*) "ifirst"
- write (irest1,*) (ifirst(i),i=1,remd_m(1))
- do il=1,nodes
- write (irest1,*) "nupa",il
- write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
- write (irest1,*) "ndowna",il
- write (irest1,*) ndowna(0,il),&
- (ndowna(i,il),i=1,ndowna(0,il))
- enddo
- if(usampl) then
- write (irest1,*) "nset"
- write (irest1,*) nset
- write (irest1,*) "mset"
- write (irest1,*) (mset(i),i=1,nset)
- write (irest1,*) "i2set"
- write (irest1,*) (i2set(i),i=0,nodes-1)
- write (irest1,*) "i_index"
- do il=1,nset
- do il1=1,mset(il)
- do i=1,nrep
- write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
- enddo
- enddo
- enddo
-
- endif
- close(irest1)
- endif
- open(irest2,file=rest2name,status='unknown')
- write(irest2,*) totT,EK,potE,totE,t_bath
- do i=1,2*nres
- write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
- enddo
- do i=1,2*nres
- write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
- enddo
- if(usampl) then
- write (irest2,*) iset
- endif
- close(irest2)
- rstcount=0
- endif
-
-! REMD - exchange
-! forced synchronization
- if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
- .and. .not. mremdsync) then
- synflag=.false.
- call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
- if (synflag) then
- call mpi_recv(itime_master, 1, MPI_INTEGER,&
- 0,101,CG_COMM, status, ierr)
- call mpi_barrier(CG_COMM, ierr)
-!deb if (out1file.or.traj1file) then
-!deb call mpi_gather(itime,1,mpi_integer,
-!deb & icache_all,1,mpi_integer,king,
-!deb & CG_COMM,ierr)
- if(traj1file) &
- call mpi_gather(ntwx_cache,1,mpi_integer,&
- icache_all,1,mpi_integer,king,&
- CG_COMM,ierr)
- if (.not.out1file) &
- write(iout,*) 'REMD synchro at',itime_master,itime
- if (itime_master.ge.n_timestep .or. ovrtim()) &
- end_of_run=.true.
-!time call flush(iout)
- endif
- endif
-
-! REMD - exchange
- if ((mod(itime,nstex).eq.0.and.me.eq.king &
- .or.end_of_run.and.me.eq.king ) &
- .and. .not. mremdsync ) then
- synflag=.true.
- do i=1,nodes-1
- call mpi_isend(itime,1,MPI_INTEGER,i,101, &
- CG_COMM, ireqi(i), ierr)
-!d write(iout,*) 'REMD synchro with',i
-!d call flush(iout)
- enddo
- call mpi_waitall(nodes-1,ireqi,statusi,ierr)
- call mpi_barrier(CG_COMM, ierr)
- time01=MPI_WTIME()
- write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
- if (out1file.or.traj1file) then
-!deb call mpi_gather(itime,1,mpi_integer,
-!deb & itime_all,1,mpi_integer,king,
-!deb & CG_COMM,ierr)
-!deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
-!deb & (itime_all(i),i=1,nodes)
- if(traj1file) then
-!deb imin_itime=itime_all(1)
-!deb do i=2,nodes
-!deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
-!deb enddo
-!deb ii_write=(imin_itime-imin_itime_old)/ntwx
-!deb imin_itime_old=int(imin_itime/ntwx)*ntwx
-!deb write(iout,*) imin_itime,imin_itime_old,ii_write
- call mpi_gather(ntwx_cache,1,mpi_integer,&
- icache_all,1,mpi_integer,king,&
- CG_COMM,ierr)
-! write(iout,'(a19,8000i8)') ' ntwx_cache',
-! & (icache_all(i),i=1,nodes)
- ii_write=icache_all(1)
- do i=2,nodes
- if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
- enddo
-! write(iout,*) "MIN ii_write=",ii_write
- endif
- endif
-!time call flush(iout)
- endif
- if(mremdsync .and. mod(itime,nstex).eq.0) then
- synflag=.true.
- if (me.eq.king .or. .not. out1file) &
- write(iout,*) 'REMD synchro at',itime
-
- if(traj1file) then
- call mpi_gather(ntwx_cache,1,mpi_integer,&
- icache_all,1,mpi_integer,king,&
- CG_COMM,ierr)
- if (me.eq.king) then
- write(iout,'(a19,8000i8)') ' ntwx_cache',&
- (icache_all(i),i=1,nodes)
- ii_write=icache_all(1)
- do i=2,nodes
- if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
- enddo
- write(iout,*) "MIN ii_write=",ii_write
- endif
- endif
- call flush(iout)
- endif
- if (synflag) then
-! Update the time safety limiy
- if (time001-time00.gt.safety) then
- safety=time001-time00+600
- write (iout,*) "****** SAFETY increased to",safety," s"
- endif
- if (ovrtim()) end_of_run=.true.
- endif
- if(synflag.and..not.end_of_run) then
- time02=MPI_WTIME()
- synflag=.false.
-
- write(iout,*) 'REMD before',me,t_bath
-
-! call mpi_gather(t_bath,1,mpi_double_precision,
-! & remd_t_bath,1,mpi_double_precision,king,
-! & CG_COMM,ierr)
- potEcomp(n_ene+1)=t_bath
- if (usampl) then
- potEcomp(n_ene+2)=iset
- if (iset.lt.nset) then
- i_set_temp=iset
- iset=iset+1
- call EconstrQ
- potEcomp(n_ene+3)=Uconst
- iset=i_set_temp
- endif
- if (iset.gt.1) then
- i_set_temp=iset
- iset=iset-1
- call EconstrQ
- potEcomp(n_ene+4)=Uconst
- iset=i_set_temp
- endif
- endif
- call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
- remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
- CG_COMM,ierr)
- if(lmuca) then
- call mpi_gather(elow,1,mpi_double_precision,&
- elowi,1,mpi_double_precision,king,&
- CG_COMM,ierr)
- call mpi_gather(ehigh,1,mpi_double_precision,&
- ehighi,1,mpi_double_precision,king,&
- CG_COMM,ierr)
- endif
-
- time03=MPI_WTIME()
- if (me.eq.king .or. .not. out1file) then
- write(iout,*) 'REMD gather times=',time03-time01 &
- ,time03-time02
- endif
-
- if (restart1file) call write1rst(i_index)
-
- time04=MPI_WTIME()
- if (me.eq.king .or. .not. out1file) then
- write(iout,*) 'REMD writing rst time=',time04-time03
- endif
-
- if (traj1file) call write1traj
-!d debugging
-!deb call mpi_gather(ntwx_cache,1,mpi_integer,
-!deb & icache_all,1,mpi_integer,king,
-!deb & CG_COMM,ierr)
-!deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
-!deb & (icache_all(i),i=1,nodes)
-!d end
-
-
- time05=MPI_WTIME()
- if (me.eq.king .or. .not. out1file) then
- write(iout,*) 'REMD writing traj time=',time05-time04
- call flush(iout)
- endif
-
-
- if (me.eq.king) then
- do i=1,nodes
- remd_t_bath(i)=remd_ene(n_ene+1,i)
- iremd_iset(i)=remd_ene(n_ene+2,i)
- enddo
-#ifdef DEBUG
- if(lmuca) then
-!o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
- do i=1,nodes
- write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
- elowi(i),ehighi(i)
- enddo
- else
- write(iout,*) 'REMD exchange temp,ene'
- do i=1,nodes
- write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
- write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
- enddo
- endif
-#endif
-!-------------------------------------
- IF(.not.usampl) THEN
- write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
- " nodes",nodes
- call flush(iout)
- write (iout,*) "remd_m(1)",remd_m(1)
- do irr=1,remd_m(1)
- i=ifirst(iran_num(1,remd_m(1)))
- write (iout,*) "i",i
- call flush(iout)
-
- do ii=1,nodes-1
-
-#ifdef DEBUG
- write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
-#endif
- if(i.gt.0.and.nupa(0,i).gt.0) then
- iex=i
-! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
-! write (iout,*)
-! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
-! call flush(iout)
-! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
-! endif
-! do while (iex.eq.i)
-! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
- iex=nupa(iran_num(1,int(nupa(0,i))),i)
-! enddo
-! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
- if (lmuca) then
- call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
- else
-! Swap temperatures between conformations i and iex with recalculating the free energies
-! following temperature changes.
- ene_iex_iex=remd_ene(0,iex)
- ene_i_i=remd_ene(0,i)
-! write (iout,*) "i",i," ene_i_i",ene_i_i,
-! & " iex",iex," ene_iex_iex",ene_iex_iex
-! write (iout,*) "rescaling weights with temperature",
-! & remd_t_bath(i)
-! call flush(iout)
- call rescale_weights(remd_t_bath(i))
-
-! write (iout,*) "0,iex",remd_t_bath(i)
-! call enerprint(remd_ene(0,iex))
-
- call sum_energy(remd_ene(0,iex),.false.)
- ene_iex_i=remd_ene(0,iex)
-! write (iout,*) "ene_iex_i",remd_ene(0,iex)
-
-! write (iout,*) "0,i",remd_t_bath(i)
-! call enerprint(remd_ene(0,i))
-
- call sum_energy(remd_ene(0,i),.false.)
-! write (iout,*) "ene_i_i",remd_ene(0,i)
-! call flush(iout)
-! write (iout,*) "rescaling weights with temperature",
-! & remd_t_bath(iex)
- if (real(ene_i_i).ne.real(remd_ene(0,i))) then
- write (iout,*) "ERROR: inconsistent energies:",i,&
- ene_i_i,remd_ene(0,i)
- endif
- call rescale_weights(remd_t_bath(iex))
-
-! write (iout,*) "0,i",remd_t_bath(iex)
-! call enerprint(remd_ene(0,i))
-
- call sum_energy(remd_ene(0,i),.false.)
-! write (iout,*) "ene_i_iex",remd_ene(0,i)
-! call flush(iout)
- ene_i_iex=remd_ene(0,i)
-
-! write (iout,*) "0,iex",remd_t_bath(iex)
-! call enerprint(remd_ene(0,iex))
-
- call sum_energy(remd_ene(0,iex),.false.)
- if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
- write (iout,*) "ERROR: inconsistent energies:",iex,&
- ene_iex_iex,remd_ene(0,iex)
- endif
-! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
-! write (iout,*) "i",i," iex",iex
-! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
-! & " ene_i_iex",ene_i_iex,
-! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
-! call flush(iout)
- delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
- (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
- delta=-delta
-! write(iout,*) 'delta',delta
-! delta=(remd_t_bath(i)-remd_t_bath(iex))*
-! & (remd_ene(i)-remd_ene(iex))/Rb/
-! & (remd_t_bath(i)*remd_t_bath(iex))
- endif
- if (delta .gt. 50.0d0) then
- delta=0.0d0
- else
-#ifdef OSF
- if(isnan(delta))then
- delta=0.0d0
- else if (delta.lt.-50.0d0) then
- delta=dexp(50.0d0)
- else
- delta=dexp(-delta)
- endif
-#else
- delta=dexp(-delta)
-#endif
- endif
- iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
- xxx=ran_number(0.0d0,1.0d0)
-! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
-! call flush(iout)
- if (delta .gt. xxx) then
- tmp=remd_t_bath(i)
- remd_t_bath(i)=remd_t_bath(iex)
- remd_t_bath(iex)=tmp
- remd_ene(0,i)=ene_i_iex
- remd_ene(0,iex)=ene_iex_i
- if(lmuca) then
- tmp=elowi(i)
- elowi(i)=elowi(iex)
- elowi(iex)=tmp
- tmp=ehighi(i)
- ehighi(i)=ehighi(iex)
- ehighi(iex)=tmp
- endif
-
-
- do k=0,nodes
- itmp=nupa(k,i)
- nupa(k,i)=nupa(k,iex)
- nupa(k,iex)=itmp
- itmp=ndowna(k,i)
- ndowna(k,i)=ndowna(k,iex)
- ndowna(k,iex)=itmp
- enddo
- do il=1,nodes
- if (ifirst(il).eq.i) ifirst(il)=iex
- do k=1,nupa(0,il)
- if (nupa(k,il).eq.i) then
- nupa(k,il)=iex
- elseif (nupa(k,il).eq.iex) then
- nupa(k,il)=i
- endif
- enddo
- do k=1,ndowna(0,il)
- if (ndowna(k,il).eq.i) then
- ndowna(k,il)=iex
- elseif (ndowna(k,il).eq.iex) then
- ndowna(k,il)=i
- endif
- enddo
- enddo
-
- iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
- itmp=i2rep(i-1)
- i2rep(i-1)=i2rep(iex-1)
- i2rep(iex-1)=itmp
-
-! write(iout,*) 'exchange',i,iex
-! write (iout,'(a8,100i4)') "@ ifirst",
-! & (ifirst(k),k=1,remd_m(1))
-! do il=1,nodes
-! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
-! & (nupa(k,il),k=1,nupa(0,il))
-! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
-! & (ndowna(k,il),k=1,ndowna(0,il))
-! enddo
-! call flush(iout)
-
- else
- remd_ene(0,iex)=ene_iex_iex
- remd_ene(0,i)=ene_i_i
- i=iex
- endif
- endif
- enddo
- enddo
-!d write (iout,*) "exchange completed"
-!d call flush(iout)
- ELSE
- do ii=1,nodes
-!d write(iout,*) "########",ii
-
- i_temp=iran_num(1,nrep)
- i_mult=iran_num(1,remd_m(i_temp))
- i_iset=iran_num(1,nset)
- i_mset=iran_num(1,mset(i_iset))
- i=i_index(i_temp,i_mult,i_iset,i_mset)
-
-!d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
-
- i_dir=iran_num(1,3)
-!d write(iout,*) "i_dir=",i_dir
-
- if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
-
- i_temp1=i_temp+1
- i_mult1=iran_num(1,remd_m(i_temp1))
- i_iset1=i_iset
- i_mset1=iran_num(1,mset(i_iset1))
- iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-
- elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
-
- i_temp1=i_temp
- i_mult1=iran_num(1,remd_m(i_temp1))
- i_iset1=i_iset+1
- i_mset1=iran_num(1,mset(i_iset1))
- iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
- econstr_temp_i=remd_ene(20,i)
- econstr_temp_iex=remd_ene(20,iex)
- remd_ene(20,i)=remd_ene(n_ene+3,i)
- remd_ene(20,iex)=remd_ene(n_ene+4,iex)
-
- elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
-
- i_temp1=i_temp+1
- i_mult1=iran_num(1,remd_m(i_temp1))
- i_iset1=i_iset+1
- i_mset1=iran_num(1,mset(i_iset1))
- iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
- econstr_temp_i=remd_ene(20,i)
- econstr_temp_iex=remd_ene(20,iex)
- remd_ene(20,i)=remd_ene(n_ene+3,i)
- remd_ene(20,iex)=remd_ene(n_ene+4,iex)
-
- else
- goto 444
- endif
-
-!d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
- call flush(iout)
-
-! Swap temperatures between conformations i and iex with recalculating the free energies
-! following temperature changes.
- ene_iex_iex=remd_ene(0,iex)
- ene_i_i=remd_ene(0,i)
-!o write (iout,*) "rescaling weights with temperature",
-!o & remd_t_bath(i)
- call rescale_weights(remd_t_bath(i))
-
- call sum_energy(remd_ene(0,iex),.false.)
- ene_iex_i=remd_ene(0,iex)
-!d write (iout,*) "ene_iex_i",remd_ene(0,iex)
-! call sum_energy(remd_ene(0,i),.false.)
-!d write (iout,*) "ene_i_i",remd_ene(0,i)
-! write (iout,*) "rescaling weights with temperature",
-! & remd_t_bath(iex)
-! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
-! write (iout,*) "ERROR: inconsistent energies:",i,
-! & ene_i_i,remd_ene(0,i)
-! endif
- call rescale_weights(remd_t_bath(iex))
- call sum_energy(remd_ene(0,i),.false.)
-!d write (iout,*) "ene_i_iex",remd_ene(0,i)
- ene_i_iex=remd_ene(0,i)
-! call sum_energy(remd_ene(0,iex),.false.)
-! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
-! write (iout,*) "ERROR: inconsistent energies:",iex,
-! & ene_iex_iex,remd_ene(0,iex)
-! endif
-!d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
-! write (iout,*) "i",i," iex",iex
-!d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
-!d & " ene_i_iex",ene_i_iex,
-!d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
- delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
- (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
- delta=-delta
-!d write(iout,*) 'delta',delta
-! delta=(remd_t_bath(i)-remd_t_bath(iex))*
-! & (remd_ene(i)-remd_ene(iex))/Rb/
-! & (remd_t_bath(i)*remd_t_bath(iex))
- if (delta .gt. 50.0d0) then
- delta=0.0d0
- else
- delta=dexp(-delta)
- endif
- if (i_dir.eq.1.or.i_dir.eq.3) &
- iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
- if (i_dir.eq.2.or.i_dir.eq.3) &
- iremd_tot_usa(int(i2set(i-1)))= &
- iremd_tot_usa(int(i2set(i-1)))+1
- xxx=ran_number(0.0d0,1.0d0)
-!d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
- if (delta .gt. xxx) then
- tmp=remd_t_bath(i)
- remd_t_bath(i)=remd_t_bath(iex)
- remd_t_bath(iex)=tmp
-
- itmp=iremd_iset(i)
- iremd_iset(i)=iremd_iset(iex)
- iremd_iset(iex)=itmp
-
- remd_ene(0,i)=ene_i_iex
- remd_ene(0,iex)=ene_iex_i
-
- if (i_dir.eq.1.or.i_dir.eq.3) &
- iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
-
- itmp=i2rep(i-1)
- i2rep(i-1)=i2rep(iex-1)
- i2rep(iex-1)=itmp
-
- if (i_dir.eq.2.or.i_dir.eq.3) &
- iremd_acc_usa(int(i2set(i-1)))= &
- iremd_acc_usa(int(i2set(i-1)))+1
-
- itmp=i2set(i-1)
- i2set(i-1)=i2set(iex-1)
- i2set(iex-1)=itmp
-
- itmp=i_index(i_temp,i_mult,i_iset,i_mset)
- i_index(i_temp,i_mult,i_iset,i_mset)= &
- i_index(i_temp1,i_mult1,i_iset1,i_mset1)
- i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
-
- else
- remd_ene(0,iex)=ene_iex_iex
- remd_ene(0,i)=ene_i_i
- remd_ene(20,iex)=econstr_temp_iex
- remd_ene(20,i)=econstr_temp_i
- endif
-
-!d do il=1,nset
-!d do il1=1,mset(il)
-!d do i=1,nrep
-!d do j=1,remd_m(i)
-!d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
-!d enddo
-!d enddo
-!d enddo
-!d enddo
-
- 444 continue
-
- enddo
-
-
- ENDIF
-
-!-------------------------------------
- write (iout,*) "NREP",nrep
- do i=1,nrep
- if(iremd_tot(i).ne.0) &
- write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
- ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
- enddo
-
- if(usampl) then
- do i=1,nset
- if(iremd_tot_usa(i).ne.0) &
- write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
- iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
- enddo
- endif
-
- call flush(iout)
-
-!d write (iout,'(a6,100i4)') "ifirst",
-!d & (ifirst(i),i=1,remd_m(1))
-!d do il=1,nodes
-!d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
-!d & (nupa(i,il),i=1,nupa(0,il))
-!d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
-!d & (ndowna(i,il),i=1,ndowna(0,il))
-!d enddo
- endif
-
- time06=MPI_WTIME()
-!d write (iout,*) "Before scatter"
-!d call flush(iout)
-#ifdef DEBUG
- if (me.eq.king) then
- write (iout,*) "t_bath before scatter",remd_t_bath
- call flush(iout)
- endif
-#endif
- call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
- t_bath,1,mpi_double_precision,king,&
- CG_COMM,ierr)
-!d write (iout,*) "After scatter"
-!d call flush(iout)
- if(usampl) &
- call mpi_scatter(iremd_iset,1,mpi_integer,&
- iset,1,mpi_integer,king,&
- CG_COMM,ierr)
-
- time07=MPI_WTIME()
- if (me.eq.king .or. .not. out1file) then
- write(iout,*) 'REMD scatter time=',time07-time06
- endif
-
- if(lmuca) then
- call mpi_scatter(elowi,1,mpi_double_precision,&
- elow,1,mpi_double_precision,king,&
- CG_COMM,ierr)
- call mpi_scatter(ehighi,1,mpi_double_precision,&
- ehigh,1,mpi_double_precision,king,&
- CG_COMM,ierr)
- endif
- call rescale_weights(t_bath)
-!o write (iout,*) "Processor",me,
-!o & " rescaling weights with temperature",t_bath
-
- stdfp=dsqrt(2*Rb*t_bath/d_time)
- do i=1,ntyp
- stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
- enddo
-
-!de write(iout,*) 'REMD after',me,t_bath
- time08=MPI_WTIME()
- if (me.eq.king .or. .not. out1file) then
- write(iout,*) 'REMD exchange time 8-0=',time08-time00
- write(iout,*) 'REMD exchange time 8-7=',time08-time07
- write(iout,*) 'REMD exchange time 7-6=',time07-time06
- write(iout,*) 'REMD exchange time 6-5=',time06-time05
- write(iout,*) 'REMD exchange time 5-4=',time05-time04
- write(iout,*) 'REMD exchange time 4-3=',time04-time03
- write(iout,*) 'REMD exchange time 3-2=',time03-time02
- write(iout,*) 'REMD exchange time 2-1=',time02-time01
- write(iout,*) 'REMD exchange time 1-0=',time01-time00
- call flush(iout)
- endif
- endif
- enddo
-
- if (restart1file) then
- if (me.eq.king .or. .not. out1file) &
- write(iout,*) 'writing restart at the end of run'
- call write1rst(i_index)
- endif
-
- if (traj1file) call write1traj
-!d debugging
-!deb call mpi_gather(ntwx_cache,1,mpi_integer,
-!deb & icache_all,1,mpi_integer,king,
-!deb & CG_COMM,ierr)
-!deb write(iout,'(a40,8000i8)')
-!deb & ' ntwx_cache after traj1file at the end',
-!deb & (icache_all(i),i=1,nodes)
-!d end
-
-
-#ifdef MPI
- t_MD=MPI_Wtime()-tt0
-#else
- t_MD=tcpu()-tt0
-#endif
- if (me.eq.king .or. .not. out1file) then
- write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
- ' Timing ',&
- 'MD calculations setup:',t_MDsetup,&
- 'Energy & gradient evaluation:',t_enegrad,&
- 'Stochastic MD setup:',t_langsetup,&
- 'Stochastic MD step setup:',t_sdsetup,&
- 'MD steps:',t_MD
- write (iout,'(/28(1h=),a25,27(1h=))') &
- ' End of MD calculation '
- endif
-!el common /przechowalnia/
-! deallocate(d_restart1)
-! deallocate(d_restart2)
-! deallocate(p_c)
-!el--------------
- return
- end subroutine MREMD
-!-----------------------------------------------------------------------------
- subroutine write1rst(i_index)
-
- use control_data
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- include 'mpif.h'
-! include 'COMMON.MD'
-! include 'COMMON.IOUNITS'
-! include 'COMMON.REMD'
-! include 'COMMON.SETUP'
-! include 'COMMON.CHAIN'
-! include 'COMMON.SBRIDGE'
-! include 'COMMON.INTERACT'
-
-!el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
-!el d_restart2(3,2*nres*maxprocs)
- real(kind=4) :: r_d(3,2*nres)
- real(kind=4) :: t5_restart1(5)
- integer :: iret,itmp
- integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
-!el common /przechowalnia/ d_restart1,d_restart2
- integer :: i,j,il,il1,ierr,ixdrf
-
- t5_restart1(1)=totT
- t5_restart1(2)=EK
- t5_restart1(3)=potE
- t5_restart1(4)=t_bath
- t5_restart1(5)=Uconst
-
- call mpi_gather(t5_restart1,5,mpi_real,&
- t_restart1,5,mpi_real,king,CG_COMM,ierr)
-
-
- do i=1,2*nres
- do j=1,3
- r_d(j,i)=d_t(j,i)
- enddo
- enddo
- call mpi_gather(r_d,3*2*nres,mpi_real,&
- d_restart1,3*2*nres,mpi_real,king,&
- CG_COMM,ierr)
-
-
- do i=1,2*nres
- do j=1,3
- r_d(j,i)=dc(j,i)
- enddo
- enddo
- call mpi_gather(r_d,3*2*nres,mpi_real,&
- d_restart2,3*2*nres,mpi_real,king,&
- CG_COMM,ierr)
-
- if(me.eq.king) then
-#ifdef AIX
- call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
- do i=0,nodes-1
- call xdrfint_(ixdrf, i2rep(i), iret)
- enddo
- do i=1,remd_m(1)
- call xdrfint_(ixdrf, ifirst(i), iret)
- enddo
- do il=1,nodes
- do i=0,nupa(0,il)
- call xdrfint_(ixdrf, nupa(i,il), iret)
- enddo
-
- do i=0,ndowna(0,il)
- call xdrfint_(ixdrf, ndowna(i,il), iret)
- enddo
- enddo
-
- do il=1,nodes
- do j=1,4
- call xdrffloat_(ixdrf, t_restart1(j,il), iret)
- enddo
- enddo
-
- do il=0,nodes-1
- do i=1,2*nres
- do j=1,3
- call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
- enddo
- enddo
- enddo
- do il=0,nodes-1
- do i=1,2*nres
- do j=1,3
- call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
- enddo
- enddo
- enddo
-
- if(usampl) then
- call xdrfint_(ixdrf, nset, iret)
- do i=1,nset
- call xdrfint_(ixdrf,mset(i), iret)
- enddo
- do i=0,nodes-1
- call xdrfint_(ixdrf,i2set(i), iret)
- enddo
- do il=1,nset
- do il1=1,mset(il)
- do i=1,nrep
- do j=1,remd_m(i)
- itmp=i_index(i,j,il,il1)
- call xdrfint_(ixdrf,itmp, iret)
- enddo
- enddo
- enddo
- enddo
-
- endif
- call xdrfclose_(ixdrf, iret)
-#else
- call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
- do i=0,nodes-1
- call xdrfint(ixdrf, i2rep(i), iret)
- enddo
- do i=1,remd_m(1)
- call xdrfint(ixdrf, ifirst(i), iret)
- enddo
- do il=1,nodes
- do i=0,nupa(0,il)
- call xdrfint(ixdrf, nupa(i,il), iret)
- enddo
-
- do i=0,ndowna(0,il)
- call xdrfint(ixdrf, ndowna(i,il), iret)
- enddo
- enddo
-
- do il=1,nodes
- do j=1,4
- call xdrffloat(ixdrf, t_restart1(j,il), iret)
- enddo
- enddo
-
- do il=0,nodes-1
- do i=1,2*nres
- do j=1,3
- call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
- enddo
- enddo
- enddo
- do il=0,nodes-1
- do i=1,2*nres
- do j=1,3
- call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
- enddo
- enddo
- enddo
-
-
- if(usampl) then
- call xdrfint(ixdrf, nset, iret)
- do i=1,nset
- call xdrfint(ixdrf,mset(i), iret)
- enddo
- do i=0,nodes-1
- call xdrfint(ixdrf,i2set(i), iret)
- enddo
- do il=1,nset
- do il1=1,mset(il)
- do i=1,nrep
- do j=1,remd_m(i)
- itmp=i_index(i,j,il,il1)
- call xdrfint(ixdrf,itmp, iret)
- enddo
- enddo
- enddo
- enddo
-
- endif
- call xdrfclose(ixdrf, iret)
-#endif
- endif
- return
- end subroutine write1rst
-!-----------------------------------------------------------------------------
- subroutine write1traj
-
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- include 'mpif.h'
-! include 'COMMON.MD'
-! include 'COMMON.IOUNITS'
-! include 'COMMON.REMD'
-! include 'COMMON.SETUP'
-! include 'COMMON.CHAIN'
-! include 'COMMON.SBRIDGE'
-! include 'COMMON.INTERACT'
-
- real(kind=4) :: t5_restart1(5)
- integer :: iret,itmp
- real(kind=4) :: xcoord(3,2*nres+2),prec
- real(kind=4) :: r_qfrag(50),r_qpair(100)
- real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
- real(kind=4) :: p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
- real(kind=4) :: p_utheta(50*maxprocs),p_ugamma(100*maxprocs),&
- p_uscdiff(100*maxprocs)
-!el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
- real(kind=4) :: r_c(3,2*nres+2)
-!el common /przechowalnia/ p_c
-
- integer :: i,j,il,ierr,ii,ixdrf
-
- call mpi_bcast(ii_write,1,mpi_integer,&
- king,CG_COMM,ierr)
-
-! debugging
- print *,'traj1file',me,ii_write,ntwx_cache
-! end debugging
-
-#ifdef AIX
- if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
-#else
- if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
-#endif
- do ii=1,ii_write
-! write (iout,*) "before gather write1traj: from node",ii
-! call flush(iout)
-! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
-! call flush(iout)
- t5_restart1(1)=totT_cache(ii)
- t5_restart1(2)=EK_cache(ii)
- t5_restart1(3)=potE_cache(ii)
- t5_restart1(4)=t_bath_cache(ii)
- t5_restart1(5)=Uconst_cache(ii)
-! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
- call flush(iout)
- call mpi_gather(t5_restart1,5,mpi_real,&
- t_restart1,5,mpi_real,king,CG_COMM,ierr)
-! do il=1,nodes
-! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
-! enddo
-
- call mpi_gather(iset_cache(ii),1,mpi_integer,&
- iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
-
- do i=1,nfrag
- r_qfrag(i)=qfrag_cache(i,ii)
- enddo
- do i=1,npair
- r_qpair(i)=qpair_cache(i,ii)
- enddo
- do i=1,nfrag_back
- r_utheta(i)=utheta_cache(i,ii)
- r_ugamma(i)=ugamma_cache(i,ii)
- r_uscdiff(i)=uscdiff_cache(i,ii)
- enddo
-
- call mpi_gather(r_qfrag,nfrag,mpi_real,&
- p_qfrag,nfrag,mpi_real,king,&
- CG_COMM,ierr)
- call mpi_gather(r_qpair,npair,mpi_real,&
- p_qpair,npair,mpi_real,king,&
- CG_COMM,ierr)
- call mpi_gather(r_utheta,nfrag_back,mpi_real,&
- p_utheta,nfrag_back,mpi_real,king,&
- CG_COMM,ierr)
- call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
- p_ugamma,nfrag_back,mpi_real,king,&
- CG_COMM,ierr)
- call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
- p_uscdiff,nfrag_back,mpi_real,king,&
- CG_COMM,ierr)
-
-#ifdef DEBUG
- write (iout,*) "p_qfrag"
- do i=1,nodes
- write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
- enddo
- write (iout,*) "p_qpair"
- do i=1,nodes
- write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
- enddo
- call flush(iout)
-#endif
- do i=1,nres*2
- do j=1,3
- r_c(j,i)=c_cache(j,i,ii)
- enddo
- enddo
-
- call mpi_gather(r_c,3*2*nres,mpi_real,&
- p_c,3*2*nres,mpi_real,king,&
- CG_COMM,ierr)
-
- if(me.eq.king) then
-#ifdef AIX
- do il=1,nodes
- call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
- call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
- call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
- call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
- call xdrfint_(ixdrf, nss, iret)
- do j=1,nss
- if (dyn_ss) then
- call xdrfint(ixdrf, idssb(j)+nres, iret)
- call xdrfint(ixdrf, jdssb(j)+nres, iret)
- else
- call xdrfint_(ixdrf, ihpb(j), iret)
- call xdrfint_(ixdrf, jhpb(j), iret)
- endif
- enddo
- call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
- call xdrfint_(ixdrf, iset_restart1(il), iret)
- do i=1,nfrag
- call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
- enddo
- do i=1,npair
- call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
- enddo
- do i=1,nfrag_back
- call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
- call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
- call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
- enddo
- prec=10000.0
- do i=1,nres
- do j=1,3
- xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
- enddo
- enddo
- do i=nnt,nct
- do j=1,3
- xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
- enddo
- enddo
- itmp=nres+nct-nnt+1
- call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
- enddo
-#else
- do il=1,nodes
- call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
- call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
- call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
- call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
-! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
- call xdrfint(ixdrf, nss, iret)
- do j=1,nss
- if (dyn_ss) then
- call xdrfint(ixdrf, idssb(j)+nres, iret)
- call xdrfint(ixdrf, jdssb(j)+nres, iret)
- else
- call xdrfint(ixdrf, ihpb(j), iret)
- call xdrfint(ixdrf, jhpb(j), iret)
- endif
- enddo
- call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
- call xdrfint(ixdrf, iset_restart1(il), iret)
- do i=1,nfrag
- call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
- enddo
- do i=1,npair
- call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
- enddo
- do i=1,nfrag_back
- call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
- call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
- call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
- enddo
- prec=10000.0
- do i=1,nres
- do j=1,3
- xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
- enddo
- enddo
- do i=nnt,nct
- do j=1,3
- xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
- enddo
- enddo
- itmp=nres+nct-nnt+1
- call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
- enddo
-#endif
- endif
- enddo
-#ifdef AIX
- if(me.eq.king) call xdrfclose_(ixdrf, iret)
-#else
- if(me.eq.king) call xdrfclose(ixdrf, iret)
-#endif
- do i=1,ntwx_cache-ii_write
-
- totT_cache(i)=totT_cache(ii_write+i)
- EK_cache(i)=EK_cache(ii_write+i)
- potE_cache(i)=potE_cache(ii_write+i)
- t_bath_cache(i)=t_bath_cache(ii_write+i)
- Uconst_cache(i)=Uconst_cache(ii_write+i)
- iset_cache(i)=iset_cache(ii_write+i)
-
- do ii=1,nfrag
- qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
- enddo
- do ii=1,npair
- qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
- enddo
- do ii=1,nfrag_back
- utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
- ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
- uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
- enddo
-
- do ii=1,nres*2
- do j=1,3
- c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
- enddo
- enddo
- enddo
- ntwx_cache=ntwx_cache-ii_write
- return
- end subroutine write1traj
-!-----------------------------------------------------------------------------
- subroutine read1restart(i_index)
-
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- include 'mpif.h'
-! include 'COMMON.MD'
-! include 'COMMON.IOUNITS'
-! include 'COMMON.REMD'
-! include 'COMMON.SETUP'
-! include 'COMMON.CHAIN'
-! include 'COMMON.SBRIDGE'
-! include 'COMMON.INTERACT'
-!el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
- real(kind=4) :: r_d(3,2*nres),t5_restart1(5)
- integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
-!el common /przechowalnia/ d_restart1
- integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
-
- write (*,*) "Processor",me," called read1restart"
-
- if(me.eq.king)then
- open(irest2,file=mremd_rst_name,status='unknown')
- read(irest2,*,err=334) i
- write(iout,*) "Reading old rst in ASCI format"
- close(irest2)
- call read1restart_old
- return
- 334 continue
-#ifdef AIX
- call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
-
- do i=0,nodes-1
- call xdrfint_(ixdrf, i2rep(i), iret)
- enddo
- do i=1,remd_m(1)
- call xdrfint_(ixdrf, ifirst(i), iret)
- enddo
- do il=1,nodes
- call xdrfint_(ixdrf, nupa(0,il), iret)
- do i=1,nupa(0,il)
- call xdrfint_(ixdrf, nupa(i,il), iret)
- enddo
-
- call xdrfint_(ixdrf, ndowna(0,il), iret)
- do i=1,ndowna(0,il)
- call xdrfint_(ixdrf, ndowna(i,il), iret)
- enddo
- enddo
- do il=1,nodes
- do j=1,4
- call xdrffloat_(ixdrf, t_restart1(j,il), iret)
- enddo
- enddo
-#else
- call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
-
- do i=0,nodes-1
- call xdrfint(ixdrf, i2rep(i), iret)
- enddo
- do i=1,remd_m(1)
- call xdrfint(ixdrf, ifirst(i), iret)
- enddo
- do il=1,nodes
- call xdrfint(ixdrf, nupa(0,il), iret)
- do i=1,nupa(0,il)
- call xdrfint(ixdrf, nupa(i,il), iret)
- enddo
-
- call xdrfint(ixdrf, ndowna(0,il), iret)
- do i=1,ndowna(0,il)
- call xdrfint(ixdrf, ndowna(i,il), iret)
- enddo
- enddo
- do il=1,nodes
- do j=1,4
- call xdrffloat(ixdrf, t_restart1(j,il), iret)
- enddo
- enddo
-#endif
- endif
- call mpi_scatter(t_restart1,5,mpi_real,&
- t5_restart1,5,mpi_real,king,CG_COMM,ierr)
- totT=t5_restart1(1)
- EK=t5_restart1(2)
- potE=t5_restart1(3)
- t_bath=t5_restart1(4)
-
- if(me.eq.king)then
- do il=0,nodes-1
- do i=1,2*nres
-! read(irest2,'(3e15.5)')
-! & (d_restart1(j,i+2*nres*il),j=1,3)
- do j=1,3
-#ifdef AIX
- call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
-#else
- call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
-#endif
- enddo
- enddo
- enddo
- endif
- call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
- r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
-
- do i=1,2*nres
- do j=1,3
- d_t(j,i)=r_d(j,i)
- enddo
- enddo
- if(me.eq.king)then
- do il=0,nodes-1
- do i=1,2*nres
-! read(irest2,'(3e15.5)')
-! & (d_restart1(j,i+2*nres*il),j=1,3)
- do j=1,3
-#ifdef AIX
- call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
-#else
- call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
-#endif
- enddo
- enddo
- enddo
- endif
- call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
- r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
- do i=1,2*nres
- do j=1,3
- dc(j,i)=r_d(j,i)
- enddo
- enddo
-
-
- if(usampl) then
-#ifdef AIX
- if(me.eq.king)then
- call xdrfint_(ixdrf, nset, iret)
- do i=1,nset
- call xdrfint_(ixdrf,mset(i), iret)
- enddo
- do i=0,nodes-1
- call xdrfint_(ixdrf,i2set(i), iret)
- enddo
- do il=1,nset
- do il1=1,mset(il)
- do i=1,nrep
- do j=1,remd_m(i)
- call xdrfint_(ixdrf,itmp, iret)
- i_index(i,j,il,il1)=itmp
- enddo
- enddo
- enddo
- enddo
- endif
-#else
- if(me.eq.king)then
- call xdrfint(ixdrf, nset, iret)
- do i=1,nset
- call xdrfint(ixdrf,mset(i), iret)
- enddo
- do i=0,nodes-1
- call xdrfint(ixdrf,i2set(i), iret)
- enddo
- do il=1,nset
- do il1=1,mset(il)
- do i=1,nrep
- do j=1,remd_m(i)
- call xdrfint(ixdrf,itmp, iret)
- i_index(i,j,il,il1)=itmp
- enddo
- enddo
- enddo
- enddo
- endif
-#endif
- call mpi_scatter(i2set,1,mpi_integer,&
- iset,1,mpi_integer,king,&
- CG_COMM,ierr)
-
- endif
-
- if(me.eq.king) close(irest2)
- return
- end subroutine read1restart
-!-----------------------------------------------------------------------------
- subroutine read1restart_old
-
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- include 'mpif.h'
-! include 'COMMON.MD'
-! include 'COMMON.IOUNITS'
-! include 'COMMON.REMD'
-! include 'COMMON.SETUP'
-! include 'COMMON.CHAIN'
-! include 'COMMON.SBRIDGE'
-! include 'COMMON.INTERACT'
-!el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
- real(kind=4) :: r_d(3,2*nres),t5_restart1(5)
-!el common /przechowalnia/ d_restart1
-
- integer :: i,j,il,ierr
-
- if(me.eq.king)then
- open(irest2,file=mremd_rst_name,status='unknown')
- read (irest2,*) (i2rep(i),i=0,nodes-1)
- read (irest2,*) (ifirst(i),i=1,remd_m(1))
- do il=1,nodes
- read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
- read (irest2,*) ndowna(0,il),&
- (ndowna(i,il),i=1,ndowna(0,il))
- enddo
- do il=1,nodes
- read(irest2,*) (t_restart1(j,il),j=1,4)
- enddo
- endif
- call mpi_scatter(t_restart1,5,mpi_real,&
- t5_restart1,5,mpi_real,king,CG_COMM,ierr)
- totT=t5_restart1(1)
- EK=t5_restart1(2)
- potE=t5_restart1(3)
- t_bath=t5_restart1(4)
-
- if(me.eq.king)then
- do il=0,nodes-1
- do i=1,2*nres
- read(irest2,'(3e15.5)') &
- (d_restart1(j,i+2*nres*il),j=1,3)
- enddo
- enddo
- endif
- call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
- r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
-
- do i=1,2*nres
- do j=1,3
- d_t(j,i)=r_d(j,i)
- enddo
- enddo
- if(me.eq.king)then
- do il=0,nodes-1
- do i=1,2*nres
- read(irest2,'(3e15.5)') &
- (d_restart1(j,i+2*nres*il),j=1,3)
- enddo
- enddo
- endif
- call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
- r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
- do i=1,2*nres
- do j=1,3
- dc(j,i)=r_d(j,i)
- enddo
- enddo
- if(me.eq.king) close(irest2)
- return
- end subroutine read1restart_old
-!----------------------------------------------------------------
- subroutine alloc_MREMD_arrays
-
-! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
- if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1)) !(ntyp1))
-! commom.remd
-! common /remdcommon/ in io: read_REMDpar
-! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
-! integer,dimension(:),allocatable :: remd_m !(maxprocs)
-! common /remdrestart/
- if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
-
- allocate(i2set(0:2*nodes)) !(0:maxprocs)
- allocate(ifirst(0:nodes)) !(maxprocs)
- allocate(nupa(0:nodes,0:2*nodes))
- allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
- allocate(t_restart1(5,nodes)) !(5,maxprocs)
- allocate(iset_restart1(nodes)) !(maxprocs)
-! common /traj1cache/
- allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
- allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
- allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
- allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
- allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
- allocate(ugamma_cache(nfrag_back,max_cache_traj))
- allocate(utheta_cache(nfrag_back,max_cache_traj))
- allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
- allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
- allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
-
- return
- end subroutine alloc_MREMD_arrays
-!-----------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
- end module MREMDyn