Introduction to the PyRETIS library

In this introduction to the PyRETIS library the main classes and functions from the PyRETIS library will be discussed. For the full description, we refer to the API documentation.

The PyRETIS library contains methods and classes that handle the different aspects of a simulation. These are grouped into sub-packages and the different sub-packages and the classes and methods defined within them will interact. As an example, assume that we are performing a RETIS simulation. In terms of objects from the PyRETIS library, this simulation can be described as follows:

  1. We first define the System we are studying. This contains information about the Particles, Box and the ForceField.
  2. Next, an OrderParameter is defined and the order parameter it is representing can be calculated using a System object as the argument.
  3. The RETIS simulation is handled by the SimulationRETIS class which will use a specific EngineBase like object in order to generate several Path objects for a collection of PathEnsemble objects. This generation is done by a set of methods defined in the modules pyretis.core.tis and pyretis.core.retis.
  4. An analysis is carried out by making use of methods from the analysis sub-package.

The core and forcefield sub-packages

The pyretis.core library defines the core method and classes and here we will introduce the

  • System class which defines the system we are investigating. This class is actually composed of several other objects, which we will also discuss here:
    • Particles: A class which represents the particles.
    • Box: A class which represents the simulation box.
    • ForceField: A class representing the force field. (Note: This is defined within pyretis.forcefield).
  • Path class which defines paths/trajectories.
  • PathEnsemble class which defines path ensembles.

Below is an illustration on how some of these classes are interacting. The classes (shown as boxes in this figure) will be discusses more in the following.

Illustration of core and force field classes

Fig. 10 Illustration of the relations between classes from the core and force field sub-packages. The classes as shown as green boxes where the class names are shown with white text. Examples of class methods are written as “method()” while and class attributes are indicated as just “attribute”. Some of the attributes are in fact references to other classes and these are highlighted with grey boxes and the reference is shown by arrows. As can be seen in this figure, the System is in practice composed of several other classes.

The System class

The System class defines the system we are investigating. It will typically contain particles, a simulation box and a force field. For PyRETIS, this is a very convenient class since it exposes important parts we can interact with, in particular the particles.

Example of creation:

from pyretis.core import System
new_system = System(temperature=0.8, units='lj')

This will create an empty system with a set temperature equal to 0.8 in lj units (lj refers to Lennard-Jones units. It is also possible to specify a box here in case that it needed:

new_system = System(temperature=0.8, units='lj', box=mybox)

where mybox can be created as described below.

Particles can be added by first creating an object as described below. A short example:

from pyretis.core import System, Particles
new_system = System(temperature=0.8, units='lj')
new_system.particles = Particles()
new_system.add_particle([0.0, 0.0, 0.0], mass=1.0, name='Ar', ptype=0)

Here, we are setting the System.particles attribute and using the System.add_particle() to add a particle with a given position, mass, name and type.

The Particles class

The Particles class represents a collection of particles and in many ways it can be viewed as a particle list. Internally in PyRETIS, the particle list is one of the most important classes. The positions, velocities and forces are accessed through an instance of this class using the class attributes Particles.pos, Particles.vel and Particles.force. Actually, there is an additional particle class within PyRETIS which is called ParticlesExt. This class is used when PyRETIS is using external engines. It is very similar to the Particles class but it has in addition a ParticlesExt.config which can be used to reference files which holds the current configuration of the particles.

Here are some examples of interacting with the Particles class, using Particles.add_particle() to add some particles. Actually, the system will make use of this method if you are calling System.add_particle().

from pyretis.core import Particles

part = Particles(dim=3)

pos = [0.0, 1.0, 0.0]
vel = [0.0, 0.0, 0.0]
force = [0.0, 0.0, 0.0]
part.add_particle(pos, vel, force, mass=1.0, name='Ar', ptype=0)
print(part.pos)
pos = [1.0, 0.0, 0.0]
part.add_particle(pos, vel, force, mass=1.0, name='Ar', ptype=0)
print(part.pos)

Here, we can add names to particles using the keyword name and we can also specify a particle type using ptype. The name can be used to identify/tag specific particles and is used for output purposes. Internally, the particle type is more important: The particle type can be used to specify parameters for pair interactions which is computed by the force field.

When we initiate an instance of Particles, we define the dimensionality using the dim keyword parameter.

The Box class

The Box class defines a simulation box. It is useful in simulations where we wish to have periodic boundaries. Typically, we do not interact much with the box beyond creating it. Boxes are created by passing a (optional) length which is a list of integers either of type [lengthx, lengthy, lengthz]. At the same time periodicity can be specified with the keyword periodic which is a list of boolean values that determine if a dimension is periodic or not. The default is periodic in all directions.

Some examples:

from pyretis.core import create_box
box1 = create_box()
print(box1)
box2 = create_box(length=[10, 10, 10])
print(box2)
box3 = create_box(low=[0, -10, 10], high=[10, 10, 20],
                  periodic=[True, True, False])
print(box3)

The ForceField class

The ForceField class is used to define force fields. A force field is just a list of functions (and parameters) which can be used to obtain the force and potential energy. In general, the force field expect that its constituent potential functions actually supports calling three functions which means that the potential functions must be slightly more complex than just simple functions — they need to be classes which subclass the PotentialFunction class.

If we, for the sake of an example, let an instance of the ForceField class have a constituent potential function named func, then PyRETIS will assume that it can call:

  1. func.potential(system) to obtain the potential energy.
  2. func.force(system) to obtain the forces and the virial.
  3. func.potential_and_force(system) to obtain the potential energy, forces and the virial. Typically, this can be done by just calling func.potential(system) and func.force(system).

Notice that all this functions should only take in a System as the only parameter.

Let’s see an example of how we can set-up a potential function (or class) and add it to a force field.

First, define the potential function using:

from matplotlib import pyplot as plt
import numpy as np
from pyretis.core import System, Particles
from pyretis.forcefield import ForceField
from pyretis.forcefield.potential import PotentialFunction


class Harmonic1D(PotentialFunction):
    """A 1D harmonic potential function."""

    def __init__(self):
        """Set up the potential."""
        super().__init__(dim=1, desc='1D harmonic potential')
        self.eq_pos = 0.0  # equilibrium position
        self.k_force = 1.0  # force constant

    def potential(self, system):
        """Calculate potential energy."""
        pos = system.particles.pos  # Get positions from the particle list.
        vpot = 0.5 * self.k_force * (pos - self.eq_pos)**2
        return vpot.sum()

    def force(self, system):
        """Calculate the force."""
        pos = system.particles.pos  # Get positions from the particle list.
        virial = None  # We are lazy and do not calculate the virial here.
        forces = -self.k_force * (pos - self.eq_pos)
        return forces, virial

    def potential_and_force(self, system):
        """Calculate force and potential."""
        pot = self.potential(system)
        force, virial = self.force(system)
        return pot, force, virial

Using what we have already discussed above about the System and Particles we can make a plot of the potential by adding:

if __name__ == '__main__':
    # Create empty force field:
    forcefield = ForceField(desc='1D Harmonic potential')
    pot = Harmonic1D()  # create potential
    # Add the potential function to the force field:
    forcefield.add_potential(pot)
    system = System()
    system.particles = Particles(dim=1)
    system.add_particle(pos=np.zeros(1))
    # Do some plotting, first calculate the potential as some locations:
    vpot = []
    pos = np.linspace(-2.5, 2.5, 100)
    for xi in pos:
        system.particles.pos = xi
        vpot.append(forcefield.evaluate_potential(system))
    # And plot it using matplotlib:
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ax1.plot(pos, vpot, lw=3)
    plt.show()

The Path and PathEnsemble classes

These two classes are representation of paths and path ensembles. Path are essentially trajectories and path ensembles are collections of such paths. The PathEnsemble stores some information about all paths in the attribute PathEnsemble.paths, but it will not store the full trajectory (positions, velocities, etc.). It will, however, keep a reference to the last accepted path in PathEnsemble.last_path. This Path object can for instance be inspected by using the Path.trajectory() method.

Illustration of the relation PathEnsemble->Path

Fig. 11 Illustration of the relation between the PathEnsemble and the Path.

We give here a short example on how you can interact with a path. A complicating factor here is that we need to create a random generator to use with the Path. This is done by using the RandomGenerator class and the reader is referred to the API documentation for information about this.

import numpy as np
from pyretis.core.path import Path
from pyretis.core.random_gen import RandomGenerator

rgen = RandomGenerator(seed=0)  # Create random generator
path = Path(rgen)  # Create empty path

# Add some phase points to the path:
# (A phase point is just the order parameter, the positions and
# velocities and energies.)
for i in range(10):
    phasepoint = {'order': [i],
                  'pos': np.zeros(3),
                  'vel': np.zeros(3),
                  'vpot': i,
                  'ekin': i}
    path.append(phasepoint)

# Loop over the path using .trajectory()
print('Looping forward:')
for phasepoint in path.trajectory():
    print(phasepoint['vpot'])

# Loop over the path using .trajectory()
print('Looping backward:')
for phasepoint in path.trajectory(reverse=True):
    print(phasepoint['vpot'])

# Get some randomly chosen shooting points:
print('Getting shooting points:')
for i in range(10):
    point, idx = path.get_shooting_point()
    print('Shooting point {}: index = {}, vpot = {}'.format(i, idx,
                                                            point['vpot']))

The simulation and engines sub-packages

The pyretis.simulation sub-package defines classes which are used to set-up and define different types of simulations. Typically, such simulations will need to interact with and change the state of a given System. This interaction is carried out by a particular EngineBase from the pyretis.engines sub-package. The interaction between these classes are illustrated in the figure below:

Illustration the simulation class.

Fig. 12 Illustration of the relations between the Simulation, EngineBase and System. A simulation object will typically contain references to a system object and to an engine object. The simulation can then use the engine in order to interact with the system. For instance can the Simulation.step() methods use the EngineBase.integration_step() in order to update positions, velocities and forces in the system.

In this section, we will not give a complete example on how to create new simulation class or a new engine class. We refer the reader to the examples, in particular:

The Simulation class

The Simulation class defines a simulation we can run. A simulation will typically act on a System and alter its state. We will here just describe the generic base class Simulation and we refer the reader to the extended pyretis API documentation for information about specific simulation classes, for instance the SimulationRETIS class. The most commonly used methods from the Simulation are:

  • run() which will run a simulation and for each step it will yield a dictionary with results.
  • step() which will run a single step from the simulation and return a dictionary with results.
  • is_finished() which will return True if the simulation has ended.
  • add_task() which can be used to add simulation tasks to a generic simulation.

Example of creation:

from pyretis.simulation import Simulation
new_simulation = Simulation(startcycle=0, endcycle=100)

for step in new_simulation.run():
    print(step['cycle']['stepno'])

The code block above will create a generic simulation object and run it. This simulation is not doing anything useful, it will only increment the current simulation step number from the given startcycle to the given endcycle. In order to something more productive, we can attach tasks to the simulation. This can be done as follows:

from pyretis.core import System
from pyretis.simulation import Simulation


def function(sys):
    """Return the set temperature for the system."""
    return sys.temperature


simulation = Simulation(steps=10)
system = System(temperature=300)

my_task = {'func': function,
           'args': [system],
           'first': True,
           'result': 'temperature'}

simulation.add_task(my_task)

for result in simulation.run():
    step = result['cycle']['step']
    temp = result['temperature']['set']
    print(step, temp)

As you can see, a new task is added by defining it as a dictionary:

my_task = {'func': function,
           'args': [system],
           'first': True,
           'result': 'temperature'}

The following keywords are used:

  • func which is a reference to a function to use.
  • args which is a list of arguments that should be given to the function.
  • first whether the task should be executed on the first simulation step (i.e. step 0).
  • result a string which labels the result in the dictionary returned by the methods Simulation.run() or Simulation.step().

Typically, when creating custom simulation, you will rewrite the methods Simulation.run() and Simulation.step() to fit the custom simulation you are going to perform, rather than adding tasks. However, for interactive work, short examples etc., the Simulation.add_task() can be useful.

The EngineBase class

The EngineBase is a base class defining a generic engine. In PyRETIS engines are either internal or external. External engines, e.g. ExternalMDEngine, interfaces external programs while internal engines, e.g. MDEngine interact directly with a System object. Creating an external engine may be somewhat involved depending on the program you wish to interface. A description on how to create new external engines can be found in the user guide.

The orderparameter sub-package

The pyretis.orderparameter package defines order parameters to use for path sampling simulation. In PyRETIS, such order parameters are assumed to be functions of a System only. Here, we give some examples of how a generic order parameter can be used. For information on how custom order parameters can be created, we refer to the detailed description in the user guide.

The OrderParameter class

The most important piece of the OrderParameter class is the actual calculation of the order parameter. This should be defined in a method like OrderParameter.calculate(). Here, it is assumed that the order parameter takes in an object like System. Since this is described elsewhere we will here just describe the usage of:

Here is an example on how these can be used. First we just set-up and define a test-system and we create an order parameter, for simplicity we use the pre-defined OrderParameterPosition class:

import numpy as np
from pyretis.core import System, Particles
from pyretis.orderparameter import OrderParameterPosition

system = System(temperature=1, units='reduced')
system.particles = Particles(dim=3)
system.add_particle([1.0, 2.0, 3.0], mass=1.0, name='Ar', ptype=0)

position = OrderParameterPosition(0, dim='x')

Next, we can calculate the order parameter as follows:

print('Order =', position.calculate(system))
print('Order =', position.calculate_all(system))

Note that both these lines return two values. This is because the OrderParameterPosition class also calculates the time derivative of the order parameter when calculate() or calculate_all() is called.

To show how we can add additional variables, try the following:

# Let's define additional collective variables:
def collective_y(syst):
    """Position along y-axis."""
    return [syst.particles.pos[0][1]]


def collective_2(syst):
    """Additional collective variable."""
    return [np.cos(np.pi * syst.particles.pos[0][2])]


position.add_orderparameter(collective_y)
position.add_orderparameter(collective_2)
print('Order =', position.calculate(system))
print('Order =', position.calculate_all(system))

As demonstrated by this example, one can add additional order parameters on the fly if such are needed. Note also the difference between calculate() and calculate_all(): The latter will also calculate order parameters added using add_orderparameter() while the former only calculate the order parameters defined within the calculate() method.

Important

The first order parameter returned from calculate() is taken as the progress coordinate used in path sampling simulations.

The analysis sub-package

The pyretis.analysis sub-package defines tools which are used in the analysis of simulation output. It defines methods for running predefined analysis tasks, e.g. retis_flux(), but also generic analysis methods, e.g. running_average(). Here, we refer the reader to the pyretis API documentation for more information about these methods.

The inout sub-package

The pyretis.inout sub-package contains methods which PyRETIS make use of in order to interact with you. These are ways to read input from you and to present output to you. This sub-package is relatively large and it is in fact made up by the following sub-packages:

  • pyretis.inout.analysisio which handles the input and output needed for analysis.
  • pyretis.inout.plotting which handles plotting of figures. It defines simple things like colors etc. for plotting. It also defines functions which can be used for specific plotting by the analysis and report tools.
  • pyretis.inout.report which is used to generate reports with resuts from different simulations.
  • pyretis.inout.setup which handles creation of objects from simulation settings.
  • pyretis.inout.writers for formatting and presenting text based output.

Again, we refer the reader to the pyretis API documentation for more information about these sub-packages.

The initiation sub-package

The pyretis.initiation sub-package contains methods to initialize path ensembles. We refer the reader to the pyretis API documentation for more information about this sub-package.

The tools sub-package

The tools library can be used to generate initial structures for a simulation. In the tools library the function generate_lattice() is defined and it supports the creation of the following lattices where the short hand keywords (sc, sq etc.) are used to select a specific lattice:

  • sc: A simple cubic lattice .
  • sq: Square lattice (2D) with one atom in the unit cell.
  • sq2: Square lattice with two atoms in the unit cell.
  • bcc: Body-centered cubic lattice.
  • fcc: Face-centered cubic lattice.
  • hcp: Hexagonal close-packed lattice.
  • diamond: Diamond structure.

A lattice is generated in the following way:

from pyretis.tools import generate_lattice
xyz, size = generate_lattice('diamond', [1, 1, 1], lcon=3.57)

Where the first parameter selects the lattice type, the second parameter selects give the number of repetitions in each direction and the optional parameter lcon is the lattice constant. The returned values are xyz – the coordinates – and size which is the bounding box of the generated structure. This variable can be used to define a simulation box. It is also possible to specify a number density:

from pyretis.tools import generate_lattice
xyz, size = generate_lattice('diamond', [1, 1, 1], density=1)

If the density is specified, the lattice constant lcon is deduced:

\text{lcon} = \left(\frac{n}{\rho}\right)^{\frac{1}{d}},

where n is the number of particles in the unit cell, \rho the specified number density and d the dimensionality. If we wish to save a generated lattice, this can be done as follows

from pyretis.tools import generate_lattice
from pyretis.inout.writers.xyzio import write_xyz_file
xyz, size = generate_lattice('diamond', [3, 3, 3], lcon=3.57)
write_xyz_file('diamond.xyz', xyz, names=['C']*len(xyz), header='Diamond')

xyz, size = generate_lattice('hcp', [3, 3, 3], lcon=2.5)
name = ['A', 'B'] * (len(xyz) // 2)
write_xyz_file('hcp.xyz', xyz, names=name, header='HCP test')

xyz, size = generate_lattice('sq2', [3, 3], lcon=1.0)
write_xyz_file('sq2.xyz', xyz, header='sq2 test')

Some examples of using the PyRETIS library

Here, we show some examples on how we can perform some common tasks using the PyRETIS library.

Reversing a trajectory

PyRETIS will not reverse backward trajectories during a simulation if you are using an external engine. For visualization purposes, it is very helpfull to reverse these trajectories before viewing them. This can be accomplished with the PyRETIS library as follows:

  • For GROMACS .trr trajectories:

    from pyretis.inout.writers.gromacsio import reverse_trr
    reverse_trr('trajB.trr', 'rev-trajB.trr')
    

    which will read the trajectory trajB.trr and store it as rev-trajB.trr.

  • For .xyz trajectories:

    from pyretis.inout.writers.xyzio import reverse_xyz_file
    reverse_xyz_file('trajB.xyz', 'rev-trajB.xyz')
    

    which will read the trajectory trajB.xyz and store it as rev-trajB.xyz.

Recalculating order parameters

If you have an existing trajectory and with to calculate order parameters for each step in the trajectory, this can be accomplished by making use of the pyretis.tools.recalculate_order module. For example, if you have a trajectory named traj.trr and an order parameter defined in an input file retis.rst this can be done as follows:

from pyretis.inout.writers import get_writer
from pyretis.inout.settings import parse_settings_file
from pyretis.inout.setup import create_orderparameter
from pyretis.tools.recalculate_order import recalculate_order

settings = parse_settings_file('retis.rst')
order_parameter = create_orderparameter(settings)
order = recalculate_order(order_parameter, 'traj.trr', reverse=False,
                          maxidx=None, minidx=None)
writer = get_writer('order')
with open('order.txt', 'w') as outfile:
    outfile.write('{}\n'.format(writer.header))
    for step, data in enumerate(order):
        txt = writer.format_data(step, data)
        outfile.write('{}\n'.format(txt))

This will create a new file order.txt with the re-calculated order parameters. The keyword reverse can be used to tell PyRETIS that you are looking at a time-reversed trajectory. The keyword minidx gives the frame number from which we start calculating (None is equal to reading from the first frame) and the keyword maxidx gives the last frame number we will read (setting it to None will read until the end of the trajectory file).