Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP]Tutorial for pair potential implementation #196

Merged
merged 6 commits into from
Oct 19, 2017
Merged
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,12 @@ script:
# Run unit tests and doc tests in debug mode
- cargo test -p lumol-core
- cargo test -p lumol-input
# Check if tutorials compile
- cd
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like a previous attempt leftover, but this might break the build, as the code is not in $HOME.

# Check that benchmarks still compile
- cargo bench --no-run
# Check if tutorials compile
- ./scripts/ci/compile-tutorials.py
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would rather explicitly do cd tutorial/lumol_potential_tutorial && cargo build here. We can craft a script for this later!

You might need to do a cd $TRAVIS_BUILD_DIR after that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure, that was my first attempt :)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is another way to do this, in one line: cargo build --manifest-path=tutorial/lumol_potential_tutorial/Cargo.toml.

But if you prefer the python script, you'll need to chmod +x it. There is a permission denied error on Travis.

# Misc style checking
- ./scripts/ci/check-whitespaces.py
# Generate and check doc
Expand Down
1 change: 1 addition & 0 deletions Cargo.toml
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ members = [
"src/core",
"src/input"
]
exclude = ["tutorials"]

[[bin]]
name = "lumol"
Expand Down
4 changes: 3 additions & 1 deletion doc/src/SUMMARY.md
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,9 @@
- [Logging configuration](input/log.md)

- [Advanced tutorials]()
- [Adding potentials]()
- [Adding potentials](advanced/custom_potential/intro.md)
- [Part I: Separate package](advanced/custom_potential/external.md)
- [Part II: Lumol core](advanced/custom_potential/core.md)
- [Adding outputs]()
- [Extending Molecular dynamics]()
- [Extending Monte Carlo]()
Expand Down
1 change: 1 addition & 0 deletions doc/src/advanced/custom_potential/core.md
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
# Part II: Adding the potential to lumol core
204 changes: 204 additions & 0 deletions doc/src/advanced/custom_potential/external.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,204 @@
# Part I: Implementation of the logic in a separate package

We start by creating a new package using `cargo`:

```
cargo new lumol_tutorial_potential
cd lumol_tutorial_potential
```

Open `Cargo.toml` and add the lines
```toml
[dependencies]
lumol = {git = "https://github.com/lumol-org/lumol"}
```

to add `Lumol` as a dependency to the package.
To test if everything works, open the `lib.rs` file in `src` and add
```rust
extern crate lumol;
```
in the first line.

Run `cargo build` or `cargo test` to check if errors occur.

## Defining the struct

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) = \frac{n}{n-m} \left(\frac{n}{m}\right)^{m/(n-m)}\epsilon \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,
- $\epsilon$ 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`:

```rust
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you think it would be possible to directly include the relevant part of the file in /tutorials ? This would make it easier to prevent the code from bit rotting.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm note sure what you mean here. How would you directly include that part of the file?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I checked and this is not possible with markdown. So let's keep it this way for now.

I do this with sphinx: you have example files somewhere, and then you include them in your documentation. (see here for example). This allow to ensure that the documentation and the actual files stays in sync.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah I see, yes i do that in sphinx too.

extern crate lumol;
use lumol::energy::{Potential, PairPotential};

#[derive(Clone, Copy)]
pub struct Mie {
/// Distance constant
sigma: f64,
/// Energy constant
epsilon: 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 structure.

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.

$\text{prefactor} = \frac{n}{n-m} \left(\frac{n}{m}\right)^{m/(n-m)}\epsilon$

In Rust we typically use `new` for the constructors' name.

```rust
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;
Mie {
sigma: sigma,
epsilon: epsilon,
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`.

## Implementing `Potential`

Add the following lines below the structs implementation.

```rust
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);
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);
-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 algortihms to compute
non-bonded energies and forces.
Therefore, we have to implement the `PairPotential` trait.

## Implementing `PairPotential`

Let's inspect the [documentation of `PairPotential`](http://lumol.org/lumol/latest/lumol/energy/trait.PairPotential.html).

```rust
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](https://engineering.ucsb.edu/~shell/che210d/Simulations_of_bulk_phases.pdf).

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

energy: $\int_{r_c}^{\infty} u(r) r^2 \mathrm{d}r$

virial: $\int_{r_c}^{\infty} \frac{\partial u(r)}{\partial r} r^3 \mathrm{d}r$

The implementation looks like so

```rust
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);
-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);
-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.

## Running a small simulation using the new potential

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 `lumol/tutorials/lumol_tutorial_potential` directory 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 part II.
64 changes: 64 additions & 0 deletions doc/src/advanced/custom_potential/intro.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
# Creating your own potential function

In this tutorial we will teach you how to use your own potential function with Lumol.
This tutorial consists of two parts.
In the first part, we will implement the logic in a separate project using Lumol as an external library.
The second part describes the neccessary steps to implement the logic into the core package of Lumol called `lumol-core`.

A potential is a function that describes the energy and force between interaction sites.
In Lumol we differentiate between two types of potentials: First, there are potential functions that need the global state of the system, *i.e.* all positions, as input. The Ewald summation which is used to compute electrostatic interactions is an example for this kind of potential.
We call them `GlobalPotential`s.
And second, we have potentials that take a single geometric parameter as input.
This geometric parameter can be for example a distance or an angle.
Typical examples are Van-der-Waals potentials (*e.g.* Lennard-Jones) and potential functions describing covalent bonds (*e.g.* harmonic potential, cosine potential, torsion, *etc.*).
There are plenty of potentials falling into this category, hence in Lumol we simply call them `Potential`s.

> There are two kind of potentials. A `GlobalPotential` takes the global systems' state as input
and a `Potential` takes a single scalar value (distance, angle, ...) as input.

In this tutorial we will focus on the implementation of the latter.
More specific, we will implement a potential to compute van-der-Waals interactions between pairs of particles.
We will make liberal use of the API documentation of both the Rust standard library as well as the Lumol API.
Please note that we will point to other references, such as the Rust book, concerning general Rust concepts to keep this tutorial brief.
If you have questions concerning Rust or Lumol, please don't hesitate to file an issue on [github](https://github.com/lumol-org/lumol/issues) or join the discussion on our [gitter](https://gitter.im/lumol-org/lumol).

The tutorial is structured as follows: First we will have a look at how `Potentials` are represented (what data structures are used) and what functionalities we have to implement. Then we describe how we add those functionalities. We will write a small simulation program that makes use of our newly created potential.
In the second part we talk about implementing the potential into Lumol's core.
Rust offers beautiful utilities to add documentation inside the code. We will have a look at the documentation of currently implemented potentials to guide us. We will then add some tests to make sure that our implementation is correct. This concludes the bulk the work, but to make our new potential function accessible to all Lumol users we will also add a parsing function. Doing so, our potential can be conveniently specified from an input file. We conclude the tutorial by adding a short documentation to the user manual.

We have a lot to do. Ready? Let's go.

## Rust prerequisites

We recommend reading these chapters in the rust book.

- [structs](https://doc.rust-lang.org/book/second-edition/ch05-00-structs.html)
- [traits](https://doc.rust-lang.org/book/second-edition/ch10-02-traits.html)

## Representation

As mentioned above, potential functions can be used to model all kinds of interactions between particles such as bonded and non-bonded interactions. In Lumol, `Potential` is a [trait](https://doc.rust-lang.org/book/second-edition/ch10-02-traits.html). To further distinguish between bonded interactions (bond lengths, angles and dihedrals) and non-bonded interactions, we use another trait (often called marker traits). Your possible options to further specialise a `Potential` are

- [`PairPotential`](http://lumol.org/lumol/latest/lumol/energy/trait.PairPotential.html) for non-bonded two body interactions;
- [`BondPotential`](http://lumol.org/lumol/latest/lumol/energy/trait.BondPotential.html) for covalent bonds interactions;
- [`AnglePotential`](http://lumol.org/lumol/latest/lumol/energy/trait.AnglePotential.html) for covalent angles interactions;
- [`DihedralPotential`](http://lumol.org/lumol/latest/lumol/energy/trait.DihedralPotential.html) for covalent dihedral angles interactions.

In this tutorial, we will implement a potential to describe non-bonded pair interactions, namely the [**Mie potential**](http://www.sklogwiki.org/SklogWiki/index.php/Mie_potential). We will have to implement both the `Potential` as well as the `PairPotential` traits. If we wanted to implement a function that can be used as non-bonded as well as bond-length potential, we'd have to implement `Potential`, `PairPotential` as well as `BondPotential`. "Implementing a trait" means that we will define a [structure](https://doc.rust-lang.org/book/second-edition/ch05-00-structs.html) (a `struct`) for which we will add functions to satisfy the traits' requirements.

Let's start by having a look at the API documentation for `Potential`: open the [API documentation](http://lumol.org/lumol/latest/lumol/energy/trait.Potential.html).

As you can see from the trait, a `Potential` defines the interface for two functions, `energy` and `force` (we ignore the `Sync + Send` statement for now):

```rust
pub trait Potential: Sync + Send {
fn energy(&self, x: f64) -> f64;
fn force(&self, x: f64) -> f64;
}
```

`energy` will compute the interaction energy between two sites as a function of `x`. The force is defined as the negative derivative of the energy function with respect to `x`.

Both functions take a single, scalar argument and return a single scalar value. In our case `x` stands for the distance between two interaction sites. Note that only the function definitions -- without a function body -- are specified. We will have to implement these functions for our potential.

Let's start the implementation.
12 changes: 12 additions & 0 deletions scripts/ci/compile-tutorials.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
#!/usr/bin/env python
# coding=utf-8
import os
import subprocess
import glob

TUTORIALS = glob.glob('{}/tutorials/*'.format(os.curdir))

if __name__ == '__main__':
for tutorial in TUTORIALS:
subprocess.call('Cargo build', shell=True, cwd=tutorial)

9 changes: 9 additions & 0 deletions tutorials/lumol_tutorial_potential/Cargo.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
[package]
name = "lumol_tutorial_potential"
version = "0.0.0"
publish = false
authors = ["Gernot Bauer <[email protected]>"]

[dependencies]
lumol = {git = "https://github.com/lumol-org/lumol"}
lumol-input = {git = "https://github.com/lumol-org/lumol"}
26 changes: 26 additions & 0 deletions tutorials/lumol_tutorial_potential/data/argon.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
[input]
version = 1

[[systems]]
file = "argon.xyz"
cell = 31.0

[systems.potentials.global]
cutoff = "14.0 A"
tail_correction = true

[[simulations]]
nsteps = 1_00000
outputs = [
{type = "Energy", file = "mie_energy.dat", frequency = 500},
{type = "Properties", file = "mie_prp.dat", frequency = 500}
]

[simulations.propagator]
type = "MonteCarlo"
temperature = "217.0 K"
update_frequency = 500

moves = [
{type = "Translate", delta = "10 A", frequency = 500, target_acceptance = 0.5},
]
Loading