Skip to content

Commit

Permalink
Updates nuclear data based on https://www-nds.iaea.org/
Browse files Browse the repository at this point in the history 
Squashed commit of the following:

commit c70721ea0d5387651439c78a9b879821c44b43c0
Author: BM Roberts <b.roberts@uq.edu.au>
Date:   Thu Nov 2 12:13:41 2023 +1000

    Updates nuclear data, and tests

commit 3c928b2
Author: BM Roberts <b.roberts@uq.edu.au>
Date:   Wed Oct 11 20:26:32 2023 +1000

    Updated - BUT WRONG???

commit 687beca
Author: BM Roberts <b.roberts@uq.edu.au>
Date:   Wed Oct 11 16:51:00 2023 +1000

    Startung to swap

commit 50ecb7a
Author: BM Roberts <b.roberts@uq.edu.au>
Date:   Wed Oct 11 16:48:12 2023 +1000

    Update Nuclear
  • Loading branch information
@benroberts999
benroberts999 committed Nov 2, 2023
1 parent 3724cd2 commit cf787a1
Show file tree
Hide file tree
Showing 13 changed files with 1,113 additions and 1,017 deletions.
14 changes: 9 additions & 5 deletions src/DiracOperator/DiracOperator.tests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,7 @@ TEST_CASE("DiracOperator", "[DiracOperator][unit]") {

//--------------------------------------------------------------------
SECTION("hfs") {
// test data generated with "old" mu = 2.582025
std::cout << "hfs\n";
const auto data0 = std::vector{std::tuple{"1s+", "1s+", 2.2989921474e+02},
{"1s+", "3d-", -1.2365855864e+00},
Expand Down Expand Up @@ -171,13 +172,16 @@ TEST_CASE("DiracOperator", "[DiracOperator][unit]") {
{"3d-", "1s+", -3.0880838988e-09},
{"3d-", "3d-", -0.0000000000e+00}};

// test data generated with "old" mu = 2.582025
const IO::InputBlock options{""};
auto h0 = DiracOperator::generate("hfs", {"hfs", "F(r)=pointlike;"}, wf);
auto hB = DiracOperator::generate("hfs", {"hfs", "F(r)=Ball;"}, wf);
auto hS =
DiracOperator::generate("hfs", {"hfs", "F(r)=SingleParticle;"}, wf);
auto h0 = DiracOperator::generate(
"hfs", {"hfs", "F(r)=pointlike; mu=2.582025;"}, wf);
auto hB =
DiracOperator::generate("hfs", {"hfs", "F(r)=Ball; mu=2.582025;"}, wf);
auto hS = DiracOperator::generate(
"hfs", {"hfs", "F(r)=SingleParticle; mu=2.582025;"}, wf);
auto h0_au = DiracOperator::generate(
"hfs", {"hfs", "F(r)=pointlike; units=au;"}, wf);
"hfs", {"hfs", "F(r)=pointlike; units=au; mu=2.582025;"}, wf);

REQUIRE(h0->get_d_order() == 0);
REQUIRE(h0->imaginaryQ() == false);
Expand Down
8 changes: 4 additions & 4 deletions src/DiracOperator/Operators/hfs.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -291,8 +291,8 @@ generate_hfs(const IO::InputBlock &input, const Wavefunction &wf) {

const auto nuc = wf.nucleus();
const auto isotope = Nuclear::findIsotopeData(nuc.z(), nuc.a());
const auto mu = input.get("mu", isotope.mu);
const auto I_nuc = input.get("I", isotope.I_N);
const auto mu = input.get("mu", isotope.mu ? *isotope.mu : 0.0);
const auto I_nuc = input.get("I", isotope.I_N ? *isotope.I_N : 0.0);
const auto print = input.get("print", true);
const auto k = input.get("k", 1);

Expand Down Expand Up @@ -372,7 +372,7 @@ generate_hfs(const IO::InputBlock &input, const Wavefunction &wf) {
} else if (distro_type == DistroType::shell) {
Fr = Hyperfine::sphericalShell_F();
} else if (distro_type == DistroType::SingleParticle) {
const auto pi = input.get("parity", isotope.parity);
const auto pi = input.get("parity", isotope.parity ? *isotope.parity : 0);
const auto l_tmp = int(I_nuc + 0.5 + 0.0001);
auto l = ((l_tmp % 2 == 0) == (pi == 1)) ? l_tmp : l_tmp - 1;
l = input.get("l", l); // can override derived 'l' (not recommended)
Expand All @@ -390,7 +390,7 @@ generate_hfs(const IO::InputBlock &input, const Wavefunction &wf) {
}
Fr = Hyperfine::VolotkaSP_F(mu, I_nuc, l, gl, print);
} else if (distro_type == DistroType::spu) {
const auto pi = input.get("parity", isotope.parity);
const auto pi = input.get("parity", isotope.parity ? *isotope.parity : 0);
const auto u_func = input.get("u", std::string{"u1"}); // u1=(R-r)^n, u2=r^n
const bool u_option = u_func == std::string{"u1"};
const auto n = input.get("n", 0.0); // u1=(R-r)^n, u2=r^n
Expand Down
5 changes: 4 additions & 1 deletion src/HF/HartreeFock.tests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,15 +48,18 @@ TEST_CASE("HartreeFock", "[HF][HartreeFock][integration]") {
for (auto &[Atom, Core, Valence, EnergyData, E1Data, HFSData] :
UnitTest::HF_test_data::regression_test_data) {

// Fr generated with old rms value 5.6951
double r_rms = Atom == "Fl" ? 6.0 : Atom == "Fr" ? 5.6951 : -1.0;
Wavefunction wf({points, r0, rmax, b, grid_type},
{Atom, A, nucleus_type});
{Atom, A, nucleus_type, r_rms});

const auto h =
DiracOperator::hfs(1, 1.0 / 0.5, 0.0, wf.grid(),
DiracOperator::Hyperfine::pointlike_F());
const auto d = DiracOperator::E1(wf.grid());

std::cout << "\n" << wf.atom() << "\n";
std::cout << wf.nucleus() << "\n";
wf.solve_core("HartreeFock", x_Breit, Core);
wf.solve_valence(Valence);

Expand Down
52 changes: 35 additions & 17 deletions src/HF/HartreeFock_test_data.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ namespace HF_test_data {
// (r0=1.0e-7, points=10000, rmax = 175)
const auto regression_test_data = std::vector{
// atom data:
std::tuple{"Fr", "[Rn]", "8sp6d5f",
std::tuple{std::string{"Fr"}, "[Rn]", "8sp6d5f",
// energy data:
std::vector{std::tuple{"7s+", -0.131072789359},
{"8s+", -0.055958640434},
Expand Down Expand Up @@ -42,6 +42,7 @@ const auto regression_test_data = std::vector{
{"6d+", "5f-", 3.0142798909e+00},
{"6d+", "5f+", 1.3486261312e+01}},
// Hyperfine data [with mu=1, I=0.5]
// Fr generated with old rms value 5.6951
std::vector{std::tuple{"7s+", 1.3314299360e+04},
{"8s+", 3.4123493482e+03},
{"7p-", 1.4126725069e+03},
Expand Down Expand Up @@ -573,22 +574,39 @@ const auto regression_test_data = std::vector{
{"Fl",
"[Fl]",
"",
{{"1s+", -7424.427179454213}, {"2s+", -1530.733479689491},
{"2p-", -1491.664930260463}, {"2p+", -1049.948280363067},
{"3s+", -416.046497753770}, {"3p-", -396.872539425381},
{"3p+", -289.040471120399}, {"4s+", -122.079390213186},
{"3d-", -259.087300477089}, {"3d+", -241.196228217148},
{"4p-", -112.589450006828}, {"4p+", -81.421609102561},
{"5s+", -33.128814925616}, {"4d-", -66.647586594190},
{"4d+", -61.691443481028}, {"5p-", -28.775952609054},
{"5p+", -19.818023259049}, {"6s+", -6.719521710057},
{"4f-", -42.354403509323}, {"4f+", -40.905530403658},
{"5d-", -13.393057434201}, {"5d+", -12.111237575920},
{"6p-", -5.045809341525}, {"6p+", -3.011951124376},
{"7s+", -0.723561460305}, {"5f-", -4.344103942735},
{"5f+", -4.063143919599}, {"6d-", -0.957665375424},
{"6d+", -0.793491352798}, {"7p-", -0.403495134507},
{"7p+", -0.180030415175}},
// updated, due to no rms. rms = 6.0
{{"1s+", -7425.710987424376}, {"2s+", -1531.014574604990},
{"2p-", -1491.716563951863}, {"2p+", -1049.940393130142},
{"3s+", -416.115787000977}, {"3p-", -396.886913410124},
{"3p+", -289.037887135880}, {"4s+", -122.099519310650},
{"3d-", -259.084893733668}, {"3d+", -241.194085500118},
{"4p-", -112.593541247698}, {"4p+", -81.420660760941},
{"5s+", -33.134775273775}, {"4d-", -66.646742300520},
{"4d+", -61.690680688550}, {"5p-", -28.777074200928},
{"5p+", -19.817663416715}, {"6s+", -6.721052939124},
{"4f-", -42.353754507762}, {"4f+", -40.904907556323},
{"5d-", -13.392761901021}, {"5d+", -12.110968269179},
{"6p-", -5.046045597848}, {"6p+", -3.011823311115},
{"7s+", -0.723817687551}, {"5f-", -4.343911933157},
{"5f+", -4.062959193553}, {"6d-", -0.957584114798},
{"6d+", -0.793423229448}, {"7p-", -0.403520771030},
{"7p+", -0.180006075662}},
// {{"1s+", -7424.427179454213}, {"2s+", -1530.733479689491},
// {"2p-", -1491.664930260463}, {"2p+", -1049.948280363067},
// {"3s+", -416.046497753770}, {"3p-", -396.872539425381},
// {"3p+", -289.040471120399}, {"4s+", -122.079390213186},
// {"3d-", -259.087300477089}, {"3d+", -241.196228217148},
// {"4p-", -112.589450006828}, {"4p+", -81.421609102561},
// {"5s+", -33.128814925616}, {"4d-", -66.647586594190},
// {"4d+", -61.691443481028}, {"5p-", -28.775952609054},
// {"5p+", -19.818023259049}, {"6s+", -6.719521710057},
// {"4f-", -42.354403509323}, {"4f+", -40.905530403658},
// {"5d-", -13.393057434201}, {"5d+", -12.111237575920},
// {"6p-", -5.045809341525}, {"6p+", -3.011951124376},
// {"7s+", -0.723561460305}, {"5f-", -4.344103942735},
// {"5f+", -4.063143919599}, {"6d-", -0.957665375424},
// {"6d+", -0.793491352798}, {"7p-", -0.403495134507},
// {"7p+", -0.180030415175}},
{},
{}},
// test Dy (V^N), has open f-shell
Expand Down
10 changes: 5 additions & 5 deletions src/Modules/HFAnomaly.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,8 +202,8 @@ void HFAnomaly(const IO::InputBlock &input, const Wavefunction &wf) {
"(eps). Two inputs, comma separated: state (in 'short "
"symbol' form), and eps (in %). E.g.: '6s+, -0.05' [optional]"},
{"A2", "Second isotope (for differential anomaly) [optional]"},
{"Nucleus2", "Nuclear (charge) parameters for isotope 2 (see -a "
"Nucleus); uses default for A2 if blank. [optional]"},
{"Nucleus2{}", "Nuclear (charge) parameters for isotope 2 (see -a "
"Nucleus); uses default for A2 if blank. [optional]"},
{"hfs2_options{}",
"Options for HFS operator for isotope 2 (see -o hfs)"},
{"1D2_target",
Expand Down Expand Up @@ -316,17 +316,17 @@ void HFAnomaly(const IO::InputBlock &input, const Wavefunction &wf) {
nucleus.t() = *t;
}
if (rrms) {
nucleus.r_rms() = *rrms;
nucleus.set_rrms(*rrms);
}
if (c_hdr) {
// this will over-ride given rms
nucleus.r_rms() = Nuclear::rrms_formula_c_t(*c_hdr, nucleus.t());
nucleus.set_rrms(Nuclear::rrms_formula_c_t(*c_hdr, nucleus.t()));
}
// If A or given rrms are zero, explicitely set to pointlike nucleus
// This isn't required, but makes output more explicit
if (nucleus.a() == 0.0 || nucleus.r_rms() == 0.0) {
nucleus.t() = 0.0;
nucleus.r_rms() = 0.0;
nucleus.set_rrms(0.0);
nucleus.type() = Nuclear::ChargeDistro::point;
}

Expand Down
2 changes: 1 addition & 1 deletion src/Modules/isotopeShift.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@ void fieldShift(const IO::InputBlock &input, const Wavefunction &wf) {
const auto dr2 = rB * rB - r0 * r0;

auto nuc_b = wf.nucleus();
nuc_b.r_rms() = rB;
nuc_b.set_rrms(rB);

wfB.update_Vnuc(Nuclear::formPotential(nuc_b, wf.grid().r()));

Expand Down
28 changes: 14 additions & 14 deletions src/Physics/AtomData_PeriodicTable.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ static const std::string filling_order =
"6p6,7s2,5f14,6d10,7p6,8s2,6f14,7d10,8p6,9s2";

//==============================================================================
//! Atomic element data: Z, symbol, A, name
struct Element {
int Z;
std::string symbol;
Expand All @@ -24,8 +25,7 @@ struct Element {
: Z(inZ), symbol(insymbol), A(inA), name(inname) {}
};

// Default values for A for each atom.
// Goes up to E120 (Z=120)
//! Default values for A for each atom. Goes up to E120 (Z=120)
static const std::vector<Element> periodic_table = {
{1, "H", 1, "hydrogen"},
{2, "He", 4, "helium"},
Expand Down Expand Up @@ -54,30 +54,30 @@ static const std::vector<Element> periodic_table = {
{25, "Mn", 55, "manganese"},
{26, "Fe", 56, "iron"},
{27, "Co", 59, "cobalt"},
{28, "Ni", 59, "nickel"},
{29, "Cu", 64, "copper"},
{30, "Zn", 65, "zinc"},
{31, "Ga", 70, "gallium"},
{28, "Ni", 58, "nickel"},
{29, "Cu", 63, "copper"},
{30, "Zn", 64, "zinc"},
{31, "Ga", 71, "gallium"},
{32, "Ge", 73, "germanium"},
{33, "As", 75, "arsenic"},
{34, "Se", 79, "selenium"},
{35, "Br", 80, "bromine"},
{34, "Se", 77, "selenium"},
{35, "Br", 79, "bromine"},
{36, "Kr", 84, "krypton"},
{37, "Rb", 85, "rubidium"},
{38, "Sr", 88, "strontium"},
{39, "Y", 89, "yttrium"},
{40, "Zr", 91, "zirconium"},
{41, "Nb", 93, "niobium"},
{42, "Mo", 96, "molybdenum"},
{43, "Tc", 97, "technetium"},
{43, "Tc", 99, "technetium"},
{44, "Ru", 101, "ruthenium"},
{45, "Rh", 103, "rhodium"},
{46, "Pd", 106, "palladium"},
{47, "Ag", 108, "silver"},
{47, "Ag", 107, "silver"},
{48, "Cd", 112, "cadmium"},
{49, "In", 115, "indium"},
{50, "Sn", 119, "tin"},
{51, "Sb", 122, "antimony"},
{51, "Sb", 121, "antimony"},
{52, "Te", 128, "tellurium"},
{53, "I", 127, "iodine"},
{54, "Xe", 131, "xenon"},
Expand All @@ -101,9 +101,9 @@ static const std::vector<Element> periodic_table = {
{72, "Hf", 178, "hafnium"},
{73, "Ta", 181, "tantalum"},
{74, "W", 184, "tungsten"},
{75, "Re", 186, "rhenium"},
{75, "Re", 185, "rhenium"},
{76, "Os", 190, "osmium"},
{77, "Ir", 192, "iridium"},
{77, "Ir", 191, "iridium"},
{78, "Pt", 195, "platinum"},
{79, "Au", 197, "gold"},
{80, "Hg", 201, "mercury"},
Expand All @@ -122,7 +122,7 @@ static const std::vector<Element> periodic_table = {
{93, "Np", 237, "neptunium"},
{94, "Pu", 244, "plutonium"},
{95, "Am", 243, "americium"},
{96, "Cm", 247, "curium"},
{96, "Cm", 246, "curium"},
{97, "Bk", 247, "berkelium"},
{98, "Cf", 251, "californium"},
{99, "Es", 252, "einsteinium"},
Expand Down
51 changes: 30 additions & 21 deletions src/Physics/NuclearData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,13 +3,13 @@
namespace Nuclear {
//==============================================================================
Isotope findIsotopeData(int z, int a) {
if (a == 0)
return Isotope{z, a, 0, 0, 0, 0};
if (a < z)
return Isotope{z, 0, 0.0, 0, {}, {}, {}};
for (const auto &nucleus : NuclearDataTable) {
if (nucleus.Z == z && nucleus.A == a)
return nucleus;
}
return Isotope{z, a, -1, 0, 0, -1};
return Isotope{z, a, {}, {}, {}, {}, {}};
}

std::vector<Isotope> findIsotopeList(int z) {
Expand All @@ -22,41 +22,50 @@ std::vector<Isotope> findIsotopeList(int z) {
}

double find_rrms(int z, int a) {
auto nuc = findIsotopeData(z, a);
if (!nuc.r_ok())
return 0;
return nuc.r_rms;
const auto nuc = findIsotopeData(z, a);
return nuc.r_rms ? *nuc.r_rms : 0.0;
}

double find_mu(int z, int a) {
auto nuc = findIsotopeData(z, a);
if (!nuc.mu_ok())
return 0;
return nuc.mu;
return nuc.mu ? *nuc.mu : 0.0;
}

int find_parity(int z, int a) {
auto nuc = findIsotopeData(z, a);
if (!nuc.parity_ok())
return 0;
return nuc.parity;
const auto nuc = findIsotopeData(z, a);
return nuc.parity ? *nuc.parity : 0;
}

double find_spin(int z, int a) {
auto nuc = findIsotopeData(z, a);
if (!nuc.I_ok())
return -1.0;
return nuc.I_N;
const auto nuc = findIsotopeData(z, a);
return nuc.I_N ? *nuc.I_N : 0.0;
}

//==============================================================================
double approximate_r_rms(int A)
double approximate_r_rms(int A, int Z)
// Returns approximate root-mean-square charge radius in fm [1.e-15 m]
// https://www.sciencedirect.com/science/article/pii/S0092640X12000265
// https://www-nds.iaea.org/radii/
{
return (A < 10) ? 1.15 * std::pow(A, 0.3333) :
0.836 * std::pow(A, 0.3333) + 0.570;
// return (A < 10) ? 1.15 * std::pow(A, 0.3333) :
// 0.836 * std::pow(A, 0.3333) + 0.570;

const auto a13 = std::pow(A, 1.0 / 3.0);
const auto z13 = std::pow(Z, 1.0 / 3.0);

// From mathematica fit to Angeli data
if (A <= 0)
return 0.0;
double rt{0.0};
if (A % 2 == 0) {
rt = (A < 20) ? 1.15875 - 0.384013 * a13 + 1.26460 * z13 :
0.344936 + 0.554395 * a13 + 0.439873 * z13;
} else {
rt = (A < 5) ? -1.90973 + 1.99164 * a13 + 0.796396 * z13 :
0.488201 + 0.60582 * a13 + 0.333936 * z13;
}
// round to 2 digits so as not to give impression of accuracy
return (static_cast<int>(rt * 100.0)) / 100.0;
}

//==============================================================================
Expand Down
Loading

0 comments on commit cf787a1

Please sign in to comment.