How to add a new representation
Write a Calculator for a representation
A calculator of a representation is an object that builds a representation of the atomic structure contained in a structure manager. This class:
inherits publicly from
CalculatorBaseto follow its interface and use some of the common utilies shared by such class.uses a
compute()function expecting one or a list of structure manager(s) to build the representation efficiently and attach the resulting features to their respective structure manager.
The behaviour of a representation calculator is defined at construction by hyperparameters which are specific to the representation. These hyperparameters can:
define a particular implementation from a set of conceptually equivalent methods, implemented using primitive enums
be scalar values or on/off features used as hyperparameters in the calculation of the represenation, implemented using primitive types
To illustrate the basic structure that a new representation that would be implemented in calculator_representation_name.hh should follow, let’s take the example of the CalculatorSortedCoulomb. A detailed discussion of the sorted coulomb matrix representation can be found in Ref. 1 and Ref. 2.
The representation starts with the definition of some useful short hands
/**
* Implementation of the Environmental Coulomb Matrix
*/
class CalculatorSortedCoulomb : public CalculatorBase {
public:
using Parent = CalculatorBase;
// type of the hyperparameters
using Hypers_t = typename Parent::Hypers_t;
// numeric type for the representation features
using Precision_t = typename Parent::Precision_t;
// type of the data structure for the representation feaures
template <class StructureManager>
using Property_t =
Property<Precision_t, 1, StructureManager, Eigen::Dynamic, 1>;
// short hand type to help the iteration over the structure manager
template <class StructureManager, size_t Order>
using ClusterRef_t = typename StructureManager::template ClusterRef<Order>;
// type of the datastructure used to register the list of valid
// hyperparameters
using ReferenceHypers_t = Parent::ReferenceHypers_t;
followed by the definition of its constructors and destructor
//! Constructor
explicit CalculatorSortedCoulomb(const Hypers_t & hyper)
: CalculatorBase{} {
this->set_default_prefix("sorted_coulomb_");
this->check_hyperparameters(this->reference_hypers, hyper);
// Extract the options and hyperparameters
this->set_hyperparameters(hyper);
}
//! Copy constructor
CalculatorSortedCoulomb(const CalculatorSortedCoulomb & other) = delete;
//! Move constructor
CalculatorSortedCoulomb(CalculatorSortedCoulomb && other) noexcept
: CalculatorBase{std::move(other)}, central_cutoff{std::move(
other.central_cutoff)},
central_decay{std::move(other.central_decay)},
interaction_cutoff{std::move(other.interaction_cutoff)},
interaction_decay{std::move(other.interaction_decay)},
size{std::move(other.size)} {}
//! Destructor
virtual ~CalculatorSortedCoulomb() = default;
//! Copy assignment operator
CalculatorSortedCoulomb &
operator=(const CalculatorSortedCoulomb & other) = delete;
//! Move assignment operator
CalculatorSortedCoulomb &
operator=(CalculatorSortedCoulomb && other) = default;
the declaration of the concrete implementation of the calculator interface
/**
* Compute representation for a given structure manager.
*
* @tparam StructureManager a (single or collection)
* of structure manager(s) (in an iterator) held in shared_ptr
*/
template <class StructureManager>
void compute(StructureManager & managers);
//! set hypers
void set_hyperparameters(const Hypers_t & /*hypers*/) override;
void update_central_cutoff(double /*cutoff*/);
/**
* check if size of the calculator is enough for current structure
* manager.
* size refers to the parameter that regulate the feature size of the
* calculator.
*/
template <class StructureManager>
void check_size_compatibility(StructureManager & manager) {
for (auto center : manager) {
auto n_neighbours{center.pairs().size()};
if (n_neighbours > this->size) {
std::cout << "size is too small for this "
"structure and has been reset to: "
<< n_neighbours << std::endl;
this->size = n_neighbours;
}
}
}
and the declaration of some functions for internal use in the protected section. The end of the class contains the different internal variables needed by the class
// list of hyperparameters specific to the coulomb matrix
// spherical cutoff for the atomic environment
double central_cutoff{};
double central_decay{};
double interaction_cutoff{};
double interaction_decay{};
// at least equal to the largest number of neighours
size_t size{};
//! reference the requiered hypers
ReferenceHypers_t reference_hypers{
{"central_cutoff", {}},
{"central_decay", {}},
{"interaction_cutoff", {}},
{"interaction_decay", {}},
{"size", {}},
{"sorting_algorithm", {"distance", "row_norm"}},
};
Implementing several types of the same representation
Often there are different ways to implement the same fundamental representation. To be able to do this in rascal, the compute function is used with the type of the implementation as template argument (or multiple if representation requires this). The switch between several implementations of conceptually equivalent parts of a representation can be implemented through several mechanisms such a virtual inheritance. We detail here how to implement such switch efficiently using the CalculatorSortedCoulomb as an example.
To make the coulomb matrix invariant with respect to permutation Ref. 2 proposes to sort the upper triangular part of the coulomb matrix according to the norm of each rows or the distance from the central atom (see 1 for details).
The implementation of these two behaviour is encapsulated in the SortCoulomMatrix class and the choice between them is done with a template parameter using template specialization. Note that in this particular case templated functions could be sufficient but to underline how to implement the most general case a class is used.
// Enum class defining the several possible sorting options of the Coulomb
// Matrix
enum class CMSortAlgorithm {
Distance, // sort according to the distance from the central atom
RowNorm, // sort according to the norm of the coulomb matrix rows
};
// Empty general template class implementing the determination of the
// sorting order for the coulomb matrix. It should never be instantiated.
template <CMSortAlgorithm Method>
struct SortCoulomMatrix {};
The specific implementation of the two options is done in with template specialization
template <>
struct SortCoulomMatrix<CMSortAlgorithm::Distance> {
/**
* Sort the coulomb matrix using the distance to the central atom
* as reference order.
*
* @param distance_mat distance matrix between all the atoms in the
* neighbourhood
*/
static std::vector<std::pair<size_t, distiter>>
get_coulomb_matrix_sorting_order(
const Eigen::Ref<const Eigen::MatrixXd> & distance_mat,
const Eigen::Ref<const Eigen::MatrixXd> &) {
// initialize the distances to be sorted. the center is always first
std::vector<double> distances_to_sort{0};
distances_to_sort.reserve(distance_mat.cols());
for (auto idx_i{1}; idx_i < distance_mat.cols(); ++idx_i) {
distances_to_sort.push_back(distance_mat(idx_i, 0));
}
// find the sorting order
std::vector<std::pair<size_t, distiter>> order_coulomb(
distances_to_sort.size());
size_t nn{0};
for (distiter it{distances_to_sort.begin()};
it != distances_to_sort.end(); ++it, ++nn) {
order_coulomb[nn] = make_pair(nn, it);
}
// use stable sort
std::stable_sort(order_coulomb.begin(), order_coulomb.end(),
ordering::ascending);
return order_coulomb;
}
};
template <>
struct SortCoulomMatrix<CMSortAlgorithm::RowNorm> {
/**
* Sort the coulomb matrix using the distance to the central atom
* as reference order.
*
* @param coulomb_mat coulomb matris between all the atoms in the
* neighbourhood
*/
static std::vector<std::pair<size_t, distiter>>
get_coulomb_matrix_sorting_order(
const Eigen::Ref<const Eigen::MatrixXd> &,
const Eigen::Ref<const Eigen::MatrixXd> & coulomb_mat) {
// initialize the distances to be sorted. the center is always first
std::vector<double> distances_to_sort{};
distances_to_sort.reserve(coulomb_mat.cols());
auto row_norms = coulomb_mat.colwise().squaredNorm().eval();
row_norms(0) = 1e200;
for (auto idx_i{0}; idx_i < coulomb_mat.cols(); ++idx_i) {
distances_to_sort.push_back(row_norms(idx_i));
}
std::vector<std::pair<size_t, distiter>> order_coulomb(
distances_to_sort.size());
size_t nn{0};
for (distiter it{distances_to_sort.begin()};
it != distances_to_sort.end(); ++it, ++nn) {
order_coulomb[nn] = make_pair(nn, it);
}
// use stable sort
std::stable_sort(order_coulomb.begin(), order_coulomb.end(),
ordering::descending);
return order_coulomb;
}
};
The switch between the two behaviours is done in the compute() function which chooses the right type of computation method to compute the representation for each structure manager. At the moment the compute_loop() function has to be copied to every new calculator to allow iterations over structure managers. It is not necessary to understand this code, if one is only interested to integegrate a new representation into rascal.
//! loop over a collection of manangers (note that maps would raise a
//! compilation error)
template <
internal::CMSortAlgorithm AlgorithmType, class StructureManager,
std::enable_if_t<internal::is_proper_iterator<StructureManager>::value,
int> = 0>
void compute_loop(StructureManager & managers) {
for (auto & manager : managers) {
this->compute_impl<AlgorithmType>(manager);
}
}
//! if it is not a list of managers
template <internal::CMSortAlgorithm AlgorithmType, class StructureManager,
std::enable_if_t<
not(internal::is_proper_iterator<StructureManager>::value),
int> = 0>
void compute_loop(StructureManager & manager) {
this->compute_impl<AlgorithmType>(manager);
}
The actual implementation of the representation is in the function compute_impl() which is templated with the computation method.
template <class StructureManager>
void CalculatorSortedCoulomb::compute(StructureManager & managers) {
auto option{this->options["sorting_algorithm"]};
if (option == "distance") {
compute_loop<internal::CMSortAlgorithm::Distance>(managers);
} else if (option == "row_norm") {
compute_loop<internal::CMSortAlgorithm::RowNorm>(managers);
} else {
auto error_message{std::string("Option '") + option +
std::string("' is not implemented.")};
throw std::invalid_argument(error_message.c_str());
}
}
Write the python bindings of a new representation
We use pybind11 to handle the generation of the python bindings. To make a new representation available to the python users, you need to explicitly register your representation calculator. Since the function compute() is templated with the type of the input structure manager, it has to be binded explicitly for every possible structure manager stack type that you want to allow as input. To do this you need to add you own binding code in the add_representation_calculators() function in bindings/bind_py_representation_calculator.cc using the available infrastructure. Here is an example on how it is done for the sorted coulomb representation.
// Defines a particular structure manager type
using TypeHolder_t =
StructureManagerTypeHolder<StructureManagerCenters,
AdaptorNeighbourList, AdaptorStrict>;
using ManagerList_t = typename TypeHolder_t::type_list;
// using Manager_t = typename TypeHolder_t::type;
// StructureManagerCenters,AdaptorNeighbourList, AdaptorStrict
// Defines the representation manager type for the particular structure
// manager
using Calc1_t = CalculatorSortedCoulomb;
// Bind the interface of this representation manager
auto rep_sorted_coulomb =
add_representation_calculator<Calc1_t>(mod, m_internal);
bind_compute_function_helper<ManagerList_t>(rep_sorted_coulomb);
Note
To be able to use a particular structure manager stack in python, it also has to be binded. In the case a valid structure manager stack for your representation is not already binded, you will have to register it too in the add_structure_managers() function in bindings/bind_py_structure_manager.cc like in this example:
// bind the root structure manager
bind_structure_manager<StructureManagerCenters>(m_strc_mng, m_internal);
bind_make_structure_manager<StructureManagerCenters>(m_nl);
// bind a structure manager stack
BindAdaptorStack<StructureManagerCenters, AdaptorNeighbourList,
AdaptorStrict>
adaptor_stack_1{m_nl, m_adp, m_internal, name_list};
// bind the corresponding manager collection
bind_structure_manager_collection<StructureManagerCenters,
AdaptorNeighbourList, AdaptorStrict>(
m_nl);
The last step is to write a python class in bindings/rascal/representations/ to simplify the use of the representation from the python side. You can use the implentation of the SortedCoulombMatrix in coulomb_matrix.py as a template.
Implement and test gradients
In principle, all libRascal representations should implement gradients with
respect to the atomic positions. Currently the only representations to do so
are the SphericalExpansion and SphericalInvariants
(tensor order 0, body order 0–1 a.k.a. “RadialSpectrum” and “PowerSpectrum”) in
the GTO and DVR radial basis. Until we come up with a general, standard way of
implementing gradients for any representation, please see those implementations
for guidance.
Alternatively, the utility compute_numerical_kernel_gradients() can
be used to compute the gradient of a the kernel for a particular representation
w.r.t. the atomic positions using finite differences allowing to build kernel
models with atomic forces when analytical gradients of the representation
are missing.
Once you’ve implemented the gradient – or derivative – of any function in libRascal, you must test that it actually corresponds to the gradient of the function that it purpots to be. A finite-difference testing function is provided for this purpose; see Testing gradients for details.
- 1(1,2)
- 2(1,2)
Rupp, M., Tkatchenko, A., Müller, K.-R., & von Lilienfeld, O. A. (2011). Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning. Physical Review Letters, 108(5), 58301. https://doi.org/10.1103/PhysRevLett.108.058301