# Molecular dynamics in PyRETIS with C or FORTRAN¶

In this example, we will implement a Lennard-Jones potential in C or FORTRAN and use this to run a MD simulation with PyRETIS for a simple system. The system we consider is a Lennard-Jones mixture starting for a specific configuration as shown below. After running for some steps (which melt the initial configuration), we reverse the velocities and continue the simulation.

Fig. 35 Snapshots from the MD simulation. To the left, the initial (and final) configuration and to the right, the configuration after running 2000 steps, just before reversing the velocities. The system consist of two types of particles, labeled Ar (colored silver) and Kr (colored purple), which are initially arranged in a particular configuration.

## Creating a new potential in C or FORTRAN¶

In this example, we will create a new potential function to use with PyRETIS in C or FORTRAN.

As discussed in the introduction to the library, we need to create a new potential class which PyRETIS will make use of. This new potential class will import routines from a C or FORTRAN library and use these routines to do the actual computation of the forces.

We will in the following show how this can be done, both with FORTRAN and with C. Essentially, we are going to complete the following steps:

1. We write the code responsible for the evaluation of the force in an external C or FORTRAN library.
2. We compile the external code. Here we do one of the following:
• For C, we create a setup.py which is used to compile.
• For FORTRAN, make use of a Makefile and the program f2py.
3. We write a Python module containing a new PotentialFunction sub-class representing the new potential. This class imports the library which we created in steps 1 and 2.

### Writing a new potential function with FORTRAN¶

We will now create a potential function writing the external library in FORTRAN. This is done in the following steps:

#### Step 1. Creating the FORTRAN code¶

The FORTRAN code for evaluating the potential and forces is relative simply and we assume that we are handed positions, velocities, and forces as double precision arrays and that we can directly make use of them.

! Copyright (c) 2021, PyRETIS Development Team.
module ljfortran

implicit none

private

public :: potential, force, potential_and_force
public :: pbc_dist

contains

function pbc_dist(v, box, ibox) result(dist)
! Function to apply periodic boundaries to a given displacement.
!
! Parameters
! ----------
! v : the displacement to pbc-wrap
! box : the box lengths
! ibox : the inverse box lengths, ibox = 1.0 / box
!
! Returns
! -------
! dist : the pbc-wrapped displacement
implicit none
double precision :: dist
double precision, intent(in) :: v, box, ibox
dist = v - nint(v * ibox) * box
end function pbc_dist

subroutine potential(pos, box, ibox, lj3, lj4, offset, rcut2, ptype, n, d, p, vpot)
! Function to evaluate the Lennard-Jones potential
!
! Parameters
! ----------
! pos : positions of the particles
! box : the box lengths
! ibox : the inverse box lengths, ibox = 1.0 / box
! lj3 : Lennard-Jones parameters, 4.0 * epsilonij * sigmaij**12
! lj4 : Lennard-Jones parameters, 4.0 * epsilonij * sigmaij**6
! rcut2 : Squared cut-offs
! offset : Potential energy shift
! ptype : identifier for the particle types
! n : number of particles
! d : dimensionality of the vectors
! p : number of particle types
!
! Returns
! -------
! vpot : the potential energy
implicit none
integer, intent(in) :: n, d, p
double precision, dimension(n, d), target, intent(in) :: pos
double precision, dimension(d), intent(in) :: box, ibox
double precision, dimension(p, p), intent(in) :: lj3, lj4, offset, rcut2
integer, dimension(n), intent(in) :: ptype
double precision, intent(out) :: vpot
double precision :: rsq, r2inv, r6inv, dx, dy, dz
double precision, dimension(:), pointer :: x, y, z
integer :: i, j, itype, jtype
vpot = 0.0D0
x => pos(1:n,1)
y => pos(1:n,2)
z => pos(1:n,3)
do i=1,n-1
itype = ptype(i) + 1
do j=i+1,n
jtype = ptype(j) + 1
dx = x(i) - x(j)
dy = y(i) - y(j)
dz = z(i) - z(j)
dx = pbc_dist(dx, box(1), ibox(1))
dy = pbc_dist(dy, box(2), ibox(2))
dz = pbc_dist(dz, box(3), ibox(3))
rsq = dx*dx + dy*dy + dz*dz
if (rsq < rcut2(jtype, itype)) then
r2inv = 1.0D0 / rsq
r6inv = r2inv**3
vpot = vpot + r6inv * (lj3(jtype, itype) * r6inv - lj4(jtype, itype))
vpot = vpot - offset(jtype, itype)
end if
end do
end do
end subroutine potential

subroutine force(pos, box, ibox, lj1, lj2, rcut2, ptype, n, d, p, forces, virial)
! Function to evaluate the Lennard-Jones force
!
! Parameters
! ----------
! pos : positions of the particles
! box : the box lengths
! ibox : the inverse box lengths, ibox = 1.0 / box
! lj1 : Lennard-Jones parameters, 48.0 * epsilonij * sigmaij**12
! lj2 : Lennard-Jones parameters, 24.0 * epsilonij * sigmaij**6
! rcut2 : Squared cut-offs
! ptype : identifier for the particle types
! n : number of particles
! d : dimensionality of the vectors
! p : number of particle types
!
! Returns
! -------
! forces : the Lennard-Jones forces on the particles
! virial : the virial matrix
implicit none
integer, intent(in) :: n, d, p
double precision, dimension(n, d), target, intent(out) :: forces
double precision, dimension(d, d), intent(out) :: virial
double precision, dimension(n, d), target, intent(in) :: pos
double precision, dimension(d), intent(in) :: box, ibox
double precision, dimension(p, p), intent(in) :: lj1, lj2, rcut2
integer, dimension(n), intent(in) :: ptype
double precision, dimension(:), pointer :: x, y, z, fx, fy, fz
double precision :: rsq, r2inv, r6inv, forcelj, dx, dy, dz
double precision :: forceij_x, forceij_y, forceij_z
integer :: i, j, itype, jtype
x => pos(1:n, 1)
y => pos(1:n,2)
z => pos(1:n,3)
fx => forces(1:n,1)
fy => forces(1:n,2)
fz => forces(1:n,3)
forces = 0.0D0
virial = 0.0D0
do i=1,n-1
itype = ptype(i) + 1
do j=i+1,n
jtype = ptype(j) + 1
dx = x(i) - x(j)
dy = y(i) - y(j)
dz = z(i) - z(j)
dx = pbc_dist(dx, box(1), ibox(1))
dy = pbc_dist(dy, box(2), ibox(2))
dz = pbc_dist(dz, box(3), ibox(3))
rsq = dx*dx + dy*dy + dz*dz
if (rsq < rcut2(jtype, itype)) then
r2inv = 1.0D0 / rsq
r6inv = r2inv*r2inv*r2inv
forcelj = r2inv * r6inv * (lj1(jtype, itype) * r6inv - lj2(jtype, itype))
forceij_x = forcelj * dx
forceij_y = forcelj * dy
forceij_z = forcelj * dz
fx(i) = fx(i) + forceij_x
fy(i) = fy(i) + forceij_y
fz(i) = fz(i) + forceij_z
fx(j) = fx(j) - forceij_x
fy(j) = fy(j) - forceij_y
fz(j) = fz(j) - forceij_z
! accumulate for the virial:
virial(1,1) = virial(1,1) + forceij_x * dx
virial(1,2) = virial(1,2) + forceij_x * dy
virial(1,3) = virial(1,3) + forceij_x * dz
virial(2,1) = virial(2,1) + forceij_y * dx
virial(2,2) = virial(2,2) + forceij_y * dy
virial(2,3) = virial(2,3) + forceij_y * dz
virial(3,1) = virial(3,1) + forceij_z * dx
virial(3,2) = virial(3,2) + forceij_z * dy
virial(3,3) = virial(3,3) + forceij_z * dz
end if
end do
end do
end subroutine force

subroutine potential_and_force(pos, box, ibox, lj1, lj2, lj3, lj4, offset, rcut2, ptype, n, d, p, forces, virial, vpot)
! Function to evaluate the Lennard-Jones force and potential
!
! Parameters
! ----------
! pos : positions of the particles
! box : the box lengths
! ibox : the inverse box lengths, ibox = 1.0 / box
! lj1 : Lennard-Jones parameters, 48.0 * epsilonij * sigmaij**12
! lj2 : Lennard-Jones parameters, 24.0 * epsilonij * sigmaij**6
! lj3 : Lennard-Jones parameters, 4.0 * epsilonij * sigmaij**12
! lj4 : Lennard-Jones parameters, 4.0 * epsilonij * sigmaij**6
! rcut2 : Squared cut-offs
! offset : Potential energy shift
! ptype : identifier for the particle types
! n : number of particles
! d : dimensionality of the vectors
! p : number of particle types
!
! Returns
! -------
! vpot : the potential energy
! forces : the Lennard-Jones forces on the particles
! virial : the virial matrix
implicit none
integer, intent(in) :: n, d, p
double precision, dimension(n, d), target, intent(out) :: forces
double precision, dimension(d, d), intent(out) :: virial
double precision, intent(out) :: vpot
double precision, dimension(n, d), target, intent(in) :: pos
double precision, dimension(d), intent(in) :: ibox, box
double precision, dimension(p, p), intent(in) :: lj1, lj2, lj3, lj4, offset, rcut2
double precision, dimension(:), pointer :: x, y, z, fx, fy, fz
integer, dimension(n), intent(in) :: ptype
double precision :: rsq, r2inv, r6inv, forcelj, dx, dy, dz
double precision :: forceij_x, forceij_y, forceij_z
integer :: i, j, itype, jtype
forces = 0.0D0
virial = 0.0D0
vpot = 0.0D0
x => pos(1:n,1)
y => pos(1:n,2)
z => pos(1:n,3)
fx => forces(1:n,1)
fy => forces(1:n,2)
fz => forces(1:n,3)
do i=1,n-1
itype = ptype(i) + 1
do j=i+1,n
jtype = ptype(j) + 1
dx = x(i) - x(j)
dy = y(i) - y(j)
dz = z(i) - z(j)
dx = pbc_dist(dx, box(1), ibox(1))
dy = pbc_dist(dy, box(2), ibox(2))
dz = pbc_dist(dz, box(3), ibox(3))
rsq = dx*dx + dy*dy + dz*dz
if (rsq < rcut2(jtype, itype)) then
r2inv = 1.0D0 / rsq
r6inv = r2inv*r2inv*r2inv
forcelj = r2inv * r6inv * (lj1(jtype, itype) * r6inv - lj2(jtype, itype))
forceij_x = forcelj * dx
forceij_y = forcelj * dy
forceij_z = forcelj * dz
fx(i) = fx(i) + forceij_x
fy(i) = fy(i) + forceij_y
fz(i) = fz(i) + forceij_z
fx(j) = fx(j) - forceij_x
fy(j) = fy(j) - forceij_y
fz(j) = fz(j) - forceij_z
vpot = vpot + r6inv * (lj3(jtype, itype) * r6inv - lj4(jtype, itype))
vpot = vpot - offset(jtype, itype)
! accumulate for the virial:
virial(1,1) = virial(1,1) + forceij_x * dx
virial(1,2) = virial(1,2) + forceij_x * dy
virial(1,3) = virial(1,3) + forceij_x * dz
virial(2,1) = virial(2,1) + forceij_y * dx
virial(2,2) = virial(2,2) + forceij_y * dy
virial(2,3) = virial(2,3) + forceij_y * dz
virial(3,1) = virial(3,1) + forceij_z * dx
virial(3,2) = virial(3,2) + forceij_z * dy
virial(3,3) = virial(3,3) + forceij_z * dz
end if
end do
end do
end subroutine potential_and_force

end module ljfortran


#### Step 2. Creating the Makefile and compiling¶

In order to compile the FORTRAN code created in the previous step, we make use of f2py.

Note that in some systems you might actually need to use f2py3 (or, for a specific version: f2py3.7 etc.) if you have several versions of Python installed on your system. This is to ensure that the FORTRAN code is compiled with a version that matches the version of Python you are using. You will then have to modify the F2PY = f2py setting in the Makefile.

F2PY = f2py
F2PY_FLAGS = --fcompiler=gfortran

MODULE=ljfortran

all: ${MODULE}.so${MODULE}.so: ${MODULE}.f90${F2PY} ${F2PY_FLAGS} -c$< -m ${MODULE} clean: find -name \*.so -delete find -name \*.pyc -delete find -name __pycache__ -delete veryclean:${MAKE} clean
find -name \*energy.dat* -delete
find -name \*thermo.txt* -delete
find -name \*traj.gro* -delete


Note that the FORTRAN compiler is specified using --fcompiler=gfortran and other choices can be seen by running:

f2py -c --help-fcompiler


Further, the Makefile assumes that the FORTRAN module is named ljfortran.f90.

#### Step 3. Creating a new Python class for the potential function¶

For the Python class representing the potential function, we import the module we just compiled and make use of the methods defined in that module.

# -*- coding: utf-8 -*-
# Copyright (c) 2021, PyRETIS Development Team.
"""Example of using a Lennard-Jones potential implemented in FORTRAN."""
import sys
import os
import logging
import numpy as np
from pyretis.forcefield.potentials import PairLennardJonesCut
from pyretis.forcefield.potentials.pairpotentials import (
generate_pair_interactions
)
# Just to handle imports of the library:
sys.path.insert(0, os.path.abspath(os.path.dirname(__file__)))
logger = logging.getLogger(__name__)  # pylint: disable=invalid-name
try:
from ljfortran import ljfortran
except ImportError:
MSG = ('Could not import external FORTRAN library.'
logger.critical(MSG)
raise ImportError(MSG)

__all__ = ['PairLennardJonesCutF']

class PairLennardJonesCutF(PairLennardJonesCut):
r"""class PairLennardJonesCutF(PairLennardJonesCut).

This class implements as simple Lennard-Jones 6-12 potential which
employs a simple cut-off and can be shifted. The potential energy
(:math:V_\text{pot}) is defined in the usual way for an
interacting pair of particles a distance :math:r apart,

.. math::

V_\text{pot} = 4 \varepsilon \left( x^{12} - x^{6} \right),

where :math:x = \sigma/r and :math:\varepsilon
and :math:\sigma are the potential parameters. The parameters are
stored as attributes of the potential and we store one set for each
kind of pair interaction. Parameters can be generated with a
specific mixing rule by the force field.

Attributes
----------
params : dict
The parameters for the potential. This dict is assumed to
contain parameters for pairs, i.e. for interactions.
_lj1 : numpy.array
Lennard-Jones parameters used for calculation of the force.
Keys are the pairs (particle types) that may interact.
Calculated as: 48.0 * epsilon * sigma**12
_lj2 : numpy.array
Lennard-Jones parameters used for calculation of the force.
Keys are the pairs (particle types) that may interact.
Calculated as: 24.0 * epsilon * sigma**6
_lj3 : numpy.array
Lennard-Jones parameters used for calculation of the potential.
Keys are the pairs (particle types) that may interact.
Calculated as: 4.0 * epsilon * sigma**12
_lj4 : numpy.array
Lennard-Jones parameters used for calculation of the potential.
Keys are the pairs (particle types) that may interact.
Calculated as: 4.0 * epsilon * sigma**6
_offset : numpy.array
Potential values for shifting the potential if requested.
This is the potential evaluated at the cutoff.
_rcut2 : numpy.array
The squared cut-off for each interaction type.
Keys are the pairs (particle types) that may interact.

"""

def __init__(self, dim=3, shift=True, mixing='geometric',
desc='Lennard-Jones pair potential (FORTRAN)'):
"""Set up the Lennard-Jones potential.

Parameters
----------
dim : int
The dimensionality to use.
shift : boolean
Determines if the potential should be shifted or not.
mixing : string
Determines how we should mix potential parameters.

"""
super().__init__(dim=dim, desc=desc, mixing=mixing)
self.ntype = 0

def set_parameters(self, parameters):
"""Update all parameters.

Here, we generate pair interactions, since that is what this
potential actually is using.

Parameters
----------
parameters : dict
The input base parameters.

"""
self.params = {}
pair_param = generate_pair_interactions(parameters, self.mixing)
self.ntype = max(int(np.sqrt(len(pair_param))), 2)
self._lj1 = np.zeros((self.ntype, self.ntype))
self._lj2 = np.zeros_like(self._lj1)
self._lj3 = np.zeros_like(self._lj1)
self._lj4 = np.zeros_like(self._lj1)
self._rcut2 = np.zeros_like(self._lj1)
self._offset = np.zeros_like(self._lj1)
for pair in pair_param:
eps_ij = pair_param[pair]['epsilon']
sig_ij = pair_param[pair]['sigma']
rcut = pair_param[pair]['rcut']
self._lj1[pair] = 48.0 * eps_ij * sig_ij**12
self._lj2[pair] = 24.0 * eps_ij * sig_ij**6
self._lj3[pair] = 4.0 * eps_ij * sig_ij**12
self._lj4[pair] = 4.0 * eps_ij * sig_ij**6
self._rcut2[pair] = rcut**2
vcut = 0.0
if self.shift:
try:
vcut = 4.0 * eps_ij * ((sig_ij / rcut)**12 -
(sig_ij / rcut)**6)
except ZeroDivisionError:
vcut = 0.0
self._offset[pair] = vcut
self.params[pair] = pair_param[pair]

def potential(self, system):
"""Calculate the potential energy for the Lennard-Jones interaction.

Parameters
----------
system : object like :py:class:.System
The system we are evaluating the potential in.

Returns
-------
v_pot : float
The potential energy.

"""
particles = system.particles
box = system.box
v_pot = ljfortran.potential(particles.pos,
box.length, box.ilength,
self._lj3, self._lj4, self._offset,
self._rcut2, particles.ptype,
particles.npart,
box.dim, self.ntype)
return v_pot

def force(self, system):
"""Calculate the force for the Lennard-Jones interaction.

We also calculate the virial here, since the force
is evaluated.

Parameters
----------
system : object like :py:class:.System
The system we are evaluating the force in.

Returns
-------
forces : numpy.array
The forces on the particles.
virial : numpy.array
The virial obtained from the forces.

"""
particles = system.particles
box = system.box
forces, virial = ljfortran.force(particles.pos,
box.length, box.ilength,
self._lj1, self._lj2, self._rcut2,
particles.ptype,
particles.npart,
box.dim, self.ntype)
return forces, virial

def potential_and_force(self, system):
"""Calculate potential and force for the Lennard-Jones interaction.

Since the force is evaluated, the virial is also calculated.

Parameters
----------
system : object like :py:class:.System
The system we are evaluating the potential and force in.

Note
----
Currently, the virial is only calculated for all the particles.
It is not calculated as a virial per atom. The virial
per atom might be useful to obtain a local pressure or stress,
however, this needs some consideration. Perhaps it's best to
fully implement this as a method of planes or something similar.

Returns
-------
out[0] : float
The potential energy as a float.
out[1] : numpy.array
The force as a numpy.array of the same shape as the
positions in particles.pos.
out[2] : numpy.array
The virial, as a symmetric matrix with dimensions
(dim, dim) where dim is given by the box/system dimensions.

"""
particles = system.particles
box = system.box
forces, virial, vpot = ljfortran.potential_and_force(particles.pos,
box.length,
box.ilength,
self._lj1,
self._lj2,
self._lj3,
self._lj4,
self._offset,
self._rcut2,
particles.ptype,
particles.npart,
box.dim,
self.ntype)
return vpot, forces, virial


### Writing a new potential function with C¶

We will now create a new potential function by writing the external library in C. This is done by the following steps:

#### Step 1. Creating the C code¶

The C code is more involved, and perhaps cumbersome compared to the FORTRAN code above. We are here explicitly assuming that our system is going to be 3D.

/* Copyright (c) 2021, PyRETIS Development Team.
#include <Python.h>
#include <numpy/arrayobject.h>
#include <math.h>

// Forward function declaration.
static PyObject *potential(PyObject *self, PyObject *args);
static PyObject *force(PyObject *self, PyObject *args);
static PyObject *force_and_pot(PyObject *self, PyObject *args);

// Boilerplate: function list.
static PyMethodDef ljmethods[] = {
{ "potential", potential, METH_VARARGS, "Calculate potential energy."},
{ "force", force, METH_VARARGS, "Calculate force."},
{ "force_and_pot", force_and_pot, METH_VARARGS, "Calculate force and potential energy."},
{ NULL, NULL, 0, NULL } /* Sentinel */
};

static struct PyModuleDef ljc =
{
"ljc",
"Compute Lennard-Jones interactions",
-1,
ljmethods
};

// Boilerplate: Module initialization.
PyMODINIT_FUNC PyInit_ljc(void) {
import_array();
return PyModule_Create(&ljc);
}

static inline npy_float64 pbc_dist(npy_float64 x,
npy_float64 box,
npy_float64 ibox) {
return x - nearbyint(x * ibox) * box;
}

static PyObject *potential(PyObject *self, PyObject *args) {
// Input variables:
npy_int64 npart, dim, ntype;
npy_int64 itype, loc;
PyArrayObject *py_pos, *py_box, *py_ibox, *py_lj3, *py_lj4, *py_rcut2, *py_offset, *py_ptype;
// Internal variables:
npy_float64 dx, dy, dz, rsq, r2inv, r6inv;
npy_int64 i, xi, yi, zi;
npy_int64 j, xj, yj, zj;

// Parse arguments.
if (!PyArg_ParseTuple(args, "O!O!O!O!O!O!O!O!lll",
&PyArray_Type, &py_pos,
&PyArray_Type, &py_box,
&PyArray_Type, &py_ibox,
&PyArray_Type, &py_lj3,
&PyArray_Type, &py_lj4,
&PyArray_Type, &py_offset,
&PyArray_Type, &py_rcut2,
&PyArray_Type, &py_ptype,
&npart, &dim, &ntype)){
return NULL;
}
// Get underlying arrays from numpy arrays.
npy_float64 *pos = (npy_float64*)PyArray_DATA(py_pos);
npy_float64 *box = (npy_float64*)PyArray_DATA(py_box);
npy_float64 *ibox = (npy_float64*)PyArray_DATA(py_ibox);
npy_float64 *lj3 = (npy_float64*)PyArray_DATA(py_lj3);
npy_float64 *lj4 = (npy_float64*)PyArray_DATA(py_lj4);
npy_float64 *rcut2 = (npy_float64*)PyArray_DATA(py_rcut2);
npy_float64 *offset = (npy_float64*)PyArray_DATA(py_offset);
npy_int64 *ptype = (npy_int64*)PyArray_DATA(py_ptype);
npy_float64 vpot = 0.0;
// Calculate potential:
for(i = 0; i < npart-1; i++){
xi = 3*i;
yi = 3*i + 1;
zi = 3*i + 2;
itype = ptype[i];
for (j = i + 1; j < npart; j++){
loc = itype + ntype * ptype[j];
xj = 3*j;
yj = 3*j + 1;
zj = 3*j + 2;
dx = pos[xi] - pos[xj];
dy = pos[yi] - pos[yj];
dz = pos[zi] - pos[zj];
dx = pbc_dist(dx, box[0], ibox[0]);
dy = pbc_dist(dy, box[1], ibox[1]);
dz = pbc_dist(dz, box[2], ibox[2]);
rsq = dx*dx + dy*dy + dz*dz;
if (rsq < rcut2[loc]){
r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv;
vpot += r6inv * (lj3[loc] * r6inv - lj4[loc]) - offset[loc];
}
}
}
/* Clean up */
/* Build the output tuple */
PyObject *ret = Py_BuildValue("d", vpot);
return ret;
}

static PyObject *force(PyObject *self, PyObject *args) {
// Input variables:
npy_int64 npart, dim, ntype;
npy_int64 itype, loc;
PyArrayObject *py_pos, *py_box, *py_ibox, *py_force, *py_virial;
PyArrayObject *py_lj1, *py_lj2, *py_rcut2, *py_ptype;
// Internal variables:
npy_float64 dx, dy, dz, rsq, r2inv, r6inv;
npy_float64 forcelj, fx, fy, fz;
npy_int64 i, xi, yi, zi;
npy_int64 j, xj, yj, zj;

// Parse arguments.
if (!PyArg_ParseTuple(args, "O!O!O!O!O!O!O!O!O!lll",
&PyArray_Type, &py_pos,
&PyArray_Type, &py_box,
&PyArray_Type, &py_ibox,
&PyArray_Type, &py_lj1,
&PyArray_Type, &py_lj2,
&PyArray_Type, &py_rcut2,
&PyArray_Type, &py_ptype,
&PyArray_Type, &py_force,
&PyArray_Type, &py_virial,
&npart, &dim, &ntype)){
return NULL;
}
// Get underlying arrays from numpy arrays.
npy_float64 *force = (npy_float64*)PyArray_DATA(py_force);
npy_float64 *virial = (npy_float64*)PyArray_DATA(py_virial);
npy_float64 *pos = (npy_float64*)PyArray_DATA(py_pos);
npy_float64 *box = (npy_float64*)PyArray_DATA(py_box);
npy_float64 *ibox = (npy_float64*)PyArray_DATA(py_ibox);
npy_float64 *lj1 = (npy_float64*)PyArray_DATA(py_lj1);
npy_float64 *lj2 = (npy_float64*)PyArray_DATA(py_lj2);
npy_float64 *rcut2 = (npy_float64*)PyArray_DATA(py_rcut2);
npy_int64 *ptype = (npy_int64*)PyArray_DATA(py_ptype);
// Create array for forces:
//PyObject *ret = PyArray_SimpleNew(dim, test, NPY_DOUBLE);
// Zero virial:
for(i = 0; i < dim*dim; i++) {
virial[i] = 0;
}
// Set all forces to zero.
for(i = 0; i < npart; i++) {
force[3*i] = 0;
force[3*i+1] = 0;
force[3*i+2] = 0;
}
// Calculate forces:
for(i = 0; i < npart-1; i++){
xi = 3*i;
yi = 3*i + 1;
zi = 3*i + 2;
itype = ptype[i];
for (j = i + 1; j < npart; j++){
loc = itype + ntype * ptype[j];
xj = 3*j;
yj = 3*j + 1;
zj = 3*j + 2;
dx = pos[xi] - pos[xj];
dy = pos[yi] - pos[yj];
dz = pos[zi] - pos[zj];
dx = pbc_dist(dx, box[0], ibox[0]);
dy = pbc_dist(dy, box[1], ibox[1]);
dz = pbc_dist(dz, box[2], ibox[2]);
rsq = dx*dx + dy*dy + dz*dz;
if (rsq < rcut2[loc]){
r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv;
forcelj = r2inv * r6inv * (lj1[loc] * r6inv - lj2[loc]);
fx = forcelj * dx;
fy = forcelj * dy;
fz = forcelj * dz;
force[xi] += fx;
force[yi] += fy;
force[zi] += fz;
force[xj] -= fx;
force[yj] -= fy;
force[zj] -= fz;
virial[0] += fx*dx;
virial[1] += fx*dy;
virial[2] += fx*dz;
virial[3] += fy*dx;
virial[4] += fy*dy;
virial[5] += fy*dz;
virial[6] += fz*dx;
virial[7] += fz*dy;
virial[8] += fz*dz;
}
}
}
Py_RETURN_NONE;
}

static PyObject *force_and_pot(PyObject *self, PyObject *args) {
// Input variables:
npy_int64 npart, dim, ntype;
npy_int64 itype, loc;
PyArrayObject *py_pos, *py_box, *py_ibox, *py_force, *py_virial;
PyArrayObject *py_lj1, *py_lj2, *py_lj3, *py_lj4, *py_rcut2, *py_offset, *py_ptype;
// Internal variables:
npy_float64 dx, dy, dz, rsq, r2inv, r6inv;
npy_float64 forcelj, fx, fy, fz;
npy_int64 i, xi, yi, zi;
npy_int64 j, xj, yj, zj;

// Parse arguments.
if (!PyArg_ParseTuple(args, "O!O!O!O!O!O!O!O!O!O!O!O!lll",
&PyArray_Type, &py_pos,
&PyArray_Type, &py_box,
&PyArray_Type, &py_ibox,
&PyArray_Type, &py_lj1,
&PyArray_Type, &py_lj2,
&PyArray_Type, &py_lj3,
&PyArray_Type, &py_lj4,
&PyArray_Type, &py_offset,
&PyArray_Type, &py_rcut2,
&PyArray_Type, &py_ptype,
&PyArray_Type, &py_force,
&PyArray_Type, &py_virial,
&npart, &dim, &ntype)){
return NULL;
}
// Get underlying arrays from numpy arrays.
npy_float64 *force = (npy_float64*)PyArray_DATA(py_force);
npy_float64 *virial = (npy_float64*)PyArray_DATA(py_virial);
npy_float64 *pos = (npy_float64*)PyArray_DATA(py_pos);
npy_float64 *box = (npy_float64*)PyArray_DATA(py_box);
npy_float64 *ibox = (npy_float64*)PyArray_DATA(py_ibox);
npy_float64 *lj1 = (npy_float64*)PyArray_DATA(py_lj1);
npy_float64 *lj2 = (npy_float64*)PyArray_DATA(py_lj2);
npy_float64 *lj3 = (npy_float64*)PyArray_DATA(py_lj3);
npy_float64 *lj4 = (npy_float64*)PyArray_DATA(py_lj4);
npy_float64 *offset = (npy_float64*)PyArray_DATA(py_offset);
npy_float64 *rcut2 = (npy_float64*)PyArray_DATA(py_rcut2);
npy_int64 *ptype = (npy_int64*)PyArray_DATA(py_ptype);
npy_float64 vpot = 0.0;
// Create array for forces:
//PyObject *ret = PyArray_SimpleNew(dim, test, NPY_DOUBLE);
// Zero virial:
for(i = 0; i < dim*dim; i++) {
virial[i] = 0;
}
// Set all forces to zero.
for(i = 0; i < npart; i++) {
force[3*i] = 0;
force[3*i+1] = 0;
force[3*i+2] = 0;
}
// Calculate forces:
for(i = 0; i < npart-1; i++){
xi = 3*i;
yi = 3*i + 1;
zi = 3*i + 2;
itype = ptype[i];
for (j = i + 1; j < npart; j++){
loc = itype + ntype * ptype[j];
xj = 3*j;
yj = 3*j + 1;
zj = 3*j + 2;
dx = pos[xi] - pos[xj];
dy = pos[yi] - pos[yj];
dz = pos[zi] - pos[zj];
dx = pbc_dist(dx, box[0], ibox[0]);
dy = pbc_dist(dy, box[1], ibox[1]);
dz = pbc_dist(dz, box[2], ibox[2]);
rsq = dx*dx + dy*dy + dz*dz;
if (rsq < rcut2[loc]){
r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv;
forcelj = r2inv * r6inv * (lj1[loc] * r6inv - lj2[loc]);
vpot += r6inv * (lj3[loc] * r6inv - lj4[loc]) - offset[loc];
fx = forcelj * dx;
fy = forcelj * dy;
fz = forcelj * dz;
force[xi] += fx;
force[yi] += fy;
force[zi] += fz;
force[xj] -= fx;
force[yj] -= fy;
force[zj] -= fz;
virial[0] += fx*dx;
virial[1] += fx*dy;
virial[2] += fx*dz;
virial[3] += fy*dx;
virial[4] += fy*dy;
virial[5] += fy*dz;
virial[6] += fz*dx;
virial[7] += fz*dy;
virial[8] += fz*dz;
}
}
}
PyObject *ret = Py_BuildValue("d", vpot);
return ret;
}


As can be seen in the C-code, there is some boilerplate code and we make use of PyArg_ParseTuple in order to parse the parameters to the function.

#### Step 2. Creating a setup.py file and compiling¶

The C-code can be compiled using a setup.py file.

# -*- coding: utf-8 -*-
# Copyright (c) 2021, PyRETIS Development Team.
"""A script for building the C extension."""
from distutils.core import (  # pylint: disable=import-error,no-name-in-module
setup,
Extension,
)
import numpy as np

LJMODULE = Extension(
'ljc',
sources=['ljc.c'],
extra_compile_args=['-Ofast', '-march=native'],
)
setup(
name='PyRETIS Lennard-Jones c extension',
description='C extension for the Lennard-Jones potential.',
ext_modules=[LJMODULE],
include_dirs=[np.get_include()],
)


The setup.py file is used to compile via the command

python setup.py build_ext --inplace


Here, build_ext is used to tell setup.py to compile the C extension and the --inplace will put the compiled extensions into the directory you have the source code in.

#### Step 3. Creating a new Python class for the potential function¶

The final step is to create a Python class which is making use of the C-code.

# -*- coding: utf-8 -*-
# Copyright (c) 2021, PyRETIS Development Team.
"""Example of using a Lennard-Jones potential implemented in C."""
import logging
import sys
import os
import numpy as np
from pyretis.forcefield.potentials import PairLennardJonesCut
from pyretis.forcefield.potentials.pairpotentials import (
generate_pair_interactions
)
# Just to handle imports of the library:
sys.path.insert(0, os.path.abspath(os.path.dirname(__file__)))
logger = logging.getLogger(__name__)  # pylint: disable=invalid-name
try:
import ljc
except ImportError:
MSG = ('Could not import external C library.'
'\nPlease compile with: "python setup.py build_ext --inplace"')
logger.critical(MSG)
raise ImportError(MSG)

__all__ = ['PairLennardJonesCutC']

class PairLennardJonesCutC(PairLennardJonesCut):
r"""class PairLennardJonesCutC(PairLennardJonesCut).

This class implements as simple Lennard-Jones 6-12 potential which
employs a simple cut-off and can be shifted. The potential energy
(:math:V_\text{pot}) is defined in the usual way for an
interacting pair of particles a distance :math:r apart,

.. math::

V_\text{pot} = 4 \varepsilon \left( x^{12} - x^{6} \right),

where :math:x = \sigma/r and :math:\varepsilon
and :math:\sigma are the potential parameters. The parameters are
stored as attributes of the potential and we store one set for each
kind of pair interaction. Parameters can be generated with a
specific mixing rule by the force field.

Attributes
----------
params : dict
The parameters for the potential. This dict is assumed to
contain parameters for pairs, i.e. for interactions.
_lj1 : numpy.array
Lennard-Jones parameters used for calculation of the force.
Keys are the pairs (particle types) that may interact.
Calculated as: 48.0 * epsilon * sigma**12
_lj2 : numpy.array
Lennard-Jones parameters used for calculation of the force.
Keys are the pairs (particle types) that may interact.
Calculated as: 24.0 * epsilon * sigma**6
_lj3 : numpy.array
Lennard-Jones parameters used for calculation of the potential.
Keys are the pairs (particle types) that may interact.
Calculated as: 4.0 * epsilon * sigma**12
_lj4 : numpy.array
Lennard-Jones parameters used for calculation of the potential.
Keys are the pairs (particle types) that may interact.
Calculated as: 4.0 * epsilon * sigma**6
_offset : numpy.array
Potential values for shifting the potential if requested.
This is the potential evaluated at the cutoff.
_rcut2 : numpy.array
The squared cut-off for each interaction type.
Keys are the pairs (particle types) that may interact.

"""

def __init__(self, dim=3, shift=True, mixing='geometric',
desc='Lennard-Jones pair potential (C)'):
"""Set up the Lennard-Jones potential.

Parameters
----------
dim : int
The dimensionality to use.
shift : boolean
Determines if the potential should be shifted or not.
mixing : string
Determines how we should mix potential parameters.

"""
super().__init__(dim=dim, desc=desc, mixing=mixing)
self.ntype = 0

def set_parameters(self, parameters):
"""Update all parameters.

Here, we generate pair interactions, since that is what this
potential actually is using.

Parameters
----------
parameters : dict
The input base parameters.

"""
self.params = {}
pair_param = generate_pair_interactions(parameters, self.mixing)
self.ntype = max(int(np.sqrt(len(pair_param))), 2)
self._lj1 = np.zeros((self.ntype, self.ntype))
self._lj2 = np.zeros_like(self._lj1)
self._lj3 = np.zeros_like(self._lj1)
self._lj4 = np.zeros_like(self._lj1)
self._rcut2 = np.zeros_like(self._lj1)
self._offset = np.zeros_like(self._lj1)
for pair in pair_param:
eps_ij = pair_param[pair]['epsilon']
sig_ij = pair_param[pair]['sigma']
rcut = pair_param[pair]['rcut']
self._lj1[pair] = 48.0 * eps_ij * sig_ij**12
self._lj2[pair] = 24.0 * eps_ij * sig_ij**6
self._lj3[pair] = 4.0 * eps_ij * sig_ij**12
self._lj4[pair] = 4.0 * eps_ij * sig_ij**6
self._rcut2[pair] = rcut**2
vcut = 0.0
if self.shift:
try:
vcut = 4.0 * eps_ij * ((sig_ij / rcut)**12 -
(sig_ij / rcut)**6)
except ZeroDivisionError:
vcut = 0.0
self._offset[pair] = vcut
self.params[pair] = pair_param[pair]

def potential(self, system):
"""Calculate the potential energy for the Lennard-Jones interaction.

Parameters
----------
system : object like :py:class:.System
The system we are evaluating the potential in.

Returns
-------
v_pot : float
The potential energy.

"""
particles = system.particles
box = system.box
v_pot = ljc.potential(particles.pos,
box.length, box.ilength,
self._lj3, self._lj4, self._offset,
self._rcut2, particles.ptype,
particles.npart,
box.dim, self.ntype)
return v_pot

def force(self, system):
"""Calculate the force for the Lennard-Jones interaction.

We also calculate the virial here, since the force
is evaluated.

Parameters
----------
system : object like :py:class:.System
The system we are evaluating the force in.

Returns
-------
forces : numpy.array
The forces on the particles.
virial : numpy.array
The virial obtained from the forces.

"""
particles = system.particles
box = system.box
forces = np.zeros((particles.npart, box.dim))
virial = np.zeros((box.dim, box.dim))
ljc.force(particles.pos,
box.length, box.ilength,
self._lj1, self._lj2, self._rcut2,
particles.ptype,
forces, virial,
particles.npart,
box.dim, self.ntype)
return forces, virial

def potential_and_force(self, system):
"""Calculate potential and force for the Lennard-Jones interaction.

Since the force is evaluated, the virial is also calculated.

Parameters
----------
system : object like :py:class:.System
The system we are evaluating the potential and force in.

Note
----
Currently, the virial is only calculated for all the particles.
It is not calculated as a virial per atom. The virial
per atom might be useful to obtain a local pressure or stress,
however, this needs some consideration. Perhaps it's best to
fully implement this as a method of planes or something similar.

Returns
-------
out[0] : float
The potential energy as a float
out[1] : numpy.array
The force as a numpy.array of the same shape as the
positions in particles.pos.
out[2] : numpy.array
The virial, as a symmetric matrix with dimensions
(dim, dim) where dim is given by the box/system dimensions.

"""
particles = system.particles
box = system.box
forces = np.zeros((particles.npart, box.dim))
virial = np.zeros((box.dim, box.dim))
vpot = ljc.force_and_pot(particles.pos,
box.length,
box.ilength,
self._lj1,
self._lj2,
self._lj3,
self._lj4,
self._offset,
self._rcut2,
particles.ptype,
forces,
virial,
particles.npart,
box.dim,
self.ntype)
return vpot, forces, virial


## Running the MD simulation using the PyRETIS library¶

Below we give a Python script which makes use of the PyRETIS library to run a MD simulation of 2000 steps in the forward and backward directions. It will generate a trajectory (named traj.xyz) and display a plot of the energies (an example is shown in the figure below).

Before running the script you will have to download the initial configuration initial.gro and place it in the same directory as the script. Further, the script makes the following assumptions:

• If you are using FORTRAN (that is if USE = fortran) the script assumes that there is a folder named fortran which contains the compiled library and the Python script named ljpotentialf.py.
• If you are using C (that is if USE = cpython3) the script assumes that there is a folder named cpython3 which contains the compiled library and the Python script named ljpotentialc.py.

# -*- coding: utf-8 -*-
# Copyright (c) 2021, PyRETIS Development Team.
"""
Example of running a MD NVE simulation.

Here we are running forward and then backward using a potential
written in FORTRAN or C. This can be selected below by the user.
"""

# pylint: disable=invalid-name
import sys
# For plotting:
from matplotlib import pyplot as plt
from matplotlib import gridspec
# PyRETIS imports:
from pyretis.core.units import create_conversion_factors
from pyretis.inout.setup import (create_simulation, create_force_field,
create_system, create_engine)
from pyretis.inout.fileio import FileIO
from pyretis.inout.formats import ThermoTableFormatter
from pyretis.inout.plotting import mpl_set_style

# Define simulation settings:

klass = {'fortran': 'PairLennardJonesCutF',
'cpython3': 'PairLennardJonesCutC'}
module = {'fortran': 'fortran/ljpotentialf.py',
'cpython3': 'cpython3/ljpotentialc.py'}

USE = 'cpython3'  # or 'fortran'

settings = {}
settings['simulation'] = {
'steps': 2000,
}
settings['system'] = {
'units': 'lj',
'temperature': 2.5,
'dimensions': 3
}
settings['engine'] = {
'class': 'velocityverlet',
'timestep': 0.002
}
settings['output'] = {
'backup': 'overwrite',
'energy-file': 1,
'screen': 10,
'trajectory-file': 1
}
settings['potential'] = [
{'class': klass[USE],
'module': module[USE],
'parameter': {0: {'sigma': 1, 'epsilon': 1, 'factor': 2.5},
1: {'sigma': 1, 'epsilon': 1, 'factor': 2.5}},
'shift': True}
]
settings['particles'] = {
'position': {'file': 'initial.gro'},
'velocity': {'generate': 'maxwell', 'momentum': True, 'seed': 0}
}

create_conversion_factors(settings['system']['units'])
print('# Creating system from settings.')
ljsystem = create_system(settings)
ljsystem.forcefield = create_force_field(settings)
kwargs = {'system': ljsystem, 'engine': create_engine(settings)}
simulation_nve = create_simulation(settings, kwargs)

# Set up extra output:
thermo_file = FileIO('thermo.txt', 'w', ThermoTableFormatter(), backup=False)
thermo_file.open()
store_results = []
# Also create some other outputs:
simulation_nve.set_up_output(settings, progress=False)
# Run the simulation forward:
for result in simulation_nve.run():
stepno = result['cycle']['stepno']
thermo_file.output(stepno, result['thermo'])
result['thermo']['stepno'] = stepno
store_results.append(result['thermo'])
# Run backward:
ljsystem.particles.vel *= -1.0
simulation_nve.extend_cycles(settings['simulation']['steps'] - 1)
for result in simulation_nve.run():
stepno = result['cycle']['stepno']
thermo_file.output(stepno, result['thermo'])
result['thermo']['stepno'] = stepno
store_results.append(result['thermo'])
thermo_file.close()

mpl_set_style()  # Load PyRETIS plotting style
step = [res['stepno'] for res in store_results]
pot_e = [res['vpot'] for res in store_results]
kin_e = [res['ekin'] for res in store_results]
tot_e = [res['etot'] for res in store_results]
pressure = [res['press'] for res in store_results]
temp = [res['temp'] for res in store_results]
# Plot some energies:
fig1 = plt.figure()
gs = gridspec.GridSpec(2, 2)
ax1.plot(step, pot_e, label='Potential')
ax1.plot(step, kin_e, label='Kinetic')
ax1.plot(step, tot_e, label='Total')
ax1.set_xlabel('Step no.')
ax1.set_ylabel('Energy per particle')
ax1.legend(loc='center left', prop={'size': 'small'})
ax2.plot(step, temp)
ax2.set_ylabel('Temperature')
ax3.plot(step, pressure)
ax3.set_xlabel('Step no.')
ax3.set_ylabel('Pressure')

Copy this script and store it in a new Python file, say md_forward_backward_py. Before running the script, you will have to set the USE variable in the script to either USE = fortran or USE = cpython3 to select the external library to use. After this has been done, you can execute the script as usual
python md_forward_backward.py

This will run the simulation, and when it is done, it should display a figure similar to the one given below. Further, you may wish to inspect the generated trajectory traj.xyz visually.