md1.f90

分子动力学程序

  ! MD1: a simple, minimal molecular dynamics program in Fortran 90
!
! Lennard-Jones potential, 'velocity' Verlet time integration algorithm.
! Computes kinetic and potential energy, density, pressure.
!
! Furio Ercolessi, SISSA, Trieste, May 1995, revised May 1997
!
! Files used by this program:
!
! Unit  I/O  Meaning
! ----  ---  ----------------------------------------------------------------
!   1    I   Input sample (coordinates, and perhaps also velocities and
!            accelerations) read at the beginning of the run
!   2    O   Output sample (coordinates, velocities, accelerations) written
!            at the end of the run
!   *    I   Standard input for the simulation control
!   *    O   Standard output containing various informations on the run
!
! Output suitable to be directly fed to 'gnuplot' to produce plots as a
! function of time.
 !!
 !!
 !!here,i have do some veifing.
 !!去掉了lis_trim(sampin) and sampout 的类似语句
 !!	其中包含他们的打印语句，读入数据的格式也作改为默认格式了
 !!	去掉了title的输入语句
 !!	增加了一些判断出错在那里的检测语句
 !!数据仍然是用文件的形式读入，自己给它另外建了文件
 !!
 !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  MODULES CONTAINING
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    DATA STRUCTURES

module Particles
!
!  Contains all the data structures containing informations
!  about atom properties
!
   integer, parameter :: DIM=3     ! change DIM to switch between 2D and 3D !
   logical :: VelAcc = .FALSE.     ! velocities and accelerations in file?
   integer :: N=0
   double precision, dimension(DIM) :: BoxSize
!     the following arrays are allocated at run time, 
!     when the number of atoms N is known.
   double precision, dimension(:,:), allocatable :: pos     ! positions
   double precision, dimension(:,:), allocatable :: vel     ! velocities
   double precision, dimension(:,:), allocatable :: acc     ! accelerations
   double precision, dimension(:), allocatable :: ene_pot   ! potential energies
   double precision, dimension(:), allocatable :: ene_kin   ! kinetic energies
!     scalars derived from the above quantities
   double precision :: volume, density                   ! these are constants 
   double precision :: virial      ! virial term, used to compute the pressure
end module Particles

module Simulation_Control
!
!  Contains the parameters controlling the simulation, which will be 
!  supplied as input by the user.
!
  
   character*80     :: SampIn        !  Name of file containing input sample
   character*80     :: SampOut       !  Name of file to contain output sample
   integer          :: Nsteps        !  Number of time steps to do
   double precision :: deltat        !  Time steps (redueced units)
   double precision :: RhoRequested  !  Desired density.  0 leaves it unchanged.
   double precision :: TRequested    !  Desired temperature, or <0 if constant E
   logical          :: ChangeRho     !  .TRUE. when user is changing the density
   logical          :: ConstantT     !  made .TRUE. when (TRequested >= 0)
end module Simulation_Control

module Potential
!     
!     Contains parameters of the Lennard-Jones potential
!
   double precision, parameter :: Rcutoff = 2.5d0   ! cutoff distance
   double precision, parameter :: phicutoff = &     ! potential at cutoff
                                4.d0/(Rcutoff**12) - 4.d0/(Rcutoff**6)
end module Potential

module Statistics
!
!  Contains statistical quantities accumulated during the run.
!  All quantities are resetted to zero.
!
   double precision :: temperature_sum = 0.d0
   double precision :: ene_kin_sum     = 0.d0
   double precision :: ene_pot_sum     = 0.d0
   double precision :: pressure_sum    = 0.d0
end module Statistics

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  CODE STARTS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     HERE

program MD1
!
!  This is the main driver for the molecular dynamics program.
!
implicit none
call Initialize

call Evolve_Sample
call Terminate
end program MD1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  INITIALIZATION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     ROUTINES

subroutine Initialize
!
!  Initialization procedure (called once at the beginning, before the
!  time evolution starts)
!
use Particles
implicit none
!
!  Read the user directives controlling the simulation
!
call Read_Input
!
!  Read the input sample containing the initial particle coordinates
!  (and perhaps also velocities and accelerations)
!
call Read_Sample
!
!  Print informations on the run on the standard output file
!
call Initial_Printout
end subroutine Initialize


subroutine Read_Sample
!
!  Reads the initial sample from file unit 1
!
use Particles
use Simulation_Control
implicit none
double precision, dimension(DIM) :: PosAtomReal,VelAtomReal,AccAtomReal
double precision, dimension(DIM) :: Mass_center
double precision :: scale
integer :: i,k


open(unit=1,file=SampIn,status='old', &
     action='read',err=700)   
read(1,*,end=600,err=900) VelAcc,N,BoxSize
!!
!! add some sentence by myself
print*,velacc,N,boxsize
!!
!!
if ( N <= 0 ) then
   print*,'Read_Sample: FATAL: N is',N
   stop
endif

 
!
!  compute volume and density once for all (they do not change in the run)
!
volume  = product(BoxSize)

density = N / volume
!
!   ... unless the user wants to change the density, in this case we do
!   it here:
!
if (ChangeRho) then
   scale = ( density / RhoRequested ) ** ( 1.d0 / DIM )
   BoxSize = scale*BoxSize
   volume  = product(BoxSize)
   density = N / volume
endif
!
!  now that we know the system size, we can dynamically allocate the
!  arrays containing atomic informations
!
allocate( pos(DIM,N) )
allocate( vel(DIM,N) )
allocate( acc(DIM,N) )
allocate( ene_pot(N) )
allocate( ene_kin(N) )
!
!  read the coordinates from the file (one line per atom), normalize
!  them to the box size along each direction and store them.
!  Energies are set initially to zero.
!
do i=1,N
    read(1,*,end=800,err=900) PosAtomReal
	
    pos(:,i) = PosAtomReal/BoxSize
   ene_pot(i) = 0.d0
   ene_kin(i) = 0.d0
enddo
!
!  For "new" samples (that is, samples just created by defining the atomic
!  coordinates and not the result of previous simulations), we have now
!  read everything, and velocities and accelerations are set to zero.
!  For sample which have been produced by previous simulations, we also
!  have to read velocities and accelerations.
!  The logical variable 'VelAcc' distinguishes between these two cases.
!
if (VelAcc) then     ! also read velocities and accelerations
   do i=1,N
       read(1,*,end=800,err=900) velAtomReal
        vel(:,i) = VelAtomReal/BoxSize
   enddo
   do i=1,N
       read(1,*,end=800,err=900) accAtomReal
	   acc(:,i) = AccAtomReal/BoxSize
   enddo
else                 ! set velocities and accelerations to zero
   vel = 0.d0
   acc = 0.d0
endif
!
!  compute center of mass coordinates
!
Mass_Center = sum( pos , dim=2 ) / N
!
!  translate atoms so that center of mass is at the origin
!
do k=1,DIM
   pos(k,:) = pos(k,:) - Mass_Center(k)
enddo
close(unit=1)
return

!
!  all coordinates read successfully if we get to this point
!
!
!  handling of various kinds of errors
!
600 continue
   print*,'Read_Sample: FATAL: ',SampIn,' is empty?'
   stop
700 continue
   print*,'Read_Sample: FATAL: ',SampIn,' not found.'
   stop

800 continue
   print*,'Read_Sample: FATAL: premature end-of-file at atom ',i
   close(unit=1)
   stop   
900 continue
   print*,'Read_Sample: FATAL: read error in ',SampIn
   close(unit=1)
   stop





end subroutine Read_Sample


subroutine Read_Input
!
!  Read the parameters controlling the simulation from the standard input.
!
use Simulation_Control
implicit none
!
!  Read the input parameters controlling the simulation from the standard
!  input.
!

read(*,'(a)',end=200,err=800) SampIn
read(*,'(a)',end=200,err=800) SampOut
read(*,  *  ,end=200,err=800) Nsteps
read(*,  *  ,end=200,err=800) deltat
read(*,  *  ,end=200,err=800) RhoRequested
read(*,  *  ,end=200,err=800) TRequested
ChangeRho = ( RhoRequested > 0 )
ConstantT = ( TRequested >= 0 )
return   !  successful exit
200 continue
   print*,'Read_Input: FATAL: premature end-of-file in standard input'
   stop
800 continue
   print*,'Read_Input: FATAL: read error in standard input'
   stop
end subroutine Read_Input


subroutine Initial_Printout
!
!  Prints informations on the run parameters on the standard output
!  Leading # are to directly use the output as a gnuplot input.
!


use Particles
use Simulation_Control
implicit none



print '(a,/,a,/,a)','#','# MD1: a minimal molecular dynamics program','#'