grid.utils module

Utility Module.

convert_cart_to_sph(points, center=None)

Convert a set of points from cartesian to spherical coordinates.

The convention that \(\theta \in [-\pi, \pi]\) and \(\phi \in [0, \pi)\) is chosen.

Parameters

• points (np.ndarray(n, 3)) – The (3-dimensional) Cartesian coordinates of points.

• center (np.ndarray(3,), list, optional) – Cartesian coordinates of the center of spherical coordinate system. If None, origin is used.

Returns

Spherical coordinates of atoms respect to the center [radius \(r\), azumuthal \(\theta\), polar \(\phi\)]

Return type

np.ndarray(N, 3)

convert_derivative_from_spherical_to_cartesian(deriv_r, deriv_theta, deriv_phi, r, theta, phi)

Convert the derivative from spherical coordinates to Cartesian.

If \(r\) is zero, then the derivative wrt to theta and phi are set to zero. If \(\phi\) is zero, then the derivative wrt to theta is set to zero.

Parameters

• deriv_r (float, float, float) – Derivative wrt to spherical coordinates

• deriv_theta (float, float, float) – Derivative wrt to spherical coordinates

• deriv_phi (float, float, float) – Derivative wrt to spherical coordinates

• r (float, float, float) – The spherical points where the derivative is taken.

• theta (float, float, float) – The spherical points where the derivative is taken.

• phi (float, float, float) – The spherical points where the derivative is taken.

Returns

The derivative wrt to (x, y, z) Cartesian coordinates.

Return type

ndarray(3,)

dipole_moment_of_molecule(grid, density, coords, charges)

Calculate the dipole of a molecule.

This is defined as the observable form a wavefunction \(\Psi\):

\[\vec{\mu} = \int \Psi \hat{\mu} \Psi \vec{r}\]

which results in

\[\vec{\mu} = \sum_{i=1}^{N_{atoms}} Z_i (\vec{R_i} - \vec{R_c}) - \int ((\vec{r} - \vec{R_c})) \rho(\vec{r}) dr,\]

where \(N_{atoms}\) is the number of atoms, \(Z_i\) is the atomic charge of the ith atom, \(\vec{R_i}\) is the ith coordinate of the atom, \(\vec{R_c}\) is the center of mass of the molecule in atomic units and \(\rho\) is the electron density of the molecule.

Parameters

• grid (Grid) – The grid used to perform the integration.

• density (ndarray) – The electron density evaluated on points on grid object.

• coords (ndarray[M, 3]) – The coordinates of the atoms.

• charges (ndarray[M]) – The atomic charge of each atom.

Returns

Returns array of size three corresponding to the dipole moment of a molecule.

Return type

ndarray

generate_derivative_real_spherical_harmonics(l_max, theta, phi)

Generate derivatives of real spherical harmonics.

If \(\phi\) is zero, then the first component of the derivative wrt to \(phi\) is set to zero.

Parameters

• l_max (int) – Largest angular degree of the spherical harmonics.

• theta (np.ndarray(N,)) – Azimuthal angle \(\theta \in [0, 2\pi]\) that are being evaluated on. If this angle is outside of bounds, then periodicity is used.

• phi (np.ndarray(N,)) – Polar angle \(\phi \in [0, \pi]\) that are being evaluated on. If this angle is outside of bounds, then periodicity is used.

Returns

Derivative of spherical harmonics, (theta first, then phi) of all degrees up to \(l_{max}\) and orders \(m\) in Horton 2 order, i.e. \(m=0, 1, -1, \cdots, l, -l\).

Return type

ndarray(2, (l_max^2 + 1)^2, M)

Notes

• The derivative with respect to \(\phi\) is

\[\frac{\partial Y_l^m(\theta, \phi)}{\partial\phi} = -m Y_l^{-m}(\theta, \phi),\]

• The derivative with respect to \(\theta\) is

\[\frac{\partial Y_l^m(\theta, \phi)}{\partial\phi} = |m| \cot(\phi) \Re (Y_l^{|m|}(\theta, \phi)) + \sqrt{(n - |m|)(n + |m| + 1)} Re(e^{-i\theta} Y_l^{|m| + 1}(\theta, \phi))\]

for positive \(m\) and for negative \(m\), the real projection \(\Re\) is replaced with the imaginary projection \(\Im\).

generate_orders_horton_order(order, type_ord, dim=3)

Generate all orders from an integer \(l\).

For Cartesian, the orders are \((n_x, n_y, n_z)\) such that they sum to order. If dim=1,2, then it generates Cartesian orders \((n_x)\), \((n_x, n_y)\), respectively, such that they sum to order.

For radial, the orders is just the order \(l\).

For pure, the orders \((l, m)\) following the order \([(l, 0), (l, 1), (l, -1), \cdots, (l, l), (l, -l)]\).

Parameters

• order (int) – The order \(l\).

• type_ord (str, optional) – The type of the order, it is either “cartesian”, “radial”, “pure” or “pure-radial”.

• dim (int, optional) – The dimension of the orders for only Cartesian.

Returns

Each row is a list of D integers (e.g. \((n_x, n_y, n_z)\) or \((l, m)\)). The output is in Horton 2 order. e.g. order=2 it’s [[2, 0, 0], [1, 1, 0], [1, 0, 1], [0, 2, 0], ….]

Return type

ndarray(3 * order , D)

generate_real_spherical_harmonics(l_max, theta, phi)

Compute the real spherical harmonics recursively up to a maximum angular degree l.

\[Y_l^m(\theta, \phi) = \sqrt{\frac{(2l + 1) (l - m)!}{4\pi (l + m)!}} f(m, \theta) P_l^m(\cos(\phi)),\]

where \(l\) is the angular degree, \(m\) is the order and \(f(m, \theta) = \sqrt{2} \cos(m \theta)\) when \(m>0\) otherwise \(f(m, \theta) = \sqrt{2} \sin(m\theta)\) when \(m<0\), and equal to one when \(m= 0\). \(P_l^m\) is the associated Legendre polynomial without the Conway phase factor. The angle \(\theta \in [0, 2\pi]\) is the azimuthal angle and \(\phi \in [0, \pi]\) is the polar angle.

Parameters

• l_max (int) – Largest angular degree of the spherical harmonics.

• theta (np.ndarray(N,)) – Azimuthal angle \(\theta \in [0, 2\pi]\) that are being evaluated on. If this angle is outside of bounds, then periodicity is used.

• phi (np.ndarray(N,)) – Polar angle \(\phi \in [0, \pi]\) that are being evaluated on. If this angle is outside of bounds, then periodicity is used.

Returns

Value of real spherical harmonics of all orders \(m\),and degree \(l\) spherical harmonics. For each degree \(l\), the orders \(m\) are in Horton 2 order, i.e. \(m=0, 1, -1, 2, -2, \cdots, l, -l\).

Return type

ndarray((l_max + 1)**2, N)

Notes

• The associated Legendre polynomials are computed using the forward recursion: \(P_l^m(\phi) = \frac{2l + 1}{l - m + 1}\cos(\phi) P_{l-1}^m(x) - \frac{(l + m)}{l - m + 1} P_{l-1}^m(x)\), and \(P_l^l(\phi) = (2l + 1) \sin(\phi) P_{l-1}^{l-1}(\phi)\).

• For higher maximum degree \(l_{max} > 1900\) with double precision the computation of spherical harmonic will underflow within the range \(20^\circ \leq \phi \leq 160^\circ\). This code does not implement the modified recursion formulas in [2] and instead opts for higher precision defined by the computer system and np.longdouble.

• The mapping from \((l, m)\) to the correct row index in the spherical harmonic array is given by \((l + 1)^2 + 2 * m - 1\) if \(m > 0\), else it is \((l+1)^2 + 2 |m|\).

References

1

Colombo, Oscar L. Numerical methods for harmonic analysis on the sphere. Ohio State Univ Columbus Dept of Geodetic Science And Surveying, 1981.

2

Holmes, Simon A., and Will E. Featherstone. “A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions.” Journal of Geodesy 76.5 (2002): 279-299.

generate_real_spherical_harmonics_scipy(l_max, theta, phi)

Generate real spherical harmonics up to degree \(l\) and for all orders \(m\).

The real spherical harmonics are defined as a function \(Y_{lm} : L^2(S^2) \rightarrow \mathbb{R}\) such that

\[\begin{split}Y_{lm} = \begin{cases} \frac{i}{\sqrt{2}} (Y^m_l - (-1)^m Y_l^{-m} & \text{if } < m < 0 \\ Y_l^0 & \text{if } m = 0 \\ \frac{1}{\sqrt{2}} (Y^{-|m|}_{l} + (-1)^m Y_l^m) & \text{if } m > 0 \end{cases},\end{split}\]

where \(l \in \mathbb{N}\), \(m \in \{-l_{max}, \cdots, l_{max} \}\) and \(Y^m_l\) is the complex spherical harmonic.

Parameters

• l_max (int) – Largest angular degree of the spherical harmonics.

• theta (np.ndarray(N,)) – Azimuthal angle \(\theta \in [0, 2\pi]\) that are being evaluated on. If this angle is outside of bounds, then periodicity is used.

• phi (np.ndarray(N,)) – Polar angle \(\phi \in [0, \pi]\) that are being evaluated on. If this angle is outside of bounds, then periodicity is used.

Returns

Value of real spherical harmonics of all orders \(m\), and degree \(l\) spherical harmonics. For each degree \(l\), the orders \(m\) are in HORTON 2 order, i.e. \(m=0, 1, -1, 2, -2, \cdots, l, -l\).

Return type

ndarray((l_max + 1)**2, N)

Notes

• SciPy spherical harmonics is known (Jan 30, 2024) to give NaN when the degree is large, for our experience, when l >= 86.

get_cov_radii(atnums, cov_type='bragg')

Get the covalent radii for given atomic number(s).

Parameters

• atnums (int or np.ndarray) – Atomic number(s) of the element(s) of interest.

• cov_type (str, optional) – Type of covalent radii for elements. Possible values are: - “bragg”: Bragg-Slater empirically measured covalent radii. - “cambridge”: Covalent radii from analysis of the Cambridge Database. - “alvarez”: Covalent radii from https://doi.org/10.1039/B801115J.

Returns

Covalent radii of the desired atom(s) in atomic units.

Return type

np.ndarray

isotopic_masses = {1: 1.007825, 2: 4.002603, 3: 7.016004, 4: 9.012182, 5: 11.009305, 6: 12.0, 7: 14.003074, 8: 15.994915, 9: 18.998403, 10: 19.99244, 11: 22.98977, 12: 23.985042, 13: 26.981538, 14: 27.976927, 15: 30.973762, 16: 31.972071, 17: 34.968853, 18: 39.962383, 19: 38.963707, 20: 39.962591, 21: 44.95591, 22: 47.947947, 23: 50.943964, 24: 51.940512, 25: 54.93805, 26: 55.934942, 27: 58.9332, 28: 57.935348, 29: 62.929601, 30: 63.929147, 31: 68.925581, 32: 73.921178, 33: 74.921596, 34: 79.916522, 35: 78.918338, 36: 83.911507, 37: 84.911789, 38: 87.905614, 39: 88.905848, 40: 89.904704, 41: 92.906378, 42: 97.905408, 43: 97.907216, 44: 101.90435, 45: 102.905504, 46: 105.903483, 47: 106.905093, 48: 113.903358, 49: 114.903878, 50: 119.902197, 51: 120.903818, 52: 129.906223, 53: 126.904468, 54: 131.904154, 55: 132.905447, 56: 137.905241, 57: 138.906348, 58: 139.905434, 59: 140.907648, 60: 141.907719, 61: 144.912744, 62: 151.919728, 63: 152.921226, 64: 157.924101, 65: 158.925343, 66: 161.926795, 67: 164.930319, 68: 165.93029, 69: 168.934211, 70: 173.938858, 71: 174.940768, 72: 179.946549, 73: 183.950933, 74: 186.955751, 75: 186.955751, 76: 191.961479, 77: 192.962924, 78: 194.964774, 79: 196.966552, 80: 201.970626, 81: 204.974412, 82: 207.976636}

Obtained from theochem/horton/data/grids/tv-13.5.txt with the following preamble

These grids were created by Toon Verstraelen in July 2013 in a second effort to generate tuned grids based on a series of diatomic benchmarks computed with the QZVP basis of Ahlrichs and coworkers. They are constructed to give on average 5 digits of precision for a variety of molecular integrals, using the Becke scheme with k=3.

solid_harmonics(l_max, sph_pts)

Generate the solid harmonics from zero to a maximum angular degree.

\[R_l^m(r, \theta, \phi) = \sqrt{\frac{4\pi}{2l + 1}} r^l Y_l^m(\theta, \phi)\]

Parameters

• l_max (int) – Largest angular degree of the spherical harmonics.

• sph_pts (ndarray(M, 3)) – Three-dimensional points in spherical coordinates \((r, \theta, \phi)\), where \(r\geq 0, \theta \in [0, 2\pi]\) and \(\phi \in [0, \pi]\).

Returns

The solid harmonics from degree \(l=0=, \cdots, l_{max}\) and for all orders \(m\), ordered as \(m=0, 1, -1, 2, -2, \cdots, l, -l\) evaluated on \(M\) points.

Return type

ndarray((l_max + 1)^2, M)