diff --git a/src/DiracOperator/TensorOperator.cpp b/src/DiracOperator/TensorOperator.cpp index b69e818dc..ed229ff79 100644 --- a/src/DiracOperator/TensorOperator.cpp +++ b/src/DiracOperator/TensorOperator.cpp @@ -14,12 +14,21 @@ namespace DiracOperator { //============================================================================== bool TensorOperator::isZero(const int ka, int kb) const { // checks m_rank and m_parity - if (m_rank < std::abs(Angular::twoj_k(ka) - Angular::twoj_k(kb)) / 2) + + // if (m_rank < std::abs(Angular::twoj_k(ka) - Angular::twoj_k(kb)) / 2) + // return true; + + if (Angular::triangle(Angular::twoj_k(ka), Angular::twoj_k(kb), 2 * m_rank) == + 0) return true; + if ((m_parity == Parity::even) != (Angular::parity_k(ka) == Angular::parity_k(kb))) return true; - return (angularF(ka, kb) == 0); + + return false; + // can remove this?? + // return (angularF(ka, kb) == 0); } bool TensorOperator::isZero(const DiracSpinor &Fa, diff --git a/src/DiracOperator/TensorOperator.hpp b/src/DiracOperator/TensorOperator.hpp index 6ccce7825..a168137af 100644 --- a/src/DiracOperator/TensorOperator.hpp +++ b/src/DiracOperator/TensorOperator.hpp @@ -118,7 +118,7 @@ class TensorOperator { opC(RorI), m_freqDependantQ(freq_dep), m_constant(constant), - m_vec(inv){}; + m_vec(inv) {}; public: virtual ~TensorOperator() = default; @@ -146,13 +146,17 @@ class TensorOperator { bool selectrion_rule(int twoJA, int piA, int twoJB, int piB) const { if (twoJA == twoJB && twoJA == 0.0) return false; - if (2 * m_rank < std::abs(twoJA - twoJB)) + + if (Angular::triangle(twoJA, twoJB, 2 * m_rank) == 0) return false; + + // if (2 * m_rank < std::abs(twoJA - twoJB)) + // return false; return (m_parity == Parity::even) == (piA == piB); } //! Update frequency for frequency-dependant operators. - virtual void updateFrequency(const double){}; + virtual void updateFrequency(const double) {}; //! Returns a const ref to vector v const std::vector &getv() const { return m_vec; } diff --git a/src/ExternalField/MixedStates.tests.cpp b/src/ExternalField/MixedStates.tests.cpp index 3a08160dd..5d98ce1ab 100644 --- a/src/ExternalField/MixedStates.tests.cpp +++ b/src/ExternalField/MixedStates.tests.cpp @@ -83,232 +83,3 @@ TEST_CASE("External Field: Mixed-states", } std::cout << "\nMixed-States summary:\n" << summary.str() << "\n"; } - -//============================================================================== -TEST_CASE("External Field: Mixed-states2", - "[ExternalField][MixedStates][integration][!mayfail]") { - std::cout << "\n----------------------------------------\n"; - std::cout << "External Field: Mixed-states2\n"; - - // This test illucidates the issue with TDHF for hyperfine: - - // Create wavefunction object - Wavefunction wf({2000, 1.0e-5, 80.0, 20.0, "loglinear", -1.0}, - {"Rb", -1, "Fermi", -1.0, -1.0}, 1.0); - wf.solve_core("HartreeFock", 0.0, "[Kr]"); - wf.solve_valence("5sp"); - - SplineBasis::Parameters bspl_param; - { - bspl_param.states = "80spdfg"; - bspl_param.n = 90; - bspl_param.k = 9; - bspl_param.r0 = 1.0e-6; - bspl_param.reps = 0.0; - bspl_param.rmax = 60.0; - bspl_param.positronQ = false; - } - - wf.formBasis(bspl_param); - - for (std::string_view oper : {"E1", "hfs"}) { - std::cout << "\n" << oper << "\n"; - auto h = DiracOperator::generate(oper, {}, wf); - - ExternalField::TDHF tdhf(h.get(), wf.vHF()); - ExternalField::TDHFbasis basis(h.get(), wf.vHF(), wf.basis()); - - const auto r = wf.grid().r(); - - for (auto omega : {0.0, 0.05}) { - - tdhf.clear(); - basis.clear(); - // 1 iteration: i.e., just solve once, without iterating - tdhf.solve_core(omega, 1, false); - basis.solve_core(omega, 1, false); - - if (oper != "E1" && omega != 0.0) - continue; - for (auto &Fc : wf.core()) { - for (const auto &Fm : wf.core()) { - if (Fc <= Fm || h->isZero(Fm, Fc) || std::abs(Fm.n() - Fc.n()) > 1) - continue; - if (omega == 0.0 && Fc == Fm) - continue; - - auto hmc1 = h->reducedME(Fm, Fc); - if (std::abs(hmc1) < 1.0e-5) - continue; - - const auto hFc = h->reduced_rhs(Fm.kappa(), Fc); - auto imag = h->imaginaryQ(); - auto de = Fc.kappa() == Fm.kappa() && !imag ? Fc * hFc : 0.0; - const auto &dF = ExternalField::solveMixedState( - Fc, omega, hFc - de * Fc, wf.vHF()); - - auto dF2 = - tdhf.get_dPsi_x(Fc, ExternalField::dPsiType::X, hFc.kappa()); - auto dF3 = - basis.get_dPsi_x(Fc, ExternalField::dPsiType::X, hFc.kappa()); - - auto hmc2 = (Fm * dF) * (Fc.en() - Fm.en() + omega); - auto eps = std::abs((hmc1 - hmc2) / hmc1); - auto del = std::abs((hmc1 - hmc2)); - eps = std::min(eps, del); - std::cout << Fm << "|" << Fc << " "; - printf("w=%.2f %+.3e %+.3e ", omega, hmc1, hmc2); - - auto hmc3 = (Fm * dF2) * (Fc.en() - Fm.en() + omega); - auto hmc4 = (Fm * dF3) * (Fc.en() - Fm.en() + omega); - - printf("%+.3e %+.3e %.1e", hmc3, hmc4, eps); - if (eps > 1.0e-2) - std::cout << " ****"; - std::cout << "\n"; - - auto n1 = dF2.norm(); - auto n2 = dF3.norm(); - auto d3 = (dF2 - dF3).norm() / dF2.norm(); - printf(" %+.3e %+.3e %.1e", n1, n2, d3); - if (d3 > 1.0e-2) - std::cout << " ******"; - std::cout << "\n"; - - auto dv1 = tdhf.dV(Fc, Fm); - auto dv2 = basis.dV(Fc, Fm); - // auto dv3 = basis.dV1(Fc, Fm); - auto eps3 = std::abs((dv1 - dv2) / dv1); - printf(" %+.3e %+.3e %.1e", dv1, dv2, eps3); - if (eps3 > 1.0e-2) - std::cout << " ******"; - std::cout << "\n"; - } - } - - double worst = 0.0; - for (auto &Fc : wf.core()) { - const auto &dF_t = tdhf.get_dPsis(Fc, ExternalField::dPsiType::X); - const auto &dF_b = basis.get_dPsis(Fc, ExternalField::dPsiType::X); - const auto &dY_t = tdhf.get_dPsis(Fc, ExternalField::dPsiType::Y); - const auto &dY_b = basis.get_dPsis(Fc, ExternalField::dPsiType::Y); - assert(dF_t.size() == dF_b.size()); - for (std::size_t ix = 0; ix < dF_t.size(); ix++) { - const auto &dF_tx = dF_t[ix]; - const auto &dF_bx = dF_b[ix]; - const auto &dY_tx = dY_t[ix]; - const auto &dY_bx = dY_b[ix]; - assert(dF_tx.kappa() == dF_bx.kappa()); - std::cout << Fc << "/" << dF_tx << ": " << dF_tx.norm() << " " - << dF_bx.norm() << " | "; - std::cout << dY_tx.norm() << " " << dY_bx.norm() << "\n"; - auto eps = std::abs(dF_tx.norm() / dF_bx.norm() - 1.0); - if (eps > worst) { - worst = eps; - } - } - } - REQUIRE(worst < 1.0e-2); - } - } -} - -//============================================================================== -TEST_CASE("External Field: Mixed-states3", - "[ExternalField][MixedStates][integration][!mayfail][XYZ]") { - std::cout << "\n----------------------------------------\n"; - std::cout << "External Field: Mixed-states3\n"; - - // This test illucidates the issue with TDHF for hyperfine: - - // Create wavefunction object - Wavefunction wf({4000, 1.0e-6, 90.0, 20.0, "loglinear", -1.0}, - {"Rb", -1, "Fermi", -1.0, -1.0}, 1.0); - wf.solve_core("HartreeFock", 0.0, "[Kr]"); - wf.solve_valence("5sp"); - - SplineBasis::Parameters bspl_param; - { - bspl_param.states = "85spdfg"; - bspl_param.n = 90; - bspl_param.k = 9; - bspl_param.r0 = 1.0e-5; - bspl_param.reps = 0.0; - bspl_param.rmax = 75.0; - bspl_param.positronQ = false; - } - - wf.formBasis(bspl_param); - - for (std::string_view oper : {"E1", "hfs"}) { - std::cout << "\n" << oper << "\n"; - auto h = DiracOperator::generate(oper, {}, wf); - - // const auto r = wf.grid().r(); - - double omega = 0.00; - - // - fmt::print("\n{:3s} {:3s} {:>11s} {:>11s} {:>11s} | {:>12s} {:>12s} eps\n", - "n", "a", "", "TDHF", "basis", "_TDHF", - "_basis"); - - const auto states = qip::merge(wf.core(), wf.valence()); - - // use basis states, - - for (auto Fm : states) { - for (auto Fc : wf.core()) { - - // Use basis states, instead of HF, for "exact" orthogonality - // Fm = *DiracSpinor::find(Fm.n(), Fm.kappa(), wf.basis()); - // Fc = *DiracSpinor::find(Fc.n(), Fc.kappa(), wf.basis()); - - if (Fc > Fm || h->isZero(Fm, Fc)) - continue; - - const auto hFc = h->reduced_rhs(Fm.kappa(), Fc); - auto imag = h->imaginaryQ(); - auto de = Fc.kappa() == Fm.kappa() && !imag ? Fc * hFc : 0.0; - - // Solves: (H - e - de)*dF = -h * F - const auto &dF_t = - ExternalField::solveMixedState(Fc, omega, hFc - de * Fc, wf.vHF()); - - // Directly finds: dF = \sum_n |n> / (ec - en + w) - const auto &dF_b = - ExternalField::solveMixedState_basis(Fc, hFc, omega, wf.basis()); - - auto h0 = Fm * hFc; - auto h1 = (Fm * dF_t) * (Fc.en() - Fm.en() + omega); - auto h2 = (Fm * dF_b) * (Fc.en() - Fm.en() + omega); - - if (Fm == Fc) { - // Fm*dF is not defined in this case! - // Should be zero exactly from sum-over-state method - isn't, because - // states aren't exactly orthogonal. Should diverge in TDHF case? - fmt::print( - "{:3s} {:3s} {:11.4e} {:^11s} {:^11s} | {:11.4e} {:11.4e}\n", - Fc.shortSymbol(), Fm.shortSymbol(), h0, "-----", "-----", - Fm * dF_t, Fm * dF_b); - } else { - const auto eps = - 0.5 * (std::abs(h1 / h0 - 1.0) + std::abs(h2 / h0 - 1.0)); - fmt::print("{:3s} {:3s} {:11.4e} {:11.4e} {:11.4e} | {:11.4e} " - "{:11.4e} {:.0e}\n", - Fc.shortSymbol(), Fm.shortSymbol(), h0, h1, h2, Fm * dF_t, - Fm * dF_b, eps); - - // base target on overlap of wavefunctions: - const auto dn = std::abs(Fc.n() - Fm.n()); - const auto target = dn == 0 ? 1.0e-4 : - dn == 1 ? 1.0e-3 : - dn == 2 ? 1.0e-3 : - 1.0e-2; - REQUIRE(eps < target); - } - } - std::cout << "\n"; - } - } -} diff --git a/src/ExternalField/TDHF.cpp b/src/ExternalField/TDHF.cpp index f75f6fe0f..1e225f5df 100644 --- a/src/ExternalField/TDHF.cpp +++ b/src/ExternalField/TDHF.cpp @@ -47,6 +47,10 @@ void TDHF::initialise_dPsi() { for (int tj = tjmin; tj <= tjmax; tj += 2) { const auto l_minus = (tj - 1) / 2; const auto pi_chla = Angular::parity_l(l_minus) * pi_ch; + + if (Angular::triangle(tj_c, tj, 2 * m_rank) == 0) + continue; + const auto l = (pi_chla == 1) ? l_minus : l_minus + 1; const auto kappa = Angular::kappa_twojl(tj, l); m_X[ic].emplace_back(0, kappa, Fc.grid_sptr()); @@ -245,12 +249,12 @@ std::pair TDHF::tdhf_core_it(double omega, const auto &oldX = m_X[ib][ibeta]; const auto t_DdF2 = (Xx - oldX).norm2(); - const auto t_dF2 = Xx.norm2(); + const auto t_dF2 = 0.5 * (Xx + oldX).norm2(); DdF2 += t_DdF2; dF2 += t_dF2; // find worst orbital - const auto t_eps = t_DdF2 / t_dF2; + const auto t_eps = t_dF2 == 0.0 ? 1.0 : t_DdF2 / t_dF2; if (t_eps > epss[ib].first) { epss[ib].first = t_eps; epss[ib].second = Fb.shortSymbol() + "," + Xx.shortSymbol(); @@ -274,7 +278,7 @@ std::pair TDHF::tdhf_core_it(double omega, //============================================================================== void TDHF::solve_core(const double omega, int max_its, const bool print) { - const double converge_targ = 1.0e-8; + const double converge_targ = 1.0e-9; const auto eta_damp0 = 0.35; auto eta_damp = eta_damp0; diff --git a/src/ExternalField/TDHF.tests.cpp b/src/ExternalField/TDHF.tests.cpp index 351e4f757..b1c9fcbf3 100644 --- a/src/ExternalField/TDHF.tests.cpp +++ b/src/ExternalField/TDHF.tests.cpp @@ -3,6 +3,7 @@ #include "IO/ChronoTimer.hpp" #include "Wavefunction/Wavefunction.hpp" #include "catch2/catch.hpp" +#include "fmt/format.hpp" #include "qip/Vector.hpp" #include #include @@ -63,15 +64,128 @@ TEST_CASE("External Field: TDHF - basic unit tests", } //============================================================================== -namespace helper { +struct TestData { + std::string a, b; + double h, rpa; +}; + +// All data with: Omega = 0.0 or 0.1 + +std::vector dzuba_e1_0{ + {"6p-", "6s+", 5.277685, 4.974416}, {"6p-", "7s+", -4.413137, -4.449363}, + {"6p+", "6s+", 7.426433, 7.013095}, {"6p+", "7s+", -6.671013, -6.712218}, + {"7p-", "6s+", 0.3717475, 0.2387373}, {"7p-", "7s+", 11.00887, 10.92106}, + {"7p+", "6s+", 0.6947538, 0.5087653}, {"7p+", "7s+", 15.34478, 15.22744}, + {"5d-", "6p-", -8.978332, -8.639903}, {"5d-", "6p+", -4.062459, -3.916489}, + {"5d-", "7p-", 4.039455, 4.147366}, {"5d-", "7p+", 1.688036, 1.736073}, + {"5d+", "6p+", -12.18643, -11.75336}, {"5d+", "7p+", 5.024624, 5.166189}, + {"4f-", "5d-", -10.65973, -10.52085}, {"4f-", "5d+", -2.840239, -2.803322}, + {"4f+", "5d+", -12.70344, -12.5383}}; + +std::vector dzuba_e1_w{ + {"6p-", "6s+", 5.277685, 4.971342}, {"6p-", "7s+", -4.413137, -4.451618}, + {"6p+", "6s+", 7.426433, 7.009477}, {"6p+", "7s+", -6.671013, -6.715528}, + {"7p-", "6s+", 0.3717475, 0.2377626}, {"7p-", "7s+", 11.00887, 10.92009}, + {"7p+", "6s+", 0.6947538, 0.5077115}, {"7p+", "7s+", 15.34478, 15.22634}, + {"5d-", "6p-", -8.978332, -8.627918}, {"5d-", "6p+", -4.062459, -3.911142}, + {"5d-", "7p-", 4.039455, 4.147084}, {"5d-", "7p+", 1.688036, 1.735976}, + {"5d+", "6p+", -12.18643, -11.73544}, {"5d+", "7p+", 5.024624, 5.164451}, + {"4f-", "5d-", -10.65973, -10.51826}, {"4f-", "5d+", -2.840239, -2.802684}, + {"4f+", "5d+", -12.70344, -12.53541}}; + +std::vector dzuba_e2_0{ + {"6p+", "6p-", 68.48282, 68.23295}, {"6p+", "6p+", -70.2262, -69.97999}, + {"7p-", "6p+", 49.3854, 49.48805}, {"7p+", "6p-", -42.52614, -42.63387}, + {"7p+", "6p+", 46.81519, 46.9199}, {"7p+", "7p-", 299.8769, 299.8039}, + {"7p+", "7p+", -305.114, -305.0414}, {"5d-", "6s+", -43.84649, -43.63158}, + {"5d-", "7s+", 80.74682, 80.78614}, {"5d-", "5d-", -67.06309, -66.83697}, + {"5d+", "6s+", -53.71202, -53.47216}, {"5d+", "7s+", 98.51161, 98.55041}, + {"5d+", "5d-", 43.84507, 43.70795}, {"5d+", "5d+", -87.57531, -87.30408}}; + +std::vector dzuba_e2_w{ + {"6p+", "6p-", 68.48282, 68.23044}, {"6p+", "6p+", -70.2262, -69.97613}, + {"7p-", "6p+", 49.3854, 49.49049}, {"7p+", "6p-", -42.52614, -42.63489}, + {"7p+", "6p+", 46.81519, 46.92167}, {"7p+", "7p-", 299.8769, 299.8032}, + {"7p+", "7p+", -305.114, -305.0403}, {"5d-", "6s+", -43.84649, -43.62903}, + {"5d-", "7s+", 80.74682, 80.7866}, {"5d-", "5d-", -67.06309, -66.83467}, + {"5d+", "6s+", -53.71202, -53.46243}, {"5d+", "7s+", 98.51161, 98.54769}, + {"5d+", "5d-", 43.84507, 43.70911}, {"5d+", "5d+", -87.57531, -87.30183}}; + +std::vector dzuba_m1_0{{"6p+", "6p-", 1.153521, 1.153505}, + {"7p-", "6p+", 0.03014502, 0.03016005}, + {"7p+", "6p-", 0.02855965, 0.02855698}, + {"7p+", "7p-", 1.153342, 1.153337}, + {"5d+", "5d-", 1.549161, 1.54951}, + {"4f+", "4f-", 1.851639, 1.85164}}; + +std::vector dzuba_m1_w{{"6p+", "6p-", 1.153511, 1.153558}, + {"7p-", "6p+", 0.03013823, 0.0301901}, + {"7p+", "6p-", 0.0285655, 0.02859815}, + {"7p+", "7p-", 1.153301, 1.153316}, + {"5d+", "5d-", 1.549149, 1.549712}, + {"4f+", "4f-", 1.851568, 1.85157}}; + +std::vector dzuba_pnc{{"6p-", "6s+", -5.72823E-4, -7.270836e-4}, + {"6p-", "7s+", -3.002691e-4, -3.794671e-4}, + {"7p-", "6s+", -3.429023e-4, -4.330309e-4}, + {"7p-", "7s+", -1.797466e-4, -2.260384e-4}}; + +// with g = 1.0; +std::vector dzuba_hfs{{"6s+", "6s+", 1943.394, 2342.449}, + {"7s+", "6s+", 1018.711, 1227.127}, + {"7s+", "7s+", 533.9993, 642.53}, + {"6p-", "6p-", 218.2665, 273.257}, + {"6p+", "6p+", 32.41904, 58.06549}, + {"7p-", "6p-", 130.5002, 162.7951}, + {"7p-", "7p-", 78.14994, 97.1154}, + {"7p+", "6p+", 19.46554, 34.75759}, + {"7p+", "7p+", 11.71131, 20.81699}, + {"5d-", "5d-", 24.71128, 21.90427}, + {"5d+", "5d+", 10.1202, -33.08472}, + {"4f-", "4f-", 0.05106033, 0.04225674}, + {"4f+", "4f+", 0.02838632, -0.01949691}}; + +//-/-/----------------------------------------------------------------/////----- +void test_RPA(const Wavefunction &wf, DiracOperator::TensorOperator &h, + double omega, const std::vector &test_data, + double target_1, double target_2) { + std::cout << "\n"; + ExternalField::TDHF dV(&h, wf.vHF()); + if (h.freqDependantQ()) { + h.updateFrequency(omega); + } + dV.solve_core(omega); + + for (const auto &[a, b, h0, rpa0] : test_data) { + const auto &Fa = *wf.getState(a); + const auto &Fb = *wf.getState(b); + + const auto f = h.name() == "hfs1" ? + DiracOperator::Hyperfine::convert_RME_to_AB( + 1, Fa.kappa(), Fb.kappa()) : + 1.0; + const auto hab = f * h.reducedME(Fa, Fb); + const auto dv = f * dV.dV(Fa, Fb); + + // don't worry about sign (tested elsewhere) + const auto s1 = hab / std::abs(hab); + const auto s2 = h0 / std::abs(h0); -// Calculates core-polarisation correction to matrix elements between all -// valence states, returns a vector of {states(string), value} -inline std::vector> -dV_result(const Wavefunction &wf, const DiracOperator::TensorOperator &h, - double ww); + const auto delta1 = (s1 * (hab + dv) - s2 * (rpa0)); + const auto eps1 = std::abs(delta1 / rpa0); -} // namespace helper + const auto delta2 = (s1 * dv - s2 * (rpa0 - h0)); + const auto eps2 = std::abs(delta2 / (rpa0 - h0)); + + fmt::print("{:3s} {:3s} {:12.5e} [{:12.5e}] {:.0e} {:12.5e} [{:12.5e}] " + "{:.0e}\n", + a, b, s1 * (hab + dv), s2 * rpa0, eps1, s1 * dv, + s2 * (rpa0 - h0), eps2); + + REQUIRE(eps1 < target_1); + REQUIRE(eps2 < target_2); + } +} //============================================================================== //============================================================================== @@ -85,110 +199,31 @@ TEST_CASE("External Field: TDHF (RPA)", // Create wavefunction object, solve HF for core+valence Wavefunction wf({6000, 1.0e-6, 175.0, 20.0, "loglinear", -1.0}, - {"Cs", 133, "Fermi", -1.0, -1.0}, 1.0); + {"Cs", 133, "Fermi", 4.8041, 2.0}, 1.0); wf.solve_core("HartreeFock", 0.0, "[Xe]"); wf.solve_valence("7sp5d4f"); - // Lambda to compare mine to test data - auto cmpr = [](const auto &ds, const auto &v) { - return (ds.second - v) / v; - }; - - { // E1, w = 0.0 - const auto ww = 0.00; - auto h = DiracOperator::E1(wf.grid()); - // Expected values from V. Dzuba code: - std::vector expected_VD = { - -0.303269160, 0.413337460, -0.133010110, 0.185988360, -0.036226022, - 0.041204439, -0.087804876, 0.117345560, -0.303269160, -0.036226022, - -0.338428260, -0.413337460, -0.041204439, 0.145969420, -0.433068210, - -0.133010110, -0.087804876, -0.107910900, -0.185988360, -0.117345560, - 0.048037122, -0.141564780, 0.338428260, 0.145969420, 0.107910900, - 0.048037122, -0.138879730, 0.433068210, 0.141564780, 0.036916812, - -0.165144850, 0.138879730, 0.036916812, 0.165144850}; - // -ve: |e|r vs er = -|e|r - qip::scale(&expected_VD, -1.0); - // sort to allow easy comparison: - std::sort(begin(expected_VD), end(expected_VD)); - - const auto result = helper::dV_result(wf, h, ww); - const auto [eps, at] = qip::compare(result, expected_VD, cmpr); - const std::string worst = at == result.end() ? "" : at->first; - // pass &= qip::check_value(&obuff, "RPA E1 w=0 " + worst, eps, - // 0.0, 5.0e-5); - std::cout << "TDHF: RPA E1 w=0 " << worst << " " << eps << "\n"; - REQUIRE(std::abs(eps) < 5.0e-5); - } - - { // E1, w = 0.05 - const auto ww = 0.05; - auto h = DiracOperator::E1(wf.grid()); - std::vector expected_VD = { - -0.303213300, 0.412942330, -0.132817610, 0.185545850, -0.037186418, - 0.042653995, -0.087825599, 0.117263410, -0.303213300, -0.037186418, - -0.342497480, -0.412942330, -0.042653995, 0.147791470, -0.439517730, - -0.132817610, -0.087825599, -0.107218100, -0.185545850, -0.117263410, - 0.047727162, -0.139992440, 0.342497480, 0.147791470, 0.107218100, - 0.047727162, -0.139387570, 0.439517730, 0.139992440, 0.037028834, - -0.165663560, 0.139387570, 0.037028834, 0.165663560}; - // -ve: |e|r vs er = -er - qip::scale(&expected_VD, -1.0); - // sort to allow easy comparison: - std::sort(begin(expected_VD), end(expected_VD)); - - const auto result = helper::dV_result(wf, h, ww); - const auto [eps, at] = qip::compare(result, expected_VD, cmpr); - const std::string worst = at == result.end() ? "" : at->first; - // pass &= - // qip::check_value(&obuff, "RPA E1 w=0.05 " + worst, eps, 0.0, 5.0e-5); - std::cout << "TDHF: RPA E1 w=0.05 " << worst << " " << eps << "\n"; - REQUIRE(std::abs(eps) < 5.0e-5); - } - - { // PNC, w = 0.0 - // Note: even zero-order PNC disagrees at ~5th digit - possibly due to c,t? - const auto ww = 0.0; - const auto c = Nuclear::c_hdr_formula_rrms_t(wf.get_rrms()); - auto h = DiracOperator::PNCnsi(c, Nuclear::default_t, wf.grid()); - std::vector expected_VD = { - 1.5428e-04, 9.0137e-05, 7.9206e-05, 4.6296e-05, -1.5428e-04, - -7.9206e-05, 4.9201e-05, -9.0137e-05, -4.6296e-05, 3.0130e-05, - -4.9201e-05, -3.0130e-05, 3.3958e-06, -3.3958e-06}; - std::sort(begin(expected_VD), end(expected_VD)); - - const auto result = helper::dV_result(wf, h, ww); - const auto [eps, at] = qip::compare(result, expected_VD, cmpr); - const std::string worst = at == result.end() ? "" : at->first; - // pass &= qip::check_value(&obuff, "RPA PNC w=0 " + worst, eps, - // 0.0, 5.0e-4); - std::cout << "TDHF: RPA PNC w=0 " << worst << " " << eps << "\n"; - REQUIRE(std::abs(eps) < 5.0e-4); - } - } -} + DiracOperator::E1 dE1(wf.grid()); + DiracOperator::Ek dE2(wf.grid(), 2); + DiracOperator::M1 m1(wf.grid(), wf.alpha()); + DiracOperator::PNCnsi dpnc(5.67073, 2.3, wf.grid()); + DiracOperator::hfs hfs(1, 1.0, 0.0, wf.grid(), + DiracOperator::Hyperfine::pointlike_F()); -//============================================================================== -inline std::vector> -helper::dV_result(const Wavefunction &wf, - const DiracOperator::TensorOperator &h, double ww) { - - // Form TDHF (RPA) object for this operator - auto dV = ExternalField::TDHF(&h, wf.vHF()); - // Solve set of TDHF equations for core, with frequency ww - const auto max_iterations = 150; - const auto print_details = true; - dV.solve_core(ww, max_iterations, print_details); - - std::vector> result; - for (const auto &Fv : wf.valence()) { - for (const auto &Fm : wf.valence()) { - if (h.isZero(Fv.kappa(), Fm.kappa())) - continue; - - result.emplace_back(Fv.shortSymbol() + Fm.shortSymbol(), dV.dV(Fv, Fm)); - } + test_RPA(wf, dE1, 0.0, dzuba_e1_0, 1.0e-4, 1.0e-4); + + test_RPA(wf, dE1, 0.1, dzuba_e1_w, 1.0e-4, 1.0e-4); + + // RPA is very small for E2, so test less stringent (simply, number of digits!) + test_RPA(wf, dE2, 0.0, dzuba_e2_0, 1.0e-5, 1.0e-3); + test_RPA(wf, dE2, 0.1, dzuba_e2_w, 1.0e-5, 1.0e-3); + + // RPA is extremely small for M1, so test less stringent (simply, number of digits!) + test_RPA(wf, m1, 0.0, dzuba_m1_0, 1.0e-4, 3.0e-1); + test_RPA(wf, m1, 0.1, dzuba_m1_w, 1.0e-4, 3.0e-1); + + test_RPA(wf, hfs, 0.0, dzuba_hfs, 1.0e-3, 1.0e-3); + + test_RPA(wf, dpnc, 0.0, dzuba_pnc, 1.0e-3, 1.0e-3); } - std::sort(begin(result), end(result), - [](auto a, auto b) { return a.second < b.second; }); - return result; }