The combination of Monte Carlo and Molecular Dynamic simulations offers many possibilities beyond using each technique solely on their own. In general, MD simulations are restricted to relatively small time scales and, depending on the modeled interactions, to rather small particle numbers as well. However, they have the advantage of computing collective particle movements with the correct dynamics. The time between two MD steps is bound by the fastest motion of the system of interest (for all-atom MD, those are atomic vibrations).

MC simulations, on the other hand, do not have such limitations. Even nonphysical moves can be considered, still leading to results with correct thermodynamic properties. However, no valuable information of the system dynamics can be obtained. The combination of MD and MC moves allows us to mimic any ensemble.

In this tutorial, you will compute the excess chemical potential of a monovalent salt solution by using a combination of MD and MC techniques. For this purpose, you will learn about the Widom particle insertion method [1] and how it can be used in ESPResSo. The chemical potential is important for simulating chemical equilibria and systems where particles can be exchanged with an external reservoir, which is too large to be simulated explicitly, or between two subsystems.

In a pure MD simulation, there is no direct way to access the chemical potential, since it is related to changes in the particle number that is typically fixed. In 1963, Widom proposed a MC scheme to measure the chemical potential in a system using trial particle insertions [1].

The chemical potential for species $\alpha$ is defined as \begin{equation} \mu_{\alpha} = \left(\frac{\partial G}{\partial N_\alpha} \right)_{p,T,N_{\beta\neq\alpha}} = \left(\frac{\partial F}{\partial N_\alpha}\right)_{V,T,N_{\beta\neq\alpha}} \quad ,\tag{1} \label{eq:chem_pot_def} \end{equation} and corresponds to the difference in the Helmholtz free energy $F$ or Gibbs free energy $G$ when adding or removing one particle in the system.

The canonical partition function is defined as \begin{equation} Z(N,V,T)=\frac{V^N}{\Lambda^{3N}N!} \int_{\left[0,1\right]^{3N}}\exp\left(-\beta U\left(\left\{\boldsymbol{s}^i\right\}\right)\right)\,\mathrm{d}\boldsymbol{s}^1\ldots\mathrm{d}\boldsymbol{s}^N \quad ,\tag{2} \label{eq:partition_function} \end{equation} where V is the volume, $\Lambda=h/\sqrt{2\pi m k_\mathrm{B}T}$ the thermal wavelength, $\beta$ the inverse temperature, $U$ the interaction potential and $\boldsymbol{s}^i$ the rescaled (dimensionless) particle coordinates.

The Helmholtz free energy $F$ can be computed from the partition function: $F=-k_BT\ln\left(Z(N,V,T)\right)$. For sufficiently large N, the following relation approximates the chemical potential: \begin{align} \mu & = -k_\mathrm{B}T \ln\left(\frac{Z(N+1,V,T)}{Z(N,V,T)}\right) \label{eq:mu_limit_largeN}\\ & = \mu_\textrm{ideal}+\mu_\textrm{ex} \nonumber\\ & = - k_\mathrm{B}T\ln\left(\frac{V}{\Lambda^{3}(N+1)}\right)+\mu_\textrm{ex}\, , \nonumber\tag{3} \end{align} where $\mu_\textrm{ideal}$ depends only on the density of the system, so it is a constant for systems with different interactions but the same particle density. The interesting part is the excess chemical potential $\mu_\textrm{ex}$, which can be expressed in a form involving a canonical average: \begin{equation} \mu-\mu_\textrm{ideal}=\mu_\textrm{ex}=-k_\mathrm{B}T\ln\left(\int_{\left[0,1\right]^{3}}\langle\exp(-\beta\Delta U)\rangle\,\mathrm{d}\boldsymbol{s}^{N+1}\right)\quad .\tag{4} \label{eq:mu_ex} \end{equation}

In the above expression, $\Delta U$ is the energy difference of the system before and after the insertion of an additional particle ("test particle") into the system, and $\langle...\rangle$ is the canonical average for a system with $N$ particles. The expression thus suggests that the excess chemical potential can be computed using the following procedure. The canonical average $\langle\exp(-\beta\Delta U)\rangle$ can be computed by inserting a test particle into many different configurations of the system (sampled according to the canonical distribution) and then averaging the associated Boltzmann factor $\exp(-\beta\Delta U)$. Note that this canonical average can be computed using either MC or canonical MD, in the following we will use Langevin MD for this task. The integral over $\boldsymbol{s}^{N+1}$ corresponds to an average over different insertion positions of the test particle and can be calculated using a brute-force MC approach by simply inserting test particles at many different (random) trial positions for each sampled state of the system.

The system under consideration here consists of $N$ monovalent ion pairs in an implicit solvent with $\varepsilon_{\mathrm{r}} = 78$, which is a simple model ("restricted primitive model") for an aqueous salt solution of e.g. NaCl. A Weeks–Chandler–Andersen potential is used to avoid particles overlapping one another. In the unit system employed in the simulation, the following parameters will be used:

\begin{align*} \lambda_b=2.0\,,\quad k_BT=1.0\,,\quad z_i=1.0\,,\quad N=N_i=50 \end{align*}Since the investigated system is a monovalent salt, the concentrations (and thus, densities) of positive and negative ions are the same. The Widom insertion method can be used to calculate the excess chemical potential for this system. However, we cannot calculate the excess chemical potential of a single ion (a quantitity that is also not accessible in experiments) but only that of an ion pair, which is given by \begin{align} \mu_\textrm{ex}^\textrm{pair} = \mu_\textrm{ex}^+ + \mu_\textrm{ex}^- = \mu-\mu_\textrm{ideal}=\mu_\textrm{ex}=-k_BT\ln\left(\int_{\left[0,1\right]^{6}}\langle\exp(-\beta\Delta U)\rangle \mathrm{d}\boldsymbol{s}_{+}^{N+1}\mathrm{d}\boldsymbol{s}_{-}^{N+1}\right).\tag{5} \end{align} Inserting ion pairs ensures the overall electroneutrality of the system and prevents a diverging electrostatic energy.

In a first step, we import the required **ESPResSo** features as well as some additional python modules:

In [1]:

```
import espressomd
import espressomd.electrostatics
import espressomd.reaction_methods
espressomd.assert_features(["ELECTROSTATICS", "P3M", "WCA"])
import tqdm
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
np.random.seed(42)
plt.rcParams.update({"font.size": 18})
```

In [2]:

```
# we assume a reduced unit system where the elementary charge,
# the thermal energy and the WCA length scale (sigma) are all unity
sigma = 3.55e-10 # sigma in SI units
avo = 6.022e+23 # Avogadro's number in SI units
PREF = 1 / (10**3 * avo * sigma**3) # prefactor to mol/L
BJERRUM_LENGTH = 2.0
# Weeks-Chandler-Andersen interaction
LJ_EPSILON = 1.0
LJ_SIGMA = 1.0
# Langevin thermostat
KT = 1.0
GAMMA = 1.0
# number of salt ion pairs
N_ION_PAIRS = 50
MASS=1.0
# particle types and charges
types = {
"Xplus": 1, # Positive ions
"Xminus": 2, # Negative ions
}
charges = {
"Xplus": 1.0,
"Xminus": -1.0,
}
```

Now we are ready to set up the system. Because we will need to rescale the simulation box, we will initially only set up the short-range interactions.

**Exercise:**

- Set up a system with box length $10\,\sigma$, integrator time step $\Delta t=0.01 \,\tau$ ($\tau$ is the Lennard-Jones timescale, i.e. equal to 1 in the employed unit system) and Verlet skin width of $\delta = 0.4\sigma$
- Add a total of
`N_ION_PAIRS`

ion pairs at random positions to the simulation box - Add a WCA interaction (purely repulsive Lennard-Jones interaction) between the particles

**Hints:**

- Refer to the documentation to find out how to set up the WCA interaction (Non-bonded Interactions)

In [3]:

```
# SOLUTION CELL
box_l_init = 10.0
system = espressomd.System(box_l=3 * [box_l_init])
system.time_step = 0.01
system.cell_system.skin = 0.4
system.part.add(pos=np.random.rand(N_ION_PAIRS, 3) * box_l_init,
type=[types["Xplus"]] * N_ION_PAIRS,
q=[charges["Xplus"]] * N_ION_PAIRS)
system.part.add(pos=np.random.rand(N_ION_PAIRS, 3) * box_l_init,
type=[types["Xminus"]] * N_ION_PAIRS,
q=[charges["Xminus"]] * N_ION_PAIRS)
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)
```

**Exercise:**

- Implement a function
`system_setup(c_salt_SI)`

that takes the desired salt concentration`c_salt_SI`

in SI-units (mol/l) as an argument and rescales the simulation box accordingly - Afterwards, the function should also add the P3M actor to account for electrostatic interactions

**Hints:**

- Use the prefactor
`PREF`

defined above to convert the concentration from SI-units to the unit system employed in the simulation - An accuracy of $10^{-3}$ for the P$^3$M should be enough for this exercise

In [4]:

```
# SOLUTION CELL
ci_params = {'mesh':(12,12,12),'cao':6}
def system_setup(c_salt_SI):
print(f"Salt concentration: {c_salt_SI:.12g} mol/l", flush=True)
c_salt_sim = c_salt_SI / PREF
box_length = np.cbrt(N_ION_PAIRS / c_salt_sim)
# rescale box
system.change_volume_and_rescale_particles(d_new=box_length)
print("Rescaled the simulation box.", flush=True)
# add P3M
p3m = espressomd.electrostatics.P3M(prefactor=BJERRUM_LENGTH, accuracy=1e-3, **ci_params)
system.electrostatics.solver = p3m
```

**Exercise:**

- Implement a function
`warmup()`

that removes potential overlaps between particles with 10000 steps of the steepest descent integrator and performs a warmup integration with 100 loops of 100 simulation steps

**Hints:**

- keep in mind that after the overlaps have been removed, the integrator should be switched
back to Velocity-Verlet with the Langevin thermostat (using
`GAMMA`

and`KT`

) in order to perform a warmup integration - to show the progression of
`for`

loops, use the syntax`for i in tqdm.tqdm(range(100))`

In [5]:

```
# SOLUTION CELL
def warmup():
FMAX = 0.01 * LJ_SIGMA * MASS / system.time_step**2
system.integrator.set_steepest_descent(
f_max=FMAX,
gamma=1e-3,
max_displacement=0.01)
print("Remove overlaps...", flush=True)
system.integrator.run(5000)
assert np.all(np.abs(system.part.all().f) < FMAX), "Overlap removal did not converge!"
system.integrator.set_vv()
system.thermostat.set_langevin(kT=KT, gamma=GAMMA, seed=42)
print("Running warmup integration...", flush=True)
for i in tqdm.tqdm(range(100)):
system.integrator.run(100)
```

**Exercise:**

- Add the functionality to perform Widom insertion moves. Make sure that you always insert an ion pair in order to conserve the system electroneutrality.

**Hints:**

- Refer to the documentation to find out how to set up Widom particle insertion (Widom Insertion)

In [6]:

```
# SOLUTION CELL
# add reaction for Widom particle insertion
widom = espressomd.reaction_methods.WidomInsertion(kT=KT, seed=42)
widom.add_reaction(reactant_types=[], reactant_coefficients=[],
product_types=[types["Xplus"], types["Xminus"]], product_coefficients=[1, 1],
default_charges={types["Xplus"]: charges["Xplus"], types["Xminus"]: charges["Xminus"]})
widom.set_non_interacting_type(type=len(types) + 1)
```

**Exercise:**

- Implement a function
`calculate_excess_chemical_potential(n_md_steps, n_mc_steps, sample_size)`

that performs`n_mc_steps`

Widom particle insertions every`n_md_steps`

MD steps in a loop that repeats`sample_size`

times and returns the chemical potential and its error as a tuple via widom.calculate_excess_chemical_potential() - Implement a function
`system_teardown()`

that removes the electrostatics solver and turns off the thermostat

In [7]:

```
# SOLUTION CELL
def calculate_excess_chemical_potential(n_md_steps, n_mc_steps, sample_size):
potential_energy_samples = []
print("Sampling...", flush=True)
for i in tqdm.tqdm(range(sample_size)):
system.integrator.run(n_md_steps)
for j in range(n_mc_steps):
e_pot = widom.calculate_particle_insertion_potential_energy(reaction_id=0)
potential_energy_samples.append(e_pot)
excess_chemical_potential, error_excess_chemical_potential = \
widom.calculate_excess_chemical_potential(
particle_insertion_potential_energy_samples=potential_energy_samples)
print(f"Excess chemical potential: {excess_chemical_potential:.4g}\n", flush=True)
return excess_chemical_potential, error_excess_chemical_potential
def system_teardown():
system.electrostatics.clear()
system.thermostat.turn_off()
```

In [8]:

```
salt_concentrations = np.array([0.003, 0.01, 0.03, 0.1, 0.3])
excess_chemical_potentials = []
sample_size = 250
import time
start = time.time()
for c_salt_SI in salt_concentrations:
system_setup(c_salt_SI)
warmup()
mu_ex = calculate_excess_chemical_potential(n_md_steps=100, n_mc_steps=100, sample_size=sample_size)
excess_chemical_potentials.append(mu_ex)
system_teardown()
excess_chemical_potentials = np.asarray(excess_chemical_potentials)
end = time.time()
print(f"Elapsed time: {(end - start) / 60:.1f} min", flush=True)
```

Salt concentration: 0.003 mol/l Rescaled the simulation box. Remove overlaps... Running warmup integration...

100%|██████████| 100/100 [00:01<00:00, 82.66it/s]

Sampling...

100%|██████████| 250/250 [00:07<00:00, 34.43it/s]

Excess chemical potential: -0.1387 Salt concentration: 0.01 mol/l

Rescaled the simulation box. Remove overlaps... Running warmup integration...

100%|██████████| 100/100 [00:01<00:00, 82.08it/s]

Sampling...

100%|██████████| 250/250 [00:07<00:00, 31.82it/s]

Excess chemical potential: -0.2151

Salt concentration: 0.03 mol/l Rescaled the simulation box. Remove overlaps... Running warmup integration...

100%|██████████| 100/100 [00:01<00:00, 76.28it/s]

Sampling...

100%|██████████| 250/250 [00:09<00:00, 26.78it/s]

Excess chemical potential: -0.3335 Salt concentration: 0.1 mol/l Rescaled the simulation box. Remove overlaps... Running warmup integration...

100%|██████████| 100/100 [00:01<00:00, 68.54it/s]

Sampling...

100%|██████████| 250/250 [00:14<00:00, 17.35it/s]

Excess chemical potential: -0.5224 Salt concentration: 0.3 mol/l Rescaled the simulation box. Remove overlaps... Running warmup integration...

100%|██████████| 100/100 [00:01<00:00, 61.32it/s]

Sampling...

100%|██████████| 250/250 [00:14<00:00, 17.10it/s]

Excess chemical potential: -0.6911 Elapsed time: 1.0 min

We will now plot the measured excess chemical potential as a function of $c_{\mathrm{salt}}$ using a logarithmically scaled x-axis.

**Exercise:**

- Explain the observed behaviour qualitatively
- How do you expect the excess chemical potential to behave in the limit $c_{\mathrm{salt}}\rightarrow 0$ mol/l?

In [9]:

```
fig = plt.figure(figsize=(10, 8))
plt.errorbar(salt_concentrations,
excess_chemical_potentials[:, 0],
yerr=excess_chemical_potentials[:, 1],
linestyle="none", marker="o", markersize=5., capsize=5.,
label=r"Widom method")
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((2e-3, 0.6))
plt.ylim((-0.8, 0.0))
plt.legend()
plt.show()
```

There is a number of analytical and phenomenological theories to describe the excess chemical potential of salt solutions using analytical expressions [2], using the property that $\mu_{\textrm{ex}} = k_{\mathrm{B}}T \ln{\gamma_{\pm}}$, with $\gamma_{\pm}$ the mean activity coefficient of the salt. Examples are the extended Debye-Hückel (DH) equation,

\begin{equation} \mu_{\textrm{ex}}^{\textrm{DH}}=-2\cdot k_{\mathrm{B}}T\cdot \ln{10} \cdot \frac{A \sqrt{c_{\mathrm{salt}}/c^{\ominus}}}{1+B\sqrt{c_{\mathrm{salt}}/c^{\ominus}}} \end{equation}for ionic strengths up to 0.01 mol/l, and the Davies equation,

\begin{equation} \mu_{\textrm{ex}}^{\textrm{Davies}}=-2\cdot k_{\mathrm{B}}T\cdot \ln{10} \cdot A \cdot\left(\frac{\sqrt{c_{\mathrm{salt}}/c^{\ominus}}}{1+\sqrt{c_{\mathrm{salt}}/c^{\ominus}}}-C\cdot c_{\mathrm{salt}}/c^{\ominus}\right) \end{equation}which results from an empirical extension of the Debye-Hückel theory for ionic strengths up to 1 mol/l, with $c^{\ominus} = 1$ mol/l the standard concentration, and parameters $A = 1.82 \cdot 10^{6} (\varepsilon_0 \varepsilon_{\mathrm{r}} T)^{-3/2} \simeq 0.509$, $B \simeq 1.0$ and $C = 0.2$ in reduced units for NaCl in an aqueous solution at 25°C.

**Exercise:**

- Compare your simulation results with the Debye-Hückel theory and the Davies equation by plotting all three curves in a single plot

**Hint:**

- re-use the code from the previous figure
- create a logarithmic concentration range from $10^{-4}$ to $10^{0}$ mol/l with
`np.logspace(-4, 0.0, num=500, base=10)`

to plot the analytical solutions

In [10]:

```
# SOLUTION CELL
def davies_equation(c_salt):
A = 0.509
C = 0.2
mu_ex = - 2 * np.log(10) * A * (np.sqrt(c_salt) / (1 + np.sqrt(c_salt)) - C * c_salt)
return mu_ex
def dh_equation(c_salt):
A = 0.509
B = 1.0
mu_ex = - 2 * np.log(10) * (A * np.sqrt(c_salt) / (1 + B * np.sqrt(c_salt)))
return mu_ex
concentration_range = np.logspace(-4, 0.0, num=500, base=10)
fig = plt.figure(figsize=(10, 8))
plt.plot(concentration_range, davies_equation(concentration_range),
label=r"Davies equation", color="r")
plt.plot(concentration_range, dh_equation(concentration_range),
label=r"Debye-Hückel equation", color="black", linestyle="dashed")
plt.errorbar(salt_concentrations,
excess_chemical_potentials[:, 0],
yerr=excess_chemical_potentials[:, 1],
linestyle="none", marker="o", markersize=5., capsize=5.,
label=r"Widom method")
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((2e-3, 0.6))
plt.ylim((-0.8, 0.0))
plt.legend()
plt.show()
```

- Derive equation (4) starting from the canonical partition function.
- What problems would you run into if you tried to use the Widom insertion method with large molecules (e.g. polymers) as test particles?

[1] B. Widom. Some topics in the theory of fluids. Journal of Chemical Physics, 39:2808, 1963, doi:10.1063/1.1734110.

[2] P. Atkins, J. de Paula and J. Keeler. Atkins' Physical Chemistry, 11^{th} edition, Volume 3: Molecular Thermodynamics and Kinetics, Chapter 5F.4: The activities of ions, pp. 205–208. Oxford University Press, 2019, ISBN: 978-0-19-882-336-0.