Alexei Abrikossov and Mikael Lund, December 2016
In this Notebook we setup both MC and MD simulation to calculate the interaction free energy between a pair of CPPM's. The layout is as follows:
In [1]:
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate
import pandas as pd
import os.path, os, sys, json, shutil
from pathlib import Path
import pickle
from sys import stdout
from math import exp, sqrt
import mdtraj as md
import simtk.unit
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
plt.rcParams.update({'font.size': 16, 'figure.figsize': [6.0, 5.0]})
try:
workdir
except NameError:
workdir=%pwd
else:
%cd $workdir
In [18]:
rho1 = 9.0/nanometer**2
(1/rho1).in_units_of(angstrom**2)
Out[18]:
In [25]:
r=15*angstrom
a=4*np.pi*r**2
rho2 = a/10
(1/rho2).in_units_of(nanometer**-2)
Out[25]:
In this section we use OpenMM to run a steered MD simulation where the two CPPM's are pulled towards each other using an external, harmonic potential. Done slowly, this enables us to sample the equilibrium force as a function of separation. We will lastly integrate the mean force to obtain the potential of mean force, PMF.
Some details:
MP.
In [37]:
dr=0.1 # distance spacing in nm
param = pd.DataFrame({
'Neutral' :
dict(
L=15*simtk.unit.nanometer, sep=4*simtk.unit.nanometer,
R=np.arange(4.2, 5.0, dr)*simtk.unit.nanometer, label=r'Neutral',
charge=0,
ions = [ ]
),
'P00' :
dict(
L=15*simtk.unit.nanometer, sep=4*simtk.unit.nanometer,
R=np.arange(4.2, 5.0, dr)*simtk.unit.nanometer, label=r'$P_0^0$',
charge=-8,
ions = [ dict(Na=100, Cl=100), dict(La=100, Cl=300) ]
),
'P18' :
dict(
L=15*simtk.unit.nanometer, sep=4*simtk.unit.nanometer,
R=np.arange(4.2, 5.0, dr)*simtk.unit.nanometer, label=r'$P_8^1$',
charge=-8,
ions = [ dict(Na=100, Cl=100), dict(La=100, Cl=300) ]
)
})
param.T
Out[37]:
In [52]:
# Generate initial system with two CPPMs, one locked in the middle of the box, the other
# some distance away from it. The final file is called "system.gro".
%cd -q $workdir
def makeInitialConfiguration(pdbin, sep, L, ionlist={}):
''' Generate initial .gro file
Two molecules of `pdbin` will be inserter, one in the box center, the
other a distance `separation` from it, along x. If an ion list is given, these
will be inserted such that they do not overlap with the two molecules.
For each ion, a single line PDB with name `ionname.pdb` must exist.
Arguments:
pdbin - input structure. Two will be inserted, one in the box center (str)
sep - COM-COM separation (float)
L - Cube side length (float)
ionlist - List of ions to insert, for example `dict(Na=10, Cl=10)`
'''
L = L.value_in_unit(simtk.unit.angstrom)
sep = sep.value_in_unit(simtk.unit.angstrom)
!echo "ATOM 1 Na NA X 1 149.529 134.999 136.634 0.00 0.00" > Na.pdb
!echo "ATOM 1 Cl CL X 1 149.529 134.999 136.634 0.00 0.00" > Cl.pdb
!echo "ATOM 1 La LA X 1 149.529 134.999 136.634 0.00 0.00" > La.pdb
topol = '''
#include "amber99.ff/forcefield.itp"
[ atomtypes ]
; name1 name2 mass charge ptype sigma epsilon
UP UP 1 0.0000 A 0.3 0.24746
NP NP 1 0.0000 A 0.3 0.24746
MP MP 1 0.0000 A 0.3 0.24746
PP PP 1 0.0000 A 0.3 0.24746
Na Na 1 1.0000 A 0.3 0.24746
Cl Cl 1 -1.0000 A 0.3 0.24746
La La 1 +3.0000 A 0.3 0.24746
;M1 M1 1 1.0000 A 0.3 0.24746
;D2 D2 1 2.0000 A 0.3 0.24746
;T3 T3 1 3.0000 A 0.3 0.24746
;AM1 AM1 1 -1.0000 A 0.3 0.24746
;AD2 AD2 1 -2.0000 A 0.3 0.24746
;AT3 AT3 1 -3.0000 A 0.3 0.24746
[ moleculetype ]
NA 1
[ atoms ]
1 Na 1 NA NA 1 1.00000
[ moleculetype ]
CL 1
[ atoms ]
1 Cl 1 CL CL 1 -1.00000
[ moleculetype ]
LA 1
[ atoms ]
1 La 1 LA LA 1 +3.00000
#include "sphere.itp"
[ system ]
Patchy particles in water
[ molecules ]
SPH 1
SPH 1
'''
PACKMOL_INPUT = """
sidemax {a}
tolerance 2.0
filetype pdb
output system.pdb
structure {pdbin}
resnumbers 2
number 1
center
fixed 0 0 0 0. 0. 0.
end structure
structure {pdbin}
resnumbers 2
number 1
center
fixed {sep} 0 0 0. 0. 0.
end structure
""".format( a=0.5*L, sep=sep, pdbin=pdbin)
# add ions
r=20 # radius
for ion, N in ionlist.items():
s = '''
structure {ion}.pdb
number {N}
outside box {x1} {y1} {z1} {x2} {y2} {z2}
end structure'''.format(ion=ion, N=N, x1=-r, y1=-r, z1=-r, x2=sep+r, y2=r, z2=r)
PACKMOL_INPUT += s
topol += ion.upper()+' '+str(N)+'\n'
!echo '$PACKMOL_INPUT' > packmol.inp
!~/Downloads/packmol/packmol sidemax $L < packmol.inp > /dev/null
with open('topol.top', 'w') as f:
f.write( topol )
cryst = "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1 1" % (L,L,L,90,90,90)
!sed '5s/.*/$cryst/' system.pdb > tmp.pdb ; mv tmp.pdb system.pdb # add box info
md.load('system.pdb').save_gro('system.gro') # convert pdb -> gro
#!rm -fR system.pdb
#L=30*simtk.unit.nanometer
#makeInitialConfiguration('sphere.pdb', L=L, sep=5*simtk.unit.nanometer, ionlist=dict(Na=1000, Cl=1000) )
In [53]:
%cd -q $workdir
for particletype, d in param.items(): # 'P00', 'P18', ...
if os.path.isdir( particletype ):
%cd $particletype
md.load('sphere.gro').save_pdb('sphere.pdb') # gro -> pdb
for ionlist in d.ions: # loop over all salt conditions
iondir = ''
for ion, N in ionlist.items():
iondir = ion+str(N) + iondir
print(iondir)
if not os.path.isdir(iondir):
%mkdir $iondir
if os.path.isdir(iondir):
%cd -q $iondir
makeInitialConfiguration('../sphere.pdb', L=d.L, sep=d.R[-1], ionlist=ionlist)
%cd -q ..
%cd -q ..
In [ ]:
def scaleCharges(system, scale=1.0):
"""
Scale charges in all appropriate force classes of system. This can
be used to simulate a system with a different background dielectric
by setting scale=1/sqrt(epsilon_r), for example.
"""
qsum=0
for force in system.getForces():
if 'getParticleParameters' in dir(force):
for i in range(force.getNumParticles()):
charge, sigma, epsilon = force.getParticleParameters(i)
force.setParticleParameters(i, scale*charge, sigma, epsilon)
qsum += charge / elementary_charge
if 'getExceptionParameters' in dir(force):
for i in range(force.getNumExceptions()):
p1, p2, chargeProd, sigma, epsilon = force.getExceptionParameters(i)
force.setExceptionParameters(i, p1, p2, chargeProd*(scale**2), sigma, epsilon)
print('Net charge before scaling =', qsum)
class MeanForceReporter(object):
''' Calculate mean force between two particle ranges
This will calculate the mean force on particle index 1 and particle
index 2 along the COM-COM distance between the two particle sets.
Only forces in given OpenMM force is used to calculate
the forces, hence the manual call to `getState()`. Periodic boundaries
and minimum image distances are *not* considered, so make sure molecules
are not wrapped, nor can the COM-COM distance be longer than half the
box length.
Keyword arguments:
cm1 -- index of particle situated on mass center 1 (int)
cm2 -- index of particle situated on mass center 1 (int)
index1 -- particle index of group1 (list)
index2 -- particle index of group1 (list)
groups -- force groups to extract forces from (set)
reportInterval -- Steps between each sample (int)
'''
def __init__(self, cm1, cm2, index1, index2, groups, reportInterval):
self._reportInterval = reportInterval
self._cm1 = cm1
self._cm2 = cm2
self._index1 = index1
self._index2 = index2
self.mf1 = 0 # mean force on cm1 (kJ/mol/nm)
self.mf2 = 0 # mean force on cm2 (kj/mol/nm)
self.mR = 0 # mean COM-COM distance
self.cnt = 0 # number of samples
self.groups = groups
def describeNextReport(self, simulation):
steps = self._reportInterval - simulation.currentStep%self._reportInterval
return (steps, False, False, False, False)
def clear(self):
'''Resets the mean force average'''
self.mf1=0
self.mf2=0
self.mR=0
self.cnt=0
def meanforce(self):
'''Returns average meanforce along COM-COM vector (kJ/mol/nm)'''
if self.cnt>0:
return 0.5*(self.mf1+self.mf2) / self.cnt, self.mf1/self.cnt, self.mf2/self.cnt, self.mR/self.cnt
else:
print('no force samples')
def report(self, simulation, state):
s = simulation.context.getState(getPositions=True, getForces=True, groups={0}) # system
f = s.getForces(asNumpy=True).value_in_unit(kilojoules/mole/nanometer) # forces
r = s.getPositions(asNumpy=True).value_in_unit(nanometer) # distances
R = r[ self._cm1 ] - r[ self._cm2 ] # COM-COM distance vector...
Rhat = R / np.linalg.norm(R) # ...and scalar
self.mR += np.linalg.norm(R) # mean COM-COM distance
self.mf1 += np.inner( f[ self._index1 ], Rhat).sum() # Mean force 1 along R
self.mf2 += np.inner( f[ self._index2 ], -Rhat).sum() # Mean force 2 along -R
self.cnt += 1 # number of samples since last call to clear()
gro = GromacsGroFile('system.gro')
top = GromacsTopFile('system.top', periodicBoxVectors=gro.getPeriodicBoxVectors(),
includeDir='/opt/local/share/gromacs/top/')
system = top.createSystem( nonbondedCutoff=1.2*nanometer, nonbondedMethod=PME)
integrator = LangevinIntegrator( 300*kelvin, 1/picosecond, 0.002*picoseconds)
platform = Platform.getPlatformByName( 'OpenCL' )
properties = {}#'OpenCLDeviceIndex': '0,1', 'OpenCLPrecision': 'mixed'}
# OpenMM does not allow to set the background dielectric. Instead we scale charges:
scaleCharges(system, scale=1/sqrt(78.44))
# modify automatically added forces
for i in range(0, system.getNumForces()):
f = system.getForce(i)
if (type(f) == CMMotionRemover): # remove COM motion remover
system.removeForce(i)
if type(f) == NonbondedForce: # remove dispersion correction
f.setUseDispersionCorrection(False)
f.setEwaldErrorTolerance(0.05) # default=0.0005
if type(f) == HarmonicBondForce: # move constraints to another
f.setForceGroup(1) # force group
# add harmonic bond between COMs
t = md.load('system.gro') # mdtraj trajectory
cm = t.top.select('name MP').tolist() # particle index of the two COMs
harmonic = HarmonicBondForce()
harmonic.addBond( cm[0], cm[1], 5.0*nanometer, 10000*kilojoule_per_mole/nanometer**2 )
harmonic.setForceGroup(2)
system.addForce( harmonic )
# fix COM of molecule 1 to middle of box
freezeforce = CustomExternalForce('8000*r^2; r=sqrt((x-L/2)^2 + (y-L/2)^2 + (z-L/2)^2); L=%f' % L.value_in_unit(nanometer) )
freezeforce.addParticle(cm[0], [])
freezeforce.setForceGroup(2)
system.addForce(freezeforce)
simulation = Simulation(top.topology, system, integrator, platform)
simulation.context.setPositions(gro.positions)
nb = system.getForce(0)
print('PME parameters: ', nb.getPMEParametersInContext(simulation.context))
for f in system.getForces():
print(type(f), 'force group =',f.getForceGroup())
In [ ]:
rho = 1000./(30*nanometer)**3 / (6.022e23/mole)
rho.in_units_of(millimolar)
In [ ]:
# (de)serialize to/from XML files. Works for system, force, integrator, state, platform
with open("system.xml", "w") as f:
f.write( XmlSerializer.serialize(system) )
with open('system.xml', 'r') as f:
system = XmlSerializer.deserialize( f.read() )
In [ ]:
print('Minimization...')
simulation.minimizeEnergy() # minization
pos = simulation.context.getState( getPositions=True ).getPositions()
PDBFile.writeFile(simulation.topology, pos, open( 'minimized.pdb', 'w'))
simulation.step(1000) # pre-equilibration
res1 = t.top.select('resid 0') # index of molecule 1
res2 = t.top.select('resid 1') # index of molecule 2
simulation.reporters.clear()
meanforce = MeanForceReporter( cm[0], cm[1], res1, res2, {0}, 50 )
simulation.reporters.append( meanforce )
#simulation.reporters.append( DCDReporter('output.dcd', 100))
#simulation.reporters.append( StateDataReporter(stdout, 100, step=True, temperature=True) )
print('Distance scan...')
rinterval=[4.0, 4.6] # distance range (nm)
dr=0.1 # distance steps (nm)
nsteps=1000 # steps per distance, 2 fs each
_r = []
_f = []
if True:
for req in np.arange(rinterval[1], rinterval[0], -dr):
# update spring eq. distance in context
harmonic.setBondParameters(0, cm[0], cm[1], req*nanometers, 10000*kilojoule_per_mole/nanometer**2)
harmonic.updateParametersInContext(simulation.context)
simulation.step( 1000 ) # equilibration for new distance
meanforce.clear() # clear mean force
simulation.step( nsteps ) # propagate and sample mean force
_r.append( req )
_f.append( meanforce.meanforce()[1] )
print(req, meanforce.meanforce())
print('done.')
In [ ]:
#rmc, g = np.loadtxt('mc/rdf.dat', unpack=True)
#rmc = 0.1*rmc
#w = -np.log( g / g[rmc>4.75].mean() )
pmf = integrate.cumtrapz( _f, _r, initial=0 ) # integrate force --> potential of mean force
plt.plot(_r, -(pmf/2.5), 'bo-', label='OpenMM (MD)')
#plt.plot(rmc, w, 'r-', label='Faunus (MC)', lw=2)
plt.legend(loc=0, frameon=False)
#plt.xlim(4.15, 4.6)
In [ ]:
t, r = np.loadtxt('pullx_old.xvg', skiprows=15, unpack=True) # distance (nm)
t, f = np.loadtxt('pullf_old.xvg', skiprows=15, unpack=True) # force (kJ/mol/nm)
sort = np.argsort(r) # sort w. respect to distance
r = r[sort]
f = f[sort]
t = t[sort]
In [ ]:
pmf = integrate.cumtrapz( f, r, initial=0 ) # integrate force
pmf = pmf - pmf[ (r>10) & (r<11)].mean() # shift pmf to zero at long separations
pmf = pmf / 2.48 # kJ/mol --> kT
plt.plot( r, pmf )
plt.xlabel('$R$ (nm)')
plt.ylabel('PMF ($k_BT$)')
plt.xlim(4.1, 5)
plt.ylim(-5,4)
print('PMF minimum =', pmf[r<6].min(), 'kT')
In [ ]:
%%bash -s "$workdir"
cd $1
if [ ! -d "faunus" ]; then
git clone https://github.com/mlund/faunus.git
cd faunus
#git checkout dfdddf3
else
cd faunus
fi
# if different, copy custom .cpp file into faunus
if ! cmp ../mc/twobody.cpp src/examples/twobody.cpp >/dev/null 2>&1
then
cp ../mc/twobody.cpp src/examples/
fi
cmake . -DCMAKE_BUILD_TYPE=RelWithDebInfo -DENABLE_APPROXMATH=on &>/dev/null
make example_twobody -j
cd $1
In [ ]:
%cd $workdir/mc
def mkinput():
d = {
"energy" : {
"nonbonded" : {
"coulomb" : { "epsr" : 78.7, "ionicstrength" : 0.1 }
},
"cmconstrain" : {
"sphere1 sphere2" : { "mindist": 0, "maxdist": 50 }
}
},
"atomlist" : {
"UP": dict(q=0, sigma=3.0, eps=0.2479, mw=1e-3),
"NP": dict(q=0, sigma=3.0, eps=0.2479, mw=1e-3),
"PP": dict(q=0, sigma=3.0, eps=0.2479, mw=1e-3),
"MP": dict(q=0, sigma=3.0, eps=0.2479, mw=1e6)
},
"moleculelist": {
"sphere1": { "structure":"sphere.xyz", "Ninit":1, "insdir":"0 0 1", "insoffset":"0 0 -20"},
"sphere2": { "structure":"sphere.xyz", "Ninit":1, "insdir":"0 0 1", "insoffset":"0 0 20"}
},
"moves" : {
"moltransrot2body" : {
"sphere1" : { "dp":30, "dprot":1, "prob":1.0 },
"sphere2" : { "dp":30, "dprot":1, "prob":1.0 }
}
},
"analysis" : {
"pqrfile" : { "file": "confout.pqr" },
"statefile" : { "file": "state" }
},
"system" : {
"temperature" : 298.15,
"cylinder" : { "length" : 400, "radius" : 200 },
"mcloop" : { "macro" : 10, "micro" : 5000 }
}
}
with open('twobody.json', 'w+') as f:
f.write(json.dumps(d, indent=4))
mkinput()
#!rm -fR state
!../faunus/src/examples/twobody
In [ ]:
%cd $workdir/mc
rmc, g = np.loadtxt('rdf.dat', unpack=True)
rmc = 0.1*rmc
w = -np.log( g / g[rmc>4.8].mean() )
plt.plot(rmc, w, label='faunus' )
plt.xlabel('$R$ (nm)')
plt.ylabel('PMF ($k_BT$)')
plt.xlim(4.1, 5)
plt.ylim(-5.2,4)
plt.legend(loc=0)
print('PMF minimum =', w[rmc<6].min(), 'kT')