ampsci
c++ program for highprecision atomic structure calculations of singlevalence systems

brief Advanced ampsci tutorial: correlations
This assumes you already have ampsci compiled and have a basic understanding of how to run and use it.
The full atomic Hamiltonian for an \(N\)electron atom:
\[ H = \sum_i^N h_0(r_i) + \sum_{i < j}\frac{1}{r_i  r_j} \]
\[ h_0(\vec{r_i}) = c \vec\alpha_i\cdot\vec{p_i} + c^2 (\beta_i1) + V_{\rm nuc}. \]
The starting approximation is the HartreeFock method:
\[ H \approx \sum_i^N [h_0(\vec r_i) + v^{\rm HF} ]. \]
We focus on case of singlevalence systems, and start with socalled \(V^{N1}\) approximation, in which HartreeFock potential is due to the \(N1\) core electrons:
\[ \hat v^{\rm HF}\phi_a(\vec{r_1}) = \sum_{i\neq a}^{N_c}\Bigg( \int \frac{\phi_i^\dagger(\vec{r_2})\phi_i(\vec{r_2})}{r_{12}}d^3\vec{r_2}\,\phi_a(\vec{r_1}) \int \frac{\phi_i^\dagger(\vec{r_2})\phi_a(\vec{r_2})}{r_{12}}d^3\vec{r_2}\,\phi_i(\vec{r_1}) \Bigg), \]
First, HartreeFock equations are solved selfconsistently for all core electrons, then the valence states are found in the Frozen HartreeFock potential due to the core. For example, to run HartreeFock for Cs with a \(N1\) (implies Xelike core),
Besides HartreeFock, other methods are available:
e.g.,
If no innputs are given, the code will assume a Fermilike distribution for the nuclear charge, and look up default values for the nuclear parameters (charge radius and skin thickness). The default values are chosen from the specified isotope in the Atom{}
block. These may also be specified specifically:
Default rms values are taken from:
It's typically best to leave these as the default parameters, but you're free to update them, e.g., the following will assume a sphericallysymmetric nucleus with rms charge radius of 3.5 fm:
Nucleus{ rrms = 3.5; type = spherical; }
The Breit Hamiltonian accounts for magnetic interactions between electrons (also known as the Gaunt interaction), and retardation effects. It leads to a correction to the electronelectron Coulomb term in the manybody Hamiltonian:
\[ \sum_{ij}\frac{1}{r_{ij}} \to \sum_{ij}\left( \frac{1}{r_{ij}} + \hat h^B_{ij}\right), \]
where, in the limit of zero frequency, the twoparticle Breit Hamiltonian is
\[ h^B_{ij} =  \frac{\vec{\alpha_i}\cdot\vec{\alpha_j} + (\vec{\alpha_i}\cdot\hat{n_{ij}})(\vec{\alpha_j}\cdot\hat{n_{ij}})}{2\, r_{ij}}. \]
This can be included at the HartreeFock level with the HartreeFock{Breit = 1.0;}
setting.
This setting is a scaling factor; i.e., setting Breit = 0.5;
will add effective factor of 0.5 in front of \(h^B\). Typically, only 0 or 1 is set; other values are useful for checking for nonphysical nonlinearinBreit effects.
Note: it's very important when including Breit effects along with correlations that the Breit Hamiltonian is included in the Dirac equation when forming the basis used to construct the correlation potential! This is to ensure basis is orthogonal to the HartreeFock core states At the moment, there is no way to include the Breit effect into the HartreeFock Green's function; as a result, we cannot include Breit at the level of allorder correlations (see below).
Radiative QED corrections can be included into the wavefunctions using the FlambaumGinges radiative potential method. An effective potential, \(V_{\rm rad}\), is added to the Hamiltonian before the equations are solved. The potential can be written as the sum of the Uehling (vacuum polarisation), WichmannKroll (higherorder vacuum polarisation) and selfenergy potentials; the selfenergy potential itself is written as the sum of the high and lowfrequency electric contributions, and the magnetic contribution:
\[ V_{\rm rad}(\vec{r}) = V_{\rm Ueh}(r) + V_{\rm WK}(r) + V_{\rm SE}^{h}(r) + V_{\rm SE}^{l}(r) + i (\vec{\gamma}\cdot\hat{n}) V^{\rm mag}(r). \]
If the RadPot{}
block is included, then the QED radiative potential will be included with the default parameters.
equivalent to:
Ueh, SE_h, SE_l, SE_m, WK
) options are scaling factors for their corresponding terms in the potential; typically 0 or 1, but can be tuned (for testing).core_qed
means radiative potential will be included into the core states  this gives an important contribution and should only be set to false for testing.scale_l
takes a list of scaling factors for each l; these will rescale \(V_{\rm rad}\) for each partial wave. e.g., scale_l = 0,1,0;
will include \(V_{\rm rad}\) for p states, but not s or d states.scale_rN;
Rescales the effective nuclear radius; used to test finitenuclear size effects on \(V_{\rm rad}\). =1 means normal, =0 means assume pointlike nucleus (when calculating \(V_{\rm rad}\)).We often require a full "compete" set of solutions to the HartreeFock equation in order to perform sumsoverstates required in perturbation theory. To form these, the set of atomic orbitals are expanded as
\[ \phi_{n\kappa}(\vec{r}) = \sum_i^{2N} p_i S_i(\vec{r}), \]
where \(S_i\) are a set of \(2N\) basis orbitals that form an approximately complete set over a subdomain of the radial grid \([0,r_{\rm max}]\) ( \(N\) is defined this way because of the duel set of positive/negative energy Dirac solutions). The \(p_i\) expansion coefficients are found by diagonalising the set of basis orbitals with respect to the Hamiltonian matrix, equivalent to solving the eigenvalue problem:
\[ \langle{S_i}h_{\rm HF} S_j \rangle p_i = \varepsilon\langle S_iS_j \rangle p_i. \]
There are \(2N\) solutions of eigenvalues \(\varepsilon\) with corresponding eigenvectors \(\vec{p}\), which correspond to the spectrum of stationary states; \(N\) of these correspond to negativeenergy ( \(\varepsilon<mc^2\)) states.
For the \(S_i\) basis orbitals, we use the DuelKineticBalence basis of Beloy and Dereviano, and \(S_i\) is built from \(N_{\rm spl}\) Bsplines of order \(k\).
r0
is the location of the first 'internal' knot (first actual knot is always placed at r=0). Typically, we instead set r0_eps
, which automatically choses r0
such that the core density \(\rho_l(r) = \sum_n\phi_{nl}(r)^2\) drops below given relative value \(\rho_l(r0)/\rho_{\rm max}<{\rm r0\_{\rm eps}}\).
e.g, to calculate basis including up to \(n=30\) for states up to \(l=6\) in a cavity of 40 \(a_0\), using 40 Bsplines of order 7:
If Breit and/or QED is set (see above), these are included into Hamiltonian and thus into the basis also.
It's important to look at the output when generating the basis. In the above case (with all other default options for Cs), we get output like:
The line Spline cavity l=0 s: (1.4e04, 40.0)aB.
tells us the actual point of the first (internal) and final Bspline knot.
The most important output is the Basis/core:
output  this tells us how orthogonal the generated basis is to the HartreeFock core. The MBPT formalism relies on the orthogonality here, so it's important to check. If the results are not good enough (the code will warn you with **
), try increasing rmax
and/or number
.
In this case, the worst normality occurred for the 3s state, where the innerproduct of the finitedifference (Hartree Fock) \(3s\) state and the corresponding basis \(3s\) state was different from 1 by parts in \(10^6\). The worst energy comparison was for \(5s\) state (parts in \(10^5\)), and the worst orthogonality comparison was for the \(3s\) finitedifference HartreeFock state, and the \(22s\) basis state, which were orthogonal to parts in \(10^4\).
The spectrum is the same as the basis (takes the same options), except that it also includes correlations (see below):
\[ \langle{S_i}\hat h_{\rm HF} + \hat \Sigma{S_j}\rangle p_i = \varepsilon\langle{S_iS_j}\rangle p_i. \]
We typically use Basis to calculate \(\Sigma\) (or other MBPT corrections), and then use spectrum to calculate atomic properties where a direct sumoverstates is required (e.g., polarisabilities). Since we often only use basis to calculate MBPT corrections, it doesn't need to be very large. We directly use the spectrum, so we typically require a larger basis set. It's now typically also crucial for the spectrum to be orthogonal to the valence states (not just the core states). e.g.,
\[ H = \sum_i h_{\rm HF}(\vec{r_i}) + \delta V_{\rm corr}, \]
where \(h_{\rm HF}(\vec{r}_i)\) is the singleparticle HF Hamiltonian, and
\[ \delta V_{\rm corr} = \sum_{i < j}\frac{1}{r_{ij}}  \sum_i V_{\rm HF}(\vec{r_i}) \]
\[ \delta \varepsilon_v = \sum_{amn} \frac{g_{vamn}\widetilde g_{nmav}}{\varepsilon_v+\varepsilon_a  \varepsilon_m\varepsilon_n} +\sum_{abn} \frac{g_{vnab}\widetilde g_{banv}}{\varepsilon_v+\varepsilon_n\varepsilon_a\varepsilon_b} , \]
We define correlation potential, \(\Sigma\):
\[ \delta \varepsilon_v =\langle v  \Sigma v \rangle \]
We then solve the HartreeFock equation, including the correlation potential. The solutions are known as Brueckner orbitals, and include correlation corrections.
\[ [h_{\rm HF} + \Sigma(\varepsilon)]\phi^{\rm Br} = \varepsilon^{\rm Br}\phi^{\rm Br} \]
There are quite a few options available for Correlation corrections (see ampsci a Correlations
) – we will focus on the most important here.
The read
and write
option are filenames the correlation potential will be written to and/or read from (this saves time). If none are given, it will read/write from the default filename. Set to 'false' to not read/write.
n_min_core
is minimum n for core states to include in polarisation loops. Very low core states contribute very little (due to large excitation energy), so this saves much time for almost no degradation in results.
If each_valence
is true, a new correlation potential will be calculated for each valence state, fully taking energy dependence into account; if false, only 1 \(\Sigma\) for \(\kappa\) will be formed (i.e. assumes the same energy for calculating \(\Sigma\) for each angular symmetry).
fitTo_cm
takes a list of experimental energies (in \({\rm cm}^{1}\)). A scaling factor \(\lambda\) will be introduced in front of the correlation potential, and will be tuned so that energies exactly match those given. This is a semiempirical method for accounting for higherorder correlations in wavefunctions. lambda_kappa
is for manually setting the scaling factors. Note that for both of these, the input list must be in the exact same order as the valence states  check the 'valence' output.
Will form the correlation potential \(\Sigma(\varepsilon)\) for the listed states at the given energy \(\varepsilon\).
Will use the correlation potential for \(\kappa=1\) and \(\varepsilon = 0.127\) for the \(6s\) state, and \(\kappa=+1\) and \(\varepsilon = 0.127\) for the \(6p_{1/2}\) state.
The most important corrections beyond the secondorder correlation potential are the screening of the residual Coulomb interaction by the core electrons and holeparticle interaction. These are taken into account using the allorders correlation potential method, also called Perturbation Theory in the Screened Coulomb Interaction (PTSCI), developed by Dzuba, Flambaum and Sushkov [13].
The starting point uses the Feynman technique, in which the direct part of the correlation potential can be expressed
\[ \Sigma_{\rm d} = \int\frac{{\rm d}\omega}{2\pi} G_{12}(\varepsilon+\omega) Q_{1i}\Pi_{ij}(\omega) Q_{j2}(\omega), \]
where \(G_{12} = G(r_1,r_2)\) in the HarteeFock Feynman Green's function, \(Q\) is the (nonrelativistic) Coulomb operator, and \(\Pi\) is the polarisation operator (subscripts are coordinate indices; integration is assumed over internal \(i\) and \(j\)). Note that \(\Pi\) represents polarisation of the atomic core.
The screening of the coulomb interaction can be represented by a series of polarisation corrections
\[ \widetilde Q \equiv Q + \ Q(i\, \Pi Q) + Q(i\, \Pi Q)^2+\ldots \]
which can be summed exactly:
\[ \widetilde Q(\omega) = Q\left[1+i\,\Pi(\omega) Q\right]^{1}. \]
The screening is accounted for in the direct diagrams via \(Q\to \widetilde Q\) for one of the Coulomb operators in \(\Sigma_{\rm d}\). (Screening for exchange diagrams is taken into account in a simpler fasion, see pdf for details).
The holeparticle interaction arises due to the deviation of the HartreeFock potential for the excited core electron in the polarisation loop from that for the nonexcited one. The potential that simultaneously describes the occupied core and excited states is
\[ \hat V = V^{N1}  (1\hat P_{\rm core})V_{\rm self} (1\hat P_{\rm core}), \]
where \(P_{\rm core}\) is the operator of projection onto the core, and \(V_{\rm self}\) is the selfinteraction part of the HartreeFock potential for the outgoing electron. Therefore, holeparticle interaction is accounted for by using this potential when forming the polarisation operator.
To include allorders correlations, the Feynman
, screening
, and holeParticle
options in the Correlations{}
block should be set to true
. Alternativel, just set AllOrder=true;
equivilant to:
More options are available, though they rarely need to be changed from the default. See ampsci.pdf for full details.
Note: currently, there is a numerical problem in calculating polarisation loop for very deep core states. This leads to no issues, unless you attempt to polarise a deep core shell, where large numerical errors are encountered. Best to set n_min_core = 2
for Rb, n_min_core = 3
for Cs, n_min_core = 4
for Fr etc. The effect of this is negligable, as can be checked by seeing impact on secondorder results. This should be addressed soon.
Note also: it is currently not possible to correctly account for Breit Hamiltonian within the allorders method (problem is in construction of Green's function). Therefore, Breit should only be included at secondorder level.
Core polarisation (RPA) is included in the matrix elements.
The best method to use is TDHF, which is numerically stable, and includes contribution from negative energy states automatically. However, it doesn't work (at the moment) for even parity operators, so for these we must use the diagram method. The diagram method uses Basis{}
for internal summations, so one should ensure the basis is complete enough for this to work.
To improve tha accuracy of matrix elements, structure radiation and normalisation corrections should be included. There is an option to do this in the MatrixElements module. There is also a ModulebriefStructureRad
module, which gives some finer control.
If a Qk filename is given, program will first calculate all required Q^k Coulomb integrals before calculating structure radiation. This speeds up the calculation, at a great memory cost.