grid.utils module

Utility Module.

grid.utils.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

pointsnp.ndarray(n, 3)

The (3-dimensional) Cartesian coordinates of points.

centernp.ndarray(3,), list, optional

Cartesian coordinates of the center of spherical coordinate system. If None, origin is used.

Returns

np.ndarray(N, 3)

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

grid.utils.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, deriv_theta, deriv_phifloat, float, float

Derivative wrt to spherical coordinates

r, theta, phifloat, float, float

The spherical points where the derivative is taken.

Returns

ndarray(3,) :

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

grid.utils.dipole_moment_of_molecule(grid, density: ndarray, coords: ndarray, charges: ndarray)

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

ndarray:

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

grid.utils.generate_derivative_real_spherical_harmonics(l_max: int, theta: ndarray, phi: ndarray)

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_maxint

Largest angular degree of the spherical harmonics.

thetanp.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.

phinp.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

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

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\).

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\).

grid.utils.generate_orders_horton_order(order: int, type_ord: str, dim: int = 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_ordstr, optional

The type of the order, it is either “cartesian”, “radial”, “pure” or “pure-radial”.

dimint, optional

The dimension of the orders for only Cartesian.

Returns

ndarray(3 * order , D)

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], ….]

grid.utils.generate_real_spherical_harmonics(l_max: int, theta: ndarray, phi: ndarray)

Compute the real spherical harmonics recursively up to a maximum angular degree \(l\). [1] [2]

\[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_maxint

Largest angular degree of the spherical harmonics.

thetanp.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.

phinp.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

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

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\).

Notes

• The associated Legendre polynomials are computed using the forward recursion: [2] \(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] (1,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.

grid.utils.generate_real_spherical_harmonics_scipy(l_max: int, theta: ndarray, phi: ndarray)

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_maxint

Largest angular degree of the spherical harmonics.

thetanp.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.

phinp.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

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

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\).

Notes

• Spherical harmonics in SciPy may return NaNs for large degrees due to numerical instability in associated Legendre functions. We tested up to \(l_{max} = 600\) without NaNs.

grid.utils.get_cov_radii(atnums, cov_type='bragg')

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

Parameters

atnumsint or np.ndarray

Atomic number(s) of the element(s) of interest.

cov_typestr, 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

np.ndarray

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

grid.utils.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.

grid.utils.solid_harmonics(l_max: int, sph_pts: ndarray)

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_maxint

Largest angular degree of the spherical harmonics.

sph_ptsndarray(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

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

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.