Descriptor convergence

docs/source/tutorials/generating_tutorial.rst

@@ -31,6 +31,15 @@ The arguments are as follows:
 - ``return_exit_code``: Whether to return the exit code of the generator.
 
 
+The accepted keys for the ``method_ratio`` dictionary are:
+
+- ``min`` (or ``minimum`` or ``global``): The minimum method.
+- ``walk``: The walk method.
+- ``grow`` or (``growth``): The grow method.
+- ``void``: The void method.
+- ``rand`` or (``random``): The random method.
+
+
 Example
 -------
 

docs/source/tutorials/parameters_tutorial.rst

@@ -29,6 +29,39 @@ It is recommended to use the Atomic Simulation Environment (ASE)~\cite{ase-paper
 Whilst RAFFLE can handle its own atomic structure object, ASE is more widely used and has a more extensive feature set.
 
 
+Convergence criteria
+--------------------
+
+RAFFLE has a convergence criterion that compares the amount of change in the RAFFLE descriptor over the last ``history_len`` update iterations.
+When the total change in the descriptor for the last ``history_len`` iterations is less than ``threshold``, the generator can be considered converged.
+Without explicitly setting the ``history_len``, the convergence is not checked.
+
+To set the convergence criteria, the user can use the ``set_history_len`` method.
+
+.. code-block:: python
+
+    # Set convergence criteria
+    generator.set_history_len(
+        history_len=10
+    )
+
+To check the convergence, the user can use the ``is_converged`` method (default ``threshold = 1e-4``).
+
+.. code-block:: python
+
+    # Check convergence
+    converged = generator.is_converged(threshold=1e-5)
+    print(f"Converged: {converged}")
+
+Note, the convergence criterion is not enforced, i.e. the generator will not stop generating structures when the convergence criterion is met and the ``update`` method can still be called.
+As such, it is recommended that the user implement their own convergence check in their code.
+
+.. note::
+    The convergence criterion is not tested without explicitly setting the ``history_len`` parameter.
+    The convergence criterion is also not enforced even when set, it is simply a check for the user to use.
+
+
+
 Energy references
 -----------------
 

docs/source/tutorials/quick_guide.rst

@@ -54,7 +54,7 @@ The generator is then run, and the structures are written to an output file.
 Iterative structure search
 --------------------------
 
-Here is an example of how to run an iterative structure search with RAFFLE.
+Here is an example of how to run an iterative structure search with RAFFLE and check for convergence.
 
 .. code-block:: python
 
@@ -65,6 +65,7 @@ Here is an example of how to run an iterative structure search with RAFFLE.
     from mace.calculators import mace_mp
 
     generator = raffle_generator()
+    generator.set_history_len(10)
     mace = mace_mp(model="medium", dispersion=False, default_dtype="float32", device='cpu')
 
     host = read("host.xyz")
@@ -92,3 +93,7 @@ Here is an example of how to run an iterative structure search with RAFFLE.
 
         generator.update(structures)
         num_structures_old = len(structures)
+        if generator.is_converged():
+            break
+
+    write("output.xyz", structures)

example/python_pkg/graphene_grain_boundary_learn/DRAFFLE/learn.py

@@ -0,0 +1,261 @@
+from mace.calculators import mace_mp
+from ase.calculators.singlepoint import SinglePointCalculator
+from raffle.generator import raffle_generator
+from ase import build
+from ase.optimize import FIRE
+from ase.io import read, write
+import numpy as np
+import os
+from joblib import Parallel, delayed
+from ase.constraints import FixAtoms
+
+import logging
+logging.basicConfig(level=logging.DEBUG)
+
+# Function to relax a structure
+def process_structure(i, atoms, num_structures_old, calc_params, optimise_structure, iteration, calc):
+    # 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
+
+    # Calculate and save the initial energy per atom
+    energy_unrlxd = atoms.get_potential_energy() / len(atoms)
+    print(f"Initial energy per atom: {energy_unrlxd}")
+
+    # Optimise the structure if requested
+    if optimise_structure:
+        optimizer = FIRE(atoms, trajectory = f"traje{inew}.traj", logfile=f"optimisation{inew}.log")
+        try:
+            optimizer.run(fmax=0.05, steps=100)
+        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)
+
+    # Get the distance matrix
+    distances = atoms.get_all_distances(mic=True)
+
+    # Set self-distances (diagonal entries) to infinity to ignore them
+    np.fill_diagonal(distances, np.inf)
+
+    # Check if the minimum non-self distance is below 1.5
+    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"):
+        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 = mace_mp(**calc_params)
+
+    # set up the hosts
+    hosts = []
+    host = read("../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")
+    h2 = build.molecule("H2")
+    graphene.calc = calc
+    C_reference_energy = graphene.get_potential_energy() / len(graphene)
+    h2.calc = calc
+    H_reference_energy = h2.get_potential_energy() / len(h2)
+    generator.distributions.set_element_energies(
+        {
+            'C': C_reference_energy,
+            'H': H_reference_energy,
+        }
+    )
+
+    # set energy scale
+    generator.distributions.set_kBT(0.2)
+
+    # set the distribution function widths (2-body, 3-body, 4-body)
+    generator.distributions.set_width([0.04, np.pi/160.0, np.pi/160.0])
+
+    # set the initial database
+    initial_database = [graphene]
+    generator.distributions.create(initial_database, deallocate_systems=False)
+
+    # set the parameters for the generator
+    seed = 0
+
+    # 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
+    else:
+        open(energies_rlxd_filename, "w").close()
+
+    if os.path.exists(energies_unrlxd_filename):
+        with open(energies_unrlxd_filename, "w") as energy_file:
+            pass
+    else:
+        open(energies_unrlxd_filename, "w").close()
+
+    # initialise the number of structures generated
+    iter = -1
+    unrlxd_structures = []
+    rlxd_structures = []
+    num_structures_old = 0
+    optimise_structure = True
+    # start the iterations, loop over the hosts, then repeat the process X times
+    for ival in range(40):
+        for host in hosts:
+            print("setting host")
+            generator.set_host(host)
+            generator.set_bounds([[0.20, 0, 0.18], [0.30, 1, 0.33]])
+
+            iter += 1
+            print(f"Iteration {iter}")
+
+            # generate the structures
+            generator.generate(
+                num_structures = 5,
+                stoichiometry = { 'C': 15 },
+                seed = seed*1000+iter,
+                method_ratio = {"void": 0.1, "rand": 0.01, "walk": 0.25, "grow": 0.25, "min": 1.0},
+                verbose = 0,
+            )
+
+            # print the number of structures generated
+            print("Total number of structures generated: ", generator.num_structures)
+            generated_structures = generator.get_structures(calc)
+            num_structures_new = len(generated_structures)
+
+            # check if directory iteration[iter] exists, if not create it
+            iterdir = f"iteration{iter}/"
+            if not os.path.exists(iterdir):
+                os.makedirs(iterdir)
+            generator.print_settings(iterdir+"generator_settings.txt")
+
+            # set up list of energies
+            energy_unrlxd = np.zeros(num_structures_new - num_structures_old)
+            energy_rlxd = np.zeros(num_structures_new - num_structures_old)
+            for i in range(num_structures_new - num_structures_old):
+                write(iterdir+f"POSCAR_unrlxd_{i}", generated_structures[num_structures_old + i])
+                print(f"Structure {i} energy per atom: {generated_structures[num_structures_old + i].get_potential_energy() / len(generated_structures[num_structures_old + i])}")
+                # print(f"Structure {i} energy per unit area: {( generated_structures[num_structures_old + i].get_potential_energy() - Si_slab.get_potential_energy() - Ge_slab.get_potential_energy() ) / (2 * ( cell[0, 0] * cell[1, 1] ) )}")
+                # fix all hydrogen atoms to stop them from moving
+                c = FixAtoms(indices=[atom.index for atom in generated_structures[num_structures_old + i] if atom.symbol == 'H'])
+                generated_structures[num_structures_old + i].set_constraint(c)
+                unrlxd_structures.append(generated_structures[num_structures_old + i].copy())
+                unrlxd_structures[-1].calc = SinglePointCalculator(
+                    generated_structures[num_structures_old + i],
+                    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)(
+                delayed(process_structure)(i, generated_structures[i].copy(), num_structures_old, calc_params, optimise_structure, iteration=seed, calc=calc)
+                for i in range(num_structures_old, num_structures_new)
+            )
+
+            # Wait for all futures to complete
+            for j, result in enumerate(results):
+                generated_structures[num_structures_old + j], energy_unrlxd[j], energy_rlxd[j] = result
+                if generated_structures[num_structures_old + j] is None:
+                    print("Structure failed the checks")
+                    continue
+                rlxd_structures.append(generated_structures[num_structures_old + j].copy())
+                rlxd_structures[-1].calc = SinglePointCalculator(
+                    generated_structures[j+num_structures_old],
+                    energy=generated_structures[num_structures_old + j].get_potential_energy(),
+                    forces=generated_structures[num_structures_old + j].get_forces()
+                )
+            print("All futures completed")
+
+            # Remove structures that failed the checks
+            for j, atoms in reversed(list(enumerate(generated_structures))):
+                if j < num_structures_old:
+                    continue
+                if atoms is None:
+                    energy_unrlxd = np.delete(energy_unrlxd, j-num_structures_old)
+                    energy_rlxd = np.delete(energy_rlxd, j-num_structures_old)
+                    del generated_structures[j]
+                    del unrlxd_structures[j]
+                    del rlxd_structures[j]
+                    generator.remove_structure(j)
+            num_structures_new = len(generated_structures)
+
+            # write the structures to files
+            for i in range(num_structures_new - num_structures_old):
+                write(iterdir+f"POSCAR_{i}", generated_structures[num_structures_old + i])
+                print(f"Structure {i} energy per atom: {energy_rlxd[i]}")
+                # append energy per atom to the 'energies_unrlxd_filename' file
+                with open(energies_unrlxd_filename, "a") as energy_file:
+                    energy_file.write(f"{i+num_structures_old} {energy_unrlxd[i]}\n")
+                # append energy per atom to the 'energies_rlxd_filename' file
+                with open(energies_rlxd_filename, "a") as energy_file:
+                    energy_file.write(f"{i+num_structures_old} {energy_rlxd[i]}\n")
+
+            # update the distribution functions
+            print("Updating distributions")
+            generator.distributions.update(generated_structures[num_structures_old:], from_host=False, deallocate_systems=False)
+            if generator.distributions.is_converged(threshold=1e-3):
+                print("Converged")
+                break
+            else:
+              print("deltas", generator.distributions.history_deltas)
+
+            # print the new distribution functions to a file
+            print("Printing distributions")
+            generator.distributions.write_dfs(iterdir+"distributions.txt")
+            generator.distributions.write_2body(iterdir+"df2.txt")
+            generator.distributions.write_3body(iterdir+"df3.txt")
+            generator.distributions.write_4body(iterdir+"df4.txt")
+            generator.distributions.deallocate_systems()
+
+            # update the number of structures generated
+            num_structures_old = num_structures_new
+
+
+    # Write the final distribution functions to a file
+    generator.distributions.write_gdfs(f"gdfs_seed{seed}.txt")
+
+    # Read energies from the file
+    with open(energies_rlxd_filename, "r") as energy_file:
+        energies = energy_file.readlines()
+
+    # Parse and sort the energies
+    energies = [line.strip().split() for line in energies]
+    energies = sorted(energies, key=lambda x: float(x[1]))
+
+    # Write the sorted energies back to the file
+    with open(f"sorted_{energies_rlxd_filename}", "w") as energy_file:
+        for entry in energies:
+            energy_file.write(f"{int(entry[0])} {float(entry[1])}\n")
+
+    # Write the structures to files
+    write(f"unrlxd_structures_seed{seed}.traj", unrlxd_structures)
+    write(f"rlxd_structures_seed{seed}.traj", rlxd_structures)
+    print("All generated and relaxed structures written")
+
+    print("Learning complete")

example/python_pkg/graphene_grain_boundary_learn/POSCAR_graphene

@@ -0,0 +1,10 @@
+C2
+1.0
+   1.2336456308015411   -2.1367369110836258    0.0000000000000000
+   1.2336456308015411    2.1367369110836258    0.0000000000000000
+   0.0000000000000000    0.0000000000000000    10.0
+C
+2
+direct
+   0.0000000000000000    0.0000000000000000    0.2500000000000000 C0+
+   0.3333333333333330    0.6666666666666661    0.2500000000000000 C0+