+++ /dev/null
- subroutine MREMD
- 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'
- double precision cm(3),L(3),vcm(3)
- double precision energia(0:n_ene)
- double precision remd_t_bath(maxprocs)
- double precision remd_ene(0:n_ene+1,maxprocs)
- integer iremd_acc(maxprocs),iremd_tot(maxprocs)
- integer ilen,rstcount
- external ilen
- character*50 tytul
- common /gucio/ cm
- integer itime
- integer nup(0:maxprocs),ndown(0:maxprocs)
- integer rep2i(0:maxprocs)
- integer itime_all(maxprocs)
- integer status(MPI_STATUS_SIZE)
- logical synflag,end_of_run,file_exist
-
- time00=MPI_WTIME()
- if(me.eq.king.or..not.out1file)
- & write (iout,*) 'MREMD',nodes,'time before',time00-walltime
-
- synflag=.false.
- mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
-
-cd print *,'MREMD',nodes
-
-cd print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
-
-cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath
- k=0
- rep2i(k)=-1
- do i=1,nrep
- iremd_acc(i)=0
- iremd_tot(i)=0
- do j=1,remd_m(i)
- i2rep(k)=i
- rep2i(i)=k
- k=k+1
- enddo
- enddo
-
-c print *,'i2rep',me,i2rep(me)
-c print *,'rep2i',(rep2i(i),i=0,nrep)
-
- if (i2rep(me).eq.nrep) then
- nup(0)=0
- else
- nup(0)=remd_m(i2rep(me)+1)
- k=rep2i(i2rep(me))+1
- do i=1,nup(0)
- nup(i)=k
- k=k+1
- enddo
- endif
-
-cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
-
- if (i2rep(me).eq.1) then
- ndown(0)=0
- else
- ndown(0)=remd_m(i2rep(me)-1)
- k=rep2i(i2rep(me)-2)+1
- do i=1,ndown(0)
- ndown(i)=k
- k=k+1
- enddo
- endif
-
-cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
-
-
- if(rest.and.restart1file) then
- inquire(file=mremd_rst_name,exist=file_exist)
- if(file_exist) call read1restart
- endif
-
- if(me.eq.king) then
- if (rest.and..not.restart1file)
- & inquire(file=mremd_rst_name,exist=file_exist)
- IF (rest.and.file_exist.and..not.restart1file) THEN
- 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
- 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
- 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(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
-
-c write (iout,'(a6,100i4)') "ifirst",
-c & (ifirst(i),i=1,remd_m(1))
-c do il=1,nodes
-c write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
-c & (nupa(i,il),i=1,nupa(0,il))
-c write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
-c & (ndowna(i,il),i=1,ndowna(0,il))
-c enddo
-
- ENDIF
- endif
-c
-c 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
-cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
- if (remd_tlist) t_bath=remd_t(i2rep(me))
- endif
-c
- stdfp=dsqrt(2*Rb*t_bath/d_time)
- do i=1,ntyp
- stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
- enddo
-
-c 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
-
-
-c------copy MD--------------
-c The driver for molecular dynamics subroutines
-c------------------------------------------------
- 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"
- tt0 = tcpu()
-c Determine the inverse of the inertia matrix.
- call setup_MD_matrices
-c 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
-
- t_MDsetup = tcpu()-tt0
- rstcount=0
-c Entering the MD loop
- tt0 = tcpu()
- if (lang.eq.2 .or. lang.eq.3) then
- call setup_fricmat
- if (lang.eq.2) then
- call sd_verlet_p_setup
- else
- call sd_verlet_ciccotti_setup
- endif
-#ifndef LANG0
- 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
-#endif
- flag_stoch(0)=.true.
- do i=1,maxflag_stoch
- flag_stoch(i)=.false.
- enddo
- 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)
- t_langsetup=tcpu()-tt0
- tt0=tcpu()
- 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
- call setup_fricmat
- if (lang.eq.2) then
- call sd_verlet_p_setup
- else
- call sd_verlet_ciccotti_setup
- endif
-#ifndef LANG0
- 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
-#endif
- flag_stoch(0)=.true.
- do i=1,maxflag_stoch
- flag_stoch(i)=.false.
- enddo
- 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/(dimen*Rb)*EK
- scalfac=dsqrt(T_bath/kinetic_T)
-cd 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
-c Time-reversible RESPA algorithm
-c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
- call RESPA_step(itime)
- else
-c Variable time step algorithm.
- call velverlet_step(itime)
- endif
- else
- call brown_step(itime)
- 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 pdbout(potE,tytul,ipdb)
- else
- call cartout(totT)
- endif
- 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
- 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
- close(irest2)
- rstcount=0
- endif
-
-c REMD - exchange
-c forced synchronization
- if (mod(itime,100).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)
- if (.not. out1file) then
- write(iout,*) 'REMD synchro at',itime_master,itime
- else
- call mpi_gather(itime,1,mpi_integer,
- & itime_all,1,mpi_integer,king,
- & CG_COMM,ierr)
- endif
- if(itime_master.ge.n_timestep) end_of_run=.true.
- call flush(iout)
- endif
- endif
-
-c 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.
- time00=MPI_WTIME()
- do i=1,nodes-1
- call mpi_isend(itime,1,MPI_INTEGER,i,101,
- & CG_COMM, ireq, ierr)
-cd write(iout,*) 'REMD synchro with',i
-cd call flush(iout)
- enddo
- time02=MPI_WTIME()
- write(iout,*) 'REMD synchro at',itime,'time=',time02-time00
- if (out1file) then
- call mpi_gather(itime,1,mpi_integer,
- & itime_all,1,mpi_integer,king,
- & CG_COMM,ierr)
- write(iout,'(a18,8000i8)') 'REMD synchro itime',
- & (itime_all(i),i=1,nodes)
- endif
- 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
- call flush(iout)
- endif
- if(synflag.and..not.end_of_run) then
- time00=MPI_WTIME()
- synflag=.false.
-
-cd write(iout,*) 'REMD before',me,t_bath
-
-c call mpi_gather(t_bath,1,mpi_double_precision,
-c & remd_t_bath,1,mpi_double_precision,king,
-c & CG_COMM,ierr)
- potEcomp(n_ene+1)=t_bath
- call mpi_gather(potEcomp(0),n_ene+2,mpi_double_precision,
- & remd_ene(0,1),n_ene+2,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
-
- if (restart1file) call write1rst
- if (traj1file) call write1traj
-
- if (me.eq.king) then
- do i=1,nodes
- remd_t_bath(i)=remd_ene(n_ene+1,i)
- enddo
- if(lmuca) then
-co 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
-cd write(iout,*) 'REMD exchange temp,ene'
-c do i=1,nodes
-co write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
-c write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
-c enddo
- endif
-c-------------------------------------
- do irr=1,remd_m(1)
- i=ifirst(iran_num(1,remd_m(1)))
- do ii=1,nodes-1
-
- if(i.gt.0.and.nupa(0,i).gt.0) then
- iex=nupa(iran_num(1,int(nupa(0,i))),i)
- if (lmuca) then
- call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
- else
-c Swap temperatures between conformations i and iex with recalculating the free energies
-c following temperature changes.
- ene_iex_iex=remd_ene(0,iex)
- ene_i_i=remd_ene(0,i)
-co write (iout,*) "rescaling weights with temperature",
-co & 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)
-cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
- call sum_energy(remd_ene(0,i),.false.)
-cd write (iout,*) "ene_i_i",remd_ene(0,i)
-c write (iout,*) "rescaling weights with temperature",
-c & 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.)
-c 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
-c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
-c write (iout,*) "i",i," iex",iex
-co write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
-co & " ene_i_iex",ene_i_iex,
-co & " 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
-c write(iout,*) 'delta',delta
-c delta=(remd_t_bath(i)-remd_t_bath(iex))*
-c & (remd_ene(i)-remd_ene(iex))/Rb/
-c & (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(i2rep(i-1))=iremd_tot(i2rep(i-1))+1
- xxx=ran_number(0.0d0,1.0d0)
-co 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
- 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(i2rep(i-1))=iremd_acc(i2rep(i-1))+1
- itmp=i2rep(i-1)
- i2rep(i-1)=i2rep(iex-1)
- i2rep(iex-1)=itmp
-
-cd write(iout,*) 'exchange',i,iex
-cd write (iout,'(a8,100i4)') "@ ifirst",
-cd & (ifirst(k),k=1,remd_m(1))
-cd do il=1,nodes
-cd write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
-cd & (nupa(k,il),k=1,nupa(0,il))
-cd write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
-cd & (ndowna(k,il),k=1,ndowna(0,il))
-cd enddo
-
- else
- remd_ene(0,iex)=ene_iex_iex
- remd_ene(0,i)=ene_i_i
- i=iex
- endif
- endif
- enddo
- enddo
-c-------------------------------------
- do i=1,nrep
- if(iremd_tot(i).ne.0)
- & write(iout,'(a3,i4,2f12.5)') 'ACC',i,remd_t(i)
- & ,iremd_acc(i)/(1.0*iremd_tot(i))
- enddo
-
-cd write (iout,'(a6,100i4)') "ifirst",
-cd & (ifirst(i),i=1,remd_m(1))
-cd do il=1,nodes
-cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
-cd & (nupa(i,il),i=1,nupa(0,il))
-cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
-cd & (ndowna(i,il),i=1,ndowna(0,il))
-cd enddo
- endif
-
- call mpi_scatter(remd_t_bath,1,mpi_double_precision,
- & t_bath,1,mpi_double_precision,king,
- & CG_COMM,ierr)
- 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)
-co write (iout,*) "Processor",me,
-co & " 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
-
-cde write(iout,*) 'REMD after',me,t_bath
- time02=MPI_WTIME()
- if (me.eq.king .or. .not. out1file) then
- write(iout,*) 'REMD exchange time=',time02-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
- endif
-
- if (traj1file) call write1traj
-
-
- t_MD=tcpu()-tt0
- 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
- return
- end
-
- subroutine write1rst
- 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 d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
- & d_restart2(3,2*maxres*maxprocs)
- real t5_restart1(5)
- integer iret,itmp
-
-
- 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
- open(irest1,file=mremd_rst_name,status='unknown')
- write (irest1,*) (i2rep(i),i=0,nodes-1)
- write (irest1,*) (ifirst(i),i=1,remd_m(1))
- do il=1,nodes
- write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
- write (irest1,*) ndowna(0,il),
- & (ndowna(i,il),i=1,ndowna(0,il))
- enddo
-
- do il=1,nodes
- write(irest1,*) (t_restart1(j,il),j=1,4)
- enddo
-
- do il=0,nodes-1
- do i=1,2*nres
- write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
- enddo
- enddo
- do il=0,nodes-1
- do i=1,2*nres
- write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
- enddo
- enddo
- close(irest1)
- endif
-
-
- return
- end
-
- 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 t5_restart1(5)
- integer iret,itmp
- real xcoord(3,maxres2+2),prec
- real r_qfrag(50),r_qpair(100)
- real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
- real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
-
- if(.not.restart1file) then
- 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)
- endif
-
- do i=1,nfrag
- r_qfrag(i)=qfrag(i)
- enddo
- do i=1,npair
- r_qpair(i)=qpair(i)
- enddo
-
- call mpi_gather(r_qfrag,nfrag,mpi_real,
- & p_qfrag,nfrag,mpi_real,king,
- & CG_COMM,ierr)
- call mpi_gather(qpair,nfrag,mpi_real,
- & p_qpair,nfrag,mpi_real,king,
- & CG_COMM,ierr)
-
- do i=1,nres*2
- do j=1,3
- r_c(j,i)=c(j,i)
- 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
- call xdrfopen(ixdrf,cartname, "a", iret)
- 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
- call xdrfint(ixdrf, ihpb(j), iret)
- call xdrfint(ixdrf, jhpb(j), iret)
- enddo
- call xdrfint(ixdrf, nfrag+npair, 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
- 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
- call xdrfclose(ixdrf, iret)
- endif
-
- return
- end
-
-
- subroutine read1restart
- 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 d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
- & t5_restart1(5)
-
- 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
-