diff --git a/iodata/basis.py b/iodata/basis.py index c9aff450a..c8e81cb62 100644 --- a/iodata/basis.py +++ b/iodata/basis.py @@ -24,6 +24,7 @@ import attr import numpy as np +from .overlap_cartpure import tfs from .attrutils import validate_shape @@ -37,11 +38,13 @@ def _alsolist(f): """Wrap a function to accepts also list as first argument and then return list.""" + @wraps(f) def wrapper(firsts, *args, **kwargs): if isinstance(firsts, (Integral, str)): return f(firsts, *args, **kwargs) return [f(first, *args, **kwargs) for first in firsts] + return wrapper @@ -88,7 +91,7 @@ def angmom_its(angmom: Union[int, List[int]]) -> Union[str, List[str]]: @attr.s(auto_attribs=True, slots=True, on_setattr=[attr.setters.validate, attr.setters.convert]) class Shell: - """A shell in a molecular basis representing (generalized) contractions with the same exponents. + """Contracted shell representing (generalized) contractions with the same exponents. Attributes ---------- @@ -102,7 +105,8 @@ class Shell: and 'p' for pure. Pure functions are only allowed for angmom>1. The length equals the number of contractions: len(angmoms)=ncon. exponents - The array containing the exponents of the primitives, with shape (nprim,). + The array containing the exponents of the primitives, with shape (nprim,), + where nprim is the number of primitive contracted shells. coeffs The array containing the coefficients of the normalized primitives in each contraction; shape = (nprim, ncon). @@ -110,6 +114,10 @@ class Shell: (densities) normalized, but contractions are not necessarily normalized. (This depends on the code which generated the contractions.) + Notes + ----- + Basis set conventions and terminology are documented in :ref:`basis_conventions`. + """ icenter: int @@ -119,7 +127,7 @@ class Shell: coeffs: np.ndarray = attr.ib(validator=validate_shape(("exponents", 0), ("kinds", 0))) @property - def nbasis(self) -> int: # noqa: D401 + def nbasis(self) -> int: # noqa: D401 """Number of basis functions (e.g. 3 for a P shell and 4 for an SP shell).""" result = 0 for angmom, kind in zip(self.angmoms, self.kinds): @@ -132,20 +140,23 @@ def nbasis(self) -> int: # noqa: D401 return result @property - def nprim(self) -> int: # noqa: D401 - """Number of primitives, also known as the contraction length.""" + def nprim(self) -> int: # noqa: D401 + """Number of primitives in the contracted shells, also known as the contraction length""" return len(self.exponents) @property - def ncon(self) -> int: # noqa: D401 - """Number of contractions. This is usually 1; e.g., it would be 2 for an SP shell.""" + def ncon(self) -> int: # noqa: D401 + """Number of contractions with distinct angular momentum numbers. + + This is usually 1; e.g., it would be 2 for an SP shell. + """ return len(self.angmoms) @attr.s(auto_attribs=True, slots=True, on_setattr=[attr.setters.validate, attr.setters.convert]) class MolecularBasis: - """A complete molecular orbital or density basis set. + """A complete molecular orbital or density basis set as a collection of contracted shells. Attributes ---------- @@ -190,6 +201,19 @@ class MolecularBasis: The normalization convention of primitives, which can be 'L2' (orbitals) or 'L1' (densities) normalized. + Methods + ------- + get_segmented + Return molecular basis object that is segmented. + get_decontracted + Return molecular basis object that is decontracted. + convert_kind + Return molecular basis object whose kinds are the same type. + + Notes + ----- + Basis set conventions and terminology are documented in :ref:`basis_conventions`. + """ shells: List[Shell] @@ -197,12 +221,12 @@ class MolecularBasis: primitive_normalization: str @property - def nbasis(self) -> int: # noqa: D401 + def nbasis(self) -> int: # noqa: D401 """Number of basis functions.""" return sum(shell.nbasis for shell in self.shells) def get_segmented(self): - """Unroll generalized contractions.""" + """Convert Generalized Molecular Basis to Segmented Molecular Basis.""" shells = [] for shell in self.shells: for angmom, kind, coeffs in zip(shell.angmoms, shell.kinds, shell.coeffs.T): @@ -211,6 +235,139 @@ def get_segmented(self): # pylint: disable=no-member return attr.evolve(self, shells=shells) + def get_decontracted(self): + r"""Convert to Decontracted Molecular Basis from a Molecular Basis.""" + shells = [] + for shell in self.shells: + for i, (angmom, kind) in enumerate(zip(shell.angmoms, shell.kinds)): + for exponent, coeff in zip(shell.exponents, shell.coeffs[:, i]): + shells.append( + Shell(shell.icenter, [angmom], [kind], np.array([exponent]), + coeff.reshape(-1, 1)) + ) + # pylint: disable=no-member + return attr.evolve(self, shells=shells) + + def convert_kind(self, to: str): + r""" + Convert all contractions from Pure (or Cartesian) kind to Cartesian (or Pure). + + Parameters + ---------- + to + Specifying which one to convert to. It is either "c" for Cartesian or "p" for Pure. + + Raises + ------ + ValueError : + If "to" is not recognized to be either "c" or "p". + + Notes + ----- + - Some of these conversion will cause the contraction coefficients to be + different across different shells. In such a scenario, + this function will average all coefficients so that the shell-type inside a contracted + shell will have the same coefficients. + + See Also + -------- + convert_primitive_kind : Converts primitive shells without averaging. + + """ + if to not in ("c", "p"): + raise ValueError("The to string was not recognized: %s" % to) + + shells = [] + for shell in self.shells: + kind_new = [] + coeffs_new = [] + + for angmom, kind, coeffs in zip(shell.angmoms, shell.kinds, shell.coeffs.T): + # Need to convert the kind of shell. S-type doesn't need to change. + if kind != to and angmom >= 2: + if kind == "c": + b = np.repeat(coeffs[np.newaxis], tfs[angmom].shape[1], axis=0) + coeff = tfs[angmom].dot(b) + # Take the mean, since the coefficients vary between primitives. + coeff = np.mean(coeff, axis=0) + elif kind == "p": + b = np.repeat(coeffs[np.newaxis], tfs[angmom].shape[0], axis=0) + # Solve it using least-squared projection. + coeff = np.linalg.lstsq(tfs[angmom], b, rcond=None)[0] + # Take the mean, since the coefficients vary between primitives. + coeff = np.mean(coeff, axis=0) + kind_new.append(to) + coeffs_new.append(coeff) + # No need to convert. + else: + kind_new.append(to) + coeffs_new.append(coeffs) + + shells.append( + Shell( + shell.icenter, shell.angmoms, kind_new, shell.exponents, np.array(coeffs_new).T + ) + ) + # pylint: disable=no-member + return attr.evolve(self, shells=shells) + + +def convert_primitive_kind(angmom: int, kind: str, coeffs: np.ndarray, to: str) -> np.ndarray: + r""" + Convert coefficients in Cartesian to Pure or vice-versa of a Primitive shell. + + Parameters + ---------- + angmom + Integer representing the total angular momentum number. + kind + String representing "p" for Pure primitives and "c" for Cartesian primitives. + coeffs + Coefficients of the contraction of type `kind`. Shape is (N,) where + N is the number of primitives with angular momentum `angmom` and type `kind`. + For example, angmom=2 and kind="c" implies N=6, + The order should match the conventions. + to + Specifying which one to convert to. It is either "c" for Cartesian or "p" for Pure. + + Returns + ------- + new_coeffs + Coefficients in the new basis (either Cartesian or Pure). + + Examples + -------- + Suppose one has the following contraction with angular momentum 2 and kind is Cartesian with + integers (2, 0, 0) and (1, 1, 0): + + .. math:: + c_1 x^2 y^0 z^0 e^{-\alpha r^2} + c_2 x y z^0 e^{-\alpha r^2} + + Then the coefficients of this in the right convention [xx, xy, xz, yy, yz, zz] is + + >> coeff = [c_1, c_2, 0, 0, 0, 0] + + The basis function conventions dictionary are documented in :ref:`basis_conventions and + in the MolecularBasis class documentation. + + To convert this Cartesian contraction to Pure type, one would do + + >> new_coeffs = convert_primitive_kind(1, "c", coeff, "p") + + """ + if to not in ("p", "c"): + raise ValueError("The to string was not recognized: %s" % to) + if kind not in ("p", "c"): + raise ValueError("The kind string was not recognized: %s" % kind) + if to != kind: + if kind == "c": + coeff = tfs[angmom].dot(coeffs) + else: + # Solve it using least-squared projection. + coeff = np.linalg.lstsq(tfs[angmom], coeffs, rcond=None)[0] + return coeff + return coeffs + def convert_convention_shell(conv1: List[str], conv2: List[str], reverse=False) \ -> Tuple[np.ndarray, np.ndarray]: diff --git a/iodata/overlap.py b/iodata/overlap.py index 75a71de3c..72e527657 100644 --- a/iodata/overlap.py +++ b/iodata/overlap.py @@ -26,7 +26,9 @@ from .basis import convert_conventions, iter_cart_alphabet, MolecularBasis from .basis import HORTON2_CONVENTIONS as OVERLAP_CONVENTIONS -__all__ = ['OVERLAP_CONVENTIONS', 'compute_overlap', 'gob_cart_normalization'] +__all__ = [ + 'OVERLAP_CONVENTIONS', 'compute_overlap', 'gob_cart_normalization', "convert_vector_basis" +] # pylint: disable=too-many-nested-blocks,too-many-statements @@ -223,3 +225,64 @@ def gob_cart_normalization(alpha: np.ndarray, n: np.ndarray) -> np.ndarray: """ return np.sqrt((4 * alpha) ** sum(n) * (2 * alpha / np.pi) ** 1.5 / np.prod(factorial2(2 * n - 1))) + + +def convert_vector_basis( + coeffs1: np.ndarray, + basis2_overlap: np.ndarray, + basis21_overlap: np.ndarray +) -> np.ndarray: + r""" + Convert a vector from basis 1 of size M to another basis 2 of size N. + + Parameters + ---------- + coeffs1: + Coefficients of the vector expanded in basis set 1. Shape is (M,). + basis2_overlap: + Symmetric matrix whose entries are the inner-product between basis vectors + inside basis set 2. Shape is (N, N). + basis21_overlap: + The overlap matrix between basis set 1 and basis set 2. Shape is (N, M). + + Returns + ------- + coeffs2 : + Coefficients of the vector expanded in basis set 2. Shape is (N,). + + Raises + ------ + ValueError : + If shapes of the matrices don't match the requirements in the docs. + LinAlgError : + If least squared solution does not converge. + + Notes + ----- + - Basis vectors are defined to be linearly independent wrt to one another + and need not be orthogonal. + + - `basis2_overlap` is the matrix with (i, j)th entries + :math:`\braket{\psi_i, \psi_j}` where :math:`\psi_i` are in basis set 2. + + - `basis21_overlap` is the matrix whose (i, j)th entries are + :math:`\braket{\psi_i, \phi_j}` where :math:`\phi_j` are in basis set 1 and + :math:`\psi_i` are in basis set 2. + + - If `basis2_overlap` is not full rank, then least-squared solution is solved instead. This + is effectively solving :math:`||b - Ax||`. + + """ + if basis2_overlap.shape[0] != basis2_overlap.shape[1]: + raise ValueError("The `basis2_overlap` should be a square matrix.") + if np.any(np.abs(basis2_overlap.T - basis2_overlap) > 1e-10): + raise ValueError("The `basis2_overlap` should be a symmetric matrix.") + + b = basis21_overlap.dot(coeffs1) + try: + # Try solving exact solution if `basis12_overlap` is full rank. + coeffs2 = np.linalg.solve(basis2_overlap, b) + except np.linalg.LinAlgError: + # Try least-squared solution. + coeffs2, _, _, _ = np.linalg.lstsq(basis2_overlap, b) + return coeffs2 diff --git a/iodata/test/test_basis.py b/iodata/test/test_basis.py index 73140fcdf..c4038099e 100644 --- a/iodata/test/test_basis.py +++ b/iodata/test/test_basis.py @@ -21,11 +21,11 @@ import attr import numpy as np -from numpy.testing import assert_equal +from numpy.testing import assert_equal, assert_allclose import pytest from ..basis import (angmom_sti, angmom_its, Shell, MolecularBasis, - convert_convention_shell, convert_conventions, + convert_convention_shell, convert_conventions, convert_primitive_kind, iter_cart_alphabet, HORTON2_CONVENTIONS, PSI4_CONVENTIONS) from ..formats.cp2klog import CONVENTIONS as CP2K_CONVENTIONS @@ -228,6 +228,146 @@ def test_convert_convention_shell(): assert_equal(vec2[permutation] * signs, vec1) +def test_get_decontracted(): + obasis = MolecularBasis( + [Shell(0, [0], ['c'], [0.01], np.array([[3.]])), + Shell(1, [0, 1], ['c', 'c'], [0.03, 0.04], np.array([[3., 4.], [1., 2.]])), + Shell(0, [2, 1, 2], ['p', 'c', 'c'], [0.05], np.array([[5., 5., 6.]]))], + {(0, 'c'): ['s'], + (1, 'c'): ['x', 'z', '-y'], + (2, 'p'): ['dc0', 'dc1', '-ds1', 'dc2', '-ds2']}, + 'L2' + ) + obasis = obasis.get_decontracted() + + ncontshells = 8 # Number of contracted shells now is equal to the number of primitives. + assert_equal(len(obasis.shells), ncontshells) + # Assert each shell has length one exponents + for _ in range(0, ncontshells): + assert_equal(len(obasis.shells[0].exponents), 1) + # Assert they have the right exponents. + assert_equal(obasis.shells[0].exponents, 0.01) + assert_equal(obasis.shells[1].exponents, 0.03) + assert_equal(obasis.shells[2].exponents, 0.04) + assert_equal(obasis.shells[3].exponents, 0.03) + assert_equal(obasis.shells[4].exponents, 0.04) + assert_equal(obasis.shells[5].exponents, 0.05) + assert_equal(obasis.shells[6].exponents, 0.05) + assert_equal(obasis.shells[7].exponents, 0.05) + # Assert they have the right coefficients. + assert_equal(obasis.shells[0].coeffs, 3.) + assert_equal(obasis.shells[1].coeffs, np.array([[3.]])) + assert_equal(obasis.shells[2].coeffs, np.array([[1.]])) + assert_equal(obasis.shells[3].coeffs, np.array([[4.]])) + assert_equal(obasis.shells[4].coeffs, np.array([[2.]])) + assert_equal(obasis.shells[5].coeffs, np.array([[5.]])) + assert_equal(obasis.shells[6].coeffs, np.array([[5.]])) + assert_equal(obasis.shells[7].coeffs, np.array([[6.]])) + # Assert right centers + assert_equal(obasis.shells[0].icenter, 0) + assert_equal(obasis.shells[1].icenter, 1) + assert_equal(obasis.shells[2].icenter, 1) + assert_equal(obasis.shells[3].icenter, 1) + assert_equal(obasis.shells[4].icenter, 1) + assert_equal(obasis.shells[5].icenter, 0) + assert_equal(obasis.shells[6].icenter, 0) + assert_equal(obasis.shells[7].icenter, 0) + # Assert right kind + assert_equal(obasis.shells[0].kinds, ['c']) + assert_equal(obasis.shells[1].kinds, ['c']) + assert_equal(obasis.shells[2].kinds, ['c']) + assert_equal(obasis.shells[3].kinds, ['c']) + assert_equal(obasis.shells[4].kinds, ['c']) + assert_equal(obasis.shells[5].kinds, ['p']) + assert_equal(obasis.shells[6].kinds, ['c']) + assert_equal(obasis.shells[7].kinds, ['c']) + + +def test_convert_kind_to_cartesian(): + obasis = MolecularBasis( + [Shell(0, [0], ['c'], [0.01], np.array([[3.]])), + Shell(1, [0, 1], ['c', 'c'], [0.03, 0.04], np.array([[3., 4.], [1., 2.]])), + Shell(0, [2, 1], ['p', 'c'], [0.05], np.array([[5., 6.]]))], + {(0, 'c'): ['s'], + (1, 'c'): ['x', 'z', '-y'], + (2, 'p'): ['dc0', 'dc1', '-ds1', 'dc2', '-ds2']}, + 'L2' + ) + obasis = obasis.convert_kind("c") + # Assert number of shells is the same. + assert_equal(len(obasis.shells), 3) + # Assert the already cartesian are still cartesian. + assert_equal(obasis.shells[0].kinds, ["c"]) + assert_equal(obasis.shells[1].kinds, ["c", "c"]) + assert_equal(obasis.shells[2].kinds, ["c", "c"]) + # Assert the exponents should not change. + assert_equal(obasis.shells[0].exponents, [0.01]) + assert_equal(obasis.shells[1].exponents, [0.03, 0.04]) + assert_equal(obasis.shells[2].exponents, [0.05]) + # Assert that the coefficients expect for the pure is the same. + assert_equal(obasis.shells[0].coeffs, np.array([[3.]])) + assert_equal(obasis.shells[1].coeffs, np.array([[3., 4.], [1., 2.]])) + + # Assert that the coefficient of the pure-shell inside the third contraction changed to cart. + # First solve for the coefficients, this was done manually. + coeff_pure = np.array([[1.22008468], [5.], [5.], [-4.55341801], [5.], [3.33333333]]) + assert_allclose(obasis.shells[2].coeffs, np.array([[np.mean(coeff_pure), 6.0]])) + + # Convert Easier Example where nothing happens. + obasis = MolecularBasis( + [Shell(0, [1], ['p'], [0.05, 0.01], np.array([[5.], [10.]]))], None, 'L2' + ) + obasis = obasis.convert_kind("c") + assert_equal(len(obasis.shells), 1) + assert_equal(obasis.shells[0].kinds, ["c"]) + assert_equal(obasis.shells[0].exponents, [0.05, 0.01]) + assert_equal(obasis.shells[0].coeffs, np.array([[5.], [10.]])) + + +def test_convert_kind_to_pure(): + obasis = MolecularBasis( + [Shell(0, [0], ['c'], [0.01], np.array([[3.]])), + Shell(1, [1, 2], ['c', 'c'], [0.03, 0.04], np.array([[3., 4.], [1., 2.]]))], + {(0, 'c'): ['s'], + (1, 'c'): ['x', 'z', '-y'], + (2, 'p'): ['dc0', 'dc1', '-ds1', 'dc2', '-ds2']}, + 'L2' + ) + obasis = obasis.convert_kind("p") + assert_equal(len(obasis.shells), 2) + assert_equal(obasis.shells[0].coeffs, np.array([[3.]])) + new_coeff = np.mean(np.array([[0.], [4.], [4.], [0.], [4.]])) + new_coeff2 = np.mean(np.array([[0.], [2.], [2.], [0.], [2.]])) + desired = np.array([[3., new_coeff], [1., new_coeff2]]) + assert_equal(obasis.shells[1].coeffs, desired) + + +def test_convert_primitive_kind(): + # Test converting P-type from kind="c" to kind="p". + coeff = np.array([1., 2., 3.]) + coeff2 = convert_primitive_kind(1, "c", coeff, "p") + desired = np.array([3., 1., 2.]) + assert_equal(coeff2, desired) + + # Test converting D-type from kind="c" to kind="p". + coeff = np.array([1., 2., 3., 4., 5., 6.]) + coeff = convert_primitive_kind(2, "c", coeff, "p") + desired = np.array([3.5, 3., 5., -2.59807621, 2.]) + assert_allclose(coeff, desired) + + # Test converting P-type from kind="p" to kind="c". + coeff = np.array([1., 2., 3.]) + coeff2 = convert_primitive_kind(1, "p", coeff, "c") + desired = np.array([2., 3., 1.]) + assert_equal(coeff2, desired) + + # Test converting D-type from kind="p" to kind="c". + coeff = np.array([1., 2., 3., 4., 5.]) + ccoeff = convert_primitive_kind(2, "p", coeff, "c") + desired = np.array([1.97606774, 5., 2., -2.64273441, 3., 0.66666667]) + assert_allclose(ccoeff, desired) + + def test_convert_convention_obasis(): obasis = MolecularBasis( [Shell(0, [0], ['c'], np.zeros(3), np.zeros((3, 1))), diff --git a/iodata/test/test_overlap.py b/iodata/test/test_overlap.py index 797d7de49..b574a5bc2 100644 --- a/iodata/test/test_overlap.py +++ b/iodata/test/test_overlap.py @@ -20,12 +20,13 @@ import attr import numpy as np -from numpy.testing import assert_allclose +from numpy.testing import assert_allclose, assert_equal from pytest import raises from ..api import load_one from ..basis import MolecularBasis, Shell -from ..overlap import compute_overlap, OVERLAP_CONVENTIONS +from ..overlap import compute_overlap, OVERLAP_CONVENTIONS, convert_vector_basis +from ..overlap_cartpure import tfs try: from importlib_resources import path @@ -85,3 +86,76 @@ def test_overlap_l1(): atcoords = np.zeros((1, 3)) with raises(ValueError): _ = compute_overlap(dbasis, atcoords) + + +def test_converting_between_orthonormal_basis_set(): + # Test converting from basis set of M Elements to N Elements where M <= N. + # All basis sets are assumed orthonormal and 1st basis set is assumed to be + # in the second basis set. + M = np.random.randint(1, 25) + N = np.random.randint(M, 50) + coeffs1 = np.random.random(M) + basis2_overlap = np.eye(N) # Since it is orthonormal, it is identity. + + basis21_overlap = np.zeros((N, M)) + for i in range(0, M): + basis21_overlap[i, i] = 1. # Since it is a subset. + + coeffs2 = convert_vector_basis(coeffs1, basis2_overlap, basis21_overlap) + assert_allclose(coeffs1, coeffs2[:M], rtol=1.e-5, atol=1.e-8) + assert_allclose(np.zeros(coeffs2[M:].shape), coeffs2[M:], rtol=1.e-5, atol=1.e-8) + + # Test converting in the reverse direction ie Large to small. + coeffs2 = np.random.random(N) + basis1_overlap = np.eye(M) + basis12_overlap = np.zeros((M, N)) + for i in range(0, M): + basis12_overlap[i, i] = 1. # Since it is a subset. + + coeffs1 = convert_vector_basis(coeffs2, basis1_overlap, basis12_overlap) + assert_equal(len(coeffs1), M) + assert_allclose(coeffs1, coeffs2[:M], rtol=1.e-5, atol=1.e-8) + + +def test_converting_from_cartesian_to_pure(): + # Test converting simple one coefficient, S-type. + overlap_cp = tfs[0] + coeff = np.array([5.]) + coeff2 = convert_vector_basis(coeff, np.eye(1), overlap_cp) + assert_allclose(coeff, coeff2, rtol=1.e-5, atol=1.e-8) + + # Test converting p-type. + overlap_cp = tfs[1] + coeff = np.array([1., 2., 3.]) + coeff2 = convert_vector_basis(coeff, np.eye(3), overlap_cp) + desired = np.array([3., 1., 2.]) + assert_allclose(coeff2, desired, rtol=1.e-5, atol=1.e-8) + + # Test converting d-type. + overlap_cp = tfs[2] + coeff = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]) + coeff2 = convert_vector_basis(coeff, np.eye(5), overlap_cp) + desired = overlap_cp.dot(coeff) + assert_allclose(coeff2, desired, rtol=1.e-5, atol=1.e-8) + + +def test_converting_from_pure_to_cartesian(): + # Test converting S-type. + overlap_cp = tfs[0] + coeff = np.array([5.]) + coeff2 = convert_vector_basis(coeff, np.eye(1), overlap_cp.T) + assert_allclose(coeff, coeff2, rtol=1.e-5, atol=1.e-8) + + # Test converting P-type. + overlap_cp = tfs[1] + coeff = np.array([1., 2., 3.]) + coeff2 = convert_vector_basis(coeff, np.eye(3), overlap_cp.T) + desired = np.linalg.lstsq(overlap_cp, coeff, rcond=None)[0] + assert_allclose(coeff2, desired, rtol=1.e-5, atol=1.e-8) + + # Test converting D-type. + overlap_cp = tfs[2] + pcoeff = np.array([1., 2., 3., 4., 5.]) + ccoeff = convert_vector_basis(pcoeff, np.eye(6), np.linalg.pinv(overlap_cp)) + desired = np.linalg.lstsq(overlap_cp, pcoeff, rcond=None)[0] + assert_allclose(ccoeff, desired, rtol=1.e-5, atol=1.e-8)