FRQ: lattice function

[kay-ro]

Proposed new feature or change:

The local density in particle simulations is calculated at the locations of the particles but in cases where there MIPS is collapse-like, this causes a sometimes unwanted underrepresentation of the gaseous phase. One option is to calculate the local density by first calculating a density field from the particle simulations but this is computationally expensive in large simulations. An alternative is to use the same method as with the local density but evaluate at evenly spaced grid points throughout the simulation box.

To do so, a lattice function should should be implemented that allows to calculate the coordinates of specific lattices in 2D (and 3D). A square grid should be the focus. The user should be able to give a number of points or a specific number of grid points in all directions.

ToDo:

• [ ] add function amep.utils.lattice (gets number of lattice points and box)
• [ ] add detailed documentation
• [ ] detailed testing of the module to verify the created lattices
• [ ] add example (also to calculation of local density distribution) where the lattice is used to calculate the local density distribution

[kay-ro]

A hexagonal lattice could be implemented as follows:

def hexagonal(a, N, f=1.0, center=np.array([0.0,0.0,0.0])):
    '''
    Returns the coordinates of the points on a 2D hexagonal
    lattice of sites with a width-to-height ratio f centered
    around the given center. The final particle number N can
    slightly differ from the given value of N.

    Notes
    -----
    See also notes from September 1, 2022.

    Parameters
    ----------
    a : float
        Length of the lattice site.
    N : float
        Number of particles.
    f : float, optional
        Width-to-height ratio. The default is 1.0.
    center : np.ndarray, optional
        Center of the lattice. The default is the origin.

    Returns
    -------
    coords : np.ndarray
        Coordinate frame of the lattice points.

    Examples
    --------
    >>> hex = hexagonal(2, 100, f=1)
    >>> fig, axs = amep.plot.new(figsize=(8,8))
    >>> axs.scatter(hex[:,0], hex[:,1])
    >>> 

    '''
    Nx = int(np.sqrt(np.sqrt(3)*N*f/2)+0.5)
    Ny = int(np.sqrt(2*N/f/np.sqrt(3))+0.5)
    coords = np.zeros((int(Nx*Ny),3))
    n = 0
    for k in range(Nx):
        for j in range(Ny):
            if j%2 == 0:
                coords[n] = np.array([(2*k)*a/2, j*np.sqrt(3)*a/2, 0.0])
            else:
                coords[n] = np.array([(1+2*k)*a/2, j*np.sqrt(3)*a/2, 0.0])
            n += 1

    # shift coordinates to the given center
    coords = coords - np.mean(coords, axis=0) + center
    return coords