Simulating toy polymers using openMM

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

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

import openmm as mm
from   openmm import app
from   openmm.unit import * 

import parmed
import mdtraj as md
import py3Dmol
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Simulation parameters
temperature = 293.15 * kelvin
pressure = 1 * bar
mass = 39.948 * amu
sig = 3.419 * angstrom
eps = 117.8 * kelvin * BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA

box_size = 150 * angstrom  
natom    = 12
cutoff   = 3 * sig

Create polymer structure and topology

• Here we will use parmed which is. apowerful tool for creating and manipulating structures and topologies of complex molecules.

from parmed import Structure, Atom, Bond, Angle, Dihedral

s  = Structure()

for i in range(natom):
  s.add_atom(atom     =  Atom(name='A',  mass=mass), 
             resname  = "LJ", 
             resnum   =  i, 
             chain    = 'A')
  
s.positions = [(0, 0, z) for z in range(natom)]

for i in range(natom-1):
  s.bonds.append(Bond(s.atoms[i], s.atoms[i+1]))

s.save('lj.pdb', overwrite=True)
s.save('lj.psf', overwrite=True)

print(s.topology)
s.to_dataframe().head()

Create openmm systema and add particles

# Add particles to the system
system = mm.System()

for i in range(natom):
  system.addParticle(s.atoms[i].mass)

box_vecs = box_size  *np.eye(3)
system.setDefaultPeriodicBoxVectors(*box_vecs )

### Create LJ force
E_lj  = '4*eps*((sig/r)^12-(sig/r)^6); sig=0.5*(sig1+sig2); eps=sqrt(eps1*eps2)'
force = mm.CustomNonbondedForce(E_lj)
force.addPerParticleParameter('sig') 
force.addPerParticleParameter('eps')

# Particles are assigned properties in the same order as they appear in the System object
for _ in range(natom): 
    force.addParticle([sig, eps])

# Set force cutoff parameters
force.setNonbondedMethod(mm.NonbondedForce.CutoffPeriodic)
force.setCutoffDistance(3.0 * sig)       # set cutoff (truncation) distance at 3*sigma
force.createExclusionsFromBonds([(i, i+1) for i in range(natom-1)], 1) # exlcude bonded particles from LJ

### Create Harmonic force
force2 = mm.HarmonicBondForce()

for i in range(natom-1):
  force2.addBond(i, i+1, 1.5*sig, 100) 

### Add a force to remove Center of Mass motion to prevent drifting of molecule
force3 = mm.CMMotionRemover()

# Added forces to system
system.addForce(force)
system.addForce(force2)
system.addForce(force3)

print('Check if force uses PBC: ', force.usesPeriodicBoundaryConditions() )

print('No particles: ', force.getNumParticles() )

print('per-particle parameters for particle-0: ', force.getParticleParameters(0))

print('Check Energy function: ', force.getEnergyFunction() )

Create OpenMM simulation object

integrator =  mm.LangevinIntegrator(temperature, 1/picosecond, 2*femtoseconds) 

simulation = app.Simulation(s.topology, 
                            system, 
                            integrator) 

simulation.context.setPositions(s.positions)

# - Minimize the energy
simulation.minimizeEnergy()

# - Initialize velocities with random values at 300K.
simulation.context.setVelocitiesToTemperature(300)


# Reporters
simulation.reporters = []
simulation.reporters.append(app.DCDReporter('ljtraj.dcd', 100))
simulation.reporters.append(app.PDBReporter('ljtraj.pdb', 100))
simulation.reporters.append(app.StateDataReporter("ljscalars.csv", 10, 
                                              time=True, 
                                              potentialEnergy=True, 
                                              totalEnergy=True, 
                                              temperature=True, 
                                              volume=True))

simulation.step(50000)

# 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()