Constraints

When performing minimizations or dynamics one may wish to keep some degrees of freedom in the system fixed. One way of doing this is by attaching constraint object(s) directly to the atoms object.

Important: setting constraints will freeze the corresponding atom positions. Changing such atom positions can be achieved:

  • by directly setting the positions attribute (see example of setting Special attributes),

  • alternatively, by removing the constraints first:

    del atoms.constraints

    or:

    atoms.set_constraint()

    and using the set_positions() method.

The FixAtoms class

This class is used for fixing some of the atoms.

class ase.constraints.FixAtoms(indices=None, mask=None)[source]

You must supply either the indices of the atoms that should be fixed or a mask. The mask is a list of booleans, one for each atom, being true if the atoms should be kept fixed.

For example, to fix the positions of all the Cu atoms in a simulation with the indices keyword:

>>> from ase.constraints import FixAtoms
>>> c = FixAtoms(indices=[atom.index for atom in atoms if atom.symbol == 'Cu'])
>>> atoms.set_constraint(c)

or with the mask keyword:

>>> c = FixAtoms(mask=[atom.symbol == 'Cu' for atom in atoms])
>>> atoms.set_constraint(c)

The FixBondLength class

This class is used to fix the distance between two atoms specified by their indices (a1 and a2)

class ase.constraints.FixBondLength(a1, a2)[source]

Example of use:

>>> c = FixBondLength(0, 1)
>>> atoms.set_constraint(c)

In this example the distance between the atoms with indices 0 and 1 will be fixed in all following dynamics and/or minimizations performed on the atoms object.

This constraint is useful for finding minimum energy barriers for reactions where the path can be described well by a single bond length (see the Dissociation of a molecule using the NEB method tutorial).

Important: If fixing multiple bond lengths, use the FixBondLengths class below, particularly if the same atom is fixed to multiple partners.

The FixBondLengths class

RATTLE-type holonomic constraints. More than one bond length can be fixed by using this class. Especially for cases in which more than one bond length constraint is applied on the same atom. It is done by specifying the indices of the two atoms forming the bond in pairs.

class ase.constraints.FixBondLengths(pairs)[source]

Example of use:

>>> c = FixBondLengths([[0, 1], [0, 2]])
  >>> atoms.set_constraint(c)

  Here the distances between atoms with indices 0 and 1 and atoms with
  indices 0 and 2 will be fixed. The constraint is for the same purpose
  as the FixBondLength class.

The FixLinearTriatomic class

This class is used to keep the geometry of linear triatomic molecules rigid in geometry optimizations or molecular dynamics runs. Rigidness of linear triatomic molecules is impossible to attain by constraining all interatomic distances using FixBondLength, as this won’t remove an adequate number of degrees of freedom. To overcome this, FixLinearTriatomic fixes the distance between the outer atoms with RATTLE and applies a linear vectorial constraint to the central atom using the RATTLE-constrained positions of the outer atoms (read more about the method here: G. Ciccotti, M. Ferrario, J.-P. Ryckaert, Molecular Physics 47, 1253 (1982)).

When setting these constraints one has to specify a list of triples of atomic indices, each triple representing a specific triatomic molecule.

class ase.constraints.FixLinearTriatomic(triples)[source]

Holonomic constraints for rigid linear triatomic molecules.

Apply RATTLE-type bond constraints between outer atoms n and m and linear vectorial constraints to the position of central atoms o to fix the geometry of linear triatomic molecules of the type:

n–o–m

Parameters:

triples: list

Indices of the atoms forming the linear molecules to constrain as triples. Sequence should be (n, o, m) or (m, o, n).

When using these constraints in molecular dynamics or structure optimizations, atomic forces need to be redistributed within a triple. The function redistribute_forces_optimization implements the redistribution of forces for structure optimization, while the function redistribute_forces_md implements the redistribution for molecular dynamics.

References:

Ciccotti et al. Molecular Physics 47 (1982) doi: 10.1080/00268978200100942

The example below shows how to fix the geometry of two carbon dioxide molecules:

>>> from ase.build import molecule
>>> from ase.constraints import FixLinearTriatomic
>>> atoms = molecule('CO2')
>>> dimer = atoms + atoms.copy()
>>> c = FixLinearTriatomic(triples=[(1, 0, 2), (4, 3, 5)])
>>> dimer.set_constraint(c)

Note

When specifying a triple of indices, the second element must correspond to the index of the central atom.

The FixedLine class

class ase.constraints.FixedLine(indices, direction)[source]

Constrain an atom index or a list of atom indices to move on a line only.

The line is defined by its vector direction

Constrain chosen atoms.

Parameters:

  • indices (int or list of int) – Index or indices for atoms that should be constrained

  • direction (list of 3 int) – Direction of the vector defining the line

Examples

Fix all Copper atoms to only move in the x-direction:

>>> from ase.constraints import FixedLine
>>> c = FixedLine(
...     indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],
...     direction=[1, 0, 0],
... )
>>> atoms.set_constraint(c)

or constrain a single atom with the index 0 to move in the z-direction:

>>> c = FixedLine(indices=0, direction=[0, 0, 1])
>>> atoms.set_constraint(c)

The FixedPlane class

class ase.constraints.FixedPlane(indices, direction)[source]

Constraint object for fixing chosen atoms to only move in a plane.

The plane is defined by its normal vector direction

Constrain chosen atoms.

Parameters:

  • indices (int or list of int) – Index or indices for atoms that should be constrained

  • direction (list of 3 int) – Direction of the normal vector

Examples

Fix all Copper atoms to only move in the yz-plane:

>>> from ase.build import bulk
>>> from ase.constraints import FixedPlane

>>> atoms = bulk('Cu', 'fcc', a=3.6)
>>> c = FixedPlane(
...     indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],
...     direction=[1, 0, 0],
... )
>>> atoms.set_constraint(c)

or constrain a single atom with the index 0 to move in the xy-plane:

>>> c = FixedPlane(indices=0, direction=[0, 0, 1])
>>> atoms.set_constraint(c)

Example of use: Surface diffusion energy barriers using ASE constraints.

The FixedMode class

class ase.constraints.FixedMode(mode)[source]

Constrain atoms to move along directions orthogonal to a given mode only. Initialize with a mode, such as one produced by ase.vibrations.Vibrations.get_mode().

A mode is a list of vectors specifying a direction for each atom. It often comes from ase.vibrations.Vibrations.get_mode().

The FixCom class

class ase.constraints.FixCom[source]

Constraint class for fixing the center of mass.

Example of use:

>>> from ase.constraints import FixCom
>>> c = FixCom()
>>> atoms.set_constraint(c)

The FixSubsetCom class

class ase.constraints.FixSubsetCom(indices)[source]

Constraint class for fixing the center of mass of a subset of atoms.

The Hookean class

This class of constraints, based on Hooke’s Law, is generally used to conserve molecular identity in optimization schemes and can be used in three different ways. In the first, it applies a Hookean restorative force between two atoms if the distance between them exceeds a threshold. This is useful to maintain the identity of molecules in quenched molecular dynamics, without changing the degrees of freedom or violating conservation of energy. When the distance between the two atoms is less than the threshold length, this constraint is completely inactive.

class ase.constraints.Hookean(a1, a2, k, rt=None)[source]

Applies a Hookean restorative force between a pair of atoms, an atom and a point, or an atom and a plane.

Forces two atoms to stay close together by applying no force if they are below a threshold length, rt, and applying a Hookean restorative force when the distance between them exceeds rt. Can also be used to tether an atom to a fixed point in space or to a distance above a plane.

a1int

Index of atom 1

a2one of three options

  1. index of atom 2

  2. a fixed point in cartesian space to which to tether a1

  3. a plane given as (A, B, C, D) in A x + B y + C z + D = 0.

kfloat

Hooke’s law (spring) constant to apply when distance exceeds threshold_length. Units of eV A^-2.

rtfloat

The threshold length below which there is no force. The length is 1) between two atoms, 2) between atom and point. This argument is not supplied in case 3. Units of A.

If a plane is specified, the Hooke’s law force is applied if the atom is on the normal side of the plane. For instance, the plane with (A, B, C, D) = (0, 0, 1, -7) defines a plane in the xy plane with a z intercept of +7 and a normal vector pointing in the +z direction. If the atom has z > 7, then a downward force would be applied of k * (atom.z - 7). The same plane with the normal vector pointing in the -z direction would be given by (A, B, C, D) = (0, 0, -1, 7).

References

Andrew A. Peterson, Topics in Catalysis volume 57, pages40–53 (2014) https://link.springer.com/article/10.1007%2Fs11244-013-0161-8

The below example tethers atoms at indices 3 and 4 together:

>>> c = Hookean(a1=3, a2=4, rt=1.79, k=5.)
>>> atoms.set_constraint(c)

Alternatively, this constraint can tether a single atom to a point in space, for example to prevent the top layer of a slab from subliming during a high-temperature MD simulation. An example of tethering atom at index 3 to its original position:

>>> from ase.constraints import Hookean
>>> c = Hookean(a1=3, a2=atoms[3].position, rt=0.94, k=2.)
>>> atoms.set_constraint(c)

Reasonable values of the threshold (rt) and spring constant (k) for some common bonds are below.

Bond

rt (Angstroms)

k (eV Angstrom^-2)

O-H

1.40

5

C-O

1.79

5

C-H

1.59

7

C=O

1.58

10

Pt sublimation

0.94

2

Cu sublimation

0.97

2

A third way this constraint can be applied is to apply a restorative force if an atom crosses a plane in space. For example:

>>> c = Hookean(a1=3, a2=(0, 0, 1, -7), k=10.)
>>> atoms.set_constraint(c)

This will apply a restorative force on atom 3 in the downward direction of magnitude k * (atom.z - 7) if the atom’s vertical position exceeds 7 Angstroms. In other words, if the atom crosses to the (positive) normal side of the plane, the force is applied and directed towards the plane. (The same plane with the normal direction pointing in the -z direction would be given by (0, 0, -1, 7).)

For an example of use, see the Constrained minima hopping (global optimization) tutorial.

Note

In previous versions of ASE, this was known as the BondSpring constraint.

The ExternalForce class

This class can be used to simulate a constant external force (e.g. the force of atomic force microscope). One can set the absolute value of the force f_ext (in eV/Ang) and two atom indices a1 and a2 to define on which atoms the force should act. If the sign of the force is positive, the two atoms will be pulled apart. The external forces which acts on both atoms are parallel to the connecting line of the two atoms.

class ase.constraints.ExternalForce(a1, a2, f_ext)[source]

Example of use:

>>> form ase.constraints import ExternalForce
>>> c = ExternalForce(0, 1, 0.5)
>>> atoms.set_constraint(c)

One can combine this constraint with FixBondLength but one has to consider the correct ordering when setting both constraints. ExternalForce must come first in the list as shown in the following example.

>>> from ase.constraints import ExternalForce, FixBondLength
>>> c1 = ExternalForce(0, 1, 0.5)
>>> c2 = FixBondLength(1, 2)
>>> atoms.set_constraint([c1, c2])

The FixInternals class

This class allows to fix an arbitrary number of bond lengths, angles and dihedral angles as well as linear combinations of bond lengths (‘bondcombos’). A fixed linear combination of bond lengths fulfils \(\sum_i \mathrm{coef}_i \times \mathrm{bond_length}_i = \mathrm{constant}\). The defined constraints are satisfied self consistently. To define the constraints one needs to specify the atoms object on which the constraint works (needed for atomic masses), a list of bond, angle and dihedral constraints. Those constraint definitions are always list objects containing the value to be set and a list of atomic indices. For the linear combination of bond lengths the list of atomic indices is a list of bond definitions with coeficients ([[a1, a2, coef],[a3, a4, coef],]). The usage of mic is supported by providing the keyword argument \(mic=True\). Using mic slows the algorithm and is probably not necessary in most cases. The epsilon value specifies the accuracy to which the constraints are fulfilled. Please specify angles and dihedrals in degrees using the keywords angles_deg and dihedrals_deg.

class ase.constraints.FixInternals(bonds=None, angles=None, dihedrals=None, angles_deg=None, dihedrals_deg=None, bondcombos=None, mic=False, epsilon=1e-07)[source]

Constraint object for fixing multiple internal coordinates.

Allows fixing bonds, angles, dihedrals as well as linear combinations of bonds (bondcombos).

Please provide angular units in degrees using \(angles_deg\) and \(dihedrals_deg\). Fixing planar angles is not supported at the moment.

A constrained internal coordinate is defined as a nested list: ‘[value, [atom indices]]’. The constraint is initialized with a list of constrained internal coordinates, i.e. ‘[[value, [atom indices]], …]’. If ‘value’ is None, the current value of the coordinate is constrained.

Parameters:

  • bonds (nested python list, optional) – List with targetvalue and atom indices defining the fixed bonds, i.e. [[targetvalue, [index0, index1]], …]

  • angles_deg (nested python list, optional) – List with targetvalue and atom indices defining the fixedangles, i.e. [[targetvalue, [index0, index1, index3]], …]

  • dihedrals_deg (nested python list, optional) – List with targetvalue and atom indices defining the fixed dihedrals, i.e. [[targetvalue, [index0, index1, index3]], …]

  • bondcombos (nested python list, optional) – List with targetvalue, atom indices and linear coefficient defining the fixed linear combination of bonds, i.e. [[targetvalue, [[index0, index1, coefficient_for_bond], [index1, index2, coefficient_for_bond]]], …]

  • mic (bool, optional, default: False) – Minimum image convention.

  • epsilon (float, optional, default: 1e-7) – Convergence criterion.

static get_bondcombo(atoms, indices, mic=False)[source]

Convenience function to return the value of the bondcombo coordinate (linear combination of bond lengths) for the given Atoms object ‘atoms’. Example: Get the current value of the linear combination of two bond lengths defined as \(bondcombo = [[0, 1, 1.0], [2, 3, -1.0]]\).

Example of use:

>>> bond1 = [1.20, [1, 2]]
>>> angle_indices1 = [2, 3, 4]
>>> dihedral_indices1 = [2, 3, 4, 5]
>>> bondcombo_indices1 = [[6, 7, 1.0], [8, 9, -1.0]]
>>> angle1 = [atoms.get_angle(*angle_indices1), angle_indices1]
>>> dihedral1 = [atoms.get_dihedral(*dihedral_indices1), dihedral_indices1]
>>> bondcombo1 = [0.0, bondcombo_indices1]
>>> c = FixInternals(bonds=[bond1], angles_deg=[angle1],
...                  dihedrals_deg=[dihedral1], bondcombos=[bondcombo1])
>>> atoms.set_constraint(c)

This example defines a bond, an angle and a dihedral angle constraint to be fixed at the same time at which also the linear combination of bond lengths \(1.0 * \text{bond}_{6-7} -1.0 * \text{bond}_{8-9}\) is fixed to the value of 0.0 Ångstrom.

Combining constraints

It is possible to supply several constraints on an atoms object. For example one may wish to keep the distance between two nitrogen atoms fixed while relaxing it on a fixed ruthenium surface:

>>> pos = [[0.00000, 0.00000,  9.17625],
...        [0.00000, 0.00000, 10.27625],
...        [1.37715, 0.79510,  5.00000],
...        [0.00000, 3.18039,  5.00000],
...        [0.00000, 0.00000,  7.17625],
...        [1.37715, 2.38529,  7.17625]]
>>> unitcell = [5.5086, 4.7706, 15.27625]

>>> atoms = Atoms(positions=pos,
...               symbols='N2Ru4',
...               cell=unitcell,
...               pbc=[True,True,False])

>>> fa = FixAtoms(mask=[a.symbol == 'Ru' for a in atoms])
>>> fb = FixBondLength(0, 1)
>>> atoms.set_constraint([fa, fb])

When applying more than one constraint they are passed as a list in the set_constraint() method, and they will be applied one after the other.

Important: If wanting to fix the length of more than one bond in the simulation, do not supply a list of FixBondLength instances; instead, use a single instance of FixBondLengths.

Making your own constraint class

A constraint class must have these two methods:

ase.constraints.adjust_positions(oldpositions, newpositions)

Adjust the newpositions array inplace.

ase.constraints.adjust_forces(positions, forces)

Adjust the forces array inplace.

A simple example:

import numpy as np
class MyConstraint:
    """Constrain an atom to move along a given direction only."""
    def __init__(self, a, direction):
        self.a = a
        self.dir = direction / sqrt(np.dot(direction, direction))

    def adjust_positions(self, atoms, newpositions):
        step = newpositions[self.a] - atoms.positions[self.a]
        step = np.dot(step, self.dir)
        newpositions[self.a] = atoms.positions[self.a] + step * self.dir

    def adjust_forces(self, atoms, forces):
        forces[self.a] = self.dir * np.dot(forces[self.a], self.dir)

A constraint can optionally have two additional methods, which will be ignored if missing:

ase.constraints.adjust_momenta(atoms, momenta)

Adjust the momenta array inplace.

ase.constraints.adjust_potential_energy(atoms, energy)

Provide the difference in the potential energy due to the constraint. (Note that inplace adjustment is not possible for energy, which is a float.)

The FixSymmetry class

class ase.constraints.FixSymmetry(atoms, symprec=0.01, adjust_positions=True, adjust_cell=True, verbose=False)[source]

Constraint to preserve spacegroup symmetry during optimisation.

Requires spglib package to be available.

The following are some utility functions to prepare symmetrized configurations and to check symmetry.

ase.spacegroup.symmetrize.refine_symmetry(atoms, symprec=0.01, verbose=False)[source]

Refine symmetry of an Atoms object

Parameters:

  • object (atoms - input Atoms)

  • precicion (symprec - symmetry)

  • True (verbose - if)

  • after (print out symmetry information before and)

Return type:

spglib dataset

ase.spacegroup.symmetrize.check_symmetry(atoms, symprec=1e-06, verbose=False)[source]

Check symmetry of \(atoms\) with precision \(symprec\) using \(spglib\)

Prints a summary and returns result of \(spglib.get_symmetry_dataset()\)

Here is an example of using these tools to demonstrate the difference between minimising a perturbed bcc Al cell with and without symmetry-preservation. Since bcc is unstable with respect to fcc with a Lennard Jones model, the unsymmetrised case relaxes to fcc, while the constraint keeps the original symmetry.

import numpy as np

from ase.build import bulk
from ase.calculators.lj import LennardJones
from ase.constraints import FixSymmetry, UnitCellFilter
from ase.optimize import BFGS
from ase.spacegroup.symmetrize import check_symmetry

# We setup a bcc Al cell - bcc is unstable with LJ potential
# so without constraint this would relax back to an fcc structure
atoms_prim = bulk('Al', 'bcc', a=2 / np.sqrt(3), cubic=True)

# Now we setup a 2x2x2 supercell, and break the symmetry slightly
atoms_init = atoms_prim * [2, 2, 2]
atoms_init.positions[0, 0] += 1.0e-7  # break symmetry by 1e-7

# We use an LJ calculator, and allow the cell and atomic positions to relax
atoms_unsym = atoms_init.copy()
atoms_unsym.calc = LennardJones()
ucf_unsym = UnitCellFilter(atoms_unsym)

dyn = BFGS(ucf_unsym)
print("Initial Energy", atoms_unsym.get_potential_energy())
dyn.run(fmax=0.001)
print("Final Energy", atoms_unsym.get_potential_energy())

# Now we repeat the optimisation with the symmetrization constraint in place
atoms_sym = atoms_init.copy()
atoms_sym.calc = LennardJones()
atoms_sym.set_constraint(FixSymmetry(atoms_sym))
ucf_sym = UnitCellFilter(atoms_sym)

dyn = BFGS(ucf_sym)
print("Initial Energy", atoms_sym.get_potential_energy())
dyn.run(fmax=0.001)
print("Final Energy", atoms_sym.get_potential_energy())

print("position difference", np.linalg.norm(atoms_unsym.get_positions() -
                                            atoms_sym.get_positions()))

# We print out the initial symmetry groups at two different precision levels
print("initial symmetry at precision 1e-6")
check_symmetry(atoms_init, 1.0e-6, verbose=True)
print("initial symmetry at precision 1e-8")
check_symmetry(atoms_init, 1.0e-8, verbose=True)

# Printing the final symmetries shows that
# the "unsym" case relaxes to a lower energy fcc structure
# with a change in spacegroup, while the "sym" case stays as bcc
print("unsym symmetry after relaxation")
d_unsym = check_symmetry(atoms_unsym, verbose=True)
print("sym symmetry after relaxation")
d_sym = check_symmetry(atoms_sym, verbose=True)