|
ampsci
c++ program for high-precision atomic structure calculations of single-valence 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_i-1) + V_{\rm nuc}. \]
The starting approximation is the Hartree-Fock method:
\[ H \approx \sum_i^N [h_0(\vec r_i) + v^{\rm HF} ]. \]
We focus on case of single-valence systems, and start with so-called \(V^{N-1}\) approximation, in which Hartree-Fock potential is due to the \(N-1\) 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, Hartree-Fock equations are solved self-consistently for all core electrons, then the valence states are found in the Frozen Hartree-Fock potential due to the core. For example, to run Hartree-Fock for Cs with a \(N-1\) (implies Xe-like core),
Besides Hartree-Fock, other methods are available:
e.g.,
If no innputs are given, the code will assume a Fermi-like 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 spherically-symmetric 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 electron-electron Coulomb term in the many-body 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 two-particle 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 Hartree-Fock 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 non-physical non-linear-in-Breit 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 Hartree-Fock core states At the moment, there is no way to include the Breit effect into the Hartree-Fock Green's function; as a result, we cannot include Breit at the level of all-order correlations (see below).
Radiative QED corrections can be included into the wavefunctions using the Flambaum-Ginges 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), Wichmann-Kroll (higher-order vacuum polarisation) and self-energy potentials; the self-energy potential itself is written as the sum of the high- and low-frequency 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; Re-scales the effective nuclear radius; used to test finite-nuclear size effects on \(V_{\rm rad}\). =1 means normal, =0 means assume point-like nucleus (when calculating \(V_{\rm rad}\)).We often require a full "compete" set of solutions to the Hartree-Fock equation in order to perform sums-over-states 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 sub-domain 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_i|S_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 negative-energy ( \(\varepsilon<-mc^2\)) states.
For the \(S_i\) basis orbitals, we use the Duel-Kinetic-Balence basis of Beloy and Dereviano, and \(S_i\) is built from \(N_{\rm spl}\) B-splines 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 B-splines 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.4e-04, 40.0)aB. tells us the actual point of the first (internal) and final B-spline knot.
The most important output is the Basis/core: output - this tells us how orthogonal the generated basis is to the Hartree-Fock 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 inner-product of the finite-difference (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\) finite-difference Hartree-Fock 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_i|S_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 sum-over-states 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 single-particle 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 Hartree-Fock 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 semi-empirical method for accounting for higher-order 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 second-order correlation potential are the screening of the residual Coulomb interaction by the core electrons and hole-particle interaction. These are taken into account using the all-orders correlation potential method, also called Perturbation Theory in the Screened Coulomb Interaction (PTSCI), developed by Dzuba, Flambaum and Sushkov [1-3].
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 Hartee-Fock Feynman Green's function, \(Q\) is the (non-relativistic) 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 hole-particle interaction arises due to the deviation of the Hartree-Fock potential for the excited core electron in the polarisation loop from that for the non-excited one. The potential that simultaneously describes the occupied core and excited states is
\[ \hat V = V^{N-1} - (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 self-interaction part of the Hartree-Fock potential for the outgoing electron. Therefore, hole-particle interaction is accounted for by using this potential when forming the polarisation operator.
To include all-orders 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 second-order results. This should be addressed soon.
Note also: it is currently not possible to correctly account for Breit Hamiltonian within the all-orders method (problem is in construction of Green's function). Therefore, Breit should only be included at second-order 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.
1.9.1