[lammps-users] Silicon

Dear all,

I want to compute the total energy of some atoms of Silicon in a cubic
box with periodic boundary conditions . but I think I am wrong because my
results are not as same as some benchmarks.
I want to reach to a specific temperature for example 500K. then I run 1e6
steps in NVT ensemble. then I run 1e6 steps in NVE ensemble for
equilibrating the system and then run 8e6 steps in NVE ensemble again to
compute the variables. Here's my input subscript. Would you please tell me
what errors are in it?
Do I simulate the silicon correctly? Are the mass and unit right?
Any help would be appreciated.

Ali Rajabpour


boundary p p p
units metal

# create atom------------------------------------------------------

atom_style atomic

# create geometry--------------------------------------------------

  lattice diamond 5.43
region box block 0 6 0 6 0 6
create_box 3 box
create_atoms 3 box

  mass 1 28
  mass 2 28
  mass 3 28

#Tersoff potential--------------------------------------------------

newton on
pair_style tersoff
pair_coeff * * ../../../potentials/Si.tersoff Si Si Si

# initial velocities-------------------------------------------------

  variable T0 equal 500
  velocity all create ${T0} 482748 mom yes rot yes dist gaussian

# Create thermostat------------------------------------------------

   fix 1 all nvt 500 500 0.01 drag 0.01
   timestep 0.1e-3
   run 1e6
   unfix 1

# Run to equilibrate -----------------------------------------------
   fix 2 all nve
   timestep 0.1e-3
   run 1e6
   unfix 2

# Compute thermophysical variables.------------------------

  compute pot all pe/atom
  compute kin all ke/atom

   variable energy atom c_kin[]+c_pot[]

# Run-----------------------------------------------------

  timestep 0.1e-3
  dump 1 all custom 100 data.txt v_energy
  run 8e6

I don't see where you define a time integration fix (e.g. fix nve)
for the final run. Other than that your script looks fine.