diff --git a/docs/source/tutorials/generating_tutorial.rst b/docs/source/tutorials/generating_tutorial.rst index 8fdaa63..6d38e72 100644 --- a/docs/source/tutorials/generating_tutorial.rst +++ b/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 ------- diff --git a/docs/source/tutorials/parameters_tutorial.rst b/docs/source/tutorials/parameters_tutorial.rst index b39f2d0..9cd0299 100644 --- a/docs/source/tutorials/parameters_tutorial.rst +++ b/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 ----------------- diff --git a/docs/source/tutorials/quick_guide.rst b/docs/source/tutorials/quick_guide.rst index 5254bfb..21d8957 100644 --- a/docs/source/tutorials/quick_guide.rst +++ b/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) diff --git a/example/python_pkg/graphene_grain_boundary_learn/DRAFFLE/learn.py b/example/python_pkg/graphene_grain_boundary_learn/DRAFFLE/learn.py new file mode 100644 index 0000000..17a6139 --- /dev/null +++ b/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") diff --git a/example/python_pkg/graphene_grain_boundary_learn/POSCAR_graphene b/example/python_pkg/graphene_grain_boundary_learn/POSCAR_graphene new file mode 100644 index 0000000..a6619dc --- /dev/null +++ b/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+ diff --git a/example/python_pkg/graphene_grain_boundary_learn/POSCAR_host_gb b/example/python_pkg/graphene_grain_boundary_learn/POSCAR_host_gb new file mode 100644 index 0000000..fac8134 --- /dev/null +++ b/example/python_pkg/graphene_grain_boundary_learn/POSCAR_host_gb @@ -0,0 +1,105 @@ +C2+C2 + 1.000000000 + 10.000000000 0.000000000 0.000000000 + 0.000000000 17.271038055 0.000000000 + 0.000000000 0.000000000 27.684305191 +C H +82 15 +Direct + 0.249999985 0.857590735 0.000000000 + 0.249999985 0.929019332 0.077182270 + 0.249999985 0.000447869 0.154364484 + 0.249999985 0.714733541 0.000000000 + 0.249999985 0.786162138 0.077182270 + 0.249999985 0.857590735 0.154364503 + 0.249999985 0.571876407 0.000000000 + 0.249999985 0.643305004 0.077182270 + 0.249999985 0.714733541 0.154364522 + 0.249999985 0.429019302 0.000000000 + 0.249999985 0.500447869 0.077182270 + 0.249999985 0.571876466 0.154364503 + 0.249999985 0.286162138 0.000000000 + 0.249999985 0.357590735 0.077182276 + 0.249999985 0.429019272 0.154364503 + 0.249999985 0.143305004 0.000000000 + 0.249999985 0.214733586 0.077182270 + 0.249999985 0.286162138 0.154364484 + 0.249999985 0.000447861 0.000000000 + 0.249999985 0.071876436 0.077182261 + 0.249999985 0.143305019 0.154364503 + 0.249999985 0.929019272 0.025727424 + 0.249999985 0.000447869 0.102909669 + 0.249999985 0.071876407 0.180091925 + 0.249999985 0.786162198 0.025727415 + 0.249999985 0.857590795 0.102909679 + 0.249999985 0.929019392 0.180091906 + 0.249999985 0.643305004 0.025727436 + 0.249999985 0.714733601 0.102909679 + 0.249999985 0.786162138 0.180091925 + 0.249999985 0.500447869 0.025727432 + 0.249999985 0.571876466 0.102909679 + 0.249999985 0.643305004 0.180091925 + 0.249999985 0.357590765 0.025727420 + 0.249999985 0.429019332 0.102909679 + 0.249999985 0.500447929 0.180091925 + 0.249999985 0.214733586 0.025727432 + 0.249999985 0.286162168 0.102909679 + 0.249999985 0.357590735 0.180091925 + 0.249999985 0.071876436 0.025727432 + 0.249999985 0.143305019 0.102909679 + 0.249999985 0.214733601 0.180091925 + 0.249999985 0.326032639 0.498560566 + 0.249999985 0.826032877 0.320315826 + 0.249999985 0.201032996 0.364877030 + 0.249999985 0.576032877 0.409438234 + 0.249999985 0.951032758 0.453999476 + 0.249999985 0.076032877 0.498560604 + 0.249999985 0.576032937 0.320315826 + 0.249999985 0.951032937 0.364877030 + 0.249999985 0.326032877 0.409438196 + 0.249999985 0.701032758 0.453999476 + 0.249999985 0.826032758 0.498560604 + 0.249999985 0.326032966 0.320315788 + 0.249999985 0.701032996 0.364877030 + 0.249999985 0.076032877 0.409438196 + 0.249999985 0.451032758 0.453999476 + 0.249999985 0.576032758 0.498560604 + 0.249999985 0.076032981 0.320315788 + 0.249999985 0.451032996 0.364877030 + 0.249999985 0.826032877 0.409438196 + 0.249999985 0.201032758 0.453999400 + 0.249999985 0.409365892 0.498560604 + 0.249999985 0.909366190 0.320315788 + 0.249999985 0.284366131 0.364876992 + 0.249999985 0.659366250 0.409438234 + 0.249999985 0.034366131 0.453999400 + 0.249999985 0.159366131 0.498560642 + 0.249999985 0.659366190 0.320315788 + 0.249999985 0.034366250 0.364876992 + 0.249999985 0.409366250 0.409438234 + 0.249999985 0.784366131 0.453999400 + 0.249999985 0.909366131 0.498560642 + 0.249999985 0.409366310 0.320315826 + 0.249999985 0.784366310 0.364876992 + 0.249999985 0.159366250 0.409438196 + 0.249999985 0.534366250 0.453999438 + 0.249999985 0.659366131 0.498560680 + 0.249999985 0.159366310 0.320315788 + 0.249999985 0.534366310 0.364877030 + 0.249999985 0.909366190 0.409438196 + 0.249999985 0.284366250 0.453999438 + 0.250000000 0.044199999 0.532874050 + 0.250000000 0.191000000 0.532874050 + 0.250000000 0.294000000 0.532874050 + 0.250000000 0.441000015 0.532874050 + 0.250000000 0.544000030 0.532874050 + 0.250000000 0.690999985 0.532874050 + 0.250000000 0.794000030 0.532874050 + 0.250000000 0.940999985 0.532874050 + 0.250000000 0.000450000 0.960395375 + 0.250000000 0.143309996 0.960395375 + 0.250000000 0.286159992 0.960395375 + 0.250000000 0.429019988 0.960395375 + 0.250000000 0.571879983 0.960395375 + 0.250000000 0.714730024 0.960395375 + 0.250000000 0.857590020 0.960395375 diff --git a/src/fortran/lib/mod_distribs.f90 b/src/fortran/lib/mod_distribs.f90 index ca21f1a..5e17c15 100644 --- a/src/fortran/lib/mod_distribs.f90 +++ b/src/fortran/lib/mod_distribs.f90 @@ -17,12 +17,12 @@ module raffle__distribs element_database, element_bond_database implicit none - + private public :: distribs_base_type, distribs_type, get_distrib - - + + type :: distribs_base_type !! Base type for distribution functions. real(real32), dimension(:,:), allocatable :: df_2body @@ -31,6 +31,9 @@ module raffle__distribs !! 3-body distribution function. real(real32), dimension(:,:), allocatable :: df_4body !! 4-body distribution function. + contains + procedure, pass(this) :: compare + !! Compare this distribution function with another. end type distribs_base_type type, extends(distribs_base_type) :: distribs_type @@ -60,6 +63,7 @@ module raffle__distribs !! Weights for the 2-body and species distribution functions. contains procedure, pass(this) :: calculate + !! Calculate the distribution functions for the structure. end type distribs_type @@ -131,7 +135,7 @@ subroutine calculate(this, basis, & nbins, width, sigma, cutoff_min, cutoff_max, radius_distance_tol) !! Calculate the distribution functions for the container. !! - !! This procedure calculates the 2-, 3-, and 4-body distribution function + !! This procedure calculates the 2-, 3-, and 4-body distribution function !! for a given atomic structure (i.e. basis). implicit none @@ -223,7 +227,7 @@ subroutine calculate(this, basis, & else radius_distance_tol_ = [1.5_real32, 2.5_real32, 3._real32, 6._real32] end if - + !--------------------------------------------------------------------------- @@ -328,12 +332,12 @@ subroutine calculate(this, basis, & ], basis_extd%lat ) & ) bondlength = modu( vector ) - + if( & bondlength .lt. cutoff_min_(1) .or. & bondlength .gt. cutoff_max_(1) & ) cycle atom_loop - + ! add 2-body bond to store if within tolerances for 3-body ! distance if( & @@ -355,7 +359,7 @@ subroutine calculate(this, basis, & ! add 2-body bond to store if within tolerances for 4-body ! distance if( & - bondlength .ge. ( & + bondlength .ge. ( & bond_info(pair_index(is, js))%radius_covalent * & radius_distance_tol_(3) & ) .and. & @@ -374,7 +378,7 @@ subroutine calculate(this, basis, & itmp1 = itmp1 + 1 bondlength_list(itmp1) = bondlength distance(itmp1) = 1._real32 - + end associate end do atom_loop @@ -391,12 +395,12 @@ subroutine calculate(this, basis, & ) bondlength = modu( vector ) - + if( & bondlength .lt. cutoff_min_(1) .or. & bondlength .gt. cutoff_max_(1) & ) cycle image_loop - + ! add 2-body bond to store if within tolerances for 3-body ! distance if( & @@ -419,7 +423,7 @@ subroutine calculate(this, basis, & ! add 2-body bond to store if within tolerances for 4-body ! distance if( & - bondlength .ge. ( & + bondlength .ge. ( & bond_info(pair_index(is, js))%radius_covalent * & radius_distance_tol_(3) & ) .and. & @@ -438,7 +442,7 @@ subroutine calculate(this, basis, & itmp1 = itmp1 + 1 bondlength_list(itmp1) = bondlength distance(itmp1) = 1._real32 - + end associate end do image_loop @@ -696,4 +700,65 @@ function get_distrib(value_list, nbins, eta, width, cutoff_min, & end function get_distrib !############################################################################### -end module raffle__distribs \ No newline at end of file + +!############################################################################### + function compare(this, input) result(output) + !! Compare this distribution function with another. + implicit none + + ! Arguments + class(distribs_base_type), intent(in) :: this + !! Parent of the procedure. Instance of distribution functions container. + class(distribs_base_type), intent(in) :: input + !! Distribution function to compare with. + + ! Local variables + integer :: num_bins_2body, num_bins_3body, num_bins_4body + !! Number of bins for the distribution functions. + real(real32) :: diff_2body, diff_3body, diff_4body + !! Difference between the two distribution functions. + real(real32) :: output + !! Output comparison value (i.e. how much the two dfs overlap). + integer :: i + !! Loop index. + + + output = 0._real32 + + !--------------------------------------------------------------------------- + ! compare the 2-body distribution functions + !--------------------------------------------------------------------------- + num_bins_2body = size(this%df_2body, dim=1) + num_bins_3body = size(this%df_3body, dim=1) + num_bins_4body = size(this%df_4body, dim=1) + do i = 1, size(this%df_2body, dim=2) + if(any(abs(this%df_2body(:,i)).gt.1.E-6))then + diff_2body = sum( & + abs( this%df_2body(:,i) - input%df_2body(:,i) ) & + ) / num_bins_2body + output = output + diff_2body + end if + end do + + !--------------------------------------------------------------------------- + ! compare the 3-body distribution functions + !--------------------------------------------------------------------------- + do i = 1, size(this%df_3body, dim=2) + if(any(abs(this%df_3body(:,i)).gt.1.E-6))then + diff_3body = sum( & + abs( this%df_3body(:,i) - input%df_3body(:,i) ) & + ) / num_bins_3body + output = output + diff_3body + end if + if(any(abs(this%df_4body(:,i)).gt.1.E-6))then + diff_4body = sum( & + abs( this%df_4body(:,i) - input%df_4body(:,i) ) & + ) / num_bins_4body + output = output + diff_4body + end if + end do + + end function compare +!############################################################################### + +end module raffle__distribs diff --git a/src/fortran/lib/mod_distribs_container.f90 b/src/fortran/lib/mod_distribs_container.f90 index d03a4b7..a2de17e 100644 --- a/src/fortran/lib/mod_distribs_container.f90 +++ b/src/fortran/lib/mod_distribs_container.f90 @@ -38,6 +38,10 @@ module raffle__distribs_container !! Number of evaluated systems still allocated. real(real32) :: kBT = 0.2_real32 !! Boltzmann constant times temperature. + integer :: history_len = 0 + !! Length of the history for the distribution functions. + real(real32), dimension(:), allocatable :: history_deltas + !! History of the changes in the distribution functions. logical :: weight_by_hull = .false. !! Boolean whether to weight the distribution functions by the energy !! above the hull. If false, the formation energy from the element @@ -104,6 +108,8 @@ module raffle__distribs_container !! Set the maximum cutoff for the 2-, 3-, and 4-body. procedure, pass(this) :: set_radius_distance_tol !! Set the tolerance for the distance between atoms for 3- and 4-body. + procedure, pass(this) :: set_history_len + !! Set the length of the history for the distribution functions. procedure, pass(this) :: create !! Create the distribution functions for all systems, and the learned one. @@ -120,6 +126,8 @@ module raffle__distribs_container !! Set the list of elements for the container. procedure, pass(this), private :: update_element_info !! Update the element information in the container. + procedure, pass(this) :: add_to_element_info + !! Add an element to the container. procedure, pass(this) :: set_element_energy !! Set the energy of an element in the container. procedure, pass(this) :: set_element_energies @@ -156,6 +164,8 @@ module raffle__distribs_container !! Set the generalised distribution function to the default value. procedure, pass(this) :: evolve !! Evolve the learned distribution function. + procedure, pass(this) :: is_converged + !! Check if the learned distribution function has converged. procedure, pass(this) :: write_gdfs @@ -178,12 +188,19 @@ module raffle__distribs_container !! Return the index for element_info given one element. procedure, pass(this) :: get_bin !! Return the bin index for a given distance. + procedure, pass(this) :: get_2body + !! Return the 2-body distribution function. + procedure, pass(this) :: get_3body + !! Return the 3-body distribution function. + procedure, pass(this) :: get_4body + !! Return the 4-body distribution function. end type distribs_container_type interface distribs_container_type !! Interface for the distribution functions container. module function init_distribs_container( & - nbins, width, sigma, cutoff_min, cutoff_max & + nbins, width, sigma, cutoff_min, cutoff_max, & + history_len & ) result(distribs_container) !! Initialise the distribution functions container. integer, dimension(3), intent(in), optional :: nbins @@ -195,6 +212,8 @@ module function init_distribs_container( & real(real32), dimension(3), intent(in), optional :: & cutoff_min, cutoff_max !! Optional. Minimum and maximum cutoff for the 2-, 3-, and 4-body. + integer, intent(in), optional :: history_len + !! Optional. Length of the history for the distribution functions. type(distribs_container_type) :: distribs_container !! Instance of the distribution functions container. end function init_distribs_container @@ -206,7 +225,8 @@ end function init_distribs_container !############################################################################### module function init_distribs_container( & nbins, width, sigma, & - cutoff_min, cutoff_max & + cutoff_min, cutoff_max, & + history_len & ) result(distribs_container) !! Initialise the distribution functions container. implicit none @@ -220,6 +240,8 @@ module function init_distribs_container( & !! 4-body. real(real32), dimension(3), intent(in), optional :: cutoff_min, cutoff_max !! Optional. Minimum and maximum cutoff for the 2-, 3-, and 4-body. + integer, intent(in), optional :: history_len + !! Optional. Length of the history for the distribution functions. type(distribs_container_type) :: distribs_container !! Instance of the distribution functions container. @@ -261,6 +283,14 @@ module function init_distribs_container( & return end if + + if(present(history_len)) distribs_container%history_len = history_len + if(distribs_container%history_len.ge.0) & + allocate( & + distribs_container%history_deltas(history_len), & + source = huge(0._real32) & + ) + end function init_distribs_container !############################################################################### @@ -354,6 +384,40 @@ end subroutine set_radius_distance_tol !############################################################################### +!############################################################################### + subroutine set_history_len(this, history_len) + !! Set the length of the history for the distribution functions. + implicit none + + ! Arguments + class(distribs_container_type), intent(inout) :: this + !! Parent. Instance of distribution functions container. + integer, intent(in) :: history_len + !! Length of the history for the distribution functions. + + ! Local variables + integer :: min_len + !! Minimum length of the history. + real(real32), dimension(:), allocatable :: history_deltas + !! History of the changes in the distribution functions. + + if(history_len.gt.0)then + allocate(history_deltas(history_len), source = huge(0._real32) ) + if(allocated(this%history_deltas))then + min_len = min(this%history_len, history_len) + history_deltas(1:min_len) = this%history_deltas(1:min_len) + deallocate(this%history_deltas) + end if + call move_alloc(history_deltas, this%history_deltas) + else + if(allocated(this%history_deltas)) deallocate(this%history_deltas) + end if + this%history_len = history_len + + end subroutine set_history_len +!############################################################################### + + !############################################################################### subroutine create( & this, basis_list, energy_above_hull_list, deallocate_systems, & @@ -502,6 +566,7 @@ subroutine update( & !! Verbosity level. logical :: suppress_warnings_store !! Boolean to store the suppress_warnings value. + type(distribs_base_type) :: gdf_old ! Set the verbosity level @@ -597,11 +662,24 @@ subroutine update( & end if end if + + ! If history_len is set, temporarily store the old descriptor + if(this%history_len.gt.0)then + gdf_old = this%gdf + this%history_deltas(2:) = this%history_deltas(1:this%history_len-1) + end if + + ! Evolve the distribution functions call this%evolve() if(deallocate_systems_) call this%deallocate_systems() if(this%host_system%defined) & call this%host_system%set_element_map(this%element_info) + ! Evaluate the change in the descriptor from the last iteration + if(this%history_len.gt.0.and.this%num_evaluated.gt.0)then + this%history_deltas(1) = this%gdf%compare(gdf_old) + end if + if(verbose_ .eq. 0)then suppress_warnings = suppress_warnings_store end if @@ -1075,6 +1153,54 @@ end subroutine write_4body !############################################################################### +!############################################################################### + function get_2body(this) result(output) + !! Get the 2-body distribution functions. + implicit none + + ! Arguments + class(distribs_container_type), intent(in) :: this + !! Parent. Instance of distribution functions container. + real(real32), dimension(this%nbins(1),size(this%bond_info)) :: output + + output = this%gdf%df_2body + + end function get_2body +!############################################################################### + + +!############################################################################### + function get_3body(this) result(output) + !! Get the 3-body distribution functions. + implicit none + + ! Arguments + class(distribs_container_type), intent(in) :: this + !! Parent. Instance of distribution functions container. + real(real32), dimension(this%nbins(2),size(this%element_info)) :: output + + output = this%gdf%df_3body + + end function get_3body +!############################################################################### + + +!############################################################################### + function get_4body(this) result(output) + !! Get the 4-body distribution functions. + implicit none + + ! Arguments + class(distribs_container_type), intent(in) :: this + !! Parent. Instance of distribution functions container. + real(real32), dimension(this%nbins(3),size(this%element_info)) :: output + + output = this%gdf%df_4body + + end function get_4body +!############################################################################### + + !############################################################################### subroutine add(this, system) !! Add a system to the container. @@ -1338,6 +1464,29 @@ end subroutine update_element_info !############################################################################### +!############################################################################### + subroutine add_to_element_info(this, element) + !! Add an element to the element_info array. + implicit none + + ! Arguments + class(distribs_container_type), intent(inout) :: this + !! Parent of the procedure. Instance of distribution functions container. + character(len=3), intent(in) :: element + !! Element name. + + + if(.not.allocated(this%element_info)) allocate(this%element_info(0)) + this%element_info = [ & + this%element_info(:), & + element_type(element) & + ] + call this%element_info(size(this%element_info))%set(element) + + end subroutine add_to_element_info +!############################################################################### + + !############################################################################### subroutine set_element_energy(this, element, energy) !! Set the energy of an element in the container. @@ -1358,6 +1507,8 @@ subroutine set_element_energy(this, element, energy) !! Element radius. character(len=3) :: element_ !! Element name without null characters. + logical :: in_element_info + !! Boolean whether the element is in the element_info array. !--------------------------------------------------------------------------- @@ -1369,9 +1520,13 @@ subroutine set_element_energy(this, element, energy) !--------------------------------------------------------------------------- ! if element_info is allocated, update energy of associated index !--------------------------------------------------------------------------- + in_element_info = .false. if(allocated(this%element_info))then idx = findloc( [ this%element_info(:)%name ], element_, dim=1 ) - if(idx.ge.1) this%element_info(idx)%energy = energy + if(idx.ge.1)then + this%element_info(idx)%energy = energy + in_element_info = .true. + end if end if @@ -1390,6 +1545,12 @@ subroutine set_element_energy(this, element, energy) element_database(idx_db)%energy = energy end if + + !--------------------------------------------------------------------------- + ! if element is not in element_info, add it to reserve_element_names + !--------------------------------------------------------------------------- + if(.not.in_element_info) call this%add_to_element_info(element_) + end subroutine set_element_energy !############################################################################### @@ -2319,7 +2480,7 @@ subroutine evolve(this, system) num_evaluated = num_evaluated + 1 if(this%weight_by_hull)then weight = exp( this%system(i)%energy_above_hull / this%kBT ) - if(weight.lt.1.E-6) cycle + if(weight.lt.1.E-6_real32) cycle end if !------------------------------------------------------------------------ ! get the list of 2-body species pairs the system @@ -2383,7 +2544,7 @@ subroutine evolve(this, system) ) & ) / this%kBT & ) - if(weight.lt.1.E-6) cycle + if(weight.lt.1.E-6_real32) cycle end if this%gdf%df_3body(:,idx1) = this%gdf%df_3body(:,idx1) + & @@ -2423,7 +2584,7 @@ subroutine evolve(this, system) ) & ) / this%kBT & ) - if(weight.lt.1.E-6) cycle + if(weight.lt.1.E-6_real32) cycle end if this%gdf%df_2body(:,j) = this%gdf%df_2body(:,j) + & @@ -2442,19 +2603,19 @@ subroutine evolve(this, system) ! if not in the dataset, set distribution functions to default !--------------------------------------------------------------------------- do j = 1, size(this%gdf%df_2body,2) - if(all(abs(this%gdf%df_2body(:,j)).lt.1.E-6))then + if(all(abs(this%gdf%df_2body(:,j)).lt.1.E-6_real32))then call this%set_gdfs_to_default(2, j) else this%in_dataset_2body(j) = .true. end if end do do is = 1, size(this%element_info) - if(all(abs(this%gdf%df_3body(:,is)).lt.1.E-6))then + if(all(abs(this%gdf%df_3body(:,is)).lt.1.E-6_real32))then call this%set_gdfs_to_default(3, is) else this%in_dataset_3body(is) = .true. end if - if(all(abs(this%gdf%df_4body(:,is)).lt.1.E-6))then + if(all(abs(this%gdf%df_4body(:,is)).lt.1.E-6_real32))then call this%set_gdfs_to_default(4, is) else this%in_dataset_4body(is) = .true. @@ -2464,7 +2625,7 @@ subroutine evolve(this, system) allocate(this%norm_2body(size(this%gdf%df_2body,2))) do j = 1, size(this%gdf%df_2body,2) this%norm_2body(j) = maxval(this%gdf%df_2body(:,j)) - if(abs(this%norm_2body(j)).lt.1.E-6)then + if(abs(this%norm_2body(j)).lt.1.E-6_real32)then call stop_program( "Zero norm for 2-body distribution function" ) return end if @@ -2486,12 +2647,12 @@ subroutine evolve(this, system) allocate(this%norm_4body(size(this%element_info))) do is = 1, size(this%element_info) this%norm_3body(is) = maxval(this%gdf%df_3body(:,is)) - if(abs(this%norm_3body(is)).lt.1.E-6)then + if(abs(this%norm_3body(is)).lt.1.E-6_real32)then call stop_program( "Zero norm for 3-body distribution function" ) return end if this%norm_4body(is) = maxval(this%gdf%df_4body(:,is)) - if(abs(this%norm_4body(is)).lt.1.E-6)then + if(abs(this%norm_4body(is)).lt.1.E-6_real32)then call stop_program( "Zero norm for 4-body distribution function" ) return end if @@ -2528,4 +2689,41 @@ subroutine evolve(this, system) end subroutine evolve !############################################################################### + +!############################################################################### + function is_converged(this, threshold) result(converged) + !! Check if the distribution functions have converged. + implicit none + + ! Arguments + class(distribs_container_type), intent(in) :: this + !! Parent of the procedure. Instance of distribution functions container. + real(real32), intent(in), optional :: threshold + !! Threshold for convergence. + logical :: converged + !! Convergence flag. + + ! Local variables + integer :: i, j + !! Loop index. + real(real32) :: threshold_ + !! Threshold for convergence. + + + threshold_ = 1.E-4_real32 + if(present(threshold)) threshold_ = threshold + + if(any(abs(this%history_deltas-huge(0._real32)).lt.1.E-6_real32))then + converged = .false. + return + end if + if(all(abs(this%history_deltas - this%history_deltas(1)).lt.threshold))then + converged = .true. + else + converged = .false. + end if + + end function is_converged +!############################################################################### + end module raffle__distribs_container diff --git a/src/fortran/lib/mod_generator.f90 b/src/fortran/lib/mod_generator.f90 index 075faea..b498b50 100644 --- a/src/fortran/lib/mod_generator.f90 +++ b/src/fortran/lib/mod_generator.f90 @@ -112,12 +112,15 @@ module raffle__generator !! Constructor for the raffle generator type. module function init_raffle_generator( & host, & - width, sigma, cutoff_min, cutoff_max) result(generator) + width, sigma, cutoff_min, cutoff_max, & + history_len & + ) result(generator) type(basis_type), intent(in), optional :: host real(real32), dimension(3), intent(in), optional :: width real(real32), dimension(3), intent(in), optional :: sigma real(real32), dimension(3), intent(in), optional :: cutoff_min real(real32), dimension(3), intent(in), optional :: cutoff_max + integer, intent(in), optional :: history_len type(raffle_generator_type) :: generator end function init_raffle_generator end interface raffle_generator_type @@ -127,7 +130,8 @@ end function init_raffle_generator !############################################################################### module function init_raffle_generator( & - host, width, sigma, cutoff_min, cutoff_max & + host, width, sigma, cutoff_min, cutoff_max, & + history_len & ) result(generator) !! Initialise an instance of the raffle generator. !! @@ -147,6 +151,8 @@ module function init_raffle_generator( & !! Minimum cutoff for the 2-, 3-, and 4-body distribution functions. real(real32), dimension(3), intent(in), optional :: cutoff_max !! Maximum cutoff for the 2-, 3-, and 4-body distribution functions. + integer, intent(in), optional :: history_len + !! Length of the history for the 2-, 3-, and 4-body distribution functions. ! Local variables type(raffle_generator_type) :: generator @@ -157,15 +163,18 @@ module function init_raffle_generator( & if(present(host)) call generator%set_host(host) ! Set up the distribution function parameters - if( present(width) ) & + if(present(width)) & call generator%distributions%set_width(width) - if( present(sigma) ) & + if(present(sigma)) & call generator%distributions%set_sigma(sigma) - if( present(cutoff_min) ) & + if(present(cutoff_min)) & call generator%distributions%set_cutoff_min(cutoff_min) - if( present(cutoff_max) ) & + if(present(cutoff_max)) & call generator%distributions%set_cutoff_max(cutoff_max) + if(present(history_len)) & + call generator%distributions%set_history_len(history_len) + end function init_raffle_generator !############################################################################### diff --git a/src/raffle/raffle.py b/src/raffle/raffle.py index 19340c8..3de40db 100644 --- a/src/raffle/raffle.py +++ b/src/raffle/raffle.py @@ -713,6 +713,17 @@ def set_radius_distance_tol(self, radius_distance_tol): _raffle.f90wrap_raffle__dc__set_radius_distance_tol__binding__dc_type(this=self._handle, \ radius_distance_tol=radius_distance_tol) + def set_history_len(self, history_len : int = None): + """ + Set the history length for the convergence check of the RAFFLE descriptor. + + Parameters: + history_len (int): + Length of the history for checking convergence of the descriptor. + """ + _raffle.f90wrap_raffle__dc__set_history_len__binding__dc_type(this=self._handle, \ + history_len=history_len) + def create(self, basis_list, energy_above_hull_list=None, deallocate_systems=True, verbose=None): """ Create the distribution functions. @@ -934,6 +945,21 @@ def evolve(self): #, system=None): # _raffle.f90wrap_raffle__dc__evolve__binding__dc_type(this=self._handle, \ # system=None if system is None else system._handle) + def is_converged(self, threshold : float = 1e-4): + """ + Parameters: + threshold (float): + Threshold for convergence. + + Returns: + converged (bool): + Boolean indicating whether the distribution functions have converged. + """ + converged = \ + _raffle.f90wrap_raffle__dc__is_converged__binding__dc_type(this=self._handle, \ + threshold=threshold) + return converged + def write_gdfs(self, file): """ Write the generalised distribution functions (GDFs) to a file. @@ -1047,8 +1073,7 @@ def _get_bin(self, value, dim): Dimension of the distribution function. 1 for 2-body, 2 for 3-body, and 3 for 4-body. - Returns - ------- + Returns: bin (int): Bin index for the value in the distribution functions. @@ -1058,6 +1083,50 @@ def _get_bin(self, value, dim): value=value, dim=dim) return bin + def get_2body(self): + """ + Get the 2-body distribution function. + + Returns: + output (float array): 2D array of shape (num_species_pairs, nbins) + 2-body distribution function. + + """ + n0 = self.nbins[0] + n1 = _raffle.f90wrap_raffle__dc__get_num_pairs__dc_type(this=self._handle) + output = \ + _raffle.f90wrap_raffle__dc__get_2body__binding__dc_type(this=self._handle, n0=n0, n1=n1).reshape((n0, n1), order='F').T + return output + + def get_3body(self): + """ + Get the 3-body distribution function. + + Returns: + output (float array): 2D array of shape (num_species, nbins) + 3-body distribution function. + + """ + n0 = self.nbins[1] + n1 = _raffle.f90wrap_raffle__dc__get_num_species__dc_type(this=self._handle) + output = \ + _raffle.f90wrap_raffle__dc__get_3body__binding__dc_type(this=self._handle, n0=n0, n1=n1).reshape((n0, n1), order='F').T + return output + + def get_4body(self): + """ + Get the 4-body distribution function. + + Returns: + output (float array): 2D array of shape (num_species, nbins) + 4-body distribution function. + """ + n0 = self.nbins[2] + n1 = _raffle.f90wrap_raffle__dc__get_num_species__dc_type(this=self._handle) + output = \ + _raffle.f90wrap_raffle__dc__get_4body__binding__dc_type(this=self._handle, n0=n0, n1=n1).reshape((n0, n1), order='F').T + return output + @property def num_evaluated(self): """ @@ -1094,6 +1163,38 @@ def kBT(self): def kBT(self, kBT): _raffle.f90wrap_distribs_container_type__set__kbt(self._handle, kBT) + @property + def history_len(self): + """ + Length of the history for convergence checking. + """ + return _raffle.f90wrap_distribs_container_type__get__history_len(self._handle) + + @history_len.setter + def history_len(self, history_len): + _raffle.f90wrap_distribs_container_type__set__history_len(self._handle, \ + history_len) + + @property + def history_deltas(self): + """ + History of the descriptor convergence. + """ + array_ndim, array_type, array_shape, array_handle = \ + _raffle.f90wrap_distribs_container_type__array__history_deltas(self._handle) + if array_handle in self._arrays: + history_deltas = self._arrays[array_handle] + else: + history_deltas = f90wrap.runtime.get_array(f90wrap.runtime.sizeof_fortran_t, + self._handle, + _raffle.f90wrap_distribs_container_type__array__history_deltas) + self._arrays[array_handle] = history_deltas + return history_deltas + + @history_deltas.setter + def history_deltas(self, history_deltas): + self.history_deltas[...] = history_deltas + @property def weight_by_hull(self): """ @@ -1245,6 +1346,10 @@ def __str__(self): ret.append(repr(self.num_evaluated_allocated)) ret.append(',\n kBT : ') ret.append(repr(self.kBT)) + ret.append(',\n history_len : ') + ret.append(repr(self.history_len)) + ret.append(',\n history_deltas : ') + ret.append(repr(self.history_deltas)) ret.append(',\n weight_by_hull : ') ret.append(repr(self.weight_by_hull)) ret.append(',\n nbins : ') @@ -1428,17 +1533,23 @@ def deallocate(self): @f90wrap.runtime.register_class("raffle.raffle_generator") class raffle_generator(f90wrap.runtime.FortranDerivedType): - def __init__(self, handle=None): + def __init__(self, handle=None, history_len=None): """ Create a ``raffle_generator`` object. This object is used to generate structures using the RAFFLE method. The object has procedures to set the parameters for the generation and to generate the structures. + + Parameters: + history_len (int): + Length of the history for checking convergence of descriptor. """ f90wrap.runtime.FortranDerivedType.__init__(self) result = _raffle.f90wrap_generator__raffle_generator_type_initialise() self._handle = result[0] if isinstance(result, tuple) else result + if history_len is not None: + self.distributions.set_history_len(history_len) def __del__(self): """ @@ -1683,7 +1794,7 @@ def generate( # exit if any other method is provided for key in method_ratio.keys(): - if key not in ["voide", "rand", "walk", "grow", "min"]: + if key not in ["void", "rand", "walk", "grow", "min"]: raise ValueError(f"Unknown method {key} provided, use only void, rand (random), walk, grow (growth), or min (minimum/global)") method_ratio_list = [] @@ -1807,6 +1918,21 @@ def read_settings(self, file): _raffle.f90wrap_generator__read_settings__binding__rgt(this=self._handle, \ file=file) + def get_descriptor(self): + """ + Get the descriptor for the generated structures. + + Returns: + descriptor (list of numpy arrays): + The descriptor for the generated structures. + """ + descriptor = [ + self.distributions.get_2body(), + self.distributions.get_3body(), + self.distributions.get_4body() + ] + return descriptor + @property def num_structures(self): """ diff --git a/src/wrapper/f90wrap_mod_distribs_container.f90 b/src/wrapper/f90wrap_mod_distribs_container.f90 index 8fa370d..2b74713 100644 --- a/src/wrapper/f90wrap_mod_distribs_container.f90 +++ b/src/wrapper/f90wrap_mod_distribs_container.f90 @@ -14,7 +14,7 @@ subroutine f90wrap_distribs_container_type__get__num_evaluated( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr integer, intent(out) :: f90wrap_num_evaluated - + this_ptr = transfer(this, this_ptr) f90wrap_num_evaluated = this_ptr%p%num_evaluated end subroutine f90wrap_distribs_container_type__get__num_evaluated @@ -30,7 +30,7 @@ subroutine f90wrap_distribs_container_type__set__num_evaluated( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr integer, intent(in) :: f90wrap_num_evaluated - + this_ptr = transfer(this, this_ptr) this_ptr%p%num_evaluated = f90wrap_num_evaluated end subroutine f90wrap_distribs_container_type__set__num_evaluated @@ -46,7 +46,7 @@ subroutine f90wrap_distribs_container_type__get__num_evaluated_allocated( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr integer, intent(out) :: f90wrap_num_evaluated_allocated - + this_ptr = transfer(this, this_ptr) f90wrap_num_evaluated_allocated = this_ptr%p%num_evaluated_allocated end subroutine f90wrap_distribs_container_type__get__num_evaluated_allocated @@ -62,7 +62,7 @@ subroutine f90wrap_distribs_container_type__set__num_evaluated_allocated( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr integer, intent(in) :: f90wrap_num_evaluated_allocated - + this_ptr = transfer(this, this_ptr) this_ptr%p%num_evaluated_allocated = f90wrap_num_evaluated_allocated end subroutine f90wrap_distribs_container_type__set__num_evaluated_allocated @@ -70,7 +70,7 @@ end subroutine f90wrap_distribs_container_type__set__num_evaluated_allocated !############################################################################### -! set energy scaling +! energy scaling !############################################################################### subroutine f90wrap_distribs_container_type__get__kBT(this, f90wrap_kBT) use raffle__distribs_container, only: distribs_container_type @@ -81,7 +81,7 @@ subroutine f90wrap_distribs_container_type__get__kBT(this, f90wrap_kBT) integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr real(4), intent(out) :: f90wrap_kBT - + this_ptr = transfer(this, this_ptr) f90wrap_kBT = this_ptr%p%kBT end subroutine f90wrap_distribs_container_type__get__kBT @@ -95,13 +95,77 @@ subroutine f90wrap_distribs_container_type__set__kBT(this, f90wrap_kBT) integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr real(4), intent(in) :: f90wrap_kBT - + this_ptr = transfer(this, this_ptr) this_ptr%p%kBT = f90wrap_kBT end subroutine f90wrap_distribs_container_type__set__kBT !############################################################################### +!############################################################################### +! convergence history +!############################################################################### +subroutine f90wrap_distribs_container_type__get__history_len( & + this, f90wrap_history_len & +) + use raffle__distribs_container, only: distribs_container_type + implicit none + type distribs_container_type_ptr_type + type(distribs_container_type), pointer :: p => NULL() + end type distribs_container_type_ptr_type + integer, intent(in) :: this(2) + type(distribs_container_type_ptr_type) :: this_ptr + integer, intent(out) :: f90wrap_history_len + + this_ptr = transfer(this, this_ptr) + f90wrap_history_len = this_ptr%p%history_len +end subroutine f90wrap_distribs_container_type__get__history_len + +subroutine f90wrap_distribs_container_type__set__history_len( & + this, f90wrap_history_len & +) + use raffle__distribs_container, only: distribs_container_type + implicit none + type distribs_container_type_ptr_type + type(distribs_container_type), pointer :: p => NULL() + end type distribs_container_type_ptr_type + integer, intent(in) :: this(2) + type(distribs_container_type_ptr_type) :: this_ptr + integer, intent(in) :: f90wrap_history_len + + this_ptr = transfer(this, this_ptr) + this_ptr%p%history_len = f90wrap_history_len +end subroutine f90wrap_distribs_container_type__set__history_len + +subroutine f90wrap_distribs_container_type__array__history_deltas( & + this, nd, dtype, dshape, dloc & +) + use raffle__distribs_container, only: distribs_container_type + use, intrinsic :: iso_c_binding, only : c_int + implicit none + type distribs_container_type_ptr_type + type(distribs_container_type), pointer :: p => NULL() + end type distribs_container_type_ptr_type + integer(c_int), intent(in) :: this(2) + type(distribs_container_type_ptr_type) :: this_ptr + integer(c_int), intent(out) :: nd + integer(c_int), intent(out) :: dtype + integer(c_int), dimension(10), intent(out) :: dshape + integer*8, intent(out) :: dloc + + nd = 1 + dtype = 11 + this_ptr = transfer(this, this_ptr) + if (allocated(this_ptr%p%history_deltas)) then + dshape(1:1) = shape(this_ptr%p%history_deltas) + dloc = loc(this_ptr%p%history_deltas) + else + dloc = 0 + end if +end subroutine f90wrap_distribs_container_type__array__history_deltas +!############################################################################### + + !############################################################################### ! boolean for using hull weighting or empirical method !############################################################################### @@ -116,7 +180,7 @@ subroutine f90wrap_distribs_container_type__get__weight_by_hull( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr logical, intent(out) :: f90wrap_weight_by_hull - + this_ptr = transfer(this, this_ptr) f90wrap_weight_by_hull = this_ptr%p%weight_by_hull end subroutine f90wrap_distribs_container_type__get__weight_by_hull @@ -132,7 +196,7 @@ subroutine f90wrap_distribs_container_type__set__weight_by_hull( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr logical, intent(in) :: f90wrap_weight_by_hull - + this_ptr = transfer(this, this_ptr) this_ptr%p%weight_by_hull = f90wrap_weight_by_hull end subroutine f90wrap_distribs_container_type__set__weight_by_hull @@ -153,7 +217,7 @@ subroutine f90wrap_distribs_container_type__get__viability_3body_default( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr real(4), intent(out) :: f90wrap_viability_3body_default - + this_ptr = transfer(this, this_ptr) f90wrap_viability_3body_default = this_ptr%p%viability_3body_default end subroutine f90wrap_distribs_container_type__get__viability_3body_default @@ -169,7 +233,7 @@ subroutine f90wrap_distribs_container_type__set__viability_3body_default( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr real(4), intent(in) :: f90wrap_viability_3body_default - + this_ptr = transfer(this, this_ptr) this_ptr%p%viability_3body_default = f90wrap_viability_3body_default end subroutine f90wrap_distribs_container_type__set__viability_3body_default @@ -185,7 +249,7 @@ subroutine f90wrap_distribs_container_type__get__viability_4body_default( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr real(4), intent(out) :: f90wrap_viability_4body_default - + this_ptr = transfer(this, this_ptr) f90wrap_viability_4body_default = this_ptr%p%viability_4body_default end subroutine f90wrap_distribs_container_type__get__viability_4body_default @@ -201,7 +265,7 @@ subroutine f90wrap_distribs_container_type__set__viability_4body_default( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr real(4), intent(in) :: f90wrap_viability_4body_default - + this_ptr = transfer(this, this_ptr) this_ptr%p%viability_4body_default = f90wrap_viability_4body_default end subroutine f90wrap_distribs_container_type__set__viability_4body_default @@ -226,7 +290,7 @@ subroutine f90wrap_distribs_container_type__array__nbins( & integer(c_int), intent(out) :: dtype integer(c_int), dimension(10), intent(out) :: dshape integer*8, intent(out) :: dloc - + nd = 1 dtype = 5 this_ptr = transfer(this, this_ptr) @@ -249,7 +313,7 @@ subroutine f90wrap_distribs_container_type__array__sigma( & integer(c_int), intent(out) :: dtype integer(c_int), dimension(10), intent(out) :: dshape integer*8, intent(out) :: dloc - + nd = 1 dtype = 11 this_ptr = transfer(this, this_ptr) @@ -272,7 +336,7 @@ subroutine f90wrap_distribs_container_type__array__width( & integer(c_int), intent(out) :: dtype integer(c_int), dimension(10), intent(out) :: dshape integer*8, intent(out) :: dloc - + nd = 1 dtype = 11 this_ptr = transfer(this, this_ptr) @@ -295,7 +359,7 @@ subroutine f90wrap_distribs_container_type__array__cutoff_min( & integer(c_int), intent(out) :: dtype integer(c_int), dimension(10), intent(out) :: dshape integer*8, intent(out) :: dloc - + nd = 1 dtype = 11 this_ptr = transfer(this, this_ptr) @@ -318,7 +382,7 @@ subroutine f90wrap_distribs_container_type__array__cutoff_max( & integer(c_int), intent(out) :: dtype integer(c_int), dimension(10), intent(out) :: dshape integer*8, intent(out) :: dloc - + nd = 1 dtype = 11 this_ptr = transfer(this, this_ptr) @@ -341,7 +405,7 @@ subroutine f90wrap_distribs_container_type__array__radius_distance_tol( & integer(c_int), intent(out) :: dtype integer(c_int), dimension(10), intent(out) :: dshape integer*8, intent(out) :: dloc - + nd = 1 dtype = 11 this_ptr = transfer(this, this_ptr) @@ -357,7 +421,7 @@ end subroutine f90wrap_distribs_container_type__array__radius_distance_tol subroutine f90wrap_raffle__dc__dc_type_initialise(this) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -370,7 +434,7 @@ end subroutine f90wrap_raffle__dc__dc_type_initialise subroutine f90wrap_raffle__dc__dc_type_finalise(this) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -390,7 +454,7 @@ subroutine f90wrap_raffle__dc__set_width__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -406,7 +470,7 @@ subroutine f90wrap_raffle__dc__set_sigma__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -422,7 +486,7 @@ subroutine f90wrap_raffle__dc__set_cutoff_min__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -438,7 +502,7 @@ subroutine f90wrap_raffle__dc__set_cutoff_max__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -454,7 +518,7 @@ subroutine f90wrap_raffle__dc__set_radius_distance_tol__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -466,6 +530,22 @@ subroutine f90wrap_raffle__dc__set_radius_distance_tol__binding__dc_type( & radius_distance_tol=radius_distance_tol & ) end subroutine f90wrap_raffle__dc__set_radius_distance_tol__binding__dc_type + +subroutine f90wrap_raffle__dc__set_history_len__binding__dc_type( & + this, history_len & +) + use raffle__distribs_container, only: distribs_container_type + implicit none + + type distribs_container_type_ptr_type + type(distribs_container_type), pointer :: p => NULL() + end type distribs_container_type_ptr_type + type(distribs_container_type_ptr_type) :: this_ptr + integer, intent(in), dimension(2) :: this + integer, intent(in) :: history_len + this_ptr = transfer(this, this_ptr) + call this_ptr%p%set_history_len(history_len=history_len) +end subroutine f90wrap_raffle__dc__set_history_len__binding__dc_type !############################################################################### @@ -478,7 +558,7 @@ subroutine f90wrap_raffle__dc__create__binding__dc_type( & use raffle__geom_rw, only: basis_type use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -517,7 +597,7 @@ subroutine f90wrap_raffle__dc__update__binding__dc_type( & use raffle__geom_rw, only: basis_type use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -560,7 +640,7 @@ subroutine f90wrap_raffle__dc__deallocate_systems__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -582,7 +662,7 @@ subroutine f90wrap_raffle__dc__add_basis__binding__dc_type( & use raffle__geom_rw, only: basis_type use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -614,7 +694,7 @@ subroutine f90wrap_raffle__dc__get__num_elements( & integer, intent(in) :: this(2) type(distribs_container_type_ptr_type) :: this_ptr integer, intent(out) :: ret_num_elements - + this_ptr = transfer(this, this_ptr) if(.not.allocated(this_ptr%p%element_info)) then ret_num_elements = 0 @@ -632,7 +712,7 @@ subroutine f90wrap_raffle__dc__set_element_energy__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -649,7 +729,7 @@ subroutine f90wrap_raffle__dc__set_element_energies__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -668,7 +748,7 @@ subroutine f90wrap_raffle__dc__get_element_energies_sm__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -688,7 +768,7 @@ subroutine f90wrap_raffle__dc__set_bond_radius__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -705,7 +785,7 @@ subroutine f90wrap_raffle__dc__set_bond_radii__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -728,7 +808,7 @@ subroutine f90wrap_raffle__dc__get_bond_radii_staticmem__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -750,7 +830,7 @@ end subroutine f90wrap_raffle__dc__get_bond_radii_staticmem__binding__dc_type subroutine f90wrap_raffle__dc__initialise_gdfs__binding__dc_type(this) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -768,7 +848,7 @@ end subroutine f90wrap_raffle__dc__initialise_gdfs__binding__dc_type subroutine f90wrap_raffle__dc__evolve__binding__dc_type(this) !, system) use raffle__distribs_container, only: distribs_container_type implicit none - + ! type distribs_type_ptr_type ! type(distribs_type), pointer :: p => NULL() ! end type distribs_type_ptr_type @@ -790,6 +870,28 @@ end subroutine f90wrap_raffle__dc__evolve__binding__dc_type !############################################################################### +!############################################################################### +! check for convergence of the descriptor +!############################################################################### +subroutine f90wrap_raffle__dc__is_converged__binding__dc_type( & + this, ret_converged, threshold & +) + use raffle__distribs_container, only: distribs_container_type + implicit none + + type distribs_container_type_ptr_type + type(distribs_container_type), pointer :: p => NULL() + end type distribs_container_type_ptr_type + type(distribs_container_type_ptr_type) :: this_ptr + integer, intent(in), dimension(2) :: this + logical, intent(out) :: ret_converged + real(4), intent(in), optional :: threshold + this_ptr = transfer(this, this_ptr) + ret_converged = this_ptr%p%is_converged(threshold=threshold) +end subroutine f90wrap_raffle__dc__is_converged__binding__dc_type +!############################################################################### + + !############################################################################### ! read and write distribution functions to file !############################################################################### @@ -798,7 +900,7 @@ subroutine f90wrap_raffle__dc__read_gdfs__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -814,7 +916,7 @@ subroutine f90wrap_raffle__dc__write_gdfs__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -830,7 +932,7 @@ subroutine f90wrap_raffle__dc__read_dfs__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -846,7 +948,7 @@ subroutine f90wrap_raffle__dc__write_dfs__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -862,7 +964,7 @@ subroutine f90wrap_raffle__dc__write_2body__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -878,7 +980,7 @@ subroutine f90wrap_raffle__dc__write_3body__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -894,7 +996,7 @@ subroutine f90wrap_raffle__dc__write_4body__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -915,7 +1017,7 @@ subroutine f90wrap_raffle__dc__get_pair_index__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -933,7 +1035,7 @@ subroutine f90wrap_raffle__dc__get_bin__binding__dc_type( & ) use raffle__distribs_container, only: distribs_container_type implicit none - + type distribs_container_type_ptr_type type(distribs_container_type), pointer :: p => NULL() end type distribs_container_type_ptr_type @@ -947,5 +1049,95 @@ subroutine f90wrap_raffle__dc__get_bin__binding__dc_type( & end subroutine f90wrap_raffle__dc__get_bin__binding__dc_type !############################################################################### -! End of module raffle__distribs_container defined in file ../src/lib/mod_distribs_container.f90 +!############################################################################### +! get the generalised descriptor n-body components +!############################################################################### +subroutine f90wrap_raffle__dc__get_2body__binding__dc_type( & + this, ret_output, n0, n1 & +) + use raffle__distribs_container, only: distribs_container_type + implicit none + + type distribs_container_type_ptr_type + type(distribs_container_type), pointer :: p => NULL() + end type distribs_container_type_ptr_type + real(4), intent(out), dimension(n0,n1) :: ret_output + type(distribs_container_type_ptr_type) :: this_ptr + integer, intent(in), dimension(2) :: this + integer :: n0 + integer :: n1 + this_ptr = transfer(this, this_ptr) + ret_output = this_ptr%p%get_2body() +end subroutine f90wrap_raffle__dc__get_2body__binding__dc_type + +subroutine f90wrap_raffle__dc__get_3body__binding__dc_type( & + this, ret_output, n0, n1 & +) + use raffle__distribs_container, only: distribs_container_type + implicit none + + type distribs_container_type_ptr_type + type(distribs_container_type), pointer :: p => NULL() + end type distribs_container_type_ptr_type + real(4), intent(out), dimension(n0,n1) :: ret_output + type(distribs_container_type_ptr_type) :: this_ptr + integer, intent(in), dimension(2) :: this + integer :: n0 + integer :: n1 + this_ptr = transfer(this, this_ptr) + ret_output = this_ptr%p%get_3body() +end subroutine f90wrap_raffle__dc__get_3body__binding__dc_type + +subroutine f90wrap_raffle__dc__get_4body__binding__dc_type( & + this, ret_output, n0, n1 & +) + use raffle__distribs_container, only: distribs_container_type + implicit none + + type distribs_container_type_ptr_type + type(distribs_container_type), pointer :: p => NULL() + end type distribs_container_type_ptr_type + real(4), intent(out), dimension(n0,n1) :: ret_output + type(distribs_container_type_ptr_type) :: this_ptr + integer, intent(in), dimension(2) :: this + integer :: n0 + integer :: n1 + this_ptr = transfer(this, this_ptr) + ret_output = this_ptr%p%get_4body() +end subroutine f90wrap_raffle__dc__get_4body__binding__dc_type + +subroutine f90wrap_raffle__dc__get_num_species__dc_type( & + this, ret_num_species & +) + use raffle__distribs_container, only: distribs_container_type + implicit none + + type distribs_container_type_ptr_type + type(distribs_container_type), pointer :: p => NULL() + end type distribs_container_type_ptr_type + type(distribs_container_type_ptr_type) :: this_ptr + integer, intent(in), dimension(2) :: this + integer, intent(out) :: ret_num_species + this_ptr = transfer(this, this_ptr) + ret_num_species = size(this_ptr%p%element_info, dim = 1) +end subroutine f90wrap_raffle__dc__get_num_species__dc_type + +subroutine f90wrap_raffle__dc__get_num_pairs__dc_type( & + this, ret_num_pairs & +) + use raffle__distribs_container, only: distribs_container_type + implicit none + + type distribs_container_type_ptr_type + type(distribs_container_type), pointer :: p => NULL() + end type distribs_container_type_ptr_type + type(distribs_container_type_ptr_type) :: this_ptr + integer, intent(in), dimension(2) :: this + integer, intent(out) :: ret_num_pairs + this_ptr = transfer(this, this_ptr) + ret_num_pairs = size(this_ptr%p%bond_info, dim = 1) +end subroutine f90wrap_raffle__dc__get_num_pairs__dc_type +!############################################################################### + +! End of module raffle__distribs_container defined in file ../src/lib/mod_distribs_container.f90