benroberts999 / ampsci Goto Github PK

• For hyperfine TDHF doesn't work
  • Even first-order TDHF isn't correct, so not a convergence issue
  • (Values are roughly correct size, but unstable)
  • Extremely sensitive to \delta e term - possibly just numerically unstable?
  • 'MixedStates' does work for HFS, so long as F(r) != 1 (i.e., 1/r^2 probably causes problems)
  • Issue with equations? Or just with numerics?
• Also doesn't work for QED radiative potential, which also has de term
• Seems to work for E2 - doesn't "converge", but gets correct answer

Note: Uses "Green" method to solve inhomogeneous equation twice for TDHF - once for mixed states, and once for exchange. Possibly this is the issue. Switching to new inhomogeneous Hartree-Fock method might solve problem?

Diagram RPA method does work well for these operators (though, Breit is not included in diagram RPA, so it's not a perfect workaround)

Possibly related to #26

Breit is not working with Sigma(2) after update of spline.

210117-cs-btest.txt

• Write a general modern c++ DE solver
• Allow non-local term (requires normalisation) (#11)
• Efficient + numerically stable

Currently use Green's function method; works, but not as numerically stable.
Update DiracODE to take (VxPsi) term, include into derivative.
Should involve just updating the Dirac derivative part.
This may also solve the TDHF issue #3

Cleaner way to construct general operators - make consistent

• Fix up 'generate operator'
  • Have function: takes (+oper name)

The output files can be very large, and there can be many of them
To make finding and cleaning/removing them easier, they should have the same extension.
e.g., "filename.rpad.amp" or similar

Option already there in the functions, but the way RadPot is stored makes it not so simple.
Should be a simple re-design of the class to add this extra feature.
Low priority.

• Update Coulomb tables, angular tables

Some basis instabilities

• Add other basis types
• Check for spurious states
• Sign of orbitals

Despite the observed instabilities (sometimes get wrong sign, sometimes get spurious states), the basis always passes all the tests (completeness, low-r behaviour etc. etc.)

• Fix raw B-spline class
• Fix B-spline basis:
• Correct the r0 issue; correct number of splines
• Implement Johnson version
• Use more efficient integration?
  • Store coefs instead of expand?

StructureRad routine:
At the moment, only uses TDHF method.
Also: can make use of meTable to speed up (particularly in case with RPA)

Hartree Fock

• Re-write class (make less inter-dependent)
• Option to use non-local DiracODE (#11)
• Two versions of Vdir (one in wf, one in HF) - BAD.

Wavefunction

• Major cleanup
• WF vs HF;

B-spline stability

• Also: sign of some basis functions wrong (not just w/ Breit)?
• Different versions working
• Sensitive to r0(spl); especially for higher l states, and when including Breit

• Construct correlation potential at specified energies
• Easy way to use correct correlation potential matrix..?
  • e.g., in SolveMixedStates.....
• Major cleanup..
• filename: label.sig2 -> CsI_label.sig2
• Option to use Sigma2 from Gold + rest from feyn

Spectrum energies do not perfectly match with valence energies when using Feynman Sigma.
When using Goldstone sigma, however, they do.
This implies it's probably a numerical error stemming from Feynman Sigma - but is not obvious.
Typically difference is small, and doesn't seem to negatively affect much. However, when scaling sigma to exactly reproduce energy intervals, the scaling is done for valence states, so the spectrum states will not match exactly!
Possibly related to #22
Possibly related to #26

MBPT/Feynman.hpp

E.g., For Cs with n<=2, becomes numerically unstable.
Not huge issue, since these states don't contribute much..but indicative of underlying numerical problem.

TDHF (RPA) is much slower than I think it could be when Breit is included.
Not parallelised well.

• Polarisation operator instability: (#8)
  • Fails, e.g., for Cs with n_core < 3
• Hole-particle; k=0 vs k=1?
• Overall numerical stability
• Exchange (2nd order): w1 and w1w2
• g-part for Feynman?
• Breit issue?

• Allow a basis string be specified for Structure Radiation (and diagram RPA!), similar to CI case
• Make consistent input block structure

• TDHF method doesn't work for frequency-dependent E1v operator, though diagram method does
• Possibly linked to #3 (possibly unrelated though)
• Neither method is correctly symmetric for E1v (probably $|\omega|$ vs. $\omega$)
• "Basis" method performs the same as TDHF - which is somewhat interesting
  • "Basis" method is equivalent to TDHF, but instead of solving TDHF equations, expands solutions using a basis
  • nb: "basis" method is extremely slow, and should only be used for testing

Length form (independent of $\omega$)

Velocity form ($\omega$ dependent)

Full code input/output: E1_LvsV.txt

e.g., (8.26)
Particularly for M1

Spuriously combines modules input blocks.
Supposed to combine other blocks, so that we can specify options in any order.
But don't want to combine Modules blocks (we should be able to run same block twice, with different options)

Input options to reproduce the error:

MatrixElements::E1 {}
MatrixElements::E1 {omega = each;}

Currently, cannot construct Green function with Breit operator.
This means Feynman method cannot be used if Breit is included.
Would be better to be able to include Breit into construction of Green's function.
Requires modification of Dyson equation to include Hartree-Fock-Breit potential.

Currently, we have

$$ \hat G = \hat G_0 + \hat G_0 V_x \hat G = \left[ 1 - \hat G_0 V_x \right]^{-1} , $$

where $\hat G_0$ is the Green's function for the Homogeneous Dirac-Hartree-Fock equation (excluding the non-local exchange potential, $V_x$), and

$$ V_x = - \sum_a^{\rm core} \left| a \right \rangle \hat Q \left \langle a \right | $$

is the non-local exchange potential (with $Q$ the Coulomb operator), which can be easily represented as a coordinate matrix.

The issue to include Breit is

$$ V_x \to V_x + V_{\rm Br} , $$

but the matrix $V_{\rm Br}$ is much more difficult to construct.

Probably, it requires keeping the $g$-parts of the Green's function intact to make sense.

It might be easier to form the matrix $V_{\rm Br} G_0$, rather than directly calculating $V_{\rm Br}$.

Matrix elements, frequency dependent, omega = each
Fails to work properly with omega=each for E2, though it works if i manually set omega = each.
The operator seems to work, but not the RPA part.

For Cs, 6s+ - 6p-, HF+RPA:
Module::matrixElements {
// omega = 0.0417522;
omega = each;
operator = E1;
options = [gauge = vform;];
rpa_diagram = true;
}

For hyperfine: possibly linked to diagram method?