% LaTeX
\documentclass[12pt]{iopart}
\begin{document}
\jl{2}
\title[Kinetic ballance]
{Kinetically-balanced Dirac equation:\\
Properties and applications}

\author{Jacek Karwowski\dag\footnote[3]{To whom 
correspondence should be addressed}, 
Monika Stanke\dag and Hiroshi Tatewaki\ddag}
 
\address{\dag\ Instytut Fizyki, Uniwersytet Miko\l{}aja Kopernika, 
Grudzi\c{a}dzka 5, 87-100 Toru\'n, Poland}

\address{\ddag\ Computation Center and Institute of Natural Sciences,
Nagoya City University, Nagoya, Aichi 467-8501, Japan} 

\begin{abstract}
Properties of the Dirac equation defined in a space restricted by imposing
the non-relativistic kinetic balance condition on the small component of the
Dirac spinor (referred to as the kinetically-balanced Dirac equation) are
analysed. In particular links between the kinetically balanced Dirac
equation and the corresponding Schr\"odinger equation are discussed. It is
proposed that the kinetically-balanced Dirac equation may be used as a
simple tool to generate basis sets for relativistic calculations from their
non-relativistic counterparts.
\end{abstract}

\pacs{}
\submitted

\section{Introduction}

With the formulation by Hall \cite{hall} and Roothaan \cite{roothaan} of the
method of solving Hartree-Fock equations based on basis-set expansions of
orbitals, emerged the problem of an optimum choice of the basis set for the
expansion. Over the years two main approaches to the basis set problem
developed. One, usually referred to as the {\it algebraic} approach is based
on a systematic generation of basis sets imposing some simple conditions
upon the basis functions (as for example selection of the exponents so that
they form a geometric progression). A very good review of this approach has
recently been given by Wilson \cite{wilson}. Another, {\it numerical},
approach is based on a generation of the basis sets from a condition
(usually derived from the variational principle) for some model systems
(usually for atoms). These basis sets are then tabulated for future
applications. There are several reviews of this approach \cite{r1, r2} and
extensive tabulations \cite{t1, t2, t3}.

Probably the most general an efficent implementation of the numerical
approach in non-relativistic calculations is based on the idea of a
simultaneous optimization of both linear and nonlinear parameters in the
basis-set-expanded orbitals \cite{tahu}. In a series of recent papers
\cite{tat1b, tat2b, tat3b} sets of the energetically optimal contracted
Gaussian orbitals have been derived by a direct minimization of the
Hartree-Fock-type Rayleigh quotient in the space of all parameters (the
linear coefficients and the exponents) defining the orbital. The
calculations performed for the ground states of all atoms from He through 
Kr resulted in very high quality contracted Gaussian basis sets.

The problem of selection of an optimum basis set for relativistic
calculations is apparently similar to that in the non-relativistic case,
however it is much more subtle. An incorrect choosing basis sets in early
attempts to solving Dirac-Fock equations \cite{malli} resulted in the so
called {\it variational collapse}. The reasons for this behaviour of the
algebraic Dirac-Fock equations have been explained in works of Grant {\it et
al.} \cite{grant1, grant2} and of Kutzelnigg \cite{kutz1}. The explanation
carried also a prescription on what to do in order to prevent this
udesirable phenomenon: the basis set has to fulfill the correct boundary
conditions. Though the conditions which are sufficient to prevent the
variational collapse in relativistic calculations are easy to formulate,
their direct implementation is not always possible. In particular, an
incorrect asymptotic behaviour at both nucleus and infinity is inherent for
the Gaussian functions, most commonly used in molecular variational
calculations. Fortunately there are several theorems which allow us for a
safe use of the variational principle also in this case. As proven by Grant
\cite{grant2}, if the Gaussian basis is kinetically balanced and if the
linear parameters in the expansions of the Dirac-Fock orbitals are the only
ones which are varied, then the variational collapse would not appear unless
the maximum exponent is larger than a certain (rather large) threshold
value. This prescription leads to a straightforward generalization of the
non-relativistic parameter-optimization-based approach: take a basis set in
which the small component is related to the large one by the kinetic balance
condition and minmize the Rayleigh quotient corresponding to a properly
selected model Dirac-Fock problem with respect to both linear and non-linear
basis set parameters. Unfortunately, the results of this procedure obtained
with several-term expansions suffered from the variational collapse already
for the medium $Z$ atoms \cite{ht}. An analysis of this phenomenon and
formulation of an alternative approach is the aim of this paper.

Atomic units are used in this paper. The velocity of light is taken to be
$c=1/\alpha=137.0359895$.

\section{Kinetically balanced Dirac equation}

Let us consider a single Dirac electron in a field of a Coulomb-type radial
potential U(r). For the time being $U(r)$ may represent any kind of a local
atomic potential as e.g. a Dirac-Slater one, any kind of a model potential
or simply a hydrogenic Coulomb potential. We assume that the Dirac spinor 
is selected as 
\begin{equation}
\Psi={{\Psi^{\rm L}}\choose{\Psi^{\rm S}}},
\end{equation}
where
\begin{equation}
\Psi^{\rm L}=\frac{1}{r}A^{\rm L}(\Omega)\,\Phi^{\rm L}(r)\;\;\mbox{and}\;\;
\Psi^{\rm S}=\frac{1}{r}A^{\rm S}(\Omega)\,\Phi^{\rm S}(r)
\end{equation}
are, respectively, large and small components of the wavefunction. The
angular and spin parts of the trial functions, $A(\Omega)$, are equal to the
exact ones.
The radial Dirac equation, for the case under consideration, reads
\begin{equation}
\left( 
\begin{array}{cc}
U(r)-E&c\,{\hat{\cal O}}^\dag\\[10pt]
c\,{\hat{\cal O}}&U(r)-2c^2-E 
\end{array}
\right) 
\left( 
\begin{array}{c} \Phi^{\rm L}\\[10pt] \Phi^{\rm S}
\end{array} 
\right) = 0, 
\label{dirac1}
\end{equation}
where
\begin{equation}
{\hat{\cal O}}=\frac{\kappa}{r}+\frac{\rmd}{\rmd r}
\label{calo}
\end{equation}
is the radial part of $(\mbox{\boldmath$\sigma$}\cdot\bi{p})$ operator
($\kappa$ is the Dirac angular momentum quantum number, $\mbox{\boldmath$
\sigma$}$ are  the $2\times 2$ Pauli matrices and and $\bi{p}$ is the
momentum operator). As one can easily see,
\begin{equation}
{\hat{\cal O}}^\dag{\hat{\cal O}}=
-\frac{\rmd^2}{\rmd r^2}+\frac{\kappa(\kappa+1)}{r^2}
\label{calo1}
\end{equation}
and in the case of $s_{1/2}$ states $(\kappa=-1)$ it becomes equal to the
negative of the second derivative with respect to $r$.

According to the kinetic balance condition
\begin{equation}
\Phi^{\rm S}(r)={\hat{\cal O}}\Phi^{\rm L}(r).
\label{kb}
\end{equation}
If the small component is determined by this condition then the the Dirac 
equation (\ref{dirac1}) becomes
\begin{equation}
\left( 
\begin{array}{cc}
U-E&c\,{\hat{\cal O}}^\dag\\[10pt]
c\,{\hat{\cal O}}&U-2c^2-E 
\end{array}
\right) 
\left( 
\begin{array}{c} C^{\rm L}\Phi^{\rm L}\\[10pt] 
C^{\rm S}{\hat{\cal O}}\Phi^{\rm L}
\end{array} 
\right) = 0.
\label{dirackb}
\end{equation}
Equation (\ref{dirackb}) is referred to as the {\it kinetically balanced}
Dirac equation. Obviously, a kinetically balanced wavefunction cannot be an
exact solution of the Dirac equation, except for the free-particle. However
the equation may be used as a generator of kinetically balanced basis sets.

Equation (\ref{dirackb}) may be rewriten as a system of two equations:
\begin{eqnarray}
\left[\frac{1}{2m}{\hat{\cal O}}^\dag{\hat{\cal O}}
+U-E_m\right]\Phi_m&=&0, \label{schr} \\[8pt]
\left[2m-2+\alpha^2U-\alpha^2E_m\right]{\hat{\cal O}}\Phi_m&=&0,
\label{cond1}
\end{eqnarray}
where 
\begin{equation}
m=\alpha\frac{C^{\rm L}}{2C^{\rm S}} 
\label{mass}
\end{equation}
and $\Phi_m$ stands for $\Phi^{\rm L}$
normalized to unity. The first of them is identical to the Schr\"odinger
equation for a particle with an effective mass $m$.  The second equation
cannot be fulfilled if $U$ is a function of $r$. Then, the best we can do is
to convert it to a relation between the expectation values which is
appropriate for determining $m$ in equation (\ref{schr}). By mutiplying
equation (\ref{cond1}) by $({\hat{\cal O}}\Phi_m)^\dag$ and integrating over
$r$ we obtain
\begin{equation}
2T_m(2m-2-\alpha^2E_m)+\alpha^2\,X_m=0,
\label{cond2}
\end{equation}  
where $\Omega_m=\langle\Phi_m|{\hat\Omega}|\Phi_m\rangle$. If we introduce
a new variable $\rho=mr$ then one can easily see that
\begin{equation}
E_m=m\,E_1,\;\;\;T_m=m\,T_1,\;\;\;X_m=m^2\,X_1
\end{equation} 
and 
\begin{equation}
m=\frac{1}{1+\delta}, \;\;\; \mbox{where} \;\;\; 
\delta=\frac{\alpha^2}{4}\left(\frac{X_1}{T_1}-2\,E_1\right).
\label{mhyd}
\end{equation}
If $U=-\frac{Z}{R}$ (the hydrogenic potential), then $\Phi_m(r)$ is the 
hydrogenic radial wavefunction, 
\begin{equation}
E_1=-\frac{Z^2}{2n^2},\;\;\;V_1=2E_1,\;\;\;T_1=-E_1,\;\;\;
X_1=-\frac{Z^4}{n^4}\left(\frac{2n}{|\kappa|}-1\right)
\end{equation}
and
\begin{equation}
\delta=\frac{\alpha^2Z^2}{n^2}
\left(\frac{3}{4}-\frac{n}{|\kappa|}\right).
\end{equation}
Consequently,
\begin{equation}
E=E_m=-\frac{Z^2}{2n^2(1+\delta)}.
\end{equation}
By expanding the last equation into a power series with respect to
$\alpha Z$ we obtain
\begin{equation}
E=-\frac{Z^2}{2n^2}\left[1-\frac{\alpha^2Z^2}{2n^2}\left(
\frac{3}{4}-\frac{n}{|\kappa|}\right)+O(Z^4\alpha^4)\right],
\end{equation}
i.e. the correct Pauli energy.
    
\section{Algebraic representations}

Let $\phi_k(r)$, $k=1,2,\ldots N$ be a primitive basis set of radial
functions which depend upon some nonlinear parameters and let
\begin{equation}
\Phi^{\rm L}(r)=\sum_{k=1}^{N}C_k^{\rm L}\phi_k(r)
\label{phil}
\end{equation}
where $C_k^{\rm L}$ are the linear variational parameters.
We shall consider two different ways of representing 
the small component of the radial part of the wavefunction. In the 
first one, referred to hereafter as the {\it uncontracted} case
\begin{equation}
\Phi_{\rm u}^{\rm S}(r)=\sum_{k=1}^{N}C_k^{\rm S}{\hat{\cal O}}\phi_k(r)
\label{phisu}
\end{equation}
where $C_k^{\rm S}$ are variational parameters independent of 
$C_k^{\rm L}$. In the second representation, referred to as the {\it
contracted} case,
\begin{equation}
\Phi_{\rm c}^{\rm S}(r)={\hat{\cal O}}\Phi_{\rm c}^{\rm L}(r).
\label{phisc}
\end{equation}
In this case the small component is determined solely by the kinetic
balance condition.

In the {\it uncontracted} case the radial Dirac equation (\ref{dirac1})
leads to the following algebraic $2N\times 2N$ eigenvalue problem:
\begin{equation}
\mbox{\boldmath$H$}(E)\mbox{\boldmath$C$}=0,
\label{dirac2}
\end{equation}
where
\begin{equation}
\mbox{\boldmath$H$}(E)=
\left(
\begin{array}{cc}
\mbox{\boldmath$V$}-E\mbox{\boldmath$S$}&
\;\;\;2c\,\mbox{\boldmath$T$}\\[5pt]
2c\,\mbox{\boldmath$T$}&
\;\;\;\mbox{\boldmath$X$}-2(E+2c^2)\mbox{\boldmath$T$}
\end{array}
\right),
\label{matrix}
\end{equation}
$E$ is the eigenvalue, 
\[
\mbox{\boldmath$C$}=\left(
\begin{array}{c}
\mbox{\boldmath$C$}^{\rm L}\\
\mbox{\boldmath$C$}^{\rm S}
\end{array}
\right)
\]
is the corresponding eigenvector and
\begin{eqnarray*} 
\left({\mbox{\boldmath$V$}}\right)_{kk^\prime}&=&\langle\phi_k|
U|\phi_{k^\prime}\rangle,\\
\left({\mbox{\boldmath$T$}}\right)_{kk^\prime}
&=&{\textstyle\frac{1}{2}}\langle\phi_k|
{\hat{\cal O}}^\dag{\hat{\cal O}}|\phi_{k^\prime}\rangle,\\
\left({\mbox{\boldmath$X$}}\right)_{kk^\prime}&=&\langle\phi_k|
{\hat{\cal O}}^\dag\,U\,{\hat{\cal O}}|\phi_{k^\prime}\rangle,\\
\left({\mbox{\boldmath$S$}}\right)_{kk^\prime}
&=&\langle\phi_k|\phi_{k^\prime}\rangle.
\end{eqnarray*}
As it is known \cite{grant1, grant2}, the lowest positive energy solutions
of equation (\ref{dirac2}) converge to the bound-state solutions of the
corresponding Dirac equations if the primitive basis set approaches
completeness. The $2N \times 2N$ matrix $\mbox{\boldmath$H$}(E)$ is the
representative of the Dirac operator of equation (\ref{dirac1})
in the large/small component space spanned by 
$\{\phi_k\}$/$\{{\hat{\cal O}}\phi_k\}$, $k=1,2,\ldots,N$ or, alternatively, 
as a representative of the kinetically balanced Dirac operator of
equation (\ref{dirackb}) in the space spanned by $\{\phi_k\}$. In both
cases the result is the same since the spinor and the orbital degrees
of freedom are in each of these cases combined differently. 

The {\it contracted} case, by definition, cannot lead to an exact solution
of the Dirac equation even in the limit of a complete basis set. Therefore
its analysis did not receive much attention in the literature. However it
may be very useful as a generator of kinetically-balanced contracted basis
sets for relativistic calculations. The algebraic representation of the
Dirac equation in the {\it contracted} case is given by a $2\times 2$ matrix
and reads
\begin{equation}
\left( 
\begin{array}{cc}
V_{\rm c}-ES_{\rm c} & 2cT_{\rm c}\\[10pt]
2cT_{\rm c} & X_{\rm c}-2(E+2c^2)T_{\rm c} 
\end{array}
\right) 
\left( 
\begin{array}{c} C^{\rm L}\\[10pt] 
C^{\rm S}
\end{array} 
\right) = 0,
\label{kbmatrix}
\end{equation}
where
\begin{eqnarray} 
V_{\rm c}&=&\langle\Phi^{\rm L}|U|\Phi^{\rm L}\rangle=
\sum_{k=1}^N\sum_{k^\prime=1}^Nc_k^*c_k^\prime
\langle\phi_k|U|\phi_{k^\prime}\rangle,
\label{defv}\\
T_{\rm c}&=&{\textstyle\frac{1}{2}}
\langle\Phi^{\rm L}|{\hat{\cal O}}^\dag\,{\hat{\cal O}}
|\Phi^{\rm L}\rangle=
{\textstyle\frac{1}{2}}\sum_{k=1}^N\sum_{k^\prime=1}^N
c_k^*c_k^\prime\langle\phi_k|
{\hat{\cal O}}^\dag{\hat{\cal O}}|\phi_{k^\prime}\rangle,
\label{deft}\\
X_{\rm c}&=&\langle\Phi^{\rm L}|{\hat{\cal O}}^\dag\,U\,
{\hat{\cal O}}|\Phi^{\rm L}\rangle=
\sum_{k=1}^N\sum_{k^\prime=1}^Nc_k^*c_k^\prime
\langle\phi_k|{\hat{\cal O}}^\dag\,U\,
{\hat{\cal O}}|\phi_{k^\prime}\rangle,
\label{defx}\\
S_{\rm c}&=&\langle\Phi^{\rm L}|\Phi^{\rm L}\rangle=
\sum_{k=1}^N\sum_{k^\prime=1}^Nc_k^*c_k^\prime
\langle\phi_k|\phi_{k^\prime}\rangle.
\label{defs}
\end{eqnarray}
The linear coefficients $c_k$ and the non-linear parameters in the basis
functions may be determined from a variational condition imposed
upon the pertinent Rayleigh quotient \cite{ht}. 

The question to be answered is the way one should treat variationally 
equation (\ref{kbmatrix}). One might be tempted to treat it as a $2 \times
2$ secular problem and to determine the remaining linear and non-lienar 
parameters from the energy minimum condition or from a mini-max principle
\cite{talman, kutzmm}. However, consulting the ``exact'' kinetically
balanced Dirac equation (ref{dirackb}) - (\ref{cond1}), we can see that
such an approach would imply that we are trying to solve a
self-contradictory (???? I have to find in a dictionary the correct 
English name ????) system of equations. Therefore, with an enlargement
of the basis set one may the results to deteriorate rather than to improve
(an illustration of this behaviour is given in the next section).
Instead we can rearrange equation (\ref{kbmatrix}):
\begin{eqnarray}
{}&\frac{1}{m}T_{\rm c}+V_{\rm c}-ES_{\rm c}=0\\
\label{alkb1}
{}&m+\frac{\alpha^2}{4}\left(\frac{X_{\rm c}}{T_{\rm c}}-2E\right)-1=0,
\label{alkb2}
\end{eqnarray}
where $m$ is given by equation (\ref{mass}).
The first of these equations is the algebraic representation of equation
(\ref{schr}) and the second one -- of the condition (\ref{cond1})
which determines $m$. The equations may be solved in the following way.
first, by the substitution $\rho=mr$ and by the subsequent application
of the Rayleigh-Ritz variational principle one transforms it into
the algebraic eigenvalue problem
\begin{equation}
\left(\mbox{\boldmath$T$}+
\mbox{\boldmath$V$}-E\mbox{\boldmath$S$}\right)
\mbox{\boldmath$c$}=0.
\label{alkb3}
\end{equation}
The energy $E$ may be minimized also with respect to the nonlinear
parameters (an alternative to this procedure may be a minimization of
\[
E=\left(\mbox{\boldmath$T$}+\mbox{\boldmath$V$}\right)
\mbox{\boldmath$S$}^{-1}
\] 
simultaneously with respect to all parameters,
as it has been done in references \cite{tat1b, tat2b, tat3b}.
Then the parameter $m$ may be determined from equation (\ref{alkb1}).
The large component of the contracted Dirac basis set function
is given by the solution of eqution (\ref{alkb3}) and the small component
by equation (\ref{kb}). The values of $C^{\rm L}$ and $C^{\rm S}$
are determined by equation (\ref{mass}) and by the normalization condition
\begin{equation}
\mbox{\boldmath$S$}|C^{\rm L}|^2+2\mbox{\boldmath$T$}|C^{\rm S}|^2=0.
\label{norm}
\end{equation}

A generalization of this procedure to the case of several contracted
basis functions is straightforward. An extension to the case of a
$\phi$-dependent potentials (as it is in the case of Dirac-Fock equations)
will be subject of a subsequent publication \cite{htjk}.

\section{Examples}
In this section we compare the ``exact'' solutions of equations
(\ref{dirackb}) and the algebraic solutions of equation (\ref{kbmatrix})
with the exact Dirac wavefunctions of the $1s_{1/2}$ states of hydrogen-like
atoms. The exact Dirac radial wavefunction reads
\begin{equation}
\Phi^{\rm{exact}}_{\rm D}=(2Z)^{\gamma+1/2}(2\Gamma(2\gamma+1))^{-1/2}
\left(
\begin{array}{c} 
(1+\gamma)^{1/2}\\[10pt] 
-(1-\gamma)^{1/2}
\end{array} 
\right)r^\gamma e^{-Zr},
\end{equation}
where $\gamma=\sqrt{1-\alpha^2Z^2}$. 
The radial part of the ``exact'' solution of the kinetically-balanced 
Dirac equation is equal to
\begin{equation}
\Phi^{\rm{exact}}_{\rm kb}=2(Zm)^{3/2}(4+\alpha^2Z^2)^{-1/2}
\left(
\begin{array}{c} 
2\\[10pt] 
-\alpha Z
\end{array} 
\right)r e^{-Zmr},
\end{equation}
where $m=1/(1-\alpha^2Z^2/4)$. 
The ground-state energy obtained as the solution of the $2\times 2$ 
secular problem corresponding to equations (\ref{dirackb}) is equal to
\begin{equation}
E = \langle V_+\rangle - c^2 +\left|\langle V_-\rangle+c^2\right|
\sqrt{1+\frac{2\langle T\rangle}{c^2}
\left(1+\frac{\langle V_-\rangle}{c^2}\right)^{-2}},
\label{energy1}
\end{equation} 
where 
\begin{equation}
\langle V_\pm\rangle =
\frac{1}{2}\left(\frac{V}{S}\pm\frac{X}{2T}\right),\;\mbox{and}\;\;
\langle T\rangle=\frac{T}{S}.
\end{equation}
The resulting wavefunction obtained in the Slater-type minimum-basis-set
is 
\begin{equation}
\Phi^{\rm{se}}_{\rm kb}=(Z\gamma)^{3/2}[(1+\gamma)/2]^{1/2}
\left(
\begin{array}{c} 
2\\[10pt] 
-2\alpha Z/(1+\gamma)
\end{array} 
\right)r e^{-Zr\gamma}.
\end{equation}
In figure 1 the behaviour of overlap integrals 
$S_1=\langle\Phi^{\rm{exact}}_{\rm D}|\Phi^{\rm{exact}}_{\rm kb}\rangle$ and 
$S_2=\langle\Phi^{\rm{exact}}_{\rm D}|\Phi^{\rm{se}}_{\rm kb}\rangle$
versus $Z$ is displayed. As one can see, in all cases the wavefunction
obtained as a solution of the Schr\"odinger-type equation gives a 
better approximation. 

In order to check the quality of the wavefunctions constructed as
expansion in a Gaussian-type basis, we have selected the large-component
basis sets as
\begin{equation}
\phi_k(r)^{\rm GL}=r^{n_k}\rm e^{-a_kr^2},\qquad, k=1,2,\ldots,N
\label{lgauss}
\end{equation}
with $n_k=1$ and $N=1,2,3,4,5$. Then the small-component basis generated
by the kinetic balance condition comprises functions
\begin{equation}
\phi_k(r)^{\rm GS}=r^{n_k+1}\rme^{-a_kr^2},\qquad, k=1,2,\ldots,N
\label{sgauss}
\end{equation}
The energy derived from the $2\times 2$ secular equation is given by
equation (\ref{energy1}). In the case of the Schr\"odinger-type equation
it is equal to
\begin{equation}
E=\langle T\rangle+\langle V_+\rangle+\langle V_-\rangle.
\end{equation}
If $\Omega$ stands for any of the matrix elements $V$, $T$, $X$, $S$ and
${\hat \Omega}$ --- for the corresponding operator, i.e., respectively, for
$U(r)$, $\frac{1}{2}{\hat{\cal O}}^\dag{\hat{\cal O}}$, ${\hat{\cal O}}^\dag
\,U(r)\,{\hat{\cal O}}$, $I$ then, in the contracted basis, $\Omega$ may be
expressed as a combination of two-orbital contributions:
\begin{equation}
\Omega=\sum_{i<j}^NR_{ij}
\end{equation}
where
\begin{equation}
R_{ij}={c_i}^2\frac{\omega_{ii}}{N-1}
\left[1+(N-1)f_{ji}r_1(q_{ji})+f_{ji}^2r_2(q_{ji})\right],
\end{equation}
$f_{ji}=\frac{c_j}{c_i}$, $q_{ji}=\frac{a_j}{a_i}$, $\omega_{ji}=\langle
\phi_j|{\hat\Omega}|\phi_i\rangle$. Explicit expressions for $\omega_{ji}$,
and for $r_p(q)$, $q=1,2$, may readily be obtained by a direct integration.
The resulting formulae, for the basis functions defined in equation
(\ref{lgauss}) with $n_k=1$ are collected in table 1.
\begin{table}
\caption{Expressions defining matrix elements of the algebraic 
form of the kinetically balanced Dirac equation for the $1s_{1/2}$-type
contracted Gaussian basis set. For explanation of symbols see text}
\begin{indented}
\item[]
\begin{tabular}{@{}cccc}
\br
$\Omega$&$\omega_{ij}$&$r_1(q)$&$r_2(q)$\\
\mr
$V$&$-\frac{Z}{2(a_i+a_j)}$&$\frac{4}{1+q}$&$\frac{1}{q}$\\[8pt] 
$T$&$\frac{3\sqrt{\pi}\;a_ia_j}{2(a_i+a_j)^{5/2}}$
&$\frac{8\sqrt{2}\;q}{(1+q)^{5/2}}$&$\frac{1}{\sqrt{q}}$\\[8pt]
$X$&$-\frac{2Z\;a_ia_j}{(a_i+a_j)^2}$&$\frac{8\;q}{(1+q)^2}$&$1$\\[8pt]
$S$&$\frac{\sqrt{\pi}}{4(a_i+a_j)^{3/2}}$
&$\frac{4\sqrt{2}}{(1+q)^{3/2}}$&$\frac{1}{q^{3/2}}$\\
\br
\end{tabular}
\end{indented}
\end{table}
Using these equations the energies may be expressed as functions 
of the parameters defining the basis sets. The energy minimization
with respect to these parameters has been performed using the Powell's
method \cite{powell}. The overlap integrals between the resulting
wavefunctions and the exact one are compared in figure 2. 

\section{Conclusions} 
We have noticed a possibility of an easy generation of basis sets
suitable for relativistic calculations from the nonrelativistic 
ones. The method consists of the following steps
\begin{itemize}
\item
Determining $\Phi(\rho)$ by minimizing the Rayleigh quotient resulting
from equation (\ref{schr}).
\item
Determining $m$ using equation (\ref{cond2}).
\item
Rescaling $\Phi$ by substitution $\rho=mr$.
\item
Generation of the small component using equation (\ref{kb}).
\item
Normalization of the wavefunction using equations (\ref{mass}) and
(\ref{norm}) 
\end{itemize}
An extension of this approach to the case of many-electron
systems described by Dirack-Fock equations is also feasible
and will be subject of a separate publication \cite{htjk}.

\ack
This work has been supported by a Grant-in-Aid from the Ministry of
Education, Science and Culture of Japan and by the Polish KBN, project No. 2
P03B 126 14. J. K. thanks to Nagoya City University for a visiting-scientist
fellowship which allowed him for a three-month stay at the Institute of
Natural Sciences, where this work has been executed, and to Professor
Hiroshi Tatewaki for his kind hospitality.

\section*{References}
\begin{thebibliography}{99}
\bibitem{hall} 
Hall G 1950 
\bibitem{roothaan} 
Roothaan C 1950
\bibitem{wilson} 
Wilson S 199x 
\bibitem{r1}
review on basis sets 1
\bibitem{r2}
review on basis sets 2
\bibitem{t1}
Huzinaga et al 198x
\bibitem{t2}
Mezey 198x
\bibitem{t3}
basis set tabulation 3
\bibitem{tahu}
Tatewaki - Huzinaga method
\bibitem{tat1b}
Tatewaki H and Koga T 1996 \JCP {\bf 104} 8493
\bibitem{tat2b}
Tatewaki H, Koga T and Takashima H 1997 {\em Theor. Chem. Acc.} 
{\bf 96} 243  
\bibitem{tat3b}
Koga T, Tatewaki H, Matsuyama H and Satoh Y 1999 {\em Theor. Chem. Acc.} 
{\bf 102} 105 
\bibitem{malli} Malli G and Oreg  
\bibitem{grant1}
Dyall K G, Grant I P and Wilson S 1984 \JPB {\bf 17} L45, 1201
\bibitem{grant2}
Grant I P 1986 \JPB {\bf 19} 3187
\bibitem{kutz1} Kutzelnigg W 1984 {\em Int. J. Quantum Chem.} {\bf 25} 107
\bibitem{ht}
Tatewaki H {\it unpublished results}
\bibitem{talman}
Talman J D 1986 \PRL {\bf 57} 1091
\bibitem{kutzmm}
Kutzelnigg W 1997 {\em Chem. Phys.} {\bf 225} 203
\bibitem{htjk}
Tatewaki H and Karwowski J {\it in preparation}
\bibitem{powell}
Powell M J D 1964 {\em Comput. J.} {\bf 7} 155. 
\end{thebibliography}
\end{document}



The properties of the uncontracted case were subject of extensive studies
reported in several papers \cite{grant1, grant2}. It has been demonstrated
that the algebraic problem defined in equation (\ref{dirac2}) gives a
correct description of the bound states located in the energy region between
$0$ and $-2c^2$ if the basis set exponents are chosen in such a way that
\begin{equation}
v_{min}=\,\mbox{ \raisebox{-2mm}{$\stackrel{\displaystyle
\inf}{\scriptscriptstyle \{\bf C\}}$}}
\frac{({\bf C}^\dag{\bf V}{\bf C})}{({\bf C}^\dag{\bf S}{\bf C})}
\ge -2c^2.
\end{equation}



