Monte Carlo Simulation Exercise: Electrolyte Solutions

Mikael Lund, Division of Theoretical Chemistry, Lund university

Introduction to Jupyter Notebooks

This document is a Notebook consisting of code, documentation, and plots, each listed in a cell.

• Double click on a cell to edit it - try for example on this cell.
• Evaluate an active cell by pressing shift+return.
• The default cell language is Python, but by prefixing a command with !, it is interpreted as a BASH command.
• For getting help on a function, place the cursor inside the () brackets and press shift+tab-tab.
• More on text formatting, equations etc. here.
• Try to evaluate the cell just below this one; it should show a short movie about Monte Carlo simulations.

Download

git clone http://github.com/mlund/labs
    cd labs/excess

Installation

To open the Notebook, install python via Miniconda and make sure all required packages are loaded by issuing the following terminal commands,

conda env create -f environment.yml
    source activate excess
    jupyter-notebook excess.ipynb

Theory

Metropolis Monte Carlo Simulation

Litterature: Frenkel and Smith, Understanding Molecular Simulation, Chapters 3 and 5. [Open Access within LU]

Metropolis Monte Carlo is an importance sampling technique where we bias the otherwise random MC method to sample points in coordinate space that contribute to properties of interest. In molecular systems, configurations (microstates) with energies much higher than the thermal energy, $k_BT$, have a very low probability due to the Boltzmann factor, $P(U) = \exp(-U/k_BT)$, and therefore do not need to be frequently explored. In MC, we perform a finite number of trial moves where the system is taken from on old state ($o$) to a new state ($n$) giving rise to an energy change $\Delta U$. The move is accepted or rejected based in the following acceptance criterion,

which will ensure that low energy states are more likely to occur than high energy states. Thus, if the energy decreases ($e^{-\Delta U/k_BT}>1$) the move is always accepted. Temperature and entropy enters the algorithm by allowing acceptance also for moves that leads to an energy increase. As $T\rightarrow\infty$, all configurations are accepted leading to equal population of all states. When $T\rightarrow 0$ the system propagates towards a (local) minimum or ground state as expected from statistical thermodynamics. It is important to realize that the Metropolis MC algorithm only provide information on equilibrium properties. Thus from a MC simulation we cannot say anything about time-dependent or dynamic properties.

In the MC program we will use in this lab, two moves are performed:

1. Translational particle moves
2. Particle insertion and deletion

The first move we simply attempt to move a particle by changing its coordinates by a random unit vector, scaled by a displacement parameter. The second move is used to maintain a constant chemical potential of the salt particles so that we sample the Grand Canonical Ensemble ($\mu V T$). Here, each particle is associated with a chemical potential, $\mu$, which contribute to the trial energy when deleting or inserting particles. The acceptance rules become,

and $$ acc(N \rightarrow N-1) = \min \left(1, \frac{N}{V}e^{-(\Delta U+\mu)/k_BT}\right ) $$

The insertion/deletion move leads to a fluctuating number of particles in the simulation box, and we can therefore calculate an average concentraion. The ratio between this and the input activity gives us the activity coefficient which is the topic of the next section.

Excess Chemical Potential

In this lab you will use Metropolis Monte Carlo (MC) simulations to study properties of strong electrolyte solutions in water. We will use the primitive model of electrolytes where solvent, water, is treated as a dielectric continuum and salt particles as hard spheres. The inter-particle pair-potential as a function is separation $r$ is given by,

Question: insert particle pair-potential here

When salt is added to water, it dissociates into cations and anions, and due to the long-range nature of electrostatics ($1/r$), these influence each other over long distances. The particles thus cannot be regarded as ideal - especially as the concentration is increased. Non-ideality in chemical systems is captured by the activity coefficient, $\gamma_i$, which is a measure of the chemical potential ($\mu_i$) change due to interactions,

Here, the first term is a reference state which can be arbitrarily defined, the second is from the translational partition function (ideal gas), and the third is the excess chemical potential, which is related to the activity coefficient,

where $a_i$ is the activity. Recall that the chemical potential, $\mu_i$, is the free energy associated with adding a single species, $i$ to the system, while keeping all others constant. In practice this is not possible for salts, as they always come as electroneutral pairs. Instead we can define a mean activity coefficient,

valid for salt compound $M_pX_q$. Hence, for NaCl, $\gamma_{\pm} = \sqrt{\gamma_{Na}\gamma_{Cl}}$.

Question:

1. Give an expression for $\gamma_{\pm}$ for sodium sulphate?
2. Give the equivalent expression for the mean excess chemical potential, i.e. the free energy change of adding a single, electroneutral salt molecule. Simplify the expression as much as possible.

Statistical Mechanical Perturbation

Literature: D. McQuirrie, Chapter 14-1. [Full text] and B. Widom, "Some Topics in the Theory of Fluids" [Full text].

Consider the Canonical Ensemble ($NVT$) where the potential energy function can be written as

$U = U^0 + U^1$.

$U^0$ is the energy of some unperturbed reference system, and $U^1$ is a small perturbation that will give us $U$. The configuational integral of the final system is

where the integral runs over all states and $\beta=1/k_BT$. Likewise, the configurational integral of the unperturbed system is $Z^0=\int e^{-\beta U^0}d\mathcal{R}^N$ and multiplying and dividing $Z$ with this we get,

That is, the configurational integral of the final system can be obtained by averaging the exponential of the perturbation using configurations from the unperturbed system (0). The Helmholtz free energy change associated with the perturbation naturally follows,

This important result is often used in MD/MC simulations to estimate free energy changes, for example the excess pressure by performing a virtual volume move or the excess chemical potential by performing a ghost particle insertion. The latter is known as Widom particle insertion - see link above - and we will use this to calculate activity coefficients. As we already covered, the excess chemical potential, $\mu_i^{ex}$, is the free energy due to interactions upon inserting species $i$. During simulation we sample an unperturbed system (0) of $N$ particles and occacionally insert a ghost or virtual particle. Calculating the energy change of this process, $U^1=U-U^0$ we can thus estimate $\mu^{ex}$. After the energy calculation, the particle is removed and it will not influence configurations of the reference system (0) which must remain unperturbed.

Questions:

1. Below you'll find a snapshot of the first page of Widom's article from 1963. Show explicitly how Eq. 1 leads to Eq. 2.

2. Discuss the effectiveness of Widom's method if the system is dense.

B. Widom, 1963:

Part I. Mean Activity Coefficients in Bulk Solution

In this first part we estimate mean activity coefficients, $\gamma_{\pm}$, of salts in aqueous solution using MC simulation in the Grand Canonical ensemble simulations. In particular we will investigate the three salts in the table below and your task is to find the best parameters to describe the salt guanidinium chloride (GdnCl).

Note: The experimental GdnCl data from the Macaskill reference is available in the file exp-gdncl-coeff.csv.

Import required python modules

Download and compile the MC software

This will download and compile the free Monte Carlo software Faunus. It requires CMake and a reasonable modern C++ compiler (clang, gcc4.9+)

Run Monte Carlo simulations

The block below defines the properties of the salts you wish to investigate. If you are happy with the MC results fo a certain salt, set run=False in order to skip simulation runs in the next step (faster).

Tasks:

1. Go through and understand the input below. For more information about Faunus input, see here. In particular, test how to change the number of MC steps, change the displacement parameter, box length.
2. Have a look at the main C++ file mc/excess.cpp and understand the program flow, energy setup etc.

Check MC output

Questions:

1. Investigate how the displacement parameter (dp) affects the translational move acceptance and mean aquare displacement (see output about "Markov Move: Single Particle Translation")

2. Is there a best value for dp?

3. How is the energy drift defined?

4. Explain the difference between the $NPT$ and $\mu VT$ ensembles.

Plot mean activity coefficients

For compound $M_pX_q$, the mean ionic activity coefficient is $ \gamma_{\pm} = \sqrt[q+p]{ \gamma_M^p \gamma_X^q}$. Note that the experimental data is in molality, i.e. moles salt per kilo solvent and should in principle be converted to the molarity scale as outlined in i.e. Robinson and Stokes classic book. We ignore this here as the intended use for the model is for sub-molar concentrations. All salt concentrations are converted to the ionic strength, $I=\frac{1}{2}\sum c_i z_i^2$ where $c$ are molar concentrations, $z$ valencies.

Questions

1. Explain what the meaning of an activity coefficient below, equal to, and above unity.
2. Below we have already tabulated simulation data for NaCl and Na3Cit using radii from the table above. For NaCl, verify that the first point ($c=0.1315$ mol/l) is correct by re-doing the simulation.
3. Adjust the cation and anion radii for guanidinium chloride (GdnCl) to match experiment and add data points to the arrays below.
4. Give a physical explanation for why hard sphere radii alters $\gamma_{\pm}.$ What physical features have we ignored?
5. Discuss and explain the differences between the different salts.
6. Why do the activity coefficients for NaCl tend to increase at high concentrations?

Ion-ion radial distribution function

We now plot the radial distribution function, $g(r)$, between ions as sampled during MC simulation. This analysis is activated via the input in analysis/atomrdf.

Questions:

1. Explain the meaning and definition of $g(r).$

2. Plot and discuss $g(r)$ for cation-anion and cation-anion. How are these influenced by ion radius?

3. Discuss the differences between the three different salts. Focus on the cation-anion interaction.

4. What happens to $g(r)$ as the concentration increases. Look i.e. at NaCl and explain.

5. How is $g(r)$ affected by the box size? What happens at long separation, when $r$ exceeds half the box length?

Widom Particle Insertion

In the above, activity coefficients or excess chemical potentials of ions were calculated simply from the average ion density in the grand canonical ensemble where activities are known and constant. An alternative method is to use Widom particle insertion which we covered in the theory section above. In the MC code this can be done by activating an insertion analysis in the input file. The appropriate input is already there (analysis->widom) but the interval to which to run it, nstep is set to zero. Setting this to a positive integers will activate the sample.

For example, if the number of MC steps, set by macro=10 and micro=1000, then an nstep=50 will lead to 200 calls to the Widom analysis. For each call you may specify how many times to insert a particle, ninsert, which will multiply the number of samples and thus improve statistics. For each insertion, a ghost particle is inserted at a random position and its energy is calculated, whereafter it is removed.

Questions:

1. For four different salt activities in the interval 50-300 mM, plot the activity coefficient of a neutral particle (add data to the arrays below).
2. Repeat the above with a larger radius of the neutral particle.
3. Discuss the results and estimate the so-called Setschenow coefficient which describes how the solubility of gasses is affected by electrolyte.
4. How is the solubility of a greenhouse gas like methane affected by the salt content of the oceans?

Comments

Please leave comments for how we may improve this lab in the future. What was good, what was bad? Did you find errors or unclear sections?