-
Notifications
You must be signed in to change notification settings - Fork 37
Adds get_atomic_property method (issue #216) #292
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: master
Are you sure you want to change the base?
Changes from 1 commit
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change | ||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
@@ -611,6 +611,49 @@ def __getitem__(self, index: int): | |||||||||||||||||||||||
| ) | ||||||||||||||||||||||||
| return self._atgrids[index] | ||||||||||||||||||||||||
|
|
||||||||||||||||||||||||
| def get_atomic_property(self, property_values): | ||||||||||||||||||||||||
| """Convert property at grid points to property per atom. | ||||||||||||||||||||||||
|
|
||||||||||||||||||||||||
| This method takes a property evaluated at all grid points and calculated how much | ||||||||||||||||||||||||
| of that propery belongs to each atom in the molecule. | ||||||||||||||||||||||||
|
|
||||||||||||||||||||||||
| Parameters | ||||||||||||||||||||||||
| ---------- | ||||||||||||||||||||||||
| property_values : np.narray | ||||||||||||||||||||||||
| Must have the same length as the number of points in the molecular grid. | ||||||||||||||||||||||||
| """ | ||||||||||||||||||||||||
|
|
||||||||||||||||||||||||
| # check if the input is valid | ||||||||||||||||||||||||
| if not isinstance(property_values, np.ndarray): | ||||||||||||||||||||||||
| raise TypeError( | ||||||||||||||||||||||||
| f"property_values must be of numpy array type, got {type(property_values)}" | ||||||||||||||||||||||||
| ) | ||||||||||||||||||||||||
|
|
||||||||||||||||||||||||
| if property_values.size != self.size: | ||||||||||||||||||||||||
| raise ValueError( | ||||||||||||||||||||||||
| "property_values is not the same size as grid. \n" | ||||||||||||||||||||||||
|
||||||||||||||||||||||||
| "property_values is not the same size as grid. \n" | |
| "property_values is not the same size as grid\n" |
Copilot
AI
Feb 8, 2026
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Input validation only checks property_values.size and not its shape. This allows e.g. shape (self.size, 1) to pass, which will broadcast against 1D weights and can produce incorrect results or massive intermediate arrays during multiplication. Please validate that property_values.shape == (self.size,) (mirroring Grid.integrate) or coerce with np.asarray(...).reshape(-1) and then enforce 1D length.
| if property_values.size != self.size: | |
| raise ValueError( | |
| "property_values is not the same size as grid. \n" | |
| f"property_values.size: {property_values.size}, self.size: {self.size}" | |
| # ensure a one-dimensional array to avoid unintended broadcasting | |
| property_values = np.asarray(property_values).reshape(-1) | |
| if property_values.size != self.size: | |
| raise ValueError( | |
| "property_values must be one-dimensional with the same length as the grid.\n" | |
| f"property_values.shape: {property_values.shape}, expected: ({self.size},)" |
Copilot
AI
Feb 8, 2026
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
self.weights already represents self.atweights * self.aim_weights (see MolGrid.__init__). Using self.weights[start:end] here would avoid recomputing the product and keep this method consistent if the weight definition changes. Optionally, you could precompute property_values * self.weights once and use a segmented reduction to avoid the Python loop.
| # finding property of each atom | |
| atomic_property[i] = np.sum( | |
| property_values[start_index:end_index] | |
| * self.aim_weights[start_index:end_index] | |
| * self.atweights[start_index:end_index] | |
| # finding property of each atom using precomputed molecular weights | |
| atomic_property[i] = np.sum( | |
| property_values[start_index:end_index] | |
| * self.weights[start_index:end_index] |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -810,6 +810,55 @@ def test_integrate_hirshfeld_weights_pair_1s(self): | |
| occupation = mg.integrate(fn) | ||
| assert_almost_equal(occupation, 2.5, decimal=5) | ||
|
|
||
| def test_get_atomic_property_constant(self): | ||
| """Test with constant property values""" | ||
| coords = np.array([[0, 0, -0.5], [0, 0, 0.5]]) | ||
| mg = MolGrid.from_preset( | ||
| atnums=np.array([1, 1]), | ||
| atcoords=coords, | ||
| preset="coarse", | ||
| aim_weights=BeckeWeights(order=3), | ||
| ) | ||
|
|
||
| property_values = np.ones(mg.size) | ||
| atomic_properties = mg.get_atomic_property(property_values) | ||
|
|
||
| assert atomic_properties.size == 2 | ||
| assert_allclose(np.sum(atomic_properties), mg.integrate(property_values)) | ||
|
|
||
| def test_get_atomic_property_wrong_size(self): | ||
| """Test error handling for wrong array size""" | ||
| coords = np.array([[0, 0, -0.5], [0, 0, 0.5]]) | ||
| mg = MolGrid.from_preset( | ||
| atnums=np.array([1, 1]), | ||
| atcoords=coords, | ||
| preset="coarse", | ||
| aim_weights=BeckeWeights(order=3), | ||
| ) | ||
|
|
||
| wrong_size = np.ones(mg.size + 1) | ||
| with pytest.raises(ValueError, match="property_values is not the same size as grid"): | ||
| mg.get_atomic_property(wrong_size) | ||
|
||
|
|
||
| def test_get_atomic_property_nonuniform(self): | ||
| """Test with non-uniform property(distance from origin)""" | ||
| coords = np.array([[0, 0, -0.5], [0, 0, 0.5]]) | ||
| mg = MolGrid.from_preset( | ||
| atnums=np.array([1, 1]), | ||
| atcoords=coords, | ||
| preset="coarse", | ||
| aim_weights=BeckeWeights(order=3), | ||
| ) | ||
|
|
||
| property_values = np.linalg.norm(mg.points, axis=1) | ||
| atomic_properties = mg.get_atomic_property(property_values) | ||
|
|
||
| # Due to symmetry, both atoms should have equal properties | ||
| assert atomic_properties.size == 2 | ||
|
|
||
| assert_allclose(atomic_properties[0], atomic_properties[1], rtol=1e-3) | ||
| assert_allclose(np.sum(atomic_properties), mg.integrate(property_values), rtol=1e-3) | ||
|
|
||
|
|
||
| def test_interpolation_with_gaussian_center(): | ||
| r"""Test interpolation with molecular grid of sum of two Gaussian examples.""" | ||
|
|
||
Uh oh!
There was an error while loading. Please reload this page.