Getting started

pele is a package for global optimization and landscape analysis. Follow Installation to install pele before running these examples.

The general sequence of steps when using these methods is

1. Find the global minimum and build up a database of other important minima

   • We use Basinhopping to do the global optimization.

2. Find connections between the minima.

   • This means finding the saddle point separating two minima such that the minimum energy path between the two minima crosses through the saddle point.

   • We use DoubleEndedConnect which is a combination of the Nudged Elastic Band (NEBDriver) and Hybrid Eigenvector Following (FindTransitionState)

3. Landscape analysis.

   • Visualize the landscape with a disconnectivity graph.

   • Compute rates between minima or groups of minima

Note

This tutorial mirrors the examples in examples/getting_started/.

First steps

pele is first and foremost a library. There is no single executable that you can call. Instead, you use python as a scripting language to explicitly state what you want to do, importing the tools from pele to do it.

As a very simple example, let’s generate a random Lennard-Jones configuration and minimize it.

import numpy as np
from pele.potentials import LJ
from pele.optimize import lbfgs_py

natoms = 5
x = np.random.uniform(-2, 2, natoms*3)
pot = LJ()
result = lbfgs_py(x, pot)
print(result)

In the above we used three imports. We used the numpy library to construct a random one dimensional array. We used the Lennard-Jones potential LJ, and we used the minimization routine lbfgs_py() which is just a wrapper for the class LBFGS.

The return value is an optimization result, which is just a container (Result) that stores the final energy, final coordinates, the number of function calls, etc.

If we want to then save the minimized coordinates in an xyz file we can use the function write_xyz()

from pele.utils.xyz import write_xyz
with open("out.xyz", "w") as fout:
    title = "energy = " + str(result.energy)
    write_xyz(fout, result.coords, title=title)

Global optimization

We primarily use the global optimization method BasinHopping because it has been shown to be very effective on the very high dimensional, smooth landscapes that we work with.

BasinHopping is iterative with each cycle composed of the following features

1. random perturbation of the coordinates

2. local minimization

3. accept or reject the new coordinates based on the minimized function value

Let’s run BasinHopping on a Lennard-Jones cluster (LJ).

Note

This example involves quite a number of steps. Normally these are done automatically in functions defined in the system class (systems). With the system class the following can be done in just a few steps:

from pele.systems import LJCluster
natoms = 17
system = LJCluster(natoms)
database = system.create_database('lj17.sqlite')
bh = system.get_basinhopping(database)
bh.run(10)

The following example shows how to do this manually

We start by defining the potential and choosing a random set of starting coordinates

import numpy as np
from pele.potentials import LJ
natoms = 17
potential = LJ()
x0 = np.random.uniform(-1, 1, 3*natoms)

We also need to set up the take-step class (takestep). We will use simple random displacements of the coordinates (RandomDisplacement) and wrap it with the class AdaptiveStepsizeTemperature that adjusts both the stepsize and temperature to achieve optimum results.

from pele.takestep import RandomDisplacement, AdaptiveStepsizeTemperature
displace = RandomDisplacement()
adaptive_displacement = AdaptiveStepsizeTemperature(displace)

Next we set up a SQLite database (storage) which we use to store the minima that we find during the basinhopping run.

from pele.storage import Database
database = Database("lj17.sqlite")

There is now an sqlite database named “lj17.sqlite” in the folder where this script was run. We use sqlalchemy to communicate with the sqlite database from python.

Now construct and run BasinHopping. Calling database.minimum_adder() returns the callback that stores each minimum found during the run:

from pele.basinhopping import BasinHopping
bh = BasinHopping(
    x0, potential, adaptive_displacement, storage=database.minimum_adder()
)
bh.run(10)

The minima in the database can be accessed simply by:

for m in database.minima():
    print(m.energy)

You can access other attributes of the Minimum as minimum.coords in the same way.

Landscape exploration

See the tutorial on finding pathways for detailed information about how to connect the minima using transition states.