Simulating Ethane#
!pip install -q condacolab
import condacolab
condacolab.install()
%%capture
!conda install -c conda-forge openmm mdtraj parmed
!pip install py3dmol
import parmed
import mdtraj as md
import py3Dmol
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import openmm as mm
from openmm import app
from openmm.unit import *
Structure + topology file#
ATOM 1 C1 ETH 1 -3.553 2.382 0.000 1.00 0.00 C
ATOM 2 H11 ETH 1 -3.940 1.922 0.912 1.00 0.00 H
ATOM 3 H12 ETH 1 -3.941 1.831 -0.859 1.00 0.00 H
ATOM 4 H13 ETH 1 -3.919 3.410 -0.053 1.00 0.00 H
ATOM 5 C2 ETH 1 -2.016 2.361 0.000 1.00 0.00 C
ATOM 6 H21 ETH 1 -1.649 1.333 0.053 1.00 0.00 H
ATOM 7 H22 ETH 1 -1.627 2.912 0.859 1.00 0.00 H
ATOM 8 H23 ETH 1 -1.629 2.821 -0.912 1.00 0.00 H
CONECT 1 3 4 5 2
CONECT 2 1
CONECT 3 1
CONECT 4 1
CONECT 5 8 1 6 7
CONECT 6 5
CONECT 7 5
CONECT 8 5
MASTER 0 0 0 0 0 0 0 0 8 0 8 0
END
Forcefield file#
<ForceField>
<AtomTypes>
<Type name="0" class="c3" element="C" mass="12.01078"/>
<Type name="1" class="hc" element="H" mass="1.007947"/>
</AtomTypes>
<Residues>
<Residue name="ETH">
<Atom name="C1" type="0"/>
<Atom name="H11" type="1"/>
<Atom name="H12" type="1"/>
<Atom name="H13" type="1"/>
<Atom name="C2" type="0"/>
<Atom name="H21" type="1"/>
<Atom name="H22" type="1"/>
<Atom name="H23" type="1"/>
<Bond atomName1="C1" atomName2="H11"/>
<Bond atomName1="C1" atomName2="H12"/>
<Bond atomName1="C1" atomName2="H13"/>
<Bond atomName1="C1" atomName2="C2"/>
<Bond atomName1="C2" atomName2="H21"/>
<Bond atomName1="C2" atomName2="H22"/>
<Bond atomName1="C2" atomName2="H23"/>
</Residue>
</Residues>
<HarmonicBondForce>
<Bond class1="c3" class2="c3" length="0.15380" k="1945727.36"/>
<Bond class1="c3" class2="hc" length="0.10970" k="3145687.56"/>
</HarmonicBondForce>
<HarmonicAngleForce>
<Angle class1="c3" class2="c3" class3="hc" angle="1.91637152" k="391.756288"/>
<Angle class1="hc" class2="c3" class3="hc" angle="1.87762521" k="326.01728"/>
</HarmonicAngleForce>
<PeriodicTorsionForce>
<Proper class1="hc" class2="c3" class3="c3" class4="hc" periodicity1="3" phase1="0.0" k1="0.50208"/>
</PeriodicTorsionForce>
<NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
<Atom type="0" charge="-0.094100" sigma="0.3397710" epsilon="0.4510352"/>
<Atom type="1" charge="0.031700" sigma="0.2600177" epsilon="0.0870272"/>
</NonbondedForce>
</ForceField>
Creating system from structure and forcefield files#
pdb = app.PDBFile('ethane.pdb')
forcefield = app.ForceField('ethane.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.NoCutoff, constraints=app.HBonds)
integrator = mm.LangevinIntegrator(298.15*kelvin, 5.0/picoseconds, 2.0*femtoseconds)
integrator.setConstraintTolerance(1e-5)
simulation = app.Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
print('Minimizing...')
st = simulation.context.getState(getPositions=True,getEnergy=True)
print(f"Potential energy before minimization is {st.getPotentialEnergy()}")
simulation.minimizeEnergy(maxIterations=100)
st = simulation.context.getState(getPositions=True,getEnergy=True)
print(f"Potential energy after minimization is {st.getPotentialEnergy()}")
print('Equilibration...')
from sys import stdout
simulation.reporters.append(app.StateDataReporter(stdout, 100, step=True,
potentialEnergy=True, temperature=True, separator=','))
simulation.context.setVelocitiesToTemperature(150.0*kelvin)
simulation.step(2500)
Production run#
import time as time
print('Running Production...')
# Begin timer
tinit=time.time()
# Clear simulation reporters
simulation.reporters.clear()
# Reinitialize simulation reporters. We do this because we want different information printed from the production run than the equilibration run.
# output basic simulation information below every 250000 steps - (which is equal to 2 fs(250,000) = 500,000 fs = 500 ps)
simulation.reporters.append(app.StateDataReporter(stdout,
250000,
step=True,
time=True,
potentialEnergy=True,
temperature=True,
speed=True,
separator=','))
# write out a trajectory (i.e., coordinates vs. time) to a DCD
# file every 100 steps - 0.2 ps
simulation.reporters.append(app.DCDReporter('ethane_sim.dcd', 100))
# run the simulation for 1.0x10^7 steps - 20 ns
simulation.step(10000000)
# End timer
tfinal=time.time()
print('Done!')
print('Time required for simulation:', tfinal-tinit, 'seconds')
Analysis#
traj = md.load('ethane_sim.dcd', top='ethane.pdb')
atoms, bonds = traj.topology.to_dataframe()
atoms
bond_indices = [0, 4] # atoms to define the bond length
bond_length = md.compute_distances(traj, [bond_indices])
bondcounts, binedges, otherstuff = plt.hist(bond_length, bins=120)
plt.title('C-C bond length histogram')
plt.xlabel('Bond length (nm)')
plt.ylabel('Counts')
plt.show()
phi_indices = [1, 0, 4, 5] # atoms to define the torsion angle
phi = md.compute_dihedrals(traj, [phi_indices])
phicounts, phi_binedges, otherstuff = plt.hist(phi, bins=90) # create a histogram with 90 bins
plt.title('H-C-C-H torsion angle')
plt.xlabel(r'$\phi$ (rad)')
plt.ylabel('Counts')
Potential of Mean Force (PMF)#
So far in our analysis, we have looked at the distribution of bond lengths and torsion angles for ethane. However, we can also use our simulations to calculate thermodynamics properties of our system. For example, we can use our calculated distributions along with Boltzmann’s constant to calculate the potential of mean force (pmf), or energy change associated with changes in the bond length or torsion angle.
The potential of mean force is defined by the expression
\[ W(x) = -k_{B}T*ln[p(x)] + C \]
where $\(p(x)\)$ is the probability, or the histogram we calculated previously.
For our torsion angle, we can calculate and plot the PMF:
kB = 8.31446/1000 # Boltzmann constant in kJ/mol
Temp = 298.15 # simulation temperature
phicounts[phicounts==0] = 0.1 # get rid of any bins with 0 counts/infinite energy
pmf = -kB*Temp*np.log(phicounts) # W(x) = -kT*ln[p(x)] = -kT*ln[n(x)] + C
pmf = pmf - np.min(pmf) # subtract off minimum value so that energies start from 0
phi_bincenters = (phi_binedges[1:] + phi_binedges[:-1])/2 # compute centers of histogram bins
plt.plot(phi_bincenters, pmf)
plt.title('H-C-C-H torsion pmf')
plt.xlabel(r'$\phi$ (rad)')
plt.ylabel('Relative free energy (kJ/mol)')