About Issues Tutorials Documentation
</span>
This notebook takes a crystal structure from the PDB and prepares it for simulation.
In [ ]:
%matplotlib inline
from matplotlib.pyplot import *
import moldesign as mdt
from moldesign import units as u
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.
In [ ]:
xtal = mdt.from_pdb('1kbu')
xtal.draw()
In [ ]:
assembly = mdt.build_assembly(xtal, 1)
assembly.draw()
By evaulating the assembly object (it's a normal instance of the moldesign.Molecule class), we can get some information about it's content:
In [ ]:
assembly
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.
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:
In [ ]:
dna_atoms = assembly.get_atoms('dna')
dna_only = mdt.Molecule(dna_atoms)
dna_only.draw3d(display=True)
dna_only
In [ ]:
ff = mdt.forcefields.DefaultAmber()
dna = ff.create_prepped_molecule(dna_only)
In [ ]:
dna.set_energy_model(mdt.models.OpenMMPotential, implicit_solvent='obc')
dna.configure_methods()
In [ ]:
minimization = dna.minimize()
In [ ]:
minimization.draw()
In [ ]:
# Freeze the backbone:
for residue in dna.residues:
for atom in residue.backbone:
dna.constrain_atom(atom)
In [ ]:
dna.set_integrator(mdt.integrators.OpenMMLangevin,
timestep=2.0*u.fs,
frame_interval=1.0*u.ps,
remove_rotation=True)
dna.integrator.configure()
And now we run it. This is may take a while, depending on your hardware.
In [ ]:
equil1 = dna.run(20.0*u.ps)
In [ ]:
equil1.draw()
Next, we'll remove the constraints and do full dynamics:
In [ ]:
dna.clear_constraints()
equil2 = dna.run(20.0*u.ps)
In [ ]:
equil = equil1 + equil2
equil.draw()
In [ ]:
plot(equil2.time, equil2.rmsd())
xlabel('time / fs'); ylabel(u'rmsd / Å'); grid()
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.
In [ ]:
trajectory = dna.run(40.0*u.ps)
In [ ]:
trajectory.draw()
In [ ]:
trajectory.write('holliday_traj.P.gz')
To load the saved object, use:
In [ ]:
traj = mdt.read('holliday_traj.P.gz')
In [ ]:
traj.draw()