From d9b18a7c20ce023018755c202f8d06cbf8bd27c5 Mon Sep 17 00:00:00 2001 From: Marnik Bercx Date: Thu, 15 Jun 2023 13:12:53 +0200 Subject: [PATCH] Add the `create_magnetic_configuration` function Here we add the `create_magnetic_configuration` calculation function, which returns a new `StructureData` with the required kinds as well as a `Dict` node with the corresponding magnetic moments for each kind, based on the input structure and magnetic moments per site. To create the new list of kinds, the algorithm loops over all the elements in the structure and makes a list of the sites with that element and their corresponding magnetic moment. Next, it splits this list in three lists: * Zero magnetic moments: Any site that has an absolute magnetic moment lower than `ztol` * Positive magnetic moments * Negative magnetic moments The algorithm then sorts the positive and negative lists from large to small absolute value, and loops over each list. New magnetic kinds will be created when the absolute difference between the magnetic moment of the current kind and the site exceeds `atol`. The positive and negative magnetic moments are handled separately to avoid assigning two sites with opposite signs in their magnetic moment to the same kind and make sure that each kind has the correct magnetic moment, i.e. the largest magnetic moment in absolute value of the sites corresponding to that kind. To make it easier to pass the magnetic moments between calculations, a new tool is added to `PwCalculation`: `get_magnetic_configuration`. This will check the final magnetic moments from the calculation, create a new `StructureData` if necessary, and return the structure and dictionary that maps the kinds to the magnetic moments. A new tutorial is added on "Magnetic configurations" that explains how to work with spin-polarised calculations in the plugin. --- docs/Makefile | 1 + docs/source/_static/aiida-custom.css | 13 + docs/source/_static/aiida-qe-custom.css | 21 +- docs/source/_templates/icon-links.html | 2 +- docs/source/conf.py | 35 +- docs/source/howto/calculations/ph.md | 6 +- docs/source/howto/calculations/pw.md | 18 +- docs/source/tutorials/first_pw.md | 295 ++++++++++++++ docs/source/tutorials/index.md | 328 ++-------------- docs/source/tutorials/magnetism.md | 366 ++++++++++++++++++ pyproject.toml | 1 + .../create_magnetic_configuration.py | 123 ++++++ .../tools/calculations/pw.py | 46 ++- .../workflows/protocols/utils.py | 6 +- .../test_create_magnetic_configuration.py | 207 ++++++++++ tests/conftest.py | 23 ++ 16 files changed, 1178 insertions(+), 313 deletions(-) create mode 100644 docs/source/tutorials/first_pw.md create mode 100644 docs/source/tutorials/magnetism.md create mode 100644 src/aiida_quantumespresso/calculations/functions/create_magnetic_configuration.py create mode 100644 tests/calculations/functions/test_create_magnetic_configuration.py diff --git a/docs/Makefile b/docs/Makefile index 20f595cd5..6b0e93dfe 100755 --- a/docs/Makefile +++ b/docs/Makefile @@ -28,6 +28,7 @@ all: clean html view clean: rm -rf $(BUILDDIR) + rm -rf $(shell find . -type d -wholename "*/reference/api/auto") html: $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html diff --git a/docs/source/_static/aiida-custom.css b/docs/source/_static/aiida-custom.css index a3ebd1c63..42d2fa2dd 100644 --- a/docs/source/_static/aiida-custom.css +++ b/docs/source/_static/aiida-custom.css @@ -112,6 +112,19 @@ div.admonition :last-child { div.highlight-bash div.highlight { background-color: aliceblue; } + div.highlight-console div.highlight { background-color: aliceblue; } + +.aiida-green { + color: #32B805; +} + +.aiida-blue { + color: #0496DE; +} + +.aiida-orange { + color: #FF7D16; +} diff --git a/docs/source/_static/aiida-qe-custom.css b/docs/source/_static/aiida-qe-custom.css index fc02e48e1..8ebb4c265 100644 --- a/docs/source/_static/aiida-qe-custom.css +++ b/docs/source/_static/aiida-qe-custom.css @@ -26,10 +26,15 @@ img.logo__image { font-size: 40px!important; } -img.aiida-logo { +img.aiida-sidebar-logo { height: 40px!important; } +img.aiida-logo { + width: 20px; + padding-bottom: 3px; +} + .fa { color: var(--pst-color-primary); } @@ -37,3 +42,17 @@ img.aiida-logo { .bd-search { padding: 0 1rem; } + +.tutor-footer { + padding-top: 0rem; + border-top: none!important; +} + +.footer-table { + margin-bottom: 0rem; + border-color: transparent; +} + +.footer-table td:last-child { + text-align: right; +} diff --git a/docs/source/_templates/icon-links.html b/docs/source/_templates/icon-links.html index 8ec32e8c0..356c4f733 100644 --- a/docs/source/_templates/icon-links.html +++ b/docs/source/_templates/icon-links.html @@ -30,7 +30,7 @@ aria-label="{{ theme_icon_links_label }}">
  • - +
    • {%- for icon_link in theme_icon_links -%} diff --git a/docs/source/conf.py b/docs/source/conf.py index b93e99c2b..7c3445c45 100755 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -30,9 +30,6 @@ # -- General configuration ------------------------------------------------ -# If your documentation needs a minimal Sphinx version, state it here. -#needs_sphinx = '1.0' - # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. @@ -55,12 +52,6 @@ 'aiida_pseudo': ('http://aiida-pseudo.readthedocs.io/en/latest/', None), } -myst_enable_extensions = [ - 'deflist', - 'colon_fence', - 'substitution' -] - # Settings for the `autoapi.extenstion` automatically generating API docs filepath_docs = pathlib.Path(__file__).parent.parent filepath_src = filepath_docs.parent / 'src' @@ -90,6 +81,32 @@ # directories to ignore when looking for source files. exclude_patterns = ['**.ipynb_checkpoints', 'reference/api/auto/aiida_quantumespresso/index.rst'] +# -- MyST options + +myst_enable_extensions = [ + 'deflist', + 'colon_fence', + 'substitution', + 'attrs_inline', + 'substitution' +] + +myst_substitutions = { + 'aiida_logo': '', + 'create_magnetic_configuration': \ + '{func}`~aiida_quantumespresso.calculations.functions.create_magnetic_configuration.create_magnetic_configuration`', + 'get_builder_from_protocol': \ + '{meth}`~aiida_quantumespresso.workflows.pw.base.PwBaseWorkChain.get_builder_from_protocol`', + 'get_magnetic_configuration': \ + '{meth}`~aiida_quantumespresso.tools.calculations.pw.PwCalculationTools.get_magnetic_configuration`', + 'nspin': '[`nspin`](https://www.quantum-espresso.org/Doc/INPUT_PW.html#idm412)', + 'PwBaseWorkChain': '{class}`~aiida_quantumespresso.workflows.pw.base.PwBaseWorkChain`', + 'PwCalculation': '{class}`~aiida_quantumespresso.calculations.pw.PwCalculation`', + 'SpinType': '{class}`~aiida_quantumespresso.common.types.SpinType`', + 'starting_magnetization': '[`starting_magnetization`](https://www.quantum-espresso.org/Doc/INPUT_PW.html#idm299)', + 'StructureData': '{class}`~aiida.orm.StructureData', +} + # -- Options for HTML output ---------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for diff --git a/docs/source/howto/calculations/ph.md b/docs/source/howto/calculations/ph.md index a601f8ffd..3bedbf106 100644 --- a/docs/source/howto/calculations/ph.md +++ b/docs/source/howto/calculations/ph.md @@ -16,18 +16,18 @@ In order to run a `ph.x` calculation, you first need to have completed a `pw.x` See the [tutorial](#tutorials-pw-through-api) or [how-to guide](#howto-calculations-pw) for more information. ::: -Once you have successfully run a `PwCalculation` you can run a `ph.x` calculation through the `PhCalculation` plugin as follows: +Once you have successfully run a {{ PwCalculation }} you can run a `ph.x` calculation through the `PhCalculation` plugin as follows: ```{literalinclude} ../../tutorials/include/scripts/run_ph_basic.py :language: python ``` -Note that you will have to replace `IDENTIFIER_PW_CALCULATION` with the identifier (pk or UUID) of the completed `PwCalculation`. +Note that you will have to replace `IDENTIFIER_PW_CALCULATION` with the identifier (pk or UUID) of the completed {{ PwCalculation }}. ## How to define input file parameters The `ph.x` code supports many parameters that can be defined through the input file, as shown on the [official documentation](https://www.quantum-espresso.org/Doc/INPUT_PH.html). -Parameters that are part of the `INPUTPH` card should be specified through the `parameters` input of the `PwCalculation` plugin. +Parameters that are part of the `INPUTPH` card should be specified through the `parameters` input of the {{ PwCalculation }} plugin. The parameters are specified using a Python dictionary, for example: ```python diff --git a/docs/source/howto/calculations/pw.md b/docs/source/howto/calculations/pw.md index cef883eee..1947a93fe 100644 --- a/docs/source/howto/calculations/pw.md +++ b/docs/source/howto/calculations/pw.md @@ -18,7 +18,7 @@ Examples of these properties include ground-state energy and one-electron (Kohn- ## How to launch a `pw.x` calculation -Below is a script with a basic example of how to run a `pw.x` calculation through the `PwCalculation` plugin that computes the electronic ground state of an fcc silicon crystal: +Below is a script with a basic example of how to run a `pw.x` calculation through the {{ PwCalculation }} plugin that computes the electronic ground state of an fcc silicon crystal: ```{literalinclude} ../../tutorials/include/scripts/run_pw_basic.py :language: python @@ -30,7 +30,7 @@ Note that you may have to change the name of the code that is loaded using `load The `pw.x` code supports many parameters that can be defined through the input file, as shown on the [official documentation](https://www.quantum-espresso.org/Doc/INPUT_PW.html). The parameters are divided into section or "cards". -Parameters that are part of cards that start with an ampersand (`&`) should be specified through the `parameters` input of the `PwCalculation` plugin. +Parameters that are part of cards that start with an ampersand (`&`) should be specified through the `parameters` input of the {{ PwCalculation }} plugin. The parameters are specified using a Python dictionary, where each card is its own sub-dictionary, for example: ```python @@ -78,13 +78,15 @@ These include: Defining them anyway will result in an exception when launching the calculation. ::: +(howto-calculations-pw-multidimensional-parameters)= + ### Multidimensional parameters The input format of `pw.x` contains various keywords that do not simply take the format of a key value pair, but rather there will some indices in the key itself. Take for example the [`starting_magnetization`](https://www.quantum-espresso.org/Doc/INPUT_PW.html#idm287) keyword of the `SYSTEM` card. The starting magnetization value needs to be applied to a specific species and therefore the index `i` is required to be able to make this distinction. -The `PwCalculation` plugin makes this easy as it will do the conversion from kind name to species index automatically. +The {{ PwCalculation }} plugin makes this easy as it will do the conversion from kind name to species index automatically. This allows you to specify a starting magnetization value by using a dictionary notation, where the key is the kind name to which it should be applied. For example, if you have a structure with the kind `Co` and want it to have a given starting magnetization, one can add the following in the parameter data dictionary: @@ -171,7 +173,7 @@ Assuming the input structure contained the kinds `Ni` and `Fe`, which would have ## How to define pseudopotentials Each `pw.x` calculation requires a pseudopotential to be specified for each kind in the structure. -These pseudopotentials can be specified in the `PwCalculation` plugin through the `pseudos` input namespace. +These pseudopotentials can be specified in the {{ PwCalculation }} plugin through the `pseudos` input namespace. This input takes a dictionary, where the keys are the kind names and the values are instances of the {class}`~aiida_pseudo.data.pseudo.upf.UpfData` data plugin of the {{ aiida_pseudo }} plugin package. For example, if the input structure is a `GaAs` crystal, the pseudopotentials could be specified as follows: @@ -263,7 +265,7 @@ builder.settings = Dict({'GAMMA_ONLY': True}) ## How to fix the coordinates of atoms Quantum ESPRESSO pw.x allows to optionally fix the coordinates of atoms during relaxation and molecular-dynamics simulations. -This functionality is enabled in `PwCalculation` through the `FIXED_COORDS` setting which is a list of length equal to the number of sites in the input structure. +This functionality is enabled in {{ PwCalculation }} through the `FIXED_COORDS` setting which is a list of length equal to the number of sites in the input structure. Each element is a list of length three containing booleans, where `True` means that the position of the site in that direction should be fixed. For example: @@ -277,7 +279,7 @@ All other coordinates are allowed to change. ## How to retrieve additional files -The `PwCalculation` plugin will retrieve the most important and common output files by default. +The {{ PwCalculation }} plugin will retrieve the most important and common output files by default. To retrieve additional output files, specify the list of files in the `ADDITIONAL_RETRIEVE_LIST` key in the `settings` input: ```python @@ -288,13 +290,13 @@ builder.settings = Dict({'ADDITIONAL_RETRIEVE_LIST': ['custom-file.txt', 'some-o ## How to analyze the results -When a `PwCalculation` is completed, there are quite a few possible analyses to perform. +When a {{ PwCalculation }} is completed, there are quite a few possible analyses to perform. ### How to check the SCF accuracy during the self-consistent cycle During the self-consistent field cycle, the difference in energy of the newly computed charge density and the starting one is computed and stored. It can easily be retrieved through the {meth}`~aiida_quantumespresso.tools.calculations.pw.PwCalculationTools.get_scf_accuracy` method. -This method can be accessed directly through the ``tools`` of a completed ``PwCalculation`` node: +This method can be accessed directly through the `tools` of a completed {{ PwCalculation }} node: ```python In [1]: node = load_node(IDENTIFIER) diff --git a/docs/source/tutorials/first_pw.md b/docs/source/tutorials/first_pw.md new file mode 100644 index 000000000..6d8657d3f --- /dev/null +++ b/docs/source/tutorials/first_pw.md @@ -0,0 +1,295 @@ + +(tutorials-running-pw)= + +# Running a `pw.x` calculation + +(tutorials-pw-through-cli)= + +## Running a `pw.x` calculation through the CLI + +The simplest way to run a first `pw.x` calculation with `aiida-quantumespresso` is through the [command line interface (CLI)](#reference-cli) that ships with the package. + +```console +aiida-quantumespresso calculation launch pw -X pw@locahost -F SSSP/1.1/PBE/efficiency +``` + +If everything is setup correctly, you should output that looks something like the following: + +```console +Running a PwCalculation... +PwCalculation<59517> terminated with state: finished [0] + +Output link Node pk and type +------------------------------------------------------------ +output_band BandsData<59520> +output_parameters Dict<59522> +output_trajectory TrajectoryData<59521> +remote_folder RemoteData<59518> +retrieved FolderData<59519> +``` + +The command ran a {{ PwCalculation }}, with the identifier `59517` and it finished correctly. + +:::{note} +We can know that the job finished correctly because it says its state is `finished [0]`. +The `0` is the exit status of the job, which means it was successful. +A non-zero exit code means the primary objective of the calculation was not (fully) achieved or an error occurred. +::: + +The outputs that were generated by the calculation are printed to screen. +To obtain more information about the calculation that ran, such as its inputs, run `verdi process show ` + +```console +verdi process show 59517 + +Property Value +----------- ------------------------------------ +type PwCalculation +state Finished [0] +pk 59517 +uuid 982e6b36-bdf6-4525-8fd3-86e546afcd7f +label +description +ctime 2022-10-09 16:39:49.026195+02:00 +mtime 2022-10-09 16:39:51.543732+02:00 +computer [1] localhost + +Inputs PK Type +---------- ----- ------------- +pseudos + Si 56785 UpfData +code 56729 Code +kpoints 59515 KpointsData +parameters 59516 Dict +structure 56731 StructureData + +Outputs PK Type +----------------- ----- -------------- +output_band 59520 BandsData +output_parameters 59522 Dict +output_trajectory 59521 TrajectoryData +remote_folder 59518 RemoteData +retrieved 59519 FolderData +``` + +The calculation received a number of inputs: `pseudos`, `code`, `kpoints`, `parameters` and a `structure`. +The `code` was specified through the `-X` option when we invoked the CLI command to launch the calculation. +Other inputs, such as `structure`, `parameters` etc. were created automatically by the CLI command. + +The inputs can be introspected through the `verdi` CLI that ships with `aiida-core`. +For example, the structure can be visualized using the command: + +```console +verdi data core.structure show 56731 +``` + +This should open a new dialog with a visualization of the silicon crystal structure that was automatically created for this calculation. +Similarly, the `parameters` can be shown using: + +```console +verdi node attributes 59516 + +PK: 59516 +{ + "CONTROL": { + "calculation": "scf" + }, + "SYSTEM": { + "ecutrho": 240.0, + "ecutwfc": 30.0 + } +} +``` + +Likewise, the output parameters (`output_parameters`), which are of the same type (`Dict`) of the input parameters (`parameters`) can also be shown: + +```console +verdi node attributes 59522 + +PK: 59522 +{ + "beta_real_space": false, + "convergence_info": { + "scf_conv": { + "convergence_achieved": true, + "n_scf_steps": 6, + "scf_error": 2.1572731407236e-07 + } + }, + "creator_name": "pwscf", + "creator_version": "6.6", + "energy": -308.19211637089, + ... + "energy_accuracy": 5.850447441879e-06, + "energy_accuracy_units": "eV", + "fermi_energy": 6.508401516556, + "fermi_energy_units": "eV", + "number_of_atoms": 2, + "number_of_bands": 4, + "number_of_bravais_symmetries": 48, + "number_of_electrons": 8.0, + "number_of_k_points": 3, + "number_of_species": 1, + "number_of_spin_components": 1, + "number_of_symmetries": 48, + "scf_iterations": 6, + "volume": 40.02575122515, + "wall_time_seconds": 1.23, +} +``` + +The output will show a number of parameters that have been automatically parsed from the output of the `pw.x` calculation. +Note that the output has been shortened for clarity and not all parsed output parameters are shown. + +(tutorials-pw-through-api)= + +## Running a `pw.x` calculation through the API + +The following minimal example shows how to run a `pw.x` calculation through AiiDA's Python API. +For the purposes of this demonstration, the electronic ground state of an fcc silicon crystal structure is computed. + +:::{tip} +The code that is shown below in snippets can be {download}`downloaded as a script `. +Save the script to a file, make it executable and then simply execute it to run the example calculation. +::: + +First, import the required classes and functions: + +```python +from aiida.engine import run +from aiida.orm import Dict, KpointsData, StructureData, load_code, load_group +``` + +Then, load the code that was setup in AiiDA for `pw.x` and get an instance of the [process builder](https://aiida.readthedocs.io/projects/aiida-core/en/latest/topics/processes/usage.html#process-builder): + +```python +# Load the code configured for ``pw.x``. Make sure to replace this string +# with the label that you used in the setup of the code. +code = load_code('pw@localhost') +builder = code.get_builder() +``` + +The process builder can be used to assign the inputs that are to be used for the calculation. +When an input is assigned, it is automatically validated. +Let's first define a structure and assign it to the builder: + +```python +# Create a silicon fcc crystal +from ase.build import bulk +structure = StructureData(ase=bulk('Si', 'fcc', 5.43)) +builder.structure = structure +``` + +To construct a crystal structure, we use the `bulk` method of the [ASE library](https://wiki.fysik.dtu.dk/ase/ase/build/build.html#ase.build.bulk). +This is then wrapped in a `StructureData` node such that it can be stored in AiiDA's provenance graph, and assigned to the builder. + +Next, we need to define the pseudopotentials. +If the pseudopotentials were installed using the `aiida-pseudo` package, as described in the [installation guide](installation-setup-pseudopotentials), that can be done as follows: + +```python +# Load the pseudopotential family. +pseudo_family = load_group('SSSP/1.1/PBE/efficiency') +builder.pseudos = pseudo_family.get_pseudos(structure=structure) +``` + +The parameters of the input file that will be read by `pw.x` are defined through the `parameters` input. +A minimal setup to run an SCF calculation is shown below: + +```python +# Request the recommended wavefunction and charge density cutoffs +# for the given structure and energy units. +cutoff_wfc, cutoff_rho = pseudo_family.get_recommended_cutoffs( + structure=structure, + unit='Ry' +) + +parameters = Dict({ + 'CONTROL': { + 'calculation': 'scf' + }, + 'SYSTEM': { + 'ecutwfc': cutoff_wfc, + 'ecutrho': cutoff_rho, + } +}) +builder.parameters = parameters +``` + +:::{note} +The `pw.x` code actually requires more input parameters than have been specified in the example above. +The other required parameters are actually automatically defined by the {{ PwCalculation }} plugin, based on for example the input structure. +To find out what other parameters can be specified, please refer to the [documentation of Quantum ESPRESSO](https://www.quantum-espresso.org/Doc/INPUT_PW.html). +Note that certain parameters cannot be manually specified as these are reserved for the plugin to be set. +::: + +The k-points to be used for the calculation are defined using the `KpointsData` class: + +```python +# Generate a 2x2x2 Monkhorst-Pack mesh +kpoints = KpointsData() +kpoints.set_kpoints_mesh([2, 2, 2]) +builder.kpoints = kpoints +``` + +Finally, specify the amount of CPUs that the calculation should use and how long it can run before it should be killed: + +```python +# Run the calculation on 1 CPU and kill it if it runs longer than 1800 seconds. +# Set ``withmpi`` to ``False`` if ``pw.x`` was compiled without MPI support. +builder.metadata.options = { + 'resources': { + 'num_machines': 1, + }, + 'max_wallclock_seconds': 1800, + 'withmpi': False, +} +``` + +All required inputs have now been defined. +The calculation can be started by launching the process builder: + +```python +results, node = run.get_node(builder) +``` + +When the calculation is finished, `run.get_node` will return a tuple of two values: `results` and `node`. +The `results` will be a dictionary of all the output nodes that were produced and the `node` will hold a reference to the node that represents the calculation in AiiDA's provenance graph. +When the calculation is finished, check whether it was successful through its exit status: + +```python +node.exit_status +``` + +If the calculation finished successfully, this will return `0`. +A non-zero exit code means the primary objective of the calculation was not (fully) achieved or an error occurred. + +The `results` variable contains a dictionary of the created output nodes. +Print it to see exactly which outputs were created and what their identifiers are: + +```python +print(results) +{ + 'output_band': , + 'output_trajectory': , + 'output_parameters': , + 'remote_folder': , + 'retrieved': +} +``` + +To show the parsed output parameters, call: + +```python +results['output_parameters'].base.attributes.all +{ + 'energy': -153.30519007432, + 'volume': 40.02575122515, + .... +} +``` + +The complete output that was written by `pw.x` to stdout, can be retrieved as follows: + +```python +results['retrieved'].base.repository.get_object_content('aiida.out') +``` diff --git a/docs/source/tutorials/index.md b/docs/source/tutorials/index.md index 53520b770..8a44c7d66 100644 --- a/docs/source/tutorials/index.md +++ b/docs/source/tutorials/index.md @@ -10,293 +10,49 @@ Before you get started, make sure that you have: - installed the SSSP pseudopotential family ([see instructions](#installation-setup-pseudopotentials)) ::: -(tutorials-pw-through-cli)= - -## Running a `pw.x` calculation through the CLI - -The simplest way to run a first `pw.x` calculation with `aiida-quantumespresso` is through the [command line interface (CLI)](#reference-cli) that ships with the package. - -```console -aiida-quantumespresso calculation launch pw -X pw@locahost -F SSSP/1.1/PBE/efficiency -``` - -If everything is setup correctly, you should output that looks something like the following: - -```console -Running a PwCalculation... -PwCalculation<59517> terminated with state: finished [0] - -Output link Node pk and type ------------------------------------------------------------- -output_band BandsData<59520> -output_parameters Dict<59522> -output_trajectory TrajectoryData<59521> -remote_folder RemoteData<59518> -retrieved FolderData<59519> -``` - -The command ran a `PwCalculation`, with the identifier `59517` and it finished correctly. - -:::{note} -We can know that the job finished correctly because it says its state is `finished [0]`. -The `0` is the exit status of the job, which means it was successful. -A non-zero exit code means the primary objective of the calculation was not (fully) achieved or an error occurred. -::: - -The outputs that were generated by the calculation are printed to screen. -To obtain more information about the calculation that ran, such as its inputs, run `verdi process show ` - -```console -verdi process show 59517 - -Property Value ------------ ------------------------------------ -type PwCalculation -state Finished [0] -pk 59517 -uuid 982e6b36-bdf6-4525-8fd3-86e546afcd7f -label -description -ctime 2022-10-09 16:39:49.026195+02:00 -mtime 2022-10-09 16:39:51.543732+02:00 -computer [1] localhost - -Inputs PK Type ----------- ----- ------------- -pseudos - Si 56785 UpfData -code 56729 Code -kpoints 59515 KpointsData -parameters 59516 Dict -structure 56731 StructureData - -Outputs PK Type ------------------ ----- -------------- -output_band 59520 BandsData -output_parameters 59522 Dict -output_trajectory 59521 TrajectoryData -remote_folder 59518 RemoteData -retrieved 59519 FolderData -``` - -The calculation received a number of inputs: `pseudos`, `code`, `kpoints`, `parameters` and a `structure`. -The `code` was specified through the `-X` option when we invoked the CLI command to launch the calculation. -Other inputs, such as `structure`, `parameters` etc. were created automatically by the CLI command. - -The inputs can be introspected through the `verdi` CLI that ships with `aiida-core`. -For example, the structure can be visualized using the command: - -```console -verdi data core.structure show 56731 -``` - -This should open a new dialog with a visualization of the silicon crystal structure that was automatically created for this calculation. -Similarly, the `parameters` can be shown using: - -```console -verdi node attributes 59516 - -PK: 59516 -{ - "CONTROL": { - "calculation": "scf" - }, - "SYSTEM": { - "ecutrho": 240.0, - "ecutwfc": 30.0 - } -} -``` - -Likewise, the output parameters (`output_parameters`), which are of the same type (`Dict`) of the input parameters (`parameters`) can also be shown: - -```console -verdi node attributes 59522 - -PK: 59522 -{ - "beta_real_space": false, - "convergence_info": { - "scf_conv": { - "convergence_achieved": true, - "n_scf_steps": 6, - "scf_error": 2.1572731407236e-07 - } - }, - "creator_name": "pwscf", - "creator_version": "6.6", - "energy": -308.19211637089, - ... - "energy_accuracy": 5.850447441879e-06, - "energy_accuracy_units": "eV", - "fermi_energy": 6.508401516556, - "fermi_energy_units": "eV", - "number_of_atoms": 2, - "number_of_bands": 4, - "number_of_bravais_symmetries": 48, - "number_of_electrons": 8.0, - "number_of_k_points": 3, - "number_of_species": 1, - "number_of_spin_components": 1, - "number_of_symmetries": 48, - "scf_iterations": 6, - "volume": 40.02575122515, - "wall_time_seconds": 1.23, -} -``` - -The output will show a number of parameters that have been automatically parsed from the output of the `pw.x` calculation. -Note that the output has been shortened for clarity and not all parsed output parameters are shown. - -(tutorials-pw-through-api)= - -## Running a `pw.x` calculation through the API - -The following minimal example shows how to run a `pw.x` calculation through AiiDA's Python API. -For the purposes of this demonstration, the electronic ground state of an fcc silicon crystal structure is computed. - -:::{tip} -The code that is shown below in snippets can be {download}`downloaded as a script `. -Save the script to a file, make it executable and then simply execute it to run the example calculation. +```{toctree} +:hidden: true +:maxdepth: 1 + +first_pw +magnetism +``` + +:::{card} +:class-header: panel-header-text +:class-footer: tutor-footer +:link: tutorials-pw-through-cli +:link-type: ref +:margin: 4 + +{fa}`fa-regular fa-rocket` Running your first `pw.x` calculation +^^^ +Learn how to run the Quantum ESPRESSO `pw.x` with AiiDA, both from the command line and using the Python API. ++++ +::::{list-table} +:class: footer-table +:widths: 50 50 +* - {fa}`fa-sharp fa-regular fa-clock` 30 min + - {{ aiida_logo }} [Beginner]{.aiida-green} +:::: ::: -First, import the required classes and functions: - -```python -from aiida.engine import run -from aiida.orm import Dict, KpointsData, StructureData, load_code, load_group -``` - -Then, load the code that was setup in AiiDA for `pw.x` and get an instance of the [process builder](https://aiida.readthedocs.io/projects/aiida-core/en/latest/topics/processes/usage.html#process-builder): - -```python -# Load the code configured for ``pw.x``. Make sure to replace this string -# with the label that you used in the setup of the code. -code = load_code('pw@localhost') -builder = code.get_builder() -``` - -The process builder can be used to assign the inputs that are to be used for the calculation. -When an input is assigned, it is automatically validated. -Let's first define a structure and assign it to the builder: - -```python -# Create a silicon fcc crystal -from ase.build import bulk -structure = StructureData(ase=bulk('Si', 'fcc', 5.43)) -builder.structure = structure -``` - -To construct a crystal structure, we use the `bulk` method of the [ASE library](https://wiki.fysik.dtu.dk/ase/ase/build/build.html#ase.build.bulk). -This is then wrapped in a `StructureData` node such that it can be stored in AiiDA's provenance graph, and assigned to the builder. - -Next, we need to define the pseudopotentials. -If the pseudopotentials were installed using the `aiida-pseudo` package, as described in the [installation guide](installation-setup-pseudopotentials), that can be done as follows: - -```python -# Load the pseudopotential family. -pseudo_family = load_group('SSSP/1.1/PBE/efficiency') -builder.pseudos = pseudo_family.get_pseudos(structure=structure) -``` - -The parameters of the input file that will be read by `pw.x` are defined through the `parameters` input. -A minimal setup to run an SCF calculation is shown below: - -```python -# Request the recommended wavefunction and charge density cutoffs -# for the given structure and energy units. -cutoff_wfc, cutoff_rho = pseudo_family.get_recommended_cutoffs( - structure=structure, - unit='Ry' -) - -parameters = Dict({ - 'CONTROL': { - 'calculation': 'scf' - }, - 'SYSTEM': { - 'ecutwfc': cutoff_wfc, - 'ecutrho': cutoff_rho, - } -}) -builder.parameters = parameters -``` -:::{note} -The `pw.x` code actually requires more input parameters than have been specified in the example above. -The other required parameters are actually automatically defined by the `PwCalculation` plugin, based on for example the input structure. -To find out what other parameters can be specified, please refer to the [documentation of Quantum ESPRESSO](https://www.quantum-espresso.org/Doc/INPUT_PW.html). -Note that certain parameters cannot be manually specified as these are reserved for the plugin to be set. +:::{card} +:class-header: panel-header-text +:class-footer: tutor-footer +:link: tutorials-magnetic-configurations +:link-type: ref +:margin: 4 + +{fa}`fa-solid fa-arrow-down-up-across-line` Magnetic configurations +^^^ +Learn how to assign magnetic configurations to your structure, and retrieve the final magnetic configuration from a `pw.x` calculation. ++++ +::::{list-table} +:class: footer-table +:widths: 50 50 +* - {fa}`fa-sharp fa-regular fa-clock` 20 min + - {{ aiida_logo }} [Beginner]{.aiida-green} +:::: ::: - -The k-points to be used for the calculation are defined using the `KpointsData` class: - -```python -# Generate a 2x2x2 Monkhorst-Pack mesh -kpoints = KpointsData() -kpoints.set_kpoints_mesh([2, 2, 2]) -builder.kpoints = kpoints -``` - -Finally, specify the amount of CPUs that the calculation should use and how long it can run before it should be killed: - -```python -# Run the calculation on 1 CPU and kill it if it runs longer than 1800 seconds. -# Set ``withmpi`` to ``False`` if ``pw.x`` was compiled without MPI support. -builder.metadata.options = { - 'resources': { - 'num_machines': 1, - }, - 'max_wallclock_seconds': 1800, - 'withmpi': False, -} -``` - -All required inputs have now been defined. -The calculation can be started by launching the process builder: - -```python -results, node = run.get_node(builder) -``` - -When the calculation is finished, `run.get_node` will return a tuple of two values: `results` and `node`. -The `results` will be a dictionary of all the output nodes that were produced and the `node` will hold a reference to the node that represents the calculation in AiiDA's provenance graph. -When the calculation is finished, check whether it was successful through its exit status: - -```python -node.exit_status -``` - -If the calculation finished successfully, this will return `0`. -A non-zero exit code means the primary objective of the calculation was not (fully) achieved or an error occurred. - -The `results` variable contains a dictionary of the created output nodes. -Print it to see exactly which outputs were created and what their identifiers are: - -```python -print(results) -{ - 'output_band': , - 'output_trajectory': , - 'output_parameters': , - 'remote_folder': , - 'retrieved': -} -``` - -To show the parsed output parameters, call: - -```python -results['output_parameters'].base.attributes.all -{ - 'energy': -153.30519007432, - 'volume': 40.02575122515, - .... -} -``` - -The complete output that was written by `pw.x` to stdout, can be retrieved as follows: - -```python -results['retrieved'].base.repository.get_object_content('aiida.out') -``` diff --git a/docs/source/tutorials/magnetism.md b/docs/source/tutorials/magnetism.md new file mode 100644 index 000000000..8ff5cf597 --- /dev/null +++ b/docs/source/tutorials/magnetism.md @@ -0,0 +1,366 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.14.5 + kernelspec: + display_name: QE-dev + language: python + name: qe-dev +--- + +(tutorials-magnetic-configurations)= + +# Magnetic configurations + +There are several ways of defining a magnetic configuration in Quantum ESPRESSO. +This tutorial focuses on using {{ starting_magnetization }} to specify the initial magnetization for each kind in our system. + +(tutorials-magnetic-configuration-parameters)= + +## Using the `parameters` + +One way to specify the initial magnetic configuration of your system is by simply using the `parameters` input of the {{ PwCalculation }}. +Start by loading your default profile: + +```python +from aiida import orm, load_profile + +load_profile() +``` + +Next, set up a basic structure for HCP cobalt, and load the code you have set up for `pw.x`: + +:::{margin} +❗️Make sure to replace the label of the `pw.x` code in the {func}`~aiida.orm.load_code` function by that of _your_ installed code. +::: + +```python +from ase.build import bulk + +structure = orm.StructureData(ase=bulk('Co', 'hcp', 2.5, 4.06)) +pw_code = orm.load_code('qe-7.2-pw@localhost') +``` + +Obtain a pre-populated builder using the {{ get_builder_from_protocol }} method: + +:::{margin} +


    +🚀 Use the `fast` protocol for testing and demonstrations! +::: + +```python +from aiida_quantumespresso.workflows.pw.base import PwBaseWorkChain + +builder = PwBaseWorkChain.get_builder_from_protocol( + structure=structure, + code=pw_code, + protocol='fast' +) +``` + +In order to set the magnetic moments, you need to change the `parameters` input of the {{ PwCalculation }}: + +```python +parameters = builder.pw.parameters.get_dict() +parameters['SYSTEM']['nspin'] = 2 +parameters['SYSTEM']['starting_magnetization'] = {'Co': 1} +builder.pw.parameters = orm.Dict(parameters) +``` + +The commands above do the following: + +1. Get the `parameters` input of the `pw` namespace as a regular Python `dict`. +2. Set {{ nspin }} to `2` so Quantum ESPRESSO does a spin-polarised calculation. +3. Set the {{ starting_magnetization }} for the `Co` atomic species or kind to 1. +4. Replace the `parameters` of the {{ PwCalculation }} by a new `Dict` node that has the {{ starting_magnetization }} set. + +:::{tip} +Note that {{ starting_magnetization }} is a "multidimensional parameter", which in the `aiida-quantumespresso` plugin you shouldn't specify using indices between brackets as you would in a regular Quantum ESPRESSO input file. +See the [corresponding how-to section](howto-calculations-pw-multidimensional-parameters) for more details. +::: + +Now we ready to submit the builder to the AiiDA daemon: + +:::{margin} +⏱ Since this cell actually runs a `pw.x` calculation, it might take a bit longer to complete. +::: + +```python +from aiida.engine import run_get_node + +results, pwbase_node = run_get_node(builder) +``` + +Once the calculation is completed, you can for example get the total magnetization in the unit cell from the output parameters: + +```python +results['output_parameters'].get_dict()['total_magnetization'] +``` + +:::{important} +The {{ starting_magnetization }} input is defined as "spin polarization" of the atom species, with inputs ranging from -1 (all spins down) to 1 (all spins up). +However, the "magnetization" outputs are always defined in units of Bohr magneton. +We will commonly refer to the values in Bohr magneton as the "magnetic moments" to differentiate them from "magnetization". +::: + +(tutorials-magnetic-configuration-spintype)= + +## Using the {{ SpinType }} + +An alternative way of performing spin-polarised calculations is by specifying a {{ SpinType }} and passing it to the {{ get_builder_from_protocol }} method: + +```python +from aiida_quantumespresso.common.types import SpinType + +builder = PwBaseWorkChain.get_builder_from_protocol( + structure=structure, + code=pw_code, + protocol='fast', + spin_type=SpinType.COLLINEAR +) +``` + +In this case, an initial {{ starting_magnetization }} will be set for each element based on the following logic: + +- If the element has a partially occupied *d* or *f* shell, calculate what the {{ starting_magnetization }} should be in case we aim for a maximally possible magnetic moment, i.e. 5 or 7 Bohr magneton for *d* or *f* shells. +- For other elements, simply set {{ starting_magnetization }} to 0.1. + + +:::{note} +The approach outlined above will always set *positive* values for {{ starting_magnetization }}, and hence can only initialise the system in a ferro or ferrimagnetic state. +::: + +You can have a look to see what {{ starting_magnetization }} was set by inspecting the contents of the `parameters` set in the builder: + +```python +builder.pw.parameters.get_dict()['SYSTEM']['starting_magnetization'] +``` + +Once again, try running the calculation: + +```python +results, pwbase_node = run_get_node(builder) +``` +The calculation should have the same final total magnetization as the one in [section using the parameters](tutorials-magnetic-configuration-parameters). + +```python +results['output_parameters'].get_dict()['total_magnetization'] +``` + +## Anti-ferromagnetic configurations + + +What if you want to initialize the system in an anti-ferromagnetic state? +In this case, Quantum ESPRESSO requires two different kinds to be specified in the inputs, but looking at the `sites` in the input structure: + +```python +structure.sites +``` + +You can see that although our structure has two _sites_, they are both of the same _kind_. +This means you'll have to create a new `` node with two _different_ kinds. +For this purpose, you can use the {{ create_magnetic_configuration }} calculation function: + +```python +from aiida_quantumespresso.calculations.functions.create_magnetic_configuration import create_magnetic_configuration + +results = create_magnetic_configuration(structure, [-2, 2]) +``` + +Note that as {{ create_magnetic_configuration }} is a calculation function, and so its execution is tracked in the provenance: + +```console +❯ verdi process list -ap1 + PK Created Process label Process State Process status +---- --------- ---------------------------- --------------- ---------------- + [...] + 695 12s ago create_magnetic_configuration ⏹ Finished [0] + +Total results: 8 + +Report: last time an entry changed state: 3m ago (at 02:06:20 on 2023-05-28) +``` + +You can see the in- and outputs of the {{ create_magnetic_configuration }} process using `verdi process show`: + +```bash +❯ verdi process show 695 +Property Value +----------- ------------------------------------ +type create_magnetic_configuration +state Finished [0] +pk 695 +uuid bb8d4026-19c2-4efd-a99d-b55336df0eb3 +label create_magnetic_configuration +description +ctime 2023-05-28 04:50:25.887472+02:00 +mtime 2023-05-28 04:50:26.098163+02:00 + +Inputs PK Type +------------------------ ---- ------------- +atol 693 Float +magnetic_moment_per_site 692 List +structure 659 StructureData +ztol 694 Float + +Outputs PK Type +---------------- ---- ------------- +magnetic_moments 697 Dict +structure 696 StructureData +``` + +:::{important} +Every time you execute the {{ create_magnetic_configuration }} calculation function, it will store another process with corresponding in- and output nodes to the database! +If you want to avoid this, you can set `store_provenance` to `False` in the `metadata`: + +```python +results = create_magnetic_configuration(structure, [-2, 2], metadata={"store_provenance": False}) +``` +::: + +Let's have a look at the `structure` output: + +```python +new_structure = results['structure'] +new_structure +``` + +As you can see, this is the new {{ StructureData }}, which this time has two different kind names for each of the `Co` sites: + +```python user_expressions=[] +results['structure'].sites +``` + +The `results` also contain the `magnetic_moments` output: + +```python +results['magnetic_moments'].get_dict() +``` + +This can be passed to the `initial_magnetic_moments` input of the {{ get_builder_from_protocol }} method, along with the `new_structure`: + +```python +builder = PwBaseWorkChain.get_builder_from_protocol( + structure=new_structure, + code=pw_code, + protocol='fast', + spin_type=SpinType.COLLINEAR, + initial_magnetic_moments=results['magnetic_moments'] +) +``` + +Let's have another look at the {{ starting_magnetization }} input: + +```python +builder.pw.parameters['SYSTEM']['starting_magnetization'] +``` + +As always, we can run the calculation: + +```python +results, pwbase_node = run_get_node(builder) +``` + +In this case, the _total_ magnetization in the unit cell is zero: + +```python +results['output_parameters'].get_dict()['total_magnetization'] +``` + +But the _absolute_ magnetization is not: + +```python +results['output_parameters'].get_dict()['absolute_magnetization'] +``` + +## Passing the magnetic configuration + +Often you'll want to pass the magnetic configuration obtained from one calculation to the next. +The final magnetic configuration of the `pw.x` run can be found in the `output_trajectory` output: + +```python +trajectory = pwbase_node.outputs.output_trajectory +``` + +This {class}`aiida.orm.TrajectoryData` node contains a lot of information for each _ionic_ step performed in the `pw.x` calculation, all stored as `numpy` arrays in the repository. +To see all the arrays stored you can use the `get_arraynames()` method: + +```python +trajectory.get_arraynames() +``` + +We're interested in the `atomic_magnetic_moments`. +Retrieve the _final_ magnetic moments using the `get_array()` method: + +```python +magnetic_moments = trajectory.get_array('atomic_magnetic_moments')[-1] +``` + +:::{note} +Here you retrieved the _final_ array of magnetic moments by specifying the `-1` index. +Of course, for a static (`scf`) calculation there is only one ionic step in the trajectory. +If you had run a `vc-relax` calculation, the arrays in `output_trajectory` would contain all the information above for _each_ ionic step. +::: + +Have a look at the contents of `magnetic_moments`. +In this case, the magnetic moments are stored as a one-dimensional array of length 2. +Note that Quantum ESPRESSO prints the magnetic configuration as a list of magnetic moments: + +```fortran + Magnetic moment per site (integrated on atomic sphere of radius R) + atom 1 (R=0.411) charge= 15.6004 magn= 1.2090 + atom 2 (R=0.411) charge= 15.6004 magn= -1.2087 +``` + +And that, same as for the total magnetization, the magnetic moments are expressed in Bohr magneton. + +Due to this discrepancy between the in- and outputs in Quantum ESPRESSO, it becomes a bit challenging to properly pass magnetic configurations between calculations. +You'd have to perform roughly the following steps: + +1. Retrieve the final structure and atomic moments. + Which structure to use will depend on the calculation you have run. +2. Check if the structure needs to have different kinds defined in order to specify the magnetic moments. +3. Properly map the kinds to the the correct magnetic moments, grouping together sites whose magnetic moments are close. +3. Check if all magnetic moments are very small. + In this case continue with a non-magnetic calculation. + +Fortunately, the {{ PwCalculation }} comes with a handy tool for obtaining the magnetic configuration. +First, obtain the final {{ PwCalculation }} from the `called` processes of the {{ PwBaseWorkChain }} node: + +```python +pw_node = pwbase_node.called[-1] +``` + +Then, you can simply use the {{ get_magnetic_configuration }} method in the `tools` name space: + +```python +configuration = pw_node.tools.get_magnetic_configuration() +``` + +:::{note} +In case new kinds are needed for the configuration, the {{ create_magnetic_configuration }} calculation function is run inside {{ get_magnetic_configuration }} to generate the new structure. +This step is added to the provenance and hence stored in the database. +::: + +The configuration contains both the new `structure` and `magnetic_moments`, which can then be easily passed to the {{ get_builder_from_protocol }} method: + +```python +builder = PwBaseWorkChain.get_builder_from_protocol( + structure=configuration.structure, + code=pw_code, + spin_type=SpinType.COLLINEAR, + initial_magnetic_moments=configuration.magnetic_moments +) +``` + +Have a look at the {{ starting_magnetization }} input that was generated: + +```python +builder.pw.parameters['SYSTEM']['starting_magnetization'] +``` + +Is it different from the one at the start of the initial run? diff --git a/pyproject.toml b/pyproject.toml index 6d11aa0a8..1f7454141 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -68,6 +68,7 @@ aiida-quantumespresso = 'aiida_quantumespresso.cli:cmd_root' [project.entry-points.'aiida.calculations'] 'quantumespresso.cp' = 'aiida_quantumespresso.calculations.cp:CpCalculation' 'quantumespresso.create_kpoints_from_distance' = 'aiida_quantumespresso.calculations.functions.create_kpoints_from_distance:create_kpoints_from_distance' +'quantumespresso.create_magnetic_configuration' = 'aiida_quantumespresso.calculations.functions.create_magnetic_configuration:create_magnetic_configuration' 'quantumespresso.merge_ph_outputs' = 'aiida_quantumespresso.calculations.functions.merge_ph_outputs:merge_ph_outputs' 'quantumespresso.dos' = 'aiida_quantumespresso.calculations.dos:DosCalculation' 'quantumespresso.epw' = 'aiida_quantumespresso.calculations.epw:EpwCalculation' diff --git a/src/aiida_quantumespresso/calculations/functions/create_magnetic_configuration.py b/src/aiida_quantumespresso/calculations/functions/create_magnetic_configuration.py new file mode 100644 index 000000000..b670123fd --- /dev/null +++ b/src/aiida_quantumespresso/calculations/functions/create_magnetic_configuration.py @@ -0,0 +1,123 @@ +# -*- coding: utf-8 -*- +"""Create a new magnetic configuration from the given structure based on the desired magnetic moments.""" +from aiida.engine import calcfunction +from aiida.orm import Float +import numpy + + +@calcfunction +def create_magnetic_configuration( + structure, magnetic_moment_per_site, atol=lambda: Float(0.5), ztol=lambda: Float(0.05) +): + """Create a new magnetic configuration from the given structure based on a list of magnetic moments per site. + + To create the new list of kinds, the algorithm loops over all the elements in the structure and makes a list of the + sites with that element and their corresponding magnetic moment. Next, it splits this list in three lists: + + * Zero magnetic moments: Any site that has an absolute magnetic moment lower than ``ztol`` + * Positive magnetic moments + * Negative magnetic moments + + The algorithm then sorts the positive and negative lists from large to small absolute value, and loops over each of + list. New magnetic kinds will be created when the absolute difference between the magnetic moment of the current + kind and the site exceeds ``atol``. + + The positive and negative magnetic moments are handled separately to avoid assigning two sites with opposite signs + in their magnetic moment to the same kind and make sure that each kind has the correct magnetic moment, i.e. the + largest magnetic moment in absolute value of the sites corresponding to that kind. + + .. important:: the function currently does not support alloys. + + :param structure: a `StructureData` instance. + :param magnetic_moment_per_site: list of magnetic moments for each site in the structure. + :param atol: the absolute tolerance on determining if two sites have the same magnetic moment. + :param ztol: threshold for considering a kind to have non-zero magnetic moment. + """ + # pylint: disable=too-many-locals,too-many-branches,too-many-statements + import string + + from aiida.orm import Dict, StructureData + + if structure.is_alloy: + raise ValueError('Alloys are currently not supported.') + + atol = atol.value + rtol = 0 # Relative tolerance used in the ``numpy.is_close()`` calls. + ztol = ztol.value + + new_structure = StructureData(cell=structure.cell, pbc=structure.pbc) + magnetic_configuration = {} + + for element in structure.get_symbols_set(): + + # Filter the sites and magnetic moments on the site element + element_sites, element_magnetic_moments = zip( + *[(site, magnetic_moment) + for site, magnetic_moment in zip(structure.sites, magnetic_moment_per_site) + if site.kind_name.rstrip(string.digits) == element] + ) + + # Split the sites and their magnetic moments by sign to filter out the sites with magnetic moment lower than + # `ztol`and deal with the positive and negative magnetic moment separately. This is important to avoid assigning + # two sites with opposite signs to the same kind and make sure that each kind has the correct magnetic moment, + # i.e. the largest magnetic moment in absolute value of the sites corresponding to that kind. + zero_sites = [] + pos_sites = [] + neg_sites = [] + + for site, magnetic_moment in zip(element_sites, element_magnetic_moments): + + if abs(magnetic_moment) <= ztol: + zero_sites.append((site, 0)) + elif magnetic_moment > 0: + pos_sites.append((site, magnetic_moment)) + else: + neg_sites.append((site, magnetic_moment)) + + kind_index = -1 + kind_names = [] + kind_sites = [] + kind_magnetic_moments = {} + + for site_list in (zero_sites, pos_sites, neg_sites): + + if not site_list: + continue + + # Sort the site list in order to build the kind lists from large to small absolute magnetic moment. + site_list = sorted(site_list, key=lambda x: abs(x[1]), reverse=True) + + sites, magnetic_moments = zip(*site_list) + + kind_index += 1 + current_kind_name = f'{element}{kind_index}' + kind_sites.append(sites[0]) + kind_names.append(current_kind_name) + kind_magnetic_moments[current_kind_name] = magnetic_moments[0] + + for site, magnetic_moment in zip(sites[1:], magnetic_moments[1:]): + + if not numpy.isclose(magnetic_moment, kind_magnetic_moments[current_kind_name], rtol, atol): + kind_index += 1 + current_kind_name = f'{element}{kind_index}' + kind_magnetic_moments[current_kind_name] = magnetic_moment + + kind_sites.append(site) + kind_names.append(current_kind_name) + + # In case there is only a single kind for the element, remove the 0 kind index + if current_kind_name == f'{element}0': + kind_names = len(element_magnetic_moments) * [element] + kind_magnetic_moments = {element: kind_magnetic_moments[current_kind_name]} + + magnetic_configuration.update(kind_magnetic_moments) + + for name, site in zip(kind_names, kind_sites): + new_structure.append_atom( + name=name, + symbols=(element,), + weights=(1.0,), + position=site.position, + ) + + return {'structure': new_structure, 'magnetic_moments': Dict(dict=magnetic_configuration)} diff --git a/src/aiida_quantumespresso/tools/calculations/pw.py b/src/aiida_quantumespresso/tools/calculations/pw.py index bf8ce43f3..498065102 100644 --- a/src/aiida_quantumespresso/tools/calculations/pw.py +++ b/src/aiida_quantumespresso/tools/calculations/pw.py @@ -1,8 +1,10 @@ # -*- coding: utf-8 -*- """Tools for nodes created by running the `PwCalculation` class.""" -from aiida.common import exceptions +from aiida.common import AttributeDict, exceptions from aiida.tools.calculations.base import CalculationTools +from aiida_quantumespresso.calculations.functions.create_magnetic_configuration import create_magnetic_configuration + class PwCalculationTools(CalculationTools): """Calculation tools for `PwCalculation`. @@ -51,3 +53,45 @@ def get_scf_accuracy(self, index=0): return scf_accuracy[scf_accuracy_index[index - 1]:scf_accuracy_index[index]] return scf_accuracy[scf_accuracy_index[index]:scf_accuracy_index[index + 1]] + + def get_magnetic_configuration(self, atol: float = 0.5, ztol: float = 0.05) -> AttributeDict: + """Get the final magnetic configuration of a ``pw.x`` calculation.""" + try: + structure = self._node.outputs.output_structure + except AttributeError: + structure = self._node.inputs.structure + + try: + magnetic_moments = self._node.outputs.output_trajectory.get_array('atomic_magnetic_moments')[-1].tolist() + except KeyError as exc: + raise KeyError('Could not find the `atomic_magnetic_moments` in the `output_trajectory`.') from exc + + results = create_magnetic_configuration( + structure, + magnetic_moments, + atol=atol, + ztol=ztol, + metadata={'store_provenance': False}, + ) + structure_kindname_position = [(site.kind_name, site.position) for site in structure.sites] + allo_kindname_position = [(site.kind_name, site.position) for site in results['structure'].sites] + structure_kindnames_sorted = sorted(structure_kindname_position, key=lambda l: l[1]) + allo_kindnames_sorted = sorted(allo_kindname_position, key=lambda l: l[1]) + + requires_new_kinds = not structure_kindnames_sorted == allo_kindnames_sorted + non_magnetic = all(abs(magn) < ztol for magn in results['magnetic_moments'].get_dict().values()) + + if requires_new_kinds: + results = create_magnetic_configuration( + structure, + magnetic_moments, + atol=atol, + ztol=ztol, + metadata={'store_provenance': False}, + ) + structure = results['results'] + + return AttributeDict({ + 'structure': structure, + 'magnetic_moments': None if non_magnetic else results['magnetic_moments'] + }) diff --git a/src/aiida_quantumespresso/workflows/protocols/utils.py b/src/aiida_quantumespresso/workflows/protocols/utils.py index 3eace037c..664fe17c9 100644 --- a/src/aiida_quantumespresso/workflows/protocols/utils.py +++ b/src/aiida_quantumespresso/workflows/protocols/utils.py @@ -3,12 +3,10 @@ import pathlib from typing import Optional, Union -from aiida.plugins import DataFactory, GroupFactory +from aiida.orm import StructureData +from aiida_pseudo.groups.family import PseudoPotentialFamily import yaml -StructureData = DataFactory('core.structure') -PseudoPotentialFamily = GroupFactory('pseudo.family') - class ProtocolMixin: """Utility class for processes to build input mappings for a given protocol based on a YAML configuration file.""" diff --git a/tests/calculations/functions/test_create_magnetic_configuration.py b/tests/calculations/functions/test_create_magnetic_configuration.py new file mode 100644 index 000000000..38f4bfe82 --- /dev/null +++ b/tests/calculations/functions/test_create_magnetic_configuration.py @@ -0,0 +1,207 @@ +# -*- coding: utf-8 -*- +"""Tests for the `create_magnetic_configuration` calculation function.""" +from aiida.orm import Float, List +import pytest + +from aiida_quantumespresso.calculations.functions.create_magnetic_configuration import create_magnetic_configuration + + +@pytest.mark.usefixtures('aiida_profile') +def test_configuration_00(generate_structure_from_kinds): + """Test `create_magnetic_configuration` calculation function. + + Case: one kind but with equal magnetic moments. + Expected result: no new kind names should be introduced. + """ + kind_names = ['Fe', 'Fe'] + magnetic_moments = List(list=[0.2, 0.2]) + + structure = generate_structure_from_kinds(kind_names) + allotrope, allotrope_magnetic_moments = create_magnetic_configuration(structure, magnetic_moments).values() + + assert set(allotrope.get_kind_names()) == {'Fe'} + assert allotrope_magnetic_moments.get_dict() == {'Fe': 0.2} + + +@pytest.mark.usefixtures('aiida_profile') +def test_configuration_01(generate_structure_from_kinds): + """Test `create_magnetic_configuration` calculation function. + + Case: two kinds all with equal magnetic moments. + Expected result: no new kind names should be introduced. + """ + kind_names = ['Fe', 'Fe', 'Ni', 'Ni', 'Ni'] + magnetic_moments = List(list=[0.2, 0.2, 0.5, 0.5, 0.5]) + + structure = generate_structure_from_kinds(kind_names) + allotrope, allotrope_magnetic_moments = create_magnetic_configuration(structure, magnetic_moments).values() + + assert set(allotrope.get_kind_names()) == {'Fe', 'Ni'} + assert allotrope_magnetic_moments.get_dict() == {'Fe': 0.2, 'Ni': 0.5} + + +@pytest.mark.usefixtures('aiida_profile') +def test_configuration_02(generate_structure_from_kinds): + """Test `create_magnetic_configuration` calculation function. + + Case: only one kind but with unequal magnetic moments. + Expected result: two new kinds introduced one for each magnetic moment. + """ + kind_names = ['Fe', 'Fe'] + magnetic_moments = List(list=[0.2, 1.0]) + + structure = generate_structure_from_kinds(kind_names) + allotrope, allotrope_magnetic_moments = create_magnetic_configuration(structure, magnetic_moments).values() + + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1'} + assert allotrope_magnetic_moments.get_dict() == {'Fe0': 1.0, 'Fe1': 0.2} + + +@pytest.mark.usefixtures('aiida_profile') +def test_configuration_03(generate_structure_from_kinds): + """Test `create_magnetic_configuration` calculation function. + + Case: only one kind but with three types of magnetic moments that are not grouped together. + Expected result: two new kinds introduced one for each magnetic moment. + """ + kind_names = ['Fe', 'Fe', 'Fe', 'Fe'] + magnetic_moments = List(list=[0.2, 0.8, 1.5, 0.8]) + + structure = generate_structure_from_kinds(kind_names) + allotrope, allotrope_magnetic_moments = create_magnetic_configuration(structure, magnetic_moments).values() + + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1', 'Fe2'} + assert allotrope_magnetic_moments.get_dict() == {'Fe0': 1.5, 'Fe1': 0.8, 'Fe2': 0.2} + + +@pytest.mark.usefixtures('aiida_profile') +def test_configuration_04(generate_structure_from_kinds): + """Test `create_magnetic_configuration` calculation function. + + Case: only one kind but with four different values of magnetic moments but middle two are within tolerance. + Expected result: two new kinds introduced one for each magnetic moment. + """ + kind_names = ['Fe', 'Fe', 'Fe', 'Fe'] + magnetic_moments = List(list=[0.0, 0.50, 0.45, 0.40]) + + structure = generate_structure_from_kinds(kind_names) + + # Default tolerances: just two different kinds + allotrope, allotrope_magnetic_moments = create_magnetic_configuration(structure, magnetic_moments).values() + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1'} + assert [site.kind_name for site in allotrope.sites] == ['Fe0', 'Fe1', 'Fe1', 'Fe1'] + assert allotrope_magnetic_moments.get_dict() == {'Fe0': 0.0, 'Fe1': 0.5} + + # Lower atol to 0.05: 0.5 & 0.45 now one kind, 0.4 new kind -> three different kinds + allotrope, allotrope_magnetic_moments = create_magnetic_configuration( + structure, magnetic_moments, atol=Float(0.05) + ).values() + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1', 'Fe2'} + assert [site.kind_name for site in allotrope.sites] == ['Fe0', 'Fe1', 'Fe1', 'Fe2'] + assert allotrope_magnetic_moments.get_dict() == {'Fe0': 0.0, 'Fe1': 0.5, 'Fe2': 0.4} + + # Increase atol to 0.1, again only two different kinds + allotrope, allotrope_magnetic_moments = create_magnetic_configuration(structure, magnetic_moments, + atol=Float(0.1)).values() + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1'} + assert [site.kind_name for site in allotrope.sites] == ['Fe0', 'Fe1', 'Fe1', 'Fe1'] + assert allotrope_magnetic_moments.get_dict() == {'Fe0': 0.0, 'Fe1': 0.5} + + # Really strict tolerance or atol = 0.01: All sites get different kinds + allotrope, allotrope_magnetic_moments = create_magnetic_configuration( + structure, magnetic_moments, atol=Float(1E-2) + ).values() + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1', 'Fe2', 'Fe3'} + assert [site.kind_name for site in allotrope.sites] == ['Fe0', 'Fe1', 'Fe2', 'Fe3'] + assert allotrope_magnetic_moments.get_dict() == {'Fe0': 0.0, 'Fe1': 0.5, 'Fe2': 0.45, 'Fe3': 0.4} + + +@pytest.mark.usefixtures('aiida_profile') +def test_configuration_05(generate_structure_from_kinds): + """Test `create_magnetic_configuration` calculation function. + + Case: One kind, only negative magnetic moments with one close to zero + Expected result: Depends on tolerance, see below + """ + kind_names = ['Fe', 'Fe', 'Fe', 'Fe'] + magnetic_moments = List(list=[-0.5, -0.6, -1.5, -0.01]) + + structure = generate_structure_from_kinds(kind_names) + + # Default tolerance values, one zero site and two magnetic + allotrope, allotrope_magnetic_moments = create_magnetic_configuration(structure, magnetic_moments).values() + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1', 'Fe2'} + assert [site.kind_name for site in allotrope.sites] == ['Fe0', 'Fe1', 'Fe2', 'Fe2'] + assert allotrope_magnetic_moments.get_dict() == {'Fe0': 0.0, 'Fe1': -1.5, 'Fe2': -0.6} + + # Strict absolute tolerance, one zero site and three magnetic + allotrope, allotrope_magnetic_moments = create_magnetic_configuration( + structure, magnetic_moments, atol=Float(0.05) + ).values() + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1', 'Fe2', 'Fe3'} + assert [site.kind_name for site in allotrope.sites] == ['Fe0', 'Fe1', 'Fe2', 'Fe3'] + assert allotrope_magnetic_moments.get_dict() == {'Fe0': 0.0, 'Fe1': -1.5, 'Fe2': -0.6, 'Fe3': -0.5} + + # Strict absolute and zero tolerance, four magnetic sites + allotrope, allotrope_magnetic_moments = create_magnetic_configuration( + structure, magnetic_moments, atol=Float(0.05), ztol=Float(1E-3) + ).values() + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1', 'Fe2', 'Fe3'} + assert [site.kind_name for site in allotrope.sites] == ['Fe0', 'Fe1', 'Fe2', 'Fe3'] + assert allotrope_magnetic_moments.get_dict() == {'Fe0': -1.5, 'Fe1': -0.6, 'Fe2': -0.5, 'Fe3': -0.01} + + +@pytest.mark.usefixtures('aiida_profile') +def test_configuration_06(generate_structure_from_kinds): + """Test `create_magnetic_configuration` calculation function. + + Case: Two kinds, magnetic moments with different signs for the first (Fe) + Expected result: Depends on tolerance, see below + """ + kind_names = ['Fe', 'Fe', 'Fe', 'Fe', 'Ni', 'Ni'] + magnetic_moments = List(list=[-0.1, 0.1, -0.2, 0.01, 0.2, 0.25]) + + structure = generate_structure_from_kinds(kind_names) + + # Default tolerance values, one zero and two magnetic sites for Fe, one magnetic site for Ni + allotrope, allotrope_magnetic_moments = create_magnetic_configuration(structure, magnetic_moments).values() + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1', 'Fe2', 'Ni'} + assert allotrope_magnetic_moments.get_dict() == {'Fe0': 0.0, 'Fe1': 0.1, 'Fe2': -0.2, 'Ni': 0.25} + + # Very strict absolute tolerance, all different sites + allotrope, allotrope_magnetic_moments = create_magnetic_configuration( + structure, magnetic_moments, atol=Float(0.02) + ).values() + assert set(allotrope.get_kind_names()) == {'Fe0', 'Fe1', 'Fe2', 'Fe3', 'Ni0', 'Ni1'} + assert allotrope_magnetic_moments.get_dict() == { + 'Fe0': 0.0, + 'Fe1': 0.1, + 'Fe2': -0.2, + 'Fe3': -0.1, + 'Ni0': 0.25, + 'Ni1': 0.2 + } + + +@pytest.mark.usefixtures('aiida_profile') +def test_configuration_07(generate_structure_from_kinds): + """Test `create_magnetic_configuration` calculation function. + + Case: Two different symbols but the same magnetic moment. + Expected result: Depends on tolerance, see below + """ + kind_names = ['Fe0', 'Fe1'] + magnetic_moments = List(list=[0.1, 0.1]) + + structure = generate_structure_from_kinds(kind_names) + + # Default tolerance values, one kind with name equal to the element symbol + allotrope, allotrope_magnetic_moments = create_magnetic_configuration(structure, magnetic_moments).values() + assert set(allotrope.get_kind_names()) == {'Fe'} + assert allotrope_magnetic_moments.get_dict() == {'Fe': 0.1} + + # Very loose zero tolerance, one kind with zero magnetic moment + allotrope, allotrope_magnetic_moments = create_magnetic_configuration(structure, magnetic_moments, + ztol=Float(0.2)).values() + assert set(allotrope.get_kind_names()) == {'Fe'} + assert allotrope_magnetic_moments.get_dict() == {'Fe': 0} diff --git a/tests/conftest.py b/tests/conftest.py index e0b08b57e..a42e5c0ca 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -410,6 +410,29 @@ def _generate_structure(structure_id='silicon'): return _generate_structure +@pytest.fixture +def generate_structure_from_kinds(): + """Return a dummy `StructureData` instance with the specified kind names.""" + + def _generate_structure_from_kinds(site_kind_names): + """Return a dummy `StructureData` instance with the specified kind names.""" + import re + + from aiida import orm + + if not isinstance(site_kind_names, (list, tuple)): + site_kind_names = (site_kind_names,) + + structure = orm.StructureData(cell=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + + for kind_name in site_kind_names: + structure.append_atom(name=kind_name, symbols=re.sub('[0-9]', '', kind_name), position=(0., 0., 0.)) + + return structure + + return _generate_structure_from_kinds + + @pytest.fixture def generate_kpoints_mesh(): """Return a `KpointsData` node."""