feat(species): added half life, spin parity, and decay modes to species database

meson.build

@@ -18,7 +18,7 @@
 #   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 #
 # *********************************************************************** #
-project('libcomposition', 'cpp', version: 'v1.0.9', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
+project('libcomposition', 'cpp', version: 'v1.1.0', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
 
 # Add default visibility for all C++ targets
 add_project_arguments('-fvisibility=default', language: 'cpp')

src/composition/include/fourdst/composition/atomicSpecies.h

@@ -15,11 +15,42 @@ namespace fourdst::atomic {
         double m_bindingEnergy; //< Binding energy
         std::string m_betaCode; //< Beta decay code
         double m_betaDecayEnergy; //< Beta decay energy
+        double m_halfLife_s; //< Half-life in seconds
+        std::string m_spinParity; //< Spin and parity
+        std::string m_decayModes; //< Decay modes
         double m_atomicMass; //< Atomic mass
         double m_atomicMassUnc; //< Atomic mass uncertainty
 
-        Species(const std::string_view name, const std::string_view el, const int nz, const int n, const int z, const int a, const double bindingEnergy, const std::string_view betaCode, const double betaDecayEnergy, const double atomicMass, const double atomicMassUnc)
-            : m_name(name), m_el(el), m_nz(nz), m_n(n), m_z(z), m_a(a), m_bindingEnergy(bindingEnergy), m_betaCode(betaCode), m_betaDecayEnergy(betaDecayEnergy), m_atomicMass(atomicMass), m_atomicMassUnc(atomicMassUnc) {};
+        Species(
+            const std::string_view name,
+            const std::string_view el,
+            const int nz,
+            const int n,
+            const int z,
+            const int a,
+            const double bindingEnergy,
+            const std::string_view betaCode,
+            const double betaDecayEnergy,
+            const double halfLife_s,
+            const std::string_view spinParity,
+            const std::string_view decayModes,
+            const double atomicMass,
+            const double atomicMassUnc
+        ) :
+        m_name(name),
+        m_el(el),
+        m_nz(nz),
+        m_n(n),
+        m_z(z),
+        m_a(a),
+        m_bindingEnergy(bindingEnergy),
+        m_betaCode(betaCode),
+        m_betaDecayEnergy(betaDecayEnergy),
+        m_halfLife_s(halfLife_s),
+        m_spinParity(spinParity),
+        m_decayModes(decayModes),
+        m_atomicMass(atomicMass),
+        m_atomicMassUnc(atomicMassUnc) {};
 
         //Copy constructor
         Species(const Species& species) {
@@ -32,6 +63,9 @@ namespace fourdst::atomic {
             m_bindingEnergy = species.m_bindingEnergy;
             m_betaCode = species.m_betaCode;
             m_betaDecayEnergy = species.m_betaDecayEnergy;
+            m_halfLife_s = species.m_halfLife_s;
+            m_spinParity = species.m_spinParity;
+            m_decayModes = species.m_decayModes;
             m_atomicMass = species.m_atomicMass;
             m_atomicMassUnc = species.m_atomicMassUnc;
         }
@@ -45,6 +79,18 @@ namespace fourdst::atomic {
             return m_atomicMassUnc;
         }
 
+        double halfLife() const {
+            return m_halfLife_s;
+        }
+
+        std::string_view spinParity() const {
+            return m_spinParity;
+        }
+
+        std::string_view decayModes() const {
+            return m_decayModes;
+        }
+
         double bindingEnergy() const {
             return m_bindingEnergy;
         }
@@ -82,7 +128,7 @@ namespace fourdst::atomic {
         }
 
         friend std::ostream& operator<<(std::ostream& os, const Species& species) {
-            os << species.m_name << " (" << species.m_atomicMass << " u)";
+            os << species.m_name;
             return os;
         }
 

src/composition/include/fourdst/composition/composition.h

@@ -61,7 +61,7 @@ namespace fourdst::composition {
      */
     struct CompositionEntry {
         std::string m_symbol; ///< The chemical symbol of the species.
-        fourdst::atomic::Species m_isotope; ///< The isotope of the species.
+        atomic::Species m_isotope; ///< The isotope of the species.
         bool m_massFracMode = true; ///< The mode of the composition entry. True if mass fraction, false if number fraction.   
 
         double m_massFraction = 0.0; ///< The mass fraction of the species.

src/composition/lib/composition.cpp

@@ -35,9 +35,10 @@
 
 namespace fourdst::composition {
 
-    CompositionEntry::CompositionEntry() : m_symbol("H-1"), m_isotope(fourdst::atomic::species.at("H-1")),
-                                           m_initialized(false) {
-    }
+    CompositionEntry::CompositionEntry() :
+    m_symbol("H-1"),
+    m_isotope(fourdst::atomic::species.at("H-1")),
+    m_initialized(false) {}
 
     CompositionEntry::CompositionEntry(const std::string& symbol, const bool massFracMode) : m_symbol(symbol), m_isotope(fourdst::atomic::species.at(symbol)), m_massFracMode(massFracMode) {
         setSpecies(symbol);

tests/composition/compositionTest.cpp

@@ -26,6 +26,12 @@ TEST_F(compositionTest, isotopeMasses) {
     EXPECT_DOUBLE_EQ(fourdst::atomic::species.at("He-4").mass(),4.0026032541300003);
 }
 
+TEST_F(compositionTest, isotopeHalfLives) {
+    EXPECT_DOUBLE_EQ(fourdst::atomic::H_1.halfLife(), std::numeric_limits<double>::infinity());
+    EXPECT_DOUBLE_EQ(fourdst::atomic::F_18.halfLife(), 6584.04);
+    EXPECT_DOUBLE_EQ(fourdst::atomic::B_20.halfLife(), 0.0);
+}
+
 TEST_F(compositionTest, constructor) {
     fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
     EXPECT_NO_THROW(fourdst::composition::Composition comp);
@@ -218,3 +224,4 @@ TEST_F(compositionTest, molarAbundance) {
     EXPECT_DOUBLE_EQ(comp1.getMolarAbundance("H-1"), 0.5/fourdst::atomic::H_1.mass());
     EXPECT_DOUBLE_EQ(comp1.getMolarAbundance("He-4"), 0.5/fourdst::atomic::He_4.mass());
 }
+

tests/composition/meson.build

@@ -11,7 +11,13 @@ foreach test_file : test_sources
   test_exe = executable(
       exe_name,
       test_file,
-      dependencies: [gtest_dep, species_weight_dep, gtest_main, composition_dep],
+      dependencies: [
+          gtest_dep,
+          species_weight_dep,
+          gtest_main,
+          composition_dep,
+          config_dep,
+      ],
       install_rpath: '@loader_path/../../src'  # Ensure runtime library path resolves correctly
   )
 

utils/atomic/convertWeightsToHeader.py

@@ -1,151 +0,0 @@
-import pandas as pd
-
-# Define fixed-width column specifications based on the format:
-# a1 (width 1), i3 (width 3), i5 (width 5), i5 (width 5), i5 (width 5),
-# 1x (skip 1), a3 (width 3), a4 (width 4), 1x (skip 1),
-# f14.6 (width 14), f12.6 (width 12), f13.5 (width 13),
-# 1x (skip 1), f10.5 (width 10), 1x (skip 1),
-# a2 (width 2), f13.5 (width 13), f11.5 (width 11),
-# 1x (skip 1), i3 (width 3), 1x (skip 1),
-# f13.6 (width 13), f12.6 (width 12)
-# Compute cumulative positions (0-indexed):
-colSpecs = [
-    (0, 1),    # control
-    (1, 4),    # NZ
-    (4, 9),    # N
-    (9, 14),   # Z
-    (14, 19),  # A
-    # skip 1 char at position 19; next field starts at 20
-    (20, 23),  # el
-    (23, 27),  # o
-    # skip 1 char at position 27; next field starts at 28
-    (28, 42),  # massExcess (f14.6)
-    (42, 54),  # massExcessUnc (f12.6)
-    (54, 67),  # bindingEnergy (f13.5)
-    # skip 1 char at position 67; next field starts at 68
-    (68, 78),  # bindingEnergyUnc (f10.5)
-    # skip 1 char at position 78; next field starts at 79
-    (79, 81),  # betaCode (a2)
-    (81, 94),  # betaDecayEnergy (f13.5)
-    (94, 105), # betaDecayEnergyUnc (f11.5)
-    # skip 1 char at position 105; next field starts at 106
-    (106, 109),# atomicMassInt (i3)
-    # skip 1 char at position 109; next field starts at 110
-    (110, 123),# atomicMassFrac (f13.6)
-    (123, 135) # atomicMassUnc (f12.6)
-]
-
-# Define column names (using camelCase for variables)
-columnNames = [
-    "control",
-    "nz",
-    "n",
-    "z",
-    "a",
-    "el",
-    "o",
-    "massExcess",
-    "massExcessUnc",
-    "bindingEnergy",
-    "bindingEnergyUnc",
-    "betaCode",
-    "betaDecayEnergy",
-    "betaDecayEnergyUnc",
-    "atomicMassInt",
-    "atomicMassFrac",
-    "atomicMassUnc"
-]
-
-def combine_atomic_mass(row):
-    """
-    Combine the integer and fractional parts of the atomic mass.
-    For example, if atomicMassInt is '1' and atomicMassFrac is '008664.91590',
-    this function returns float('1008664.91590').
-    """
-    intPart = str(row["atomicMassInt"]).strip()
-    fracPart = str(row["atomicMassFrac"]).strip()
-    try:
-        combined = int(intPart) + float(fracPart)/1e6
-        return combined
-    except ValueError:
-        return None
-
-def mkInstanceName(row):
-    """
-    Make a c++ instance name from the element and atomic number.
-    """
-    speciesName = f"{row['el'].strip()}-{row['a']}"
-    return speciesName.replace("-", "_")
-
-def formatSpecies(row):
-    """
-    Format c++ instantiation of Species struct from row data.
-    """
-    name = f"{row['el'].strip()}-{row['a']}"
-    instanceName = name.replace("-", "_")
-    nz = int(row['nz'])
-    n = int(row['n'])
-    z = int(row['z'])
-    a = int(row['a'])
-    bindingEnergy = float(row['bindingEnergy'])
-    atomicMass = float(row['atomicMass'])
-    atomicMassUnc = float(row['atomicMassUnc'])
-    NaN = "std::numeric_limits<double>::quiet_NaN()"
-    try:
-        betaDecayEnergy = float(row['betaDecayEnergy'].replace("#", "").replace("*", ""))
-    except ValueError:
-        betaDecayEnergy = NaN
-    instantiation = f"static const Species {instanceName}(\"{name}\", \"{row['el']}\", {nz}, {n}, {z}, {a}, {bindingEnergy}, \"{row['betaCode']}\", {betaDecayEnergy}, {atomicMass}, {atomicMassUnc});"
-    return instantiation, instanceName
-
-def formatSpeciesDefines(row):
-    instanceName = f"SERIF_SPECIES_{formatSpecies(row)[1]}"
-    define = f"""#ifndef {instanceName.upper()}
-    #define {instanceName.upper()}
-#endif // {instanceName.upper()}"""
-    return define
-
-def formatHeader(dataFrame):
-    """
-    Format c++ header file from DataFrame.
-    """
-    header = f"""#pragma once
-#include <unordered_map>
-#include <string_view>
-#include <string>
-#include "atomicSpecies.h"
-
-namespace fourdst::atomic {{
-    {'\n    '.join([formatSpecies(row)[0] for index, row in dataFrame.iterrows()])}
-    static const std::unordered_map<std::string, Species> species = {{
-        {'\n        '.join([f'{{"{row["el"].strip()}-{row["a"]}", {mkInstanceName(row)}}},' for index, row in dataFrame.iterrows()])}
-    }};
-}}; // namespace fourdst::atomic
-
-{'\n'.join([formatSpeciesDefines(row) for index, row in dataFrame.iterrows()])}
-"""
-    return header
-
-if __name__ == "__main__":
-    import argparse 
-    import os
-    parser = argparse.ArgumentParser(description="Convert mass data to c++ header file.")
-    parser.add_argument("input", help="Input file path.")
-    parser.add_argument("-o", "--output", help="Output file path.", default="../../assets/static/atomic/include/species.h")
-    args = parser.parse_args()
-
-
-    if not os.path.exists(args.input):
-        raise FileNotFoundError(f"File not found: {args.input}")
-
-    # Read the file (adjust the skiprows value if your header differs)
-    dataFrame = pd.read_fwf(args.input, colspecs=colSpecs, names=columnNames, skiprows=36)
-
-    # Combine the two atomic mass fields into one float column
-    dataFrame["atomicMass"] = dataFrame.apply(combine_atomic_mass, axis=1)
-    dataFrame.drop(columns=["atomicMassInt", "atomicMassFrac"], inplace=True)
-
-    # Format the header
-    header = formatHeader(dataFrame)
-    with open(args.output, "w") as f:
-        f.write(header)

utils/atomic/format.py

@@ -0,0 +1,207 @@
+import pandas as pd
+import numpy as np
+import argparse
+import os
+
+
+nubase_col_specs = [
+    (0, 3),    
+    (4, 8),    
+    (11, 16),  
+    (69, 78),  
+    (79, 81),  
+    (89, 102), 
+    (120, 209) 
+]
+nubase_column_names = [
+    "a",             
+    "ZZZi",
+    "A_El",          
+    "halfLife",
+    "halfLifeUnit",
+    "spinParity",
+    "decayModes"
+]
+
+ame_col_specs = [
+    (0, 1),
+    (1, 4),
+    (4, 9),
+    (9, 14),
+    (14, 19),
+    (20, 23),
+    (23, 27),
+    (28, 42),
+    (42, 54),
+    (54, 67),
+    (68, 78),
+    (79, 81),
+    (81, 94),
+    (94, 105),
+    (106, 109),
+    (110, 123),
+    (123, 135)
+]
+ame_column_names = [
+    "control",
+    "nz",
+    "n",
+    "z",
+    "a",
+    "el",
+    "o",
+    "massExcess",
+    "massExcessUnc",
+    "bindingEnergy",
+    "bindingEnergyUnc",
+    "betaCode",
+    "betaDecayEnergy",
+    "betaDecayEnergyUnc",
+    "atomicMassInt",
+    "atomicMassFrac",
+    "atomicMassUnc"
+]
+
+
+def combine_atomic_mass(row):
+    intPart = str(row["atomicMassInt"]).strip()
+    fracPart = str(row["atomicMassFrac"]).strip()
+    try:
+        combined = int(intPart) + float(fracPart)/1e6
+        return combined
+    except ValueError:
+        return None
+
+def mkInstanceName(row):
+    speciesName = f"{row['el'].strip()}-{row['a']}"
+    return speciesName.replace("-", "_")
+
+def convert_half_life_to_seconds(row):
+    value_str = str(row['halfLife']).strip()
+    unit = str(row['halfLifeUnit']).strip()
+
+    if value_str == 'stbl':
+        return np.inf
+
+    try:
+        value = float(value_str.replace('#', ''))
+    except ValueError:
+        return 0.0 
+
+    factors = {
+        'ys': 1e-24, 'zs': 1e-21, 'as': 1e-18, 'fs': 1e-15,
+        'ps': 1e-12, 'ns': 1e-9,  'us': 1e-6,  'ms': 1e-3,
+        's': 1.0, 'm': 60.0, 'h': 3600.0, 'd': 86400.0,
+        'y': 3.15576e7,  
+        'ky': 3.15576e10, 'My': 3.15576e13, 'Gy': 3.15576e16,
+        'Ty': 3.15576e19, 'Py': 3.15576e22, 'Ey': 3.15576e25
+    }
+    return value * factors.get(unit, 0.0)
+
+
+def formatSpecies(row):
+    name = f"{row['el'].strip()}-{row['a']}"
+    instanceName = name.replace("-", "_")
+    nz = int(row['nz'])
+    n = int(row['n'])
+    z = int(row['z'])
+    a = int(row['a'])
+    bindingEnergy = float(row['bindingEnergy'])
+    atomicMass = float(row['atomicMass'])
+    atomicMassUnc = float(row['atomicMassUnc'])
+    NaN = "std::numeric_limits<double>::quiet_NaN()"
+    
+    try:
+        betaDecayEnergy = float(row['betaDecayEnergy'].replace("#", "").replace("*", ""))
+    except (ValueError, AttributeError):
+        betaDecayEnergy = NaN
+
+    ### --- Add new NUBASE fields ---
+    halfLife_s = row.get('halfLife_s', 'std::numeric_limits<double>::infinity()')
+    if halfLife_s == np.inf:
+        halfLife_s = 'std::numeric_limits<double>::infinity()'
+        
+    spinParity = str(row.get('spinParity', '')).strip().replace('"', '\\"')
+    decayModes = str(row.get('decayModes', '')).strip().replace('"', '\\"')
+
+    instantiation = (
+        f"static const Species {instanceName}(\"{name}\", \"{row['el'].strip()}\", {nz}, {n}, {z}, {a}, "
+        f"{bindingEnergy}, \"{row['betaCode']}\", {betaDecayEnergy}, "
+        f"{halfLife_s}, \"{spinParity}\", \"{decayModes}\", "
+        f"{atomicMass}, {atomicMassUnc});"
+    )
+    return instantiation, instanceName
+
+def formatSpeciesDefines(row):
+    instanceName = f"SERIF_SPECIES_{formatSpecies(row)[1]}"
+    define = f"""#ifndef {instanceName.upper()}
+    #define {instanceName.upper()}
+#endif // {instanceName.upper()}"""
+    return define
+
+
+def formatHeader(dataFrame):
+    header = f"""#pragma once
+#include <unordered_map>
+#include <string_view>
+#include <string>
+#include <limits> // Required for std::numeric_limits
+#include "fourdst/composition/atomicSpecies.h"
+
+namespace fourdst::atomic {{
+    {'\n    '.join([formatSpecies(row)[0] for index, row in dataFrame.iterrows()])}
+    
+    static const std::unordered_map<std::string, const Species&> species = {{
+        {'\n        '.join([f'{{"{row["el"].strip()}-{row["a"]}", {mkInstanceName(row)}}},' for index, row in dataFrame.iterrows()])}
+    }};
+}}; // namespace fourdst::atomic
+{'\n'.join([formatSpeciesDefines(row) for index, row in dataFrame.iterrows()])}
+"""
+    return header
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Convert AME2020 and NUBASE2020 data to a C++ header file.")
+    parser.add_argument("ame_input", help="Input file path for AME2020 (mass.mas20).")
+    parser.add_argument("nubase_input", help="Input file path for NUBASE2020.")
+    parser.add_argument("-o", "--output", help="Output file path.", default="species.h")
+    args = parser.parse_args()
+
+    for path in [args.ame_input, args.nubase_input]:
+        if not os.path.exists(path):
+            raise FileNotFoundError(f"File not found: {path}")
+
+    ame_df = pd.read_fwf(args.ame_input, colspecs=ame_col_specs, names=ame_column_names, skiprows=36, keep_default_na=False)
+    ame_df["atomicMass"] = ame_df.apply(combine_atomic_mass, axis=1)
+    ame_df.drop(columns=["atomicMassInt", "atomicMassFrac"], inplace=True)
+    ame_df['el'] = ame_df['el'].str.strip()
+
+    nubase_df = pd.read_fwf(args.nubase_input, colspecs=nubase_col_specs, names=nubase_column_names, skiprows=1, keep_default_na=False)
+    nubase_df['a'] = pd.to_numeric(nubase_df['a'], errors='coerce')
+    nubase_df.dropna(subset=['a'], inplace=True)
+    nubase_df['a'] = nubase_df['a'].astype(int)
+    nubase_df['el'] = nubase_df['A_El'].str.extract(r'([a-zA-Z]+)')[0].str.strip()
+    nubase_df = nubase_df[nubase_df['ZZZi'].str[3:4].isin([' ', '0'])]
+    nubase_df.drop_duplicates(subset=['a', 'el'], keep='first', inplace=True)
+    nubase_df['halfLife_s'] = nubase_df.apply(convert_half_life_to_seconds, axis=1)
+    
+
+    print("--- AME DataFrame ---")
+    print(ame_df[['a', 'el']].head())
+    print("\n--- NUBASE DataFrame ---")
+    print(nubase_df[['a', 'el', 'halfLife_s']].head())
+    print("\n")
+
+    merged_df = pd.merge(ame_df, nubase_df[['a', 'el', 'halfLife_s', 'spinParity', 'decayModes']], on=['a', 'el'], how='left')
+
+    print("--- Merged DataFrame ---")
+    print(merged_df[['a', 'el', 'halfLife_s']].head(10))
+    print("\n--- Merge Stats ---")
+    print(f"Total rows in final data: {len(merged_df)}")
+    print(f"Number of rows with a valid half-life after merge: {merged_df['halfLife_s'].notna().sum()}")
+    print("-" * 20)
+
+    header = formatHeader(merged_df)
+    with open(args.output, "w") as f:
+        f.write(header)
+
+    print(f"Successfully generated C++ header at {args.output}")

utils/atomic/readme.md

@@ -9,5 +9,5 @@ In order to use this utility you will need
 
 ## Usage
 ```bash
-python convertWeightsToHeader.py <path/to/weights.dat> -o atomicWeights.h
-```
+python format.py <path/to/AME.txt> <path/to/nubase.asc> -o speciesData.h
+```