<jupyter_start><jupyter_text>About&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Issues&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Tutorials&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Documentation![Molecular Design Toolkit](img/Top.png)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_: [<NAME>](https://github.com/avirshup), Autodesk Research - _Created on_: July 1, 2016 - _Tags_: DNA, molecular dynamics<jupyter_code>import moldesign as mdt
from moldesign import units as u

%matplotlib inline
from matplotlib.pyplot import *

# seaborn is optional -- it makes plots nicer
try: import seaborn  
except ImportError: pass<jupyter_output><empty_output><jupyter_text>Contents=======---   - [1. Create a DNA helix](1.-Create-a-DNA-helix)   - [2. Forcefield](2.-Forcefield)   - [3. Constraints](3.-Constraints)   - [4. MD Setup](4.-MD-Setup)   - [5. Minimization](5.-Minimization)   - [6. Dynamics](6.-Dynamics)   - [7. Analysis](7.-Analysis)  1. Create a DNA helix<jupyter_code>dna_structure = mdt.build_dna_helix('ACTGACTG', helix_type='b')
dna_structure.draw()
dna_structure<jupyter_output><empty_output><jupyter_text>2. ForcefieldThe cell below adds forcefield parameters to the molecule.**NOTE:** This molecule because is not missing expected atoms. If your molecule _is_ missing atoms (e.g., it's missing its hydrogens), use the `Forcefield.create_prepped_molecule` method instead of `Forcefield.assign`.**Click on the ERRORS/WARNING tab** to see any warnings raised during assignment.<jupyter_code>ff = mdt.forcefields.DefaultAmber()
ff.assign(dna_structure)<jupyter_output><empty_output><jupyter_text>3. ConstraintsThis section uses an interactive selection to constrain parts of the DNA.After executing the following cells, **click on the 3' and 5' bases:**<jupyter_code>rs = mdt.widgets.ResidueSelector(dna_structure)
rs
if len(rs.selected_residues) == 0:
    raise ValueError("You didn't click on anything!")
    
rs.selected_residues
for residue in rs.selected_residues:
    print('Constraining position for residue %s' % residue)
    
    for atom in residue.atoms:
        dna_structure.constrain_atom(atom)<jupyter_output><empty_output><jupyter_text>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 SetupThis section adds an OpenMM energy model and a Langevin integrator to the DNA.<jupyter_code>dna_structure.set_energy_model(mdt.models.OpenMMPotential,
                               implicit_solvent='obc')

dna_structure.set_integrator(mdt.integrators.OpenMMLangevin,
                             timestep=2.0*u.fs,
                             temperature=300.0*u.kelvin,
                             frame_interval=1.0*u.ps)<jupyter_output><empty_output><jupyter_text>You can interactively configure these methods:<jupyter_code>dna_structure.configure_methods()<jupyter_output><empty_output><jupyter_text>5. MinimizationNearly 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.<jupyter_code>trajectory = dna_structure.minimize(nsteps=200)
trajectory.draw()
plot(trajectory.potential_energy)

xlabel('steps');ylabel('energy / %s' % trajectory.unit_system.energy)
title('Energy relaxation'); grid('on')<jupyter_output><empty_output><jupyter_text>6. DynamicsWe'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.<jupyter_code>traj = dna_structure.run(run_for=25.0*u.ps)
traj.draw()<jupyter_output><empty_output><jupyter_text>7. AnalysisThe trajectory object (named `traj`) gives direct access to the timeseries data:<jupyter_code>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')
# Using the trajectory's 'plot' method will autogenerate axes labels with the appropriate units
traj.plot('time','kinetic_temperature')
title('Temperature'); grid('on')<jupyter_output><empty_output><jupyter_text>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**<jupyter_code>from ipywidgets import interact_manual
from IPython.display import display

rs = mdt.widgets.ResidueSelector(dna_structure)
def plot_rmsd(): 
    plot(traj.time, traj.rmsd(rs.selected_atoms))
    xlabel('time / fs'); ylabel(u'RMSD / Å')
interact_manual(plot_rmsd, description='plot rmsd')
rs<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10815568159593193, 'fraction_numerical': 0.009778994719342852, 'mean_word_length': 5.539641943734015, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 42, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '298936', 'license': 'Apache-2.0', 'n_tokens_mistral': 1666, 'n_tokens_neox': 1567, 'n_words': 433, 'path': 'moldesign/_notebooks/Example 1. Build and simulate DNA.ipynb', 'repo_name': 'Autodesk/molecular-design-toolkit', 'size': '5120'}