Guide to Lammps - Choosing a Timestep

As previously mentioned, choosing a suitable timestep is an important part while setting up a simulation. Luckily, the potential we use would be of use here. Tersoff's potential is essentially a pair potential, with a harmonic behaviour. We can use this fact to experiment with different timesteps and find a value large enough to avoid unneeded computations, and small enough to avoid numeric errors.

The idea is to record a short simulation with a high sampling rate, each step for example, and to plot a trajectory for an atom on a single axis. If the trajectory is still harmonic, the sampling rate is high enough, meaning the timestep is suitable.

Below is an example of a simple simulation doing that.

# --------------------- Initialization -----------------
units metal
atom_style atomic

boundary p p s
newton on

processors * * *

# --------------------- Create Atoms -------------------
lattice diamond 3.57
region diamondBox block 0 5 0 5 0 5 units lattice
create_box 1 diamondBox
create_atoms 1 region diamondBox
mass * 12.0107

# --------------------- Potential Defs -----------------
pair_style tersoff
pair_coeff * * SiC.tersoff C

# --------------------- Simulation Settings ------------
neighbor 2.0 bin
neigh_modify delay 3

timestep 0.00001
thermo_style custom step pe etotal temp
thermo 1

dump diamond all xyz 1 out.diamond.*.xyz

# --------------------- Relaxation ---------------------
min_style cg
minimize 1e-6 1e-6 5000 10000

# --------------------- Dynamics -----------------------
fix dynamics all nve
fix tempControl all temp/rescale 100 300.0 300.0 0.02 1.0
run 1000
unfix tempControl

minimize 1e-6 1e-6 5000 10000

# --------------------- Simulation End -----------------
unfix dynamics
undump diamond

notice we did not use the command velocity to assign initial velocities, and that the thermostat used is temp/rescale. Using Matlab to read the files and produce a plot, we receive the following trajectory for the first atom on the X axis:

This is obviously an odd result, which suggests the atom has no velocity on the x axis. Careful observation of other trajectories suggests this is the case for other atoms, on other axes, as well. To fix this error, we add the following line in the 'Create Atoms' section:

velocity all create 3000.0 4928459 dist gaussian

Notice we assign velocities compatible to 3000 K, so we shall also change the fix line:

fix tempControl all temp/rescale 100 3000.0 3000.0 0.02 1.0
run 10000

This time the result is:

This plot features a proper behaviour under the Tersoff potential. It also features another interesting trait - we see the atom oscillating around two different values, with a change of DC value between the 6000th step and the 7000th step. There is a simple explanation to this behaviour - the atom moved between two different sites on the lattice. Atoms 'travelling' through different sites in the lattice are a clear sign of a phase change of the material - from solid to liquid. This signifies a temperature high enough to break the SP3 bonds.

Next, let us try and run the simulation again, with a longer timestep and a shorter simulation:

timestep 0.0001
run 1000

The result this time:

Which is once again a proper behaviour. Notice the jittering at the last steps - this is the result of the minimization algorithm.

To conclude, let us try to increase the timestep even further:

timestep 0.01
run 100

And the result we get:

This is clearly not a proper behaviour for the Tersoff potential. We see that this timestep is too long to be of use. The sampling rate is too low to allow accurate measurements, leading to inaccurate computations of the system's dynamics.

This demonstration exhibits the importance of choosing a proper timestep, and the method used here can be of use when trying to determine a good timestep, as well as a good temperature for a phase change in a crystal.

To continue to useful scripts, press here.
To return backwards in the LAMMPS guide, press here.