From 137abc3b1a8ccfa804b36628ff496e8d9efd1bdf Mon Sep 17 00:00:00 2001 From: Ali-Tehrani Date: Wed, 4 Aug 2021 12:04:45 -0400 Subject: [PATCH 1/9] Fix calculating neighbors in critcal.py Moving a vector to a new-direction requires the formula: "vector + direction". The old way wasn't calculating the direction but rather just adding two vectors --- chemtools/topology/critical.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/chemtools/topology/critical.py b/chemtools/topology/critical.py index 8afbf6d0..7e2ae7ce 100644 --- a/chemtools/topology/critical.py +++ b/chemtools/topology/critical.py @@ -124,7 +124,9 @@ def find_critical_points(self): # compute distance to 4 closest grid points dists, _ = self._kdtree.query(point, 4) # store coordinates of neighbouring polyhedron vertices surrounding the point - neighs[4 * index: 4 * index + 4, :] = point + np.max(dists) * self._neighbours + neighs[4 * index: 4 * index + 4, :] = point + np.max(dists) * ( + self._neighbours - point + ) # compute the gradient norm of points & surrounding vertices points_norm = np.linalg.norm(self.grad(self._kdtree.data), axis=-1) From d6dbf0f84c278b2e58ce58c7559ab1f533c701be Mon Sep 17 00:00:00 2001 From: Ali-Tehrani Date: Wed, 4 Aug 2021 12:07:41 -0400 Subject: [PATCH 2/9] Add nan check to critical point finder The New-ton raphson algorithm for finding the roots/critical points sometimes returns nan. I added a check here to the if statement. --- chemtools/topology/critical.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/chemtools/topology/critical.py b/chemtools/topology/critical.py index 7e2ae7ce..862d8f2f 100644 --- a/chemtools/topology/critical.py +++ b/chemtools/topology/critical.py @@ -135,8 +135,11 @@ def find_critical_points(self): for index, point in enumerate(self._kdtree.data): point_norm = points_norm[index] neigh_norm = neighs_norm[4 * index: 4 * index + 4] - # use central point as initial guess for critical point finding - if index < len(self._coords) or np.all(point_norm < neigh_norm): + # If the gradient at the point is less than the neighbours then run/refine it further + # with root_vector_func or enumerate all the points that are the centers (if centers + # aren't None). + if np.all(point_norm < neigh_norm) or \ + (self._coords is not None and index < len(self._coords)): try: coord = self._root_vector_func(point.copy()) except np.linalg.LinAlgError as _: From 4ee8e958901092575d1ae8b0bf5ea81deb28f570 Mon Sep 17 00:00:00 2001 From: Ali-Tehrani Date: Wed, 4 Aug 2021 12:08:26 -0400 Subject: [PATCH 3/9] Change the order of if-statements in critical algo First change is to make sure, if there is LinAlg error then to discontinue. If the density is zero, then discontinue if the gradient isn't zero, then discontinue. --- chemtools/topology/critical.py | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/chemtools/topology/critical.py b/chemtools/topology/critical.py index 862d8f2f..a6f889de 100644 --- a/chemtools/topology/critical.py +++ b/chemtools/topology/critical.py @@ -142,19 +142,23 @@ def find_critical_points(self): (self._coords is not None and index < len(self._coords)): try: coord = self._root_vector_func(point.copy()) + # add critical point if it is new and it doesn't contain any + if not ( + np.any([np.linalg.norm(coord - cp.coordinate) < 1.e-3 for cp in self.cps]) + and not np.any(np.isnan(coord)) + ): + dens = self.func(coord) + # skip critical point if its density value is zero + if np.abs(dens) > 1.e-4: + grad = self.grad(coord) + # Make sure gradient is zero, as it's a critical point. + if np.all(np.abs(grad) < 1e-4): + # compute rank & signature of critical point + eigenvals, eigenvecs = np.linalg.eigh(self.hess(coord)) + cp = CriticalPoint(coord, eigenvals, eigenvecs, 1e-4) + self._cps.setdefault((cp.rank[0], cp.signature[0]), []).append(cp) except np.linalg.LinAlgError as _: continue - # add critical point if it is new - if not np.any([np.linalg.norm(coord - cp.coordinate) < 1.e-3 for cp in self.cps]): - dens = self.func(coord) - grad = self.grad(coord) - # skip critical point if its dens & grad are zero - if abs(dens) < 1.e-4 and np.all(abs(grad) < 1.e-4): - continue - # compute rank & signature of critical point - eigenvals, eigenvecs = np.linalg.eigh(self.hess(coord)) - cp = CriticalPoint(coord, eigenvals, eigenvecs, 1e-4) - self._cps.setdefault((cp.rank[0], cp.signature[0]), []).append(cp) # check Poincare–Hopf equation if not self.poincare_hopf_equation: warnings.warn("Poincare–Hopf equation is not satisfied.", RuntimeWarning) @@ -178,6 +182,7 @@ def _root_vector_func(self, guess, maxiter=5000): guess = guess - np.dot(np.linalg.inv(hess), grad[:, np.newaxis]).flatten() norm = np.linalg.norm(grad, axis=-1) niter += 1 + # Sometimes guess returns [np.nan, np.nan, np.nan] return guess @staticmethod From cf933dfe5b733e30cc16580fe8ee343e0fc5fea1 Mon Sep 17 00:00:00 2001 From: Ali-Tehrani Date: Wed, 4 Aug 2021 12:09:21 -0400 Subject: [PATCH 4/9] Add tests to critical point finder --- chemtools/topology/test/test_critical.py | 163 +++++++++++------------ 1 file changed, 75 insertions(+), 88 deletions(-) diff --git a/chemtools/topology/test/test_critical.py b/chemtools/topology/test/test_critical.py index 9ff49461..a86cd0bc 100644 --- a/chemtools/topology/test/test_critical.py +++ b/chemtools/topology/test/test_critical.py @@ -1,6 +1,6 @@ """Test critical point finder.""" -""" + from unittest import TestCase from chemtools.topology.critical import Topology, CriticalPoint @@ -11,7 +11,6 @@ class TestCriticalPoints(TestCase): # Test critical point finder class. - def gauss_func(self, coors, centers=np.zeros(3), alphas=1): # Generate gaussian function value. return np.prod(np.exp((-alphas * (coors - centers) ** 2)), axis=-1) @@ -41,77 +40,36 @@ def gauss_deriv2(self, coors, centers=np.zeros(3), alphas=1): hm[j][i] = hm[i][j] return hm - def test_topo(self): + def test_instantiation_of_topology_class(self): # Test properly initiate topo instance. coors = np.array([[1, 1, 1], [-1, -1, -1]]) pts = np.random.rand(4, 3) - topo = Topology(coors, self.gauss_func, self.gauss_deriv, self.gauss_deriv2, pts) + topo = Topology(self.gauss_func, self.gauss_deriv, self.gauss_deriv2, pts, coords=coors) assert isinstance(topo, Topology) - def test_default_cube(self): - # Test default cube for points. - coors = np.array([[1, 1, 1], [-1, -1, -1]]) - topo = Topology(coors, self.gauss_func, self.gauss_deriv, self.gauss_deriv2) - assert topo._kdtree.data.shape == (60 * 60 * 60, 3) - - def test_construct_cage(self): - # Test construct cage among target points. - pts = Topology._construct_cage(np.array([0, 0, 0]), 1) - assert len(pts) == 4 - dis = np.linalg.norm(pts[:] - np.array([0, 0, 0]), axis=-1) - assert_allclose(dis, np.ones(4) * 4.89898, rtol=1e-5) - - def test_get_gradietn(self): - # Test get proper gradient. - # initiate topo obj. - coors = np.array([[1, 1, 1], [-1, -1, -1]]) + def test_center_of_gaussian_is_a_nuclear_critical_point(self): + r"""Test the center of a single Gaussian is a nuclear critical point.""" + coords = np.array([[0., 0., 0.]]) pts = np.random.rand(4, 3) - topo = Topology(coors, self.gauss_func, self.gauss_deriv, self.gauss_deriv2, pts) - # get cage points - pts = Topology._construct_cage(np.array([0.5, 0.5, 0.5]), 0.1) - g_pts = topo.get_gradient(pts) - # assert len(g_pts) == 4 - ref_g = self.gauss_deriv(pts) - assert_allclose(ref_g, g_pts) - - def test_add_critical_points(self): - # Test add critical points to class. - coors = np.array([[1, 1, 1], [-1, -1, -1]]) - pts = np.random.rand(4, 3) - topo = Topology(coors, self.gauss_func, self.gauss_deriv, self.gauss_deriv2, pts) - pt = CriticalPoint(np.random.rand(3), None, None) - # bond - ct_type = -1 - topo._add_critical_point(pt, ct_type) - assert_allclose(topo._bcp[0].point, pt.point, atol=1e-10) - # maxima - ct_type = -3 - topo._add_critical_point(pt, ct_type) - assert_allclose(topo._nna[0].point, pt.point, atol=1e-10) - # ring - ct_type = 1 - topo._add_critical_point(pt, ct_type) - assert_allclose(topo._rcp[0].point, pt.point, atol=1e-10) - # cage - ct_type = 3 - topo._add_critical_point(pt, ct_type) - assert_allclose(topo._ccp[0].point, pt.point, atol=1e-10) - - def test_is_coors_pt(self): - # Test same pts as nuclear position. - coors = np.array([[1, 1, 1], [-1, -1, -1]]) - pts = np.random.rand(4, 3) - topo = Topology(coors, self.gauss_func, self.gauss_deriv, self.gauss_deriv2, pts) - for i in coors: - result = topo._is_coors_pt(i) - assert result - pt = (np.random.rand(4, 3) + 0.1) / 2 - for i in pt: - result = topo._is_coors_pt(i) - assert not result - - def test_find_one_critical_pt(self): - # Test two atoms critical pts. + topo = Topology(self.gauss_func, self.gauss_deriv, self.gauss_deriv2, pts, coords) + topo.find_critical_points() + # Test that only one critical point is found. + assert len(topo.cps) == 1 + # Test critical point is center + assert_allclose(np.array([0., 0., 0.]), topo.cps[0].coordinate) + # Test it is a nuclear critical point. + assert len(topo.nna) == 1 + assert_allclose(np.array([0., 0., 0.]), topo.nna[0].coordinate) + + def test_find_critical_pts_of_sum_of_two_gaussians(self): + r""" + Test finding critical points of a sum of two Gaussian functions. + + There should be three critical points, two of them are the center of the Gaussians. + The third should be in-between them. + + """ + # The center of the two Gaussians and the alpha parameters atoms = np.array([[-2, -2, -2], [2, 2, 2]]) alf = 0.3 @@ -137,19 +95,32 @@ def fun_d2(coors): pts = np.vstack( (pts, np.array([0.05, 0.05, 0.05]), np.array([-0.05, -0.05, -0.05])) ) - tp_ins = Topology(atoms, fun_v, fun_d, fun_d2, pts) - tp_ins.find_critical_pts() - # one critical pt - assert len(tp_ins._found_ct) == 1 - assert len(tp_ins._bcp) == 1 - # one critical type bond - assert tp_ins._found_ct_type[0] == -1 - # one critical point at origin - assert_allclose(tp_ins._found_ct[0], [0, 0, 0], atol=1e-10) - - def test_find_triangle_critical_pt(self): - # Test three atom ring critical pts. - atoms = np.array([[-2, -2, 0], [2, -2, 0], [0, 1, 0]]) + tp_ins = Topology(fun_v, fun_d, fun_d2, pts, atoms) + tp_ins.find_critical_points() + # Test that three critical points are found + assert len(tp_ins.cps) == 3 + # Test that two of them are the centers + assert len(tp_ins.nna) == 2 + assert_allclose(tp_ins.nna[0].coordinate, atoms[0, :], rtol=1e-5) + assert_allclose(tp_ins.nna[1].coordinate, atoms[1, :], rtol=1e-5) + assert_allclose(tp_ins.cps[-2].coordinate, atoms[0, :], rtol=1e-5) + assert_allclose(tp_ins.cps[-1].coordinate, atoms[1, :], rtol=1e-5) + # Test the bond critical-point inbetween the two centers + assert len(tp_ins.bcp) == 1 + assert_allclose(tp_ins.bcp[0].coordinate, np.array([0., 0., 0.]), atol=1e-4) + # Poincare-Hopf is satisfied + assert tp_ins.poincare_hopf_equation + # Test all other types of critical points weren't found + assert len(tp_ins.ccp) == 0 + assert len(tp_ins.rcp) == 0 + + def test_find_critical_points_of_three_gaussians_in_triangle_format(self): + # Test three Gaussians centered at "atoms" with hyper-parameter alpha + # The three centers of the Gaussians should form a equaliteral triangle. + # Since it is a equaliteral triangle, there should 9 critical points, + # First three are the centers, the next three are inbetween the the tree Gaussians + # the final critical point is the center of the equilateral triangle. + atoms = np.array([[-2, -2, 0], [2, -2, 0], [0, 2.0 * (np.sqrt(3.0) - 1.0), 0]]) alf = 1 def fun_v(coors): @@ -172,11 +143,27 @@ def fun_d2(coors): + self.gauss_deriv2(coors, atoms[1], alphas=alf) + self.gauss_deriv2(coors, atoms[2], alphas=alf) ) - - tp_ins = Topology(atoms, fun_v, fun_d, fun_d2, extra=1) - tp_ins.find_critical_pts() - assert len(tp_ins._bcp) == 3 - assert len(tp_ins._rcp) == 1 - assert len(tp_ins._nna) == 0 - assert len(tp_ins._ccp) == 0 -""" + # Construct a Three-Dimensional Grid + one_dgrid = np.arange(-2.0, 2.0, 0.1) + one_dgridz = np.arange(-0.5, 0.5, 0.1) # Reduce computation time. + pts = np.vstack(np.meshgrid(one_dgrid, one_dgrid, one_dgridz)).reshape(3,-1).T + tp_ins = Topology(fun_v, fun_d, fun_d2, pts, atoms) + tp_ins.find_critical_points() + # Test that there are 7 critical points (c.p.) in total + assert len(tp_ins.cps) == 7 + # Test that there are 3 nuclear c.p., 3 bond c.p., and 1 ring cp.p at the center. + assert len(tp_ins.nna) == 3 + assert len(tp_ins.bcp) == 3 + assert len(tp_ins.rcp) == 1 + assert len(tp_ins.ccp) == 0 + # Test the Nuclear c.p. is the same as the center of Gaussians. + assert_allclose(tp_ins.nna[0].coordinate, atoms[0], rtol=1e-6) + assert_allclose(tp_ins.nna[1].coordinate, atoms[1], rtol=1e-6) + assert_allclose(tp_ins.nna[2].coordinate, atoms[2], rtol=1e-6) + # Test the bond c.p. is the same as the center between the center of Gaussians. + assert_allclose(tp_ins.bcp[0].coordinate, (atoms[0] + atoms[1]) / 2.0, atol=1e-3) + assert_allclose(tp_ins.bcp[1].coordinate, (atoms[1] + atoms[2]) / 2.0, atol=1e-3) + assert_allclose(tp_ins.bcp[2].coordinate, (atoms[2] + atoms[0]) / 2.0, atol=1e-3) + # Test the ring c.p. is the as the center of the equilateral triangle. + # Calculated using centroid of equilateral triangle online calculator. + assert_allclose(tp_ins.rcp[0].coordinate, np.array([0 , -0.845, 0]), atol=1e-3) From 40c44816751c409106ed316aba7e5ac352e25024 Mon Sep 17 00:00:00 2001 From: Ali-Tehrani Date: Wed, 4 Aug 2021 12:39:57 -0400 Subject: [PATCH 5/9] Fix comment error in critical.py --- chemtools/topology/critical.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/chemtools/topology/critical.py b/chemtools/topology/critical.py index a6f889de..023e629c 100644 --- a/chemtools/topology/critical.py +++ b/chemtools/topology/critical.py @@ -142,7 +142,7 @@ def find_critical_points(self): (self._coords is not None and index < len(self._coords)): try: coord = self._root_vector_func(point.copy()) - # add critical point if it is new and it doesn't contain any + # add critical point if it is new and it doesn't contain any nans if not ( np.any([np.linalg.norm(coord - cp.coordinate) < 1.e-3 for cp in self.cps]) and not np.any(np.isnan(coord)) From 448fcaef1aff0b6b160db0d38b452c8c1e33ec7e Mon Sep 17 00:00:00 2001 From: Ali-Tehrani Date: Wed, 4 Aug 2021 12:57:40 -0400 Subject: [PATCH 6/9] Clean up if statements in critical.py --- chemtools/topology/critical.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/chemtools/topology/critical.py b/chemtools/topology/critical.py index 023e629c..ba1132e5 100644 --- a/chemtools/topology/critical.py +++ b/chemtools/topology/critical.py @@ -124,9 +124,8 @@ def find_critical_points(self): # compute distance to 4 closest grid points dists, _ = self._kdtree.query(point, 4) # store coordinates of neighbouring polyhedron vertices surrounding the point - neighs[4 * index: 4 * index + 4, :] = point + np.max(dists) * ( - self._neighbours - point - ) + neighs[4 * index: 4 * index + 4, :] = point + np.max(dists) * \ + (self._neighbours - point) # compute the gradient norm of points & surrounding vertices points_norm = np.linalg.norm(self.grad(self._kdtree.data), axis=-1) @@ -143,15 +142,15 @@ def find_critical_points(self): try: coord = self._root_vector_func(point.copy()) # add critical point if it is new and it doesn't contain any nans - if not ( - np.any([np.linalg.norm(coord - cp.coordinate) < 1.e-3 for cp in self.cps]) - and not np.any(np.isnan(coord)) - ): + new_pt = not np.any([ + np.linalg.norm(coord - cp.coordinate) < 1.e-3 for cp in self.cps + ]) + if new_pt and not np.any(np.isnan(coord)): dens = self.func(coord) # skip critical point if its density value is zero if np.abs(dens) > 1.e-4: grad = self.grad(coord) - # Make sure gradient is zero, as it's a critical point. + # make sure gradient is zero, as it's a critical point. if np.all(np.abs(grad) < 1e-4): # compute rank & signature of critical point eigenvals, eigenvecs = np.linalg.eigh(self.hess(coord)) From ca3b08a133854afa118d95e8cb80ba28eabcc262 Mon Sep 17 00:00:00 2001 From: Ali-Tehrani Date: Wed, 4 Aug 2021 12:58:02 -0400 Subject: [PATCH 7/9] Add header to test critical --- chemtools/topology/test/test_critical.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/chemtools/topology/test/test_critical.py b/chemtools/topology/test/test_critical.py index a86cd0bc..038ca056 100644 --- a/chemtools/topology/test/test_critical.py +++ b/chemtools/topology/test/test_critical.py @@ -1,3 +1,25 @@ +# -*- coding: utf-8 -*- +# ChemTools is a collection of interpretive chemical tools for +# analyzing outputs of the quantum chemistry calculations. +# +# Copyright (C) 2016-2019 The ChemTools Development Team +# +# This file is part of ChemTools. +# +# ChemTools is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# ChemTools is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, see +# +# -- """Test critical point finder.""" From 32aa9d472ca5853250682d7614d96033f6d9788f Mon Sep 17 00:00:00 2001 From: Ali-Tehrani Date: Wed, 4 Aug 2021 13:25:37 -0400 Subject: [PATCH 8/9] Revert "Fix calculating neighbors in critcal.py" This reverts commit 137abc3b1a8ccfa804b36628ff496e8d9efd1bdf. --- chemtools/topology/critical.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/chemtools/topology/critical.py b/chemtools/topology/critical.py index ba1132e5..75c71f2d 100644 --- a/chemtools/topology/critical.py +++ b/chemtools/topology/critical.py @@ -124,8 +124,7 @@ def find_critical_points(self): # compute distance to 4 closest grid points dists, _ = self._kdtree.query(point, 4) # store coordinates of neighbouring polyhedron vertices surrounding the point - neighs[4 * index: 4 * index + 4, :] = point + np.max(dists) * \ - (self._neighbours - point) + neighs[4 * index: 4 * index + 4, :] = point + np.max(dists) * self._neighbours # compute the gradient norm of points & surrounding vertices points_norm = np.linalg.norm(self.grad(self._kdtree.data), axis=-1) From 64f04050dfc91a1d04cc041a909f47632ff6932d Mon Sep 17 00:00:00 2001 From: Ali-Tehrani Date: Wed, 4 Aug 2021 16:23:13 -0400 Subject: [PATCH 9/9] Revert "Change the order of if-statements in critical algo" This reverts commit 4ee8e958901092575d1ae8b0bf5ea81deb28f570. --- chemtools/topology/critical.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/chemtools/topology/critical.py b/chemtools/topology/critical.py index 75c71f2d..8ad00fdf 100644 --- a/chemtools/topology/critical.py +++ b/chemtools/topology/critical.py @@ -140,23 +140,24 @@ def find_critical_points(self): (self._coords is not None and index < len(self._coords)): try: coord = self._root_vector_func(point.copy()) - # add critical point if it is new and it doesn't contain any nans - new_pt = not np.any([ - np.linalg.norm(coord - cp.coordinate) < 1.e-3 for cp in self.cps - ]) - if new_pt and not np.any(np.isnan(coord)): - dens = self.func(coord) - # skip critical point if its density value is zero - if np.abs(dens) > 1.e-4: - grad = self.grad(coord) - # make sure gradient is zero, as it's a critical point. - if np.all(np.abs(grad) < 1e-4): - # compute rank & signature of critical point - eigenvals, eigenvecs = np.linalg.eigh(self.hess(coord)) - cp = CriticalPoint(coord, eigenvals, eigenvecs, 1e-4) - self._cps.setdefault((cp.rank[0], cp.signature[0]), []).append(cp) except np.linalg.LinAlgError as _: + # Go to the next iteration. continue + # add critical point if it is new + new_pt = not np.any([ + np.linalg.norm(coord - cp.coordinate) < 1.e-3 for cp in self.cps + ]) + if new_pt and not np.any(np.isnan(coord)): + dens = self.func(coord) + # skip critical point if its dens is greater than zero + if abs(dens) > 1.e-4: + # make sure gradient is zero, as it's a critical point. + grad = self.grad(coord) + if np.all(abs(grad) < 1.e-4): + # compute rank & signature of critical point + eigenvals, eigenvecs = np.linalg.eigh(self.hess(coord)) + cp = CriticalPoint(coord, eigenvals, eigenvecs, 1e-4) + self._cps.setdefault((cp.rank[0], cp.signature[0]), []).append(cp) # check Poincare–Hopf equation if not self.poincare_hopf_equation: warnings.warn("Poincare–Hopf equation is not satisfied.", RuntimeWarning) @@ -180,7 +181,6 @@ def _root_vector_func(self, guess, maxiter=5000): guess = guess - np.dot(np.linalg.inv(hess), grad[:, np.newaxis]).flatten() norm = np.linalg.norm(grad, axis=-1) niter += 1 - # Sometimes guess returns [np.nan, np.nan, np.nan] return guess @staticmethod