Dihedral rotator function

mbuild/compound.py

@@ -2433,6 +2433,50 @@ def spin(self, theta, around):
         self.rotate(theta, around)
         self.translate(center_pos)
 
+    def rotate_dihedral(self, bond, phi):
+        """Rotate a dihedral about a central bond.
+
+        Parameters
+        ----------
+        bond : indexable object, length=2, dtype=mb.Compound
+            The pair of bonded Particles in the central bond of the dihedral
+        phi : float
+            The angle by which to rotate the dihedral, in radians.
+        """
+        nx = import_("networkx")
+
+        # Generate a bond graph and convert to networkX
+        mb_bondgraph = self.bond_graph
+        G = nx.Graph(mb_bondgraph.edges())
+
+        # Remove separate the compound in to two pieces by removing the bond
+        G.remove_edge(*bond)
+        assert len([i for i in nx.connected_components(G)]) == 2
+        components = [G.subgraph(c).copy() for c in nx.connected_components(G)]
+        component1 = components[1]  # One piece of the compound
+
+        # Get original coordinates
+        original_bond_positions = [bond[0].pos, bond[1].pos]
+
+        # Get the vector along the bond
+        bond_vec = bond[1].pos - bond[0].pos
+
+        # Rotate the coordinates of the piece by phi about the bond vector
+        xyz = np.array([p.pos for p in component1.nodes])
+        transformed_xyz = _rotate(xyz, phi, bond_vec)
+        for atom, coord in zip(component1.nodes, transformed_xyz):
+            atom.translate_to(coord)
+
+        # Move atoms involved in the bond to original positions
+        # This is neccessary since the piece is rotated about its center
+        if bond[0] in set(component1.nodes):
+            trans_vec = original_bond_positions[0] - bond[0].pos
+        elif bond[1] in set(component1.nodes):
+            trans_vec = original_bond_positions[1] - bond[1].pos
+
+        for atom in component1.nodes:
+            atom.translate(trans_vec)
+
     # Interface to GMSO Topology for reading/writing mol2 files
     def from_gmso(self, topology, coords_only=False, infer_hierarchy=True):
         """Convert a GMSO Topology to mBuild Compound.

mbuild/tests/test_coordinate_transform.py

@@ -346,6 +346,50 @@ def test_rotate_around_z_away_from_origin(self, sixpoints):
             before[:, 1], -1 * after[:, 1], atol=1e-16
         ) and np.allclose(before[:, 0], -1 * after[:, 0], atol=1e-16)
 
+    def test_rotate_dihedral(self, ethane):
+        bond = (ethane[0], ethane[4])
+        rotate_angle = np.deg2rad(60)
+        ethane.rotate_dihedral(bond, rotate_angle)
+
+        CH_vec1 = ethane[1].pos - ethane[0].pos
+        CH_vec2 = ethane[5].pos - ethane[4].pos
+        cos_dihedral = np.dot(CH_vec1, CH_vec2) / (
+            np.linalg.norm(CH_vec1) * np.linalg.norm(CH_vec2)
+        )
+        dihedral = np.rad2deg(np.arccos(cos_dihedral))
+        assert np.allclose(dihedral, 120, atol=1e-15)
+
+        # Extra test on asymmetric molecule
+        compound = mb.load("FCCO", smiles=True)
+
+        # Making sure this is the dihedral we are expecting
+        assert compound[0].element.symbol == "F"
+        assert compound[1].element.symbol == "C"
+        assert compound[2].element.symbol == "C"
+        assert compound[3].element.symbol == "O"
+
+        original_dihedral = [
+            compound[0].pos,
+            compound[1].pos,
+            compound[2].pos,
+            compound[3].pos,
+        ]
+        original_angle = dihedral_angle(original_dihedral)
+
+        bond = (compound[1], compound[2])
+        rotate_angle = np.deg2rad(68)
+        compound.rotate_dihedral(bond, rotate_angle)
+
+        new_dihedral = [
+            compound[0].pos,
+            compound[1].pos,
+            compound[2].pos,
+            compound[3].pos,
+        ]
+        new_angle = dihedral_angle(new_dihedral)
+
+        assert np.allclose(new_angle - original_angle, 68, atol=1e-15)
+
     def test_equivalence_transform(self, ch2, ch3, methane):
         ch2_atoms = list(ch2.particles())
         methane_atoms = list(methane.particles())
@@ -461,3 +505,18 @@ def test_bondgraph(self, ch3):
         for node in bgraph.nodes():
             x = map(lambda node: node.name, bgraph._adj[node])
             assert neighbors[node.name] == len(list(x))
+
+
+# Code to calculate angle of dihedral (stole from https://stackoverflow.com/questions/20305272/dihedral-torsion-angle-from-four-points-in-cartesian-coordinates-in-python#:~:text=The%20angle%20is%20given%20by,Similar%20to%20your%20dihedral2.)
+def dihedral_angle(p):
+    b1 = p[2] - p[1]
+    b0, b1, b2 = -(p[1] - p[0]), b1 / np.sqrt((b1 * b1).sum()), p[3] - p[2]
+    v = b0 - (b0[0] * b1[0] + b0[1] * b1[1] + b0[2] * b1[2]) * b1
+    w = b2 - (b2[0] * b1[0] + b2[1] * b1[1] + b2[2] * b1[2]) * b1
+    x = v[0] * w[0] + v[1] * w[1] + v[2] * w[2]
+    y = (
+        (b1[1] * v[2] - b1[2] * v[1]) * w[0]
+        + (b1[2] * v[0] - b1[0] * v[2]) * w[1]
+        + (b1[0] * v[1] - b1[1] * v[0]) * w[2]
+    )
+    return 180 * np.arctan2(y, x) / np.pi