Skip to content

Commit

Permalink
Updates TDHF tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@benroberts999
benroberts999 committed Jul 12, 2024
1 parent fefa5bc commit 16c2fe5
Show file tree
Hide file tree
Showing 5 changed files with 168 additions and 345 deletions.
13 changes: 11 additions & 2 deletions src/DiracOperator/TensorOperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
10 changes: 7 additions & 3 deletions src/DiracOperator/TensorOperator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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<double> &getv() const { return m_vec; }
Expand Down
229 changes: 0 additions & 229 deletions src/ExternalField/MixedStates.tests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

// <a|Xn>
fmt::print("\n{:3s} {:3s} {:>11s} {:>11s} {:>11s} | {:>12s} {:>12s} eps\n",
"n", "a", "<a|h|n>", "TDHF", "basis", "<a|Xn>_TDHF",
"<a|Xn>_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><n||h|c> / (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";
}
}
}
10 changes: 7 additions & 3 deletions src/ExternalField/TDHF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down Expand Up @@ -245,12 +249,12 @@ std::pair<double, std::string> 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();
Expand All @@ -274,7 +278,7 @@ std::pair<double, std::string> 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;
Expand Down
Loading

0 comments on commit 16c2fe5

Please sign in to comment.