diff --git a/chemtools/topology/critical.py b/chemtools/topology/critical.py index 8afbf6d0..8ad00fdf 100644 --- a/chemtools/topology/critical.py +++ b/chemtools/topology/critical.py @@ -133,23 +133,31 @@ 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 _: + # Go to the next iteration. 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]): + 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) - 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) + # 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) diff --git a/chemtools/topology/test/test_critical.py b/chemtools/topology/test/test_critical.py index 9ff49461..038ca056 100644 --- a/chemtools/topology/test/test_critical.py +++ b/chemtools/topology/test/test_critical.py @@ -1,6 +1,28 @@ +# -*- 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 unittest import TestCase from chemtools.topology.critical import Topology, CriticalPoint @@ -11,7 +33,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 +62,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 +117,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 +165,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)