This post details how to rotate around a dihedral and save the various rotations to a file. For this we will be using the Python packages MDAnalysis, networkx, and nglview.

First, import these packages then create the Universe which contains your molecular structure. We will be working on an ethanol molecule. As we are working on an xyz file, there is no connectivity information given in this file, however we can use the method guess_bonds to guess this based on the positions of the atoms.

# for reading and writing our molecular simulation files
import MDAnalysis as mda
# for working with connectivity
import networkx as nx
# for visualising our result
import nglview

u = mda.Universe('./ethanol.xyz')
u.atoms.guess_bonds()

For this example, we want to rotate around the C-O bond. We can select this as so:

dih = u.dihedrals.select_bonds(('H', 'O', 'C', 'C'))[0]
bond = u.atoms.bonds.select_bonds(('C', 'O'))[0]

The [0] are necessary to convert from a TopologyGroup to the singular objects. The selections will normally return a Group which is similar to an array of bonds.

Next we need to split our molecule into two parts, the stationary and mobile parts. To do this, we can treat our molecule as a Graph and use networkX to identify the individual Graphs within.

g = nx.Graph()
g.add_edges_from(u.atoms.bonds.to_indices())
# remove the bond from the graph
# at this point the graph can be separated into two chunks
g.remove_edge(*bond.indices)

# unpack the two unconnected graphs
a, b = nx.connected_component_subgraphs(g)
# call the graph with 0 as a node the 'head'
head = a if 0 in a.nodes() else b

We can then return back into dealing with AtomGroups by slicing u.atoms with the indices of the nodes in our head group. The newly added set operations in MDAnalysis make it possible to invert the selection of head to give us the tail AtomGroup.

# use the node indices to slice u.atoms to form our first AtomGroup
head = u.atoms[head.nodes()]
# then use set logic to select all other atoms
tail = u.atoms ^ head

For our rotation, we will need to provide a vector to rotate around, and the center of rotation. The vector to rotate around can be found through considering the positions of the two atoms in the C-O bond. We want the center of rotation to be the atom which is in both the head group and the C-O bond, again we can use a set operation (&) to find the atom which is in both AtomGroups.

# calculate the vector which represents our bond
bvec = bond[1].position - bond[0].position

# find the center of our rotation.
# we are rotating `head` so want the atom in both head AND bond
center = (head & bond.atoms)[0]

We can set the initial rotation of the dihedral to 0, and then check its value:

head.rotateby(dih.value(), bvec, point=center.position)
dih.value()

With all our reference points defined, we can now write our trajectory. We can repeatedly apply the AtomGroup.rotateby method and write a new frame after each call.

with mda.Writer('out.pdb', n_atoms=len(u.atoms)) as w:
    for i in range(36):
        # rotateby takes an angle in degrees
        head.rotateby(10, bvec, point=center.position)
        w.write(u.atoms)

To check our results, we can visualise the trajectory we have just made inside our notebook using nglview:

new = mda.Universe('out.pdb')

v = nglview.show_mdanalysis(new)
v.add_ball_and_stick()
v

This will give something like this in the notebook itself: