"""Build your own shell model -- step 3: read the interaction.

The physics of the valence space lives in one small text file: usdb.int, the
USDB interaction (Brown & Richter 2006), fit to hundreds of measured sd-shell
levels. The format (oxbash convention) is:

  header line:   -63  2.1117  -3.9257  -3.2079   16.0  18.0  0.3
                  |     |________|________|        |     |    |
                  |     single-particle energies   |     |    |
                  |     for labels 1,2,3 (MeV):    A_core A_ref power
                  |     1=d3/2  2=d5/2  3=s1/2
                  '-- line count (negative flags the mass scaling below)

  then one line per two-body matrix element (TBME):
      a  b  c  d   J  T   V
  meaning <ab; JT| V |cd; JT> in MeV: the interaction energy between a PAIR
  of particles put in orbitals (a,b) and a pair in (c,d), with the pair
  coupled to total angular momentum J and isospin T. Isospin is the
  proton/neutron label treated like a spin-1/2: T=1 pairs are the
  symmetric combination (includes nn and pp), T=0 pairs antisymmetric
  (proton-neutron only).

The mass scaling: USDB was fit at A = 18, but orbitals shrink as more valence
particles pile in, so every TBME is scaled by (18/A)^0.3 for nucleus mass A.
For 18O the factor is exactly 1; for 24Mg it is 0.917315 -- and BIGSTICK's
input file carries the same "1.0 18.0 24.0 0.3" line, which is how the
convention was cross-checked.
"""

def read_int(path):
    spe, tbme, tail, n = {}, [], None, None
    with open(path) as fh:
        for raw in fh:
            line = raw.strip()
            if not line or line.startswith(('!', '#')):
                continue
            w = line.split()
            if n is None:
                n = int(w[0])
                spe = {1: float(w[1]), 2: float(w[2]), 3: float(w[3])}
                tail = [float(x) for x in w[4:]]
                continue
            a, b, c, d, J, T = (int(x) for x in w[:6])
            tbme.append(((a, b, c, d, J, T), float(w[6])))
    return spe, tbme, tail

NAMES = {1: 'd3/2', 2: 'd5/2', 3: 's1/2'}
spe, tbme, tail = read_int('usdb.int')

print("single-particle energies (the seat prices, MeV):")
for label in (2, 3, 1):   # printed in energy order
    print(f"  {NAMES[label]:>5}: {spe[label]:+8.4f}")
print("d5/2 lowest -- filling it first is the cheap move; the interaction")
print("decides whether that naive filling survives.")

print(f"\n{len(tbme)} two-body matrix elements. The first few:")
print(f"{'a':>5} {'b':>5} {'c':>5} {'d':>5} {'J':>2} {'T':>2} {'V (MeV)':>9}")
for (key, V) in tbme[:5]:
    a, b, c, d, J, T = key
    print(f"{NAMES[a]:>5} {NAMES[b]:>5} {NAMES[c]:>5} {NAMES[d]:>5} {J:>2} {T:>2} {V:>9.4f}")

a_core, a_ref, power = tail[-3], tail[-2], tail[-1]
for A in (18, 20, 24):
    scale = (a_ref / A) ** power
    print(f"mass scaling for A = {A}:  ({a_ref:.0f}/{A})^{power} = {scale:.6f}")
print("\nEvery number the final calculation produces is built from these")
print(f"{len(tbme)} matrix elements, 3 seat prices, and antisymmetry. Nothing else.")