Add smooth cutoffs to distance tolerances (#44) and improve documentation

## Added
- Added smooth cutoffs to distance tolerances
- Added smooth cutoffs to 3- and 4-body contributions to `evaluate_point()`
  - Viability function is now a more continuous surface (abrupt when transitioning below cutoff_min)
- Added descriptor generator procedure
- Added unit tests for distribution convergence
- Added probability density generation procedure
- Added `set_method_ratio_default` procedure

## Changed
- Changed Python docstrings style to match [NumPy Style Guide](https://numpydoc.readthedocs.io/en/latest/format.html)
- Specified expected types for Python arguments
- Converted all custom `modu` calls to `norm2`
- Optimised `evaluate` and `viability` modules
- Optimised array indexing
- Moved `atom_ignore_index` array to boolean inside `basis_type`
  - Acts as a mask for atoms

## Fixed
- Fixed handling of optional method ratio

Author: nedtaylor

.github/PULL_REQUEST_TEMPLATE/pull_request_template.md

@@ -0,0 +1,14 @@
+# Description
+
+> _Briefly describe the purpose and content of this merge request._
+> _E.g. "Add structure generation constraints", or "Fix bug in energy parsing"._
+
+# Checklist
+
+Mark with `x` when complete, or `~` if not applicable.
+
+- [ ] [I have read and followed the **RAFFLE's contribution guidelines.**](https://github.com/ExeQuantCode/RAFFLE/blob/main/CONTRIBUTING.md)
+- [ ] **Code is commented** appropriately, and API docstrings follow NumPy or FORD style.
+- [ ] **Read*the*Docs** documentation is added/updated (if applicable).
+- [ ] **Unit tests** are added/updated (if applicable).
+- [ ] **Linked issue** is resolved with a `closes #XXXX` reference (if applicable).

.pre-commit-config.yaml

@@ -5,7 +5,7 @@ repos:
     -   id: trailing-whitespace
     -   id: end-of-file-fixer
 -   repo: https://github.com/nedtaylor/fortran-format-hooks
-    rev: 1aad19af1f0c87027829f0f872d4bdc6ed9251e2
+    rev: 749c61fce3a5f29f6d5d3b236ed9b693c6b2bf4e
     hooks:
     -   id: check-fortran-indentation
         args: [--line-length=80, --ignore-directories=src/wrapper]

CMakeLists.txt

@@ -74,6 +74,7 @@ set(LIB_DIR ${FORTRAN_SRC_DIR}/lib)
 set(LIB_FILES
 	mod_io_utils.F90
 	mod_constants.f90
+    mod_cache.f90
 	mod_misc.f90
     mod_tools_infile.f90
 	mod_misc_maths.f90

CONTRIBUTING.md

@@ -92,7 +92,7 @@ When submitting your contributions, please ensure the following:
 - Reference any related issues or pull requests, if applicable.
 - Write unit tests for your contributions
 - Ensure all existing tests pass before submitting your changes.
-- Update the documentation to reflect your changes, if necessary (i.e. through FORD style commenting).
+- Update the documentation to reflect your changes, if necessary (i.e. through FORD style commenting for Fortran and NumPy docstrings style for Python).
 - Provide examples and usage instructions, if applicable.
 
 Follow the [Code Style](#code-style) when contributing code to this project to ensure compatibility and a uniform format to the project.
@@ -101,7 +101,7 @@ Follow the [Code Style](#code-style) when contributing code to this project to e
 ### Code Style
 - Follow the existing code style and conventions.
 - Use meaningful variable and function names.
-- Write clear and concise comments. For the Fortran library, use comments compatible with the [FORD Fortran Documenter](https://forddocs.readthedocs.io/en/stable/). For the Python wrapper, use comments compatible with [pandoc](https://pandoc.org).
+- Write clear and concise comments. For the Fortran library, use comments compatible with the [FORD Fortran Documenter](https://forddocs.readthedocs.io/en/stable/). For the Python wrapper, use comments compatible with [pandoc](https://pandoc.org), following the [NumPy style guide](https://numpydoc.readthedocs.io/en/latest/format.html).
 
 
 

docs/source/tutorials/graphene_grain_boundary_tutorial.rst

@@ -30,6 +30,9 @@ First, we must import the required packages:
     from raffle.generator import raffle_generator
     from mace.calculators import mace_mp
     import numpy as np
+    from pathlib import Path
+
+    script_dir = Path(__file__).resolve().parent
 
 
 Next, we need to set up the RAFFLE generator and the calculator to calculate the energies of the structures.
@@ -55,14 +58,14 @@ The host is read in from a file, but it can also be generated using the ARTEMIS
 
   .. code-block:: python
 
-    host = read("../POSCAR_host_gb")
+    host = read(script_dir  / ".." / "POSCAR_host_gb")
     generator.set_host(host)
 
 We then need to set up the RAFFLE generator by creating the descriptor.
 
 .. code-block:: python
 
-    graphene = read("../POSCAR_graphene")
+    graphene = read(script_dir / ".." / "POSCAR_graphene")
     h2 = build.molecule("H2")
     graphene.calc = calc
     C_reference_energy = graphene.get_potential_energy() / len(graphene)

docs/source/tutorials/index.rst

@@ -35,3 +35,9 @@ Whilst RAFFLE is a random sturcture search package designed primarlily for inter
    aluminium_tutorial
    Si-Ge_tutorial
    graphene_grain_boundary_tutorial
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Visualisation:
+
+   visualisation

docs/source/tutorials/parameters_tutorial.rst

@@ -14,6 +14,7 @@ For a guide on how to build a database, see the :doc:`Databases tutorial </tutor
 
 Initialisation
 --------------
+
 RAFFLE is initialised by importing the generator object.
 This object is the main interface for the user to interact with the RAFFLE package.
 
@@ -61,6 +62,54 @@ As such, it is recommended that the user implement their own convergence check i
     The convergence criterion is also not enforced even when set, it is simply a check for the user to use.
 
 
+Placement methods
+-----------------
+
+RAFFLE uses a random sampling of five placement methods to place atoms in the host structure.
+These methods are:
+- ``rand``: Randomly place atoms in the host structure.
+- ``void``: Place atoms at the point furthest from the nearest atom in the host structure.
+- ``walk``: Start from a random position and walk to a local minimum in the RAFFLE-calculated probability.
+- ``grow``: Start from the previously placed atom and walk to a local minimum in the RAFFLE-calculated probability.
+- ``min``: Place atoms at the point of highest RAFFLE-calculated probability.
+
+The ratio of these methods can be set by the user.
+This can either be done on a per-generate method basis, or as a global setting.
+
+To set the default placement method ratios, the user can use the ``set_method_ratio_default`` method.
+
+.. code-block:: python
+
+    # Set default placement method ratios
+    generator.set_method_ratio_default(
+        method_ratio = {
+            'rand': 0.1,
+            'void': 0.1,
+            'walk': 0.25,
+            'grow': 0.25,
+            'min': 1.0
+        }
+    )
+
+These defaults will then be used for all generation methods where the user does not specify a method ratio.
+
+To set the placement method ratios for a specific generation method, the user can define the method ratio ``dict`` in the ``generate`` method.
+
+.. code-block:: python
+
+    # Set placement method ratios for a specific generation method
+    generator.generate(
+        method_ratio = {
+            'rand': 0.1,
+            'void': 0.1,
+            'walk': 0.25,
+            'grow': 0.25,
+            'min': 1.0
+        },
+        # Other generation parameters
+        ...
+    )
+
 
 Energy references
 -----------------

docs/source/tutorials/visualisation.rst

@@ -0,0 +1,229 @@
+.. visualisation:
+
+=============
+Visualisation
+=============
+
+Here we detail the methods available for outputting the learned RAFFLE descriptor and RAFFLE fingerprints (distribution functions) for individual structures.
+
+
+Visualising the learned RAFFLE descriptor
+-----------------------------------------
+
+RAFFLE generates 2-body, 3-body, and 4-body distribution functions for each atomic species in the system (element pairs for the 2-body function).
+These are generated by combining formation energy-weighted (or convex hull-weighted) `n`-body distribution functions for each structure provided in the learning database.
+
+The ``raffle_generator`` object has a method ``get_descriptor()`` that returns the 2-, 3-, and 4-body forms of the learned RAFFLE generalised descriptor.
+The output is a list of three numpy arrays, with the first array containing the 2-body descriptor, the second array containing the 3-body descriptor, and the third array containing the 4-body descriptor.
+Each `n`-body descriptor is a 2D array, with the first column containing the species index (or element pair index for the 2-body descriptor) and the second column containing the binned descriptor value.
+The bin lengths are (``nbins`` component) set either explicitly or determined by the ``cutoff_min``, ``cutoff_max``, and ``width`` components of the generator.
+
+Here is an example of how to use the ``get_descriptor()`` method:
+
+.. code-block:: python
+
+    # Initialise RAFFLE generator
+    from raffle.generator import raffle_generator
+
+    generator = raffle_generator()
+
+    # Set the host structure
+    host = Atoms(
+        # Host structure for the generator
+    )
+    generator.set_host(host)
+
+    # Set the reference energies (i.e. chemical potential references)
+    generator.distributions.set_element_energies(
+        {
+            # reference energies for all elements in the systems
+        }
+    )
+
+    # Optional parameters
+    generator.distributions.set_kBT(0.2)
+    generator.distributions.set_width([0.04, np.pi/160.0, np.pi/160.0])
+    generator.distributions.set_cutoff_min([0.5, 0.0, 0.0])
+    generator.distributions.set_cutoff_max([6.0, np.pi, np.pi])
+
+    # Set and learn from the initial database
+    database = [
+        # List of structures in the learning database
+    ]
+    generator.distributions.create(database)
+
+    # Retrieve the descriptor
+    descriptor_init = generator.get_descriptor()
+
+    # Print the 2-body descriptor
+    print("2-body descriptor:")
+    print(descriptor_init[0])
+
+    # Print the 3-body descriptor
+    print("3-body descriptor:")
+    print(descriptor_init[1])
+
+    # Print the 4-body descriptor
+    print("4-body descriptor:")
+    print(descriptor_init[2])
+
+With this, we can now plot the descriptor using any plotting library of your choice.
+
+.. code-block:: python
+
+    import matplotlib.pyplot as plt
+    import numpy as np
+
+    # Create a figure with 3 subplots side by side
+    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
+
+    # Plot for each n-body descriptor (2-body, 3-body, 4-body)
+    for j in range(3):
+        # Calculate x-axis values
+        x = np.arange(generator.distributions.cutoff_min[j],
+                    generator.distributions.cutoff_max[j] + generator.distributions.width[j],
+                    generator.distributions.width[j])
+
+        # Plot on the respective subplot
+        for idx in range(len(descriptor_init[j])):
+            axes[j].plot(x, descriptor_init[j][idx,:])
+
+        # Set labels and title for each subplot
+        axes[j].set_ylabel('Descriptor value')
+        axes[j].set_title(f'{j+2}-body descriptor')
+
+    axes[0].set_xlabel('Distance (Å)')
+    axes[1].set_xlabel('3-body angle (radians)')
+    axes[2].set_xlabel('Improper dihedral angle (radians)')
+    plt.tight_layout()
+    plt.show()
+
+An example python notebook is provided in :git:`examples/python_pkg/visualisation/descriptor.ipynb <examples/python_pkg/visualisation/descriptor.ipynb>`
+
+We can now use this to compare the initial descriptor with the updated descriptor after generating new structures.
+
+.. code-block:: python
+
+    # Generate new structures and update the descriptor
+    structures = [
+        # List of structures to be generated
+    ]
+    generator.distributions.update(structures)
+
+    # Retrieve the updated descriptor
+    descriptor_new = generator.get_descriptor()
+
+    # Print the updated descriptor on the plots and compare
+    ...
+
+
+Visualising a RAFFLE fingerprint
+--------------------------------
+
+RAFFLE fingerprints are the distribution functions for each structure in the learning database.
+These are then weighted by energy (formation or convex hull) to form the RAFFLE descriptor.
+
+However, the individual fingerprints can also be extracted and visualised.
+
+The `raffle_generator` object has a method `get_fingerprint()` that returns the distribution functions for a provided structure.
+The output is a list of three numpy arrays, with the first array containing the 2-body fingerprint, the second array containing the 3-body fingerprint, and the third array containing the 4-body fingerprint.
+Each `n`-body fingerprint is a 2D array, with the first column containing the species index (or element pair index for the 2-body fingerprint) and the second column containing the binned fingerprint value.
+Like above, the bin lengths are set either explicitly or determined by the `cutoff_min`, `cutoff_max`, and `width` components of the generator.
+Here is an example of how to use the `get_fingerprint()` method:
+
+.. code-block:: python
+
+    # Initialise RAFFLE generator
+    from raffle.generator import raffle_generator
+
+    generator = raffle_generator()
+
+    # Optional parameters
+    generator.distributions.set_width([0.04, np.pi/160.0, np.pi/160.0])
+    generator.distributions.set_cutoff_min([0.5, 0.0, 0.0])
+    generator.distributions.set_cutoff_max([6.0, np.pi, np.pi])
+
+    # Structure to obtain the fingerprint for
+    structure = Atoms(
+        # Structure to be used for the fingerprint
+    )
+
+    fingerprint = generator.distributions.generate_fingerprint(structure)
+
+This can then be visualised in a similar way to the descriptor.
+
+.. code-block:: python
+
+    # Create a figure with 3 subplots side by side
+    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
+
+    # Plot for each n-body function (2-body, 3-body, 4-body)
+    for j in range(3):
+        # Calculate x-axis values
+        x = np.arange(generator.distributions.cutoff_min[j],
+                    generator.distributions.cutoff_max[j] + generator.distributions.width[j],
+                    generator.distributions.width[j])
+
+        # Plot on the respective subplot
+        for idx in range(len(fingerprint[j])):
+            axes[j].plot(x, fingerprint[j][idx,:])
+
+        # Set labels and title for each subplot
+        axes[j].set_ylabel('Fingerprint value')
+        axes[j].set_title(f'{j+2}-body fingerprint')
+
+    axes[0].set_xlabel('Distance (Å)')
+    axes[1].set_xlabel('3-body angle (radians)')
+    axes[2].set_xlabel('Improper dihedral angle (radians)')
+    plt.tight_layout()
+    plt.show()
+
+An example python notebook is provided in  :git:`examples/python_pkg/visualisation/fingerprint.ipynb <examples/python_pkg/visualisation/fingerprint.ipynb>`.
+
+
+Visualising RAFFLE probability density
+--------------------------------------
+
+RAFFLE probability density is the probability of finding a given element in a given position in the system.
+This is calculated by the RAFFLE generator and can be visualised using the `get_probability_density()` method.
+The output is a 2D array, with the first column containing the coordinates (spatial and species) the second column containing the binned probability density value.
+
+A structure is provided to the `get_probability_density()`, along with a list of elements to calculate the probability density for.
+These elements must be present in the RAFFLE descriptor.
+
+An example of how to use the `get_probability_density()` method is shown below:
+
+.. code-block:: python
+
+    # Initialise RAFFLE generator
+    from raffle.generator import raffle_generator
+
+    generator = raffle_generator()
+
+    generator.distributions.set_element_energies(
+        {
+            # reference energies for all elements in the systems
+        }
+    )
+
+    database = [
+        # List of structures in the learning database
+    ]
+    generator.distributions.create(database)
+
+    # Structure to obtain the probability density for
+    structure = Atoms(
+        # Structure to be used for the probability density
+    )
+
+    species = 'SiGe'
+
+    probability_density, grid = generator.get_probability_density(structure, species, return_grid=True)
+
+The first index of the first column of `probability_density` is the x-coordinate, the second index is the y-coordinate, and the third index is the z-coordinate.
+The fourth index is the distance between the position and the nearest atom (i.e. the void value).
+The fifth index onwards is the species index (in order of the species list provided).
+The second column is the site index.
+
+For a more extensive example, see the `examples/python_pkg/visualisation/probability_density.ipynb` notebook.
+This also provides a visualisation of the probability density using the `matplotlib` library.