Hello Sodium Chloride

In this example, we will simulate a Sodium Chloride crystal using molecular dynamics to propagate the position throughout time. Sodium Chloride adds a challenge in simulations because each atom carries a charge. These charges interact with a Coulomb potential which goes to zero as \(1 / r\). The problem is that the cutoff scheme used for pair potentials in most molecular simulations can not be used for anything that goes to zero slower than \(1 / r^3\). So we need to use alternative methods to compute the potential for the interactions between pairs of charges.

For this simulation, you will need the following files:

  • the initial configuration nacl.xyz`
  • the input file nacl.toml

You can download both files here.

Again, you can run the simulation which should complete in a minute with lumol nacl.toml. This will perform a molecular dynamics simulation of a NaCl crystal using electrostatic interactions between atomic charges.

The input file commented

We start with the input version again:

[input]
version = 1

Then we load the system from the nacl.xyz file and define the unit cell.

[[systems]]
file = "nacl.xyz"
cell = 22.5608

Next, we define our potential. This section is way bigger than the one for our previous Lennard-Jones example:

[systems.potentials.global]
cutoff = "8 A"

[systems.potentials.charges]
Na = 1.0
Cl = -1.0

[[systems.potentials.pairs]]
atoms = ["Na", "Cl"]
lj = {sigma = "3.5545 A", epsilon = "0.04425 kcal/mol"}

[[systems.potentials.pairs]]
atoms = ["Na", "Na"]
lj = {sigma = "2.497 A", epsilon = "0.07826 kcal/mol"}

[[systems.potentials.pairs]]
atoms = ["Cl", "Cl"]
lj = {sigma = "4.612 A", epsilon = "0.02502 kcal/mol"}

[systems.potentials.coulomb]
wolf = {cutoff = "8 A"}

Let’s break it down. First, we define some global values for the interactions: setting systems.potentials.global.cutoff will use the given cutoff for all pair interactions in the system. We set the charges of atoms in the systems.potentials.charges section.

[systems.potentials.global]
cutoff = "8 A"

[systems.potentials.charges]
Na = 1.0
Cl = -1.0

Then, we need to define the pair interactions for all possible pair combinations in the system, i.e. (Na, Na), (Cl, Cl), and (Na, Cl).

[[systems.potentials.pairs]]
atoms = ["Na", "Cl"]
lj = {sigma = "3.5545 A", epsilon = "0.04425 kcal/mol"}

[[systems.potentials.pairs]]
atoms = ["Na", "Na"]
lj = {sigma = "2.497 A", epsilon = "0.07826 kcal/mol"}

[[systems.potentials.pairs]]
atoms = ["Cl", "Cl"]
lj = {sigma = "4.612 A", epsilon = "0.02502 kcal/mol"}

Because our system contains charges, we need to use an electrostatic potential solver. Here we are going for the Wolf solver, with a cutoff of 8 A. Note that if we’d chose a cutoff different from the global one defined above, we would overwrite the global one for the coulombic interactions.

[systems.potentials.coulomb]
wolf = {cutoff = "8 A"}

We can now define the simulation and the outputs for this simulation.

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

We are using a molecular dynamics simulation of the NaCl crystal with a timestep of 1 fs for integration (this will produce a NVE ensemble).

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

As before, run the simulation via

lumol nacl.toml