benroberts999/ampsci

A c++ program for high-precision atomic structure calculations of one and two valence systems. Uses Hartree-Fock + correlation potential method. Can calculate ionisation cross sections with large of energy/momentum transfer.

ampsci

Atomic Many-body Perturbation theory in the Screened Coulomb Interaction

ampsci is a c++ program for high-precision atomic structure calculations of one and two valence atomic systems, developed and maintained by Benjamin M. Roberts, University of Queensland, Australia

It solves the correlated Dirac equation using the Hartree-Fock + correlation potential method (based on Dzuba-Flambaum-Sushkov method) to produce a set of atomic wavefunctions and energies. For two-valence systems, uses the CI+MBPT (Configuration Interaction with many-body perturbation theory) method. The method is fully relativistic, includes electron correlations, all-orders screening and hole-particle interaction, finite-nuclear size, Breit interaction, radiative QED effects, RPA for matrix elements, and structure radiation/renormalisation. QED is included via the Flambaum-Ginges radiative potential method. Can solve for continuum states with high energy, and calculate ionisation cross sections with large energy/momentum transfer.

Designed to be fast, accurate, and easy to use. The "modules" system (see doc/modules.md) makes it simple to add your own routines to use the atomic wavefunctions to calculate whatever properties you may be interested in.

• The code is on GitHub: github.com/benroberts999/ampsci
• See ampsci.dev/ for full documentation
• A full description of the physics methods and approximations, including references, is given in the physics documentation: ampsci.pdf.
• Important: this is a pre-release version of the code: not fully tested or documented, and should not be used for publishable calculations (without consultation)

---

Compilation and usage

Quick start

• The setup.sh bash script should compile and build ampsci (uses Make)
  • It uses only defaults, and may not work on all systems. See documentation for full guide to compilation.
  • It assumes dependencies have already been installed. If not, see next:
• The install-dependencies.sh bash script should install all required dependencies.
  • It uses only defaults, and may not work on all systems.
• If there are issue with compilation
  • See full compilation instructions: doc/compilation.md
• Check out the example input files to get running:
  • doc/examples/ampsci.in -- an example/template input file
  • More: in doc/examples/ there are several example input files, along with the expected output; use these to test if everything is working

Tutorials/examples

• The fastest way to get familiar with ampsci is to follow the tutorials
• A basic step-by-step tutorial: doc/tutorial.md
• More advanced tutorials follow: doc/tutorial_advanced.md, doc/tutorial_CI.md
• See which ampsci input options are available: ./ampsci -a
  • See available input options for each input block by following with its name
  • e.g., ./ampsci -a HartreeFock
• Check which Modules are aviable: ./ampsci -m
  • See available input options for each module by following with its name
  • e.g., ./ampsci -m MatrixElements
• Check which operators are aviable: ./ampsci -o
  • See available input options for each operator by following with its name
  • e.g., ./ampsci -o hfs

Looking for atomic ionisation form-factors for dark-matter-electron scattering?

• See Kionisation module (./ampsci -m Kionisation)

Documentation

Full documentation available online: ampsci.dev/. Divided into sections:

1. Usage intructions, input options

   • Compilation instructions (for linux/mac/windows): doc/compilation.md

2. Physics documentation: ampsci.dev/ampsci.pdf

   • Description of physics/methods used in the code
   • Includes many references to the works where the methods implemented here were developed.

3. Modules

   • The modules system allows the easy calculation of any atomic properties after the wavefunction has been calculated. See doc/modules.md for description
   • The code is designed so that you can easily create your own modules. See doc/writing_modules.md for details

4. Code documentation -- details on classes/functions in the code

   • Available online: ampsci.dev/
   • This should only be required if you plan to edit the code or add new modules

---

A selection of publications resulting from ampsci

• Accurate electron-recoil ionization factors for dark matter direct detection in xenon, krypton and argon, A. R. Caddell, V. V. Flambaum, B. M. Roberts, Physical Review D 108, 083030 (2023), [arXiv:2305.05125]
• The neutrino fog for dark matter-electron scattering experiments, B. Carew, A. R. Caddell, T. N. Maity, C. A. J. O'Hare, arXiv:2312.04303
• Electric dipole transition amplitudes for atoms and ions with one valence electron, B. M. Roberts, C. J. Fairhall, J. S. M. Ginges, Physical Review A 107, 052812 (2023), [arXiv:2211.11134]
• Experimental and Theoretical Study of Dynamic Polarizabilities in the $5S - 5D_{5/2}$ Clock Transition in Rubidium-87 and Determination of Electric Dipole Matrix Elements, R. Hamilton et al., Physical Review Applied 19, 054059 (2023), [arXiv:2212.10743]
• QED radiative corrections to electric dipole amplitudes in heavy atoms, C. J. Fairhall, B. M. Roberts, J. S. M. Ginges, Physical Review A 107, 022813 (2023), [arXiv:2212.11490]

---

ampsci's Issues

Generalised input options for RPA (+SR)

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

Breit, s2

Possibly related to #26

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

210117-cs-btest.txt

TDHF for even operators

• 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)

B-spline stability

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

Breit/Sigma memory issue

• Including Breit into correlations uses absurd memory
• Maybe can be made more efficient?
• At the very least, allow a "basis string" input to limit number Breit integrals required to be stored

TDHF for frequency-dependent operators

• 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$)

E1(L)	$\omega$	HF	TDHF	Diagram	Basis
6p- - 6s+	0.0418	-5.2777	-4.9747	-4.9747	-4.9747
6p+ - 6s+	0.0436	7.4264	7.0137	7.0137	7.0137
6s+ - 6p-	-0.0418	-5.2777	-4.9747	-4.9747	-4.9747
6s+ - 6p+	-0.0436	-7.4264	-7.0137	-7.0137	-7.0137

Velocity form ($\omega$ dependent)

E1(v)	$\omega$	HF	TDHF	Diagram	Basis
6p- - 6s+	0.0418	-5.0371	1.9758	-4.9747	1.9755
6p+ - 6s+	0.0436	7.0662	-1.9837	7.0137	-1.9834
6s+ - 6p-	-0.0418	5.0371	12.0498	5.0994	12.0495
6s+ - 6p+	-0.0436	7.0662	16.1161	7.1187	16.1157

Full code input/output: E1_LvsV.txt

Matrix elements, frequency dependent, omega = each

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;
}

Correlations

• 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

Breit into Green's function

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}$.

Add self-energy vertex and magnetic loop Uehling (VP vertex)

DiracOperator

Cleaner way to construct general operators - make consistent

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

Catch2 error with g++ on mac

With g++ on macOS, the "mayfail" MixedStates test causes a SIGABRT failure (output below).
Doesn't impact actual code - but limits usefulness of tests.

It works fine macOS with clang, and on linux with g++ and clang.
Fails both with and without openMP.
It happens on a catch2 'REQUIRE' clause that fails - this test is known to fail, and the program should mark it as 'failed as expected' and continue. Instead, the Catch::TestFailureException exception is thrown.

I think this is a maybe bug with Catch2?
Not sure best workaround?

Output:

AMPSCI v: 0.0 [main/0e3a2d69]
Libraries:
  GSL (GNU Scientific Libraries): 2.7.1
  OpenMP: 201511
Compiled: g++-13 [Homebrew GCC 13.2.0] 13.2.0 2024-05-23 15:20 AEST

./tests \[MixedStates\]

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
tests is a Catch v2.13.9 host application.
Run with -? for options

-------------------------------------------------------------------------------
External Field: Mixed-states2
-------------------------------------------------------------------------------
src/ExternalField/MixedStates.tests.cpp:87
...............................................................................

src/ExternalField/MixedStates.tests.cpp:209: FAILED:
  REQUIRE( worst < 1.0e-2 )
with expansion:
  0.1529438228 < 0.01

terminate called after throwing an instance of 'Catch::TestFailureException'
src/ExternalField/MixedStates.tests.cpp:209: FAILED:
  {Unknown expression after the reported line}
due to a fatal error condition:
  SIGABRT - Abort (abnormal termination) signal

===============================================================================
test cases:  2 | 1 passed | 1 failed
assertions: 10 | 8 passed | 2 failed as expected

Expected output

AMPSCI v: 0.0 [main/0e3a2d69]
Libraries:
  GSL (GNU Scientific Libraries): 2.7.1
  OpenMP: 201511
Compiled: g++-13 [Ubuntu 13.1.0-8ubuntu1~22.04] 13.1.0 2024-05-23 15:20 AEST

./tests [MixedStates]

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
tests is a Catch v2.13.9 host application.
Run with -? for options

-------------------------------------------------------------------------------
External Field: Mixed-states2
-------------------------------------------------------------------------------
src/ExternalField/MixedStates.tests.cpp:87
...............................................................................

src/ExternalField/MixedStates.tests.cpp:209: FAILED:
  REQUIRE( worst < 1.0e-2 )
with expansion:
  0.1529438228 < 0.01

===============================================================================
test cases: 2 | 1 passed | 1 failed as expected
assertions: 9 | 8 passed | 1 failed as expected

Feynman

• 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?

Add third-order matrix element method

e.g., (8.26)
Particularly for M1

Coulomb and Angular Tables

• Update Coulomb tables, angular tables

Add l-dependence to Radiative Potential

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.

Add ".amp" or similar to all binary outputs to help cleanup

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

Polarisation operator for low core n

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.

B-splines

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?

Spectrum stability with Feynman Sigma

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

Hartree-Fock update: in-homogeneous method

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

Hartree Fock + Wavefunction

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;

DiracODE

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

Breit with RPA [probably no issue]

For hyperfine: possibly linked to diagram method?

TDHF is slow - parallelised poorly

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

Spuriously combines modules input blocks

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;}

StrucRad to use full RPA methods

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