ReaDDy - A particle-based
reaction-diffusion simulator

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

import readdy
print(readdy.__version__)
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.):
    system = readdy.ReactionDiffusionSystem(
        box_size=[edge_length, edge_length, edge_length],
        unit_system=None
    )
    system.kbt = temperature
    system.add_species("A", diffusion_constant=1.)
    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.add_particles("A", pos)

    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)

    traj = readdy.Trajectory(sim.output_file)
    traj.convert_to_xyz(particle_radii={"A": 0.5})
    
    # measure
    sim = system.simulation(kernel="CPU")
    sim.add_particles("A", pos)
    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
    traj = readdy.Trajectory(sim.output_file)
    _, energy = traj.read_observable_energy()
    _, bin_centers, rdf = traj.read_observable_rdf()
    _, pressure = traj.read_observable_pressure()
    
    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
© Copyright 2018 CMB group