Deprecate top-level cutoff, always deduce it from bonds::rcut

.github/workflows/test.yml

@@ -42,6 +42,16 @@ jobs:
         if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
         #python setup.py install
         pip install . 
+    - name: Install maxvolpy
+      # `pip install .` builds a wheel and never runs setup.py's custom
+      # InstallMaxVolPyLocalPackage(install) hook, so lib/maxvolpy is not
+      # installed and the active-learning tests fail to import. Build it
+      # explicitly; --no-build-isolation is required because
+      # lib/maxvolpy/setup.py imports Cython and numpy at module level
+      # without declaring them in a pyproject.toml build-system.
+      run: |
+        pip install Cython scipy
+        pip install --no-build-isolation ./lib/maxvolpy
     - name: Test with python-ace with tensorpotential
       run: |
         python -m pytest tests/ --runtensorpot

bin/pacemaker

@@ -159,10 +159,9 @@ def main(args):
 
     assert isinstance(args_yaml, dict)
     if "cutoff" in args_yaml:
-        cutoff = args_yaml["cutoff"]
-    else:
-        log.warning("No 'cutoff' provided in YAML file, please specify it")
-        raise ValueError("No 'cutoff' provided in YAML file, please specify it")
+        log.warning("Top-level 'cutoff' is deprecated and is ignored. "
+                    "The cutoff for the neighbour list construction is deduced from the potential "
+                    "as the maximum over `potential::bonds::rcut`.")
 
     if "seed" in args_yaml:
         seed = args_yaml["seed"]
@@ -325,8 +324,8 @@ def generate_template_input():
         """Enter number of functions per element ([ENTER] - default 700): """) or 700)
     print("Number of functions per element: ", number_of_functions_per_element)
 
-    cutoff = float(input("Enter cutoff (Angstrom, default:7.0): ") or 7.0)
-    print("Cutoff: ", cutoff)
+    rcut = float(input("Enter outer cutoff `rcut` (Angstrom, default:7.0): ") or 7.0)
+    print("Outer cutoff rcut: ", rcut)
 
     # weighting scheme
     default_energy_based_weighting = """{ type: EnergyBasedWeightingPolicy, DElow: 1.0, DEup: 10.0, DFup: 50.0, DE: 1.0, DF: 1.0, wlow: 0.75, energy: convex_hull, reftype: all,seed: 42}"""
@@ -349,7 +348,7 @@ def generate_template_input():
         input_yaml_text = f.read()
 
     input_yaml_text = input_yaml_text.replace("{{ELEMENTS}}", str(elements))
-    input_yaml_text = input_yaml_text.replace("{{CUTOFF}}", str(cutoff))
+    input_yaml_text = input_yaml_text.replace("{{RCUT}}", str(rcut))
     input_yaml_text = input_yaml_text.replace("{{DATAFILENAME}}", train_filename)
     input_yaml_text = input_yaml_text.replace("{{number_of_functions_per_element}}",
                                               "number_of_functions_per_element: {}".format(

docs/pacemaker/faq.md

@@ -163,7 +163,8 @@ data:
 ## I want to change the cutoff, what should I do ?
 
 If you change cutoff, i.e. from `rcut: 7` to `rcut: 6.5`, then potential should be refitted from the scratch.
-`pacemaker` will recompute neighbourlists on every run, so, no need to extra options except for specifying cutoff. 
+Change `potential::bonds::rcut` only - the cutoff for the neighbour lists is deduced from it automatically.
+`pacemaker` recomputes the neighbour lists on every run, so no extra options are needed.
 
 
 

docs/pacemaker/inputfile.md

@@ -5,11 +5,13 @@ particular settings of `pacemaker`. The sections are listed below.
 
 ## Cutoff and (optional) metadata
 
-* Global cutoff for the neighborlist constructor is setup as:
+* The cutoff for the neighbour list construction is **not** specified separately. It is deduced from the potential
+  as the maximum of the bonds' outer cutoffs `potential::bonds::rcut` (see the `potential` section below),
+  which is exactly the cutoff that the ACE calculators (ASE, LAMMPS) use. The atomic environments seen during
+  fitting are therefore identical to those seen in production, and changing `rcut` is enough - nothing else
+  has to be kept in sync.
 
-```YAML
-cutoff: 10.0
-```
+  A top-level `cutoff` key is deprecated: it is ignored and only produces a warning.
 
 * Metadata (optional)
 

docs/pacemaker/quickstart.md

@@ -34,8 +34,8 @@ In order to interactively generate default `pacemaker` input file `input.yaml`,
 ```
 pacemaker -t
 ```
-and enter requested information, such as dataset filename, test set size (optional), list of elements, cutoff,
-number of functions.  Doing so will produce an `input.yaml` file with the most general
+and enter requested information, such as dataset filename, test set size (optional), list of elements,
+outer cutoff `rcut`, number of functions.  Doing so will produce an `input.yaml` file with the most general
 settings that can be adjusted for a particular task. Detailed overview of the input file parameters can be found in the [Input file](inputfile.md) section below.
 
 ## (optional) Manual fitting dataset preparation

examples/Cu-I/input.yaml

@@ -1,4 +1,3 @@
-cutoff: 7 # cutoff for neighbour list construction
 seed: 42  # random seed
 
 

examples/Cu-II/input.yaml

@@ -1,4 +1,3 @@
-cutoff: 7.5  # cutoff for neighbour list construction
 seed: 42  # random seed
 
 

examples/Ethanol/input.yaml

@@ -1,4 +1,3 @@
-cutoff: 4 # cutoff for neighbour list construction
 seed: 42  # random seed
 
 

examples/HEA/input.yaml

@@ -1,4 +1,3 @@
-cutoff: 7 # cutoff for neighbour list construction
 seed: 42  # random seed
 
 

src/pyace/data/input_template.yaml

@@ -1,4 +1,3 @@
-cutoff: {{CUTOFF}} # cutoff for neighbour list construction
 seed: 42  # random seed
 
 #################################################################
@@ -25,7 +24,7 @@ potential:
     ALL: {
       radbase: SBessel,
       radparameters: [ 5.25 ],
-      rcut: {{CUTOFF}},
+      rcut: {{RCUT}},
       dcut: 0.01,
       NameOfCutoffFunction: cos,
     }

src/pyace/generalfit.py

@@ -113,6 +113,19 @@ def get_initial_potential(start_potential):
     return initial_bbasisconfig
 
 
+def get_maximal_rcut(bbasisconfig: BBasisConfiguration) -> float:
+    """
+    Maximal outer cutoff `rcut` over all bond (i.e. one- and two-species) blocks.
+
+    This is the same quantity as `ACEBBasisSet.cutoffmax`, which the ACE calculators
+    (ASE, LAMMPS) use to construct their neighbour lists. Deducing the fitting cutoff
+    from it guarantees that the atomic environments seen during the fit are identical
+    to those seen in production.
+    """
+    return max(block.rcutij for block in bbasisconfig.funcspecs_blocks
+               if len(block.elements_vec) <= 2)
+
+
 def save_dataset(dataframe, fname):
     # columns to save: w_energy, w_forces, NUMBER_OF_ATOMS, PROTOTYPE_NAME, prop_id,structure_id, gen_id, if any
     # columns_to_save = ["PROTOTYPE_NAME", "NUMBER_OF_ATOMS", "prop_id", "structure_id", "gen_id", "pbc"] + \
@@ -143,6 +156,11 @@ class GeneralACEFit:
     :param data_config:  training data specification
     :param backend_config: specification of potential evaluation and fitting backend (pyace / tensorpot)
                     - dict ['evaluator']
+    :param cutoff: cutoff for the neighbour lists construction. If None (default), it is deduced from
+                    `potential_config` as the maximum over the bonds' outer cutoffs `rcut`, which is what
+                    the ACE calculators use in production. Provide it explicitly only if the neighbour
+                    lists have to extend beyond the potential's own cutoff; note that neighbours
+                    further than the corresponding bond's `rcut` contribute exactly zero to the fit.
     """
 
     def __init__(self,
@@ -265,11 +283,13 @@ def __init__(self,
             log.info("Ladder-scheme fitting is ON")
 
         if cutoff is None:
-            rcut = max(
-                [block.rcutij for block in self.target_bbasisconfig.funcspecs_blocks if len(block.elements_vec) <= 2])
-            self.cutoff = rcut
+            self.cutoff = get_maximal_rcut(self.target_bbasisconfig)
+            log.info("Cutoff for neighbour lists construction is deduced from the potential: "
+                     "{:.3f} A (maximum over bonds::rcut)".format(self.cutoff))
         else:
             self.cutoff = cutoff
+            log.info("Cutoff for neighbour lists construction is provided explicitly: "
+                     "{:.3f} A".format(self.cutoff))
 
         if self.ladder_scheme:
             if FIT_LADDER_TYPE_KW in fit_config:

tests/test-CLI/Cu-I/input.yaml

@@ -1,4 +1,3 @@
-cutoff: 7 # cutoff for neighbour list construction
 seed: 42  # random seed
 
 

tests/test-CLI/Cu-I/input_aug.yaml

@@ -1,4 +1,3 @@
-cutoff: 7 # cutoff for neighbour list construction
 seed: 42  # random seed
 
 

tests/test-CLI/Ethanol/input.yaml

@@ -1,4 +1,3 @@
-cutoff: 5.0 # cutoff for neighbour list construction
 seed: 42  # random seed
 
 #################################################################

tests/test-CLI/Ethanol/input_aug.yaml

@@ -1,4 +1,3 @@
-cutoff: 5.0 # cutoff for neighbour list construction
 seed: 42  # random seed
 
 #################################################################

tests/test_PyACECalculator.py

@@ -188,7 +188,11 @@ def test_fcc_stress():
 
 
 def test_relaxation():
-    from ase.constraints import UnitCellFilter
+    try:
+        # ASE >= 3.23
+        from ase.filters import UnitCellFilter
+    except ImportError:
+        from ase.constraints import UnitCellFilter
     from ase.optimize import QuasiNewton
 
     calc = PyACECalculator(basis_set="tests/Al.pbe.rhocore-v2.ace")

tests/test_cutoff_deduction.py

@@ -0,0 +1,67 @@
+import pytest
+
+from pyace import ACEBBasisSet
+from pyace.basisextension import construct_bbasisconfiguration
+from pyace.generalfit import get_maximal_rcut
+
+
+def build_potential_config(bonds):
+    return {
+        'deltaSplineBins': 0.001,
+        'elements': ['Al', 'Ni'],
+
+        'embeddings': {'ALL': {'drho_core_cut': 250,
+                               'fs_parameters': [1, 1],
+                               'ndensity': 1,
+                               'npot': 'FinnisSinclairShiftedScaled',
+                               'rho_core_cut': 200000}},
+
+        'bonds': bonds,
+
+        'functions': {
+            'UNARY': {
+                'nradmax_by_orders': [5, 2],
+                'lmax_by_orders': [0, 0]
+            },
+            'BINARY': {
+                'nradmax_by_orders': [5, 2, 2],
+                'lmax_by_orders': [0, 2, 2],
+            },
+        }
+    }
+
+
+DEFAULT_BOND = {'NameOfCutoffFunction': 'cos',
+                'core-repulsion': [10000.0, 5.0],
+                'dcut': 0.01,
+                'radbase': 'ChebPow',
+                'radparameters': [2.0],
+                'lmax': 2,
+                'rcut': 3.9}
+
+
+def test_get_maximal_rcut_uniform():
+    config = construct_bbasisconfiguration(build_potential_config({'ALL': DEFAULT_BOND}))
+    assert get_maximal_rcut(config) == pytest.approx(3.9)
+
+
+def test_get_maximal_rcut_per_bond():
+    bonds = {'ALL': DEFAULT_BOND,
+             'AlNi': dict(DEFAULT_BOND, rcut=5.5)}
+    config = construct_bbasisconfiguration(build_potential_config(bonds))
+    assert get_maximal_rcut(config) == pytest.approx(5.5)
+
+
+@pytest.mark.parametrize("binary_rcut", [3.9, 5.5])
+def test_get_maximal_rcut_equals_cutoffmax(binary_rcut):
+    """
+    The cutoff deduced for the neighbour lists must coincide with `ACEBBasisSet.cutoffmax`,
+    which is what the ACE calculators (ASE, LAMMPS) use in production. `get_maximal_rcut`
+    only inspects one- and two-species blocks, while `cutoffmax` is a maximum over all
+    blocks - these agree because the radial basis (`rcut` included) is copied from the
+    pair block into the higher-rank blocks.
+    """
+    bonds = {'ALL': DEFAULT_BOND,
+             'AlNi': dict(DEFAULT_BOND, rcut=binary_rcut)}
+    config = construct_bbasisconfiguration(build_potential_config(bonds))
+    assert get_maximal_rcut(config) == pytest.approx(ACEBBasisSet(config).cutoffmax)