How can I get analytical form of GTO primitives for a atom?
Main differenece betweet STO's and GTO's is a $r$-expotetnt. Normalizing factor of GTO one can get from common procedure:
\begin{equation}
N^2 \int_0^\infty \left(r^{n-1}e^{-\zeta r^2}\right)^2 r^2 dr =1 
\end{equation}
Of course, computation is not so simple for $n &gt; 1$, but we can use online WolframAlpha (clickable calculation) and we get the formulas from the question post.
For a Lithium atom we have contraction (6s,3p) -&gt; [2s,1p], which means
the first 3s of 6s GTO's cotract to one s &quot;1s&quot;-STO
$$\mathrm{STO}(1s) = \sqrt{\frac{\zeta^3}{\pi}} e^{-\zeta_1 r},$$
second 3s of 6s GTO's contract to another s &quot;2s&quot;-STO
$$\mathrm{STO}(2s) = \sqrt{\frac{\zeta^5}{3\pi}} r e^{-\zeta_2 r},$$
3p GTO's contract to one p - &quot;2p&quot;-STO
$$\mathrm{STO}(2p) = \sqrt{\frac{\zeta_2^5}{\pi}} x e^{-\zeta_2 r}$$
(notation in orbitals in quotation marks &quot;1s&quot;, &quot;2s&quot; means the same as for the hydrogen atom). So, as we can see, every atomic STO contract by only 3 GTO's (which justifies their name).

For a Litium atom we have for &quot;1s&quot; atomic orbital
\begin{equation}
\sqrt{\frac{\zeta_1^3}{\pi}} e^{-\zeta_1 r} \approx 
C_1 \left(\frac{2\alpha_1}{\pi}\right)^{3/4}e^{-\alpha_1 r^2} + 
C_2 \left(\frac{2\alpha_2}{\pi}\right)^{3/4}e^{-\alpha_2 r^2} +
C_3 \left(\frac{2\alpha_3}{\pi}\right)^{3/4}e^{-\alpha_3 r^2}
\end{equation}
where
\begin{align}
\alpha_1 = 0.1611957475E+02; C_1 = 0.1543289673E+00;\\
\alpha_2 = 0.2936200663E+01; C_2 = 0.5353281423E+00;\\
\alpha_3 = 0.7946504870E+00; C_3 = 0.4446345422E+00;
\end{align}

For &quot;2s&quot; atomic orbital
\begin{equation}
\sqrt{\frac{\zeta_2^5}{3\pi}} r e^{-\zeta_2 r} \approx 
C_4 \left(\frac{2\alpha_1}{\pi}\right)^{3/4}e^{-\alpha_4 r^2} + 
C_5 \left(\frac{2\alpha_2}{\pi}\right)^{3/4}e^{-\alpha_5 r^2} +
C_6 \left(\frac{2\alpha_3}{\pi}\right)^{3/4}e^{-\alpha_6 r^2}
\end{equation}
where
\begin{align}
\alpha_4 = 0.6362897469E+00; C_4 = -0.9996722919E-01;\\ 
\alpha_5 = 0.1478600533E+00; C_5 = 0.3995128261E+00; \\ 
\alpha_6 = 0.4808867840E-01; C_6 = 0.7001154689E+00;  
\end{align}

For &quot;2p_x&quot; atomic orbital
\begin{equation}
\sqrt{\frac{\zeta_2^5}{\pi}} x e^{-\zeta_2 r} \approx 
D_1 \left(\frac{128\alpha_4^5}{\pi^3}\right)^{1/4} xe^{-\alpha_4 r^2} + 
D_2 \left(\frac{128\alpha_5^5}{\pi^3}\right)^{1/4} xe^{-\alpha_5 r^2} +
D_3 \left(\frac{128\alpha_6^5}{\pi^3}\right)^{1/4} xe^{-\alpha_6 r^2}
\end{equation}
where
\begin{align}
\alpha_4 = 0.6362897469E+00; D_1 = 0.1559162750E+00;;\\ 
\alpha_5 = 0.1478600533E+00; D_2 = 0.6076837186E+00; \\ 
\alpha_6 = 0.4808867840E-01; D_3 = 0 0.3919573931E+00;  
\end{align}

(For $p_y$ and $p_z$ --- do the same)
Also, for STO's best fit $\zeta_1 = 2.69$ and $\zeta_2 = 0.81$ are consistent with those given in the article Hehre, W. J., Ditchfield, R., Stewart, R. F., &amp; Pople, J. A. (1970). Self‐Consistent Molecular Orbital Methods. IV. Use of Gaussian Expansions of Slater‐Type Orbitals. Extension to Second‐Row Molecules. The Journal of Chemical Physics, 52(5), 2769–2773. doi:10.1063/1.1673374