diff --git a/bfit/density.py b/bfit/density.py index ad92da5..04b0c05 100644 --- a/bfit/density.py +++ b/bfit/density.py @@ -204,7 +204,7 @@ def orbitals_cusp(self): return self._orbitals_cusp @staticmethod - def radial_slater_orbital(exponent, number, points): + def radial_slater_orbital(exponent, number, points, normalized=True): r""" Compute the radial component of Slater-type orbitals on the given points. @@ -227,6 +227,8 @@ def radial_slater_orbital(exponent, number, points): The principal quantum numbers :math:`n` of :math:`M` Slater orbitals. points : ndarray, (N,) The radial :math:`r` grid points. + normalized : bool + If true, adds the normalization constant :math:`N`. Returns ------- @@ -242,14 +244,19 @@ def radial_slater_orbital(exponent, number, points): """ if points.ndim != 1: raise ValueError("The argument point should be a 1D array.") - # compute norm & pre-factor - norm = np.power(2. * exponent, number) * np.sqrt((2. * exponent) / factorial(2. * number)) - pref = np.power(points, number - 1).T + # compute pre-factor + with np.errstate(divide='ignore'): + pref = np.power(points, number - 1, dtype=np.float64).T # compute slater orbitals - slater = norm.T * pref * np.exp(-exponent * points).T + slater = pref * np.exp(-exponent * points).T + # compute normalization + if normalized: + norm = np.power(2. * exponent, number) * \ + np.sqrt((2. * exponent) / factorial(2. * number)) + slater *= norm.T return slater - def phi_matrix(self, points, deriv=False): + def phi_matrix(self, points, deriv=None): r""" Compute the linear combination of Slater-type orbitals on the given points. @@ -271,8 +278,9 @@ def phi_matrix(self, points, deriv=False): ---------- points : ndarray, (N,) The radial grid points. - deriv : bool - If true, use the derivative of the slater-orbitals. + deriv : int, optional + If it is one, return the derivative of the slater-orbitals. + If it is two, return the second derivative of the slater-orbitals. Returns ------- @@ -287,12 +295,21 @@ def phi_matrix(self, points, deriv=False): zero instead. See "derivative_radial_slater_type_orbital". """ + if deriv is not None: + if deriv not in [1, 2]: + raise ValueError( + f"The deriv parameter {deriv} should be either one or two. Higher-order" + f"is not supported." + ) + # compute orbital composed of a linear combination of Slater phi_matrix = np.zeros((len(points), len(self.orbitals))) for index, orbital in enumerate(self.orbitals): exps, number = self.orbitals_exp[orbital[1]], self.basis_numbers[orbital[1]] - if deriv: - slater = self.derivative_radial_slater_type_orbital(exps, number, points) + if deriv == 1: + slater = self.first_derivative_radial_slater_type_orbital(exps, number, points) + elif deriv == 2: + slater = self.second_derivative_radial_slater_type_orbital(exps, number, points) else: slater = self.radial_slater_orbital(exps, number, points) phi_matrix[:, index] = np.dot(slater, self.orbitals_coeff[orbital]).ravel() @@ -359,9 +376,9 @@ def atomic_density(self, points, mode="total"): return dens @staticmethod - def derivative_radial_slater_type_orbital(exponent, number, points): + def first_derivative_radial_slater_type_orbital(exponent, number, points): r""" - Compute the derivative of the radial component of Slater-type orbital on the given points. + Compute the first derivative of the radial component of Slater-type orbital. The derivative of the Slater-type orbital is defined as: @@ -387,30 +404,120 @@ def derivative_radial_slater_type_orbital(exponent, number, points): Returns ------- slater : ndarray, (N, M) - The Slater-type orbitals evaluated on the :math:`N` grid points. + First derivative of Slater-type orbitals evaluated on the :math:`N` grid points. - Notes - ----- - - At r = 0, the derivative is undefined and this function returns zero instead. + """ + norm = np.power(2. * exponent, number) * np.sqrt((2. * exponent) / factorial(2. * number)) + # compute R(r) = N r^{n - 1} e^{-\alpha r} + slater = SlaterAtoms.radial_slater_orbital(exponent, number, points, normalized=True) + + # Derivative of R(r) w.r.t. r: + # n=1: R(r) = N e^{-\alpha r}, so dR(r)/dr = N (-\alpha) e^{-\alpha r} = -\alpha R(r) + # n!=1: dR(r)/dr = N (-\alpha) r^{n - 1} e^{-\alpha r} + N (n - 1) r^{n - 2} e^{-\alpha r} + + # compute N (-\alpha) e^{-\alpha r} part of derivative which exists for for all n + # ------------------------------------------------------------------------------- + deriv = np.zeros((len(points), number.shape[0])) + deriv -= exponent.T * slater + + # compute part of the derivative which only exists for n != 1 + # ----------------------------------------------------------- + # calculate the un-normalized Slater with n - 1; i.e., r^{n - 2} e^{-\alpha r} + slater_minus_one = SlaterAtoms.radial_slater_orbital( + exponent, number - 1, points, normalized=False + ) + # compute N (n - 1) r^{n - 2} e^{-\alpha r} for n != 1 + deriv_pref = norm.T * slater_minus_one * (number.T - 1.0) + i_numb_one = np.where(number == 1)[0] + deriv_pref[:, i_numb_one] = 0.0 + deriv += deriv_pref + return deriv + + @staticmethod + def second_derivative_radial_slater_type_orbital(exponent, number, points): + r""" + Compute the second derivative of the radial component of Slater-type orbital. + + The derivative of the Slater-type orbital is defined as: + + .. math:: + \frac{d^2 R(r)}{dr^2} = \bigg[ + \frac{(n-1)(n-2)}{r^2} - \frac{2\alpha (n-1)}{r} + \alpha^2 \bigg] + N r^{n-1} e^{- \alpha r}, + + where + :math:`n` is the principal quantum number of that orbital, + :math:`N` is the normalizing constant, + :math:`r` is the radial point, distance to the origin, and + :math:`\alpha` is the zeta exponent of that orbital. - References + Parameters ---------- - See wikipedia page on "Slater-Type orbitals". + exponent : ndarray, (M, 1) + The zeta exponents of Slater orbitals. + number : ndarray, (M, 1) + The principle quantum numbers of Slater orbitals. + points : ndarray, (N,) + The radial grid points. If points contain zero, then it is undefined at those + points and set to zero. + + Returns + ------- + slater : ndarray, (N, M) + Second derivative of Slater-type orbitals evaluated on the :math:`N` grid points. """ - slater = SlaterAtoms.radial_slater_orbital(exponent, number, points) - # Consider the case when dividing by zero. + norm = np.power(2. * exponent, number) * np.sqrt((2. * exponent) / factorial(2. * number)) + # compute R(r) = N r^{n - 1} e^{-\alpha r} + slater = SlaterAtoms.radial_slater_orbital(exponent, number, points, normalized=True) + # calculate the un-normalized Slater with n - 1 and n - 2 with np.errstate(divide='ignore'): - # derivative - deriv_pref = (number.T - 1.) / np.reshape(points, (points.shape[0], 1)) - exponent.T - deriv_pref[np.abs(points) < 1e-10, :] = 0.0 - deriv = deriv_pref * slater + slater_minus_one = SlaterAtoms.radial_slater_orbital( + exponent, number - 1, points, normalized=False + ) + slater_minus_two = SlaterAtoms.radial_slater_orbital( + exponent, number - 2, points, normalized=False + ) + + # General formula for the first derivative: + # dR(r)/dr = N (-\alpha) r^{n - 1} e^{-\alpha r} + N (n - 1) r^{n - 2} e^{-\alpha r} + + # compute N \alpha^2 e^{-\alpha r} part of derivative which exists for for all n + # ------------------------------------------------------------------------------ + # when n=1, R(r) = N e^{-\alpha r}, so d^2R(r)/dr^2 = N \alpha^2 e^{-\alpha r} + deriv = exponent.T**2.0 * slater + + # compute part of the derivative which only exists for n != 1 + # ----------------------------------------------------------- + # calculate 2 * N (n - 1) (-\alpha) r^{n - 2} e^{-\alpha r} + deriv_pref = 2.0 * exponent.T * norm.T * slater_minus_one * (number.T - 1.0) + i_numb_one = np.where(number == 1)[0] + deriv_pref[:, i_numb_one] = 0.0 + deriv -= deriv_pref + + # compute part of the derivative which only exists for n != 1 & n != 2 + # -------------------------------------------------------------------- + # calculate N (n - 1) (n - 2) r^{n - 3} e^{-\alpha r} + deriv_pref = norm.T * slater_minus_two * (number.T - 1.0) * (number.T - 2.0) + i_numb_one = np.where(number == 1)[0] + deriv_pref[:, i_numb_one] = 0.0 + deriv_pref[:, np.where(number == 2)[0]] = 0.0 + deriv += deriv_pref return deriv def positive_definite_kinetic_energy(self, points): r""" Positive definite or Lagrangian kinetic energy density. + .. math:: + K(r) = \sum_i n_i \bigg[ \bigg(\frac{d \sum_j N_j c^i_j R_{n_j}}{dr} \bigg)^2 + + \frac{l_i (l_i + 1}{r^2} \bigg(\sum_j N_j c^i_j R_{n_j} \bigg)^2 \bigg], + + where + :math:`n` is the principal quantum number of that orbital, + :math:`N` is the normalizing constant, andx + :math:`r` is the radial point, distance to the origin + Parameters ---------- points : ndarray,(N,) @@ -421,6 +528,10 @@ def positive_definite_kinetic_energy(self, points): energy : ndarray, (N,) The kinetic energy on the grid points. + Notes + ----- + - When :math:`n=1` and :math:`r=0`, then kinetic energy is either nan or infinity. + """ phi_matrix = np.zeros((len(points), len(self.orbitals))) angular = { @@ -430,18 +541,27 @@ def positive_definite_kinetic_energy(self, points): exps, numbers = self.orbitals_exp[orbital[1]], self.basis_numbers[orbital[1]] # Calculate del^2 of radial component # derivative of the radial component - deriv_radial = self.derivative_radial_slater_type_orbital(exps, numbers, points) + deriv_radial = self.first_derivative_radial_slater_type_orbital(exps, numbers, points) phi_matrix[:, index] = np.ravel(np.dot(deriv_radial, self.orbitals_coeff[orbital])**2.0) # Calculate del^2 of spherical component - slater = SlaterAtoms.radial_slater_orbital(exps, numbers, points) - with np.errstate(divide='ignore'): - gradient_sph = ( - angular[orbital[1]] * - np.ravel(np.dot(slater, self.orbitals_coeff[orbital]))**2.0 / - points**2.0 - ) - gradient_sph[np.abs(points) < 1e-10] = 0.0 + norm = np.power(2. * exps, numbers) * np.sqrt((2. * exps) / factorial(2. * numbers)) + # Take unnormalized slater with number n-1, this is needed to remove divide by r^2 + slater_minus_one = SlaterAtoms.radial_slater_orbital( + exps, numbers - 1, points, normalized=False + ) + deriv_pref = norm.T * slater_minus_one + # When r=0 and n = 1, then the derivative is undefined and this returns infinity or nan + i_r_zero = np.where(np.abs(points) == 0.0)[0] + i_numb_one = np.where(numbers[0] == 1)[0] + indices = np.array([[x, y] for x in i_r_zero for y in i_numb_one]) + if len(indices) != 0: # if-statement needed to remove numpy warning using list + deriv_pref[indices] = np.inf + # Sum the slater orbital multipled with their coefficients + gradient_sph = ( + angular[orbital[1]] * + np.ravel(np.dot(deriv_pref, self.orbitals_coeff[orbital]))**2.0 + ) phi_matrix[:, index] += gradient_sph orb_occs = self.orbitals_occupation @@ -463,6 +583,70 @@ def derivative_density(self, points): The derivative of atomic density on the grid points. """ - factor = self.phi_matrix(points) * self.phi_matrix(points, deriv=True) + factor = self.phi_matrix(points) * self.phi_matrix(points, deriv=1) derivative = np.dot(2. * factor, self.orbitals_occupation).ravel() / (4 * np.pi) return derivative + + def laplacian_of_atomic_density(self, points): + r""" + Return the Laplacian of the atomic density on a set of points. + + The Laplacian in radial coordinates only is: + + .. math:: + \Delta f = \frac{1}{r^2}\frac{\partial }{\partial r}\bigg(r^2 \frac{df}{dr} \bigg). + + Letting f be the atomic density we have the following equation: + + .. math:: + \Delta f = 2 \bigg[\sum n_i \big( \frac{d \phi_i}{dr}^2 + + \frac{2 \phi \frac{d\phi_i}{dr}}{r} + + \phi_i \frac{d^2 \phi_i}{dr^2} \big) + \bigg], + + where + :math:`\phi_i` is the ith molecular orbital with occupation number :math:`n_i`. + + Parameters + ---------- + points : ndarray(N,) + The :math:`N` radial grid points. + + Returns + ------- + laplacian : ndarray(N,) + The Laplacian of the atomic density on the grid points. + + """ + molecular_orb = self.phi_matrix(points) + molecular_orb_deriv = self.phi_matrix(points, deriv=1) + molecular_orb_sec_deriv = self.phi_matrix(points, deriv=2) + + # phi_i * phi_i^\prime / r + with np.errstate(divide='ignore'): + # Absorb phi_i / r together + phi_i_r = np.zeros((len(points), len(self.orbitals))) + for index, orbital in enumerate(self.orbitals): + exps, numbers = self.orbitals_exp[orbital[1]], self.basis_numbers[orbital[1]] + # Calculate slater divided by r + # Unnormalized slater with number n-1, this is needed to remove divide by r + norm = np.power(2. * exps, numbers) * np.sqrt((2. * exps) / factorial(2. * numbers)) + slater_minus_one = SlaterAtoms.radial_slater_orbital( + exps, numbers - 1, points, normalized=False + ) + slater_r = norm.T * slater_minus_one + # When r=0 and n = 1, then slater/r is infinity. + i_r_zero = np.where(np.abs(points) == 0.0)[0] + i_numb_one = np.where(numbers[0] == 1)[0] + indices = np.array([[x, y] for x in i_r_zero for y in i_numb_one]) + if len(indices) != 0: # if-statement needed to remove numpy warning using list + slater_r[indices] = np.inf + phi_i_r[:, index] += np.ravel(np.dot(slater_r, self.orbitals_coeff[orbital])) + + factor = 2.0 * phi_i_r * molecular_orb_deriv + factor = ( + molecular_orb_deriv**2.0 + factor + molecular_orb * molecular_orb_sec_deriv + ) + # Multiply by occupation number to get molecular orbitals. + laplacian = np.dot(2.0 * factor, self.orbitals_occupation).ravel() / (4 * np.pi) + return laplacian diff --git a/bfit/test/test_density.py b/bfit/test/test_density.py index 1b76a7d..6790216 100644 --- a/bfit/test/test_density.py +++ b/bfit/test/test_density.py @@ -25,7 +25,7 @@ from bfit.density import SlaterAtoms from bfit.grid import ClenshawRadialGrid import numpy as np -from numpy.testing import assert_almost_equal, assert_equal, assert_raises +from numpy.testing import assert_allclose, assert_almost_equal, assert_equal, assert_raises import scipy @@ -35,6 +35,8 @@ def slater(e, n, r, derivative=False): slater = norm * np.power(r, n - 1) * np.exp(-e * r) if derivative: + if n == 1: + return - e * slater return (((n - 1) / r) - e) * slater return slater @@ -49,6 +51,12 @@ def test_slater_type_orbital_be(): # check values of a single orbital at r=2.0 orbital = be.radial_slater_orbital(np.array([[0.821620]]), np.array([[2]]), np.array([2])) assert_almost_equal(orbital, slater(0.821620, 2, 2.0), decimal=6) + # check values of a single orbital at r=0.0 and n=2 + orbital = be.radial_slater_orbital(np.array([[0.821620]]), np.array([[2]]), np.array([0.0])) + assert_almost_equal(orbital, slater(0.821620, 2, 0.0), decimal=6) + # check values of a single orbital at r=0.0 and n=1 + orbital = be.radial_slater_orbital(np.array([[0.821620]]), np.array([[1]]), np.array([0.0])) + assert_almost_equal(orbital, slater(0.821620, 1, 0.0), decimal=6) # check value of tow orbitals at r=1.0 & r=2.0 exps, nums = np.array([[12.683501], [0.821620]]), np.array([[1], [2]]) orbitals = be.radial_slater_orbital(exps, nums, np.array([1., 2.])) @@ -62,20 +70,21 @@ def test_derivative_slater_type_orbital_be(): # load Be atomic wave function be = SlaterAtoms("Be") # check values of a single orbital at r=1.0 - orbital = be.derivative_radial_slater_type_orbital(np.array([[12.683501]]), - np.array([[1]]), np.array([1])) + orbital = be.first_derivative_radial_slater_type_orbital(np.array([[12.683501]]), + np.array([[1]]), np.array([1])) assert_almost_equal(orbital, slater(12.683501, 1, 1.0, derivative=True), decimal=6) # check values of a single orbital at r=2.0 - orbital = be.derivative_radial_slater_type_orbital(np.array([[0.821620]]), np.array([[2]]), - np.array([2])) + orbital = be.first_derivative_radial_slater_type_orbital(np.array([[0.821620]]), np.array([[2]]), + np.array([2])) assert_almost_equal(orbital, slater(0.821620, 2, 2.0, derivative=True), decimal=6) # check value of single orbital at r = 0.0 - orbital = be.derivative_radial_slater_type_orbital(np.array([[0.821620]]), np.array([[2]]), - np.array([0.])) - assert_almost_equal(orbital, 0.0, decimal=6) + orbital = be.first_derivative_radial_slater_type_orbital(np.array([[0.821620]]), np.array([[2]]), + np.array([0.])) + actual = np.power(2. * 0.821620, 2) * np.sqrt((2. * 0.821620) / (4 * 3 * 2)) + assert_almost_equal(orbital, actual, decimal=6) # check value of tow orbitals at r=1.0 & r=2.0 exps, nums = np.array([[12.683501], [0.821620]]), np.array([[1], [2]]) - orbitals = be.derivative_radial_slater_type_orbital(exps, nums, np.array([1., 2.])) + orbitals = be.first_derivative_radial_slater_type_orbital(exps, nums, np.array([1., 2.])) expected = np.array([[slater(exps[0, 0], nums[0, 0], 1., True), slater(exps[1, 0], nums[1, 0], 1., True)], [slater(exps[0, 0], nums[0, 0], 2., True), @@ -83,26 +92,47 @@ def test_derivative_slater_type_orbital_be(): assert_almost_equal(orbitals, expected, decimal=6) +def test_second_derivative_of_slater_type_orbital_against_finite_difference(): + r"""Test second derivative of radial Slater-type orbital against finite-difference.""" + # load Be atomic wave function + be = SlaterAtoms("Be") + + # check values of a single orbital at various points + for r in np.arange(0.0, 10., 0.1): + exponent, number, pt = np.array([[12.683501]]), np.array([[1]]), np.array([r]) + actual = be.second_derivative_radial_slater_type_orbital(exponent, number, pt) + + first_deriv = lambda x: slater(exponent[0, 0], number[0, 0], x, derivative=True) + desired = (first_deriv(r + 1e-8) - first_deriv(r)) / 1e-8 + assert_almost_equal(actual, desired, decimal=3) + + def test_positive_definite_kinetic_energy_he(): r"""Test integral of kinetic energy density of helium against actual value.""" # load he atomic wave function he = SlaterAtoms("he") # compute density on an equally distant grid - grid = ClenshawRadialGrid(4, 30000, 35000) + grid = ClenshawRadialGrid(4, 30000, 35000, include_origin=False) energ = he.positive_definite_kinetic_energy(grid.points) integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0) assert np.all(np.abs(integral - he.kinetic_energy) < 1e-5) + # Test at r=0, that it returns nan + energ = he.positive_definite_kinetic_energy(np.array([0.])) + assert np.all(np.isnan(energ[0])) def test_positive_definite_kinetic_energy_li(): r"""Test integral of kinetic energy density of lithium against actual value.""" # load be atomic wave function - be = SlaterAtoms("li") + li = SlaterAtoms("li") # compute density on an equally distant grid - grid = ClenshawRadialGrid(3, 20000, 35000) - energ = be.positive_definite_kinetic_energy(grid.points) + grid = ClenshawRadialGrid(3, 20000, 35000, include_origin=False) + energ = li.positive_definite_kinetic_energy(grid.points) integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0) - assert np.all(np.abs(integral - be.kinetic_energy) < 1e-5) + assert np.all(np.abs(integral - li.kinetic_energy) < 1e-5) + # Test at r=0, that it returns nan + energ = li.positive_definite_kinetic_energy(np.array([0.])) + assert np.all(np.isnan(energ[0])) def test_positive_definite_kinetic_energy_c(): @@ -110,7 +140,7 @@ def test_positive_definite_kinetic_energy_c(): # load be atomic wave function be = SlaterAtoms("c") # compute density on an equally distant grid - grid = ClenshawRadialGrid(6, 20000, 35000) + grid = ClenshawRadialGrid(6, 20000, 35000, include_origin=False, dtype=np.float64) energ = be.positive_definite_kinetic_energy(grid.points) integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0) assert np.all(np.abs(integral - be.kinetic_energy) < 1e-5) @@ -121,7 +151,7 @@ def test_positive_definite_kinetic_energy_p(): # load be atomic wave function be = SlaterAtoms("p") # compute density on ClenshawCurtis Grid - grid = ClenshawRadialGrid(15, 20000, 35000) + grid = ClenshawRadialGrid(15, 20000, 35000, include_origin=False, dtype=np.longlong) energ = be.positive_definite_kinetic_energy(grid.points) integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0) assert np.all(np.abs(integral - be.kinetic_energy) < 1e-3) @@ -132,7 +162,7 @@ def test_positive_definite_kinetic_energy_ag(): # load c atomic wave function adens = SlaterAtoms("ag") # compute density on ClenshawCurtis Grid - grid = ClenshawRadialGrid(47, 20000, 35000) + grid = ClenshawRadialGrid(47, 20000, 35000, include_origin=False, dtype=np.longlong) energ = adens.positive_definite_kinetic_energy(grid.points) integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0) assert np.all(np.abs(integral - adens.kinetic_energy) < 1e-3) @@ -143,7 +173,7 @@ def test_phi_derivative_lcao_b(): # load Be atomic wave function b = SlaterAtoms("b") # check the values of the phi_matrix at point 1.0 - phi_matrix = b.phi_matrix(np.array([1]), deriv=True) + phi_matrix = b.phi_matrix(np.array([1]), deriv=1) def _slater_deriv(r): # compute expected value of 1S @@ -173,6 +203,17 @@ def _slater_deriv(r): assert_almost_equal(phi_matrix[1, :], _slater_deriv(2.), decimal=4) assert_almost_equal(phi_matrix[2, :], _slater_deriv(3.), decimal=4) + # Test second derivative of phi at point 1.0, 2.0, 3.0 + pts = np.arange(0.5, 5., 0.1) + actual = b.phi_matrix(pts, deriv=2) + for i, r in enumerate(pts): + desired = (np.array(_slater_deriv(r + 1e-8)) - np.array(_slater_deriv(r))) / 1e-8 + assert_almost_equal(actual[i], desired, decimal=3) + + # Test raises error if value is not correct. + assert_raises(ValueError, b.phi_matrix, points=pts, deriv=3) + assert_raises(ValueError, b.phi_matrix, points=pts, deriv=b) + def test_coeff_matrix_be(): r"""Test the coefficients of Beryllium is correctly parsed.""" @@ -371,51 +412,68 @@ def test_atomic_density_heavy_rn(): def test_kinetic_energy_cation_anion_c(): r"""Test integral of kinetic energy density of cation/anion of carbon.""" c = SlaterAtoms("c", cation=True) - grid = ClenshawRadialGrid(6, 20000, 35000) + grid = ClenshawRadialGrid(6, 20000, 35000, include_origin=False) energ = c.positive_definite_kinetic_energy(grid.points) integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0) assert_almost_equal(integral, c.kinetic_energy, decimal=6) c = SlaterAtoms("c", anion=True) - grid = ClenshawRadialGrid(7, 20000, 35000) + grid = ClenshawRadialGrid(7, 20000, 35000, include_origin=False) energ = c.positive_definite_kinetic_energy(grid.points) integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0) assert_almost_equal(integral, c.kinetic_energy, decimal=5) +def test_derivative_electron_density_h(): + r"""Test derivative of atomic density of hydrogen.""" + h = SlaterAtoms("h") + pts = np.random.uniform(0, 1, size=(100,)) + actual = h.derivative_density(pts) + finite_diff = (h.atomic_density(pts + 1e-8) - h.atomic_density(pts)) / 1e-8 + assert_allclose(actual, finite_diff) + + def test_derivative_electron_density_c(): r"""Test derivative of atomic density of cation of carbon.""" c = SlaterAtoms("c", cation=True) - eps = 1e-10 - grid = np.array([0.1, 0.1 + eps, 0.5, 0.5 + eps]) - dens = c.atomic_density(grid) - actual = c.derivative_density(np.array([0.1, 0.5])) - desired_0 = (dens[1] - dens[0]) / eps - desired_1 = (dens[3] - dens[2]) / eps - assert_almost_equal(actual, np.array([desired_0, desired_1]), decimal=4) + pts = np.random.uniform(0, 2, size=(100,)) + actual = c.derivative_density(pts) + finite_diff = (c.atomic_density(pts + 1e-8) - c.atomic_density(pts)) / 1e-8 + assert_allclose(actual, finite_diff, atol=1e-4) def test_derivative_electron_density_cr(): r"""Test derivative of atomic density of chromium.""" cr = SlaterAtoms("cr") - eps = 2.0e-8 - grid = np.array([0.1 - eps, 0.1, 0.1 + eps, 1. - eps, 1., 1. + eps]) - dens = cr.atomic_density(grid) - actual = cr.derivative_density(np.array([0.1, 1.])) - desired_0 = (dens[2] - dens[0]) / (2. * eps) - desired_1 = (dens[5] - dens[3]) / (2. * eps) - assert_almost_equal(actual, np.array([desired_0, desired_1]), decimal=4) + pts = np.random.uniform(0, 5, size=(100,)) + actual = cr.derivative_density(pts) + finite_diff = (cr.atomic_density(pts + 1e-8) - cr.atomic_density(pts)) / 1e-8 + assert_allclose(actual, finite_diff, rtol=1e-6, atol=1e-4) def test_kinetic_energy_heavy_element_ce(): r"""Test integral of kinetic energy of cesium.""" c = SlaterAtoms("ce") - grid = ClenshawRadialGrid(55, 10000, 30000) + grid = ClenshawRadialGrid(55, 10000, 30000, include_origin=False) energ = c.positive_definite_kinetic_energy(grid.points) integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0) assert_almost_equal(integral, c.kinetic_energy, decimal=2) +def test_laplacian_of_hydrogen_slater(): + r"""Test laplacian of hydrogen Slater.""" + # Slater is N e^(-r) + h = SlaterAtoms("h") + + pts = np.arange(0.0, 10., 0.1) + lap = h.laplacian_of_atomic_density(pts) + + # Analytic derive it + solution = np.exp(-2.0 * pts) * (pts - 1) * 16.0 / pts + solution /= (4.0 * np.pi) + assert_almost_equal(lap, solution, decimal=7) + + def test_raises(): r"""Test raises error of SlaterAtoms.""" assert_raises(TypeError, SlaterAtoms, 25)