Ethane: Reading from files, Ports and coordinate transforms

Note: mBuild expects all distance units to be in nanometers.

In this example, we’ll cover reading molecular components from files, introduce the concept of Ports and start using some coordinate transforms.

First, we need to import the mbuild package:

import mbuild as mb

As you probably noticed while creating your methane molecule in the last tutorial, manually adding Particles and Bonds to a Compound is a bit cumbersome. The easiest way to create small, reusable components, such as methyls, amines or monomers, is to hand draw them using software like Avogadro and export them as either a .pdb or .mol2 file (the file should contain connectivity information).

Let’s start by reading a methyl group from a .pdb file:

import mbuild as mb

ch3 = mb.load('ch3.pdb')
ch3.visualize()

Now let’s use our first coordinate transform to center the methyl at its carbon atom:

import mbuild as mb

ch3 = mb.load('ch3.pdb')
ch3.translate(-ch3[0].pos)  # Move carbon to origin.

Now we have a methyl group loaded up and centered. In order to connect Compounds in mBuild, we make use of a special type of Compound: the Port. A Port is a Compound with two sets of four “ghost” Particles that assist in bond creation. In addition, Ports have an anchor attribute which typically points to a particle that the Port should be associated with. In our methyl group, the Port should be anchored to the carbon atom so that we can now form bonds to this carbon:

import mbuild as mb

ch3 = mb.load('ch3.pdb')
ch3.translate(-ch3[0].pos)  # Move carbon to origin.

port = mb.Port(anchor=ch3[0])
ch3.add(port, label='up')

# Place the port at approximately half a C-C bond length.
ch3['up'].translate([0, -0.07, 0])

By default, Ports are never output from the mBuild structure. However, it can be useful to look at a molecule with the Ports to check your work as you go:

ch3.visualize(show_ports=True)

Now we wrap the methyl group into a python class, so that we can reuse it as a component to build more complex molecules later.

import mbuild as mb

class CH3(mb.Compound):
    def __init__(self):
        super(CH3, self).__init__()

        mb.load('ch3.pdb', compound=self)
        self.translate(-self[0].pos)  # Move carbon to origin.

        port = mb.Port(anchor=self[0])
        self.add(port, label='up')
        # Place the port at approximately half a C-C bond length.
        self['up'].translate([0, -0.07, 0])

When two Ports are connected, they are forced to overlap in space and their parent Compounds are rotated and translated by the same amount.

Note: If we tried to connect two of our methyls right now using only one set of four ghost particles, not only would the Ports overlap perfectly, but the carbons and hydrogens would also perfectly overlap - the 4 ghost atoms in the Port are arranged identically with respect to the other atoms. For example, if a Port and its direction is indicated by “<-”, forcing the port in <-CH3 to overlap with <-CH3 would just look like <-CH3 (perfectly overlapping atoms).

To solve this problem, every port contains a second set of 4 ghost atoms pointing in the opposite direction. When two Compounds are connected, the port that places the anchor atoms the farthest away from each other is chosen automatically to prevent this overlap scenario.

When <->CH3 and <->CH3 are forced to overlap, the CH3<->CH3 is automatically chosen.

Now the fun part: stick ’em together to create an ethane:

ethane = mb.Compound()

ethane.add(CH3(), label="methyl_1")
ethane.add(CH3(), label="methyl_2")
mb.force_overlap(move_this=ethane['methyl_1'],
                         from_positions=ethane['methyl_1']['up'],
                         to_positions=ethane['methyl_2']['up'])

Above, the force_overlap() function takes a Compound and then rotates and translates it such that two other Compounds overlap. Typically, as in this case, those two other Compounds are Ports - in our case, methyl1['up'] and methyl2['up'].

ethane.visualize()

ethane.visualize(show_ports=True)

Similarly, if we want to make ethane a reusable component, we need to wrap it into a python class.

import mbuild as mb

class Ethane(mb.Compound):
    def __init__(self):
        super(Ethane, self).__init__()

        self.add(CH3(), label="methyl_1")
        self.add(CH3(), label="methyl_2")
        mb.force_overlap(move_this=self['methyl_1'],
                         from_positions=self['methyl_1']['up'],
                         to_positions=self['methyl_2']['up'])

ethane = Ethane()
ethane.visualize()

# Save to .mol2
ethane.save('ethane.mol2', overwrite=True)