Matching a Slater-type wavefunction with a minimal (STO-nG) Gaussian basis set
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}) &amp;= 
  (\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}) &amp;= 
  (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}) &amp;= 
  \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}) &amp;= 
  \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 &amp;=
  \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 &amp; \alpha_{\mathrm{1s},k} &amp; d_{\mathrm{1s},k}      &amp; I\\\hline
  &amp;                        &amp;                        &amp; 3.31\times10^{-4}\\
1 &amp; 1.09818 \times10^{-1}  &amp; 4.44635 \times10^{-1}  &amp;\\
2 &amp; 4.05771 \times10^{-1}  &amp; 5.35328 \times10^{-1}  &amp;\\
3 &amp; 2.22766                &amp; 1.54329 \times10^{-1}  &amp;\\\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 &amp; \alpha_{\mathrm{1s},k} &amp; d_{\mathrm{1s},k}      \\\hline
1 &amp;       0.16885540       &amp;      0.44463454        \\
2 &amp;       0.62391373       &amp;      0.53532814        \\
3 &amp;       3.42525091       &amp;      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} &amp; 0\\   
L &amp; K &amp; \text{scale}\\
α_1  &amp; C_1^{L_\mathrm{min}} &amp; C_1^{L_\mathrm{min}+1} &amp;\dots&amp; C_1^{L_\mathrm{max}}&amp;\\
α_2  &amp; C_2^{L_\mathrm{min}} &amp; C_2^{L_\mathrm{min}+1} &amp;\dots&amp; C_2^{L_\mathrm{max}}&amp;\\
\vdots    &amp; \vdots &amp; \vdots &amp;\ddots &amp; \vdots \\
α_K  &amp; C_K^{L_\mathrm{min}} &amp; C_K^{L_\mathrm{min}+1} &amp;\dots&amp; C_K^{L_\mathrm{max}}&amp;\\
\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. 

Addendum:

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)