!pip install -q condacolab
import condacolab
condacolab.install()

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 2
      1 get_ipython().system('pip install -q condacolab')
----> 2 import condacolab
      3 condacolab.install()

File /opt/hostedtoolcache/Python/3.12.12/x64/lib/python3.12/site-packages/condacolab.py:22
     20 from typing import Dict, AnyStr
     21 from urllib.request import urlopen
---> 22 from distutils.spawn import find_executable
     24 from IPython import get_ipython
     26 try:

ModuleNotFoundError: No module named 'distutils'

%%capture
!conda install -c conda-forge openmm mdtraj parmed
!pip install py3dmol

import openmm as mm
import openmm.app as app
import openmm.unit as unit
import numpy as np
from sys import stdout
from parmed import Structure, Atom, Bond

Create Structure and Topology objects¶

# System parameters
M = 100   # Number of chains
N = 10    # Atoms per chain
mass = 120.0  # Mass in amu

# Create an empty structure
structure = Structure()

# Add atoms and residues
for m in range(M):
    for n in range(N):
        atom = Atom(name='CG', type='CG', atomic_number=0, mass=mass, charge=0.0)
        resname = f"CG{m}"
        resnum = m + 1
        structure.add_atom(atom=atom, resname=resname, resnum=resnum, chain='A')

# Assign simple positions: spaced along z-axis
positions = []
for m in range(M):
    x0 = np.array(((m%10)*1.0, (m//10)*1.0, 0))
    positions.append(x0)
    for i in range(1,N):
        xi = positions[-1] + np.array((0, 0, 0.38))
        positions.append(xi)

# Convert the list into OpenMM Quantity with units of nm
structure.positions = positions * unit.nanometer

# Add bonds between consecutive atoms in each chain
for m in range(M):
    offset = m * N
    for n in range(N - 1):
        i = offset + n
        structure.bonds.append(Bond(structure.atoms[i], structure.atoms[i + 1]))

# Save the files
structure.save('cg_polymer.pdb', overwrite=True)
structure.save('cg_polymer.psf', overwrite=True)

# Check basic info
print(structure.topology)
structure.to_dataframe().head()

Defining the force field¶

We will cover two methods for defining the forcefield:

Using the Python API.
Creating a forcefield xml file.

Both methods achieve the same result of defining the forcefield. There are pros and cons of both ways. The key points are that using the Python API is more flexible, while creating a forcefield xml file makes it easier to share your forcefield and helps reproducibility. This is demonstrated in the other OpenMM tutorials --- the standard forcefields such as Amber are distributed with OpenMM as xml files, while the custom forces used to do free energy calculations are defined using the Python API.

Option-1: Defining a force-field using the Python API¶

We can create instances of the OpenMM Force classes, assign parameters, and add them to the system. First we must create a System.

# create the system and add the particles to it
system = mm.System()

for i in range(len(structure.atoms)):
  system.addParticle(structure.atoms[i].mass)

box_vecs = 11 * np.eye(3) *unit.nanometer
system.setDefaultPeriodicBoxVectors(*box_vecs )

harmonic_bond_force = mm.HarmonicBondForce()

# Add each bond to the force from the topology
for bond in structure.topology.bonds():
    harmonic_bond_force.addBond(bond.atom1.index, bond.atom2.index, 0.38, 1000)

# Define a Lennard-Jones potential
expression = '4*epsilon*((sigma/r)^12-(sigma/r)^6);'\
            + ' sigma=0.5*(sigma1+sigma2);'\
            + ' epsilon=sqrt(epsilon1*epsilon2)'

custom_nb_force = mm.CustomNonbondedForce(expression)

custom_nb_force.addPerParticleParameter('sigma')
custom_nb_force.addPerParticleParameter('epsilon')

# Add the parameters for each particle
for atom in structure.topology.atoms():
    custom_nb_force.addParticle([0.5, 1.0])

# Create exclusions for directly bonded atoms
custom_nb_force.createExclusionsFromBonds([(bond[0].index, bond[1].index) for bond in structure.topology.bonds()], 1)

# Set a cutoff of 1.5nm
custom_nb_force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic)
custom_nb_force.setCutoffDistance(1.5*unit.nanometers)

# Add the forces to the system
system.addForce(harmonic_bond_force)
system.addForce(custom_nb_force)

Option-2" Creating a force-field xml file¶

Alternatively, we can create a custom force-field xml file for our system. You should look here first for the full information: http://docs.openmm.org/latest/userguide/application/05_creating_ffs.html.

You will need to create a new file called “cg_ff.xml” using a text editor of your choice. Then paste the following lines into it:

<!-- cg_ff.xml 
Coarse-grained forcefield for a bead-spring polymer. -->
<ForceField>
	
    <AtomTypes>
	    <Type name="CG-bead" class="CG-bead" element="CG" mass="120.0"/>
	</AtomTypes>

	<Residues>
        <!-- Each bead is a single residue.
        Need a template for the different Residue types. 
        First type is the end. This only has one bond. -->
        <Residue name="CG-residue-end">
            <Atom name="CG-bead" type="CG-bead"/>
            <ExternalBond atomName="CG-bead"/>
        </Residue>

        <!-- Second type is in the middle of the chain. This has two bonds. -->
        <Residue name="CG-residue-middle">
            <Atom name="CG-bead" type="CG-bead"/>
            <ExternalBond atomName="CG-bead"/>
            <ExternalBond atomName="CG-bead"/>
        </Residue>
    </Residues>

    <!-- Define a harmonic bond between the CA-beads -->
    <HarmonicBondForce>
        <Bond class1="CG-bead" class2="CG-bead" length="0.38" k="1000.0"/>
    </HarmonicBondForce>

    <!-- Use a custom non-bonded force for maximum flexibility.
    The bondCutoff=1 tells it to only exclude interactions between directly bonded atoms. -->
    <CustomNonbondedForce energy="4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); epsilon=sqrt(epsilon1*epsilon2)" bondCutoff="1">
        <PerParticleParameter name="sigma"/>
        <PerParticleParameter name="epsilon"/>
		<Atom type="CG-bead" sigma="0.5" epsilon="1.0"/>
	</CustomNonbondedForce>
</ForceField>

(If you have cloned the cookbook repo then this file will already exist) The comments in the above code explain what the difference sections are for.

Create the system¶

We can now load in our previously created custom ForceField and use the createSystem method with the topology.

# load in the ForceField we created
ff = app.ForceField('cg_ff.xml')
system2 = ff.createSystem(topology, nonbondedMethod=app.CutoffPeriodic, nonbondedCutoff=1.5*unit.nanometers)

We can compare the two systems to check if they are the same. A simple way to do this is to serialize the systems as save them as xml files.

with open('system1.xml', 'w') as output:
    output.write(mm.XmlSerializer.serialize(system))

with open('system2.xml', 'w') as output:
    output.write(mm.XmlSerializer.serialize(system2))

!diff system1.xml system2.xml

We actually find they are slightly different. system2 created by ForceField.createSystem has a CMMotionRemover that was added by default when using the createSystem method. If we wanted to we could add this to the first system.

Setup and run the simulation¶

We will use a LangevinMiddleIntegrator with a friction term of 0.01ps^-1 and a timestep of 0.01ps as used in similar coarse-grained polymer models [1]. For your own models these parameters will be important and you will need to choose them carefully!

integrator = mm.LangevinMiddleIntegrator(300*unit.kelvin, 0.01/unit.picosecond, 0.010*unit.picoseconds)
simulation = app.Simulation(structure.topology, system, integrator)
simulation.context.setPositions(positions)

# setup simulation reporters
# Write the trajectory to a file called 'traj.dcd'
simulation.reporters.append(app.DCDReporter('traj.dcd', 1000, enforcePeriodicBox=False))

# Report information to the screen as the simulation runs
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True,
        potentialEnergy=True, temperature=True, volume=True, speed=True))


# NVT equilibration
simulation.step(10000)

# Add a barostat
barostatIndex=system.addForce(mm.MonteCarloBarostat(1.0*unit.bar, 300*unit.kelvin))
simulation.context.reinitialize(preserveState=True)

# Run NPT equilibration
simulation.step(100000)

Visualize Simulation¶

# Animation
view = py3Dmol.view(width=800, height=800)

view.addModelsAsFrames(open('ljtraj.pdb', 'r').read(),'pdb')
view.setBackgroundColor('white')
view.setStyle({'sphere': {'color':'green'}}) 

#
view.zoomTo()
view.animate({'loop': "forward"})
view.show()

import mdtraj as md
import nglview as nv
from IPython.display import display

# Load trajectory with PSF as topology
traj = md.load('traj.dcd', top='equilibrated_config.pdb')

# Visualize with connectivity
view = nv.show_mdtraj(traj)
view.add_representation('ball+stick', selection='all')  # or 'licorice'

# Display the viewer
display(view)