threshold for the convergence of the observable of interest
+
str/float
+
+
+
+
variables
+
The observable of interest will be converged with respect to this (these) simulation variable(s)
+
list/str
+
['ecutwfc']
+
+
\ No newline at end of file
diff --git a/docs/_static/keywords/generate_keywords.py b/docs/_static/keywords/generate_keywords.py
index 688a8bd53..66248003e 100644
--- a/docs/_static/keywords/generate_keywords.py
+++ b/docs/_static/keywords/generate_keywords.py
@@ -1,4 +1,5 @@
-from koopmans.settings import (MLSettingsDict, PlotSettingsDict,
+from koopmans.settings import (ConvergenceSettingsDict, MLSettingsDict,
+ PlotSettingsDict,
UnfoldAndInterpolateSettingsDict,
WorkflowSettingsDict)
@@ -53,3 +54,6 @@ def print_settings_as_html(settings, table_id, fd):
with open('ml_keywords.html', 'w') as fd:
print_settings_as_html(MLSettingsDict(), 'mlTable', fd)
+
+ with open('convergence_keywords.html', 'w') as fd:
+ print_settings_as_html(ConvergenceSettingsDict(), 'convTable', fd)
diff --git a/docs/_static/keywords/workflow_keywords.html b/docs/_static/keywords/workflow_keywords.html
index dca792272..b56079ddb 100644
--- a/docs/_static/keywords/workflow_keywords.html
+++ b/docs/_static/keywords/workflow_keywords.html
@@ -38,7 +38,7 @@
pseudo_directory
the folder containing the pseudopotentials to use (mutually exclusive with "pseudo_library")
-
str/Path
+
Path
@@ -131,6 +131,12 @@
float/list
0.6
+
+
alpha_mixing
+
mixing parameter for updating alpha
+
float
+
1.0
+
alpha_from_file
if True, uses the file_alpharef.txt from the base directory as a starting guess
@@ -162,22 +168,22 @@
-
convergence_observable
-
System observable of interest which we converge
-
str
-
total energy
+
orbital_groups_spread_tol
+
when calculating alpha parameters, the code will group orbitals together only if their spread is within this threshold
+
float
+
-
convergence_threshold
-
Convergence threshold for the system observable of interest
-
str/float
-
+
converge
+
If True, repeat the workflow increasing the convergence_parameters until the convergence_observable converges within the convergence_threshold
+
bool
+
False
-
convergence_parameters
-
The observable of interest will be converged with respect to this/these simulation parameter(s)
-
list/str
-
['ecutwfc']
+
dfpt_coarse_grid
+
The coarse k-point grid on which to perform the DFPT calculations
+
list
+
eps_cavity
diff --git a/docs/input_file.rst b/docs/input_file.rst
index 33fe7c021..df62cd21e 100644
--- a/docs/input_file.rst
+++ b/docs/input_file.rst
@@ -163,6 +163,15 @@ This subblock controls the unfolding and interpolation procedure for generating
List of valid keywords
+The convergence block
+^^^^^^^^^^^^^^^^^^^^^
+This block can be used to customize a convergence calculation.
+
+.. toctree::
+
+ List of valid keywords
+
+See :ref:`here for a more detailed tutorial on performing convergence calculations `.
The plotting block
^^^^^^^^^^^^^^^^^^
diff --git a/docs/input_file/convergence_keywords.rst b/docs/input_file/convergence_keywords.rst
new file mode 100644
index 000000000..f6a1b312c
--- /dev/null
+++ b/docs/input_file/convergence_keywords.rst
@@ -0,0 +1,15 @@
+Valid keywords
+==============
+
+.. raw:: html
+
+
+
+
+
+
+.. raw:: html
+ :file: ../_static/keywords/convergence_keywords.html
+
+.. raw:: html
+ :file: ../_static/keywords/lookup_table.html
diff --git a/docs/tutorials/tutorial_4.rst b/docs/tutorials/tutorial_4.rst
index 4f37c3af9..8d5c3a236 100644
--- a/docs/tutorials/tutorial_4.rst
+++ b/docs/tutorials/tutorial_4.rst
@@ -1,30 +1,40 @@
Tutorial 4: running convergence tests
=====================================
-While ``koopmans`` is a package primarily oriented towards performing Koopmans spectral functional calculations, it does have a couple of other useful functionalities. In this tutorial, we will make use of its ``convergence`` task to determine how large a cell size and energy cutoff is required to converge the PBE energy of the highest occupied molecular orbital (HOMO) of a water molecule.
+While ``koopmans`` is a package primarily oriented towards performing Koopmans functional calculations, it does have a couple of other useful functionalities. Among these functionalities is the ability to perform arbitrary covnergence tests.
+
+
+A simple convergence test
+-------------------------
+
+In this tutorial, we will make use of its ``convergence`` task to determine how large a cell size and energy cutoff is required to converge the PBE energy of the highest occupied molecular orbital (HOMO) of a water molecule. This workflow was chosen for its simplicity; it is possible to run convergence tests on any workflow implemented in the ``koopmans`` package.
The input file
---------------
-In order to run this calculation, our ``workflow`` block needs a few particular keywords:
+''''''''''''''
+In order to run this calculation, in the ``workflow`` block we need to set the ``converge`` parameter to true:
.. literalinclude:: ../../tutorials/tutorial_4/h2o_conv.json
- :lines: 1-11
- :linenos:
- :emphasize-lines: 4, 6-10
+ :lines: 1-9
+ :emphasize-lines: 6
-The important lines are highlighted. ``"task": "convergence"`` means that we will be performing a convergence test. The other three highlighted keywords specifying that we are going to converge the HOMO energy to within 0.01 eV, with respect to *both* the energy cutoff ``ecutwfc`` and the cell size. The full input file can be found :download:`here <../../tutorials/tutorial_4/h2o_conv.json>`.
+and then provide the convergence information in the ``convergence`` block:
+
+.. literalinclude:: ../../tutorials/tutorial_4/h2o_conv.json
+ :lines: 10-17
+
+These settings state that we are going to converge the HOMO energy to within 0.01 eV, with respect to *both* the energy cutoff ``ecutwfc`` and the cell size. The full input file can be found :download:`here <../../tutorials/tutorial_4/h2o_conv.json>`.
The output file
----------------
+'''''''''''''''
When you run the calculation, you should see something like this after the header:
-.. literalinclude:: ../../tutorials/tutorial_4/h2o_conv.out
- :lines: 14-32
+.. literalinclude:: ../../tutorials/tutorial_4/h2o_conv.log
+ :lines: 15-34
:language: text
-Here, the code is attempting to use progressively larger energy cutoffs and cell sizes. It will ultimately arrive at a converged solution, with a ``ecutwfc`` of 50.0 Ha and a cell 1.3 times larger than that provided in the ``.json`` input file.
+Here, the code is attempting to use progressively larger energy cutoffs and cell sizes. It will ultimately arrive at a converged solution, with a ``ecutwfc`` of 50.0 Ha and a cell slightly larger than that provided in the ``.json`` input file.
Plotting
---------
+''''''''
The individual ``Quantum ESPRESSO`` calculations reside in nested subdirectories. If you plot the HOMO energies from each of these, you should get something like this:
@@ -35,4 +45,30 @@ The individual ``Quantum ESPRESSO`` calculations reside in nested subdirectories
Plot of the HOMO energy of water with respect to the energy cutoff and cell size (generated using :download:`this script <../../tutorials/tutorial_4/plot.py>`)
-We can see that indeed the calculation with ``ecutwfc = 50.0`` and ``cell_size = 1.3`` is the one where the energy of the HOMO goes within (and stays within) 0.01 eV of the most accurate calculation.
+We can see that indeed the calculation with ``ecutwfc = 50.0`` and ``celldm(1) = 13.3`` is the one where the energy of the HOMO goes within (and stays within) 0.01 eV of the most accurate calculation.
+
+A custom convergence test
+-------------------------
+
+In the previous example, we performed a convergence test with respect to ``ecutwfc`` and ``celldm1``. A full list of supported convergence variables can be found :ref:`here `. You will see that only a couple of variables are implemented by default. But don't let that limit you! It is possible to perform a convergence test on arbitrary keywords using factories.
+
+First, try taking the input file from the first part of the tutorial, and switch the pseudopotential library to ``pseudo_dojo_standard``. What do you notice?
+
+Hopefully, the first thing you will see is that there are now some warnings about the small box parameters ``nrb`` like this,
+
+.. code-block:: text
+
+ UserWarning: Small box parameters "nrb" not provided in input: these will be automatically set to safe default values. These values can probably be decreased, but this would require convergence tests.
+ Estimated real mesh dimension (nr1, nr2, nr3) = ...
+ Small box mesh dimension (nr1b, nr2b, nr3b) = ...
+
+These parameters are associated with the way ``Quantum ESPRESSO`` handles non-local core corrections in pseudopotentials, and corrections are present in the new set of pseudopotentials but absent in the SG15 pseudopotentials.
+
+So, let's perform a convergence test! The ``nrb`` parameters are not directly implemented as a convergence variable in ``koopmans``, but we can use a factory to perform a convergence test on them, by making use of the ``ConvergenceVariable`` and ``ConvergenceWorkflowfactory`` classes.
+
+.. literalinclude:: ../../tutorials/tutorial_4/custom_convergence.py
+
+Running this script will perform a convergence test with respect to ``nrb 1-3``.
+
+.. warning::
+ This tutorial performs convergence tests in a slightly incorrect way, To see this, add the keyword ``length = 10`` to the ``ConvergenceVarialbe`` in the above script. How is the behaviour different? Which behaviour is correct? Why?
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index 8df72c280..38539dfd8 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
[project]
name = "koopmans"
-version = "1.0.0"
+version = "1.0.1"
description = "Koopmans spectral functional calculations with python and Quantum ESPRESSO"
readme = "README_pypi.rst"
requires-python = ">=3.8"
diff --git a/src/koopmans/cell.py b/src/koopmans/cell.py
index 8de5c9d61..8159d37ae 100644
--- a/src/koopmans/cell.py
+++ b/src/koopmans/cell.py
@@ -1,5 +1,5 @@
import copy
-from typing import Dict, Union
+from typing import Dict, TypedDict, Union
import numpy as np
from ase.cell import Cell
@@ -8,6 +8,11 @@
from ase.units import Bohr
+class CellParams(TypedDict):
+ ibrav: int
+ celldms: Dict[int, float]
+
+
def parameters_to_cell(ibrav: int, celldms: Dict[int, float]) -> Cell:
'''
Converts QE cell parameters to the corresponding ASE `Cell` object.
@@ -129,7 +134,7 @@ def parameters_to_cell(ibrav: int, celldms: Dict[int, float]) -> Cell:
return Cell(new_array)
-def cell_to_parameters(cell: Cell) -> Dict[str, Union[int, Dict[int, float]]]:
+def cell_to_parameters(cell: Cell) -> CellParams:
'''
Identifies a cell in the language of Quantum ESPRESSO
@@ -140,8 +145,8 @@ def cell_to_parameters(cell: Cell) -> Dict[str, Union[int, Dict[int, float]]]:
Returns
-------
- Dict
- a dictionary containing the ibrav and celldms
+ CellParams
+ a typed dictionary containing the ibrav and celldms
'''
@@ -254,7 +259,7 @@ def cell_to_parameters(cell: Cell) -> Dict[str, Union[int, Dict[int, float]]]:
# Convert to Bohr radii
celldms[1] /= Bohr
- return {'ibrav': ibrav, 'celldms': celldms}
+ return CellParams(ibrav=ibrav, celldms=celldms)
def cell_follows_qe_conventions(cell: Cell) -> bool:
diff --git a/src/koopmans/io/_json.py b/src/koopmans/io/_json.py
index a4d292b85..392657bf9 100644
--- a/src/koopmans/io/_json.py
+++ b/src/koopmans/io/_json.py
@@ -34,8 +34,6 @@ def read_json(fd: TextIO, override: Dict[str, Any] = {}):
workflow_cls: Type[workflows.Workflow]
if parameters.task == 'singlepoint':
workflow_cls = workflows.SinglepointWorkflow
- elif parameters.task == 'convergence':
- workflow_cls = workflows.ConvergenceWorkflow
elif parameters.task == 'wannierize':
workflow_cls = workflows.WannierizeWorkflow
elif parameters.task == 'environ_dscf':
diff --git a/src/koopmans/io/_kwf.py b/src/koopmans/io/_kwf.py
index 4067a3105..f8bad1b8d 100644
--- a/src/koopmans/io/_kwf.py
+++ b/src/koopmans/io/_kwf.py
@@ -6,6 +6,7 @@
"""
+import functools
import inspect
import json
import os
@@ -16,6 +17,7 @@
from ase.io import jsonio as ase_json
import koopmans.workflows as workflows
+from koopmans import utils
class KoopmansEncoder(ase_json.MyEncoder):
@@ -26,6 +28,14 @@ def default(self, obj) -> dict:
return {'__path__': os.path.relpath(obj, '.')}
elif inspect.isclass(obj):
return {'__class__': {'__name__': obj.__name__, '__module__': obj.__module__}}
+ elif isinstance(obj, functools.partial):
+ return {'__partial__': {'func': obj.func, 'args': obj.args, 'keywords': obj.keywords}}
+ elif inspect.isfunction(obj):
+ module = import_module(obj.__module__)
+ if hasattr(module, obj.__name__):
+ return {'__function__': {'__name__': obj.__name__, '__module__': obj.__module__}}
+ else:
+ return {'__function__': {}}
elif hasattr(obj, 'todict'):
d = obj.todict()
if '__koopmans_name__' in d:
@@ -49,6 +59,17 @@ def object_hook(dct):
subdct = dct['__class__']
module = import_module(subdct['__module__'])
return getattr(module, subdct['__name__'])
+ elif '__partial__' in dct:
+ subdct = dct['__partial__']
+ return functools.partial(subdct['func'], **subdct['keywords'])
+ elif '__function__' in dct:
+ subdct = dct['__function__']
+ if '__module__' in subdct:
+ module = import_module(subdct['__module__'])
+ if hasattr(module, subdct['__name__']):
+ return getattr(module, subdct['__name__'])
+ utils.warn('Function was not able to be serialized')
+ return f''
else:
# Patching bug in ASE where allocating an np.empty(dtype=str) will assume a particular length for each
# string. dtype=object allows for individual strings to be different lengths
diff --git a/src/koopmans/settings/__init__.py b/src/koopmans/settings/__init__.py
index 33471f24c..297f8a82e 100644
--- a/src/koopmans/settings/__init__.py
+++ b/src/koopmans/settings/__init__.py
@@ -6,6 +6,7 @@
'''
+from ._convergence import ConvergenceSettingsDict
from ._koopmans_cp import KoopmansCPSettingsDict
from ._koopmans_ham import KoopmansHamSettingsDict
from ._koopmans_screen import KoopmansScreenSettingsDict
diff --git a/src/koopmans/settings/_convergence.py b/src/koopmans/settings/_convergence.py
new file mode 100644
index 000000000..d5b8dc03b
--- /dev/null
+++ b/src/koopmans/settings/_convergence.py
@@ -0,0 +1,36 @@
+import os
+from pathlib import Path
+from typing import Any
+
+from koopmans import pseudopotentials
+
+from ._utils import Setting, SettingsDictWithChecks
+
+
+class ConvergenceSettingsDict(SettingsDictWithChecks):
+
+ def __init__(self, **kwargs) -> None:
+ settings = [
+ Setting('observable',
+ 'System observable of interest which we converge',
+ str, 'total energy', None),
+ Setting('threshold',
+ 'threshold for the convergence of the observable of interest',
+ (str, float), None, None),
+ Setting('variables',
+ 'The observable of interest will be converged with respect to this (these) '
+ 'simulation variable(s)',
+ (list, str), ['ecutwfc'], None)]
+
+ super().__init__(settings=settings, physicals=['threshold'], **kwargs)
+
+ @property
+ def _other_valid_keywords(self):
+ return []
+
+ def __setitem__(self, key: str, value: Any):
+ # Convert parameters to a list
+ if key == 'parameters' and isinstance(value, str):
+ value = [value]
+
+ return super().__setitem__(key, value)
diff --git a/src/koopmans/settings/_workflow.py b/src/koopmans/settings/_workflow.py
index 01505772d..80da2eb47 100644
--- a/src/koopmans/settings/_workflow.py
+++ b/src/koopmans/settings/_workflow.py
@@ -108,16 +108,10 @@ def __init__(self, **kwargs) -> None:
'when calculating alpha parameters, the code will group orbitals '
'together only if their spread is within this threshold',
float, None, None),
- Setting('convergence_observable',
- 'System observable of interest which we converge',
- str, 'total energy', None),
- Setting('convergence_threshold',
- 'Convergence threshold for the system observable of interest',
- (str, float), None, None),
- Setting('convergence_parameters',
- 'The observable of interest will be converged with respect to this/these '
- 'simulation parameter(s)',
- (list, str), ['ecutwfc'], None),
+ Setting('converge',
+ 'If True, repeat the workflow increasing the convergence_parameters until the '
+ 'convergence_observable converges within the convergence_threshold',
+ bool, False, (True, False)),
Setting('dfpt_coarse_grid',
'The coarse k-point grid on which to perform the DFPT calculations',
list, None, None),
@@ -128,7 +122,7 @@ def __init__(self, **kwargs) -> None:
# Defer storing init_empty_orbitals...
init_empty_orbitals = kwargs.pop('init_empty_orbitals', 'same')
- super().__init__(settings=settings, physicals=['alpha_conv_thr', 'convergence_threshold'], **kwargs)
+ super().__init__(settings=settings, physicals=['alpha_conv_thr'], **kwargs)
# ... until we are sure that init_orbitals has been defined
self.init_empty_orbitals = init_empty_orbitals
diff --git a/src/koopmans/utils/_io.py b/src/koopmans/utils/_io.py
index a9cdd9712..557ab5a6f 100644
--- a/src/koopmans/utils/_io.py
+++ b/src/koopmans/utils/_io.py
@@ -49,7 +49,7 @@ def parse_dict(dct: Dict[str, Any]) -> Dict[str, Any]:
def construct_cell_parameters_block(atoms: Atoms) -> Dict[str, Any]:
params: Dict[str, Any]
if cell_follows_qe_conventions(atoms.cell):
- params = cell_to_parameters(atoms.cell)
+ params = dict(**cell_to_parameters(atoms.cell))
else:
params = {'vectors': [list(row) for row in atoms.cell[:]], 'units': 'angstrom'}
params['periodic'] = all(atoms.pbc)
diff --git a/src/koopmans/workflows/__init__.py b/src/koopmans/workflows/__init__.py
index 4f6f4e120..e87425865 100644
--- a/src/koopmans/workflows/__init__.py
+++ b/src/koopmans/workflows/__init__.py
@@ -1,7 +1,10 @@
'Workflows for use with koopmans'
from ._anion_dscf import DeltaSCFWorkflow
-from ._convergence import ConvergenceWorkflow
+from ._convergence import (ConvergenceVariable, ConvergenceVariableFactory,
+ ConvergenceWorkflow, ConvergenceWorkflowFactory,
+ ObservableFactory, get_calculator_parameter,
+ set_calculator_parameter)
from ._convergence_ml import ConvergenceMLWorkflow
from ._dft import DFTBandsWorkflow, DFTCPWorkflow, DFTPhWorkflow, DFTPWWorkflow
from ._folding import FoldToSupercellWorkflow
diff --git a/src/koopmans/workflows/_convergence.py b/src/koopmans/workflows/_convergence.py
index da42917c2..d052611de 100644
--- a/src/koopmans/workflows/_convergence.py
+++ b/src/koopmans/workflows/_convergence.py
@@ -7,40 +7,77 @@
"""
+from __future__ import annotations
+
import copy
import itertools
import shutil
-from dataclasses import dataclass, field
+from dataclasses import dataclass
from functools import partial
from pathlib import Path
-from typing import Any, Callable, Generic, List, Optional, TypeVar, Union, cast
+from typing import (Any, Callable, Dict, Generic, List, Optional, Type,
+ TypeVar, Union, cast)
import numpy as np
-from koopmans import utils
+from koopmans import cell, utils
from ._workflow import Workflow
-T = TypeVar('T')
+# T is either an int or a float or a list of ints/floats
+T = TypeVar('T', int, float, List[int], List[float])
+
+
+def generate_values_via_addition(initial_value: T, increment: T, length: int) -> List[T]:
+ if isinstance(initial_value, list):
+ return [[v + i * delta for v, delta in zip(initial_value, increment)] for i in range(length)]
+ else:
+ return [initial_value + i * increment for i in range(length)]
@dataclass
-class ConvergenceParameter(Generic[T]):
+class ConvergenceVariable(Generic[T]):
name: str
increment: T
get_value: Callable[[Workflow], T]
set_value: Callable[[Workflow, T], None]
- values: List[T] = field(default_factory=list)
+ initial_value: Optional[T] = None
+ length: int = 3
+ generate_values: Callable[[T, T, int], List[T]] = generate_values_via_addition
converged_value: Optional[T] = None
def extend(self):
- if isinstance(self.values[-1], list):
- self.values.append([x + self.increment for x in self.values[-1]])
- else:
- self.values.append(self.values[-1] + self.increment)
+ self.length += 1
+
+ def __len__(self) -> int:
+ return self.length
+
+ @property
+ def values(self) -> List[T]:
+ if self.initial_value is None:
+ raise ValueError(
+ f'{self.__class__.__name__}.initial_value has not been set. Use {self.__class__.__name__}.get_initial_value() to set it.')
+
+ return self.generate_values(self.initial_value, self.increment, self.length)
+
+ def get_initial_value(self, workflow: Workflow):
+ if self.initial_value is not None:
+ raise ValueError(
+ f'Do not call {self.__class__.__name__}.get_initial_value() once {self.__class__.__name__}.initial_value has been initialized')
+ self.initial_value = self.get_value(workflow)
+
+ def todict(self) -> Dict[str, Any]:
+ dct = dict(self.__dict__)
+ dct['__koopmans_name__'] = self.__class__.__name__
+ dct['__koopmans_module__'] = self.__class__.__module__
+ return dct
+ @classmethod
+ def fromdict(cls, dct: Dict[str, Any]):
+ return cls(**dct)
-def get_calc_value(wf: Workflow, key: str) -> None:
+
+def get_calculator_parameter(wf: Workflow, key: str) -> None:
'''
Gets a particular calculator setting from a workflow
'''
@@ -49,7 +86,8 @@ def get_calc_value(wf: Workflow, key: str) -> None:
if settings.is_valid(key):
val = settings[key]
if val is None:
- raise AttributeError(f'In order to converge wrt {key}, specify a baseline value for it')
+ raise AttributeError(
+ f'In order to converge wrt {key}, specify a baseline value for it')
values.add(val)
if len(values) > 1:
raise ValueError(f'{key} has different values for different calculators. This is not compatible with the'
@@ -57,7 +95,7 @@ def get_calc_value(wf: Workflow, key: str) -> None:
return values.pop()
-def set_calc_value(wf: Workflow, value: Union[int, float], key: str) -> None:
+def set_calculator_parameter(wf: Workflow, key: str, value: Union[int, float]) -> None:
'''
Sets a particular calculator setting for every calculator in the workflow
'''
@@ -80,124 +118,204 @@ def fetch_eps_inf(wf: Workflow) -> float:
return np.mean(np.diag(calc.results['dielectric tensor']))
-class ConvergenceWorkflow(Workflow):
+def ObservableFactory(obs_str: str) -> Callable[[Workflow], float]:
+ if obs_str == 'eps_inf':
+ return fetch_eps_inf
+ else:
+ return partial(fetch_result_default, observable=obs_str)
- def _run(self, initial_depth: int = 3) -> None:
- # Deferred import to allow for monkeypatchings
- from koopmans import workflows
+def _get_celldm1(wf: Workflow) -> float:
- convergence_parameters = []
+ params = cell.cell_to_parameters(wf.atoms.cell)
+ assert 1 in params['celldms']
+ return params['celldms'][1]
- for conv_param in self.parameters.convergence_parameters:
- get_value: Callable[[Workflow], Any]
- set_value: Callable[[Workflow, Any], None]
- values: List[Any]
+def _set_celldm1(wf: Workflow, value: float) -> None:
+ params = cell.cell_to_parameters(wf.atoms.cell)
+ params['celldms'][1] = value
+ wf.atoms.cell = cell.parameters_to_cell(**params)
- if conv_param == 'cell_size':
- increment = 0.1
- def get_value(wf: Workflow) -> float:
- mask = self.atoms.cell != 0
- cell_ratio = wf.atoms.cell[mask] / self.atoms.cell[mask]
- return list(set(cell_ratio))[0]
+def conv_var_celldm1(increment: float = 1.0, **kwargs) -> ConvergenceVariable:
- def set_value(wf: Workflow, value: float) -> None:
- wf.atoms.cell = self.atoms.cell * value
- values = [get_value(self) + i * increment for i in range(initial_depth)]
+ return ConvergenceVariable('celldm1', increment, _get_celldm1, _set_celldm1, **kwargs)
- elif conv_param == 'ecutwfc':
- increment = 10
- def set_value(wf: Workflow, value: float) -> None:
- set_calc_value(wf, value, key='ecutwfc')
- set_calc_value(wf, 4 * value, key='ecutrho')
- get_value = cast(Callable[[Workflow], float], partial(get_calc_value, key='ecutwfc'))
- values = [get_value(self) + i * increment for i in range(initial_depth)]
- elif conv_param == 'nbnd':
- increment = 1
- set_value = cast(Callable[[Workflow, int], None], partial(set_calc_value, key='nbnd'))
- get_value = cast(Callable[[Workflow], int], partial(get_calc_value, key='nbnd'))
- values = [get_value(self) + i * increment for i in range(initial_depth)]
+def _set_ecutwfc(wf: Workflow, value: float) -> None:
+ set_calculator_parameter(wf, 'ecutwfc', value)
+ set_calculator_parameter(wf, 'ecutrho', 4 * value)
- elif conv_param == 'kgrid':
- increment = 2
- def get_value(wf: Workflow) -> List[int]:
- assert wf.kpoints.grid
- return wf.kpoints.grid
+_get_ecutwfc = cast(Callable[[Workflow], float], partial(get_calculator_parameter, key='ecutwfc'))
- def set_value(wf: Workflow, value: List[int]) -> None:
- wf.kpoints.grid = value
- values = [[x + i * increment for x in get_value(self)] for i in range(initial_depth)]
- else:
- raise NotImplementedError(f'Convergence with respect to {conv_param} has not yet been implemented')
- convergence_parameters.append(ConvergenceParameter(conv_param, increment, get_value, set_value, values))
+def conv_var_ecutwfc(increment: float = 10.0, **kwargs) -> ConvergenceVariable:
+
+ return ConvergenceVariable('ecutwfc', increment, _get_ecutwfc, _set_ecutwfc, **kwargs)
+
+
+_set_nbnd = cast(Callable[[Workflow, int], None], partial(set_calculator_parameter, key='nbnd'))
+_get_nbnd = cast(Callable[[Workflow], int], partial(get_calculator_parameter, key='nbnd'))
+
+
+def conv_var_nbnd(increment: int = 1, **kwargs) -> ConvergenceVariable:
+
+ return ConvergenceVariable('nbnd', increment, _get_nbnd, _set_nbnd, **kwargs)
+
+
+def _get_kgrid(wf: Workflow) -> List[int]:
+ assert wf.kpoints.grid
+ return wf.kpoints.grid
+
+
+def _set_kgrid(wf: Workflow, value: List[int]) -> None:
+ wf.kpoints.grid = value
+
+
+def _generate_kgrid_values(initial_value: List[int], increment: List[int], length: int) -> List[List[int]]:
+ values: List[List[int]] = [initial_value]
+
+ while len(values) < length:
+ candidates = []
+ for delta in itertools.product([0, 1], repeat=len(initial_value)):
+ if 1 not in delta:
+ continue
+ candidates.append([x + d * i for x, d, i in zip(values[-1], delta, increment)])
+
+ def cos_angle(c):
+ return np.dot(c, initial_value) / np.linalg.norm(c) / np.linalg.norm(initial_value)
+ best_candidate = sorted(candidates, key=cos_angle)[-1]
+ values.append(best_candidate)
+
+ return values
+
+
+def conv_var_kgrid(increment: List[int] = [1, 1, 1], **kwargs) -> ConvergenceVariable:
- if self.parameters.convergence_observable == 'eps_inf':
- fetch_result = fetch_eps_inf
- else:
- fetch_result = partial(fetch_result_default, observable=self.parameters.convergence_observable)
+ return ConvergenceVariable('kgrid', increment, _get_kgrid, _set_kgrid,
+ generate_values=_generate_kgrid_values, **kwargs)
+
+
+def ConvergenceVariableFactory(conv_var, **kwargs) -> ConvergenceVariable:
+
+ if conv_var == 'celldm1':
+ return conv_var_celldm1(**kwargs)
+ elif conv_var == 'ecutwfc':
+ return conv_var_ecutwfc(**kwargs)
+ elif conv_var == 'nbnd':
+ return conv_var_nbnd(**kwargs)
+ elif conv_var == 'kgrid':
+ return conv_var_kgrid(**kwargs)
+ else:
+ raise NotImplementedError(f'Convergence with respect to {conv_var} has not been directly implemented. You '
+ 'can still perform a convergence calculation with respect to this variable, but '
+ 'you must first create an appropriate ConvergenceVariable object and then '
+ 'construct your ConvergenceWorkflow using the ConvergenceWorkflowFactory')
+
+
+class ConvergenceWorkflow(Workflow):
+
+ '''
+ A Workflow class that wraps another workflow in a convergence procedure in order to converge the observable within the specified tolerance with respect to the variables
+ '''
+
+ def __init__(self, subworkflow_class: Type[Workflow], observable: Optional[Callable[[Workflow], float]] = None, threshold: Optional[float] = None, variables: List[ConvergenceVariable] = [], *args, **kwargs):
+ super().__init__(*args, **kwargs)
+
+ self._subworkflow_class = subworkflow_class
+ self.observable = observable
+ self.threshold = threshold
+ self.variables = variables
+
+ @classmethod
+ def fromdict(cls, dct: Dict[str, Any], **kwargs) -> Workflow:
+ kwargs.update(subworkflow_class=dct.pop('_subworkflow_class'),
+ observable=dct.pop('observable'),
+ threshold=dct.pop('threshold'),
+ variables=dct.pop('variables'))
+ return super(ConvergenceWorkflow, cls).fromdict(dct, **kwargs)
+
+ def _run(self) -> None:
+
+ # Deferred import to allow for monkeypatching
+ from koopmans import workflows
+
+ # Set the initial value for each of the convergence variables
+ for c in self.variables:
+ if c.initial_value is None:
+ c.get_initial_value(self)
+
+ # Check that everything has been initialized
+ if self.observable is None:
+ raise ValueError(
+ f'{self.__class__.__name__} has not been provided with an observable to converge')
+
+ if len(self.variables) == 0:
+ raise ValueError(
+ f'{self.__class__.__name__} has not been provided with any variables with which to perform convergence')
+
+ if self.threshold is None:
+ raise ValueError(
+ f'{self.__class__.__name__} has not been provided with a threshold with which to perform convergence')
if self.parameters.from_scratch:
- for c in convergence_parameters:
+ for c in self.variables:
for path in Path().glob(c.name + '*'):
shutil.rmtree(str(path))
# Create array for storing calculation results
- results = np.empty([initial_depth for _ in convergence_parameters])
+ results = np.empty([len(v) for v in self.variables])
results[:] = np.nan
- # Continue to increment the convergence parameters until convergence in the convergence observable
+ # Continue to increment the convergence variables until convergence in the convergence observable
# is achieved. A lot of the following code is quite obtuse, in order to make it work for an arbitrary
- # number of convergence parameters
+ # number of convergence variables
# Record the alpha values for the original calculation
- provide_alpha = any([p.name == 'nbnd' for p in convergence_parameters]) \
+ provide_alpha = any([p.name == 'nbnd' for p in self.variables]) \
and self.parameters.functional in ['ki', 'kipz', 'pkipz', 'all'] \
and self.parameters.alpha_from_file
if provide_alpha:
master_alphas = utils.read_alpha_file(directory=Path())
if self.parameters.orbital_groups is None:
- self.parameters.orbital_groups = list(range(self.calculator_parameters['kcp'].nbnd))
- master_orbital_groups = copy.deepcopy(self.parameters.orbital_groups)
+ self.parameters.orbital_groups = list(
+ range(self.calculator_parameters['kcp'].nbnd))
+ master_orbital_groups = copy.deepcopy(
+ self.parameters.orbital_groups)
while True:
- # Loop over all possible permutations of convergence parameter settings
- for indices in itertools.product(*[range(len(x.values)) for x in convergence_parameters]):
+ # Loop over all possible permutations of convergence variables
+ for indices in itertools.product(*[range(len(x)) for x in self.variables]):
# If we already have results for this permutation, don't recalculate it
if not np.isnan(results[tuple(indices)]):
continue
- # Create a new workflow
- subwf: workflows.Workflow
- if self.parameters.convergence_observable == 'eps_inf':
- subwf = workflows.DFTPhWorkflow.fromparent(self)
- else:
- subwf = workflows.SinglepointWorkflow.fromparent(self)
+ # Create a new subworkflow
+ subwf = self._subworkflow_class.fromparent(self)
# For each parameter we're converging wrt...
header = ''
subdir = Path()
- for index, parameter in zip(indices, convergence_parameters):
- value = parameter.values[index]
+ for index, variable in zip(indices, self.variables):
+ value = variable.values[index]
if isinstance(value, float):
value_str = f'{value:.1f}'
else:
value_str = str(value)
- header += f'{parameter.name} = {value_str}, '
+ header += f'{variable.name} = {value_str}, '
if isinstance(value, list):
value_str = ''.join([str(x) for x in value])
# Create new working directory
- subdir /= f'{parameter.name}_{value_str}'.replace(' ', '_').replace('.', 'd')
+ subdir /= f'{variable.name}_{value_str}'.replace(
+ ' ', '_').replace('.', 'd')
# Set the value
- parameter.set_value(subwf, value)
+ variable.set_value(subwf, value)
if provide_alpha:
# Update alpha files and orbitals
@@ -217,63 +335,104 @@ def set_value(wf: Workflow, value: List[int]) -> None:
subwf.run(subdirectory=subdir)
# Store the result
- results[indices] = fetch_result(subwf)
+ results[indices] = self.observable(subwf)
# Check convergence
- converged = np.empty([len(p.values) for p in convergence_parameters], dtype=bool)
+ converged = np.empty([len(p) for p in self.variables], dtype=bool)
converged[:] = False
- converged_result = results[tuple([-1 for _ in convergence_parameters])]
- threshold = self.parameters.convergence_threshold
+ converged_result = results[tuple([-1 for _ in self.variables])]
subarray_slice: List[Union[int, slice]]
- for indices in itertools.product(*[range(len(p.values)) for p in convergence_parameters]):
+ for indices in itertools.product(*[range(len(p)) for p in self.variables]):
subarray_slice = [slice(None) for _ in range(len(indices))]
for i_index, index in enumerate(indices):
subarray_slice[i_index] = slice(index, None)
- if np.all(np.abs(results[tuple(subarray_slice)] - converged_result) < threshold):
+ if np.all(np.abs(results[tuple(subarray_slice)] - converged_result) <