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

.. only:: html

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

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

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

.. _sphx_glr_gallery_monte_carlo_vs_analytical.py:


Monte Carlo RSA vs Percus Yevick
================================

This example compares the partial pair correlation functions :math:`g_{ij}(r)` obtained from:

1. A Monte Carlo Random Sequential Addition (RSA) packing in 3D (PackLab.monte_carlo)
2. The analytical Percus Yevick solution for a polydisperse hard sphere mixture (PackLab.analytical)

We use a two class discrete radius distribution with radii :math:`r = 1` and :math:`r = 2`,
sampled with equal number fractions. The Monte Carlo simulation is configured to
enforce the target radii distribution during packing, which is useful when larger
particles are rejected more often.

The workflow is:

1. Configure a periodic cubic domain, a discrete radius sampler, and RSA options
2. Run RSA and compute partial :math:`g_{ij}(r)` from the resulting configuration
3. Build the matching analytical polydisperse domain and solve Percus Yevick
4. Plot Monte Carlo curves and overlay the analytical solution

.. GENERATED FROM PYTHON SOURCE LINES 22-30

.. code-block:: Python


    import numpy as np
    import matplotlib.pyplot as plt

    from PackLab import monte_carlo, analytical, samplers
    from PackLab import ureg









.. GENERATED FROM PYTHON SOURCE LINES 31-34

Monte Carlo RSA setup
---------------------
We use a periodic cubic domain and a two radius discrete sampler.

.. GENERATED FROM PYTHON SOURCE LINES 34-69

.. code-block:: Python


    domain = monte_carlo.MCDomain(
        length_x=40.0 * ureg.micrometer,
        length_y=40.0 * ureg.micrometer,
        length_z=40.0 * ureg.micrometer,
        use_periodic_boundaries=True,
    )

    sampler = samplers.Discrete(
        radii=[1.0, 2.0] * ureg.micrometer,
        weights=[0.5, 0.5],
    )

    options = monte_carlo.Options()
    options.maximum_attempts = 2_500_000
    options.maximum_consecutive_rejections = 500_000
    options.target_packing_fraction = 0.24
    options.minimum_center_separation_addition = 0.0
    options.enforce_radii_distribution = True

    rsa_simulator = monte_carlo.Simulator(
        domain=domain,
        radius_sampler=sampler,
        options=options,
    )

    result = rsa_simulator.run()
    result.statistics.print()


    mc_centers, mc_g_ij = result.compute_partial_pair_correlation_function(
        n_bins=1_000,
        maximum_pairs=0,
    )








.. GENERATED FROM PYTHON SOURCE LINES 70-74

Analytical Percus Yevick setup
------------------------------
We construct an analytical polydisperse domain matching the Monte Carlo mixture.
The analytical domain uses Pint quantities.

.. GENERATED FROM PYTHON SOURCE LINES 74-103

.. code-block:: Python

    particle_radii, number_fractions = sampler.to_bins()

    py_domain = analytical.PYDomain(
        size=100_000 * ureg.micrometer,
        radii=particle_radii,
        volume_fraction=0.24,
        number_fractions=number_fractions,
    )

    # Percus Yevick solver radial frequency grid
    # Because we want to plot g we need to have a large p_max to capture the oscillations at small r.
    # The p_max should be at least 10 times 2 * pi/r_min, where r_min is the smallest particle radius.
    p_max = 1e3 / py_domain.radii.min()
    p = np.linspace(0, p_max * 5, 30_000)

    solver = analytical.Solver(
        densities=py_domain.particle_densities_per_radius,
        radii=py_domain.radii,
        p=p,
    )

    distances = np.linspace(
        py_domain.radii.min() * 2,
        py_domain.radii.max() * 10,
        400,
    )

    py_result = solver.compute(distances=distances)








.. GENERATED FROM PYTHON SOURCE LINES 104-107

Compare Monte Carlo and analytical g_ij(r)
------------------------------------------
We plot all partial curves on the same axes and overlay the Percus Yevick result in black.

.. GENERATED FROM PYTHON SOURCE LINES 107-120

.. code-block:: Python

    fig, ax = plt.subplots(1, 1)

    K = 2
    for i in range(K):
        for j in range(K):
            ax.plot(mc_centers, mc_g_ij[i, j], label=rf"RSA $g_{{{i}{j}}}(r)$")
            ax.plot(py_result.distances, py_result.g[i, j], color="black", linewidth=1.5)

    ax.set_xlabel("r")
    ax.set_ylabel(r"$g_{ij}(r)$")
    ax.set_title("Partial pair correlation: RSA vs Percus Yevick")
    ax.legend()
    plt.show()



.. image-sg:: /gallery/images/sphx_glr_monte_carlo_vs_analytical_001.png
   :alt: Partial pair correlation: RSA vs Percus Yevick
   :srcset: /gallery/images/sphx_glr_monte_carlo_vs_analytical_001.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_gallery_monte_carlo_vs_analytical.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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