Add new entries in Physical Species properties (+ refactoring of the SpeciesPhysicalProperties file)

[lucafedeli88]

This PR does mainly three things:

• It adds new species to the possible PhysicalSpecies (muon, antimuon, protium, deuterium, tritium, helium3, helium4, lithium, lithium6, lithium7, beryllium, carbon12, carbon13, nitrogen14, nitrogen15, oxygen16, oxygen17, oxygen18, fluorine, neon, neon20, neon21, neon22, aluminum, gold). All the stable isotopes with Z<=10 and selected heavier elements are provided.

• It consistently uses atomic mass units for the ion masses. Note that the mass of the proton is now slightly different from the mass of protium, which includes the mass of the bound electron. Similarly, the mass of alpha is not exactly equal to the mass of helium4

• It refactors the SpeciesPhysicalProperties source file. All the properties are now in a .cpp file and are stored using std::map

The documentation has been updated to include the new options for the physical species.

⚠️ Some tests are failing because now the hydrogen mass and the helium mass are not equal to the proton mass and the alpha mass (atomic masses include the electron contributions). We may want to discuss this point.

[ax3l]

Hi @lucafedeli88 , splendid :)

• I added the mass question to tomorrows agenda for discussion
• please use https://en.cppreference.com/w/cpp/container/unordered_map over map, since we don't need the elements to be sorted here by key

[ax3l]

@lucafedeli88 please feel free to remove the [WIP] once ready.

I think we would use this in #3153

[RemiLehe]

Could you also add neutron?

[RemiLehe]

Thanks a lot for this PR! This is very useful.
I added a few inline comments.

[RemiLehe]

I have updated the values according to this NIST page

Here is the Python script that I used in order to parse and re-format this data:

import re
from scipy.constants import m_e, m_u
import periodictable

# Extract the element names from the package `periodictable`
element_names = {}
for el in periodictable.elements:
    element_names[el.number] = el.name

# Read file (downloaded from
# https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=&ascii=ascii2&isotype=some
with open('./data.html') as f:
    text = f.read()

# Parse the file to extract the numbers
atomic_data = re.findall(r'Atomic Number = (.*)\nAtomic Symbol = (.*)\nMass Number = (.*)\nRelative Atomic Mass = (\d+.\d*).*\n.*\nStandard Atomic Weight = (.*)', text)

# Go through the isotopes and print the information in the right format.
already_printed_Z_avg = []

for line in atomic_data:
    Z = int(line[0])
    N = line[2]

    # Average property across isotopes
    if Z <= 10 or element_names[Z] in ["aluminium", "argon", "copper", "xenon", "gold"]:
        if Z not in already_printed_Z_avg:
            already_printed_Z_avg.append(Z)
            values = re.findall('\d+\.\d+', line[4])
            n_digits = len(values[0])
            m_avg = sum([ float(v) for v in values])/len(values)
            m_avg = str(m_avg)[:n_digits]
            print('{PhysicalSpecies::' + element_names[Z]
                  + ', Properties{\n     amrex::Real(' + str(m_avg)
                  +') * PhysConst::m_u,\n     amrex::Real(' + str(Z) + ') * PhysConst::q_e}},')

    # Properties specific to that isotope
    m = line[3]
    if Z <= 10:
        print('{PhysicalSpecies::' + element_names[Z] + N
              + ', Properties{\n     amrex::Real(' + str(m)
              +') * PhysConst::m_u,\n     amrex::Real(' + str(Z) + ') * PhysConst::q_e}},')

[RemiLehe]

Latest update:
The changes in this PR cause a few checksum changes (which is to be expected) but also causes the physics check of the automated test for the fusion module to fail (which is not really expected).
This will probably need to be investigated.

[NeilZaim]

Hi @RemiLehe. I've looked into the fusion test failing with this PR. Basically this is because the masses of helium and of boron11 have been changed (at the keV level, which is small but sufficient to trigger the CI test). There are a few places where this change needs to be taken into account:

• https://github.com/ECP-WarpX/WarpX/blob/542ef3c09b467974b5bd9d1e3b25a208e3296562/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H#L84
  The alpha mass needs to be replaced here with PhysConst::m_u * 4.002602_prt. (it might be a good idea to pass the mass of the product species as an argument to this function rather than to redefine it here, but this can potentially be done in a separate PR)

• https://github.com/ECP-WarpX/WarpX/blob/542ef3c09b467974b5bd9d1e3b25a208e3296562/Examples/Modules/nuclear_fusion/inputs_3d#L35
  The boron mass needs to be replaced here with 11.00930536*m_u.

• https://github.com/ECP-WarpX/WarpX/blob/542ef3c09b467974b5bd9d1e3b25a208e3296562/Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py#L69
  The boron mass needs to be replaced here with 11.00930536*scc.m_u.

• https://github.com/ECP-WarpX/WarpX/blob/542ef3c09b467974b5bd9d1e3b25a208e3296562/Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py#L71
  The alpha mass needs to be replaced here with 4.002602*scc.m_u.

I have tested with these changes locally and the CI tests almost passed, I just needed to make one more small modification here:

https://github.com/ECP-WarpX/WarpX/blob/542ef3c09b467974b5bd9d1e3b25a208e3296562/Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py#L384
I've replaced that line with return E_com_kin*(p_proton_sq>0.) (because of machine precision errors, that function was returning a small but negative number instead of 0 when p_proton_sq was 0, which made the test fail later on).

[RemiLehe]

@NeilZaim Thanks a lot for looking into this.
Do you want me to make the above-mentioned changes on this PR, or are you also able to push those changes here?

[NeilZaim]

@RemiLehe I don't think I have write access, so it's probably easier if you or @lucafedeli88 do the changes.

[RemiLehe]

OK, I have made the above-mentioned changes and have reset the benchmarks.