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

brief Advanced ampsci tutorial: CI+MBPT for twovalence 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 oneparticle the HartreeFock Hamiltonian (with HF potential due to the \(NM\) core electrons), \(\Sigma^1\) accounts for the corevalence correlations, and \(\Sigma^2\) accounts for the screening of the valencevalence Coulomb interaction by the core electrons.
The CI routines in ampsci are only for twovalence systems: \(M=2\).
In the CI method, approximate valencespace wavefunctions, \(\Psi\), are expanded over \(M\)particle wavefunctions called ConfigurationState 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 Slaterdeterminants formed from singleparticle 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}\rightH_{\rm eff}\left{J}\right\rangle = E c_I, \]
For the singleparticle basis, we use eigenfunctions of the same \(h^{\rm HF}\) Hamiltonian from the CI Hamiltonian. This is known as the \(V^{NM}\) approximation, with simplifies the MBPT part of the calculation, and is very accurate for twovalence 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 HartreeFock, so, this can be left blank. We don't include the two valence states in the core, hence use \(V^{N2}\) approximation.
To begin the CI calculation
The ci_basis
options means we will use twoparticle CSFs formed from combinations of singleparticle 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 precomputes all the twobody 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 gfactor. 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 nonrelativistic notation), and the %
column shows the combined constructions from each relativistic configuration with the same nonrelativistic 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 twobody \(\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 valencespace 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 frequencydependence 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.