This will simulate:
This Jupyter Notebook was originally run in MacOS 10.11 with GCC 4.8, Python2, matplotlib, pandas within the Anaconda environment. Contemporary Linux distributions such as Ubuntu 14.04 should work as well.
We use a custom Metropolis Monte Carlo (MC) program build within the Faunus framework. The sections below will fetch the complete faunus project and compile the program.
In [ ]:
from __future__ import division, unicode_literals, print_function
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np, pandas as pd
import os.path, os, sys, json
plt.rcParams.update({'font.size': 16, 'figure.figsize': [8.0, 6.0]})
try:
workdir
except NameError:
workdir=%pwd
else:
%cd $workdir
In [ ]:
%%bash -s "$workdir"
cd $1
if [ ! -d "faunus/" ]; then
echo 'fau_example(surf "./" surf.cpp)' > mc/CMakeLists.txt
git clone https://github.com/mlund/faunus.git
cd faunus
git checkout 4da3195977b44012bfacbcd63ec11f0f82e3b836
else
cd faunus
fi
pwd
CXX=c++ CC=cc cmake . -DCMAKE_BUILD_TYPE=Release -DENABLE_APPROXMATH=on -DMYPLAYGROUND=`pwd`/../mc &>/dev/null
make surf -j4
In [ ]:
def mkinput():
js = {
# customize LJ interaction between pairs (Ångstrom & kJ/mol)
"customlj" : {
"Na COOH" : { "sigma":4.0, "eps":0.7 },
"Na COO" : { "sigma":4.0, "eps":0.7 }
},
# initial degree of protonation of the lipids
"initial_degree_of_protonation" : 0.5,
# to turn of lipid titration, just rename to i.e. "_processes"
"processes": {
"acid": { "pKd": 4.75, "pX": pH, "bound": "COOH", "free": "COO" }
},
"energy": {
"eqstate": { "processfile": "gctit.json" },
"nonbonded": {
"coulomb": { "epsr": 80 },
"ljsimple": { "eps": 0.05 }
}
},
"system": {
"temperature": 298.15,
"cuboid": { "xyzlen": str(xylen)+' '+str(xylen)+' '+str(zlen) },
"mcloop": { "macro": 10, "micro": micro }
},
"moves": {
"gctit" : { "molecule": "salt", "prob": 0.5 },
"atomtranslate" : {
"lipid": { "prob": 1.0, "dir":"1 1 0" },
"salt": { "prob": 0.5 }
}
},
"moleculelist": {
"lipid": { "Ninit": Nlipid, "atomic": True, "atoms": "COO", "insdir": "1 1 0",
"insoffset":"0 0 "+str(surfpos) },
"salt": { "Ninit": 50, "atomic": True, "atoms": "Na Cl" }
},
"atomlist": {
"Na": { "q": 1, "r": 2, "dp": 50, "eps":0.05, "activity": salt },
"Li": { "q": 1, "r": 2, "dp": 50, "eps":0.05, "activity": 0.0 },
"Cl": { "q": -1, "r": 2, "dp": 50, "eps":0.05, "activity": salt },
"COOH": { "q": 0, "r": 2, "dp": 10, "eps":0.05 },
"COO": { "q": -1, "r": 2, "dp": 10, "eps":0.05 }
}
}
with open('surf.json', 'w+') as f:
f.write(json.dumps(js, indent=4))
%cd $workdir'/mc'
micro=250000
xylen=40
zlen=60
surfpos=-0.5*zlen+1e-6
for pH in [7.0]:
for salt in [1.0]: # mol/l
for rho in [25]: # angstrom^2 per lipid
Nlipid = int(xylen**2/rho)
print('Number of head groups = ', Nlipid)
pfx='rho'+str(rho)+'-pH'+str(pH)+'-I'+str(salt)
if True: #not os.path.isdir(pfx):
%mkdir -p $pfx
%cd $pfx
!rm -fR state
mkinput()
!../surf
%cd ..
In [ ]:
%cd $workdir'/mc'
for pH in [7.0]:
for salt in [1.0]: # mol/l
for rho in [25]: # angstrom^2 per lipid
pfx='rho'+str(rho)+'-pH'+str(pH)+'-I'+str(salt)
if os.path.isdir(pfx):
%cd $pfx
r1, P1 = np.loadtxt("zhist_cation.dat", unpack=True)
r2, P2 = np.loadtxt("zhist_anion.dat", unpack=True)
P1 = P1 / P1[ r1>15 ].mean() # normalize to unity at long sep.
P2 = P2 / P2[ r2>15 ].mean() # normalize to unity at long sep.
plt.plot(r1-r1.min(), P1, 'r-', label='cation')
plt.plot(r2-r1.min(), P2, 'b-', label='anion')
%cd ..
plt.legend(loc=0, frameon=False)
plt.xlabel(r'$z$-distance from lipids (Å)')
plt.ylabel(r'Relative probability')
plt.ylim( (0, 20) )
plt.xlim( (0, 20) )
plt.savefig('fig.pdf', bbox_inches='tight')
In [ ]: