Hydrogenic Computation
using JenaAtomicCalculator
Perhaps, the simplest atomic computations can be made for hydrogenic ions. If we consider, for example, hydrogen-like argon $(Z=18)$, we can first determine and compare the energies of the $1s$ and $2p$ levels from the (non-relativistic) Schroedinger equation with those from Dirac's relativistic equation by
Z = 18.0
e1s = HydrogenicIon.energy(Shell("1s"), Z)
-162.0
e2p = HydrogenicIon.energy(Shell("2p"), Z)
e2p_1 = HydrogenicIon.energy(Subshell("2p_1/2"), Z)
e2p_3 = HydrogenicIon.energy(Subshell("2p_3/2"), Z)
-40.54376720995463
While the (one-electron) energies are displayed by the function HydrogenicIon.energy()
in the default units (presently eV here and as could be overwritten by the user), all computations are internally performed and returned always in atomic units. This applies to all computations if not stated otherwise in the description of some particular function; indeed, the user-defined units mainly in the tabulation of results but are not returned. Most of these tabulations are generated by display method which print some table to screen but return nothing
otherwise.
This clear distinction between the use of atomic and (user-specified) default units can be seen easily from the output of the variables e2p
, e2p_1
, e2p_3
above:
(e2p, e2p_1, e2p_3)
(-40.5, -40.720363843040175, -40.54376720995463)
The units of energies, rates and several other physical properties can be however quite easily converted among each other, for instance, by:
e1s_eV = convertUnits("energy: from atomic to eV", e1s)
e1s_au = convertUnits("energy: from eV to atomic", e1s_eV)
(e1s, e1s_eV, e1s_au)
(-162.0, -4408.24453524, -162.0)
... and similarly also for other energy units as well as for other physical entities; cf. ? convertUnits()
.
From these energies, we can easily compute the fine-structure splitting of the $2p$ level into the $2p_{1/2}$ and $2p_{3/2}$ (fine-structure) levels for hydrogen-like argon from above:
e13 = e2p_1 - e2p_3
e13_eV = convertUnits("energy: from atomic to eV", e13)
-4.805439152723064
Apart from the (single-electron) energies, we can generate also the radial orbitals, i.e. the $P(r)$ in the non-relativistic theory or $P(r)$ and $Q(r)$ in the relativistic theory. Both, the relativistic and non-relativistic orbitals, can be calculated by using the general (and analytically well-known) solutions of either the Dirac's or Schrödinger's equation.
However, since all (radial) orbital functions are represented on some Radial.Grid
, special care has to be taken in order to define a proper grid, and which affects also all subsequent numerical computations as well as the accuracy that can be obtained. To better understand the definition of the grid, let us ask for:
? Radial.Grid
JenaAtomicCalculator.Radial.Grid
— Typestruct Radial.Grid
... defines a type for the radial grid which contains all information about the grid parameters, the genration of the B-spline basis as well as for performing radial integrations.
** Physical grid parameter **
+ rnt ::Float64 ... smalles grid point > 0.
+ h ::Float64 ... stepsize in the construction of the exponential grid.
+ hp ::Float64 ... asymptotic stepsize of the log-lin grid.
+ NoPoints ::Int64 ... No. of grid points so that r[NoPoints] coincides also
with the largest break point of the B-spline knot.
** B-spline grid parameters and break points **
+ tL ::Array{Float64,1} ... radial break points for the B-splines of the large c.
+ tS ::Array{Float64,1} ... radial break points for the B-splines of the small c.
+ ntL ::Int64 ... number of break points in the t-grid of the large c.
+ ntS ::Int64 ... number of break points in the t-grid of the small c.
+ orderL ::Int64 ... B-spline order of large components.
+ orderS ::Int64 ... B-spline order of small components.
+ nsL ::Int64 ... number of B-splines for large components.
+ nsS ::Int64 ... number of B-splines for small components.
+ orderGL ::Int64
... order of the Gauss-Lengedre integration if mesh == Radial.MeshGL(); this order also determines
the (number of) break points by taking the orderGL-th point from the physical grid points.
** Radial mesh points **
+ meshType ::Radial.AbstractMesh
+ r ::Array{Float64,1} ... radial grid points
+ rp ::Array{Float64,1} ... derivative of the radial grid at the grid points
+ rpor ::Array{Float64,1} ... rp over r
+ wr ::Array{Float64,1}
... integration weights for all grid points, for instance, GL weights.
which shows the internal definition and the constructors to define such a grid. In practice, there are three (logical) levels that are considered and realized for each grid:
(i) The phenomenological grid definition that is based on our physics understanding and intuition of which grid is appropriate to decribe a particular atomic property or process. Two currently implemented grid refers to an exponential grid (with exponentially increasing stepsize) as appropriate for many bound-state computations and a linear-logarithmic grid that start exponentially but becomes linear in its stepsize for large values of $r$. Such a linear-logarithmic grid is typically needed to describe electron ionization and capture processes as the incoming or outgoing electron waves are sinusoidal and should be represented with a proper number of grid points, say 20-30, per period of the electron wave. In JAC, the phenomenologically part of the grid is chosen (analogue as in GRASP) by the parameters rnt
, h
, hp
and NoPoints
in the definition above.
(ii) The phenomenological definition of the grid is internally translated into a sequence of knots upon which the B-spline primitives are defined. Here, each nth
mesh point from the phenomenological is chosen to keep the number of B-splines and the corresponding size of matrices (that need to be diagonalized) moderate. These knots are used to define the B-splines and to determine the eigenvectors of all orbitals (within some given potential) but they are not applied in order evaluate matrix elements or radial integrals.
(iii) The phenomenological grid and the definition of the B-splines (knots) are eventually combined into a physical grid upon which all radial functions are represented. This physical grid resembles the phenomenological grid but with modified grid points in between the knots of the B-spline grid (t-grid). Here the radial points and corresponding weights are chosen due to a Gauss-Legendre distribution and goal to determine all radial integrals exact up to a given Gauss-Legendre order.
The physical grid is thus defined by the three arrays r
(the mesh points along $r$), rp
(the derivatives d r
/ dr), rpor
(the values rp / r) as well as w
(the corresponding weights). All these arrays are of length nr
$\approx$ NoPoints
but not necessarely equal as the grid points nr
are coupled to the underlying integration, respectively, interpolation scheme.
In principle, this physical grid could chosen also on other interpolation/integration formulas, such as Gauss-Laguerre or others, that have been utilized in atomic physics. All what is needed would be to adapt the intermediate grid point and weights accordingly.
The clear distinction between the phenomenological grid and the physical grid help avoid that every new (radial) operator as well as particular boundary condition of the radial orbital functions need to be treated independently within the B-spline basis. Instead, a proper interpolation/integration formula should guarantee that all results are integrated sufficiently accurate, and this is first of all tested by enlarging the number of grid point NoPoint
–> nr
Here, we first apply an exponential grid, and which is appropriate below in order to compute various expectation values:
grid = Radial.Grid(true, printout= true)
Radial grid: rnt = 2.0e-6, h = 0.05, hp = 0.0, NoPoints = 392, ntL = 69, ntS = 71, orderL = 7, orderS = 8, nsL = 62, nsS = 63, mesh = JenaAtomicCalculator.Radial.MeshGL(), ...
r: [1.780504510281598e-8, 9.042759120177923e-8, 2.0787030685460475e-7] ... [564.4302351890452, 595.097239730825, 614.0606632064549]
wr: [4.530145586285996e-8, 9.785739473901763e-8, 1.3358660616640612e-7] ... [34.882547080539766, 25.552825070730517, 11.829256033233948]
tS: [0.0, 0.0, 0.0] ... [618.7099715607404, 618.7099715607404, 618.7099715607404]
With this grid, the non-relativistic radial orbitals $P(r)$ can be obtained either for a single $r$-value, for a list of $r$-values as well as for all $r$-values on a given radial grid:
Pnr_1s = HydrogenicIon.radialOrbital(Shell("1s"), Z, grid)
Pnr_2p = HydrogenicIon.radialOrbital(Shell("2p"), Z, grid)
392-element Vector{Float64}:
8.895323125542194e-14
2.2944426826828984e-12
1.2124385701573003e-11
3.434471347946261e-11
6.787882523004707e-11
1.0416489853450738e-10
1.3047592028166753e-10
1.4891535552172264e-10
2.0078900763041305e-10
3.010512642484143e-10
⋮
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
Here, the exponential tails of the radial orbitals are simply set to zero if $|P| < 1.0e-15$ (and similarly for $|Q|$).
Of course, we could plot the two functions Pnr_1s
and Pnr_2p
directly by some proper call of PyPlot
, and if we give the correct arrays of radial mesh points from above.
using Plots
plot(grid.r, [Pnr_1s, Pnr_2p],xlims=(0,1), xlabel = "r (a.u.)", ylabel="Amplitude (a.u.)", label=["1s" "2p"])
WARNING: using Plots.grid in module Main conflicts with an existing identifier.
GKS: cannot open display - headless operation mode active
We can also make use of this function to visualise the relativistic orbitals if we specify the total angular momentum and generate the relativistic orbitals, which contain both the large and the small component.
setDefaults("standard grid", grid)
PQr_1s = HydrogenicIon.radialOrbital(Subshell("1s_1/2"), Z, grid)
PQr_2p1 = HydrogenicIon.radialOrbital(Subshell("2p_1/2"), Z, grid)
PQr_2p3 = HydrogenicIon.radialOrbital(Subshell("2p_3/2"), Z, grid)
plot(grid.r, [PQr_1s,PQr_2p1, PQr_2p3],xlims=(0,1), xlabel = "r (a.u.)", ylabel="Amplitude (a.u.)",
label=["1s - Large component" "1s - Small component" "2p_1/2 - Large component" "2p_1/2 - Small component" "2p_3/2 - Large component" "2p_3/2 - Small component"])
(Re-) Define the standard grid with 392 grid points.
In practice, using Pyplot is often quite tedious since all the specification of arrays, axes, labels, etc. need to be given explicitly. To display the shape of such orbitals, a special plotting method is prepared in JAC that displays one or several orbitals together.
? Basics.plot
JenaAtomicCalculator.Basics.plot
— FunctionBasics.plot()
... plots various quantities, often in a new window.
("radial potentials", potentials::Array{Radial.Potential,1}, grid::Radial.Grid; N::Int64 = 0)
... to plot one or more radial potentials, and where N::Int64 describes the number of grid points to be considered. call: using Plots; pyplot() ... to access this method by plot(...)
("radial orbitals: large", orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N::Int64 = 0)
... to plot the large component of one or more radial orbitals, and where N::Int64 describes the number of grid points to be considered.("radial orbitals: small", orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N::Int64 = 0)
... to plot the small component of one or more radial orbitals.("radial orbitals: both", orbitals::Array{Radial.Orbital,1}, grid::Radial.Grid; N::Int64 = 0)
... to plot the large and small component of one or more radial orbitals.call: using Plots; pyplot() ... to access this method by plot(...)
("spectrum: transition rates over energy", lines::Array{PhotoEmission.Line,1})
... to plot the transition rates of all lines as function of their transition energies. The plot is shown in a new window but nothing is returned otherwise. Not yet implemented !("spectrum: oscillator strength over energy, emission", lines::Array{PhotoEmission.Line,1})
or("spectrum: oscillator strength over energy, absorption", lines::Array{PhotoEmission.Line,1})
... to plot the absorption oscillator strength of all lines as function of their transition energies. Again, a new window is opened but nothing returned by this method. Not yet implemented !
("spectrum: transition rates over energy, Gaussian", lines::Array{PhotoEmission.Line,1}; widths=value::Float64)
or("spectrum: transition rates over energy, Lorentzian", lines::Array{PhotoEmission.Line,1}; widths=value::Float64)
... to plot the transition rates of all lines as function of their transition energies but with a Gaussian or Lorentzian distribution. Again, a new window is opened but nothing returned by this method. It still need to be decided how the widths (and, perhaps, other parameters) will be communicated to the method. Not yet implemented !
This Basics.plot
function (from the Basics module) accepts for instance a list of (relativistic) radial Orbitals
:
? Orbital
JenaAtomicCalculator.Radial.Orbital
— Typestruct Radial.Orbital
... defines a type for a single-electron radial orbital function with a large and small component, and which can refer to either the standard or an explicitly given grid due to the logical flag useStandardGrid. Bound-state orbitals with energy < 0 are distinguished from free-electron orbitals by the flag isBound.
+ subshell ::Subshell ... Relativistic subshell.
+ isBound ::Bool ... Logical flag to distinguish between bound (true) and free-electron orbitals (false).
+ useStandardGrid ::Bool ... Logical flag for using the standard grid (true) or an explicitly given grid (false).
+ energy ::Float64 ... Single-electron energies of bound orbitals are always negative.
+ P ::Array{Float64,1} ... Large and ..
+ Q ::Array{Float64,1} ... small component of the radial orbital.
+ Pprime ::Array{Float64,1} ... dP/dr.
+ Qprime ::Array{Float64,1} ... dQ/dr.
+ grid ::Array{Float64,1} ... explic. defined radial grid array for P, Q, if StandardGrid = false.
Of course, the main purpose of this plot
method is to display and compare a number of generated orbitals, either for getting physical insight or just for test purposes. We can make use of the same plot
method also for the non-relativistic Pnr_1s
$(r)$ orbital if we assign this array as large component of a relativistic orbital and specify the additional quantum numbers. Moreover, we here first internally to define a standard grid as well as such an (instance of the) Orbital by using its standard constructor, and by setting the small $Q(r)$ component as well as $dP/dr$ and $dQ/dr$ simply to zero:
setDefaults("standard grid", grid)
Qnr_1s = Pprime = Qprime = zeros( length(Pnr_1s) )
nrOrb_1s = Orbital( Subshell("1s_1/2"), true, true, e1s, Pnr_1s, Qnr_1s, Pprime, Qprime, Radial.Grid())
typeof(nrOrb_1s)
Orbital
We could also compute a relativistic $1s_{1/2}$ orbital on the given grid and plot the two $1s$-orbitals together; while the large components of these orbitals coincide of course, due to its special construction via the non-relativistic orbitals, the small compoments differ. Note that one need to use the standard REPL in order to display these functions in some separate window.
Orb_1s = HydrogenicIon.radialOrbital(Subshell("1s_1/2"), Z, grid)
e1s_1 = HydrogenicIon.energy(Subshell("1s_1/2"), Z)
rOrb_1s = Orbital( Subshell("1s_1/2"), true, true, e1s_1, Orb_1s[1], Orb_1s[2], Pprime, Qprime, Radial.Grid())
Bound-state orbital 1s_1/2 with energy -162.70485811496013 a.u. is defined with 392 (grid) points on the standard grid:
Large component P: [3.094976712530428e-6, 1.5498862894779233e-5, 3.537189759693442e-5, 5.926504419660713e-5, 8.30717237172039e-5, 0.00010271653265137339, 0.00011484727355055586, 0.0001226242231422872, 0.00014220454983176329, 0.00017382034714231485 … 0.0003400832988479155, 0.0003846215960535446, 0.0004384081093685234, 0.0004921362282859167, 0.0005365367822203975, 0.0005639761006343588, 0.0005802444216048941, 0.000619152996620868, 0.0006820292546595882, 0.0007579783496550785] ... [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Small component Q: [2.0415063278116693e-7, 1.0223348868983198e-6, 2.3331985820277163e-6, 3.9092366109110816e-6, 5.479571104502363e-6, 6.7753805878302265e-6, 7.5755476523015185e-6, 8.088530245621624e-6, 9.380086355732164e-6, 1.1465525319036404e-5 … 2.243255026023918e-5, 2.5370382238332118e-5, 2.8918244386660844e-5, 3.246225472790821e-5, 3.5391000894191195e-5, 3.720095125492078e-5, 3.827404108043568e-5, 4.084052572568854e-5, 4.498796496603807e-5, 4.999771374369833e-5] ... [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Pprime: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ... [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Qprime: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] ... [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Defined on Grid: [1.780504510281598e-8, 9.042759120177923e-8, 2.0787030685460475e-7] ... [527.353759459587, 564.4302351890452, 595.097239730825, 614.0606632064549]
plot("radial orbitals: both", Orbital[nrOrb_1s, rOrb_1s], grid; N = 230)
We can use these orbitals to quickly evaluate some overlap integrals and $\langle r^k \rangle$ expectation values; for this, we make use of the non-relativistic $2p$ radial orbital Pnr_2p
from above and compute the normalization and $\langle r \rangle$ expectation values:
N_1s = RadialIntegrals.overlap(Pnr_1s, Pnr_1s, grid)
N_2p = RadialIntegrals.overlap(Pnr_2p, Pnr_2p, grid)
r_exp_1s = RadialIntegrals.rkDiagonal(1, Pnr_1s, Pnr_1s,grid)
r_exp_2p = RadialIntegrals.rkDiagonal(1, Pnr_2p, Pnr_2p,grid)
println("N_1s = $(N_1s)","\n","N_2p = $(N_2p)","\n","r_exp_1s = $(r_exp_1s)","\n","r_exp_2p = $(r_exp_2p)")
N_1s = 0.999999999999998
N_2p = 0.9999999999999959
r_exp_1s = 0.08333333333333318
r_exp_2p = 0.2777777777777762
As we can see the radial orbitals are already normalized by construction and the expectation values match with the analytic values of $\langle r \rangle_{nl} = \frac{1}{2Z}[3n²-l(l+1)]$.
For the sake of simplicity, all hydrogenic computations are performed for a point-like nucleus. For many-electron computations, in constrast, a more realistic nucleus can be taken into accout by selecting a proper nuclear model. See the tutorial: Define the nuclear model parameters or ? Nuclear.Model
... for further details.
? Nuclear.Model
JenaAtomicCalculator.Nuclear.Model
— Typestruct Nuclear.Model
... defines a type for the nuclear model, i.e. for its form and parameters.
+ Z ::Float64 ... nuclear charge
+ model ::String ... identifier of the nuclear model: {"Fermi", "Point", "Uniform"}
+ mass ::Float64 ... atomic mass
+ radius ::Float64 ... (root-mean square) radius of a uniform or Fermi-distributed nucleus
+ spinI ::AngularJ64 ... nuclear spin I, must be >= 0
+ mu ::Float64 ... magnetic dipole moment in Bohr magnetons
+ Q ::Float64 ... electric quadrupole moment
If we would like to compute a relativistic orbital by taking into account the type of the nucleus, we can simply define a nuclear model:
Nucleus = Model(Z, "Fermi")
Fermi nuclear model for Z = 18.0 with mass = 37.62, radius R = 3.37121133160794 fm and nuclear spin I = 0, dipole moment mu = 0.0 and quadrupole moment Q = 0.0.
and set it as an argument value of the radialOrbital-function. By comparing the analytic orbitals which are based on a point like nucleus and the orbitals taking into account a fermi distributed nucleus one can see, that the influence of the nucleus is negligible for large r. However by looking at the functions at a closer angle the influence of the nucleus is well visible for small r.
nuc_orb_1s = HydrogenicIon.radialOrbital(Subshell("1s_1/2"),Nucleus, grid)
plot("radial orbitals: both", [nuc_orb_1s,rOrb_1s], grid; N = 230)
savefig("hdy-4.svg"); nothing # hide
plot("radial orbitals: both", [nuc_orb_1s,rOrb_1s], grid; N = 75)
savefig("hdy-2.svg"); nothing # hide
Compute the low-lying levels of $C^{2+} 1s^2 (2s^2 + 2s2p + 2p^2)$: $\;$ SCF and configuration interaction calculations
using JenaAtomicCalculator
The low-lying levels (level structure) of beryllium-like ions, and especially of C$^{2+}$, has been calculated in many case studies in the literature. While the level structure of these ions is still quite simple, it exhibits a considerable admixture of the $2s^22p^2$ configuration already for the $1s^{2}2s^{2}$ $^{1}S_{0}$ ground level.
We here show how the low-lying levels of C$^{2+}$ can be readily calculated in JAC by either following the default settings or by specifying further details for both, the SCF and configuration-interaction (CI) computations. As usual, we first need to specify a radial grid as well as the nuclear model for the subsequent computations:
grid = Radial.Grid(true)
nucModel = Nuclear.Model(6., "Fermi")
Fermi nuclear model for Z = 6.0 with mass = 12.18, radius R = 2.4934845703596404 fm and nuclear spin I = 0, dipole moment mu = 0.0 and quadrupole moment Q = 0.0.
For a quick computation of the ground level of C$^{2+}$ ions, we can simply use the standard settings as given by AsfSettings()
:
multiplet = SelfConsistent.performSCF([Configuration("1s^2 2s^2")], nucModel, grid, AsfSettings());
>>> include Configuration: 1s_1/2^2 2s_1/2^2
(Re-) Define a new standard subshell list.
> Start SCF process with hydrogenic orbitals.
>> (Re-) Define a storage array for dealing with single-electron TTp B-spline matrices:
-----------------------------------------------------------------------------
Index Subshell Energies [a.u.] Dirac-E [a.u.] Delta-E / |E|
-----------------------------------------------------------------------------
1 1s_1/2 -1.80086331e+01 -1.80086350e+01 +1.03994916e-07
2 2s_1/2 -4.50269964e+00 -4.50269857e+00 -2.38959499e-07
3 3s_1/2 -2.00098006e+00 -2.00095940e+00 -1.03240049e-05
4 4s_1/2 -1.12562939e+00 -1.12543844e+00 -1.69632831e-04
5 5s_1/2 -7.21116908e-01 -7.20234830e-01 -1.22321194e-03
6 6s_1/2 -5.29936446e-01 -5.00139888e-01 -5.62266634e-02
7 7s_1/2 -4.90049428e-01 -3.67436826e-01 -2.50204560e-01
: :
61 61s_1/2 +1.07252829e+09 -4.83755953e-03 +1.00000000e+00
62 62s_1/2 +2.82418557e+09 -4.68276538e-03 +1.00000000e+00
-----------------------------------------------------------------------------
>> (Re-) Define a storage array for dealing with single-electron TTp B-spline matrices:
Iteration 1 for symmetries ...
1s_1/2:: en [a.u.] = -1.0316230e+01; self-cons'cy = 2.7158e-01 [1.0000e+02 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.2238995e+00; self-cons'cy = 5.7256e-01 [1.0000e+02 for sym-block kappa = -1]
Iteration 2 for symmetries ...
1s_1/2:: en [a.u.] = -1.1298727e+01; self-cons'cy = 4.5454e-02 [1.4446e-01 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4578099e+00; self-cons'cy = 8.7224e-02 [1.4446e-01 for sym-block kappa = -1]
Iteration 3 for symmetries ...
1s_1/2:: en [a.u.] = -1.1135297e+01; self-cons'cy = 7.2849e-03 [2.9588e-02 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4168210e+00; self-cons'cy = 1.4259e-02 [2.9588e-02 for sym-block kappa = -1]
Iteration 4 for symmetries ...
1s_1/2:: en [a.u.] = -1.1162291e+01; self-cons'cy = 1.2106e-03 [4.6048e-03 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4235023e+00; self-cons'cy = 2.3523e-03 [4.6048e-03 for sym-block kappa = -1]
Iteration 5 for symmetries ...
1s_1/2:: en [a.u.] = -1.1157826e+01; self-cons'cy = 2.0004e-04 [7.6877e-04 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4223956e+00; self-cons'cy = 3.8887e-04 [7.6877e-04 for sym-block kappa = -1]
Iteration 6 for symmetries ...
1s_1/2:: en [a.u.] = -1.1158565e+01; self-cons'cy = 3.3083e-05 [1.2687e-04 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4225785e+00; self-cons'cy = 6.4293e-05 [1.2687e-04 for sym-block kappa = -1]
Iteration 7 for symmetries ...
1s_1/2:: en [a.u.] = -1.1158443e+01; self-cons'cy = 5.4705e-06 [2.0985e-05 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4225483e+00; self-cons'cy = 1.0631e-05 [2.0985e-05 for sym-block kappa = -1]
Iteration 8 for symmetries ...
1s_1/2:: en [a.u.] = -1.1158463e+01; self-cons'cy = 9.0461e-07 [3.4698e-06 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4225533e+00; self-cons'cy = 1.7579e-06 [3.4698e-06 for sym-block kappa = -1]
Iteration 9 for symmetries ...
1s_1/2:: en [a.u.] = -1.1158459e+01; self-cons'cy = 1.4959e-07 [5.7373e-07 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4225525e+00; self-cons'cy = 2.9069e-07 [5.7373e-07 for sym-block kappa = -1]
> Compute CI matrix of dimension 1 x 1 for the symmetry 0^+ ... ... done.
Eigenenergies:
Level J Parity Hartrees eV [eV]
1 0 + -3.641932007141357e+01 -9.910201770491686e+02 -9.910201770491686e+02
Energy of each level relative to immediately lower level:
Level J Parity Hartrees eV [eV]
Energy of each level relative to lowest level:
Level J Parity Hartrees eV [eV]
Eigenenergies:
Level J Parity Hartrees eV [eV]
1 0 + -3.641932007141357e+01 -9.910201770491686e+02 -9.910201770491686e+02
Energy of each level relative to lowest level:
Level J Parity Hartrees eV [eV]
... and, similarly, also for the 10 lowest levels of C$^{2+}$ ions as they arise from the $C^{2+} 1s^2 (2s^2 + 2s2p + 2p^2)$ configurations:
multiplet = SelfConsistent.performSCF([Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p"), Configuration("1s^2 2p^2")],
nucModel, grid, AsfSettings());
>>> include Configuration: 1s_1/2^2 2s_1/2^2
>>> include Configuration: 1s_1/2^2 2s_1/2^1 2p_1/2^0 2p_3/2^1
>>> include Configuration: 1s_1/2^2 2s_1/2^1 2p_1/2^1 2p_3/2^0
>>> include Configuration: 1s_1/2^2 2p_1/2^0 2p_3/2^2
>>> include Configuration: 1s_1/2^2 2p_1/2^1 2p_3/2^1
>>> include Configuration: 1s_1/2^2 2p_1/2^2 2p_3/2^0
(Re-) Define a new standard subshell list.
> Start SCF process with hydrogenic orbitals.
>> (Re-) Define a storage array for dealing with single-electron TTp B-spline matrices:
-----------------------------------------------------------------------------
Index Subshell Energies [a.u.] Dirac-E [a.u.] Delta-E / |E|
-----------------------------------------------------------------------------
1 1s_1/2 -1.80086331e+01 -1.80086350e+01 +1.03994916e-07
2 2s_1/2 -4.50269964e+00 -4.50269857e+00 -2.38959499e-07
3 3s_1/2 -2.00098006e+00 -2.00095940e+00 -1.03240049e-05
4 4s_1/2 -1.12562939e+00 -1.12543844e+00 -1.69632831e-04
5 5s_1/2 -7.21116908e-01 -7.20234830e-01 -1.22321194e-03
6 6s_1/2 -5.29936446e-01 -5.00139888e-01 -5.62266634e-02
7 7s_1/2 -4.90049428e-01 -3.67436826e-01 -2.50204560e-01
: :
61 61s_1/2 +1.07252829e+09 -4.83755953e-03 +1.00000000e+00
62 62s_1/2 +2.82418557e+09 -4.68276538e-03 +1.00000000e+00
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
Index Subshell Energies [a.u.] Dirac-E [a.u.] Delta-E / |E|
-----------------------------------------------------------------------------
1 2p_1/2 -4.50269888e+00 -4.50269857e+00 -6.91701387e-08
2 3p_1/2 -2.00096714e+00 -2.00095940e+00 -3.86903735e-06
3 4p_1/2 -1.12549933e+00 -1.12543844e+00 -5.41023530e-05
4 5p_1/2 -7.20400487e-01 -7.20234830e-01 -2.29951315e-04
5 6p_1/2 -5.00841849e-01 -5.00139888e-01 -1.40156319e-03
6 7p_1/2 -3.91479761e-01 -3.67436826e-01 -6.14155231e-02
7 8p_1/2 -3.49594545e-01 -2.81311119e-01 -1.95321771e-01
: :
61 62p_1/2 +1.31159793e+09 -4.68276538e-03 +1.00000000e+00
62 63p_1/2 +6.75311123e+09 -4.53528382e-03 +1.00000000e+00
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
Index Subshell Energies [a.u.] Dirac-E [a.u.] Delta-E / |E|
-----------------------------------------------------------------------------
1 2p_3/2 -4.50053984e+00 -4.50053930e+00 -1.19777072e-07
2 3p_3/2 -2.00033044e+00 -2.00031959e+00 -5.42476571e-06
3 4p_3/2 -1.12535222e+00 -1.12516853e+00 -1.63226680e-04
4 5p_3/2 -7.21810564e-01 -7.20096642e-01 -2.37447501e-03
5 6p_3/2 -5.88521709e-01 -5.00059921e-01 -1.50311851e-01
6 7p_3/2 -4.94086009e-01 -3.67386470e-01 -2.56432153e-01
7 8p_3/2 -3.59118075e-01 -2.81277385e-01 -2.16755142e-01
: :
61 62p_3/2 +1.75386123e+09 -4.68269292e-03 +1.00000000e+00
62 63p_3/2 +6.75822628e+09 -4.53521476e-03 +1.00000000e+00
-----------------------------------------------------------------------------
>> (Re-) Define a storage array for dealing with single-electron TTp B-spline matrices:
Iteration 1 for symmetries ...
1s_1/2:: en [a.u.] = -1.0025691e+01; self-cons'cy = 2.8476e-01 [1.0000e+02 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.1510995e+00; self-cons'cy = 5.9280e-01 [1.0000e+02 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -9.0244503e-01; self-cons'cy = 6.6608e-01 [1.0000e+02 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -9.7438984e-01; self-cons'cy = 6.4405e-01 [1.0000e+02 for sym-block kappa = -2]
Iteration 2 for symmetries ...
1s_1/2:: en [a.u.] = -1.1312283e+01; self-cons'cy = 6.0296e-02 [1.9039e-01 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4565873e+00; self-cons'cy = 1.1715e-01 [1.9039e-01 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.1346471e+00; self-cons'cy = 1.1399e-01 [3.6941e+00 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.1175834e+00; self-cons'cy = 6.8449e-02 [7.6589e-02 for sym-block kappa = -2]
Iteration 3 for symmetries ...
1s_1/2:: en [a.u.] = -1.1095187e+01; self-cons'cy = 9.6886e-03 [4.1453e-02 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4000354e+00; self-cons'cy = 1.9797e-02 [4.1453e-02 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.0939225e+00; self-cons'cy = 1.8274e-02 [1.1376e-01 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.0954890e+00; self-cons'cy = 9.9836e-03 [1.7762e-02 for sym-block kappa = -2]
Iteration 4 for symmetries ...
1s_1/2:: en [a.u.] = -1.1125825e+01; self-cons'cy = 1.3788e-03 [5.1060e-03 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4073480e+00; self-cons'cy = 2.6048e-03 [5.1060e-03 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.0985915e+00; self-cons'cy = 2.1295e-03 [1.7140e-02 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.0978026e+00; self-cons'cy = 1.0549e-03 [1.7598e-03 for sym-block kappa = -2]
Iteration 5 for symmetries ...
1s_1/2:: en [a.u.] = -1.1121783e+01; self-cons'cy = 1.8166e-04 [6.2051e-04 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4064550e+00; self-cons'cy = 3.1739e-04 [6.2051e-04 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.0980772e+00; self-cons'cy = 2.3413e-04 [2.0654e-03 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.0975647e+00; self-cons'cy = 1.0834e-04 [1.8483e-04 for sym-block kappa = -2]
Iteration 6 for symmetries ...
1s_1/2:: en [a.u.] = -1.1122289e+01; self-cons'cy = 2.2733e-05 [7.0320e-05 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4065576e+00; self-cons'cy = 3.6480e-05 [7.0320e-05 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.0981296e+00; self-cons'cy = 2.3863e-05 [2.1074e-04 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.0975873e+00; self-cons'cy = 1.0299e-05 [1.7791e-05 for sym-block kappa = -2]
Iteration 7 for symmetries ...
1s_1/2:: en [a.u.] = -1.1122228e+01; self-cons'cy = 2.7232e-06 [7.5613e-06 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4065464e+00; self-cons'cy = 3.9804e-06 [7.5613e-06 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.0981247e+00; self-cons'cy = 2.2246e-06 [1.9880e-05 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.0975854e+00; self-cons'cy = 8.7245e-07 [1.5454e-06 for sym-block kappa = -2]
Iteration 8 for symmetries ...
1s_1/2:: en [a.u.] = -1.1122235e+01; self-cons'cy = 3.1302e-07 [7.6455e-07 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4065475e+00; self-cons'cy = 4.1032e-07 [7.6455e-07 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.0981251e+00; self-cons'cy = 1.7992e-07 [1.6377e-06 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.0975855e+00; self-cons'cy = 5.8204e-08 [1.0933e-07 for sym-block kappa = -2]
Iteration 9 for symmetries ...
1s_1/2:: en [a.u.] = -1.1122235e+01; self-cons'cy = 3.4462e-08 [7.1250e-08 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.4065474e+00; self-cons'cy = 3.9328e-08 [7.1250e-08 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.0981251e+00; self-cons'cy = 1.0489e-08 [1.0104e-07 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.0975855e+00; self-cons'cy = 1.2988e-09 [6.4886e-09 for sym-block kappa = -2]
>>> Sign changed for orbital 2p_3/2
> Compute CI matrix of dimension 1 x 1 for the symmetry 0^- ... ... done.
> Compute CI matrix of dimension 1 x 1 for the symmetry 1^+ ... ... done.
> Compute CI matrix of dimension 2 x 2 for the symmetry 2^+ ... ... done.
> Compute CI matrix of dimension 1 x 1 for the symmetry 2^- ... ... done.
> Compute CI matrix of dimension 2 x 2 for the symmetry 1^- ... ... done.
> Compute CI matrix of dimension 3 x 3 for the symmetry 0^+ ... ... done.
Eigenenergies:
Level J Parity Hartrees eV [eV]
1 0 + -3.648596893065901e+01 -9.928337848858889e+02 -9.928337848858889e+02
2 0 - -3.624303134670400e+01 -9.862231165101230e+02 -9.862231165101230e+02
3 1 - -3.624287306106216e+01 -9.862188093384214e+02 -9.862188093384214e+02
4 2 - -3.624255595204814e+01 -9.862101803626306e+02 -9.862101803626306e+02
5 1 - -3.597168309684437e+01 -9.788393545373413e+02 -9.788393545373413e+02
6 0 + -3.585249805874401e+01 -9.755961644577839e+02 -9.755961644577839e+02
7 1 + -3.585234658730297e+01 -9.755920427099306e+02 -9.755920427099306e+02
8 2 + -3.585204525341215e+01 -9.755838429971067e+02 -9.755838429971067e+02
9 2 + -3.577601908774491e+01 -9.735150656555151e+02 -9.735150656555151e+02
10 0 + -3.559491143203561e+01 -9.685868753248319e+02 -9.685868753248319e+02
Energy of each level relative to immediately lower level:
Level J Parity Hartrees eV [eV]
2 0 - 2.429375839550119e-01 6.610668375765986e+00 6.610668375765986e+00
3 1 - 1.582856418451684e-04 4.307171701672341e-03 4.307171701672341e-03
4 2 - 3.171090140128285e-04 8.628975790724667e-03 8.628975790724667e-03
5 1 - 2.708728552037698e-01 7.370825825289345e+00 7.370825825289345e+00
6 0 + 1.191850381003619e-01 3.243190079557354e+00 3.243190079557354e+00
7 1 + 1.514714410433271e-04 4.121747853235646e-03 4.121747853235646e-03
8 2 + 3.013338908175456e-04 8.199712823944766e-03 8.199712823944766e-03
9 2 + 7.602616566723697e-02 2.068777341591656e+00 2.068777341591656e+00
10 0 + 1.811076557093045e-01 4.928190330683143e+00 4.928190330683143e+00
Energy of each level relative to lowest level:
Level J Parity Hartrees eV [eV]
2 0 - 2.429375839550119e-01 6.610668375765986e+00 6.610668375765986e+00
3 1 - 2.430958695968570e-01 6.614975547467658e+00 6.614975547467658e+00
4 2 - 2.434129786108699e-01 6.623604523258383e+00 6.623604523258383e+00
5 1 - 5.142858338146397e-01 1.399443034854773e+01 1.399443034854773e+01
6 0 + 6.334708719150015e-01 1.723762042810508e+01 1.723762042810508e+01
7 1 + 6.336223433560448e-01 1.724174217595832e+01 1.724174217595832e+01
8 2 + 6.339236772468624e-01 1.724994188878226e+01 1.724994188878226e+01
9 2 + 7.099498429140993e-01 1.931871923037392e+01 1.931871923037392e+01
10 0 + 8.910574986234039e-01 2.424690956105706e+01 2.424690956105706e+01
Eigenenergies:
Level J Parity Hartrees eV [eV]
1 0 + -3.648596893065901e+01 -9.928337848858889e+02 -9.928337848858889e+02
2 0 - -3.624303134670400e+01 -9.862231165101230e+02 -9.862231165101230e+02
3 1 - -3.624287306106216e+01 -9.862188093384214e+02 -9.862188093384214e+02
4 2 - -3.624255595204814e+01 -9.862101803626306e+02 -9.862101803626306e+02
5 1 - -3.597168309684437e+01 -9.788393545373413e+02 -9.788393545373413e+02
6 0 + -3.585249805874401e+01 -9.755961644577839e+02 -9.755961644577839e+02
7 1 + -3.585234658730297e+01 -9.755920427099306e+02 -9.755920427099306e+02
8 2 + -3.585204525341215e+01 -9.755838429971067e+02 -9.755838429971067e+02
9 2 + -3.577601908774491e+01 -9.735150656555151e+02 -9.735150656555151e+02
10 0 + -3.559491143203561e+01 -9.685868753248319e+02 -9.685868753248319e+02
Energy of each level relative to lowest level:
Level J Parity Hartrees eV [eV]
2 0 - 2.429375839550119e-01 6.610668375765986e+00 6.610668375765986e+00
3 1 - 2.430958695968570e-01 6.614975547467658e+00 6.614975547467658e+00
4 2 - 2.434129786108699e-01 6.623604523258383e+00 6.623604523258383e+00
5 1 - 5.142858338146397e-01 1.399443034854773e+01 1.399443034854773e+01
6 0 + 6.334708719150015e-01 1.723762042810508e+01 1.723762042810508e+01
7 1 + 6.336223433560448e-01 1.724174217595832e+01 1.724174217595832e+01
8 2 + 6.339236772468624e-01 1.724994188878226e+01 1.724994188878226e+01
9 2 + 7.099498429140993e-01 1.931871923037392e+01 1.931871923037392e+01
10 0 + 8.910574986234039e-01 2.424690956105706e+01 2.424690956105706e+01
From the comparison of the two ground-state energies, we see that the admixture of the $2p^2$ configuration has lowered the (total) ground state energy by about 1.8 eV, a rather remarkable admixture, as the $\: ^3P_0$ is just 6.6 eV above of the ground level.
Further control about these electronic computations can be obtained by modifying the (so-called) settings. In general, all computations of the electronic structure, properties and processes as well as all more advanced computations can be controlled quite in details by various settings that are associated to the different computational requests. The SCF and configuration interaction calculations are controlled by AsfSettings
that specify all details for the generation of the ASF. We can first have a look at the internal representation of these settings:
? AsfSettings
JenaAtomicCalculator.ManyElectron.AsfSettings
— Typestruct ManyElectron.AsfSettings
... a struct for defining the settings for the atomic state functions, i.e. the self-consistent-field (SCF) and CI computations
+ generateScf ::Bool ... True, if a SCF need to be generated, and false otherwise
(frozen orbitals).
+ eeInteraction ::AbstractEeInteraction ... Specify the e-e interaction to be included into the SCF
computations.
+ scField ::AbstractScField ... Specify the self-consistent field, for instance,
Basics.ALField(), etc.
+ startScfFrom ::AbstractStartOrbitals ... Specify the orbitals to start the SCF computations
+ maxIterationsScf ::Int64 ... maximum number of SCF iterations
+ accuracyScf ::Float64 ... convergence criterion for the SCF field.
+ shellSequenceScf ::Array{Subshell,1} ... Sequence of subshells to be optimized.
+ frozenSubshells ::Array{Subshell,1} ... Sequence of subshells to be kept frozen.
+ eeInteractionCI ::AbstractEeInteraction ... Specify the e-e interaction to be included into the
CI computations.
+ qedModel ::AbstractQedModel ... model for estimating QED corrections {NoneQed(),
QedPetersburg(), QedSydney()}.
+ jjLS ::LSjjSettings ... settings to control a jj-LS transformation of atomic
level, if requested.
+ levelSelectionCI ::LevelSelection ... Specifies the selected levels, if any.
... as well as how the standard settings are presently defined:
defaultAsfSettings = AsfSettings()
generateScf: true
eeInteraction: CoulombInteraction()
scField: JenaAtomicCalculator.Basics.DFSField(1.0)
startScfFrom: StartFromHydrogenic()
maxIterationsScf: 24
accuracyScf: 1.0e-6
shellSequenceScf: Subshell[]
frozenSubshells: Subshell[]
eeInteractionCI: CoulombInteraction()
qedModel : NoneQed()
jjLS: false
levelSelectionCI: Inactive LevelSelection.
From this list, we easily see that the self-consistent field is by default based on a (mean) Dirac-Fock-Slater potential, a choice which we could overwrite by meanCH
or any pre-defined potential. At present, however, no full treatment of the exchange interaction has yet been implemented in this first release of the program. The standard settings also show that the SCF is usually based on just the lowest level (cf. levelScf
) as well as on the accuracy 1.0e-6
to terminate the SCF computations. Moreover, the individual orbitals are improved due to the standard subshell order (cf. shellSequenceScf
).
For the CI parameters, the treatment of the Breit and QED interaction is of particular interest. At present, the defaults does not included neither Breit interactions nor QED. Such QED estimates can either be neglected (NoneQed()
) or estimated by using an effective Hamiltonian approach due to Shabaev and coworkers (QedPetersburg()
) or effective potential approach (QedSydney()
; Flambaum et al.) However, further tests need to be done to better understand the reliability of these QED estimates to the level structure and state represetation of the ASF.
As seen from the settings above, moreover, there are special features in order to select individual levels for the CI computations, either in terms of their (relative) level No within the given multiplet or in terms of their level symmetry, i.e. their total angular momentum and parity, respectively. The levelSelectionCI::LevelSelection
here tells whether (and which) selections were made; apparently, no selection of level numbers of symmetries is made by default though this can be overwritten. The selection of individual symmetries, in particular, may considerably reduced the computational effort as the Hamiltonian matrix need then to be calculated only for the selected symmetries.
In principle, these standard settings can be easily re-defined within the code by simply modifying the constructor AsfSetings()
with no additional arguments. Alternatively, we can easily overwrite those parameters in some given (instance of) AsfSetting
which we just wish to modify. This is achieved by
asfSettings = AsfSettings(defaultAsfSettings; generateScf=true, jjLS=LSjjSettings(true),
levelSelectionCI=LevelSelection(true, symmetries=[LevelSymmetry(0,"+"), LevelSymmetry(1,"-")]) )
generateScf: true
eeInteraction: CoulombInteraction()
scField: JenaAtomicCalculator.Basics.DFSField(1.0)
startScfFrom: StartFromHydrogenic()
maxIterationsScf: 24
accuracyScf: 1.0e-6
shellSequenceScf: Subshell[]
frozenSubshells: Subshell[]
eeInteractionCI: CoulombInteraction()
qedModel : NoneQed()
jjLS: true
levelSelectionCI: LevelSelection: indices = Int64[]; symmetries = LevelSymmetry[0 +, 1 -].
multiplet = SelfConsistent.performSCF([Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p"),
Configuration("1s^2 2p^2"), Configuration("1s^2 3s^2"),
Configuration("1s^2 3p^2")], nucModel, grid, asfSettings);
>>> include Configuration: 1s_1/2^2 2s_1/2^2
>>> include Configuration: 1s_1/2^2 2s_1/2^1 2p_1/2^0 2p_3/2^1
>>> include Configuration: 1s_1/2^2 2s_1/2^1 2p_1/2^1 2p_3/2^0
>>> include Configuration: 1s_1/2^2 2p_1/2^0 2p_3/2^2
>>> include Configuration: 1s_1/2^2 2p_1/2^1 2p_3/2^1
>>> include Configuration: 1s_1/2^2 2p_1/2^2 2p_3/2^0
>>> include Configuration: 1s_1/2^2 3s_1/2^2
>>> include Configuration: 1s_1/2^2 3p_1/2^0 3p_3/2^2
>>> include Configuration: 1s_1/2^2 3p_1/2^1 3p_3/2^1
>>> include Configuration: 1s_1/2^2 3p_1/2^2 3p_3/2^0
(Re-) Define a new standard subshell list.
> Start SCF process with hydrogenic orbitals.
>> (Re-) Define a storage array for dealing with single-electron TTp B-spline matrices:
-----------------------------------------------------------------------------
Index Subshell Energies [a.u.] Dirac-E [a.u.] Delta-E / |E|
-----------------------------------------------------------------------------
1 1s_1/2 -1.80086331e+01 -1.80086350e+01 +1.03994916e-07
2 2s_1/2 -4.50269964e+00 -4.50269857e+00 -2.38959499e-07
3 3s_1/2 -2.00098006e+00 -2.00095940e+00 -1.03240049e-05
4 4s_1/2 -1.12562939e+00 -1.12543844e+00 -1.69632831e-04
5 5s_1/2 -7.21116908e-01 -7.20234830e-01 -1.22321194e-03
6 6s_1/2 -5.29936446e-01 -5.00139888e-01 -5.62266634e-02
7 7s_1/2 -4.90049428e-01 -3.67436826e-01 -2.50204560e-01
: :
61 61s_1/2 +1.07252829e+09 -4.83755953e-03 +1.00000000e+00
62 62s_1/2 +2.82418557e+09 -4.68276538e-03 +1.00000000e+00
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
Index Subshell Energies [a.u.] Dirac-E [a.u.] Delta-E / |E|
-----------------------------------------------------------------------------
1 2p_1/2 -4.50269888e+00 -4.50269857e+00 -6.91701387e-08
2 3p_1/2 -2.00096714e+00 -2.00095940e+00 -3.86903735e-06
3 4p_1/2 -1.12549933e+00 -1.12543844e+00 -5.41023530e-05
4 5p_1/2 -7.20400487e-01 -7.20234830e-01 -2.29951315e-04
5 6p_1/2 -5.00841849e-01 -5.00139888e-01 -1.40156319e-03
6 7p_1/2 -3.91479761e-01 -3.67436826e-01 -6.14155231e-02
7 8p_1/2 -3.49594545e-01 -2.81311119e-01 -1.95321771e-01
: :
61 62p_1/2 +1.31159793e+09 -4.68276538e-03 +1.00000000e+00
62 63p_1/2 +6.75311123e+09 -4.53528382e-03 +1.00000000e+00
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
Index Subshell Energies [a.u.] Dirac-E [a.u.] Delta-E / |E|
-----------------------------------------------------------------------------
1 2p_3/2 -4.50053984e+00 -4.50053930e+00 -1.19777072e-07
2 3p_3/2 -2.00033044e+00 -2.00031959e+00 -5.42476571e-06
3 4p_3/2 -1.12535222e+00 -1.12516853e+00 -1.63226680e-04
4 5p_3/2 -7.21810564e-01 -7.20096642e-01 -2.37447501e-03
5 6p_3/2 -5.88521709e-01 -5.00059921e-01 -1.50311851e-01
6 7p_3/2 -4.94086009e-01 -3.67386470e-01 -2.56432153e-01
7 8p_3/2 -3.59118075e-01 -2.81277385e-01 -2.16755142e-01
: :
61 62p_3/2 +1.75386123e+09 -4.68269292e-03 +1.00000000e+00
62 63p_3/2 +6.75822628e+09 -4.53521476e-03 +1.00000000e+00
-----------------------------------------------------------------------------
>> (Re-) Define a storage array for dealing with single-electron TTp B-spline matrices:
Iteration 1 for symmetries ...
1s_1/2:: en [a.u.] = -1.0558212e+01; self-cons'cy = 2.6081e-01 [1.0000e+02 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.3423610e+00; self-cons'cy = 5.4069e-01 [1.0000e+02 for sym-block kappa = -1]
3s_1/2:: en [a.u.] = -4.1672621e-01; self-cons'cy = 6.5527e-01 [1.0000e+02 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.0764473e+00; self-cons'cy = 6.1412e-01 [1.0000e+02 for sym-block kappa = 1]
3p_1/2:: en [a.u.] = -3.5366557e-01; self-cons'cy = 6.9960e-01 [1.0000e+02 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.1545030e+00; self-cons'cy = 5.9169e-01 [1.0000e+02 for sym-block kappa = -2]
3p_3/2:: en [a.u.] = -3.8676450e-01; self-cons'cy = 6.7595e-01 [1.0000e+02 for sym-block kappa = -2]
Iteration 2 for symmetries ...
1s_1/2:: en [a.u.] = -1.1609809e+01; self-cons'cy = 4.7438e-02 [3.3975e-01 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.6448839e+00; self-cons'cy = 1.0127e-01 [3.3975e-01 for sym-block kappa = -1]
3s_1/2:: en [a.u.] = -5.1720554e-01; self-cons'cy = 1.0759e-01 [3.3975e-01 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.3383100e+00; self-cons'cy = 1.0844e-01 [2.0010e-01 for sym-block kappa = 1]
3p_1/2:: en [a.u.] = -4.4197623e-01; self-cons'cy = 1.1099e-01 [2.0010e-01 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.3252820e+00; self-cons'cy = 6.8868e-02 [1.7677e-01 for sym-block kappa = -2]
3p_3/2:: en [a.u.] = -4.3842354e-01; self-cons'cy = 6.2603e-02 [1.7677e-01 for sym-block kappa = -2]
Iteration 3 for symmetries ...
1s_1/2:: en [a.u.] = -1.1462046e+01; self-cons'cy = 6.4045e-03 [5.7217e-02 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.6025415e+00; self-cons'cy = 1.3039e-02 [5.7217e-02 for sym-block kappa = -1]
3s_1/2:: en [a.u.] = -5.0586302e-01; self-cons'cy = 1.1087e-02 [5.7217e-02 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.3061648e+00; self-cons'cy = 1.2156e-02 [2.8380e-02 for sym-block kappa = 1]
3p_1/2:: en [a.u.] = -4.3292972e-01; self-cons'cy = 1.0340e-02 [2.8380e-02 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.3066949e+00; self-cons'cy = 7.0621e-03 [2.6977e-02 for sym-block kappa = -2]
3p_3/2:: en [a.u.] = -4.3313053e-01; self-cons'cy = 6.0731e-03 [2.6977e-02 for sym-block kappa = -2]
Iteration 4 for symmetries ...
1s_1/2:: en [a.u.] = -1.1478600e+01; self-cons'cy = 7.2161e-04 [4.5128e-03 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.6064829e+00; self-cons'cy = 1.2282e-03 [4.5128e-03 for sym-block kappa = -1]
3s_1/2:: en [a.u.] = -5.0688854e-01; self-cons'cy = 1.0126e-03 [4.5128e-03 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.3087889e+00; self-cons'cy = 1.0035e-03 [2.0972e-03 for sym-block kappa = 1]
3p_1/2:: en [a.u.] = -4.3367610e-01; self-cons'cy = 8.6126e-04 [2.0972e-03 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.3080689e+00; self-cons'cy = 5.2547e-04 [1.7483e-03 for sym-block kappa = -2]
3p_3/2:: en [a.u.] = -4.3352376e-01; self-cons'cy = 4.5373e-04 [1.7483e-03 for sym-block kappa = -2]
Iteration 5 for symmetries ...
1s_1/2:: en [a.u.] = -1.1476785e+01; self-cons'cy = 7.9067e-05 [4.1574e-04 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.6061178e+00; self-cons'cy = 1.1364e-04 [4.1574e-04 for sym-block kappa = -1]
3s_1/2:: en [a.u.] = -5.0679662e-01; self-cons'cy = 9.0679e-05 [4.1574e-04 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.3085843e+00; self-cons'cy = 7.8167e-05 [1.7869e-04 for sym-block kappa = 1]
3p_1/2:: en [a.u.] = -4.3361701e-01; self-cons'cy = 6.8134e-05 [1.7869e-04 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.3079682e+00; self-cons'cy = 3.8486e-05 [1.4050e-04 for sym-block kappa = -2]
3p_3/2:: en [a.u.] = -4.3349401e-01; self-cons'cy = 3.4304e-05 [1.4050e-04 for sym-block kappa = -2]
Iteration 6 for symmetries ...
1s_1/2:: en [a.u.] = -1.1476978e+01; self-cons'cy = 8.3848e-06 [3.3502e-05 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.6061497e+00; self-cons'cy = 9.9264e-06 [3.3502e-05 for sym-block kappa = -1]
3s_1/2:: en [a.u.] = -5.0680430e-01; self-cons'cy = 7.5788e-06 [3.3502e-05 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.3085979e+00; self-cons'cy = 5.2040e-06 [1.2077e-05 for sym-block kappa = 1]
3p_1/2:: en [a.u.] = -4.3362104e-01; self-cons'cy = 4.6510e-06 [1.2077e-05 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.3079742e+00; self-cons'cy = 2.2956e-06 [8.6883e-06 for sym-block kappa = -2]
3p_3/2:: en [a.u.] = -4.3349589e-01; self-cons'cy = 2.1603e-06 [8.6883e-06 for sym-block kappa = -2]
Iteration 7 for symmetries ...
1s_1/2:: en [a.u.] = -1.1476958e+01; self-cons'cy = 8.7366e-07 [2.6315e-06 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.6061470e+00; self-cons'cy = 8.4117e-07 [2.6315e-06 for sym-block kappa = -1]
3s_1/2:: en [a.u.] = -5.0680369e-01; self-cons'cy = 6.0715e-07 [2.6315e-06 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.3085972e+00; self-cons'cy = 2.6300e-07 [7.0594e-07 for sym-block kappa = 1]
3p_1/2:: en [a.u.] = -4.3362082e-01; self-cons'cy = 2.5590e-07 [7.0594e-07 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.3079740e+00; self-cons'cy = 8.2561e-08 [4.1323e-07 for sym-block kappa = -2]
3p_3/2:: en [a.u.] = -4.3349580e-01; self-cons'cy = 9.8340e-08 [4.1323e-07 for sym-block kappa = -2]
Iteration 8 for symmetries ...
1s_1/2:: en [a.u.] = -1.1476960e+01; self-cons'cy = 8.9535e-08 [1.8383e-07 for sym-block kappa = -1]
2s_1/2:: en [a.u.] = -1.6061472e+00; self-cons'cy = 6.7721e-08 [1.8383e-07 for sym-block kappa = -1]
3s_1/2:: en [a.u.] = -5.0680373e-01; self-cons'cy = 4.4683e-08 [1.8383e-07 for sym-block kappa = -1]
2p_1/2:: en [a.u.] = -1.3085972e+00; self-cons'cy = 3.6359e-10 [1.8259e-08 for sym-block kappa = 1]
3p_1/2:: en [a.u.] = -4.3362082e-01; self-cons'cy = 4.2555e-09 [1.8259e-08 for sym-block kappa = 1]
2p_3/2:: en [a.u.] = -1.3079740e+00; self-cons'cy = 6.3222e-09 [6.3222e-09 for sym-block kappa = -2]
3p_3/2:: en [a.u.] = -4.3349580e-01; self-cons'cy = 2.0040e-09 [6.3222e-09 for sym-block kappa = -2]
> Compute CI matrix of dimension 2 x 2 for the symmetry 1^- ... ... done.
> Compute CI matrix of dimension 6 x 6 for the symmetry 0^+ ... ... done.
Eigenenergies:
Level J Parity Hartrees eV [eV]
1 0 + -3.648798247813691e+01 -9.928885763035797e+02 -9.928885763035797e+02
2 1 - -3.624410615291826e+01 -9.862523634769160e+02 -9.862523634769160e+02
3 1 - -3.596230472421729e+01 -9.785841560195463e+02 -9.785841560195463e+02
4 0 + -3.584833202046239e+01 -9.754828007819285e+02 -9.754828007819285e+02
5 0 + -3.557520093359942e+01 -9.680505253432382e+02 -9.680505253432382e+02
6 0 + -3.407634141250472e+01 -9.272644803249779e+02 -9.272644803249779e+02
7 0 + -3.390434882134372e+01 -9.225843235343160e+02 -9.225843235343160e+02
8 0 + -3.379320402360886e+01 -9.195599195390378e+02 -9.195599195390378e+02
Energy of each level relative to immediately lower level:
Level J Parity Hartrees eV [eV]
2 1 - 2.438763252186504e-01 6.636212826663757e+00 6.636212826663757e+00
3 1 - 2.818014287009731e-01 7.668207457369687e+00 7.668207457369687e+00
4 0 + 1.139727037549036e-01 3.101355237617787e+00 3.101355237617787e+00
5 0 + 2.731310868629677e-01 7.432275438690366e+00 7.432275438690366e+00
6 0 + 1.498859521094701e+00 4.078604501826025e+01 4.078604501826025e+01
7 0 + 1.719925911609934e-01 4.680156790661831e+00 4.680156790661831e+00
8 0 + 1.111447977348661e-01 3.024403995278262e+00 3.024403995278262e+00
Energy of each level relative to lowest level:
Level J Parity Hartrees eV [eV]
2 1 - 2.438763252186504e-01 6.636212826663757e+00 6.636212826663757e+00
3 1 - 5.256777539196236e-01 1.430442028403344e+01 1.430442028403344e+01
4 0 + 6.396504576745272e-01 1.740577552165123e+01 1.740577552165123e+01
5 0 + 9.127815445374949e-01 2.483805096034160e+01 2.483805096034160e+01
6 0 + 2.411641065632196e+00 6.562409597860184e+01 6.562409597860184e+01
7 0 + 2.583633656793189e+00 7.030425276926367e+01 7.030425276926367e+01
8 0 + 2.694778454528056e+00 7.332865676454193e+01 7.332865676454193e+01
Eigenenergies:
Level J Parity Hartrees eV [eV]
1 0 + -3.648798247813691e+01 -9.928885763035797e+02 -9.928885763035797e+02
2 1 - -3.624410615291826e+01 -9.862523634769160e+02 -9.862523634769160e+02
3 1 - -3.596230472421729e+01 -9.785841560195463e+02 -9.785841560195463e+02
4 0 + -3.584833202046239e+01 -9.754828007819285e+02 -9.754828007819285e+02
5 0 + -3.557520093359942e+01 -9.680505253432382e+02 -9.680505253432382e+02
6 0 + -3.407634141250472e+01 -9.272644803249779e+02 -9.272644803249779e+02
7 0 + -3.390434882134372e+01 -9.225843235343160e+02 -9.225843235343160e+02
8 0 + -3.379320402360886e+01 -9.195599195390378e+02 -9.195599195390378e+02
Energy of each level relative to lowest level:
Level J Parity Hartrees eV [eV]
2 1 - 2.438763252186504e-01 6.636212826663757e+00 6.636212826663757e+00
3 1 - 5.256777539196236e-01 1.430442028403344e+01 1.430442028403344e+01
4 0 + 6.396504576745272e-01 1.740577552165123e+01 1.740577552165123e+01
5 0 + 9.127815445374949e-01 2.483805096034160e+01 2.483805096034160e+01
6 0 + 2.411641065632196e+00 6.562409597860184e+01 6.562409597860184e+01
7 0 + 2.583633656793189e+00 7.030425276926367e+01 7.030425276926367e+01
8 0 + 2.694778454528056e+00 7.332865676454193e+01 7.332865676454193e+01
We finish this (simple) tutorial by enlarging the configuration basis for the low-lying levels but by restricting the CI computations to the level symmetries $J^P = 0^+$ and $1^-$. This is achieved by specifying the settings to: