<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <span style="float
Step1:  Contents
A. View the crystal structure
B. Build the biomolecular assembly
C. Isolate the DNA
D. Prep for simulation
E. Dynamics - equilibration
F. Dynamics - production
A. View the crystal structure
We start by downloading the 1KBU crystal structure.
It will generate several warnings. Especially note that it contains biomolecular "assembly" information. This means that the file from PDB doesn't contain the complete structure, but we can generate the missing parts using symmetry operations.
Step2:  B. Build the biomolecular assembly
As you can read in the warning, 1KBU only has one biomolecular assembly, conveniently named '1'. This cell builds and views it
Step3:  By evaulating the assembly object (it's a normal instance of the moldesign.Molecule class), we can get some information about it's content
Step4:  Because we're only interested in DNA, we'll create a new molecule using only the DNA residues, and then assign a forcefield to it.
C. Isolate the DNA
This example will focus only on the DNA components of this structure, so we'll isolate the DNA atoms and create a new molecule from them.
We could do this with a list comprehension, e.g.
mdt.Molecule([atom for atom in assembly.atoms if atom.residue.type == 'dna'])
Here, however we'll use a shortcut for this - the molecule.get_atoms method, which allows you to run queries on the atoms
Step5:  D. Prep for simulation
Next, we'll assign a forcefield and energy model, then minimize the structure.
Step6:  E. Dynamics - equilibration
The structure is ready. We'll associate an integrator with the molecule, then do a 2 step equilibration - first freezing the peptide backbone and running 300K dynamics, then unfreezing and continuing dyanmics.
Step7:  And now we run it. This is may take a while, depending on your hardware.
Step8:  Next, we'll remove the constraints and do full dynamics
Step9:  NOTE
Step10:  G. Save your results
Any MDT object can be saved to disk. We recommend saving objects with the "Pickle" format to make sure that all the data is preserved.
This cell saves the final trajectory to disk as a compressed pickle file
Step11:  To load the saved object, use
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from matplotlib.pyplot import *
import moldesign as mdt
from moldesign import units as u
Explanation: <span style="float:right"><a href="http://moldesign.bionano.autodesk.com/" target="_blank" title="About">About</a>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="https://github.com/autodesk/molecular-design-toolkit/issues" target="_blank" title="Issues">Issues</a>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="http://bionano.autodesk.com/MolecularDesignToolkit/explore.html" target="_blank" title="Tutorials">Tutorials</a>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="http://autodesk.github.io/molecular-design-toolkit/" target="_blank" title="Documentation">Documentation</a></span>
</span>
<br>
<center><h1>Example 3: Simulating a Holliday Junction PDB assembly </h1> </center>
This notebook takes a crystal structure from the PDB and prepares it for simulation.
Author: Aaron Virshup, Autodesk Research
Created on: July 1, 2016
Tags: DNA, holliday junction, assembly, PDB, MD
End of explanation
xtal = mdt.from_pdb('1kbu')
xtal.draw()
Explanation: Contents
A. View the crystal structure
B. Build the biomolecular assembly
C. Isolate the DNA
D. Prep for simulation
E. Dynamics - equilibration
F. Dynamics - production
A. View the crystal structure
We start by downloading the 1KBU crystal structure.
It will generate several warnings. Especially note that it contains biomolecular "assembly" information. This means that the file from PDB doesn't contain the complete structure, but we can generate the missing parts using symmetry operations.
End of explanation
assembly = mdt.build_assembly(xtal, 1)
assembly.draw()
Explanation: B. Build the biomolecular assembly
As you can read in the warning, 1KBU only has one biomolecular assembly, conveniently named '1'. This cell builds and views it:
End of explanation
assembly
Explanation: By evaulating the assembly object (it's a normal instance of the moldesign.Molecule class), we can get some information about it's content:
End of explanation
dna_atoms = assembly.get_atoms('dna')
dna_only = mdt.Molecule(dna_atoms)
dna_only.draw3d(display=True)
dna_only
Explanation: Because we're only interested in DNA, we'll create a new molecule using only the DNA residues, and then assign a forcefield to it.
C. Isolate the DNA
This example will focus only on the DNA components of this structure, so we'll isolate the DNA atoms and create a new molecule from them.
We could do this with a list comprehension, e.g.
mdt.Molecule([atom for atom in assembly.atoms if atom.residue.type == 'dna'])
Here, however we'll use a shortcut for this - the molecule.get_atoms method, which allows you to run queries on the atoms:
End of explanation
ff = mdt.forcefields.DefaultAmber()
dna = ff.create_prepped_molecule(dna_only)
dna.set_energy_model(mdt.models.OpenMMPotential, implicit_solvent='obc')
dna.configure_methods()
minimization = dna.minimize()
minimization.draw()
Explanation: D. Prep for simulation
Next, we'll assign a forcefield and energy model, then minimize the structure.
End of explanation
# Freeze the backbone:
for residue in dna.residues:
    for atom in residue.backbone:
        dna.constrain_atom(atom)
dna.set_integrator(mdt.integrators.OpenMMLangevin,
                   timestep=2.0*u.fs,
                   frame_interval=1.0*u.ps,
                   remove_rotation=True)
dna.integrator.configure()
Explanation: E. Dynamics - equilibration
The structure is ready. We'll associate an integrator with the molecule, then do a 2 step equilibration - first freezing the peptide backbone and running 300K dynamics, then unfreezing and continuing dyanmics.
End of explanation
equil1 = dna.run(20.0*u.ps)
equil1.draw()
Explanation: And now we run it. This is may take a while, depending on your hardware.
End of explanation
dna.clear_constraints()
equil2 = dna.run(20.0*u.ps)
equil = equil1 + equil2
equil.draw()
plot(equil2.time, equil2.rmsd())
xlabel('time / fs'); ylabel(u'rmsd / Å'); grid()
Explanation: Next, we'll remove the constraints and do full dynamics:
End of explanation
trajectory = dna.run(40.0*u.ps)
trajectory.draw()
Explanation: NOTE: THIS IS NOT A SUFFICIENT EQUILIBRATION FOR PRODUCTION MOLECULAR DYNAMICS! 
In practice, before going to "production", we would at least want to run dynamics until the RMSD and thermodynamic observabled have converged. A variety of equilibration protocols are used in practice, including slow heating, reduced coupling, multiple constraints, etc.
F. Dynamics - production
Assuming that we're satisfied with our system's equilibration, we now gather data for "production". This will take a while.
End of explanation
trajectory.write('holliday_traj.P.gz')
Explanation: G. Save your results
Any MDT object can be saved to disk. We recommend saving objects with the "Pickle" format to make sure that all the data is preserved.
This cell saves the final trajectory to disk as a compressed pickle file:
End of explanation
traj = mdt.read('holliday_traj.P.gz')
traj.draw()
Explanation: To load the saved object, use:
End of explanation
<END_TASK>