“Know how to solve every problem that has been solved.” “What I cannot create, I do not understand.” — Richard Feynman
python 50 lines · 2.1 KB
"""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.")