Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
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
48 changes: 40 additions & 8 deletions deepmd/utils/pair_tab.py
Original file line number Diff line number Diff line change
Expand Up @@ -57,17 +57,49 @@ def reinit(self, filename: str, rcut: float | None = None) -> None:
if filename is None:
self.tab_info, self.tab_data = None, None
return
self.vdata = np.loadtxt(filename, dtype=self.data_type)
self.rmin = self.vdata[0][0]
self.rmax = self.vdata[-1][0]
self.hh = self.vdata[1][0] - self.vdata[0][0]
ncol = self.vdata.shape[1] - 1
vdata = np.loadtxt(filename, dtype=self.data_type)
rmin = vdata[0][0]
rmax = vdata[-1][0]
dx = np.diff(vdata[:, 0])
if not np.all(dx > 0):
raise ValueError(
f"The distance grid in the pairwise table {filename} is not "
"strictly increasing. The tabulated potential must be provided "
"on a uniform grid with distances sorted in ascending order and "
"without duplicated rows. Please regrid the table."
)
# validate against absolute node positions rather than per-interval
# spacing: consumers (the C++ kernel and _make_data) index by
# rmin + i * hh, so that is what must stay accurate, not each dx.
n = vdata.shape[0]
hh = (rmax - rmin) / (n - 1)
deviation = np.abs(
vdata[:, 0] - (rmin + hh * np.arange(n, dtype=self.data_type))
)
tol = 1e-2 * abs(hh)
if np.any(deviation > tol):
bad_row = int(np.argmax(deviation > tol))
raise ValueError(
f"The distance grid in the pairwise table {filename} is not "
"evenly spaced. The tabulated potential must be provided on a "
f"uniform grid, but row {bad_row} (distance "
f"{vdata[bad_row, 0]}) does not match the constant step "
f"inferred from rmin and rmax ({hh}). Please regrid the "
"table to use a constant distance step."
)
ncol = vdata.shape[1] - 1
n0 = (-1 + np.sqrt(1 + 8 * ncol)) * 0.5
self.ntypes = int(n0 + 0.1)
assert self.ntypes * (self.ntypes + 1) // 2 == ncol, (
f"number of volumes provided in {filename} does not match guessed number of types {self.ntypes}"
ntypes = int(n0 + 0.1)
assert ntypes * (ntypes + 1) // 2 == ncol, (
f"number of volumes provided in {filename} does not match guessed number of types {ntypes}"
)

self.vdata = vdata
self.rmin = rmin
self.rmax = rmax
self.hh = hh

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the line that keeps the fix from closing the bug it targets, and the reason it is written this way is my comment last round, which I am withdrawing.

The check above proves every node sits within one percent of a cell of rmin + i * hh_ref. But what gets stored, and what every consumer then uses, is the first interval: tab_info[1] feeds uu = (rr - rmin) / hh in pair_tab.cc and the dpmodel/pt _pair_tabulated_inter, _make_data scales derivatives by it, rcut_idx and the padding linspace in _check_table_upper_boundary are computed from it. Nothing bounds how far this value may sit from hh_ref; the check only bounds node positions, so the first interval can be off by nearly a full percent while the table passes. A relative error in hh accumulates linearly in the index.

I reproduced it against this head: linspace(0, 1, 1001) with only row 1 moved by 9.9e-6, i.e. 0.99 of a cell, is accepted, self.hh comes out as 0.0010099, and at r = 0.995 the consumer's index arithmetic lands on cell 985 while the node is 995. Ten cells off, silently, which is precisely the failure mode in the PR title. The neighbouring comment says the consumers' rmin + i * hh is "what must stay accurate", and it is, but it is hh_ref that was made accurate, not hh.

The fix is one line: self.hh = hh_ref, and drop the first-interval hh. I asked last round to keep them separate to avoid changing tab_info for existing models; having worked through it, that caution was misplaced. For a table with a round stride the two values are bit-identical, and for a printed table hh_ref averages out the rounding of row 1 and is the better estimate, so the change to existing models is at the printing-precision level and in the right direction. The error message already reports hh_ref as the step, so this also makes what the object stores agree with what it tells the user.

A test that would have caught it: a uniform grid with a single node perturbed inside the tolerance, asserting tab.hh against the true step with rtol rather than assertAlmostEqual.

self.ntypes = ntypes

# check table data against rcut and update tab_file if needed, table upper boundary is used as rcut if not provided.
self.rcut = rcut if rcut is not None else self.rmax
self._check_table_upper_boundary()
Expand Down
4 changes: 4 additions & 0 deletions doc/model/pairtab.md
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,10 @@ DeePMD-kit also supports combination with a pairwise potential {{ tensorflow_ico

The table file should be a text file that can be read by {py:meth}`numpy.loadtxt`.
The first column is the distance between two atoms, where upper range should be larger than the cutoff radius.
It must be strictly increasing and evenly spaced: every distance has to sit within one percent of a grid step of `rmin + i * hh`,
where `rmin` is the first distance and `hh` is the constant step inferred from the first and last distances.
A table that violates this raises a `ValueError` when the model is constructed, because the spline coefficients and both
evaluators index the table by that constant step and cannot represent a non-uniform grid.
Other columns are two-body interaction energies for pairs of certain types,
in the order of Type_0-Type_0, Type_0-Type_1, ..., Type_0-Type_N, Type_1-Type_1, ..., Type_1-Type_N, ..., and Type_N-Type_N.

Expand Down
139 changes: 139 additions & 0 deletions source/tests/common/dpmodel/test_pairtab_preprocess.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,7 @@
# SPDX-License-Identifier: LGPL-3.0-or-later
import copy
import os
import tempfile
import unittest
from unittest.mock import (
patch,
Expand Down Expand Up @@ -275,5 +278,141 @@ def test_preprocess(self) -> None:
)


class TestPairTabGridSpacing(unittest.TestCase):
@patch("numpy.loadtxt")
def test_non_uniform_grid(self, mock_loadtxt) -> None:
mock_loadtxt.return_value = np.array(
[
[0.00, 1.0],
[0.01, 0.8],
[0.02, 0.6],
[0.09, 0.3],
[0.16, 0.0],
]
)
with self.assertRaisesRegex(ValueError, "evenly spaced"):
PairTab(filename="dummy_path", rcut=0.16)

@patch("numpy.loadtxt")
def test_duplicate_distances(self, mock_loadtxt) -> None:
mock_loadtxt.return_value = np.array(
[
[0.01, 1.0],
[0.01, 0.8],
[0.01, 0.6],
[0.01, 0.0],
]
)
with self.assertRaisesRegex(ValueError, "strictly increasing"):
PairTab(filename="dummy_path", rcut=0.04)

@patch("numpy.loadtxt")
def test_descending_grid(self, mock_loadtxt) -> None:
mock_loadtxt.return_value = np.array(
[
[0.04, 1.0],
[0.03, 0.8],
[0.02, 0.6],
[0.01, 0.0],
]
)
with self.assertRaisesRegex(ValueError, "strictly increasing"):
PairTab(filename="dummy_path", rcut=0.04)

@patch("numpy.loadtxt")
def test_non_uniform_fine_grid(self, mock_loadtxt) -> None:
mock_loadtxt.return_value = np.array(
[
[0.0, 1.0],
[1e-10, 0.8],
[1.1e-9, 0.6],
[2.1e-9, 0.3],
[3.1e-9, 0.0],
]
)
with self.assertRaisesRegex(ValueError, "evenly spaced"):
PairTab(filename="dummy_path", rcut=3.1e-9)

@patch("numpy.loadtxt")
def test_uniform_fine_grid(self, mock_loadtxt) -> None:
mock_loadtxt.return_value = np.array(
[
[0.0, 1.0],
[1e-9, 0.8],
[2e-9, 0.6],
[3e-9, 0.3],
[4e-9, 0.0],
]
)
tab = PairTab(filename="dummy_path", rcut=4e-9)
np.testing.assert_allclose(tab.hh, 1e-9, rtol=1e-6)

@patch("numpy.loadtxt")
def test_hh_from_node_positions_not_first_interval(self, mock_loadtxt) -> None:
rr = np.linspace(0.0, 1.0, 1001)
rr[1] += 9.9e-6
mock_loadtxt.return_value = np.stack((rr, np.zeros_like(rr)), axis=1)
tab = PairTab(filename="dummy_path", rcut=1.0)
np.testing.assert_allclose(tab.hh, 1.0 / 1000, rtol=1e-6)

@patch("numpy.loadtxt")
def test_failed_reinit_keeps_state(self, mock_loadtxt) -> None:
uniform = np.array(
[
[0.00, 1.0],
[0.01, 0.8],
[0.02, 0.6],
[0.03, 0.3],
[0.04, 0.0],
]
)
mock_loadtxt.return_value = uniform
tab = PairTab(filename="dummy_path", rcut=0.04)
expected = copy.deepcopy(tab.serialize())

mock_loadtxt.return_value = np.array(
[
[0.00, 1.0],
[0.01, 0.8],
[0.02, 0.6],
[0.09, 0.3],
[0.16, 0.0],
]
)
with self.assertRaisesRegex(ValueError, "evenly spaced"):
tab.reinit(filename="dummy_path", rcut=0.16)

actual = tab.serialize()
for key in ("rmin", "rmax", "hh", "ntypes", "rcut", "nspline"):
self.assertEqual(actual[key], expected[key])
for key in ("vdata", "tab_info", "tab_data"):
np.testing.assert_array_equal(
actual["@variables"][key], expected["@variables"][key]
)

@patch("numpy.loadtxt")
def test_uniform_grid(self, mock_loadtxt) -> None:
mock_loadtxt.return_value = np.array(
[
[0.00, 1.0],
[0.01, 0.8],
[0.02, 0.6],
[0.03, 0.3],
[0.04, 0.0],
]
)
tab = PairTab(filename="dummy_path", rcut=0.04)
np.testing.assert_allclose(tab.hh, 0.01)

def test_uniform_grid_rounded_text_precision(self) -> None:
rr = np.linspace(0.0, 6.0, 1000)
ee = np.exp(-rr)
with tempfile.TemporaryDirectory() as tmpdir:
path = os.path.join(tmpdir, "table.txt")
np.savetxt(path, np.stack((rr, ee), axis=1), fmt="%.6f")
tab = PairTab(filename=path)
np.testing.assert_allclose(tab.hh, 6.0 / 999, rtol=1e-6)


if __name__ == "__main__":
unittest.main(warnings="ignore")
Loading