About Forum Issues Tutorials Documentation </span>

Example 1: Build and simulate DNA

This notebook builds a small DNA double helix, assigns a forcefield to it, and runs a molecular dynamics simulation.

Author: Aaron Virshup, Autodesk Research
Created on: July 1, 2016
Tags: DNA, molecular dynamics



In [ ]:

    
%matplotlib inline
from matplotlib.pyplot import *
try: import seaborn  #optional, makes plots nicer
except ImportError: pass

import moldesign as mdt
from moldesign import units as u

1. Create a DNA helix
2. Forcefield
3. Constraints
4. MD Setup
5. Minimization
6. Dynamics
7. Analysis

1. Create a DNA helix



In [ ]:

    
dna_structure = mdt.build_dna_helix('ACTGACTG', helix_type='b')
dna_structure.draw()



In [ ]:

    
dna_structure

2. Forcefield

Next, we'll assign forcefield parameters to our system. This will also add any missing atoms (typically hydrogens) to the molecule.

Click on the ERRORS/WARNING tab to see any warnings raised during assignment.



In [ ]:

    
mol = mdt.assign_forcefield(dna_structure)

3. Constraints

This section uses an interactive selection to constrain parts of the DNA.

After executing the following cells, click on the 3' and 5' bases:



In [ ]:

    
rs = mdt.widgets.ResidueSelector(mol)
rs



In [ ]:

    
if len(rs.selected_atoms) == 0:
    raise ValueError("You didn't click on anything!")
    
rs.selected_residues



In [ ]:

    
for residue in rs.selected_residues:
    print 'Constraining position for residue %s' % residue
    
    for atom in residue.atoms:
        mol.constrain_atom(atom)

Of course, fixing the positions of the terminal base pairs is a fairly extreme step. For extra credit, see if you can find a less heavy-handed keep the terminal base pairs bonded. (Try using tab-completion to see what other constraint methods are available)

4. MD Setup

This section adds an OpenMM energy model and a Langevin integrator to the DNA.



In [ ]:

    
mol.set_energy_model(mdt.models.OpenMMPotential,
                     implicit_solvent='obc',
                     cutoff=7.0 * u.angstrom)

mol.set_integrator(mdt.integrators.OpenMMLangevin,
                   timestep=2.0*u.fs,
                   temperature=300.0*u.kelvin,
                   frame_interval=1.0*u.ps)

You can interactively configure these methods:



In [ ]:

    
mol.configure_methods()

5. Minimization

Nearly every MD simulation should be preceded by an energy minimization, especially for crystal structure data. This will remove any energetically catastrophic clashes between atoms and prevent our simulation from blowing up.



In [ ]:

    
trajectory = mol.minimize(nsteps=200)
trajectory.draw()



In [ ]:

    
plot(trajectory.potential_energy)

xlabel('steps');ylabel('energy / %s' % trajectory.unit_system.energy)
title('Energy relaxation'); grid('on')

6. Dynamics

We're ready to run 25 picoseconds of dynamics at room temperature (that's 300º Kelvin). This will probably take a few minutes - if you're on an especially pokey computer, you might want to reduce the length of the simulation.



In [ ]:

    
traj = mol.run(run_for=25.0*u.ps)
traj.draw()

7. Analysis

The trajectory object (named traj) gives direct access to the timeseries data:



In [ ]:

    
plot(traj.time, traj.kinetic_energy, label='kinetic energy')
plot(traj.time, traj.potential_energy - traj.potential_energy[0], label='potential_energy')
xlabel('time / {time.units}'.format(time=traj.time))
ylabel('energy / {energy.units}'.format(energy=traj.kinetic_energy))
title('Energy vs. time'); legend(); grid('on')



In [ ]:

    
# Using the trajectory's 'plot' method will autogenerate axes labels with the appropriate units
traj.plot('time','kinetic_temperature')
title('Temperature'); grid('on')

This cell sets up an widget that plots the RMSDs of any selected group of atoms. Select a group of atoms, then click "Run plot_rmsd" to generate a plot



In [ ]:

    
from ipywidgets import interact_manual
from IPython.display import display

rs = mdt.widgets.ResidueSelector(mol)
def plot_rmsd(): 
    plot(traj.time, traj.rmsd(rs.selected_atoms))
    xlabel('time / fs'); ylabel(u'RMSD / Å')
interact_manual(plot_rmsd)
rs



In [ ]: