Jastrow

Jastrow factor is the sum of homogeneous, isotropic electron–electron terms \(u(r_{ij})\), isotropic electron–nucleus terms \(\chi(r_{iI})\) centered on the nuclei, isotropic electron–electron–nucleus terms \(f(r_{ij}, r_{iI}, r_{jI})\), also centered on the nuclei:

\[J = \sum_{i>j}^{N_e} u(r_{ij}) + \sum_{I=1}^{N_I}\sum_{i=1}^{N_e} \chi(r_{iI}) + \sum_{I=1}^{N_I}\sum_{i>j}^{N_e} f(r_{ij}, r_{iI}, r_{jI})\]

Jastrow part of wavefunction is represented by the casino.Jastrow class.

It must be initialized from the configuration files:

from casino.readers import CasinoConfig
from casino.jastrow import Jastrow

config_path = <path to a directory containing input file>
config = CasinoConfig(config_path)
config.read()
jastrow = Jastrow(config)

To prevent code duplication, we need to prepare the necessary intermediate data:

import numpy as np

neu, ned = config.input.neu, config.input.ned
ne = neu + ned
r_e = np.random.uniform(-1, 1, ne * 3).reshape((ne, 3))
atom_positions = config.wfn.atom_positions
e_vectors = np.expand_dims(r_e, 1) - np.expand_dims(r_e, 0)
n_vectors = np.expand_dims(r_e, 0) - np.expand_dims(atom_positions, 1)
e_powers = jastrow.ee_powers(e_vectors)
n_powers = jastrow.en_powers(n_vectors)
r_ij = np.linalg.norm(e_vectors, axis=2)
r_iI = np.linalg.norm(n_vectors, axis=2)
ee_spin_mask = np.ones(shape=(ne, ne), dtype=int)
ee_spin_mask[:neu, :neu] = 0
ee_spin_mask[neu:, neu:] = 2
en_spin_mask = np.zeros(shape=(ne,), dtype=int)
en_spin_mask[neu:] = 1
C = jastrow.trunc

Summary of Methods

Jastrow class has a following methods:

Method	Output	Shape
u_term	\(u(r_{ij})\)	scalar
chi_term	\(\chi(r_{iI})\)	scalar
f_term	\(f(r_{ij}, r_{iI}, r_{jI})\)	scalar
u_term_gradient	\(\nabla u(r_{ij})\)	\((3N_e,)\)
chi_term_gradient	\(\nabla \chi(r_{iI})\)	\((3N_e,)\)
f_term_gradient	\(\nabla f(r_{ij}, r_{iI}, r_{jI})\)	\((3N_e,)\)
u_term_laplacian	\(\Delta u(r_{ij})\)	scalar
chi_term_laplacian	\(\Delta \chi(r_{iI})\)	scalar
f_term_laplacian	\(\Delta f(r_{ij}, r_{iI}, r_{jI})\)	scalar

u-term

\(u(r_{ij})\) term consists of a complete power expansion in electron-electron distances \(r_{ij}\):

\[u(r_{ij}) = (r_{ij} - L_u)^C\Theta(L_u - r_{ij})\sum_{l=0}^{N_u}\alpha_lr^l_{ij}\]

where \(\Theta\) is the Heaviside function. This term goes to zero at the cutoff length \(L_u\) with \(C - 1\) continuous derivatives at. In this expression C determines the behavior at the cutoff length. If \(C = 2\), the gradient of this term is continuous but the second derivative and hence the local energy is discontinuous, and if \(C = 3\) then both the gradient of this term and the local energy are continuous. Electron-electron Kato cusp conditions at \(r_{ij} = 0\) satisfied by constraint:

\[\alpha_1 = \left[\frac{\Gamma_{ij}}{(-L_u)^C} + \frac{\alpha_0C}{L_u}\right]\]

where \(\Gamma_{ij} = 1/2\) if electrons \(i\) and \(j\) have opposite spins and \(\Gamma_{ij} = 1/4\) if \(i\) and \(j\) have the same spin. For certain electron coordinates, \(u\) term can be obtained with casino.Jastrow.u_term() method:

jastrow.u_term(e_powers)

this is equivalent to (continues from):

from numpy.polynomial.polynomial import polyval
poly = polyval(r_ij, jastrow.u_parameters.T)
cutoff = np.minimum(r_ij - jastrow.u_cutoff, 0) ** C
np.sum(np.triu(cutoff * np.choose(ee_spin_mask, poly, mode='wrap'), 1))

chi-term

\(\chi(r_{iI})\) term consists of a complete power expansion in electron-nucleus distances \(r_{iI}\):

\[\chi(r_{iI}) = (r_{iI} - L_{\chi I})^C\Theta(L_{\chi I} - r_{iI})\sum_{m=0}^{N_\chi}\beta_mr^m_{iI}\]

where \(\Theta\) is the Heaviside function. This term goes to zero at the cutoff length \(L_{\chi I}\). The term involving the ionic charge \(Z_I\) enforces the electron–nucleus Kato cusp condition:

\[\beta_1 = \left[\frac{-Z_I}{(-L_{\chi I})^C} + \frac{\beta_{0I}C}{L_{\chi I}}\right]\]

if the Slater part of the wave function satisfies the electron–nucleus cusp condition, or if a non-divergent pseudopotential is used, then the Jastrow factor is required to be cuspless at the nuclei, i.e it should satisfy the above equation with \(Z_I = 0\)

For certain electron coordinates, \(\chi\) term can be obtained with casino.Jastrow.chi_term() method:

jastrow.chi_term(n_powers)

this is equivalent to (continues from):

from numpy.polynomial.polynomial import polyval
poly = polyval(r_iI, jastrow.chi_parameters[0].T)
cutoff = np.minimum(r_iI - jastrow.chi_cutoff, 0) ** C
np.sum(cutoff[0] * np.choose(en_spin_mask, poly, mode='wrap'))

f-term

\(f(r_{ij}, r_{iI}, r_{jI})\) term is the most general expansion of a function of \(r_{ij}\) , \(r_{iI}\) , and \(r_{jI}\) that is cuspless at the coalescence point and goes smoothly to zero when either \(r_{iI}\) or \(r_{jI}\) reach cutoff lengths:

\[f(r_{ij}, r_{iI}, r_{jI}) = (r_{iI} - L_{fI})^C(r_{jI} - L_{fI})^C \Theta(L_{fI} - r_{iI})\Theta(L_{fI} - r_{jI}) \sum_{l=0}^{N_{fI}^{eN}}\sum_{m=0}^{N_{fI}^{eN}}\sum_{n=0}^{N_{fI}^{ee}}\gamma_{lmnI}r_{iI}^lr_{jI}^mr_{ij}^n\]

To ensure electron–electron Kato cusp conditions folowing \(2N_{fI}^{eN} + 1\) constraints is applied:

\[\sum_{l,m}^{l+m=k}\gamma_{lm1I} = 0\]

and to ensure electron–nucleus Kato cusp conditions folowing \(N_{fI}^{eN} + N_{fI}^{ee} + 1\) constraints is applied:

\[\sum_{l,m}^{l+m=k'}(C\gamma_{0mnI} - L_{fI}\gamma_{1mnI}) = 0\]

If desired, there are \(N_{fI}^{ee}\) constraints imposed to prevent duplication of \(u\) term \((γ_{00nI} = 0 \ \forall n)\) and there are \(N_{fI}^{eI}\) constraints imposed to prevent duplication of \(\chi\) term \((γ_{l00I} = 0 \ \forall l)\) also the Jastrow factor to be symmetric under electron exchanges it is required that \(\gamma_{lmnI} = \gamma_{mlnI} \ \forall I, m, l, n\).

For certain electron coordinates, \(f\) term can be obtained with casino.Jastrow.f_term() method:

jastrow.f_term(e_powers, n_powers)

this is equivalent to (continues from):

from numpy.polynomial.polynomial import polyval3d
r_ijI = np.tile(r_iI[0], (ne, 1))
cutoff = np.minimum(r_iI - jastrow.f_cutoff, 0) ** C
poly = polyval3d(r_ijI, r_ijI.T, r_ij, jastrow.f_parameters[0].T)
np.sum(np.triu(np.outer(cutoff[0], cutoff[0]) * np.choose(ee_spin_mask, poly, mode='wrap'), 1))

u-term gradient

Considering that gradient of spherically symmetric function (in 3-D space) is:

\[\nabla f(r) = f'(r) \mathbf{\hat r}\]

There is only two non-zero terms of \(u(r_{ij})\) gradient, i.e. by \(i\)-th or \(j\)-th electron coordinates:

\[\nabla_{e_i} u(r_{ij}) = -\nabla_{e_j} u(r_{ij}) = (r_{ij} - L_u)^C\Theta(L_u - r_{ij})\mathbf{\hat r}_{ij}\sum_{l=0}^{N_u}\left(\frac{C}{r_{ij} - L_u} + \frac{l}{r_{ij}}\right)\alpha_lr^l_{ij}\]

where \(\mathbf{\hat r}_{ij}\) is the unit vector in the direction of the \(\mathbf{r}_{ij}\)

For certain electron coordinates, \(u\) term gradient can be obtained with casino.Jastrow.u_term_gradient() method:

jastrow.u_term_gradient(e_powers, e_vectors)

this is equivalent to (continues from):

import numpy as np
from numpy.polynomial.polynomial import polyval
L = jastrow.u_cutoff
l = np.arange(jastrow.u_parameters.shape[1])
cutoff = np.minimum(r_ij - L, 0) ** C
poly = polyval(r_ij, jastrow.u_parameters.T) * C / (r_ij - L)
poly += polyval(r_ij, (l * jastrow.u_parameters).T) / r_ij
g_ij = np.nan_to_num(cutoff * np.choose(ee_spin_mask, poly, mode='wrap') * e_vectors.T / r_ij)
np.sum(g_ij, axis=1).T.ravel()

chi-term gradient

Considering that gradient of spherically symmetric function (in 3-D space) is:

\[\nabla f(r) = f'(r) \mathbf{\hat r}\]

There is only one non-zero term of \(\chi(r_{iI})\) gradient, i.e. by \(i\)-th electron coordinates:

\[\nabla_{e_i} \chi(r_{iI}) = (r_{iI} - L_{\chi I})^C\Theta(L_{\chi I} - r_{iI})\mathbf{\hat r}_{iI}\sum_{m=0}^{N_\chi}\left(\frac{C}{r_{iI} - L_{\chi I}} + \frac{m}{r_{iI}}\right)\beta_mr^m_{iI}\]

where \(\mathbf{\hat r}_{iI}\) is the unit vector in the direction of the \(\mathbf{r}_{iI}\)

For certain electron coordinates, \(\chi\) term gradient can be obtained with casino.Jastrow.chi_term_gradient() method:

jastrow.chi_term_gradient(n_powers, n_vectors)

this is equivalent to (continues from):

from numpy.polynomial.polynomial import polyval
L = jastrow.chi_cutoff
cutoff = np.minimum(r_iI - L, 0) ** C
r_iI = np.linalg.norm(n_vectors, axis=2)
m = np.arange(jastrow.chi_parameters[0].shape[1])
poly = polyval(r_iI, jastrow.chi_parameters[0].T) * C / (r_iI[0] - L[0])
poly += polyval(r_iI, (m * jastrow.chi_parameters[0]).T) / r_iI[0]
(cutoff[0] * np.choose(en_spin_mask, poly, mode='wrap') * n_vectors[0].T / r_iI[0]).T.ravel()

f-term gradient

Considering that gradient of spherically symmetric function (in 3-D space) is:

\[\nabla f(r) = f'(r) \mathbf{\hat r}\]

There is only two non-zero terms of \(f(r_{ij}, r_{iI}, r_{jI})\) gradient, i.e. by \(i\)-th or \(j\)-th electron coordinates:

\[g_{ij} = \mathbf{\hat r}_{ij} \sum_{l=0}^{N_{fI}^{eN}}\sum_{m=0}^{N_{fI}^{eN}}\sum_{n=0}^{N_{fI}^{ee}}\left(\frac{n}{r_{ij}}\right)\gamma_{lmnI}r_{iI}^lr_{jI}^mr_{ij}^n\]

\[g_{iI} = \mathbf{\hat r}_{iI} \sum_{l=0}^{N_{fI}^{eN}}\sum_{m=0}^{N_{fI}^{eN}}\sum_{n=0}^{N_{fI}^{ee}}\left(\frac{C}{r_{iI} - L_{fI}} + \frac{l}{r_{iI}}\right)\gamma_{lmnI}r_{iI}^lr_{jI}^mr_{ij}^n\]

\[g_{jI} = \mathbf{\hat r}_{jI} \sum_{l=0}^{N_{fI}^{eN}}\sum_{m=0}^{N_{fI}^{eN}}\sum_{n=0}^{N_{fI}^{ee}}\left(\frac{C}{r_{jI} - L_{fI}} + \frac{m}{r_{jI}}\right)\gamma_{lmnI}r_{iI}^lr_{jI}^mr_{ij}^n\]

\[\nabla_{e_i} f(r_{ij}, r_{iI}, r_{jI}) = (r_{iI} - L_{fI})^C(r_{jI} - L_{fI})^C \Theta(L_{fI} - r_{iI})\Theta(L_{fI} - r_{jI})(g_{iI} + g_{ij})\]

\[\nabla_{e_j} f(r_{ij}, r_{iI}, r_{jI}) = (r_{iI} - L_{fI})^C(r_{jI} - L_{fI})^C \Theta(L_{fI} - r_{iI})\Theta(L_{fI} - r_{jI})(g_{jI} - g_{ij})\]

For certain electron coordinates, \(f\) term gradient can be obtained with casino.Jastrow.f_term_gradient() method:

jastrow.f_term_gradient(e_powers, n_powers, e_vectors, n_vectors)

this is equivalent to (continues from):

from numpy.polynomial.polynomial import polyval3d
n = np.expand_dims(np.arange(jastrow.f_parameters[0].shape[1]), axis=(1, 2))
m = np.expand_dims(np.arange(jastrow.f_parameters[0].shape[2]), axis=1)
l = np.arange(jastrow.f_parameters[0].shape[3])
L = jastrow.f_cutoff
cutoff = np.minimum(r_iI - L, 0) ** C
r_ijI = np.tile(r_iI[0], (ne, 1))
poly = polyval3d(r_ijI, r_ijI.T, r_ij, jastrow.f_parameters[0].T)
poly_l = polyval3d(r_ijI, r_ijI.T, r_ij, (l * jastrow.f_parameters[0]).T)
poly_m = polyval3d(r_ijI, r_ijI.T, r_ij, (m * jastrow.f_parameters[0]).T)
poly_n = polyval3d(r_ijI, r_ijI.T, r_ij, (n * jastrow.f_parameters[0]).T)

g_ijI = np.choose(ee_spin_mask, poly, mode='wrap') * C / (r_iI[0] - L[0])
g_ijI += np.choose(ee_spin_mask, poly_l, mode='wrap') / r_iI[0]
g_ijI = np.triu(g_ijI, 1) * np.expand_dims(n_vectors[0].T / r_iI[0], 1)

g_jiI = np.choose(ee_spin_mask, poly, mode='wrap').T * C / (r_iI[0] - L[0])
g_jiI += np.choose(ee_spin_mask, poly_m, mode='wrap').T / r_iI[0]
g_jiI = np.tril(g_jiI, -1) * np.expand_dims(n_vectors[0].T / r_iI[0], 1)

g_ij = np.nan_to_num(np.choose(ee_spin_mask, poly_n / r_ij, mode='wrap') * e_vectors.T / r_ij)

np.sum(np.outer(cutoff[0], cutoff[0]) * (g_ijI + g_jiI + g_ij), axis=1).T

u-term laplacian

Considering that Laplace operator of spherically symmetric function (in 3-D space) is:

\[\Delta f(r) = f''(r) + \frac{2}{r} f'(r)\]

There is only two non-zero terms of \(u(r_{ij})\) laplacian, i.e. by \(i\)-th or \(j\)-th electron coordinates:

\[\Delta_{e_i} u(r_{ij}) = \Delta_{e_j} u(r_{ij}) = (r_{ij} - L_u)^C\Theta(L_u - r_{ij}) \sum_{l=0}^{N_u}\left(\frac{C(C-1)}{(r_{ij} - L_u)^2} + \frac{2C(l+1)}{r_{ij}(r_{ij} + L_u)} + \frac{l(l+1)}{r_{ij}^2}\right)\alpha_lr^l_{ij}\]

For certain electron coordinates, \(u\) term laplacian can be obtained with casino.Jastrow.u_term_laplacian() method:

jastrow.u_term_laplacian(e_powers)

this is equivalent to (continues from):

from numpy.polynomial.polynomial import polyval
L = jastrow.u_cutoff
l = np.arange(jastrow.u_parameters.shape[1])
cutoff = np.minimum(r_ij - jastrow.u_cutoff, 0) ** C
poly =  polyval(r_ij, jastrow.u_parameters.T) * C * (C - 1) / (r_ij - L) ** 2
poly += 2 * polyval(r_ij, ((l + 1) * jastrow.u_parameters).T) * C / r_ij / (r_ij - L)
poly += polyval(r_ij, (l * (l + 1) * jastrow.u_parameters).T) / r_ij ** 2
2 * np.sum(np.triu(cutoff * np.choose(ee_spin_mask, poly, mode='wrap'), 1))

chi-term laplacian

Considering that Laplace operator of spherically symmetric function (in 3-D space) is:

\[\Delta f(r) = f''(r) + \frac{2}{r} f'(r)\]

then \(\chi(r_{iI})\) term laplacian:

\[\Delta_{e_i} \chi(r_{iI}) = (r_{iI} - L_{\chi I})^C\Theta(L_{\chi I} - r_{iI}) \sum_{l=0}^{N_\chi}\left(\frac{C(C-1)}{(r_{iI} - L_{\chi I})^2} + \frac{2C(m+1)}{r_{iI}(r_{iI} - L_{\chi I})} + \frac{m(m+1)}{r_{iI}^2}\right)\beta_mr^m_{iI}\]

For certain electron coordinates, \(\chi\) term laplacian can be obtained with casino.Jastrow.chi_term_laplacian() method:

jastrow.chi_term_laplacian(n_powers)

this is equivalent to (continues from):

from numpy.polynomial.polynomial import polyval
L = jastrow.chi_cutoff
m = np.arange(jastrow.chi_parameters[0].shape[1])
cutoff = np.minimum(r_iI - L, 0) ** C
poly = polyval(r_iI, jastrow.chi_parameters[0].T) * C * (C - 1) / (r_iI[0] - L[0]) ** 2
poly += 2 * polyval(r_iI, ((m + 1) * jastrow.chi_parameters[0]).T) * C / r_iI / (r_iI[0] - L[0])
poly += polyval(r_iI, (m * (m + 1) * jastrow.chi_parameters[0]).T) / r_iI ** 2
np.sum(cutoff[0] * np.choose(en_spin_mask, poly, mode='wrap'))

f-term laplacian

Considering that Laplace operator of spherically symmetric function (in 3-D space) is:

\[\Delta f(r) = f''(r) + \frac{2}{r} f'(r)\]

and \(f\) term is a product of two spherically symmetric functions \(g(r_{ij})\) and \(h(r_{iI})\) so using:

\[\Delta_{e_i}(g(r_{ij})h(r_{iI})) = g(r_{ij})\Delta_{e_i}h(r_{iI}) + 2\nabla_{e_i}g(r_{ij})\nabla_{e_i}h(r_{iI}) + g(r_{ij})\Delta_{e_i}h(r_{iI})\]

There is only two non-zero terms of \(f(r_{ij}, r_{iI}, r_{jI})\) laplacian, i.e. by \(i\)-th or \(j\)-th electron coordinates:

\[l_{iI} = \sum_{l=0}^{N_{fI}^{eN}}\sum_{m=0}^{N_{fI}^{eN}}\sum_{n=0}^{N_{fI}^{ee}} \left(\frac{C(C-1)}{(r_{iI} - L_{fI})^2} + \frac{2C(l+1)}{r_{iI}(r_{iI} - L_{fI})} + \frac{l(l+1)}{r_{iI}^2}\right) \gamma_{lmnI}r_{iI}^lr_{jI}^mr_{ij}^n\]

\[l_{jI} = \sum_{l=0}^{N_{fI}^{eN}}\sum_{m=0}^{N_{fI}^{eN}}\sum_{n=0}^{N_{fI}^{ee}} \left(\frac{C(C-1)}{(r_{jI} - L_{fI})^2} + \frac{2C(l+1)}{r_{jI}(r_{jI} - L_{fI})} + \frac{l(l+1)}{r_{jI}^2}\right) \gamma_{lmnI}r_{iI}^lr_{jI}^mr_{ij}^n\]

\[l_{dot,i} = \mathbf{\hat r}_{ij} \cdot \mathbf{\hat r}_{iI} \sum_{l=0}^{N_{fI}^{eN}}\sum_{m=0}^{N_{fI}^{eN}}\sum_{n=0}^{N_{fI}^{ee}} \left(\frac{C}{(r_{iI} - L_{fI})} + \frac{l}{r_{iI}}\right) \frac{n}{r_{ij}} \gamma_{lmnI}r_{iI}^lr_{jI}^mr_{ij}^n\]

\[l_{dot,j} = \mathbf{\hat r}_{ij} \cdot \mathbf{\hat r}_{jI} \sum_{l=0}^{N_{fI}^{eN}}\sum_{m=0}^{N_{fI}^{eN}}\sum_{n=0}^{N_{fI}^{ee}} \left(\frac{C}{(r_{jI} - L_{fI})} + \frac{l}{r_{jI}}\right) \frac{n}{r_{ij}} \gamma_{lmnI}r_{iI}^lr_{jI}^mr_{ij}^n\]

\[l_{ij} = \sum_{l=0}^{N_{fI}^{eN}}\sum_{m=0}^{N_{fI}^{eN}}\sum_{n=0}^{N_{fI}^{ee}} \frac{n(n+1)}{r_{ij}^2} \gamma_{lmnI}r_{iI}^lr_{jI}^mr_{ij}^n\]

\[\Delta_{e_i} f(r_{ij}, r_{iI}, r_{jI}) = (r_{iI} - L_{fI})^C(r_{jI} - L_{fI})^C \Theta(L_{fI} - r_{iI})\Theta(L_{fI} - r_{jI}) (l_{iI} + 2l_{dot,i} + l_{ij})\]

\[\Delta_{e_j} f(r_{ij}, r_{iI}, r_{jI}) = (r_{iI} - L_{fI})^C(r_{jI} - L_{fI})^C \Theta(L_{fI} - r_{iI})\Theta(L_{fI} - r_{jI}) (l_{jI} - 2l_{dot,j} + l_{ij})\]

For certain electron coordinates, \(f\) term laplacian can be obtained with casino.Jastrow.f_term_laplacian() method:

jastrow.f_term_laplacian(e_powers, n_powers, e_vectors, n_vectors)

this is equivalent to (continues from):

from numpy.polynomial.polynomial import polyval3d
n = np.expand_dims(np.arange(jastrow.f_parameters[0].shape[1]), axis=(1, 2))
m = np.expand_dims(np.arange(jastrow.f_parameters[0].shape[2]), axis=1)
l = np.arange(jastrow.f_parameters[0].shape[3])
L = jastrow.f_cutoff
cutoff = np.minimum(r_iI - L, 0) ** C
r_ijI = np.tile(r_iI[0], (ne, 1))
poly = polyval3d(r_ijI, r_ijI.T, r_ij, jastrow.f_parameters[0].T)
poly_l = polyval3d(r_ijI, r_ijI.T, r_ij, ((l + 1) * jastrow.f_parameters[0]).T)
poly_ll = polyval3d(r_ijI, r_ijI.T, r_ij, (l * (l + 1) * jastrow.f_parameters[0]).T)
poly_n = polyval3d(r_ijI, r_ijI.T, r_ij, (n * jastrow.f_parameters[0]).T)
poly_nn = polyval3d(r_ijI, r_ijI.T, r_ij, (n * (n + 1) * jastrow.f_parameters[0]).T)
poly_lm = polyval3d(r_ijI, r_ijI.T, r_ij, (l * m * jastrow.f_parameters[0]).T)
poly_ln = polyval3d(r_ijI, r_ijI.T, r_ij, (l * n * jastrow.f_parameters[0]).T)

l_1 = np.choose(ee_spin_mask, poly, mode='wrap') * C * (C - 1) / (r_iI[0] - L[0]) ** 2
l_1 += 2 * np.choose(ee_spin_mask, poly_l, mode='wrap') * C / (r_iI[0] - L[0]) / r_iI[0]
l_1 += np.choose(ee_spin_mask, poly_ll, mode='wrap') / r_iI[0] ** 2

l_dot = np.choose(ee_spin_mask, poly_n, mode='wrap') * C / (r_iI[0] - L[0]) / r_ij
l_dot += np.choose(ee_spin_mask, poly_ln, mode='wrap') / r_iI[0] / r_ij
l_dot *= -np.einsum('aij,ai->ji', (e_vectors.T / r_ij), (n_vectors[0].T / r_iI[0]))

l_2 = np.choose(ee_spin_mask, poly_nn, mode='wrap') / r_ij ** 2

np.sum(np.outer(cutoff[0], cutoff[0]) * np.nan_to_num(l_1 + 2 * l_dot + l_2))