From 1a840d4d071bb28fc7a7634028d17543ad7ecb70 Mon Sep 17 00:00:00 2001 From: Ned Taylor <71959356+nedtaylor@users.noreply.github.com> Date: Sat, 14 Jun 2025 13:26:01 +0100 Subject: [PATCH 01/11] Fix typo --- example/python_pkg/agox_runs/n-body_run/agox_run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/python_pkg/agox_runs/n-body_run/agox_run.py b/example/python_pkg/agox_runs/n-body_run/agox_run.py index 886df74..af97e9c 100644 --- a/example/python_pkg/agox_runs/n-body_run/agox_run.py +++ b/example/python_pkg/agox_runs/n-body_run/agox_run.py @@ -113,7 +113,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: From fed446860496abb786b2c2108f1f43eec4df8c35 Mon Sep 17 00:00:00 2001 From: Ned Taylor <71959356+nedtaylor@users.noreply.github.com> Date: Sat, 14 Jun 2025 14:09:25 +0100 Subject: [PATCH 02/11] Improve AGOX compatibility --- .../agox_runs/n-body_run/agox_run.py | 2 -- src/raffle/raffle.py | 34 ++++++++++++++++++- src/wrapper/f90wrap_mod_generator.f90 | 21 ++++++++++++ 3 files changed, 54 insertions(+), 3 deletions(-) diff --git a/example/python_pkg/agox_runs/n-body_run/agox_run.py b/example/python_pkg/agox_runs/n-body_run/agox_run.py index af97e9c..953b1e4 100644 --- a/example/python_pkg/agox_runs/n-body_run/agox_run.py +++ b/example/python_pkg/agox_runs/n-body_run/agox_run.py @@ -101,10 +101,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, diff --git a/src/raffle/raffle.py b/src/raffle/raffle.py index d4b72bb..d0a14d3 100644 --- a/src/raffle/raffle.py +++ b/src/raffle/raffle.py @@ -6,6 +6,8 @@ import warnings from ase import Atoms import re +from collections import Counter +from ase.data import chemical_symbols class Geom_Rw(f90wrap.runtime.FortranModule): @@ -2169,6 +2171,9 @@ def get_host(self): ase.Atoms The host structure for the generation. """ + is_defined = _raffle.f90wrap_raffle_generator_type__get__host_defined(self._handle) + if not is_defined: + return None output = \ _raffle.f90wrap_generator__get_host__binding__rgt(this=self._handle) output = f90wrap.runtime.lookup_class("raffle.basis").from_handle(output, \ @@ -2268,6 +2273,29 @@ def reset_grid(self): """ _raffle.f90wrap_generator__reset_grid__binding__raffle_generator_type(this=self._handle) + def get_bounds(self): + """ + Get the bounding box for the generation. + + Returns + ------- + bounds : list[list[float]] + The bounding box within which to constrain placement of atoms. + In the form [[a_min, b_min, c_min], [a_max, b_max, c_max]]. + Values given in direct (crystal) coordinates, ranging from 0 to 1. + """ + + array_ndim, array_type, array_shape, array_handle = \ + _raffle.f90wrap_raffle_generator_type__array__bounds(this=self._handle) + if array_handle in self._arrays: + bounds = self._arrays[array_handle] + else: + bounds = f90wrap.runtime.get_array(f90wrap.runtime.sizeof_fortran_t, + self._handle, + _raffle.f90wrap_raffle_generator_type__array__bounds) + self._arrays[array_handle] = bounds + return bounds + def set_bounds(self, bounds = None): """ Set the bounding box for the generation. @@ -2276,7 +2304,7 @@ def set_bounds(self, bounds = None): ---------- bounds : list[list[float]]: The bounding box within which to constrain placement of atoms. - In the form [[a_min, a_max], [b_min, b_max], [c_min, c_max]]. + In the form [[a_min, b_min, c_min], [a_max, b_max, c_max]]. Values given in direct (crystal) coordinates, ranging from 0 to 1. """ # check if bounds is a list of length 2 @@ -2403,6 +2431,10 @@ def generate(self, num = re.findall(r'\d+', s)[0] stoichiometry[element] = int(num) stoichiometry = Generator.stoichiometry_array(dict=stoichiometry) + elif isinstance(stoichiometry, list) and all(isinstance(i, int) for i in stoichiometry): + stoichiometry = [chemical_symbols[n] for n in stoichiometry] + stoichiometry = dict(Counter(stoichiometry)) + stoichiometry = Generator.stoichiometry_array(dict=stoichiometry) exit_code = _raffle.f90wrap_generator__generate__binding__rgt( this=self._handle, diff --git a/src/wrapper/f90wrap_mod_generator.f90 b/src/wrapper/f90wrap_mod_generator.f90 index 8b2d808..a4fddc8 100644 --- a/src/wrapper/f90wrap_mod_generator.f90 +++ b/src/wrapper/f90wrap_mod_generator.f90 @@ -283,6 +283,27 @@ end subroutine f90wrap_raffle_generator_type__set__num_structures !############################################################################### ! host handling !############################################################################### +subroutine f90wrap_raffle_generator_type__get__host_defined( & + this, f90wrap_host_is_defined & +) + use raffle__generator, only: raffle_generator_type + use raffle__geom_rw, only: basis_type + implicit none + type raffle_generator_type_ptr_type + type(raffle_generator_type), pointer :: p => NULL() + end type raffle_generator_type_ptr_type + type basis_type_ptr_type + type(basis_type), pointer :: p => NULL() + end type basis_type_ptr_type + integer, intent(in) :: this(2) + type(raffle_generator_type_ptr_type) :: this_ptr + logical, intent(out) :: f90wrap_host_is_defined + type(basis_type_ptr_type) :: host_ptr + + this_ptr = transfer(this, this_ptr) + f90wrap_host_is_defined = this_ptr%p%distributions%host_system%defined +end subroutine f90wrap_raffle_generator_type__get__host_defined + subroutine f90wrap_raffle_generator_type__get__host(this, f90wrap_host) use raffle__generator, only: raffle_generator_type use raffle__geom_rw, only: basis_type From 13b4a0f28ab9593e32370b3e220990d64191784e Mon Sep 17 00:00:00 2001 From: Ned Taylor <71959356+nedtaylor@users.noreply.github.com> Date: Sun, 15 Jun 2025 06:12:31 +0100 Subject: [PATCH 03/11] Handle zero volume bounds --- .../agox_runs/n-body_run/agox_run.py | 23 +++++++------------ src/fortran/lib/mod_generator.f90 | 6 +++++ src/raffle/raffle.py | 10 ++++++++ 3 files changed, 24 insertions(+), 15 deletions(-) diff --git a/example/python_pkg/agox_runs/n-body_run/agox_run.py b/example/python_pkg/agox_runs/n-body_run/agox_run.py index 953b1e4..0265cc2 100644 --- a/example/python_pkg/agox_runs/n-body_run/agox_run.py +++ b/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,6 +26,7 @@ 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 @@ -47,6 +34,12 @@ def raffle_agox(train=True, n_dist=4, index=42): 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} diff --git a/src/fortran/lib/mod_generator.f90 b/src/fortran/lib/mod_generator.f90 index 78accad..10abcf7 100644 --- a/src/fortran/lib/mod_generator.f90 +++ b/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() diff --git a/src/raffle/raffle.py b/src/raffle/raffle.py index d0a14d3..dac6813 100644 --- a/src/raffle/raffle.py +++ b/src/raffle/raffle.py @@ -2315,6 +2315,16 @@ def set_bounds(self, bounds = None): if not ( isinstance(bound, list) or isinstance(bound, numpy.ndarray) ) or len(bound) != 3: raise ValueError("each element of bounds must be a list of length 3") + # assert that bounds are in the range [0, 1] + for bound in bounds: + for value in bound: + if not (0.0 <= value <= 1.0): + raise ValueError("bounds must be in the range [0, 1]") + + # assert that a_max > a_min, b_max > b_min, c_max > c_min + if not (bounds[0][0] < bounds[1][0] and bounds[0][1] < bounds[1][1] and bounds[0][2] < bounds[1][2]): + raise ValueError("bounds must be in the form [[a_min, b_min, c_min], [a_max, b_max, c_max]] with a_max > a_min, b_max > b_min, c_max > c_min") + _raffle.f90wrap_generator__set_bounds__binding__rgt(this=self._handle, \ bounds=bounds) From 009a93b2f55786b4b2aa405c7ecb08ea589b4905 Mon Sep 17 00:00:00 2001 From: Ned Taylor <71959356+nedtaylor@users.noreply.github.com> Date: Sun, 15 Jun 2025 06:32:21 +0100 Subject: [PATCH 04/11] Fix segmentation fault --- src/fortran/lib/mod_cache.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fortran/lib/mod_cache.f90 b/src/fortran/lib/mod_cache.f90 index 579e9a9..fba52a1 100644 --- a/src/fortran/lib/mod_cache.f90 +++ b/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 From 5121d03a59e07d80fa681a723cfc01f9fcea445c Mon Sep 17 00:00:00 2001 From: Ned Taylor <71959356+nedtaylor@users.noreply.github.com> Date: Sun, 15 Jun 2025 07:36:42 +0100 Subject: [PATCH 05/11] Fix warnings --- src/fortran/lib/mod_misc_maths.f90 | 2 +- src/wrapper/f90wrap_mod_generator.f90 | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/fortran/lib/mod_misc_maths.f90 b/src/fortran/lib/mod_misc_maths.f90 index b410e36..5ca75df 100644 --- a/src/fortran/lib/mod_misc_maths.f90 +++ b/src/fortran/lib/mod_misc_maths.f90 @@ -47,7 +47,7 @@ pure function triangular_number(n) result(output) integer, intent(in) :: n !! The index of the triangular number to return. - real(real32) :: output + integer :: output !! The nth triangular number. output = n * ( n + 1 ) / 2 diff --git a/src/wrapper/f90wrap_mod_generator.f90 b/src/wrapper/f90wrap_mod_generator.f90 index a4fddc8..fb68de8 100644 --- a/src/wrapper/f90wrap_mod_generator.f90 +++ b/src/wrapper/f90wrap_mod_generator.f90 @@ -298,7 +298,6 @@ subroutine f90wrap_raffle_generator_type__get__host_defined( & integer, intent(in) :: this(2) type(raffle_generator_type_ptr_type) :: this_ptr logical, intent(out) :: f90wrap_host_is_defined - type(basis_type_ptr_type) :: host_ptr this_ptr = transfer(this, this_ptr) f90wrap_host_is_defined = this_ptr%p%distributions%host_system%defined From b9f76e6c56ce9b42fa00a56e4c67916fc83c2fcc Mon Sep 17 00:00:00 2001 From: Ned Taylor <71959356+nedtaylor@users.noreply.github.com> Date: Sun, 15 Jun 2025 07:36:52 +0100 Subject: [PATCH 06/11] Improve file reading and writing --- src/fortran/lib/mod_distribs_container.f90 | 140 ++++++++++++++------- 1 file changed, 94 insertions(+), 46 deletions(-) diff --git a/src/fortran/lib/mod_distribs_container.f90 b/src/fortran/lib/mod_distribs_container.f90 index fafc8c8..ef9575d 100644 --- a/src/fortran/lib/mod_distribs_container.f90 +++ b/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,6 +840,8 @@ subroutine read_gdfs(this, file) ! Local variables integer :: unit !! File unit. + integer :: iostat + !! I/O status. integer :: i !! Loop index. @@ -843,7 +849,7 @@ subroutine read_gdfs(this, file) !! 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. From c66a5683a2a42c0e658a5a0d541284613febb5eb Mon Sep 17 00:00:00 2001 From: Ned Taylor <71959356+nedtaylor@users.noreply.github.com> Date: Sun, 15 Jun 2025 07:37:17 +0100 Subject: [PATCH 07/11] Add unit tests --- test/CMakeLists.txt | 3 +- test/test_cache.f90 | 103 +++++++++++++++ test/test_distribs_container.f90 | 215 +++++++++++++++++++++++++++++++ test/test_tools_infile.f90 | 3 +- 4 files changed, 321 insertions(+), 3 deletions(-) create mode 100644 test/test_cache.f90 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7ef750b..02401fe 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -3,6 +3,7 @@ foreach(execid misc misc_maths misc_linalg + cache element_utils geom_rw geom_utils @@ -29,4 +30,4 @@ foreach(execid endif() add_test(NAME test_${execid} COMMAND test_${execid}) -endforeach() \ No newline at end of file +endforeach() diff --git a/test/test_cache.f90 b/test/test_cache.f90 new file mode 100644 index 0000000..3fddccb --- /dev/null +++ b/test/test_cache.f90 @@ -0,0 +1,103 @@ +program test_cache + use raffle__io_utils + use raffle__cache + use raffle__constants, only: real32 + implicit none + + logical :: success = .true. + + call test_retrieve_unallocated(success) + call test_store_retrieve_probability_density(success) + + !----------------------------------------------------------------------------- + ! check for any failed tests + !----------------------------------------------------------------------------- + write(*,*) "----------------------------------------" + if(success)then + write(*,*) 'test_cache passed all tests' + else + write(0,*) 'test_cache failed one or more tests' + stop 1 + end if + +contains + + subroutine test_store_retrieve_probability_density(success) + implicit none + logical, intent(inout) :: success + real(real32), allocatable :: test_data(:,:), retrieved_data(:,:) + integer :: i, j + + ! Create test data - a 3x3 matrix with values + allocate(test_data(3,3)) + do i = 1, 3 + do j = 1, 3 + test_data(i,j) = real(i*10 + j, real32) + end do + end do + + ! Store the test data in the cache + call store_probability_density(test_data) + + ! Retrieve the data from the cache + retrieved_data = retrieve_probability_density() + + ! Check dimensions match + if (.not. all(shape(retrieved_data) == shape(test_data))) then + write(0,*) 'test_store_retrieve_probability_density: array shapes do not match' + success = .false. + return + end if + + ! Check all values match + do i = 1, 3 + do j = 1, 3 + if (abs(retrieved_data(i,j) - test_data(i,j)) .gt. 1.e-6_real32) then + write(0,*) & + 'test_store_retrieve_probability_density: values do not match at', & + i, j, retrieved_data(i,j), test_data(i,j) + success = .false. + return + end if + end do + end do + + ! Clean up + if (allocated(test_data)) deallocate(test_data) + if (allocated(retrieved_data)) deallocate(retrieved_data) + + write(*,*) 'test_store_retrieve_probability_density: PASSED' + end subroutine test_store_retrieve_probability_density + + subroutine test_retrieve_unallocated(success) + implicit none + logical, intent(inout) :: success + real(real32), allocatable :: retrieved_data(:,:) + + ! Try to retrieve data when cache is not allocated + retrieved_data = retrieve_probability_density() + + ! Verify returned array is a scalar with value 0 + if (size(retrieved_data) .ne. 1) then + write(0,*) & + 'test_retrieve_unallocated: returned array should be scalar but has size', & + size(retrieved_data) + success = .false. + return + end if + + if (abs(retrieved_data(1,1)) .gt. 1.e-6_real32) then + write(0,*) & + 'test_retrieve_unallocated: returned value should be 0 but is', & + retrieved_data(1,1) + success = .false. + return + end if + + ! Clean up + if (allocated(retrieved_data)) deallocate(retrieved_data) + + write(*,*) 'test_retrieve_unallocated: PASSED' + end subroutine test_retrieve_unallocated + +end program test_cache diff --git a/test/test_distribs_container.f90 b/test/test_distribs_container.f90 index 45aafea..0829073 100644 --- a/test/test_distribs_container.f90 +++ b/test/test_distribs_container.f90 @@ -83,9 +83,11 @@ program test_distribs_container call test_set_radius_distance_tol(success) call test_create([basis_graphite], success) call test_update([basis_graphite, basis_diamond], success) + call test_read_write_gdfs([basis_graphite, basis_diamond], success) call test_create([basis_diamond, basis_mgo], success) call test_update([basis_diamond, basis_mgo], success) + call test_read_write_gdfs([basis_diamond, basis_mgo], success) ! call test_write_read(success) ! call test_write_2body(success) ! call test_write_3body(success) @@ -232,6 +234,41 @@ subroutine test_set_history_len(success) success & ) + ! Call the subroutine to set the history length when already allocated + history_len = 3 + call distribs_container%set_history_len(history_len) + + ! Check if the history length was set correctly + call assert( & + distribs_container%history_len .eq. history_len, & + "History length was not set correctly after reallocation", & + success & + ) + + call assert( & + allocated(distribs_container%history_deltas) .and. & + size(distribs_container%history_deltas, 1) .eq. history_len, & + "History delta list was not reallocated", & + success & + ) + + ! Call the subroutine to set the history length to zero + history_len = 0 + call distribs_container%set_history_len(history_len) + + ! Check if the history length was set correctly + call assert( & + distribs_container%history_len .eq. history_len, & + "History length was not set to zero", & + success & + ) + + call assert( & + .not.allocated(distribs_container%history_deltas), & + "History deltas was not deallocated", & + success & + ) + end subroutine test_set_history_len subroutine test_is_converged(success) @@ -392,6 +429,7 @@ subroutine test_create(basis, success) type(distribs_container_type) :: distribs_container type(basis_type), dimension(size(basis,1)) :: basis_list character(len=3), dimension(:), allocatable :: elements + real(real32), dimension(:), allocatable :: energy_above_hull ! Initialise basis_list do i = 1, size(basis,1) @@ -584,6 +622,29 @@ subroutine test_create(basis, success) success & ) + allocate(energy_above_hull(10)) + energy_above_hull = 0._real32 + + write(*,*) "Testing energy_above_hull error handling" + call distribs_container%create( & + basis_list, deallocate_systems=.false., & + energy_above_hull_list=energy_above_hull & + ) + write(*,*) "Handled error: energy_above_hull provided but weight_by_hull is false" + distribs_container%weight_by_hull = .true. + call distribs_container%create( & + basis_list, deallocate_systems=.false., & + energy_above_hull_list=energy_above_hull & + ) + write(*,*) "Handled error: energy_above_hull and basis_list sizes do not match" + deallocate(energy_above_hull) + allocate(energy_above_hull(size(basis_list,1))) + energy_above_hull = 0._real32 + call distribs_container%create( & + basis_list, deallocate_systems=.false., & + energy_above_hull_list=energy_above_hull & + ) + end subroutine test_create subroutine test_update(basis, success) @@ -597,6 +658,7 @@ subroutine test_update(basis, success) type(distribs_container_type) :: distribs_container type(basis_type), dimension(size(basis,1)) :: basis_list character(len=3), dimension(:), allocatable :: elements + real(real32), dimension(:), allocatable :: energy_above_hull ! Initialise basis_list do i = 1, size(basis,1) @@ -652,6 +714,42 @@ subroutine test_update(basis, success) success & ) + ! Check history length + call distribs_container%set_history_len(10) + call distribs_container%update( & + basis_list, deallocate_systems=.false. & + ) + call assert( & + distribs_container%history_deltas(1) .ne. huge(0._real32), & + "History delta not set on update", & + success & + ) + + allocate(energy_above_hull(10)) + energy_above_hull = 0._real32 + + write(*,*) "Testing energy_above_hull error handling" + call distribs_container%update( & + basis_list, deallocate_systems=.false., & + energy_above_hull_list=energy_above_hull & + ) + write(*,*) "Handled error: energy_above_hull provided but weight_by_hull is false" + distribs_container%weight_by_hull = .true. + call distribs_container%update( & + basis_list, deallocate_systems=.false., & + energy_above_hull_list=energy_above_hull & + ) + write(*,*) "Handled error: energy_above_hull and basis_list sizes do not match" + deallocate(energy_above_hull) + allocate(energy_above_hull(size(basis_list,1))) + energy_above_hull = 0._real32 + call distribs_container%update( & + basis_list, deallocate_systems=.false., & + energy_above_hull_list=energy_above_hull, & + from_host = .true. & + ) + write(*,*) "Handled error: from_host and energy_above_hull provided" + end subroutine test_update subroutine test_add(basis, success) @@ -909,6 +1007,123 @@ subroutine test_get_bin(success) end subroutine test_get_bin + subroutine test_read_write_gdfs(basis, success) + implicit none + type(basis_type), dimension(:), intent(in) :: basis + logical, intent(inout) :: success + + type(distribs_container_type) :: distribs_container, distribs_container_read + character(len=100) :: filename + + ! Set up a test filename + filename = 'test_gdfs.dat' + + ! Create a dummy distribs_container + call distribs_container%set_history_len(10) + call distribs_container%set_width([0.1_real32, 0.1_real32, 0.1_real32]) + call distribs_container%set_sigma([0.2_real32, 0.2_real32, 0.2_real32]) + call distribs_container%create([basis], deallocate_systems=.false.) + + ! Write the GDFs to a file + call distribs_container%write_gdfs(filename) + + ! Check if the file exists + inquire(file=filename, exist=success) + call assert(success, "GDF file was not created", success) + + ! Read the GDFs from the file + call distribs_container_read%read_gdfs(filename) + + ! Check if the GDFs were read correctly + call assert( & + allocated(distribs_container_read%gdf%df_2body), & + "2-body GDF was not read correctly", & + success & + ) + call assert( & + allocated(distribs_container_read%gdf%df_3body), & + "3-body GDF was not read correctly", & + success & + ) + call assert( & + allocated(distribs_container_read%gdf%df_4body), & + "4-body GDF was not read correctly", & + success & + ) + + ! Check if the read width, sigma, and history length match + call assert( & + all( & + abs( & + distribs_container_read%width - & + distribs_container%width & + ) .lt. 1.E-6_real32 & + ), & + "Width was not read correctly", & + success & + ) + call assert( & + all( & + abs( & + distribs_container_read%sigma - & + distribs_container%sigma & + ) .lt. 1.E-6_real32 & + ), & + "Sigma was not read correctly", & + success & + ) + call assert( & + distribs_container_read%history_len .eq. distribs_container%history_len, & + "History length was not read correctly", & + success & + ) + call assert( & + all( & + abs( 1._real32 - & + distribs_container_read%history_deltas / & + distribs_container%history_deltas & + ) .lt. 1.E-4_real32 & + ), & + "History deltas were not read correctly", & + success & + ) + call assert( & + all( & + abs( & + distribs_container_read%cutoff_min - & + distribs_container%cutoff_min & + ) .lt. 1.E-6_real32 & + ), & + "Cutoff min was not read correctly", & + success & + ) + call assert( & + all( & + abs( & + distribs_container_read%cutoff_max - & + distribs_container%cutoff_max & + ) .lt. 1.E-6_real32 & + ), & + "Cutoff max was not read correctly", & + success & + ) + call assert( & + all( & + abs( & + distribs_container_read%radius_distance_tol - & + distribs_container%radius_distance_tol & + ) .lt. 1.E-6_real32 & + ), & + "Radius distance tolerance was not read correctly", & + success & + ) + call assert( & + all( distribs_container_read%nbins .eq. distribs_container%nbins ), & + "Number of bins was not read correctly", & + success & + ) + end subroutine test_read_write_gdfs + !############################################################################### subroutine assert(condition, message, success) diff --git a/test/test_tools_infile.f90 b/test/test_tools_infile.f90 index 12c062c..6ac923e 100644 --- a/test/test_tools_infile.f90 +++ b/test/test_tools_infile.f90 @@ -6,7 +6,6 @@ program test_tools_infile integer :: ival = 0 logical :: ltmp1 - real(real32) :: rtmp1 character(256) :: stmp1 character(256) :: line @@ -80,4 +79,4 @@ program test_tools_infile -end program test_tools_infile \ No newline at end of file +end program test_tools_infile From 14eea1fe76464c064185c7050917c37fc4bed4c9 Mon Sep 17 00:00:00 2001 From: Ned Taylor <71959356+nedtaylor@users.noreply.github.com> Date: Sun, 15 Jun 2025 08:00:02 +0100 Subject: [PATCH 08/11] Improve AGOX documentation --- docs/source/tutorials/agox_tutorial.rst | 55 +++++++++++++++++-- docs/source/tutorials/parameters_tutorial.rst | 2 + 2 files changed, 53 insertions(+), 4 deletions(-) diff --git a/docs/source/tutorials/agox_tutorial.rst b/docs/source/tutorials/agox_tutorial.rst index c253766..c83e801 100644 --- a/docs/source/tutorials/agox_tutorial.rst +++ b/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/ diff --git a/docs/source/tutorials/parameters_tutorial.rst b/docs/source/tutorials/parameters_tutorial.rst index 2537ac0..0ebfdca 100644 --- a/docs/source/tutorials/parameters_tutorial.rst +++ b/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 ----------------- From bf62b7170c2fbd0025a0651ab30361ee007a0b43 Mon Sep 17 00:00:00 2001 From: Ned Taylor <71959356+nedtaylor@users.noreply.github.com> Date: Sun, 15 Jun 2025 08:10:48 +0100 Subject: [PATCH 09/11] Update AGOX examples --- .../python_pkg/agox_runs/Si-Ge_run/agox_run.py | 1 - example/python_pkg/agox_runs/Si-Ge_run/host.traj | Bin 0 -> 5480 bytes .../graphene_grain_boundary_run/agox_run.py | 1 - 3 files changed, 2 deletions(-) create mode 100644 example/python_pkg/agox_runs/Si-Ge_run/host.traj diff --git a/example/python_pkg/agox_runs/Si-Ge_run/agox_run.py b/example/python_pkg/agox_runs/Si-Ge_run/agox_run.py index b49b5d5..c94dd97 100644 --- a/example/python_pkg/agox_runs/Si-Ge_run/agox_run.py +++ b/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, diff --git a/example/python_pkg/agox_runs/Si-Ge_run/host.traj b/example/python_pkg/agox_runs/Si-Ge_run/host.traj new file mode 100644 index 0000000000000000000000000000000000000000..8564c9885c5936d8a6fea88e042db1e99a74b061 GIT binary patch literal 5480 zcmeH~?Q2ta7=~*Vwdnj11ubeh#yTdZhm$l(qjt?vSHXSQ@C}7pG>xgmrHv%9F?0pR zp_{hg6btGEsrXt2k&38DaI>OS+=nvR+|;3^D88bMiAvosFINXo{10M3VJTwS&#;#(W=$HI|>rVTB~{0+%YBRf6gGH=FY|Nl7ivkxtM|3A)5 z5&71`Vdv?h%1Qm+*$o@Y&Z)WWPU}BcK7E?3n%(u_hkXU6|JotjF*R=_-~4=T*I(X6 zkKSG$c50Lvk3apSR89D<>r;BYbz{vdJNL%FR#WoiAGtYig?ZkRW?^Mc6}e=0CLr**bZQ&UeC z&tJAuH_uN0abUN4?k9U+#t*0o|L~h9;qTRR|FOCMM4OtDC*MdMuB_f;nDFHHSj<~@ z3mqmt5m^doul+wYA&*y(lW?0?$xSVhkM!?Qo+OP+r4fu|qoLqG8N(~soI#|NH# z=|es|KIFs8{l|y<#eU)=_aB~ppdZQ0{^7}&KI}jK_>eDo`oRaDexwimz~e(dk|!S@ zc=DwW`SAFV4_}!5{g&(~t&iBsO`_A|NEc>R}+7yUHgW*Q&EX{7ytiEVq z<^M;l<+{VE8&-8 Date: Sun, 15 Jun 2025 08:11:19 +0100 Subject: [PATCH 10/11] Bump version number --- fpm.toml | 2 +- src/fortran/lib/mod_io_utils.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/fpm.toml b/fpm.toml index 788e755..84267b3 100644 --- a/fpm.toml +++ b/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" diff --git a/src/fortran/lib/mod_io_utils.F90 b/src/fortran/lib/mod_io_utils.F90 index 119f093..e1dab98 100644 --- a/src/fortran/lib/mod_io_utils.F90 +++ b/src/fortran/lib/mod_io_utils.F90 @@ -7,7 +7,7 @@ module raffle__io_utils logical :: test_error_handling = .false. logical :: suppress_warnings = .false. - character(len=*), parameter :: raffle__version__ = "1.1.0" + character(len=*), parameter :: raffle__version__ = "1.1.1" private From bfeffd9ed56854923c41d0f48f08e3363fedae05 Mon Sep 17 00:00:00 2001 From: Ned Taylor <71959356+nedtaylor@users.noreply.github.com> Date: Sun, 15 Jun 2025 08:13:31 +0100 Subject: [PATCH 11/11] Fix link --- docs/source/tutorials/agox_tutorial.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/tutorials/agox_tutorial.rst b/docs/source/tutorials/agox_tutorial.rst index c83e801..1352841 100644 --- a/docs/source/tutorials/agox_tutorial.rst +++ b/docs/source/tutorials/agox_tutorial.rst @@ -189,4 +189,4 @@ Other arguments exist that are not specific to the RAFFLE generator, but are use - ``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/ +Documentation of these arguments can be found in the AGOX documentation: https://agox.gitlab.io/agox/generators/generators.html