!pip install -q condacolab
import condacolab
condacolab.install()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/condacolab.py:26
25 try:
---> 26 import google.colab
27 except ImportError:
ModuleNotFoundError: No module named 'google'
During handling of the above exception, another exception occurred:
RuntimeError 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.8.18/x64/lib/python3.8/site-packages/condacolab.py:28
26 import google.colab
27 except ImportError:
---> 28 raise RuntimeError("This module must ONLY run as part of a Colab notebook!")
31 __version__ = "0.1.6"
32 __author__ = "Jaime Rodríguez-Guerra <jaimergp@users.noreply.github.com>"
RuntimeError: This module must ONLY run as part of a Colab notebook!
%%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()
<Topology; 1 chains, 100 residues, 1000 atoms, 900 bonds>
| number | name | type | atomic_number | charge | mass | nb_idx | solvent_radius | screen | occupancy | ... | rmin_14 | epsilon_14 | resname | resid | resnum | chain | segid | xx | xy | xz | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | -1 | CG | CG | 0 | 0.0 | 120.0 | 0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | CG0 | 0 | 1 | A | 0.0 | 0.0 | 0.0 | |
| 1 | -1 | CG | CG | 0 | 0.0 | 120.0 | 0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | CG0 | 0 | 1 | A | 0.0 | 0.0 | 3.8 | |
| 2 | -1 | CG | CG | 0 | 0.0 | 120.0 | 0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | CG0 | 0 | 1 | A | 0.0 | 0.0 | 7.6 | |
| 3 | -1 | CG | CG | 0 | 0.0 | 120.0 | 0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | CG0 | 0 | 1 | A | 0.0 | 0.0 | 11.4 | |
| 4 | -1 | CG | CG | 0 | 0.0 | 120.0 | 0 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | CG0 | 0 | 1 | A | 0.0 | 0.0 | 15.2 |
5 rows × 27 columns
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)
1
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)
#"Step","Potential Energy (kJ/mole)","Temperature (K)","Box Volume (nm^3)","Speed (ns/day)"
1000,-2504.836181640625,191.5762888669726,1331.0,0
2000,-2815.38623046875,228.05667202079903,1331.0,1e+03
3000,-2720.246337890625,237.99379315778745,1331.0,968
4000,-2717.5224609375,246.43023451709598,1331.0,942
5000,-2659.786865234375,253.65382630735783,1331.0,900
6000,-2649.3740234375,262.1551793714623,1331.0,877
7000,-2439.745361328125,257.93077347642696,1331.0,879
8000,-2385.37646484375,261.5250756318168,1331.0,894
9000,-2507.65869140625,277.1601449505948,1331.0,894
10000,-2277.37890625,263.72359415211446,1331.0,888
11000,-2352.224853515625,277.5869879541712,1224.1066608916633,881
12000,-2151.421875,272.3363540797704,1217.1759308796754,850
13000,-2169.27490234375,274.05078020934116,1126.4046429984535,849
14000,-2102.018310546875,272.95167374749195,1059.0785394736092,831
15000,-2159.517578125,280.71576208676345,1019.7380219560032,810
16000,-2180.33447265625,283.8484832712969,916.285214431259,796
17000,-2051.6181640625,272.3146143999463,904.4493001807455,781
18000,-2094.2294921875,284.0290534109619,899.8734092637802,770
19000,-2165.8671875,283.0528967485372,851.9844830974769,760
20000,-2205.318115234375,281.1227509728958,819.3514649819555,750
21000,-2240.06640625,277.7291262216464,744.7108329580016,743
22000,-2270.175048828125,289.7231251012086,690.54029952904,735
23000,-2166.90283203125,288.02704420530813,714.4953652595874,730
24000,-2042.0665283203125,282.2783723471896,657.17060647423,724
25000,-2157.97265625,289.81985838884697,568.0055136847603,726
26000,-2064.509033203125,279.2794794212596,528.1222343854555,720
27000,-2097.9609375,278.8585177835825,476.022732738361,715
28000,-2136.488525390625,285.4858725279452,454.2461928375308,717
29000,-2275.820068359375,297.4297078624183,411.57515560451657,717
30000,-2302.36328125,296.8539829212256,382.9434600057386,713
31000,-2322.546630859375,293.28292644135826,351.7996888497053,708
32000,-2323.501708984375,301.54835140926133,326.3492824758923,708
33000,-2344.83203125,296.83241114510923,305.0296352469343,704
34000,-2371.13818359375,292.0429706275744,272.28412890464665,700
35000,-2395.903564453125,299.9779724446933,255.96621936747422,697
36000,-2426.533203125,300.50611639627925,245.95579266983378,694
37000,-2463.66845703125,301.7052921624469,231.4306650768564,690
38000,-2523.1689453125,303.3020174178514,222.77053028028178,687
39000,-2468.69091796875,305.6769871541834,215.18216250964574,683
40000,-2439.684814453125,296.2453221377333,202.64792525833616,683
41000,-2619.2255859375,304.21391337578046,182.74158676803813,684
42000,-2653.97216796875,303.5549828120358,170.45828251946833,682
43000,-2828.11474609375,322.2711930312677,164.31983469702786,683
44000,-2836.41455078125,313.82151050042694,162.49357145363723,681
45000,-2784.0439453125,317.83260171624175,161.80449110324318,679
46000,-2966.98291015625,326.5362778178724,154.9536375508575,680
47000,-2886.562744140625,317.998989801577,155.0226264149907,682
48000,-3015.0478515625,311.83503874485706,157.29521573680645,683
49000,-3177.96435546875,321.1500196360962,155.53068279557974,684
50000,-3078.247802734375,304.8727344661126,153.33806897954847,684
51000,-3134.77392578125,316.3601262254108,154.3395998209216,683
52000,-3081.000732421875,310.44632688484296,155.36769816752778,684
53000,-3122.693359375,302.79666128507097,153.3511198045237,682
54000,-3178.40673828125,309.7085792320148,152.94065454850175,681
55000,-3078.410888671875,306.5856397639294,153.7837242107676,679
56000,-3115.92236328125,313.06766546869636,151.96553426237358,678
57000,-3112.86669921875,307.88280624752576,151.66837533343082,678
58000,-3350.757080078125,317.4928445252218,149.93140945617438,676
59000,-3119.8427734375,301.5432656423632,153.63852149589025,674
60000,-3099.111328125,304.16718443278126,157.00929518595655,670
61000,-3029.989013671875,305.5084019305884,156.908485836683,670
62000,-3159.145263671875,310.9742098333687,153.90360261610758,670
63000,-3122.4853515625,313.8566009005241,154.16403424292062,669
64000,-3102.88623046875,318.05293753609993,152.6325878428767,670
65000,-3174.193115234375,313.31148772433824,153.2000406220117,671
66000,-3150.25146484375,307.8944321803099,156.35209583434252,672
67000,-3111.649169921875,310.2161733009433,153.3931228026114,674
68000,-3064.90869140625,308.1185115335678,154.47906490680788,674
69000,-3107.314453125,314.2226441187978,158.47411143390255,674
70000,-3020.36669921875,312.8071934702927,157.91053181231752,674
71000,-2915.658935546875,314.3003198591404,158.6898914328477,676
72000,-2982.724609375,321.3533961423947,157.8584424606822,677
73000,-2993.283935546875,321.6774788462723,155.41351035702309,678
74000,-2958.010986328125,318.93175837874475,157.91788898881128,678
75000,-3051.702880859375,318.42951885216286,157.48571267645266,679
76000,-2923.89306640625,305.11767557968307,156.4514785054374,680
77000,-2974.4990234375,306.5367243491593,160.3551099114306,679
78000,-3047.22021484375,306.0771864505819,157.94235644004854,678
79000,-3165.8349609375,305.71471397570383,152.0875388951903,677
80000,-3296.365478515625,309.95971729495164,148.9298681726276,677
81000,-3283.8798828125,309.73738877238293,151.35860837380065,676
82000,-3100.02490234375,292.8003349870861,157.51971469678924,674
83000,-3257.4453125,304.68531198358164,153.81058752207832,673
84000,-3231.6474609375,310.0404361112034,150.90510791929688,672
85000,-3145.6240234375,307.5033965192968,153.2782580280336,671
86000,-3214.52294921875,319.09915644248554,151.83007815204505,669
87000,-3218.677734375,317.28295016042006,150.82525814455423,668
88000,-3302.59228515625,319.6713106386491,151.79410226860819,667
89000,-3147.118408203125,310.6033133316036,152.1655106684356,666
90000,-3246.8740234375,310.14530159966756,151.2346444030287,667
91000,-3209.808837890625,306.74162854886487,150.17269011389885,668
92000,-3243.70068359375,306.5829723467821,146.6921216035493,668
93000,-3276.2587890625,307.24945198497886,147.8624787966208,669
94000,-3287.388916015625,317.2650315294048,148.32022445617525,669
95000,-3130.58935546875,302.3975100257787,151.06547142013616,668
96000,-3104.04248046875,300.41556443930267,154.17320423925096,668
97000,-3246.49267578125,306.18753227541663,147.45068295808045,669
98000,-3310.83642578125,317.5655163341493,150.5140459515414,669
99000,-3089.177734375,298.15170496697965,156.55254032343626,670
100000,-3332.8876953125,313.1646673624205,153.56280187500843,670
101000,-3228.376708984375,299.9204040863644,156.1136791113453,670
102000,-3369.59033203125,308.1169559210054,149.28130652753856,670
103000,-3327.7177734375,288.62873353467984,148.83696334740617,671
104000,-3511.173828125,296.1635038932168,146.49862081363338,671
105000,-3465.52587890625,297.8163413025938,145.32428207096052,670
106000,-3553.03857421875,312.52290924248274,144.98710882745706,670
107000,-3371.2021484375,296.193638173217,150.1696419692844,670
108000,-3450.92626953125,294.6664424358834,147.3674255651175,670
109000,-3472.06787109375,295.36016825585796,144.40580488044853,671
110000,-3406.15771484375,295.70452258307296,145.63809854730596,671
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
/var/folders/3h/w4n54dz57dvbhgvz68gdv0f0f3hkmm/T/ipykernel_64115/607502326.py in <cell line: 26>()
26 with open('equilibrated_config.pdb','w') as f:
27 state = simulation.context.getState(getPositions=True, enforcePeriodicBox=True)
---> 28 topology.setPeriodicBoxVectors(state.getPeriodicBoxVectors())
29 app.PDBFile.writeFile(topology, state.getPositions(), f)
NameError: name 'topology' is not defined
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)