Add smooth cutoffs to distance tolerances

.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).

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.
 

docs/source/tutorials/visualisation.rst

@@ -0,0 +1,181 @@
+.. 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>`.

example/python_pkg/Si-Ge_learn/DRAFFLE/learn.py

@@ -7,6 +7,9 @@
 import numpy as np
 import os
 from joblib import Parallel, delayed
+from pathlib import Path
+
+script_dir = Path(__file__).resolve().parent
 
 import logging
 logging.basicConfig(level=logging.DEBUG)
@@ -16,7 +19,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
     # Check if the structure has already been processed
     if i < num_structures_old:
         return None, None, None
-    
+
     # calc = Vasp(**calc_params, label=f"struct{i}", directory=f"iteration{iteration}/struct{i}/", txt=f"stdout{i}.o")
     inew = i - num_structures_old
     atoms.calc = calc
@@ -33,7 +36,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
         except Exception as e:
             print(f"Optimisation failed: {e}")
             return None, None, None
-    
+
     # Save the optimised structure and its energy per atom
     energy_rlxd = atoms.get_potential_energy() / len(atoms)
 
@@ -47,24 +50,24 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
     if distances.min() < 1.0:
         print(f"Distance too small: {atoms.get_all_distances(mic=True).min()}")
         return None, None, None
-    
+
     if abs(energy_rlxd - energy_unrlxd) > 10.0:
         print(f"Energy difference too large: {energy_rlxd} vs {energy_unrlxd}")
         return None, None, None
-    
+
     return atoms, energy_unrlxd, energy_rlxd
 
 
 if __name__ == "__main__":
 
     # check if mace file exists
-    if not os.path.exists("../mace-mpa-0-medium.model"):
+    if not os.path.exists(script_dir  / ".." / ".." / "mace-mpa-0-medium.model"):
         print("MACE-MPA-0 model file not found. Please download the model from the MACE website.")
         print("https://github.com/ACEsuit/mace-foundations/releases/tag/mace_mpa_0")
         exit(1)
 
     # set up the calculator
-    calc_params = { 'model': "../mace-mpa-0-medium.model" }
+    calc_params = { 'model':  script_dir/ ".." / ".." / "mace-mpa-0-medium.model" }
     calc = mace_mp(**calc_params)
 
     # set up the hosts
@@ -123,7 +126,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
     # check if the energies file exists, if not create it
     energies_rlxd_filename = f"energies_rlxd_seed{seed}.txt"
     energies_unrlxd_filename = f"energies_unrlxd_seed{seed}.txt"
-    
+
     if os.path.exists(energies_rlxd_filename):
         with open(energies_rlxd_filename, "w") as energy_file:
             pass
@@ -135,7 +138,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
             pass
     else:
         open(energies_unrlxd_filename, "w").close()
-        
+
     # initialise the number of structures generated
     iter = -1
     unrlxd_structures = []
@@ -185,7 +188,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
                     energy=generated_structures[num_structures_old + i].get_potential_energy(),
                     forces=generated_structures[num_structures_old + i].get_forces()
                 )
-            
+
             # Start parallel execution
             print("Starting parallel execution")
             results = Parallel(n_jobs=5)(
@@ -218,7 +221,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
                     del unrlxd_structures[j]
                     del rlxd_structures[j]
                     generator.remove_structure(j)
-            num_structures_new = len(generated_structures) 
+            num_structures_new = len(generated_structures)
 
             # write the structures to files
             for i in range(num_structures_new - num_structures_old):
@@ -268,4 +271,4 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
     write(f"rlxd_structures_seed{seed}.traj", rlxd_structures)
     print("All generated and relaxed structures written")
 
-    print("Learning complete")
+    print("Learning complete")

example/python_pkg/graphene_grain_boundary_learn/DRAFFLE/learn.py

@@ -8,6 +8,9 @@
 import os
 from joblib import Parallel, delayed
 from ase.constraints import FixAtoms
+from pathlib import Path
+
+script_dir = Path(__file__).resolve().parent
 
 import logging
 logging.basicConfig(level=logging.DEBUG)
@@ -59,24 +62,24 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
 if __name__ == "__main__":
 
     # check if mace file exists
-    if not os.path.exists("../mace-mpa-0-medium.model"):
+    if not os.path.exists(script_dir  / ".." / ".." / "mace-mpa-0-medium.model"):
         print("MACE-MPA-0 model file not found. Please download the model from the MACE website.")
         print("https://github.com/ACEsuit/mace-foundations/releases/tag/mace_mpa_0")
         exit(1)
 
     # set up the calculator
-    calc_params = { 'model': "../mace-mpa-0-medium.model" }
+    calc_params = { 'model':  script_dir/ ".." / ".." / "mace-mpa-0-medium.model" }
     calc = mace_mp(**calc_params)
 
     # set up the hosts
     hosts = []
-    host = read("../POSCAR_host_gb")
+    host = read(script_dir  / ".." / "POSCAR_host_gb")
     hosts.append(host)
 
     # set the parameters for the generator
     generator = raffle_generator()
     generator.distributions.set_history_len(10)
-    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)