Skip to content

Commit

Permalink
Bigmap/minor fixes (#306)
Browse files Browse the repository at this point in the history
* Fix minor issues that I noticed when going through the QE section
* Remove the Fermi level align with zero example, which uses bands.show_mpl(), but this doesn't work from the AiiDAlab machines since we don't have X forwarding.
* Replace the query demo script by the one from the intro week, which shows off the high-throughput elements more.
  • Loading branch information
mbercx committed Dec 2, 2020
1 parent 21ef1a8 commit 6c03bf1
Show file tree
Hide file tree
Showing 3 changed files with 117 additions and 56 deletions.
Binary file modified docs/pages/2020_BIGMAP/sections/include/images/demo_query.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
112 changes: 81 additions & 31 deletions docs/pages/2020_BIGMAP/sections/include/snippets/demo_query.py
Original file line number Diff line number Diff line change
@@ -1,44 +1,94 @@
from aiida.orm import Group
from aiida.plugins import DataFactory, CalculationFactory
from IPython.display import Image
from datetime import datetime, timedelta
import numpy as np
from matplotlib import gridspec, pyplot as plt

PwCalculation = CalculationFactory('quantumespresso.pw')
StructureData = DataFactory('structure')
KpointsData = DataFactory('array.kpoints')
Dict = DataFactory('dict')
UpfData = DataFactory('upf')

import matplotlib.pyplot as plt
def plot_results(query_res):
"""
:param query_res: The result of an instance of the QueryBuilder
"""
smearing_unit_set,magnetization_unit_set,pseudo_family_set = set(), set(), set()
# Storing results:
results_dict = {}
for pseudo_family, formula, smearing, smearing_units, mag, mag_units in query_res:
if formula not in results_dict:
results_dict[formula] = {}
# Storing the results:
results_dict[formula][pseudo_family] = (smearing, mag)
# Adding to the unit set:
smearing_unit_set.add(smearing_units)
magnetization_unit_set.add(mag_units)
pseudo_family_set.add(pseudo_family)

qb = QueryBuilder().append(
Group, filters={'label': {'like': 'tutorial_%'}}, tag='group',
project='label'
).append(
PwCalculation, with_group='group', tag='pw'
).append(
StructureData, with_outgoing='pw',
project='extras.formula'
).append(
Dict, with_incoming='pw', filters={'attributes.absolute_magnetization': {'>': 0.0}},
project='attributes.absolute_magnetization'
)
# Sorting by formula:
sorted_results = sorted(results_dict.items())
formula_list = next(zip(*sorted_results))
nr_of_results = len(formula_list)

results_dict = dict()
formulas = set()
# Checks that I have not more than 3 pseudo families.
# If more are needed, define more colors
#pseudo_list = list(pseudo_family_set)
if len(pseudo_family_set) > 3:
raise Exception('I was expecting 3 or less pseudo families')

for group_label, formula, abs_magnetization in qb.iterall():
functional = group_label.split('_')[1].upper()
results_dict.setdefault(functional, {})
results_dict[functional][formula] = abs_magnetization
formulas.add(formula)
colors = ['b', 'r', 'g']

formulas = list(formulas)
# Plotting:
plt.clf()
fig=plt.figure(figsize=(16, 9), facecolor='w', edgecolor=None)
gs = gridspec.GridSpec(2,1, hspace=0.01, left=0.1, right=0.94)

for functional, results in results_dict.items():
# Defining barwidth
barwidth = 1. / (len(pseudo_family_set)+1)
offset = [-0.5+(0.5+n)*barwidth for n in range(len(pseudo_family_set))]
# Axing labels with units:
yaxis = ("Smearing energy [{}]".format(smearing_unit_set.pop()),
"Total magnetization [{}]".format(magnetization_unit_set.pop()))
# If more than one unit was specified, I will exit:
if smearing_unit_set:
raise ValueError('Found different units for smearing')
if magnetization_unit_set:
raise ValueError('Found different units for magnetization')

abs_magnetizations = [results[formula] for formula in formulas]
# Making two plots, the top one for the smearing, the bottom one for the magnetization
for index in range(2):
ax=fig.add_subplot(gs[index])
for i,pseudo_family in enumerate(pseudo_family_set):
X = np.arange(nr_of_results)+offset[i]
Y = np.array([thisres[1][pseudo_family][index] for thisres in sorted_results])
ax.bar(X, Y, width=0.2, facecolor=colors[i], edgecolor=colors[i], label=pseudo_family)
ax.set_ylabel(yaxis[index], fontsize=14, labelpad=15*index+5)
ax.set_xlim(-0.5, nr_of_results-0.5)
ax.set_xticks(np.arange(nr_of_results))
if index == 0:
plt.setp(ax.get_yticklabels()[0], visible=False)
ax.xaxis.tick_top()
ax.legend(loc=3, prop={'size': 18})
else:
plt.setp(ax.get_yticklabels()[-1], visible=False)
for i in range(0, nr_of_results, 2):
ax.axvspan(i-0.5, i+0.5, facecolor='y', alpha=0.2)
ax.set_xticklabels(list(formula_list),rotation=90, size=14, ha='center')
plt.savefig('demo_query.pdf')

plt.plot(abs_magnetizations, 's')
plt.xticks(range(len(formulas)), formulas, rotation=90)
plt.ylabel('Magnetization [Bohrmag / cell]')

plt.legend(results_dict.keys())
plt.tight_layout()
plt.savefig('demo_query.pdf')
qb = QueryBuilder().append(
Group, filters={'label':{'like':'tutorial_%'}}, project='label', tag='group'
).append(
PwCalculation, tag='calculation', with_group='group'
).append(
StructureData, project=['extras.formula'], tag='structure', with_outgoing='calculation'
).append(
Dict, tag='results',
project=['attributes.energy_smearing', 'attributes.energy_smearing_units',
'attributes.total_magnetization', 'attributes.total_magnetization_units',
], with_incoming='calculation'
)

plot_results(qb.all())
61 changes: 36 additions & 25 deletions docs/pages/2020_BIGMAP/sections/qe.rst
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ It's a good idea to mark it down, but should you forget, you can always have a l
Total results: 4
The first column (marked `Id`) are the PK's of the ``StructureData`` nodes.
The first column (marked ``Id``) are the PK's of the ``StructureData`` nodes.

.. important::

Expand Down Expand Up @@ -166,12 +166,12 @@ You can get more information on an input by adding a question mark ``?``:
Here you can see that the ``structure`` input is required, needs to be of the ``StructureData`` type and is stored in the database (``"non_db": "False"``).

Next, we'll set up a dictionary with the pseudopotentials.
This can be done easily with a little utility function
This can be done easily with the ``get_pseudos_from_structure`` utility function.

.. code-block:: ipython
In [5]: from aiida.orm.nodes.data.upf import get_pseudos_from_structure
...: pseudos = get_pseudos_from_structure(structure, '<PSEUDO_FAMILY>')
...: pseudos = get_pseudos_from_structure(structure, 'SSSP_1.1_efficiency')
If we check the content of the ``pseudos`` variable:

Expand Down Expand Up @@ -283,7 +283,7 @@ To see *all* processes, use the ``--all`` option:
Info: last time an entry changed state: 28s ago (at 16:20:43 on 2020-11-29)
Notice how the band structure workflow (``PwBandsWorkChain``) you ran in the `Quantum ESPRESSO`_ app of `AiiDAlab`_ is also in the process list!
Use the PK of the most recent `PwCalculation` (the one you just sent) to get more information on it:
Use the PK of the most recent ``PwCalculation`` (the one you just sent) to get more information on it:

.. code-block:: console
Expand Down Expand Up @@ -375,8 +375,8 @@ To see all currently available workflows in your installation, you can run the f
$ verdi plugin list aiida.workflows
We are going to choose the ``PwBandStructureWorkChain `` workflow of the ``aiida-quantumespresso`` plugin (you can see it on the list as ``quantumespresso.pw.band_structure``).
We are going to run the ``PwBandStructureWorkChain`` workflow of the ``aiida-quantumespresso`` plugin.
You can see it on the list as ``quantumespresso.pw.band_structure``, which is the *entry point* of this workflow.
This is a fully automated workflow that will:

#. Determine the primitive cell of a given input structure.
Expand All @@ -387,8 +387,8 @@ This is a fully automated workflow that will:

The workflow uses the PBE exchange-correlation functional with suitable pseudopotentials and energy cutoffs from the `SSSP library version 1.1 <https://www.materialscloud.org/discover/sssp/table/efficiency>`_.

In order to run it, we will open again the ``verdi shell``.
We will then load the workflow plugin using the previously identified label and get a builder for the workflow:
In order to run it, we will again open the ``verdi shell``.
We will then load the workflow plugin using the previously identified entry point and get a builder for the workflow:

.. code-block:: ipython
Expand All @@ -405,7 +405,7 @@ Replace the following ``<CODE_LABEL>`` and ``<PK>`` with the corresponding value
...: builder.structure = load_node(<PK>) # REPLACE <PK>
Finally, we just need to submit the builder in the same as we did before for the calculation:
Finally, we just need to submit the builder in the same way as we did before for the calculation:

.. code-block:: ipython
Expand Down Expand Up @@ -553,13 +553,12 @@ As you can see, the explore tool of the `Materials Cloud <https://www.materialsc
However, you might already imagine that doing a more intensive kind of data mining of specific results this way can quickly become tedious.
For this use cases, AiiDA has a more versatile tool: the ``QueryBuilder``.


Finishing the workchain
-----------------------

Let's stop ``ngrok`` using ``Ctrl+C`` and close its terminal, as well as stop the REST API (also using ``Ctrl+C``).
The workchain we started earlier should be finished by now, let's use ``verdi process show <PK>`` to inspect the ``PwBandsWorkChain`` and find the PK of its ``band_structure`` output.
Use this to produce a PDF of the band structure:
Let's use ``verdi process show <PK>`` to inspect the ``PwBandsWorkChain`` and find the PK of its ``band_structure`` output.
Instead of relying on the explore tool, we can also plot the band structure using the ``verdi shell``:

.. code-block:: console
Expand All @@ -570,17 +569,29 @@ Use this to produce a PDF of the band structure:

Band structure computed by the ``PwBandStructureWorkChain``.

.. note::
The ``BandsData`` node does contain information about the Fermi energy, so the energy zero in your plot will be arbitrary.
You can produce a plot with the Fermi energy set to zero (as above) using the following steps in the ``verdi shell``.
Just look for the ``scf_parameters`` and ``band_structure`` output nodes of the ``PwBandStructureWorkChain`` using ``verdi process show`` and replace them in the following code:
Finally, the ``verdi process status`` command prints a *hierarchical* overview of the processes called by the work chain:

.. code-block:: ipython
.. code-block:: console
In [1]: scf_params = load_node(<PK>) # PK of the `scf_parameters` node
...: fermi_energy = scf_params.dict.fermi_energy
...: bands = load_node(<PK>) # PK of the `band_structure` node
...: bands.show_mpl(y_origin=fermi_energy, plot_zero_axis=True)
$ verdi process status 186
PwBandStructureWorkChain<186> Finished [0] [3:results]
└── PwBandsWorkChain<201> Finished [0] [7:results]
├── PwRelaxWorkChain<203> Finished [0] [3:results]
│ ├── PwBaseWorkChain<206> Finished [0] [7:results]
│ │ ├── create_kpoints_from_distance<208> Finished [0]
│ │ └── PwCalculation<212> Finished [0]
│ └── PwBaseWorkChain<223> Finished [0] [7:results]
│ ├── create_kpoints_from_distance<225> Finished [0]
│ └── PwCalculation<229> Finished [0]
├── seekpath_structure_analysis<236> Finished [0]
├── PwBaseWorkChain<243> Finished [0] [7:results]
│ ├── create_kpoints_from_distance<245> Finished [0]
│ └── PwCalculation<249> Finished [0]
└── PwBaseWorkChain<257> Finished [0] [7:results]
└── PwCalculation<260> Finished [0]
The bracket ``[3:result]`` indicates the current step in the outline of the ``PwBandStructureWorkChain`` (step 3, with name ``result``).
The ``process status`` is particularly useful for debugging complex work chains, since it helps pinpoint where a problem occurred.

Querying the database
---------------------
Expand Down Expand Up @@ -708,7 +719,7 @@ Let's have a look at the results:
['FeO3Sr', 3.38]]
You can see that we've found 7 magnetic structures for the calculations in the ``tutorial_pbesol`` group, along with their formulas and magnetizations.
We've set up a little script (:download:`demo_query.py <include/snippets/demo_query.py>`) that performs a similar query to obtain these results for all three groups, and then postprocess the data to make a simple plot.
We've set up a script (:download:`demo_query.py <include/snippets/demo_query.py>`) that performs a similar query to obtain the magnetization and smearing energy for all results in the three groups, and then postprocess the data to visualize it.
You can find it in the dropdown panel below:

.. dropdown:: **Query demo script**
Expand All @@ -727,14 +738,14 @@ and use ``verdi run`` to execute it:
$ verdi run demo_query.py
The resulting plot should look something like the one shown in :numref:`BIGMAP_2020_Query_demo`.
The resulting plot should look like the one shown in :numref:`BIGMAP_2020_Query_demo`.

.. _BIGMAP_2020_Query_demo:
.. figure:: include/images/demo_query.png
:width: 80%
:width: 100%
:align: center

Comparison of the absolute magnetization of the cell of the perovskite structures, calculated with different functionals.
Comparison of the absolute magnetization and smearing energy of the cell of the perovskite structures, calculated with different functionals.

What next?
----------
Expand Down

0 comments on commit 6c03bf1

Please sign in to comment.