Convergence parser for properties

src/aiida_sssp_workflow/calculations/calculate_bands_distance.py

@@ -53,7 +53,7 @@ def retrieve_bands(
     collect the bands of certain number with setting the start band.
     In order to make sure that when comparing two bands distance the number of bands is the same.
 
-    aligh bands to fermi level
+    align bands to fermi level
 
     The bands calculation of magnetic elements will giving a three dimensional bands where the
     first dimension is for the up and down spin.
@@ -236,7 +236,8 @@ def get_bands_distance(
         do_smearing,
     )
 
-    # after cut and aligh in retrive band, the shapes are same now
+    # after cut and align in retrive band, the shapes are same now
+    # import ipdb; ipdb.set_trace()
     assert np.shape(bandsdata_a["bands"]) == np.shape(
         bandsdata_b["bands"]
     ), f'{np.shape(bandsdata_a["bands"])} != {np.shape(bandsdata_b["bands"])}'
@@ -275,7 +276,7 @@ def get_bands_distance(
         "eta_c": eta_c,
         "shift_c": shift_c,
         "max_diff_c": max_diff_c,
-        "units": "meV",
+        "unit": "meV",
     }
 
     return out

src/aiida_sssp_workflow/protocol/control.yml

@@ -31,6 +31,14 @@ quick:
     nonnc_dual_scan: [6.0, 7.0, 8.0]
     nonnc_high_dual_scan: [8.0, 10.0, 12.0, 18.0]
 
+experiment:
+    description: dense grid for experiments which goes to SI of the paper
+
+    wfc_scan: [20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170, 175, 180, 185, 190, 195, 200] # at fix dual
+    nc_dual_scan: [2.0, 2.5, 3.0, 3.5, 4.0] # at fix wfc
+    nonnc_dual_scan: [6.0, 6.5, 7.0, 7.5, 8.0]
+    nonnc_high_dual_scan: [8.0, 9.0, 10.0, 12.0, 16.0, 18.0]
+
 test:
     description: test only
 

src/aiida_sssp_workflow/workflows/convergence/bands.py

@@ -4,15 +4,18 @@
 """
 
 from pathlib import Path
-from typing import Union
+from typing import Union, Any
+import copy
 
 from aiida import orm
 from aiida.engine import ProcessBuilder
 from aiida_pseudo.data.pseudo import UpfData
 
 from aiida_sssp_workflow.utils import get_default_mpi_options
+from aiida_sssp_workflow.calculations.calculate_bands_distance import get_bands_distance
 from aiida_sssp_workflow.workflows.convergence._base import _BaseConvergenceWorkChain
 from aiida_sssp_workflow.workflows.evaluate._bands import BandsWorkChain
+from aiida_sssp_workflow.workflows.convergence.report import ConvergenceReport
 
 
 class ConvergenceBandsWorkChain(_BaseConvergenceWorkChain):
@@ -147,3 +150,85 @@ def prepare_evaluate_builder(self, ecutwfc, ecutrho):
         builder.run_band_structure = orm.Bool(False)
 
         return builder
+
+
+def compute_xy(
+    node: orm.Node,
+) -> dict[str, Any]:
+    """From report calculate the xy data, xs are cutoffs and ys are band distance from reference"""
+    report_dict = node.outputs.report.get_dict()
+    report = ConvergenceReport.construct(**report_dict)
+
+    reference_node = orm.load_node(report.reference.uuid)
+    band_structure_r: orm.BandsData = reference_node.outputs.bands.band_structure 
+    band_parameters_r: orm.Dict = reference_node.outputs.bands.band_parameters
+
+    bandsdata_r = {
+        "number_of_electrons": band_parameters_r["number_of_electrons"],
+        "number_of_bands": band_parameters_r["number_of_bands"],
+        "fermi_level": band_parameters_r["fermi_energy"],
+        "bands": band_structure_r.get_bands(),
+        "kpoints": band_structure_r.get_kpoints(),
+        "weights": band_structure_r.get_array("weights"),
+    }
+
+    # smearing width is from degauss
+    smearing = reference_node.inputs.bands.pw.parameters.get_dict()['SYSTEM']['degauss']
+    fermi_shift = reference_node.inputs.fermi_shift.value
+
+    # always do smearing on high bands and not include the spin since we didn't turn on the spin for all
+    # convergence test, but this may change in the future.
+    # The `get_bands_distance` function can deal with mag bands with spin_up and spin_down
+    spin = False
+    do_smearing = True
+
+    xs = []
+    ys = []
+    for node_point in report.convergence_list:
+        if node_point.exit_status != 0:
+            # TODO: log to a warning file for where the node is not finished_okay
+            continue
+
+        x = node_point.wavefunction_cutoff
+        xs.append(x)
+
+        node = orm.load_node(node_point.uuid) 
+        band_structure_p: orm.BandsData = node.outputs.bands.band_structure
+        band_parameters_p: orm.Dict = node.outputs.bands.band_parameters
+
+
+        # The raw implementation of `get_bands_distance` is in `aiida_sssp_workflow/calculations/bands_distance.py`
+        bandsdata_p = {
+            "number_of_electrons": band_parameters_p["number_of_electrons"],
+            "number_of_bands": band_parameters_p["number_of_bands"],
+            "fermi_level": band_parameters_p["fermi_energy"],
+            "bands": band_structure_p.get_bands(),
+            "kpoints": band_structure_p.get_kpoints(),
+            "weights": band_structure_p.get_array("weights"),
+        }
+        res = get_bands_distance(
+            copy.deepcopy(bandsdata_r),
+            copy.deepcopy(bandsdata_p),
+            smearing=smearing,
+            fermi_shift=fermi_shift,
+            do_smearing=do_smearing,
+            spin=spin,
+        )
+        eta_c = res.get("eta_c", None)
+        shift_c = res.get("shift_c", None)
+        max_diff_c = res.get("max_diff_c", None)
+        unit = res.get("unit", None)
+
+        # eta_c is the y, others are write into as metadata
+        ys.append(eta_c)
+
+
+    return {
+        'x': xs,
+        'y': ys,
+        'metadata': {
+            'shift_c': shift_c,
+            'max_diff_c': max_diff_c,
+            'unit': unit,
+        }
+    }

src/aiida_sssp_workflow/workflows/convergence/cohesive_energy.py

@@ -4,13 +4,14 @@
 """
 
 from pathlib import Path
-from typing import Union
+from typing import Union, Any
 
 from aiida import orm
 from aiida.engine import ProcessBuilder
 from aiida_pseudo.data.pseudo import UpfData
 
 from aiida_sssp_workflow.utils import get_default_mpi_options
+from aiida_sssp_workflow.workflows.convergence.report import ConvergenceReport
 from aiida_sssp_workflow.workflows.convergence._base import _BaseConvergenceWorkChain
 from aiida_sssp_workflow.workflows.evaluate._cohesive_energy import (
     CohesiveEnergyWorkChain,
@@ -170,3 +171,39 @@ def prepare_evaluate_builder(self, ecutwfc, ecutrho) -> ProcessBuilder:
         builder.atom.pw["metadata"]["options"] = self.inputs.atom_mpi_options.get_dict()
 
         return builder
+
+def compute_xy(
+    node: orm.Node,
+) -> dict[str, Any]:
+    """From report calculate the xy data, xs are cutoffs and ys are cohesive energy diff from reference"""
+    report_dict = node.outputs.report.get_dict()
+    report = ConvergenceReport.construct(**report_dict)
+
+    reference_node = orm.load_node(report.reference.uuid)
+    output_parameters_r: orm.Dict = reference_node.outputs.output_parameters
+    y_ref = output_parameters_r['cohesive_energy_per_atom']
+
+    xs = []
+    ys = []
+    for node_point in report.convergence_list:
+        if node_point.exit_status != 0:
+            # TODO: log to a warning file for where the node is not finished_okay
+            continue
+
+        x = node_point.wavefunction_cutoff
+        xs.append(x)
+
+        node = orm.load_node(node_point.uuid)
+        output_parameters_p: orm.Dict = node.outputs.output_parameters
+
+        y = abs((output_parameters_p['cohesive_energy_per_atom'] - y_ref) / y_ref) * 100
+        ys.append(y)
+
+    return {
+        'x': xs,
+        'y': ys,
+        'metadata': {
+            'unit': '%',
+        }
+    }
+

src/aiida_sssp_workflow/workflows/convergence/eos.py

@@ -4,14 +4,16 @@
 """
 
 from pathlib import Path
-from typing import Union
+from typing import Union, Any
 
 from aiida import orm
 from aiida.engine import ProcessBuilder
 from aiida_pseudo.data.pseudo import UpfData
 
+from aiida_sssp_workflow.calculations.calculate_metric import rel_errors_vec_length
 from aiida_sssp_workflow.utils import get_default_mpi_options
 from aiida_sssp_workflow.workflows.convergence._base import _BaseConvergenceWorkChain
+from aiida_sssp_workflow.workflows.convergence.report import ConvergenceReport
 from aiida_sssp_workflow.workflows.evaluate._metric import MetricWorkChain
 
 
@@ -127,3 +129,59 @@ def prepare_evaluate_builder(self, ecutwfc, ecutrho) -> ProcessBuilder:
         builder.eos.pw["metadata"]["options"] = self.inputs.mpi_options.get_dict()
 
         return builder
+
+def compute_xy(
+    report: ConvergenceReport,
+) -> dict[str, Any]:
+    """From report calculate the xy data, xs are cutoffs and ys are band distance from reference"""
+    reference_node = orm.load_node(report.reference.uuid)
+    output_parameters_r: orm.Dict = reference_node.outputs.output_parameters
+    ref_V0, ref_B0, ref_B1 = output_parameters_r['birch_murnaghan_results']
+
+
+
+    xs = []
+    ys = []
+    for node_point in report.convergence_list:
+        if node_point.exit_status != 0:
+            # TODO: log to a warning file for where the node is not finished_okay
+            continue
+
+        x = node_point.wavefunction_cutoff
+        xs.append(x)
+
+        node = orm.load_node(node_point.uuid)
+        output_parameters_p: orm.Dict = node.outputs.output_parameters
+
+        V0, B0, B1 = output_parameters_p['birch_murnaghan_results']
+
+        y_nu = rel_errors_vec_length(ref_V0, ref_B0, ref_B1, V0, B0, B1)
+
+        ys.append(y_nu)
+
+    return {
+        'x': xs,
+        'y': ys,
+        'metadata': {
+            'unit': 'n/a',
+        }
+    }
+
+# def compute_xy_epsilon(
+#     report: ConvergenceReport,
+# ) -> dict[str, Any]:
+#     sample_node = orm.load_node(sample_uuid)
+#     ref_node = orm.load_node(ref_uuid)
+#
+#     arr_sample = np.array(sample_node.outputs.eos.output_volume_energy.get_dict()["energies"])
+#     arr_ref = np.array(ref_node.outputs.eos.output_volume_energy.get_dict()["energies"])
+#
+#     avg_sample = np.average(arr_sample)
+#     avg_ref = np.average(arr_ref)
+#
+#     # eq.6 of nat.rev.phys
+#     A = np.sum(np.square(arr_sample - arr_ref))
+#     B = np.sum(np.square(arr_sample - avg_sample))
+#     C = np.sum(np.square(arr_ref - avg_ref))
+#
+#     epsilon = np.sqrt(A / (np.sqrt(B * C)))

src/aiida_sssp_workflow/workflows/convergence/phonon_frequencies.py

@@ -3,7 +3,7 @@
 Convergence test on phonon frequencies of a given pseudopotential
 """
 
-from typing import Union
+from typing import Union, Any
 from pathlib import Path
 
 from aiida import orm
@@ -13,6 +13,7 @@
 from aiida_sssp_workflow.workflows.evaluate._phonon_frequencies import (
     PhononFrequenciesWorkChain,
 )
+from aiida_sssp_workflow.workflows.convergence.report import ConvergenceReport
 from aiida_sssp_workflow.utils import get_default_mpi_options
 
 
@@ -167,3 +168,72 @@ def prepare_evaluate_builder(self, ecutwfc, ecutrho):
         builder.phonon.ph["metadata"]["options"] = self.inputs.ph_mpi_options.get_dict()
 
         return builder
+
+def compute_xy(
+    node: orm.Node,
+) -> dict[str, Any]:
+    """From report calculate the xy data, xs are cutoffs and ys are phonon frequencies diff from reference"""
+    import numpy as np
+
+    report_dict = node.outputs.report.get_dict()
+    report = ConvergenceReport.construct(**report_dict)
+
+    reference_node = orm.load_node(report.reference.uuid)
+    output_parameters_r: orm.Dict = reference_node.outputs.ph.output_parameters
+    y_ref = output_parameters_r["dynamical_matrix_1"]["frequencies"]
+
+    xs = []
+    ys = []
+    for node_point in report.convergence_list:
+        if node_point.exit_status != 0:
+            # TODO: log to a warning file for where the node is not finished_okay
+            continue
+
+        x = node_point.wavefunction_cutoff
+        xs.append(x)
+
+        node = orm.load_node(node_point.uuid)
+        output_parameters_p: orm.Dict = node.outputs.ph.output_parameters
+
+        y_p = output_parameters_p["dynamical_matrix_1"]["frequencies"]
+
+        # calculate the diff
+        diffs = np.array(y_p) - np.array(y_ref)
+        weights = np.array(y_ref)
+
+        relative_diff = np.sqrt(np.mean((diffs / weights) ** 2))
+        y = relative_diff
+
+        # XXX: properties that should write down for GUI rendering
+        # omega_max = np.amax(y_p)
+        # absolute_diff = np.mean(diffs)
+        # absolute_max_diff = np.amax(diffs)
+        #
+        # relative_max_diff = np.amax(diffs / weights)
+
+        # Legacy modification required when configuration is `GS` which is not run for convergence anymore
+        # Keep it here just for reference.
+        # if configuration == "GS":
+        #     # leftover setting from SSSP v1
+        #     # Otherwise the phonon frequencies calculated at BZ boundary qpoint (1/2, 1/2, 1/2) are not converged.
+        #     if element == "N" or element == "Cl":
+        #         start_idx = 12
+        #     elif element == "H" or element == "I":
+        #         start_idx = 4
+        #     elif element == "O":
+        #         start_idx = 6
+        #     else:
+        #         start_idx = 0
+        # else:
+        #     start_idx = 0
+        #
+        ys.append(y)
+
+    return {
+        'x': xs,
+        'y': ys,
+        'metadata': {
+            'unit': '%',
+        }
+    }
+