<jupyter_start><jupyter_text>Elastic network (harmonic springs) in OpenMMThis network shows how to add harmonic restraints between different atoms in OpenMM. In this case, we'll tie backbone atoms to the "$i+2$" backbone atom (i.e., next-nearest neighbor; nearest backbone atom not already bonded to this one).<jupyter_code>import openmm_scaled_md as scaled
import simtk.openmm as mm
from simtk.openmm import app
from simtk import unit

import mdtraj as md
pdb_file = '../resources/AD_initial_frame.pdb'
pdb = app.PDBFile(pdb_file)
forcefield = app.ForceField('amber96.xml', 'tip3p.xml')

pdb_system = forcefield.createSystem(pdb.topology,
                                     nonbondedMethod=app.PME, 
                                     nonbondedCutoff=1.0*unit.nanometers,
                                     constraints=app.HBonds,
                                     rigidWater=True,
                                     ewaldErrorTolerance=0.0005)<jupyter_output><empty_output><jupyter_text>The first two cells are identical to the normal setup.The next cells add restraints on the absolute positions. We'll use [MDTraj's atom selection language](http://mdtraj.org/latest/atom_selection.html) to select atoms; this *does* assume that MDTraj does not change the atom ordering scheme when creating an `mdtraj.Topology` from an `openmm.app.Topology`. We'll also assume that we want to restrain to the distances as given in the PDB.<jupyter_code>topology = md.Topology.from_openmm(pdb.topology)
pos_restrained_atoms = topology.select("backbone")
atom_pairs = [(pos_restrained_atoms[i], pos_restrained_atoms[i+2])
              for i in range(len(pos_restrained_atoms)-2)]
# calculate the distances based on the PDB file
traj = md.load(pdb_file)
default_distances = md.compute_distances(traj, atom_pairs=atom_pairs)[0]  # first frame in traj
elastic_network = mm.HarmonicBondForce()
length_unit = unit.nanometer
energy_unit = unit.kilojoule_per_mole
for ((atom_a, atom_b), r_0) in zip(atom_pairs, default_distances):
    elastic_network.addBond(atom_a, atom_b, 
                            length=float(r_0)*length_unit, 
                            k=5.0*energy_unit/length_unit**2)
pdb_system.addForce(elastic_network)<jupyter_output><empty_output><jupyter_text>Running with an elastic networkNow we'll create a trajectory that uses this system for dynamics. Getting the `simulation` object and using it to run MD is the same as is normally done.<jupyter_code># this is equivalent to the standard BAOAB integrator, since force_scaling is 1
integrator = scaled.integrators.BAOABIntegrator(
    temperature=300.0*unit.kelvin,
    collision_rate=1.0/unit.picosecond,
    timestep=2.0*unit.femtosecond,
    force_scaling=1.0
)
sim = app.Simulation(pdb.topology, pdb_system, integrator)
sim.context.setPositions(pdb.positions)
sim.reporters.append(app.PDBReporter('elastic_network.pdb', 10))
sim.step(1000)<jupyter_output><empty_output><jupyter_text>Visualizing with NGLViewIf you also have [NGLView](http://nglviewer.org/nglview/latest/) installed, you can visualize the trajectory in-notebook.<jupyter_code>import nglview as nv
traj = md.load("./elastic_network.pdb")
view = nv.show_mdtraj(traj)
view.add_ball_and_stick("ACE ALA NME")
view.add_point("water and .O")
view<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup