Molecular dynamics simulation of water

In this example, we will simulate a molecular liquid of high scientific interest: water. At the same time we will learn how we can reuse a potential definition between multiple simulations.

You will need three input files for this simulation:

• the initial configuration water.xyz;
• the simulation input water.toml;
• the force field (potential definitions) water-fSCP.toml.

You can download them here.

As usual, you can run the simulation with lumol water.toml. The simulation should finish in a few minutes.

The input files commented

water.toml file

The main input file is pretty similar to the previous examples, with two novelties:

• The guess_bonds = true entry tells lumol to try to guess the bonds from the distances in the XYZ file. This is needed because we want to simulate a molecule but there is no bonding information inside the XYZ format. If we were to use a PDB file with connectivity instead, this would not be needed;
• The potentials = "water-fSCP.toml" tells lumol to read the potentials from the water-fSCP.toml file. Using this type of input for the potentials we can reuse the same potential for multiple simulations and also easily change the potential used in the simulation.

[input]
version = 1

[[systems]]
file = "water.xyz"
guess_bonds = true
cell = 28.0
potentials = "water-fSCP.toml"

[[simulations]]
nsteps = 5000
outputs = [
    {type = "Trajectory", file = "trajectory.xyz", frequency = 10}
]

[simulations.propagator]
type = "MolecularDynamics"
timestep = "1 fs"

water-fSCP.toml file

water-fSCP.toml is a standalone potential input file. It contains the same data as a potential definition inside a [[systems]] section, but without the system.potential prefix on all section names.

In this input, we start by defining which version of the input format we are using.

[input]
version = 1

Then, we can define some global input data: the pair potential cutoff and the atomic charges.

[global]
cutoff = "14 A"

[charges]
O = -0.82
H = 0.41

The pair potential section contains the usual declarations for pairs, with a few additional options.

[pairs]
O-O = {type = "lj", sigma = "3.16 A", epsilon = "0.155 kcal/mol"}

We can add a restriction to restrict a specific pair interaction to some kind of pairs. Here, we will only account for H-H pairs inside the same molecule ("intra-molecular" interactions).

H-H = {type = "harmonic", k = "79.8 kcal/mol/A^2", x0 = "1.633 A", restriction = "intra-molecular"}

We can also define a non-interacting pair interaction! This is useful when some atoms do not interact in the model we use.

H-O = {type = "null"}

Next, the bonds and angles are defined. These are interactions used between bonded particles (or angles formed by two bonds and dihedral angles formed by three bonds). This section follows the same pattern as the [pairs] section.

[bonds]
O-H = {type = "harmonic", k = "1054.2 kcal/mol/A^2", x0 = "1.0 A"}

[angles]
H-O-H = {type = "harmonic", k = "75.9 kcal/mol/rad^2", x0 = "109.5 deg"}

Finally we specify how to compute the electrostatic interaction, this time using the Ewald summation method. We can restrict the coulombic interactions to only apply between particles not in the same molecule by using restriction = "inter-molecular".

[coulomb]
ewald = {cutoff = "8.5 A", kmax = 3}
restriction = "inter-molecular"

The simulation is run via

lumol water.toml