.. 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_04_wca_repulsive_reference.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_04_wca_repulsive_reference.py>`
        to download the full example code.

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

.. _sphx_glr_api_gallery_md_lj_fluid_plot_04_wca_repulsive_reference.py:


Weeks-Chandler-Andersen: repulsive forces set the structure of a dense liquid
================================================================================

Weeks, Chandler, and Andersen split the Lennard-Jones potential at its
minimum :math:`r_m=2^{1/6}\sigma` into a purely repulsive reference part
:math:`u_0` and a smooth attractive perturbation :math:`w`:

.. math::

   u_0(r) = \begin{cases} u(r) + \epsilon, & r < r_m \\ 0, & r \ge r_m \end{cases}
   \qquad
   w(r) = \begin{cases} -\epsilon, & r < r_m \\ u(r), & r \ge r_m \end{cases}

and showed that at liquid densities the structure of the full fluid is
almost the same as that of the repulsive reference fluid alone. This
example draws the split, then simulates both fluids with
:class:`~chemistrykit.md.systems.lj_fluid.LJFluid` at a dense liquid
state point near the triple point (:math:`\rho^*=0.84`,
:math:`T^*=0.75`) -- the full potential with the usual cutoff
:math:`2.5\sigma`, the WCA reference by cutting at
:attr:`~chemistrykit.md.systems.lj_fluid.LennardJones.r_min` -- and
compares their radial distribution functions.

.. GENERATED FROM PYTHON SOURCE LINES 27-51

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from chemistrykit.md.systems.lj_fluid import LennardJones, LJFluid
    from chemistrykit.md.systems.thermostats import VelocityRescalingThermostat

    lj = LennardJones()
    r_m = lj.r_min
    r = np.linspace(0.9, 2.5, 400)
    u = lj.energy(r)
    u0 = np.where(r < r_m, u + 1.0, 0.0)
    w = np.where(r < r_m, -1.0, u)

    fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
    axes[0].plot(r, u, color="black", label="full LJ  u = u0 + w")
    axes[0].plot(r, u0, color="darkorange", linestyle="--", label="repulsive reference u0")
    axes[0].plot(r, w, color="steelblue", linestyle="--", label="attractive perturbation w")
    axes[0].axvline(r_m, color="gray", linestyle=":", linewidth=0.8)
    axes[0].set_ylim(-1.3, 2.0)
    axes[0].set_xlabel(r"r / $\sigma$")
    axes[0].set_ylabel(r"energy / $\epsilon$")
    axes[0].set_title("The WCA split at the potential minimum")
    axes[0].legend()




.. image-sg:: /api/gallery/md/lj_fluid/images/sphx_glr_plot_04_wca_repulsive_reference_001.png
   :alt: The WCA split at the potential minimum
   :srcset: /api/gallery/md/lj_fluid/images/sphx_glr_plot_04_wca_repulsive_reference_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x11ccfcd70>



.. GENERATED FROM PYTHON SOURCE LINES 52-55

Now the two fluids at the same dense state point. Each is equilibrated
with a velocity-rescaling thermostat and g(r) is averaged over an
NVE production run.

.. GENERATED FROM PYTHON SOURCE LINES 55-74

.. code-block:: Python


    density, T = 0.84, 0.75
    g_curves = {}
    for label, cutoff, color in [("full LJ (cutoff 2.5 sigma)", 2.5, "black"), ("WCA repulsive reference", r_m, "darkorange")]:
        fluid = LJFluid.from_lattice(n_per_side=7, density=density, temperature=T, cutoff=cutoff, rng=0)
        fluid.run(dt=0.004, n_steps=3000, thermostat=VelocityRescalingThermostat(T, interval=10), sample_every=3000)
        result = fluid.run(dt=0.004, n_steps=1500, thermostat=VelocityRescalingThermostat(T, interval=10), sample_every=100)
        rr, _ = fluid.radial_distribution_function(n_bins=120)
        g = np.mean([fluid.radial_distribution_function(n_bins=120, positions=p)[1] for p in result.positions], axis=0)
        g_curves[label] = g
        axes[1].plot(rr, g, color=color, label=label)

    axes[1].axhline(1.0, color="gray", linestyle=":", linewidth=0.8)
    axes[1].set_xlabel(r"r / $\sigma$")
    axes[1].set_ylabel("g(r)")
    axes[1].set_title(rf"Dense liquid, $\rho^*$={density}, $T^*$={T}")
    axes[1].legend()
    fig.tight_layout()








.. GENERATED FROM PYTHON SOURCE LINES 75-79

The two radial distribution functions nearly coincide: removing the
entire attractive tail barely changes where the neighbors sit, because
in a dense liquid the attractions from all sides roughly cancel and the
packing is dictated by the repulsive cores.

.. GENERATED FROM PYTHON SOURCE LINES 79-85

.. code-block:: Python


    g_full, g_wca = g_curves.values()
    print(f"first-peak height: full LJ {g_full.max():.2f}, WCA reference {g_wca.max():.2f}")
    print(f"largest pointwise difference in g(r): {np.max(np.abs(g_full - g_wca)):.2f}")

    plt.show()




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

 .. code-block:: none

    first-peak height: full LJ 2.94, WCA reference 2.97
    largest pointwise difference in g(r): 0.27





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

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


.. _sphx_glr_download_api_gallery_md_lj_fluid_plot_04_wca_repulsive_reference.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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