MDT Quickstart
or, A Hitchhiker's Guide to the Chemical Modeling Universe

This notebook gives you the Molecular Design Toolkit crash course.

Click on the following cell and press Shift + Return to execute it.



In [ ]:

    
# This cell sets up both the python and notebook environments
%matplotlib inline  
import moldesign as mdt  # import the buckyball package
from moldesign import units as u  # import the buckyball unit system

Molecules

Let's start with something relaxing.

Execute each cell by selecting it and pressing Shift + Return



In [ ]:

    
my_first_molecule = mdt.from_name('morphine')
my_first_molecule.draw()

You've just created your first molecule in buckyball - the bb.from_name method can translate many common chemical names, as well as IUPAC nomenclature, into 2D and 3D chemical structures.

If morphine is a little intimdating, you can try making caffeine, or even just water, by editing the code in the above cell. You can edit any code in this notebook at any time, and press Shift + Return to execute the new commands.

Let's take a look at your Molecule's attributes. Try using Jupyter's live autocompletion - typing my_first_molecule.[tab], for instance, will show you a list of the object's attributes.



In [ ]:

    
my_first_molecule

Every chemical system in buckyball is called a Molecule. These objects store information about atoms, chemical bonds, 3D positions, and chemical properties. A Molecule object is more than just a data store, though - they also have methods for simulating, modeling, and visualizing this information.

A list of the molecule's atoms are stored at my_first_molecule.atoms. Let's print out some custom information about the first five atoms:



In [ ]:

    
print 'This molecule has %d atoms' % my_first_molecule.num_atoms
for atom in my_first_molecule.atoms[:5]:
    print 'Atom {atom} (atomic number {atom.atnum}) weighs {atom.mass}, and ' \
    'is bonded to {atom.num_bonds} other atoms.'.format(atom=atom)



In [ ]:

    
my_first_molecule.atoms[-1]

Manipulating 3D structures

Sometimes, a 3D structure needs to be fixed - perhaps the wrong rotamer was generated, or someone gave you a non-planar geometry of a planar molecule.

For instance, the OpenBabel SMILES converter will generate a gauche-conformation alkane, when we really want an anti conformation. Let's use the GeometryBuilder to fix it - click on the central dihedral bond, and adjust the angle to 180º, either with the slider or by typing "180".



In [ ]:

    
mol = mdt.from_smiles('CCCC')
mdt.ui.GeometryBuilder(mol)

Readin', Writin', and Introspectin'

Of course, you'll probably want to import and export your own molecular structures. Buckyball can read and write molecules from a variety of 3D formats, including pdb, sdf, mol2, cif, and xyz.

Let's start by reading in a file that was downloaded from the PDB database (you can also use bb.from_pdb(3AID) to download it directly from the PDB site). We can get some quick information about the molecule in the file just by its name on the last line of the cell.



In [ ]:

    
hivprotease = mdt.read('data/3AID.pdb')
hivprotease

This is a crystal structure of a protein (HIV-1 protease) that's bound to a small organic drug molecule (that's the single unknown residue type). In addition to the protein and drug, there are also several water molecules hanging around the structure.



In [ ]:

    
hivprotease.draw()

Let's pull the drug molecule out of the structure and save it to disk. First, click on one of the atoms in the drug. You'll see that it's Residue ARQ401 in Chain A. We can grab this residue by going into the molecule's secondary structure:



In [ ]:

    
drugresidue = hivprotease.chains['A'].residues['ARQ401']



In [ ]:

    
drug = mdt.Molecule(drugresidue)
drug.draw3d()

Now, we'll write it to disk:



In [ ]:

    
drug.write(filename='/tmp/drugmol.xyz')



In [ ]:

    
!cat /tmp/drugmol.xyz

Simulation and modeling

As a quick first example, let's take a look at a quantum mechanical model of benzene.

This time, we'll construct the molecule in yet another way - by using a SMILES string.



In [ ]:

    
mol = mdt.from_smiles('C1=CC=CC=C1')
mol.draw()

Next, we associate the quantum mechanical model with the molecule, and run a calculation:



In [ ]:

    
mol.set_energy_model(mdt.models.PySCFPotential(theory='rhf',basis='sto-3g'))
mol.calculate()

There's a bunch of data here, which we'll dive into in other notebooks. For now, we can visualize the results of the calculation:



In [ ]:

    
mol.draw_orbitals()

You can visualize two different kinds of orbitals here: the canonical molecular orbitals, often just referred to as "molecular orbitals" (MOs); and atomic orbitals, which were combined to create the MOs.

Visualization API

You can access and create molecule visualizations in the notebook as well.

First, let's draw some more information into our visualization. We'll draw all the protein's potential hydrogen bond donors (basically Oxygens and Nitrogens) as spheres, and highlight those close to the ligand.



In [ ]:

    
viewer = hivprotease.draw3d()
hdonors = [atom for atom in hivprotease.atoms
           if atom.elem in ('O','N') and atom.residue.type == 'protein']

viewer.add_style(style='vdw', atoms=hdonors, radius=0.5)
viewer.highlight_atoms(atoms=[a for a in hdonors if a.distance(drugresidue) <= 4.0*u.angstrom])

viewer

The cell below constructs a quick visualization of the HIV drug's contacts with the surrounding protein residues, with 2D rendering, a 3D rendering, and an inspector that provides information about atoms that the user clicks on: