Slater determinant

The Slater determinant component of the wavefunction is implemented in the casino.slater.Slater class. This class provides methods to compute the value, gradient, Laplacian, Hessian, and Tressian of a multi-determinant Slater wavefunction.

It must be initialized from the configuration files:

from casino.readers import CasinoConfig
from casino.slater import Slater

config_path = <path to a directory containing input file>
config = CasinoConfig(config_path)
config.read()
slater = Slater(config, cusp=None)

Summary of Methods

The Slater class has the following methods:

Method	Output	Shape
value_matrix	\(A^\uparrow, A^\downarrow\)	\((MO^\uparrow, N^\uparrow_e), (MO^\downarrow, N^\downarrow_e)\)
gradient_matrix	\(G^\uparrow, G^\downarrow\)	\((MO^\uparrow, N^\uparrow_e, 3), (MO^\downarrow, N^\downarrow_e, 3)\)
laplacian_matrix	\(L^\uparrow, L^\downarrow\)	\((MO^\uparrow, N^\uparrow_e), (MO^\downarrow, N^\downarrow_e)\)
hessian_matrix	\(H^\uparrow, H^\downarrow\)	\((MO^\uparrow, N^\uparrow_e, 3, 3), (MO^\downarrow, N^\downarrow_e, 3, 3)\)
tressian_matrix	\(T^\uparrow, T^\downarrow\)	\((MO^\uparrow, N^\uparrow_e, 3, 3, 3), (MO^\downarrow, N^\downarrow_e, 3, 3, 3)\)
value	\(\Psi(r)\)	\(scalar\)
gradient	\(\nabla \Psi(r)/\Psi(r)\)	\((3N_e,)\)
laplacian	\(\Delta \Psi(r)/\Psi(r)\)	\(scalar\)
hessian	\(\nabla^2 \Psi(r)/\Psi(r), \nabla \Psi(r)/\Psi(r)\)	\((3N_e, 3N_e), (3N_e,)\)
tressian	\(\nabla^3 \Psi(r)/\Psi(r), \nabla^2 \Psi(r)/\Psi(r), \nabla \Psi(r)/\Psi(r)\)	\((3N_e, 3N_e, 3N_e), (3N_e, 3N_e), (3N_e,)\)
tressian_dot	\(T_B, \nabla^2 \Psi(r)/\Psi(r), \nabla \Psi(r)/\Psi(r)\)	\((3N_e,), (3N_e, 3N_e), (3N_e,)\)

value matrix

In quantum chemistry, molecular orbitals (MOs) are normally expanded in a set of atom-centered basis functions, or localized atomic orbitals (AOs):

\[\phi_p(\mathbf{r}) = \sum_{\alpha}c_{\alpha p}\chi_\alpha(\mathbf{r}-\mathbf{R}_\alpha)\]

where \(\mathbf{r}=\{r_{1}...r_{N}\}\) are the coordinates of the N spin-up and spin-down electrons, \(\mathbf{R}_\alpha\) denotes the atomic position center of basis function \(\chi_\alpha\), and the expansion coefficients \(c_{\alpha p}\) are known as molecular orbital (MO) coefficients, also to avoid overflows and underflows a normalization coefficient is multiplied:

\[A_{ip} = \frac{1}{\sqrt[2N]{N^\uparrow! \, N^\downarrow!}} \phi_p(r_i)\]

In the multi-determinant case the \(p\) index should include occupied plus virtual MOs to span excited states. Therefore it is reasonable to define electrons \(\times\) (occupied + virtual MOs) matrix \(\mathcal{A}\). For a system described by a spin-independent Hamiltonian, the spatial and spin degrees of freedom are separable and we can split \(\mathcal{A}_{ip}\) into two: \(\mathcal{A}^\uparrow_{ip}\) for spin-up and \(\mathcal{A}^\downarrow_{ip}\) for spin-down electrons. For certain electron coordinates, the values of these matrices can be obtained with casino.slater.Slater.value_matrix() method:

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
n_vectors = np.expand_dims(r_e, 0) - np.expand_dims(atom_positions, 1)
A_up, A_down = slater.value_matrix(n_vectors)

note that the returned arrays are indexed as (orbital, electron), i.e. the transpose of \(A_{ip}\) defined above; since all expressions below are traces, they are invariant under this transposition.

the inverse matrix will be needed to calculate the gradient, laplacian, hessian and tressian (in the single-determinant case the matrices are square; in the multi-determinant case the rows must first be selected with permutation_up/permutation_down):

inv_A_up = np.linalg.inv(A_up)
inv_A_down = np.linalg.inv(A_down)

gradient matrix

Consider the gradient operator for \(i\)-th electron:

\[\nabla_{e_i} = \left[\frac{\partial}{\partial{x_i}}, \frac{\partial}{\partial{y_i}}, \frac{\partial}{\partial{z_i}}\right]\]

It is easy to check that:

\[\nabla_{e_i} A_{jp} = 0 \quad \text{if} \quad i \neq j\]

hence all non-zero values compose the matrix of vectors: \((x, y, z)\) indexed by \(a \in (x, y, z)\):

\[G_{ipa} = \nabla_{e_i} A_{ip}\]

In the multi-determinant case the \(p\) index should include occupied plus virtual MOs to span excited states. Therefore it is reasonable to define electrons \(\times\) (occupied + virtual MOs) matrix \(\mathcal{G}_{ip}\). For a system described by a spin-independent Hamiltonian, the spatial and spin degrees of freedom are separable and we can split \(\mathcal{G}_{ip}\) into two: \(\mathcal{G}^\uparrow_{ip}\) for spin-up and \(\mathcal{G}^\downarrow_{ip}\) for spin-down electrons. For certain electron coordinates, the values of these matrices can be obtained with casino.slater.Slater.gradient_matrix() method:

G_up, G_down = slater.gradient_matrix(n_vectors)

laplacian matrix

Consider the laplacian operator for \(i\)-th electron:

\[\Delta_{e_i} = \frac{\partial^2}{\partial{x_i}^2} + \frac{\partial^2}{\partial{y_i}^2} + \frac{\partial^2}{\partial{z_i}^2}\]

It is easy to check that:

\[\Delta_{e_i} A_{jp} = 0 \quad \text{if} \quad i \neq j\]

hence all non-zero values compose the matrix of scalars:

\[L_{ip} = \Delta_{e_i} A_{ip}\]

In the multi-determinant case the \(p\) index should include occupied plus virtual MOs to span excited states. Therefore it is reasonable to define electrons \(\times\) (occupied + virtual MOs) matrix \(\mathcal{L}_{ip}\). For a system described by a spin-independent Hamiltonian, the spatial and spin degrees of freedom are separable and we can split \(\mathcal{L}_{ip}\) into two: \(\mathcal{L}^\uparrow_{ip}\) for spin-up and \(\mathcal{L}^\downarrow_{ip}\) for spin-down electrons. For certain electron coordinates, the values of these matrices can be obtained with casino.slater.Slater.laplacian_matrix() method:

L_up, L_down = slater.laplacian_matrix(n_vectors)

hessian matrix

Consider the hessian operator for \(i\)-th electron:

\[\nabla_{e_i} \otimes \nabla_{e_i}\]

It is easy to check that:

\[(\nabla_{e_i} \otimes \nabla_{e_i}) A_{jp} = 0 \quad \text{if} \quad i \neq j\]

hence all non-zero values compose the matrix of hessians: \((x, y, z) \otimes (x, y, z)\) indexed by \(a,b \in (x, y, z)\):

\[H_{ipab} = (\nabla_{e_i} \otimes \nabla_{e_i}) A_{ip}\]

In the multi-determinant case the \(p\) index should include occupied plus virtual MOs to span excited states. Therefore it is reasonable to define electrons \(\times\) (occupied + virtual MOs) matrix \(\mathcal{H}_{ip}\). For a system described by a spin-independent Hamiltonian, the spatial and spin degrees of freedom are separable and we can split \(\mathcal{H}_{ip}\) into two: \(\mathcal{H}^\uparrow_{ip}\) for spin-up and \(\mathcal{H}^\downarrow_{ip}\) for spin-down electrons. For certain electron coordinates, the values of these matrices can be obtained with casino.slater.Slater.hessian_matrix() method:

H_up, H_down = slater.hessian_matrix(n_vectors)

tressian matrix

Consider the tressian operator for \(i\)-th electron:

\[\nabla_{e_i} \otimes \nabla_{e_i} \otimes \nabla_{e_i}\]

It is easy to check that:

\[(\nabla_{e_i} \otimes \nabla_{e_i} \otimes \nabla_{e_i}) A_{jp} = 0 \quad \text{if} \quad i \neq j\]

hence all non-zero values compose the matrix of tressians: \((x, y, z) \otimes (x, y, z) \otimes (x, y, z)\) indexed by \(a,b,c \in (x, y, z)\):

\[T_{ipabc} = (\nabla_{e_i} \otimes \nabla_{e_i} \otimes \nabla_{e_i}) A_{ip}\]

In the multi-determinant case the \(p\) index should include occupied plus virtual MOs to span excited states. Therefore it is reasonable to define electrons \(\times\) (occupied + virtual MOs) matrix \(\mathcal{T}_{ip}\). For a system described by a spin-independent Hamiltonian, the spatial and spin degrees of freedom are separable and we can split \(\mathcal{T}_{ip}\) into two: \(\mathcal{T}^\uparrow_{ip}\) for spin-up and \(\mathcal{T}^\downarrow_{ip}\) for spin-down electrons. For certain electron coordinates, the values of these matrices can be obtained with casino.slater.Slater.tressian_matrix() method:

T_up, T_down = slater.tressian_matrix(n_vectors)

value

Consider contribution of single Slater determinant:

\[\psi(\mathbf{r}) = \det(A)\]

we can get the value of multideterminant wavefunction:

\[\Psi(\mathbf{r}) = \sum_n c_n \psi(\mathbf{r})_n\]

and \(\mathbf{r}=\{r_{1}...r_{N}\}\) are the coordinates of the N spin-up and spin-down electrons.

For certain electron coordinates, the value can be obtained with casino.slater.Slater.value() method:

value = slater.value(n_vectors)

gradient

Consider Slater determinant gradient by \(i\)-th electron coordinates:

\[\frac{\nabla_{e_i} \psi(\mathbf{r})}{\psi(\mathbf{r})} = \left[ tr\left(A^{-1}\frac{\partial{A}}{\partial{x_i}}\right), tr\left(A^{-1}\frac{\partial{A}}{\partial{y_i}}\right), tr\left(A^{-1}\frac{\partial{A}}{\partial{z_i}}\right) \right] = tr(A^{-1} \nabla_{e_i} A)\]

to express the trace through sum using equality:

\[tr(AB) = \sum_{ij} a_{ij}b_{ji} = {a_i}^j {b_j}^i\]

notice that the \(\nabla_{e_i} A\) has the only one non-zero \(row_i(\nabla_{e_i} A) = row_i(G)\):

\[tr(A^{-1} \nabla_{e_i} A) = {(A^{-1})_i}^j {(\nabla_{e_i} A)_j}^{ia}\]

expand gradient vector over \(i\):

\[\frac{\nabla \psi(\mathbf{r})}{\psi(\mathbf{r})} = {(A^{-1})_i}^j G_{jia}\]

and get gradient of multideterminant wavefunction:

\[\nabla \Psi(\mathbf{r}) / \Psi(\mathbf{r}) = \sum_n c_n \nabla \psi(\mathbf{r})_n / \sum_n c_n \psi(\mathbf{r})_n\]

where \(\mathbf{r}=\{r_{1}...r_{N}\}\) are the coordinates of the N spin-up and spin-down electrons

For certain electron coordinates, the gradient vector can be obtained with casino.slater.Slater.gradient() method:

slater.gradient(n_vectors)

this is equivalent to (continues from):

G_up, G_down = slater.gradient_matrix(n_vectors)
tr_grad_u = np.einsum('ij,jia->ia', inv_A_up, G_up).reshape(neu * 3)
tr_grad_d = np.einsum('ij,jia->ia', inv_A_down, G_down).reshape(ned * 3)
np.concatenate((tr_grad_u, tr_grad_d))

laplacian

Consider Slater determinant laplacian by \(i\)-th electron coordinates:

\[\frac{\Delta_{e_i} \psi(\mathbf{r})}{\psi(\mathbf{r})} = tr\left(A^{-1}\frac{\partial^2{A}}{\partial{x_i}^2}\right) + tr\left(A^{-1}\frac{\partial^2{A}}{\partial{y_i}^2}\right) + tr\left(A^{-1}\frac{\partial^2{A}}{\partial{z_i}^2}\right) = tr(A^{-1} \Delta_{e_i} A)\]

to express the trace through sum using equality:

\[tr(AB) = \sum_{ij} a_{ij}b_{ji} = {a_i}^j {b_j}^i\]

notice that the \(\Delta_{e_i} A\) has the only one non-zero \(row_i(\Delta_{e_i} A) = row_i(L)\):

\[tr(A^{-1} \Delta_{e_i} A) = {(A^{-1})_i}^j {(\Delta_{e_i} A)_j}^i\]

sum laplacian over \(i\):

\[\frac{\Delta \psi(\mathbf{r})}{\psi(\mathbf{r})} = (A^{-1})_{ij} L^{ji}\]

and get laplacian of multideterminant wavefunction:

\[\Delta \Psi(\mathbf{r}) / \Psi(\mathbf{r}) = \sum_n c_n \Delta \psi(\mathbf{r})_n / \sum_n c_n \psi(\mathbf{r})_n\]

where \(\mathbf{r}=\{r_{1}...r_{N}\}\) are the coordinates of the N spin-up and spin-down electrons

For certain electron coordinates, the laplacian can be obtained with casino.slater.Slater.laplacian() method:

slater.laplacian(n_vectors)

this is equivalent to (continues from):

L_up, L_down = slater.laplacian_matrix(n_vectors)
lap_up = np.einsum('ij,ji', inv_A_up, L_up)
lap_down = np.einsum('ij,ji', inv_A_down, L_down)
lap_up + lap_down

hessian

Consider Slater determinant hessian by \(i\)-th and \(j\)-th electrons coordinates:

\[\frac{\nabla^2_{{e_i}{e_j}} \psi(\mathbf{r})}{\psi(\mathbf{r})} = tr(A^{-1} \nabla_{e_i} \nabla_{e_j} A - (A^{-1} \nabla_{e_i} A)(A^{-1} \nabla_{e_j} A)) + \frac{\nabla_{e_i} \psi(\mathbf{r})}{\psi(\mathbf{r})} \otimes \frac{\nabla_{e_j} \psi(\mathbf{r})}{\psi(\mathbf{r})}\]

to express the trace through sum using equality:

\[tr(AB) = \sum_{ij} a_{ij}b_{ji} = {a_i}^j {b_j}^i\]

notice that the \(\nabla_{e_i} A\) has the only one non-zero \(row_i(\nabla_{e_i} A) = row_i(G)\) and the \(\nabla_{e_i} \nabla_{e_i} A\) has only non-zero \(row_i(\nabla_{e_i} \nabla_{e_i} A) = row_i(H)\):

\[tr(A^{-1} \nabla_{e_i} \nabla_{e_j} A - (A^{-1} \nabla_{e_i} A)(A^{-1} \nabla_{e_j} A)) = {(A^{-1})_i}^j (\nabla_{e_i} {\nabla_{e_j} A)_j}^{iab} - {(A^{-1} \nabla_{e_i} A)_j}^{ia} {(A^{-1} \nabla_{e_j} A)_i}^{jb}\]

expand gradient vectors and hessian tensor over \(i\) and \(j\) (with Kronecker delta \(\delta_{ij}\)):

\[\begin{split}\frac{\nabla^2 \psi(\mathbf{r})}{\psi(\mathbf{r})} = \delta_{ij}{(A^{-1})_i}^j H_{jiab} - (A^{-1} G)_{jia} (A^{-1} G)_{ijb} + \frac{\nabla \psi(\mathbf{r})}{\psi(\mathbf{r})} \otimes \frac{\nabla \psi(\mathbf{r})}{\psi(\mathbf{r})} \\\end{split}\]

we can get hessian of multideterminant wavefunction:

\[\nabla^2 \Psi(\mathbf{r}) / \Psi(\mathbf{r}) = \sum_n c_n \nabla^2 \psi(\mathbf{r})_n / \sum_n c_n \psi(\mathbf{r})_n\]

where \(\mathbf{r}=\{r_{1}...r_{N}\}\) are the coordinates of the N spin-up and spin-down electrons

For certain electron coordinates, the hessian matrix can be obtained with casino.slater.Slater.hessian() method:

slater.hessian(n_vectors)[0]

this is equivalent to (continues from):

G_up, G_down = slater.gradient_matrix(n_vectors)
tr_grad_u = np.einsum('ij,jia->ia', inv_A_up, G_up).reshape(neu * 3)
tr_grad_d = np.einsum('ij,jib->ib', inv_A_down, G_down).reshape(ned * 3)
mul_grad_u = np.einsum('ij,jka->ika', inv_A_up, G_up)
mul_grad_d = np.einsum('ij,jkb->ikb', inv_A_down, G_down)
grad = np.concatenate((tr_grad_u, tr_grad_d))

H_up, H_down = slater.hessian_matrix(n_vectors)
tr_hess_u = np.einsum('ij,jiab->iab', inv_A_up, H_up)
tr_hess_d = np.einsum('ij,jiab->iab', inv_A_down, H_down)
hess_u = np.einsum('ij,iab->iajb', np.eye(neu), tr_hess_u)
hess_d = np.einsum('ij,iab->iajb', np.eye(ned), tr_hess_d)
hess_u -= np.einsum('ijb,jia->iajb', mul_grad_u, mul_grad_u)
hess_d -= np.einsum('ijb,jia->iajb', mul_grad_d, mul_grad_d)
hess = np.zeros((ne * 3, ne * 3))
hess[:neu * 3, :neu * 3] = hess_u.reshape(neu * 3, neu * 3)
hess[neu * 3:, neu * 3:] = hess_d.reshape(ned * 3, ned * 3)
hess += np.outer(grad, grad)

tressian

Consider Slater determinant tressian by \(i\)-th, \(j\)-th and \(k\)-th electrons coordinates:

\[\begin{split}\begin{align} & \frac{\nabla^3_{{e_i}{e_j}{e_k}} \psi(\mathbf{r})}{\psi(\mathbf{r})} = tr(A^{-1} \nabla_{e_i} \nabla_{e_j} \nabla_{e_k} A) - 2 \cdot \frac{\nabla_{e_i} \psi(\mathbf{r})}{\psi(\mathbf{r})} \otimes \frac{\nabla_{e_j} \psi(\mathbf{r})}{\psi(\mathbf{r})} \otimes \frac{\nabla_{e_k} \psi(\mathbf{r})}{\psi(\mathbf{r})} \\ & + \frac{\nabla^2_{{e_i}{e_j}} \psi(\mathbf{r})}{\psi(\mathbf{r})} \otimes \frac{\nabla_{e_k} \psi(\mathbf{r})}{\psi(\mathbf{r})} + \frac{\nabla^2_{{e_i}{e_k}} \psi(\mathbf{r})}{\psi(\mathbf{r})} \otimes \frac{\nabla_{e_j} \psi(\mathbf{r})}{\psi(\mathbf{r})} + \frac{\nabla^2_{{e_j}{e_k}} \psi(\mathbf{r})}{\psi(\mathbf{r})} \otimes \frac{\nabla_{e_i} \psi(\mathbf{r})}{\psi(\mathbf{r})} \\ & - tr((A^{-1} \nabla_{e_i} \nabla_{e_j} A)(A^{-1} \nabla_{e_k} A) + (A^{-1} \nabla_{e_i} \nabla_{e_k} A)(A^{-1} \nabla_{e_j} A) + (A^{-1} \nabla_{e_j} \nabla_{e_k} A)(A^{-1} \nabla_{e_i} A)) \\ & + tr((A^{-1} \nabla_{e_i} A)(A^{-1} \nabla_{e_j} A)(A^{-1} \nabla_{e_k} A)) + tr((A^{-1} \nabla_{e_k} A)(A^{-1} \nabla_{e_j} A)(A^{-1} \nabla_{e_i} A)) \end{align}\end{split}\]

to express the trace through sum using equalities:

\[tr(AB) = \sum_{ij} a_{ij}b_{ji} = {a_i}^j {b_j}^i\]

\[tr(ABC) = \sum_{ijk} a_{ij}b_{jk}c_{ki} = {a_i}^j {b_j}^k {c_k}^i\]

\[\begin{split}\begin{align} & tr(A^{-1} \nabla_{e_i} \nabla_{e_j} \nabla_{e_k} A) \\ & - tr((A^{-1} \nabla_{e_i} \nabla_{e_j} A)(A^{-1} \nabla_{e_k} A) + (A^{-1} \nabla_{e_i} \nabla_{e_k} A)(A^{-1} \nabla_{e_j} A) + (A^{-1} \nabla_{e_j} \nabla_{e_k} A)(A^{-1} \nabla_{e_i} A)) \\ & + tr((A^{-1} \nabla_{e_i} A)(A^{-1} \nabla_{e_j} A)(A^{-1} \nabla_{e_k} A) + (A^{-1} \nabla_{e_k} A)(A^{-1} \nabla_{e_j} A)(A^{-1} \nabla_{e_i} A)) \\ & = {(A^{-1})_i}^j {(\nabla_{e_i} \nabla_{e_j} \nabla_{e_k} A)_j}^{iabc} - {(A^{-1} \nabla_{e_i} \nabla_{e_j} A)_i}^{jab}{(A^{-1} \nabla_{e_k} A)_j}^{ic} \\ & - {(A^{-1} \nabla_{e_i} \nabla_{e_k} A)_i}^{jac}{(A^{-1} \nabla_{e_j} A)_j}^{ib} - {(A^{-1} \nabla_{e_j} \nabla_{e_k} A)_i}^{jbc}{(A^{-1} \nabla_{e_i} A)_j}^{ia} \\ & + {(A^{-1} \nabla_{e_i} A)_j}^{ia}{(A^{-1} \nabla_{e_j} A)_k}^{jb}{(A^{-1} \nabla_{e_k} A)_i}^{kc} + {(A^{-1} \nabla_{e_i} A)_k}^{ia}{(A^{-1} \nabla_{e_j} A)_i}^{jb}{(A^{-1} \nabla_{e_k} A)_j}^{kc} \end{align}\end{split}\]

notice that the \(\nabla_i A\) has only non-zero \(row_i(\nabla_i A) = row_i(G)\) and the \(\nabla_i \nabla_i A\) has only non-zero \(row_i(\nabla_i \nabla_i A) = row_i(H)\) and the \(\nabla_i \nabla_i \nabla_i A\) has only non-zero \(row_i(\nabla_i \nabla_i \nabla_i A) = row_i(T)\) and expand gradient vectors, hessian and tressian tensors over \(i\), \(j\), \(k\):

\[\begin{split}\begin{align} & \frac{\nabla^3 \psi(\mathbf{r})}{\psi(\mathbf{r})} = \delta_{ijk}{(A^{-1})_i}^jT_{jiabc} - 2 \cdot \frac{\nabla \psi(\mathbf{r})}{\psi(\mathbf{r})} \otimes \frac{\nabla \psi(\mathbf{r})}{\psi(\mathbf{r})} \otimes \frac{\nabla \psi(\mathbf{r})}{\psi(\mathbf{r})} \\ & + \left(\frac{\nabla^2 \psi(\mathbf{r})}{\psi(\mathbf{r})}\right)_{iajb} \left(\frac{\nabla \psi(\mathbf{r})}{\psi(\mathbf{r})}\right)_{kc} + \left(\frac{\nabla^2 \psi(\mathbf{r})}{\psi(\mathbf{r})}\right)_{iakc} \left(\frac{\nabla \psi(\mathbf{r})}{\psi(\mathbf{r})}\right)_{jb} + \left(\frac{\nabla^2 \psi(\mathbf{r})}{\psi(\mathbf{r})}\right)_{jbkc} \left(\frac{\nabla \psi(\mathbf{r})}{\psi(\mathbf{r})}\right)_{ia} \\ & - \delta_{ij}(A^{-1} H)_{ijab}(A^{-1} G)_{ijc} - \delta_{jk}(A^{-1} H)_{jkac}(A^{-1} G)_{jkb} - \delta_{ki}(A^{-1} G)_{kia}(A^{-1} H)_{kibc} \\ & + (A^{-1} G)_{jia}(A^{-1} G)_{kjb}(A^{-1} G)_{ikc} + (A^{-1} G)_{kia}(A^{-1} G)_{ijb}(A^{-1} G)_{jkc} \end{align}\end{split}\]

we can get tressian of multideterminant wavefunction:

\[\nabla^3 \Psi(\mathbf{r}) / \Psi(\mathbf{r}) = \sum_n c_n \nabla^3 \psi(\mathbf{r})_n / \sum_n c_n \psi(\mathbf{r})_n\]

where \(\mathbf{r}=\{r_{1}...r_{N}\}\) are the coordinates of the N spin-up and spin-down electrons

For certain electron coordinates, the tressian matrix can be obtained with casino.slater.Slater.tressian() method:

slater.tressian(n_vectors)[0]

this is equivalent to (continues from):

G_up, G_down = slater.gradient_matrix(n_vectors)
tr_grad_u = np.einsum('ij,jia->ia', inv_A_up, G_up).reshape(neu * 3)
tr_grad_d = np.einsum('ij,jib->ib', inv_A_down, G_down).reshape(ned * 3)
grad = np.concatenate((tr_grad_u, tr_grad_d))

H_up, H_down = slater.hessian_matrix(n_vectors)
tr_hess_u = np.einsum('ij,jiab->iab', inv_A_up, H_up)
tr_hess_d = np.einsum('ij,jiab->iab', inv_A_down, H_down)
mul_grad_u = np.einsum('ik,kja->ija', inv_A_up, G_up)
mul_grad_d = np.einsum('ik,kjb->ijb', inv_A_down, G_down)
hess_u = np.einsum('ij,iab->iajb', np.eye(neu), tr_hess_u)
hess_d = np.einsum('ij,iab->iajb', np.eye(ned), tr_hess_d)
hess_u -= np.einsum('ijb,jia->iajb', mul_grad_u, mul_grad_u)
hess_d -= np.einsum('ijb,jia->iajb', mul_grad_d, mul_grad_d)
hess = np.zeros((ne * 3, ne * 3))
hess[:neu * 3, :neu * 3] = hess_u.reshape(neu * 3, neu * 3)
hess[neu * 3:, neu * 3:] = hess_d.reshape(ned * 3, ned * 3)
hess += np.outer(grad, grad)

T_up, T_down = slater.tressian_matrix(n_vectors)
tr_tress_u = np.einsum('ij,jiabc->iabc', inv_A_up, T_up)
tr_tress_d = np.einsum('ij,jiabc->iabc', inv_A_down, T_down)
mul_hess_u = np.einsum('ik,kjab->iajb', inv_A_up, H_up)
mul_hess_d = np.einsum('ik,kjab->iajb', inv_A_down, H_down)
tress_u = np.einsum('ij,jk,iabc->iajbkc', np.eye(neu), np.eye(neu), tr_tress_u)
tress_d = np.einsum('ij,jk,iabc->iajbkc', np.eye(ned), np.eye(ned), tr_tress_d)
tress_u -= np.einsum('ij,kajb,jkc->iajbkc', np.eye(neu), mul_hess_u, mul_grad_u)
tress_u -= np.einsum('ki,jaic,ijb->iajbkc', np.eye(neu), mul_hess_u, mul_grad_u)
tress_u -= np.einsum('jk,ibkc,kia->iajbkc', np.eye(neu), mul_hess_u, mul_grad_u)
tress_d -= np.einsum('ij,kajb,jkc->iajbkc', np.eye(ned), mul_hess_d, mul_grad_d)
tress_d -= np.einsum('ki,jaic,ijb->iajbkc', np.eye(ned), mul_hess_d, mul_grad_d)
tress_d -= np.einsum('jk,ibkc,kia->iajbkc', np.eye(ned), mul_hess_d, mul_grad_d)
tress_u += np.einsum('jia,kjb,ikc->iajbkc', mul_grad_u, mul_grad_u, mul_grad_u)
tress_d += np.einsum('jia,kjb,ikc->iajbkc', mul_grad_d, mul_grad_d, mul_grad_d)
tress_u += np.einsum('kia,ijb,jkc->iajbkc', mul_grad_u, mul_grad_u, mul_grad_u)
tress_d += np.einsum('kia,ijb,jkc->iajbkc', mul_grad_d, mul_grad_d, mul_grad_d)
tress = np.zeros((ne * 3, ne * 3, ne * 3))
tress[:neu * 3, :neu * 3, :neu * 3] = tress_u.reshape(neu * 3, neu * 3, neu * 3)
tress[neu * 3:, neu * 3:, neu * 3:] = tress_d.reshape(ned * 3, ned * 3, ned * 3)
tress += (
    np.einsum('i,jk->ijk', grad, hess) +
    np.einsum('k,ij->ijk', grad, hess) +
    np.einsum('j,ki->ijk', grad, hess) -
    2 * np.einsum('i,j,k->ijk', grad, grad, grad)
)

Implementation

The tressian tensor \(T[a, b, c]\) has shape \((3N_e, 3N_e, 3N_e)\) where \(N_e\) is the number of electrons of a given spin. A naïve implementation iterates over all triples \((e_1, e_2, e_3)\) with conditional checks inside the loop body:

for e1 in range(neu):
    for e2 in range(neu):
        for e3 in range(neu):
            res = 0
            if e1 == e2 == e3:
                res += tr_tress[e1, ...]
            if e1 == e2:
                res -= matrix_hess[e3, e1, ...] * matrix_grad[e1, e3, ...]
            # ... etc
            tress[e1*3+r1, e2*3+r2, e3*3+r3] += c * res

This costs \(O(N_e^6 \cdot 27)\) iterations and the runtime branches prevent LLVM from auto-vectorising the inner loops.

The key structural observation is that the determinant-specific contributions to \(T[a, b, c]\) are sparse in the electron indices: each term is non-zero only when at least two of \(e_1, e_2, e_3\) coincide. Specifically:

• \(\mathrm{tr}(A^{-1} \nabla^3 A)\) — non-zero only when \(e_1 = e_2 = e_3\)

• \(\mathrm{tr}(A^{-1} \nabla^2_{e_i e_j} A \cdot A^{-1} \nabla_{e_k} A)\) — non-zero only when \(e_i = e_j\) (one constrained pair, \(e_k\) free)

• \((A^{-1}G)_{e_3 e_2}(A^{-1}G)_{e_1 e_3}(A^{-1}G)_{e_2 e_1}\) — non-zero for all \((e_1, e_2, e_3)\)

This motivates decomposing the computation into five branch-free loop nests:

Loop nest	Complexity	Replaces
\(e_1 = e_2 = e_3 = e\) (tr_tress diagonal)	\(O(N_e \cdot 27)\)	if e1 == e2 == e3
\(e_1 = e_2\), free \(e_3\) (hess×grad, pair 12)	\(O(N_e^2 \cdot 27)\)	if e1 == e2
\(e_1 = e_3\), free \(e_2\) (hess×grad, pair 13)	\(O(N_e^2 \cdot 27)\)	if e1 == e3
\(e_2 = e_3\), free \(e_1\) (hess×grad, pair 23)	\(O(N_e^2 \cdot 27)\)	if e2 == e3
all \((e_1, e_2, e_3)\) free (triple product)	\(O(N_e^3 \cdot 27)\)	always executed

Because no loop nest carries runtime conditionals, LLVM auto-vectorises all five loops via SIMD instructions. The outer-product term \(\nabla\phi \otimes H + H \otimes \nabla\phi + \nabla\phi \otimes \nabla\phi \otimes \nabla\phi\) involves all \((N_e \cdot 3)^3\) elements and is unchanged.

# tr_tress: e1 == e2 == e3
for e in range(neu):
    for r1, r2, r3 in product(range(3), repeat=3):
        tress[e*3+r1, e*3+r2, e*3+r3] += c * tr_tress_u[e, r1, r2, r3]

# hess × grad: e1 == e2, free e3
for e12 in range(neu):
    for e3 in range(neu):
        for r1, r2, r3 in product(range(3), repeat=3):
            tress[e12*3+r1, e12*3+r2, e3*3+r3] -= (
                c * matrix_hess_u[e3, e12, r1, r2] * matrix_grad_u[e12, e3, r3]
            )

# ... similarly for (e1==e3, free e2) and (e2==e3, free e1) ...

# triple product: all (e1, e2, e3) free
for e1 in range(neu):
    for e2 in range(neu):
        for e3 in range(neu):
            for r1, r2, r3 in product(range(3), repeat=3):
                tress[e1*3+r1, e2*3+r2, e3*3+r3] += c * (
                    matrix_grad_u[e3, e2, r2]
                    * matrix_grad_u[e1, e3, r3]
                    * matrix_grad_u[e2, e1, r1]
                    + matrix_grad_u[e2, e3, r3]
                    * matrix_grad_u[e3, e1, r1]
                    * matrix_grad_u[e1, e2, r2]
                )

On the He atom (\(N_e = 1\) per spin, ne3 = 6) the refactored implementation is 5.5× faster than the original (33 µs vs 184 µs per call).

The actual wall-clock speedup depends on system size. The total cost of casino.slater.Slater.tressian() is:

\[t_\text{total} = t_\text{matrix calls} + t_\text{index loops}\]

where \(t_\text{matrix calls}\) is the combined time of value_matrix + gradient_matrix + hessian_matrix + tressian_matrix (all scale with system size and basis set), and \(t_\text{index loops}\) is the time of the electron-index loop nests described above.

For small systems (He, \(N_e = 1\)) the index loops dominate and the 5.5× speedup is realised. For larger systems such as Ar (\(N_e = 9\) per spin) the matrix calls dominate — in particular tressian_matrix grows with the number of orbitals and primitives — and the loop optimisation contributes only ~20% of the total runtime. The next bottleneck to address for large-\(N_e\) systems is therefore tressian_matrix itself.

tressian_dot

The tressian enters the local-energy parameter derivatives only through its contraction with a symmetric matrix \(B = b_g\,b_g^\top\) over the last two axes:

\[(T_B)_a = \sum_{b, c} T_{abc}\, B_{bc}\]

casino.slater.Slater.tressian_dot() returns this vector of shape (ne3,) without ever forming the (ne3, ne3, ne3) tensor, so the working set stays \(O(\text{ne3}^2)\) and cache-resident.

Contracting the dense outer part with \(B\) collapses to a closed form (using \(B = B^\top\)):

\[\sum_{bc} \big(g_c\,\mathrm{PH}_{ab} + g_b\,\mathrm{PH}_{ac} + g_a\,\mathrm{PH}_{bc}\big) B_{bc} = 2\,\big(\mathrm{PH}\,(B\,g)\big)_a + g_a\,\langle \mathrm{PH}, B\rangle\]

with \(\langle \mathrm{PH}, B\rangle = \sum_{bc}\mathrm{PH}_{bc} B_{bc}\). This costs \(O(\text{ne3}^2)\).

The sparse determinant part is contracted inside the same per-spin six-fold loop, accumulating directly into the scalar \((T_B)_a\) instead of storing the tensor. The flop count stays \(O(\text{ne3}^3)\) but the memory footprint drops to \(O(\text{ne3}^2)\).

For certain electron coordinates:

bb = b_g @ b_g.T
t_bb, hess, grad = slater.tressian_dot(n_vectors, bb)

is equivalent to:

tress, hess, grad = slater.tressian(n_vectors)
t_bb = np.tensordot(tress, bb, axes=([1, 2], [0, 1]))

Because the \(\text{ne3}^3\) tensor is never materialised, the speedup over casino.slater.Slater.tressian() grows with system size as the tensor outgrows the cache (wall-clock seconds for an equal number of calls):

System	ne3	tressian	tressian_dot	speedup
Ne	30	0.8	0.6	1.3×
Ar	54	3.0	1.6	1.9×
O3	72	9.9	4.2	2.4×
Kr	108	28.4	7.8	3.6×