|
ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
|
brief Advanced ampsci tutorial: CI+MBPT 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 -a CI
./ampsci -a 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 'guessed' based on the g-factor. 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\) | -179539 | -182939 | -1.9% |
| \(3s4s\) | \({}^1S_0\) | 42665 | 43503 | -1.9% |
| \(3s6s\) | \({}^1S_0\) | 51644 | 52556 | -1.7% |
| \(3s4s\) | \({}^3S_1\) | 40403 | 41197 | -1.9% |
| \(3s3d\) | \({}^3D_1\) | 46971 | 47957 | -2.1% |
| \(3s5s\) | \({}^3S_1\) | 50973 | 51873 | -1.7% |
| \(3s3d\) | \({}^1D_2\) | 45121 | 46403 | -2.8% |
| \(3s3d\) | \({}^3D_2\) | 46971 | 47957 | -2.1% |
| \(3s4d\) | \({}^1D_2\) | 52038 | 53135 | -2.1% |
| \(3s3p\) | \({}^3P^o_0\) | 20908 | 21850 | -4.3% |
| \(3s4p\) | \({}^3P^o_0\) | 46917 | 47841 | -1.9% |
| \(3s6p\) | \({}^3P^o_0\) | 53311 | 54249 | -1.7% |
| \(3s3p\) | \({}^3P^o_1\) | 20928 | 21870 | -4.3% |
| \(3s3p\) | \({}^1P^o_1\) | 34491 | 35051 | -1.6% |
| \(3s4p\) | \({}^3P^o_1\) | 46920 | 47844 | -1.9% |
| \(3s3p\) | \({}^3P^o_2\) | 20969 | 21911 | -4.3% |
| \(3s4p\) | \({}^3P^o_2\) | 46927 | 47851 | -1.9% |
| \(3s6p\) | \({}^3P^o_2\) | 53315 | 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.
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 leads to a significant improvement in the accuracy:
| Level | AMPSCI | Exp. | \(\Delta\) | |
|---|---|---|---|---|
| \(3s^2\) | \({}^1S_0\) | -182804 | -182939 | -0.07% |
| \(3s4s\) | \({}^1S_0\) | 43490 | 43503 | -0.03% |
| \(3s6s\) | \({}^1S_0\) | 52526 | 52556 | -0.06% |
| \(3s4s\) | \({}^3S_1\) | 41208 | 41197 | 0.03% |
| \(3s3d\) | \({}^3D_1\) | 47957 | 47957 | 0.00% |
| \(3s5s\) | \({}^3S_1\) | 51875 | 51873 | 0.00% |
| \(3s3d\) | \({}^1D_2\) | 46385 | 46403 | -0.04% |
| \(3s3d\) | \({}^3D_2\) | 47960 | 47957 | 0.01% |
| \(3s4d\) | \({}^1D_2\) | 53113 | 53135 | -0.04% |
| \(3s3p\) | \({}^3P^o_0\) | 21822 | 21850 | -0.13% |
| \(3s4p\) | \({}^3P^o_0\) | 47837 | 47841 | -0.01% |
| \(3s6p\) | \({}^3P^o_0\) | 54271 | 54249 | 0.04% |
| \(3s3p\) | \({}^3P^o_1\) | 21844 | 21870 | -0.12% |
| \(3s3p\) | \({}^1P^o_1\) | 35090 | 35051 | 0.11% |
| \(3s4p\) | \({}^3P^o_1\) | 47840 | 47844 | -0.01% |
| \(3s3p\) | \({}^3P^o_2\) | 21876 | 21911 | -0.16% |
| \(3s4p\) | \({}^3P^o_2\) | 47845 | 47851 | -0.01% |
| \(3s6p\) | \({}^3P^o_2\) | 54272 | 54253 | 0.03% |
The discrepancies are now at the level of 0.1% or below.
On my pc, this entire calculation took less than two minutes.
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.
1.9.1