Skip to content

Rigid bodies

ASE constraint that keeps nodes and linkers rigid during optimisation: cage_isomer_builder.utils.rigid_body.

FixRigidBodies 

FixRigidBodies(groups: Sequence[Sequence[int]], group_masses, group_ref_rel)

Bases: FixConstraint

Constrain groups of atoms to move only as rigid bodies.

Each group's reference shape (atom positions relative to its mass-weighted centroid) never changes; the constrained optimiser may only translate and rotate each group as a whole. Build one from a structure with :meth:FixRigidBodies.from_atoms. The constructor takes the reference data directly (not an Atoms object) so that todict() and dict2constraint round-trip, as ASE's Trajectory writer requires.

Parameters:

Name Type Description Default
groups sequence of sequence of int

Atom indices per rigid body.

 required
group_masses sequence of sequence of float

Per-atom mass, aligned with groups.

 required
group_ref_rel sequence of (N, 3) array-like

Reference positions of each group's atoms relative to its mass-weighted centroid, aligned with groups.

 required

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
>>> a, b = molecule("H2O"), molecule("H2O")
>>> b.translate([3.0, 0.0, 0.0])
>>> pair = a + b                               # two waters, atoms 0-2 and 3-5
>>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
>>> pair.set_constraint(constraint)
>>> new = pair.positions.copy()
>>> new[1] += [0.3, 0.0, 0.0]                  # try to stretch an O-H bond
>>> pair.set_positions(new)                    # the constraint keeps each water rigid
>>> round(float(pair.get_distance(0, 1)), 4) == round(float(a.get_distance(0, 1)), 4)
True
Source code in cage_isomer_builder/utils/rigid_body.py 
60
61
62
63
def __init__(self, groups: Sequence[Sequence[int]], group_masses, group_ref_rel):
    self.group_indices = [np.asarray(g, dtype=int) for g in groups]
    self.group_masses = [np.asarray(m, dtype=float) for m in group_masses]
    self.group_ref_rel = [np.asarray(r, dtype=float) for r in group_ref_rel]

from_atoms classmethod 

from_atoms(atoms: Atoms, groups: Sequence[Sequence[int]]) -> FixRigidBodies

Build from a structure, keeping each group's current shape.

Groups with fewer than 2 atoms are dropped: a single atom has no shape to keep.

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
>>> a, b = molecule("H2O"), molecule("H2O")
>>> b.translate([3.0, 0.0, 0.0])
>>> pair = a + b                               # two waters, atoms 0-2 and 3-5
>>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
>>> [g.tolist() for g in constraint.group_indices]
[[0, 1, 2], [3, 4, 5]]
Source code in cage_isomer_builder/utils/rigid_body.py 
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
@classmethod
def from_atoms(
    cls, atoms: Atoms, groups: Sequence[Sequence[int]]
) -> "FixRigidBodies":
    """
    Build from a structure, keeping each group's current shape.

    Groups with fewer than 2 atoms are dropped: a single atom has no
    shape to keep.

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
    >>> a, b = molecule("H2O"), molecule("H2O")
    >>> b.translate([3.0, 0.0, 0.0])
    >>> pair = a + b                               # two waters, atoms 0-2 and 3-5
    >>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
    >>> [g.tolist() for g in constraint.group_indices]
    [[0, 1, 2], [3, 4, 5]]
    """
    masses = atoms.get_masses()
    positions = atoms.get_positions()

    kept_groups: list[np.ndarray] = []
    kept_masses: list[np.ndarray] = []
    kept_ref_rel: list[np.ndarray] = []

    for group in groups:
        indices = np.asarray(group, dtype=int)
        if len(indices) < 2:
            continue
        m = masses[indices]
        pos = positions[indices]
        centroid = (m[:, None] * pos).sum(axis=0) / m.sum()
        kept_groups.append(indices)
        kept_masses.append(m)
        kept_ref_rel.append(pos - centroid)

    return cls(kept_groups, kept_masses, kept_ref_rel)

get_removed_dof 

get_removed_dof(atoms)

Degrees of freedom removed: 3N - 6 per rigid group of N atoms (ASE constraint interface).

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
>>> a, b = molecule("H2O"), molecule("H2O")
>>> b.translate([3.0, 0.0, 0.0])
>>> pair = a + b                               # two waters, atoms 0-2 and 3-5
>>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
>>> constraint.get_removed_dof(pair)           # 2 waters x (9 - 6)
6
Source code in cage_isomer_builder/utils/rigid_body.py 
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def get_removed_dof(self, atoms):
    """
    Degrees of freedom removed: 3N - 6 per rigid group of N atoms
    (ASE constraint interface).

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
    >>> a, b = molecule("H2O"), molecule("H2O")
    >>> b.translate([3.0, 0.0, 0.0])
    >>> pair = a + b                               # two waters, atoms 0-2 and 3-5
    >>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
    >>> constraint.get_removed_dof(pair)           # 2 waters x (9 - 6)
    6
    """
    return sum(max(0, 3 * len(g) - 6) for g in self.group_indices)

adjust_positions 

adjust_positions(atoms, new)

Replace new positions in place by the closest rigid-body move of each group (ASE constraint interface).

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
>>> a, b = molecule("H2O"), molecule("H2O")
>>> b.translate([3.0, 0.0, 0.0])
>>> pair = a + b                               # two waters, atoms 0-2 and 3-5
>>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
>>> new = pair.positions.copy()
>>> new[0] += [0.0, 0.0, 0.5]                  # move one atom only
>>> constraint.adjust_positions(pair, new)
>>> bool(np.isclose(np.linalg.norm(new[0] - new[1]), a.get_distance(0, 1)))
True
Source code in cage_isomer_builder/utils/rigid_body.py 
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def adjust_positions(self, atoms, new):
    """
    Replace ``new`` positions in place by the closest rigid-body move of
    each group (ASE constraint interface).

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
    >>> a, b = molecule("H2O"), molecule("H2O")
    >>> b.translate([3.0, 0.0, 0.0])
    >>> pair = a + b                               # two waters, atoms 0-2 and 3-5
    >>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
    >>> new = pair.positions.copy()
    >>> new[0] += [0.0, 0.0, 0.5]                  # move one atom only
    >>> constraint.adjust_positions(pair, new)
    >>> bool(np.isclose(np.linalg.norm(new[0] - new[1]), a.get_distance(0, 1)))
    True
    """
    for indices, masses, ref_rel in zip(
        self.group_indices, self.group_masses, self.group_ref_rel
    ):
        trial = new[indices]
        total_mass = masses.sum()
        trial_centroid = (masses[:, None] * trial).sum(axis=0) / total_mass
        trial_rel = trial - trial_centroid

        rotation, _ = Rotation.align_vectors(trial_rel, ref_rel, weights=masses)
        new[indices] = trial_centroid + rotation.apply(ref_rel)

adjust_forces 

adjust_forces(atoms, forces)

Replace forces in place by each group's net force and torque, spread back over its atoms as a rigid-body motion (ASE constraint interface).

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
>>> a, b = molecule("H2O"), molecule("H2O")
>>> b.translate([3.0, 0.0, 0.0])
>>> pair = a + b                               # two waters, atoms 0-2 and 3-5
>>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
>>> forces = np.zeros((6, 3))
>>> forces[1] = [1.0, 0.0, 0.0]                # push one H only
>>> constraint.adjust_forces(pair, forces)
>>> round(float(forces[:3].sum(axis=0)[0]), 6)   # the whole water feels the net force
1.0
Source code in cage_isomer_builder/utils/rigid_body.py 
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
def adjust_forces(self, atoms, forces):
    # Newton-Euler rigid-body force projection: keep only the force
    # components that produce net translation/rotation of the group as a
    # whole, discarding anything that would deform it internally.
    """
    Replace forces in place by each group's net force and torque, spread
    back over its atoms as a rigid-body motion (ASE constraint
    interface).

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
    >>> a, b = molecule("H2O"), molecule("H2O")
    >>> b.translate([3.0, 0.0, 0.0])
    >>> pair = a + b                               # two waters, atoms 0-2 and 3-5
    >>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
    >>> forces = np.zeros((6, 3))
    >>> forces[1] = [1.0, 0.0, 0.0]                # push one H only
    >>> constraint.adjust_forces(pair, forces)
    >>> round(float(forces[:3].sum(axis=0)[0]), 6)   # the whole water feels the net force
    1.0
    """
    positions = atoms.get_positions()
    for indices, masses in zip(self.group_indices, self.group_masses):
        pos = positions[indices]
        total_mass = masses.sum()
        centroid = (masses[:, None] * pos).sum(axis=0) / total_mass
        d = pos - centroid

        f = forces[indices]
        f_net = f.sum(axis=0)
        torque = np.cross(d, f).sum(axis=0)

        inertia = np.zeros((3, 3))
        for m_i, d_i in zip(masses, d):
            inertia += m_i * (np.dot(d_i, d_i) * np.eye(3) - np.outer(d_i, d_i))
        # pinv handles the singular case (e.g. a linear/2-atom group,
        # where rotation about the inter-atom axis is undefined).
        alpha = np.linalg.pinv(inertia) @ torque

        f_trans = (masses[:, None] / total_mass) * f_net[None, :]
        f_rot = masses[:, None] * np.cross(alpha[None, :], d)
        forces[indices] = f_trans + f_rot

todict 

todict()

Dictionary form, as used by ASE to store the constraint in a trajectory.

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
>>> a, b = molecule("H2O"), molecule("H2O")
>>> b.translate([3.0, 0.0, 0.0])
>>> pair = a + b                               # two waters, atoms 0-2 and 3-5
>>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
>>> constraint.todict()["name"]
'FixRigidBodies'
Source code in cage_isomer_builder/utils/rigid_body.py 
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
def todict(self):
    """
    Dictionary form, as used by ASE to store the constraint in a trajectory.

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> from cage_isomer_builder.utils.rigid_body import FixRigidBodies
    >>> a, b = molecule("H2O"), molecule("H2O")
    >>> b.translate([3.0, 0.0, 0.0])
    >>> pair = a + b                               # two waters, atoms 0-2 and 3-5
    >>> constraint = FixRigidBodies.from_atoms(pair, [[0, 1, 2], [3, 4, 5]])
    >>> constraint.todict()["name"]
    'FixRigidBodies'
    """
    return {
        "name": "FixRigidBodies",
        "kwargs": {
            "groups": [g.tolist() for g in self.group_indices],
            "group_masses": [m.tolist() for m in self.group_masses],
            "group_ref_rel": [r.tolist() for r in self.group_ref_rel],
        },
    }