Step 2 — The basis as bitstrings
Nuclear Physics
Step 1's determinants were "seat 3 and seat 7 occupied." Now give each one a name a computer loves: a bitstring — a single integer in which bit i set means seat i occupied. Two neutrons in twelve seats is a 12-bit number with exactly two bits set, and the whole ¹⁸O basis is fourteen small integers.
The payoff is that the three operations a shell-model code performs billions of times are each one machine instruction on this representation: is seat i occupied (shift and mask), add or remove a particle (set or clear a bit), and count the occupied seats below position i — a popcount, which looks like bookkeeping now but becomes physics in step 4, where that count sets the ± sign of a matrix element. Production codes are built on exactly this idea: CENS stores every determinant in one 64-bit integer, OpenFCI packs the same scheme into 512-bit registers, and this file is the idea at 12 bits, where you can still read every state by eye.
The program
"""Build your own shell model -- step 2: the basis as bitstrings.
Step 1 listed 12 seats and counted 14 determinants. Now give each determinant
a NAME a computer loves: one integer, one bit per seat, bit i set = seat i
occupied. Two neutrons in twelve seats is a 12-bit number with exactly two
bits set.
Why bits and not lists? Three operations dominate a shell-model code, and all
three are single machine instructions on a bitstring:
occupied? det >> i & 1
create/destroy det | (1 << i), det & ~(1 << i)
count below i popcount(det & ((1 << i) - 1)) <- the fermionic sign,
coming in step 4
(Production codes live this way: CENS stores each determinant in an
`unsigned long long`; BIGSTICK and OpenFCI pack the same idea into bigger
registers. This file is that idea at 12 bits.)
The enumeration below is deliberately dumb -- try all pairs, keep Jz = 0.
For two particles that is instant. Step 6 will need the smarter version
(prune while enumerating, an idea borrowed from cosmo), because 24Mg has
28,503 determinants hiding in a space of C(24,8) = 735,471 candidates.
"""
from itertools import combinations
ORBITS = {1: ('d3/2', 2, 3), 2: ('d5/2', 2, 5), 3: ('s1/2', 0, 1)}
SEATS = []
for label, (name, l, j2) in ORBITS.items():
for m2 in range(-j2, j2 + 1, 2):
SEATS.append((name, j2, m2))
def seat_tag(i):
name, j2, m2 = SEATS[i]
return f"{name}[m={m2:+d}/2]"
basis = []
for (i, j) in combinations(range(12), 2):
if SEATS[i][2] + SEATS[j][2] == 0:
basis.append((1 << i) | (1 << j))
basis.sort()
print("the 14 determinants of 18O (2 neutrons, Jz = 0), as 12-bit integers:")
print(f"{'k':>2} {'integer':>7} {'bits (seat 11 ... seat 0)':>26} occupied seats")
for k, det in enumerate(basis):
occ = [i for i in range(12) if det >> i & 1]
tags = ", ".join(seat_tag(i) for i in occ)
print(f"{k:>2} {det:>7} {det:012b} {tags}")
print(f"\ndimension: {len(basis)}")
print("Each row IS a many-body basis state. The Hamiltonian will be a")
print("14 x 14 matrix connecting these integers to each other.")000000000110 (both d₃/₂ partners near m = 0) to
110000000000 (the s₁/₂ pair) — each with its occupied seats
named. Read one row and you have read a many-body quantum state off an
integer.