import numpy as np
from ase.constraints.constraint import FixConstraint, slice2enlist
from ase.geometry.geometry import find_mic
[docs]
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
