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

# TODO: Use another base class to hide internals from user
[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=np.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(, 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