Introductory Tutorial: Lennard-Jones Liquid

Table of Contents

1. Introduction
2. Background
3. The Lennard-Jones Potential
4. Units
5. First steps
6. Overview of a simulation script
   1. System setup
   2. Placing and accessing particles
   3. Setting up non-bonded interactions
   4. Energy minimization
   5. Choosing the thermodynamic ensemble, thermostat
   6. Integrating equations of motion and taking manual measurements
   7. Automated data collection
7. Further Exercises
   1. Binary Lennard-Jones Liquid
8. References

Introduction

Welcome to the basic ESPResSo tutorial!

In this tutorial, you will learn, how to use the ESPResSo package for your research. We will cover the basics of ESPResSo, i.e., how to set up and modify a physical system, how to run a simulation, and how to load, save and analyze the produced simulation data.

More advanced features and algorithms available in the ESPResSo package are described in additional tutorials.

Background

Research on Soft Condensed Matter has brought the needs for having a flexible, extensible, reliable, and efficient (parallel) molecular simulation package. For this reason ESPResSo (Extensible Simulation Package for Research on Soft Matter Systems) [1] has been developed at the Max Planck Institute for Polymer Research, Mainz, and at the Institute for Computational Physics at the University of Stuttgart in the group of Prof. Dr. Christian Holm [2,3]. The ESPResSo package is a flexible and extensible simulation package. It is specifically developed for coarse-grained molecular dynamics (MD) simulation of polyelectrolytes but is not necessarily limited to this. For example, it could also be used to simulate granular media. ESPResSo has been nominated for the Heinz-Billing-Preis for Scientific Computing in 2003 [4].

Prior knowledge

• Basics on thermodynamics and statistical mechanics.
• Basics on coarse-grained simulation methods: molecular dynamics and Langevin dynamics.
• Correlation functions.
• Pair structure and radial distribution function.

The Lennard-Jones Potential

A pair of neutral atoms or molecules is subject to two distinct forces in the limit of large separation and small separation: an attractive force at long ranges (van der Waals force, or dispersion force) and a repulsive force at short ranges (the result of overlapping electron orbitals, referred to as Pauli repulsion from the Pauli exclusion principle). The Lennard-Jones potential (also referred to as the LJ potential, 6-12 potential or, less commonly, 12-6 potential) is a simple mathematical model that represents this behavior. It was proposed in 1924 by John Lennard-Jones. The LJ potential is of the form

\begin{equation} V(r) = 4\epsilon \left[ \left( \dfrac{\sigma}{r} \right)^{12} - \left( \dfrac{\sigma}{r} \right)^{6} \right] \end{equation}

where $\epsilon$ is the depth of the potential well, providing a measure of the attraction strength, $\sigma$ is the (finite) distance at which the inter-particle potential is zero and $r$ is the distance between the particles. The $\left(\frac{1}{r}\right)^{12}$ term describes repulsion and the $\left(\frac{1}{r}\right)^{6}$ term describes attraction. The Lennard-Jones potential is an approximation. The form of the repulsion term has no theoretical justification; the repulsion force should depend exponentially on the distance, but the inverse power form of the repulsion term of the LJ formula is more convenient due to the ease and efficiency of computing $r^{12}$ as the square of $r^6$.

In simulations, the LJ potential is typically cut off beyond a specified distance $r_{\mathrm{cut}}$, i.e. the potential is zero for distances larger than the cutoff distance. Also an additional potential shift $c_{\mathrm{shift}}$ can be added:

\begin{equation} V(r) = \begin{cases}4\epsilon \left[ \left( \dfrac{\sigma}{r} \right)^{12} - \left( \dfrac{\sigma}{r} \right)^{6} \right] + c_{\mathrm{shift}} && 0<r\leq r_{\mathrm{cut}} \\ 0 && r > r_{\mathrm{cut}} \end{cases} \end{equation}

For example, the LJ potential could be cut off at $r_{\mathrm{cut}}=2^{1/6}\sigma$, which corresponds to the location of the potential minimum. One can then shift the potential up by $c_{\mathrm{shift}}=\epsilon$ to make it continuous, resulting in a purely repulsive version of the LJ potential, commonly known as the Weeks–Chandler–Andersen potential. Another common choice for the cutoff is $r_{\mathrm{cut}}=2.5\sigma$, beyond which the potential is consistently $|V(r>r_{\mathrm{cut}})| < 0.02 \epsilon$. In this case, $c_{\mathrm{shift}}$ is chosen such that the potential is continuous at $r_{\mathrm{cut}}$. The main advantage of this truncated version of the LJ potential is that it is significantly cheaper computationally, while the error due to the missing tail is negligible.

Originally LJ has been introduced to describe Van der Waals interactions in molecular systems, but nowadays it is often used to represent generic particles in coarse-grained simulations, such as polymer systems.

import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt

plt.rcParams.update({'font.size': 18})

def lj_pot(x, epsilon, sigma, r_cut, c_shift=0.0):
    pot = 4.0 * epsilon * ((sigma / x)**12 - (sigma / x)**6) + c_shift
    pot[x > r_cut] = 0.
    return pot

epsilon = 1.0  # energy in units of k_BT
sigma = 1.0    # distance in units of sigma

xs = np.linspace(0.5, 3, 100)
ys_lj = lj_pot(xs, epsilon, sigma, r_cut=xs[-1])
ys_WCA = lj_pot(xs, epsilon, sigma, r_cut=2**(1/6) * sigma, c_shift=epsilon)

fig = plt.figure(figsize=(10, 6))
plt.plot(xs, ys_lj, label='LJ')
plt.plot(xs, ys_WCA, label='WCA')
plt.axhline(y=0, color='grey')
plt.xlabel(r"$r/\sigma$")
plt.ylabel(r"$V(r)/(k_{\mathrm{B}}T)$")
plt.legend()
plt.ylim(-1.5, 2.5)
plt.show()

Units

Novice users must understand that ESPResSo has no fixed unit system. The unit system is set by the user. Conventionally, reduced units are employed. For instance, the potential can be expressed in units of $k_{\mathrm{B}}T$: $$V^*(r^*)=\beta V(r^*) = 4\epsilon^* \left[ \left( \dfrac{1}{r^*} \right)^{12} - \left( \dfrac{1}{r^*} \right)^{6} \right]$$ with $\beta=(k_{\mathrm{B}}T)^{-1}$, $\epsilon^*=\beta\epsilon$ and $r^*=r/\sigma$. We have used $\sigma$ as unit of distance.

For more information on units, please check the documentation on units.

Lennard-Jones fluid

In the following simulations, we are going to model a single-species Lennard-Jones fluid at constant number of particles $N$, density $\rho=\frac{N}{V}$ and temperature $T$. Depending on the chosen parameters, the system might present different phase behavior, leading to solid, liquid and gas phases. Simulations will allow us to characterize these different phases structurally by means of the radial distribution function.

For a more detailed discussion on simulation results about phase diagrams of Lennard-Jones fluids, see [5].

First steps

What is ESPResSo? It is an extensible, efficient Molecular Dynamics package specially powerful on simulating charged systems. In depth information about the package can be found in the relevant sources [1,4,2,3].

ESPResSo consists of two components. The simulation engine is written in C++ for the sake of computational efficiency. The steering or control level is interfaced to the kernel via an interpreter of the Python scripting languages.

The kernel performs all computationally demanding tasks. Before all, integration of Newton's equations of motion, including calculation of energies and forces. It also takes care of internal organization of data, storing the data about particles, communication between different processors or cells of the cell-system.

The scripting interface (Python) is used to setup the system (particles, boundary conditions, interactions etc.), control the simulation, run analysis, and store and load results. The user has at hand the full reliability and functionality of the scripting language. For instance, it is possible to use the SciPy package for analysis and PyPlot for plotting. With a certain overhead in efficiency, it can also be used to reject/accept new configurations in combined MD/MC schemes. In principle, any parameter which is accessible from the scripting level can be changed at any moment of runtime. In this way methods like thermodynamic integration become readily accessible.

Note: This tutorial assumes that you already have a working ESPResSo installation on your system. If this is not the case, please consult the first chapters of the user's guide for installation instructions.

Overview of a simulation script

Typically, a simulation script consists of the following parts:

• System setup (box geometry, thermodynamic ensemble, integrator parameters)
• Placing the particles
• Setup of interactions between particles
• Equilibration (bringing the system into a state suitable for measurements)
• Integration loop (propagate the system in time and record measurements)

System setup

The functionality of ESPResSo for python is provided via a python module called espressomd. At the beginning of the simulation script, it has to be imported.

import espressomd
import espressomd.observables
import espressomd.accumulators
import espressomd.analyze
required_features = ["LENNARD_JONES"]
espressomd.assert_features(required_features)

The function espressomd.assert_features() expects a list of features as argument and checks they are available in the ESPResSo executable. If a required feature is missing, the program will print an error message and halt. To compile ESPResSo with a different set of features, see the documentation on features.

import scipy.optimize
np.random.seed(42)

# System parameters
N_PART = 200
DENSITY = 0.75

BOX_L = np.power(N_PART / DENSITY, 1.0 / 3.0) * np.ones(3)

The next step would be to create an instance of the System class. This instance is used as a handle to the simulation system. At any time, only one instance of the System class can exist.

Exercise:

• Create an instance of an espresso system and store it in a variable called system; use BOX_L as box length.

See ESPResSo documentation and module documentation.

system = espressomd.System(box_l=BOX_L)

# Test solution of Exercise 1
assert isinstance(system, espressomd.System)

It can be used to store and manipulate the crucial system parameters like the time step and the size of the simulation box (time_step, and box_l).

SKIN = 0.4
TIME_STEP = 0.01

system.time_step = TIME_STEP
system.cell_system.skin = SKIN

The parameter SKIN affects how often the Verlet lists will be updated. This parameter does not influence the physics of the simulation. It can however have a significant impact on the performance of the simulation. Depending on system parameters such as density and temperature, the optimal SKIN parameter can vary considerably. Please be aware that ESPResSo implements the function tune_skin() that automatically tunes SKIN for optimal performance.

Placing and accessing particles

Particles in the simulation can be added and accessed via the part property of the System class. Individual particles should be referred to by the particle handle returned upon creation. You can also retrieve them by an integer id, e.g. system.part.by_id(0). If id is unspecified when creating a particle, an unused particle id is automatically assigned. It is also possible to use common python iterators and slicing operations to add or access several particles at once.

Particles can be grouped into several types, so that, e.g., a binary fluid can be simulated. Particle types are identified by integer ids, which are set via the particles' type attribute. If it is not specified, zero is implied.

• Create N_PART particles at random positions, store their handles in a variable called particles.

  Use system.part.add().

  Use an (N_PART x 3) array for positions. Use np.random.random() to generate random numbers.

particles = system.part.add(type=[0] * N_PART, pos=np.random.random((N_PART, 3)) * system.box_l)

# Test that now we have indeed N_PART particles in the system
assert len(system.part) == N_PART

The particle properties can be accessed using standard numpy slicing syntax:

# Access position of a single particle
print(f"position of particle with id 0: {system.part.by_id(0).pos}")

# Iterate over the first five particles for the purpose of demonstration.
first_five = system.part.by_ids(range(5))
for p in first_five:
    print(f"id {p.id} position: {p.pos}")
    print(f"id {p.id} velocity: {p.v}")

# Obtain particle positions for the particles created until now
cur_pos = particles.pos

position of particle with id 0: [2.41076339 6.1193638  4.7115492 ]
id 0 position: [2.41076339 6.1193638  4.7115492 ]
id 0 velocity: [0. 0. 0.]
id 1 position: [3.85332274 1.00422894 1.00407369]
id 1 velocity: [0. 0. 0.]
id 2 position: [0.37386074 5.57522583 3.86913442]
id 2 velocity: [0. 0. 0.]
id 3 position: [4.55757705 0.13249407 6.24291778]
id 3 velocity: [0. 0. 0.]
id 4 position: [5.35809689 1.36674105 1.17033384]
id 4 velocity: [0. 0. 0.]

You can also get all particles using system.part.all(), but particles already contains all particles that are in the simulation so far.

Many objects in ESPResSo have a string representation, and thus can be displayed via python's print function:

print(system.part.by_id(0))

ParticleHandle([('id', 0), ('pos', (2.4107633923737293, 6.119363804209853, 4.711549202763617)), ('_id', 0), ('bonds', ()), ('dip', (0.0, 0.0, 0.0)), ('dipm', 0.0), ('director', (0.0, 0.0, 1.0)), ('exclusions', ()), ('ext_force', (0.0, 0.0, 0.0)), ('ext_torque', (0.0, 0.0, 0.0)), ('f', (0.0, 0.0, 0.0)), ('fix', (0, 0, 0)), ('gamma', (-1.0, -1.0, -1.0)), ('gamma_rot', (-1.0, -1.0, -1.0)), ('image_box', (0, 0, 0)), ('lees_edwards_flag', 0), ('lees_edwards_offset', 0.0), ('mass', 1.0), ('mol_id', 0), ('node', 0), ('omega_body', (0.0, 0.0, 0.0)), ('omega_lab', (0.0, 0.0, 0.0)), ('q', 0.0), ('quat', (1.0, 0.0, 0.0, 0.0)), ('rinertia', (1.0, 1.0, 1.0)), ('rotation', (0, 0, 0)), ('swimming', {'v_swim': 0.0, 'f_swim': 0.0, 'mode': 'N/A', 'dipole_length': 0.0}), ('torque_lab', (0.0, 0.0, 0.0)), ('type', 0), ('v', (0.0, 0.0, 0.0)), ('virtual', False), ('vs_quat', (1.0, 0.0, 0.0, 0.0)), ('vs_relative', (-1, 0.0, array([1., 0., 0., 0.])))])

Setting up non-bonded interactions

Non-bonded interactions act between all particles of a given combination of particle types. In this tutorial, we use the Lennard-Jones non-bonded interaction. First we define the LJ parameters

# use LJ units: EPS=SIG=1
LJ_EPS = 1.0
LJ_SIG = 1.0
LJ_CUT = 2.5 * LJ_SIG

In a periodic system, it is in general not straightforward to calculate all non-bonded interactions. As mentioned earlier in the text, usually a cutoff distance $r_{\mathrm{cut}}$ is applied for infinite-range potentials like Lennard-Jones, such that $V(r>r_{\mathrm{cut}}) = 0$. The potential can be shifted to zero at the cutoff value to ensure continuity using the shift='auto' option of espressomd.interactions.LennardJonesInteraction. For comparison with the fundamental work on MD simulations of LJ systems [7], we do NOT shift the potential to zero at the cutoff in this tutorial and instead correct for the long-range error $V_{\mathrm{lr}}$ later. However, it is strongly recommended to implement shift='auto' in any other case to ensure continuity of the potential!

To avoid spurious self-interactions of particles with their periodic images one usually forces that the shortest box length is at least twice the cutoff distance:

assert (BOX_L - 2 * SKIN > LJ_CUT).all()

Exercise:

• Setup a Lennard-Jones interaction with $\epsilon=$LJ_EPS and $\sigma=$LJ_SIG that is cut at $r_{\mathrm{cut}}=$LJ_CUT$\times\sigma$ and not shifted.

Hint:

• Have a look at the documentation on non-bonded interactions

system.non_bonded_inter[0, 0].lennard_jones.set_params(
    epsilon=LJ_EPS, sigma=LJ_SIG, cutoff=LJ_CUT, shift=0)

Energy minimization

In many cases, including this tutorial, particles are initially placed randomly in the simulation box. It is therefore possible that particles overlap, resulting in a huge repulsive force between them. In this case, integrating the equations of motion would not be numerically stable. Hence, it is necessary to remove possible overlaps. This is typically done by performing a steepest descent minimization of the potential energy until a maximal force criterion is reached.

Note: Making sure a system is properly equilibrated highly depends on the system details. In most cases a relative convergence criterion on the forces and/or energies works well but you might have to make sure that the total force is smaller than a threshold value f_max at the end of the minimization. Depending on the simulated system, other strategies might be necessary, such as simulating with a small time step or with capped forces.

# suitable minimization parameters for this LJ system
F_TOL = 1e-2
DAMPING = 30
MAX_STEPS = 10000
MAX_DISPLACEMENT = 0.01 * LJ_SIG
EM_STEP = 10

Exercise:

• Use espressomd.integrate.set_steepest_descent to relax the initial configuration. Use a maximal displacement $\vec{r}_{\mathrm{max}}$ of MAX_DISPLACEMENT. The particle displacement is related to the particle force via a damping constant $\gamma$, such that $\vec{r}(t + dt) = \vec{r}(t) + \min(\gamma \vec{F}(t), \vec{r}_{\mathrm{max}})$. Use a damping constant gamma = DAMPING.
• Use the relative change of the system maximal force as a convergence criterion. Check the documentation espressomd.particle_data module to obtain the forces. The steepest descent has converged if the relative force change between two rounds of minimizations is less than the threshold value F_TOL. Note that by default espressomd.integrate.set_steepest_descent will halt when the system maximal force is less than some value f_max. When a custom convergence criterion is implemented, as it is the case here, the default convergence criterion needs to be disabled by setting f_max=0.
• Break the minimization loop after a maximal number of MAX_STEPS steps or if convergence is achieved. Check for convergence every EMSTEP steps.

Hint: To obtain the initial forces one has to initialize the integrator using integ_steps=0, i.e. call system.integrator.run(0) before accessing the force array.

# Set up steepest descent integration
system.integrator.set_steepest_descent(f_max=0,  # use a relative convergence criterion only
                                       gamma=DAMPING,
                                       max_displacement=MAX_DISPLACEMENT)

# Initialize integrator to obtain initial forces
system.integrator.run(0)
old_force = np.max(np.linalg.norm(system.part.all().f, axis=1))


while system.time / system.time_step < MAX_STEPS:
    system.integrator.run(EM_STEP)
    force = np.max(np.linalg.norm(system.part.all().f, axis=1))
    rel_force = np.abs((force - old_force) / old_force)
    print(f'rel. force change: {rel_force:.2e}')
    if rel_force < F_TOL:
        break
    old_force = force

rel. force change: 1.00e+00
rel. force change: 9.99e-01
rel. force change: 9.38e-01
rel. force change: 6.00e-01
rel. force change: 6.81e-01
rel. force change: 5.53e-01
rel. force change: 3.49e-01
rel. force change: 2.62e-01
rel. force change: 2.26e-03

# reset clock
system.time = 0.

Choosing the thermodynamic ensemble, thermostat

Simulations can be carried out in different thermodynamic ensembles such as NVE (particle Number, Volume, Energy), NVT (particle Number, Volume, Temperature) or isotropic NpT (particle Number, pressure, Temperature).

In this tutorial, we use a NVT ensemble with the Langevin thermostat.

# Parameters for the Langevin thermostat
# reduced temperature T* = k_B T / LJ_EPS
TEMPERATURE = 0.827  # value from Tab. 1 in [7]
GAMMA = 1.0

Exercise:

• Use system.integrator.set_vv() to use a Velocity Verlet integration scheme and system.thermostat.set_langevin() to turn on the Langevin thermostat. Set the reduced temperature to TEMPERATURE and the solvent friction coefficient to GAMMA. Initialize the random number generator of the thermostat with an integer seed.

The Langevin thermostat maintains the system temperature at a constant value by mimicking the effect of an implicit solvent with large heat capacity that acts as a heat reservoir. This is achieved by introducing both a friction term and a stochastic term in the Newton's equations of motion (for details see the documentation on thermostats).

system.integrator.set_vv()
system.thermostat.set_langevin(kT=TEMPERATURE, gamma=GAMMA, seed=42)

Integrating equations of motion and taking manual measurements

Now, we integrate the equations of motion and take measurements of relevant quantities.

# Integration parameters
STEPS_PER_SAMPLE = 20
N_SAMPLES = 1000

times = np.zeros(N_SAMPLES)
e_total = np.zeros_like(times)
e_kin = np.zeros_like(times)
T_inst = np.zeros_like(times)

Exercise:

• Integrate the system and measure the total and kinetic energy. Take N_SAMPLES measurements every STEPS_PER_SAMPLE integration steps. Notice that the total simulated time in LJ units is given by the product N_SAMPLES $\times$ STEPS_PER_SAMPLE $\times$ TIME_STEP.
• Calculate the total and kinetic energies using the analysis method system.analysis.energy().
• Use the containers times, e_total and e_kin from the cell above to store the time series.
• From the simulation results, calculate the kinetic temperature $T_{\mathrm{inst}} = 2/3 \times E_{\mathrm{kin}}$ / N_PART.

for i in range(N_SAMPLES):
    times[i] = system.time
    energy = system.analysis.energy()
    e_total[i] = energy['total']
    e_kin[i] = energy['kinetic']
    system.integrator.run(STEPS_PER_SAMPLE)
T_inst = 2. / 3. * e_kin / N_PART

plt.figure(figsize=(10, 6))
plt.plot(times, T_inst, label='$T_{\\mathrm{inst}}$')
plt.plot(times, [TEMPERATURE] * len(times), label='$T$ set by thermostat')
plt.legend()
plt.xlabel('t')
plt.ylabel('T')
plt.show()

Since the ensemble average $\langle E_{\mathrm{kin}}\rangle=3/2 N k_B T$ is related to the temperature, we may compute the actual temperature of the system via $k_B T= 2/(3N) \langle E_{\mathrm{kin}}\rangle$. Note the thermodynamic temperature $T$ is fixed and does not fluctuate in the NVT ensemble! The kinetic temperature calculated via $T_{\mathrm{inst}} = 2/(3N) E_{\mathrm{kin}}$ (without ensemble averaging) is allowed to fluctuate, but it is not the temperature of the system ($T = \langle T_{\mathrm{inst}} \rangle$).

In the first simulation run, we picked STEPS_PER_SAMPLE arbitrarily. To ensure proper statistics, we will subsample the time series to reduce correlation between consecutive measurements. In order to do this, we first have to remove the beginning of the time series, because the system was out of equilibrium. The time to reach equilibrium depends on the thermostat friction coefficient gamma and can be determined visually from the plot.

# Use only the data where the system is at equilibrium
equilibration_time = 15
e_total = e_total[times > equilibration_time]
e_kin = e_kin[times > equilibration_time]
times = times[times > equilibration_time]
times -= times[0]

Notice that equilibration_time is not known a priori. It has to be found by trial and error for the given set of input parameters.

Autocorrelation function and correlation time

In simulations, the knowledge of the time evolution of the particle positions and velocities allows for the calculation of autocorrelation functions and autocorrelation time.

Correlation functions are useful to quantify how microscopic observables at different positions or times of a system are related one another. One particular example is the equilibrium velocity autocorrelation function $$C(\tau)=\langle {\bf v}(0)\cdot{\bf v}(\tau)\rangle,$$ that relates the velocities at different times along an equilibrium trajectory [8]. Time correlation functions are typically decaying functions in fluid systems, whose decay is characterized by the so-called correlation time $\xi$. We say that two samples of a given system taken at time $t'$ and $t''$, respectively, are uncorrelated for a certain observable $X$, if $t'-t''\gg \xi_X $.

The ESPResSo module espressomd.analyze provides a function to compute the autocorrelation of time series. By means of espressomd.analyze.autocorrelation(X), it is possible to calculate the unnormalized autocorrelation function of an observable $X$ measured at time $t$ with constant time step for lag times $\tau$. This method neither subtracts the mean value nor normalizes by the variance of the provided time series.

def autocor(x):
    x = np.asarray(x)
    mean = np.mean(x)
    var = np.var(x)
    xp = x - mean
    corr = espressomd.analyze.autocorrelation(xp) / var
    return corr


def fit_correlation_time(data, ts):
    data = np.asarray(data)
    data /= data[0]

    def fitfn(t, t_corr): return np.exp(-t / t_corr)
    popt, pcov = scipy.optimize.curve_fit(fitfn, ts, data)
    return popt[0]

Exercise

• Calculate the autocorrelation of the total energy (store in e_total_autocor). Calculate the correlation time $\xi$ (corr_time). Calculate a quantity steps_per_subsample that represents the number of integration steps necessary to advance the simulation time by $3\xi$. This is a conservative quantity that will help us subsample the time series in such a way that the correlation between two consecutive samples is small ($e^{-3} \simeq 5\%$).

The value steps_per_subsample is somewhat arbitrary and was chosen as a trade-off between accuracy and simulation runtime for the purpose of this tutorial. In a research project, you would run simulations for a lot longer and increase the time between subsamples to reduce the residual correlation even further. Please refer to the Error Analysis tutorial for an in-depth discussion of time series autocorrelation.

e_total_autocor = autocor(e_total)
corr_time = fit_correlation_time(e_total_autocor[:100], times[:100])
steps_per_subsample = int(np.ceil(3 * corr_time / system.time_step))

print(f'steps_per_subsample = {steps_per_subsample}')

steps_per_subsample = 201

We plot the autocorrelation function and the fit to visually confirm a roughly exponential decay

plt.figure(figsize=(10, 6))
plt.plot(times[1:], e_total_autocor[1:], label='data')
plt.plot(times[1:], np.exp(-times[1:] / corr_time), label='exponential fit')
plt.plot(2 * [steps_per_subsample * system.time_step],
         [min(e_total_autocor), 1], label='subsampling interval')
plt.ylim(top=1.2, bottom=-0.15)
plt.legend()
plt.xscale('log')
plt.xlabel('t')
plt.ylabel('total energy autocorrelation')
plt.show()

In order to obtain equilibrium properties, we need to consider ensemble-averaged quantities. Assuming that the simulation is an ergodic process, averaging over uncorrelated samples would provide equilibrium results.

Exercise:

• Calculate the mean and standard error of the mean potential energy per particle using the formula for uncorrelated samples (define mean_pot_energy and SEM_pot_energy).

Hint

• You know how many steps are between samples in e_total and how many steps are between subsamples. So you have to figure out how many samples to skip.

subsample_step = int(np.ceil(steps_per_subsample / STEPS_PER_SAMPLE))
pot_energies = (e_total - e_kin)[::subsample_step] / N_PART
mean_pot_energy = np.mean(pot_energies)
SEM_pot_energy = np.std(pot_energies) / np.sqrt(len(pot_energies))

print(f'mean potential energy = {mean_pot_energy:.2f} +/- {SEM_pot_energy:.2f}')

mean potential energy = -4.96 +/- 0.01

For comparison to literature values we need to account for the error made by the LJ truncation. For an isotropic system one can assume that the density is homogeneous behind the cutoff, which allows to calculate the so-called long-range corrections to the energy and pressure: $$V_{\mathrm{lr}} = \frac{1}{2} \rho \int_{r_{\mathrm{cut}}}^\infty 4 \pi r^2 g(r) V(r) \,\mathrm{d}r$$ Using that the radial distribution function $g(r)=1$ for $r>r_{\mathrm{cut}}$ one obtains $$V_{\mathrm{lr}} = -\frac{8}{3}\pi \rho \epsilon \sigma^3 \left[\frac{1}{3} \left(\frac{\sigma}{r_{\mathrm{cut}}}\right)^9 - \left(\frac{\sigma}{r_{\mathrm{cut}}}\right)^3 \right].$$ Similarly, a long-range contribution to the pressure can be derived [6].

tail_energy_per_particle = 8. / 3. * np.pi * DENSITY * LJ_EPS * \
    LJ_SIG**3 * (1. / 3. * (LJ_SIG / LJ_CUT)**9 - (LJ_SIG / LJ_CUT)**3)
mean_pot_energy_corrected = mean_pot_energy + tail_energy_per_particle
print(f'corrected mean potential energy = {mean_pot_energy_corrected:.2f}')

corrected mean potential energy = -5.36

This value differs quite strongly from the uncorrected one but agrees well with the literature value $U^i = -5.38$ given in Table 1 of Ref. [7].

Automated data collection

As we have seen, it is easy to manually extract information from an ESPResSo simulation, but it can get quite tedious. Therefore, ESPResSo provides a number of data collection tools to make life easier (and less error-prone). We will now demonstrate those with the calculation of the radial distribution function.

Observables extract properties from the particles and calculate some quantity with it, e.g. the center of mass, the total energy or a histogram.

Accumulators allow the calculation of observables while running the system and then doing further analysis. Examples are a simple time series or more advanced methods like correlators. For our purposes we need an accumulator that calculates the average of the RDF samples.

# Parameters for the radial distribution function
BIN_WIDTH = 0.05
R_MIN = 0.0
R_MAX = system.box_l[0] / 2.0
N_BINS = int((R_MAX-R_MIN)/BIN_WIDTH)

Exercise

• Instantiate a RDF observable
• Instantiate a MeanVarianceCalculator accumulator to track the RDF over time. Samples should be taken every steps_per_subsample steps.
• Add the accumulator to the auto_update_accumulators of the system for automatic updates

rdf_obs = espressomd.observables.RDF(ids1=system.part.all().id, min_r=R_MIN, max_r=R_MAX, n_r_bins=N_BINS)
rdf_acc = espressomd.accumulators.MeanVarianceCalculator(obs=rdf_obs, delta_N=steps_per_subsample)
system.auto_update_accumulators.add(rdf_acc)

Now we don't need an elaborate integration loop anymore, instead the RDFs are calculated and accumulated automatically.

system.integrator.run(N_SAMPLES * steps_per_subsample)

Exercise

• Get the mean RDF (define rdf) from the accumulator
• Get the histogram bin centers (define rs) from the observable

rdf = rdf_acc.mean()
rs = rdf_obs.bin_centers()

fig, ax = plt.subplots(figsize=(10, 7))
ax.plot(rs, rdf, label='simulated')
plt.legend()
plt.xlabel('r')
plt.ylabel('RDF')
plt.show()

We now plot the experimental radial distribution. Empirical radial distribution functions have been determined for pure fluids [9], mixtures [10] and confined fluids [11]. We will compare our distribution $g(r)$ to the theoretical distribution $g(r^*, \rho^*, T^*)$ of a pure fluid [9].

# comparison to literature
def calc_literature_rdf(rs, temperature, density, LJ_eps, LJ_sig):
    T_star = temperature / LJ_eps
    rho_star = density * LJ_sig**3

    # expression of the factors Pi from Equations 2-8 with coefficients qi from Table 1
    # expression for a,g
    def P(q1, q2, q3, q4, q5, q6, q7, q8, q9): return \
        q1 + q2 * np.exp(-q3 * T_star) + q4 * np.exp(-q5 * T_star) + q6 / rho_star + q7 / rho_star**2 \
        + q8 * np.exp(-q3 * T_star) / rho_star**3 + q9 * \
        np.exp(-q5 * T_star) / rho_star**4
    a = P(9.24792, -2.64281, 0.133386, -1.35932, 1.25338,
          0.45602, -0.326422, 0.045708, -0.0287681)
    g = P(0.663161, -0.243089, 1.24749, -2.059, 0.04261,
          1.65041, -0.343652, -0.037698, 0.008899)

    # expression for c,k
    def P(q1, q2, q3, q4, q5, q6, q7, q8): return \
        q1 + q2 * np.exp(-q3 * T_star) + q4 * rho_star + q5 * rho_star**2 + q6 * rho_star**3 \
        + q7 * rho_star**4 + q8 * rho_star**5
    c = P(-0.0677912, -1.39505, 0.512625, 36.9323, -
          36.8061, 21.7353, -7.76671, 1.36342)
    k = P(16.4821, -0.300612, 0.0937844, -61.744,
          145.285, -168.087, 98.2181, -23.0583)

    # expression for b,h
    def P(q1, q2, q3): return q1 + q2 * np.exp(-q3 * rho_star)
    b = P(-8.33289, 2.1714, 1.00063)
    h = P(0.0325039, -1.28792, 2.5487)

    # expression for d,l
    def P(q1, q2, q3, q4): return q1 + q2 * \
        np.exp(-q3 * rho_star) + q4 * rho_star
    d = P(-26.1615, 27.4846, 1.68124, 6.74296)
    l = P(-6.7293, -59.5002, 10.2466, -0.43596)

    # expression for s
    def P(q1, q2, q3, q4, q5, q6, q7, q8): return \
        (q1 + q2 * rho_star + q3 / T_star + q4 / T_star**2 + q5 / T_star**3) \
        / (q6 + q7 * rho_star + q8 * rho_star**2)
    s = P(1.25225, -1.0179, 0.358564, -0.18533,
          0.0482119, 1.27592, -1.78785, 0.634741)

    # expression for m
    def P(q1, q2, q3, q4, q5, q6): return \
        q1 + q2 * np.exp(-q3 * T_star) + q4 / T_star + \
        q5 * rho_star + q6 * rho_star**2
    m = P(-5.668, -3.62671, 0.680654, 0.294481, 0.186395, -0.286954)

    # expression for n
    def P(q1, q2, q3): return q1 + q2 * np.exp(-q3 * T_star)
    n = P(6.01325, 3.84098, 0.60793)

    # fitted expression (=theoretical curve)
    # slightly more than 1 to smooth out the discontinuity in the range [1.0, 1.02]
    theo_rdf_cutoff = 1.02

    theo_rdf = 1 + 1 / rs**2 * (np.exp(-(a * rs + b)) * np.sin(c * rs + d)
                                + np.exp(-(g * rs + h)) * np.cos(k * rs + l))
    theo_rdf[np.nonzero(rs <= theo_rdf_cutoff)] = \
        s * np.exp(-(m * rs + n)**4)[np.nonzero(rs <= theo_rdf_cutoff)]
    return theo_rdf

theo_rdf = calc_literature_rdf(rs, TEMPERATURE, DENSITY, LJ_EPS, LJ_SIG)

ax.plot(rs, theo_rdf, label='literature')
ax.legend()
fig

Further Exercises

Binary Lennard-Jones Liquid

A two-component Lennard-Jones liquid can be simulated by placing particles of two types (0 and 1) into the system. Depending on the Lennard-Jones parameters, the two components either mix or separate.

1. Modify the code such that half of the particles are of type=1. Type 0 is implied for the remaining particles.
2. Specify Lennard-Jones interactions between type 0 particles with other type 0 particles, type 1 particles with other type 1 particles, and type 0 particles with type 1 particles (set parameters for system.non_bonded_inter[i,j].lennard_jones where {i,j} can be {0,0}, {1,1}, and {0,1}. Use the same Lennard-Jones parameters for interactions within a component, but use a different lj_cut_mixed parameter for the cutoff of the Lennard-Jones interaction between particles of type 0 and particles of type 1. Set this parameter to $2^{\frac{1}{6}}\sigma$ to get de-mixing or to $2.5\sigma$ to get mixing between the two components.
3. Record the radial distribution functions separately for particles of type 0 around particles of type 0, type 1 around particles of type 1, and type 0 around particles of type 1. This can be done by changing the ids1/ids2 arguments of the espressomd.observables.RDF command. You can record all three radial distribution functions in a single simulation. It is also possible to write them as several columns into a single file.
4. Plot the radial distribution functions for all three combinations of particle types. The mixed case will differ significantly, depending on your choice of lj_cut_mixed. Explain these differences.

References

[1] https://espressomd.org
[2] H. J. Limbach, A. Arnold, B. Mann and C. Holm. ESPResSo: An extensible simulation package for research on soft matter systems. Computer Physics Communications, 174(9):704–727, 2006. DOI:10.1016/j.cpc.2005.10.005
[3] A. Arnold, O. Lenz, S. Kesselheim, R. Weeber, F. Fahrenberger, D. Rohm, P. Košovan, and C. Holm. ESPResSo 3.1 — molecular dynamics software for coarse-grained models. In M. Griebel and M. A. Schweitzer, editors, Meshfree Methods for Partial Differential Equations VI, volume 89 of Lecture Notes in Computational Science and Engineering, pages 1–23. Springer Berlin Heidelberg, 2013. DOI:10.1007/978-3-642-32979-1_1
[4] A. Arnold, BA Mann, HJ Limbach, and C. Holm. ESPResSo–An Extensible Simulation Package for Research on Soft Matter Systems. In K. Kremer and V. Macho, editors, Forschung und wissenschaftliches Rechnen, Beiträge zum Heinz-Billing-Preis 2003, volume 63, pages 43–59. Gesellschaft für wissenschaftliche Datenverarbeitung mbH Göttingen, 2004. https://www.gwdg.de/documents/20182/31188/gwdg-bericht-63.pdf
[5] B. Smit. Phase diagrams of Lennard‐Jones fluids. Journal of Chemical Physics, 96:8639, 1992. DOI:10.1063/1.462271
[6] W. P. Allen and D. J. Tildesley. Computer Simulation of Liquids. Oxford University Press, 2017.
[7] L. Verlet, Computer “Experiments” on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules, Physical Review, 159(1):98–103, 1967. DOI:10.1103/PhysRev.159.98
[8] D. Frenkel and B. Smit. Understanding Molecular Simulation: From Algorithms to Applications. Second edition. Academic Press, 2002.
[9] Morsali, Goharshadi, Mansoori, Abbaspour. An accurate expression for radial distribution function of the Lennard-Jones fluid. Chemical Physics, 310(1–3):11–15, 2005. DOI:10.1016/j.chemphys.2004.09.027
[10] Matteoli. A simple expression for radial distribution functions of pure fluids and mixtures. The Journal of Chemical Physics, 103(11):4672, 1995. DOI:10.1063/1.470654
[11] Abbaspour, Akbarzadeha, Abroodia. A new and accurate expression for the radial distribution function of confined Lennard-Jones fluid in carbon nanotubes. RSC Advances, 5(116): 95781–95787, 2015. DOI:10.1039/C5RA16151G