diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index a72772ec..0cdd96e1 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -69,7 +69,7 @@ def evaluate_density( points, transform=None, threshold=1.0e-8, - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the density of the given basis set at the given points. @@ -104,8 +104,7 @@ def evaluate_density( The absolute value below which negative density values are acceptable. Any negative density value with an absolute value smaller than this threshold will be set to zero. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -137,7 +136,7 @@ def evaluate_deriv_reduced_density_matrix( points, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the derivative of the first-order reduced density matrix at the given points. @@ -192,8 +191,7 @@ def evaluate_deriv_reduced_density_matrix( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -241,7 +239,7 @@ def evaluate_deriv_density( points, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the derivative of density of the given transformed basis set at the given points. @@ -291,8 +289,7 @@ def evaluate_deriv_density( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -357,7 +354,7 @@ def evaluate_density_gradient( points, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the gradient of the density evaluated at the given points. @@ -397,8 +394,7 @@ def evaluate_density_gradient( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -448,7 +444,7 @@ def evaluate_density_laplacian( points, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the Laplacian of the density evaluated at the given points. @@ -486,8 +482,7 @@ def evaluate_density_laplacian( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -565,7 +560,7 @@ def evaluate_density_hessian( points, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the Hessian of the density evaluated at the given points. @@ -611,8 +606,7 @@ def evaluate_density_hessian( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -719,7 +713,7 @@ def evaluate_posdef_kinetic_energy_density( transform=None, deriv_type="general", threshold=1.0e-8, - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return evaluations of positive definite kinetic energy density at the given points. @@ -765,8 +759,7 @@ def evaluate_posdef_kinetic_energy_density( The absolute value below which negative density values are acceptable. Any negative density value with an absolute value smaller than this threshold will be set to zero. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -808,7 +801,7 @@ def evaluate_general_kinetic_energy_density( alpha, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return evaluations of general form of the kinetic energy density at the given points. @@ -846,8 +839,7 @@ def evaluate_general_kinetic_energy_density( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that diff --git a/gbasis/evals/electrostatic_potential.py b/gbasis/evals/electrostatic_potential.py index 4d850975..6cdd4209 100644 --- a/gbasis/evals/electrostatic_potential.py +++ b/gbasis/evals/electrostatic_potential.py @@ -136,6 +136,7 @@ def electrostatic_potential( raise TypeError( "`coord_type` must be a list/tuple of the strings 'spherical' or 'cartesian'." ) + # compute Hartree potential (no screening since it is cumulative for all points) hartree_potential = point_charge_integral( basis, points, -np.ones(points.shape[0]), transform=transform ) diff --git a/gbasis/evals/eval.py b/gbasis/evals/eval.py index de06e235..6c2fa110 100644 --- a/gbasis/evals/eval.py +++ b/gbasis/evals/eval.py @@ -55,7 +55,7 @@ class Eval(BaseOneIndex): """ @staticmethod - def construct_array_contraction(contractions, points, screen_basis=True, tol_screen=1e-8): + def construct_array_contraction(contractions, points, screen_basis=False, tol_screen=1e-8): r"""Return the evaluations of the given contractions at the given coordinates. Parameters @@ -69,8 +69,7 @@ def construct_array_contraction(contractions, points, screen_basis=True, tol_scr Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` components. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -145,7 +144,7 @@ def construct_array_contraction(contractions, points, screen_basis=True, tol_scr return output -def evaluate_basis(basis, points, transform=None, screen_basis=True, tol_screen=1e-8): +def evaluate_basis(basis, points, transform=None, screen_basis=False, tol_screen=1e-8): r"""Evaluate the basis set in the given coordinate system at the given points. Parameters @@ -164,8 +163,7 @@ def evaluate_basis(basis, points, transform=None, screen_basis=True, tol_screen= and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that diff --git a/gbasis/evals/eval_deriv.py b/gbasis/evals/eval_deriv.py index 5bd8edcd..d524acfa 100644 --- a/gbasis/evals/eval_deriv.py +++ b/gbasis/evals/eval_deriv.py @@ -58,7 +58,7 @@ class EvalDeriv(BaseOneIndex): @staticmethod def construct_array_contraction( - contractions, points, orders, deriv_type="general", screen_basis=True, tol_screen=1e-8 + contractions, points, orders, deriv_type="general", screen_basis=False, tol_screen=1e-8 ): r"""Return the array associated with a set of contracted Cartesian Gaussians. @@ -82,7 +82,7 @@ def construct_array_contraction( contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional Whether to screen out points that are too far from the contraction center. Default value - is True (enable screening). + is False. tol_screen : float Screening tolerance for excluding points. This value, together with the minimum contraction parameters, determines a cutoff distance. Points @@ -177,7 +177,7 @@ def evaluate_deriv_basis( orders, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Evaluate the derivative of the basis set in the given coordinate system at the given points. @@ -220,7 +220,7 @@ def evaluate_deriv_basis( derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional Whether to screen out points that are too far from the contractions center to contribute to - the result above the given tolerance. Default value is True (enable screening). + the result above the given tolerance. Default value is False. tol_screen : float Screening tolerance for excluding points. This value, together with the minimum contraction parameters, determines a cutoff distance. Points diff --git a/gbasis/evals/stress_tensor.py b/gbasis/evals/stress_tensor.py index 6a33ac51..76994c07 100644 --- a/gbasis/evals/stress_tensor.py +++ b/gbasis/evals/stress_tensor.py @@ -1,4 +1,5 @@ """Module for computing properties related to the stress tensor.""" + from gbasis.evals.density import ( evaluate_density_laplacian, evaluate_deriv_density, @@ -8,7 +9,16 @@ # TODO: need to be tested against reference -def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, transform=None): +def evaluate_stress_tensor( + one_density_matrix, + basis, + points, + alpha=1, + beta=0, + transform=None, + screen_basis=False, + tol_screen=1e-8, +): r"""Return the stress tensor evaluated at the given coordinates. Stress tensor is defined here as: @@ -70,6 +80,13 @@ def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, t beta : {int, float} Second parameter of the stress tensor. Default value is 0. + screen_basis : bool, optional + Whether to screen out points with negligible contributions. Default value is False. + tol_screen : float + Screening tolerance for excluding evaluations. Points with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -99,6 +116,8 @@ def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, t basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if alpha != 1: output[i, j] += (1 - alpha) * evaluate_deriv_reduced_density_matrix( @@ -108,6 +127,8 @@ def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, t basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if i == j and beta != 0: output[i, j] -= ( @@ -118,6 +139,8 @@ def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, t basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) ) output[j, i] = output[i, j] @@ -125,7 +148,16 @@ def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, t # TODO: need to be tested against reference -def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, transform=None): +def evaluate_ehrenfest_force( + one_density_matrix, + basis, + points, + alpha=1, + beta=0, + transform=None, + screen_basis=False, + tol_screen=1e-8, +): r"""Return the Ehrenfest force. Ehrenfest force is the negative of the divergence of the stress tensor: @@ -183,6 +215,13 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, beta : {int, float} Second parameter of the stress tensor. Default value is 0. + screen_basis : bool, optional + Whether to screen out points with negligible contributions. Default value is False. + tol_screen : float + Screening tolerance for excluding evaluations. Points with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -211,6 +250,8 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if alpha != 1: output[i] -= (1 - alpha) * evaluate_deriv_reduced_density_matrix( @@ -220,6 +261,8 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if alpha != 0.5: output[i] -= (1 - 2 * alpha) * evaluate_deriv_reduced_density_matrix( @@ -229,6 +272,8 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if beta != 0: output[i] += ( @@ -240,6 +285,8 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) ) return output.T @@ -254,6 +301,8 @@ def evaluate_ehrenfest_hessian( beta=0, transform=None, symmetric=False, + screen_basis=False, + tol_screen=1e-8, ): r"""Return the Ehrenfest Hessian. @@ -326,6 +375,13 @@ def evaluate_ehrenfest_hessian( Flag for symmetrizing the Hessian. If True, then the Hessian is symmetrized by averaging it with its transpose. Default value is False. + screen_basis : bool, optional + Whether to screen out points with negligible contributions. Default value is False. + tol_screen : float + Screening tolerance for excluding evaluations. Points with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -356,6 +412,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) output[i, j] += alpha * evaluate_deriv_reduced_density_matrix( 2 * orders_one, @@ -364,6 +422,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if alpha != 1: output[i, j] -= (1 - alpha) * evaluate_deriv_reduced_density_matrix( @@ -373,6 +433,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) output[i, j] -= (1 - alpha) * evaluate_deriv_reduced_density_matrix( 2 * orders_one + orders_two, @@ -381,6 +443,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if alpha != 0.5: output[i, j] -= (1 - 2 * alpha) * evaluate_deriv_reduced_density_matrix( @@ -390,6 +454,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) output[i, j] -= (1 - 2 * alpha) * evaluate_deriv_reduced_density_matrix( orders_one + orders_two, @@ -398,6 +464,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if beta != 0: output[i, j] += ( @@ -409,6 +477,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) ) if symmetric: diff --git a/gbasis/integrals/angular_momentum.py b/gbasis/integrals/angular_momentum.py index 9f6503bc..ff01dc72 100644 --- a/gbasis/integrals/angular_momentum.py +++ b/gbasis/integrals/angular_momentum.py @@ -60,7 +60,7 @@ class AngularMomentumIntegral(BaseTwoIndexSymmetric): @staticmethod def construct_array_contraction( - contractions_one, contractions_two, screen_basis=True, tol_screen=1e-8 + contractions_one, contractions_two, screen_basis=False, tol_screen=1e-8 ): """Return the integrals over the angular momentum operator of the given contractions. @@ -73,7 +73,7 @@ def construct_array_contraction( Contracted Cartesian Gaussians (of the same shell) associated with the second index of the angular momentum integral. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance @@ -183,7 +183,7 @@ def construct_array_contraction( return -1j * np.transpose(output, (3, 2, 4, 1, 0)) -def angular_momentum_integral(basis, transform=None, screen_basis=True, tol_screen=1e-8): +def angular_momentum_integral(basis, transform=None, screen_basis=False, tol_screen=1e-8): r"""Return the integral over :math:`hat{L}` of the given basis set. .. math:: @@ -218,7 +218,7 @@ def angular_momentum_integral(basis, transform=None, screen_basis=True, tol_scre and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/kinetic_energy.py b/gbasis/integrals/kinetic_energy.py index 3276abee..eb016633 100644 --- a/gbasis/integrals/kinetic_energy.py +++ b/gbasis/integrals/kinetic_energy.py @@ -55,7 +55,7 @@ class KineticEnergyIntegral(BaseTwoIndexSymmetric): @staticmethod def construct_array_contraction( - contractions_one, contractions_two, screen_basis=True, tol_screen=1e-8 + contractions_one, contractions_two, screen_basis=False, tol_screen=1e-8 ): """Return the evaluations of the given contractions at the given coordinates. @@ -68,7 +68,7 @@ def construct_array_contraction( Contracted Cartesian Gaussians (of the same shell) associated with the second index of the kinetic energy integral. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance @@ -146,7 +146,7 @@ def construct_array_contraction( return -0.5 * np.sum(output, axis=0) -def kinetic_energy_integral(basis, transform=None, screen_basis=True, tol_screen=1e-8): +def kinetic_energy_integral(basis, transform=None, screen_basis=False, tol_screen=1e-8): r"""Return kinetic energy integral of the given basis set. .. math:: @@ -174,7 +174,7 @@ def kinetic_energy_integral(basis, transform=None, screen_basis=True, tol_screen and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/moment.py b/gbasis/integrals/moment.py index 682c95d3..e2f9a6d3 100644 --- a/gbasis/integrals/moment.py +++ b/gbasis/integrals/moment.py @@ -56,7 +56,7 @@ def construct_array_contraction( contractions_two, moment_coord, moment_orders, - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): """Return the evaluations of the given contractions at the given coordinates. @@ -76,7 +76,7 @@ def construct_array_contraction( Note that a two dimensional array must be given, even if there is only one set of orders of the moment. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance @@ -188,7 +188,7 @@ class if the first four indices correspond to the segmented contraction and the def moment_integral( - basis, moment_coord, moment_orders, transform=None, screen_basis=True, tol_screen=1e-8 + basis, moment_coord, moment_orders, transform=None, screen_basis=False, tol_screen=1e-8 ): r"""Return moment integral of the given basis set. @@ -217,7 +217,7 @@ def moment_integral( and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/momentum.py b/gbasis/integrals/momentum.py index 86a2fc37..b505a0b8 100644 --- a/gbasis/integrals/momentum.py +++ b/gbasis/integrals/momentum.py @@ -57,7 +57,7 @@ class MomentumIntegral(BaseTwoIndexSymmetric): @staticmethod def construct_array_contraction( - contractions_one, contractions_two, screen_basis=True, tol_screen=1e-8 + contractions_one, contractions_two, screen_basis=False, tol_screen=1e-8 ): """Return the integrals over the momentum operator of the given contractions. @@ -70,7 +70,7 @@ def construct_array_contraction( Contracted Cartesian Gaussians (of the same shell) associated with the second index of the momentum energy integral. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance @@ -140,7 +140,7 @@ def construct_array_contraction( return -1j * np.transpose(output, (1, 2, 3, 4, 0)) -def momentum_integral(basis, transform=None, screen_basis=True, tol_screen=1e-8): +def momentum_integral(basis, transform=None, screen_basis=False, tol_screen=1e-8): r"""Return integral over momentum operator of the given basis set. .. math:: @@ -166,7 +166,7 @@ def momentum_integral(basis, transform=None, screen_basis=True, tol_screen=1e-8) Default is no transformation. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/nuclear_electron_attraction.py b/gbasis/integrals/nuclear_electron_attraction.py index 00e2241e..4efd5180 100644 --- a/gbasis/integrals/nuclear_electron_attraction.py +++ b/gbasis/integrals/nuclear_electron_attraction.py @@ -1,9 +1,12 @@ """Module for computing the nuclear electron attraction.""" + from gbasis.integrals.point_charge import point_charge_integral import numpy as np -def nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, transform=None): +def nuclear_electron_attraction_integral( + basis, nuclear_coords, nuclear_charges, transform=None, screen_basis=False, tol_screen=1e-8 +): """Return the nuclear electron attraction integrals of the basis set in the Cartesian form. .. math:: @@ -27,6 +30,13 @@ def nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` and index 0 of the array for contractions. Default is no transformation. + screen_basis : bool, optional + Whether to enable screening of the basis functions. Default value is False. + tol_screen : float, optional + Screening tolerance for excluding evaluations. Integrals with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -37,6 +47,13 @@ def nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, """ return np.sum( - point_charge_integral(basis, nuclear_coords, nuclear_charges, transform=transform), + point_charge_integral( + basis, + nuclear_coords, + nuclear_charges, + transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, + ), axis=2, ) diff --git a/gbasis/integrals/overlap.py b/gbasis/integrals/overlap.py index 7ee3c9ce..d5b7b69f 100644 --- a/gbasis/integrals/overlap.py +++ b/gbasis/integrals/overlap.py @@ -53,7 +53,7 @@ class Overlap(BaseTwoIndexSymmetric): @staticmethod def construct_array_contraction( - contractions_one, contractions_two, screen_basis=True, tol_screen=1e-8 + contractions_one, contractions_two, screen_basis=False, tol_screen=1e-8 ): """Return the evaluations of the given contractions at the given coordinates. @@ -66,7 +66,7 @@ def construct_array_contraction( Contracted Cartesian Gaussians (of the same shell) associated with the second index of the overlap. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance @@ -134,7 +134,7 @@ def construct_array_contraction( )[0] -def overlap_integral(basis, transform=None, screen_basis=True, tol_screen=1e-8): +def overlap_integral(basis, transform=None, screen_basis=False, tol_screen=1e-8): r"""Return overlap integral of the given basis set. .. math:: @@ -154,7 +154,7 @@ def overlap_integral(basis, transform=None, screen_basis=True, tol_screen=1e-8): and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/overlap_asymm.py b/gbasis/integrals/overlap_asymm.py index 0a4c80a6..28f61f9d 100644 --- a/gbasis/integrals/overlap_asymm.py +++ b/gbasis/integrals/overlap_asymm.py @@ -65,7 +65,7 @@ class OverlapAsymmetric(BaseTwoIndexAsymmetric): def overlap_integral_asymmetric( - basis_one, basis_two, transform_one=None, transform_two=None, screen_basis=True, tol_screen=1e-8 + basis_one, basis_two, transform_one=None, transform_two=None, screen_basis=False, tol_screen=1e-8 ): r"""Return overlap integrals between two basis sets. @@ -94,7 +94,7 @@ def overlap_integral_asymmetric( and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/point_charge.py b/gbasis/integrals/point_charge.py index 294f1bd2..e8ce88a7 100644 --- a/gbasis/integrals/point_charge.py +++ b/gbasis/integrals/point_charge.py @@ -1,7 +1,9 @@ """Module for computing point charge integrals.""" + from gbasis.base_two_symm import BaseTwoIndexSymmetric from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals._one_elec_int import _compute_one_elec_integrals +from gbasis.screening import is_two_index_overlap_screened, get_points_mask_for_contraction import numpy as np from scipy.special import hyp1f1 # pylint: disable=E0611 @@ -117,7 +119,13 @@ def boys_func(orders, weighted_dist): @classmethod def construct_array_contraction( - cls, contractions_one, contractions_two, points_coords, points_charge + cls, + contractions_one, + contractions_two, + points_coords, + points_charge, + screen_basis=False, + tol_screen=1e-8, ): r"""Return point charge interaction integral for the given contractions and point charges. @@ -135,6 +143,13 @@ def construct_array_contraction( :math:`x, y, \text{and} z` components. points_charge : np.ndarray(N) Charge of each point charge. + screen_basis : bool, optional + Whether to enable screening of the basis functions. Default value is False. + tol_screen : float, optional + Screening tolerance for excluding evaluations. Integrals with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -198,7 +213,26 @@ def construct_array_contraction( "`points_charge`." ) - # TODO: Overlap screening + # Overlap screening + if screen_basis and is_two_index_overlap_screened( + contractions_one, contractions_two, tol_screen + ): + # Shape of output array: + # axis 0 : index for segmented contractions of contraction one + # axis 1 : angular momentum vector of contraction one (in the same order as angmoms_a) + # axis 2 : index for segmented contractions of contraction two + # axis 3 : angular momentum vector of contraction two (in the same order as angmoms_b) + # axis 4 : point charge + return np.zeros( + ( + contractions_one.num_seg_cont, + len(contractions_one.angmom_components_cart), + contractions_two.num_seg_cont, + len(contractions_two.angmom_components_cart), + points_coords.shape[0], + ), + dtype=np.float64, + ) coord_a = contractions_one.coord angmom_a = contractions_one.angmom @@ -221,18 +255,66 @@ def construct_array_contraction( coeffs_a, coeffs_b = coeffs_b, coeffs_a ab_swapped = True - integrals = _compute_one_elec_integrals( - points_coords, - cls.boys_func, - coord_a, - angmom_a, - exps_a, - coeffs_a, - coord_b, - angmom_b, - exps_b, - coeffs_b, - ) + if screen_basis: + pts_mask_one = get_points_mask_for_contraction( + contractions_one, points_coords, deriv_order=0, tol_screen=tol_screen + ) + pts_mask_two = get_points_mask_for_contraction( + contractions_two, points_coords, deriv_order=0, tol_screen=tol_screen + ) + pts_mask = np.logical_and(pts_mask_one, pts_mask_two) + + # NOTE: Ordering convention for horizontal recursion of integrals + # axis 0 : a_x (size: angmom_a + 1) + # axis 1 : a_y (size: angmom_a + 1) + # axis 2 : a_z (size: angmom_a + 1) + # axis 3 : b_x (size: angmom_b + 1) + # axis 4 : b_y (size: angmom_b + 1) + # axis 5 : b_z (size: angmom_b + 1) + # axis 6 : point charge (size: N) + # axis 7 : index for segmented contractions of contraction a (size: M_a) + # axis 8 : index for segmented contractions of contraction b (size: M_b) + integrals = np.zeros( + ( + angmom_a + 1, + angmom_a + 1, + angmom_a + 1, + angmom_b + 1, + angmom_b + 1, + angmom_b + 1, + points_coords.shape[0], + coeffs_a.shape[1], + coeffs_b.shape[1], + ), + dtype=np.float64, + ) + if np.any(pts_mask): + integrals[:, :, :, :, :, :, pts_mask, :, :] = _compute_one_elec_integrals( + points_coords[pts_mask], + cls.boys_func, + coord_a, + angmom_a, + exps_a, + coeffs_a, + coord_b, + angmom_b, + exps_b, + coeffs_b, + ) + + else: + integrals = _compute_one_elec_integrals( + points_coords, + cls.boys_func, + coord_a, + angmom_a, + exps_a, + coeffs_a, + coord_b, + angmom_b, + exps_b, + coeffs_b, + ) integrals = np.transpose(integrals, (7, 0, 1, 2, 8, 3, 4, 5, 6)) angmoms_a_x, angmoms_a_y, angmoms_a_z = angmoms_a.T @@ -266,7 +348,9 @@ def construct_array_contraction( return output -def point_charge_integral(basis, points_coords, points_charge, transform=None): +def point_charge_integral( + basis, points_coords, points_charge, transform=None, screen_basis=False, tol_screen=1e-8 +): r"""Return the point-charge interaction integrals of basis set in the given coordinate systems. .. math:: @@ -293,6 +377,13 @@ def point_charge_integral(basis, points_coords, points_charge, transform=None): Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` and index 0 of the array for contractions. Default is no transformation. + screen_basis : bool, optional + Whether to enable screening of the basis functions. Default value is False. + tol_screen : float, optional + Screening tolerance for excluding evaluations. Integrals with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -305,19 +396,24 @@ def point_charge_integral(basis, points_coords, points_charge, transform=None): """ coord_type = [ct for ct in [shell.coord_type for shell in basis]] + kwargs = {"tol_screen": tol_screen, "screen_basis": screen_basis} if transform is not None: return PointChargeIntegral(basis).construct_array_lincomb( - transform, coord_type, points_coords=points_coords, points_charge=points_charge + transform, + coord_type, + points_coords=points_coords, + points_charge=points_charge, + **kwargs, ) if all(ct == "cartesian" for ct in coord_type): return PointChargeIntegral(basis).construct_array_cartesian( - points_coords=points_coords, points_charge=points_charge + points_coords=points_coords, points_charge=points_charge, **kwargs ) if all(ct == "spherical" for ct in coord_type): return PointChargeIntegral(basis).construct_array_spherical( - points_coords=points_coords, points_charge=points_charge + points_coords=points_coords, points_charge=points_charge, **kwargs ) return PointChargeIntegral(basis).construct_array_mix( - coord_type, points_coords=points_coords, points_charge=points_charge + coord_type, points_coords=points_coords, points_charge=points_charge, **kwargs ) diff --git a/tests/test_angular_momentum.py b/tests/test_angular_momentum.py index f915e26c..b23fd688 100644 --- a/tests/test_angular_momentum.py +++ b/tests/test_angular_momentum.py @@ -1,4 +1,5 @@ """Test gbasis.integrals.angular_momentum.""" + from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals._diff_operator_int import ( _compute_differential_operator_integrals_intermediate, @@ -11,7 +12,8 @@ from utils import find_datafile -def test_angular_momentum_construct_array_contraction(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_angular_momentum_construct_array_contraction(screen_basis, tol_screen): """Test integrals.angular_momentum.angular_momentumIntegral.construct_array_contraction.""" test_one = GeneralizedContractionShell( 1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), "spherical" @@ -346,7 +348,9 @@ def test_angular_momentum_construct_array_contraction(): int_pz_dyz = int_pz_dyz.squeeze() int_pz_dzz = int_pz_dzz.squeeze() - test = AngularMomentumIntegral.construct_array_contraction(test_one, test_two) + test = AngularMomentumIntegral.construct_array_contraction( + test_one, test_two, screen_basis=screen_basis, tol_screen=tol_screen + ) assert test.shape == (1, 3, 1, 6, 3) assert np.allclose( test.squeeze(), @@ -358,6 +362,7 @@ def test_angular_momentum_construct_array_contraction(): [int_pz_dxx, int_pz_dxy, int_pz_dxz, int_pz_dyy, int_pz_dyz, int_pz_dzz], ] ), + atol=tol_screen, ) with pytest.raises(TypeError): @@ -366,50 +371,66 @@ def test_angular_momentum_construct_array_contraction(): AngularMomentumIntegral.construct_array_contraction(None, test_two) -def test_angular_momentum_integral_cartesian(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_angular_momentum_integral_cartesian(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_cartesian.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "cartesian") angular_momentum_integral_obj = AngularMomentumIntegral(basis) assert np.allclose( - angular_momentum_integral_obj.construct_array_cartesian(), - angular_momentum_integral(basis), + angular_momentum_integral_obj.construct_array_cartesian( + screen_basis=screen_basis, tol_screen=tol_screen + ), + angular_momentum_integral(basis, screen_basis=False), + atol=tol_screen, ) -def test_angular_momentum_integral_spherical(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_angular_momentum_integral_spherical(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_spherical.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") angular_momentum_integral_obj = AngularMomentumIntegral(basis) assert np.allclose( - angular_momentum_integral_obj.construct_array_spherical(), - angular_momentum_integral(basis), + angular_momentum_integral_obj.construct_array_spherical( + screen_basis=screen_basis, tol_screen=tol_screen + ), + angular_momentum_integral(basis, screen_basis=False), + atol=tol_screen, ) -def test_angular_momentum_integral_mix(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_angular_momentum_integral_mix(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_mix.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["spherical"] * 8) angular_momentum_integral_obj = AngularMomentumIntegral(basis) assert np.allclose( - angular_momentum_integral_obj.construct_array_mix(["spherical"] * 8), - angular_momentum_integral(basis), + angular_momentum_integral_obj.construct_array_mix( + ["spherical"] * 8, screen_basis=screen_basis, tol_screen=tol_screen + ), + angular_momentum_integral(basis, screen_basis=False), + atol=tol_screen, ) -def test_angular_momentum_integral_lincomb(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_angular_momentum_integral_lincomb(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_lincomb.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") angular_momentum_integral_obj = AngularMomentumIntegral(basis) transform = np.random.rand(14, 18) assert np.allclose( - angular_momentum_integral_obj.construct_array_lincomb(transform, ["spherical"]), - angular_momentum_integral(basis, transform), + angular_momentum_integral_obj.construct_array_lincomb( + transform, ["spherical"], screen_basis=screen_basis, tol_screen=tol_screen + ), + angular_momentum_integral(basis, transform, screen_basis=False), + atol=tol_screen, ) @@ -422,8 +443,8 @@ def test_angular_momentum_screening_accuracy(precision): atcoords = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]]) contraction = make_contractions(basis_dict, atsymbols, atcoords, "cartesian") - # the screening tolerance needs to be 1e-4 times the desired precision - tol_screen = precision * 1e-4 + # the screening tolerance needs to be 1e-2 times the desired precision + tol_screen = precision * 1e-2 angular_momentum = angular_momentum_integral(contraction, tol_screen=tol_screen) angular_momentum_no_screen = angular_momentum_integral(contraction, screen_basis=False) assert np.allclose(angular_momentum, angular_momentum_no_screen, atol=precision) diff --git a/tests/test_density.py b/tests/test_density.py index 844d2e0c..8a69f648 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -60,8 +60,7 @@ def test_evaluate_density_using_evaluated_orbs(): evaluate_density_using_evaluated_orbs(density_mat, orb_eval) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) def test_evaluate_density(screen_basis, tol_screen): """Test gbasis.evals.density.evaluate_density.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) @@ -83,8 +82,7 @@ def test_evaluate_density(screen_basis, tol_screen): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_deriv_density(screen_basis, tol_screen, deriv_type): """Test gbasis.evals.density.evaluate_deriv_density.""" @@ -272,8 +270,7 @@ def test_evaluate_deriv_density(screen_basis, tol_screen, deriv_type): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_density_gradient(screen_basis, tol_screen, deriv_type): """Test gbasis.evals.density.evaluate_density_gradient.""" @@ -338,8 +335,7 @@ def test_evaluate_density_gradient(screen_basis, tol_screen, deriv_type): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) def test_evaluate_density_horton(screen_basis, tol_screen): """Test gbasis.evals.density.evaluate_density against result from HORTON. @@ -370,8 +366,7 @@ def test_evaluate_density_horton(screen_basis, tol_screen): assert np.allclose(dens, horton_density, atol=tol_screen) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_density_gradient_horton(screen_basis, tol_screen, deriv_type): """Test gbasis.evals.density.evaluate_density_gradient against result from HORTON. @@ -405,8 +400,7 @@ def test_evaluate_density_gradient_horton(screen_basis, tol_screen, deriv_type): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_hessian_deriv_horton(screen_basis, tol_screen, deriv_type): """Test gbasis.evals.density.evaluate_density_hessian against result from HORTON. @@ -447,8 +441,7 @@ def test_evaluate_hessian_deriv_horton(screen_basis, tol_screen, deriv_type): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_laplacian_deriv_horton(screen_basis, tol_screen, deriv_type): """Test gbasis.evals.density.evaluate_density_laplacian against result from HORTON. @@ -483,8 +476,7 @@ def test_evaluate_laplacian_deriv_horton(screen_basis, tol_screen, deriv_type): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_posdef_kinetic_energy_density(screen_basis, tol_screen, deriv_type): """Test evaluate_posdef_kinetic_energy_density against results from HORTON. @@ -519,8 +511,7 @@ def test_evaluate_posdef_kinetic_energy_density(screen_basis, tol_screen, deriv_ assert np.allclose(dens, horton_density_kinetic_density, atol=tol_screen) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_general_kinetic_energy_density_horton(screen_basis, tol_screen, deriv_type): """Test evaluate_general_kinetic_energy_density against results from HORTON. @@ -560,8 +551,7 @@ def test_evaluate_general_kinetic_energy_density_horton(screen_basis, tol_screen ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_general_kinetic_energy_density(screen_basis, tol_screen, deriv_type): """Test density.evaluate_general_kinetic_energy_density.""" diff --git a/tests/test_nuclear_electron_attraction.py b/tests/test_nuclear_electron_attraction.py index 0e67871b..1840ed3e 100644 --- a/tests/test_nuclear_electron_attraction.py +++ b/tests/test_nuclear_electron_attraction.py @@ -1,17 +1,21 @@ """Test gbasis.integrals.nuclear_electron_attraction.""" + from gbasis.integrals.nuclear_electron_attraction import nuclear_electron_attraction_integral from gbasis.integrals.point_charge import point_charge_integral from gbasis.parsers import make_contractions, parse_nwchem import numpy as np from utils import find_datafile, HortonContractions +import pytest -def test_nuclear_electron_attraction_horton_anorcc_hhe(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_horton_anorcc_hhe(screen_basis, tol_screen): """Test nuclear_electron_attraciton.nuclear_electron_attraction_cartesian with HORTON. The test case is diatomic with H and He separated by 0.8 angstroms with basis set ANO-RCC. """ + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem")) # NOTE: used HORTON's conversion factor for angstroms to bohr coords = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]]) @@ -20,17 +24,20 @@ def test_nuclear_electron_attraction_horton_anorcc_hhe(): horton_nucattract = np.load(find_datafile("data_horton_hhe_cart_nucattract.npy")) assert np.allclose( - nuclear_electron_attraction_integral(basis, coords, np.array([1, 2])), + nuclear_electron_attraction_integral(basis, coords, np.array([1, 2]), **kwargs), horton_nucattract, + atol=tol_screen, ) -def test_nuclear_electron_attraction_horton_anorcc_bec(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_horton_anorcc_bec(screen_basis, tol_screen): """Test nuclear_electron_attraciton.nuclear_electron_attraction_cartesian with HORTON. The test case is diatomic with B and C separated by 1.0 angstroms with basis set ANO-RCC. """ + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem")) # NOTE: used HORTON's conversion factor for angstroms to bohr coords = np.array([[0, 0, 0], [1.0 * 1.0 / 0.5291772083, 0, 0]]) @@ -39,65 +46,80 @@ def test_nuclear_electron_attraction_horton_anorcc_bec(): horton_nucattract = np.load(find_datafile("data_horton_bec_cart_nucattract.npy")) assert np.allclose( - nuclear_electron_attraction_integral(basis, coords, np.array([4, 6])), + nuclear_electron_attraction_integral(basis, coords, np.array([4, 6]), **kwargs), horton_nucattract, + atol=tol_screen, ) -def test_nuclear_electron_attraction_cartesian(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_cartesian(screen_basis, tol_screen): """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_cartesian.""" + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "cartesian") nuclear_coords = np.random.rand(5, 3) nuclear_charges = np.random.rand(5) - ref = point_charge_integral(basis, nuclear_coords, nuclear_charges) + ref = point_charge_integral(basis, nuclear_coords, nuclear_charges, **kwargs) assert np.allclose( ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4], - nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges), + nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, **kwargs), + atol=1e-10, ) -def test_nuclear_electron_attraction_spherical(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_spherical(screen_basis, tol_screen): """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_spherical.""" + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") nuclear_coords = np.random.rand(5, 3) nuclear_charges = np.random.rand(5) - ref = point_charge_integral(basis, nuclear_coords, nuclear_charges) + ref = point_charge_integral(basis, nuclear_coords, nuclear_charges, **kwargs) assert np.allclose( ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4], - nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges), + nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, **kwargs), + atol=1e-10, ) -def test_nuclear_electron_attraction_mix(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_mix(screen_basis, tol_screen): """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_mix.""" + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["spherical"] * 8) nuclear_coords = np.random.rand(5, 3) nuclear_charges = np.random.rand(5) - ref = point_charge_integral(basis, nuclear_coords, nuclear_charges) + ref = point_charge_integral(basis, nuclear_coords, nuclear_charges, **kwargs) assert np.allclose( ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4], - nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges), + nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, **kwargs), + atol=1e-10, ) -def test_nuclear_electron_attraction_lincomb(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_lincomb(screen_basis, tol_screen): """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_lincomb.""" + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") nuclear_coords = np.random.rand(5, 3) nuclear_charges = np.random.rand(5) transform = np.random.rand(14, 18) - ref = point_charge_integral(basis, nuclear_coords, nuclear_charges, transform=transform) + ref = point_charge_integral( + basis, nuclear_coords, nuclear_charges, transform=transform, **kwargs + ) assert np.allclose( ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4], nuclear_electron_attraction_integral( - basis, nuclear_coords, nuclear_charges, transform=transform + basis, nuclear_coords, nuclear_charges, transform=transform, **kwargs ), + atol=1e-10, ) diff --git a/tests/test_point_charge.py b/tests/test_point_charge.py index 0b4fd6a7..7d0819ca 100644 --- a/tests/test_point_charge.py +++ b/tests/test_point_charge.py @@ -1,4 +1,5 @@ """Test gbasis.integrals.point_charge.""" + from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals.point_charge import point_charge_integral, PointChargeIntegral from gbasis.parsers import make_contractions, parse_nwchem @@ -30,7 +31,8 @@ def test_boys_func(): assert np.allclose(test, ref.reshape(10, 20, 30)) -def test_construct_array_contraction(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_construct_array_contraction(screen_basis, tol_screen): """Test gbasis.integrals.point_charge.PointChargeIntegral.construct_array_contraction.""" coord_one = np.array([0.5, 1, 1.5]) test_one = GeneralizedContractionShell( @@ -58,7 +60,9 @@ def test_construct_array_contraction(): PointChargeIntegral.construct_array_contraction(test_one, test_two, coord, np.array([0, 1])) assert np.allclose( - PointChargeIntegral.construct_array_contraction(test_one, test_two, coord, charge), + PointChargeIntegral.construct_array_contraction( + test_one, test_two, coord, charge, screen_basis=screen_basis, tol_screen=tol_screen + ), ( 2 * np.pi @@ -71,6 +75,7 @@ def test_construct_array_contraction(): * 3 * (-1) ), + atol=tol_screen, ) test_one = GeneralizedContractionShell( @@ -88,7 +93,9 @@ def test_construct_array_contraction(): for i in range(2) ] assert np.allclose( - PointChargeIntegral.construct_array_contraction(test_one, test_two, coord, charge).ravel(), + PointChargeIntegral.construct_array_contraction( + test_one, test_two, coord, charge, screen_basis=screen_basis, tol_screen=tol_screen + ).ravel(), ( ((coord_wac - coord_one) * v_000_000[0] - (coord_wac - coord[0]) * v_000_000[1]) * (2 * 0.1 / np.pi) ** (3 / 4) @@ -98,6 +105,7 @@ def test_construct_array_contraction(): * 3 * (-1) ), + atol=tol_screen, ) test_one = GeneralizedContractionShell( @@ -107,7 +115,9 @@ def test_construct_array_contraction(): 1, np.array([1.5, 2, 3]), np.array([3.0]), np.array([0.2]), "spherical" ) assert np.allclose( - PointChargeIntegral.construct_array_contraction(test_one, test_two, coord, charge).ravel(), + PointChargeIntegral.construct_array_contraction( + test_one, test_two, coord, charge, screen_basis=screen_basis, tol_screen=tol_screen + ).ravel(), ( ((coord_wac - coord_two) * v_000_000[0] - (coord_wac - coord[0]) * v_000_000[1]) * (2 * 0.1 / np.pi) ** (3 / 4) @@ -117,6 +127,7 @@ def test_construct_array_contraction(): * 3 * (-1) ), + atol=tol_screen, ) @@ -130,7 +141,8 @@ def test_point_charge_cartesian(): points_charge = np.random.rand(5) assert np.allclose( point_charge_obj.construct_array_cartesian( - points_coords=points_coords, points_charge=points_charge + points_coords=points_coords, + points_charge=points_charge, ), point_charge_integral(basis, points_coords=points_coords, points_charge=points_charge), ) diff --git a/tests/test_stress_tensor.py b/tests/test_stress_tensor.py index a9ef3f56..da45c979 100644 --- a/tests/test_stress_tensor.py +++ b/tests/test_stress_tensor.py @@ -1,4 +1,5 @@ """Test gbasis.evals.stress_tensor.""" + from gbasis.evals.density import ( evaluate_density_laplacian, evaluate_deriv_density, @@ -15,7 +16,9 @@ from utils import find_datafile, HortonContractions -def test_evaluate_stress_tensor(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_evaluate_stress_tensor(screen_basis, tol_screen): """Test gbasis.evals.stress_tensor.evaluate_stress_tensor.""" basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem")) coords = np.array([[0, 0, 0]]) @@ -29,10 +32,47 @@ def test_evaluate_stress_tensor(): with pytest.raises(TypeError): evaluate_stress_tensor(np.identity(40), basis, points, 1, 0j, np.identity(40)) - test_a = evaluate_stress_tensor(np.identity(40), basis, points, 0, 0, np.identity(40)) - test_b = evaluate_stress_tensor(np.identity(40), basis, points, 1, 0, np.identity(40)) - test_c = evaluate_stress_tensor(np.identity(40), basis, points, 1, 2, np.identity(40)) - test_d = evaluate_stress_tensor(np.identity(40), basis, points, 0.5, 2, np.identity(40)) + test_a = evaluate_stress_tensor( + np.identity(40), + basis, + points, + 0, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_b = evaluate_stress_tensor( + np.identity(40), + basis, + points, + 1, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_c = evaluate_stress_tensor( + np.identity(40), + basis, + points, + 1, + 2, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_d = evaluate_stress_tensor( + np.identity(40), + basis, + points, + 0.5, + 2, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + # compute reference without screening for i in range(3): for j in range(3): orders_i = np.array([0, 0, 0]) @@ -41,10 +81,22 @@ def test_evaluate_stress_tensor(): orders_j[j] += 1 temp1 = evaluate_deriv_reduced_density_matrix( - orders_i, orders_j, np.identity(40), basis, points, np.identity(40) + orders_i, + orders_j, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) temp2 = evaluate_deriv_reduced_density_matrix( - orders_j, orders_i, np.identity(40), basis, points, np.identity(40) + orders_j, + orders_i, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) temp3 = evaluate_deriv_reduced_density_matrix( orders_i + orders_j, @@ -53,6 +105,7 @@ def test_evaluate_stress_tensor(): basis, points, np.identity(40), + screen_basis=False, ) temp4 = evaluate_deriv_reduced_density_matrix( orders_i + orders_j, @@ -61,20 +114,28 @@ def test_evaluate_stress_tensor(): basis, points, np.identity(40), + screen_basis=False, ) if i == j: - temp5 = evaluate_density_laplacian(np.identity(40), basis, points, np.identity(40)) + temp5 = evaluate_density_laplacian( + np.identity(40), basis, points, np.identity(40), screen_basis=False + ) else: temp5 = 0 - assert np.allclose(test_a[:, i, j], 0.5 * temp3 + 0.5 * temp4) - assert np.allclose(test_b[:, i, j], -0.5 * temp1 - 0.5 * temp2) - assert np.allclose(test_c[:, i, j], -0.5 * temp1 - 0.5 * temp2 - temp5) + # check that non screened reference matches screened result within tol_screen + assert np.allclose(test_a[:, i, j], 0.5 * temp3 + 0.5 * temp4, atol=tol_screen) + assert np.allclose(test_b[:, i, j], -0.5 * temp1 - 0.5 * temp2, atol=tol_screen) + assert np.allclose(test_c[:, i, j], -0.5 * temp1 - 0.5 * temp2 - temp5, atol=tol_screen) assert np.allclose( - test_d[:, i, j], -0.25 * temp1 - 0.25 * temp2 + 0.25 * temp3 + 0.25 * temp4 - temp5 + test_d[:, i, j], + -0.25 * temp1 - 0.25 * temp2 + 0.25 * temp3 + 0.25 * temp4 - temp5, + atol=tol_screen, ) -def test_evaluate_ehrenfest_force(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_evaluate_ehrenfest_force(screen_basis, tol_screen): """Test gbasis.evals.stress_tensor.evaluate_ehrenfest_force.""" basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem")) coords = np.array([[0, 0, 0]]) @@ -88,16 +149,53 @@ def test_evaluate_ehrenfest_force(): with pytest.raises(TypeError): evaluate_ehrenfest_force(np.identity(40), basis, points, 1, 0j, np.identity(40)) - test_a = evaluate_ehrenfest_force(np.identity(40), basis, points, 0, 0, np.identity(40)) - test_b = evaluate_ehrenfest_force(np.identity(40), basis, points, 1, 0, np.identity(40)) - test_c = evaluate_ehrenfest_force(np.identity(40), basis, points, 0.5, 0, np.identity(40)) - test_d = evaluate_ehrenfest_force(np.identity(40), basis, points, 0, 2, np.identity(40)) + test_a = evaluate_ehrenfest_force( + np.identity(40), + basis, + points, + 0, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_b = evaluate_ehrenfest_force( + np.identity(40), + basis, + points, + 1, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_c = evaluate_ehrenfest_force( + np.identity(40), + basis, + points, + 0.5, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_d = evaluate_ehrenfest_force( + np.identity(40), + basis, + points, + 0, + 2, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) for j in range(3): ref_a = np.zeros(points.shape[0]) ref_b = np.zeros(points.shape[0]) ref_c = np.zeros(points.shape[0]) ref_d = np.zeros(points.shape[0]) + # compute reference without screening orders_j = np.array([0, 0, 0]) orders_j[j] += 1 for i in range(3): @@ -105,10 +203,22 @@ def test_evaluate_ehrenfest_force(): orders_i[i] += 1 temp1 = evaluate_deriv_reduced_density_matrix( - 2 * orders_i, orders_j, np.identity(40), basis, points, np.identity(40) + 2 * orders_i, + orders_j, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) temp2 = evaluate_deriv_reduced_density_matrix( - orders_i, orders_i + orders_j, np.identity(40), basis, points, np.identity(40) + orders_i, + orders_i + orders_j, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) temp3 = evaluate_deriv_reduced_density_matrix( 2 * orders_i + orders_j, @@ -117,25 +227,40 @@ def test_evaluate_ehrenfest_force(): basis, points, np.identity(40), + screen_basis=False, ) temp4 = evaluate_deriv_reduced_density_matrix( - orders_i + orders_j, orders_i, np.identity(40), basis, points, np.identity(40) + orders_i + orders_j, + orders_i, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) temp5 = evaluate_deriv_density( - 2 * orders_i + orders_j, np.identity(40), basis, points, np.identity(40) + 2 * orders_i + orders_j, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) ref_a += temp3 + temp4 ref_b += -temp1 - temp2 ref_c += -0.5 * temp1 - 0.5 * temp2 + 0.5 * temp3 + 0.5 * temp4 ref_d += temp3 + temp4 - temp5 - assert np.allclose(test_a[:, j], -ref_a) - assert np.allclose(test_b[:, j], -ref_b) - assert np.allclose(test_c[:, j], -ref_c) - assert np.allclose(test_d[:, j], -ref_d) + # check that non screened reference matches screened result within tol_screen + assert np.allclose(test_a[:, j], -ref_a, atol=tol_screen) + assert np.allclose(test_b[:, j], -ref_b, atol=tol_screen) + assert np.allclose(test_c[:, j], -ref_c, atol=tol_screen) + assert np.allclose(test_d[:, j], -ref_d, atol=tol_screen) -def test_evaluate_ehrenfest_hessian(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_evaluate_ehrenfest_hessian(screen_basis, tol_screen): """Test gbasis.evals.stress_tensor.evaluate_ehrenfest_hessian.""" basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem")) coords = np.array([[0, 0, 0]]) @@ -151,10 +276,48 @@ def test_evaluate_ehrenfest_hessian(): with pytest.raises(TypeError): evaluate_ehrenfest_hessian(np.identity(40), basis, points, 1, 0j, np.identity(40)) - test_a = evaluate_ehrenfest_hessian(np.identity(40), basis, points, 0, 0, np.identity(40)) - test_b = evaluate_ehrenfest_hessian(np.identity(40), basis, points, 1, 0, np.identity(40)) - test_c = evaluate_ehrenfest_hessian(np.identity(40), basis, points, 0.5, 0, np.identity(40)) - test_d = evaluate_ehrenfest_hessian(np.identity(40), basis, points, 0, 2, np.identity(40)) + test_a = evaluate_ehrenfest_hessian( + np.identity(40), + basis, + points, + 0, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_b = evaluate_ehrenfest_hessian( + np.identity(40), + basis, + points, + 1, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_c = evaluate_ehrenfest_hessian( + np.identity(40), + basis, + points, + 0.5, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_d = evaluate_ehrenfest_hessian( + np.identity(40), + basis, + points, + 0, + 2, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + + # compute reference without screening for j in range(3): for k in range(3): ref_a = np.zeros(points.shape[0]) @@ -177,6 +340,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp2 = evaluate_deriv_reduced_density_matrix( 2 * orders_i, @@ -185,6 +349,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp3 = evaluate_deriv_reduced_density_matrix( orders_i + orders_k, @@ -193,6 +358,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp4 = evaluate_deriv_reduced_density_matrix( orders_i, @@ -201,6 +367,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp5 = evaluate_deriv_reduced_density_matrix( 2 * orders_i + orders_j + orders_k, @@ -209,6 +376,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp6 = evaluate_deriv_reduced_density_matrix( 2 * orders_i + orders_j, @@ -217,6 +385,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp7 = evaluate_deriv_reduced_density_matrix( orders_i + orders_j + orders_k, @@ -225,6 +394,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp8 = evaluate_deriv_reduced_density_matrix( orders_i + orders_j, @@ -233,6 +403,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp9 = evaluate_deriv_density( 2 * orders_i + orders_j + orders_k, @@ -240,6 +411,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) ref_a += temp5 + temp6 + temp7 + temp8 @@ -255,18 +427,30 @@ def test_evaluate_ehrenfest_hessian(): + 0.5 * temp8 ) ref_d += temp3 + temp4 + temp5 + temp6 - temp9 - assert np.allclose(test_a[:, j, k], -ref_a) - assert np.allclose(test_b[:, j, k], -ref_b) - assert np.allclose(test_c[:, j, k], -ref_c) - assert np.allclose(test_d[:, j, k], -ref_d) + assert np.allclose(test_a[:, j, k], -ref_a, atol=tol_screen) + assert np.allclose(test_b[:, j, k], -ref_b, atol=tol_screen) + assert np.allclose(test_c[:, j, k], -ref_c, atol=tol_screen) + assert np.allclose(test_d[:, j, k], -ref_d, atol=tol_screen) + # check symmetry if requested assert np.allclose( evaluate_ehrenfest_hessian( - np.identity(40), basis, points, 0, 0, np.identity(40), symmetric=True + np.identity(40), + basis, + points, + 0, + 0, + np.identity(40), + symmetric=True, + screen_basis=False, ), ( - evaluate_ehrenfest_hessian(np.identity(40), basis, points, 0, 0, np.identity(40)) + evaluate_ehrenfest_hessian( + np.identity(40), basis, points, 0, 0, np.identity(40), screen_basis=False + ) + np.swapaxes( - evaluate_ehrenfest_hessian(np.identity(40), basis, points, 0, 0, np.identity(40)), + evaluate_ehrenfest_hessian( + np.identity(40), basis, points, 0, 0, np.identity(40), screen_basis=False + ), 1, 2, )