 |
|
High-precision calculations for one- and two-valence atomic systems
|
|
|
High-precision calculations for one- and two-valence atomic systems
|
CI+MBPT calculations for two-valence atoms.
This assumes you already have ampsci compiled and have a basic understanding of how to run and use it.
For an \(M\)-valence atomic system, the effective Hamiltonian is
\[ H_{\rm CI} = H_{\rm CI}^1 + H_{\rm CI}^2 = \sum_i^M \left(h^{\rm HF}(r_i) + \Sigma^1(r_i)\right) +\sum_{i < j}\left(r^{-1}_{ij} +\Sigma^2(r_i,r_j)\right), \]
where \(h^{\rm HF}\) is one-particle the Hartree-Fock Hamiltonian (with HF potential due to the \(N-M\) core electrons), \(\Sigma^1\) accounts for the core-valence correlations, and \(\Sigma^2\) accounts for the screening of the valence-valence Coulomb interaction by the core electrons.
The CI routines in ampsci are only for two-valence systems: \(M=2\).
In the CI method, approximate valence-space wavefunctions, \(\Psi\), are expanded over \(M\)-particle wavefunctions called Configuration-State Functions (CSFs), \(\psi_I\):
\[ \left|{\Psi,J^\pi J_z}\right\rangle = \sum_I c_I \left|{I,J^\pi J_z}\right\rangle. \]
The CSFs are combinations of Slater-determinants formed from single-particle eigenfunctions. The CSFs are eigenfunctions of \(J^2\), \(J_z\), and parity ( \(\pi\)). For each \(J^\pi\) symmetry, the energies and wavefunctions (expansion coefficients) are found by solving the Schr"odinger equation, which for a finite set of \(N_{CSF}\) CSFs, is cast to an \(N_{CSF}^2\) eigenvalue problem:
\[ \sum_J c_J\left\langle{I}\right|H_{\rm eff}\left|{J}\right\rangle = E c_I, \]
For the single-particle basis, we use eigenfunctions of the same \(h^{\rm HF}\) Hamiltonian from the CI Hamiltonian. This is known as the \(V^{N-M}\) approximation, with simplifies the MBPT part of the calculation, and is very accurate for two-valence systems.
As always, we check the available input options: ./ampsci -i CI
./ampsci -i CI
For a simple example, we'll consider neutral Mg.
The first part of the input file will be familiar from previous examples. We don't need valence states in Hartree-Fock, so, this can be left blank. We don't include the two valence states in the core, hence use \(V^{N-2}\) approximation.
To begin the CI calculation
The ci_basis options means we will use two-particle CSFs formed from combinations of single-particle basis states up to \(n=20\) for \(s\), \(p\), \(d\), and \(f\) states. We will find the lowest 3 solutions for each \(J^\pi\) symmetry up to \(J=2\).
The first part of the output:
The code pre-computes all the two-body Coulomb integrals \(Q^k_{ijkl}\). This uses a significant memory footprint, but makes the subsequent calculations much faster. It only calculates the required integrals. Even though our total basis was up to 30spdfghi, we only used 20spdf in the CI expansion (the remaining basis states will be used for MBPT in the next example).
Then, the code runs CI routine for each symmetry. The output for each symmetry will look something like this:
After each line, the code will list all CSFs that contribute to the solution at above the 1% level (i.e., with \(c_I^2>0.01\)). The code also calculates the \(g\)-factor (without RPAd), which are useful for level identification. The term symbol ( \({}^{2S+1}L_J\) e.g., 1^D_2) and leading configuration are given – the term symbol is determined from the expectation values \(\langle L^2\rangle\) and \(\langle S^2\rangle\) (evaluated in the non-relativistic limit, and printed as \(L_{\rm eff}\) and \(S_{\rm eff}\), where \(L_{\rm eff}(L_{\rm eff}+1)=\langle L^2\rangle\)). It is not well defined relativistically, so should be seen as indicative only.
The final output shows a summary of the calculation:
Where conf is the leading configuration (in non-relativistic notation), and the % column shows the combined constructions from each relativistic configuration with the same non-relativistic configuration (e.g., 3s3d may have contributions from 3s3d- and 3s3d+).
| Level | AMPSCI | Exp. | \(\Delta\) | |
|---|---|---|---|---|
| \(3s^2\) | \({}^1S_0\) | -179533 | -182939 | -1.9% |
| \(3s4s\) | \({}^1S_0\) | 42665 | 43503 | -1.9% |
| \(3s6s\) | \({}^1S_0\) | 51644 | 52556 | -1.7% |
| \(3s4s\) | \({}^3S_1\) | 40405 | 41197 | -1.9% |
| \(3s3d\) | \({}^3D_1\) | 46971 | 47957 | -2.1% |
| \(3s5s\) | \({}^3S_1\) | 50974 | 51873 | -1.7% |
| \(3s3d\) | \({}^1D_2\) | 45119 | 46403 | -2.8% |
| \(3s3d\) | \({}^3D_2\) | 46971 | 47957 | -2.1% |
| \(3s4d\) | \({}^1D_2\) | 52038 | 53135 | -2.1% |
| \(3s3p\) | \({}^3P^o_0\) | 20905 | 21850 | -4.3% |
| \(3s4p\) | \({}^3P^o_0\) | 46917 | 47841 | -1.9% |
| \(3s6p\) | \({}^3P^o_0\) | 53312 | 54249 | -1.7% |
| \(3s3p\) | \({}^3P^o_1\) | 20925 | 21870 | -4.3% |
| \(3s3p\) | \({}^1P^o_1\) | 34491 | 35051 | -1.6% |
| \(3s4p\) | \({}^3P^o_1\) | 46921 | 47844 | -1.9% |
| \(3s3p\) | \({}^3P^o_2\) | 20966 | 21911 | -4.3% |
| \(3s4p\) | \({}^3P^o_2\) | 46927 | 47851 | -1.9% |
| \(3s6p\) | \({}^3P^o_2\) | 53316 | 54253 | -1.7% |
The table shows the results of the calculation, and comparison to experimental excitation energies, in units of \({\rm cm}^{-1}\) (for the ground state, the ionisation potential is instead shown). The agreement is at the ~few % level.
CI solutions (energies and eigenvectors) can be saved to disk and read back on subsequent runs, which avoids re-diagonalising the Hamiltonian. By default, a filename is generated automatically based on the atomic identity, CI basis, and a hash of the other settings that change the CI Hamiltonian (e.g., Mg2_20spdf_u05ndb.ci.abf). This can be overridden with the ci_file option; set it to false to disable read/write entirely.
On the first run the solutions are computed and written to file. On subsequent runs ampsci reads them back and skips diagonalisation. If more solutions are requested than are stored, ampsci re-solves (from scratch) and overwrites those objects in the file.
If the CI solutions already exist on disk and you do not need to recalculate (e.g., you are only computing matrix elements), set read_only = true. This skips all integral calculations and reads directly from the file:
The same CI output file stores solutions for each J/parity independently. For example, if you required 20 J=1- levels, but only a few from the other states, you could run twice like:
J=1- solutions will be read in; there are no existing solutions for the other states, so the requested 3 will be calculated.Here, we will improve the accuracy of the calculation by including the MBPT corrections. We do this by setting sigma1 and sigma2 options to true.
This will use the full basis (from Basis{}) for the MBPT part. In this case, it was 30spdfghi. (Different subsets of the full basis can be used for the internal lines of the \(\Sigma^1\) and \(\Sigma^2\) diagrams using the s1_basis and s2_basis options, respectively).
The code will now calculate more \(Q^k\) Coulomb integrals, since they are required in the MBPT calculations. It will read in the existing file (if there is one), so it doesn't need to start from scratch, and only calculate the missing integrals:
This uses close to 5Gb of memory, so may already be difficulty on some laptops. Then, is will calculate the matrix elements of the two-body \(\Sigma^2\) operator. This is the slow part of the calculation. Fortunately, the \(\Sigma^2\) correction is rather small, and good accuracy can be obtained by only including matrix elements between the lowest few valence-space basis states. By default, it will include up to \(n=n_{\rm core}+3\), where \(n_{\rm core}\) is the largest \(n\) in the core, for each \(l\) in the ci_basis. This can be controlled manually with the cis2_basis option.
The For: 5spdf, using 30spdfghi means the external lines in \(\Sigma_{ijkl}\) include up to 5spdf, while the internal lines use the full 30spdfghi basis. This is usually the bottle-neck; these integrals take a very long time to compute, though typically make only a modest difference to the calculations.
The combined MBPT corrections lead to a significant improvement in the accuracy:
| Level | AMPSCI | Exp. | \(\Delta\) | |
|---|---|---|---|---|
| \(3s^2\) | \({}^1S_0\) | -182703 | -182939 | -0.13% |
| \(3s4s\) | \({}^1S_0\) | 43440 | 43503 | -0.14% |
| \(3s6s\) | \({}^1S_0\) | 52473 | 52556 | -0.16% |
| \(3s4s\) | \({}^3S_1\) | 41146 | 41197 | -0.12% |
| \(3s3d\) | \({}^3D_1\) | 47895 | 47957 | -0.13% |
| \(3s5s\) | \({}^3S_1\) | 51810 | 51873 | -0.12% |
| \(3s3d\) | \({}^1D_2\) | 46326 | 46403 | -0.17% |
| \(3s3d\) | \({}^3D_2\) | 47895 | 47957 | -0.13% |
| \(3s4d\) | \({}^1D_2\) | 53052 | 53135 | -0.16% |
| \(3s3p\) | \({}^3P^o_0\) | 21787 | 21850 | -0.29% |
| \(3s4p\) | \({}^3P^o_0\) | 47777 | 47841 | -0.13% |
| \(3s6p\) | \({}^3P^o_0\) | 54204 | 54249 | -0.08% |
| \(3s3p\) | \({}^3P^o_1\) | 21809 | 21870 | -0.28% |
| \(3s3p\) | \({}^1P^o_1\) | 35041 | 35051 | -0.03% |
| \(3s4p\) | \({}^3P^o_1\) | 47781 | 47844 | -0.13% |
| \(3s3p\) | \({}^3P^o_2\) | 21852 | 21911 | -0.27% |
| \(3s4p\) | \({}^3P^o_2\) | 47788 | 47851 | -0.13% |
| \(3s6p\) | \({}^3P^o_2\) | 54208 | 54253 | -0.08% |
The discrepancies are now at the level of ~0.1%.
On my PC, this entire calculation took less than two minutes.
For two-valence CI, the memory bottleneck is calculating the two-body Coulomb ( \(Q^k\)) integrals; the time bottleneck is usually calculating the two-body Correlation.
You can estimate the memory required using the ampsci -z command-line option. In the above example, we used:
20spdf for CI30spdfghi for MBPT, with a multipolarity cut-off of 7 (default)The CI calculation requires only excited-excited Coulomb integrals, while the MBPT only requires Coulomb integrals involving 1 or 2 core states.
We can therefore check the memory requirements by running:
Should return something like:
and
Combining, we see that we should require 15896631 for the CI part, and 29763502 for MBPT. This exactly matched the calculated Qk integral summaries we saw above. In this case, the total memory requirement is very small (~2 Gb). Notice, however, if we had tried to run CI with the full MBPT basis, we would have required nearly 50Gb!
This takes very similar options to the regular MatrixElements{} module:
It's typically not recommended to use omega=each, as this solves the RPA equations for each transition, and the frequency-dependence is usually small.
The output is something like:
(this is for a reduced set of levels).
This implies the lifetime of the first \({}^1P^o_1\) level is 2.12 ns, in excellent agreement with experiment.
Second-order (in the external field) amplitudes between CI states – polarisabilities (scalar, tensor, vector/transition) and PNC amplitudes – are calculated with the CI_secondOrder module. The sums over intermediate states are evaluated with CI mixed states, so they are complete (no sum over individual CI solutions, and no truncation of the spectrum); each sum is formed two independent ways as an internal check.
For two one-body operators \( t \) (rank \( k_t \), frequency \( \omega \)) and \( s \) (rank \( k_s \), frequency \( \omega_s \)), the module reports the reduced amplitude of rank \( K \),
\[ A^K = \sum_n \left[ c_1 \frac{\langle B||t||n\rangle\langle n||s||A\rangle}{E_A+\omega_s-E_n} + c_2 \frac{\langle B||s||n\rangle\langle n||t||A\rangle}{E_A+\omega-E_n} \right], \]
and the specific cases ( \( [J]\equiv 2J+1 \)):
\[ \alpha_0 = \frac{A^0}{\sqrt{3[J_A]}}, \qquad \alpha_2 = -\sqrt{\frac{2J(2J-1)}{3(J+1)(2J+1)(2J+3)}}\;A^2, \qquad \beta = \frac{A^1}{\sqrt{2}\,\langle B||\boldsymbol\sigma||A\rangle}, \]
for \( t=s=E1 \) and \( K = 0, 2, 1 \), and \( E_{\rm PNC} = A^1_0 \) for \( t=E1 \), \( s=h_W \). The definitions of \( A^{k_tk_s}_{q_1q_2} \), \( A^K_Q \) and \( A^K \), and the sign convention of the coupling, are on the Second-Order Amplitudes page (identical, with \( j \to J \)) and in CI::A_K_coefs().
For example, the static scalar polarisability of the ground state:
States are given as J{+,-}:index (so 0+ is the lowest even J=0 level, and 1-:1 the second-lowest odd J=1 level). Run ./ampsci -m CI_secondOrder for the full list of options; see CI::A_K() for the formulas and conventions. The single-valence analogue is the Second-Order Amplitudes module.
1.9.8