We start by creating a new package using `cargo`

:

```
cargo new potential
cd potential
```

Open `Cargo.toml`

and add the lines

```
[dependencies]
lumol = {git = "https://github.com/lumol-org/lumol"}
```

to add the `lumol`

crate as a dependency to the package. To test if everything
works, run `cargo build`

and check if an error occurs.

For the first part of the tutorial, the complete code will be written into the
`lib.rs`

file.

The energy function of the Mie potential reads

\[u(x) = \varepsilon \frac{n}{n-m} \left(\frac{n}{m}\right)^{\frac{m}{n-m}} \left[ \left( \frac{\sigma}{x}\right)^n - \left( \frac{\sigma}{x}\right)^m \right]\]

where \(x\) denotes the distance between two interaction sites \(i, j\), with \(x = x_{ij} = | \mathbf{r}_j - \mathbf{r}_i |\). The parameters of the potential are

- \(n, m\) the repulsive and attractive exponents, respectively,
- \(\varepsilon\) the energetic paramater,
- \(\sigma\) the particle diameter or structural parameter.

We start by defining the `struct`

for our potential. Add the following lines
to `lib.rs`

:

```
use lumol::energy::{Potential, PairPotential};
#[derive(Clone, Copy)]
pub struct Mie {
/// Distance constant
sigma: f64,
/// Exponent of repulsive contribution
n: f64,
/// Exponent of attractive contribution
m: f64,
/// Energetic prefactor computed from the exponents and epsilon
prefactor: f64,
}
```

In the first two lines we define our imports from `Lumol`

, following with our
`Mie`

structure. Notice that we don’t store the `epsilon`

value, instead we
store an energetic prefactor that will make it easier to compute the potential.

\[\text{prefactor} = \varepsilon \frac{n}{n-m} \left(\frac{n}{m}\right)^{m/(n-m)}\]

Next, we implement a constructor function. That’s usefull in this case since we want to compute the prefactor of the potential once before we start our simulation.

In Rust we typically use `new`

for the constructors’ name.

```
impl Mie {
pub fn new(sigma: f64, epsilon: f64, n: f64, m: f64) -> Mie {
if m >= n {
panic!("The repulsive exponent n has to be larger than the attractive exponent m")
};
let prefactor = n / (n - m) * (n / m).powf(m / (n - m)) * epsilon;
return Mie {
sigma: sigma,
n: n,
m: m,
prefactor: prefactor,
}
}
}
```

Our function takes the parameter set as input, computes the prefactor and
returns a `Mie`

struct. Notice that it panics, for `n`

smaller than or equal
to `m`

. The next step is to implement the `Potential`

trait for `Mie`

.

`Potential`

¶Add the following lines below the structs implementation.

```
impl Potential for Mie {
fn energy(&self, r: f64) -> f64 {
let sigma_r = self.sigma / r;
let repulsive = f64::powf(sigma_r, self.n);
let attractive = f64::powf(sigma_r, self.m);
return self.prefactor * (repulsive - attractive);
}
fn force(&self, r: f64) -> f64 {
let sigma_r = self.sigma / r;
let repulsive = f64::powf(sigma_r, self.n);
let attractive = f64::powf(sigma_r, self.m);
return -self.prefactor * (self.n * repulsive - self.m * attractive) / r;
}
}
```

`energy`

is the implementation of the Mie potential equation shown above.
`force`

is the negative derivative of the energy with respect to the distance,
`r`

. To be more precise, the vectorial force can readily be computed by
multiplying the result of `force`

with the connection vector \(\vec{r}\).

The next step is to make our `Potential`

usable in Lumol’s algorithms to
compute non-bonded energies and forces. Therefore, we have to implement the
`PairPotential`

trait.

`PairPotential`

¶Let’s inspect the documentation for `PairPotential`

.

```
pub trait PairPotential: Potential + BoxClonePair {
fn tail_energy(&self, cutoff: f64) -> f64;
fn tail_virial(&self, cutoff: f64) -> f64;
fn virial(&self, r: &Vector3D) -> Matrix3 { ... }
}
```

First, we can see that `PairPotential`

enforces the implementation of
`Potential`

which is denoted by `pub trait PairPotential: Potential ...`

(we
ignore `BoxClonePair`

for now, as it is automatically implemented for us if we
implement `PairPotential`

manually). That makes sense from a didactive point
of view since we said that `PairPotential`

is a “specialization” of
`Potential`

and furthermore, we can make use of all functions that we had to
implement for `Potential`

.

There are three functions in the `PairPotential`

trait. The first two
functions start with `tail_`

. These are functions to compute long range or
tail corrections. Often, we introduce a cutoff distance into our potential
beyond which we set the energy to zero. Doing so we intoduce an error which we
can account for using a tail correction. We need two of these corrections, one
for the energy, `tail_energy`

, and one for the pressure (which uses
`tail_virial`

under the hood). For a beautiful derivation of tail corrections
for truncated potentials, see here.

The third function, `virial`

, already has its body implemented – we don’t
have to write an implementation for our potential.

We will omit the derivation of the formulae for tail corrections here but they are computed by solving these equations

\[\text{tail energy} = \int_{r_c}^{\infty} u(r) r^2 \mathrm{d}r\]

\[\text{tail virial} = \int_{r_c}^{\infty} \frac{\partial u(r)}{\partial r} r^3 \mathrm{d}r\]

The implementation looks like so

```
impl PairPotential for Mie {
fn tail_energy(&self, cutoff: f64) -> f64 {
if self.m < 3.0 {
return 0.0;
};
let sigma_rc = self.sigma / cutoff;
let n_3 = self.n - 3.0;
let m_3 = self.m - 3.0;
let repulsive = f64::powf(sigma_rc, n_3);
let attractive = f64::powf(sigma_rc, m_3);
return -self.prefactor * self.sigma.powi(3) * (repulsive / n_3 - attractive / m_3);
}
fn tail_virial(&self, cutoff: f64) -> f64 {
if self.m < 3.0 {
return 0.0;
};
let sigma_rc = self.sigma / cutoff;
let n_3 = self.n - 3.0;
let m_3 = self.m - 3.0;
let repulsive = f64::powf(sigma_rc, n_3);
let attractive = f64::powf(sigma_rc, m_3);
return -self.prefactor * self.sigma.powi(3) * (repulsive * self.n / n_3 - attractive * self.m / m_3);
}
}
```

Note that we cannot correct every kind of energy function. In fact, the
potential has to be a *short ranged* potential. For our Mie potential, both the
exponents have to be larger than 3.0 else our potential will be *long ranged*
and the integral that has to be solved to compute the tail corrections diverges.
We return zero in that case.

That concludes the first part. To test your new and shiny potential, you can run
a small simulation. You’ll find a minimal Monte Carlo simulation example in the
`tutorials/potential`

directory of the main lumol repository
where you will also find the `src/lib.rs`

file we created in this tutorial.
You can then run the simulation via

```
cargo run --release
```

Fantastic! You implemented a new potential and ran a simulation with it!

If you want to share your implementation with other Lumol users only some small additional steps are neccessary. We will talk about them in the next part of this tutorial (which is not yet written).