# RETIS in a 2D potential¶

In this example, you will explore a rare event with the Replica Exchange Transition Interface Sampling (RETIS) algorithm.

In this example, we will implement a new potential for PyRETIS and we will create three new order parameters which we will use to investigate a transition between two stable states in the potential. Using the RETIS algorithm, we will compute the rate constant for the transition between the two stable states.

This example is structured as follows:

## Creating the potential function¶

We will consider a relatively simple 2D potential where a single particle is moving. The potential energy () is given by

where and gives the positions. The potential energy is shown in Fig. 25 below.

Fig. 25 The potential energy as a function of the position. We have two stable states, at (x, y) = (-0.2, -0.4) and (x, y) = (0.2, 0.4) (indicated with the circles), separated by a saddle point at (x, y) = (0, 0).

We begin by creating a new potential function to use with PyRETIS, see the user guide for a short introduction on how custom potential functions can be created. In short, we will now have to define a new class which calculates the potential energy and the corresponding force.

We begin by simply creating a new class for the new potential function. Here, we define the parameters our new potential should accept. We implement the potential slightly more general than the potential given above so that we can change the parameters more easily, if we wish to do so, later. The potential we will be implementing is given by:

where , , , , , , , and are potential parameters.

Next, we define a new method as part of our new potential class which will be responsible for calculating the potential energy and a method responsible for calculating the force. Finally, we add a method that will calculate both the potential and the force at the same time. The full Python code for the new potential class can be found below.

# -*- coding: utf-8 -*-
# Copyright (c) 2019, PyRETIS Development Team.
"""This is a 2D example potential."""
import logging
import numpy as np
from pyretis.forcefield.potential import PotentialFunction
logger = logging.getLogger(__name__)  # pylint: disable=invalid-name

class Hyst2D(PotentialFunction):
r"""Hyst2D(PotentialFunction).

This class defines a 2D dimensional potential with two
stable states.
The potential energy (:math:V_\text{pot}) is given by

.. math::

V_\text{pot}(x, y) = \gamma_1 (x^2 + y^2)^2 +
\gamma_2 \exp(\alpha_1 (x - x_0)^2 + \alpha_2 (y - y_0)^2) +
\gamma_3 \exp(\beta_1 (x + x_0)^2 + \beta_2(y + y_0)^2)

where :math:x and :math:y gives the positions and :math:\gamma_1,
:math:\gamma_2, :math:\gamma_3, :math:\alpha_1, :math:\alpha_2,
:math:\beta_1, :math:\beta_2, :math:x_0 and :math:y_0 are
potential parameters.

Attributes
----------
params : dict
Contains the parameters. The keys are:

* gamma1: The :math:\gamma_1 parameter for the potential.
* gamma2: The :math:\gamma_2 parameter for the potential.
* gamma3: The :math:\gamma_3 parameter for the potential.
* alpha1: The :math:\alpha_1 parameter for the potential.
* alpha2: The :math:\alpha_2 parameter for the potential.
* beta1: The :math:\beta_1 parameter for the potential.
* beta2: The :math:\beta_2 parameter for the potential.
* x0: The :math:x_0 parameter for the potential.
* y0: The :math:y_0 parameter for the potential.

"""

def __init__(self, desc='2D hysteresis'):
"""Set up the potential.

Parameters
----------
a : float, optional
Parameter for the potential.
b : float, optional
Parameter for the potential.
c : float, optional
Parameter for the potential.
desc : string, optional
Description of the force field.

"""
super().__init__(dim=2, desc=desc)
self.params = {'gamma1': 0.0, 'gamma2': 0.0, 'gamma3': 0.0,
'alpha1': 0.0, 'alpha2': 0.0, 'beta1': 0.0,
'beta2': 0.0, 'x0': 0.0, 'y0': 0.0}

def potential(self, system):
"""Evaluate the potential.

Parameters
----------
system : object like System
The system we evaluate the potential for. Here, we
make use of the positions only.

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

"""
x = system.particles.pos[:, 0]  # pylint: disable=invalid-name
y = system.particles.pos[:, 1]  # pylint: disable=invalid-name
gam1 = self.params['gamma1']
gam2 = self.params['gamma2']
gam3 = self.params['gamma3']
alf1 = self.params['alpha1']
alf2 = self.params['alpha2']
bet1 = self.params['beta1']
bet2 = self.params['beta2']
x_0 = self.params['x0']
y_0 = self.params['y0']
v_pot = (gam1 * (x**2 + y**2)**2 +
gam2 * np.exp(alf1 * (x - x_0)**2 + alf2 * (y - y_0)**2) +
gam3 * np.exp(bet1 * (x + x_0)**2 + bet2 * (y + y_0)**2))
return v_pot.sum()

def force(self, system):
"""Evaluate forces.

Parameters
----------
system : object like System
The system we evaluate the potential for. Here, we
make use of the positions only.

Returns
-------
out[0] : numpy.array
The calculated force.
out[1] : numpy.array
The virial, currently not implemented for this potential.

"""
x = system.particles.pos[:, 0]  # pylint: disable=invalid-name
y = system.particles.pos[:, 1]  # pylint: disable=invalid-name
gam1 = self.params['gamma1']
gam2 = self.params['gamma2']
gam3 = self.params['gamma3']
alf1 = self.params['alpha1']
alf2 = self.params['alpha2']
bet1 = self.params['beta1']
bet2 = self.params['beta2']
x_0 = self.params['x0']
y_0 = self.params['y0']
term = 4.0 * gam1 * (x**2 + y**2)
exp1 = gam2 * np.exp(alf1 * (x - x_0)**2 + alf2 * (y - y_0)**2)
exp2 = gam3 * np.exp(bet1 * (x + x_0)**2 + bet2 * (y + y_0)**2)
forces = np.zeros_like(system.particles.pos)
forces[:, 0] = -(x * term +
2.0 * alf1 * (x - x_0) * exp1 +
2.0 * bet1 * (x + x_0) * exp2)
forces[:, 1] = -(y * term +
2.0 * alf2 * (y - y_0) * exp1 +
2.0 * bet2 * (y + y_0) * exp2)
virial = np.zeros((self.dim, self.dim))  # just return zeros here
return forces, virial

def potential_and_force(self, system):
"""Evaluate the potential and the force.

Parameters
----------
system : object like System
The system we evaluate the potential for. Here, we
make use of the positions only.

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, currently not implemented for this potential.

"""
virial = np.zeros((self.dim, self.dim))  # just return zeros here
x = system.particles.pos[:, 0]  # pylint: disable=invalid-name
y = system.particles.pos[:, 1]  # pylint: disable=invalid-name
gam1 = self.params['gamma1']
gam2 = self.params['gamma2']
gam3 = self.params['gamma3']
alf1 = self.params['alpha1']
alf2 = self.params['alpha2']
bet1 = self.params['beta1']
bet2 = self.params['beta2']
x_0 = self.params['x0']
y_0 = self.params['y0']
term0 = (x**2 + y**2)
exp1 = gam2 * np.exp(alf1 * (x - x_0)**2 + alf2 * (y - y_0)**2)
exp2 = gam3 * np.exp(bet1 * (x + x_0)**2 + bet2 * (y + y_0)**2)
v_pot = gam1 * term0**2 + exp1 + exp2
term = 4.0 * gam1 * term0
forces = np.zeros_like(system.particles.pos)
forces[:, 0] = -(x * term +
2.0 * alf1 * (x - x_0) * exp1 +
2.0 * bet1 * (x + x_0) * exp2)
forces[:, 1] = -(y * term +
2.0 * alf2 * (y - y_0) * exp1 +
2.0 * bet2 * (y + y_0) * exp2)
return v_pot.sum(), forces, virial


This new potential function can be included in a PyRETIS simulation by referencing it in the input file as follows (assuming that it has been stored in a file named potential.py):

Forcefield settings
-------------------
description = 2D hysteresis

Potential
---------
class = Hyst2D
module = potential.py
parameter gamma1 = 1
parameter gamma2 = -10
parameter gamma3 = -10
parameter alpha1 = -30
parameter alpha2 = -3
parameter beta1 = -30
parameter beta2 = -3
parameter x0 = 0.2
parameter y0 = 0.4


## Creating the order parameters¶

In this case, it is not so obvious how to define the order parameter. Three simple possibilities are (see Fig. 26 for an illustration):

1. The x coordinate of the particle, where we include the potential energy in order to define the two stable states.
2. The y coordinate of the particle, where we include the potential energy in order to define the two stable states.
3. The projection of the distance vector from the position of the particle (x, y) to (-0.2, -0.4) onto the straight line between (-0.2, -0.4) and (0.2, 0.4). Here, we also include the potential energy in order to define the two stable states.

Fig. 26 Illustration of the three different order parameters with location of interfaces (dotted lines). From left to right: using the x position as the order parameter, using the y position as the order parameter and, finally, in the rightmost figure, using a combination of x and y as the order parameter. The white line in the rightmost figure show the line which we project the order parameter onto. The blue paths are valid paths for the final path ensembles and in the bottom row, the corresponding order parameters are shown for these paths. The light dotted interface lines represent the locations where we include the potential energy in the definition of the order parameters. The two marked points (cross and triangle) in the paths in the bottom row show the second last and last points in the paths, respectively. The two marked points (cross and circle) in the paths in the top row show the first and last points in the paths, respectively.

Since we include the potential energy in the order parameter definition, we will have to create customized order parameters. There is a short introduction to how the order parameters are implemented in PyRETIS in the user guide and there is also a more lengthy recipe for creating custom order parameters available.

Here, we will first consider the x (or y) order parameter and then the more involved combination.

### Creating the x (or y) positional order parameter¶

The order parameter we are now creating will just be the position of the particle with two extra conditions:

1. If the position is below a certain value (inter_a) close to the first interface and if the potential energy is above a certain threshold value (energy_a), then the order parameter will not be allowed to come closer to the first interface.
2. If the position is closer to the last interface than a certain value (inter_b) and if the potential energy is above a certain threshold value (energy_b), then the order parameter will not be allowed to come closer to the last interface.

Thus, in addition to the position of the particle, we need to consider the 4 additional parameters inter_a, inter_b, energy_a, energy_b. The Python code for this new order parameter is given below:

# -*- coding: utf-8 -*-
# Copyright (c) 2019, PyRETIS Development Team.
"""The order parameter for the hysteresis example."""
import logging
from pyretis.orderparameter.orderparameter import OrderParameter
logger = logging.getLogger(__name__)  # pylint: disable=invalid-name

class OrderX(OrderParameter):
"""A positional order parameter.

Order parameter for the hysteresis example. In addition to using
the position, we also use the energy to tell if we are in states A/B.

Attributes
----------
index : integer
This is the index of the atom which will be used, i.e.
system.particles.pos[index] will be used.
inter_a : float
An interface such that we are in state A for postions < inter_a.
inter_b : float
An interface such that we are in state B for postions > inter_b.
energy_a : float
An energy such that we are in state A for potential energy < energy_a.
energy_b : float
An energy such that we are in state A for potential energy < energy_b.
dim : integer
This is the dimension of the coordinate to use.
0, 1 or 2 for 'x', 'y' or 'z'.
periodic : boolean
This determines if periodic boundaries should be applied to
the position or not.

"""

def __init__(self, index, inter_a, inter_b, energy_a, energy_b,
dim='x', periodic=False):
"""Initialise the order parameter.

Parameters
----------
index : int
This is the index of the atom we will use the position of.
inter_a : float
An interface such that we are in state A for postions < inter_a.
inter_b : float
An interface such that we are in state B for postions > inter_b.
energy_a : float
An energy such that we are in state A for
potential energy < energy_a.
energy_b : float
An energy such that we are in state A for
potential energy < energy_b.
dim : string
This select what dimension we should consider,
it should equal 'x', 'y' or 'z'.
periodic : boolean, optional
This determines if periodic boundary conditions should be
applied to the position.

"""
txt = 'Position of particle {} (dim: {})'.format(index, dim)
super().__init__(description=txt)
self.inter_a = inter_a
self.inter_b = inter_b
self.energy_a = energy_a
self.energy_b = energy_b
self.periodic = periodic
self.index = index
dims = {'x': 0, 'y': 1, 'z': 2}
try:
self.dim = dims[dim]
except KeyError:
msg = 'Unknown dimension {} requested'.format(dim)
logger.critical(msg)
raise

def calculate(self, system):
"""Calculate the order parameter.

Parameters
----------
system : object like :py:class:.System
This object is used for the actual calculation, typically
only system.particles.pos and/or system.particles.vel
will be used. In some cases system.forcefield can also be
used to include specific energies for the order parameter.

Returns
-------
out : float
The order parameter.

"""
particles = system.particles
pos = particles.pos[self.index]
lmb = pos[self.dim]
if lmb < self.inter_a:
if particles.vpot > self.energy_a:
lmb = self.inter_a
elif lmb > self.inter_b:
if particles.vpot > self.energy_b:
lmb = self.inter_b
if self.periodic:
return [system.box.pbc_coordinate_dim(lmb, self.dim)]
return [lmb]


This new order parameter can be included in a PyRETIS simulation by adding the following order parameter section (assuming that the order parameter is stored in a file order.py):

Orderparameter
--------------
class = OrderX
module = order.py
dim = x
index = 0
periodic = False
inter_a = -0.15
inter_b = 0.18
energy_a = -9.0
energy_b = -9.0


### Creating the projection-of-positions order parameter¶

For this order parameter, we will consider the distance from the particle to the stable A state and project the distance vector onto the line joining states A and B. In addition, as in the previously defined order parameters, we will include the energy in the definition of the states.

The Python code for this new order parameter is given below:

# -*- coding: utf-8 -*-
# Copyright (c) 2019, PyRETIS Development Team.
"""This file defines the order parameter used for the WCA example."""
import logging
import numpy as np
from pyretis.orderparameter import OrderParameter
logger = logging.getLogger(__name__)  # pylint: disable=invalid-name

class OrderXY(OrderParameter):
"""OrderXY(OrderParameter).

This class defines a 2D order parameter such that x + y = constant.

Attributes
----------
index : integer
This selects the particle to use for the order parameter.
inter_a : float
An interface such that we are in state A for postions < inter_a.
inter_b : float
An interface such that we are in state B for postions > inter_b.
energy_a : float
An energy such that we are in state A for potential energy < energy_a.
energy_b : float
An energy such that we are in state A for potential energy < energy_b.

"""

def __init__(self, index, inter_a, inter_b, energy_a, energy_b):
"""Set up the order parameter.

Parameters
----------
index : tuple of ints
This is the indices of the atom we will use the position of.
inter_a : float
An interface such that we are in state A for postions < inter_a.
inter_b : float
An interface such that we are in state B for postions > inter_b.
energy_a : float
An energy such that we are in state A for
potential energy < energy_a.
energy_b : float
An energy such that we are in state A for
potential energy < energy_b.

"""
super().__init__(description='2D->1D projection')
self.index = index
x_1 = 0.2
y_1 = 0.4
x_0 = -0.2
y_0 = -0.4
self.x_0 = x_0
self.y_0 = y_0
self.origin = np.array([x_0, y_0])
self.vec = np.array([x_1 - x_0, y_1 - y_0])
self.vec /= np.sqrt(np.dot(self.vec, self.vec))
self.inter_a = inter_a
self.inter_b = inter_b
self.energy_a = energy_a
self.energy_b = energy_b

def calculate(self, system):
"""Calculate the order parameter.

Here, the order parameter is just the distance between two
particles.

Parameters
----------
system : object like :py:class:.System
This object is used for the actual calculation, typically
only system.particles.pos and/or system.particles.vel
will be used. In some cases system.forcefield can also be
used to include specific energies for the order parameter.

Returns
-------
out : float
The order parameter.

"""
pos = system.particles.pos[self.index]
vec = np.array([pos[0] - self.x_0, pos[1] - self.y_0])
proj = np.dot(vec, self.vec)
if proj < self.inter_a:
if system.particles.vpot > self.energy_a:
proj = self.inter_a
elif proj > self.inter_b:
if system.particles.vpot > self.energy_b:
proj = self.inter_b
return [proj]


Note that this order parameter is less general than the previous ones, for instance, we make explicit use of the location of the two minimums. If you are feeling adventurous, please try to add these locations as parameters for the order parameter.

This new order parameter can be included in a PyRETIS simulation by adding the following order parameter section (assuming that the order parameter is stored in a file order.py):

Orderparameter
--------------
class = OrderXY
module = orderp.py
index = 0
inter_a = 0.15
inter_b = 0.75
energy_a = -9.0
energy_b = -9.0


## Running the RETIS simulation with PyRETIS¶

We have now defined the potential and the order parameters we are going to use and we are ready to run the RETIS simulations.

### Running the RETIS simulation using x as the order parameter¶

In order to run the RETIS simulation, the following steps are suggested:

1. Create a new folder, for instance 2D-hysteresis.

2. Place the order parameter script in this folder and name it order.py.

3. Download the initial configuration initial.xyz and place it in the same directory.

4. Create a new sub-directory, retis-x.

5. Add the following input script retis-x.rst to the retis-x folder.

6. Run the simulation using (in the retis-x folder):

pyretisrun -i retis-x.rst -p


Some examples from the analysis can be found below:

Fig. 27 Example of paths harvested for the RETIS simulation. All paths shown are accepted paths for the last path ensemble.

Fig. 28 The crossing probability after 1 000 000 RETIS cycles. The numerical value is . The initial flux is .

### Running the RETIS simulation using y as the order parameter¶

In order to run the RETIS simulation, the following steps are suggested:

1. Create a new folder, for instance 2D-hysteresis.

2. Place the order parameter script in this folder and name it order.py.

3. Download the initial configuration initial.xyz and place it in the same directory.

4. Create a new sub-directory, retis-y.

5. Add the following input script retis-y.rst to the retis-y folder.

6. Run the simulation using (in the retis-y folder):

pyretisrun -i retis-y.rst -p


Some examples from the analysis can be found below:

Fig. 29 Example of paths harvested for the RETIS simulation. All paths shown are accepted paths for the last path ensemble.

Fig. 30 The crossing probability after 1 000 000 RETIS cycles. The numerical value is . The initial flux is .

### Running the RETIS simulation using the projection-of-position as the order parameter¶

In order to run the RETIS simulation, the following steps are suggested:

1. Create a new folder, for instance 2D-hysteresis.

2. Create a new sub-directory, retis-xy and move into this directory.

3. Place the order parameter script in this folder and name it order.py.

4. Download the initial configuration initial2.xyz and place it in the same directory and rename it to initial.xyz.

5. Add the following input script retis-xy.rst to the retis-xy folder.

6. Run the simulation using (still in the retis-xy folder):

pyretisrun -i retis-xy.rst -p


Some examples from the analysis can be found below:

Fig. 31 Example of paths harvested for the RETIS simulation. All paths shown are accepted paths for the last path ensemble.

Fig. 32 The crossing probability after 1 000 000 RETIS cycles. The numerical value is . The initial flux is .

## Comparing results from the analysis¶

By using the three different order parameters, we obtain different estimates for the rate constant as summarized in the table below.

Table 50 Comparison of the rate constant for the 2D potential
Order parameter Rate constant / 1e-7
Lower bound Upper bound
x 7.3 8.3
y 7.7 8.6
xy 7.3 8.3