# Matching a Slater-type wavefunction with a minimal (STO-nG) Gaussian basis set

I have on hand a high-precision wavefunction expressed in Slater orbitals. I need to express it as accurately as possible in a minimal Gaussian basis set.

For background: I am currently working with STO-6G, but it is not good enough to describe the cusp condition with the accuracy I need. How can I get to a basis like STO-7G, or even beyond to something like STO-10G, preferably in a form that is compatible with QChem?

## TL;DR

The procedure to represent a Slater-type orbital (STO) as a linear combination of Gaussian-type orbitals (GTO) is outlined in W. J. Hehre, R. F. Stewart and J. A. Pople, J. Chem. Phys. 1969, 51, 2657 and is a quite straight forward least square fitting.

## Derivation

I'll be quoting/ rephrasing quite liberally and take the 1s orbital as an example. The STO is defined as \begin{align} \phi_\mathrm{1s}(\zeta_1, \mathbf{r}) &= (\zeta_1^3/\pi)^{1/2}\exp(-\zeta_1 \mathbf{r}). \tag1\label{STO-1s} \end{align}

The general form of a Gaussian-type 1s orbital is defined as \begin{align} g_\mathrm{1s}(\alpha, \mathbf{r}) &= (2\alpha/\pi)^{3/4}\exp(-\alpha \mathbf{r}^2). \tag2\label{GTO-1s-general} \end{align}

The STO is expressed as a sum of $K$ GTOs; the coefficients and exponents are obtained for a STO with $\zeta=1.0$ and the scaled uniformly. Hence we can write \begin{align} \phi_\mathrm{1s}'(\zeta, \mathbf{r}) &= \zeta^{3/2}\phi_\mathrm{1s}'(1.0,\zeta\mathbf{r}), \tag3\label{GTO-1s-zeta1} \end{align} where \begin{align} \phi_\mathrm{1s}'(1.0, \mathbf{r}) &= \sum_k^K d_{\mathrm{1s},k}\,g_\mathrm{1s}(\alpha_k,\mathbf{r}). \tag4\label{GTO-1s-lc} \end{align}

In order to obtain the coefficients $d$ and exponents $\alpha$ the contracted GTO is fitted to the form of the STO via a least square method. Note that this is a quite unique approach, since most other basis sets are optimised to produce the minimal energy. It all comes down to minimising the integral \begin{align} I &= \int \left( \phi_\mathrm{1s}(1.0,\mathbf{r})-\phi_\mathrm{1s}'(1.0,\mathbf{r}) \right)^2 \mathrm{d}\mathbf{r}. \tag5\label{min-int} \end{align}

For $K=3$, i.e. STO-3G, the following values are obtained: \begin{array}{llll}\hline k & \alpha_{\mathrm{1s},k} & d_{\mathrm{1s},k} & I\\\hline & & & 3.31\times10^{-4}\\ 1 & 1.09818 \times10^{-1} & 4.44635 \times10^{-1} &\\ 2 & 4.05771 \times10^{-1} & 5.35328 \times10^{-1} &\\ 3 & 2.22766 & 1.54329 \times10^{-1} &\\\hline \end{array}

## Implementation

Finally you need to scale the exponents to use the right exponent. From $\eqref{GTO-1s-zeta1}$ we see $$\mathbf{r}\mapsto\zeta\mathbf{r}_{\zeta=1.0}$$ and from $\eqref{GTO-1s-general}$ we know $\mathbf{r}\mapsto\mathbf{r}^2$ and therefore the exponent $$\alpha_\zeta\mapsto\alpha_{\zeta=1.0}\zeta^2.$$

Hydrogen for example has $\zeta=1.24$ (Tab. IX) and therefore the exponents that need to be implemented are: \begin{array}{lll}\hline k & \alpha_{\mathrm{1s},k} & d_{\mathrm{1s},k} \\\hline 1 & 0.16885540 & 0.44463454 \\ 2 & 0.62391373 & 0.53532814 \\ 3 & 3.42525091 & 0.15432897 \\\hline \end{array}

For Q-Chem this is explained in the manual under 7.4.3 and needs to have the form \begin{array}{llll} \text{\$basis}\\ \ce{X} & 0\\ L & K & \text{scale}\\ α_1 & C_1^{L_\mathrm{min}} & C_1^{L_\mathrm{min}+1} &\dots& C_1^{L_\mathrm{max}}&\\ α_2 & C_2^{L_\mathrm{min}} & C_2^{L_\mathrm{min}+1} &\dots& C_2^{L_\mathrm{max}}&\\ \vdots & \vdots & \vdots &\ddots & \vdots \\ α_K & C_K^{L_\mathrm{min}} & C_K^{L_\mathrm{min}+1} &\dots& C_K^{L_\mathrm{max}}&\\ \text{****}\\ \text{\$end}\\ \end{array} where $\ce{X}$ is the atomic symbol, $L$ is the angular momentum, $K$ is the degree of contraction, $\text{scale}$ is the scaling applied to all exponents, $\alpha_i$ are the exponents, and $C_i^L$ are the contraction coefficients.

As a practical example, the STO-3G basis set shall look like this:

$basis H 0 S 3 1.00 3.42525091 0.15432897 0.62391373 0.53532814 0.16885540 0.44463454 ****$end


For the actual fitting procedure, i.e. minimising $\eqref{min-int}$ you should consult a more technical oriented network page. It might be on topic at Mathematics, Computational Science, or even Mathematica depending on the question that troubles you. Please check their guidelines before asking.

## General remarks

The most principle shortcoming of GTOs is that they cannot describe the cusp correctly. It's a trade-off you just make for faster integral evaluation. With STO-6G you are getting the orbital quite right already and there is a reason why we usually use STO-3G and not some higher contraction. It just offers the best cost : accuracy ratio.
This can be quite neatly seen in the graphical representation. The black curve is the original STO, purple is STO-6G, cyan is STO-4G, and red is STO-2G plotted online at the wonderful FooPlot. Hack it here (I hope that works).

Going from STO-6G to STO-7G won't give you much more accuracy, maybe not even going to STO-10G. You probably have to go more into the direction of STO-20/30G. At that point it is probably more economic to switch to a program that can handle STO natively. From Wikipedia I see two, that handle those: ADF, and MolDS. The latter seems young, but is free.

If the FooPlot link does not work, here are the used formulae:

(1/pi)^(1/2)*exp(-abs(x))

0.679*(2*0.152/pi)^(3/4)*exp(-0.152*x^2)+
0.430*(2*0.852/pi)^(3/4)*exp(-0.852*x^2)

0.292*(2*0.088/pi)^(3/4)*exp(-0.088*x^2)+
0.533*(2*0.265/pi)^(3/4)*exp(-0.265*x^2)+
0.260*(2*0.955/pi)^(3/4)*exp(-0.955*x^2)+
0.0568*(2*5.22/pi)^(3/4)*exp(-5.22*x^2)

0.130*(2*0.0651/pi)^(3/4)*exp(-0.0651*x^2)+
0.416*(2*0.158/pi)^(3/4)*exp(-0.158*x^2)+
0.371*(2*0.407/pi)^(3/4)*exp(-0.407*x^2)+
0.169*(2*1.19/pi)^(3/4)*exp(-1.19*x^2)+
0.0494*(2*4.24/pi)^(3/4)*exp(-4.24*x^2)+
0.00916*(2*23.1/pi)^(3/4)*exp(-23.1*x^2)

• Thanks a lot for your contribution! Yes, we can find the STO-6G from textbooks or the paper you mentioned. After searching for a couple of days, I really cannot find the coefficients and exponents for STO-NG with N>6. I will see if I can work it out and share here. – James LT Dec 5 '16 at 1:11