LAMMPS from scratch · 20 min

read_data files

Real systems are rarely built with lattice commands: they are read from a data file that lists every atom, bond and angle. You write one by hand for a single water molecule, let LAMMPS read it, and check the bond and angle energies, the excluded Coulomb energy and the dipole moment by hand, then edit the file and check again.

The layout of a data file: header, Masses, Atoms, Bonds, Angles, and the columns of atom_style full
Predict bond, angle and Coulomb energies and the dipole of a water molecule by hand
Why the Coulomb energy inside a molecule is zero by default (special_bonds)

1 · Understand: a data file is the system, written out

Instead of lattice, region, create_box and create_atoms, block 2 of the script (page A.3) can be one line: read_data water.data. The file has a header with counts and the box, then sections:

Part	Contains	Water
header	numbers of atoms, bonds, angles, and of their types; the box	3 atoms, 2 bonds, 1 angle; 2 atom types
`Masses`	one mass per atom type	O 15.9994, H 1.008 g/mol
`Atoms # full`	one line per atom; the columns depend on atom_style	id, molecule, type, charge, x, y, z
`Bonds`	bond id, bond type, the two atom ids	O–H1, O–H2
`Angles`	angle id, angle type, three atom ids (the middle one is the vertex)	H1–O–H2

atom_style	columns of an Atoms line
atomic	id type x y z
charge	id type q x y z
full	id molecule type q x y z

Three optional integers at the end of an Atoms line are image flags (page D.2).

2 · Modify: the data file and the script

Edit any number in the data file (for example a coordinate). The hand calculation below updates at once; then run LAMMPS.

water.data

the input script

count the H–H Coulomb energy (special_bonds lj/coul 0.0 1.0 0.0)

3 · Predict by hand

LAMMPS's harmonic bond and angle are E = K(r − r₀)² and E = K(θ − θ₀)²: there is no ½ in front, unlike the textbook ½k x². Force fields that publish k must be divided by 2 before use.

What the molecule shows

Every energy is a few lines of arithmetic once the coordinates are known; the file is the whole system.
The Coulomb energy is zero by default even though the atoms carry charges: atoms joined by a bond (1–2) or by an angle (1–3) are excluded (special_bonds 0 0 0). The two H atoms, 1.61 Å apart, would otherwise add about +37 kcal/mol, far more than the bond and angle energies. Each force field states its own exclusion rule; the data file and the script must follow it.
The dipole μ = Σ q r comes out as 0.529 e·Å = 2.54 D for this geometry (1 e·Å = 4.8032 D).

Your turn

1. You change the header to 4 atoms but leave three Atoms lines. What happens?

LAMMPS reads the counts first and reads exactly that many Atoms lines, so the count and the file disagree. This version stops with ERROR: Did not assign all atoms correctly. Counts in the header must match the sections.

2. A paper gives the O–H bond as k = 1059.2 kcal/mol/Å² in E = ½k(r − r₀)². What do you write in bond_coeff?

3. Move H2 so that the angle is exactly 109.47°. What is E_angle, and does the dipole get larger or smaller?