# Lennard-Jones fluid - thermodynamic state variables

We consider a simple fluid of particles that interact due to the Lennard Jones potential

$$U(r) = 4\varepsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]$$

To validate the integration of the diffusion process and correct implementation of energies and forces, we compare observables to other known results. The observables are: the mean potential energy of the system per particle $U$, the pressure $P$ and the radial distribution function $g(r)$.

The thermodynamic quantities of a Lennard Jones fluid are typically expressed in terms of the rescaled density $\rho^* = (N/V)\sigma^3$ and the rescaled temperature $T^* = k_BT/\varepsilon$, where $N$ is the number of particles in the system, which is constant over time, and $V$ is the available volume. In simulation practice we set $\sigma=1$ and $\varepsilon=1$ to achieve the reduced units. The time units are $\sigma^2 / 4D$, where $D$ is the self diffusion coefficient of the particles. Hence in practice we set the diffusion coefficient to 1 and have a fully rescaled unit system.

We use an Euler scheme to integrate the positions of particles according to the overdamped Langevin equation of motion, in contrast to a full integration of positions and momenta as in the underdamped limit.

The pressure can be measured from the acting forces according to [4].

$$PV = Nk_BT + \langle \mathscr{W} \rangle$$

where

$$\mathscr{W} = \frac{1}{3} \sum_i \sum_{j>i} \mathbf{r}_{ij} \mathbf{f}_{ij},$$

This is implemented by ReaDDy's pressure observable.

Results

origin cutoff radius $r_c$ density $\rho$ temperature $T$ pressure $P$ potential energy per particle $U$
[1] 4 0.3 3 1.023(2) -1.673(2)
[2] 4 0.3 3 1.0245 -1.6717
HALMD [3] 4 0.3 3 1.0234(3) -1.6731(4)
ReaDDy 4 0.3 3 1.035(7) -1.679(3)
[1] 4 0.6 3 3.69(1) -3.212(3)
[2] 4 0.6 3 3.7165 -3.2065
HALMD [3] 4 0.6 3 3.6976(8) -3.2121(2)
ReaDDy 4 0.6 3 3.70(2) -3.208(7)

[1] Molecular dynamics simulations, J. K. Johnson, J. A. Zollweg, and K. E. Gubbins, The Lennard-Jones equation of state revisited, Mol. Phys. 78, 591 (1993)

[2] Integral equations theory, A. Ayadim, M. Oettel, and S Amokrane, Optimum free energy in the reference functional approach for the integral equations theory, J. Phys.: Condens. Matter 21, 115103 (2009).

[3] HAL's MD package, http://halmd.org/validation.html

[4] Allen, M. P., & Tildesley, D. J. (1987). Computer Simulation of Liquids. New York: Oxford University Press.

In [2]:
import os
import numpy as np


87-1234567


Utility methods

In [3]:
def average_across_first_axis(values):
n_values = len(values)
mean = np.sum(values, axis=0) / n_values  # shape = n_bins
difference = values - mean  # broadcasting starts with last axis
std_dev = np.sqrt(np.sum(difference * difference, axis=0) / n_values)
std_err = np.sqrt(np.sum(difference * difference, axis=0) / n_values ** 2)
return mean, std_dev, std_err

def lj_system(edge_length, temperature=1.):
box_size=[edge_length, edge_length, edge_length],
unit_system=None
)
system.kbt = temperature
system.potentials.add_lennard_jones("A", "A", m=12, n=6, epsilon=1., sigma=1., cutoff=4., shift=True)
return system


Wrap the whole simulation and analysis in a function and perform it for the two densities 0.3 and 0.6.

In [4]:
def equilibrate_and_measure(density=0.3):
n_per_axis = 12
n_particles = n_per_axis ** 3
edge_length = (n_particles / density) ** (1. / 3.)
pos_x = np.linspace(-edge_length/2., edge_length/2.-1., n_per_axis)
pos = []
for x in pos_x:
for y in pos_x:
for z in pos_x:
pos.append([x,y,z])
pos = np.array(pos)
print("n_particles", len(pos), "box edge_length", edge_length)
assert len(pos)==n_particles

def pos_callback(x):
nonlocal pos
n = len(x)
pos = np.zeros((n,3))
for i in range(n):
pos[i][0] = x[i][0]
pos[i][1] = x[i][1]
pos[i][2] = x[i][2]
print("saved positions")

# create system
system = lj_system(edge_length, temperature=3.)

# equilibrate
sim = system.simulation(kernel="CPU")

sim.observe.particle_positions(2000, callback=pos_callback, save=None)
sim.observe.energy(500, callback=lambda x: print(x), save=None)

sim.record_trajectory(stride=1)
sim.output_file = "lj_eq.h5"
if os.path.exists(sim.output_file):
os.remove(sim.output_file)

sim.run(n_steps=10000, timestep=1e-4)

# measure
sim = system.simulation(kernel="CPU")
sim.observe.energy(200)
sim.observe.pressure(200)
sim.observe.rdf(
200, bin_borders=np.linspace(0.5, 4., 50),
types_count_from="A", types_count_to="A", particle_to_density=density)

sim.output_file = "lj_measure.h5"
if os.path.exists(sim.output_file):
os.remove(sim.output_file)

sim.run(n_steps=10000, timestep=1e-4)

# obtain results

energy_mean, _, energy_err = average_across_first_axis(energy) # time average
energy_mean /= n_particles
energy_err /= n_particles

pressure_mean, _, pressure_err = average_across_first_axis(pressure) # time average

rdf_mean, _, rdf_err = average_across_first_axis(rdf) # time average

return {
"energy_mean": energy_mean, "energy_err": energy_err,
"pressure_mean": pressure_mean, "pressure_err": pressure_err,
"rdf_mean": rdf_mean, "rdf_err": rdf_err, "rdf_bin_centers": bin_centers
}

In [5]:
result_low_dens = equilibrate_and_measure(density=0.3)

n_particles 1728 box edge_length 17.925618986228656
Configured kernel context with:
--------------------------------
- kBT = 3
- periodic b.c. = (true, true, true)
- box size = (17.9256, 17.9256, 17.9256)
- particle types:
*  particle type "A" with D=1
- potentials of order 2:
* for types "A" and "A"
* 12-6-Lennard-Jones potential with cutoff=4, epsilon=1, k=4, and with energy shift


saved positions
-2006.0926382760335
-2558.484698188715
-2683.608719842995
-2789.801886566924
saved positions
-2847.0129868701574
-2827.8083868918084
-2838.9676917594834
-2766.19525249713
saved positions
-2844.0582093567423
-2923.5162765016416
-2901.8219351147673
-2798.1972554259396
saved positions
-2882.6260913979854
-2894.1440317159036
-2810.6748812845053
-2862.09364697972
saved positions
-2946.395565470109
-2914.778299999032
-2891.783581769918
-2981.4438207366065
saved positions
-2887.27813830342
Configured kernel context with:
--------------------------------
- kBT = 3
- periodic b.c. = (true, true, true)
- box size = (17.9256, 17.9256, 17.9256)
- particle types:
*  particle type "A" with D=1
- potentials of order 2:
* for types "A" and "A"
* 12-6-Lennard-Jones potential with cutoff=4, epsilon=1, k=4, and with energy shift



In [ ]:
result_low_dens

Out[ ]:
{'energy_err': 0.003865080409371339,
'energy_mean': -1.6789505081645193,
'pressure_err': 0.007208254634738095,
'pressure_mean': 1.0356710396566735,
'rdf_bin_centers': array([0.53571429, 0.60714286, 0.67857143, 0.75      , 0.82142857,
0.89285714, 0.96428571, 1.03571429, 1.10714286, 1.17857143,
1.25      , 1.32142857, 1.39285714, 1.46428571, 1.53571429,
1.60714286, 1.67857143, 1.75      , 1.82142857, 1.89285714,
1.96428571, 2.03571429, 2.10714286, 2.17857143, 2.25      ,
2.32142857, 2.39285714, 2.46428571, 2.53571429, 2.60714286,
2.67857143, 2.75      , 2.82142857, 2.89285714, 2.96428571,
3.03571429, 3.10714286, 3.17857143, 3.25      , 3.32142857,
3.39285714, 3.46428571, 3.53571429, 3.60714286, 3.67857143,
3.75      , 3.82142857, 3.89285714, 3.96428571]),
'rdf_err': array([0.        , 0.        , 0.        , 0.        , 0.00051733,
0.00389936, 0.00758116, 0.00938435, 0.01152146, 0.00730478,
0.00639527, 0.00782449, 0.00638545, 0.00624519, 0.00497622,
0.00506881, 0.00451621, 0.00434441, 0.00469345, 0.00464019,
0.00408783, 0.00449142, 0.0043387 , 0.00377135, 0.00411156,
0.00372055, 0.0034688 , 0.00382219, 0.00327985, 0.00303205,
0.00316754, 0.00312204, 0.00272095, 0.00319668, 0.00277032,
0.00340234, 0.0030637 , 0.0029004 , 0.00296017, 0.00246223,
0.00234491, 0.0023808 , 0.0026758 , 0.00238175, 0.00250338,
0.00213579, 0.00181433, 0.00226536, 0.00239341]),
'rdf_mean': array([0.        , 0.        , 0.        , 0.        , 0.00149789,
0.14211495, 0.80149169, 1.39163288, 1.48260458, 1.42187995,
1.3107626 , 1.1942978 , 1.12439772, 1.04761764, 1.00864915,
0.98014494, 0.96831762, 0.96513191, 0.96171873, 0.97111502,
0.98073241, 1.00448663, 1.0144512 , 1.01665651, 1.01685795,
1.01207641, 1.00869444, 1.00878599, 1.01079231, 1.00434769,
1.00111387, 1.00042249, 1.00358795, 1.00785319, 1.00263455,
0.99972083, 1.00176477, 1.00691315, 1.0043989 , 1.0010145 ,
1.00141435, 1.00459604, 0.99894663, 1.00164453, 1.00343149,
1.00086893, 1.00622066, 1.00326314, 1.00458794])}
In [ ]:
result_hi_dens = equilibrate_and_measure(density=0.6)

n_particles 1728 box edge_length 14.227573217960249
Configured kernel context with:
--------------------------------
- kBT = 3
- periodic b.c. = (true, true, true)
- box size = (14.2276, 14.2276, 14.2276)
- particle types:
*  particle type "A" with D=1
- potentials of order 2:
* for types "A" and "A"
* 12-6-Lennard-Jones potential with cutoff=4, epsilon=1, k=4, and with energy shift


saved positions
-6881.589189751805
-5578.631408286779
-5553.045271029449
-5598.347635152243
saved positions
-5510.905288239784
-5494.26918009796
-5497.542802090885
-5650.728669069902
saved positions
-5547.803465031961
-5664.062193482018
-5470.78844587492
-5412.782056757059
saved positions
-5633.598773743057
-5585.178042184006
-5398.992296398446
-5651.0359006886965
saved positions
-5503.306269838077
-5527.3704879582365
-5522.8550998741
-5447.638729437222
saved positions
-5638.783772161413
Configured kernel context with:
--------------------------------
- kBT = 3
- periodic b.c. = (true, true, true)
- box size = (14.2276, 14.2276, 14.2276)
- particle types:
*  particle type "A" with D=1
- potentials of order 2:
* for types "A" and "A"
* 12-6-Lennard-Jones potential with cutoff=4, epsilon=1, k=4, and with energy shift


In [10]:
result_hi_dens

Out[10]:
{'energy_err': 0.007639845277153105,
'energy_mean': -3.2083370227235433,
'pressure_err': 0.025471798740872232,
'pressure_mean': 3.70051682353479,
'rdf_bin_centers': array([0.53571429, 0.60714286, 0.67857143, 0.75      , 0.82142857,
0.89285714, 0.96428571, 1.03571429, 1.10714286, 1.17857143,
1.25      , 1.32142857, 1.39285714, 1.46428571, 1.53571429,
1.60714286, 1.67857143, 1.75      , 1.82142857, 1.89285714,
1.96428571, 2.03571429, 2.10714286, 2.17857143, 2.25      ,
2.32142857, 2.39285714, 2.46428571, 2.53571429, 2.60714286,
2.67857143, 2.75      , 2.82142857, 2.89285714, 2.96428571,
3.03571429, 3.10714286, 3.17857143, 3.25      , 3.32142857,
3.39285714, 3.46428571, 3.53571429, 3.60714286, 3.67857143,
3.75      , 3.82142857, 3.89285714, 3.96428571]),
'rdf_err': array([0.        , 0.        , 0.        , 0.        , 0.00036706,
0.00360158, 0.00659086, 0.00752906, 0.00701379, 0.00766467,
0.00537993, 0.0052711 , 0.00339262, 0.00353886, 0.00351387,
0.00389974, 0.00361009, 0.00340888, 0.00325845, 0.00329971,
0.00322622, 0.00288786, 0.0029233 , 0.00311885, 0.002958  ,
0.00248027, 0.00225245, 0.00228648, 0.00220771, 0.00237958,
0.00256574, 0.00215778, 0.00196384, 0.00196462, 0.00209791,
0.00188725, 0.00200307, 0.00180375, 0.00230419, 0.00186438,
0.00190152, 0.00159949, 0.00191004, 0.00187251, 0.00183949,
0.00166731, 0.00167584, 0.00165395, 0.0014815 ]),
'rdf_mean': array([0.        , 0.        , 0.        , 0.        , 0.00224684,
0.21333092, 1.08795227, 1.71876871, 1.71941855, 1.49818423,
1.29725486, 1.11868612, 1.00248554, 0.93143182, 0.88929835,
0.86376007, 0.86582752, 0.88624511, 0.93118798, 0.97239683,
1.02267644, 1.04902949, 1.05419416, 1.05665055, 1.03252189,
1.02100544, 1.00978356, 0.9927093 , 0.98594912, 0.98351264,
0.98009485, 0.98725083, 0.99319732, 0.99637319, 1.00368474,
1.00537685, 1.00781406, 1.00839372, 1.00684037, 0.99703831,
0.99873489, 0.9962783 , 0.9965838 , 0.99828945, 0.99931174,
0.9978545 , 1.00089406, 1.00142241, 0.99987158])}
In [11]:
%matplotlib inline
import matplotlib.pyplot as plt

print("density 0.3:")
print("mean energy per particle {}\nerr energy per particle {}".format(
result_low_dens["energy_mean"], result_low_dens["energy_err"]))
print("pressure {}\nerr pressure {}".format(
result_low_dens["pressure_mean"], result_low_dens["pressure_err"]))
print("density 0.6:")
print("mean energy per particle {}\nerr energy per particle {}".format(
result_hi_dens["energy_mean"], result_hi_dens["energy_err"]))
print("pressure {}\nerr pressure {}".format(
result_hi_dens["pressure_mean"], result_hi_dens["pressure_err"]))

plt.plot(result_low_dens["rdf_bin_centers"], result_low_dens["rdf_mean"], label=r"density $\rho=0.3$")
plt.plot(result_hi_dens["rdf_bin_centers"], result_hi_dens["rdf_mean"], label=r"density $\rho=0.6$")
plt.xlabel(r"distance $r/\sigma$")
plt.ylabel(r"radial distribution $g(r)$")
plt.legend(loc="best")
plt.show()

density 0.3:
mean energy per particle -1.6789505081645193
err energy per particle 0.003865080409371339
pressure 1.0356710396566735
err pressure 0.007208254634738095
density 0.6:
mean energy per particle -3.2083370227235433
err energy per particle 0.007639845277153105
pressure 3.70051682353479
err pressure 0.025471798740872232