Often, in soft-matter physics and chemistry we are interested in systems where the particle numbers are not fixed. This is for instance the case in systems where chemical reactions occur or when particles can be exchanged with a reservoir. The case of changing particle numbers is best described by the grand canonical ensemble.
Canonical Monte-Carlo or molecular dynamics simulations are not suitable for this task, since they conserve the particle numbers. However, the Metropolis Monte Carlo method can be generalized to the Grand-Canonical ensemble in a straightforward manner, yielding the so-called Grand-Canonical Monte Carlo method (GCMC), which is introduced in this tutorial.
This tutorial builds on the tutorial Widom insertion, in which you measured the chemical potential of a monovalent salt solution by using a combination of MD and MC techniques. In this part you will make use of the results (i. e. the chemical potentials) of the previous tutorial to simulate the partitioning of salt ions between a reservoir and a polymer solution. Therefore we introduce and use the grand-canonical Monte Carlo method.
The case of changing particle numbers is best described by the grand canonical ensemble [1]. In the grand canonical ensemble, the volume $V$, the temperature $T$ and the chemical potentials $\mu_i$ of the exchangeable species $i$ rather than the particle numbers $N_i$ are fixed. The chemical potentials $\mu_i$ correspond to the free energy (Helmholtz free energy $F$, will be defined later) cost of adding a particle of type $i$ while keeping everything else fixed:
$$ \mu_i = \left(\frac{\partial F}{\partial N_i}\right)_{T, V, \left\{N_j\right\}_{j\neq i}}. $$The grand canonical ensemble is thus also referred to as the $(\mu, V, T)$-ensemble. In the grand canonical ensemble, the probability for the system to contain $N_i$ particles of type $i$, i.e. a total of $N=\sum_{i=1}^{N_\mathrm{s}}N_i$ ($N_\mathrm{s}$ is the number of different types) and to be in a specific microstate $\xi\in\Gamma_{N}$ is given by the probability distribution
$$ p^\mathrm{G}\left(\left\{N_i\right\},\xi\right) = \frac{\exp\left(-\beta \left(H(\xi)-\sum_{i=1}^{N_\mathrm{s}}\mu_iN_i\right)\right)}{Z^\mathrm{G}h^{3N}\prod_{i=1}^{N_\mathrm{s}} N_i!}. \quad \tag{1} $$Here, $H$ denotes the Hamiltonian of the system with $N_i$ particles of type $i$, $h$ is the planck constant, $\beta=1/k_\mathrm{B}T$ the inverse temperature and $Z^\mathrm{G}$ is the grand-canonical partition function which is given by
$$ Z^\mathrm{G}\left(\left\{\mu_i\right\},V,T\right) = \sum_{N_1=0}^{\infty}...\sum_{N_\mathrm{s}=0}^{\infty}Z\left(\left\{N_i\right\},V,T\right)\exp\left(\beta \sum_{i=1}^{N_\mathrm{s}}\mu_iN_i\right), $$with the canonical partition function
$$ Z\left(\left\{N_i\right\},V,T\right) = \frac{1}{h^{3N} \prod_{i=1}^{N_s}N_i!}\int_{\Gamma_N} \mathrm{d}\xi \ \exp(-\beta H(\xi)). $$The partition function is connected to a thermodynamic potential, called the grand potential (Landau free energy):
$$ \Omega\left(\left\{\mu_i\right\},V,T\right) = -k_\mathrm{B}T\ln\left(Z^\mathrm{G}\left(\left\{\mu_i\right\},V,T\right)\right). $$In a similiar way, the Helmholtz free energy is defined as
$$ F\left(\left\{N_i\right\},V,T\right) = -k_\mathrm{B}T\ln\left(Z\left(\left\{N_i\right\},V,T\right)\right). $$When considering a system that can exchange particles with a reservoir, the chemical potential $\mu_i$ for each exchangeable species $i$ has to be equal in the reservoir and in the system: \begin{equation} \mu_i^{\mathrm{sys}} = \mu_i^{\mathrm{res}}. \end{equation} This condition is referred to as a chemical equilibrium between the two phases and is equivalent to a minimization of the free energy of the total system (system + reservoir). In the case where one requires electroneutrality of the phases, due to the additional constraint, the electrochemical potential $\overline{\mu}_i=\mu_i+z_ie\psi$ rather than the chemical potential has to be equal. Here, $z_i$ is the valency of type $i$, $e$ is the elementary charge and $\psi$ is the local electrostatic potential.
The grand-canonical distribution (see eq. (1)) can be sampled using the Metropolis Monte Carlo method, similiar to the canonical ensemble. Before we derive the acceptance criterion for the GCMC method, let us shortly recall the generic Metropolis-Hastings algorithm.
We want to calculate averages of observables $A$ which are distributed according to some distribution $p^\mathrm{eq}_s$ (e.g. the canonical or grand-canonical distribution)
$$ \langle A\rangle = \frac{\sum_{s\in\Gamma} A_s p^\mathrm{eq}_s}{\sum_{s\in\Gamma} p^\mathrm{eq}_s}. \quad \tag{2} \label{eq:can_average_mc} $$Here, $s$ labels the different states of the systems, $A_s$ is the value of the observable $A$ in the state $s$ and the sum runs over the set of all states $\Gamma$. The sum can also be an integral for systems with continuous degrees of freedom, for this derivation it does not matter. For almost all systems, the ensemble average (2) can not be evaluated exactly and only numerical approximations can be obtained. The basic idea of all MC methods is to randomly sample only a subset of states. In the most naive approach (simple sampling), states are simply selected with a uniform probability from the space $\Gamma$ of states $s$, resulting in a sample $\gamma\subset\Gamma$. The canonical average can then be approximated by a sum over $\gamma$:
$$ \langle A\rangle \approx \frac{\sum_{s\in\gamma} A_s p^\mathrm{eq}_s}{\sum_{s\in\gamma} p^\mathrm{eq}_s} \quad \tag{3} $$This method is highly inefficient, as most states will typically not contribute significantly to the average. To get a more efficient sampling method, let us select the states $s$ according to a non-uniform probability distribution $p_s$. Then, the average can be written as
$$ \langle A\rangle^{\mathrm{can}} \approx \frac{\sum_{s\in\gamma} A_s p^\mathrm{eq}_s/p_s}{\sum_{s\in\gamma} p^\mathrm{eq}_s/p_s}, \quad \tag{4} \label{eq:can_average_mc_prob} $$where we have to divide in the numerator and in the denominator by $p_s$ to compensate for the fact that the samples are distributed non-uniformly. It is obvious that we recover (3) from (4) as the special case of a uniform distribution. To sample the averages more efficiently, we simply pick
$$ p_s \propto p^\mathrm{eq}_s \quad \tag{5} \label{eq:mc_prob} $$such that the selected samples are representative of the desired equilibrium distribution. In this case, the estimators of the averages are thus simply given by
$$ \langle A\rangle \approx \frac{1}{N}\sum_{s\in\gamma}A_s. $$The remaining problem is now to find a method which allows us to sample states according to (5). This can be done in terms of so-called Markov chains. A Markov chain is a stochastic process without memory. In our case, the Markov chain is simply a sequence $\{s_i\}_{i=1,...N}$ of $N$ states $s_i$ of the system and the fact that the Markov chain has no memory means that the transition probability from a state $s_i$ to a state $s_j$ can be given in terms of a transition matrix $W_{i\rightarrow j}$ which is constant. The equilibrium probability distribution has to be invariant under $W$. Mathematically this condition can be expressed in the form
$$ \sum_{i}W_{i\rightarrow j}p_i^\mathrm{eq} = p_j^\mathrm{eq} $$A sufficient constraint on $W_{i\rightarrow j}$ which ensures the validity of this condition is the so-called criterion of detailed balance:
$$ \frac{W_{i\rightarrow j}}{W_{j\rightarrow i}} = \frac{p_j^\mathrm{eq}} {p_i^\mathrm{eq}}. $$Typically, the transition matrix elements are factorized into a proposal probability $g_{i\rightarrow j}$ and acceptance probability $P_{i\rightarrow j}^\mathrm{acc}$:
$$ W_{i\rightarrow j} = g_{i\rightarrow j} P_{i\rightarrow j}^\mathrm{acc}. $$For the symmetric choice $g_{i\rightarrow j}=g_{j\rightarrow i}$, the condition of detailed balance reduces to
$$ \frac{P_{i\rightarrow j}^\mathrm{acc}}{P_{j\rightarrow i}^\mathrm{acc}} = \frac{p_j^\mathrm{eq}}{p_i^\mathrm{eq}}. \quad \tag{6} $$By checking the cases $p_j^\mathrm{eq}>p_i^\mathrm{eq}$ and $p_j^\mathrm{eq} \leq p_i^\mathrm{eq}$ one can see that the acceptance probability
$$ P_{i\rightarrow j}^\mathrm{acc,\text{ }Metropolis} = \min\left(1,\frac{p_j^\mathrm{eq}}{p_i^\mathrm{eq}}\right) \quad \tag{7} $$fulfills equation (6). With this, we can now formulate the Metropolis-Hastings algorithm [2]. We start with an arbitrary initial state $s_1$ and generate a sequence $\{s_i\}_{i=1,...N}$ of $N$ states. Each new state (n) is generated from the previous state (o) in the following way. First, we propose a new state according to the symmetric proposal probability $g_{\mathrm{o}\rightarrow\mathrm{n}}$. Then, we accept the new configuration with a probability given by equation (7). If the proposed new state is rejected, the old state is retained. Applying this algorithm many times results in the desired Markov chain of representative configurations.
To simulate a system in a grand-canonical ensemble the described algorithm can be applied the grand-canonical distribution [1]. When considering the possible states of the system, it becomes clear that every state can be reached by a combination of two kinds of moves:
Displacement moves that leave the particle numbers fixed and only change positions. This is equivalent to sampling in the canonical ensemble with the current particle number. Due to the ergodic hypothesis this operation can be performed either using canonical Monte Carlo methods or canonical molecular dynamic simulations (e. g. Langevin dynamics). In this tutorial we choose the latter.
Insertion and Deletion moves that add or delete a single particle while keeping all other particle positions fixed.
In the GCMC method the insertion and deletion moves are implemented in the following fashion. First of all it is determined randomly with equal probability if an insertion or deletion move is performed. In the case of an insertion the additional particle is placed at a uniformly chosen random position in the simulation box. Then the proposed new state is accepted according to the criterion
$$ \begin{split} P_{N_i\rightarrow N_i+1}^\mathrm{acc,\text{ }GCMC} &= \min\left(1,\frac{1}{N_i+1}\left(\frac{V}{\Lambda_i^3}\right)\times \exp\left(-\beta \left(U^{N_i+1}(\mathbf{q})-U^{N_i}(\mathbf{q})\right)+\beta \mu_i\right)\right). \end{split} $$In the case of a deletion a randomly selected particle is removed from the system. We have the analogous result
$$ \begin{split} P_{N_i\rightarrow N_i-1}^\mathrm{acc,\text{ }GCMC} &= \min\left(1,N_i\left(\frac{\Lambda_i^3}{V}\right)\times \exp\left(-\beta \left(U^{N_i-1}(\mathbf{q})-U^{N_i}(\mathbf{q})\right)-\beta \mu_i\right)\right). \end{split} $$Exercise
Hint
For the case ($N_i\rightarrow N_i + 1$) one gets
\begin{align} P_{N_i\rightarrow N+1}^\mathrm{acc,\text{ }GCMC} &= \min\left(1,\frac{p^\mathrm{eq}_{N_i+1}}{p^\mathrm{eq}_{N_i}}\right) = \min\left(1,\frac{\frac{1}{(N_i+1)!}\left(\frac{V}{\Lambda_i^3}\right)^{N_i+1}\times \exp\left(-\beta U^{N_i+1}(\mathbf{q})+\beta \mu (N_i+1)\right)}{\frac{1}{N_i!}\left(\frac{V}{\Lambda_i^3}\right)^{N_i}\times \exp\left(-\beta U^{N_i}(\mathbf{q})+\beta \mu N_i\right)}\right)\\ &= \min\left(1,\frac{1}{N_i+1}\left(\frac{V}{\Lambda_i^3}\right)\times \exp\left(-\beta \left(U^{N_i+1}(\mathbf{q})-U^{N_i}(\mathbf{q})\right)+\beta \mu_i\right)\right). \end{align}For the other case ($N_i\rightarrow N_i-1$) we have the analogous result
\begin{align} P_{N_i\rightarrow N_i-1}^\mathrm{acc,\text{ }GCMC} &= \min\left(1,\frac{p^\mathrm{eq}_{N_i-1}}{p^\mathrm{eq}_{N_i}}\right) = \min\left(1,\frac{\frac{1}{(N_i-1)!}\left(\frac{V}{\Lambda_i^3}\right)^{N_i-1}\times \exp\left(-\beta U^{N_i-1}(\mathbf{q})+\beta \mu_i (N_i-1)\right)}{\frac{1}{N_i!}\left(\frac{V}{\Lambda_i^3}\right)^{N_i}\times \exp\left(-\beta U^{N_i}(\mathbf{q})+\beta \mu_i N_i\right)}\right)\\ &= \min\left(1,N_i\left(\frac{\Lambda_i^3}{V}\right)\times \exp\left(-\beta \left(U^{N_i-1}(\mathbf{q})-U^{N_i}(\mathbf{q})\right)-\beta \mu_i\right)\right). \end{align}Here we want to investigate a polyelectrolyte solution that is coupled to a reservoir containing a monovalent salt, as shown schematically in figure 1. Similar systems are useful to investigate for example the swelling behavior of hydrogels. The surrounding of the gel can be considered as the reservoir and the inner of the gel is the system. The small salt ions (X$^+$ and X$^-$) can be exchanged between the two phases while the polyelectrolyte is confined to the system. Experimentally such a setup could be realized by seperating the two phases with a semi-permeable membrane. Similiar setups also occur in the study of polyelectrolyte hydrogels.
Due to the macroscopic electroneutrality of both phases the salt ions obey an electrochemical equilibrium, i. e. there is a difference in the electrostatic potential between the two phases
$$ \mu_i^{\mathrm{sys}}+z_ie\psi^{\mathrm{Donnan}} = \mu_i^{\mathrm{res}}. $$This effect is called the Donnan effect [3] and leads to an unequal partioning of salt between the reservoir and the system. The partitioning can be quantified using the partition coefficients $\xi_i \equiv c^{\mathrm{sys}}_i/c^{\mathrm{res}}_i$, where $c^{\mathrm{sys}}_i$ and $c^{\mathrm{res}}_i$ are the particle concentrations of type $i$ in the system and reservoir, respectively. Defining the "universal" partition coefficient,
\begin{align*} \xi\equiv\left\{ \begin{array}{ll} \displaystyle\xi_+ - \frac{c_{\text{M}^-}^{\text{sys}}}{c_{\text{X}^+,\text{X}^-}^{\text{res}}}& \text{for positive ions}\\ \displaystyle\xi_- & \text{for negative ions,}\\ \end{array}\right. \end{align*}one can obtain the following solution for an ideal system without any interaction between the particles:
$$ \xi^{\text{Donnan}} = -\frac{c_{\text{M}^-}^{\text{sys}}}{2c_{\text{X}^+,\text{X}^-}^{\text{res}}} + \sqrt{\left(\frac{c_{\text{M}^-}^{\text{sys}}}{2c_{\text{X}^+,\text{X}^-}^{\text{res}}}\right)^2+1}. $$This solution shows that salt is rejected by the polyelctrolyte solution. In the following we want to investigate this partitioning for a system with electrostatic interactions.
Because of the electroneutrality condition, the two definitions of the "universal" partition coefficients for positive and negative ions are equal. The additional term in the definition for positive ions is a correction due to the counterions in the system, which neutralize the charged polyelectrolytes.
In the folllowing, we plot the universal partition coefficient over the ratio of the concentrations.
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(42)
plt.rcParams.update({"font.size": 18})
ratios = np.logspace(-1, 2, 1000)
plt.figure(figsize=(10, 6))
plt.loglog(ratios, -ratios + np.sqrt((ratios)**2 + 1), label=r"$\xi^{\mathrm{Donnan}}$")
plt.xlabel(r"$\dfrac{c_{\mathrm{M}^-}^{\mathrm{sys}}}{2c_{\mathrm{X}^+,\mathrm{X}^-}^{\mathrm{res}}}$")
plt.ylabel(r"$\xi$")
plt.legend()
plt.show()
Let us first add required modules.
import tqdm # module for progress bar
import scipy.interpolate
import pint # module for unit conversions
import espressomd
import espressomd.interactions
import espressomd.electrostatics
import espressomd.reaction_methods
import espressomd.polymer
espressomd.assert_features(["ELECTROSTATICS", "P3M", "WCA"])
In the next step we define the Lennard-Jones units. We use the module pint in order to make unit conversions between reduced simulation units and SI units easier.
##### Units
ureg = pint.UnitRegistry()
sigma = 0.355 * ureg.nm # Sigma in SI units.
# This is corresponds to the half of the Bjerrum length in water (eps_r approx 80) and room temperature (293 K)
T = 298.15 * ureg.K # 25 degree celcius corresponds to approx room temperature
CREF_IN_MOL_PER_L = 1.0 * ureg.molar #in mol/l
# Define the reduced units
ureg.define(f"reduced_length = {sigma}")
ureg.define(f"reduced_charge = 1* e")
ureg.define(f"reduced_energy = {T} * boltzmann_constant")
We load the resulting excess chemical potentials for different salt concentrations of the reservoir from the Widom insertion tutorial. The loaded excess chemical potentials are in units of $k_\mathrm{B}T$ (i.e. in reduced units).
salt_concentration_magnitudes_si = [0.001, 0.003, 0.01, 0.03, 0.1]
salt_concentration_si = np.array(salt_concentration_magnitudes_si) * ureg.molar
salt_concentration_sim = salt_concentration_si.to("molecules/reduced_length^3")
excess_chemical_potential_data = [
-0.0672291585811938, -0.123732052028611, -0.218687259792629,
-0.326207586236404, -0.510091568808112
]
excess_chemical_potential_data_error = [
0.000994798843587, 0.00152160176511698, 0.00220162667953136,
0.0029299206854553, 0.00422309102221265
]
excess_chemical_potential_monovalent_pairs_in_bulk = scipy.interpolate.interp1d(
salt_concentration_sim.magnitude, excess_chemical_potential_data)
fig = plt.figure(figsize=(10, 6))
plt.errorbar(salt_concentration_si.magnitude,
excess_chemical_potential_data,
yerr=excess_chemical_potential_data_error,
linestyle="-",
marker="o",
markersize=5.,
capsize=5.)
plt.xlabel(r"salt concentration $c_{\mathrm{salt}}$ (mol/l)")
plt.ylabel(r"excess chemical potential $\mu_{\mathrm{ex}}$ ($k_{\mathrm{B}}T$)")
plt.xscale("log")
plt.xlim((5e-4, 0.5))
plt.ylim((-0.8, 0));
Now we are ready to set up the system. In the first step we define the system parameters such as charges and system size as well as simulation paramters such as the time step. For small salt concentrations we need to produce more samples to get good statistics for $\xi$ due to the finite system size. For a real production run it is recommended to use larger sample sizes.
##### Particle types and charges
types = {
"Xplus": 1, # Positive ions
"Xminus": 2, # Negative ions
"Monomer": 3, # Negative monomers
}
charges = {
"Xplus": 1.0,
"Xminus": -1.0,
"Monomer": -1.0,
}
# Initialize concentrations
c_salt_res_sim = salt_concentration_sim[0]
c_monomer = 0.1 * ureg.molar
C_MONOMERS_SIM = c_monomer.to("molecules/reduced_length^3")
##### Number of monomers
N_CHAINS = 10
N_MONOMERS_PER_CHAIN = 5
##### System size
BOX_LENGTH = np.power(
N_CHAINS * N_MONOMERS_PER_CHAIN / C_MONOMERS_SIM.magnitude, 1.0 / 3.0)
##### Initial number of salt ion pairs
n_salt = int(c_salt_res_sim.magnitude * BOX_LENGTH**3)
##### Weeks-Chandler-Andersen interaction
LJ_EPSILON = 1.0 # in reduced units, i.e. 1 sigma
LJ_SIGMA = 1.0 # in reduced units, i.e. 1 k_B T
##### Electrostatic interaction
L_BJERRUM = 2.0 # in reduced units, i.e. 2 sigma
##### Langevin-Thermostat
KT = 1.0 # in reduced units
GAMMA = 1.0 # in reduced units
##### Integrator parameters
DT = 0.01 # Integrator time step
SKIN = 0.4
warmup_loops = 50
number_of_loops_for_each_concentration = [
200, 200, 150, 150, 100
] # set these higher for better accuracy
steps_per_loop = 200
samples_per_loop = 25
##### Seeding
langevin_seed = 42
reaction_seed = 42
Now we define our box and add our particles. We create the polymer chains as in the corresponding tutorial. In addition, we enable the electrostatic interactions.
##### Create an instance of the system class
system = espressomd.System(box_l=[BOX_LENGTH] * 3)
system.time_step = DT
system.cell_system.skin = SKIN
print("Created system class")
#### Add the FENE interaction
fene = espressomd.interactions.FeneBond(k=30, d_r_max=1.5)
system.bonded_inter.add(fene)
#### Add the polymer chains
polymer_positions = espressomd.polymer.linear_polymer_positions(
n_polymers=N_CHAINS,
beads_per_chain=N_MONOMERS_PER_CHAIN,
bond_length=0.9,
seed=42)
for positions in polymer_positions:
monomers = system.part.add(pos=positions,
type=[types["Monomer"]] * N_MONOMERS_PER_CHAIN,
q=[charges["Monomer"]] * N_MONOMERS_PER_CHAIN)
previous_part = None
for part in monomers:
if not previous_part is None:
part.add_bond((fene, previous_part))
previous_part = part
##### Add the particles
N_plus = n_salt + N_CHAINS * N_MONOMERS_PER_CHAIN
system.part.add(type=[types["Xplus"]] * N_plus,
pos=np.random.rand(N_plus, 3) * BOX_LENGTH,
q=[charges["Xplus"]] * N_plus)
system.part.add(type=[types["Xminus"]] * n_salt,
pos=np.random.rand(n_salt, 3) * BOX_LENGTH,
q=[charges["Xminus"]] * n_salt)
print("Added particles")
##### Add non-bonded interactions
for i in types:
for j in types:
system.non_bonded_inter[types[i],
types[j]].wca.set_params(epsilon=LJ_EPSILON,
sigma=LJ_SIGMA)
print("Added WCA interaction")
##### Minimize forces
system.integrator.set_steepest_descent(f_max=N_plus, gamma=10.0,
max_displacement=LJ_SIGMA / 100.)
n_steps = system.integrator.run(steps_per_loop)
assert n_steps < steps_per_loop
print("Removed overlaps with steepest descent...")
##### Add thermostat and long-range solvers
system.thermostat.set_langevin(kT=KT, gamma=GAMMA, seed=langevin_seed)
system.integrator.set_vv()
p3m_params = {'mesh':10,'cao':6,'r_cut':8.22}
p3m = espressomd.electrostatics.P3M(prefactor=L_BJERRUM * KT, accuracy=1e-3, **p3m_params)
system.electrostatics.solver = p3m
system.integrator.run(2 * steps_per_loop)
print("Added electrostatics")
Created system class Added particles Added WCA interaction Removed overlaps with steepest descent... Added electrostatics
Now we want to add a coupling to the reservoir. We can formally represent the insertion and deletion of ion pairs as a chemical reaction:
$$ \emptyset \rightleftharpoons \mathrm{X}^- + \mathrm{X}^+ $$It is important to note that the reaction can run in both directions, i. e. insertion and deletion of salt ions. These reactions are described by a concentration based equilibrium constant:
$$ K = c_{\mathrm{salt,res}}^2 \exp \left ( \mu_{\mathrm{salt,res}}^{\mathrm{ex}} \right ) $$Here $\mu_{\mathrm{salt,res}}^{\mathrm{ex}}$ denotes the excess chemical potential of a salt ion pair as we measured in the Widom insertion tutorial.
Exercise
Hint
# SOLUTION CELL
K_XX = c_salt_res_sim.magnitude**2 * np.exp(
excess_chemical_potential_monovalent_pairs_in_bulk(
c_salt_res_sim.magnitude))
RE = espressomd.reaction_methods.ReactionEnsemble(kT=KT,
exclusion_range=1.0,
seed=reaction_seed)
RE.add_reaction(gamma=K_XX,
reactant_types=[],
reactant_coefficients=[],
product_types=[types["Xplus"], types["Xminus"]],
product_coefficients=[1, 1],
default_charges={
types["Xplus"]: charges["Xplus"],
types["Xminus"]: charges["Xminus"]
})
# Set the non interacting type at the smallest integer.
# This may speed up the simulation (see Espresso docummentation section 19. Reaction methods)
RE.set_non_interacting_type(type=len(types) + 1)
Now we are ready to perform the actual sampling.
Exercise
perform_sampling(salt_concentration)
that takes as input the salt concentration of the reservoir.Hint
# SOLUTION CELL
def perform_sampling(salt_concentration, number_of_loops):
K_XX = salt_concentration**2 * np.exp(
excess_chemical_potential_monovalent_pairs_in_bulk(salt_concentration))
RE.change_reaction_constant(gamma=K_XX, reaction_id=0)
for i in tqdm.trange(warmup_loops, bar_format="{l_bar}{bar:60}{r_bar}{bar:-60b}"):
system.integrator.run(steps_per_loop)
RE.reaction(samples_per_loop)
particle_numbers_positive = []
particle_numbers_negative = []
for i in tqdm.trange(number_of_loops, bar_format="{l_bar}{bar:60}{r_bar}{bar:-60b}"):
system.integrator.run(steps_per_loop)
RE.reaction(samples_per_loop)
particle_numbers_positive.append(
system.number_of_particles(type=types["Xplus"]))
particle_numbers_negative.append(
system.number_of_particles(type=types["Xminus"]))
partition_coefficient_positive = np.mean(np.asarray(particle_numbers_positive))\
/ (BOX_LENGTH**3 * salt_concentration)
partition_coefficient_negative = np.mean(np.asarray(particle_numbers_negative))\
/ (BOX_LENGTH**3 * salt_concentration)
return partition_coefficient_positive, partition_coefficient_negative
Now we can perform the actual simulations for the different salt concentrations and measure the partition coefficients.
partition_coefficients_positives_array = np.zeros_like(salt_concentration_sim)
partition_coefficients_negatives_array = np.zeros_like(salt_concentration_sim)
for i, salt_concentration in enumerate(salt_concentration_sim.magnitude):
print(f"Salt concentration {i+1}/{len(salt_concentration_sim)}")
partition_coefficients_positives_array[
i], partition_coefficients_negatives_array[i] = perform_sampling(
salt_concentration, number_of_loops_for_each_concentration[i])
Salt concentration 1/5
100%|████████████████████████████████████████████████████████████| 50/50 [00:01<00:00, 25.49it/s] 100%|████████████████████████████████████████████████████████████| 200/200 [00:06<00:00, 32.64it/s]
Salt concentration 2/5
100%|████████████████████████████████████████████████████████████| 50/50 [00:01<00:00, 31.30it/s] 100%|████████████████████████████████████████████████████████████| 200/200 [00:06<00:00, 31.63it/s]
Salt concentration 3/5
100%|████████████████████████████████████████████████████████████| 50/50 [00:01<00:00, 30.02it/s] 100%|████████████████████████████████████████████████████████████| 150/150 [00:05<00:00, 29.78it/s]
Salt concentration 4/5
100%|████████████████████████████████████████████████████████████| 50/50 [00:01<00:00, 26.24it/s] 100%|████████████████████████████████████████████████████████████| 150/150 [00:05<00:00, 26.17it/s]
Salt concentration 5/5
100%|████████████████████████████████████████████████████████████| 50/50 [00:03<00:00, 16.27it/s] 100%|████████████████████████████████████████████████████████████| 100/100 [00:06<00:00, 15.53it/s]
To compare the results of our simulations we define a function for the analytical solution.
def analytical_solution(ratio):
return -ratio + np.sqrt(ratio**2 + 1)
Then we can measure the partition coefficients derived from the simulations and compare them to the analytical results for an ideal system.
universal_partion_coefficient_positive = \
partition_coefficients_positives_array - c_monomer.magnitude / salt_concentration_si.magnitude
fig = plt.figure(figsize=(10, 6))
ratios = np.logspace(-1, 2, num=200)
plt.loglog(ratios, analytical_solution(ratios), color="black", label=r"$\xi^{\mathrm{Donnan}}$")
plt.loglog(c_monomer.magnitude / (2 * salt_concentration_si.magnitude), partition_coefficients_negatives_array,
marker="o", markersize=10, markeredgewidth=1.5, markeredgecolor="k", markerfacecolor="none",
linestyle="none", label=r"$\xi_-$")
plt.loglog(c_monomer.magnitude / (2 * salt_concentration_si.magnitude), universal_partion_coefficient_positive,
marker="o", markersize=4.5, markeredgewidth=0., markeredgecolor="k", markerfacecolor="k",
linestyle="none", label=r"$\xi_+ - c_{\mathrm{M}^-}^{\mathrm{sys}}/c_{\mathrm{salt}}^{\mathrm{res}}$")
plt.xlabel(r"$\dfrac{c_{\mathrm{M}^-}^{\mathrm{sys}}}{2c_{\mathrm{salt}}^{\mathrm{res}}}$")
plt.ylabel(r"$\xi$")
plt.legend()
plt.show()
Exercise
First, we note that the data points for negative and positive ions always collapse onto the same value for ξ, i.e. ξ is indeed a universal value that is always the same for negative and positive ions. When comparing these results to the analytical prediction for an ideal gas, it is easy to see that the partition coefficient is higher for the case with electrostatic interactions. This is essentially due to the fact that the overall charge density inside the system is higher than in the reservoir, and thus the activity coefficient of an ion pair inside the system is smaller than in the reservoir, i.e. the free energy cost of putting additional ion pairs into the system is overall lowered as compared to an ideal gas. Still, we see that ξ < 1 even for the interacting case, i.e. the salt concentration is still smaller inside the system than in the reservoir.
[1] Daan Frenkel, Berend Smit. Understanding Molecular Simulation: From Algorithms to Applications. 2nd edition, chapter 5: Monte Carlo Simulations in Various Ensembles, section 5.6: Grand-Canonical Ensemble, pp. 126–135. Academic Press, 2002, ISBN: 978-0-12-267351-1, doi:10.1016/B978-012267351-1/50007-9.
[2] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller. Equation of state calculations by fast computing machines. J. Chem. Phys. 21, pp. 1087–1092 (1953), doi:10.1063/1.1699114.
[3] F. Donnan. The theory of membrane equilibria. Chemical Reviews 1(1), pp. 73–90 (1924), doi:10.1021/cr60001a003.