Data stored in TREXIO

We consider here basis functions centered on nuclei. Hence, it is possibile to define dummy atoms to place basis functions in arbitrary positions.

The atomic basis set is defined as a list of shells. Each shell \(s\) is centered on a center \(A\), possesses a given angular momentum \(l\) and a radial function \(R_s\). The radial function is a linear combination of \(N_{\text{prim}}\) primitive functions that can be of type Slater (\(p=1\)) or Gaussian (\(p=2\)), parameterized by exponents \(\gamma_{ks}\) and coefficients \(a_{ks}\): \[ R_s(\mathbf{r}) = \mathcal{N}_s \vert\mathbf{r}-\mathbf{R}_A\vert^{n_s} \sum_{k=1}^{N_{\text{prim}}} a_{ks}\, f_{ks}(\gamma_{ks},p)\, \exp \left( - \gamma_{ks} \vert \mathbf{r}-\mathbf{R}_A \vert ^p \right). \]

In the case of Gaussian functions, \(n_s\) is always zero.

Different codes normalize functions at different levels. Computing normalization factors requires the ability to compute overlap integrals, so the normalization factors should be written in the file to ensure that the file is self-contained and does not need the client program to have the ability to compute such integrals.

Some codes assume that the contraction coefficients are for a linear combination of normalized primitives. This implies that a normalization constant for the primitive \(ks\) needs to be computed and stored. If this normalization factor is not required, \(f_{ks}=1\).

Some codes assume that the basis functions are normalized. This implies the computation of an extra normalization factor, \(\mathcal{N}_s\). If the the basis function is not considered normalized, \(\mathcal{N}_s=1\).

All the basis set parameters are stored in one-dimensional arrays.

Trexio supports numerical atom centered orbitals. The implementation is based on the approach of FHI-aims [Blum, V. et al; Ab initio molecular simulations with numeric atom-centered orbitals; Computer Physics Communications 2009]. These orbitals are defined by the atom they are centered on, their angular momentum and a radial function \(R_s\), which is of the form \[ R_s(\mathbf{r}) = \mathcal{N}_s \frac{u_i(\mathbf{r})}{r^{-n_s}}. \] Where \(u_i(\mathbf{r})\) is numerically tabulated on a dense logarithmic grid. It is constructed to vanish for any \(\mathbf{r}\) outside of the grid. The reference points are stored in nao_grid_r and nao_grid_phi. Additionaly, a separate spline for the first and second derivative of \(u(\mathbf{r})\) can be stored in nao_grid_grad and nao_grid_lap. Storing them in this form allows to calculate the actual first and second derivatives easily as follows:

\[ \frac{\partial \phi}{\partial x} = \frac{x}{r^2}\left( u^\prime\left(r\right) - \frac{u\left(r\right)}{r}\right) \] \[ \frac{\partial^2 \phi}{\partial x^2} = \frac{1}{r^3}\left(x^2 u^{\prime\prime}(r) + \left( 3x^2-r^2\right) \left( \frac{u(r)}{r^2} - \frac{u'(r)}{r}\right) \right) \]

The index of the first data point for each shell is stored in nao_grid_start, the number of data points per spline is stored in nao_grid_size for convenience.

What kind of spline is used can be provided in the interpolator_kind field. For example, FHI-aims uses a cubic spline, so the interpolator_kind is \"Polynomial\" and the interp_coeff_cnt is \(4\). In this case, the first interpolation coefficient per data point is the absolute term, the second is for the linear term etc. The interpolation coefficients for the wave function are given in the interpolator_phi array. The interpolator_grad and interpolator_lap arrays provide a spline for the gradient and Laplacian, respectively. The argument passed to the interpolants is on the logarithmic scale of the reference points: If the argument is an integer \(i\), the interpolant will return the value of \(u(\mathbf{r})\) at the $i$th reference point. A radius is converted to this scale by (note the zero-indexing) \[ i_{\log} = \frac{1}{c} \cdot \log \left( \frac{r}{r_0} \right) \] where \[ c = \log\left(\frac{r_1}{r_0}\right) \] For convenience, this conversion and functions to evaluate the splines are provided with trexio. Since these implementations are not adapted to a specific software architecture, a programm using these orbitals should reimplement them with consideration for its specific needs.

A plane wave is defined as

\[ \chi_j(\mathbf{r}) = \exp \left( -i \mathbf{G}_j \cdot \mathbf{r} \right) \]

The basis set is defined as the array of $k$-points in the reciprocal space \(\mathbf{G}_j\), defined in the pbc group. The kinetic energy cutoff e_cut is the only input data relevant to plane waves.

Basis functions can be made oscillating as \[ R_s(\mathbf{r}) = \mathcal{N}_s \vert\mathbf{r}-\mathbf{R}_A\vert^{n_s} \sum_{k=1}^{N_{\text{prim}}} a_{ks}\, f_{ks}(\gamma_{ks},p)\, \exp \left( - \gamma_{ks} \vert \mathbf{r}-\mathbf{R}_A \vert ^p \right)\, \cos \left( \beta_{ks} \vert \mathbf{r}-\mathbf{R}_A \vert ^q \right) \]

Oscillation kind can be:

• Cos1 for \(q=1\)
• Cos2 for \(q=2\)

For example, consider H₂ with the following basis set (in GAMESS format), where both the AOs and primitives are considered normalized:

HYDROGEN
S   5
1         3.387000E+01           6.068000E-03
2         5.095000E+00           4.530800E-02
3         1.159000E+00           2.028220E-01
4         3.258000E-01           5.039030E-01
5         1.027000E-01           3.834210E-01
S   1
1         3.258000E-01           1.000000E+00
S   1
1         1.027000E-01           1.000000E+00
P   1
1         1.407000E+00           1.000000E+00
P   1
1         3.880000E-01           1.000000E+00
D   1
1         1.057000E+00           1.000000E+00

In TREXIO representaion we have:

type  = "Gaussian"
prim_num   = 20
shell_num   = 12

# 6 shells per H atom
nucleus_index =
[ 0, 0, 0, 0, 0, 0,
  1, 1, 1, 1, 1, 1 ]

# 3 shells in S (l=0), 2 in P (l=1), 1 in D (l=2)
shell_ang_mom =
[ 0, 0, 0, 1, 1, 2,
  0, 0, 0, 1, 1, 2 ]

# no need to renormalize shells
shell_factor =
[ 1., 1., 1., 1., 1., 1.,
  1., 1., 1., 1., 1., 1. ]

# 5 primitives for the first S shell and then 1 primitive per remaining shells in each H atom
shell_index =
[ 0, 0, 0, 0, 0, 1, 2, 3, 4, 5,
  6, 6, 6, 6, 6, 7, 8, 9, 10, 11 ]

# parameters of the primitives (10 per H atom)
exponent =
[ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407, 0.388, 1.057,
  33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407, 0.388, 1.057 ]

coefficient =
[ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0,
  0.006068, 0.045308, 0.202822, 0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0 ]

prim_factor =
[ 1.0006253235944540e+01, 2.4169531573445120e+00, 7.9610924849766440e-01
  3.0734305383061117e-01, 1.2929684417481876e-01, 3.0734305383061117e-01,
  1.2929684417481876e-01, 2.1842769845268308e+00, 4.3649547399719840e-01,
  1.8135965626177861e+00, 1.0006253235944540e+01, 2.4169531573445120e+00,
  7.9610924849766440e-01, 3.0734305383061117e-01, 1.2929684417481876e-01,
  3.0734305383061117e-01, 1.2929684417481876e-01, 2.1842769845268308e+00,
  4.3649547399719840e-01, 1.8135965626177861e+00 ]

For example, consider H₂ molecule with the following effective core potential (in GAMESS input format for the H atom):

H-ccECP GEN 0 1
3
1.00000000000000    1 21.24359508259891
21.24359508259891   3 21.24359508259891
-10.85192405303825  2 21.77696655044365
1
0.00000000000000    2 1.000000000000000

In TREXIO representation this would be:

num = 8

# lmax+1 per atom
max_ang_mom_plus_1 = [ 1, 1 ]

# number of core electrons to remove per atom
zcore = [ 0, 0 ]

# first 4 ECP elements correspond to the first H atom ; the remaining 4 elements are for the second H atom
nucleus_index = [
  0, 0, 0, 0,
  1, 1, 1, 1
  ]

# 3 first ECP elements correspond to potential of the P orbital (l=1), then 1 element for the S orbital (l=0) ; similar for the second H atom
ang_mom = [
  1, 1, 1, 0,
  1, 1, 1, 0
  ]

# ECP quantities that can be attributed to atoms and/or angular momenta based on the aforementioned ecp_nucleus and ecp_ang_mom arrays
coefficient = [
  1.00000000000000, 21.24359508259891, -10.85192405303825, 0.00000000000000,
  1.00000000000000, 21.24359508259891, -10.85192405303825, 0.00000000000000
  ]

exponent = [
  21.24359508259891, 21.24359508259891, 21.77696655044365, 1.000000000000000,
  21.24359508259891, 21.24359508259891, 21.77696655044365, 1.000000000000000
  ]

power = [
  -1, 1, 0, 0,
  -1, 1, 0, 0
  ]

• \[ \hat{V}_{\text{ne}} = \sum_{A=1}^{N_\text{nucl}} \sum_{i=1}^{N_\text{elec}} \frac{-Z_A }{\vert \mathbf{R}_A - \mathbf{r}_i \vert} \] : electron-nucleus attractive potential,
• \[ \hat{T}_{\text{e}} = \sum_{i=1}^{N_\text{elec}} -\frac{1}{2}\hat{\Delta}_i \] : electronic kinetic energy
• \(\hat{h} = \hat{T}_\text{e} + \hat{V}_\text{ne} + \hat{V}_\text{ECP}\) : core electronic Hamiltonian

The one-electron integrals for a one-electron operator \(\hat{O}\) are \[ \langle p \vert \hat{O} \vert q \rangle \], returned as a matrix over atomic orbitals.

The two-electron integrals for a two-electron operator \(\hat{O}\) are \[ \langle p q \vert \hat{O} \vert r s \rangle \] in physicists notation, where \(p,q,r,s\) are indices over atomic orbitals.

• \[ \hat{W}_{\text{ee}} = \sum_{i=2}^{N_\text{elec}} \sum_{j=1}^{i-1} \frac{1}{\vert \mathbf{r}_i - \mathbf{r}_j \vert} \] : electron-electron repulsive potential operator.
• \[ \hat{W}^{lr}_{\text{ee}} = \sum_{i=2}^{N_\text{elec}} \sum_{j=1}^{i-1} \frac{\text{erf}(\mu\, \vert \mathbf{r}_i - \mathbf{r}_j \vert)}{\vert \mathbf{r}_i - \mathbf{r}_j \vert} \] : electron-electron long range potential

The Cholesky decomposition of the integrals can also be stored:

\[ \langle ij | kl \rangle = \sum_{\alpha} G_{ik\alpha} G_{jl\alpha} \]

The operators as the same as those defined in the AO one-electron integrals section. Here, the integrals are given in the basis of molecular orbitals.

The operators are the same as those defined in the AO two-electron integrals section. Here, the integrals are given in the basis of molecular orbitals.

The first form of Jastrow factor is the one used in the CHAMP program:

\[ J(\mathbf{r},\mathbf{R}) = J_{\text{eN}}(\mathbf{r},\mathbf{R}) + J_{\text{ee}}(\mathbf{r}) + J_{\text{eeN}}(\mathbf{r},\mathbf{R}) \]

\(J_{\text{eN}}\) contains electron-nucleus terms:

\[ J_{\text{eN}}(\mathbf{r},\mathbf{R}) = \sum_{i=1}^{N_\text{elec}} \sum_{\alpha=1}^{N_\text{nucl}}\left[ \frac{a_{1,\alpha}\, f_\alpha(R_{i\alpha})}{1+a_{2,\alpha}\, f_\alpha(R_{i\alpha})} + \sum_{p=2}^{N_\text{ord}^a} a_{p+1,\alpha}\, [f_\alpha(R_{i\alpha})]^p - J_{\text{eN}}^\infty \right] \]

\(J_{\text{ee}}\) contains electron-electron terms:

\[ J_{\text{ee}}(\mathbf{r}) = \sum_{i=1}^{N_\text{elec}} \sum_{j=1}^{i-1} \left[ \frac{\frac{1}{2}\big(1 + \delta^{\uparrow\downarrow}_{ij}\big)\,b_1\, f_{\text{ee}}(r_{ij})}{1+b_2\, f_{\text{ee}}(r_{ij})} + \sum_{p=2}^{N_\text{ord}^b} b_{p+1}\, [f_{\text{ee}}(r_{ij})]^p - J_{\text{ee},ij}^\infty \right] \]

\(\delta^{\uparrow\downarrow}_{ij}\) is zero when the electrons \(i\) and \(j\) have the same spin, and one otherwise.

\(J_{\text{eeN}}\) contains electron-electron-Nucleus terms:

\[ J_{\text{eeN}}(\mathbf{r},\mathbf{R}) = \sum_{\alpha=1}^{N_{\text{nucl}}} \sum_{i=1}^{N_{\text{elec}}} \sum_{j=1}^{i-1} \sum_{p=2}^{N_{\text{ord}}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}} c_{lkp\alpha} \left[ g_{\text{ee}}({r}_{ij}) \right]^k \nonumber \\ \left[ \left[ g_\alpha({R}_{i\alpha}) \right]^l + \left[ g_\alpha({R}_{j\alpha}) \right]^l \right] \left[ g_\alpha({R}_{i\,\alpha}) \, g_\alpha({R}_{j\alpha}) \right]^{(p-k-l)/2} \] \(c_{lkp\alpha}\) are non-zero only when \(p-k-l\) is even.

The terms \(J_{\text{ee},ij}^\infty\) and \(J_{\text{eN}}^\infty\) are shifts to ensure that \(J_{\text{eN}}\) and \(J_{\text{ee}}\) have an asymptotic value of zero:

\[ J_{\text{eN}}^{\infty} = \frac{a_{1,\alpha}\, \kappa_\alpha^{-1}}{1+a_{2,\alpha}\, \kappa_\alpha^{-1}} + \sum_{p=2}^{N_\text{ord}^a} a_{p+1,\alpha}\, \kappa_\alpha^{-p} \] \[ J_{\text{ee},ij}^{\infty} = \frac{\frac{1}{2}\big(1 + \delta^{\uparrow\downarrow}_{ij}\big)\,b_1\, \kappa_{\text{ee}}^{-1}}{1+b_2\, \kappa_{\text{ee}}^{-1}} + \sum_{p=2}^{N_\text{ord}^b} b_{p+1}\, \kappa_{\text{ee}}^{-p} \]

\(f\) and \(g\) are scaling function defined as

\[ f_\alpha(r) = \frac{1-e^{-\kappa_\alpha\, r}}{\kappa_\alpha} \text{ and } g_\alpha(r) = e^{-\kappa_\alpha\, r}, \]

Mu-Jastrow is based on a one-parameter correlation factor that has been introduced in the context of transcorrelated methods. This correlation factor imposes the electron-electron cusp, and it is built such that the leading order in \(1/r_{12}\) of the effective two-electron potential reproduces the long-range interaction of the range-separated density functional theory. Its analytical expression reads

\[ J(\mathbf{r}, \mathbf{R}) = J_{\text{eeN}}(\mathbf{r}, \mathbf{R}) + J_{\text{eN}}(\mathbf{r}, \mathbf{R}) \].

The electron-electron cusp is incorporated in the three-body term

\[ J_\text{eeN} (\mathbf{r}, \mathbf{R}) = \sum_{i=1}^{N_\text{elec}} \sum_{j=1}^{i-1} \, u\left(\mu, r_{ij}\right) \, \Pi_{\alpha=1}^{N_{\text{nucl}}} \, E_\alpha({R}_{i\alpha}) \, E_\alpha({R}_{j\alpha}), \]

where ww\(u\) is an electron-electron function

\[ u\left(\mu, r\right) = \frac{r}{2} \, \left[ 1 - \text{erf}(\mu\, r) \right] - \frac{1}{2 \, \mu \, \sqrt{\pi}} \exp \left[ -(\mu \, r)^2 \right]. \]

This electron-electron term is tuned by the parameter \(\mu\) which controls the depth and the range of the Coulomb hole between electrons.

An envelope function has been introduced to cancel out the Jastrow effects between two-electrons when at least one is close to a nucleus (to perform a frozen-core calculation). The envelope function is given by

\[ E_\alpha(R) = 1 - \exp\left( - \gamma_{\alpha} \, R^2 \right). \]

In particular, if the parameters \(\gamma_\alpha\) tend to zero, the Mu-Jastrow factor becomes a two-body Jastrow factor:

\[ J_{\text{ee}}(\mathbf{r}) = \sum_{i=1}^{N_\text{elec}} \sum_{j=1}^{i-1} \, u\left(\mu, r_{ij}\right) \]

and for large \(\gamma_\alpha\) it becomes zero.

To increase the flexibility of the Jastrow and improve the electron density the following electron-nucleus term is added

\[ J_{\text{eN}}(\mathbf{r},\mathbf{R}) = \sum_{i=1}^{N_\text{elec}} \sum_{\alpha=1}^{N_\text{nucl}} \, \left[ \exp\left( a_{\alpha} R_{i \alpha}^2 \right) - 1\right]. \]

The parameter \(\mu\) is stored in the ee array, the parameters \(\gamma_\alpha\) are stored in the een array, and the parameters \(a_\alpha\) are stored in the en array.