index | modules | gitlab | page source

# 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)

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)

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)

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)

  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) http://dx.doi.org/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(pairs=[(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`(a, direction)

  Constrain an atom index a to move on a given line only.

  The line is defined by its vector direction.

## The FixedPlane class

class `ase.constraints.``FixedPlane`(a, direction)

  Constrain an atom index a to move in a given plane only.

  The plane is defined by its normal vector direction.

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

## The FixedMode class

class `ase.constraints.``FixedMode`(mode)

  Constrain atoms to move along directions orthogonal to a given mode only.

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`

  Constraint class for fixing the center of mass.

  References

  https://pubs.acs.org/doi/abs/10.1021/jp9722824

Example of use:

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

## 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.

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)

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. 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. The epsilon value specifies the accuracy to which the constraints are fulfilled.

class `ase.constraints.``FixInternals`(bonds=None, angles=None, dihedrals=None, epsilon=1e-07)

  Constraint object for fixing multiple internal coordinates.

  Allows fixing bonds, angles, and dihedrals.

Note

The `FixInternals` class use radians for angles! Most other places in ASE degrees are used.

Example of use:

    >>> from math import pi
    >>> bond1 = [1.20, [1, 2]]
    >>> angle_indices1 = [2, 3, 4]
    >>> dihedral_indices1 = [2, 3, 4, 5]
    >>> angle1 = [atoms.get_angle(*angle_indices1) * pi / 180,
                  angle_indices1]
    >>> dihedral1 = [atoms.get_dihedral(*dihedral_indices1) * pi / 180,
    ...              dihedral_indices1]
    >>> c = FixInternals(bonds=[bond1], angles=[angle1],
    ...                  dihedrals=[dihedral1])
    >>> atoms.set_constraint(c)

This example defines a bond, an angle and a dihedral angle constraint to be fixed at the same time.

## 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 Filter class

Constraints can also be applied via filters, which acts as a wrapper around an atoms object. A typical use case will look like this:

     -------       --------       ----------
    |       |     |        |     |          |
    | Atoms |<----| Filter |<----| Dynamics |
    |       |     |        |     |          |
     -------       --------       ----------

and in Python this would be:

    >>> atoms = Atoms(...)
    >>> filter = Filter(atoms, ...)
    >>> dyn = Dynamics(filter, ...)

This class hides some of the atoms in an Atoms object.

class `ase.constraints.``Filter`(atoms, indices=None, mask=None)

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

Example of use:

    >>> from ase import Atoms, Filter
    >>> atoms=Atoms(positions=[[ 0    , 0    , 0],
    ...                        [ 0.773, 0.600, 0],
    ...                        [-0.773, 0.600, 0]],
    ...             symbols='OH2')
    >>> f1 = Filter(atoms, indices=[1, 2])
    >>> f2 = Filter(atoms, mask=[0, 1, 1])
    >>> f3 = Filter(atoms, mask=[a.Z == 1 for a in atoms])
    >>> f1.get_positions()
    [[ 0.773  0.6    0.   ]
     [-0.773  0.6    0.   ]]

In all three filters only the hydrogen atoms are made visible. When asking for the positions only the positions of the hydrogen atoms are returned.

## The UnitCellFilter class

The unit cell filter is for optimizing positions and unit cell simultaneously. Note that `ExpCellFilter` will probably perform better.

class `ase.constraints.``UnitCellFilter`(atoms, mask=None, cell_factor=None, hydrostatic_strain=False, constant_volume=False, scalar_pressure=0.0)

  Modify the supercell and the atom positions.

  Create a filter that returns the atomic forces and unit cell stresses together, so they can simultaneously be minimized.

  The first argument, atoms, is the atoms object. The optional second argument, mask, is a list of booleans, indicating which of the six independent components of the strain are relaxed.

  - True = relax to zero
  - False = fixed, ignore this component

  Degrees of freedom are the positions in the original undeformed cell, plus the deformation tensor (extra 3 “atoms”). This gives forces consistent with numerical derivatives of the potential energy with respect to the cell degreees of freedom.

  For full details see:
    E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras, Phys. Rev. B 59, 235 (1999)

  You can still use constraints on the atoms, e.g. FixAtoms, to control the relaxation of the atoms.

      >>> # this should be equivalent to the StrainFilter
      >>> atoms = Atoms(...)
      >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
      >>> ucf = UnitCellFilter(atoms)

  You should not attach this UnitCellFilter object to a trajectory. Instead, create a trajectory for the atoms, and attach it to an optimizer like this:

      >>> atoms = Atoms(...)
      >>> ucf = UnitCellFilter(atoms)
      >>> qn = QuasiNewton(ucf)
      >>> traj = Trajectory('TiO2.traj', 'w', atoms)
      >>> qn.attach(traj)
      >>> qn.run(fmax=0.05)

  Helpful conversion table:

  - 0.05 eV/A^3 = 8 GPA
  - 0.003 eV/A^3 = 0.48 GPa
  - 0.0006 eV/A^3 = 0.096 GPa
  - 0.0003 eV/A^3 = 0.048 GPa
  - 0.0001 eV/A^3 = 0.02 GPa

  Additional optional arguments:

  cell_factor: float (default float(len(atoms)))
    Factor by which deformation gradient is multiplied to put it on the same scale as the positions when assembling the combined position/cell vector. The stress contribution to the forces is scaled down by the same factor. This can be thought of as a very simple preconditioners. Default is number of atoms which gives approximately the correct scaling.
  hydrostatic_strain: bool (default False)
    Constrain the cell by only allowing hydrostatic deformation. The virial tensor is replaced by np.diag([np.trace(virial)]*3).
  constant_volume: bool (default False)
    Project out the diagonal elements of the virial tensor to allow relaxations at constant volume, e.g. for mapping out an energy-volume curve. Note: this only approximately conserves the volume and breaks energy/force consistency so can only be used with optimizers that do require do a line minimisation (e.g. FIRE).
  scalar_pressure: float (default 0.0)
    Applied pressure to use for enthalpy pV term. As above, this breaks energy/force consistency.

## The StrainFilter class

The strain filter is for optimizing the unit cell while keeping scaled positions fixed.

class `ase.constraints.``StrainFilter`(atoms, mask=None, include_ideal_gas=False)

  Modify the supercell while keeping the scaled positions fixed.

  Presents the strain of the supercell as the generalized positions, and the global stress tensor (times the volume) as the generalized force.

  This filter can be used to relax the unit cell until the stress is zero. If MDMin is used for this, the timestep (dt) to be used depends on the system size. 0.01/x where x is a typical dimension seems like a good choice.

  The stress and strain are presented as 6-vectors, the order of the components follow the standard engingeering practice: xx, yy, zz, yz, xz, xy.

  Create a filter applying a homogeneous strain to a list of atoms.

  The first argument, atoms, is the atoms object.

  The optional second argument, mask, is a list of six booleans, indicating which of the six independent components of the strain that are allowed to become non-zero. It defaults to [1,1,1,1,1,1].

## The ExpCellFilter class

The exponential cell filter is an improved `UnitCellFilter` which is parameter free.

class `ase.constraints.``ExpCellFilter`(atoms, mask=None, cell_factor=None, hydrostatic_strain=False, constant_volume=False, scalar_pressure=0.0)

  Modify the supercell and the atom positions.

  Create a filter that returns the atomic forces and unit cell stresses together, so they can simultaneously be minimized.

  The first argument, atoms, is the atoms object. The optional second argument, mask, is a list of booleans, indicating which of the six independent components of the strain are relaxed.

  - True = relax to zero
  - False = fixed, ignore this component

  Degrees of freedom are the positions in the original undeformed cell, plus the log of the deformation tensor (extra 3 “atoms”). This gives forces consistent with numerical derivatives of the potential energy with respect to the cell degrees of freedom.

  For full details see:
    E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras, Phys. Rev. B 59, 235 (1999)

  You can still use constraints on the atoms, e.g. FixAtoms, to control the relaxation of the atoms.

      >>> # this should be equivalent to the StrainFilter
      >>> atoms = Atoms(...)
      >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
      >>> ecf = ExpCellFilter(atoms)

  You should not attach this ExpCellFilter object to a trajectory. Instead, create a trajectory for the atoms, and attach it to an optimizer like this:

      >>> atoms = Atoms(...)
      >>> ecf = ExpCellFilter(atoms)
      >>> qn = QuasiNewton(ecf)
      >>> traj = Trajectory('TiO2.traj', 'w', atoms)
      >>> qn.attach(traj)
      >>> qn.run(fmax=0.05)

  Helpful conversion table:

  - 0.05 eV/A^3 = 8 GPA
  - 0.003 eV/A^3 = 0.48 GPa
  - 0.0006 eV/A^3 = 0.096 GPa
  - 0.0003 eV/A^3 = 0.048 GPa
  - 0.0001 eV/A^3 = 0.02 GPa

  Additional optional arguments:

  cell_factor: (DEPRECATED)
    Retained for backwards compatibility, but no longer used.
  hydrostatic_strain: bool (default False)
    Constrain the cell by only allowing hydrostatic deformation. The virial tensor is replaced by np.diag([np.trace(virial)]*3).
  constant_volume: bool (default False)
    Project out the diagonal elements of the virial tensor to allow relaxations at constant volume, e.g. for mapping out an energy-volume curve.
  scalar_pressure: float (default 0.0)
    Applied pressure to use for enthalpy pV term. As above, this breaks energy/force consistency.

  Implementation details:

  The implementation is based on that of Christoph Ortner in JuLIP.jl: https://github.com/libAtoms/JuLIP.jl/blob/expcell/src/Constraints.jl#L244

  We decompose the deformation gradient as

  > F = exp(U) F0 x = F * F0^{-1} z = exp(U) z

  If we write the energy as a function of U we can transform the stress associated with a perturbation V into a derivative using a linear map V -> L(U, V).

  phi( exp(U+tV) (z+tv) ) ~ phi’(x) . (exp(U) v) + phi’(x) . ( L(U, V) exp(-U) exp(U) z )

        >>> \nabla E(U) : V  =  [S exp(-U)'] : L(U,V)
                             =  L'(U, S exp(-U)') : V
                             =  L(U', S exp(-U)') : V
                             =  L(U, S exp(-U)) : V     (provided U = U')

  where the : operator represents double contraction, i.e. A:B = trace(A’B), and

  > F = deformation tensor - 3x3 matrix F0 = reference deformation tensor - 3x3 matrix, np.eye(3) here U = cell degrees of freedom used here - 3x3 matrix V = perturbation to cell DoFs - 3x3 matrix v = perturbation to position DoFs x = atomic positions in deformed cell z = atomic positions in original cell phi = potential energy S = stress tensor [3x3 matrix] L(U, V) = directional derivative of exp at U in direction V, i.e d/dt exp(U + t V)|_{t=0} = L(U, V)

  This means we can write

  > d/dt E(U + t V)|_{t=0} = L(U, S exp (-U)) : V

  and therefore the contribution to the gradient of the energy is

  > nabla E(U) / nabla U_ij = [L(U, S exp(-U))]_ij

Documentation v: 3.19.0

Versions
  3.20.0b1 (development)
  3.19.0 (latest stable)