Here we define a custom force class where particles interact through a Yukawa potential and a soft repulsion,
\begin{equation} w(r) / k_BT = \frac{\lambda_Bz_iz_j}{r}e^{-r/\lambda_D} + 4\beta\epsilon_{ij} \left ( \frac{\sigma_{ij}}{r}\right )^{12} \end{equation}where $\lambda_B=e^2/4\pi\epsilon_0\epsilon_rk_BT$ and $\lambda_D=(4\pi\lambda_B\sum \rho_iz_i^2)^{-1/2}$ are the Bjerrum and Debye lengths, respectively. $\rho_i$ is the number density of the $i$th ion. In this example we also create two rigid bodies using harmonic bonds to constrain the positions.
Some comments:
CustomNonbonded is defined in cg.zml and must return energy in kJ/mol.
In [7]:
%matplotlib inline
import numpy as np
from __future__ import print_function
from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
from sys import stdout, exit
import math
import mdtraj as mdtraj
from itertools import combinations
In [8]:
cutoff = 50*unit.angstrom
useMinimize = True
epsilon_r = 80.
temperature = 300*unit.kelvin
kT = unit.BOLTZMANN_CONSTANT_kB*temperature
timestep = 10*unit.femtoseconds;
steps_eq = 5000
steps_production = 2e4
steps_total = steps_eq + steps_production
In [9]:
def findForce(system, forcetype, add=True):
""" Finds a specific force in the system force list - added if not found."""
for force in system.getForces():
if isinstance(force, forcetype):
return force
if add==True:
system.addForce(forcetype())
return findForce(system, forcetype)
return None
def setGlobalForceParameter(force, key, value):
for i in range(force.getNumGlobalParameters()):
if force.getGlobalParameterName(i)==key:
print('setting force parameter', key, '=', value)
force.setGlobalParameterDefaultValue(i, value);
def atomIndexInResidue(residue):
""" list of atom index in residue """
index=[]
for a in list(residue.atoms()):
index.append(a.index)
return index
def getResiduePositions(residue, positions):
""" Returns array w. atomic positions of residue """
ndx = atomIndexInResidue(residue)
return np.array(positions)[ndx]
def uniquePairs(index):
""" list of unique, internal pairs """
return list(combinations( range(index[0],index[-1]+1),2 ) )
def addHarmonicConstraint(harmonicforce, pairlist, positions, threshold, k):
""" add harmonic bonds between pairs if distance is smaller than threshold """
print('Constraint force constant =', k)
for i,j in pairlist:
distance = unit.norm( positions[i]-positions[j] )
if distance<threshold:
harmonicforce.addBond( i,j,
distance.value_in_unit(unit.nanometer),
k.value_in_unit( unit.kilojoule/unit.nanometer**2/unit.mole ))
print("added harmonic bond between", i, j, 'with distance',distance)
def addExclusions(nonbondedforce, pairlist):
""" add nonbonded exclusions between pairs """
for i,j in pairlist:
nonbondedforce.addExclusion(i,j)
def rigidifyResidue(residue, harmonicforce, positions, nonbondedforce=None,
threshold=6.0*unit.angstrom, k=2500*unit.kilojoule/unit.nanometer**2/unit.mole):
""" make residue rigid by adding constraints and nonbonded exclusions """
index = atomIndexInResidue(residue)
pairlist = uniquePairs(index)
addHarmonicConstraint(harmonic, pairlist, pdb.positions, threshold, k)
if nonbondedforce is not None:
for i,j in pairlist:
print('added nonbonded exclusion between', i, j)
nonbonded.addExclusion(i,j)
def centerOfMass(positions, box):
""" Calculates the geometric center taking into account periodic boundaries
More here: https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions
"""
theta=np.divide(positions, box).astype(np.float) * 2*np.pi
x1=np.array( [np.cos(theta[:,0]).mean(), np.cos(theta[:,1]).mean(), np.cos(theta[:,2]).mean()] )
x2=np.array( [np.sin(theta[:,0]).mean(), np.sin(theta[:,1]).mean(), np.sin(theta[:,2]).mean()] )
return box * (np.arctan2(-x1,-x2)+np.pi) / (2*np.pi)
In [10]:
pdb = app.PDBFile('squares.pdb')
forcefield = app.ForceField('yukawa.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffPeriodic, nonbondedCutoff=cutoff )
box = np.array(pdb.topology.getPeriodicBoxVectors()).diagonal()
harmonic = findForce(system, mm.HarmonicBondForce)
nonbonded = findForce(system, mm.CustomNonbondedForce)
setGlobalForceParameter(nonbonded, 'lB', 0.7*unit.nanometer)
setGlobalForceParameter(nonbonded, 'kappa', 0.0)
for residue in pdb.topology.residues():
p = getResiduePositions(residue, pdb.positions)
print(centerOfMass(p, box))
rigidifyResidue(residue, harmonicforce=harmonic, nonbondedforce=nonbonded, positions=pdb.positions)
integrator = mm.LangevinIntegrator(temperature, 1.0/unit.picoseconds, timestep)
integrator.setConstraintTolerance(0.0001)
In [11]:
simulation = app.Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
if useMinimize:
print('Minimizing...')
simulation.minimizeEnergy()
print('Equilibrating...')
simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
simulation.step(steps_eq)
simulation.reporters.append(mdtraj.reporters.HDF5Reporter('trajectory.h5', 100))
simulation.reporters.append(
app.StateDataReporter(stdout, int(steps_total/10), step=True,
potentialEnergy=True, temperature=True, progress=True, remainingTime=False,
speed=True, totalSteps=steps_total, volume=True, separator='\t'))
print('Production...')
simulation.step(steps_production)
print('Done!')
In [ ]: