Skip to content
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
Add screening and test to evaluate_stress_tensor
  • Loading branch information
marco-2023 committed Oct 25, 2025
commit 64f58cec244b6b4ed6cdea4714641ff83307d094
26 changes: 25 additions & 1 deletion gbasis/evals/stress_tensor.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
"""Module for computing properties related to the stress tensor."""

from gbasis.evals.density import (
evaluate_density_laplacian,
evaluate_deriv_density,
Expand All @@ -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=True,
tol_screen=1e-8,
):
r"""Return the stress tensor evaluated at the given coordinates.

Stress tensor is defined here as:
Expand Down Expand Up @@ -70,6 +80,14 @@ 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 True
(enable screening).
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
-------
Expand Down Expand Up @@ -99,6 +117,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(
Expand All @@ -108,6 +128,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] -= (
Expand All @@ -118,6 +140,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]
Expand Down
83 changes: 71 additions & 12 deletions tests/test_stress_tensor.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
"""Test gbasis.evals.stress_tensor."""

from gbasis.evals.density import (
evaluate_density_laplacian,
evaluate_deriv_density,
Expand All @@ -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]])
Expand All @@ -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])
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -61,16 +114,22 @@ 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,
)


Expand Down