Common format/API for hessians

[prs513rosewood]

The current API to compute/get a hessian in matscipy is the following function:

def get_hessian(self,
                atoms: ase.Atoms,
                format: str = 'sparse',
                divide_by_masses: bool = False) -> Union[scipy.sparse.spmatrix, np.ndarray]:
    ...

This has the following disatvantages:

1. interoperability with ASE's property system is limited since there are function parameters. I've worked around that in MatscipyCalculator by defining different properties for different parameter combinations (e.g. "dynamical_matrix" calls with divide_by_masses=True.
2. Not every calculator returns a sparse matrix (e.g. Ewald), which means conversions must sometimes happen
3. The sparse format should be bsr_matrix for all calculators: this makes it harder to change for specific calculators in case it's sub-optimal, and makes us less interoperable with other codes or future calculators that might use a different format, e.g. hierarchical matrices, or matrix-free formulations (other kspace implementations, fast-multipole, etc.)
4. There are two accepted formats: "sparse" and "neighbor-list", which are not interoperable

To solve these issues, I propose a new API:

def get_hessian(self,
                atoms: ase.Atoms) -> scipy.sparse.linalg.LinearOperator:
    ...

The abstract class LinearOperator is designed to abstract away details of how matrix-vector (or matrix-matrix) products are done and are fully composable, i.e. P = A * B and S = A + B can be defined for A and B that have different representation without needing conversion. They are also accepted as parameters to the sparse solve/diagonalization routines in scipy. Most calculators in matscipy would just need to call scipy.sparse.linalg.aslinearoperator() on the sparse matrix that they assemble, which returns a concrete instance of MatrixLinearOperator, from which the matrix representation can be retrieved, if need be (e.g. for tests). This addresses points 2 & 3 above.

For point 1, the divide by masses operation should be done in a separate function. For point 4, after discussing with @griessej , we figured that the neighbor-list format is only really used in the pair and polydisperse calculators for the calculation of the Born constants. It then makes sense to move whatever code in MatscipyCalculator that uses that format to these calculators (and actually have the polydisperse calculator inherit from the pair calculator).

LMK what you think of this change. The main burden of making this change would probably be to change the tests of hessian values.

[prs513rosewood]

One more thing to consider for sparse hessians: Scipy is switching to an array interface for sparse matrices to be compatible with Numpy (see note here).

@pastewka @jameskermode any opinions on the suggested changes?

[pastewka]

I think we should switch to the array format. The distinction between an array and a matrix is just super confusing. I tend to avoid matrix in all my codes, and I think we should do the same with the sparse formats.

[jameskermode]

Agreed

[pastewka]

@prs513rosewood and myself just discussed this:

• Returning a linear operator generically seems more portable. The method should then reflect this and be called get_hessian_operator or similar
• The internal block coordinate structure of the sparse Hessian (that is based on the neighbor list) should somehow be kept, as this appears the most convenient way to actually compute the Hessian. The idea is so have a method get_hessian_block_format (or similar) that returns it. The get_hessian could then be in the superclass and be generic among calculators.

[pastewka]

-> get_block_sparse_hessian