From a0002d7dc209773611753f5d8c22453ac81ab4b9 Mon Sep 17 00:00:00 2001 From: Johannes Steinmetzer Date: Mon, 22 Jan 2024 16:32:08 +0100 Subject: [PATCH] add: make (parts) of sympleint testable add: basic Boys-function implementation for testing add: int2c-driver for testing add: tests against pyscf/libcint add: cart2sph module from pysisyphus to remove dependecy on it add: improved generation of spherical integrals by using more precise transformation factors generated using mpmath add: sympleints.defs.overlap.gen_prefactor_shell for GDMA --- ressources/cgto_normalization.ipynb | 2 +- sympleints/Functions.py | 10 +- sympleints/PythonRenderer.py | 4 +- sympleints/Renderer.py | 16 +- sympleints/boys.py | 147 ++++++++++++++ sympleints/cart2sph.py | 270 +++++++++++++++++++++++++ sympleints/config.py | 2 +- sympleints/defs/__init__.py | 3 +- sympleints/defs/multipole.py | 2 +- sympleints/defs/overlap.py | 28 +++ sympleints/defs/twocenter1d.py | 43 ++++ sympleints/defs/twocenter3d.py | 68 +++++++ sympleints/main.py | 276 +++++++++++++++++--------- sympleints/symbols.py | 16 ++ sympleints/templates/py_equi_func.tpl | 7 +- sympleints/templates/py_func.tpl | 10 +- sympleints/templates/py_module.tpl | 4 +- sympleints/testing/__init__.py | 0 sympleints/testing/intor.py | 76 +++++++ sympleints/testing/normalization.py | 43 ++++ sympleints/testing/pyscf_interface.py | 101 ++++++++++ sympleints/testing/shells.py | 26 +++ tests/test_boys/test_boys.py | 39 ++++ tests/test_intor/test_intor.py | 79 ++++++++ 24 files changed, 1147 insertions(+), 125 deletions(-) create mode 100644 sympleints/boys.py create mode 100755 sympleints/cart2sph.py create mode 100644 sympleints/defs/twocenter3d.py create mode 100644 sympleints/testing/__init__.py create mode 100644 sympleints/testing/intor.py create mode 100644 sympleints/testing/normalization.py create mode 100644 sympleints/testing/pyscf_interface.py create mode 100644 sympleints/testing/shells.py create mode 100644 tests/test_boys/test_boys.py create mode 100644 tests/test_intor/test_intor.py diff --git a/ressources/cgto_normalization.ipynb b/ressources/cgto_normalization.ipynb index ece223c..8e40b63 100644 --- a/ressources/cgto_normalization.ipynb +++ b/ressources/cgto_normalization.ipynb @@ -661,7 +661,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.10" + "version": "3.11.5" } }, "nbformat": 4, diff --git a/sympleints/Functions.py b/sympleints/Functions.py index 601200f..b7a30a3 100644 --- a/sympleints/Functions.py +++ b/sympleints/Functions.py @@ -26,7 +26,7 @@ class Functions: header: str = "" doc_func: Optional[Callable] = None comment: str = "" - boys: bool = False + boys_func: Optional[str] = None ncomponents: int = 0 with_ref_center: bool = True full_name: Optional["str"] = None @@ -141,7 +141,8 @@ def ls(self): @property def nbfs(self): - return len(self.coeffs) + # Return number of total angular momenta + return len(self.ls_exprs[0][0]) @property def ls_name_map(self): @@ -170,8 +171,3 @@ def ndim_act(self): @property def cartesian(self): return not self.spherical - - # def numba_func_type(self): - # pass - - # def numba_driver_type(self): diff --git a/sympleints/PythonRenderer.py b/sympleints/PythonRenderer.py index f7f462e..62639b1 100644 --- a/sympleints/PythonRenderer.py +++ b/sympleints/PythonRenderer.py @@ -81,6 +81,7 @@ def render_equi_function( name, equi_name, equi_inds, + reshape, from_axes, to_axes, ): @@ -96,6 +97,7 @@ def render_equi_function( args=args, from_axes=from_axes, to_axes=to_axes, + reshape=reshape, ) return rendered @@ -110,7 +112,7 @@ def render_module(self, functions, rendered_funcs, **tpl_kwargs): _tpl_kwargs = { "header": functions.header, "comment": functions.comment, - "boys": functions.boys, + "boys": functions.boys_func, "funcs": rendered_funcs, "func_dict": func_dict, } diff --git a/sympleints/Renderer.py b/sympleints/Renderer.py index 49b618e..a45cdca 100644 --- a/sympleints/Renderer.py +++ b/sympleints/Renderer.py @@ -64,6 +64,7 @@ def render_equi_function( name, act_name, equi_inds, + L_tots, from_axes, to_axes, ): @@ -133,16 +134,21 @@ def render_functions(self, functions: Functions): from_axes = tuple(from_axes) more_components = functions.ncomponents > 1 - print( - f"{more_components=}, {functions.ncomponents=}, {self._drop_dim=}" - ) add_axis = more_components or (not self._drop_dim and nbfs == 2) - print(f"{add_axis=}") # , {inds_missing=}") if add_axis: from_axes = tuple([0] + [fa + 1 for fa in from_axes]) to_axes = tuple([0] + [ta + 1 for ta in to_axes]) + + reshape = shape if functions.ncomponents > 1 else shape[1:] + func_equi = self.render_equi_function( - functions, name, name_equi, equi_inds, from_axes, to_axes + functions, + name, + name_equi, + equi_inds, + reshape, + from_axes, + to_axes, ) rendered_funcs.append( RenderedFunction(name=name_equi, Ls=L_tots_equi, text=func_equi) diff --git a/sympleints/boys.py b/sympleints/boys.py new file mode 100644 index 0000000..f2a0717 --- /dev/null +++ b/sympleints/boys.py @@ -0,0 +1,147 @@ +from typing import Callable + +import numpy as np +import scipy as sp +from scipy.special import factorial, factorial2 + + +def boys_quad(n, x, err=1e-13): + """Boys function from quadrature.""" + + def inner(t): + return t ** (2 * n) * np.exp(-x * t**2) + + y, err = sp.integrate.quad(inner, 0.0, 1.0, epsabs=err, epsrel=err) + return y + + +def get_table(nmax, order, dx, xmax): + nsteps = (xmax / dx) + 1 + # assert that mantissa of nsteps is 0.0 + assert ( + float(int(nsteps)) == nsteps + ), f"{xmax=} must be an integer multiple of {dx=:.8e} but it isn't!" + nsteps = int(nsteps) + xs = np.linspace(0.0, xmax, num=nsteps) + table = np.empty_like(xs) + for i, x in enumerate(xs): + table[i] = boys_quad(nmax + order, x) + return xs, table + + +def boys_x0(n): + """Boys-function when x equals 0.0.""" + return 1 / (2 * n + 1) + + +def boys_xlarge(n, x): + """Boys-function for large x values (x >= 30 is sensible).""" + _2n = 2 * n + f2 = factorial2(_2n - 1) if n > 0 else 1.0 + return f2 / 2 ** (n + 1) * np.sqrt(np.pi / x ** (_2n + 1)) + + +def get_boys_func( + nmax: int, order: int = 6, dx: float = 0.1, xlarge: float = 30.0 +) -> Callable[[int, float], float]: + """Wrapper that returns Boys function object. + + 1 + / + | + | 2 + F_n(x) = | 2*n -t x + | t e dt + | + / + 0 + + Parameters + ---------- + nmax + Positive integer Maximum n value for which the Boys function can be evaluated. + order + Order of the Tayor-expansion that is used to evaluate the Boys-function, defaults + to 6. + dx + Positive float; Distance between points on the grid. + xlarge + Cutoff value. For 0 < x < xlarge a Taylor expansion will be used to evaluate + the boys function. For x >= xlarge the Boys-function is evaluated using the + incomplete gamma function. For all n and x = 0.0 the integral is calculated + analytically. + + Returns + ------- + boys + Callable F(n, x) that accepts two arguments: an integer 0 <= n <= nmax and a + float 0.0 <= x. It returns the values of the Boys-function F_n at the given x. + """ + # Pretabulated Boys-values on a grid for n = nmax + order + xs, table = get_table(nmax, order, dx, xlarge) + nxs = len(xs) + + # Build up full table for all n = {0, 1, ... nmax + order} via Downward recursion + table_full = np.empty((nmax + order + 1, nxs)) + table_full[-1] = table + exp_minus_xs = np.exp(-xs) + _2xs = 2 * xs + # Loop over all n-values + for n in range(nmax + order - 1, -1, -1): + table_full[n] = (_2xs * table_full[n + 1] + exp_minus_xs) / (2 * n + 1) + # 1 / k! + _1_factorials = 1 / factorial(np.arange(order, dtype=int)) + + # Prefactors for Boys-values for large x > xlarge + xlarge_prefact = np.empty(nmax + 1) + for n in range(nmax + 1): + xlarge_prefact[n] = ( + (factorial2(2 * n - 1) if n > 0 else 1.0) / 2 ** (n + 1) * np.sqrt(np.pi) + ) + + def boys_xlarge_precalc(n, x): + return xlarge_prefact[n] * x ** (-n - 0.5) + + def boys_taylor(n, x): + # Rounding to the nearest integer always yields a grid point before (or on) x + # so we interpolate and don't extrapolate + factor = int(1 / dx) + ind0 = int(x * factor) + # Closest grid point (on/before) x + xg = ind0 * dx + + # Step length for Taylor expansion; depends on closest grid point + dx_taylor = x - xg + result = 0.0 + # The loop is somehow faster than the vectorized expression + for k in range(order): + result += table_full[n + k, ind0] * (-dx_taylor) ** k * _1_factorials[k] + return result + + def boys(n: int, xs: float) -> float: + def func(n, x): + if x == 0.0: + # Take values directly from table; should correspond to 1 / (2n + 1) + return table_full[n, 0] + elif x < xlarge: + return boys_taylor(n, x) + else: + return boys_xlarge_precalc(n, x) + + if isinstance(xs, np.ndarray): + boys_list = list() + + with np.nditer(xs) as it: + for x in it: + boys_list.append(func(n, x)) + boys_table = np.reshape(boys_list, xs.shape) + return boys_table + else: + return func(n, xs) + + return boys + + +# Calculation should take ~ 20 ms +# python -m timeit -n 25 "from sympleints.boys import boys" +boys = get_boys_func(nmax=8) diff --git a/sympleints/cart2sph.py b/sympleints/cart2sph.py new file mode 100755 index 0000000..b25d714 --- /dev/null +++ b/sympleints/cart2sph.py @@ -0,0 +1,270 @@ +#!/usr/bin/env python + +# [1] https://doi.org/10.1002/qua.560540202 +# Transformation between Cartesian and pure spherical harmonic Gaussians +# Schlegel, Frisch, 1995 + +import argparse +import functools +from cmath import sqrt +from math import floor +import sys +from typing import Dict, List, Tuple + +import mpmath +from mpmath import factorial as mp_fact, sqrt as mp_sqrt, binomial as mp_binom +import numpy as np +from scipy.special import factorial as fact +from scipy.special import binom + +from pysisyphus.wavefunction.helpers import canonical_order +from pysisyphus.wavefunction.normalization import get_lmn_factors + + +ZERO_THRESH = 1e-14 + + +@functools.cache +def cart2sph_coeff(lx: int, ly: int, lz: int, m: int) -> complex: + """Coefficient to convert Cartesian to spherical harmonics. + + Eq. (15) in [1]. + + Paramters + --------- + lx + Cartesian quantum number. + ly + Cartesian quantum number. + lz + Cartesian quantum number. + m + Magnetic quantum number. + + Returns + ------- + coeff + Contribution of the given Cartesian basis function to the spherical + harmonic Y^(lx + ly + lz)_m. + """ + + l = lx + ly + lz + assert -l <= m <= l + + j = (lx + ly - abs(m)) / 2 + if (j % 1) != 0: # Return when j is not integer + return 0 + j = int(j) + + lfact = fact(l) # l! + lmm = l - abs(m) # (l - |m|) + sign = 1 if m >= 0 else -1 + + coeff = 0.0j + for i in range(floor(lmm / 2) + 1): + i_fact = ( + binom(l, i) + * binom(i, j) + * (-1) ** i + * fact(2 * l - 2 * i) + / fact(lmm - 2 * i) + ) + k_sum = 0.0 + for k in range(j + 1): + prefact = (-1) ** (sign * (abs(m) - lx + 2 * k) / 2) + k_sum += prefact * binom(j, k) * binom(abs(m), lx - 2 * k) + coeff += i_fact * k_sum + coeff *= sqrt( + fact(2 * lx) + * fact(2 * ly) + * fact(2 * lz) + * lfact + * fact(lmm) + / (fact(2 * l) * fact(lx) * fact(ly) * fact(lz) * fact(l + abs(m))) + ) + coeff *= 1 / (2**l * lfact) + + # print(f"{lx=}, {ly=}, {lz=}, {m=}, {coeff}") + return coeff + + +@functools.cache +@mpmath.workdps(24) +def mp_cart2sph_coeff(lx: int, ly: int, lz: int, m: int) -> mpmath.mpc: + """Higher precision Cartesian to Spherical transformation coefficients. + + Eq. (15) in [1]. + + Paramters + --------- + lx + Cartesian quantum number. + ly + Cartesian quantum number. + lz + Cartesian quantum number. + m + Magnetic quantum number. + + Returns + ------- + coeff + Contribution of the given Cartesian basis function to the spherical + harmonic Y^(lx + ly + lz)_m. + """ + + l = lx + ly + lz + assert -l <= m <= l + + j = (lx + ly - abs(m)) / 2 + if (j % 1) != 0: # Return when j is not integer + return mpmath.mp.mpc(0) + j = int(j) + + lfact = mp_fact(l) # l! + lmm = l - abs(m) # (l - |m|) + sign = 1 if m >= 0 else -1 + + coeff = mpmath.mp.mpc(0.0j) + for i in range(floor(lmm / 2) + 1): + i_fact = ( + mp_binom(l, i) + * mp_binom(i, j) + * (-1) ** i + * mp_fact(2 * l - 2 * i) + / mp_fact(lmm - 2 * i) + ) + k_sum = 0.0 + for k in range(j + 1): + prefact = (-1) ** (sign * (abs(m) - lx + 2 * k) / 2) + k_sum += prefact * mp_binom(j, k) * mp_binom(abs(m), lx - 2 * k) + coeff += i_fact * k_sum + coeff *= mp_sqrt( + mp_fact(2 * lx) + * mp_fact(2 * ly) + * mp_fact(2 * lz) + * lfact + * mp_fact(lmm) + / ( + mp_fact(2 * l) + * mp_fact(lx) + * mp_fact(ly) + * mp_fact(lz) + * mp_fact(l + abs(m)) + ) + ) + coeff *= 1 / (2**l * lfact) + + # print(f"{lx=}, {ly=}, {lz=}, {m=}, {coeff}") + return coeff + + +@functools.cache +def cart2sph_coeffs_for( + l: int, + real: bool = True, + zero_small: bool = True, + zero_thresh: float = ZERO_THRESH, + use_mp=False, +) -> np.ndarray: + if use_mp: + cart2sph_coeff_func = mp_cart2sph_coeff + sqrt_func = mp_sqrt + # dtype = np.longcomplex + else: + cart2sph_coeff_func = cart2sph_coeff + sqrt_func = sqrt + # dtype = complex + all_coeffs = list() + for lx, ly, lz in canonical_order(l): + coeffs = list() + for m in range(-l, l + 1): + coeff = cart2sph_coeff_func(lx, ly, lz, m) + # Form linear combination for m != 0 to get real spherical orbitals. + if real and m != 0: + coeff_minus = cart2sph_coeff_func(lx, ly, lz, -m) + sign = 1 if m < 0 else -1 + coeff = (coeff + sign * coeff_minus) / sqrt_func(sign * 2) + coeffs.append(coeff) + all_coeffs.append(coeffs) + C = np.array(all_coeffs, dtype=complex) # or use dtype=dtype + if real: + assert real == ~np.iscomplex(C).all() # C should be real + C = np.real(C) + if zero_small: + mask = np.abs(C) <= zero_thresh + C[mask] = 0.0 + # Final shape will be (spherical, Cartesian) + # Return read-only view; otherwise it would be too easy to mess up the + # cached result by inplace operations on C. + C = C.T.view() + C.flags.writeable = False + return C + + +@functools.cache +def expand_sph_quantum_numbers( + Lm, zero_thresh=ZERO_THRESH, with_lmn_factors=False +) -> Tuple[np.ndarray, List[Tuple[int, int, int]]]: + """Factors and Cart. angular momentum vectors for given sph. quantum numbers.""" + L, m = Lm + m += L # Map m from [-L, L] onto [0, 2*L] + C = cart2sph_coeffs_for(L) # shape (spherical, Cartesian) + # Include lmn-factor that appears when contracted functions are normalized + # according to 'normalization.norm_cgto_lmn(L). + if with_lmn_factors: + lmn_factors = get_lmn_factors(L) + C = C * lmn_factors[None, :] + indices, *_ = np.where(np.abs(C[m]) > zero_thresh) + factors = C[m, indices].copy() + + cart_lmn = [lmn for i, lmn in enumerate(canonical_order(L)) if i in indices] + return factors, cart_lmn + + +def cart2sph_coeffs( + l_max: int, + **kwargs, +) -> Dict[int, np.ndarray]: + coeffs = {l: cart2sph_coeffs_for(l, **kwargs) for l in range(l_max + 1)} + return coeffs + + +def cart2sph_nlms(l_max: int) -> Dict[int, Tuple[Tuple[int, int, int]]]: + nlms = dict() + for l in range(l_max + 1): + n = l + nlms[l] = tuple([(n, l, m) for m in range(-l, l + 1)]) + return nlms + + +def print_coeffs(l: int, C: np.ndarray): + print("".join([f"{m: >13d}" for m in range(-l, l + 1)])) + for (lx, ly, lz), line in zip(canonical_order(l), C): + fmt = "{:>12.4f}" * len(line) + verb = "".join(["x"] * lx + ["y"] * ly + ["z"] * lz) + print(f"{verb}: {fmt.format(*line)}") + + +def parse_args(args): + parser = argparse.ArgumentParser() + + parser.add_argument("l_max", type=int, default=3) + + return parser.parse_args(args) + + +def run(): + args = parse_args(sys.argv[1:]) + + l_max = args.l_max + + Cs = cart2sph_coeffs(l_max) + for l, C in Cs.items(): + print(f"{l=}") + print_coeffs(l, C.T) + print() + + +if __name__ == "__main__": + run() diff --git a/sympleints/config.py b/sympleints/config.py index 299fcad..20a2f8a 100644 --- a/sympleints/config.py +++ b/sympleints/config.py @@ -3,7 +3,7 @@ L_MAX = 4 # Up to G by default L_AUX_MAX = 6 # Up to I by default -PREC = 16 +PREC = 18 # Base name for arrays in generated code ARR_BASE_NAME = "work" CACHE_DIR = Path(".sympleints") diff --git a/sympleints/defs/__init__.py b/sympleints/defs/__init__.py index 6b2c003..9a4187f 100644 --- a/sympleints/defs/__init__.py +++ b/sympleints/defs/__init__.py @@ -1,4 +1,5 @@ from sympleints.defs.strategy import RecurStrategy, Strategy -from sympleints.defs.twocenter1d import TwoCenter1d +from sympleints.defs.twocenter1d import TwoCenter1d, TwoCenterFunc1d +from sympleints.defs.twocenter3d import TwoCenter3dFunc from sympleints.defs.fourcenter1d import FourCenter1d from sympleints.defs.multipole import Multipole1d diff --git a/sympleints/defs/multipole.py b/sympleints/defs/multipole.py index a69b38f..b8068dd 100644 --- a/sympleints/defs/multipole.py +++ b/sympleints/defs/multipole.py @@ -108,7 +108,7 @@ def gen_multipole_sph_shell(La_tot, Lb_tot, a, b, A, B, Le_tot=2): exprs = list() all_lmns = list() - # Loop over all operators that we wan't to generate. + # Loop over all operators that we want to generate. for pterms in poly_terms: # Every spherical multipole integral is the sum of multiple basic multipole integrals. # Depending on the angular momenta of the involved basis functions some integrals diff --git a/sympleints/defs/overlap.py b/sympleints/defs/overlap.py index af090d8..b38d4aa 100644 --- a/sympleints/defs/overlap.py +++ b/sympleints/defs/overlap.py @@ -1,6 +1,34 @@ +import functools + +import sympy as sym + +from sympleints import shell_iter from sympleints.defs.multipole import gen_multipole_shell +from sympleints.symbols import center_rP, center_rP2 def gen_overlap_shell(La_tot, Lb_tot, a, b, A, B): + """Overlap integral between two shells.""" exprs, lmns = gen_multipole_shell(La_tot, Lb_tot, a, b, A, B) return exprs, lmns + + +class Prefactor: + rP = sym.Matrix(center_rP).reshape(1, 3) + rP2 = sym.Matrix(center_rP2).reshape(1, 3) + + @functools.cache + def __call__(self, i, k, m, j, l, n): + rPx, rPy, rPz = self.rP + rPx2, rPy2, rPz2 = self.rP2 + expr = rPx ** (i + j) * rPy ** (k + l) * rPz ** (m + n) + # Substitute squares in, as these are also available + expr = expr.subs(rPx**2, rPx2).subs(rPy**2, rPy2).subs(rPz**2, rPz2) + return expr + + +def gen_prefactor_shell(La_tot, Lb_tot): + lmns = list(shell_iter((La_tot, Lb_tot))) + pref = Prefactor() + exprs = [pref(*La, *Lb) for La, Lb in lmns] + return exprs, lmns diff --git a/sympleints/defs/twocenter1d.py b/sympleints/defs/twocenter1d.py index 4ffa9da..780d841 100644 --- a/sympleints/defs/twocenter1d.py +++ b/sympleints/defs/twocenter1d.py @@ -49,3 +49,46 @@ def PR(self): @property def K(self): return exp(-self.mu * self.AB**2) + + +class TwoCenterFunc1d: + def __init__(self, ax, A, bx, B, R, r): + self.ax = ax + self.A = A + self.bx = bx + self.B = B + self.r = r + self.R = R + + @property + def p(self): + """Total exponent p.""" + return self.ax + self.bx + + @property + def mu(self): + """Reduced exponent mu.""" + return self.ax * self.bx / self.p + + @property + def AB(self): + """Relative coordinate/Gaussian separation X_AB.""" + return self.A - self.B + + @property + def P(self): + """Center-of-charge coordinate.""" + return (self.ax * self.A + self.bx * self.B) / self.p + + @property + def K(self): + return exp(-self.mu * self.AB**2) + + @property + def Pr(self): + """Relative coordinate/Gaussian separation X_Pr.""" + return self.P - self.r + + @property + def product(self): + return self.K * exp(-self.p * self.Pr**2) diff --git a/sympleints/defs/twocenter3d.py b/sympleints/defs/twocenter3d.py new file mode 100644 index 0000000..59b65dd --- /dev/null +++ b/sympleints/defs/twocenter3d.py @@ -0,0 +1,68 @@ +import sympy as sym + + +class TwoCenter3dFunc: + def __init__(self, ax, A, bx, B, R, r): + self.ax = ax + self.A = A + self.bx = bx + self.B = B + self.R = sym.Matrix(R).reshape(1, 3) + self.r = r + + @property + def p(self): + """Total exponent p.""" + return self.ax + self.bx + + @property + def mu(self): + """Reduced exponent mu.""" + return self.ax * self.bx / self.p + + @property + def AB(self): + """Relative coordinate/Gaussian separation X_AB.""" + return self.A - self.B + + @property + def P(self): + """Center-of-charge coordinate.""" + return (self.ax * self.A + self.bx * self.B) / self.p + + @property + def rA(self): + return self.r - self.A + + @property + def rB(self): + return self.r - self.B + + @property + def rR(self): + return self.r - self.R + + def prefact(self, coords, i, j, k): + return coords[0] ** i * coords[1] ** j * coords[2] ** k + + def prefact_rA(self, i, j, k): + return self.prefact(self.rA, i, j, k) + + def prefact_rB(self, i, j, k): + return self.prefact(self.rB, i, j, k) + + @property + def K(self): + return sym.exp(-self.mu * self.AB.dot(self.AB)) + + @property + def Pr(self): + """Relative coordinate/Gaussian separation X_Pr.""" + return self.P - self.r + + @property + def product(self): + return self.K * sym.exp(-self.p * self.Pr.dot(self.Pr)) + + def full_product(self, i, k, m, j, l, n): + return self.prefact_rA(i, k, m) * self.prefact_rB(j, l, n) * self.product diff --git a/sympleints/main.py b/sympleints/main.py index d7df1d3..00b0b08 100644 --- a/sympleints/main.py +++ b/sympleints/main.py @@ -29,16 +29,15 @@ import argparse from enum import Enum - from datetime import datetime import functools import itertools as it from math import prod -import os from pathlib import Path import sys import textwrap import time +from typing import Optional from jinja2 import Template from sympy import __version__ as sympy_version @@ -64,6 +63,7 @@ shell_iter, get_timer_getter, ) +from sympleints.cart2sph import cart2sph_coeffs, cart2sph_nlms from sympleints.config import L_MAX, L_AUX_MAX, PREC # Integral definitions @@ -82,7 +82,7 @@ gen_multipole_shell, gen_multipole_sph_shell, ) -from sympleints.defs.overlap import gen_overlap_shell +from sympleints.defs.overlap import gen_overlap_shell, gen_prefactor_shell from sympleints.Functions import Functions from sympleints.helpers import L_MAP from sympleints.l_iters import ll_iter, lllaux_iter @@ -91,17 +91,11 @@ from sympleints.FortranRenderer import FortranRenderer from sympleints.NumbaRenderer import NumbaRenderer from sympleints.PythonRenderer import PythonRenderer -from sympleints.symbols import center_R, R, R_map - -try: - from pysisyphus.wavefunction.cart2sph import cart2sph_coeffs, cart2sph_nlms - - can_sph = True -except ModuleNotFoundError: - can_sph = False +from sympleints.symbols import center_R, R, R_map, r, r_map, center_r, rP_map, rP2_map KEYS = ( + "prefactor", "gto", "ovlp", "dpm", @@ -122,10 +116,9 @@ } -if can_sph: - # Pregenerate coefficients - CART2SPH = cart2sph_coeffs(max(L_MAX, L_AUX_MAX) + 2, zero_small=True) - NLMS = cart2sph_nlms(max(L_MAX, L_AUX_MAX) + 2) +# Pregenerate coefficients; now with mpmath enabled +CART2SPH = cart2sph_coeffs(max(L_MAX, L_AUX_MAX) + 2, zero_small=True, use_mp=True) +NLMS = cart2sph_nlms(max(L_MAX, L_AUX_MAX) + 2) def cart2spherical(L_tots, exprs): @@ -377,7 +370,10 @@ def inner(Ls): return integral_gen -def make_header(args): +def make_header(): + argv = sys.argv + args = parse_args(argv[1:]) + cmd = " ".join(sys.argv) tpl = Template( """ @@ -403,80 +399,40 @@ def make_header(args): return header -def parse_args(args): - parser = argparse.ArgumentParser() - - parser.add_argument( - "--lmax", - default=L_MAX, - type=int, - help="Generate 1e-integrals up to this maximum angular momentum.", - ) - parser.add_argument( - "--lauxmax", - default=L_AUX_MAX, - type=int, - help="Maximum angular moment for integrals using auxiliary functions.", - ) - parser.add_argument( - "--write", - action="store_true", - help="Write out generated integrals to the current directory, potentially " - "overwriting the present modules.", - ) - parser.add_argument( - "--out-dir", - default="ints", - help="Directory, where the generated integrals are written.", - ) - keys_str = f"({', '.join(KEYS)})" - parser.add_argument( - "--keys", - nargs="+", - help=f"Generate only certain expressions. Possible keys are: {keys_str}. " - "If not given, all expressions are generated.", - ) - parser.add_argument("--sph", action="store_true") - parser.add_argument( - "--opt-basic", action="store_true", help="Turn on basic optimizations in CSE." - ) - parser.add_argument("--normalize", choices=normalization_map.keys(), default="none") - - return parser.parse_args(args) - - -def run(args): - l_max = args.lmax - l_aux_max = args.lauxmax - sph = args.sph - normalization = normalization_map[args.normalize] - out_dir = Path(args.out_dir if not args.write else ".").absolute() - keys = args.keys - - cse_kwargs = None - if args.opt_basic: - cse_kwargs = { - "optimizations": "basic", - } +def run( + l_max: int, + l_aux_max: int, + sph: bool = False, + normalization=Normalization.NONE, + out_dir=".", + prefix="ints", + keys=None, + cse_kwargs: Optional[dict] = None, + boys_func: Optional[str] = None, + cli=False, +): + out_dir = Path(out_dir).absolute() if keys is None: keys = list() - header = make_header(args) + if cse_kwargs is None: + cse_kwargs = {} - INT_KIND = "Spherical" if sph else "Cartesian" + if cli: + header = make_header() + else: + header = """TODO: fix headers for API call!""" + + int_kind = "Spherical" if sph else "Cartesian" # Cartesian basis function centers A, B, C and D. center_A = get_center("A") center_B = get_center("B") center_C = get_center("C") center_D = get_center("D") - # Multipole origin or nuclear position - # center_R = get_center("R") - # Orbital exponents ax, bx, cx, dx. ax, bx, cx, dx = symbols("ax bx cx dx", real=True) - # Contraction coefficients da, db, dc, dd. contr_coeffs = symbols("da db dc dd", real=True) da, db, dc, dd = contr_coeffs @@ -497,24 +453,64 @@ def run(args): ) renderers = [ PythonRenderer(), - NumbaRenderer(), + # NumbaRenderer(), # FortranRenderer(), ] + results = dict() def render_write(funcs): nonlocal out_dir fns = list() - # Force a trailing slash, as pathlib strips it out - out_dir = f"{out_dir}{os.sep}" for renderer in renderers: lang = renderer.language.lower() - rend_out_dir = Path(str(out_dir) + f"_{lang}") + lang_results = results.setdefault(lang, dict()) + rend_out_dir = out_dir / f"{prefix}_{lang}" if not rend_out_dir.exists(): rend_out_dir.mkdir() - fns.append(renderer.render_write(funcs, rend_out_dir)) + fn = renderer.render_write(funcs, rend_out_dir) + lang_results[funcs.name] = fn + fns.append(fn) return fns + ############# + # Prefactor # + ############ + + def prefactor(): + def doc_func(L_tots): + La_tot, Lb_tot = L_tots + shell_a = L_MAP[La_tot] + shell_b = L_MAP[Lb_tot] + prefix = "Spherical" if sph else "Cartesian" + return f"{prefix} prefactors for 3D ({shell_a}{shell_b}) overlap at distance RP." + + ls_exprs = integral_gen( + lambda La_tot, Lb_tot: gen_prefactor_shell(La_tot, Lb_tot), + (l_max, l_max), + (ax, bx), + "prefactor", + (rP_map, rP2_map), + L_iter=ll_iter(l_max), + ) + + prefactor_funcs = Functions( + name="prefactors", + l_max=l_max, + coeffs=[da, db], + exponents=[ax, bx], + centers=[A, B], + ls_exprs=ls_exprs, + doc_func=doc_func, + header=header, + spherical=sph, + primitive=True, + hermitian=[0, 1], + parallel=False, + with_ref_center=False, + ) + render_write(prefactor_funcs) + ################# # Cartesian GTO # ################# @@ -524,7 +520,7 @@ def doc_func(L_tot): (La_tot,) = L_tot shell_a = L_MAP[La_tot] return ( - f"3D {INT_KIND} {shell_a}-Gaussian shell.\n" + f"3D {int_kind} {shell_a}-Gaussian shell.\n" "Exponent ax, contraction coeff. da, centered at A, evaluated at R." ) @@ -559,7 +555,7 @@ def doc_func(L_tots): La_tot, Lb_tot = L_tots shell_a = L_MAP[La_tot] shell_b = L_MAP[Lb_tot] - return f"{INT_KIND} 3D ({shell_a}{shell_b}) overlap integral." + return f"{int_kind} 3D ({shell_a}{shell_b}) overlap integral." ls_exprs = integral_gen( lambda La_tot, Lb_tot: gen_overlap_shell(La_tot, Lb_tot, ax, bx, A, B), @@ -594,7 +590,7 @@ def doc_func(L_tots): shell_a = L_MAP[La_tot] shell_b = L_MAP[Lb_tot] return ( - f"{INT_KIND} 3D ({shell_a}{shell_b}) dipole moment integrals.\n" + f"{int_kind} 3D ({shell_a}{shell_b}) dipole moment integrals.\n" "The origin is at R." ) @@ -649,7 +645,7 @@ def doc_func(L_tots): shell_a = L_MAP[La_tot] shell_b = L_MAP[Lb_tot] return ( - f"{INT_KIND} 3D ({shell_a}{shell_b}) quadrupole moment integrals\n" + f"{int_kind} 3D ({shell_a}{shell_b}) quadrupole moment integrals\n" "for operators x², y² and z². The origin is at R." ) @@ -699,7 +695,7 @@ def doc_func(L_tots): shell_a = L_MAP[La_tot] shell_b = L_MAP[Lb_tot] return ( - f"{INT_KIND} 3D ({shell_a}{shell_b}) quadrupole moment integrals.\n" + f"{int_kind} 3D ({shell_a}{shell_b}) quadrupole moment integrals.\n" "The origin is at R." ) @@ -757,7 +753,7 @@ def doc_func(L_tots): shell_a = L_MAP[La_tot] shell_b = L_MAP[Lb_tot] return ( - f"Primitive {INT_KIND} 3D ({shell_a}{shell_b}) spherical multipole integrals.\n" + f"Primitive {int_kind} 3D ({shell_a}{shell_b}) spherical multipole integrals.\n" "In contrast to the other multipole integrals, the origin R is calculated\n" "inside the function and is (possibly) unique for all primitive pairs." ) @@ -809,7 +805,7 @@ def doc_func(L_tots): La_tot, Lb_tot = L_tots shell_a = L_MAP[La_tot] shell_b = L_MAP[Lb_tot] - return f"{INT_KIND} 3D ({shell_a}{shell_b}) kinetic energy integral." + return f"{int_kind} 3D ({shell_a}{shell_b}) kinetic energy integral." ls_exprs = integral_gen( lambda La_tot, Lb_tot: gen_kinetic_shell( @@ -846,7 +842,7 @@ def doc_func(L_tots): La_tot, Lb_tot = L_tots shell_a = L_MAP[La_tot] shell_b = L_MAP[Lb_tot] - return f"{INT_KIND} ({shell_a}{shell_b}) 1-electron Coulomb integral." + return f"{int_kind} ({shell_a}{shell_b}) 1-electron Coulomb integral." ls_exprs = integral_gen( lambda La_tot, Lb_tot: CoulombShell( @@ -866,7 +862,7 @@ def doc_func(L_tots): centers=[A, B], ls_exprs=ls_exprs, ncomponents=1, - boys=True, + boys_func=boys_func, doc_func=doc_func, header=header, spherical=sph, @@ -884,7 +880,7 @@ def doc_func(L_tots): shell_a = L_MAP[La_tot] shell_b = L_MAP[Lb_tot] return ( - f"{INT_KIND} ({shell_a}|{shell_b}) " + f"{int_kind} ({shell_a}|{shell_b}) " "two-center two-electron repulsion integral." ) @@ -911,7 +907,7 @@ def doc_func(L_tots): centers=[A, B], ls_exprs=ls_exprs, ncomponents=1, - boys=True, + boys_func=boys_func, doc_func=doc_func, header=header, spherical=sph, @@ -962,10 +958,10 @@ def doc_func(L_tots): shell_b = L_MAP[Lb_tot] shell_c = L_MAP[Lc_tot] doc_str = ( - f"{INT_KIND} ({shell_a}{shell_b}|{shell_c}) three-center " + f"{int_kind} ({shell_a}{shell_b}|{shell_c}) three-center " "two-electron repulsion integral." ) - if INT_KIND == "Cartesian": + if int_kind == "Cartesian": doc_str += ( "\nThese integrals MUST BE converted to spherical harmonics!\n" "\nIntegral generation utilized Ahlrichs (truncated) vertical " @@ -991,7 +987,7 @@ def doc_func(L_tots): exponents=[ax, bx, cx], centers=[A, B, C], ls_exprs=ls_exprs, - boys=True, + boys_func=boys_func, doc_func=doc_func, header=header, spherical=sph, @@ -1049,6 +1045,7 @@ def doc_func(L_tots): """ funcs = { + "prefactor": prefactor, "gto": gto, # Cartesian Gaussian-type-orbital for density evaluation "ovlp": overlap, # Overlap integrals "dpm": dipole, # Linear moment (dipole) integrals @@ -1062,7 +1059,9 @@ def doc_func(L_tots): "3c2e_sph": _3center2electron_sph, # Sph. 3-center-2-electron DF integrals # "4covlp": fourcenter_overlap, # Four center overlap integral } - assert set(funcs.keys()) == set(KEYS) + assert set(funcs.keys()) == set( + KEYS + ), "Did you just add a key to 'keys' but forgot to also add it to 'KEYS'?" # Generate all possible integrals, when no 'keys' were supplied negate_keys = list() @@ -1081,12 +1080,91 @@ def doc_func(L_tots): duration_hms = str(duration).split(".")[0] # Only keep hh:mm:ss print(f"sympleint run took {duration_hms} h.") - return 0 + return results + + +def parse_args(args): + parser = argparse.ArgumentParser() + + parser.add_argument( + "--lmax", + default=L_MAX, + type=int, + help="Generate 1e-integrals up to this maximum angular momentum.", + ) + parser.add_argument( + "--lauxmax", + default=L_AUX_MAX, + type=int, + help="Maximum angular moment for integrals using auxiliary functions.", + ) + """ + parser.add_argument( + "--write", + action="store_true", + help="Write out generated integrals to the current directory, potentially " + "overwriting the present modules.", + ) + """ + parser.add_argument( + "--out-dir", + default=".", + help="Directory, where the generated integrals are written.", + ) + parser.add_argument( + "--prefix", default="ints", help="Prefix for the integral directories." + ) + keys_str = f"({', '.join(KEYS)})" + parser.add_argument( + "--keys", + nargs="+", + help=f"Generate only certain expressions. Possible keys are: {keys_str}. " + "If not given, all expressions are generated.", + ) + parser.add_argument("--sph", action="store_true") + parser.add_argument( + "--opt-basic", action="store_true", help="Turn on basic optimizations in CSE." + ) + parser.add_argument( + "--boys-func", + default="pysisyphus.wavefunction.ints.boys", + help="Which Boys-function to use.", + ) + parser.add_argument("--normalize", choices=normalization_map.keys(), default="none") + + return parser.parse_args(args) def run_cli(): args = parse_args(sys.argv[1:]) - return run(args) + + l_max = args.lmax + l_aux_max = args.lauxmax + sph = args.sph + normalization = normalization_map[args.normalize] + out_dir = args.out_dir + prefix = args.prefix + keys = args.keys + boys_func = args.boys_func + + cse_kwargs = None + if args.opt_basic: + cse_kwargs = { + "optimizations": "basic", + } + + return run( + l_max=l_max, + l_aux_max=l_aux_max, + sph=sph, + normalization=normalization, + out_dir=out_dir, + prefix=prefix, + keys=keys, + cse_kwargs=cse_kwargs, + boys_func=boys_func, + cli=True, + ) if __name__ == "__main__": diff --git a/sympleints/symbols.py b/sympleints/symbols.py index 99d7f91..7c0becd 100644 --- a/sympleints/symbols.py +++ b/sympleints/symbols.py @@ -4,3 +4,19 @@ # Reference center/multipole origin center_R = get_center("R") R, R_map = get_map("R", center_R) + +# Coordinates at which a function is evaluated +center_r = get_center("r") +r, r_map = get_map("r", center_r) + +# Center-of-charge coordinate +center_P = get_center("P") +P, P_map = get_map("P", center_P) + + +center_rP = get_center("rP") +rP, rP_map = get_map("rP", center_rP) + + +center_rP2 = get_center("rP2") +rP2, rP2_map = get_map("rP2", center_rP2) diff --git a/sympleints/templates/py_equi_func.tpl b/sympleints/templates/py_equi_func.tpl index 3aeeb90..a4361d9 100644 --- a/sympleints/templates/py_equi_func.tpl +++ b/sympleints/templates/py_equi_func.tpl @@ -1,6 +1,7 @@ -def {{ equi_name }}({{ args }}): +def {{ equi_name }}({{ args }}, result): """See docstring of {{ name }}.""" + # Calculate values w/ swapped arguments + {{ name }}({{ equi_args }}, result) # Swap two axes - result = numpy.moveaxis({{ name }}({{ equi_args }}), {{ from_axes }}, {{ to_axes }}) - return result + result[:] = numpy.moveaxis(result.reshape({{ reshape|join(",") }}), {{ from_axes }}, {{ to_axes }}).flatten() diff --git a/sympleints/templates/py_func.tpl b/sympleints/templates/py_func.tpl index cc3349c..2594b2b 100644 --- a/sympleints/templates/py_func.tpl +++ b/sympleints/templates/py_func.tpl @@ -1,16 +1,16 @@ -def {{ name }}({{ args }}): +def {{ name }}({{ args }}, result): {% if doc_str %} """{{ doc_str }}""" {% endif %} - result = numpy.zeros({{ shape }}, dtype=float) + {# result = numpy.zeros({{ shape }}, dtype=float) #} {% for line in py_lines %} {{ line }} {% endfor %} # {{ n_return_vals }} item(s) - {% for inds, res_line in results_iter %} - result[{{ inds|join(", ")}}] = {{ res_line }} + # TODO: check if removing the += gives problems w/ primitive=True + {% for _, res_line in results_iter %} + result[{{ loop.index0 }}] = {{ res_line }} {% endfor %} - return result diff --git a/sympleints/templates/py_module.tpl b/sympleints/templates/py_module.tpl index ee85a89..4f899fb 100644 --- a/sympleints/templates/py_module.tpl +++ b/sympleints/templates/py_module.tpl @@ -11,7 +11,9 @@ import numpy {% if boys %} -from pysisyphus.wavefunction.ints.boys import boys +# The default boys-import can be changed via the --boys-func argument. +# sympleints also includes a basic implementation in 'sympleints.boys' +from {{ boys }} import boys {% endif %} {% for ai in add_imports %} diff --git a/sympleints/testing/__init__.py b/sympleints/testing/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/sympleints/testing/intor.py b/sympleints/testing/intor.py new file mode 100644 index 0000000..da5d15b --- /dev/null +++ b/sympleints/testing/intor.py @@ -0,0 +1,76 @@ +import numpy as np + + +def get_size_func(spherical): + if spherical: + + def size_func(L): + return 2 * L + 1 + + else: + + def size_func(L): + return (L + 2) * (L + 1) // 2 + + return size_func + + +def intor_2c( + shells, + func_dict, + ncomponents=1, + spherical=False, + **kwargs, +): + """Basic driver for integrals over hermitian operators.""" + if spherical: + nbfs = sum([shell.sph_size for shell in shells]) + else: + nbfs = sum([shell.cart_size for shell in shells]) + + L_max = max([shell.L for shell in shells]) + size_func = get_size_func(spherical) + size_max = size_func(L_max) + + # Preallocate empty matrices and directly assign the calculated values + integrals = np.zeros((ncomponents, nbfs, nbfs)) + + a_start = 0 + result = np.empty(ncomponents * size_max**2) + for i, shell_a in enumerate(shells): + L_a = shell_a.L + coeffs_a = shell_a.coeffs + exps_a = shell_a.exps + center_a = shell_a.center + a_size = size_func(L_a) + a_slice = slice(a_start, a_start + a_size) + b_start = a_start + for shell_b in shells[i:]: + b_size = size_func(shell_b.L) + batch_size = ncomponents * a_size * b_size + b_slice = slice(b_start, b_start + b_size) + func_dict[(L_a, shell_b.L)]( + exps_a[:, None], + coeffs_a[:, None], + center_a, + shell_b.exps[None, :], + shell_b.coeffs[None, :], + shell_b.center, + result=result[:batch_size], + **kwargs, + ) + integrals[:, a_slice, b_slice] = result[:batch_size].reshape( + ncomponents, a_size, b_size + ) + # Take symmetry into account + if a_start != b_start: + integrals[:, b_slice, a_slice] = np.transpose( + integrals[:, a_slice, b_slice], axes=(0, 2, 1) + ) + b_start += b_size + a_start += a_size + + # Return plain 2d array if components is set to 0, i.e., remove first axis. + if ncomponents == 1: + integrals = np.squeeze(integrals, axis=0) + return integrals diff --git a/sympleints/testing/normalization.py b/sympleints/testing/normalization.py new file mode 100644 index 0000000..8ea3349 --- /dev/null +++ b/sympleints/testing/normalization.py @@ -0,0 +1,43 @@ +from math import prod + +import numpy as np +from scipy.special import factorial2 as sp_factorial2 + +from sympleints.helpers import canonical_order + + +def factorial2(n: int) -> int: + """Scipy 1.11 decided that (-1)!! is not 1 anymore! + + Please see https://github.com/scipy/scipy/issues/18813.""" + if n == -1: + return 1 + elif n < -1: + raise Exception(f"Only supported negative argument is -1, but got {n}!") + return sp_factorial2(n) + + +def get_lmn_factors(L: int) -> np.ndarray: + lmns = canonical_order(L) + lmn_factors = np.zeros(len(lmns)) + for i, lmn in enumerate(lmns): + lmn_factors[i] = prod([factorial2(2 * am - 1) for am in lmn]) + lmn_factors = 1 / np.sqrt(lmn_factors) + return lmn_factors + + +def norm_cgto_lmn( + coeffs: np.ndarray, exps: np.ndarray, L: int +) -> tuple[np.ndarray, np.ndarray]: + N = 0.0 + for i, expi in enumerate(exps): + for j, expj in enumerate(exps): + tmp = coeffs[i] * coeffs[j] / (expi + expj) ** (L + 1.5) + tmp *= np.sqrt(expi * expj) ** (L + 1.5) + N += tmp + N = np.sqrt(exps ** (L + 1.5) / (np.pi**1.5 / 2**L * N)) + + mod_coeffs = N * coeffs + lmn_factors = get_lmn_factors(L) + + return mod_coeffs, lmn_factors diff --git a/sympleints/testing/pyscf_interface.py b/sympleints/testing/pyscf_interface.py new file mode 100644 index 0000000..b819fe0 --- /dev/null +++ b/sympleints/testing/pyscf_interface.py @@ -0,0 +1,101 @@ +import numpy as np + +from pyscf import gto +import scipy as sp + +from sympleints.testing.shells import Shell + + +def get_mol(atoms: str, gbs_basis: dict) -> gto.Mole: + mol = gto.Mole() + mol.atom = atoms + basis = {atom: gto.basis.parse(atom_bas) for atom, atom_bas in gbs_basis.items()} + mol.basis = basis + mol.build() + return mol + + +_L_MAP = { + 0: "S", + 1: "P", + 2: "D", + 3: "F", + 4: "G", + 5: "H", + 6: "I", +} + + +def get_he_chain(length=2, L_max=4): + atoms = "; ".join([f"He 0.0 0.0 {i:.1f}" for i in range(length)]) + cgtos = list() + for L in range(L_max + 1): + L_char = _L_MAP[L] + cgtos.append( + f"""He {L_char} + 1.0 1.0 + 2.0 1.0""" + ) + gbs_basis = { + "He": "\n".join(cgtos), + } + return get_mol(atoms, gbs_basis) + + +def get_cart_norms(mol: gto.Mole) -> np.ndarray: + S = mol.intor("int1e_ovlp_cart") + N = 1 / np.diag(S) ** 0.5 + NN = N[:, None] * N[None, :] + return NN + + +def shells_from_pyscf_mol(mol): + shells = list() + for bas_id in range(mol.nbas): + L = mol.bas_angular(bas_id) + center = mol.bas_coord(bas_id) + coeffs = mol.bas_ctr_coeff(bas_id).flatten() + exps = mol.bas_exp(bas_id) + assert coeffs.size == exps.size, "General contractions are not yet supported." + center_ind = mol.bas_atom(bas_id) + atomic_num = mol.atom_charge(center_ind) + shell = Shell(L, center, coeffs, exps, center_ind, atomic_num) + shells.append(shell) + return shells + + +PYSCF_SPH_PS = { + 0: [[1]], # s + 1: [[1, 0, 0], [0, 0, 1], [0, 1, 0]], # px py pz + 2: [ + [0, 0, 0, 0, 1], # dxy + [0, 0, 0, 1, 0], # dyz + [0, 0, 1, 0, 0], # dz² + [0, 1, 0, 0, 0], # dxz + [1, 0, 0, 0, 0], # dx² - y² + ], + 3: [ + [0, 0, 0, 0, 0, 0, 1], + [0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 1, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0], + [1, 0, 0, 0, 0, 0, 0], + ], + 4: [ + [0, 0, 0, 0, 0, 0, 0, 0, 1], + [0, 0, 0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 0, 0, 0, 0, 0, 0], + [1, 0, 0, 0, 0, 0, 0, 0, 0], + ], +} + + +def get_sph_permutation_matrix(Ls): + return sp.linalg.block_diag(*[PYSCF_SPH_PS[L] for L in Ls]) diff --git a/sympleints/testing/shells.py b/sympleints/testing/shells.py new file mode 100644 index 0000000..5bb68aa --- /dev/null +++ b/sympleints/testing/shells.py @@ -0,0 +1,26 @@ +from dataclasses import dataclass + +import numpy as np + +from sympleints.testing.normalization import norm_cgto_lmn + + +@dataclass +class Shell: + L: int + center: np.ndarray + coeffs: np.ndarray + exps: np.ndarray + center_ind: int + atomic_num: int + + def __post_init__(self): + self.coeffs, _ = norm_cgto_lmn(self.coeffs, self.exps, self.L) + + @property + def cart_size(self) -> int: + return (self.L + 2) * (self.L + 1) // 2 + + @property + def sph_size(self) -> int: + return 2 * self.L + 1 diff --git a/tests/test_boys/test_boys.py b/tests/test_boys/test_boys.py new file mode 100644 index 0000000..2e05f4b --- /dev/null +++ b/tests/test_boys/test_boys.py @@ -0,0 +1,39 @@ +import time + +import numpy as np + +from sympleints.boys import boys_quad, get_boys_func + + +def test_boys(atol=2e-10): + nmax = 10 + boys = get_boys_func(nmax, order=6) + + ntaylor = 5000 + nlarge = 10 + tot_dur = 0.0 + quad_dur = 0.0 + for n in range(nmax + 1): + xs = np.concatenate( + ( + [0.0], + np.random.rand(ntaylor) * 25, + np.full(nlarge, 30.0) + 5 * np.random.rand(nlarge), + ) + ) + for x in xs: + dur = time.time() + bc = boys(n, x) + tot_dur += time.time() - dur + dur = time.time() + br = boys_quad(n, x) + quad_dur += time.time() - dur + assert abs(bc - br) <= atol, f"{n=}, {x=:.4f}, {bc=:.8e}, {br=:.8e}" + # print(f"Checked {n=} at {x=:.8e}") + # print() + nevals = (nlarge + 1) * ntaylor + per_eval = tot_dur / nevals + per_quad_eval = quad_dur / nevals + print( + f"Eval duration: {tot_dur:.4f} s, {per_eval:.4e} s/eval, {per_quad_eval:.4e} s/quad eval" + ) diff --git a/tests/test_intor/test_intor.py b/tests/test_intor/test_intor.py new file mode 100644 index 0000000..00c16a2 --- /dev/null +++ b/tests/test_intor/test_intor.py @@ -0,0 +1,79 @@ +import importlib + +import numpy as np +import pytest +import tempfile + +import sympleints.testing.intor as intor +import sympleints.testing.pyscf_interface as pi +import sympleints.main + + +def gen_and_compare_to_pyscf(key, name, ncomponents, pyscf_name, spherical, L_max=2): + with tempfile.TemporaryDirectory(prefix="sympleints_test_") as out_dir: + results = sympleints.main.run( + l_max=L_max, + l_aux_max=L_max, + sph=spherical, + keys=(key,), + out_dir=out_dir, + normalization=sympleints.main.Normalization.CGTO, + boys_func="sympleints.boys", + ) + # Get path to generated integrals + int_fn = results["python"][name] + spec = importlib.util.spec_from_file_location(name, int_fn) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + func_dict = getattr(module, name) + + # Compare against Helium-dimer w/ basis functions up to L_max + mol = pi.get_he_chain(2, L_max=L_max) + origin = np.full(3, 0.5) + + shells = pi.shells_from_pyscf_mol(mol) + if key == "coul": + integrals = list() + for R, Z in zip(mol.atom_coords(), mol.atom_charges()): + integrals.append( + -Z * intor.intor_2c(shells, func_dict, ncomponents, spherical, R=R) + ) + integrals = np.sum(integrals, axis=0) + else: + integrals = intor.intor_2c( + shells, func_dict, ncomponents, spherical, R=origin + ) + + # Compare against (renormalized) PySCF integrals + with mol.with_common_orig(origin): + integrals_ref = mol.intor(pyscf_name) + # Bring our integrals in proper order for PySCF + if spherical: + Ls = [shell.L for shell in shells] + P_sph = pi.get_sph_permutation_matrix(Ls) + integrals = np.einsum( + "ij,...jk,kl->...il", P_sph, integrals, P_sph.T, optimize="greedy" + ) + # Ensure self-overlaps of 1.0 for Cartesian integrals + else: + integrals_ref *= pi.get_cart_norms(mol) + + np.testing.assert_allclose(integrals, integrals_ref, atol=1e-14) + + +@pytest.mark.parametrize("spherical", (False, True)) +@pytest.mark.parametrize( + "key, ncomponents, name, pyscf_name", + ( + ("ovlp", 1, "ovlp3d", "int1e_ovlp_cart"), + ("kin", 1, "kinetic3d", "int1e_kin_cart"), + ("dpm", 3, "dipole3d", "int1e_r_cart"), + # ("qpm", 6, "quadrupole3d", "int1e_rr_cart"), + ("coul", 1, "coulomb3d", "int1e_nuc_cart"), + ("2c2e", 1, "int2c2e3d", "int2c2e_cart"), + ), +) +def test_integrals(key, name, ncomponents, pyscf_name, spherical): + if spherical: + pyscf_name = pyscf_name[:-4] + "sph" + gen_and_compare_to_pyscf(key, name, ncomponents, pyscf_name, spherical, L_max=2)