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.

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
```