Merge pull request #157 from theochem/resolve_issue_117

Check the magnitude of negative density values

.github/workflows/main.yml

@@ -6,7 +6,7 @@ on:
     branches:
       - master
       - website
-      
+
   pull_request:
     types: [opened, synchronize, reopened, closed]
     branches:

.github/workflows/pytest.yaml

@@ -20,7 +20,7 @@ jobs:
 
     steps:
       - uses: "actions/checkout@v3"
-        
+
       - name: Setup python for test ${{ matrix.py }}
         uses: actions/setup-python@v4
         with:

CODE_OF_CONDUCT.md

@@ -3,7 +3,7 @@
 - Communication is a two-party activity. If you’re not sure what someone means, then ask them; don’t assume.
 - If someone asks you to stop, then stop.
 - When critiquing the work of others, focus on ways to improve that work and keep in mind that your perspective is unique: what is obvious to you may not be obvious to others.
-- If a dispute can’t be resolved directly, contact the qc.devs administrators ([email protected]). If you do not feel comfortable contacting the qc.devs administrators directly, you can contact any of the individual administrators ([email protected], [email protected], [email protected], or [email protected]). 
+- If a dispute can’t be resolved directly, contact the qc.devs administrators ([email protected]). If you do not feel comfortable contacting the qc.devs administrators directly, you can contact any of the individual administrators ([email protected], [email protected], [email protected], or [email protected]).
 
 # Statement of Principles:
 The QCDevs team embraces the diversity of its past, current, and future members, recognizing that the diversity of our backgrounds, preferences, personalities, ambitions, priorities, knowledge, and skills provide unique and useful perspectives that enrich our collaboration and enhance the quality of the software we build together and helps us grow, together and separately, as developers and as people. This certainly includes diversity of ethnicity, race, body shape, personal appearance, age, gender identity and expression, sexual identity and orientation, culture, religion, political affiliation, socioeconomic status, physical and mental capabilities, and national origin, but also diversity of perspectives, skills, knowledge, educational background, and priorities. It is expected that members of QCDevs are not only open-minded towards diverse backgrounds and perspectives, but actively encourage the participation of others. We, both individually and collectively, have zero tolerance for disrespectful or inconsiderate behavior of any kind.
@@ -22,11 +22,11 @@ The QCDevs team strives to be welcoming and respectful, and we want to ensure th
 - be supportive,
 - be open to new ideas and approaches,
 - accept responsibility for their actions and the effects they have on others, sincerely, apologize without hesitation, and learn from their mistakes to become better,
-- communicate clearly, constructively, and politely during disagreements and, after any tense discussions are resolved, revisit them for lessons on how future interactions can be improved, 
+- communicate clearly, constructively, and politely during disagreements and, after any tense discussions are resolved, revisit them for lessons on how future interactions can be improved,
 - give others’ intentions the benefit of the doubt. Assume that others are trying to be constructive and helpful, because they probably are.
 - think about the long term. Just as you are building on others’ contributions, they will build on your contributions. Try to make it easier for them; think about the implications of your actions and decisions.
 - steadfastly strive to become a better colleague and developer,
- 
+
 Similarly, the following activities are always unacceptable:
 - sexualized language or imagery, and sexual attention or advances of any kind
 - trolling, insulting or derogatory comments, and personal or political attacks
@@ -37,13 +37,13 @@ Similarly, the following activities are always unacceptable:
 
 This isn't an exhaustive list of things that you can and cannot do. Rather, take it in the spirit in which it's intended - a guide to make it easier to communicate and participate in QCDevs. As a single overarching rule: do not focus on what is best for you, but focus on what is best for the QCDevs team as a whole. This does not mean that your personal success does not matter! It means that everyone’s success matters equally, because our team is strong only when all of us are thriving personally and professionally.
 
-## Enforcement: 
-If you believe someone is violating the code of conduct we ask that you report it to the QCDevs administrative team at [email protected]. However, if you believe anyone is in physical danger, please notify appropriate law enforcement first. If you do not feel comfortable contacting the QCDevs administrators collectively, you can contact any of the individual administrators ([email protected], [email protected], [email protected], or [email protected]). 
+## Enforcement:
+If you believe someone is violating the code of conduct we ask that you report it to the QCDevs administrative team at [email protected]. However, if you believe anyone is in physical danger, please notify appropriate law enforcement first. If you do not feel comfortable contacting the QCDevs administrators collectively, you can contact any of the individual administrators ([email protected], [email protected], [email protected], or [email protected]).
 
 All reports will be kept confidential whenever possible. The QCDevs Administrators will strive to protect the identity and safety of reporters. In some cases we may need to make a public statement of some form, in which case we will use the minimum details and identifying information necessary to protect our community. In rare cases, we may need to identify some of the people involved to comply with the law or protect other potential victims. In these cases, we will consult with the reporter to find out what their wishes are and take them into account. In all cases, we aim not to directly or indirectly identify reporters without their consent.
 
-In your report please include 
-- your name and contact information, 
+In your report please include
+- your name and contact information,
 - names (real, nicknames, and/or pseudonyms) of any individuals involved. If there were other witnesses besides you, please try to include them as well.
 - when and where the incident occurred. Please be as specific as possible.
 - your account of what occurred. If there is a publicly available record please include a link and/or screen shots.
@@ -54,11 +54,11 @@ In your report please include
 As quickly as possible, the QCDevs administrators will review the incident and determine whether an investigation is needed, including interviewing additional parties or witnesses. After we have determined what has happened, we’ll determine whether the behaviour contravenes the code of conduct. Actions taken may include one or more of the following:
 - nothing (for example, if we determine that no violation occurred).
 - a private reprimand to the individual(s) involved, explaining why the conduct was inappropriate. This is appropriate when the behavior is minor, potentially accidental, and there is no damage to the QCDevs team as a whole, nor to its principles.
-- a public reprimand. This is appropriate when the action compromises the cohesiveness and integrity of the QCDevs team as a whole. 
+- a public reprimand. This is appropriate when the action compromises the cohesiveness and integrity of the QCDevs team as a whole.
 - a request for a private or public apology.
 - a request to compose an accountability plan or to engage in mediation.
 - a temporary or permanent suspension from QCDevs.
 
 Once the QCDevs administrators have determined our final action, we will contact the original reporter to let them know what action (if any) we will be taking. We welcome feedback from the reporter on the appropriateness of our response, and we may (or may not) modulate our response based on that feedback.
 
-**Credits:** This document was written by Wil Adams, Paul Ayers, Fanwang Meng, Lian Pharoah, and Michael Richer based on previous codes of conduct in QCDevs projects (primarily from Farnaz Heidar-Zadeh and Toon Verstraelen) and the freeBSD code of conduct, the Python code of conduct, the Open-Source Contributor Code of Conduct, and resources provided by the Office of Equity and Inclusion at McMaster University. The current version of the code of conduct was edited, then approved, by the QCDevs administrators. 
+**Credits:** This document was written by Wil Adams, Paul Ayers, Fanwang Meng, Lian Pharoah, and Michael Richer based on previous codes of conduct in QCDevs projects (primarily from Farnaz Heidar-Zadeh and Toon Verstraelen) and the freeBSD code of conduct, the Python code of conduct, the Open-Source Contributor Code of Conduct, and resources provided by the Office of Equity and Inclusion at McMaster University. The current version of the code of conduct was edited, then approved, by the QCDevs administrators.

gbasis/contractions.py

@@ -519,14 +519,17 @@ def coord_type(self, coord_type):
                 - 'spherical'
         """
 
-
         if not isinstance(coord_type, str):
             raise TypeError("Coordinate type must be given as a string.")
         if coord_type not in ["c", "cartesian", "p", "spherical"]:
-            raise ValueError("`coord_type` is incorrectly specified. It must be either 'c' "
-                             "or 'cartesian' for Cartesian coordinates, or 'p' or 'spherical' "
-                             "for spherical coordinates.")
-        self._coord_type = {"c": "cartesian",
-                            "cartesian": "cartesian",
-                            "spherical": "spherical",
-                            "p": "spherical"}[coord_type]
+            raise ValueError(
+                "`coord_type` is incorrectly specified. It must be either 'c' "
+                "or 'cartesian' for Cartesian coordinates, or 'p' or 'spherical' "
+                "for spherical coordinates."
+            )
+        self._coord_type = {
+            "c": "cartesian",
+            "cartesian": "cartesian",
+            "spherical": "spherical",
+            "p": "spherical",
+        }[coord_type]

gbasis/evals/_deriv.py

@@ -143,7 +143,7 @@ def _eval_deriv_contractions(coords, orders, center, angmom_comps, alphas, prim_
 
 
 def _eval_first_second_order_deriv_contractions(
-        coords, orders, center, angmom_comps, alphas, prim_coeffs, norm
+    coords, orders, center, angmom_comps, alphas, prim_coeffs, norm
 ):
     """Return the evaluation of direct 1st and 2nd derivative orders of a Cartesian contraction.
 
@@ -195,7 +195,7 @@ def _eval_first_second_order_deriv_contractions(
 
     # Useful variables
     new_coords = coords.T - center[None, :].T
-    gauss = np.exp(-alphas[:, None, None] * (new_coords ** 2))
+    gauss = np.exp(-alphas[:, None, None] * (new_coords**2))
 
     # Filters derivative orders
     indices_noderiv = orders <= 0
@@ -216,17 +216,21 @@ def _eval_first_second_order_deriv_contractions(
 
     # Calling 1st and 2nd derivatives functions for different combination of orders
     if indices_first_deriv.any():
-        first_deriv = _first_derivative(new_coords, gauss, indices_first_deriv,
-                                        angmom_comps, alphas)
+        first_deriv = _first_derivative(
+            new_coords, gauss, indices_first_deriv, angmom_comps, alphas
+        )
         if indices_second_deriv.any():
-            second_deriv = _second_derivative(new_coords, gauss, indices_second_deriv,
-                                              angmom_comps, alphas)
+            second_deriv = _second_derivative(
+                new_coords, gauss, indices_second_deriv, angmom_comps, alphas
+            )
     elif indices_second_deriv.any():
-        second_deriv = _second_derivative(new_coords, gauss, indices_second_deriv,
-                                          angmom_comps, alphas)
+        second_deriv = _second_derivative(
+            new_coords, gauss, indices_second_deriv, angmom_comps, alphas
+        )
         if indices_first_deriv.any():
-            first_deriv = _first_derivative(new_coords, gauss, indices_first_deriv,
-                                            angmom_comps, alphas)
+            first_deriv = _first_derivative(
+                new_coords, gauss, indices_first_deriv, angmom_comps, alphas
+            )
     # Combining all the derivatives
     norm = norm.T[:, :, np.newaxis]
     output = np.tensordot(prim_coeffs, norm * zeroth_deriv * first_deriv * second_deriv, (0, 0))
@@ -279,14 +283,15 @@ def _first_derivative(center_coords, gauss, indices_first_deriv, angmom_comps, a
     # zero (and negative power is undefined).
     power_part_1[power_part_1 < 0] = 0
     part1 = first_coords[:, None, :] ** power_part_1[:, :, None]
-    part2 = (2 * alphas[:, None, None]) * (first_coords ** 2)
+    part2 = (2 * alphas[:, None, None]) * (first_coords**2)
     part2 = first_ang_comp[None, :, :, None] - part2[:, :, None, :]
     # NOTE: Using an array of ones with same shape as first_ang_comp to power part2_zero_ang_mom
     # variable in order to get the same shape as part2. This is done in order to make easier
     # to filter at the end for the angular components corresponding to n=0
     array_ones = np.ones(first_ang_comp.shape)
-    part2_n_0 = -2 * alphas[:, None, None, None] * (
-            first_coords[:, None, :] ** array_ones[:, :, None])
+    part2_n_0 = (
+        -2 * alphas[:, None, None, None] * (first_coords[:, None, :] ** array_ones[:, :, None])
+    )
     raw_first_deriv = part1 * part2
     # Substitute angular components n=0 with correct derivative
     raw_first_deriv[:, n_0_indices, :] = part2_n_0[:, n_0_indices, :]
@@ -342,13 +347,15 @@ def _second_derivative(center_coords, gauss, indices_second_deriv, angmom_comps,
     n_2_indices = second_ang_comp >= 2
 
     # angular momentum == 0
-    total_n_0 = ((4 * alphas[:, None, None] ** 2) * (second_coords ** 2)) - \
-                (2 * alphas[:, None, None])
+    total_n_0 = ((4 * alphas[:, None, None] ** 2) * (second_coords**2)) - (
+        2 * alphas[:, None, None]
+    )
     raw_second_deriv = total_n_0[:, :, None, :] ** array_ones[None, :, :, None]
     # angular momentum == 1
     if any(second_ang_comp[0] == 1):
-        total_n_1 = ((4 * alphas[:, None, None] ** 2) * (second_coords ** 3)) \
-                    - ((6 * alphas[:, None, None]) * second_coords)
+        total_n_1 = ((4 * alphas[:, None, None] ** 2) * (second_coords**3)) - (
+            (6 * alphas[:, None, None]) * second_coords
+        )
         total_n_1 = total_n_1[:, :, None, :] ** array_ones[None, :, :, None]
         # Substitute angular components n=1 with correct derivative
         raw_second_deriv[:, n_1_indices, :] = total_n_1[:, n_1_indices, :]
@@ -363,9 +370,12 @@ def _second_derivative(center_coords, gauss, indices_second_deriv, angmom_comps,
             # Calculating
             # ..math:: 4 \alpha^{2}\left(x-X_{A}\right)^{4}-
             # \alpha(4 n+2)\left(x-X_{A}\right)^{2}+n(n-1)
-            part2_1_n_2 = (4 * alphas[:, None, None] ** 2) * (second_coords ** 4)
-            part2_2_n_2 = alphas[:, None, None, None] * (4 * second_ang_comp[:, :, None] + 2) \
-                                                      * second_coords[:, None, :] ** 2
+            part2_1_n_2 = (4 * alphas[:, None, None] ** 2) * (second_coords**4)
+            part2_2_n_2 = (
+                alphas[:, None, None, None]
+                * (4 * second_ang_comp[:, :, None] + 2)
+                * second_coords[:, None, :] ** 2
+            )
             part2_3_n_2 = second_ang_comp * (second_ang_comp - 1)
             part2_n_2 = part2_1_n_2[:, :, None, :] - part2_2_n_2 + part2_3_n_2[None, :, :, None]
             total_n_2 = part1_n_2[None, :, :, :] * part2_n_2

gbasis/evals/density.py

@@ -62,7 +62,7 @@ def evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval):
     return np.sum(density, axis=0)
 
 
-def evaluate_density(one_density_matrix, basis, points, transform=None):
+def evaluate_density(one_density_matrix, basis, points, transform=None, threshold=1.0e-8):
     r"""Return the density of the given basis set at the given points.
 
     Parameters
@@ -84,6 +84,9 @@ def evaluate_density(one_density_matrix, basis, points, 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.
+    threshold : float, optional
+        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.
 
     Returns
     -------
@@ -92,8 +95,12 @@ def evaluate_density(one_density_matrix, basis, points, transform=None):
 
     """
     orb_eval = evaluate_basis(basis, points, transform=transform)
-    # Fix: # 117; to avoid small negative density values, the array is clipped
-    return evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval).clip(min=0.0)
+    output = evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval)
+    # Fix #117: check magnitude of small negative density values, then use clip to remove them
+    min_output = np.min(output)
+    if min_output < 0.0 and abs(min_output) > threshold:
+        raise ValueError(f"Found negative density <= {-threshold}, got {min_output}.")
+    return output.clip(min=0.0)
 
 
 def evaluate_deriv_reduced_density_matrix(
@@ -553,6 +560,7 @@ def evaluate_posdef_kinetic_energy_density(
     points,
     transform=None,
     deriv_type="general",
+    threshold=1.0e-8,
 ):
     r"""Return evaluations of positive definite kinetic energy density at the given points.
 
@@ -581,17 +589,19 @@ def evaluate_posdef_kinetic_energy_density(
         are evaluated.
         Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z`
         components.
-    transform : np.ndarray(K_orbs, K_cont)
+    transform : np.ndarray(K_orbs, K_cont), optional
         Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
         combinations of contractions (e.g. MO).
         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.
-    deriv_type : "general" or "direct"
+    deriv_type : "general" or "direct", optional
         Specification of derivative of contraction function in _deriv.py. "general" makes reference
         to general implementation of any order derivative function (_eval_deriv_contractions())
         and "direct" makes reference to specific implementation of first and second order
         derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()).
+    threshold : float, optional
+        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.
 
     Returns
     -------
@@ -611,7 +621,10 @@ def evaluate_posdef_kinetic_energy_density(
             transform=transform,
             deriv_type=deriv_type,
         )
-    # Fix: #117; to avoid small negative values, the array is clipped
+    # Fix #117: check magnitude of small negative density values, then use clip to remove them
+    min_output = np.min(output)
+    if min_output < 0.0 and abs(min_output) > threshold:
+        raise ValueError(f"Found negative density <= {-threshold}, got {min_output}.")
     return (0.5 * output).clip(min=0.0)