TinyFF

This is a minimalistic force-field engine written in pure Python, using vectorized NumPy code. It has minimal dependencies (NumPy, SciPy, attrs and npy-append-array), so all the force-field specific code is self-contained.

This little library is geared towards teaching and favors simplicity and conciseness over fancy features and top-notch performance. TinyFF implements the nitty-gritty of linear-scaling pairwise potentials, so students can build their own molecular dynamics implementation, skipping the technicalities of implementing the correct potential energy, pressure and forces acting on atoms.

TinyFF is written by Toon Verstraelen for students of the Computational Physics course (C004504) in the Physics and Astronomy program at Ghent University. TinyFF is distributed under the conditions of the GPL-v3 license.

Installation

TinyFF is available on PyPI. In a properly configured Python virtual environment, you can install TinyFF with:

pip install tinyff

Migration from older versions

Code written for TinyFF 1 can be easily updated to work with TinyFF 2. All required changes are documented in the Migration Guide.

A complete list of changes can be found in CHANGELOG.md.

Features

TinyFF is a Python package with the following modules:

• tinyff.analysis: helper functions to compute a radial distribution function and an autocorrelation function.
• tinyff.atomsmithy: functions for creating initial atomic positions.
• tinyff.forcefield: implements a general force field interface: energy, atomic forces and pressure
• tinyff.pairwise: pairwise potentials to be used in force fields.
• tinyff.neighborlist: used by the forcefield module to compute pairwise interactions with a real-space cut-off.
• tinyff.trajectory: tools for writing molecular dynamics trajectories to file(s).
• tinyff.utils: utility functions used in tinyff.

Most relevant functions and classes can be imported directly from the top-level tinyff package. None of these modules implement any molecular dynamics integrators. This is the part that you are expected to write.

Basic usage

Computing energy, forces and pressure

The evaluation of the force field energy and its derivatives requires the following:

import numpy as np
from tinyff import CutOffWrapper, LennardJones, ForceField, NBuildSimple

# Define a pairwise potential, with energy and force shift
rcut = 5.0
lj = CutOffWrapper(LennardJones(2.5, 2.0), rcut)

# Define a force field
ff = ForceField([lj], nbuild=NBuildSimple(rcut))

# You need atomic positions and the length of a periodic cell edge.
# The following line defines just two atomic positions.
atpos = np.array([[0.0, 0.0, 1.0], [1.0, 2.0, 0.0]])

# Note that the cell must be large enough to contain the cutoff sphere.
cell_length = 20.0

# Compute a selection of results with ff.compute.
#   The ff.compute method has an `nderiv` arguments to compute only some results:
#   - `nderiv=0` (default): compute only the energy
#   - `nderiv=1`: compute the energy, forces and force contribution to the pressure.
#   Requested results are put in a list, even in only one result is requested.
potential_energy, forces, press = ff.compute(atpos, cell_length, nderiv=1)

This basic recipe can be extended by passing additional options into the ForceField constructor:

• Linear-scaling neighbor lists with the cell lists method:

  from tinyff import NBuildCellLists
  
  # Construct your force field object as follows:
  ff = ForceField([lj], nbuild=NBuildCellLists(rcut, nbin_approx=natom / 100))

  For about 430 atoms, this becomes more efficient than the simple neighborlist build.

• Verlet lists (cut-off radius + buffer):

  rmax = 6.0  # > rcut, so buffer of 1.0
  ff = ForceField([lj], nbuild=NBuildSimple(rmax, nlist_reuse=16))

Forging initial positions

The atomsmithy defines functions to generate a cubic box with standard lattices or randomized atomic positions:

from tinyff import (
    build_bcc_lattice,
    build_cubic_lattice,
    build_fcc_lattice,
    build_random_cell,
)

# Atoms on a regular lattice. args:
# - primitive cell edge length
# - number of repetitions in X, Y and Z directions.
atpos = build_cubic_lattice(2.5, 2)
atpos = build_bcc_lattice(2.5, 3)
atpos = build_fcc_lattice(2.5, 4)

# Randomize positions. args:
# - cell edge length
# - number of atoms
atpos = build_random_cell(10.0, 32, 2.5)

Writing trajectories to disk

For visualization with nglview, TinyFF provides a PDBWriter, to be used as follows:

from tinyff import PDBWriter

# Initialization of the writer: specify a file and a conversion factor to angstrom.
# If the PDB file exists, it is overwritten!
# This example shows the conversion factor when your program works in nanometer.
# Through `atnums` you can specify the chemical elements, here 50 argon atoms (Z=18).
pdb_writer = PDBWriter("trajectory.pdb", to_angstrom=10.0, atnums=[18] * 50)

# Somewhere in your code, typically inside some loop.
# cell_length(s) can be a float or an array of 3 floats.
pdb_writer.dump(atpos, cell_length)

# If you are using Jupyter, you can visualize the trajectory with nglview as follows:
import mdtraj
import nglview
traj = mdtraj.load(f"trajectory.pdb")
view = nglview.show_mdtraj(traj)
view.clear()
view.add_hyperball()
view.add_unitcell()
view

For numerical post-processing, TinyFF provides a more flexible trajectory writer, which writes NPY files. It is implemented with the npy-append-array library, which makes it possible to extend an NPY file without having to rewrite from scratch. (That would not be possible with NumPy alone.) The NPYWriter can be used as follows:

from tinyff import NPYWriter

# Initialization, will create (and possibly clean up an existing) a `traj` directory.
npy_writer = NPYWriter("traj")

# Somewhere in your production code, normally inside some loop.
# You can specify any array or float you like,
# as long as the shape and type is the same upon every call.
# This will result in files `traj/atpos.npy`, `traj/pressure.npy`, etc.
# These files will contain arrays with data passed into all `dump` calls.
npy_writer.dump(atpos=atpos, pressure=pressure, temperature=temperature, ...)

# In your post-processing code
pressure = np.load("traj/pressure.npz")
print(np.mean(pressure))

Irrespective of the file format, it is recommended to write trajectory data only every so many steps. The constructors of both trajectory writers have an optional stride argument to control the frequency of the output. For example:

npy_writer = NPYWriter("traj", stride=100)

Monte Carlo: partial updates

TinyFF supports partial updates in which the change in energy is computed due to the displacement of a single atom with the methods ForceField.try_move and ForceField.accept_move. These methods assume that a neighborlist was previously built. They will only recompute results for known pairs of atoms, without adding or removing pairs. This implies that your Monte Carlo implementation needs to take care of the following:

• Build an initial neighborlist before starting the Monte Carlo loop.
• Construct the neighborlist up to a radius rmax, which must extend beyond the cutoff radius.
• Rebuild the neighborlist at regular intervals, e.g. by calling ForceField.compute, which will also compute the energy by default. At this stage it is useful to check whether the new energy is still consistent with the sum of all energy changes due to accepted Monte Carlo moves.

The following example just shows how to call the functions. You can use these function calls in your Monte Carlo loop:

import numpy as np
from tinyff import build_fcc_lattice, CutOffWrapper, LennardJones, ForceField, NBuildSimple

# System configuration, a simple (inflated) FCC lattice of Argon atoms.
atpos = build_fcc_lattice(2.5, 4)
cell_lengths = np.array([10.0, 10.0, 10.0])
rmax = 3.0
rcut = 2.5

# Define the force field and compute the initial energy.
lj = CutOffWrapper(LennardJones(2.5, 2.0), rcut)
ff = ForceField([lj], nbuild=NBuildSimple(rmax))
energy0, = ff.compute(atpos, cell_length)

# Try and accept a move of atom 3.
iatom = 3
delta = np.array([0.1, 0.2, -0.1])
energy_change, move = ff.try_move(iatom, delta, cell_lengths)
ff.accept_move(move)
atpos[iatom] += delta

# Verify the change in energy.
energy1, = ff.compute(atpos, cell_length)
assert abs(energy_change - (energy1 - energy0)) < 1e-10