Use NUBASE2020 for half lives

[pkienzle]

10.1088/1674-1137/abddae

[pkienzle]

Tables available here: https://www-nds.iaea.org/amdc/

[pkienzle]

Some parsing code for the nubase_4.mas20 file in the link above.

Can't use it until isomers are supported.

Note: the triple (AAA, ZZZ, isomer) is unique in the dataset.

from numpy import inf, nan

# Parse the nubase_4.mas20 file
NUBASE_FILE = 'nubase_4.mas20'
NUBASE_FILE = 'periodictable/nubase_4.mas20.txt'
YEAR = 31556926 # seconds = 365.2422 days
UNIT_CONVERSION = {
    'y' : YEAR,
    'd' : 24*3600,
    'h' : 3600,
    'm' : 60,
    's' : 1,
    'ms' : 1e-3,  'ky' : 1e3 * YEAR,
    'us' : 1e-6,  'My' : 1e6 * YEAR,
    'ns' : 1e-9,  'Gy' : 1e9 * YEAR,
    'ps' : 1e-12, 'Ty' : 1e12 * YEAR,
    'fs' : 1e-15, 'Py' : 1e15 * YEAR,
    'as' : 1e-18, 'Ey' : 1e18 * YEAR,
    'zs' : 1e-21, 'Zy' : 1e21 * YEAR,
    'ys' : 1e-24, 'Yy' : 1e24 * YEAR,
}

def nubase_float(s):
    # Note contains codes regarding the value: >, <, ~, >>
    # If the value is estimated, then note will end with #.
    # Currently not returning it because it isn't used.
    note = None
    s = s.strip()
    if not s:
        return None
    #if s[-1] == '#' and s[0] not in '-0123456789':
    #    print('parsing', s)
    tail = ''
    if s[-1] == '#':
        tail = s[-1]
        s = s[:-1]
    if s.startswith('>>'):
        note = '>>' + tail
        s = s[2:]
    elif s[0] in '><~':
        note = s[0] + tail
        s = s[1:]
    if s == 'non-exist':
        value = nan
    else:
        value = float(s)
    return value

def load_halflife():
    special = {
        'stbl': inf,
        'p-unst': 0,
        'non-exist': nan,
    }
    data = []
    with open(NUBASE_FILE) as fd:
        for line in fd:
          #try:
            # Skip the comment lines.
            if line.startswith('#'):
                continue
    
            # Make sure empty fields weren't trimmed from the fixed format lines.
            line = line[:-1] + ' '*100
    
            # Split the line into fields
            # isomer code 0:gs 1:m 2:n 3:p 4:q 5:r 6:x 7:unused 8:i 9:j
            AAA, ZZZ, code, isomer = int(line[0:3]), int(line[4:7]), int(line[7]), line[16].strip()
            symbol = line[11:16].strip()
            mass_excess = nubase_float(line[18:31]) # keV
            mass_excess_unc = nubase_float(line[31:42]) # keV
            isomer_excitation_energy = nubase_float(line[42:54]) # keV
            isomer_excitation_energy_unc = nubase_float(line[54:65])# keV
            origin = line[65:67].strip()
            isomer_unc, isomer_inv = line[67], line[68]
            halflife_str = line[69:78].strip()
            halflife_units = line[78:80].strip()
            halflife_unc_str = line[81:88].strip()
            spin_parity, isospin = line[88:102], ''
            update = line[102:103]
            discovery = line[114:118].strip()
            decay_modes = line[119:].strip()

            # Fix i11B and i35S which stuffs T=3/2 into the halflife uncertainty column
            if 'T=' in halflife_unc_str:
                spin_parity = f"{spin_parity} {halflife_unc_str}"
                halflife_unc_str = ''
                
            # Process "{halflife} {units} {unc}" fields 
            label = f"{isomer:1s}{symbol:<5s}"
            scale = UNIT_CONVERSION[halflife_units] if halflife_units else 1
            if not halflife_str:
                halflife = halflife_seconds = nan
            elif halflife_str == 'stbl':
                halflife = halflife_seconds = inf
                halflife_units = 's'
            elif halflife_str == 'p-unst':
                halflife = halflife_seconds = 0
                halflife_units = 's'
            else:
                #assert halflife_unc_str != '', f"{label}: {line[69:88]}"
                halflife = nubase_float(halflife_str)
                halflife_seconds = halflife*scale

            halflife_unc = halflife_seconds_unc = 0
            if not halflife_unc_str:
                ...
            elif halflife_unc_str[-1] in 'ys': # units on halflife uncertainty
                if halflife_unc_str[-2] in '0123456789':
                    halflife_unc = nubase_float(halflife_unc_str[:-1])
                    units = halflife_unc_str[-1]
                else:
                    halflife_unc = nubase_float(halflife_unc_str[:-2])
                    units = halflife_unc_str[-2:]
                unc_scale = UNIT_CONVERSION[units]
                value = nubase_float(halflife_unc_str[:-2])
                halflife_seconds_unc = value*unc_scale
                if halflife_units:
                    # Convert uncertainty into halflife units rather than uncertainty units
                    halflife_unc = halflife_seconds_unc/scale
                else:
                    # Note: may have "stbl" elements with lower limit estimate on halflife
                    #if halflife_str != '': print(f"{label}: {line[69:88]}")
                    #if halflife_unc_str[0] not in '<>~': print(f"{label}: {line[69:88]}")
                    halflife_units = units
                    halflife_unc = value
                    # If halflife is missing then use uncertainty as halflife
                    if not halflife_str:
                        halflife, halflife_seconds = halflife_unc, halflife_seconds_unc
            else:
                halflife_unc = nubase_float(halflife_unc_str)
                halflife_seconds_unc = halflife_unc * scale

            # Various halflife parsing integrity views
            empty = line[69:88].strip() == '' or "T=" in line[69:88]
            stable = halflife_str in ('stbl', 'p-unst')
            usual = not stable and halflife_str and halflife_unc_str and halflife_unc_str[-1] not in 'ys'
            unusual = not stable and halflife_str and halflife_unc_str and halflife_unc_str[-1] in 'ys'
            partial = not stable and not empty and not (halflife_str and halflife_unc_str)
            #if unusual:
            #    print(f"{label}: {line[69:88]} => {halflife} ± {halflife_unc} {halflife_units} => {halflife_seconds:.4g} ± {halflife_seconds_unc:.4g} s")
            
            # Sort out spin, parity and isospin
            if 'T=' in spin_parity:
                spin_parity, isospin = spin_parity.split('T=')
            spin_parity, isospin = spin_parity.strip(), isospin.strip()

            # Build the data structure representing the line.
            data.append(dict(
                AAA=AAA, Z=ZZZ, code=code, isomer=isomer, symbol=symbol,
                mass_excess=mass_excess, mass_excess_unc=mass_excess_unc,
                isomer_excitation_energy=isomer_excitation_energy,
                isomer_excitation_energy_unc=isomer_excitation_energy_unc,
                origin=origin,
                isomer_unc=isomer_unc, isomer_inv=isomer_inv,
                halflife_str=halflife_str,
                halflife=halflife, halflife_units=halflife_units, halflife_unc=halflife_unc,
                halflife_seconds=halflife_seconds, halflife_seconds_unc=halflife_seconds_unc,
                spin_parity=spin_parity, isospin=isospin,
                discovery=discovery, update=update,
                decay_modes=decay_modes,                
            ))
          #except Exception as exc: print(f"{exc} {len(line)}\n{line}")
    return data

[pkienzle]

Maybe a list of excitations with dot access for each item (gs for ground state, m, n, p, q, r, x, i, j). Then every isotope would have isotope.gs, with attributes for halflife, etc., but only those with an "m" excitation would have isotope.m.

Maybe parse decay modes into a list of (mode, %, uncertainty). Note that stable isotopes have "decay mode" IS=% giving the natural abundance for that isotope. Other modes are n, p, 2n, 2p, B+, B- (β+ and β-), A (α), IT (γ internal transition), EC (ε electron capture), and rarely 3H (H-3) and 3He (He-3). Subsequent decays (e.g., B-2n) may also be listed.

Maybe parse spin_parity into a tuple of spin-parity strings, removing ()#* and frg flags.

I don't think we need to specify excitation in the chemical formula.

[pkienzle]

Note: isotopic abundance for Tc in iupac tables is 0 (no natural abundance for any isotope), but the activation calculator uses Tc-98 = 100. That is, if Tc is given in a formula, assume it is Tc-98.

[pkienzle]

See also #105. Fixed by #107.