Skip to main content
`Ctrl`+`K`

#

# Source code for ase.constraints.hookean

    import numpy as np

    from ase.constraints.constraint import FixConstraint, slice2enlist
    from ase.geometry.geometry import find_mic

    class Hookean(FixConstraint):
        """Applies a Hookean restorative force between a pair of atoms, an atom
        and a point, or an atom and a plane."""

        def __init__(self, a1, a2, k, rt=None):
            """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.

            a1 : int
               Index of atom 1
            a2 : one 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.
            k : float
               Hooke's law (spring) constant to apply when distance
               exceeds threshold_length. Units of eV A^-2.
            rt : float
               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
            """

            if isinstance(a2, int):
                self._type = 'two atoms'
                self.indices = [a1, a2]
            elif len(a2) == 3:
                self._type = 'point'
                self.index = a1
                self.origin = np.array(a2)
            elif len(a2) == 4:
                self._type = 'plane'
                self.index = a1
                self.plane = a2
            else:
                raise RuntimeError('Unknown type for a2')
            self.threshold = rt
            self.spring = k

        def get_removed_dof(self, atoms):
            return 0

        def todict(self):
            dct = {'name': 'Hookean'}
            dct['kwargs'] = {'rt': self.threshold, 'k': self.spring}
            if self._type == 'two atoms':
                dct['kwargs']['a1'] = self.indices[0]
                dct['kwargs']['a2'] = self.indices[1]
            elif self._type == 'point':
                dct['kwargs']['a1'] = self.index
                dct['kwargs']['a2'] = self.origin
            elif self._type == 'plane':
                dct['kwargs']['a1'] = self.index
                dct['kwargs']['a2'] = self.plane
            else:
                raise NotImplementedError(f'Bad type: {self._type}')
            return dct

        def adjust_positions(self, atoms, newpositions):
            pass

        def adjust_momenta(self, atoms, momenta):
            pass

        def adjust_forces(self, atoms, forces):
            positions = atoms.positions
            if self._type == 'plane':
                A, B, C, D = self.plane
                x, y, z = positions[self.index]
                d = (A * x + B * y + C * z + D) / np.sqrt(A**2 + B**2 + C**2)
                if d < 0:
                    return
                magnitude = self.spring * d
                direction = -np.array((A, B, C)) / np.linalg.norm((A, B, C))
                forces[self.index] += direction * magnitude
                return
            if self._type == 'two atoms':
                p1, p2 = positions[self.indices]
            elif self._type == 'point':
                p1 = positions[self.index]
                p2 = self.origin
            displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
            bondlength = np.linalg.norm(displace)
            if bondlength > self.threshold:
                magnitude = self.spring * (bondlength - self.threshold)
                direction = displace / np.linalg.norm(displace)
                if self._type == 'two atoms':
                    forces[self.indices[0]] += direction * magnitude
                    forces[self.indices[1]] -= direction * magnitude
                else:
                    forces[self.index] += direction * magnitude

        def adjust_potential_energy(self, atoms):
            """Returns the difference to the potential energy due to an active
            constraint. (That is, the quantity returned is to be added to the
            potential energy.)"""
            positions = atoms.positions
            if self._type == 'plane':
                A, B, C, D = self.plane
                x, y, z = positions[self.index]
                d = (A * x + B * y + C * z + D) / np.sqrt(A**2 + B**2 + C**2)
                if d > 0:
                    return 0.5 * self.spring * d**2
                else:
                    return 0.0
            if self._type == 'two atoms':
                p1, p2 = positions[self.indices]
            elif self._type == 'point':
                p1 = positions[self.index]
                p2 = self.origin
            displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
            bondlength = np.linalg.norm(displace)
            if bondlength > self.threshold:
                return 0.5 * self.spring * (bondlength - self.threshold) ** 2
            else:
                return 0.0

        def get_indices(self):
            if self._type == 'two atoms':
                return self.indices
            elif self._type == 'point':
                return self.index
            elif self._type == 'plane':
                return self.index

        def index_shuffle(self, atoms, ind):
            # See docstring of superclass
            if self._type == 'two atoms':
                newa = [-1, -1]  # Signal error
                for new, old in slice2enlist(ind, len(atoms)):
                    for i, a in enumerate(self.indices):
                        if old == a:
                            newa[i] = new
                if newa[0] == -1 or newa[1] == -1:
                    raise IndexError('Constraint not part of slice')
                self.indices = newa
            elif (self._type == 'point') or (self._type == 'plane'):
                newa = -1  # Signal error
                for new, old in slice2enlist(ind, len(atoms)):
                    if old == self.index:
                        newa = new
                        break
                if newa == -1:
                    raise IndexError('Constraint not part of slice')
                self.index = newa

        def __repr__(self):
            if self._type == 'two atoms':
                return 'Hookean(%d, %d)' % tuple(self.indices)
            elif self._type == 'point':
                return 'Hookean(%d) to cartesian' % self.index
            else:
                return 'Hookean(%d) to plane' % self.index