Notebooks · lib/skyrme/spherical

Skyrme Hartree-Fock of any spherical closed-shell nucleus

The solver of lib/skyrme/skyrme for any nucleus whose protons and neutrons each fill whole $(n, l, j)$ shells: oxygen-16 through lead-208, with the forces SLy4, SLy5, SIII and SkM*.

from skyrme.spherical import SLY5, nucleus_z, nucleus_n, solve

res = solve(0.1, 200, SLY5, nucleus_z("208Pb"), nucleus_n("208Pb"), 0.5, 1.0e-9, 300)
show res.e_per_a       # E/A, MeV

Two things change from lib/skyrme/skyrme, which solves oxygen-16 only.

The occupied orbitals are chosen by energy. Each iteration diagonalizes every $(l, j)$ block up to $l = 7$ , and the protons (or neutrons) go into the lowest orbitals across all blocks, $2 j + 1$ to an orbital, until all $Z$ (or $N$ ) are placed. A nucleus whose last orbital would be only partly filled is open-shell and is refused: a spherical solution with a partly filled shell needs pairing or a filling approximation, which this solver does not have. Inside a block the $n$ -th orbital is the $n$ -th eigenvalue, so calcium-40's $2 s_{1/2}$ is the second state of the $(0, 1/2)$ block.

The force and the center-of-mass factor are arguments. A force is a Force record (SLY4, SLY5, SIII and SKMS are defined below), with a flag for the $J^{2}$ terms of $t_{1}$ and $t_{2}$ , which SLy5 and SIII keep and SLy4 and SkM* drop. The kinetic term carries $ℏ^{2} /2 m (1 - 1/ A)$ for every force, as each was fitted.

For oxygen-16 with SLy5 the arithmetic is lib/skyrme/skyrme's, operation for operation, and gives its bits (the test at the end).

Used with import skyrme.spherical.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

This program declares a record type (Force), and records run only in compiled C for now, so the Run buttons here are off. Open it in the workbench and use its C tab to compile it.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 36–692

from scf.loop import scf
from scf.mix import MixHistory, MixStep, broyden, linear, mix_densities, pulay
from skyrme.skyrme import E2, EIG_TOL, FOLD, HBARC, MIN_ITER, M_NUC, NEAR, NEUTRON, NKIND, PI, PROTON, WS_DIFFUSE, WS_R0, coulomb, count_below, grad1, integrate_r2, inverse_iteration, lap_r, particle_number, proton_coulomb, skyrme_o16

#| A Skyrme force: the ten parameters of the functional and whether its
#| $J^2$ terms are kept. t0 is in MeV fm^3, t1, t2 and w0 in MeV fm^5, t3 in
#| MeV fm^(3 + 3 alpha); the x's and alpha have no unit.
record Force {
    t0: Num
    t1: Num
    t2: Num
    t3: Num
    x0: Num
    x1: Num
    x2: Num
    x3: Num
    w0: Num
    alpha: Num
    j2: Bool
}

# SLy4 and SLy5: Chabanat, Bonche, Haensel, Meyer and Schaeffer, Nucl. Phys.
# A 635, 231 (1998), Table 1 (SLy5's numbers as skyrme_rpa's skyrme.in lists
# them, the same to the digits printed there)
SLY4 = Force { t0: -2488.91, t1: 486.82, t2: -546.39, t3: 13777.0, x0: 0.834, x1: -0.344, x2: -1.0, x3: 1.354, w0: 123.0, alpha: 1.0 / 6.0, j2: false }
SLY5 = Force { t0: -2484.88, t1: 483.13, t2: -549.40, t3: 13763.0, x0: 0.778, x1: -0.328, x2: -1.0, x3: 1.267, w0: 126.0, alpha: 1.0 / 6.0, j2: true }
# SIII: Beiner, Flocard, Van Giai and Quentin, Nucl. Phys. A 238, 29 (1975);
# its three-body t3 term is the density-dependent one with x3 = 1, alpha = 1
SIII = Force { t0: -1128.75, t1: 395.0, t2: -95.0, t3: 14000.0, x0: 0.45, x1: 0.0, x2: 0.0, x3: 1.0, w0: 120.0, alpha: 1.0, j2: true }
# SkM*: Bartel, Quentin, Brack, Guet and Hakansson, Nucl. Phys. A 386, 79 (1982)
SKMS = Force { t0: -2645.0, t1: 410.0, t2: -135.0, t3: 15595.0, x0: 0.09, x1: 0.0, x2: 0.0, x3: 0.0, w0: 130.0, alpha: 1.0 / 6.0, j2: false }

# ---- the blocks: l = 0 .. L_MAX, j = l + 1/2 and then (for l > 0) l - 1/2.
# L_MAX = 7 reaches 1j15/2, the first neutron shell above lead-208's 126.
L_MAX = 7
NBLOCK = 2 * L_MAX + 1
# at most this many orbitals of one block are occupied (lead-208's 3s1/2 is the third)
S_MAX = 4
# spectroscopic letters for l = 0 .. 7 (nuclear usage: j for l = 7)
L_LETTERS = "spdfghij"
# the kinds, for messages
KIND_NAMES = ["neutrons", "protons"]
# a state counts as bound below this energy, MeV
BOUND = 0.0
# the margin of the Sturm-count checks on an orbital energy, MeV
CHECK_TOL = 1.0e-7
# the element symbols, Z = 1 .. 118
ELEMENTS = split("H He Li Be B C N O F Ne Na Mg Al Si P S Cl Ar K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe Cs Ba La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Db Sg Bh Hs Mt Ds Rg Cn Nh Fl Mc Lv Ts Og", " ")

#| The l of block b; the blocks run s1/2, p3/2, p1/2, d5/2, d3/2, ...
def block_l(b) {
    if b == 0 {
        return 0
    }
    return (b + 1) // 2
}

#| The j of block b: l + 1/2 for odd b and b = 0, l - 1/2 for even b.
def block_j(b) {
    l = block_l(b)

from scf.loop import scf
from scf.mix import MixHistory, MixStep, broyden, linear, mix_densities, pulay
from skyrme.skyrme import E2, EIG_TOL, FOLD, HBARC, MIN_ITER, M_NUC, NEAR, NEUTRON, NKIND, PI, PROTON, WS_DIFFUSE, WS_R0, coulomb, count_below, grad1, integrate_r2, inverse_iteration, lap_r, particle_number, proton_coulomb, skyrme_o16

#| A Skyrme force: the ten parameters of the functional and whether its
#| $J^2$ terms are kept. t0 is in MeV fm^3, t1, t2 and w0 in MeV fm^5, t3 in
#| MeV fm^(3 + 3 alpha); the x's and alpha have no unit.
record Force {
    t0: Num
    t1: Num
    t2: Num
    t3: Num
    x0: Num
    x1: Num
    x2: Num
    x3: Num
    w0: Num
    alpha: Num
    j2: Bool
}

# SLy4 and SLy5: Chabanat, Bonche, Haensel, Meyer and Schaeffer, Nucl. Phys.
# A 635, 231 (1998), Table 1 (SLy5's numbers as skyrme_rpa's skyrme.in lists
# them, the same to the digits printed there)
SLY4 = Force { t0: -2488.91, t1: 486.82, t2: -546.39, t3: 13777.0, x0: 0.834, x1: -0.344, x2: -1.0, x3: 1.354, w0: 123.0, alpha: 1.0 / 6.0, j2: false }
SLY5 = Force { t0: -2484.88, t1: 483.13, t2: -549.40, t3: 13763.0, x0: 0.778, x1: -0.328, x2: -1.0, x3: 1.267, w0: 126.0, alpha: 1.0 / 6.0, j2: true }
# SIII: Beiner, Flocard, Van Giai and Quentin, Nucl. Phys. A 238, 29 (1975);
# its three-body t3 term is the density-dependent one with x3 = 1, alpha = 1
SIII = Force { t0: -1128.75, t1: 395.0, t2: -95.0, t3: 14000.0, x0: 0.45, x1: 0.0, x2: 0.0, x3: 1.0, w0: 120.0, alpha: 1.0, j2: true }
# SkM*: Bartel, Quentin, Brack, Guet and Hakansson, Nucl. Phys. A 386, 79 (1982)
SKMS = Force { t0: -2645.0, t1: 410.0, t2: -135.0, t3: 15595.0, x0: 0.09, x1: 0.0, x2: 0.0, x3: 0.0, w0: 130.0, alpha: 1.0 / 6.0, j2: false }

# ---- the blocks: l = 0 .. L_MAX, j = l + 1/2 and then (for l > 0) l - 1/2.
# L_MAX = 7 reaches 1j15/2, the first neutron shell above lead-208's 126.
L_MAX = 7
NBLOCK = 2 * L_MAX + 1
# at most this many orbitals of one block are occupied (lead-208's 3s1/2 is the third)
S_MAX = 4
# spectroscopic letters for l = 0 .. 7 (nuclear usage: j for l = 7)
L_LETTERS = "spdfghij"
# the kinds, for messages
KIND_NAMES = ["neutrons", "protons"]
# a state counts as bound below this energy, MeV
BOUND = 0.0
# the margin of the Sturm-count checks on an orbital energy, MeV
CHECK_TOL = 1.0e-7
# the element symbols, Z = 1 .. 118
ELEMENTS = split("H He Li Be B C N O F Ne Na Mg Al Si P S Cl Ar K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe Cs Ba La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Db Sg Bh Hs Mt Ds Rg Cn Nh Fl Mc Lv Ts Og", " ")

#| The l of block b; the blocks run s1/2, p3/2, p1/2, d5/2, d3/2, ...
def block_l(b) {
    if b == 0 {
        return 0
    }
    return (b + 1) // 2
}

#| The j of block b: l + 1/2 for odd b and b = 0, l - 1/2 for even b.
def block_j(b) {
    l = block_l(b)