Pack a number of non-intersecting spheres into a periodic system.

    Can specify packing by number of spheres or packing fraction.

    This implementation uses the Metropolis-Hastings algorithm for an
    NPT system.

    Parameters
    ----------
    R: float
        Sphere radius.
    L: float array, shape (d,)
        System lengths.
    pf: float or None
        Packing fraction
    n: integer or None
        Number of spheres.
    rng: RandomState or None
        Random number generator. If None, use inbuilt numpy state.
    periodic: bool
        Whether or not the system is periodic.


    Metropolis-Hastings parameters
    ------------------------------
    Playing with these parameters may improve packing speed.

    beta_max: float, greater than zero.
        Inverse temperature which controls how little noiseis in the system.
    dL_max: float, 0 < dL_max < 1
        Maximum fraction by which to perturb the system size.
    dr_max: float, 0 < dr_max < 1
        Maximum system fraction by which to perturb sphere positions.

    Returns
    -------
    r: float array, shape (n, d)
        Coordinates of the centres of the spheres for a valid configuration.
    R_actual: float
        Actual sphere radius used in the packing.

def pack(R, L, pf=None, n=None, rng=None, periodic=False,
         beta_max=1e4, dL_max=0.02, dr_max=0.02):
    """Pack a number of non-intersecting spheres into a periodic system.

    Can specify packing by number of spheres or packing fraction.

    This implementation uses the Metropolis-Hastings algorithm for an
    NPT system.

    Parameters
    ----------
    R: float
        Sphere radius.
    L: float array, shape (d,)
        System lengths.
    pf: float or None
        Packing fraction
    n: integer or None
        Number of spheres.
    rng: RandomState or None
        Random number generator. If None, use inbuilt numpy state.
    periodic: bool
        Whether or not the system is periodic.


    Metropolis-Hastings parameters
    ------------------------------
    Playing with these parameters may improve packing speed.

    beta_max: float, greater than zero.
        Inverse temperature which controls how little noiseis in the system.
    dL_max: float, 0 < dL_max < 1
        Maximum fraction by which to perturb the system size.
    dr_max: float, 0 < dr_max < 1
        Maximum system fraction by which to perturb sphere positions.

    Returns
    -------
    r: float array, shape (n, d)
        Coordinates of the centres of the spheres for a valid configuration.
    R_actual: float
        Actual sphere radius used in the packing.
    """
    if pf is not None:
        if pf == 0.0:
            return np.array([]), R
        # If packing fraction is specified, find required number of spheres
        # and the actual packing fraction this will produce
        n, pf_actual = pf_to_n(L, pf, R)
    elif n is not None:
        if n == 0:
            return np.array([]), R
        # If n is specified, find packing fraction
        pf_actual = n_to_pf(L, n, R)

    # Calculate an initial packing fraction and system size
    # Start at at most 0.5%; lower if the desired packing fraction is very low
    pf_initial = min(0.005, pf_actual / 2.0)
    # Find system size that will create this packing fraction
    dim = L.shape[0]
    increase_initial_ratio = (pf_actual / pf_initial) ** (1.0 / dim)
    L_0 = L * increase_initial_ratio

    # Pack naively into this system
    r_0, R = pack_simple(R, L_0, n=n, rng=rng, periodic=periodic)

    mg = metro_rcp_factory(periodic, r_0, L_0, R, dr_max, dL_max, rng=rng)

    print('Initial packing done, Initial packing: {:g}'.format(mg.pf))

    t = 0
    while mg.pf < pf_actual:
        t += 1
        beta = beta_max * mg.pf
        mg.iterate(beta)

        if not t % every:
            print('Packing: {:g}%'.format(100.0 * mg.pf))

    print('Final packing: {:g}%'.format(100.0 * mg.pf))

    return mg.r, mg.R