.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/md/lj_fluid/plot_01_rahman_liquid_argon.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_md_lj_fluid_plot_01_rahman_liquid_argon.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_md_lj_fluid_plot_01_rahman_liquid_argon.py:


Rahman's liquid argon: g(r) of a Lennard-Jones liquid in real units
======================================================================

Rahman's 1964 simulation followed 864 argon atoms interacting through a
Lennard-Jones potential with :math:`\epsilon/k_B = 120\ \text{K}` and
:math:`\sigma = 3.4\ \text{\AA}`, at :math:`T = 94.4\ \text{K}` and a
mass density of :math:`1.374\ \text{g/cm}^3` -- liquid argon just above
its triple point. Here the same state point is simulated with
:class:`~chemistrykit.md.systems.lj_fluid.LJFluid` (512 atoms, reduced
units), and the radial distribution function
(:meth:`~chemistrykit.md.systems.lj_fluid.LJFluid.radial_distribution_function`)
is converted back to ångströms, as Rahman compared it with neutron- and
X-ray-scattering data on real liquid argon.

.. GENERATED FROM PYTHON SOURCE LINES 18-50

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from chemistrykit.constants import NA
    from chemistrykit.md.systems.lj_fluid import LJFluid
    from chemistrykit.md.systems.thermostats import VelocityRescalingThermostat

    epsilon_over_k, sigma_A, molar_mass = 120.0, 3.4, 39.948  # K, angstrom, g/mol (Rahman's argon)
    T_kelvin, rho_mass = 94.4, 1.374  # K, g/cm^3

    rho_reduced = rho_mass / molar_mass * NA * (sigma_A * 1e-8) ** 3
    T_reduced = T_kelvin / epsilon_over_k
    print(f"reduced density rho* = {rho_reduced:.3f}, reduced temperature T* = {T_reduced:.3f}")

    fluid = LJFluid.from_lattice(n_per_side=8, density=rho_reduced, temperature=T_reduced, cutoff=2.5, rng=0)
    # melt the starting lattice and equilibrate at Rahman's temperature, then run NVE
    fluid.run(dt=0.005, n_steps=3000, thermostat=VelocityRescalingThermostat(T_reduced, interval=10), sample_every=3000)
    result = fluid.run(dt=0.005, n_steps=2000, sample_every=100)

    g_frames = [fluid.radial_distribution_function(n_bins=150, positions=p)[1] for p in result.positions]
    r, _ = fluid.radial_distribution_function(n_bins=150)
    g = np.mean(g_frames, axis=0)
    r_A = r * sigma_A

    fig, ax = plt.subplots(figsize=(7, 4.5))
    ax.plot(r_A, g, color="steelblue")
    ax.axhline(1.0, color="gray", linestyle=":", linewidth=0.8)
    ax.set_xlabel(r"r ($\mathrm{\AA}$)")
    ax.set_ylabel("g(r)")
    ax.set_title(f"Liquid argon at {T_kelvin} K, {rho_mass} g/cm$^3$ (Rahman's state point)")
    fig.tight_layout()




.. image-sg:: /api/gallery/md/lj_fluid/images/sphx_glr_plot_01_rahman_liquid_argon_001.png
   :alt: Liquid argon at 94.4 K, 1.374 g/cm$^3$ (Rahman's state point)
   :srcset: /api/gallery/md/lj_fluid/images/sphx_glr_plot_01_rahman_liquid_argon_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    reduced density rho* = 0.814, reduced temperature T* = 0.787




.. GENERATED FROM PYTHON SOURCE LINES 51-57

The first peak sits near 3.7 Å and the structure dies away after two or
three coordination shells -- the liquid signature Rahman compared with
experiment. Integrating :math:`4\pi\rho r^2 g(r)` out to the first
minimum counts the nearest neighbors of an atom: a little more than a
dozen, a nearly close-packed first shell like that of the crystal the
liquid freezes into:

.. GENERATED FROM PYTHON SOURCE LINES 57-66

.. code-block:: Python


    i_peak = int(np.argmax(g))
    i_min = i_peak + int(np.argmin(g[i_peak : i_peak + 40]))
    dr = r[1] - r[0]
    coordination = float(np.sum(4.0 * np.pi * rho_reduced * r[:i_min] ** 2 * g[:i_min]) * dr)
    print(f"first peak at r = {r_A[i_peak]:.2f} angstrom, first minimum at {r_A[i_min]:.2f} angstrom")
    print(f"coordination number (first shell): {coordination:.1f}")

    plt.show()




.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    first peak at r = 3.74 angstrom, first minimum at 5.58 angstrom
    coordination number (first shell): 13.6





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 2.789 seconds)


.. _sphx_glr_download_api_gallery_md_lj_fluid_plot_01_rahman_liquid_argon.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_01_rahman_liquid_argon.ipynb <plot_01_rahman_liquid_argon.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_01_rahman_liquid_argon.py <plot_01_rahman_liquid_argon.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_01_rahman_liquid_argon.zip <plot_01_rahman_liquid_argon.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_