v1.1.1 Improve AGOX compatibility

docs/source/tutorials/agox_tutorial.rst

@@ -14,6 +14,9 @@ The example script can be found in the following directory:
 
     raffle/example/python_pkg/agox_runs/Si-Ge_run/agox_run.py
 
+Requirements
+------------
+
 This script utilises the `hex` slurm experiment management package.
 The script is designed to run on a slurm cluster, but can also be run locally by removing the slurm-specific code.
 `hex` can be installed via pip:
@@ -35,7 +38,9 @@ To use the RAFFLE generator, you need to install the `raffle_generator` branch o
     cd agox
     pip install -e .
 
-The RAFFLE generator differs slightly from the standard AGOX generator, as it requires a host structure to be provided directly to the generator.
+The RAFFLE generator differs slightly from the standard AGOX generator.
+If using `from_host` in RAFFLE, then the host structure must be provided to the generator directly in addition to the environment.
+Otherwise, the host structure (without information regarding the system's energy) is provided automatically by the AGOX framework.
 
 An environment must be set up for the search, this defines the host structure, the stoichiometry of atoms to be added, and the bounding box for the search.
 The bounding box must first be defined in terms of a lower left and upper right corner, which can be done using the `bounds_to_confinement` function.
@@ -93,15 +98,17 @@ The RAFFLE generator can then be set up using the `RaffleGenerator` class in AGO
 
     generator = RaffleGenerator(
         **environment.get_confinement(),
-        host = template,
+        element_energies =  {
+            'Si': -4.0,
+            'Ge': -3.5
+        }, # example element energies
         database = database,
-        environment = environment,
-        species = symbols,
         n_structures = 5,
         ...
     )
 
 This sets up the RAFFLE generator to generate 5 structures each iteration, using the host structure and the environment defined earlier.
+A more extensive list of arguments specific to the RAFFLE generator can be found at the end of this tutorial in the section :ref:`raffle_generator_arguments`.
 
 Evaluators and structure filters can be set up as usual in AGOX.
 For example, to set up an evaluator to perform structural optimisation, and a pre- and post-process filter that removes structures with bondlengths less than a certain value, you can use the following code:
@@ -143,3 +150,43 @@ Finally, the AGOX search can be set up and run using the `AGOX` class in AGOX:
 
     ## Run the AGOX search for N_iterations
     agox.run(N_iterations=40)
+
+
+.. _raffle_generator_arguments:
+
+RAFFLE generator specific arguments
+-----------------------------------
+
+The RAFFLE generator has several specific arguments that can be set to control the generation of structures.
+The required argument for the RAFFLE generator is:
+
+- ``element_energies``: A dictionary of element energies, taking the form ``{'Si': -4.0, 'Ge': -3.5}``. These are reference energies for the elements in the system, similar to chemical potentials.
+
+More information on this can be found in :ref:`element-energies`.
+
+Other optional arguments that can be set include:
+
+- ``n_structures``: The number of structures to generate in each iteration.
+- ``host``: The host structure to use for the generation. This can be an ASE Atoms object.
+- ``kBT``: The weighting factor for scaling the importance of different atomic features based on their system's relative energy.
+- ``history_len``: The length of the history for tracking change in the generalised descriptor.
+- ``width``: The width of the Gaussian functions used in the distribution functions (list of three floats, for 2-, 3-, and 4-body interactions).
+- ``sigma``: The standard deviation of the Gaussian functions used in the distribution functions (list of three floats, for 2-, 3-, and 4-body interactions).
+- ``cutoff_min``: The minimum cutoff for the Gaussian functions (list of three floats, for 2-, 3-, and 4-body interactions).
+- ``cutoff_max``: The maximum cutoff for the Gaussian functions (list of three floats, for 2-, 3-, and 4-body interactions).
+- ``radius_distance_tol``: The radius distance tolerance for the element-pair covalent radii (list of four floats, for 3-body min, 3-body max, 4-body min, and 4-body max interactions).
+- ``transfer_data``: A list of ASE Atoms objects used to initialise the generalised descriptor.
+- ``method_ratio``: The ratio of placement methods to use in the generation (dictionary with keys `void`, `rand`, `walk`, `grow`, and `min`, and values as floats representing the relative importance of each method).
+- ``deallocate_systems``: A boolean flag to indicate whether to deallocate the individual distribution functions of systems after they have been combined into the generalised descriptor.
+- ``from_host``: A boolean flag to indicate whether to represent the energies with respect to the host structure.
+- ``max_walk_attempts``: The maximum number of attempts to place an atom using the `walk` method before giving up.
+- ``walk_step_size_coarse``: The initial/coarse step size for the `walk` method.
+- ``walk_step_size_fine``: The final/fine step size for the `walk` method.
+- ``grid_spacing``: The spacing of the grid used for the `void` and `min` methods.
+- ``seed``: A random seed for reproducibility.
+
+Other arguments exist that are not specific to the RAFFLE generator, but are used in the AGOX framework, such as:
+
+- ``database``: The database object to use for storing the generated structures.
+
+Documentation of these arguments can be found in the AGOX documentation: https://agox.gitlab.io/agox/generators/generators.html

docs/source/tutorials/parameters_tutorial.rst

@@ -111,6 +111,8 @@ To set the placement method ratios for a specific generation method, the user ca
     )
 
 
+.. _element-energies:
+
 Energy references
 -----------------
 

example/python_pkg/agox_runs/Si-Ge_run/agox_run.py

@@ -97,7 +97,6 @@ def raffle_agox(index=42, void=0.1, rand=0.01, walk=0.25, grow=0.25, min=1.0):
         host=template,
         database=database,
         environment=environment,
-        species=species_to_place,
         element_energies=element_energies,
         transfer_data=[Si_bulk,Ge_bulk],
         order=1,

example/python_pkg/agox_runs/graphene_grain_boundary_run/agox_run.py

@@ -101,7 +101,6 @@ def raffle_agox(index=42, void=0.1, rand=0.01, walk=0.25, grow=0.25, min=1.0):
         host=template,
         database=database,
         environment=environment,
-        species=species_to_place,
         element_energies=element_energies,
         transfer_data=[graphene],
         order=1,

example/python_pkg/agox_runs/n-body_run/agox_run.py

@@ -2,30 +2,16 @@
 
 matplotlib.use("Agg")
 
-import numpy as np
-from ase import Atoms
-from ase.optimize import BFGS
 
 from agox import AGOX
 from agox.databases import Database
 from agox.environments import Environment
 from agox.evaluators import LocalOptimizationEvaluator
 from agox.generators import RattleGenerator
-from agox.models.descriptors import SOAP
-from agox.models.GPR import SparseGPR
-from agox.models.GPR.kernels import RBF
-from agox.models.GPR.kernels import Constant as C
-from agox.postprocessors.ray_relax import ParallelRelaxPostprocess
-from agox.samplers import MetropolisSampler
 import shephex
-from ase import build
 
 @shephex.chain()
 def raffle_agox(train=True, n_dist=4, index=42):
-    # Manually set seed and database-index
-    seed = index
-    database_index = index
-    np.random.seed(seed)
 
     # Using argparse if e.g. using array-jobs on Slurm to do several independent searches.
     # from argparse import ArgumentParser
@@ -40,13 +26,20 @@ def raffle_agox(train=True, n_dist=4, index=42):
     # System & general settings:
     ##############################################################################
 
+    import numpy as np
     from ase import Atoms
     from ase.io import write
     from ase.calculators.singlepoint import SinglePointCalculator
     from agox.generators.raffle import RaffleGenerator, bounds_to_confinement
     from mace.calculators import mace_mp
     from agox.utils.replica_exchange.priors import get_prior
 
+
+    # Manually set seed and database-index
+    seed = index
+    database_index = index
+    np.random.seed(seed)
+
     chg = get_prior("chg")
     mace = get_prior("mace")
 
@@ -77,7 +70,7 @@ def raffle_agox(train=True, n_dist=4, index=42):
     #write("template.traj", template)
     print(indices)
 
-    confinement_corner, confinement_cell = bounds_to_confinement([[0, 0, 0.0], [1, 1, 1]],template)
+    confinement_corner, confinement_cell = bounds_to_confinement([[0, 0, 0], [1, 1, 1]],template)
 
     element_energies = {'C': C_reference_energy}
     species_to_place={"C" : indices}
@@ -101,10 +94,8 @@ def raffle_agox(train=True, n_dist=4, index=42):
 
 
     generator = RaffleGenerator(**environment.get_confinement(),
-        host=template,
         database=database,
         environment=environment,
-        species=species_to_place,
         element_energies=element_energies,
         transfer_data=[template_copy],
         order=1,
@@ -113,7 +104,7 @@ def raffle_agox(train=True, n_dist=4, index=42):
         sampler=None,
         train=train,
         seed=seed,
-        gen_params = {"void": 0.0, "rand": 0.1, "walk": 0.25, "grow": 0.25, "min": 1.0})
+        method_ratio = {"void": 0.0, "rand": 0.1, "walk": 0.25, "grow": 0.25, "min": 1.0})
 
     ##############################################################################
     # Search Settings:

fpm.toml

@@ -1,5 +1,5 @@
 name = "raffle"
-version = "1.1.0"
+version = "1.1.1"
 author = "Ned Thaddeus Taylor"
 maintainer = "n.t.taylor@exeter.ac.uk"
 description = "A Fortran library and executable for structure prediction at material interfaces"

src/fortran/lib/mod_cache.f90

@@ -23,7 +23,7 @@ function retrieve_probability_density() result(probability_density)
     real(real32), allocatable :: probability_density(:,:)
     if(.not.allocated(cached_probability_density)) then
        write(0,*) "Probability density not allocated. Returning zero array."
-       probability_density = 0._real32
+       allocate(probability_density(1,1), source = 0._real32)
     else
        allocate(probability_density, source = cached_probability_density)
     end if

src/fortran/lib/mod_distribs_container.f90

@@ -754,11 +754,15 @@ subroutine write_gdfs(this, file)
        return
     end if
     open(newunit=unit, file=file)
+    write(unit, '("# history_len ",I0)') this%history_len
+    write(fmt, '("(""# history_deltas""",I0,"(1X,ES0.8))")') &
+         size(this%history_deltas)
+    write(unit, fmt) this%history_deltas
     write(unit, '("# nbins",3(1X,I0))') this%nbins
-    write(unit, '("# width",3(1X,ES0.4))') this%width
-    write(unit, '("# sigma",3(1X,ES0.4))') this%sigma
-    write(unit, '("# cutoff_min",3(1X,ES0.4))') this%cutoff_min
-    write(unit, '("# cutoff_max",3(1X,ES0.4))') this%cutoff_max
+    write(unit, '("# width",3(1X,ES0.8))') this%width
+    write(unit, '("# sigma",3(1X,ES0.8))') this%sigma
+    write(unit, '("# cutoff_min",3(1X,ES0.8))') this%cutoff_min
+    write(unit, '("# cutoff_max",3(1X,ES0.8))') this%cutoff_max
     write(unit, '("# radius_distance_tol",4(1X,ES0.4))') &
          this%radius_distance_tol
     write(fmt, '("(""# "",A,",I0,"(1X,A))")') size(this%element_info)
@@ -836,14 +840,16 @@ subroutine read_gdfs(this, file)
     ! Local variables
     integer :: unit
     !! File unit.
+    integer :: iostat
+    !! I/O status.
 
     integer :: i
     !! Loop index.
     integer :: nspec
     !! Number of species.
     logical :: exist
     !! Boolean whether the file exists.
-    character(256) :: buffer, buffer1, buffer2
+    character(256) :: buffer, buffer1
     !! Buffer for reading lines.
 
     ! check if file exists
@@ -855,46 +861,90 @@ subroutine read_gdfs(this, file)
 
     ! read the file
     open(newunit=unit, file=file)
-    read(unit, *) buffer1, buffer2, this%nbins
-    read(unit, *) buffer1, buffer2, this%width
-    read(unit, *) buffer1, buffer2, this%sigma
-    read(unit, *) buffer1, buffer2, this%cutoff_min
-    read(unit, *) buffer1, buffer2, this%cutoff_max
-    read(unit, *) buffer1, buffer2, this%radius_distance_tol
-    read(unit, '(A)') buffer
-    nspec = icount(buffer(index(buffer,"elements")+8:))
-    if(allocated(this%element_info)) deallocate(this%element_info)
-    allocate(this%element_info(nspec))
-    read(buffer, *) buffer1, buffer2, this%element_info(:)%name
-    read(unit, *) buffer1, buffer2, this%element_info(:)%energy
-    do i = 1, nspec
-       call this%set_element_energy( &
-            this%element_info(i)%name, &
-            this%element_info(i)%energy &
-       )
-       call this%element_info(i)%set(this%element_info(i)%name)
+    do
+       read(unit, '(A)', iostat=iostat) buffer
+       if(iostat.ne.0) exit
+       buffer = trim(adjustl(buffer))
+       if(trim(buffer) .eq. "") cycle
+       if(buffer(1:1) .ne. "#") cycle
+       ! get the header
+       buffer = trim(adjustl(buffer(2:)))
+       if(trim(adjustl(buffer)) .eq. "2-body") exit
+       if(index(buffer, "history_len") .ne. 0) then
+          read(buffer, *) buffer1, this%history_len
+          call this%set_history_len(this%history_len)
+       else if(index(buffer, "history_deltas") .ne. 0) then
+          read(buffer, *) buffer1, this%history_deltas
+       else if(index(buffer, "nbins") .ne. 0) then
+          read(buffer, *) buffer1, this%nbins
+       else if(index(buffer, "width") .ne. 0) then
+          read(buffer, *) buffer1, this%width
+       else if(index(buffer, "sigma") .ne. 0) then
+          read(buffer, *) buffer1, this%sigma
+       else if(index(buffer, "cutoff_min") .ne. 0) then
+          read(buffer, *) buffer1, this%cutoff_min
+       else if(index(buffer, "cutoff_max") .ne. 0) then
+          read(buffer, *) buffer1, this%cutoff_max
+       else if(index(buffer, "radius_distance_tol") .ne. 0) then
+          read(buffer, *) buffer1, this%radius_distance_tol
+       else if(index(buffer, "elements") .ne. 0) then
+          nspec = icount(buffer(index(buffer,"elements")+8:))
+          if(allocated(this%element_info)) deallocate(this%element_info)
+          allocate(this%element_info(nspec))
+          read(buffer, *) buffer1, this%element_info(:)%name
+       else if(index(buffer, "energies") .ne. 0) then
+          read(buffer, *) buffer1, this%element_info(:)%energy
+          do i = 1, nspec
+             call this%set_element_energy( &
+                  this%element_info(i)%name, &
+                  this%element_info(i)%energy &
+             )
+             call this%element_info(i)%set(this%element_info(i)%name)
+          end do
+          call this%update_bond_info()
+          allocate(this%best_energy_per_species(nspec))
+          allocate(this%norm_3body(nspec))
+          allocate(this%norm_4body(nspec))
+          allocate(this%in_dataset_3body(nspec))
+          allocate(this%in_dataset_4body(nspec))
+          allocate(this%best_energy_pair(size(this%bond_info)))
+          allocate(this%norm_2body(size(this%bond_info)))
+          allocate(this%in_dataset_2body(size(this%bond_info)))
+       else if(index(buffer, "best_energy_per_element") .ne. 0) then
+          read(buffer, *) buffer1, this%best_energy_per_species
+       else if(index(buffer, "3-body_norm") .ne. 0) then
+          read(buffer, *) buffer1, this%norm_3body
+       else if(index(buffer, "4-body_norm") .ne. 0) then
+          read(buffer, *) buffer1, this%norm_4body
+       else if(index(buffer, "in_dataset_3body") .ne. 0) then
+          read(buffer, *) buffer1, this%in_dataset_3body
+       else if(index(buffer, "in_dataset_4body") .ne. 0) then
+          read(buffer, *) buffer1, this%in_dataset_4body
+       else if(index(buffer, "element_pairs") .ne. 0) then
+          read(buffer, *) buffer1, this%bond_info(:)%element(1)
+          read(buffer, *) buffer1, this%bond_info(:)%element(2)
+       else if(index(buffer, "radii") .ne. 0) then
+          read(buffer, *) buffer1, this%bond_info(:)%radius_covalent
+       else if(index(buffer, "best_energy_per_pair") .ne. 0) then
+          read(buffer, *) buffer1, this%best_energy_pair
+       else if(index(buffer, "3-body_norm") .ne. 0) then
+          read(buffer, *) buffer1, this%norm_3body
+       else if(index(buffer, "4-body_norm") .ne. 0) then
+          read(buffer, *) buffer1, this%norm_4body
+       else if(index(buffer, "in_dataset_3body") .ne. 0) then
+          read(buffer, *) buffer1, this%in_dataset_3body
+       else if(index(buffer, "in_dataset_4body") .ne. 0) then
+          read(buffer, *) buffer1, this%in_dataset_4body
+       else if(index(buffer, "2-body_norm") .ne. 0) then
+          read(buffer, *) buffer1, this%norm_2body
+       else if(index(buffer, "in_dataset_2body") .ne. 0) then
+          read(buffer, *) buffer1, this%in_dataset_2body
+       else
+          write(0,*) "Unknown header: ", trim(buffer)
+          cycle
+       end if
     end do
-    call this%update_bond_info()
-    allocate(this%best_energy_per_species(nspec))
-    allocate(this%norm_3body(nspec))
-    allocate(this%norm_4body(nspec))
-    allocate(this%in_dataset_3body(nspec))
-    allocate(this%in_dataset_4body(nspec))
-    read(unit, *) buffer1, buffer2, this%best_energy_per_species
-    read(unit, *) buffer1, buffer2, this%norm_3body
-    read(unit, *) buffer1, buffer2, this%norm_4body
-    read(unit, *) buffer1, buffer2, this%in_dataset_3body
-    read(unit, *) buffer1, buffer2, this%in_dataset_4body
-    read(unit, *)
-    allocate(this%best_energy_pair(size(this%bond_info)))
-    allocate(this%norm_2body(size(this%bond_info)))
-    allocate(this%in_dataset_2body(size(this%bond_info)))
-    read(unit, *) buffer1, buffer2, this%bond_info(:)%radius_covalent
-    read(unit, *) buffer1, buffer2, this%best_energy_pair
-    read(unit, *) buffer1, buffer2, this%norm_2body
-    read(unit, *) buffer1, buffer2, this%in_dataset_2body
-    read(unit, *)
-    read(unit, *)
+
     read(unit, *)
     allocate(this%gdf%df_2body(this%nbins(1),size(this%bond_info)))
     do i = 1, this%nbins(1)
@@ -2803,8 +2853,6 @@ function is_converged(this, threshold) result(converged)
     !! Convergence flag.
 
     ! Local variables
-    integer :: i, j
-    !! Loop index.
     real(real32) :: threshold_
     !! Threshold for convergence.
 

src/fortran/lib/mod_generator.f90

@@ -502,6 +502,12 @@ subroutine set_bounds(this, bounds)
     real(real32), dimension(2,3), intent(in) :: bounds
     !! Bounds for atom placement.
 
+    ! check if bounds has zero volume, if so, return
+    if( any(bounds(2,:) .le. bounds(1,:)) ) then
+       call stop_program("Bounds have zero volume")
+       return
+    end if
+
     this%bounds = bounds
     call this%set_grid()