# Source code for hiphive.cutoffs

```import pickle
import numpy as np
from ase.neighborlist import NeighborList
from hiphive.input_output.pretty_table_prints import table_array_to_string

[docs]class Cutoffs:
""" This class maintains information about the cutoff configuration,
i.e. which clusters will be included (="inside cutoff"). It also
encapsulates functionality that is used e.g., during cluster space
construction.

Here, `n-body` refers to number of atoms in a cluster. For example
the cluster (0011) is a two-body cluster of fourth order and the
cluster (123) is a three-body cluster of third order.

Parameters
----------
cutoff_matrix : numpy.ndarray
the matrix element `ij` provides to the cutoff for order `j+2`
and nbody `i+2`; elements with `i>j` will be ignored
"""

def __init__(self, cutoff_matrix):

self._cutoff_matrix = np.array(cutoff_matrix, dtype=float)

if len(self._cutoff_matrix.shape) != 2:
raise ValueError('Please specify cutoff matrix as a 2D array')
for i, row in enumerate(self._cutoff_matrix):
if np.any(row[i:] < 0):
raise ValueError('Negative number as cutoff')
row[:i] = np.NaN
self._cutoff_matrix = self._cutoff_matrix[:(self.max_nbody-1), :(self.max_order-1)]

@property
def cutoff_matrix(self):
""" numpy.ndarray : copy of cutoff matrix """
return self._cutoff_matrix.copy()

@property
def orders(self):
""" list(int) : allowed orders """
return list(range(2, self.max_order + 1))

@property
def nbodies(self):
""" list(int) : allowed bodies """
return list(range(2, self.max_nbody + 1))

[docs]    def get_cutoff(self, order, nbody):
"""
Returns cutoff for a given body and order.

Parameters
----------
order : int
nbody : int

Raises
------
ValueError
if order is not in orders
ValueError
if nbody is not in nbodies
ValueError
if nbody is larger than order

Returns
-------
float
"""
if order not in self.orders:
raise ValueError('order not in orders')
if nbody not in self.nbodies:
raise ValueError('nbody not in nbodies')
if nbody > order:
raise ValueError('nbody can not be larger than order')
return self._cutoff_matrix[nbody - 2, order - 2]

@property
def max_cutoff(self):
""" float : maximum cutoff """
max_cutoff = 0
for i, row in enumerate(self._cutoff_matrix):
max_cutoff = max(max_cutoff, np.max(row[i:]))
return max_cutoff

@property
def max_nbody(self):
""" int : maximum body """
nbody = 1
for i, row in enumerate(self._cutoff_matrix):
if np.any(row[i:]):
nbody = i + 2
return nbody

@property
def max_order(self):
""" int : maximum order """
order = None
for col in range(self._cutoff_matrix.shape[1]):
if np.any(self._cutoff_matrix[:col + 1, col]):
order = col + 2
return order

[docs]    def max_nbody_cutoff(self, nbody):
""" Return maximum cutoff for a given body. """
if nbody not in self.nbodies:
raise ValueError('nbody not in nbodies')
return np.max(self._cutoff_matrix[nbody - 2, max(0, nbody - 2):])

[docs]    def max_nbody_order(self, nbody):
""" Returns maximum order for a given body """
if nbody not in self.nbodies:
raise ValueError('nbody not in nbodies')
row = self._cutoff_matrix[nbody - 2]
max_order = None
for order, cutoff in enumerate(row[nbody-2:], start=nbody):
if cutoff:
max_order = order
return max_order

[docs]    def write(self, fileobj):
""" Writes instance to file.

Parameters
----------
fileobj : file-like object
file-like object to which the cutoffs will be written to
"""
pickle.dump(self._cutoff_matrix, fileobj)

[docs]    def read(fileobj):
""" Reads an instance from file.

Parameters
----------
fileobj : file-like object
input file to read from
"""
data = pickle.load(fileobj)
return Cutoffs(data)

[docs]    def to_filename_tag(self):
""" Simple function turning cutoffs into a string to be used in e.g.
filenames. """
s = []
for i, c in enumerate(self._cutoff_matrix.tolist(), start=2):
s.append('{}body-{}'.format(i, '_'.join(map(str, c))))
return '_'.join(s)

def __str__(self):
cutoff_matrix = self._cutoff_matrix.copy()
cutoff_matrix = np.vstack(([[None] * len(self.orders)], cutoff_matrix))
s = table_array_to_string(cutoff_matrix)

width = max(len(c) for c in s.split('\n'))
header = ' Cutoffs '.center(width, '=') + '\n'
bottom = '\n' + ''.center(width, '=')
s = header + s + bottom
return s

def __repr__(self):
return 'Cutoffs({!r})'.format(self._cutoff_matrix)

[docs]class CutoffMaximumBody(Cutoffs):
""" Specify cutoff-list plus maximum body

Usefull when creating e.g. 6th order expansions but with only 3-body
interactions.

Parameters
----------
cutoff_list : list
list of cutoffs for order 2, 3, etc. Must be in decresing order
max_nbody : int
No clusters containing more than max_nbody atoms will be generated
"""

def __init__(self, cutoff_list, max_nbody):
cutoff_matrix = np.zeros((max_nbody - 1, len(cutoff_list)))
for order, cutoff in enumerate(cutoff_list, start=2):
cutoff_matrix[:, order - 2] = cutoff
super().__init__(cutoff_matrix)

[docs]def is_cutoff_allowed(atoms, cutoff):
""" Checks if atoms is compatible with cutoff

Parameters
----------
atoms : ase.Atoms
structure used for checking compatibility with cutoff
cutoff : float
cutoff to be tested

Returns
-------
bool
True if cutoff compatible with atoms object, else False
"""
nbrlist = NeighborList(cutoffs=[cutoff / 2] * len(atoms), skin=0,
self_interaction=False, bothways=True)
nbrlist.update(atoms)

for i in range(len(atoms)):
neighbors, _ = nbrlist.get_neighbors(i)
if i in neighbors:
return False
if len(neighbors) != len(set(neighbors)):
return False
return True

[docs]def estimate_maximum_cutoff(atoms, max_iter=11):
""" Estimates the maximum possible cutoff given the atoms object

Parameters
----------
atoms : ase.Atoms
structure used for checking compatibility with cutoff
max_iter : int
number of iterations in binary search
"""

# First upper boundary of cutoff
upper_cutoff = min(np.linalg.norm(atoms.cell, axis=1))

# generate all possible offsets given upper_cutoff
nbrlist = NeighborList(cutoffs=[upper_cutoff / 2] * len(atoms), skin=0,
self_interaction=False, bothways=True)
nbrlist.update(atoms)
all_offsets = []
for i in range(len(atoms)):
_, offsets = nbrlist.get_neighbors(i)
all_offsets.extend([tuple(offset) for offset in offsets])

# find lower boundary and new upper boundary
unique_offsets = set(all_offsets)
unique_offsets.discard((0, 0, 0))
upper = min(np.linalg.norm(np.dot(offset, atoms.cell))
for offset in unique_offsets)
lower = upper / 2.0

# run binary search between the upper and lower bounds
for _ in range(max_iter):
cutoff = (upper + lower) / 2
if is_cutoff_allowed(atoms, cutoff):
lower = cutoff
else:
upper = cutoff
return lower

[docs]class BaseClusterFilter:
"""Base cluster filter class.

This filter simply accepts all proposed clusters.  A proper
subclass must implement the same methods.
"""
[docs]    def setup(self, atoms):
""" The filter is passed the environment of the primitive cell.

Parameters
----------
atoms : ase.Atoms
non-pbc primitive cell plus neighboring atoms
"""
self._atoms = atoms

def __call__(self, cluster):
""" Returns True or False when a cluster is proposed.

Parameters
----------
cluster : tuple(int)
indices of proposed cluster referenced to the internal
atoms object
"""
return True
```