diff --git a/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md b/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md new file mode 100644 index 0000000..06e41b2 --- /dev/null +++ b/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md @@ -0,0 +1,14 @@ +# Description + +> _Briefly describe the purpose and content of this merge request._ +> _E.g. "Add structure generation constraints", or "Fix bug in energy parsing"._ + +# Checklist + +Mark with `x` when complete, or `~` if not applicable. + +- [ ] [I have read and followed the **RAFFLE's contribution guidelines.**](https://github.com/ExeQuantCode/RAFFLE/blob/main/CONTRIBUTING.md) +- [ ] **Code is commented** appropriately, and API docstrings follow NumPy or FORD style. +- [ ] **Read*the*Docs** documentation is added/updated (if applicable). +- [ ] **Unit tests** are added/updated (if applicable). +- [ ] **Linked issue** is resolved with a `closes #XXXX` reference (if applicable). diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 4c6267a..0922aeb 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -5,7 +5,7 @@ repos: - id: trailing-whitespace - id: end-of-file-fixer - repo: https://github.com/nedtaylor/fortran-format-hooks - rev: 1aad19af1f0c87027829f0f872d4bdc6ed9251e2 + rev: 749c61fce3a5f29f6d5d3b236ed9b693c6b2bf4e hooks: - id: check-fortran-indentation args: [--line-length=80, --ignore-directories=src/wrapper] diff --git a/CMakeLists.txt b/CMakeLists.txt index 107a7e1..b75dc66 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -74,6 +74,7 @@ set(LIB_DIR ${FORTRAN_SRC_DIR}/lib) set(LIB_FILES mod_io_utils.F90 mod_constants.f90 + mod_cache.f90 mod_misc.f90 mod_tools_infile.f90 mod_misc_maths.f90 diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 62238b6..eb6cc68 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -92,7 +92,7 @@ When submitting your contributions, please ensure the following: - Reference any related issues or pull requests, if applicable. - Write unit tests for your contributions - Ensure all existing tests pass before submitting your changes. -- Update the documentation to reflect your changes, if necessary (i.e. through FORD style commenting). +- Update the documentation to reflect your changes, if necessary (i.e. through FORD style commenting for Fortran and NumPy docstrings style for Python). - Provide examples and usage instructions, if applicable. Follow the [Code Style](#code-style) when contributing code to this project to ensure compatibility and a uniform format to the project. @@ -101,7 +101,7 @@ Follow the [Code Style](#code-style) when contributing code to this project to e ### Code Style - Follow the existing code style and conventions. - Use meaningful variable and function names. -- Write clear and concise comments. For the Fortran library, use comments compatible with the [FORD Fortran Documenter](https://forddocs.readthedocs.io/en/stable/). For the Python wrapper, use comments compatible with [pandoc](https://pandoc.org). +- Write clear and concise comments. For the Fortran library, use comments compatible with the [FORD Fortran Documenter](https://forddocs.readthedocs.io/en/stable/). For the Python wrapper, use comments compatible with [pandoc](https://pandoc.org), following the [NumPy style guide](https://numpydoc.readthedocs.io/en/latest/format.html). diff --git a/docs/source/tutorials/graphene_grain_boundary_tutorial.rst b/docs/source/tutorials/graphene_grain_boundary_tutorial.rst index a8c3f9f..9b8f273 100644 --- a/docs/source/tutorials/graphene_grain_boundary_tutorial.rst +++ b/docs/source/tutorials/graphene_grain_boundary_tutorial.rst @@ -30,6 +30,9 @@ First, we must import the required packages: from raffle.generator import raffle_generator from mace.calculators import mace_mp import numpy as np + from pathlib import Path + + script_dir = Path(__file__).resolve().parent Next, we need to set up the RAFFLE generator and the calculator to calculate the energies of the structures. @@ -55,14 +58,14 @@ The host is read in from a file, but it can also be generated using the ARTEMIS .. code-block:: python - host = read("../POSCAR_host_gb") + host = read(script_dir / ".." / "POSCAR_host_gb") generator.set_host(host) We then need to set up the RAFFLE generator by creating the descriptor. .. code-block:: python - graphene = read("../POSCAR_graphene") + graphene = read(script_dir / ".." / "POSCAR_graphene") h2 = build.molecule("H2") graphene.calc = calc C_reference_energy = graphene.get_potential_energy() / len(graphene) diff --git a/docs/source/tutorials/index.rst b/docs/source/tutorials/index.rst index a4e562e..ee4ed6e 100644 --- a/docs/source/tutorials/index.rst +++ b/docs/source/tutorials/index.rst @@ -35,3 +35,9 @@ Whilst RAFFLE is a random sturcture search package designed primarlily for inter aluminium_tutorial Si-Ge_tutorial graphene_grain_boundary_tutorial + +.. toctree:: + :maxdepth: 2 + :caption: Visualisation: + + visualisation diff --git a/docs/source/tutorials/parameters_tutorial.rst b/docs/source/tutorials/parameters_tutorial.rst index 9cd0299..2537ac0 100644 --- a/docs/source/tutorials/parameters_tutorial.rst +++ b/docs/source/tutorials/parameters_tutorial.rst @@ -14,6 +14,7 @@ For a guide on how to build a database, see the :doc:`Databases tutorial ` + +We can now use this to compare the initial descriptor with the updated descriptor after generating new structures. + +.. code-block:: python + + # Generate new structures and update the descriptor + structures = [ + # List of structures to be generated + ] + generator.distributions.update(structures) + + # Retrieve the updated descriptor + descriptor_new = generator.get_descriptor() + + # Print the updated descriptor on the plots and compare + ... + + +Visualising a RAFFLE fingerprint +-------------------------------- + +RAFFLE fingerprints are the distribution functions for each structure in the learning database. +These are then weighted by energy (formation or convex hull) to form the RAFFLE descriptor. + +However, the individual fingerprints can also be extracted and visualised. + +The `raffle_generator` object has a method `get_fingerprint()` that returns the distribution functions for a provided structure. +The output is a list of three numpy arrays, with the first array containing the 2-body fingerprint, the second array containing the 3-body fingerprint, and the third array containing the 4-body fingerprint. +Each `n`-body fingerprint is a 2D array, with the first column containing the species index (or element pair index for the 2-body fingerprint) and the second column containing the binned fingerprint value. +Like above, the bin lengths are set either explicitly or determined by the `cutoff_min`, `cutoff_max`, and `width` components of the generator. +Here is an example of how to use the `get_fingerprint()` method: + +.. code-block:: python + + # Initialise RAFFLE generator + from raffle.generator import raffle_generator + + generator = raffle_generator() + + # Optional parameters + generator.distributions.set_width([0.04, np.pi/160.0, np.pi/160.0]) + generator.distributions.set_cutoff_min([0.5, 0.0, 0.0]) + generator.distributions.set_cutoff_max([6.0, np.pi, np.pi]) + + # Structure to obtain the fingerprint for + structure = Atoms( + # Structure to be used for the fingerprint + ) + + fingerprint = generator.distributions.generate_fingerprint(structure) + +This can then be visualised in a similar way to the descriptor. + +.. code-block:: python + + # Create a figure with 3 subplots side by side + fig, axes = plt.subplots(1, 3, figsize=(18, 5)) + + # Plot for each n-body function (2-body, 3-body, 4-body) + for j in range(3): + # Calculate x-axis values + x = np.arange(generator.distributions.cutoff_min[j], + generator.distributions.cutoff_max[j] + generator.distributions.width[j], + generator.distributions.width[j]) + + # Plot on the respective subplot + for idx in range(len(fingerprint[j])): + axes[j].plot(x, fingerprint[j][idx,:]) + + # Set labels and title for each subplot + axes[j].set_ylabel('Fingerprint value') + axes[j].set_title(f'{j+2}-body fingerprint') + + axes[0].set_xlabel('Distance (Å)') + axes[1].set_xlabel('3-body angle (radians)') + axes[2].set_xlabel('Improper dihedral angle (radians)') + plt.tight_layout() + plt.show() + +An example python notebook is provided in :git:`examples/python_pkg/visualisation/fingerprint.ipynb `. + + +Visualising RAFFLE probability density +-------------------------------------- + +RAFFLE probability density is the probability of finding a given element in a given position in the system. +This is calculated by the RAFFLE generator and can be visualised using the `get_probability_density()` method. +The output is a 2D array, with the first column containing the coordinates (spatial and species) the second column containing the binned probability density value. + +A structure is provided to the `get_probability_density()`, along with a list of elements to calculate the probability density for. +These elements must be present in the RAFFLE descriptor. + +An example of how to use the `get_probability_density()` method is shown below: + +.. code-block:: python + + # Initialise RAFFLE generator + from raffle.generator import raffle_generator + + generator = raffle_generator() + + generator.distributions.set_element_energies( + { + # reference energies for all elements in the systems + } + ) + + database = [ + # List of structures in the learning database + ] + generator.distributions.create(database) + + # Structure to obtain the probability density for + structure = Atoms( + # Structure to be used for the probability density + ) + + species = 'SiGe' + + probability_density, grid = generator.get_probability_density(structure, species, return_grid=True) + +The first index of the first column of `probability_density` is the x-coordinate, the second index is the y-coordinate, and the third index is the z-coordinate. +The fourth index is the distance between the position and the nearest atom (i.e. the void value). +The fifth index onwards is the species index (in order of the species list provided). +The second column is the site index. + +For a more extensive example, see the `examples/python_pkg/visualisation/probability_density.ipynb` notebook. +This also provides a visualisation of the probability density using the `matplotlib` library. diff --git a/example/python_pkg/README.md b/example/python_pkg/README.md index c9b2b29..4adc65d 100644 --- a/example/python_pkg/README.md +++ b/example/python_pkg/README.md @@ -28,7 +28,7 @@ python_pkg/ To run a RAFFLE example and analyse the results, the following order must be performed. 1. Move to the desired `SYSTEM_learn/DRAFFLE` directory and create a directory to work in: ```bash - cd SYSTEM_learn + cd SYSTEM_learn/DRAFFLE mkdir DOutput cd DOutput ``` @@ -41,7 +41,7 @@ To run a RAFFLE example and analyse the results, the following order must be per 3. Go to parent directory and open the notebook `pca.ipynb`, run all the cells. NOTE: for the `C_learn/DRSS/` example, the `pca.ipynb` notebook expects that the `C_learn/DRAFFLE/` example script and notebook have been fully run first. -This is because it attempts to use the PCA fit to the RAFFLE results to transform its data. +This is because it attempts to use the PCA of the RAFFLE results to transform its data. Doing so enables the two techniques to be properly compared. The `[un]rlxd_structures_seed0.traj` files are provided (as are the `energies_[un]rlxd+seed0.txt`) as the key outputs from the `learn.py` runs. @@ -53,4 +53,4 @@ Whilst other output files are generated during the RAFFLE runs, these are option - Python 3.11 or higher - raffle package installed (`pip install .`) - `ase` for handling atomic structure data -- `CHGNet`, `MACE-0`, `VASP`, or some other `ase`-compatible calculator +- `CHGNet`, `MACE`, `VASP`, or some other `ase`-compatible calculator diff --git a/example/python_pkg/Si-Ge_learn/DRAFFLE/learn.py b/example/python_pkg/Si-Ge_learn/DRAFFLE/learn.py index a943831..3fb78e9 100644 --- a/example/python_pkg/Si-Ge_learn/DRAFFLE/learn.py +++ b/example/python_pkg/Si-Ge_learn/DRAFFLE/learn.py @@ -7,6 +7,9 @@ import numpy as np import os from joblib import Parallel, delayed +from pathlib import Path + +script_dir = Path(__file__).resolve().parent import logging logging.basicConfig(level=logging.DEBUG) @@ -16,7 +19,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct # Check if the structure has already been processed if i < num_structures_old: return None, None, None - + # calc = Vasp(**calc_params, label=f"struct{i}", directory=f"iteration{iteration}/struct{i}/", txt=f"stdout{i}.o") inew = i - num_structures_old atoms.calc = calc @@ -33,7 +36,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct except Exception as e: print(f"Optimisation failed: {e}") return None, None, None - + # Save the optimised structure and its energy per atom energy_rlxd = atoms.get_potential_energy() / len(atoms) @@ -47,24 +50,24 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct if distances.min() < 1.0: print(f"Distance too small: {atoms.get_all_distances(mic=True).min()}") return None, None, None - + if abs(energy_rlxd - energy_unrlxd) > 10.0: print(f"Energy difference too large: {energy_rlxd} vs {energy_unrlxd}") return None, None, None - + return atoms, energy_unrlxd, energy_rlxd if __name__ == "__main__": # check if mace file exists - if not os.path.exists("../mace-mpa-0-medium.model"): + if not os.path.exists(script_dir / ".." / ".." / "mace-mpa-0-medium.model"): print("MACE-MPA-0 model file not found. Please download the model from the MACE website.") print("https://github.com/ACEsuit/mace-foundations/releases/tag/mace_mpa_0") exit(1) # set up the calculator - calc_params = { 'model': "../mace-mpa-0-medium.model" } + calc_params = { 'model': script_dir/ ".." / ".." / "mace-mpa-0-medium.model" } calc = mace_mp(**calc_params) # set up the hosts @@ -123,7 +126,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct # check if the energies file exists, if not create it energies_rlxd_filename = f"energies_rlxd_seed{seed}.txt" energies_unrlxd_filename = f"energies_unrlxd_seed{seed}.txt" - + if os.path.exists(energies_rlxd_filename): with open(energies_rlxd_filename, "w") as energy_file: pass @@ -135,7 +138,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct pass else: open(energies_unrlxd_filename, "w").close() - + # initialise the number of structures generated iter = -1 unrlxd_structures = [] @@ -185,7 +188,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct energy=generated_structures[num_structures_old + i].get_potential_energy(), forces=generated_structures[num_structures_old + i].get_forces() ) - + # Start parallel execution print("Starting parallel execution") results = Parallel(n_jobs=5)( @@ -218,7 +221,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct del unrlxd_structures[j] del rlxd_structures[j] generator.remove_structure(j) - num_structures_new = len(generated_structures) + num_structures_new = len(generated_structures) # write the structures to files for i in range(num_structures_new - num_structures_old): @@ -268,4 +271,4 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct write(f"rlxd_structures_seed{seed}.traj", rlxd_structures) print("All generated and relaxed structures written") - print("Learning complete") \ No newline at end of file + print("Learning complete") diff --git a/example/python_pkg/graphene_grain_boundary_learn/DRAFFLE/learn.py b/example/python_pkg/graphene_grain_boundary_learn/DRAFFLE/learn.py index 6d15356..1e65854 100644 --- a/example/python_pkg/graphene_grain_boundary_learn/DRAFFLE/learn.py +++ b/example/python_pkg/graphene_grain_boundary_learn/DRAFFLE/learn.py @@ -8,6 +8,9 @@ import os from joblib import Parallel, delayed from ase.constraints import FixAtoms +from pathlib import Path + +script_dir = Path(__file__).resolve().parent import logging logging.basicConfig(level=logging.DEBUG) @@ -59,24 +62,24 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct if __name__ == "__main__": # check if mace file exists - if not os.path.exists("../mace-mpa-0-medium.model"): + if not os.path.exists(script_dir / ".." / ".." / "mace-mpa-0-medium.model"): print("MACE-MPA-0 model file not found. Please download the model from the MACE website.") print("https://github.com/ACEsuit/mace-foundations/releases/tag/mace_mpa_0") exit(1) # set up the calculator - calc_params = { 'model': "../mace-mpa-0-medium.model" } + calc_params = { 'model': script_dir/ ".." / ".." / "mace-mpa-0-medium.model" } calc = mace_mp(**calc_params) # set up the hosts hosts = [] - host = read("../POSCAR_host_gb") + host = read(script_dir / ".." / "POSCAR_host_gb") hosts.append(host) # set the parameters for the generator generator = raffle_generator() generator.distributions.set_history_len(10) - graphene = read("../POSCAR_graphene") + graphene = read(script_dir / ".." / "POSCAR_graphene") h2 = build.molecule("H2") graphene.calc = calc C_reference_energy = graphene.get_potential_energy() / len(graphene) diff --git a/example/python_pkg/visualisation/descriptor.ipynb b/example/python_pkg/visualisation/descriptor.ipynb new file mode 100644 index 0000000..48b4e0e --- /dev/null +++ b/example/python_pkg/visualisation/descriptor.ipynb @@ -0,0 +1,289 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "2010c09a", + "metadata": {}, + "outputs": [], + "source": [ + "# Iterative RAFFLE structure search\n", + "from ase.io import read\n", + "from ase import Atoms\n", + "from raffle.generator import raffle_generator\n", + "from mace.calculators import mace_mp\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pandas as pd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac87b758", + "metadata": {}, + "outputs": [], + "source": [ + "generator = raffle_generator()\n", + "generator.distributions.set_history_len(10)\n", + "calc = mace_mp(model=\"medium\", dispersion=False, default_dtype=\"float32\", device='cpu')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb5834ef", + "metadata": {}, + "outputs": [], + "source": [ + "host = Atoms('C', positions=[(0, 0, 0)], cell=[10, 10, 10], pbc=True)\n", + "host.calc = calc\n", + "generator.set_host(host)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9c4c302", + "metadata": {}, + "outputs": [], + "source": [ + "aa_stack = read(\"POSCAR_AA_stack\")\n", + "ab_stack = read(\"POSCAR_AB_stack\")\n", + "aabbcc_stack = read(\"POSCAR_AABBCC_stack\")\n", + "aba_stack = read(\"POSCAR_ABA_stack\")\n", + "abab_stack = read(\"POSCAR_ABAB_stack\")\n", + "lonsdaleite = read(\"POSCAR_lonsdaleite\")\n", + "diamond = read(\"POSCAR_diamond\")\n", + "\n", + "# 0 = Unknown mp entry\n", + "# 1 = mp-169 = 0.001 OR mp-3347313 = 0.000\n", + "# 2 = mp-2516584 = 0.002\n", + "# 3 = mp-606949 = 0.006\n", + "# 4 = mp-569416 = 0.002\n", + "# 5 = mp-47 = 0.139\n", + "# 6 = mp-66 = 0.112" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5a5d183b", + "metadata": {}, + "outputs": [], + "source": [ + "aa_stack.calc = calc\n", + "ab_stack.calc = calc\n", + "aabbcc_stack.calc = calc\n", + "aba_stack.calc = calc\n", + "abab_stack.calc = calc\n", + "lonsdaleite.calc = calc\n", + "diamond.calc = calc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7fe5192", + "metadata": {}, + "outputs": [], + "source": [ + "generator.distributions.set_element_energies(\n", + " {\n", + " \"C\": ab_stack.get_potential_energy() / len(ab_stack),\n", + " }\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "139ae70d", + "metadata": {}, + "outputs": [], + "source": [ + "database = [\n", + " aa_stack,\n", + " ab_stack,\n", + " aabbcc_stack,\n", + " aba_stack,\n", + " abab_stack,\n", + " lonsdaleite,\n", + " diamond,\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e711155", + "metadata": {}, + "outputs": [], + "source": [ + "### Optional parameters\n", + "generator.distributions.set_kBT(0.00001)\n", + "# generator.distributions.set_width([0.02, np.pi/200.0, np.pi/200.0])\n", + "# generator.distributions.set_cutoff_min([0.5, 0.0, 0.0])\n", + "# generator.distributions.set_cutoff_max([6.0, np.pi, np.pi])\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff7ef8c2", + "metadata": {}, + "outputs": [], + "source": [ + "### Explore the effect of excluding 3- and 4-body terms\n", + "# generator.distributions.set_radius_distance_tol([0.0, 0.0, 0.0, 0.0]) # 2-body\n", + "# generator.distributions.set_radius_distance_tol([1.5, 2.5, 0.0, 0.0]) # 2+3-body\n", + "# generator.distributions.set_radius_distance_tol([1.5, 2.5, 3.0, 6.0]) # 2+3+4-body" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2677f485", + "metadata": {}, + "outputs": [], + "source": [ + "generator.distributions.create(database)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7579711c", + "metadata": {}, + "outputs": [], + "source": [ + "descriptor = generator.get_descriptor()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7bf54a55", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a figure with 3 subplots side by side\n", + "fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n", + "\n", + "# Plot for each n-body descriptor (2-body, 3-body, 4-body)\n", + "for j in range(3):\n", + " # Calculate x-axis values\n", + " x = np.arange(generator.distributions.cutoff_min[j],\n", + " generator.distributions.cutoff_max[j] + generator.distributions.width[j],\n", + " generator.distributions.width[j])\n", + "\n", + " # Plot on the respective subplot\n", + " for idx in range(len(descriptor[j])):\n", + " axes[j].plot(x, descriptor[j][idx,:])\n", + "\n", + " # Set labels and title for each subplot\n", + " axes[j].set_ylabel('Descriptor value')\n", + " axes[j].set_title(f'{j+2}-body descriptor')\n", + "\n", + "axes[0].set_xlabel('Distance (Å)')\n", + "axes[1].set_xlabel('3-body angle (radians)')\n", + "axes[2].set_xlabel('Improper dihedral angle (radians)')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36d87c69", + "metadata": {}, + "outputs": [], + "source": [ + "structures_dict = {\n", + " \"AA stacked\": aa_stack,\n", + " \"AB stacked\": ab_stack,\n", + " \"AABBCC stacked\": aabbcc_stack,\n", + " \"ABA stacked\": aba_stack,\n", + " \"ABAB stacked\": abab_stack,\n", + " \"lonsdaleite\": lonsdaleite,\n", + " \"diamond\": diamond,\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27d97cbe", + "metadata": {}, + "outputs": [], + "source": [ + "n_body_dict = {\n", + " \"2-body\": [0.0, 0.0, 0.0, 0.0],\n", + " \"2+3-body\": [1.5, 2.5, 0.0, 0.0],\n", + " \"2+3+4-body\": [1.5, 2.5, 3.0, 6.0],\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d42a7a66", + "metadata": {}, + "outputs": [], + "source": [ + "viability_dict = {}\n", + "for struc_name, structure in structures_dict.items():\n", + " viability_dict[struc_name] = {\n", + " \"total energy\": structure.get_potential_energy()/len(structure),\n", + " \"formation energy\": structure.get_potential_energy()/len(structure) - ab_stack.get_potential_energy()/len(ab_stack),\n", + " }\n", + "for key, value in n_body_dict.items():\n", + " generator.distributions.set_radius_distance_tol(value)\n", + " generator.distributions.create(database)\n", + "\n", + " for struc_name, structure in structures_dict.items():\n", + " viability_dict[struc_name][\"viability \"+key] = generator.evaluate(structure)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bfc1172b", + "metadata": {}, + "outputs": [], + "source": [ + "viability_df = pd.DataFrame(viability_dict).T\n", + "viability_df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23ddc5c7", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "raffle_env", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/example/python_pkg/visualisation/fingerprint.ipynb b/example/python_pkg/visualisation/fingerprint.ipynb new file mode 100644 index 0000000..a9d3f40 --- /dev/null +++ b/example/python_pkg/visualisation/fingerprint.ipynb @@ -0,0 +1,170 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "2010c09a", + "metadata": {}, + "outputs": [], + "source": [ + "# Iterative RAFFLE structure search\n", + "from ase.io import read\n", + "from raffle.generator import raffle_generator\n", + "from mace.calculators import mace_mp\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac87b758", + "metadata": {}, + "outputs": [], + "source": [ + "generator = raffle_generator()\n", + "calc = mace_mp(model=\"medium\", dispersion=False, default_dtype=\"float32\", device='cpu')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9c4c302", + "metadata": {}, + "outputs": [], + "source": [ + "aa_stack = read(\"POSCAR_AA_stack\")\n", + "ab_stack = read(\"POSCAR_AB_stack\")\n", + "aabbcc_stack = read(\"POSCAR_AABBCC_stack\")\n", + "aba_stack = read(\"POSCAR_ABA_stack\")\n", + "abab_stack = read(\"POSCAR_ABAB_stack\")\n", + "lonsdaleite = read(\"POSCAR_lonsdaleite\")\n", + "diamond = read(\"POSCAR_diamond\")\n", + "\n", + "# 0 = Unknown mp entry\n", + "# 1 = mp-169 = 0.001 OR mp-3347313 = 0.000\n", + "# 2 = mp-2516584 = 0.002\n", + "# 3 = mp-606949 = 0.006\n", + "# 4 = mp-569416 = 0.002\n", + "# 5 = mp-47 = 0.139\n", + "# 6 = mp-66 = 0.112" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5a5d183b", + "metadata": {}, + "outputs": [], + "source": [ + "aa_stack.calc = calc\n", + "ab_stack.calc = calc\n", + "aabbcc_stack.calc = calc\n", + "aba_stack.calc = calc\n", + "abab_stack.calc = calc\n", + "lonsdaleite.calc = calc\n", + "diamond.calc = calc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40882bbe", + "metadata": {}, + "outputs": [], + "source": [ + "# change aa_stack to same height as ab_stack\n", + "# this is to confirm that the differences in descriptor are not a result of different interlayer distances\n", + "ab_height = ab_stack.get_cell()[2, 2]\n", + "cell = aa_stack.get_cell()\n", + "aa_height = cell[2, 2]\n", + "cell[2, 2] = ab_height\n", + "aa_stack.set_cell(cell, scale_atoms=True)\n", + "aa_stack.calc = calc" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2677f485", + "metadata": {}, + "outputs": [], + "source": [ + "aa_descriptor = generator.distributions.generate_fingerprint(aa_stack)\n", + "ab_descriptor = generator.distributions.generate_fingerprint(ab_stack)\n", + "aba_descriptor = generator.distributions.generate_fingerprint(aba_stack)\n", + "abab_descriptor = generator.distributions.generate_fingerprint(abab_stack)\n", + "aabbcc_descriptor = generator.distributions.generate_fingerprint(aabbcc_stack)\n", + "lonsdaleite_descriptor = generator.distributions.generate_fingerprint(lonsdaleite)\n", + "diamond_descriptor = generator.distributions.generate_fingerprint(diamond)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7bf54a55", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a figure with 3 subplots side by side\n", + "fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n", + "\n", + "# Plot for each n-body function (2-body, 3-body, 4-body)\n", + "for j in range(3):\n", + " # Calculate x-axis values\n", + " x = np.arange(generator.distributions.cutoff_min[j],\n", + " generator.distributions.cutoff_max[j] + generator.distributions.width[j],\n", + " generator.distributions.width[j])\n", + "\n", + " # Plot on the respective subplot\n", + " # axes[j].plot(x, aa_descriptor[j], label='AA stack')\n", + " axes[j].plot(x, ab_descriptor[j], label='AB stack')\n", + " # axes[j].plot(x, aabbcc_descriptor[j], label='AABBCC stack')\n", + " axes[j].plot(x, aba_descriptor[j], label='ABA stack')\n", + " axes[j].plot(x, abab_descriptor[j], label='ABAB stack')\n", + " # axes[j].plot(x, lonsdaleite_descriptor[j], label='Lonsdaleite')\n", + " # axes[j].plot(x, diamond_descriptor[j], label='Diamond')\n", + "\n", + " # Set labels and title for each subplot\n", + " axes[j].set_ylabel('Fingerprint value')\n", + " axes[j].set_title(f'{j+2}-body fingerprint')\n", + " axes[j].legend()\n", + "\n", + "axes[0].set_xlabel('Distance (Å)')\n", + "axes[1].set_xlabel('3-body angle (radians)')\n", + "axes[2].set_xlabel('Improper dihedral angle (radians)')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7a215eb5", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "raffle_env", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/example/python_pkg/visualisation/probability_density.ipynb b/example/python_pkg/visualisation/probability_density.ipynb new file mode 100644 index 0000000..39b3ab4 --- /dev/null +++ b/example/python_pkg/visualisation/probability_density.ipynb @@ -0,0 +1,463 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "8f4453d0", + "metadata": {}, + "outputs": [], + "source": [ + "from raffle.generator import raffle_generator\n", + "from ase.io import read\n", + "from ase import build\n", + "import numpy as np\n", + "import re\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import colors\n", + "from scipy.interpolate import griddata\n", + "from scipy.spatial import cKDTree\n", + "import ipywidgets as widgets\n", + "from ipywidgets import interactive\n", + "%matplotlib widget\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.ticker import AutoMinorLocator\n", + "\n", + "from IPython.display import display\n", + "import ipywidgets as widgets\n", + "import tkinter as tk\n", + "from tkinter import filedialog" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9c2804ee", + "metadata": {}, + "outputs": [], + "source": [ + "# Set the font family and size to use for Matplotlib figures\n", + "plt.rc('text', usetex=True)\n", + "plt.rc('font', family='Computer Modern')\n", + "plt.rcParams.update({\n", + " \"figure.facecolor\": (1.0, 1.0, 1.0, 1.0), # red with alpha = 30%\n", + " \"axes.facecolor\": (1.0, 1.0, 1.0, 1.0), # green with alpha = 50%\n", + " \"savefig.facecolor\": (1.0, 1.0, 1.0, 0.0), # blue with alpha = 20%\n", + "})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa42e5ab", + "metadata": {}, + "outputs": [], + "source": [ + "# Load the MACE model\n", + "from pathlib import Path\n", + "from mace.calculators import mace_mp\n", + "\n", + "calc_params = { 'model': \"../mace-mpa-0-medium.model\" }\n", + "calc = mace_mp(**calc_params)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8ac96ea9", + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the generator\n", + "generator = raffle_generator()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e292ed4a", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# set up the host\n", + "Si_bulk = build.bulk(\"Si\", crystalstructure=\"diamond\", a=5.43)\n", + "Si_bulk.calc = calc\n", + "Si_reference_energy = Si_bulk.get_potential_energy() / len(Si_bulk)\n", + "Si_cubic = build.make_supercell(Si_bulk, [[-1, 1, 1], [1, -1, 1], [1, 1, -1]])\n", + "Ge_bulk = build.bulk(\"Ge\", crystalstructure=\"diamond\", a=5.65)\n", + "Ge_bulk.calc = calc\n", + "Ge_cubic = build.make_supercell(Ge_bulk, [[-1, 1, 1], [1, -1, 1], [1, 1, -1]])\n", + "Ge_reference_energy = Ge_bulk.get_potential_energy() / len(Ge_bulk)\n", + "\n", + "Si_supercell = build.make_supercell(Si_cubic, [[2, 0, 0], [0, 2, 0], [0, 0, 1]])\n", + "Ge_supercell = build.make_supercell(Ge_cubic, [[2, 0, 0], [0, 2, 0], [0, 0, 1]])\n", + "\n", + "Si_surface = build.surface(Si_supercell, indices=(0, 0, 1), layers=2)\n", + "Ge_surface = build.surface(Ge_supercell, indices=(0, 0, 1), layers=2)\n", + "\n", + "Si_slab = build.surface(Si_supercell, indices=(0, 0, 1), layers=2, vacuum=12, periodic=True)\n", + "Si_slab.calc = calc\n", + "Ge_slab = build.surface(Ge_supercell, indices=(0, 0, 1), layers=2, vacuum=12, periodic=True)\n", + "Ge_slab.calc = calc\n", + "\n", + "host = build.stack(Si_surface, Ge_surface, axis=2, distance= 5.43/2 + 5.65/2)\n", + "cell = host.get_cell()\n", + "cell[2, 2] -= 3.8865 # (5.43 + 5.65) / 2 * 3/4\n", + "host.set_cell(cell, scale_atoms=False)\n", + "\n", + "# Set up the database and the reference energies\n", + "generator.distributions.set_element_energies(\n", + " {\n", + " 'Si': Si_reference_energy,\n", + " 'Ge': Ge_reference_energy,\n", + " }\n", + ")\n", + "database = [ Si_bulk, Ge_bulk ]\n", + "database = read(\"../Si-Ge_learn/DRAFFLE/DOutput/rlxd_structures_seed0.traj\", index=\":\")\n", + "species = \"SiGe\"\n", + "# bounds = [[0, 0, 0.34], [1, 1, 0.52]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9f07a454", + "metadata": {}, + "outputs": [], + "source": [ + "# # Set up the database and the reference energies\n", + "# host = read(\"../graphene_grain_boundary_learn/POSCAR_host_gb\")\n", + "# graphene = read(\"../graphene_grain_boundary_learn/POSCAR_graphene\")\n", + "# h2 = build.molecule(\"H2\")\n", + "# database = [host]\n", + "# generator.distributions.set_element_energies(\n", + "# {\n", + "# 'C': 0.0,\n", + "# }\n", + "# )\n", + "# species = \"C\"\n", + "# bounds = [[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ed3573e", + "metadata": {}, + "outputs": [], + "source": [ + "# Learng the RAFFLE descriptor\n", + "generator.distributions.create(database, deallocate_systems=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f0b573a", + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the species list\n", + "species_name_list = re.findall(r'[A-Z][a-z]?\\d*', species)\n", + "species_name_list = [re.sub(r'\\d+', '', s) for s in species_name_list]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a85baf0", + "metadata": {}, + "outputs": [], + "source": [ + "# Return the probability density\n", + "probability_density, grid = generator.get_probability_density(host, species=species, grid_spacing=0.1, return_grid=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "65c7b982", + "metadata": {}, + "outputs": [], + "source": [ + "descriptor = generator.get_descriptor()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e192283a", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a figure with 3 subplots side by side\n", + "fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n", + "\n", + "# Plot for each n-body descriptor (2-body, 3-body, 4-body)\n", + "for j in range(3):\n", + " # Calculate x-axis values\n", + " x = np.arange(generator.distributions.cutoff_min[j],\n", + " generator.distributions.cutoff_max[j] + generator.distributions.width[j],\n", + " generator.distributions.width[j])\n", + "\n", + " # Plot on the respective subplot\n", + " for idx in range(len(descriptor[j])):\n", + " axes[j].plot(x, descriptor[j][idx,:])\n", + "\n", + " # Set labels and title for each subplot\n", + " axes[j].set_ylabel('Descriptor value')\n", + " axes[j].set_title(f'{j+2}-body descriptor')\n", + "\n", + "axes[0].set_xlabel('Distance (Å)')\n", + "axes[1].set_xlabel('3-body angle (radians)')\n", + "axes[2].set_xlabel('Improper dihedral angle (radians)')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9ca61bc", + "metadata": {}, + "outputs": [], + "source": [ + "# Conver the probability density to a meshgrid\n", + "species_list = range(probability_density.shape[0]-4)\n", + "\n", + "cell = host.get_cell()\n", + "a = np.linalg.norm(cell[0])\n", + "b = np.linalg.norm(cell[1])\n", + "c = np.linalg.norm(cell[2])\n", + "\n", + "x_min = a * ( 0.0 )\n", + "x_max = a * ( 1.0 - 1.0/grid[0] )\n", + "y_min = b * ( 0.0 )\n", + "y_max = b * ( 1.0 - 1.0/grid[1] )\n", + "z_min = c * ( 0.0 )\n", + "z_max = c * ( 1.0 - 1.0/grid[2] )\n", + "\n", + "# Create a 3D grid\n", + "grid_x, grid_y, grid_z = np.mgrid[x_min:x_max:complex(grid[0]), y_min:y_max:complex(grid[1]), z_min:z_max:complex(grid[2])]\n", + "grid_points = np.vstack((grid_x.ravel(), grid_y.ravel(), grid_z.ravel())).T" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "81708d84", + "metadata": {}, + "outputs": [], + "source": [ + "# Extract the coordinates and values\n", + "x = probability_density[0]\n", + "y = probability_density[1]\n", + "z = probability_density[2]\n", + "# scale the data positions\n", + "x = x * a\n", + "y = y * b\n", + "z = z * c\n", + "\n", + "void = probability_density[3]\n", + "\n", + "# Populate the grid with the values\n", + "grid_values = []\n", + "distance_threshold = 1e-3\n", + "for spec in sorted(set(species_list)):\n", + " values = probability_density[spec+4]\n", + "\n", + " # Calculate distances to the nearest known data point\n", + " tree = cKDTree(np.c_[x, y, z])\n", + " distances, _ = tree.query(grid_points, k=1)\n", + " \n", + " # Interpolate data onto the 3D grid\n", + " grid_values.append(griddata((x, y, z), values, (grid_x, grid_y, grid_z), method='nearest', fill_value=0.0))\n", + " # grid_values[spec] = np.nan_to_num(grid_values)\n", + "\n", + " # Reshape distances to match the grid shape\n", + " distances = distances.reshape(grid_values[spec-1].shape)\n", + "\n", + " # Set threshold for distance (e.g., 1 unit)\n", + " grid_values[spec-1][distances > distance_threshold] = 0 # Set to zero for points beyond the threshold" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11b8644c", + "metadata": {}, + "outputs": [], + "source": [ + "# Set up interactive sliders\n", + "slice_slider = widgets.FloatSlider(value=0.5, min=0, max=1-0.01, step=0.01, description='Slice Position:')\n", + "axis_selector = widgets.RadioButtons(options=['a', 'b', 'c', 'ab'], description='Axis:', value='c')\n", + "species_selector = widgets.Dropdown(options=sorted(set(species_list)), description='Species:')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5a6c43a8", + "metadata": {}, + "outputs": [], + "source": [ + "# Define the 2D plotting function\n", + "def plot_2d_slice(slice_level, species, axis):\n", + " plt.clf() # Clear the current figure\n", + "\n", + " norm = colors.PowerNorm(gamma=0.5)\n", + " plane_label = ''\n", + " # Normalize slice level to the grid range\n", + " if axis == 'c':\n", + " slice_idx = int(slice_level * grid[2])\n", + " slice_data = grid_values[species-1][:, :, slice_idx]\n", + " plt.imshow(slice_data.T, extent=(x_min, x_max, y_min, y_max), origin='lower', cmap='viridis', norm=norm, aspect='auto')\n", + " plt.xlabel('$x$ (\\AA)')\n", + " plt.ylabel('$y$ (\\AA)')\n", + " plane_label = '$p_{3}$'\n", + " elif axis == 'b':\n", + " slice_idx = int(slice_level * grid[1])\n", + " slice_data = grid_values[species-1][:, slice_idx, :]\n", + " plt.imshow(slice_data.T, extent=(x_min, x_max, z_min, z_max), origin='lower', cmap='viridis', norm=norm, aspect='auto')\n", + " plt.xlabel('$x$ (\\AA)')\n", + " plt.ylabel('$z$ (\\AA)')\n", + " plane_label = '$p_{4}$'\n", + " elif axis == 'a':\n", + " slice_idx = int(slice_level * grid[0])\n", + " slice_data = grid_values[species-1][slice_idx, :, :]\n", + " plt.imshow(slice_data.T, extent=(y_min, y_max, z_min, z_max), origin='lower', cmap='viridis', norm=norm, aspect='auto')\n", + " plt.xlabel('$y$ (\\AA)')\n", + " plt.ylabel('$z$ (\\AA)')\n", + " plane_label = '$p_{1}$'\n", + " elif axis == 'ab':\n", + " # get diagonal list of idx\n", + " slice_idx_xy = [\n", + " ( \n", + " round( ( -0.5 + slice_level ) * (grid[0] + 1) ) + i,\n", + " round( ( -0.5 + slice_level ) * (grid[1] + 1) ) + j\n", + " ) for i, j in zip(range(grid[0],-1,-1), range(grid[1]+1))\n", + " if ( round( ( -0.5 + slice_level ) * (grid[0] + 1) ) + i ) >= 0 and ( round( ( -0.5 + slice_level ) * (grid[0] + 1) ) + i ) < grid[0] and\n", + " ( round( ( -0.5 + slice_level ) * (grid[1] + 1) ) + j ) >= 0 and ( round( ( -0.5 + slice_level ) * (grid[1] + 1) ) + j ) < grid[1]\n", + " ]\n", + " slice_idx_x = [ i for i, j in slice_idx_xy ]\n", + " slice_idx_y = [ j for i, j in slice_idx_xy ]\n", + " slice_data = grid_values[species-1][slice_idx_x, slice_idx_y, :]\n", + " min_loc = 0.0\n", + " max_loc = np.sqrt( (len(slice_idx_x) * a / grid[0])**2 + (len(slice_idx_y) * b / grid[1])**2 )\n", + " plt.imshow(slice_data.T, extent=(min_loc, max_loc, z_min, z_max), origin='lower', cmap='viridis', norm=norm, aspect='auto')\n", + " plt.xlabel(r'$xy$ (\\AA)')\n", + " plt.ylabel('$z$ (\\AA)')\n", + " plane_label = '$p_{2}$'\n", + "\n", + " plt.title(f'Plane {plane_label}', fontsize=25, y=1.0, pad=10)\n", + " min_val = min(grid_value.min() for grid_value in grid_values)\n", + " max_val = max(grid_value.max() for grid_value in grid_values)\n", + " plt.clim(min_val, max_val)\n", + " cbar = plt.colorbar(label='Viability', orientation='vertical', fraction=0.085, pad=0.02)\n", + " ax = plt.gca()\n", + " # ax.set_aspect('equal')\n", + " # Set tick and label font size\n", + " ax.tick_params(axis='both', which='major', labelsize=20)\n", + " \n", + " # Create a custom colormap where zero values are colored grey\n", + " cmap = plt.cm.viridis.copy()\n", + " cmap.set_bad('grey')\n", + " slice_data = np.ma.masked_where(slice_data < 1e-8, slice_data)\n", + " plt.imshow(slice_data.T, extent=ax.get_images()[0].get_extent(), origin='lower', cmap=cmap, norm=norm, aspect='auto')\n", + " ax.set_aspect('equal')\n", + " \n", + " plt.xlabel(plt.gca().get_xlabel(), fontsize=25)\n", + " plt.ylabel(plt.gca().get_ylabel(), fontsize=25)\n", + " ax.xaxis.set_major_locator(plt.MaxNLocator(4))\n", + " ax.yaxis.set_major_locator(plt.MaxNLocator(4))\n", + " ax.xaxis.set_minor_locator(AutoMinorLocator(2))\n", + " ax.yaxis.set_minor_locator(AutoMinorLocator(2))\n", + " cbar.ax.tick_params(labelsize=20)\n", + " cbar.set_label(r'$P_{\\mathrm{'+species_name_list[species-1]+r'}}$', fontsize=25)\n", + " # cbar.set_ticks(np.linspace(min_val, max_val, 4))\n", + " cbar.ax.yaxis.set_major_locator(plt.MaxNLocator(nbins=4))\n", + " cbar.ax.minorticks_on()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b567adce", + "metadata": {}, + "outputs": [], + "source": [ + "# Save the current plot as a PDF\n", + "def save_plot(b):\n", + " slice_level = slice_slider_2d.value\n", + " species = species_selector_2d.value\n", + " axis = axis_selector_2d.value\n", + " plot_2d_slice(slice_level, species, axis)\n", + " # Open file dialog to select file path and name\n", + " root = tk.Tk()\n", + " root.withdraw() # Hide the root window\n", + " file_path = filedialog.asksaveasfilename(defaultextension=\".pdf\", filetypes=[(\"PDF files\", \"*.pdf\")])\n", + " if file_path:\n", + " fig = plt.gcf()\n", + " fig.patch.set_facecolor('none') # Set the face color to white\n", + " ax = plt.gca()\n", + " ax.patch.set_facecolor('none') # Set the axes face color to white\n", + " plt.savefig(file_path, bbox_inches='tight', pad_inches=0, facecolor=fig.get_facecolor(), edgecolor='none')\n", + " print(f'Plot saved as {file_path}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec553454", + "metadata": {}, + "outputs": [], + "source": [ + "# Create the figure for the 2D plot\n", + "plt.figure(figsize=(8, 6))\n", + "\n", + "save_button = widgets.Button(description=\"Save Plot\")\n", + "save_button.on_click(save_plot)\n", + "display(save_button)\n", + "\n", + "# Set up interactive sliders for 2D plot\n", + "slice_slider_2d = widgets.FloatSlider(value=0.5, min=0, max=1-0.01, step=0.01, description='Slice Position:')\n", + "axis_selector_2d = widgets.RadioButtons(options=['a', 'b', 'c', 'ab'], description='Axis:', value='c')\n", + "species_selector_2d = widgets.Dropdown(options=sorted(set(species_list)), description='Species:')\n", + "\n", + "# Set up interactive plot for 2D slices\n", + "interactive_plot_2d = interactive(plot_2d_slice, slice_level=slice_slider_2d, axis=axis_selector_2d, species=species_selector_2d)\n", + "interactive_plot_2d" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a3cb8a3b", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "raffle_env", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/fortran/lib/mod_cache.f90 b/src/fortran/lib/mod_cache.f90 new file mode 100644 index 0000000..579e9a9 --- /dev/null +++ b/src/fortran/lib/mod_cache.f90 @@ -0,0 +1,32 @@ +module raffle__cache + use raffle__constants, only: real32 + implicit none + + private + public :: store_probability_density, retrieve_probability_density + + real(real32), allocatable, dimension(:,:), save :: cached_probability_density + +contains + + subroutine store_probability_density(probability_density) + implicit none + real(real32), intent(in) :: probability_density(:,:) + if (allocated(cached_probability_density)) & + deallocate(cached_probability_density) + allocate(cached_probability_density, source = probability_density) + + end subroutine store_probability_density + + function retrieve_probability_density() result(probability_density) + implicit none + 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 + else + allocate(probability_density, source = cached_probability_density) + end if + end function retrieve_probability_density + +end module raffle__cache diff --git a/src/fortran/lib/mod_constants.f90 b/src/fortran/lib/mod_constants.f90 index 78ca0a3..660a0e4 100644 --- a/src/fortran/lib/mod_constants.f90 +++ b/src/fortran/lib/mod_constants.f90 @@ -5,6 +5,7 @@ module raffle__constants !! library. implicit none integer, parameter, public :: real32 = Selected_real_kind(6,37)!(15,307) + real(real32), parameter, public :: tau = 8._real32 * atan(1._real32) real(real32), parameter, public :: pi = 4._real32 * atan(1._real32) real(real32), parameter, public :: c = 0.26246582250210965422_real32 real(real32), parameter, public :: c_vasp = 0.262465831_real32 diff --git a/src/fortran/lib/mod_dist_calcs.f90 b/src/fortran/lib/mod_dist_calcs.f90 index 6ea8c93..4ad1db0 100644 --- a/src/fortran/lib/mod_dist_calcs.f90 +++ b/src/fortran/lib/mod_dist_calcs.f90 @@ -5,7 +5,7 @@ module raffle__dist_calcs !! and other points in the system. use raffle__constants, only: pi,real32 use raffle__geom_rw, only: basis_type - use raffle__misc_linalg, only: modu, get_angle + use raffle__misc_linalg, only: get_angle implicit none @@ -20,8 +20,9 @@ module raffle__dist_calcs contains !############################################################################### - pure function get_min_dist(basis,loc,lignore_close,axis,labove,lreal,tol, & - ignore_list) result(output) + pure function get_min_dist( & + basis, loc, lignore_close, axis, labove, lreal, tol & + ) result(output) !! Return the minimum distance between a point and the nearest atom !! in a cell. !! @@ -44,8 +45,6 @@ pure function get_min_dist(basis,loc,lignore_close,axis,labove,lreal,tol, & !! The tolerance for the distance. logical, intent(in), optional :: labove, lreal !! If true, return the real distance, otherwise return the vector. - integer, dimension(:,:), intent(in), optional :: ignore_list - !! List of atoms to ignore. real(real32), dimension(3) :: output !! The minimum distance between the point and the nearest atom. @@ -84,13 +83,9 @@ pure function get_min_dist(basis,loc,lignore_close,axis,labove,lreal,tol, & output = 0._real32 do js = 1, basis%nspec atmloop: do ja = 1, basis%spec(js)%num - if(present(ignore_list))then - do i = 1, size(ignore_list,2), 1 - if(all(ignore_list(:,i).eq.[js,ja])) cycle atmloop - end do - end if + if(.not.basis%spec(js)%atom_mask(ja)) cycle atmloop vdtmp1 = basis%spec(js)%atom(ja,:3) - loc - if(lignore_close.and.modu(vdtmp1).lt.tol_) cycle atmloop + if(lignore_close.and.norm2(vdtmp1).lt.tol_) cycle atmloop if(axis_.gt.0)then if(abs(vdtmp1(axis_)).lt.tol_) cycle atmloop if(labove_)then @@ -102,7 +97,7 @@ pure function get_min_dist(basis,loc,lignore_close,axis,labove,lreal,tol, & vdtmp1 = vdtmp1 - ceiling(vdtmp1 - 0.5_real32) end if vdtmp2 = matmul(vdtmp1,basis%lat) - dtmp1 = modu(vdtmp2) + dtmp1 = norm2(vdtmp2) if(dtmp1.lt.min_bond)then min_bond = dtmp1 if(lreal_)then @@ -143,7 +138,7 @@ pure function get_min_dist_between_point_and_atom(basis,loc,atom) result(dist) vec = loc - basis%spec(atom(1))%atom(atom(2),:3) vec = vec - ceiling(vec - 0.5_real32) vec = matmul(vec,basis%lat) - dist = modu(vec) + dist = norm2(vec) end function get_min_dist_between_point_and_atom !############################################################################### @@ -151,7 +146,7 @@ end function get_min_dist_between_point_and_atom !############################################################################### pure function get_min_dist_between_point_and_species( & - basis, loc, species, ignore_list) result(dist) + basis, loc, species) result(dist) !! Return the minimum distance between a point and a species in a cell. !! !! This function returns the minimum distance between a point and any @@ -165,8 +160,6 @@ pure function get_min_dist_between_point_and_species( & !! The index of the species in the cell. real(real32), dimension(3), intent(in) :: loc !! The location of the point (in crystal coordinates). - integer, dimension(:,:), intent(in), optional :: ignore_list - !! List of atoms to ignore. real(real32) :: dist !! The minimum distance between the point and the species. @@ -181,15 +174,11 @@ pure function get_min_dist_between_point_and_species( & dist = huge(0._real32) atom_loop: do ia = 1,basis%spec(species)%num - if(present(ignore_list))then - do i = 1, size(ignore_list,2), 1 - if(all(ignore_list(:,i).eq.[species,ia])) cycle atom_loop - end do - end if + if(.not.basis%spec(species)%atom_mask(ia)) cycle atom_loop vec = loc - basis%spec(species)%atom(ia,:3) vec = vec - ceiling(vec - 0.5_real32) vec = matmul(vec, basis%lat) - rtmp1 = modu(vec) + rtmp1 = norm2(vec) if( rtmp1 .lt. dist ) dist = rtmp1 end do atom_loop @@ -220,7 +209,7 @@ pure function get_dist_between_point_and_atom(basis,loc,atom) result(dist) vec = loc - basis%spec(atom(1))%atom(atom(2),:3) vec = matmul(vec,basis%lat) - dist = modu(vec) + dist = norm2(vec) end function get_dist_between_point_and_atom !############################################################################### diff --git a/src/fortran/lib/mod_distribs.f90 b/src/fortran/lib/mod_distribs.f90 index 5e17c15..9a7882d 100644 --- a/src/fortran/lib/mod_distribs.f90 +++ b/src/fortran/lib/mod_distribs.f90 @@ -5,11 +5,11 @@ module raffle__distribs !! fucntions for individual materials. !! The distribution functions are used as fingerprints for atomic structures !! to identify similarities and differences between structures. - use raffle__constants, only: real32, pi + use raffle__constants, only: real32, pi, tau use raffle__io_utils, only: stop_program, print_warning use raffle__misc, only: strip_null, sort_str use raffle__misc_maths, only: triangular_number - use raffle__misc_linalg, only: get_angle, get_improper_dihedral_angle, modu + use raffle__misc_linalg, only: get_angle, get_improper_dihedral_angle use raffle__geom_rw, only: basis_type, get_element_properties use raffle__geom_extd, only: extended_basis_type use raffle__element_utils, only: & @@ -176,7 +176,7 @@ subroutine calculate(this, basis, & !! Loop index. integer :: num_pairs !! Number of pairs and angles. - real(real32) :: bondlength + real(real32) :: bondlength, rtmp1, dist_max_smooth, dist_min_smooth !! Temporary real variables. logical :: success !! Boolean for success. @@ -186,6 +186,8 @@ subroutine calculate(this, basis, & !! Basis for storing neighbour data. real(real32), dimension(3) :: eta !! Parameters for the distribution functions. + real(real32), dimension(4) :: tolerances + !! Tolerance for the distance between atoms for 3- and 4-body. real(real32), allocatable, dimension(:) :: angle_list, bondlength_list, & distance !! Temporary real arrays. @@ -298,10 +300,10 @@ subroutine calculate(this, basis, & allocate(neighbour_basis%spec(1)) allocate(neighbour_basis%image_spec(1)) allocate(neighbour_basis%spec(1)%atom( & - sum(basis_extd%spec(:)%num)+sum(basis_extd%image_spec(:)%num), 3 & + sum(basis_extd%spec(:)%num)+sum(basis_extd%image_spec(:)%num), 4 & ) ) allocate(neighbour_basis%image_spec(1)%atom( & - sum(basis_extd%spec(:)%num)+sum(basis_extd%image_spec(:)%num), 3 & + sum(basis_extd%spec(:)%num)+sum(basis_extd%image_spec(:)%num), 4 & ) ) neighbour_basis%nspec = basis%nspec neighbour_basis%natom = 0 @@ -319,6 +321,12 @@ subroutine calculate(this, basis, & neighbour_basis%image_spec(1)%num = 0 do js = 1, basis%nspec itmp1 = 0 + tolerances(:) = radius_distance_tol_(:) * & + bond_info(pair_index(is, js))%radius_covalent + tolerances(1) = max( cutoff_min_(1), tolerances(1) ) + tolerances(3) = max( cutoff_min_(1), tolerances(3) ) + tolerances(2) = min( cutoff_max_(1), tolerances(2) ) + tolerances(4) = min( cutoff_max_(1), tolerances(4) ) !------------------------------------------------------------------ ! loop over all atoms inside the unit cell @@ -331,7 +339,7 @@ subroutine calculate(this, basis, & basis_extd%spec(is)%atom(ia,1:3) & ], basis_extd%lat ) & ) - bondlength = modu( vector ) + bondlength = norm2( vector ) if( & bondlength .lt. cutoff_min_(1) .or. & @@ -341,37 +349,45 @@ subroutine calculate(this, basis, & ! add 2-body bond to store if within tolerances for 3-body ! distance if( & - bondlength .ge. ( & - bond_info(pair_index(is, js))%radius_covalent * & - radius_distance_tol_(1) & - ) .and. & - bondlength .le. ( & - bond_info(pair_index(is, js))%radius_covalent * & - radius_distance_tol_(2) & - ) & + bondlength .ge. tolerances(1) .and. & + bondlength .le. tolerances(2) & ) then neighbour_basis%spec(1)%num = & neighbour_basis%spec(1)%num + 1 neighbour_basis%spec(1)%atom( & - neighbour_basis%spec(1)%num,1:3) = vector + neighbour_basis%spec(1)%num,1:3 & + ) = vector + neighbour_basis%spec(1)%atom( & + neighbour_basis%spec(1)%num,4 & + ) = -0.5_real32 * ( & + cos( tau * ( bondlength - tolerances(1) ) / & + ( & + min(cutoff_max_(1), tolerances(2)) - & + tolerances(1) & + ) & + ) - 1._real32 ) end if ! add 2-body bond to store if within tolerances for 4-body ! distance if( & - bondlength .ge. ( & - bond_info(pair_index(is, js))%radius_covalent * & - radius_distance_tol_(3) & - ) .and. & - bondlength .le. ( & - bond_info(pair_index(is, js))%radius_covalent * & - radius_distance_tol_(4) & - ) & + bondlength .ge. tolerances(3) .and. & + bondlength .le. tolerances(4) & ) then neighbour_basis%image_spec(1)%num = & neighbour_basis%image_spec(1)%num + 1 neighbour_basis%image_spec(1)%atom( & - neighbour_basis%image_spec(1)%num,1:3) = vector + neighbour_basis%image_spec(1)%num,1:3 & + ) = vector + neighbour_basis%image_spec(1)%atom( & + neighbour_basis%image_spec(1)%num,4 & + ) = -0.5_real32 * ( & + cos( tau * ( bondlength - tolerances(3) ) / & + ( & + min(cutoff_max_(1), tolerances(4)) - & + tolerances(3) & + ) & + ) - 1._real32 ) end if !if(js.lt.js.or.(is.eq.js.and.ja.le.ia)) cycle @@ -394,7 +410,7 @@ subroutine calculate(this, basis, & ], basis_extd%lat ) & ) - bondlength = modu( vector ) + bondlength = norm2( vector ) if( & bondlength .lt. cutoff_min_(1) .or. & @@ -404,39 +420,45 @@ subroutine calculate(this, basis, & ! add 2-body bond to store if within tolerances for 3-body ! distance if( & - bondlength .ge. ( & - bond_info(pair_index(is, js))%radius_covalent * & - radius_distance_tol_(1) & - ) .and. & - bondlength .le. ( & - bond_info(pair_index(is, js))%radius_covalent * & - radius_distance_tol_(2) & - ) & + bondlength .ge. tolerances(1) .and. & + bondlength .le. tolerances(2) & ) then neighbour_basis%spec(1)%num = & neighbour_basis%spec(1)%num + 1 neighbour_basis%spec(1)%atom( & neighbour_basis%spec(1)%num,1:3 & ) = vector + neighbour_basis%spec(1)%atom( & + neighbour_basis%spec(1)%num,4 & + ) = -0.5_real32 * ( & + cos( tau * ( bondlength - tolerances(1) ) / & + ( & + min(cutoff_max_(1), tolerances(2)) - & + tolerances(1) & + ) & + ) - 1._real32 ) end if ! add 2-body bond to store if within tolerances for 4-body ! distance if( & - bondlength .ge. ( & - bond_info(pair_index(is, js))%radius_covalent * & - radius_distance_tol_(3) & - ) .and. & - bondlength .le. ( & - bond_info(pair_index(is, js))%radius_covalent * & - radius_distance_tol_(4) & - ) & + bondlength .ge. tolerances(3) .and. & + bondlength .le. tolerances(4) & ) then neighbour_basis%image_spec(1)%num = & neighbour_basis%image_spec(1)%num + 1 neighbour_basis%image_spec(1)%atom( & neighbour_basis%image_spec(1)%num,1:3 & ) = vector + neighbour_basis%image_spec(1)%atom( & + neighbour_basis%image_spec(1)%num,4 & + ) = -0.5_real32 * ( & + cos( tau * ( bondlength - tolerances(3) ) / & + ( & + min(cutoff_max_(1), tolerances(4)) - & + tolerances(3) & + ) & + ) - 1._real32 ) end if itmp1 = itmp1 + 1 @@ -504,8 +526,11 @@ subroutine calculate(this, basis, & ) distance(idx) = & ( & - modu(neighbour_basis%spec(1)%atom(ja,:3)) ** 2 * & - modu(neighbour_basis%spec(1)%atom(ka,:3)) ** 2 & + neighbour_basis%spec(1)%atom(ja,4) * & + neighbour_basis%spec(1)%atom(ka,4) & + ) / ( & + norm2(neighbour_basis%spec(1)%atom(ja,:3)) ** 2 * & + norm2(neighbour_basis%spec(1)%atom(ka,:3)) ** 2 & ) end do end do @@ -555,9 +580,15 @@ subroutine calculate(this, basis, & [ neighbour_basis%image_spec(1)%atom(la,:3) ] & ) distance(idx) = & - modu(neighbour_basis%spec(1)%atom(ja,:3)) ** 2 * & - modu(neighbour_basis%spec(1)%atom(ka,:3)) ** 2 * & - modu(neighbour_basis%image_spec(1)%atom(la,:3)) ** 2 + ( & + neighbour_basis%spec(1)%atom(ja,4) * & + neighbour_basis%spec(1)%atom(ka,4) * & + neighbour_basis%image_spec(1)%atom(la,4) & + ) / ( & + norm2(neighbour_basis%spec(1)%atom(ja,:3)) ** 2 * & + norm2(neighbour_basis%spec(1)%atom(ka,:3)) ** 2 * & + norm2(neighbour_basis%image_spec(1)%atom(la,:3)) ** 2 & + ) end do end do ! a NaN in the angle refers to one where two of the vectors are @@ -583,9 +614,24 @@ subroutine calculate(this, basis, & !--------------------------------------------------------------------------- ! apply the cutoff function to the 2-body distribution function !--------------------------------------------------------------------------- + dist_max_smooth = cutoff_max_(1) - 0.25_real32 + dist_min_smooth = cutoff_min_(1) + 0.25_real32 do b = 1, nbins_(1) - this%df_2body(b,:) = this%df_2body(b,:) / ( cutoff_min_(1) + & - width_(1) * real(b-1, real32) ) ** 2 + rtmp1 = cutoff_min_(1) + width_(1) * real(b-1, real32) + this%df_2body(b,:) = this%df_2body(b,:) / rtmp1 ** 2 + if( rtmp1 .gt. dist_max_smooth )then + this%df_2body(b,:) = this%df_2body(b,:) * 0.5_real32 * & + ( 1._real32 + cos( pi * & + ( rtmp1 - dist_max_smooth ) / & + ( cutoff_max_(1) - dist_max_smooth ) & + ) ) + elseif( rtmp1 .lt. dist_min_smooth )then + this%df_2body(b,:) = this%df_2body(b,:) * 0.5_real32 * & + ( 1._real32 + cos( pi * & + ( rtmp1 - dist_min_smooth ) / & + ( dist_min_smooth - cutoff_min_(1) ) & + ) ) + end if end do @@ -691,7 +737,7 @@ function get_distrib(value_list, nbins, eta, width, cutoff_min, & value_list(i) - & ( width * real(b-1, real32) + cutoff_min ) & ) ** 2._real32 & - ) / scale_list(i) + ) * scale_list(i) end do end do end do diff --git a/src/fortran/lib/mod_distribs_container.f90 b/src/fortran/lib/mod_distribs_container.f90 index 593aa03..5ac3e74 100644 --- a/src/fortran/lib/mod_distribs_container.f90 +++ b/src/fortran/lib/mod_distribs_container.f90 @@ -47,9 +47,12 @@ module raffle__distribs_container !! above the hull. If false, the formation energy from the element !! reference energies is used. real(real32) :: & + viability_2body_default = 0.1_real32, & viability_3body_default = 0.1_real32, & viability_4body_default = 0.1_real32 - !! Default viability for the 3- and 4-body distribution functions. + !! Default viability for the 2-, 3-, and 4-body distribution functions. + logical :: smooth_viability = .true. + !! DEV FEATURE. Boolean whether to smooth the viability evaluation. logical, dimension(:), allocatable :: & in_dataset_2body, in_dataset_3body, in_dataset_4body !! Whether the 2-, 3-, and 4-body distribution functions are in @@ -67,6 +70,8 @@ module raffle__distribs_container real(real32), dimension(3) :: & width = [ 0.025_real32, pi/64._real32, pi/64._real32 ] !! Width of the bins used in the 2-, 3-, and 4-body distribution functions. + real(real32), dimension(3) :: width_inv + !! Inverse of the width of the bins used in the 2-, 3-, and 4-body real(real32), dimension(3) :: & cutoff_min = [ 0.5_real32, 0._real32, 0._real32 ] !! Minimum cutoff for the 2-, 3-, and 4-body distribution functions. @@ -138,7 +143,6 @@ module raffle__distribs_container !! Return the energies of elements in the container. !! Used in Python interface. - procedure, pass(this) :: set_element_map !! Set the mapping of elements to distribution function elements. procedure, pass(this), private :: set_bond_info @@ -155,7 +159,6 @@ module raffle__distribs_container !! Return the radii of all bonds in the container. !! Used in Python interface. - procedure, pass(this) :: set_best_energy !! Set the best energy and system in the container. procedure, pass(this) :: initialise_gdfs @@ -167,7 +170,6 @@ module raffle__distribs_container procedure, pass(this) :: is_converged !! Check if the learned distribution function has converged. - procedure, pass(this) :: write_gdfs !! Write the generalised distribution functions to a file. procedure, pass(this) :: read_gdfs @@ -186,6 +188,8 @@ module raffle__distribs_container !! Return the index for bond_info given two elements. procedure, pass(this) :: get_element_index !! Return the index for element_info given one element. + procedure, pass(this) :: set_num_bins + !! Set the number of bins for the n-body distribution functions. procedure, pass(this) :: get_bin !! Return the bin index for a given distance. procedure, pass(this) :: get_2body @@ -194,6 +198,11 @@ module raffle__distribs_container !! Return the 3-body distribution function. procedure, pass(this) :: get_4body !! Return the 4-body distribution function. + + procedure, pass(this) :: generate_fingerprint + !! Calculate the distribution functions for a given system. + procedure, pass(this) :: generate_fingerprint_python + !! Calculate the distribution functions for a given system. end type distribs_container_type interface distribs_container_type @@ -285,9 +294,13 @@ module function init_distribs_container( & if(present(history_len)) distribs_container%history_len = history_len + if(allocated(distribs_container%history_deltas)) & + deallocate(distribs_container%history_deltas) if(distribs_container%history_len.ge.0) & allocate( & - distribs_container%history_deltas(history_len), & + distribs_container%history_deltas( & + distribs_container%history_len & + ), & source = huge(0._real32) & ) @@ -1201,6 +1214,61 @@ end function get_4body !############################################################################### +!############################################################################### + function generate_fingerprint(this, structure) result(output) + !! Generate a descriptor for the structure. + implicit none + + ! Arguments + class(distribs_container_type), intent(inout) :: this + !! Parent. Instance of distribution functions container. + type(basis_type), intent(in) :: structure + !! Structure to generate the descriptor for. + type(distribs_type) :: output + !! Descriptor for the structure. + + call output%calculate( & + structure, & + width = this%width, & + sigma = this%sigma, & + cutoff_min = this%cutoff_min, & + cutoff_max = this%cutoff_max, & + radius_distance_tol = this%radius_distance_tol & + ) + + end function generate_fingerprint +!------------------------------------------------------------------------------- + subroutine generate_fingerprint_python( & + this, structure, output_2body, output_3body, output_4body & + ) + !! Generate a descriptor for the structure. + implicit none + + ! Arguments + class(distribs_container_type), intent(inout) :: this + !! Parent. Instance of distribution functions container. + type(basis_type), intent(in) :: structure + !! Structure to generate the descriptor for. + real(real32), dimension(:,:), intent(out) :: output_2body + !! 2-body descriptor for the structure. + real(real32), dimension(:,:), intent(out) :: output_3body + !! 3-body descriptor for the structure. + real(real32), dimension(:,:), intent(out) :: output_4body + !! 4-body descriptor for the structure. + + ! Local variables + type(distribs_type) :: distrib + !! Descriptor for the structure. + + distrib = this%generate_fingerprint(structure) + output_2body = distrib%df_2body + output_3body = distrib%df_3body + output_4body = distrib%df_4body + + end subroutine generate_fingerprint_python +!############################################################################### + + !############################################################################### subroutine add(this, system) !! Add a system to the container. @@ -2195,14 +2263,41 @@ pure function get_element_index(this, species) result(idx) integer :: idx !! Index of the element in the element_info array. - ! Local variables - idx = findloc([ this%element_info(:)%name ], species, dim=1) end function get_element_index !############################################################################### +!############################################################################### + subroutine set_num_bins(this) + !! Set the number of bins for the n-body distribution functions. + implicit none + + ! Arguments + class(distribs_container_type), intent(inout) :: this + !! Parent of the procedure. Instance of distribution functions container. + + ! Local variables + integer :: i + !! Loop index. + + do i = 1, 3 + if(this%nbins(i).eq.-1)then + this%nbins(i) = 1 + & + nint( & + ( this%cutoff_max(i) - this%cutoff_min(i) ) / & + this%width(i) & + ) + end if + this%width_inv(i) = ( this%nbins(i) - 1 ) / & + ( this%cutoff_max(i) - this%cutoff_min(i) ) + end do + + end subroutine set_num_bins +!############################################################################### + + !############################################################################### pure function get_bin(this, value, dim) result(bin) !! Get the bin index for a value in a dimension. @@ -2218,16 +2313,19 @@ pure function get_bin(this, value, dim) result(bin) integer :: bin !! Bin index for the value. - if(value .lt. this%cutoff_min(dim) - this%width(dim) .or. & - value .gt. this%cutoff_max(dim) + this%width(dim))then - bin = 0 - else - bin = nint( & - ( this%nbins(dim) - 1 ) * & - ( value - this%cutoff_min(dim) ) / & - ( this%cutoff_max(dim) - this%cutoff_min(dim) ) & - ) + 1 - end if + ! Local variables + real(real32) :: min_val, width_inv + !! Temporary variables. + + ! Prefetch frequently accessed values + min_val = this%cutoff_min(dim) + width_inv = this%width_inv(dim) + + ! Calculate bin using optimized operations + bin = nint((value - min_val) * width_inv) + 1 + + ! Ensure bin stays within bounds (floating point safety) + bin = min(max(bin, 1), this%nbins(dim)) end function get_bin !############################################################################### @@ -2352,15 +2450,7 @@ subroutine evolve(this, system) ! if present, add the system to the container !--------------------------------------------------------------------------- if(present(system)) call this%add(system) - do i = 1, 3 - if(this%nbins(i).eq.-1)then - this%nbins(i) = 1 + & - nint( & - ( this%cutoff_max(i) - this%cutoff_min(i) ) / & - this%width(i) & - ) - end if - end do + call this%set_num_bins() !--------------------------------------------------------------------------- @@ -2681,6 +2771,8 @@ subroutine evolve(this, system) this%num_evaluated_allocated = size(this%system) this%num_evaluated = this%num_evaluated + num_evaluated + this%viability_2body_default = sum( this%gdf%df_2body ) / & + real( size( this%gdf%df_2body ), real32 ) this%viability_3body_default = sum( this%gdf%df_3body ) / & real( size( this%gdf%df_3body ), real32 ) this%viability_4body_default = sum( this%gdf%df_4body ) / & diff --git a/src/fortran/lib/mod_evaluator.f90 b/src/fortran/lib/mod_evaluator.f90 index 93fa64d..e0c0cda 100644 --- a/src/fortran/lib/mod_evaluator.f90 +++ b/src/fortran/lib/mod_evaluator.f90 @@ -5,7 +5,7 @@ module raffle__evaluator !! the system with each point in the map representing the suitability of !! that point for a new atom. The map is built by checking the bond lengths, !! bond angles and dihedral angles between the test point and all atoms. - use raffle__constants, only: real32 + use raffle__constants, only: real32, tau, pi use raffle__misc_linalg, only: modu, get_angle, get_improper_dihedral_angle use raffle__geom_extd, only: extended_basis_type use raffle__distribs_container, only: distribs_container_type @@ -20,7 +20,7 @@ module raffle__evaluator !############################################################################### pure function evaluate_point( distribs_container, & - position, species, basis, atom_ignore_list, & + position, species, basis, & radius_list & ) result(output) !! Return the viability of a point in a basis for a specified species @@ -39,23 +39,34 @@ pure function evaluate_point( distribs_container, & !! Basis of the system. real(real32), dimension(3), intent(in) :: position !! Position of the test point. - integer, dimension(:,:), intent(in) :: atom_ignore_list - !! List of atoms to ignore (i.e. indices of atoms not yet placed). real(real32), dimension(:), intent(in) :: radius_list !! List of radii for each pair of elements. real(real32) :: output !! Suitability of the test point. ! Local variables - integer :: i, is, js, ia, ja + integer :: i, is, js, ia, ja, ia_end, ja_start, is_end !! Loop counters. + integer :: element_idx + !! Index of the query element. integer :: num_2body, num_3body, num_4body !! Number of 2-, 3- and 4-body interactions. real(real32) :: viability_2body !! Viability of the test point for 2-body interactions. real(real32) :: viability_3body, viability_4body !! Viability of the test point for 3- and 4-body interactions. - real(real32) :: bondlength + real(real32) :: bondlength, rtmp1, min_distance, min_from_max_cutoff + !! Temporary variables. + logical :: has_4body + !! Boolean whether the system has 4-body interactions. + real(real32) :: cos_scale1, cos_scale2 + !! Cosine scales for the 3- and 4-body interactions. + real(real32), dimension(3) :: position_1 + !! Cartesian coordinates of the test point. + real(real32), dimension(4) :: position_2 + !! Cartesian coordinates of the second atom and its cutoff weighting. + real(real32), dimension(4) :: tolerances + !! Tolerance for the distance between atoms for 3- and 4-body. integer, dimension(:,:), allocatable :: pair_index !! Index of element pairs. type(extended_basis_type) :: neighbour_basis @@ -65,6 +76,8 @@ pure function evaluate_point( distribs_container, & ! Initialisation output = 0._real32 viability_2body = 0._real32 + min_distance = distribs_container%cutoff_max(1) + min_from_max_cutoff = 0._real32 !--------------------------------------------------------------------------- @@ -91,15 +104,21 @@ pure function evaluate_point( distribs_container, & num_2body = 0 species_loop: do is = 1, basis%nspec allocate(neighbour_basis%spec(is)%atom( & - basis%spec(is)%num+basis%image_spec(is)%num, & - size(basis%spec(is)%atom,2) & + basis%spec(is)%num+basis%image_spec(is)%num, 4 & ) ) allocate(neighbour_basis%image_spec(is)%atom( & - basis%spec(is)%num+basis%image_spec(is)%num, & - size(basis%spec(is)%atom,2) & + basis%spec(is)%num+basis%image_spec(is)%num, 4 & ) ) neighbour_basis%spec(is)%num = 0 neighbour_basis%image_spec(is)%num = 0 + tolerances = distribs_container%radius_distance_tol * & + radius_list(pair_index(species,is)) + tolerances(1) = max( distribs_container%cutoff_min(1), tolerances(1) ) + tolerances(3) = max( distribs_container%cutoff_min(1), tolerances(3) ) + tolerances(2) = min( distribs_container%cutoff_max(1), tolerances(2) ) + tolerances(4) = min( distribs_container%cutoff_max(1), tolerances(4) ) + cos_scale1 = tau / ( tolerances(2) - tolerances(1) ) + cos_scale2 = tau / ( tolerances(4) - tolerances(3) ) !------------------------------------------------------------------------ ! 2-body map ! check bondlength between test point and all other atoms @@ -107,25 +126,16 @@ pure function evaluate_point( distribs_container, & atom_loop: do ia = 1, basis%spec(is)%num ! Check if the atom is in the ignore list ! If it is, skip the atom. - do i = 1, size(atom_ignore_list,dim=2), 1 - if(all(atom_ignore_list(:,i).eq.[is,ia])) cycle atom_loop - end do + if(.not.basis%spec(is)%atom_mask(ia)) cycle atom_loop associate( position_store => [ basis%spec(is)%atom(ia,1:3) ] ) - bondlength = modu( matmul(position - position_store, basis%lat) ) + bondlength = norm2( matmul(position - position_store, basis%lat) ) if( bondlength .gt. distribs_container%cutoff_max(1) ) & cycle atom_loop - if( bondlength .lt. max( & - distribs_container%cutoff_min(1), & - radius_list(pair_index(species,is)) * & - distribs_container%radius_distance_tol(1) & - ) )then + if( bondlength .lt. tolerances(1) )then ! If the bond length is less than the minimum allowed bond, ! return 0 (i.e. the point is not viable). return - elseif( bondlength .le. ( & - radius_list(pair_index(species,is)) * & - distribs_container%radius_distance_tol(2) & - ) )then + elseif( bondlength .le. tolerances(2) )then ! If the bond length is within the tolerance of the covalent ! radius of the pair, add the atom to the list of ! atoms to consider for 3-body interactions. @@ -133,20 +143,14 @@ pure function evaluate_point( distribs_container, & neighbour_basis%spec(is)%atom( & neighbour_basis%spec(is)%num,:3 & ) = matmul(position_store, basis%lat) + neighbour_basis%spec(is)%atom( & + neighbour_basis%spec(is)%num,4 & + ) = 0.5_real32 * abs( 1._real32 - & + cos( cos_scale1 * ( bondlength - tolerances(1) ) ) ) end if - if( bondlength .ge. & - ( & - radius_list(pair_index(species,is)) * & - distribs_container%radius_distance_tol(3) & - ) .and. & - bondlength .le. min( & - distribs_container%cutoff_max(1), & - ( & - radius_list(pair_index(species,is)) * & - distribs_container%radius_distance_tol(4) & - ) & - ) & + if( bondlength .ge. tolerances(3) .and. & + bondlength .le. tolerances(4) & )then ! If the bond length is within the min and max allowed size, ! add the atom to the list of atoms to consider for 4-body. @@ -155,11 +159,23 @@ pure function evaluate_point( distribs_container, & neighbour_basis%image_spec(is)%atom( & neighbour_basis%image_spec(is)%num,:3 & ) = matmul(position_store, basis%lat) + neighbour_basis%image_spec(is)%atom( & + neighbour_basis%image_spec(is)%num,4 & + ) = 0.5_real32 * abs( 1._real32 - & + cos( cos_scale2 * ( bondlength - tolerances(3) ) ) ) end if !------------------------------------------------------------------ ! Add the contribution of the bond length to the viability !------------------------------------------------------------------ + if(bondlength - tolerances(1) .lt. min_distance)then + min_distance = bondlength - tolerances(1) + end if + if(distribs_container%cutoff_max(1) - bondlength .gt. & + min_from_max_cutoff)then + min_from_max_cutoff = distribs_container%cutoff_max(1) - & + bondlength + end if viability_2body = viability_2body + & evaluate_2body_contributions( & distribs_container, bondlength, pair_index(species,is) & @@ -175,48 +191,47 @@ pure function evaluate_point( distribs_container, & !------------------------------------------------------------------------ image_loop: do ia = 1, basis%image_spec(is)%num, 1 associate( position_store => [ basis%image_spec(is)%atom(ia,1:3) ] ) - bondlength = modu( matmul(position - position_store, basis%lat) ) + bondlength = norm2( matmul(position - position_store, basis%lat) ) if( bondlength .gt. distribs_container%cutoff_max(1) ) & cycle image_loop - if( bondlength .lt. max( & - distribs_container%cutoff_min(1), & - radius_list(pair_index(species,is)) * & - distribs_container%radius_distance_tol(1) & - ) )then + if( bondlength .lt. tolerances(1) )then return - elseif( bondlength .le. ( & - radius_list(pair_index(species,is)) * & - distribs_container%radius_distance_tol(2) & - ) )then + elseif( bondlength .le. tolerances(2) )then neighbour_basis%spec(is)%num = neighbour_basis%spec(is)%num + 1 neighbour_basis%spec(is)%atom( & neighbour_basis%spec(is)%num,:3 & ) = matmul(position_store, basis%lat) + neighbour_basis%spec(is)%atom( & + neighbour_basis%spec(is)%num,4 & + ) = 0.5_real32 * ( 1._real32 - & + cos( cos_scale1 * ( bondlength - tolerances(1) ) ) ) end if - if( bondlength .ge. & - ( & - radius_list(pair_index(species,is)) * & - distribs_container%radius_distance_tol(3) & - ) .and. & - bondlength .le. min( & - distribs_container%cutoff_max(1), & - ( & - radius_list(pair_index(species,is)) * & - distribs_container%radius_distance_tol(4) & - ) & - ) & + if( bondlength .ge. tolerances(3) .and. & + bondlength .le. tolerances(4) & )then neighbour_basis%image_spec(is)%num = & neighbour_basis%image_spec(is)%num + 1 neighbour_basis%image_spec(is)%atom( & neighbour_basis%image_spec(is)%num,:3 & ) = matmul(position_store, basis%lat) + neighbour_basis%image_spec(is)%atom( & + neighbour_basis%image_spec(is)%num,4 & + ) = 0.5_real32 * abs( 1._real32 - & + cos( cos_scale2 * ( bondlength - tolerances(3) ) ) ) end if !------------------------------------------------------------------ ! Add the contribution of the bond length to the viability !------------------------------------------------------------------ + if(bondlength - tolerances(1) .lt. min_distance)then + min_distance = bondlength - tolerances(1) + end if + if(distribs_container%cutoff_max(1) - bondlength .gt. & + min_from_max_cutoff)then + min_from_max_cutoff = distribs_container%cutoff_max(1) - & + bondlength + end if viability_2body = viability_2body + & evaluate_2body_contributions( & distribs_container, bondlength, pair_index(species,is) & @@ -224,6 +239,19 @@ pure function evaluate_point( distribs_container, & num_2body = num_2body + 1 end associate end do image_loop + !------------------------------------------------------------------------ + ! DEVELOPER TOOL: This conditional allows the user to retrieve + ! results closer to first arXiv paper release. + ! It is not recommended to use this option for production runs and is + ! only here for testing purposes. + !------------------------------------------------------------------------ + if(.not.distribs_container%smooth_viability)then + neighbour_basis%spec(is)%atom(1:neighbour_basis%spec(is)%num,4) = & + 1._real32 + neighbour_basis%image_spec(is)%atom( & + 1:neighbour_basis%image_spec(is)%num,4 & + ) = 1._real32 + end if end do species_loop neighbour_basis%natom = sum(neighbour_basis%spec(:)%num) @@ -235,13 +263,24 @@ pure function evaluate_point( distribs_container, & ! This does not matter as, if there are no 2-body bonds, the point is ! not meant to be included in the viability set. ! The evaluator cannot comment on the viability of the point. - viability_2body = 0.5_real32 + viability_2body = distribs_container%viability_2body_default else viability_2body = viability_2body / real( num_2body, real32 ) + if(min_distance .lt. 0.25_real32 )then + viability_2body = viability_2body * 0.5_real32 * & + ( 1._real32 - cos( pi * min_distance / 0.25_real32 ) ) + end if + if( min_from_max_cutoff .lt. 0.5_real32 )then + rtmp1 = 0.5_real32 * ( 1._real32 - & + cos( pi * min_from_max_cutoff / 0.5_real32 ) ) + viability_2body = & + distribs_container%viability_2body_default * & + abs( 1._real32 - rtmp1 ) + & + rtmp1 * viability_2body + end if end if - ! store 3-body viable atoms in neighbour_basis%spec ! store 4-body viable atoms in neighbour_basis%image_spec ! for 3-bdoy, just cycle over neighbour_basis%spec @@ -251,44 +290,51 @@ pure function evaluate_point( distribs_container, & num_4body = 0 viability_3body = 1._real32 viability_4body = 1._real32 - do is = 1, neighbour_basis%nspec - do ia = 1, neighbour_basis%spec(is)%num + position_1 = matmul(position, basis%lat) + element_idx = distribs_container%element_map(species) + has_4body = any(neighbour_basis%image_spec(:)%num .gt. 0) + is_end = 0 + do is = 1, neighbour_basis%nspec, 1 + if(neighbour_basis%spec(is)%num .gt. 0) is_end = is + end do + do is = 1, is_end, 1 + ia_end = neighbour_basis%spec(is)%num - merge( 1, 0, is .eq. is_end ) + do ia = 1, ia_end, 1 + position_2 = neighbour_basis%spec(is)%atom(ia,1:4) !--------------------------------------------------------------------- ! 3-body map ! check bondangle between test point and all other atoms !--------------------------------------------------------------------- - associate( & - position_1 => matmul(position, basis%lat), & - position_2 => [neighbour_basis%spec(is)%atom(ia,1:3)] & - ) - if(sum(basis%spec(is:)%num).eq.ia) cycle - num_3body = num_3body + 1 - viability_3body = viability_3body * & - evaluate_3body_contributions( distribs_container, & - position_1, & - position_2, & - neighbour_basis, species, [is, ia] & - ) + viability_3body = viability_3body * max( & + 0._real32, & + evaluate_3body_contributions( distribs_container, & + position_1, & + position_2, & + neighbour_basis, element_idx, [is, ia] & + ) ) + num_3body = num_3body + 1 + if( has_4body )then do js = is, neighbour_basis%nspec - do ja = 1, neighbour_basis%spec(js)%num - if(js.eq.is .and. ja.le.ia) cycle - if(all(neighbour_basis%image_spec(:)%num.eq.0))cycle - num_4body = num_4body + 1 + ja_start = merge(ia + 1, 1, js .eq. is) + do ja = ja_start, neighbour_basis%spec(js)%num, 1 !------------------------------------------------------------ ! 4-body map ! check improperdihedral angle between test point and all ! other atoms !------------------------------------------------------------ - viability_4body = viability_4body * & + rtmp1 = max( & + 0._real32, & evaluate_4body_contributions( distribs_container, & position_1, & position_2, & - [neighbour_basis%spec(js)%atom(ja,1:3)], & - neighbour_basis, species & - ) + [ neighbour_basis%spec(js)%atom(ja,1:4) ], & + neighbour_basis, element_idx & + ) ) + viability_4body = viability_4body * rtmp1 + num_4body = num_4body + 1 end do end do - end associate + end if end do end do @@ -296,19 +342,19 @@ pure function evaluate_point( distribs_container, & !--------------------------------------------------------------------------- ! Normalise the angular viabilities !--------------------------------------------------------------------------- - if(num_3body.eq.0)then - viability_3body = distribs_container%viability_3body_default - else + if(num_3body.gt.0)then viability_3body = viability_3body ** ( & 1._real32 / real(num_3body,real32) & ) - end if - if(num_4body.eq.0)then - viability_4body = distribs_container%viability_4body_default else + viability_3body = distribs_container%viability_3body_default + end if + if(num_4body.gt.0)then viability_4body = viability_4body ** ( & 1._real32 / real(num_4body,real32) & ) + else + viability_4body = distribs_container%viability_4body_default end if @@ -350,7 +396,7 @@ end function evaluate_2body_contributions !############################################################################### pure function evaluate_3body_contributions( distribs_container, & - position_1, position_2, basis, species, current_idx & + position_1, position_2, basis, element_idx, current_idx & ) result(output) !! Return the contribution to the viability from 3-body interactions implicit none @@ -358,11 +404,13 @@ pure function evaluate_3body_contributions( distribs_container, & ! Arguments type(distribs_container_type), intent(in) :: distribs_container !! Distribution function (gvector) container. - real(real32), dimension(3), intent(in) :: position_1, position_2 - !! Positions of the atoms. + real(real32), dimension(3), intent(in) :: position_1 + !! Positions of the atom. + real(real32), dimension(4), intent(in) :: position_2 + !! Positions of the second atom and its cutoff weighting. type(extended_basis_type), intent(in) :: basis !! Basis of the system. - integer, intent(in) :: species + integer, intent(in) :: element_idx !! Index of the query element. integer, dimension(2), intent(in) :: current_idx !! Index of the 1st-atom query element. @@ -376,28 +424,34 @@ pure function evaluate_3body_contributions( distribs_container, & !! Bin for the distribution function. integer :: num_3body_local !! Number of 3-body interactions local to the current atom pair. + real(real32) :: power + !! Power for the contribution to the viability. + integer :: start_idx + !! Start index for the loop over the atoms. + real(real32) :: rtmp1, weight_2 + !! Temporary variables. - output = 1._real32 num_3body_local = sum(basis%spec(current_idx(1):)%num) - current_idx(2) - if(num_3body_local.eq.0) return + output = 1._real32 + power = 1._real32 / real( num_3body_local, real32 ) + weight_2 = sqrt( position_2(4) ) species_loop: do js = current_idx(1), basis%nspec, 1 - atom_loop: do ja = 1, basis%spec(js)%num - if(js.eq.current_idx(1) .and. ja.le.current_idx(2))cycle - associate( position_store => [ basis%spec(js)%atom(ja,1:3) ] ) - bin = distribs_container%get_bin( & - get_angle( & - position_2, & - position_1, & - position_store & - ), dim = 2 & - ) - output = output * & - distribs_container%gdf%df_3body( & - bin, & - distribs_container%element_map(species) & - ) ** ( 1._real32 / real( num_3body_local, real32 ) ) - end associate + start_idx = merge(current_idx(2) + 1, 1, js .eq. current_idx(1)) + atom_loop: do ja = start_idx, basis%spec(js)%num + bin = distribs_container%get_bin( & + get_angle( & + position_2(1:3), & + position_1, & + [ basis%spec(js)%atom(ja,1:3) ] & + ), dim = 2 & + ) + rtmp1 = weight_2 * sqrt( basis%spec(js)%atom(ja,4) ) + output = output * ( & + distribs_container%viability_3body_default * & + abs( 1._real32 - rtmp1 ) + & + rtmp1 * distribs_container%gdf%df_3body( bin, element_idx ) & + ) ** power end do atom_loop end do species_loop @@ -407,18 +461,19 @@ end function evaluate_3body_contributions !############################################################################### pure function evaluate_4body_contributions( distribs_container, & - position_1, position_2, position_3, basis, species ) result(output) + position_1, position_2, position_3, basis, element_idx ) result(output) !! Return the contribution to the viability from 4-body interactions implicit none ! Arguments type(distribs_container_type), intent(in) :: distribs_container !! Distribution function (gvector) container. - real(real32), dimension(3), intent(in) :: position_1, position_2, position_3 + real(real32), dimension(3), intent(in) :: position_1 !! Positions of the atoms. + real(real32), dimension(4), intent(in) :: position_2, position_3 type(extended_basis_type), intent(in) :: basis !! Basis of the system. - integer, intent(in) :: species + integer, intent(in) :: element_idx !! Index of the query element. real(real32) :: output !! Contribution to the viability. @@ -430,28 +485,33 @@ pure function evaluate_4body_contributions( distribs_container, & !! Bin for the distribution function. integer :: num_4body_local !! Number of 4-body interactions local to the current atom triplet. + real(real32) :: power + !! Power for the contribution to the viability. + real(real32) :: rtmp1, third, weight_2_3 + !! Temporary variables. - output = 1._real32 num_4body_local = sum(basis%image_spec(:)%num) - if(num_4body_local.eq.0) return + output = 1._real32 + power = 1._real32 / real( num_4body_local, real32 ) + third = 1._real32 / 3._real32 + weight_2_3 = ( position_2(4) * position_3(4) ) ** third species_loop: do ks = 1, basis%nspec, 1 atom_loop: do ka = 1, basis%image_spec(ks)%num - associate( position_store => [ basis%image_spec(ks)%atom(ka,1:3) ] ) - bin = distribs_container%get_bin( & - get_improper_dihedral_angle( & - position_1, & - position_2, & - position_3, & - position_store & - ), dim = 3 & - ) - output = output * & - distribs_container%gdf%df_4body( & - bin, & - distribs_container%element_map(species) & - ) ** ( 1._real32 / real( num_4body_local, real32 ) ) - end associate + bin = distribs_container%get_bin( & + get_improper_dihedral_angle( & + position_1, & + position_2(1:3), & + position_3(1:3), & + [ basis%image_spec(ks)%atom(ka,1:3) ] & + ), dim = 3 & + ) + rtmp1 = weight_2_3 * ( basis%image_spec(ks)%atom(ka,4) ) ** third + output = output * ( & + distribs_container%viability_4body_default * & + abs( 1._real32 - rtmp1 ) + & + rtmp1 * distribs_container%gdf%df_4body( bin, element_idx ) & + ) ** power end do atom_loop end do species_loop diff --git a/src/fortran/lib/mod_generator.f90 b/src/fortran/lib/mod_generator.f90 index ee15e2e..1172247 100644 --- a/src/fortran/lib/mod_generator.f90 +++ b/src/fortran/lib/mod_generator.f90 @@ -8,7 +8,6 @@ module raffle__generator use raffle__io_utils, only: stop_program, print_warning, suppress_warnings use raffle__constants, only: real32 use raffle__tools_infile, only: assign_val, assign_vec - use raffle__misc_linalg, only: modu use raffle__misc, only: strip_null, set, shuffle, sort1D, sort2D, to_upper use raffle__geom_rw, only: basis_type use raffle__geom_extd, only: extended_basis_type @@ -74,13 +73,18 @@ module raffle__generator walk_step_size_coarse = 1._real32, & walk_step_size_fine = 0.1_real32 !! Step size for the walk and grow methods. + real(real32), dimension(5) :: method_ratio_default = & + [1.0_real32, 0.1_real32, 0.5_real32, 0.5_real32, 1.0_real32] + !! Default ratio of each placement method. real(real32), dimension(5) :: method_ratio - !! Ratio of each placement method. + !! Last used ratio of each placement method. type(basis_type), dimension(:), allocatable :: structures !! Generated structures. contains procedure, pass(this) :: init_seed !! Procedure to set the seed for the random number generator. + procedure, pass(this) :: set_method_ratio_default + !! Procedure to set the ratio of each placement method. procedure, pass(this) :: set_host !! Procedure to set the host structure. @@ -108,6 +112,8 @@ module raffle__generator !! Procedure to remove a structure from the array of generated structures. procedure, pass(this) :: evaluate !! Procedure to evaluate the viability of a structure. + procedure, pass(this) :: get_probability_density + !! Procedure to get the probability density of a structure. procedure, pass(this) :: print_settings => print_generator_settings !! Procedure to print the raffle generator settings. @@ -241,6 +247,23 @@ end subroutine init_seed !############################################################################### +!############################################################################### + subroutine set_method_ratio_default(this, method_ratio) + !! Set the ratio of each placement method. + implicit none + + ! Arguments + class(raffle_generator_type), intent(inout) :: this + !! Instance of the raffle generator. + real(real32), dimension(5), intent(in) :: method_ratio + !! Ratio of each placement method. + + this%method_ratio_default = method_ratio + + end subroutine set_method_ratio_default +!############################################################################### + + !############################################################################### subroutine set_host(this, host) !! Set the host structure. @@ -337,7 +360,7 @@ function prepare_host( & depth_ = 3._real32 if(present(depth)) depth_ = depth call host%copy(this%host) - lattice_const = modu(host%lat(axis,:)) + lattice_const = norm2(host%lat(axis,:)) location_as_fractional_ = .false. if(present(location_as_fractional)) & location_as_fractional_ = location_as_fractional @@ -441,7 +464,7 @@ subroutine set_grid(this, grid, grid_spacing, grid_offset) do i = 1, 3 this%grid(i) = nint( & ( this%bounds(2,i) - this%bounds(1,i) ) * & - modu(this%host%lat(i,:)) / this%grid_spacing & + norm2(this%host%lat(i,:)) / this%grid_spacing & ) end do end if @@ -551,9 +574,7 @@ subroutine generate(this, num_structures, & !! Boolean to store the suppress_warnings value. type(basis_type) :: basis_template !! Basis of the structure to generate (i.e. allocated species and atoms). - real(real32), dimension(5) :: & - method_rand_limit = & - [1.0_real32, 0.1_real32, 0.5_real32, 0.5_real32, 1.0_real32] + real(real32), dimension(5) :: method_rand_limit !! Default ratio of each placement method. integer, dimension(:), allocatable :: seed_arr @@ -578,23 +599,14 @@ subroutine generate(this, num_structures, & !--------------------------------------------------------------------------- - ! Set the placement method selection limit numbers + ! Handle placement method optional argument !--------------------------------------------------------------------------- - if(verbose_.gt.0) write(*,*) "Setting method ratios" if(present(method_ratio))then method_rand_limit = method_ratio else - method_rand_limit = this%method_ratio + method_rand_limit = this%method_ratio_default end if this%method_ratio = method_rand_limit - ratio_norm = real(sum(method_rand_limit), real32) - method_rand_limit = method_rand_limit / ratio_norm - do i = 2, 5, 1 - method_rand_limit(i) = method_rand_limit(i) + method_rand_limit(i-1) - end do - if(verbose_.gt.0) write(*,*) & - "Method random limits (void, rand, walk, grow, min): ", & - method_rand_limit !--------------------------------------------------------------------------- @@ -606,6 +618,21 @@ subroutine generate(this, num_structures, & end if end if + + !--------------------------------------------------------------------------- + ! Set the placement method selection limit numbers + !--------------------------------------------------------------------------- + if(verbose_.gt.0) write(*,*) "Setting method ratio limits" + ratio_norm = real(sum(method_rand_limit), real32) + method_rand_limit = method_rand_limit / ratio_norm + do i = 2, 5, 1 + method_rand_limit(i) = method_rand_limit(i) + method_rand_limit(i-1) + end do + if(verbose_.gt.0) write(*,*) & + "Method random limits (void, rand, walk, grow, min): ", & + method_rand_limit + + !--------------------------------------------------------------------------- ! Set the random seed !--------------------------------------------------------------------------- @@ -647,6 +674,8 @@ subroutine generate(this, num_structures, & j = 0 do i = 1, basis_template%nspec + basis_template%spec(i)%atom_mask = & + [ ( .false., k = 1, basis_template%spec(i)%num, 1 ) ] basis_template%spec(i)%atom_idx = & [ ( k, k = j + 1, j + basis_template%spec(i)%num, 1 ) ] j = j + basis_template%spec(i)%num @@ -659,7 +688,10 @@ subroutine generate(this, num_structures, & call stop_program("Host structure not set") return end if - basis_template = basis_merge(this%host,basis_template) + basis_template = basis_merge( & + this%host, basis_template, & + mask1 = .true., mask2 = .false. & + ) basis_template%lat = this%host%lat @@ -781,7 +813,7 @@ function generate_structure( & !! Exit code. ! Local variables - integer :: iplaced, void_ticker + integer :: iplaced, void_ticker, i !! Loop counters. integer :: num_insert_atoms !! Number of atoms to insert. @@ -815,8 +847,7 @@ function generate_structure( & !--------------------------------------------------------------------------- call basis%copy(basis_initial) call basis%create_images( & - max_bondlength = this%distributions%cutoff_max(1), & - atom_ignore_list = placement_list & + max_bondlength = this%distributions%cutoff_max(1) & ) num_insert_atoms = basis%natom - this%host%natom @@ -846,7 +877,6 @@ function generate_structure( & basis, & species_index_list, & [ this%distributions%bond_info(:)%radius_covalent ], & - placement_list_shuffled, & this%grid_offset & ) @@ -871,8 +901,7 @@ function generate_structure( & basis, & species_index_list, & [ placement_list_shuffled(:,iplaced) ], & - [ this%distributions%bond_info(:)%radius_covalent ], & - placement_list_shuffled(:,iplaced+1:) & + [ this%distributions%bond_info(:)%radius_covalent ] & ) end if if(.not.allocated(gridpoint_viability))then @@ -892,18 +921,14 @@ function generate_structure( & call random_number(rtmp1) if(rtmp1.le.method_rand_limit_(1)) then if(verbose.gt.0) write(*,*) "Add Atom Void" - point = place_method_void( & - gridpoint_viability, & - basis, & - placement_list_shuffled(:,iplaced+1:), viable & - ) + point = place_method_void( gridpoint_viability, basis, viable ) elseif(rtmp1.le.method_rand_limit_(2)) then if(verbose.gt.0) write(*,*) "Add Atom Random" point = place_method_rand( & this%distributions, & this%bounds, & basis, & - placement_list_shuffled(:,iplaced+1:), & + placement_list_shuffled(1,iplaced+1), & [ this%distributions%bond_info(:)%radius_covalent ], & this%max_attempts, & viable & @@ -914,7 +939,7 @@ function generate_structure( & this%distributions, & this%bounds, & basis, & - placement_list_shuffled(:,iplaced+1:), & + placement_list_shuffled(1,iplaced+1), & [ this%distributions%bond_info(:)%radius_covalent ], & this%max_attempts, & this%walk_step_size_coarse, this%walk_step_size_fine, & @@ -927,7 +952,7 @@ function generate_structure( & this%distributions, & this%bounds, & basis, & - placement_list_shuffled(:,iplaced+1:), & + placement_list_shuffled(1,iplaced+1), & [ this%distributions%bond_info(:)%radius_covalent ], & this%max_attempts, & viable & @@ -942,7 +967,7 @@ function generate_structure( & placement_list_shuffled(1,iplaced), & this%bounds, & basis, & - placement_list_shuffled(:,iplaced+1:), & + placement_list_shuffled(1,iplaced+1), & [ this%distributions%bond_info(:)%radius_covalent ], & this%max_attempts, & this%walk_step_size_coarse, this%walk_step_size_fine, & @@ -988,10 +1013,7 @@ function generate_structure( & placement_aborted = .true. exit placement_loop elseif(void_ticker.gt.10)then - point = place_method_void( & - gridpoint_viability, basis, & - placement_list_shuffled(:,iplaced+1:), viable & - ) + point = place_method_void( gridpoint_viability, basis, viable ) void_ticker = 0 end if if(.not.viable) cycle placement_loop @@ -1006,6 +1028,8 @@ function generate_structure( & iplaced = iplaced + 1 basis%spec(placement_list_shuffled(1,iplaced))%atom( & placement_list_shuffled(2,iplaced),:3) = point(:3) + basis%spec(placement_list_shuffled(1,iplaced))%atom_mask( & + placement_list_shuffled(2,iplaced)) = .true. call basis%update_images( & max_bondlength = this%distributions%cutoff_max(1), & is = placement_list_shuffled(1,iplaced), & @@ -1156,10 +1180,8 @@ function evaluate(this, basis) result(viability) !! Viability of the generated structures. ! Local variables - integer :: is, ia, species + integer :: is, ia, idx !! Loop indices. - integer, dimension(2,1) :: atom_ignore - !! Atom to ignore. type(extended_basis_type) :: basis_extd !! Extended basis for the structure to evaluate. @@ -1172,11 +1194,12 @@ function evaluate(this, basis) result(viability) call this%distributions%set_element_map( & [ basis_extd%spec(:)%name ] & ) + call this%distributions%set_num_bins() do is = 1, basis%nspec - species = this%distributions%get_element_index( & + idx = this%distributions%get_element_index( & basis_extd%spec(is)%name & ) - if(species.eq.0)then + if(idx.eq.0)then call stop_program( & "Species "//& trim(basis_extd%spec(is)%name)//& @@ -1185,14 +1208,14 @@ function evaluate(this, basis) result(viability) return end if do ia = 1, basis%spec(is)%num - atom_ignore(:,1) = [is,ia] + basis_extd%spec(is)%atom_mask(ia) = .false. viability = viability + & evaluate_point( this%distributions, & [ basis%spec(is)%atom(ia,1:3) ], & - species, basis_extd, & - atom_ignore, & + is, basis_extd, & [ this%distributions%bond_info(:)%radius_covalent ] & ) + basis_extd%spec(is)%atom_mask(ia) = .true. end do end do @@ -1201,6 +1224,123 @@ end function evaluate !############################################################################### +!############################################################################### + function get_probability_density(this, basis, species_list, & + grid, grid_offset, grid_spacing, bounds, & + grid_output & + ) result(probability) + !! Get the probability density of the generated structures. + implicit none + + ! Arguments + class(raffle_generator_type), intent(inout) :: this + !! Instance of the raffle generator. + type(basis_type), intent(in) :: basis + !! Structure to evaluate. + character(len=3), dimension(:), intent(in) :: species_list + !! List of species to evaluate. + integer, dimension(3), intent(in), optional :: grid + !! Number of bins to divide the host structure into along each axis. + real(real32), dimension(3), intent(in), optional :: grid_offset + !! Offset of the gridpoints. + real(real32), intent(in), optional :: grid_spacing + !! Spacing of the bins. + real(real32), dimension(2,3), intent(in), optional :: bounds + !! Bounds for atom placement. + integer, dimension(3), intent(out), optional :: grid_output + + real(real32), dimension(:,:), allocatable :: probability + + ! Local variables + integer :: i, is, ia, idx + !! Loop indices. + real(real32) :: grid_spacing_ + !! Spacing of the bins. + integer, dimension(3) :: grid_ + !! Number of bins to divide the host structure into along each axis. + integer, dimension(size(species_list,1)) :: species_idx_list + real(real32), dimension(3) :: grid_offset_ + !! Offset of the gridpoints. + real(real32), dimension(2,3) :: bounds_ + !! Bounds for atom placement. + type(extended_basis_type) :: basis_extd + !! Extended basis for the structure to evaluate. + + + !--------------------------------------------------------------------------- + ! Set the grid and bounds + !--------------------------------------------------------------------------- + grid_ = -1 + grid_spacing_ = -1._real32 + grid_offset_ = 0._real32 + bounds_(1,:) = 0._real32 + bounds_(2,:) = 1._real32 + if(present(grid)) grid_ = grid + if(present(grid_offset)) grid_offset_ = grid_offset + if(present(grid_spacing)) grid_spacing_ = grid_spacing + if(present(bounds)) bounds_ = bounds + + if(any(grid_.eq.-1))then + if(grid_spacing_.lt.0._real32)then + call stop_program("Grid or grid spacing not set. One must be set") + return + end if + do i = 1, 3 + grid_(i) = nint( & + ( bounds_(2,i) - bounds_(1,i) ) * & + norm2(basis%lat(i,:)) / grid_spacing_ & + ) + end do + end if + if(present(grid_output)) grid_output = grid_ + call this%distributions%set_num_bins() + + + call basis_extd%copy(basis) + do i = 1, size(species_list) + call basis_extd%add_atom( & + species_list(i), & + position = [0._real32, 0._real32, 0._real32], & + mask = .false. & + ) + species_idx_list(i) = & + findloc(basis_extd%spec(:)%name, strip_null(species_list(i)), dim=1) + end do + do is = 1, basis_extd%nspec + idx = this%distributions%get_element_index( & + basis_extd%spec(is)%name & + ) + if(idx.eq.0)then + call stop_program( & + "Species "//& + trim(basis_extd%spec(is)%name)//& + " not found in distribution functions" & + ) + return + end if + end do + call basis_extd%create_images( & + max_bondlength = this%distributions%cutoff_max(1) & + ) + + + call this%distributions%set_element_map( & + [ basis_extd%spec(:)%name ] & + ) + probability = get_gridpoints_and_viability( & + this%distributions, & + grid_, & + bounds_, & + basis_extd, & + species_idx_list, & + [ this%distributions%bond_info(:)%radius_covalent ], & + grid_offset = grid_offset_ & + ) + + end function get_probability_density +!############################################################################### + + !############################################################################### subroutine print_generator_settings(this, file) !! Print the raffle generator settings. diff --git a/src/fortran/lib/mod_geom_extd.f90 b/src/fortran/lib/mod_geom_extd.f90 index 3b1d4dd..849ae7a 100644 --- a/src/fortran/lib/mod_geom_extd.f90 +++ b/src/fortran/lib/mod_geom_extd.f90 @@ -5,7 +5,7 @@ module raffle__geom_extd !! within a specified distance of the unit cell. This is useful for !! calculating interactions between atoms that are not within the unit cell. use raffle__constants, only: real32, pi - use raffle__misc_linalg, only: modu, cross, inverse_3x3 + use raffle__misc_linalg, only: cross, inverse_3x3 use raffle__geom_rw, only: basis_type, species_type implicit none @@ -33,7 +33,7 @@ module raffle__geom_extd contains !############################################################################### - subroutine create_images(this, max_bondlength, atom_ignore_list) + subroutine create_images(this, max_bondlength) !! Create the images for the basis set. implicit none @@ -42,8 +42,6 @@ subroutine create_images(this, max_bondlength, atom_ignore_list) !! Parent of the procedure. Instance of the extended basis. real(real32), intent(in) :: max_bondlength !! Maximum distance to extend the basis set. - integer, dimension(:,:), intent(in), optional :: atom_ignore_list - !! List of atoms to ignore when creating images. ! Local variables integer :: is, ia, i, j, k @@ -54,18 +52,6 @@ subroutine create_images(this, max_bondlength, atom_ignore_list) !! Temporary vector for storing atom positions. type(species_type), dimension(this%nspec) :: image_species !! Temporary store for the images. - integer, dimension(:,:), allocatable :: atom_ignore_list_ - !! List of atoms to ignore when creating images. - - - !--------------------------------------------------------------------------- - ! initialise the atom_ignore_list - !--------------------------------------------------------------------------- - if(present(atom_ignore_list)) then - atom_ignore_list_ = atom_ignore_list(:,:) - else - allocate(atom_ignore_list_(0,0)) - end if !--------------------------------------------------------------------------- @@ -74,9 +60,9 @@ subroutine create_images(this, max_bondlength, atom_ignore_list) ! won't work for extremely acute/obtuse angle cells ! (due to diagonal path being shorter than individual lattice vectors) !--------------------------------------------------------------------------- - amax = ceiling(max_bondlength/modu(this%lat(1,:))) - bmax = ceiling(max_bondlength/modu(this%lat(2,:))) - cmax = ceiling(max_bondlength/modu(this%lat(3,:))) + amax = ceiling(max_bondlength/norm2(this%lat(1,:))) + bmax = ceiling(max_bondlength/norm2(this%lat(2,:))) + cmax = ceiling(max_bondlength/norm2(this%lat(3,:))) spec_loop: do is = 1, this%nspec @@ -92,9 +78,7 @@ subroutine create_images(this, max_bondlength, atom_ignore_list) image_species(is)%radius = this%spec(is)%radius image_species(is)%name = this%spec(is)%name atom_loop: do ia = 1, this%spec(is)%num - do i = 1, size(atom_ignore_list_,2), 1 - if(all(atom_ignore_list_(:,i).eq.[is,ia])) cycle atom_loop - end do + if(.not.this%spec(is)%atom_mask(ia)) cycle atom_loop do i=-amax,amax+1,1 vtmp1(1) = this%spec(is)%atom(ia,1) + real(i, real32) do j=-bmax,bmax+1,1 @@ -168,9 +152,9 @@ subroutine update_images(this, max_bondlength, is, ia) ! (due to diagonal path being shorter than individual lattice vectors) !--------------------------------------------------------------------------- num_images = this%image_spec(is)%num - amax = ceiling(max_bondlength/modu(this%lat(1,:))) - bmax = ceiling(max_bondlength/modu(this%lat(2,:))) - cmax = ceiling(max_bondlength/modu(this%lat(3,:))) + amax = ceiling(max_bondlength/norm2(this%lat(1,:))) + bmax = ceiling(max_bondlength/norm2(this%lat(2,:))) + cmax = ceiling(max_bondlength/norm2(this%lat(3,:))) dim = 3 do i = 1, this%nspec if ( size(this%spec(i)%atom,2) .gt. dim) dim = size(this%spec(i)%atom,2) diff --git a/src/fortran/lib/mod_geom_rw.f90 b/src/fortran/lib/mod_geom_rw.f90 index c681105..d817cd5 100644 --- a/src/fortran/lib/mod_geom_rw.f90 +++ b/src/fortran/lib/mod_geom_rw.f90 @@ -6,7 +6,7 @@ module raffle__geom_rw use raffle__constants, only: pi,real32 use raffle__io_utils, only: stop_program, print_warning use raffle__misc, only: to_upper, to_lower, jump, icount, strip_null - use raffle__misc_linalg, only: modu, inverse_3x3 + use raffle__misc_linalg, only: inverse_3x3 implicit none @@ -34,6 +34,8 @@ module raffle__geom_rw integer, allocatable, dimension(:) :: atom_idx !! The indices of the atoms of this species in the basis. !! For ASE compatibility. + logical, dimension(:), allocatable :: atom_mask + !! The mask of the atoms of this species. real(real32), allocatable ,dimension(:,:) :: atom !! The atomic positions of the species. real(real32) :: mass @@ -72,8 +74,12 @@ module raffle__geom_rw !! Procedure to convert the basis to cartesian coordinates. procedure, pass(this) :: copy !! Procedure to copy the basis. + procedure, pass(this) :: set_atom_mask + !! Procedure to set the atom mask of the basis. procedure, pass(this) :: get_lattice_constants !! Procedure to get the lattice constants of the basis. + procedure, pass(this) :: add_atom + !! Procedure to add an atom to the basis. procedure, pass(this) :: remove_atom !! Procedure to remove an atom from the basis. procedure, pass(this) :: remove_atoms @@ -153,6 +159,7 @@ subroutine allocate_species( & istart = 1 do i = 1, this%nspec iend = istart + this%spec(i)%num - 1 + allocate(this%spec(i)%atom_mask(this%spec(i)%num), source = .true.) allocate(this%spec(i)%atom_idx(this%spec(i)%num)) allocate(this%spec(i)%atom(this%spec(i)%num,3)) if(present(atoms))then @@ -388,6 +395,7 @@ subroutine VASP_geom_read(UNIT, basis, length, iostat) natom = 0 do i = 1, basis%nspec allocate(basis%spec(i)%atom_idx(basis%spec(i)%num)) + allocate(basis%spec(i)%atom_mask(basis%spec(i)%num), source = .true.) allocate(basis%spec(i)%atom(basis%spec(i)%num,length_)) basis%spec(i)%atom(:,:)=0._real32 do j = 1, basis%spec(i)%num @@ -587,6 +595,7 @@ subroutine QE_geom_read(UNIT,basis,length) do i = 1, basis%nspec basis%spec(i)%num = 0 allocate(basis%spec(i)%atom_idx(tmp_natom(i))) + allocate(basis%spec(i)%atom_mask(tmp_natom(i)), source = .true.) allocate(basis%spec(i)%atom(tmp_natom(i),length_)) end do @@ -1302,7 +1311,7 @@ function convert_lat_to_abc(lattice, radians) result(abc_angle) do i = 1, 3 - abc_angle(1,i)=modu(lattice(i,:)) + abc_angle(1,i)=norm2(lattice(i,:)) end do do i = 1, 3 end do @@ -1363,8 +1372,8 @@ subroutine copy(this, basis, length) ! Local variables - integer :: i - !! Loop index. + integer :: i, j + !! Loop indices. integer :: length_, length_input !! The dimension of the basis atom positions. @@ -1385,6 +1394,8 @@ subroutine copy(this, basis, length) !--------------------------------------------------------------------------- if(allocated(this%spec))then do i = 1, this%nspec + if(allocated(this%spec(i)%atom_mask)) & + deallocate(this%spec(i)%atom_mask) if(allocated(this%spec(i)%atom_idx)) deallocate(this%spec(i)%atom_idx) if(allocated(this%spec(i)%atom)) deallocate(this%spec(i)%atom) end do @@ -1397,11 +1408,19 @@ subroutine copy(this, basis, length) !--------------------------------------------------------------------------- allocate(this%spec(basis%nspec)) do i = 1, basis%nspec + allocate(this%spec(i)%atom_mask(basis%spec(i)%num), source = .true.) allocate(this%spec(i)%atom_idx(basis%spec(i)%num)) - allocate(this%spec(i)%atom(& - basis%spec(i)%num,length_)) + allocate(this%spec(i)%atom(basis%spec(i)%num,length_)) + + if(allocated(basis%spec(i)%atom_mask)) & + this%spec(i)%atom_mask = basis%spec(i)%atom_mask - this%spec(i)%atom_idx = basis%spec(i)%atom_idx + if(allocated(basis%spec(i)%atom_idx))then + this%spec(i)%atom_idx = basis%spec(i)%atom_idx + else + this%spec(i)%atom_idx = [ ( j, j = sum(basis%spec(1:i-1:1)%num) + 1, & + sum(basis%spec(1:i)%num) ) ] + end if if(length_input.eq.length_)then this%spec(i)%atom(:,:length_) = basis%spec(i)%atom(:,:length_) elseif(length_input.gt.length_)then @@ -1429,6 +1448,134 @@ end subroutine copy !############################################################################### +!############################################################################### + subroutine set_atom_mask(this, index_list) + !! Set the mask for the atoms in the basis. + implicit none + + ! Arguments + class(basis_type), intent(inout) :: this + !! Parent. The basis. + integer, dimension(:,:), intent(in), optional :: index_list + !! The list of indices to set the mask for. + + ! Local variables + integer :: i + !! Loop index. + + + do i = 1, this%nspec + if(.not.allocated(this%spec(i)%atom_mask))then + allocate( & + this%spec(i)%atom_mask(this%spec(i)%num), source = .true. & + ) + end if + end do + + if(present(index_list))then + do i = 1, size(index_list,2) + this%spec(index_list(1,i))%atom_mask(index_list(2,i)) = & + .not.this%spec(index_list(1,i))%atom_mask(index_list(2,i)) + end do + end if + + end subroutine set_atom_mask +!############################################################################### + + +!############################################################################### + subroutine add_atom(this, species, position, is_cartesian, mask) + !! Add an atom to the basis. + implicit none + + ! Arguments + class(basis_type), intent(inout) :: this + !! Parent. The basis. + character(len=3), intent(in) :: species + !! The species of the atom to add. + real(real32), dimension(3), intent(in) :: position + !! The position of the atom to add. + logical, intent(in), optional :: is_cartesian + !! Optional. Boolean whether the position is in cartesian coordinates. + !! NOT YET IMPLEMENTED. + logical, intent(in), optional :: mask + !! Optional. Boolean whether to add a mask for the atom. + + ! Local variables + integer :: j + !! Loop index. + integer :: idx + !! The index of the species in the basis. + integer :: length + !! The dimension of the basis atom positions. + logical :: mask_ + !! Boolean mask for the atom. + integer, dimension(:), allocatable :: atom_idx + !! Temporary array. + logical, dimension(:), allocatable :: atom_mask + !! Temporary array. + real(real32), dimension(:,:), allocatable :: positions + !! Temporary array. + type(species_type), dimension(:), allocatable :: species_list + !! Temporary array. + + + mask_ = .true. + if(present(mask)) mask_ = mask + + this%natom = this%natom + 1 + length = size(this%spec(1)%atom,dim=2) + idx = findloc(this%spec(:)%name, strip_null(species), dim=1) + if(idx.eq.0)then + this%nspec = this%nspec + 1 + allocate(species_list(this%nspec)) + species_list(1:this%nspec-1) = this%spec(1:this%nspec-1) + deallocate(this%spec) + species_list(this%nspec)%name = strip_null(species) + species_list(this%nspec)%num = 1 + call get_element_properties(species_list(this%nspec)%name, & + species_list(this%nspec)%mass, & + species_list(this%nspec)%charge, & + species_list(this%nspec)%radius & + ) + allocate(species_list(this%nspec)%atom_idx(1)) + allocate(species_list(this%nspec)%atom_mask(1), source = mask_) + species_list(this%nspec)%atom_idx(1) = this%natom + allocate(species_list(this%nspec)%atom(1,length)) + species_list(this%nspec)%atom(1,:) = 0._real32 + species_list(this%nspec)%atom(1,:3) = position + this%spec = species_list + deallocate(species_list) + else + allocate(atom_mask(this%spec(idx)%num+1), source = .true.) + if(allocated(this%spec(idx)%atom_mask))then + atom_mask(1:this%spec(idx)%num) = this%spec(idx)%atom_mask + end if + atom_mask(this%spec(idx)%num+1) = mask_ + allocate(atom_idx(this%spec(idx)%num+1)) + if(allocated(this%spec(idx)%atom_idx))then + atom_idx(1:this%spec(idx)%num) = this%spec(idx)%atom_idx + else + atom_idx(1:this%spec(idx)%num) = [ ( j, j = 1, this%spec(idx)%num ) ] + end if + atom_idx(this%spec(idx)%num+1) = this%natom + allocate(positions(this%spec(idx)%num+1,length)) + positions = 0._real32 + positions(1:this%spec(idx)%num,:) = this%spec(idx)%atom + positions(this%spec(idx)%num+1,:3) = position + this%spec(idx)%num = this%spec(idx)%num + 1 + this%spec(idx)%atom_mask = atom_mask + this%spec(idx)%atom_idx = atom_idx + this%spec(idx)%atom = positions + deallocate(atom_mask) + deallocate(atom_idx) + deallocate(positions) + end if + + end subroutine add_atom +!############################################################################### + + !############################################################################### subroutine remove_atom(this, ispec, iatom) !! Remove an atom from the basis. @@ -1447,6 +1594,8 @@ subroutine remove_atom(this, ispec, iatom) !! The index associated with the atom to remove. integer, dimension(:), allocatable :: atom_idx !! Temporary array to store the atomic indices. + logical, dimension(:), allocatable :: atom_mask + !! Temporary array to store the atomic masks. real(real32), dimension(:,:), allocatable :: atom !! Temporary array to store the atomic positions. @@ -1461,19 +1610,25 @@ subroutine remove_atom(this, ispec, iatom) call stop_program("Atom to remove does not exist") return end if + allocate(atom_mask(this%spec(i)%num-1), source = .true.) allocate(atom_idx(this%spec(i)%num-1)) allocate(atom(this%spec(i)%num-1,size(this%spec(i)%atom,2))) if(iatom.eq.1)then + atom_mask(1:this%spec(i)%num-1) = & + this%spec(i)%atom_mask(2:this%spec(i)%num:1) atom_idx(1:this%spec(i)%num-1) = & this%spec(i)%atom_idx(2:this%spec(i)%num:1) atom(1:this%spec(i)%num-1:1,:) = & this%spec(i)%atom(2:this%spec(i)%num:1,:) elseif(iatom.eq.this%spec(i)%num)then + atom_mask(1:this%spec(i)%num-1) = & + this%spec(i)%atom_mask(1:this%spec(i)%num-1:1) atom_idx(1:this%spec(i)%num-1) = & this%spec(i)%atom_idx(1:this%spec(i)%num-1:1) atom(1:this%spec(i)%num-1:1,:) = & this%spec(i)%atom(1:this%spec(i)%num-1:1,:) else + atom_mask(1:iatom-1:1) = this%spec(i)%atom_mask(1:iatom-1:1) atom_idx(1:iatom-1:1) = this%spec(i)%atom_idx(1:iatom-1:1) atom_idx(iatom:this%spec(i)%num-1:1) = & this%spec(i)%atom_idx(iatom+1:this%spec(i)%num:1) @@ -1482,12 +1637,14 @@ subroutine remove_atom(this, ispec, iatom) this%spec(i)%atom(iatom+1:this%spec(i)%num:1,:) end if where(atom_idx(1:this%spec(i)%num-1:1).gt.remove_idx) - atom_idx(1:this%spec(i)%num-1:1) = & - atom_idx(1:this%spec(i)%num-1:1) - 1 + atom_idx(1:this%spec(i)%num-1:1) = & + atom_idx(1:this%spec(i)%num-1:1) - 1 end where + this%spec(i)%atom_mask = atom_mask this%spec(i)%atom_idx = atom_idx this%spec(i)%atom = atom - deallocate(this%spec(i)%atom_idx) + deallocate(atom_mask) + deallocate(atom_idx) deallocate(atom) this%spec(i)%num = this%spec(i)%num - 1 this%natom = this%natom - 1 diff --git a/src/fortran/lib/mod_geom_utils.f90 b/src/fortran/lib/mod_geom_utils.f90 index a3d6d6f..797a5e4 100644 --- a/src/fortran/lib/mod_geom_utils.f90 +++ b/src/fortran/lib/mod_geom_utils.f90 @@ -17,7 +17,8 @@ module raffle__geom_utils contains !############################################################################### - function basis_merge(basis1,basis2,length,map1,map2) result(output) + function basis_merge(basis1,basis2,length,map1,map2, mask1, mask2) & + result(output) !! Merge two supplied bases !! !! Merge two bases assuming that the lattice is the same @@ -32,6 +33,8 @@ function basis_merge(basis1,basis2,length,map1,map2) result(output) !! Number of dimensions for atomic positions (default 3). integer, allocatable, dimension(:,:,:), optional, intent(inout) :: map1,map2 !! Maps for atoms in the two bases. + logical, intent(in), optional :: mask1, mask2 + !! Mask for atoms in the two bases. ! Local variables integer :: i, j, k, itmp, dim @@ -100,31 +103,61 @@ function basis_merge(basis1,basis2,length,map1,map2) result(output) ! set up atoms in merged basis !--------------------------------------------------------------------------- do i = 1, basis1%nspec + allocate(output%spec(i)%atom_mask(output%spec(i)%num), source = .true.) allocate(output%spec(i)%atom_idx(output%spec(i)%num)) allocate(output%spec(i)%atom(output%spec(i)%num,dim)) + if(allocated(basis1%spec(i)%atom_mask)) & + output%spec(i)%atom_mask(1:basis1%spec(i)%num) = basis1%spec(i)%atom_mask + if(allocated(basis1%spec(i)%atom_idx))then + output%spec(i)%atom_idx(1:basis1%spec(i)%num) = basis1%spec(i)%atom_idx + else + output%spec(i)%atom_idx(1:basis1%spec(i)%num) = [(i,i=1,basis1%spec(i)%num)] + end if output%spec(i)%atom(:,:)=0._real32 output%spec(i)%atom(1:basis1%spec(i)%num,:3) = basis1%spec(i)%atom(:,:3) - output%spec(i)%atom_idx(1:basis1%spec(i)%num) = basis1%spec(i)%atom_idx if(lmap) new_map(i,:basis1%spec(i)%num,:) = map1(i,:basis1%spec(i)%num,:) + if(present(mask1)) output%spec(i)%atom_mask(1:basis1%spec(i)%num) = mask1 end do do i = 1, basis2%nspec if(match(i).gt.basis1%nspec)then + allocate(output%spec(match(i))%atom_mask(output%spec(match(i))%num), & + source = .true.) + if(allocated(basis2%spec(i)%atom_mask)) & + output%spec(match(i))%atom_mask(:) = basis2%spec(i)%atom_mask(:) allocate(output%spec(match(i))%atom_idx(output%spec(match(i))%num)) - output%spec(match(i))%atom_idx(:) = & - basis2%spec(i)%atom_idx(:) + basis1%natom + if(allocated(basis2%spec(i)%atom_idx))then + output%spec(match(i))%atom_idx(:) = & + basis2%spec(i)%atom_idx(:) + basis1%natom + else + output%spec(match(i))%atom_idx(:) = & + [(i,i=1,basis2%spec(i)%num)] + end if allocate(output%spec(match(i))%atom(output%spec(match(i))%num,dim)) output%spec(match(i))%atom(:,:)=0._real32 output%spec(match(i))%atom(:,:3)=basis2%spec(i)%atom(:,:3) if(lmap) new_map(match(i),:basis2%spec(i)%num,:) = & map2(i,:basis2%spec(i)%num,:) + if(present(mask2)) output%spec(match(i))%atom_mask(:) = mask2 else itmp=basis1%spec(match(i))%num - output%spec(match(i))%atom_idx(itmp+1:basis2%spec(i)%num+itmp) = & - basis2%spec(i)%atom_idx(:) + basis1%natom + if(allocated(basis2%spec(i)%atom_mask)) & + output%spec(match(i))%atom_mask(itmp+1:basis2%spec(i)%num+itmp) = & + basis2%spec(i)%atom_mask(:) + if(allocated(basis2%spec(i)%atom_idx))then + output%spec(match(i))%atom_idx(itmp+1:basis2%spec(i)%num+itmp) = & + basis2%spec(i)%atom_idx(:) + basis1%natom + else + output%spec(match(i))%atom_idx(itmp+1:basis2%spec(i)%num+itmp) = & + [(i,i=1,basis2%spec(i)%num)] + end if output%spec(match(i))%atom(itmp+1:basis2%spec(i)%num+itmp,:3) = & - basis2%spec(i)%atom(:,:3) + basis2%spec(i)%atom(:,:3) if(lmap) new_map(match(i),itmp+1:basis2%spec(i)%num+itmp,:) = & - map2(i,:basis2%spec(i)%num,:) + map2(i,:basis2%spec(i)%num,:) + if(present(mask2)) & + output%spec(match(i))%atom_mask( & + itmp+1:basis2%spec(i)%num+itmp & + ) = mask2 end if end do output%natom=sum(output%spec(:)%num) diff --git a/src/fortran/lib/mod_misc_linalg.f90 b/src/fortran/lib/mod_misc_linalg.f90 index b78249d..00f2713 100644 --- a/src/fortran/lib/mod_misc_linalg.f90 +++ b/src/fortran/lib/mod_misc_linalg.f90 @@ -38,7 +38,7 @@ pure function modu(vector) real(real32) :: modu !! Output magnitude. - modu = abs( sqrt( sum(vector(:)**2) ) ) + modu = sqrt( sum( vector ** 2 ) ) end function modu !############################################################################### @@ -73,7 +73,7 @@ pure function get_distance(point1, point2) result(distance) real(real32) :: distance !! Output distance. - distance = modu( point1 - point2 ) + distance = norm2( point1 - point2 ) return end function get_distance @@ -92,7 +92,7 @@ pure function get_angle_from_vectors(vector1, vector2) result(angle) !! Output angle. angle = dot_product(vector1,vector2) / & - ( modu(vector1) * modu(vector2) ) + ( norm2(vector1) * norm2(vector2) ) if(angle .ge. 1._real32)then angle = 0._real32 elseif(angle .le. -1._real32)then @@ -118,7 +118,7 @@ pure function get_angle_from_points(point1, point2, point3) result(angle) !! Output angle. angle = dot_product( point1 - point2, point3 - point2 ) / & - ( modu( point1 - point2 ) * modu( point3 - point2 ) ) + ( norm2( point1 - point2 ) * norm2( point3 - point2 ) ) if(angle .ge. 1._real32)then angle = 0._real32 elseif(angle .le. -1._real32)then @@ -164,9 +164,9 @@ pure function get_dihedral_angle_from_points( & implicit none real(real32), dimension(3), intent(in) :: point1, point2, point3, point4 real(real32) :: angle - + angle = get_angle(cross(point2 - point1, point3 - point2), point4 - point2) - + end function get_dihedral_angle_from_points !############################################################################### @@ -204,7 +204,7 @@ pure function get_improper_dihedral_angle_from_points( & !! Return the improper dihedral angle between two planes. !! !! The dihedral angle is the angle between the plane defined by four points. - !! i.e. ( point2 - point1 ) x ( point3 - point1 ) . + !! i.e. ( point2 - point1 ) x ( point3 - point1 ) . !! ( point4 - point2 ) x ( point3 - point1 ) !! alt. angle between plane point1point2point3 and point1point3point4 implicit none diff --git a/src/fortran/lib/mod_misc_maths.f90 b/src/fortran/lib/mod_misc_maths.f90 index 594ea2f..b410e36 100644 --- a/src/fortran/lib/mod_misc_maths.f90 +++ b/src/fortran/lib/mod_misc_maths.f90 @@ -39,14 +39,18 @@ end function lnsum !############################################################################### - integer function triangular_number(n) + pure function triangular_number(n) result(output) !! Return the nth triangular number. implicit none ! Arguments - integer :: n + integer, intent(in) :: n !! The index of the triangular number to return. - triangular_number = n * ( n + 1 ) / 2 + + real(real32) :: output + !! The nth triangular number. + + output = n * ( n + 1 ) / 2 end function triangular_number !############################################################################### diff --git a/src/fortran/lib/mod_place_methods.f90 b/src/fortran/lib/mod_place_methods.f90 index d2e262a..7bc1b02 100644 --- a/src/fortran/lib/mod_place_methods.f90 +++ b/src/fortran/lib/mod_place_methods.f90 @@ -10,7 +10,7 @@ module raffle__place_methods !! as the starting point !! - min: place the atom at the gridpoint with the highest viability use raffle__constants, only: real32, pi - use raffle__misc_linalg, only: modu, inverse_3x3 + use raffle__misc_linalg, only: inverse_3x3 use raffle__geom_extd, only: extended_basis_type use raffle__dist_calcs, only: & get_min_dist, & @@ -32,7 +32,7 @@ module raffle__place_methods !############################################################################### function place_method_void( & - points, basis, atom_ignore_list, viable & + points, basis, viable & ) result(point) !! VOID placement method. !! @@ -45,8 +45,6 @@ function place_method_void( & !! List of gridpoints to consider. type(extended_basis_type), intent(inout) :: basis !! Structure to add atom to. - integer, dimension(:,:), intent(in) :: atom_ignore_list - !! List of atoms to ignore (i.e. indices of atoms not yet placed). logical, intent(out) :: viable !! Boolean to indicate if point is viable. real(real32), dimension(3) :: point @@ -79,7 +77,7 @@ end function place_method_void !############################################################################### function place_method_rand( distribs_container, & bounds, & - basis, atom_ignore_list, radius_list, max_attempts, viable & + basis, species, radius_list, max_attempts, viable & ) result(point) !! Random placement method. !! @@ -93,8 +91,8 @@ function place_method_rand( distribs_container, & !! Bounds of the unit cell. type(extended_basis_type), intent(inout) :: basis !! Structure to add atom to. - integer, dimension(:,:), intent(in) :: atom_ignore_list - !! List of atoms to ignore (i.e. indices of atoms not yet placed). + integer, intent(in) :: species + !! Species index to add atom to. real(real32), dimension(:), intent(in) :: radius_list !! List of radii for each pair of elements. integer, intent(in) :: max_attempts @@ -140,11 +138,10 @@ function place_method_rand( distribs_container, & if( & get_min_dist_between_point_and_species( & basis, point, & - species = js, & - ignore_list = atom_ignore_list & + species = js & ) .lt. max( & distribs_container%cutoff_min(1), & - radius_list(pair_index(atom_ignore_list(1,1),js)) * & + radius_list(pair_index(species,js)) * & distribs_container%radius_distance_tol(1) & ) & )then @@ -162,7 +159,7 @@ end function place_method_rand !############################################################################### function place_method_walk( distribs_container, & bounds, & - basis, atom_ignore_list, & + basis, species, & radius_list, max_attempts, & step_size_coarse, step_size_fine, & viable & @@ -185,14 +182,14 @@ function place_method_walk( distribs_container, & !! Bounds of the unit cell. type(extended_basis_type), intent(inout) :: basis !! Structure to add atom to. + integer, intent(in) :: species + !! Species index to add atom to. integer, intent(in) :: max_attempts !! Limit on number of attempts. real(real32), intent(in) :: step_size_coarse, step_size_fine !! Step sizes for random walk. logical, intent(out) :: viable !! Boolean to indicate if point is viable. - integer, dimension(:,:), intent(in) :: atom_ignore_list - !! List of atoms to ignore (i.e. indices of atoms not yet placed). real(real32), dimension(:), intent(in) :: radius_list !! List of radii for each pair of elements. real(real32), dimension(3) :: point @@ -221,7 +218,7 @@ function place_method_walk( distribs_container, & ! test a random point in the unit cell !--------------------------------------------------------------------------- do i = 1, 3 - abc(i) = modu(basis%lat(i,:)) + abc(i) = norm2(basis%lat(i,:)) end do i = 0 random_loop : do @@ -231,8 +228,7 @@ function place_method_walk( distribs_container, & site_vector = bounds(1,:) + ( bounds(2,:) - bounds(1,:) ) * site_vector site_value = evaluate_point( distribs_container, & - site_vector, atom_ignore_list(1,1), basis, & - atom_ignore_list, radius_list & + site_vector, species, basis, radius_list & ) call random_number(rtmp1) if(rtmp1.lt.site_value) exit random_loop @@ -272,8 +268,7 @@ function place_method_walk( distribs_container, & ! evaluate the test point !------------------------------------------------------------------------ test_value = evaluate_point( distribs_container, & - test_vector, atom_ignore_list(1,1), basis, & - atom_ignore_list, radius_list & + test_vector, species, basis, radius_list & ) if(test_value.lt.1.E-6) cycle walk_loop !------------------------------------------------------------------------ @@ -318,7 +313,7 @@ end function place_method_walk function place_method_growth( distribs_container, & prior_point, prior_species, & bounds, & - basis, atom_ignore_list, & + basis, species, & radius_list, max_attempts, & step_size_coarse, step_size_fine, & viable & @@ -345,14 +340,14 @@ function place_method_growth( distribs_container, & !! Bounds of the unit cell. type(extended_basis_type), intent(inout) :: basis !! Structure to add atom to. + integer, intent(in) :: species + !! Species index to add atom to. integer, intent(in) :: max_attempts !! Limit on number of attempts. real(real32), intent(in) :: step_size_coarse, step_size_fine !! Step sizes for random walk. logical, intent(out) :: viable !! Boolean to indicate if point is viable. - integer, dimension(:,:), intent(in) :: atom_ignore_list - !! List of atoms to ignore (i.e. indices of atoms not yet placed). real(real32), dimension(:), intent(in) :: radius_list !! List of radii for each pair of elements. real(real32), dimension(3) :: point @@ -382,7 +377,7 @@ function place_method_growth( distribs_container, & ! get the lattice constants and the inverse lattice !--------------------------------------------------------------------------- do i = 1, 3 - abc(i) = modu(basis%lat(i,:)) + abc(i) = norm2(basis%lat(i,:)) end do inverse_lattice = inverse_3x3(basis%lat) @@ -392,7 +387,7 @@ function place_method_growth( distribs_container, & !--------------------------------------------------------------------------- idx = distribs_container%get_pair_index( & basis%spec(prior_species)%name, & - basis%spec(atom_ignore_list(1,1))%name & + basis%spec(species)%name & ) min_radius = & max( & @@ -428,8 +423,7 @@ function place_method_growth( distribs_container, & end do ! now evaluate the point and check if it passes the initial criteria site_value = evaluate_point( distribs_container, & - site_vector, atom_ignore_list(1,1), basis, & - atom_ignore_list, radius_list & + site_vector, species, basis, radius_list & ) call random_number(rtmp1) if(rtmp1.lt.site_value) exit shell_loop @@ -469,8 +463,7 @@ function place_method_growth( distribs_container, & ! evaluate the test point !------------------------------------------------------------------------ test_value = evaluate_point( distribs_container, & - test_vector, atom_ignore_list(1,1), basis, & - atom_ignore_list, radius_list & + test_vector, species, basis, radius_list & ) if(test_value.lt.1.E-6) cycle walk_loop !------------------------------------------------------------------------ diff --git a/src/fortran/lib/mod_viability.f90 b/src/fortran/lib/mod_viability.f90 index d18af41..f67019f 100644 --- a/src/fortran/lib/mod_viability.f90 +++ b/src/fortran/lib/mod_viability.f90 @@ -7,7 +7,7 @@ module raffle__viability use omp_lib #endif use raffle__constants, only: real32 - use raffle__misc_linalg, only: modu, inverse_3x3 + use raffle__misc_linalg, only: inverse_3x3 use raffle__geom_extd, only: extended_basis_type use raffle__dist_calcs, only: & get_min_dist_between_point_and_atom, get_min_dist @@ -27,7 +27,7 @@ module raffle__viability function get_gridpoints_and_viability(distribs_container, grid, bounds, & basis, & species_index_list, & - radius_list, atom_ignore_list, grid_offset) result(points) + radius_list, grid_offset) result(points) !! Return a list of viable gridpoints and their viability for each species. !! !! This function returns the viability of all viable gridpoints. @@ -46,8 +46,6 @@ function get_gridpoints_and_viability(distribs_container, grid, bounds, & !! List of radii for each pair of elements. integer, dimension(:), intent(in) :: species_index_list !! List of species indices to add atoms to. - integer, dimension(:,:), intent(in) :: atom_ignore_list - !! List of atoms to ignore (i.e. indices of atoms not yet placed). real(real32), dimension(3), intent(in) :: grid_offset !! Offset for gridpoints. real(real32), dimension(:,:), allocatable :: points @@ -126,13 +124,11 @@ function get_gridpoints_and_viability(distribs_container, grid, bounds, & !$omp parallel do default(shared) private(i,is,ia,l,atom_idx,idx) do is = 1, basis%nspec atom_loop: do ia = 1, basis%spec(is)%num - do l = 1, size(atom_ignore_list,dim=2), 1 - if(all(atom_ignore_list(:,l).eq.[is,ia])) cycle atom_loop - end do + if(.not. basis%spec(is)%atom_mask(ia)) cycle atom_loop ! get the atom position in terms of the grid indices atom_idx = & - nint( ( basis%spec(is)%atom(ia,1:3) - offset )/ grid_scale ) + nint( ( basis%spec(is)%atom(ia,1:3) - offset ) / grid_scale ) ! if any one of the indicies is always outside of the grid, skip idx_lw = atom_idx - extent @@ -184,14 +180,10 @@ function get_gridpoints_and_viability(distribs_container, grid, bounds, & do i = 1, num_points do concurrent ( is = 1 : size(species_index_list,1) ) points(4,i) = & - modu( get_min_dist(& - basis, [ points(1:3,i) ], .false., & - ignore_list = atom_ignore_list & - ) ) + norm2( get_min_dist(basis, [ points(1:3,i) ], .false. ) ) points(4+is,i) = & evaluate_point( distribs_container, & - points(1:3,i), species_index_list(is), basis, & - atom_ignore_list, radius_list & + points(1:3,i), species_index_list(is), basis, radius_list & ) end do end do @@ -205,7 +197,7 @@ end function get_gridpoints_and_viability subroutine update_gridpoints_and_viability( & points, distribs_container, basis, & species_index_list, & - atom, radius_list, atom_ignore_list & + atom, radius_list & ) !! Update the list of viable gridpoints and their viability for each !! species. @@ -226,14 +218,14 @@ subroutine update_gridpoints_and_viability( & !! List of radii for each pair of elements. integer, dimension(:), intent(in) :: species_index_list !! List of species indices to add atoms to. - integer, dimension(:,:), intent(in) :: atom_ignore_list - !! List of atoms to ignore (i.e. indices of atoms not yet placed). ! Local variables integer :: i, is !! Loop indices. integer :: num_points !! Number of gridpoints. + integer :: num_species + !! Number of species. real(real32) :: min_radius !! Minimum radius. real(real32) :: distance @@ -244,6 +236,8 @@ subroutine update_gridpoints_and_viability( & !! Temporary list of gridpoints. real(real32), dimension(3) :: diff !! Difference between atom and gridpoint (direct coorindates). + real(real32), dimension(3) :: atom_pos + !! Position of atom in direct coordinates. real(real32), dimension(:,:), allocatable :: points_tmp !! Temporary list of gridpoints. @@ -254,29 +248,33 @@ subroutine update_gridpoints_and_viability( & if(.not.allocated(points)) return num_points = size(points,dim=2) viable = .true. - min_radius = minval(radius_list) * distribs_container%radius_distance_tol(1) - associate( atom_pos => [ basis%spec(atom(1))%atom(atom(2),1:3) ] ) + min_radius = max( & + distribs_container%cutoff_min(1), & + minval(radius_list) * distribs_container%radius_distance_tol(1) & + ) + atom_pos = basis%spec(atom(1))%atom(atom(2),1:3) + num_species = size(species_index_list,1) !$omp parallel do default(shared) private(i,is,diff,distance) - do i = 1, num_points - diff = atom_pos - points(1:3,i) - diff = diff - ceiling(diff - 0.5_real32) - distance = modu( matmul( diff, basis%lat ) ) - if( distance .lt. min_radius )then - viable(i) = .false. - cycle - elseif( distance .le. distribs_container%cutoff_max(1) )then - do concurrent( is = 1 : size(species_index_list,1) ) + do i = 1, num_points + diff = atom_pos - points(1:3,i) + diff = diff - anint(diff) + distance = norm2( matmul( diff, basis%lat ) ) + if( distance .lt. min_radius )then + viable(i) = .false. + else + if( distance .le. distribs_container%cutoff_max(1) )then + do concurrent( is = 1 : num_species ) points(4+is,i) = & evaluate_point( distribs_container, & points(1:3,i), species_index_list(is), basis, & - atom_ignore_list, radius_list & + radius_list & ) end do end if points(4,i) = min( points(4,i), distance ) - end do + end if + end do !$omp end parallel do - end associate num_points = count(viable) if(num_points.lt.1)then diff --git a/src/fortran/raffle.f90 b/src/fortran/raffle.f90 index d0c6f8a..42a7271 100644 --- a/src/fortran/raffle.f90 +++ b/src/fortran/raffle.f90 @@ -3,6 +3,8 @@ module raffle use raffle__io_utils, only: raffle__version__ use raffle__generator, only: raffle_generator_type use raffle__distribs_container, only: distribs_container_type + use raffle__cache, only: & + store_probability_density, retrieve_probability_density implicit none @@ -12,4 +14,4 @@ module raffle public :: raffle_generator_type -end module raffle \ No newline at end of file +end module raffle diff --git a/src/raffle/raffle.py b/src/raffle/raffle.py index ebb6732..9deb45e 100644 --- a/src/raffle/raffle.py +++ b/src/raffle/raffle.py @@ -31,9 +31,10 @@ def __init__(self, handle=None): """ Create a ``species_type`` object. - Returns: - species (species_type): - Object to be constructed + Returns + ------- + Geom_Rw.species_type + Object to be constructed """ f90wrap.runtime.FortranDerivedType.__init__(self) result = _raffle.f90wrap_geom_rw__species_type_initialise() @@ -41,17 +42,7 @@ def __init__(self, handle=None): def __del__(self): """ - Destructor for class species_type - - - Defined at ../src/lib/mod_geom_rw.f90 lines \ - 26-32 - - Parameters - ---------- - this : species_type - Object to be destructed - + Destructor for ``species_type`` object. Automatically generated destructor for species_type """ @@ -269,9 +260,10 @@ def __init__(self, atoms=None, handle=None): including lattice and basis information. This is confusingly named as a crystal = lattice + basis. - Returns: - basis (basis): - Object to be constructed + Returns + ------- + Geom_Rw.basis + Object to be constructed """ f90wrap.runtime.FortranDerivedType.__init__(self) result = _raffle.f90wrap_geom_rw__basis_type_initialise() @@ -282,17 +274,7 @@ def __init__(self, atoms=None, handle=None): def __del__(self): """ - Destructor for class basis - - - Defined at ../src/lib/mod_geom_rw.f90 lines \ - 34-42 - - Parameters - ---------- - this : basis - Object to be destructed - + Destructor for ``basis`` object. Automatically generated destructor for basis """ @@ -304,17 +286,18 @@ def allocate_species(self, num_species=None, species_symbols=None, species_count """ Allocate memory for the species list. - Parameters: - num_species (int): - Number of species - species_symbols (list of str): - List of species symbols - species_count (list of int): - List of species counts - atoms (list of float): - List of atomic positions - atom_idx_list (list of int): - List of atomic indices + Parameters + ---------- + num_species : int + Number of species + species_symbols : list[str] + List of species symbols + species_count : list[int] + List of species counts + atoms : list[float] + List of atomic positions + atom_idx_list : list[int] + List of atomic indices """ _raffle.f90wrap_geom_rw__allocate_species__binding__basis_type(this=self._handle, \ num_species=num_species, species_symbols=species_symbols, species_count=species_count, \ @@ -341,9 +324,10 @@ def copy(self): """ Create a copy of the basis object. - Returns: - basis (basis): - Copy of the basis object + Returns + ------- + Geom_Rw.basis + Copy of the basis object """ # Create a new basis object @@ -376,13 +360,14 @@ def toase(self, calculator=None, return_none_on_failure : bool = True): """ Convert the basis object to an ASE Atoms object. - Parameters: - calculator (ASE Calculator): - ASE calculator object to be assigned to the Atoms object. - return_none_on_failure (bool): - Boolean whether to return None on failure. - If True, return None instead of raising an exception. - If False, raise an exception on failure. + Parameters + ---------- + calculator : ase.calculators.calculator.BaseCalculator + ASE calculator object to be assigned to the Atoms object. + return_none_on_failure: bool + Boolean whether to return None on failure. + If True, return None instead of raising an exception. + If False, raise an exception on failure. """ # Set the species list @@ -521,11 +506,12 @@ def fromase(self, atoms, verbose=False): """ Convert the ASE Atoms object to a basis object. - Parameters: - atoms (ASE Atoms): - ASE Atoms object to be converted. - verbose (bool): - Boolean whether to print warnings. + Parameters + ---------- + atoms : ase.Atoms + ASE Atoms object to be converted. + verbose : bool + Boolean whether to print warnings. """ from ase.calculators.singlepoint import SinglePointCalculator @@ -710,10 +696,10 @@ def __init__(self, atoms=None, handle=None): """ Create a ``basis_array`` object. - - Returns: - basis_array (basis_array): - Object to be constructed + Returns + ------- + Geom_Rw.basis_array + Object to be constructed """ f90wrap.runtime.FortranDerivedType.__init__(self) @@ -735,16 +721,6 @@ def __del__(self): """ Destructor for class basis_array - - Defined at ../src/lib/mod_generator.f90 lines \ - 19-21 - - Parameters - ---------- - this : basis_array - Object to be destructed - - Automatically generated destructor for basis_array """ if self._alloc: @@ -782,9 +758,10 @@ def allocate(self, size): """ Allocate the items array with the given size. - Parameters: - size (int): - Size of the items array + Parameters + ---------- + size : int + Size of the items array """ _raffle.f90wrap_basis_type_xnum_array__array_alloc__items(self._handle, num=size) @@ -821,9 +798,10 @@ def __init__(self, handle=None): """ Create a ``Distribs_Container_Type`` object. - Returns: - distribution_container (Distribs_Container_Type): - Object to be constructed + Returns + ------- + Distribs_Container_Type + Object to be constructed """ f90wrap.runtime.FortranDerivedType.__init__(self) @@ -833,45 +811,37 @@ def __init__(self, handle=None): def __del__(self): """ - Destructor for class Distribs_Container_Type - - - Defined at ../fortran/lib/mod_distribs_container.f90 \ - lines 25-162 - - Parameters - ---------- - this : Distribs_Container_Type - Object to be destructed - + Destructor for ``Distribs_Container_Type`` object. Automatically generated destructor for distribs_container_type """ if self._alloc: _raffle.f90wrap_raffle__dc__dc_type_finalise(this=self._handle) - def set_kBT(self, kBT): + def set_kBT(self, kBT : float): """ Set the energy scale for the distribution functions. - Parameters: - kBT (float): - Energy scale for the distribution functions. + Parameters + ---------- + kBT : float + Energy scale for the distribution functions. """ self.kBT = kBT - def set_weight_method(self, method): + def set_weight_method(self, method : str): """ Set the weight method for combining the the distribution functions to form the generalised distribution function (GDF). - Parameters: - method (str): - Method to be used for weighting the distribution functions. - Allowed values are: - - 'formation_energy' or 'formation' or 'form' or 'e_form' - - 'energy_above_hull' or 'hull_distance' or 'hull' or 'distance' or 'convex_hull' + Parameters + ---------- + method : str + Method to be used for weighting the distribution functions. + Allowed values are: + - 'formation_energy' or 'formation' or 'form' or 'e_form' + - 'energy_above_hull' or 'hull_distance' or 'hull' or 'distance' or 'convex_hull' """ @@ -882,78 +852,98 @@ def set_weight_method(self, method): else: raise ValueError("Invalid weight method: {}".format(method)) - def set_width(self, width): + def set_width(self, width : list[float]): """ Set the distribution function widths. - Parameters: - width (list of float): - List of distribution function widths. - The first element is the 2-body distribution function width, - the second element is the 3-body distribution function width, - and the third element is the 4-body distribution function width. + Parameters + ---------- + width : list[float] + List of distribution function widths. + The first element is the 2-body distribution function width, + the second element is the 3-body distribution function width, + and the third element is the 4-body distribution function width. """ + if not len(width) == 3: + raise ValueError("width must be a list of exactly three floats") + _raffle.f90wrap_raffle__dc__set_width__binding__dc_type(this=self._handle, \ width=width) - def set_sigma(self, sigma): + def set_sigma(self, sigma : list[float]): """ Set the sigma values of the Gaussians used to build the distribution functions. - Parameters: - sigma (list of float): - List of sigma values. - The first element is the 2-body distribution function sigma, - the second element is the 3-body distribution function sigma, - and the third element is the 4-body distribution function sigma. + Parameters + ---------- + sigma : list[float] + List of sigma values. + The first element is the 2-body distribution function sigma, + the second element is the 3-body distribution function sigma, + and the third element is the 4-body distribution function sigma. """ + if not len(sigma) == 3: + raise ValueError("sigma must be a list of exactly three floats") + _raffle.f90wrap_raffle__dc__set_sigma__binding__dc_type(this=self._handle, \ sigma=sigma) - def set_cutoff_min(self, cutoff_min): + def set_cutoff_min(self, cutoff_min : list[float]): """ Set the minimum cutoff values for the distribution functions. - Parameters: - cutoff_min (list of float): - List of minimum cutoff values. - The first element is the 2-body distribution function minimum cutoff, - the second element is the 3-body distribution function minimum cutoff, - and the third element is the 4-body distribution function minimum cutoff. + Parameters + ---------- + cutoff_min : list[float] + List of minimum cutoff values. + The first element is the 2-body distribution function minimum cutoff, + the second element is the 3-body distribution function minimum cutoff, + and the third element is the 4-body distribution function minimum cutoff. """ + if not len(cutoff_min) == 3: + raise ValueError("cutoff_max must be a list of exactly three floats") + _raffle.f90wrap_raffle__dc__set_cutoff_min__binding__dc_type(this=self._handle, \ cutoff_min=cutoff_min) - def set_cutoff_max(self, cutoff_max): + def set_cutoff_max(self, cutoff_max : list[float]): """ Set the maximum cutoff values for the distribution functions. - Parameters: - cutoff_min (list of float): - List of maximum cutoff values. - The first element is the 2-body distribution function maximum cutoff, - the second element is the 3-body distribution function maximum cutoff, - and the third element is the 4-body distribution function maximum cutoff. + Parameters + ---------- + cutoff_min : list[float] + List of maximum cutoff values. + The first element is the 2-body distribution function maximum cutoff, + the second element is the 3-body distribution function maximum cutoff, + and the third element is the 4-body distribution function maximum cutoff. """ + if not len(cutoff_max) == 3: + raise ValueError("cutoff_max must be a list of exactly three floats") + _raffle.f90wrap_raffle__dc__set_cutoff_max__binding__dc_type(this=self._handle, \ cutoff_max=cutoff_max) - def set_radius_distance_tol(self, radius_distance_tol): + def set_radius_distance_tol(self, radius_distance_tol : list[float]): """ Set the radius distance tolerance for the distribution functions. The radius distance tolerance represents a multiplier to the bond radii to determine the cutoff distance for the distribution functions. - Parameters: - radius_distance_tol (list of float): - List of radius distance tolerance values. - The first two values are the lower and upper bounds for the - 3-body distribution function radius distance tolerance. - The third and fourth values are the lower and upper bounds for the - 4-body distribution function radius distance tolerance. + Parameters + ---------- + radius_distance_tol : list[float] + List of radius distance tolerance values. + The first two values are the lower and upper bounds for the + 3-body distribution function radius distance tolerance. + The third and fourth values are the lower and upper bounds for the + 4-body distribution function radius distance tolerance. """ + if not len(radius_distance_tol) == 4: + raise ValueError("radius_distance_tol must be a list of exactly four floats") + _raffle.f90wrap_raffle__dc__set_radius_distance_tol__binding__dc_type(this=self._handle, \ radius_distance_tol=radius_distance_tol) @@ -961,26 +951,33 @@ 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. + 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): + def create(self, + basis_list, + energy_above_hull_list : list[float] = None, + deallocate_systems : bool = True, + verbose : int = None + ): """ Create the distribution functions. - Parameters: - basis_list (basis_array or Atoms or list of Atoms): - List of atomic configurations to be used to create the distribution functions. - energy_above_hull_list (list of float): - List of energy above hull values for the atomic configurations. - deallocate_systems (bool): - Boolean whether to deallocate the atomic configurations after creating the distribution functions. - verbose (int): - Verbosity level. + Parameters + ---------- + basis_list : geom_rw.basis_array or ase.Atoms or list[ase.Atoms] + List of atomic configurations to be used to create the distribution functions. + energy_above_hull_list : list[float] + List of energy above hull values for the atomic configurations. + deallocate_systems : bool + Boolean whether to deallocate the atomic configurations after creating the distribution functions. + verbose : int + Verbosity level. """ if isinstance(basis_list, Atoms): basis_list = geom_rw.basis_array(basis_list) @@ -995,21 +992,28 @@ def create(self, basis_list, energy_above_hull_list=None, deallocate_systems=Tru verbose=verbose \ ) - def update(self, basis_list, energy_above_hull_list=None, from_host=True, deallocate_systems=True, verbose=None): + def update(self, + basis_list, + energy_above_hull_list : list[float] = None, + from_host : bool = True, + deallocate_systems : bool = True, + verbose : int = None + ): """ Update the distribution functions. - Parameters: - basis_list (basis_array or Atoms or list of Atoms): - List of atomic configurations to be used to create the distribution functions. - energy_above_hull_list (list of float): - List of energy above hull values for the atomic configurations. - deallocate_systems (bool): - Boolean whether to deallocate the atomic configurations after creating the distribution functions. - from_host (bool): - Boolean whether the provided basis_list is based on the host. - verbose (int): - Verbosity level. + Parameters + ---------- + basis_list : geom_rw.basis_array or ase.Atoms or list[ase.Atoms] + List of atomic configurations to be used to create the distribution functions. + energy_above_hull_list : list[float] + List of energy above hull values for the atomic configurations. + deallocate_systems : bool + Boolean whether to deallocate the atomic configurations after creating the distribution functions. + from_host : bool + Boolean whether the provided basis_list is based on the host. + verbose : int + Verbosity level. """ if isinstance(basis_list, Atoms): basis_list = geom_rw.basis_array(basis_list) @@ -1039,25 +1043,27 @@ def add_basis(self, basis): It is not recommended to use this function directly, but to use the create or update functions instead. - Parameters: - basis (geom_rw.basis): - Basis object to be added to the distribution functions. + Parameters + ---------- + basis : Geom_Rw.basis + Basis object to be added to the distribution functions. """ _raffle.f90wrap_raffle__dc__add_basis__binding__dc_type(this=self._handle, \ basis=basis._handle) - def set_element_energies(self, element_energies): + def set_element_energies(self, element_energies : dict): """ Set the element reference energies for the distribution functions. These energies are used to calculate the formation energies of the atomic configurations. - Parameters: - element_energies (dict): - Dictionary of element reference energies. - The keys are the element symbols and the values are the reference energies. + Parameters + ---------- + element_energies : dict + Dictionary of element reference energies. + The keys are the element symbols and the values are the reference energies. """ element_list = list(element_energies.keys()) @@ -1069,10 +1075,11 @@ def get_element_energies(self): """ Get the element reference energies for the distribution functions. - Returns: - element_energies (dict): - Dictionary of element reference energies. - The keys are the element symbols and the values are the reference energies + Returns + ------- + dict + Dictionary of element reference energies. + The keys are the element symbols and the values are the reference energies """ num_elements = _raffle.f90wrap_raffle__dc__get__num_elements(self._handle) @@ -1099,14 +1106,15 @@ def set_bond_info(self): """ _raffle.f90wrap_raffle__dc__set_bond_info__binding__dc_type(this=self._handle) - def set_bond_radius(self, radius_dict): + def set_bond_radius(self, radius_dict : dict): """ Set the bond radius for the distribution functions. - Parameters: - radius_dict (dict): - Dictionary of bond radii. - The keys are a tuple of the two element symbols and the values are the bond radii. + Parameters + ---------- + radius_dict : dict + Dictionary of bond radii. + The keys are a tuple of the two element symbols and the values are the bond radii. """ # convert radius_dict to elements and radius @@ -1117,14 +1125,15 @@ def set_bond_radius(self, radius_dict): _raffle.f90wrap_raffle__dc__set_bond_radius__binding__dc_type(this=self._handle, \ elements=elements, radius=radius) - def set_bond_radii(self, radius_dict): + def set_bond_radii(self, radius_dict : dict): """ Set the bond radii for the distribution functions. - Parameters: - radius_dict (dict): - Dictionary of bond radii. - The keys are a tuple of the two element symbols and the values are the bond radii. + Parameters + ---------- + radius_dict : dict + Dictionary of bond radii. + The keys are a tuple of the two element symbols and the values are the bond radii. """ # convert radius_list to elements and radii @@ -1143,10 +1152,11 @@ def get_bond_radii(self): """ Get the bond radii for the distribution functions. - Returns: - bond_radii (dict): - Dictionary of bond radii. - The keys are a tuple of the two element symbols and the values are the bond radii. + Returns + ------- + dict + Dictionary of bond radii. + The keys are a tuple of the two element symbols and the values are the bond radii. """ num_elements = _raffle.f90wrap_raffle__dc__get__num_elements(self._handle) @@ -1191,135 +1201,148 @@ def evolve(self): #, system=None): def is_converged(self, threshold : float = 1e-4): """ - Parameters: - threshold (float): - Threshold for convergence. + Parameters + ---------- + threshold : float + Threshold for convergence. - Returns: - converged (bool): - Boolean indicating whether the distribution functions have converged. + Returns + ------- + 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): + def write_gdfs(self, file : str): """ Write the generalised distribution functions (GDFs) to a file. - Parameters: - file (str): - Name of file to write the GDFs to. + Parameters + ---------- + file : str + Name of file to write the GDFs to. """ _raffle.f90wrap_raffle__dc__write_gdfs__binding__dc_type(this=self._handle, \ file=file) - def read_gdfs(self, file): + def read_gdfs(self, file : str): """ Read the generalised distribution functions (GDFs) from a file. - Parameters: - file (str): - Name of file to read the GDFs from. + Parameters + ---------- + file : str + Name of file to read the GDFs from. """ _raffle.f90wrap_raffle__dc__read_gdfs__binding__dc_type(this=self._handle, \ file=file) - def write_dfs(self, file): + def write_dfs(self, file : str): """ Write the distribution functions (DFs) associated with all allocated systems to a file. - Parameters: - file (str): - Name of file to write the DFs to. + Parameters + ---------- + file : str + Name of file to write the DFs to. """ _raffle.f90wrap_raffle__dc__write_dfs__binding__dc_type(this=self._handle, \ file=file) - def read_dfs(self, file): + def read_dfs(self, file : str): """ Read the distribution functions (DFs) associated with a set of systems from a file. - Parameters: - file (str): - Name of file to read the DFs from. + Parameters + ---------- + file : str + Name of file to read the DFs from. """ _raffle.f90wrap_raffle__dc__read_dfs__binding__dc_type(this=self._handle, \ file=file) - def write_2body(self, file): + def write_2body(self, file: str): """ Write the 2-body generalised distribution functions (GDFs) to a file. - Parameters: - file (str): - Name of file to write the 2-body GDFs to. + Parameters + ---------- + file : str + Name of file to write the 2-body GDFs to. """ _raffle.f90wrap_raffle__dc__write_2body__binding__dc_type(this=self._handle, \ file=file) - def write_3body(self, file): + def write_3body(self, file: str): """ Write the 3-body generalised distribution functions (GDFs) to a file. - Parameters: - file (str): - Name of file to write the 3-body GDFs to. + Parameters + ---------- + file : str + Name of file to write the 3-body GDFs to. """ _raffle.f90wrap_raffle__dc__write_3body__binding__dc_type(this=self._handle, \ file=file) - def write_4body(self, file): + def write_4body(self, file: str): """ Write the 4-body generalised distribution functions (GDFs) to a file. - Parameters: - file (str): - Name of file to write the 4-body GDFs to. + Parameters + ---------- + file : str + Name of file to write the 4-body GDFs to. """ _raffle.f90wrap_raffle__dc__write_4body__binding__dc_type(this=self._handle, \ file=file) - def _get_pair_index(self, species1, species2): + def _get_pair_index(self, species1 : str, species2 : str): """ Get the index of the pair of species in the distribution functions. This is meant as an internal function and not likely to be used directly. - Parameters: - species1 (str): - Name of the first species - species2 (str): - Name of the second species + Parameters + ---------- + species1 : str + Name of the first species + species2 : str + Name of the second species - Returns: - idx (int): - Index of the pair of species in the distribution functions. + Returns + ------- + int + Index of the pair of species in the distribution functions. """ idx = \ _raffle.f90wrap_raffle__dc__get_pair_index__binding__dc_type(this=self._handle, \ species1=species1, species2=species2) return idx - def _get_bin(self, value, dim): + def _get_bin(self, value : float, dim : int): """ Get the bin index for a value in the distribution functions. This is meant as an internal function and not likely to be used directly. - Parameters: - value (float): - Value to get the bin index for. - dim (int): - Dimension of the distribution function. - 1 for 2-body, 2 for 3-body, and 3 for 4-body. + Parameters + ---------- + value : float + Value to get the bin index for. + dim : int + Dimension of the distribution function. + 1 for 2-body, 2 for 3-body, and 3 for 4-body. - Returns: - bin (int): - Bin index for the value in the distribution functions. + Returns + ------- + int + Bin index for the value in the distribution functions. """ bin = \ @@ -1331,9 +1354,10 @@ 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. + Returns + ------- + float array: 2D array of shape (num_species_pairs, nbins) + 2-body distribution function. """ n0 = self.nbins[0] @@ -1346,9 +1370,10 @@ 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. + Returns + ------- + float array: 2D array of shape (num_species, nbins) + 3-body distribution function. """ n0 = self.nbins[1] @@ -1361,9 +1386,10 @@ 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. + Returns + ------- + 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) @@ -1371,6 +1397,42 @@ def get_4body(self): _raffle.f90wrap_raffle__dc__get_4body__binding__dc_type(this=self._handle, n0=n0, n1=n1).reshape((n0, n1), order='F').T return output + def generate_fingerprint(self, + structure: Atoms | Geom_Rw.basis = None, + ): + """ + Generate the fingerprint for a given structure. + + Parameters + ---------- + structure : ase.Atoms + Atomic structure to generate the fingerprint for. + + Returns + ------- + output : list[arrays] + 2D arrays + """ + if isinstance(structure, Atoms): + structure = geom_rw.basis(structure) + + num_species = structure.nspec + num_pairs = num_species * (num_species + 1) // 2 + nbins = [] + for i in range(len(self.nbins)): + nbins.append(self.nbins[i]) + if nbins[i] < 0: + nbins[i] = 1 + round( ( self.cutoff_max[i] - self.cutoff_min[i] ) / self.width[i] ) + + output_2body = numpy.asfortranarray(numpy.zeros((nbins[0], num_pairs), dtype=numpy.float32)) + output_3body = numpy.asfortranarray(numpy.zeros((nbins[1], num_species), dtype=numpy.float32)) + output_4body = numpy.asfortranarray(numpy.zeros((nbins[2], num_species), dtype=numpy.float32)) + _raffle.f90wrap_raffle__dc__generate_fingerprint_python__dc_type(this=self._handle, \ + structure=structure._handle, output_2body=output_2body, \ + output_3body=output_3body, output_4body=output_4body) + + return output_2body, output_3body, output_4body + @property def num_evaluated(self): """ @@ -1452,6 +1514,19 @@ def weight_by_hull(self, weight_by_hull): _raffle.f90wrap_distribs_container_type__set__weight_by_hull(self._handle, \ weight_by_hull) + @property + def smooth_viability(self): + """ + DEVELOPER TOOL. Boolean whether to smooth the viability evaluation of points. + """ + return \ + _raffle.f90wrap_distribs_container_type__get__smooth_viability(self._handle) + + @smooth_viability.setter + def smooth_viability(self, smooth_viability): + _raffle.f90wrap_distribs_container_type__set__smooth_viability(self._handle, \ + smooth_viability) + @property def nbins(self): """ @@ -1636,9 +1711,10 @@ def __init__(self, dict=None, element=None, num=None, handle=None): """ Object to store the stoichiometry of a structure. - Returns: - stoichiometry (stoichiometry_type): - Stoichiometry object + Returns + ------- + stoichiometry_type + Stoichiometry object """ f90wrap.runtime.FortranDerivedType.__init__(self) result = _raffle.f90wrap_stoichiometry_type_initialise() @@ -1654,10 +1730,6 @@ def __del__(self): """ Destructor for class Stoichiometry_Type - - Defined at ../src/lib/mod_generator.f90 lines \ - 19-21 - Automatically generated destructor for stoichiometry_type """ if self._alloc: @@ -1703,9 +1775,10 @@ def __init__(self, dict=None, handle=None): """ Array or list of stoichiometry objects. - Returns: - stoichiometry_array (stoichiometry_array): - Stoichiometry array object + Returns + ------- + Generator.stoichiometry_array + Stoichiometry array object """ f90wrap.runtime.FortranDerivedType.__init__(self) @@ -1724,16 +1797,6 @@ def __del__(self): """ Destructor for class Stoichiometry_Type - - Defined at ../src/lib/mod_generator.f90 lines \ - 19-21 - - Parameters - ---------- - this : Stoichiometry_Type - Object to be destructed - - Automatically generated destructor for stoichiometry_type """ if self._alloc: @@ -1757,9 +1820,10 @@ def allocate(self, size): """ Allocate the items array with the given size. - Parameters: - size (int): - Size of the items array + Parameters + ---------- + size : int + Size of the items array """ _raffle.f90wrap_stoich_type_xnum_array__array_alloc__items(self._handle, num=size) @@ -1777,7 +1841,7 @@ def deallocate(self): @f90wrap.runtime.register_class("raffle.raffle_generator") class raffle_generator(f90wrap.runtime.FortranDerivedType): - def __init__(self, handle=None, history_len=None): + def __init__(self, handle=None, history_len : int = None): """ Create a ``raffle_generator`` object. @@ -1785,9 +1849,10 @@ def __init__(self, handle=None, history_len=None): 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. + 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() @@ -1811,13 +1876,14 @@ def init_seed(self, put: int|list[int] = None, get: list[int] = None, num_thread """ Initialise the random number generator seed. - Parameters: - put (list of ints): - List of integers to be used as the seed for the random number generator. - get (list of ints): - List of integers to be used as the seed for the random number generator. - num_threads (int): - Number of threads, one for each random number generator. + Parameters + ---------- + put : list[int] + List of integers to be used as the seed for the random number generator. + get : list[int] + List of integers to be used as the seed for the random number generator. + num_threads : int + Number of threads, one for each random number generator. """ # check if put is an integer @@ -1827,15 +1893,61 @@ def init_seed(self, put: int|list[int] = None, get: list[int] = None, num_thread _raffle.f90wrap_generator__init_seed__binding__rgt(this=self._handle, \ put=put, get=get, num_threads=num_threads) - def set_max_attempts(self, max_attempts): + def set_method_ratio_default(self, method_ratio : dict[str, float]): + """ + Set the method of the five placement methods to be used. + + This will be used as the default ratio if one is not provided as an input argument to the generate() method. + + Parameters + ---------- + method_ratio : dict[str, float] + The ratio of the walk method to the void method. + """ # check if method_ratio is provided, if so, use it only if method_ratio is not provided + + # handle other names for the methods, such as random, growth, minimum, etc. + if "random" in method_ratio: + if "rand" in method_ratio: + raise ValueError("Both random and rand are provided, use only one (random or rand)") + method_ratio["rand"] = method_ratio.pop("random") + if "growth" in method_ratio: + if "grow" in method_ratio: + raise ValueError("Both growth and grow are provided, use only one (growth or grow)") + method_ratio["grow"] = method_ratio.pop("growth") + if "minimum" in method_ratio: + if "min" in method_ratio: + raise ValueError("Both minimum and min are provided, use only one (minimum or min)") + method_ratio["min"] = method_ratio.pop("minimum") + if "global" in method_ratio: + if "min" in method_ratio: + raise ValueError("Both global and min are provided, use only one (global or min)") + method_ratio["min"] = method_ratio.pop("global") + + # exit if any other method is provided + for key in method_ratio.keys(): + 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 = [] + method_ratio_list.append(method_ratio.get("void", 0.0)) + method_ratio_list.append(method_ratio.get("rand", 0.0)) # or method_ratio.get("random", 0.0)) + method_ratio_list.append(method_ratio.get("walk", 0.0)) + method_ratio_list.append(method_ratio.get("grow", 0.0)) # or method_ratio.get("growth", 0.0)) + method_ratio_list.append(method_ratio.get("min", 0.0)) # or method_ratio.get("minimum", 0.0) or method_ratio.get("global", 0.0)) + + _raffle.f90wrap_generator__set_method_ratio_default__binding__rgt(this=self._handle, \ + method_ratio=method_ratio_list) + + def set_max_attempts(self, max_attempts : int): """ Set the walk-method maximum attempts parameter. This parameter determines the maximum number of attempts to generate a structure using the walk method before reverting to the void method. - Parameters: - max_attempts (int): - The maximum number of attempts to generate a structure using the walk method. + Parameters + ---------- + max_attempts : int + The maximum number of attempts to generate a structure using the walk method. """ self.max_attempts = max_attempts @@ -1843,24 +1955,26 @@ def set_walk_step_size(self, coarse=None, fine=None): """ Set the walk-method step size parameters. - Parameters: - coarse (float): - The coarse step size for the walk method. - fine (float): - The fine step size for the walk method. + Parameters + ---------- + coarse : float + The coarse step size for the walk method. + fine : float + The fine step size for the walk method. """ if coarse is not None: self.walk_step_size_coarse = coarse if fine is not None: self.walk_step_size_fine = fine - def set_host(self, host): + def set_host(self, host : Atoms | Geom_Rw.basis): """ Set the host structure for the generation. - Parameters: - host (ase.Atoms or geom_rw.basis): - The host structure for the generation. + Parameters + ---------- + host : ase.Atoms or Geom_Rw.basis + The host structure for the generation. """ # check if host is ase.Atoms object or a Fortran derived type basis_type if isinstance(host, Atoms): @@ -1868,7 +1982,7 @@ def set_host(self, host): elif isinstance(host, geom_rw.basis): host_copy = host.copy() else: - raise TypeError("Expected host to be ASE Atoms or geom_rw.basis") + raise TypeError("Expected host to be ASE ase.Atoms or Geom_Rw.basis") _raffle.f90wrap_generator__set_host__binding__rgt(this=self._handle, \ host=host_copy._handle) @@ -1877,9 +1991,10 @@ def get_host(self): """ Get the host structure for the generation. - Returns: - host (ase.Atoms or geom_rw.basis): - The host structure for the generation. + Returns + ------- + ase.Atoms + The host structure for the generation. """ output = \ _raffle.f90wrap_generator__get_host__binding__rgt(this=self._handle) @@ -1892,28 +2007,34 @@ def get_host(self): ret_host = None return ret_host - def prepare_host(self, interface_location=None, interface_axis=3, depth=3.0, location_as_fractional=False): + def prepare_host(self, + interface_location : list[float] = None, + interface_axis : int | str = 3, + depth : float = 3.0, + location_as_fractional : bool = False + ): """ Prepare the host by removing atoms depth away from the interface and returning an associated missing stoichiometry dictionary. - Parameters: - interface_location (list of float): - The location of the interface in direct coordinates. - interface_axis (int or str): - The axis along which to remove atoms from the host. - 0 for a, 1 for b, 2 for c. - Alternatively, 'a', 'b', or 'c' can be used. - depth (float): - The depth to remove atoms from the host. - location_as_fractional (bool): - Whether the interface location is given in fractional coordinates. - - Returns: - missing_stoichiometry (dict): - A dictionary of the missing stoichiometry. - The keys are the element symbols and the values are the number of atoms removed. + Parameters + ---------- + interface_location : list[float] + The location of the interface in direct coordinates. + interface_axis (int or str): + The axis along which to remove atoms from the host. + 0 for a, 1 for b, 2 for c. + Alternatively, 'a', 'b', or 'c' can be used. + depth : float + The depth to remove atoms from the host. + location_as_fractional : bool + Whether the interface location is given in fractional coordinates. + + Returns + ------- + dict + A dictionary of the missing stoichiometry. + The keys are the element symbols and the values are the number of atoms removed. """ - # check if interface_axis is a string interface_axis_int = interface_axis if isinstance(interface_axis_int, str): @@ -1952,17 +2073,18 @@ def prepare_host(self, interface_location=None, interface_axis=3, depth=3.0, loc return missing_stoichiometry - def set_grid(self, grid=None, grid_spacing=None, grid_offset=None): + def set_grid(self, grid : list[int] = None, grid_spacing : float = None, grid_offset : list[float] = None): """ Set the grid parameters for the generation. - Parameters: - grid (list of int): - The number of grid points along each axis of the host. - grid_spacing (float) - The spacing between grid points. - grid_offset (list of float): - The offset of the grid from the origin. + Parameters + ---------- + grid : list[int] + The number of grid points along each axis of the host. + grid_spacing : float + The spacing between grid points. + grid_offset : list[float] + The offset of the grid from the origin. """ _raffle.f90wrap_generator__set_grid__binding__raffle_generator_type(this=self._handle, \ grid=grid, grid_spacing=grid_spacing, grid_offset=grid_offset) @@ -1977,12 +2099,21 @@ def set_bounds(self, bounds=None): """ Set the bounding box for the generation. - Parameters: - bounds (list of list of 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]]. - Values given in direct (crystal) coordinates, ranging from 0 to 1. - """ + Parameters + ---------- + 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]]. + Values given in direct (crystal) coordinates, ranging from 0 to 1. + """ + # check if bounds is a list of length 2 + if not isinstance(bounds, list) or len(bounds) != 2: + raise ValueError("bounds must be a list of length 2") + # check each element of bounds is a list of length 3 + for bound in bounds: + 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") + _raffle.f90wrap_generator__set_bounds__binding__rgt(this=self._handle, \ bounds=bounds) @@ -1992,41 +2123,49 @@ def reset_bounds(self): """ _raffle.f90wrap_generator__reset_bounds__binding__rgt(this=self._handle) - def generate( - self, num_structures, stoichiometry, - method_ratio: dict = {"void": 0.0, "rand": 0.0, "walk": 0.0, "grow": 0.0, "min": 0.0}, - method_probab=None, - seed=None, settings_out_file=None, verbose=0, return_exit_code=False, - calc=None + def generate(self, + num_structures : int, + stoichiometry, + method_ratio: dict[str, float] = {"void": 0.0, "rand": 0.0, "walk": 0.0, "grow": 0.0, "min": 0.0}, + method_probab : dict = None, + seed : int = None, + settings_out_file : str = None, + verbose : int = 0, + return_exit_code : bool = False, + calc = None ): """ Generate structures using the RAFFLE method. - Parameters: - num_structures (int): - The number of structures to generate. - stoichiometry (stoichiometry_array, str, or dict): - The stoichiometry of the structures to generate. - method_ratio (dict): - The ratio of using each method to generate a structure. - method_probab (dict): - DEPRECATED - The ratio of using each method to generate a structure. - seed (int): - The seed for the random number generator. - settings_out_file (str): - The file to write the settings to. - verbose (int): - The verbosity level for the generation. - return_exit_code (bool): - Whether to return the exit code from the generation. - calc (ASE calculator): - The calculator to use for the generated structures. - - Returns: - structures (ase.Atoms): - The generated structures. - exit_code (int): - The exit code from the generation. + + + Parameters + ---------- + num_structures : int + The number of structures to generate. + stoichiometry : Generator.stoichiometry_array, str, or dict + The stoichiometry of the structures to generate. + method_ratio : dict + The ratio of using each method to generate a structure. + method_probab : dict + DEPRECATED. Use `method_ratio` instead. Will be removed in a future version. + seed : int + The seed for the random number generator. + settings_out_file : str + The file to write the settings to. + verbose : int + The verbosity level for the generation. + return_exit_code : bool + Whether to return the exit code from the generation. + calc : ase.calculators.calculator.BaseCalculator + The calculator to use for the generated structures. + + Returns + ------- + structures : ase.Atoms + The generated structures. + exit_code : int + The exit code from the generation. """ exit_code = 0 @@ -2073,7 +2212,7 @@ def generate( # check if all values are 0.0, if so, set them to the default of all 1.0 if all([val < 1E-6 for val in method_ratio_list]): - method_ratio_list = [1.0, 0.1, 0.5, 0.5, 1.0] + method_ratio_list = None # if stoichiometry is a dictionary, convert it to a stoichiometry_array if isinstance(stoichiometry, dict): @@ -2111,9 +2250,10 @@ def get_structures(self, calculator=None): """ Get the generated structures as a list of ASE Atoms objects. - Parameters: - calculator (ASE calculator): - The calculator to use for the generated structures. + Parameters + ---------- + calculator : ase.calculators.calculator.BaseCalculator + The calculator to use for the generated structures. """ atoms = [] for structure in self.structures: @@ -2124,21 +2264,23 @@ def set_structures(self, structures): """ Set the list of generated structures. - Parameters: - structures (list of geom_rw.basis or list of ase.Atoms): - The list of structures to set. + Parameters + ---------- + structures : list[geom_rw.basis] or list[ase.Atoms] + The list of structures to set. """ structures = geom_rw.basis_array(atoms=structures) _raffle.f90wrap_generator__set_structures__binding__rgt(this=self._handle, \ structures=structures._handle) - def remove_structure(self, index): + def remove_structure(self, index : int | list[int]): """ Remove the structure at the given indices from the generator. - Parameters: - index (int or list of int): - The indices of the structure to remove. + Parameters + ---------- + index : int or list[int] + The indices of the structure to remove. """ index_list = [index] if isinstance(index, int) else index index_list = [ i + 1 for i in index_list ] @@ -2149,15 +2291,15 @@ def evaluate(self, basis): """ Evaluate the viability of the structures. - WARNING: This function is not implemented yet. - - Parameters: - basis (geom_rw.basis or Atoms): - The basis to use for the evaluation. + Parameters + ---------- + basis : Geom_Rw.basis or ase.Atoms + The basis to use for the evaluation. - Returns: - viability (float): - The viability of the structures. + Returns + ------- + float + The viability of the structures. """ if isinstance(basis, Atoms): basis = geom_rw.basis(atoms=basis) @@ -2167,24 +2309,94 @@ def evaluate(self, basis): basis=basis._handle) return viability - def print_settings(self, file): + def get_probability_density(self, + basis, + species : list[str] | str, + grid : list[int] = None, + grid_offset : list[float] = None, + grid_spacing : float = None, + bounds : list[list[float]] = None, + return_grid : bool = False + ): + """ + Get the probability density of the generated structures. + + Parameters + ---------- + basis : Geom_Rw.basis or ase.Atoms + The basis to use for the evaluation. + species : list[str] or str + The species to use for the evaluation. + grid : list[int] + The number of grid points along each axis of the host. + grid_offset : list[float] + The offset of the grid from the origin. + grid_spacing : float + The spacing between grid points. + 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]]. + Values given in direct (crystal) coordinates, ranging from 0 to 1. + return_grid : bool + Whether to return the number of grid points along each axis of the host. + + Returns + ------- + numpy.ndarray (4 + num_species, n_points) + The probability density of the generated structures. + list[int] + The grid parameters used for the generation. + """ + probability_density = None + ret_grid = numpy.zeros(3, dtype=numpy.int32) + if isinstance(basis, Atoms): + basis = geom_rw.basis(atoms=basis) + + if isinstance(species, str): + species = re.findall(r'[A-Z][a-z]?\d*', species) + species = [re.sub(r'\d+', '', s) for s in species] + + n_coords, n_points, ret_grid = _raffle.f90wrap_generator__get_probability_density__rgt( + this = self._handle, + basis = basis._handle, + species_list = species, + grid = grid, + grid_offset = grid_offset, + grid_spacing = grid_spacing, + bounds = bounds, + ) + ret_grid = ret_grid.tolist() + + + # allocate the structures + probability_density = numpy.asfortranarray(numpy.zeros((n_coords, n_points), dtype=numpy.float32)) + _raffle.f90wrap_retrieve_probability_density(probability_density) + + if return_grid: + return probability_density, ret_grid + else: + return probability_density + + def print_settings(self, file : str): """ Print the settings for the generation to a file. - Parameters: - file (str): - Name of the file to write the settings to. + Parameters + ---------- + file : str + Name of the file to write the settings to. """ _raffle.f90wrap_generator__print_settings__binding__rgt(this=self._handle, \ file=file) - def read_settings(self, file): + def read_settings(self, file : str): """ Read the settings for the generation from a file. - Parameters: - file (str): - Name of the file to read the settings from. + Parameters + ---------- + file : str + Name of the file to read the settings from. """ _raffle.f90wrap_generator__read_settings__binding__rgt(this=self._handle, \ file=file) @@ -2193,9 +2405,10 @@ def get_descriptor(self): """ Get the descriptor for the generated structures. - Returns: - descriptor (list of numpy arrays): - The descriptor for the generated structures. + Returns + ------- + list[numpy arrays] + The descriptor for the generated structures. """ descriptor = [ self.distributions.get_2body(), @@ -2368,7 +2581,7 @@ def walk_step_size_fine(self, walk_step_size_fine): @property def method_ratio(self): """ - The ratio of methods to employ to generate a structure. + The last used ratio of methods to employ to generate a structure. """ array_ndim, array_type, array_shape, array_handle = \ _raffle.f90wrap_raffle_generator_type__array__method_ratio(self._handle) @@ -2385,6 +2598,26 @@ def method_ratio(self): def method_ratio(self, method_ratio): self.method_ratio[...] = method_ratio + @property + def method_ratio_default(self): + """ + The default ratio of methods to employ to generate a structure. + """ + array_ndim, array_type, array_shape, array_handle = \ + _raffle.f90wrap_raffle_generator_type__array__method_ratio_default(self._handle) + if array_handle in self._arrays: + method_ratio = self._arrays[array_handle] + else: + method_ratio = f90wrap.runtime.get_array(f90wrap.runtime.sizeof_fortran_t, + self._handle, + _raffle.f90wrap_raffle_generator_type__array__method_ratio) + self._arrays[array_handle] = method_ratio + return method_ratio + + @method_ratio_default.setter + def method_ratio_default(self, method_ratio): + self.method_ratio_default[...] = method_ratio + def _init_array_structures(self): """ Initialise the structures array. @@ -2398,10 +2631,6 @@ def _init_array_structures(self): """ Element items ftype=type(basis_type) pytype=basis - - Defined at ../src/lib/mod_generator.f90 line \ - 29 - """, Geom_Rw.basis) return self.structures @@ -2441,10 +2670,6 @@ def __str__(self): class Raffle(f90wrap.runtime.FortranModule): """ Module raffle - - - Defined at ../src/raffle.f90 lines 1-4 - """ pass _dt_array_initialisers = [] diff --git a/src/wrapper/f90wrap_mod_distribs_container.f90 b/src/wrapper/f90wrap_mod_distribs_container.f90 index 2b74713..93811a1 100644 --- a/src/wrapper/f90wrap_mod_distribs_container.f90 +++ b/src/wrapper/f90wrap_mod_distribs_container.f90 @@ -269,6 +269,38 @@ subroutine f90wrap_distribs_container_type__set__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 + +subroutine f90wrap_distribs_container_type__get__smooth_viability( & + this, f90wrap_smooth_viability & +) + 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 + logical, intent(out) :: f90wrap_smooth_viability + + this_ptr = transfer(this, this_ptr) + f90wrap_smooth_viability = this_ptr%p%smooth_viability +end subroutine f90wrap_distribs_container_type__get__smooth_viability + +subroutine f90wrap_distribs_container_type__set__smooth_viability( & + this, f90wrap_smooth_viability & +) + 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 + logical, intent(in) :: f90wrap_smooth_viability + + this_ptr = transfer(this, this_ptr) + this_ptr%p%smooth_viability = f90wrap_smooth_viability +end subroutine f90wrap_distribs_container_type__set__smooth_viability !############################################################################### @@ -1107,6 +1139,47 @@ subroutine f90wrap_raffle__dc__get_4body__binding__dc_type( & ret_output = this_ptr%p%get_4body() end subroutine f90wrap_raffle__dc__get_4body__binding__dc_type +subroutine f90wrap_raffle__dc__generate_fingerprint_python__dc_type( & + this, structure, output_2body, output_3body, output_4body, & + n0, n1, n2, n3, n4, n5 & +) + use raffle__distribs_container, only: distribs_container_type + use raffle__geom_rw, only: basis_type + implicit none + + type basis_type_ptr_type + type(basis_type), pointer :: p => NULL() + end type basis_type_ptr_type + 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 + type(basis_type_ptr_type) :: structure_ptr + integer, intent(in), dimension(2) :: structure + real(4), intent(inout), dimension(n0,n1) :: output_2body + real(4), intent(inout), dimension(n2,n3) :: output_3body + real(4), intent(inout), dimension(n4,n5) :: output_4body + integer :: n0 + !f2py intent(hide), depend(output_2body) :: n0 = shape(output_2body,0) + integer :: n1 + !f2py intent(hide), depend(output_2body) :: n1 = shape(output_2body,1) + integer :: n2 + !f2py intent(hide), depend(output_3body) :: n2 = shape(output_3body,0) + integer :: n3 + !f2py intent(hide), depend(output_3body) :: n3 = shape(output_3body,1) + integer :: n4 + !f2py intent(hide), depend(output_4body) :: n4 = shape(output_4body,0) + integer :: n5 + !f2py intent(hide), depend(output_4body) :: n5 = shape(output_4body,1) + this_ptr = transfer(this, this_ptr) + structure_ptr = transfer(structure, structure_ptr) + call this_ptr%p%generate_fingerprint_python( & + structure=structure_ptr%p, output_2body=output_2body, & + output_3body=output_3body, output_4body=output_4body & + ) +end subroutine f90wrap_raffle__dc__generate_fingerprint_python__dc_type + subroutine f90wrap_raffle__dc__get_num_species__dc_type( & this, ret_num_species & ) diff --git a/src/wrapper/f90wrap_mod_generator.f90 b/src/wrapper/f90wrap_mod_generator.f90 index 033b31d..8b2d808 100644 --- a/src/wrapper/f90wrap_mod_generator.f90 +++ b/src/wrapper/f90wrap_mod_generator.f90 @@ -606,6 +606,29 @@ subroutine f90wrap_raffle_generator_type__array__method_ratio( & dshape(1:1) = shape(this_ptr%p%method_ratio) dloc = loc(this_ptr%p%method_ratio) end subroutine f90wrap_raffle_generator_type__array__method_ratio + +subroutine f90wrap_raffle_generator_type__array__method_ratio_default( & + this, nd, dtype, dshape, dloc & +) + use raffle__generator, only: raffle_generator_type + use, intrinsic :: iso_c_binding, only : c_int + implicit none + type raffle_generator_type_ptr_type + type(raffle_generator_type), pointer :: p => NULL() + end type raffle_generator_type_ptr_type + integer(c_int), intent(in) :: this(2) + type(raffle_generator_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) + dshape(1:1) = shape(this_ptr%p%method_ratio_default) + dloc = loc(this_ptr%p%method_ratio_default) +end subroutine f90wrap_raffle_generator_type__array__method_ratio_default !############################################################################### @@ -734,6 +757,28 @@ end subroutine f90wrap_generator__raffle_generator_type_finalise !############################################################################### +!############################################################################### +! set placement method ratio +!############################################################################### +subroutine f90wrap_generator__set_method_ratio_default__binding__rgt( & + this, method_ratio & +) + use raffle__generator, only: raffle_generator_type + implicit none + + type raffle_generator_type_ptr_type + type(raffle_generator_type), pointer :: p => NULL() + end type raffle_generator_type_ptr_type + type(raffle_generator_type_ptr_type) :: this_ptr + integer, intent(in), dimension(2) :: this + real(4), dimension(5), intent(in) :: method_ratio + + this_ptr = transfer(this, this_ptr) + call this_ptr%p%set_method_ratio_default(method_ratio=method_ratio) +end subroutine f90wrap_generator__set_method_ratio_default__binding__rgt +!############################################################################### + + !############################################################################### ! initialise random seed !############################################################################### @@ -1031,6 +1076,74 @@ subroutine f90wrap_generator__evaluate__binding__rgt(this, ret_viability, basis) end subroutine f90wrap_generator__evaluate__binding__rgt +subroutine f90wrap_generator__get_probability_density__rgt( & + this, basis, species_list, grid, grid_offset, grid_spacing, bounds, & + n_ret_coords, n_ret_points, ret_grid, n0 & +) + use raffle__generator, only: raffle_generator_type + use raffle__geom_rw, only: basis_type + use raffle__cache, only: store_probability_density + 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 + type(raffle_generator_type_ptr_type) :: this_ptr + integer, intent(in), dimension(2) :: this + type(basis_type_ptr_type) :: basis_ptr + integer, intent(in), dimension(2) :: basis + character(len=3), intent(in), dimension(n0) :: species_list + integer, dimension(3), intent(in), optional :: grid + real(4), dimension(3), intent(in), optional :: grid_offset + real(4), intent(in), optional :: grid_spacing + real(4), dimension(2,3), intent(in), optional :: bounds + integer, intent(out) :: n_ret_coords + integer, intent(out) :: n_ret_points + integer, dimension(3), intent(out) :: ret_grid + integer :: n0 + !f2py intent(hide), depend(species_list) :: n0 = shape(species_list,0) + + ! Local temporary array + real(4), allocatable, dimension(:,:) :: local_probability_density + + ! Call the actual function + this_ptr = transfer(this, this_ptr) + basis_ptr = transfer(basis, basis_ptr) + local_probability_density = this_ptr%p%get_probability_density( & + basis = basis_ptr%p, & + species_list = species_list, & + grid = grid, & + grid_offset = grid_offset, & + grid_spacing = grid_spacing, & + bounds = bounds, & + grid_output = ret_grid & + ) + + n_ret_coords = size(local_probability_density, 1) + n_ret_points = size(local_probability_density, 2) + + ! Store local_probability_density in the cache so Python can retrieve it + call store_probability_density( local_probability_density ) +end subroutine f90wrap_generator__get_probability_density__rgt + + +subroutine f90wrap_retrieve_probability_density(probability_density, n0, n1) + use raffle__cache, only: retrieve_probability_density + implicit none + + real(4), dimension(n0,n1) :: probability_density + integer :: n0 + !f2py intent(hide), depend(probability_density) :: n0 = shape(probability_density,0) + integer :: n1 + !f2py intent(hide), depend(probability_density) :: n1 = shape(probability_density,1) + + probability_density = retrieve_probability_density() +end subroutine f90wrap_retrieve_probability_density + + subroutine f90wrap_generator__print_settings__binding__rgt( & this, file & ) diff --git a/test/test_dist_calcs.f90 b/test/test_dist_calcs.f90 index 025e092..9842aad 100644 --- a/test/test_dist_calcs.f90 +++ b/test/test_dist_calcs.f90 @@ -33,6 +33,7 @@ program test_edit_geom bas%lat(1,:) = [0.0, 2.14, 2.14] bas%lat(2,:) = [2.14, 0.0, 2.14] bas%lat(3,:) = [2.14, 2.14, 0.0] + call bas%set_atom_mask() !----------------------------------------------------------------------------- @@ -86,6 +87,7 @@ program test_edit_geom bas2%spec(2)%name = 'Ge' allocate(bas2%spec(2)%atom(bas2%spec(2)%num, 3)) bas2%spec(2)%atom(1, :) = [0.25, 0.25, 0.25] + call bas2%set_atom_mask() rtmp1 = get_min_dist_between_point_and_species( & bas2, & @@ -98,8 +100,8 @@ program test_edit_geom rtmp2 = modu(loc) if ( abs(rtmp1 - rtmp2) .gt. 1.E-6 ) then - write(0,*) 'get_min_dist_between_point_and_species failed' - success = .false. + write(0,*) 'get_min_dist_between_point_and_species failed' + success = .false. end if rtmp1 = get_min_dist_between_point_and_species( & @@ -111,10 +113,10 @@ program test_edit_geom loc = loc - ceiling(loc - 0.5) loc = matmul(loc, bas2%lat) rtmp2 = modu(loc) - + if ( abs(rtmp1 - rtmp2) .gt. 1.E-6 ) then - write(0,*) 'get_min_dist_between_point_and_species failed' - success = .false. + write(0,*) 'get_min_dist_between_point_and_species failed' + success = .false. end if @@ -122,9 +124,9 @@ program test_edit_geom ! Test get_dist_between_point_and_atom !----------------------------------------------------------------------------- rtmp1 = get_dist_between_point_and_atom( & - bas, & - loc=[0.9, 0.9, 0.9], & - atom=[1, 1] & + bas, & + loc=[0.9, 0.9, 0.9], & + atom=[1, 1] & ) loc = bas%spec(1)%atom(1,:3) - [0.9, 0.9, 0.9] loc = matmul(loc, bas%lat) @@ -135,15 +137,15 @@ program test_edit_geom success = .false. end if rtmp1 = get_dist_between_point_and_atom( & - bas, & - loc=[0.9, 0.9, 0.9], & - atom=[1, 2] & + bas, & + loc=[0.9, 0.9, 0.9], & + atom=[1, 2] & ) loc = bas%spec(1)%atom(2,:3) - [0.9, 0.9, 0.9] loc = matmul(loc, bas%lat) rtmp2 = modu(loc) if ( abs(rtmp1 - rtmp2) .gt. 1.E-6 ) then - write(*,*) rtmp1, rtmp2 + write(*,*) rtmp1, rtmp2 write(0,*) 'get_dist_between_point_and_atom failed' success = .false. end if @@ -160,4 +162,4 @@ program test_edit_geom stop 1 end if -end program test_edit_geom \ No newline at end of file +end program test_edit_geom diff --git a/test/test_distribs_container.f90 b/test/test_distribs_container.f90 index a28adea..45aafea 100644 --- a/test/test_distribs_container.f90 +++ b/test/test_distribs_container.f90 @@ -74,6 +74,8 @@ program test_distribs_container basis_mgo%energy = -20.0 call test_init_distribs_container(success) + call test_set_history_len(success) + call test_is_converged(success) call test_set_width(success) call test_set_sigma(success) call test_set_cutoff_min(success) @@ -115,11 +117,13 @@ subroutine test_init_distribs_container(success) class(distribs_container_type), allocatable :: distribs_container character(len=10) :: test_name + integer :: history_len integer, dimension(3) :: nbins real(real32), dimension(3) :: width, sigma, cutoff_min, cutoff_max ! Test case 1: Default initialisation distribs_container = distribs_container_type() + history_len = 0 nbins = [-1, -1, -1] width = [0.025_real32, pi/64._real32, pi/64._real32] sigma = [0.1_real32, 0.1_real32, 0.1_real32] @@ -128,6 +132,19 @@ subroutine test_init_distribs_container(success) test_name = "Default" do i = 1, 2 + call assert( & + distribs_container%history_len .eq. history_len, & + trim(test_name)//" history_len initialisation failed", & + success & + ) + if(history_len.gt.0)then + call assert( & + allocated( distribs_container%history_deltas) .and. & + size( distribs_container%history_deltas, 1) .eq. history_len, & + trim(test_name)//" history_deltas initialisation failed", & + success & + ) + end if call assert( & all( distribs_container%nbins .eq. nbins ), & trim(test_name)//" nbins initialisation failed", & @@ -160,12 +177,14 @@ subroutine test_init_distribs_container(success) if(i.eq.2) exit ! Test case 2: Custom initialisation + history_len = 15 nbins = [10, 20, 30] width = [0.05_real32, 0.1_real32, 0.15_real32] sigma = [0.2_real32, 0.3_real32, 0.4_real32] cutoff_min = [1.0_real32, 2.0_real32, 3.0_real32] cutoff_max = [5.0_real32, 6.0_real32, 7.0_real32] distribs_container = distribs_container_type( & + history_len=history_len, & nbins=nbins, & width=width, & sigma=sigma, & @@ -183,9 +202,63 @@ subroutine test_init_distribs_container(success) cutoff_max=cutoff_max & ) write(*,*) "Handled error: cutoff_min > cutoff_max" - + end subroutine test_init_distribs_container + subroutine test_set_history_len(success) + implicit none + logical, intent(inout) :: success + + type(distribs_container_type) :: distribs_container + integer :: history_len + + ! Initialise test data + history_len = 10 + + ! Call the subroutine to set the history length + 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", & + success & + ) + + call assert( & + allocated(distribs_container%history_deltas) .and. & + size(distribs_container%history_deltas, 1) .eq. history_len, & + "History delta list was not allocated", & + success & + ) + + end subroutine test_set_history_len + + subroutine test_is_converged(success) + implicit none + logical, intent(inout) :: success + + type(distribs_container_type) :: distribs_container + + ! Call the subroutine to check if the system is converged + call distribs_container%set_history_len(10) + + ! Check if the system is converged + call assert( & + .not. distribs_container%is_converged( threshold = 1.E-2_real32 ), & + "System should not be converged on initialisation", & + success & + ) + + distribs_container%history_deltas(:) = 1.E-3_real32 + call assert( & + distribs_container%is_converged( threshold = 1.E-2_real32 ), & + "System should be converged", & + success & + ) + + end subroutine test_is_converged + subroutine test_set_width(success) implicit none logical, intent(inout) :: success @@ -295,7 +368,8 @@ subroutine test_set_radius_distance_tol(success) call assert( & all( & abs( & - distribs_container%radius_distance_tol - radius_distance_tol & + distribs_container%radius_distance_tol - & + radius_distance_tol & ) .lt. 1.E-6_real32 & ), & "Radius_distance_tol was not set correctly", & @@ -509,7 +583,7 @@ subroutine test_create(basis, success) "system not correctly deallocated", & success & ) - + end subroutine test_create subroutine test_update(basis, success) @@ -577,7 +651,7 @@ subroutine test_update(basis, success) "system not correctly deallocated", & success & ) - + end subroutine test_update subroutine test_add(basis, success) @@ -690,7 +764,7 @@ subroutine test_get_element_energies(basis, success) "Element energy is incorrect", & success & ) - + end subroutine test_get_element_energies subroutine test_get_element_energies_staticmem(basis, success) @@ -729,7 +803,7 @@ subroutine test_get_element_energies_staticmem(basis, success) "Element energy is incorrect", & success & ) - + end subroutine test_get_element_energies_staticmem subroutine test_set_bond_radii(basis, success) @@ -792,8 +866,8 @@ subroutine test_set_bond_radii(basis, success) "Bond radius was not set correctly", & success & ) - - + + end subroutine test_set_bond_radii @@ -806,11 +880,12 @@ subroutine test_get_bin(success) type(distribs_container_type) :: distribs_container distribs_container%nbins(1) = 10 + call distribs_container%set_num_bins() ! Check lower bound correct handling bin = distribs_container%get_bin(0._real32, 1) call assert( & - bin .eq. 0, & + bin .eq. 1, & "Bin is incorrect", & success & ) @@ -818,7 +893,7 @@ subroutine test_get_bin(success) ! Check upper bound correct handling bin = distribs_container%get_bin(100._real32, 1) call assert( & - bin .eq. 0, & + bin .eq. distribs_container%nbins(1), & "Bin is incorrect", & success & ) @@ -847,4 +922,4 @@ subroutine assert(condition, message, success) end if end subroutine assert -end program test_distribs_container \ No newline at end of file +end program test_distribs_container diff --git a/test/test_evaluator_BTO.f90 b/test/test_evaluator_BTO.f90 index 3d07760..c4c9dc9 100644 --- a/test/test_evaluator_BTO.f90 +++ b/test/test_evaluator_BTO.f90 @@ -171,10 +171,8 @@ program test_evaluator_BTO atom_ignore_list(:,4) = [3, 6] atom_ignore_list(:,5) = [2, 2] - call basis_host%create_images( & - max_bondlength = max_bondlength, & - atom_ignore_list = atom_ignore_list & - ) + call basis_host%set_atom_mask( atom_ignore_list ) + call basis_host%create_images( max_bondlength = max_bondlength ) @@ -211,7 +209,6 @@ program test_evaluator_BTO basis_host, & [ 1, 2, 3 ], & [ generator%distributions%bond_info(:)%radius_covalent ], & - atom_ignore_list, & grid_offset = generator%grid_offset & ) do i = 1, 3 @@ -236,9 +233,7 @@ program test_evaluator_BTO write(unit,fmt) basis_host%spec(:)%name do is = 1, basis_host%nspec atom_loop: do ia = 1, basis_host%spec(is)%num - do i = 1, size(atom_ignore_list,2) - if( all(atom_ignore_list(:,i).eq.[is,ia]) ) cycle atom_loop - end do + if(.not.basis_host%spec(is)%atom_mask(ia)) cycle atom_loop write(unit,*) basis_host%spec(is)%atom(ia,:3) end do atom_loop end do @@ -264,7 +259,6 @@ program test_evaluator_BTO viability_grid(is,i) = evaluate_point( generator%distributions, & gridpoints(1:3,i), is, & basis_host, & - atom_ignore_list(:,ia:), & [ generator%distributions%bond_info(:)%radius_covalent ] & ) end do @@ -298,6 +292,7 @@ program test_evaluator_BTO viability_grid(atom_ignore_list(1,ia),best_loc) write(*,*) "Gridpoint: ", gridpoints(1:3,best_loc) end if + call basis_host%set_atom_mask( atom_ignore_list(:,ia:ia) ) call basis_host%update_images( & max_bondlength = max_bondlength, & is = atom_ignore_list(1,ia), ia = atom_ignore_list(2,ia) & diff --git a/test/test_evaluator_C.f90 b/test/test_evaluator_C.f90 index 2b0f6e0..f619ab4 100644 --- a/test/test_evaluator_C.f90 +++ b/test/test_evaluator_C.f90 @@ -156,10 +156,8 @@ program test_evaluator_C basis_host%spec(1)%num - size(atom_ignore_list,2) + i end do - call basis_host%create_images( & - max_bondlength = max_bondlength, & - atom_ignore_list = atom_ignore_list & - ) + call basis_host%set_atom_mask( atom_ignore_list ) + call basis_host%create_images( max_bondlength = max_bondlength) generator%distributions%kBT = 0.2 @@ -195,7 +193,6 @@ program test_evaluator_C basis_host, & [ 1 ], & [ generator%distributions%bond_info(:)%radius_covalent ], & - atom_ignore_list, & grid_offset = generator%grid_offset & ) do i = 1, 3 @@ -220,9 +217,7 @@ program test_evaluator_C write(unit,fmt) basis_host%spec(:)%name do is = 1, basis_host%nspec atom_loop: do ia = 1, basis_host%spec(is)%num - do i = 1, size(atom_ignore_list,2) - if( all(atom_ignore_list(:,i).eq.[is,ia]) ) cycle atom_loop - end do + if(.not.basis_host%spec(is)%atom_mask(ia)) cycle atom_loop write(unit,*) basis_host%spec(is)%atom(ia,:3) end do atom_loop end do @@ -246,7 +241,6 @@ program test_evaluator_C do i = 1, size(gridpoints,dim=2) viability_grid(1,i) = evaluate_point( generator%distributions, & gridpoints(1:3,i), atom_ignore_list(1,ia), basis_host, & - atom_ignore_list(:,ia:), & [ generator%distributions%bond_info(:)%radius_covalent ] & ) end do @@ -268,6 +262,7 @@ program test_evaluator_C "Incorrect gridpoint found.", & success & ) + call basis_host%set_atom_mask( atom_ignore_list(:,ia:ia) ) call basis_host%update_images( & max_bondlength = max_bondlength, & is = 1, ia = atom_ignore_list(2,ia) & diff --git a/test/test_generator.f90 b/test/test_generator.f90 index 3ab0da9..34bcde9 100644 --- a/test/test_generator.f90 +++ b/test/test_generator.f90 @@ -12,7 +12,7 @@ program test_generator character(len=100) :: filename type(raffle_generator_type) :: generator, generator2 class(raffle_generator_type), allocatable :: generator_var - type(basis_type) :: basis_host, basis_host_expected + type(basis_type) :: basis_host, basis_expected type(basis_type), dimension(1) :: database type(basis_type), dimension(:), allocatable :: structures, structures_store integer, dimension(3) :: grid @@ -144,34 +144,34 @@ program test_generator !----------------------------------------------------------------------------- - ! set up expected generated host structure - !----------------------------------------------------------------------------- - basis_host_expected%sysname = 'diamond+inserts' - basis_host_expected%nspec = 1 - allocate(basis_host_expected%spec(basis_host_expected%nspec)) - basis_host_expected%spec(1)%num = 16 - basis_host_expected%spec(1)%name = 'C' - basis_host_expected%natom = sum(basis_host_expected%spec(:)%num) - allocate(basis_host_expected%spec(1)%atom(basis_host_expected%spec(1)%num, 3)) - basis_host_expected%spec(1)%atom(1, :3) = [0.0, 0.0, 0.0] - basis_host_expected%spec(1)%atom(2, :3) = [0.5, 0.5, 0.0] - basis_host_expected%spec(1)%atom(3, :3) = [0.5, 0.0, 0.25] - basis_host_expected%spec(1)%atom(4, :3) = [0.0, 0.5, 0.25] - basis_host_expected%spec(1)%atom(5, :3) = [0.25, 0.25, 0.125] - basis_host_expected%spec(1)%atom(6, :3) = [0.75, 0.75, 0.125] - basis_host_expected%spec(1)%atom(7, :3) = [0.75, 0.25, 0.375] - basis_host_expected%spec(1)%atom(8, :3) = [0.25, 0.75, 0.375] - basis_host_expected%spec(1)%atom(9, :3) = [0.75, 0.25, 0.875] - basis_host_expected%spec(1)%atom(10, :3) = [0.75, 0.75, 0.625] - basis_host_expected%spec(1)%atom(11, :3) = [0.5, 0.0, 0.75] - basis_host_expected%spec(1)%atom(12, :3) = [0.25, 0.25, 0.625] - basis_host_expected%spec(1)%atom(13, :3) = [0.25, 0.75, 0.875] - basis_host_expected%spec(1)%atom(14, :3) = [0.5, 0.5, 0.5] - basis_host_expected%spec(1)%atom(15, :3) = [0.0, 0.5, 0.75] - basis_host_expected%spec(1)%atom(16, :3) = [0.0, 0.0, 0.5] - basis_host_expected%lat(1,:) = [3.560745109, 0.0, 0.0] - basis_host_expected%lat(2,:) = [0.0, 3.560745109, 0.0] - basis_host_expected%lat(3,:) = [0.0, 0.0, 7.121490218] + ! set up expected generated structure + !----------------------------------------------------------------------------- + basis_expected%sysname = 'diamond+inserts' + basis_expected%nspec = 1 + allocate(basis_expected%spec(basis_expected%nspec)) + basis_expected%spec(1)%num = 16 + basis_expected%spec(1)%name = 'C' + basis_expected%natom = sum(basis_expected%spec(:)%num) + allocate(basis_expected%spec(1)%atom(basis_expected%spec(1)%num, 3)) + basis_expected%spec(1)%atom(1, :3) = [0.0, 0.0, 0.0] + basis_expected%spec(1)%atom(2, :3) = [0.5, 0.5, 0.0] + basis_expected%spec(1)%atom(3, :3) = [0.5, 0.0, 0.25] + basis_expected%spec(1)%atom(4, :3) = [0.0, 0.5, 0.25] + basis_expected%spec(1)%atom(5, :3) = [0.25, 0.25, 0.125] + basis_expected%spec(1)%atom(6, :3) = [0.75, 0.75, 0.125] + basis_expected%spec(1)%atom(7, :3) = [0.75, 0.25, 0.375] + basis_expected%spec(1)%atom(8, :3) = [0.25, 0.75, 0.375] + basis_expected%spec(1)%atom(9, :3) = [0.75, 0.25, 0.875] + basis_expected%spec(1)%atom(10, :3) = [0.75, 0.75, 0.625] + basis_expected%spec(1)%atom(11, :3) = [0.5, 0.0, 0.75] + basis_expected%spec(1)%atom(12, :3) = [0.25, 0.25, 0.625] + basis_expected%spec(1)%atom(13, :3) = [0.25, 0.75, 0.875] + basis_expected%spec(1)%atom(14, :3) = [0.5, 0.5, 0.5] + basis_expected%spec(1)%atom(15, :3) = [0.0, 0.5, 0.75] + basis_expected%spec(1)%atom(16, :3) = [0.0, 0.0, 0.5] + basis_expected%lat(1,:) = [3.560745109, 0.0, 0.0] + basis_expected%lat(2,:) = [0.0, 3.560745109, 0.0] + basis_expected%lat(3,:) = [0.0, 0.0, 7.121490218] !----------------------------------------------------------------------------- @@ -293,9 +293,9 @@ program test_generator grid(2) = nint( 0.75 * generator%host%lat(2,2) / 0.2 ) grid(3) = nint( 0.5 * generator%host%lat(3,3) / 0.2 ) call assert( & - all( generator%grid .eq. grid ), & - 'Generator failed to handle grid_spacing with bounds', & - success & + all( generator%grid .eq. grid ), & + 'Generator failed to handle grid_spacing with bounds', & + success & ) call generator%reset_bounds() @@ -306,7 +306,7 @@ program test_generator ! set up generator !----------------------------------------------------------------------------- generator%distributions%kBT = 0.2 - call generator%set_grid( grid_spacing = 0.2, grid_offset = [0.0, 0.0, 0.0] ) + call generator%set_grid( grid_spacing = 0.15, grid_offset = [0.0, 0.0, 0.0] ) generator%distributions%radius_distance_tol = [1.5, 2.5, 3.0, 6.0] do i = 1, 3 tolerance(i) = 1._real32 / real(generator%grid(i),real32) / 2._real32 @@ -341,11 +341,11 @@ program test_generator !----------------------------------------------------------------------------- - ! compare the generated host structure + ! compare the generated structure !----------------------------------------------------------------------------- call assert(& - compare_bas(generator%structures(1), basis_host_expected, tolerance), & - 'Generated host structure does not match expected host structure', & + compare_bas(generator%structures(1), basis_expected, tolerance), & + 'Generated structure does not match expected structure', & success & ) diff --git a/test/test_geom_extd.f90 b/test/test_geom_extd.f90 index f3e4319..5e7fab0 100644 --- a/test/test_geom_extd.f90 +++ b/test/test_geom_extd.f90 @@ -92,12 +92,10 @@ subroutine test_update_iamges(basis, success) call basis_copy%copy(basis) atom_ignore_list(:,1) = [1, 1] + call basis_copy%set_atom_mask( atom_ignore_list ) ! Create images - call basis_copy%create_images( & - max_bondlength = 0._real32, & - atom_ignore_list = atom_ignore_list & - ) + call basis_copy%create_images( max_bondlength = 0._real32 ) ! Check if the number of images is correct call assert( & diff --git a/test/test_geom_utils.f90 b/test/test_geom_utils.f90 index c9b9cff..6e6bb2d 100644 --- a/test/test_geom_utils.f90 +++ b/test/test_geom_utils.f90 @@ -53,7 +53,7 @@ program test_geom_utils bas2%lat(2,:) = [2.14, 0.0, 2.14] bas2%lat(3,:) = [2.14, 2.14, 0.0] - + basis_merged = basis_merge(bas, bas2) if ( basis_merged%nspec .ne. 2 ) then @@ -92,7 +92,44 @@ program test_geom_utils basis_merged%spec(2)%name success = .false. end if + if( any(.not.basis_merged%spec(1)%atom_mask) ) then + write(0,*) & + 'basis_merge failed, atom mask for species 1 not all true' + success = .false. + end if + if( any(.not.basis_merged%spec(2)%atom_mask) ) then + write(0,*) & + 'basis_merge failed, atom mask for species 2 not all true' + success = .false. + end if + + basis_merged = basis_merge(bas, bas2, mask1 = .true., mask2 = .false.) + if( any(.not.basis_merged%spec(1)%atom_mask(:bas%spec(1)%num)) ) then + write(0,*) & + 'basis_merge failed, atom mask for basis 1 not all true' + success = .false. + end if + if( any(basis_merged%spec(1)%atom_mask(bas%spec(1)%num+1:)) .or. & + any(basis_merged%spec(2)%atom_mask) & + ) then + write(0,*) & + 'basis_merge failed, atom mask for basis 2 not all false' + success = .false. + end if + basis_merged = basis_merge(bas, bas2, mask1 = .false., mask2 = .true.) + if( any(basis_merged%spec(1)%atom_mask(:bas%spec(1)%num)) ) then + write(0,*) & + 'basis_merge failed, atom mask for basis 1 not all false' + success = .false. + end if + if( any(.not.basis_merged%spec(1)%atom_mask(bas%spec(1)%num+1:)) .or. & + any(.not.basis_merged%spec(2)%atom_mask) & + ) then + write(0,*) & + 'basis_merge failed, atom mask for basis 2 not all true' + success = .false. + end if !----------------------------------------------------------------------------- ! check for any failed tests @@ -105,4 +142,4 @@ program test_geom_utils stop 1 end if -end program test_geom_utils \ No newline at end of file +end program test_geom_utils diff --git a/test/test_place_methods.f90 b/test/test_place_methods.f90 index 1098752..b67ea1d 100644 --- a/test/test_place_methods.f90 +++ b/test/test_place_methods.f90 @@ -97,10 +97,8 @@ program test_place_methods seed_arr = seed call random_seed(put=seed_arr) call basis_extd%copy(basis) - call basis_extd%create_images( & - max_bondlength = 6._real32, & - atom_ignore_list = atom_ignore_list & - ) + call basis_extd%set_atom_mask( atom_ignore_list ) + call basis_extd%create_images( max_bondlength = 6._real32 ) !----------------------------------------------------------------------------- @@ -117,7 +115,7 @@ program test_place_methods generator%distributions, & bounds, & basis_extd, & - atom_ignore_list, & + atom_ignore_list(1,1), & radius_list = [ 0.5_real32 ], & max_attempts = 1000, & viable = viable & @@ -140,7 +138,7 @@ program test_place_methods generator%distributions, & bounds, & basis_extd, & - atom_ignore_list, & + atom_ignore_list(1,1), & radius_list = [ 0.5_real32 ], & max_attempts = 1000, & step_size_fine = 0.1_real32, & @@ -167,7 +165,7 @@ program test_place_methods prior_species = 1, & bounds = bounds, & basis = basis_extd, & - atom_ignore_list = atom_ignore_list, & + species = atom_ignore_list(1,1), & radius_list = [ 0.5_real32 ], & max_attempts = 1000, & step_size_fine = 0.1_real32, & @@ -227,10 +225,8 @@ subroutine test_place_method_void(distributions, basis, success) ! Initialise basis call basis_copy%copy(basis) - call basis_copy%create_images( & - max_bondlength = 6._real32, & - atom_ignore_list = atom_ignore_list & - ) + call basis_copy%set_atom_mask( atom_ignore_list ) + call basis_copy%create_images( max_bondlength = 6._real32) !--------------------------------------------------------------------------- ! set up gridpoints @@ -242,16 +238,11 @@ subroutine test_place_method_void(distributions, basis, success) basis_copy, & [ 1 ], & [ distributions%bond_info(:)%radius_covalent ], & - atom_ignore_list, & grid_offset = grid_offset & ) ! Call the void subroutine - point = place_method_void( & - gridpoints, basis_copy, & - atom_ignore_list, & - viable & - ) + point = place_method_void( gridpoints, basis_copy, viable ) ! Check if viable call assert(viable, "No viable gridpoints found.", success) diff --git a/test/test_viability.f90 b/test/test_viability.f90 index 3a699cd..c8130c4 100644 --- a/test/test_viability.f90 +++ b/test/test_viability.f90 @@ -72,9 +72,9 @@ subroutine test_get_gridpoints_and_viability(basis, success) ! Initialise basis call basis_copy%copy(basis) + call basis_copy%set_atom_mask( atom_ignore_list ) call basis_copy%create_images( & - max_bondlength = distribs_container%cutoff_max(1), & - atom_ignore_list = atom_ignore_list & + max_bondlength = distribs_container%cutoff_max(1) & ) ! Initialise gvector container @@ -98,7 +98,6 @@ subroutine test_get_gridpoints_and_viability(basis, success) basis_copy, & [ 1 ], & radius_list, & - atom_ignore_list, & grid_offset & ) @@ -149,9 +148,9 @@ subroutine test_update_gridpoints_and_viability(basis, success) ! Initialise basis call basis_copy%copy(basis) + call basis_copy%set_atom_mask( atom_ignore_list ) call basis_copy%create_images( & - max_bondlength = distribs_container%cutoff_max(1), & - atom_ignore_list = atom_ignore_list & + max_bondlength = distribs_container%cutoff_max(1) & ) ! Initialise gvector container @@ -175,7 +174,6 @@ subroutine test_update_gridpoints_and_viability(basis, success) basis_copy, & [ 1 ], & radius_list, & - atom_ignore_list, & grid_offset & ) @@ -184,8 +182,7 @@ subroutine test_update_gridpoints_and_viability(basis, success) points, distribs_container, basis_copy, & [1], & [1,2], & - radius_list, & - atom_ignore_list & + radius_list & ) ! Check points exist @@ -211,8 +208,7 @@ subroutine test_update_gridpoints_and_viability(basis, success) points, distribs_container, basis_copy, & [1], & [1,2], & - radius_list, & - atom_ignore_list & + radius_list & ) ! Check all points have been removed