`vp.nuclear.skyrme_hf`

The Skyrme-Hartree-Fock ground state of a spherical closed-shell nucleus: binding energy, radii, single-particle levels and densities, with the forces SLy4, SLy5, SIII or SkM*.

Part of the distribution vipyr-nuclear; it runs from a checkout of the repository. Planned on PyPI as pip install "vipyr[nuclear]".

Give it a nucleus whose protons and neutrons both fill whole shells, such as "16O", "40Ca", "48Ca", "56Ni", "90Zr", "132Sn" or "208Pb", and a Skyrme force. It returns the self-consistent mean-field ground state: the total energy and the energy per nucleon, the proton, neutron and charge radii, the energy of every occupied orbital with its name ("n 3p1/2"), the other bound orbitals of the final mean field, and the densities on the radial mesh.

This is the textbook calculation of nuclear structure (Vautherin and Brink, 1972): each nucleon moves in a field made by all the others, and the equations are solved until the field and the orbitals agree. Lead-208 takes about 0.09 s on 186 mesh points.

It does not do open-shell nuclei. A nucleus whose last orbital would be partly filled, such as calcium-44, is refused with the name of that orbital: treating it needs pairing (BCS or Hartree-Fock-Bogoliubov), which this solver does not have.

Signature

vp.nuclear.skyrme_hf(*, nucleus="16O", force="sly5", t0=preset.t0, t1=preset.t1, t2=preset.t2, t3=preset.t3, x0=preset.x0, x1=preset.x1, x2=preset.x2, x3=preset.x3, w0=preset.w0, alpha=preset.alpha, j2=preset.j2, spacing=0.1, box=..., mixing="pulay", mix=0.5, history=6, e_tol=1.0e-9, max_iter=300)

Parameters

nucleus : str, default "16O": The nucleus, as its mass number and element symbol ("208Pb"; "Pb208" and "pb-208" also work).
force : str, default "sly5": The Skyrme force; t0 ... j2 start from its values. Range: one of sly4, sly5, siii, skm*. Case does not matter.
t0 : float, default preset.t0: The force's t0, the contact term.
t1 : float, default preset.t1: The force's t1, the momentum-dependent term even in the relative momentum.
t2 : float, default preset.t2: The force's t2, the term odd in the relative momentum.
t3 : float, default preset.t3: The force's t3, the density-dependent term.
x0 : float, default preset.x0: The spin exchange of the t0 term.
x1 : float, default preset.x1: The spin exchange of the t1 term.
x2 : float, default preset.x2: The spin exchange of the t2 term.
x3 : float, default preset.x3: The spin exchange of the t3 term.
w0 : float, default preset.w0: The spin-orbit strength.
alpha : float, default preset.alpha: The power of the density in the t3 term. Range: >= 0, <= 1.
j2 : bool, default preset.j2: Whether the J^2 terms of t1 and t2 are kept (SLy5 and SIII keep them, SLy4 and SkM* drop them).
spacing : float, default 0.1: The mesh spacing. Range: > 0, <= 1.
box : float, default max(15.0, 1.2 * pow(Z_NUC + N_NUC, 1.0 / 3.0) + 11.5): The width of the box; the orbitals are zero past it. Range: >= 5, <= 60.
mixing : str, default "pulay": How each pass's densities are mixed into the next pass's: linear keeps MIX of the new ones; pulay and broyden combine the last HISTORY passes (lib/scf/mix.viv). Range: one of linear, pulay, broyden, every, every simple, every accelerated. Case does not matter.
mix : float, default 0.5: The fraction of each pass's new densities mixed into the old ones; for pulay and broyden, the fraction of their combined residual. Range: > 0, <= 1.
history : int, default 6: How many passes pulay and broyden combine, the current one included (pulay and broyden only; linear ignores it). Range: >= 1, <= 50.
e_tol : float, default 1.0e-9: The iteration stops when the total energy changes by less than this between passes. Range: >= 0.
max_iter : int, default 300: At most this many passes; a run that does not settle returns converged = False. Range: >= 1, <= 10000.

Returns

A frozen SkyrmeResult. Iterating it gives the fields below in order; it also has params, tests, log and check().

z : int: The protons.
n : int: The neutrons.
e_per_a : float, MeV: The energy per nucleon, MeV.
energy : float, MeV: The total energy, MeV.
r_p : float, fm: The protons' root mean square radius, fm.
r_n : float, fm: The neutrons' root mean square radius, fm.
r_ch : float, fm: The charge radius: the proton density folded with a Gaussian of range 0.65 fm, fm.
levels : ndarray[float64], MeV: The occupied orbitals' energies, MeV: the neutrons' lowest first, then the protons'.
labels : str: The names of LEVELS, comma-separated: "n 1s1/2,n 1p3/2,...,p 1s1/2,...".
spectrum : ndarray[float64], MeV: Every bound orbital of the converged fields, MeV: the neutrons' lowest first, then the protons'.
spectrum_labels : str: The names of SPECTRUM, comma-separated.
spectrum_occ : ndarray[float64]: The nucleons in each SPECTRUM orbital: 2j + 1 if occupied, else 0.
iterations : int: The passes run.
converged : bool: Whether the energy settled to E_TOL within MAX_ITER passes.
refills : int: The passes after the first in which the occupation was chosen again, because an orbital changed places.
mesh : ndarray[float64], fm: The mesh points, fm.
rho : ndarray[float64], fm^-3: The number densities, fm^-3, [2, round(BOX / SPACING)]: row 0 neutrons, row 1 protons.
tau : ndarray[float64], fm^-5: The kinetic-energy densities, fm^-5, [2, round(BOX / SPACING)].
jso : ndarray[float64], fm^-4: The spin-orbit densities, fm^-4, [2, round(BOX / SPACING)].

Raises

vipyr.ParameterError: A ValueError: an argument outside its range, choices, rank or shape, or one of a group given without the rest; its .param names the argument.
TypeError: An argument of the wrong type, or an unknown keyword.
vipyr.VivBuildError: An ImportError: no C compiler, or the program's C does not build.
RuntimeError: The program stops with its own message (a combination it rejects, or no convergence where the program treats that as an error).

Notes

The energy functional

The energy is the integral of a local density $E = \int H (r) d^{3} r$ built from three densities of each kind of nucleon $q \in {n, p}$ : the number density $ρ_{q}$ , the kinetic-energy density $τ_{q}$ and the spin-orbit density $J_{q}$ . Without a subscript, $ρ = ρ_{n} + ρ_{p}$ and likewise for $τ$ and $J$ :

H = \frac{ℏ ^{2}}{2 m} (1 - \frac{1}{A}) τ + \frac{t _{0}}{2} [(1 + \frac{x _{0}}{2}) ρ^{2} - (x_{0} + \frac{1}{2}) q \sum ρ_{q}^{2}] + \frac{t _{1}}{4} [(1 + \frac{x _{1}}{2}) (ρ τ + \frac{3}{4} (\nabla ρ)^{2}) - (x_{1} + \frac{1}{2}) q \sum (ρ_{q} τ_{q} + \frac{3}{4} (\nabla ρ_{q})^{2})] + \frac{t _{2}}{4} [(1 + \frac{x _{2}}{2}) (ρ τ - \frac{1}{4} (\nabla ρ)^{2}) + (x_{2} + \frac{1}{2}) q \sum (ρ_{q} τ_{q} - \frac{1}{4} (\nabla ρ_{q})^{2})] + \frac{t _{3}}{12} ρ^{α} [(1 + \frac{x _{3}}{2}) ρ^{2} - (x_{3} + \frac{1}{2}) q \sum ρ_{q}^{2}] + \frac{W _{0}}{2} (J \nabla ρ + q \sum J_{q} \nabla ρ_{q}) + \frac{t _{1} - t _{2}}{16} q \sum J_{q}^{2} - \frac{t _{1} x _{1} + t _{2} x _{2}}{16} J^{2} + \frac{e ^{2}}{2} ρ_{p} (r) \int \frac{ρ _{p} ( r ^{'} )}{∣ r - r ^{'} ∣} d^{3} r^{'} - \frac{3}{4} e^{2} (\frac{3}{π})^{1/3} ρ_{p}^{4/3} .

The two $J^{2}$ terms are kept for SLy5 and SIII and dropped for SLy4 and SkM*, as each force was fitted (the j2 parameter). The last two terms are the direct Coulomb energy and its exchange in the Slater approximation. The factor $1 - 1/ A$ is the one-body center-of-mass correction, which all four forces use; it multiplies $ℏ^{2} /2 m$ in the energy and in the equations below.

The Hartree-Fock equations

For a spherical nucleus each orbital is $R_{n l j} (r)$ times a spinor spherical harmonic, and making $E$ stationary under changes of the orbitals gives, for $u = r R$ ,

- (B_{q} u^{'})^{'} + [B_{q} \frac{l ( l + 1 )}{r ^{2}} + \frac{B _{q}^{'}}{r} + U_{q} + δ_{q p} V_{C} + ⟨ σ \cdot l ⟩ \frac{W _{q}}{r}] u = ε u, ⟨ σ \cdot l ⟩ = j (j + 1) - l (l + 1) - \frac{3}{4},

with the effective mass $B_{q} = ℏ^{2} /2 m_{q}^{*}$ , the central field $U_{q}$ , the spin-orbit field $W_{q}$ and the Coulomb field $V_{C}$ :

B_{q} = U_{q} = W_{q} = V_{C} = \frac{ℏ ^{2}}{2 m} (1 - \frac{1}{A}) + \frac{1}{8} [t_{1} (2 + x_{1}) + t_{2} (2 + x_{2})] ρ + \frac{1}{8} [t_{2} (1 + 2 x_{2}) - t_{1} (1 + 2 x_{1})] ρ_{q}, t_{0} [(1 + \frac{x _{0}}{2}) ρ - (x_{0} + \frac{1}{2}) ρ_{q}] + \frac{t _{1}}{4} [(1 + \frac{x _{1}}{2}) (τ - \frac{3}{2} Δ ρ) - (x_{1} + \frac{1}{2}) (τ_{q} - \frac{3}{2} Δ ρ_{q})] + \frac{t _{2}}{4} [(1 + \frac{x _{2}}{2}) (τ + \frac{1}{2} Δ ρ) + (x_{2} + \frac{1}{2}) (τ_{q} + \frac{1}{2} Δ ρ_{q})] + \frac{t _{3}}{12} [(1 + \frac{x _{3}}{2}) (2 + α) ρ^{α + 1} - (x_{3} + \frac{1}{2}) (2 ρ^{α} ρ_{q} + α ρ^{α - 1} (ρ_{n}^{2} + ρ_{p}^{2}))] - \frac{W _{0}}{2} \nabla \cdot (J + J_{q}), \frac{W _{0}}{2} (ρ^{'} + ρ_{q}^{'}) + \frac{t _{1} - t _{2}}{8} J_{q} - \frac{t _{1} x _{1} + t _{2} x _{2}}{8} J (the J terms only with j 2), e^{2} \int \frac{ρ _{p} ( r ^{'} )}{∣ r - r ^{'} ∣} d^{3} r^{'} - e^{2} (\frac{3}{π})^{1/3} ρ_{p}^{1/3} .

The occupied orbitals make the densities, $ρ_{q} = \sum \frac{2 j + 1}{4 π} R^{2}$ , $τ_{q} = \sum \frac{2 j + 1}{4 π} [R^{'2} + \frac{l ( l + 1 )}{r ^{2}} R^{2}]$ and $J_{q} = \sum \frac{2 j + 1}{4 π} ⟨ σ \cdot l ⟩ \frac{R ^{2}}{r}$ , the sums running over the occupied $(n, l, j)$ of kind $q$ .

Discretization

The mesh is $r_{i} = (i + 1) h$ for $i < round (box / h)$ , with $u = 0$ at the origin and one step past the last point. The kinetic term is written in conservative form with $B_{q}$ at the midpoints (past the last point it is the free value), so each $(l, j)$ block is a symmetric tridiagonal matrix. Derivatives of the densities and fields are centered differences inside and second-order one-sided differences at the ends; integrals are the trapezoid rule; the direct Coulomb potential sums the charge inside and outside each mesh point.

Everything is second order in $h$ . Halving the spacing from 0.1 fm to 0.05 fm raises E/A by 0.014 to 0.019 MeV across the seven nuclei and four forces. For lead-208 with SLy4 the changes from 0.1 to 0.05 fm and from 0.05 to 0.025 fm are in the ratio 4.005, and the Richardson value $(4 E (h /2) - E (h)) /3$ from either pair agrees to $7 \times 1 0^{- 6}$ MeV per nucleon. Run two spacings and extrapolate when the third decimal of E/A matters.

The default box is 15 fm, or $1.2 A^{1/3} + 11.5$ fm when that is wider (18.6 fm for lead-208). Widening it to 30 fm changes E/A by less than $2 \times 1 0^{- 7}$ MeV for calcium-48, tin-132 and lead-208. Orbitals within about 1 MeV of zero in the spectrum feel the box more than that.

Solving

The iteration starts from Woods-Saxon densities (radius $1.2 A^{1/3}$ fm, diffuseness 0.6 fm) and $τ = J = 0$ . Each pass computes the fields from the densities, builds the matrix of every block with $l \leq 7$ , finds the occupied orbitals, builds new densities and mixes them with the old ones (Pulay's mixing by default, below). It stops when the energy changes by less than e_tol between passes (after at least four). With the defaults the seven nuclei above take 22 to 32 passes with each of the four forces (37 to 66 with mixing="linear").

The loop is lib/scf/loop.viv, with this mesh as its basis, and the mixing is the contract mix_densities of lib/scf/mix.viv. mixing="pulay" (Pulay's mixing, the default) and mixing="broyden" (Johnson's modified Broyden) combine the last history passes; history applies to those two only. mixing="linear" keeps a fraction mix of the new densities. For calcium-40 with SLy4 at e_tol=1e-10 linear mixing takes 40 passes, Pulay 32 and Broyden 24 (gallery/science/plot_scf_mixers.py). Stopped there, the three differ in $r_{c h}$ by up to $2 \times 1 0^{- 9}$ fm, because the energy settles before the densities do; at e_tol=1e-13 they agree to $2 \times 1 0^{- 13}$ fm and $2 \times 1 0^{- 14}$ MeV in E/A. At the defaults (e_tol=1e-9), oxygen-16 with SLy5 by Pulay's mixing is within $7 \times 1 0^{- 11}$ MeV in E/A and $1 \times 1 0^{- 12}$ fm in $r_{c h}$ of the run to e_tol=1e-13; by linear mixing, $5 \times 1 0^{- 11}$ MeV and $1 \times 1 0^{- 9}$ fm.

The occupation is chosen by energy. In pass 1 every bound state of every block is found by Sturm-sequence bisection, and the $Z$ protons (and $N$ neutrons) go into the lowest orbitals, $2 j + 1$ to each. Later passes refine each occupied orbital by inverse iteration from its previous vector, then use Sturm counts to confirm that it is still the same state of its block and that no empty orbital has dropped below the highest occupied one. When either check fails the occupation is chosen again; refills counts those passes. Calcium-40's $2 s_{1/2}$ , the second orbital of its block, is found the same way as the first.

Limits

Spherical symmetry and closed shells only: no pairing and no deformation. A nucleus whose last proton or neutron orbital would be partly filled is refused with the name of that orbital. A closed $j$ subshell counts as closed: carbon-12 (six of each in $1 s_{1/2}$ and $1 p_{3/2}$ ) runs and converges, as a spherical nucleus.
Bound orbitals only: a nucleus whose occupied orbitals reach zero energy (past a drip line) is refused.
Blocks up to $l = 7$ and at most four occupied orbitals per block, which covers lead-208 (its $3 s_{1/2}$ protons are the third of their block). spectrum lists bound orbitals within those bounds.
The charge radius folds the proton density with a Gaussian of range 0.65 fm, $r_{c h}^{2} = r_{p}^{2} + \frac{3}{2} (0.65)^{2}$ , as skyrme_rpa does; it has no neutron or spin-orbit contributions to the charge density.
The constants are skyrme_rpa's: $ℏ c = 197.327053$ MeV fm, the nucleon mass the average of 938.27231 and 939.56563 MeV, and $e^{2} = 1.43986$ MeV fm.

Verification

Against an independent numpy solver (the reference below, its constants set to these) for all 28 pairs of nucleus and force on meshes of 0.1 and 0.05 fm: E/A agrees to $4 \times 1 0^{- 14}$ MeV, the radii to $9 \times 1 0^{- 13}$ fm and the occupied levels to $2 \times 1 0^{- 11}$ MeV, in the same number of passes, both mixing linearly. That solver diagonalizes each block densely and fills shells in a prescribed order, so the agreement also checks the filling by energy.
Oxygen-16 with SLy5 and mixing="linear" gives the bits of the oxygen-16 solver in lib/skyrme/skyrme.viv. Extrapolated to $h \to 0$ it gives E/A = −8.0250 MeV and $r_{c h}$ = 2.7993 fm, against −8.0261 MeV and 2.7991 fm in the sample output of skyrme_rpa (Colò et al. 2013).
The exported library, a native build and Viv's interpreter give the same result bits for every case in skyrmehf.cases.json.

References

D. Vautherin and D. M. Brink, Hartree-Fock calculations with Skyrme's interaction. I. Spherical nuclei, Phys. Rev. C 5, 626 (1972). doi:10.1103/PhysRevC.5.626
E. Chabanat, P. Bonche, P. Haensel, J. Meyer and R. Schaeffer, A Skyrme parametrization from subnuclear to neutron star densities. Part II. Nuclei far from stabilities, Nucl. Phys. A 635, 231 (1998); erratum Nucl. Phys. A 643, 441 (1998). SLy4 and SLy5. doi:10.1016/S0375-9474(98)00180-8
M. Beiner, H. Flocard, Nguyen Van Giai and P. Quentin, Nuclear ground-state properties and self-consistent calculations with the Skyrme interaction. (I). Spherical description, Nucl. Phys. A 238, 29 (1975). SIII. doi:10.1016/0375-9474(75)90338-3
J. Bartel, P. Quentin, M. Brack, C. Guet and H.-B. Håkansson, Towards a better parametrisation of Skyrme-like effective forces: a critical study of the SkM force, Nucl. Phys. A 386, 79 (1982). SkM*. doi:10.1016/0375-9474(82)90403-1
G. Colò, L. Cao, N. Van Giai and L. Capelli, Self-consistent RPA calculations with Skyrme-type interactions: The skyrme_rpa program, Comput. Phys. Commun. 184, 142 (2013). The constants and the charge folding used here, and the oxygen-16 check. doi:10.1016/j.cpc.2012.07.016

Examples

import vipyr as vp
r = vp.nuclear.skyrme_hf()             # oxygen-16 with SLy5
print(r.e_per_a, r.r_ch, r.levels)
assert r.check()

-8.046538617871763 2.7966837276110867 [-36.18659981 -20.63141616 -14.42364525 -32.38991873 -17.15296823
 -11.074404  ]

Calcium-40 with SLy4: the energy per nucleon, the charge radius and the occupied orbitals. Its $2 s_{1/2}$ orbitals are the second state of the $(l, j) = (0, 1/2)$ block.

import vipyr as vp

r = vp.nuclear.skyrme_hf(nucleus="40Ca", force="SLy4")
print(f"E/A = {r.e_per_a:.4f} MeV, r_ch = {r.r_ch:.4f} fm, {r.iterations} passes")
for label, e in zip(r.labels.split(","), r.levels):
    print(f"{label:9s} {e:8.3f}")
# E/A = -8.6293 MeV, r_ch = 3.5087 fm, 27 passes
# n 1s1/2    -48.337
# ...
# n 2s1/2    -17.272
# n 1d3/2    -15.313
# p 1s1/2    -40.323
# ...
# p 1d3/2     -8.359

E/A = -8.6293 MeV, r_ch = 3.5087 fm, 27 passes
n 1s1/2    -48.337
n 1p3/2    -35.048
n 1p1/2    -31.127
n 1d5/2    -22.035
n 2s1/2    -17.272
n 1d3/2    -15.313
p 1s1/2    -40.323
p 1p3/2    -27.518
p 1p1/2    -23.683
p 1d5/2    -14.903
p 2s1/2    -10.172
p 1d3/2     -8.359

Lead-208 with the four forces, extrapolated to zero mesh spacing. The error of each run goes as $h^{2}$ , so $(4 E (h /2) - E (h)) /3$ removes it. The measured values are E/A = −7.867 MeV and $r_{c h}$ = 5.501 fm.

import vipyr as vp

for force in ["SLy4", "SLy5", "SIII", "SkM*"]:
    a = vp.nuclear.skyrme_hf(nucleus="208Pb", force=force, spacing=0.1)
    b = vp.nuclear.skyrme_hf(nucleus="208Pb", force=force, spacing=0.05)
    e = (4 * b.e_per_a - a.e_per_a) / 3
    rch = (4 * b.r_ch - a.r_ch) / 3
    print(f"{force:5s} E/A {e:8.4f} MeV   r_ch {rch:.4f} fm   r_n - r_p {b.r_n - b.r_p:.3f} fm")
# SLy4  E/A  -7.8640 MeV   r_ch 5.5154 fm   r_n - r_p 0.160 fm
# SLy5  E/A  -7.8662 MeV   r_ch 5.5073 fm   r_n - r_p 0.162 fm
# SIII  E/A  -7.8395 MeV   r_ch 5.5813 fm   r_n - r_p 0.125 fm
# SkM*  E/A  -7.8667 MeV   r_ch 5.5093 fm   r_n - r_p 0.169 fm

SLy4  E/A  -7.8640 MeV   r_ch 5.5154 fm   r_n - r_p 0.160 fm
SLy5  E/A  -7.8662 MeV   r_ch 5.5073 fm   r_n - r_p 0.162 fm
SIII  E/A  -7.8395 MeV   r_ch 5.5813 fm   r_n - r_p 0.125 fm
SkM*  E/A  -7.8667 MeV   r_ch 5.5093 fm   r_n - r_p 0.169 fm

The neutron orbitals of lead-208 around the gap at N = 126, occupied and empty, from spectrum.

import vipyr as vp

r = vp.nuclear.skyrme_hf(nucleus="208Pb", force="SkM*")
for label, e, occ in zip(r.spectrum_labels.split(","), r.spectrum, r.spectrum_occ):
    if label.startswith("n") and -12.0 < e < -2.0:
        print(f"{label[2:]:7s} {e:8.3f} MeV  {'occupied' if occ else 'empty'}")
# 2f7/2    -11.913 MeV  occupied
# 1h9/2    -11.414 MeV  occupied
# 1i13/2   -10.233 MeV  occupied
# 3p3/2     -9.195 MeV  occupied
# 2f5/2     -8.892 MeV  occupied
# 3p1/2     -8.082 MeV  occupied
# 2g9/2     -4.186 MeV  empty
# 1j15/2    -2.332 MeV  empty

2f7/2    -11.913 MeV  occupied
1h9/2    -11.414 MeV  occupied
1i13/2   -10.233 MeV  occupied
3p3/2     -9.195 MeV  occupied
2f5/2     -8.892 MeV  occupied
3p1/2     -8.082 MeV  occupied
2g9/2     -4.186 MeV  empty
1j15/2    -2.332 MeV  empty

Each example ran as written when these pages were generated, and the output shown is what it printed.

Checked against

Reference: lib/skyrme/refs/skyrme_hf.py (the personal site's numpy solver, constants set to skyrme_rpa's) for every nucleus and force on two meshes (lib/skyrme/refs/compare.py, lib/skyrme/COMPARISON.md); the post's numpy twin of oxygen-16 and skyrme_rpa (Colò, Cao, Van Giai and Capelli, CPC 184, 142 (2013)) for oxygen-16 with SLy5 (tools/export/test_skyrmehf.py)
Cases: oxygen-16 with SLy5 at the defaults, and six more nuclei with the four forces at E_TOL 10⁻¹⁰ with linear mixing, as the numpy solver mixes (skyrmehf.cases.json)
Tolerance: 10⁻⁹ against the numpy solver at the same mesh (observed 2 × 10⁻¹¹ or less: the two eigensolvers' round-off)

Source

Manifest: lib/skyrme/skyrmehf.export.toml
Viv: lib/skyrme/skyrmehf.viv
The defs it reaches and the calls between them: map of the code
Watch it run, line by line: /explore/walk/skyrme-o16

vp.nuclear.skyrme_hf