Measure precision workflow accept ecutwfc/ecutrho as inputs (#214)

* hotfix, report exception in bands

* WIP: measure precision workflow accept ecutwfc/ecutrho as inputs

* CLI: support error throw and warnings when inputs are not proper

Author: unkcpz

README.md

@@ -38,7 +38,12 @@ This what we call it residual volume is therefore a stiffness-agnostic value tha
 In convergence delta factor calculation, the GS configurations from Cottiner's paper are used.
 To have a uniform sence of precision through looking at delta factor, the value is defined as delta factor per atom.
 However, since it is hard to define delta per/atom for oxides, ACWF does not use per/atom value to represent results.
-To compatible with it, the precision verification also use the structure delta.
+It uses the nu factor which was defined in the ACWF paper.
+
+For the oxides, it needs oxygen pseudopotential first, same as nitrides.
+The verification needs to run for all know oxygen/nitrigen pseudopotentials filst. However, the cutoffs of them are not known and the cutoffs are input for the precision measure.
+So the very first run is the convergence test on all oxygen/nitrigen pseudopotentials.
+Once the recommended cutoffs are known, the precision measure can be run for all oxygen/nitrigen elements.
 
 In the convergence verification, there are many options to choose from for what configuration to use. But it is redundant to use all of them, we what the configuration that gives the largest recommended cutoff energy for the pseudopotential.
 After tests (show test results, which in `sssp-verify-scripts`), we found that the for non-magnetic elements, the Diamond configuration that gives the largest cutoff energy. For magnetic elements, the GS configuration calculation with magnetization turned on for the bulk (in cohesive energy which usually give the largest recommended cutoff, where the atomic calculation still don't have magnetization turned on and it is not needed because the pseudopotential is generated in the close shell approximiation (??)) that gives the largest cutoff energy.
@@ -190,4 +195,4 @@ MIT
 
 ## Contact
 
-📧 email: jusong.yu@epfl.ch
+📧 email: jusong.yu@psi.ch

aiida_sssp_workflow/calculations/calculate_bands_distance.py

@@ -171,7 +171,7 @@ def fun_shift(occ, bands_diff, shift):
 def get_bands_distance(
     bandsdata_a: dict,
     bandsdata_b: dict,
-    smearing: float,
+    smearing: float,  # from degauss
     fermi_shift: float,
     do_smearing: bool,
     spin: bool,
@@ -189,6 +189,9 @@ def get_bands_distance(
 
     First aligh the number of two bands, e.g tranctrate the overceed nubmer of bands
     """
+    _RY_TO_EV = 13.6056980659
+    smearing = smearing * _RY_TO_EV
+
     # post process to deserial list to numpy arrar
     for key in ["bands", "kpoints", "weights"]:
         bandsdata_a[key] = np.asarray(bandsdata_a[key])

aiida_sssp_workflow/cli/run.py

@@ -9,6 +9,7 @@
 import click
 from aiida import orm
 from aiida.cmdline.params import options, types
+from aiida.cmdline.utils import echo
 from aiida.engine import run_get_node, submit
 from aiida.plugins import DataFactory, WorkflowFactory
 
@@ -27,12 +28,13 @@
 
 
 @cmd_root.command("launch")
+@click.argument("pseudo", type=click.Path(exists=True))
 @options.OverridableOption(
     "--pw-code", "pw_code", type=types.CodeParamType(entry_point="quantumespresso.pw")
 )(required=True)
 @options.OverridableOption(
     "--ph-code", "ph_code", type=types.CodeParamType(entry_point="quantumespresso.ph")
-)(required=True)
+)(required=False)
 @click.option(
     "--property",
     multiple=True,
@@ -42,22 +44,8 @@
 @click.option(
     "protocol",
     "--protocol",
-    default="standard",
-    help="Protocol to use for the verification.",
-)
-@click.option(
-    "cutoff_control",
-    "--cutoff-control",
-    default="standard",
-    help="Control of convergence.",
-)
-@click.option(
-    "criteria", "--criteria", default="efficiency", help="Criteria of convergence."
-)
-@click.option(
-    "configuration",
-    "--configuration",
-    help="(convergence only) Configuration of structure, can be: SC, FCC, BCC, Diamond, XO, XO2, XO3, X2O, X2O3, X2O5, GS.",
+    default="acwf",
+    help="Protocol to use for the verification, (acwf, test).",
 )
 @click.option("withmpi", "--withmpi", default=True, help="Run with mpi.")
 @click.option("npool", "--npool", default=1, help="Number of pool.")
@@ -78,12 +66,39 @@
     "--comment",
     help="Add a comment word as the prefix the workchain description.",
 )
-@click.argument("pseudo", type=click.Path(exists=True))
+# ecutwfc and ecutrho are for measure workflows only
+@click.option(
+    "ecutwfc",
+    "--ecutwfc",
+    help="Cutoff energy for wavefunctions in Rydberg.",
+)
+@click.option(
+    "ecutrho",
+    "--ecutrho",
+    help="Cutoff energy for charge density in Rydberg.",
+)
+# cutoff_control, criteria, configuration is for convergence workflows only
+@click.option(
+    "cutoff_control",
+    "--cutoff-control",
+    help="Cutoff control for convergence workflow, (standard, quick, opsp).",
+)
+@click.option(
+    "criteria", "--criteria", help="Criteria for convergence (efficiency, precision)."
+)
+# configuration is hard coded for convergence workflow, but here is an interface for experiment purpose
+@click.option(
+    "configuration",
+    "--configuration",
+    help="(convergence test only) Configuration of structure, can be: SC, FCC, BCC, Diamond, XO, XO2, XO3, X2O, X2O3, X2O5, GS, RE",
+)
 def launch(
     pw_code,
     ph_code,
     property,
     protocol,
+    ecutwfc,
+    ecutrho,
     cutoff_control,
     criteria,
     configuration,
@@ -112,8 +127,41 @@ def launch(
         properties_list = list(property)
         extra_desc = f"{properties_list}"
 
+    # validate the options are all provide for the property
+    is_convergence = False
+    is_measure = False
+    is_ph = False
+    for prop in properties_list:
+        if prop.startswith("convergence"):
+            is_convergence = True
+        if prop.startswith("measure"):
+            is_measure = True
+        if "phonon_frequencies" in prop:
+            is_ph = True
+
+    # raise error if the options are not provided
+    if is_convergence and not (cutoff_control and criteria):
+        echo.echo_critical(
+            "cutoff_control, criteria must be provided for convergence workflow."
+        )
+
+    if is_measure and not (ecutwfc and ecutrho):
+        echo.echo_critical("ecutwfc and ecutrho must be provided for measure workflow.")
+
+    if is_ph and not ph_code:
+        echo.echo_critical("ph_code must be provided for phonon frequencies.")
+
+    # raise warning if the options are over provided, e.g. cutoff_control is provided for measure workflow
+    if is_measure and (cutoff_control or criteria or configuration):
+        echo.echo_warning(
+            "cutoff_control, criteria, configuration are not used for measure workflow."
+        )
+
+    if is_convergence and (ecutwfc or ecutrho):
+        echo.echo_warning("ecutwfc and ecutrho are not used for convergence workflow.")
+
     _profile = aiida.load_profile()
-    click.echo(f"Current profile: {_profile.name}")
+    echo.echo_info(f"Current profile: {_profile.name}")
 
     basename = os.path.basename(pseudo)
 
@@ -130,15 +178,15 @@ def launch(
     inputs = {
         "measure": {
             "protocol": orm.Str(protocol),
-            "cutoff_control": orm.Str(cutoff_control),
+            "wavefunction_cutoff": orm.Float(ecutwfc),
+            "charge_density_cutoff": orm.Float(ecutrho),
         },
         "convergence": {
             "protocol": orm.Str(protocol),
             "cutoff_control": orm.Str(cutoff_control),
             "criteria": orm.Str(criteria),
         },
         "pw_code": pw_code,
-        "ph_code": ph_code,
         "pseudo": pseudo,
         "label": label,
         "properties_list": orm.List(properties_list),
@@ -156,6 +204,9 @@ def launch(
         "clean_workdir": orm.Bool(clean_workdir),
     }
 
+    if is_ph:
+        inputs["ph_code"] = ph_code
+
     if configuration is not None:
         inputs["convergence"]["configuration"] = orm.Str(configuration)
 
@@ -167,8 +218,8 @@ def launch(
     description = f"{label.value} ({extra_desc})"
     node.description = f"({comment}) {description}" if comment else description
 
-    click.echo(node)
-    click.echo(f"calculated on property: {'/'.join(properties_list)}")
+    echo.echo_info(node)
+    echo.echo_success(f"calculated on property: {'/'.join(properties_list)}")
 
 
 if __name__ == "__main__":

aiida_sssp_workflow/protocol/control.yml

@@ -2,7 +2,6 @@
 precheck:
     description: run on 150 v.s 200 v.s 300 Ry wavefunction cutoff with fix dual to give sense whether 200 Ry enough
 
-    # max_wfc: 0 # Not used if run precision measurement veri with this protocol, just raise exception.
     wfc_scan: [150, 200, 300] # at fix dual
     nc_dual_scan: [] # set empty so rho test not run
     nonnc_dual_scan: []
@@ -11,7 +10,6 @@ precheck:
 standard:
     description: high wavefunction cutoff set and cutoffs dense interval therefore time consuming
 
-    max_wfc: 200 # The max wfc cut for precision measurement
     wfc_scan: [30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 90, 100, 120, 150, 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]
@@ -20,7 +18,6 @@ standard:
 quick:
     description: low wavefunction cutoff set and cutoffs sparse interval therefore can run on local
 
-    max_wfc: 200 # The max wfc cut for precision measurement
     wfc_scan: [30, 40, 50, 60, 75, 100, 150, 200] # at fix dual
     nc_dual_scan: [3.0, 3.5, 4.0] # at fix wfc
     nonnc_dual_scan: [6.0, 7.0, 8.0]
@@ -29,7 +26,6 @@ quick:
 test:
     description: test only
 
-    max_wfc: 35 # The max wfc cut for precision measurement
     wfc_scan: [30, 35] # at fix dual
     nc_dual_scan: [] # at fix wfc
     nonnc_dual_scan: [] # at fix rho
@@ -38,7 +34,6 @@ test:
 opsp:
     description: To running opsp verification, make the calculations finished fast.
 
-    max_wfc: 40 # The max wfc cut for precision measurement
     wfc_scan: [30, 35] # at fix dual
     nc_dual_scan: [2.0, 4.0] # at fix wfc
     nonnc_dual_scan: [6.0, 8.0]

aiida_sssp_workflow/protocol/converge.yml

@@ -37,42 +37,42 @@ acwf:
         mixing_beta: 0.4
 
 
-moderate:
-    name: moderate
-    description: The protocol where input parameters bring from aiidaqe moderate protocol. Only for QE >= 6.8
-
-    base:   # base parameters is inherit by other process
-        occupations: smearing
-        degauss: 0.01
-        smearing: cold
-        conv_thr_per_atom: 1.0e-10
-        kpoints_distance: 0.15
-
-    cohesive_energy:
-        atom_smearing: gaussian
-        vacuum_length: 12.0
-
-    phonon_frequencies:
-        qpoints_list:
-            - [0.5, 0.5, 0.5]
-        epsilon: false
-        tr2_ph: 1.0e-16
-
-    pressure:
-        scale_count: 7
-        scale_increment: 0.02
-        mixing_beta: 0.4
-
-    bands:
-        init_nbands_factor: 3.0
-        fermi_shift: 10.0
-        kpoints_distance_scf: 0.15
-        kpoints_distance_bands: 0.15
-
-    delta:
-        scale_count: 7
-        scale_increment: 0.02
-        mixing_beta: 0.4
+#moderate:
+#    name: moderate
+#    description: The protocol where input parameters bring from aiidaqe moderate protocol. Only for QE >= 6.8
+#
+#    base:   # base parameters is inherit by other process
+#        occupations: smearing
+#        degauss: 0.01
+#        smearing: cold
+#        conv_thr_per_atom: 1.0e-10
+#        kpoints_distance: 0.15
+#
+#    cohesive_energy:
+#        atom_smearing: gaussian
+#        vacuum_length: 12.0
+#
+#    phonon_frequencies:
+#        qpoints_list:
+#            - [0.5, 0.5, 0.5]
+#        epsilon: false
+#        tr2_ph: 1.0e-16
+#
+#    pressure:
+#        scale_count: 7
+#        scale_increment: 0.02
+#        mixing_beta: 0.4
+#
+#    bands:
+#        init_nbands_factor: 3.0
+#        fermi_shift: 10.0
+#        kpoints_distance_scf: 0.15
+#        kpoints_distance_bands: 0.15
+#
+#    delta:
+#        scale_count: 7
+#        scale_increment: 0.02
+#        mixing_beta: 0.4
 
 test:
     name: test-only

aiida_sssp_workflow/workflows/convergence/bands.py

@@ -74,8 +74,6 @@ class ConvergenceBandsWorkChain(_BaseConvergenceWorkChain):
 
     # pylint: disable=too-many-instance-attributes
 
-    _RY_TO_EV = 13.6056980659
-
     _PROPERTY_NAME = "bands"
     _EVALUATE_WORKCHAIN = BandsWorkChain
     _MEASURE_OUT_PROPERTY = "eta_c"
@@ -186,7 +184,7 @@ def helper_compare_result_extract_fun(self, sample_node, reference_node, **kwarg
             sample_band_parameters,
             reference_band_structure,
             reference_band_parameters,
-            smearing=orm.Float(self.ctx.degauss * self._RY_TO_EV),
+            smearing=orm.Float(self.ctx.degauss),
             fermi_shift=orm.Float(self.ctx.fermi_shift),
             do_smearing=orm.Bool(True),
             spin=orm.Bool(spin),

aiida_sssp_workflow/workflows/measure/__init__.py

@@ -9,8 +9,17 @@
 
 UpfData = DataFactory("pseudo.upf")
 
-
 class _BaseMeasureWorkChain(SelfCleanWorkChain):
+    """Base Measure Workchain since bands measure and precision measure share same input ports
+    """
+    # ECUT for oxygen, remember to update this if the oxygen pseudo is changed
+    _O_ECUTWFC = 75.0
+    _O_ECUTRHO = 600.0
+
+    # ECUT for nitrogen, remember to update this if the nitrogen pseudo is changed
+    _N_ECUTWFC = 80.0
+    _N_ECUTRHO = 320.0
+
     @classmethod
     def define(cls, spec):
         """Define the process specification."""
@@ -22,9 +31,24 @@ def define(cls, spec):
                     help='Pseudopotential to be verified')
         spec.input('protocol', valid_type=orm.Str, required=True,
                     help='The protocol which define input calculation parameters.')
-        spec.input('cutoff_control', valid_type=orm.Str, default=lambda: orm.Str('standard'),
-                    help='The control protocol where define max_wfc.')
+        spec.input('wavefunction_cutoff', valid_type=orm.Float, required=True, help='The wavefunction cutoff.')
+        spec.input('charge_density_cutoff', valid_type=orm.Float, required=True, help='The charge density cutoff.')
         spec.input('options', valid_type=orm.Dict, required=True,
                     help='Optional `options` to use for the `PwCalculations`.')
         spec.input('parallelization', valid_type=orm.Dict, required=True,
                     help='Parallelization options for the `PwCalculations`.')
+
+    def _get_pw_cutoff(self, structure: orm.StructureData, ecutwfc: float, ecutrho: float):
+        """Get cutoff pair, if strcture contains oxygen or nitrogen, need
+        to use the max between pseudo cutoff and the O/N cutoff.
+        """
+        elements = set(structure.get_symbols_set())
+        if "O" in elements:
+            ecutwfc = max(ecutwfc, self._O_ECUTWFC)
+            ecutrho = max(ecutrho, self._O_ECUTRHO)
+
+        if "N" in elements:
+            ecutwfc = max(ecutwfc, self._N_ECUTWFC)
+            ecutrho = max(ecutrho, self._N_ECUTRHO)
+
+        return ecutwfc, ecutrho