From 60e294c0da0378c2d99295b50848e9e55c7c7b0d Mon Sep 17 00:00:00 2001 From: Lucas Ortengren <31080445+ortengren@users.noreply.github.com> Date: Thu, 6 Jun 2024 15:00:23 -0500 Subject: [PATCH] Merge documentation branch (#34) * Configured Sphinx and reformatted docstrings * Adding subtract self to the neighborlist generator (#19) * Switching list comp (#24) Switching the code to do the list comp once instead of within each loop * Adding ability to override number of radial bases * manual override n and added test cases * removed extraneous arg * Linter * Linter * Pointing to cutoff radius properly * I messed with the history so now I'm restoring it * pass linter * removed error check restricting independent specification of max_radial and radial_gaussian_width * added tests to ensure constant n functionality is correct * linter * num_radial -> max_radial in edp.py * removed changes pertaining to coupling sigma and n to keep PR clean * restored original compute_gaussian_parameters function header * incorporated necessary changes to make constant n work for both orthonormalization and moments generation * pass linter * removed deprecated comment * made definition of maxradial consistent with num_n and updated tests * Update lint.yml (#25) add diff to linter so we can see how files are failing the linter. * Resolved merge conflict by incorporating both suggestions * formatted docstrings, adjusted Sphinx config * Fixed docstrings, adjusted sphinx config, manually rewrote modules.rst * Fixed docstring formatting and changed sphinx theme * Added custom CSS, replaced modules.rst with api.rst, adjusted sphinx conf and index.rst * added ReadTheDocs configuration * Fixed requirements.txt for RTD building * Fixed requirements.txt * Removed anisoap from requirements.txt * removed unneccessary packages from requirements.txt * changed config to not fail upon warning * edited RTD config and requirements.txt * changed requirements.txt to install metatensor and rascaline * Fixing requirements.txt * Fixed requirements.txt * fixed requirements.txt * fixed requirements.txt * Changed RTD config * fixed rust version in RTD config * fixed rust version in RTD config * changed pytorch requirement to cpu version * fixed torch dependency * edited metatensor dependency * added tqdm to requirements.txt * deleted modules.rst * fixed api.rst * fixed autosummary in api.rst * fixed autodocs in api.rst * Changed api.rst to have more detailed description * further formatted docstrings * fixed mathjax formatting * fixed mathjax formatting * restructuring documentation * Fixing LaTeX and docstring formatting * Changing copyright info, author info * Removing .idea folder * Building docs * Built HTML and fixed docstrings * linter * isort * linter * deleted anisoap/reference/projection_coefficients.py * Removed documentation for deleted module * Fixed RST syntax errors in docstrings * Formatted according to Black code style * Delete anisoap/.idea directory These seem ide specific, don't want them being on the anisoap main branch. * equistore --> metatensor * minor correction --------- Co-authored-by: Rose K. Cersonsky <47536110+rosecers@users.noreply.github.com> Co-authored-by: Arthur Lin Co-authored-by: arthur-lin1027 <35580059+arthur-lin1027@users.noreply.github.com> --- .gitignore | 3 - .readthedocs.yaml | 18 + anisoap/reference/projection_coefficients.py | 199 --------- .../ellipsoidal_density_projection.py | 156 ++++--- anisoap/representations/radial_basis.py | 379 ++++++++++++------ anisoap/utils/metatensor_utils.py | 66 ++- anisoap/utils/moment_generator.py | 112 ++++-- anisoap/utils/monomial_iterator.py | 40 +- anisoap/utils/spherical_to_cartesian.py | 66 ++- docs/Makefile | 20 + docs/make.bat | 35 ++ docs/requirements.txt | 16 + docs/source/_static/custom.css | 3 + docs/source/api.rst | 37 ++ docs/source/conf.py | 48 +++ docs/source/getting_started.rst | 17 + docs/source/index.rst | 28 ++ docs/source/installation.rst | 44 ++ docs/source/tutorials.rst | 8 + docs/source/usage.rst | 4 + 20 files changed, 819 insertions(+), 480 deletions(-) create mode 100644 .readthedocs.yaml delete mode 100644 anisoap/reference/projection_coefficients.py create mode 100644 docs/Makefile create mode 100644 docs/make.bat create mode 100644 docs/requirements.txt create mode 100644 docs/source/_static/custom.css create mode 100644 docs/source/api.rst create mode 100644 docs/source/conf.py create mode 100644 docs/source/getting_started.rst create mode 100644 docs/source/index.rst create mode 100644 docs/source/installation.rst create mode 100644 docs/source/tutorials.rst create mode 100644 docs/source/usage.rst diff --git a/.gitignore b/.gitignore index 6d39677..fc2d359 100644 --- a/.gitignore +++ b/.gitignore @@ -71,9 +71,6 @@ instance/ # Scrapy stuff: .scrapy -# Sphinx documentation -docs/_build/ - # PyBuilder target/ diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 0000000..383f647 --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,18 @@ +version: "2" + +build: + os: "ubuntu-22.04" + tools: + python: "3.10" + rust: "1.70" + +python: + install: + - requirements: docs/requirements.txt + # Install our python package before building the docs + - method: pip + path: . + +sphinx: + configuration: docs/source/conf.py + fail_on_warning: true diff --git a/anisoap/reference/projection_coefficients.py b/anisoap/reference/projection_coefficients.py deleted file mode 100644 index 0ec09d7..0000000 --- a/anisoap/reference/projection_coefficients.py +++ /dev/null @@ -1,199 +0,0 @@ -import numpy as np - -from anisoap.utils.spherical_to_cartesian import spherical_to_cartesian - -try: - from tqdm import tqdm -except ImportError: - tqdm = lambda i, **kwargs: i - -from ..utils import ( - compute_moments_diagonal_inefficient_implementation, - compute_moments_inefficient_implementation, - quaternion_to_rotation_matrix, -) -from .radial_basis import RadialBasis - - -class DensityProjectionCalculator: - """ - Compute the spherical projection coefficients. - Initialize the calculator using the hyperparameters. - ---------- - max_angular : int - Number of angular functions - radial_basis : str - The radial basis. Currently implemented are - 'GTO_primitive', 'GTO', 'monomial'. - compute_gradients : bool - Compute gradients - subtract_center_contribution : bool - Subtract contribution from the central atom. - Attributes - ---------- - features : numpy.ndarray - feature_gradients : numpy.ndarray - """ - - def __init__( - self, - max_angular, - radial_basis_name, - cutoff_radius, - compute_gradients=False, - subtract_center_contribution=False, - radial_gaussian_width=None, - ): - # Store the input variables - self.max_angular = max_angular - self.cutoff_radius = cutoff_radius - self.compute_gradients = compute_gradients - self.subtract_center_contribution = subtract_center_contribution - self.radial_basis_name = radial_basis_name - - # Currently, gradients are not supported - if compute_gradients: - raise NotImplementedError("Sorry! Gradients have not yet been implemented") - - # Precompute the spherical to Cartesian transformation - # coefficients. - num_ns = [] - for l in range(max_angular + 1): - num_ns.append(max_angular + 1 - l) - self.sph_to_cart = spherical_to_cartesian(max_angular, num_ns) - - # Initialize the radial basis class - if radial_basis_name not in ["monomial", "gto"]: - raise ValueError( - f"{self.radial_basis} is not an implemented basis" - ". Try 'monomial' or 'gto'" - ) - if radial_gaussian_width != None and radial_basis_name != "gto": - raise ValueError("Gaussian width can only be provided with GTO basis") - radial_hypers = {} - radial_hypers["radial_basis"] = radial_basis_name.lower() # lower case - radial_hypers["radial_gaussian_width"] = radial_gaussian_width - radial_hypers["max_angular"] = max_angular - self.radial_basis = RadialBasis(**radial_hypers) - self.num_ns = self.radial_basis.get_num_radial_functions() - - def transform(self, frames, show_progress=False): - """ - Computes the features and (if compute_gradients == True) gradients - for all the provided frames. The features and gradients are stored in - features and feature_gradients attribute. - Parameters - ---------- - frames : ase.Atoms - List containing all ase.Atoms structures - show_progress : bool - Show progress bar for frame analysis - Returns - ------- - None, but stores the projection coefficients and (if desired) - gradients as arrays as `features` and `features_gradients`. - """ - self.frames = frames - - # Generate a dictionary to map atomic types to array indices - # In general, the types are sorted according to atomic number - # and assigned the array indices 0, 1, 2,... - # Example: for H2O: H is mapped to 0 and O is mapped to 1. - types = set() - for frame in frames: - for atom in frame: - types.add(atom.number) - types = sorted(types) - self.types_dict = {} - for frame in frames: - # Get atomic types in dataset - self.types_dict.update( - {atom.symbol: types.index(atom.number) for atom in frame} - ) - - # Define variables determining size of feature vector coming from frames - self.num_atoms_per_frame = np.array([len(frame) for frame in frames]) - num_atoms_total = np.sum(self.num_atoms_per_frame) - num_particle_types = len(self.types_dict) - num_features_total = (self.max_angular + 1) ** 2 - - # Initialize arrays in which to store all features - self.features = np.zeros( - num_atoms_total, num_particle_types, num_features_total - ) - self.feature_gradients = 0 - - if show_progress: - frame_generator = tqdm(self.frames) - else: - frame_generator = self.frames - - for i_frame, frame in enumerate(frame_generator): - number_of_atoms = self.num_atoms_per_frame[i_frame] - results = self._transform_single_frame(frame) - - def _transform_single_frame(self, frame): - """ - Compute features for single frame and return to the transform() - method which loops over all structures to obtain the complete - vector for all environments. - """ - ### - # Initialization - ### - # Define useful shortcuts - lmax = self.max_angular - num_atoms = len(frame) - num_chem_types = len(self.types_dict) - iterator_types = np.zeros(num_atoms, dtype=int) - for i, symbol in enumerate(frame.get_chemical_symbols()): - iterator_types[i] = self.types_dict[symbol] - - # Get the arrays with all - # TODO: update with correct expressions - positions = np.zeros((num_atoms, 3)) - quaternions = np.zeros((num_atoms, 4)) - ellipsoid_lengths = np.zeros((num_atoms, 3)) - - # Convert quaternions to rotation matrices - rotation_matrices = np.zeros((num_atoms, 3, 3)) - for i, quat in enumerate(quaternions): - rotation_matrices[i] = quaternion_to_rotation_matrix(quat) - - # Generate neighbor list - # TODO: change this to proper neighbor list - neighbors = [] - for i in range(num_atoms): - # for now, just treat every atom as a neighbor - # of itself + the first two atoms in the structure - neighbors.append([0, 1, i]) - - # Compute the features for a single frame - features = [] - for l in range(lmax + 1): - features.append(np.zeros((3, 3))) - - for i in range(num_atoms): - pos_i = positions[i] - for j in neighbors[i]: - # Obtain the position and orientation defining - # the neighbor particle j - r_ij = pos_i - positions[j] - rot = rotation_matrices[j] - lengths = ellipsoid_lengths[j] - - # Compute the moments - # The moments have shape ((maxdeg+1, maxdeg+1, maxdeg+1)) - precision, center = self.radial_basis.compute_gaussian_parameters( - r_ij, lengths, rot - ) - moments = compute_moments_inefficient_implementation( - precision, center, maxdeg=lmax - ) - - for l in range(lmax + 1): - features[l] = np.einsum( - "mnpqr, pqr->mn", self.sph_to_cart[l], moments - ) - - return features diff --git a/anisoap/representations/ellipsoidal_density_projection.py b/anisoap/representations/ellipsoidal_density_projection.py index 7968e81..8f46093 100644 --- a/anisoap/representations/ellipsoidal_density_projection.py +++ b/anisoap/representations/ellipsoidal_density_projection.py @@ -31,48 +31,47 @@ def pairwise_ellip_expansion( radial_basis, show_progress=False, ): - """ - Function to compute the pairwise expansion by combining the moments and the spherical to Cartesian - transformation - -------------------------------------------------------- - Parameters: - lmax : int - Maximum angular + r"""Computes pairwise expansion - neighbor_list : Equistore TensorMap - Full neighborlist with keys (types_1, types_2) enumerating the possible types pairs. - Each block contains as samples, the atom indices of (first_atom, second_atom) that correspond to the key, - and block.value is a 3D array of the form (num_samples, num_components,properties), with num_components=3 - corresponding to the (x,y,z) components of the vector from first_atom to second_atom. - Depending on the cutoff some types pairs may not appear. Self pairs are not present but in PBC, - pairs between copies of the same atom are accounted for. + Function to compute the pairwise expansion :math:`\langle anlm|\rho_{ij} \rangle` + by combining the moments and the spherical to Cartesian transformation. - types: list of ints + Parameters + ---------- + lmax : int + Maximum angular + neighbor_list : metatensor.TensorMap + Full neighborlist with keys (types_1, types_2) enumerating the possible + species pairs. Each block contains as samples, the atom indices of + (first_atom, second_atom) that correspond to the key, and block.value is + a 3D array of the form (num_samples, num_components,properties), with + num_components=3 corresponding to the (x,y,z) components of the vector + from first_atom to second_atom. Depending on the cutoff some species + pairs may not appear. Self pairs are not present but in PBC, pairs between + copies of the same atom are accounted for. + types : list of ints List of atomic numbers present across the data frames - - frame_to_global_atom_idx: list of ints - The length of the list equals the number of frames, each entry enumerating the number of atoms in - corresponding frame - - rotation_matrices: np.array of dimension ((num_atoms,3,3)) - - ellipsoid_lengths: np.array of dimension ((num_atoms,3)) - + frame_to_global_atom_idx : list of ints + The length of the list equals the number of frames, each entry enumerating + the number of atoms in corresponding frame + rotation_matrices : np.array of dimension ((num_atoms,3,3)) + ellipsoid_lengths : np.array of dimension ((num_atoms,3)) radial_basis : Instance of the RadialBasis Class - anisoap.representations.radial_basis.RadialBasis that has been instantiated appropriately with the - cutoff radius, radial basis type. - + anisoap.representations.radial_basis.RadialBasis that has been instantiated + appropriately with the cutoff radius, radial basis type. show_progress : bool Show progress bar for frame analysis and feature generation - ----------------------------------------------------------- - Returns: - An Equistore TensorMap with keys (types_1, types_2, l) where ("types_1", "types_2") is key in the - neighbor_list and "l" is the angular channel. - Each block of this tensormap has the same samples as the corresponding block of the neighbor_list. - block.value is a 3D array of the form (num_samples, num_components, properties) where num_components - form the 2*l+1 values for the corresponding angular channel and the properties dimension corresponds to - the radial channel. + Returns + ------- + TensorMap + A metatensor TensorMap with keys (types_1, types_2, l) where + ("types_1", "types_2") is key in the neighbor_list and "l" is the + angular channel. Each block of this tensormap has the same samples as the + corresponding block of the neighbor_list. block.value is a 3D array of + the form (num_samples, num_components, properties) where num_components + form the 2*l+1 values for the corresponding angular channel and the + properties dimension corresponds to the radial channel. """ tensorblock_list = [] keys = np.asarray(neighbor_list.keys, dtype=int) @@ -190,31 +189,35 @@ def pairwise_ellip_expansion( def contract_pairwise_feat(pair_ellip_feat, types, show_progress=False): - """ - Function to sum over the pairwise expansion \sum_{j in a} = - -------------------------------------------------------- - Parameters: + """Function to sum over the pairwise expansion + + .. math:: - pair_ellip_feat : Equistore TensorMap - TensorMap returned from "pairwise_ellip_expansion()" with keys (types_1, types_2,l) enumerating - the possible types pairs and the angular channels. + \\sum_{j \\in a} \\langle anlm|\\rho_{ij} \\rangle + = \\langle anlm|\\rho_i \\rangle + Parameters + ---------- + pair_ellip_feat : metatensor.TensorMap + TensorMap returned from "pairwise_ellip_expansion()" with keys + (types_1, types_2,l) enumerating the possible species pairs and the + angular channels. types: list of ints List of atomic numbers present across the data frames - show_progress : bool Show progress bar for frame analysis and feature generation - ----------------------------------------------------------- - Returns: - An Equistore TensorMap with keys (types, l) "types" takes the value of the atomic numbers present - in the dataset and "l" is the angular channel. - Each block of this tensormap has as samples ("type", "center") yielding the indices of the frames - and atoms that correspond to "types" are present. - block.value is a 3D array of the form (num_samples, num_components, properties) where num_components - take on the same values as in the pair_ellip_feat_feat.block . block.properties now has an additional index - for neighbor_types that corresponds to "a" in - + Returns + ------- + TensorMap + A metatensor TensorMap with keys (types, l) "types" takes the value + of the atomic numbers present in the dataset and "l" is the angular + channel. Each block of this tensormap has as samples ("type", "center"), + yielding the indices of the frames and atoms that correspond to "species" + are present. block.value is a 3D array of the form (num_samples, num_components, properties) + where num_components take on the same values as in the pair_ellip_feat_feat.block. + block.properties now has an additional index for neighbor_species that + corresponds to "a" in :math:`\\langle anlm|rho_i \\rangle` """ ellip_keys = list( set([tuple(list(x)[:1] + list(x)[2:]) for x in pair_ellip_feat.keys]) @@ -402,14 +405,18 @@ def contract_pairwise_feat(pair_ellip_feat, types, show_progress=False): class EllipsoidalDensityProjection: - """ + """Class for computing spherical projection coefficients. + Compute the spherical projection coefficients for a system of ellipsoids assuming a multivariate Gaussian density. + Initialize the calculator using the hyperparameters. + + Parameters ---------- max_angular : int Number of angular functions - radial_basis : str + radial_basis : _RadialBasis The radial basis. Currently implemented are 'gto', 'monomial'. compute_gradients : bool @@ -431,6 +438,7 @@ class EllipsoidalDensityProjection: ---------- features : numpy.ndarray feature_gradients : numpy.ndarray + """ def __init__( @@ -447,6 +455,35 @@ def __init__( basis_rcond=0, basis_tol=1e-8, ): + """Instantiates an object of type EllipsoidalDensityProjection. + + Parameters + ---------- + max_angular : int + Number of angular functions + radial_basis_name : str + The radial basis. Currently implemented are 'GTO_primitive', 'GTO', + and 'monomial'. + cutoff_radius + Cutoff radius of the projection + compute_gradients : bool, optional + Compute gradients; defaults to 'False' + subtract_center_contribution : bool, optional + Subtract contribution from the central atom. Defaults to 'False' + radial_gaussian_width : float, optional + Width of the Gaussian + max_radial : int, list of int + Number of radial bases to use. Can either correspond to number of + bases per spherical harmonic or a value to use with every harmonic. + If `None`, then for every `l`, `(max_angular - l) // 2 + 1` will + be used. + rotation_key : string + Key under which rotations are stored in ase frames arrays + rotation_type : string + Type of rotation object being passed. Currently implemented + are 'quaternion' and 'matrix' + + """ # Store the input variables self.max_angular = max_angular self.cutoff_radius = cutoff_radius @@ -509,10 +546,12 @@ def __init__( self.rotation_key = rotation_key def transform(self, frames, show_progress=False, normalize=True): - """ + """Computes features and gradients for frames + Computes the features and (if compute_gradients == True) gradients for all the provided frames. The features and gradients are stored in features and feature_gradients attribute. + Parameters ---------- frames : ase.Atoms @@ -520,12 +559,13 @@ def transform(self, frames, show_progress=False, normalize=True): show_progress : bool Show progress bar for frame analysis and feature generation normalize: bool - Whether to perform Lowdin Symmetric Orthonormalization or not. Orthonormalization generally - leads to better performance. Default: True. + Whether to perform Lowdin Symmetric Orthonormalization or not. + Returns ------- None, but stores the projection coefficients and (if desired) gradients as arrays as `features` and `features_gradients`. + """ self.frames = frames diff --git a/anisoap/representations/radial_basis.py b/anisoap/representations/radial_basis.py index 7b2f952..d2b1b27 100644 --- a/anisoap/representations/radial_basis.py +++ b/anisoap/representations/radial_basis.py @@ -6,23 +6,25 @@ def inverse_matrix_sqrt(matrix: np.array, rcond=1e-8, tol=1e-3): - """ - Returns the inverse matrix square root. - The inverse square root of the overlap matrix (or slices of the overlap matrix) yields the - orthonormalization matrix - Args: - matrix: np.array - Symmetric square matrix to find the inverse square root of - - rcond: float - Lower bound for eigenvalues for inverse square root - - tol: float - Tolerance for differences between original matrix and reconstruction via - inverse square root - - Returns: - inverse_sqrt_matrix: S^{-1/2} + r"""Returns the inverse matrix square root. + + The inverse square root of the overlap matrix (or slices of the overlap matrix) + yields the orthonormalization matrix + + Parameters + ---------- + matrix : np.array + Symmetric square matrix to find the inverse square root of + rcond: float + Lower bound for eigenvalues for inverse square root + tol: float + Tolerance for differences between original matrix and reconstruction via + inverse square root + + Returns + ------- + inverse_sqrt_matrix + :math:`S^{-1/2}` """ if not np.allclose(matrix, matrix.conjugate().T): @@ -44,53 +46,87 @@ def inverse_matrix_sqrt(matrix: np.array, rcond=1e-8, tol=1e-3): def gto_square_norm(n, sigma): - """ - Compute the square norm of GTOs (inner product of itself over R^3). - An unnormalized GTO of order n is \phi_n = r^n * e^{-r^2/(2*\sigma^2)} + r"""Compute the square norm of GTOs (inner product of itself over :math:`R^3`). + + An unnormalized GTO of order n is :math:`\phi_n = r^n e^{-r^2/(2*\sigma^2)}` The square norm of the unnormalized GTO has an analytic solution: - <\phi_n | \phi_n> = \int_0^\infty dr r^2 |\phi_n|^2 = 1/2 * \sigma^{2n+3} * \Gamma(n+3/2) - Args: - n: order of the GTO - sigma: width of the GTO - Returns: - square norm: The square norm of the unnormalized GTO + .. math:: + + \braket{\phi_n | \phi_n} &= \int_0^\infty dr \, r^2 \lvert\phi_n\rvert^2 \\ + &= \frac{1}{2} \sigma^{2n+3} \Gamma(n+\frac{3}{2}) + + This function uses the above expression. + + Parameters + ---------- + n + order of the GTO + sigma + width of the GTO + + Returns + ------- + float + The square norm of the unnormalized GTO + """ return 0.5 * sigma ** (2 * n + 3) * gamma(n + 1.5) def gto_prefactor(n, sigma): - """ - Computes the normalization prefactor of an unnormalized GTO. - This prefactor is simply 1/sqrt(square_norm_area). - Scaling a GTO by this prefactor will ensure that the GTO has square norm equal to 1. - Args: - n: order of GTO - sigma: width of GTO + """Computes the normalization prefactor of an unnormalized GTO. - Returns: - N: normalization constant + This prefactor is simply :math:`\\frac{1}{\\sqrt(\\text{square_norm_area)}}`. + Scaling a GTO by this prefactor will ensure that the GTO has square norm + equal to 1. + + Parameters + ---------- + n + order of GTO + sigma + width of GTO + + Returns + ------- + float + The normalization constant """ return np.sqrt(1 / gto_square_norm(n, sigma)) def gto_overlap(n, m, sigma_n, sigma_m): - """ - Compute overlap of two *normalized* GTOs - Note that the overlap of two GTOs can be modeled as the square norm of one GTO, with an effective - n and sigma. All we need to do is to calculate those effective parameters, then compute the normalization. - <\phi_n, \phi_m> = \int_0^\infty dr r^2 r^n * e^{-r^2/(2*\sigma_n^2) * r^m * e^{-r^2/(2*\sigma_m^2) - = \int_0^\infty dr r^2 |r^{(n+m)/2} * e^{-r^2/4 * (1/\sigma_n^2 + 1/\sigma_m^2)}|^2 - = \int_0^\infty dr r^2 |r^{n_{eff}} * e^{-r^2/(2*\sigma_{eff}^2)|^2 - ---Arguments--- - n: order of the first GTO - m: order of the second GTO - sigma_n: sigma parameter of the first GTO - sigma_m: sigma parameter of the second GTO - - ---Returns--- - S: overlap of the two normalized GTOs + r"""Compute overlap of two *normalized* GTOs + + Note that the overlap of two GTOs can be modeled as the square norm of one + GTO, with an effective :math:`n` and :math:`\sigma`. All we need to do is to + calculate those effective parameters, then compute the normalization prefactor. + + .. math:: + + \langle \phi_n, \phi_m \rangle &= \int_0^\infty dr \, r^2 r^n e^{-r^2/(2\sigma_n^2)} \, r^m e^{-r^2/(2\sigma_m^2)} \\ + &= \int_0^\infty dr \, r^2 \lvert r^{(n+m)/2} e^{-r^2/4 (1/\sigma_n^2 + 1/\sigma_m^2)}\rvert^2 \\ + &= \int_0^\infty dr \, r^2 r^n_\text{eff} e^{-r^2/(2\sigma_\text{eff}^2)} + + + Parameters + ---------- + n + order of the first GTO + m + order of the second GTO + sigma_n + sigma parameter of the first GTO + sigma_m + sigma parameter of the second GTO + + Returns + ------- + float + overlap of the two normalized GTOs + """ N_n = gto_prefactor(n, sigma_n) N_m = gto_prefactor(m, sigma_m) @@ -103,11 +139,15 @@ def monomial_square_norm(n, r_cut): """ Compute the square norm of monomials (inner product of itself over R^3). - Args: - n: order of the basis + Parameters + ---------- + n + order of the basis - Returns: - square norm: The square norm of the unnormalized basis + Returns + ------- + float + The square norm of the unnormalized basis """ return r_cut ** (2 * n + 3) / (2 * n + 3) @@ -115,33 +155,50 @@ def monomial_square_norm(n, r_cut): def monomial_prefactor(n, r_cut): """ Computes the normalization prefactor of an unnormalized monomial basis. - This prefactor is simply 1/sqrt(square_norm_area). + This prefactor is simply :math:`1/sqrt{square_norm_area}`. Scaling a basis by this prefactor will ensure that the basis has square norm equal to 1. - Args: - n: order of basis - Returns: - N: normalization constant + Parameters + ---------- + n + order of the basis + + Returns + ------- + float + normalization constant """ return np.sqrt(1 / monomial_square_norm(n, r_cut)) def monomial_overlap(n, m, r_cut): - """ + r""" Compute overlap of two *normalized* monomials - Note that the overlap of two monomials can be modeled as the square norm of one monomial, with an effective - n. All we need to do is to calculate those effective parameters, then compute the normalization. - <\phi_n, \phi_m> = \int_0^\infty dr r^2 r^n r^m - = \int_0^r_{cut} dr r^2 |r^{(n+m)/2}|^2 - = \int_0^r_{cut} dr r^2 |r^{n_{eff}}|^2 - - ---Arguments--- - n: order of the first monomial - m: order of the second monomial - r_cut: float. cutoff radius - ---Returns--- - S: overlap of the two normalized monomial + + Note that the overlap of two monomials can be modeled as the square norm of + one monomial, with an effective n. All we need to do is to calculate those + effective parameters, then compute the normalization: + + .. math:: + + \langle \phi_n, \phi_m \rangle &= \int_0^\infty dr r^2 r^n r^m \\ + &= \int_0^r_{cut} dr r^2 |r^{(n+m)/2}|^2 \\ + &= \int_0^r_{cut} dr r^2 |r^{n_{eff}}|^2 + + Parameters + ---------- + n + order of the first monomial + m + order of the second monomial + r_cut: float + cutoff radius + Returns + ------- + float + overlap of the two normalized monomial + """ N_n = monomial_prefactor(n, r_cut) N_m = monomial_prefactor(m, r_cut) @@ -155,8 +212,9 @@ class _RadialBasis: This helps to keep a cleaner main code by avoiding if-else clauses related to the radial basis. - Code relating to GTO orthonormalization is heavily inspired by work done in librascal, specifically this - codebase here: https://github.com/lab-cosmo/librascal/blob/8405cbdc0b5c72a5f0b0c93593100dde348bb95f/bindings/rascal/utils/radial_basis.py + Code relating to GTO orthonormalization is heavily inspired by work done in + librascal, specifically the codebase found + `here `_ """ @@ -218,20 +276,25 @@ def get_num_radial_functions(self): Otherwise, the outputted list will specify the number of radial basis functions per l, which may be automatically calculated if max_radial=None. If a custom list of max_radial is specified when initializing, then it will return the same inputted list. - Returns: - self.num_radial_functions: list containing the number of radial basis functions considered per l. + + Returns + ------- + array_like + The attribute `self.num_radial_functions`, which is a list containing + the number of radial basis functions considered per `l`. """ return self.num_radial_functions def plot_basis(self, n_r=100): """ - Plot the normalized basis functions used in calculating the expansion coefficients - Args: - n_r: int - number of mesh points. Default: 100 + Plot the normalized basis functions used in calculating the expansion + coefficients + + Parameters + ---------- + n_r: int + number of mesh points. Default: 100 - Returns: - None """ from matplotlib import pyplot as plt @@ -240,12 +303,22 @@ def plot_basis(self, n_r=100): class MonomialBasis(_RadialBasis): - """ - A subclass of _RadialBasis that contains attributes and methods required for the Monomial basis. - The monomial basis of order n is defined to be R(r) = r^n. - For monomial basis with defined nmax and lmax, our radial basis set consists of monomials of order n=0 to n=lmax + 2nmax. - For monomial basis with coupled nmax and lmax, our radial basis set consists of monomials of order n=0 to n=max(lmax + 2nmax) - Monomials are not square-integrable from [0, ∞], so we orthonormalize by integrating up to the cutoff radius + r""" + A subclass of _RadialBasis that contains attributes and methods required for + the Monomial basis. + + The monomial basis of order n is defined to be :math:`R(r) = r^n`. + + For monomial basis with defined :math:`n_{\text{max}}` and :math:`l_{\text{max}}`, + our radial basis set consists of monomials of order :math:`n=0` to + :math:`n=l_{\text{max}} + 2n_{\text{max}}`. + + For monomial basis with coupled :math:`n_{\text{max}}` and :math:`l_{\text{max}}`, + our radial basis set consists of monomials of order :math:`n=0` to + :math:`n=\text{max}{l_{\text{max}} + 2n_{\text{max}}}` + + Monomials are not square-integrable from :math:`[0, \infty]`, so we orthonormalize + by integrating up to the cutoff radius """ def __init__( @@ -285,16 +358,24 @@ def compute_gaussian_parameters(self, r_ij, lengths, rotation_matrix): def calc_overlap_matrix(self): """ Computes the overlap matrix for Monomnials over a fixed interval. - The overlap matrix is a Gram matrix whose entries are the overlap: S_{ij} = \int_0^r_cut dr r^2 phi_i phi_j + + The overlap matrix is a Gram matrix whose entries are the overlap: + + .. math:: + S_{ij} = \int_{0}^{r_{cut} dr r^2 phi_i phi_j + The overlap has an analytic solution (see above functions). + The overlap matrix is the first step to generating an orthonormal basis set of functions (Lodwin Symmetric Orthonormalization). The actual orthonormalization matrix cannot be fully precomputed because each tensor block use a different set of bases. Hence, we precompute the full overlap matrix of dim l_max, and while orthonormalizing each tensor block, we generate the respective orthonormal matrices from slices of the full overlap matrix. - Returns: - S: 2D array. The overlap matrix + Returns + ------- + 2D array: + The overlap matrix """ max_deg = np.max( np.arange(self.max_angular + 1) + 2 * np.array(self.num_radial_functions) @@ -313,13 +394,18 @@ def orthonormalize_basis(self, features: TensorMap): An instructive example of Lodwin Symmetric Orthonormalization of a 2-element basis set is found here: https://booksite.elsevier.com/9780444594365/downloads/16755_10030.pdf - Parameters: - features: A TensorMap whose blocks' values we wish to orthonormalize. Note that features is modified in place, so a - copy of features must be made before the function if you wish to retain the unnormalized values. - radial_basis: An instance of _RadialBasis - - Returns: - normalized_features: features containing values multiplied by proper normalization factors. + Parameters + ---------- + features: TensorMap + A TensorMap whose blocks' values we wish to orthonormalize. Note that + `features` is modified in place, so a copy of `features` must be made + before the function if you wish to retain the unnormalized values. + radial_basis: _RadialBasis + + Returns + ------- + TensorMap + features containing values multiplied by proper normalization factors. """ # In-place modification. @@ -353,14 +439,20 @@ def orthonormalize_basis(self, features: TensorMap): def get_basis(self, rs): """ Evaluate orthonormalized monomial basis functions on mesh rs. + If lmax and nmax defined, then the number of functions outputted is lmax*(nmax+1) + If lmax and nmax coupled, then the number of functions outputted is \sum_{l=0}^{lmax} (number_of_radial_functions_per_l) - Args: + + Parameters + ---------- rs: np.array a mesh to evaluate the basis functions - Returns: - S: a matrix containing orthonormalized monomial basis functions evaluated on rs + Returns + ------- + np.array + a matrix containing orthonormalized monomial basis functions evaluated on rs """ all_gs = np.empty(shape=(len(rs), 1)) for l in range(0, self.max_angular): @@ -392,8 +484,11 @@ def get_basis(self, rs): class GTORadialBasis(_RadialBasis): """ A subclass of _RadialBasis that contains attributes and methods required for the GTO basis. - The GTO basis of order n is defined to be R(r) = r^n * e^(-r^2/(2*sigma^2)). + + The GTO basis of order n is defined to be :math:`R(r) = r^n e^{(-r^2/(2\sigma^2))}`. + For GTO basis with defined nmax and lmax, our radial basis set consists of GTO of order n=0 to n=lmax + 2nmax. + For GTO basis with coupled nmax and lmax, our radial basis set consists of GTO of order n=0 to n=max(lmax + 2nmax) """ @@ -443,18 +538,29 @@ def compute_gaussian_parameters(self, r_ij, lengths, rotation_matrix): return new_precision, new_center, constant def calc_overlap_matrix(self): - """ - Computes the overlap matrix for GTOs. - The overlap matrix is a Gram matrix whose entries are the overlap: S_{ij} = \int_0^\infty dr r^2 phi_i phi_j + """Computes the overlap matrix for GTOs. + + The overlap matrix is a Gram matrix whose entries are the overlap: + + .. math:: + + S_{ij} = \\int_0^\\infty dr \\, r^2 \\phi_i \\phi_j + The overlap has an analytic solution (see above functions). - The overlap matrix is the first step to generating an orthonormal basis set of functions (Lodwin Symmetric - Orthonormalization). The actual orthonormalization matrix cannot be fully precomputed because each tensor - block use a different set of GTOs. Hence, we precompute the full overlap matrix of dim l_max, and while - orthonormalizing each tensor block, we generate the respective orthonormal matrices from slices of the full + + The overlap matrix is the first step to generating an orthonormal basis + set of functions (Lodwin Symmetric Orthonormalization). The actual + orthonormalization matrix cannot be fully precomputed because each tensor + block uses a different set of GTOs. Hence, we precompute the full overlap + matrix of dim l_max, and while orthonormalizing each tensor block, we + generate the respective orthonormal matrices from slices of the full overlap matrix. - Returns: - S: 2D array. The overlap matrix + Returns + ------- + 2D array + The overlap matrix + """ max_deg = np.max( np.arange(self.max_angular + 1) + 2 * np.array(self.num_radial_functions) @@ -471,20 +577,28 @@ def calc_overlap_matrix(self): return S def orthonormalize_basis(self, features: TensorMap): - """ - Apply an in-place orthonormalization on the features, using Lodwin Symmetric Orthonormalization. - Each block in the features TensorMap uses a GTO set of l + 2n, so we must take the appropriate slices of - the overlap matrix to compute the orthonormalization matrix. - An instructive example of Lodwin Symmetric Orthonormalization of a 2-element basis set is found here: + """Applies in-place orthonormalization on the features. + + Apply an in-place orthonormalization on the features, using Lodwin Symmetric + Orthonormalization. Each block in the features TensorMap uses a GTO set + of l + 2n, so we must take the appropriate slices of the overlap matrix + to compute the orthonormalization matrix. An instructive example of Lodwin + Symmetric Orthonormalization of a 2-element basis set is found here: https://booksite.elsevier.com/9780444594365/downloads/16755_10030.pdf - Parameters: - features: A TensorMap whose blocks' values we wish to orthonormalize. Note that features is modified in place, so a - copy of features must be made before the function if you wish to retain the unnormalized values. - radial_basis: An instance of _RadialBasis + Parameters + ---------- + features : TensorMap + contains blocks whose values we wish to orthonormalize. Note that + features is modified in place, so a copy of features must be made + before the function if you wish to retain the unnormalized values. + radial_basis : RadialBasis + + Returns + ------- + TensorMap + features containing values multiplied by normalization factors. - Returns: - normalized_features: features containing values multiplied by proper normalization factors. """ # In-place modification. radial_basis_name = self.radial_basis @@ -517,16 +631,23 @@ def orthonormalize_basis(self, features: TensorMap): return features def get_basis(self, rs): - """ + r""" Evaluate orthonormalized GTO basis functions on mesh rs. + If lmax and nmax defined, then the number of functions outputted is lmax*(nmax+1) - If lmax and nmax coupled, then the number of functions outputted is \sum_{l=0}^{lmax} (number_of_radial_functions_per_l) - Args: - rs: np.array - a mesh to evaluate the basis functions - Returns: - S: a matrix containing orthonormalized GTO basis functions evaluated on rs + If lmax and nmax coupled, then the number of functions outputted is + :math:`\sum_{l=0}^{\text{lmax}} (\text{number_of_radial_functions_per_l})` + + Parameters + ---------- + rs: np.array + a mesh to evaluate the basis functions + + Returns + ------- + np.array + a matrix containing orthonormalized GTO basis functions evaluated on rs """ from matplotlib import pyplot as plt diff --git a/anisoap/utils/metatensor_utils.py b/anisoap/utils/metatensor_utils.py index ec20468..e87ebb5 100644 --- a/anisoap/utils/metatensor_utils.py +++ b/anisoap/utils/metatensor_utils.py @@ -116,11 +116,31 @@ def combine_einsum(self, rho1, rho2, L, combination_string): def _real2complex(L): - """ + r"""Computes a matrix that converts from real to complex coefficients. + Computes a matrix that can be used to convert from real to complex-valued - spherical harmonics(coefficients) of order L. + spherical harmonics (coefficients) of order `L`. + + Parameters + ---------- + L : int + The order of the spherical harmonics for which the matrix will be computed. + + Returns + ------- + np.ndarray + Matrix that can be used to convert from real to complex-valued spherical + harmonics of order `L`. + + Note + ---- + The matrix generated is meant to be applied from the left; that is, if :math:`\mathbf{R}` + is the matrix returned from the function, it should be used like this: + + .. math:: + + \mathbf{R} \dot \left[-L..L\right] - It's meant to be applied to the left, ``real2complex @ [-L..L]``. """ result = np.zeros((2 * L + 1, 2 * L + 1), dtype=np.complex128) @@ -207,15 +227,37 @@ def cg_combine( lcut=None, other_keys_match=None, ): - """ - Performs a CG product of two sets of equivariants. Only requirement is that - sparse indices are labeled as ("inversion_sigma", "spherical_harmonics_l", "order_nu"). The automatically-determined - naming of output features can be overridden by giving a list of "feature_names". - By defaults, all other key labels are combined in an "outer product" mode, i.e. if there is a key-side - neighbor_types in both x_a and x_b, the returned keys will have two neighbor_types labels, - corresponding to the parent features. By providing a list `other_keys_match` of keys that should match, these are - not outer-producted, but combined together. for instance, passing `["types center"]` means that the keys with the - same types center will be combined together, but yield a single key with the same types_center in the results. + """Performs a CG product of two sets of equivariants. + + The only requirement is that sparse indices are labeled as + ("inversion_sigma", "spherical_harmonics_l", "order_nu"). + + Parameters + ---------- + x_a + First set of equivariants + x_b + Second set of equivariants + feature_names : list, optional + Overrides automatically-generated names of output features. By default, + all other key labels are combined via their outer product, i.e. if there + is a key-side `neighbor-species` in both `x_a` and `x_b`, the returned + keys will have two `neighbor_species` labels, corresponding to the parent + features. + clebsch_gordan : ClebschGordanReal, optional + lcut : np.ndarray() + Cutoff in new features + other_keys_match : list, optional + List of keys that should match. These will not need to have their outer + product taken, but will instead be merged into a new key. For instance, + passing `["types center"]` will combine the keys with the same type + center, yielding a single key with the same types_center in the results. + + Returns + ------- + TensorMap + The Clebsch-Gordan product of `x_a` and `x_b` + """ # determines the cutoff in the new features diff --git a/anisoap/utils/moment_generator.py b/anisoap/utils/moment_generator.py index 95a29bd..252e00b 100644 --- a/anisoap/utils/moment_generator.py +++ b/anisoap/utils/moment_generator.py @@ -9,21 +9,32 @@ # Define function to compute all moments for a single # variable Gaussian. def compute_moments_single_variable(A, a, maxdeg): - """ - Parameters: - - A: float + r"""Computes all moments for a single variable Gaussian. + + Parameters + ---------- + A : float inverse of variance (see below for exact mathematical form) - - a: float + a : float Center of Gaussian function - - maxdeg: int + maxdeg : int Maximum degree for which the moments need to be computed. - Returns: - - A numpy array of size (maxdeg+1,) containing the moments defined as - = integral x^n exp(-A(x-a)^2/2) dx - Note that the Gaussian is not normalized, meaning that we - need to multiply all results by a global factor if we wish - to interpret these as moments of a probability density. + Returns + ------- + np.array + A numpy array of size (`maxdeg` + 1,) containing the moments defined as + + .. math:: + + \langle x^n \rangle = \int x^n e^{-A(x-a)^2/2} dx + + Note + ---- + The Gaussian is not normalized, meaning that we need to multiply + all results by a global factor if we wish to interpret these as moments of + a probability density. + """ assert maxdeg > 0 @@ -49,26 +60,43 @@ def compute_moments_single_variable(A, a, maxdeg): def compute_moments_diagonal_inefficient_implementation( principal_components, a, maxdeg ): - """ - Parameters: - - principal_components: np.ndarray of shape (3,) + r"""Computes all moments for a diagonal dilation matrix. + + The implementation focuses on conceptual simplicity, while sacrificing memory + efficiency. To be specific, the `moments` array allows access to the value + of the moment :math:`\langle x^{n_0} y^{n_1} z^{n_2} \rangle` simply as + `moments[n0, n1, n2]`. This leads to more intuitive code, at the cost of + wasting around a third of the memory in the array to store zeros. + + Parameters + ---------- + principal_components : np.ndarray of shape (3,) Array containing the three principal components - - a: np.ndarray of shape (3,) + a : np.ndarray of shape (3,) Vectorial center of the trivariate Gaussian - - maxdeg: int + maxdeg : int Maximum degree for which the moments need to be computed. - Returns: - - moments: np.ndarray of shape (3, maxdeg+1) - moments[n0,n1,n2] is the (n0,n1,n2)-th moment of the Gaussian defined as + Returns + ------- + np.ndarray of shape (3, `maxdeg` + 1) + Moments calculated. `moments[n0,n1,n2]` is the :math:`\left(n_0,n_1,n_2\right)^\text{th}` + moment of the Gaussian defined as .. math:: - = \int(x^{n_0} * y^{n_1} * z^{n_2}) * \exp(-0.5*(r-a).T@cov@(r-a)) dxdydz - \sum_{i=1}^{\\infty} x_{i} - Note that the term "moments" in probability theory are defined for normalized Gaussian distributions. - Here, we take the Gaussian without prefactor, meaning that all moments are scaled - by a global factor. + \langle x^{n_0} y^{n_1} z^{n_2} \rangle = + \int(x^{n_0} y^{n_1} z^{n_2}) e^{-\frac{1}{2}(r-a)^T \Sigma (r-a)} dx\,dy\,dz\, + \sum_{i=1}^{\infty} x_{i} + + where :math:`\Sigma` is the covariance matrix. + + Note + ---- + The term "moments" in probability theory are defined for normalized Gaussian + distributions. Here, we take the Gaussian without prefactor, meaning that + all moments are scaled by a global factor. + """ # Initialize the array in which to store the moments # moments[n0, n1, n2] will be set to @@ -110,25 +138,35 @@ def compute_moments_diagonal_inefficient_implementation( # This leads to more intuitive code, at the cost of wasting around # a third of the memory in the array to store zeros. def compute_moments_inefficient_implementation(A, a, maxdeg): - """ - Parameters: - - A: symmetric 3x3 matrix (np.ndarray of shape (3,3)) - Dilation matrix of the Gaussian that determines its shape. - It can be written as cov = RDR^T, where R is a rotation matrix that specifies - the orientation of the three principal axes, while D is a diagonal matrix - whose three diagonal elements are the lengths of the principal axes. - - a: np.ndarray of shape (3,) + r"""Computes all moments for a general dilation matrix. + + Parameters + ---------- + A : np.ndarray of shape (3,3) + Dilation matrix of the Gaussian that determines its shape. It can be + written as :math:`\mathbf{A} = \mathbf{R}\mathbf{D}\mathbf{R}^T`, where + R is a rotation matrix that specifies the orientation of the three principal + axes, while :math:`\mathbf{D}` is a diagonal matrix whose three diagonal + elements are the lengths of the principal axes. + a : np.ndarray of shape (3,) Vectorial center of the trivariate Gaussian. - - maxdeg: int + maxdeg : int Maximum degree for which the moments need to be computed. - Returns: - - The list of moments defined as + Returns + ------- + np.ndarray + The list of moments defined as .. math:: - = \int(x^{n_0} * y^{n_1} * z^{n_2}) * \exp(-0.5*(r-a).T@cov@(r-a)) dxdydz - \sum_{i=1}^{\\infty} x_{i} + \langle x^{n_0} y^{n_1} z^{n_2} \rangle = + \int(x^{n_0} y^{n_1} z^{n_2}) \exp(-0.5(r-a)^T \Sigma (r-a)) dx\,dy\,dz\, + \sum_{i=1}^{\infty} x_{i} + + where :math:`\Sigma` is the covariance matrix. + Note + ---- Note that the term "moments" in probability theory are defined for normalized Gaussian distributions. These recursive relations find the moments of a normalized Gaussian, but we actually want to find the moments of the unnormalized underlying gaussian (as seen in the equation above) because calculating expansion diff --git a/anisoap/utils/monomial_iterator.py b/anisoap/utils/monomial_iterator.py index fd797dc..e8b027f 100644 --- a/anisoap/utils/monomial_iterator.py +++ b/anisoap/utils/monomial_iterator.py @@ -2,27 +2,28 @@ class TrivariateMonomialIndices: - """ - Class for generating an iterator object over all trivariate - monomials of the form f(x,y,z) = x^n0 * y^n1 * z^n2 - sorted in the lexicographical order. + """Class for generating an iterator object over trivariate monomials. + + Generates an iterator object over trivariate monomials of the form + :math:`f(x,y,z) = x^{n_0} y^{n_1} z^{n_2}` sorted in lexicographical order. Without this class, iterating over all monomials at some fixed degree - requires the use of a double loop of the form: + requires the use of a double loop of the form:: + + idx = 0 + for n0 in range(deg, -1, -1): + for n1 in range(deg-n0, -1, -1): + n2 = deg - n0 - n1 - idx = 0 - for n0 in range(deg, -1, -1): - for n1 in range(deg-n0, -1, -1): - n2 = deg - n0 - n1 + ... # do something with exponents (n0,n1,n2) - ... # do something with exponents (n0,n1,n2) + idx += 1 - idx += 1 + Instead, with this class, these lines can be reduced to:: - Instead, with this class, these lines can be reduced to - myiter = iter(TrivariateMonomialIndices(deg=2)) - for idx, n0, n1, n2 in myiter: - ... # do something with exponents (n0, n1, n2) + myiter = iter(TrivariateMonomialIndices(deg=2)) + for idx, n0, n1, n2 in myiter: + ... # do something with exponents (n0, n1, n2) """ @@ -52,17 +53,16 @@ def __next__(self): raise StopIteration def find_idx(self, exponents): - """ - + """Finds the index of a given tuple of exponents. Parameters ---------- - exponents : 3-tuple (n0, n1, n2) - The exponents of the monomial x^n0 * y^n1 * z^n2 + exponents : 3-tuple of the form (n0, n1, n2) + The exponents of the monomial :math:`x^{n_0} y^{n_1} z^{n_2}` Returns ------- - The index of the tuple in the lexicographical order + The index of the tuple in lexicographical order """ assert n0 + n1 + n2 == self.deg diff --git a/anisoap/utils/spherical_to_cartesian.py b/anisoap/utils/spherical_to_cartesian.py index 8884e9d..bd437ed 100644 --- a/anisoap/utils/spherical_to_cartesian.py +++ b/anisoap/utils/spherical_to_cartesian.py @@ -15,15 +15,25 @@ def prefact_minus1(l): - """ - Parameters: - - l: int + r"""Computes the prefactor that multiplies the :math:`T_{l-1}^\text{th}` term in the iteration. + + For :math:`m \in \left[-l, -l+2, ..., l \right]`, compute the factor as + + .. math:: - Returns: - - A list of size (2*l +1) corresponding to the prefactor that multiplies the T_{l-1} term in the iteration + \left( \frac{\sqrt{(l+1-m)!}}{(l+1+m)!} \right) \left( \frac{\sqrt{(l+m)!}}{(l-m)!} \right) + \left( \frac{2l+1}{l+1-m} \right) - For m in [-l, -l+2, ..., l], compute the factor as : - sqrt(factorial(l+1-m)/ factorial(l+1+m)) sqrt(factorial(l+m)/ factorial(l-m)) (2*l+1)/(l+1-m) + Parameters + ---------- + l : int + Term immediately proceeding the term for which the prefactor is computed. + + Returns + ------- + list of size (2l + 1) + corresponds to the prefactor that multiplies the :math:`T_{l-1}^\text{th}` + term in the iteration """ m = np.arange(-l, l + 1) @@ -36,15 +46,25 @@ def prefact_minus1(l): def prefact_minus2(l): - """ - Parameters: - - l: int + r"""Computes the prefactor that multiplies the :math:`T_{l-2}^\text{th}` term in the iteration. + + For :math:`m \in \left[-l+1, -l+2, ..., l-1\right]`, compute the factor as - Returns: - - A list of size (2*l -1) corresponding to the prefactor that multiplies the T_{l-2} term in the iteration + .. math:: + + \left( \frac{\sqrt{(l+1-m)!}}{(l+1+m)!} \right) \left(\frac{\sqrt{(l-1+m)!}}{(l-1-m)!} \right) + \left( \frac{l+m}{l-m+1} \right) + + Parameters + ---------- + l : int + Term two places after the term for which the prefactor is computed + + Returns + ------- + list of size (2l - 1) + Corresponds to the prefactor that multiplies the term in question - For m in [-l+1, -l+2, ..., l-1], compute the factor as : - sqrt(factorial(l+1-m)/ factorial(l+1+m)) sqrt(factorial(l-1+m)/ factorial(l-1-m)) (l+m)/(l-m+1) """ m = np.arange(-l + 1, l) return ( @@ -65,14 +85,16 @@ def binom(n, k): def spherical_to_cartesian(lmax, num_ns): """ - Finds the coefficients for the cartesian polynomial form of solid harmonics R_{lm} = sqrt((4pi)/(2l+1))*r^l*Y_{lm}. - Note that our AniSOAP expansion does not contain the sqrt((4pi)/(2l+1)), so in calculating expansion coefficients, - we need to divide by that coefficient. - Args: - lmax: - num_ns: - - Returns: + Finds the coefficients for the cartesian polynomial form of solid harmonics + :math:`R_{lm} = sqrt((4pi)/(2l+1))*r^l*Y_{lm}`. Note that our AniSOAP + expansion does not contain the sqrt((4pi)/(2l+1)), so in calculating + expansion coefficients, we need to divide by that coefficient. + + Parameters + ---------- + lmax : int + + num_ns : int """ assert len(num_ns) == lmax + 1 diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d0c3cbf --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = source +BUILDDIR = build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..747ffb7 --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=source +set BUILDDIR=build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..68e0911 --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1,16 @@ +# sphinx dependencies +Sphinx==7.2.6 +sphinx-rtd-theme==2.0.0 + +--find-links https://download.pytorch.org/whl/torch_stable.html +torch==1.11.0+cpu +--find-links https://download.pytorch.org/whl/torch_stable.html +torchvision==0.12.0+cpu +metatensor >=0.1.0,<0.2.0 + +tqdm +chemfiles +matplotlib +skmatter +ase +rascaline @ git+https://github.com/Luthaf/rascaline diff --git a/docs/source/_static/custom.css b/docs/source/_static/custom.css new file mode 100644 index 0000000..319be8b --- /dev/null +++ b/docs/source/_static/custom.css @@ -0,0 +1,3 @@ +.wy-nav-content { + max-width: 90%; +} diff --git a/docs/source/api.rst b/docs/source/api.rst new file mode 100644 index 0000000..bcc56fc --- /dev/null +++ b/docs/source/api.rst @@ -0,0 +1,37 @@ +API Reference +============= + +Representations +--------------- +.. automodule:: anisoap.representations.ellipsoidal_density_projection + :members: + :show-inheritance: + :inherited-members: +.. automodule:: anisoap.representations.radial_basis + :members: + :show-inheritance: + :inherited-members: + +Utilities +--------- +.. automodule:: anisoap.utils.metatensor_utils + :members: + :show-inheritance: + :inherited-members: +.. automodule:: anisoap.utils.moment_generator + :members: + :show-inheritance: + :inherited-members: +.. automodule:: anisoap.utils.monomial_iterator + :members: + :show-inheritance: + :inherited-members: +.. automodule:: anisoap.utils.spherical_to_cartesian + :members: + :show-inheritance: + :inherited-members: +.. automodule:: anisoap.utils.shortcuts + :members: + :show-inheritance: + :inherited-members: + diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100644 index 0000000..6213f5f --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,48 @@ +import os, sys + +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Project information ----------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information +root_doc = "index" + +project = "AniSOAP" +copyright = "2023, Cersonsky Lab" +author = "The Cersonsky Lab" +release = "0.0.1" + +# -- General configuration --------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration + +extensions = [ + "sphinx.ext.duration", + "sphinx.ext.doctest", + "sphinx.ext.autosummary", + "sphinx.ext.autodoc", + "sphinx.ext.napoleon", + "sphinx.ext.mathjax", + "sphinx.ext.coverage", +] + +sys.path.insert(0, os.path.abspath("../../")) + +templates_path = ["_templates"] +exclude_patterns = [] + +# -- Napoleon settings ------------------------------------------------------- +napoleon_google_docstring = False +napoleon_numpy_docstring = True + +# -- Options for HTML output ------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output + +html_theme = "sphinx_rtd_theme" + +html_static_path = ["_static"] + +html_css_files = [ + "custom.css", +] diff --git a/docs/source/getting_started.rst b/docs/source/getting_started.rst new file mode 100644 index 0000000..b0976cc --- /dev/null +++ b/docs/source/getting_started.rst @@ -0,0 +1,17 @@ +=============== +Getting Started +=============== + +Introduction +------------ +Todo + +First Steps +----------- +AniSOAP depends extensively on the packages `rascaline +`_ +and `metatensor `_. +It may be helpful to familiarize yourself with these packages before diving into +AniSOAP. + + diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 0000000..5a69110 --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,28 @@ +=================================== +Welcome to AniSOAP's documentation! +=================================== + +**AniSOAP** is a Python library for creating descriptors of chemical systems +suitable for machine learning use. This project aims to extend SOAP +descriptors by computing anisotrophy extensions. + +.. note:: + + This project is under active development. + +.. note:: + + AniSOAP documentation is a work in progress; please contact us with questions + not yet addressed here. + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + installation + getting_started + tutorials + usage + api + + diff --git a/docs/source/installation.rst b/docs/source/installation.rst new file mode 100644 index 0000000..f054f92 --- /dev/null +++ b/docs/source/installation.rst @@ -0,0 +1,44 @@ +============ +Installation +============ + +Dependencies +------------ + +Before installing AniSOAP, please make sure you have the following installed\: + +* `Python: 3.6 or higher `_ +* `numPy: 1.13 or higher `_ +* `sciPy: 1.4.0 or higher `_ +* `Atomic Simulation Environment (ASE): 3.18 or higher `_ +* `Metatensor `_ +* `Rascaline `_ +* `Rust -- We reccommend using 'rustup `_ + + +Installing AniSOAP +------------------ + +Navigate to the directory where you would like the AniSOAP package to be located, then copy and paste the +following into your shell:: + + git clone https://github.com/cersonsky-lab/AniSOAP + +Then navigate to the AniSOAP directory with:: + + cd AniSOAP + +Now use pip to install the library:: + + pip install . + + +Testing +------- + +AniSOAP is still under active development, so you may want to run some tests to ensure that your installation is working properly. From the main directory you can run the internal tests with:: + + pytest tests/. + + + diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst new file mode 100644 index 0000000..fde3cf7 --- /dev/null +++ b/docs/source/tutorials.rst @@ -0,0 +1,8 @@ +Tutorials +========= + +Here is a collection of tutorials that may prove useful, including some interactive Jupyter Notebooks. + +.. Note:: + + More tutorials are currently being developed. Stay tuned! diff --git a/docs/source/usage.rst b/docs/source/usage.rst new file mode 100644 index 0000000..6f3336c --- /dev/null +++ b/docs/source/usage.rst @@ -0,0 +1,4 @@ +Usage +===== + +Check back soon!