Duo: a general program for calculating spectra of diatomic molecules 1footnote 11footnote 1This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect

Duo: a general program for calculating spectra of diatomic molecules 111This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect

Sergei N. Yurchenko s.yurchenko@ucl.ac.uk Lorenzo Lodi l.lodi@ucl.ac.uk Jonathan Tennyson j.tennyson@ucl.ac.uk Andrey V. Stolyarov avstol@phys.chem.msu.ru Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, United Kingdom Department of Chemistry, Lomonosov Moscow State University, Leninskiye gory 1/3, 119992 Moscow, Russia

Duo is a general, user-friendly program for computing rotational, rovibrational and rovibronic spectra of diatomic molecules. Duo solves the Schrödinger equation for the motion of the nuclei not only for the simple case of uncoupled, isolated electronic states (typical for the ground state of closed-shell diatomics) but also for the general case of an arbitrary number and type of couplings between electronic states (typical for open-shell diatomics and excited states). Possible couplings include spin-orbit, angular momenta, spin-rotational and spin-spin. Corrections due to non-adiabatic effects can be accounted for by introducing the relevant couplings using so-called Born-Oppenheimer breakdown curves.

Duo requires user-specified potential energy curves and, if relevant, dipole moment, coupling and correction curves. From these it computes energy levels, line positions and line intensities. Several analytic forms plus interpolation and extrapolation options are available for representation of the curves. Duo can refine potential energy and coupling curves to best reproduce reference data such as experimental energy levels or line positions. Duo is provided as a Fortran 2003 program and has been tested under a variety of operating systems.

diatomics, spectroscopy, one-dimensional Schrödinger equation, excited electronic states, intramolecular perturbation, coupled-channel radial equations, transition probabilities, intensities
journal: Computer Physics Communications

Program summary

Program title: Duo
Catalogue number:
Program summary URL:
Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland
Licensing provisions: Standard CPC licence.
No. of lines in distributed program, including test data, etc.: 160 049
No. of bytes in distributed program, including test data, etc.: 13 957 785
Distribution format: tar.gz
Programming language: Fortran 2003.
Computer: Any personal computer.
Operating system: Linux, Windows, Mac OS.
Has the code been vectorized or parallelized?: Parallelized.
Memory required to execute: case dependent, typically  10 MB
Nature of physical problem: Solving the Schrödinger equation for the nuclear motion of a diatomic molecule with an arbitrary number and type of couplings between electronic states.
Solution method: Solution of the uncoupled problem first, then basis set truncation and solution of the coupled problem. A line list can be computed if a dipole moment function is provided. The potential energy and other curves can be empirically refined by fitting to experimental energies or frequencies, when provided.
Restrictions on the complexity of the problem: The current version is restricted to bound states of the system.
Unusual features of the program: User supplied curves for all objects (potential energies, spin-orbit and other couplings, dipole moment etc) as analytic functions or tabulated on a grid is a program requirement.
Typical running time: Case dependent. The test runs provided take seconds or a few minutes on a normal PC.

1 Introduction

Within the Born-Oppenheimer or adiabatic approximation Cederbaum (2004) the rotational-vibrational (rovibrational) energy levels of a diatomic molecule with nuclei and and in a electronic state are given by the solution of the one-dimensional Schrödinger equation:


where is the reduced mass of the molecule and and are the (nuclear) masses of atoms and , respectively. is the potential energy curve (PEC) for the electronic state under study, is the total angular momentum of the molecule and is the vibrational quantum number. The solution of this one-dimensional Schrödinger equation is a well-studied mathematical problem Berezin and Shubin (1991); Simon (2000) for which many efficient numerical methods are available Blatt (382); Shore (1973); Johnson (1977, 1978); Korsch and Laurent (1981); Guardiola and Ros (1982a, b); Stolyarov and Kuz’menko (1987); Lindberg (1988); Marston and Balint-Kurti (1989); Abarenov and Stolyarov (1990); Garza and Vela (1996, 1997); Ishikawa (2002); Utsumi et al. (2004); Wang et al. (2004); the most popular of them is probably the iterative “shooting” Cooley-Numerov Cooley (1961); Cashion (1963); Noumerov (1924) method which is notably used in the program Level due to Le Roy Le Roy (2007).

As well as the ‘direct’ problem of solving the Schrödinger equation for a given PEC, also of great interest is the corresponding inverse problem Karkowski (2009); Weymuth and Reiher (2014), that is the task of determining the potential which leads to a given set of energy levels , typically obtained from experiment. A traditional way of performing this task approximately is to use the semi-classical Rydberg-Klein-Rees (RKR) method Karkowski (2009); a more precise strategy called inverse perturbation analysis (IPA) has been suggested by Kosman and Hinze Kosman and Hinze (1975); Weymuth and Reiher (2014) and a program implementing this approach was presented by Pashov et al Pashov et al. (2000). A different, grid-based fitting strategy has been recently suggested by Szidarovszky and Császár Szidarovszky and Császár (2014). The program DPotFit Le Roy (2006), a companion to Le Roy’s Level, also provides this functionality for isolated states of closed shell diatomics. Indeed, for single potential problems there is an extensive literature on the determination of potential curves from experimental data; in this context we particularly note the work of Coxon and Hajigeorgiou Coxon and Hajigeorgiou (1992, 1999, 2004) and Le Roy and co-workers Le Roy (2004); Le Roy et al. (2011); Meshkov et al. (2014); Walji et al. (2015).

When the diatomic molecule has a more complex electronic structure (i.e., the electronic term is not ) the situation is more complicated, as interactions between the various terms are present and it is not possible to treat each electronic state in isolation. Although there are a growing number of studies treating coupled electronic states, for example see Refs. Carrington et al. (1995); Tamanis et al. (2002); Bergeman et al. (2003); Meshkov et al. (2005); Hutson et al. (2008); Zhang et al. (2010); Gopakumar et al. (2013); Brooke et al. (2014), there appears to be no general program available for solving the coupled problem, the closest being a general coupled-state program due to Hutson Hutson (1994). We have therefore developed a new computer program, Duo, particularly to deal with such complex cases.

Duo is a flexible, user-friendly program written in Fortran 2003 and capable of solving both the direct and the inverse problem for a general diatomic molecule with an arbitrary number and type of couplings between electronic states, including spin-orbit, electronic-rotational, spin-rotational and spin-spin couplings. Duo also has auxiliary capabilities such as interpolating and extrapolating curves and calculating lists of line positions and line intensities (so-called line lists). Duo is currently being used as part of the ExoMol project (Tennyson and Yurchenko, 2012), whose aim is to generate high-temperature spectra for all molecules likely to be observable in exoplanet atmospheres in the foreseeable future. Completed studies based on the use of Duo include ones on AlO Patrascu et al. (2014, 2015), ScH Lodi et al. (2015), CaO Yurchenko et al. (2015) and VO McKemmish et al. (2015). Our methodology is the subject of a very recent topical review Tennyson et al. (2016a).

This paper is organised as follows. In Section 2 we review the theory and the basic equations used by Duo to solve the coupled nuclear motion problem for diatomics. In Section 3 we discuss the calculation of molecular line intensities and line lists. Section 4 is devoted to the inverse problem, i.e. to the refinement (‘fitting’) of potential and coupling curves so that they reproduce a set of reference energy levels or line positions. Section 5 reviews the functional forms implemented for the various curves. In Section 6 the program structure is explained. Finally, we draw our conclusions in Section 7. Technical details on program usage such as detailed explanations of the program options and sample inputs are reported in a separate user’s manual.

2 Method of solution

After separating out the centre-of-mass motion and having introduced a body-fixed set of Cartesian axes with origin at the centre of nuclear mass and with the axis along the internuclear direction the non-relativistic Hamiltonian of a diatomic molecule can be written as Islampour and Miralinaghi (2015); Sutcliffe (2007); Kato (1993); Pack and Hirschfelder (1968); Bunker (1968):


where the meaning of the various terms is as follows. is the electronic Hamiltonian and is given by


where is the Coulomb electrostatic interactions between all particles (electrons and nuclei) and we indicated with the internuclear coordinate and collectively with the full set of electron coordinates; is the mass-polarisation term given by


where is the total nuclear mass; is the vibrational kinetic energy operator and is given by


where is the reduced mass of the molecule. is the rotational Hamiltonian and can be expressed in terms of the body-fixed rotational angular momentum (AM) operator as


In turn, the rotational AM can be expressed as where is the total AM, is the electron orbital AM and is the electron spin AM. The total AM operator acts on the Euler angles relating the laboratory-fixed and the body-fixed Cartesian frame and its expression can be found, e.g., in Ref. Islampour and Miralinaghi (2015). Introducing the ladder operators , and we can express the rotational Hamiltonian as


The approach used by Duo to solve the total rovibronic Schrödinger equation with the Hamiltonian (2) follows closely the standard coupled-surface Born-Oppenheimer treatment Cederbaum (2004); Kutzelnigg (1997); Hutson and Howard (1980). It is assumed that one has preliminary solved the electronic motion problem with clamped nuclei


for all electronic states of interest. The electronic wave functions depend on the electron coordinates and parametrically on the internuclear distance and can be labelled by total spin , projection of along the body fixed axis , projection of along the body fixed axis and by a further label ‘state’= which counts over the electronic curves. For the spacial part of the electronic wave functions is doubly degenerate; we choose the degenerate components so that they satisfy the following conditions Kato (1993):


where is the symmetry operator corresponding to a reflection through the body-fixed -plane (parity operator) and for states and for all other states.

Once the potential energy curves have been obtained, for example using an ab initio quantum chemistry program, Duo solves the rotationless () one-dimensional Schrödinger given by Eq. (1) separately for each electronic curve , producing a set of vibrational eigenvalues and vibrational wave functions , where is the vibrational quantum number assigned on the basis of the energy ordering; technical details on this step are given in Section 2.1. A subset of vibrational functions are selected to form a basis set of rovibronic basis functions defined by


where is a symmetric-top eigenfunction Islampour and Miralinaghi (2015) (a function of the Euler angles) and describes the overall rotation of the molecule as a whole, and is the projection of the total angular momentum along the laboratory axis . Only combinations of and which satisfy are selected in the rovibronic basis set (11). The selection of vibrational basis functions to retain can be made either by specifying an energy threshold (all vibrational states below the threshold are retained) or by specifying a maximum vibrational quantum number .

The rovibrational basis set (11) is used to solve the complete rovibronic Hamiltonian given by Eq. (2); this amounts to using an expansion in Hund’s case (a) functions to solve the coupled problem. In particular, the ladder operators appearing in couple rovibrational states belonging to different electronic states; specifically, the non-vanishing matrix elements of the angular momentum operators in the rotational Hamiltonian (7) are given by the standard rigid-rotor expressions Schwenke (2015):


while matrix elements of the spin operators between electronic wave functions (omitting the ‘state’ label for simplicity) are given by


The coupling rules for the Hamiltonian (7) are as follows; the first line in Eq. (7) is the diagonal part of the rotational Hamiltonian, i.e. is non-zero only for . The term containing is called S-uncoupling and is non-zero for . The term containing is called L-uncoupling and is non-zero for . Finally, the term containing is called spin-electronic and is non-zero for .

Matrix elements of the orbital AM operators and when averaged over the electronic wave functions give rise to -dependent curves; these can be computed by ab initio methods Colbourn and Wayne (1979) or estimated semi-empirically, for example using quantum defect theory Stolyarov et al. (1997); Stolyarov and Child (2001).

The expectation value of the sum of the vibrational and the mass-polarisation Hamiltonian using the electronic wave functions gives rise Herman and Asgharian (1968); Kutzelnigg (1997) to the so-called Born-Oppenheimer diagonal correction (also called adiabatic correction), which can be added to the Born-Oppenheimer PEC if desired.

At this stage Duo builds the full Hamiltonian matrix in the basis of Eq. (11) and using the Hamiltonian operator (2), possibly complemented by supplementary terms such as spin-orbit coupling (see section 2.6 for a list of possible additional terms to the Hamiltonian). The vibrational matrix elements


for all operators including couplings, dipole moments, corrections etc. between different vibrational basis set functions are computed and stored; note that in the equation above and indicate different electronic states if .

At this point a basis set transformation is carried out, from the basis given by Eq. (11) to a symmetrized one in which the basis functions have well-defined parity; parity (even or odd) is defined with respect to inversion of all laboratory-fixed coordinates Kato (1993); Barr et al. (1975); Røeggen (1971); Pack and Hirschfelder (1968) and is equivalent to the reflection operation through the molecule-fixed plane, . The parity properties of the basis functions of Eq. (11) are given by Kato Kato (1993)


where for states and for all other states. The symmetrized basis functions are symmetric () or antisymmetric () with respect to . Use of the symmetrized basis set leads to two separate Hamiltonian blocks with defined parities.

The two parity blocks are then diagonalized (see Section 2.8 for technical details), to obtain the final rovibronic eigenvalues and corresponding eigenfunctions , where is the parity quantum number, is a simple counting index. The corresponding rovibronic wave function can be written as an expansion in the basis set (11)


where the are expansion coefficients and here is a shorthand index for the basis set labels ‘state’, , , , , , and :


As the notation above indicates, in the general case the only good quantum numbers (i.e. labels associated with the eigenvalues of symmetry operators) are the total angular momentum value and the parity . Nevertheless, Duo analyzes the eigenvectors and assigns energy levels with the approximate quantum numbers ‘state’, , , and on the basis of the largest coefficient in the basis set expansion (21). It should be noted that the absolute signs of and are not well defined, only their relative signs are. This is related to the symmetry properties of the eigenfunctions of the Hamiltonian (2), which are 50/50 symmetric and antisymmetric mixtures of the and contributions. Therefore the absolute value of the quantum number is required additionally in order to fully describe the spin-electronic-rotational contribution. In situations where some couplings are absent some approximate quantum numbers can become exact; for example, in the absence of spin-orbit and spin-spin interactions the basis functions (11) with different values of spin do not interact and, hence, becomes a “good” quantum number. As another example, without the presence of or states there is no mechanism for the ro-vibrational functions of a state to interact with other electronic states and therefore the corresponding eigenfunctions will have well defined values of .

Table 1 gives an example of a Duo output with the energy term values computed for the case of the first three electronic states, , , and , of AlO Patrascu et al. (2014).


       J      N        Energy/cm  State   v  lambda spin   sigma   omega  parity
       0.5    1          0.000000   1     0   0     0.5     0.5     0.5   +    ||X2SIGMA+
       0.5    2        965.435497   1     1   0     0.5     0.5     0.5   +    ||X2SIGMA+
       0.5    3       1916.845371   1     2   0     0.5     0.5     0.5   +    ||X2SIGMA+
       0.5    4       2854.206196   1     3   0     0.5     0.5     0.5   +    ||X2SIGMA+
       0.5    5       3777.503929   1     4   0     0.5     0.5     0.5   +    ||X2SIGMA+
       0.5    6       4686.660386   1     5   0     0.5     0.5     0.5   +    ||X2SIGMA+
       0.5    7       5346.116382   2     0   1     0.5    -0.5     0.5   +    ||A2PI
       0.5    8       5581.906844   1     6   0     0.5     0.5     0.5   +    ||X2SIGMA+
       0.5    9       6066.934830   2     1   1     0.5    -0.5     0.5   +    ||A2PI
       0.5   10       6463.039443   1     7   0     0.5     0.5     0.5   +    ||X2SIGMA+
       0.5   11       6778.997803   2     2   1     0.5    -0.5     0.5   +    ||A2PI
       0.5   12       7329.427637   1     8   0     0.5     0.5     0.5   +    ||X2SIGMA+
       0.5   13       7483.145675   2     3   1     0.5    -0.5     0.5   +    ||A2PI
       0.5   14       8159.170405   2     4   1     0.5    -0.5     0.5   +    ||A2PI
       0.5   15       8201.467744   1     9   0     0.5     0.5     0.5   +    ||X2SIGMA+
       0.5   16       8857.266385   2     5   1     0.5    -0.5     0.5   +    ||A2PI


Table 1: Sample Duo energy output for AlO Patrascu et al. (2014). The energy is given in cm, and the exact (J, n, parity) and approximate (state, v (), lambda (), spin (), sigma (), and omega ()) quantum numbers. The final column contains labels of the electronic states as given by the user and the separator || is to facilitate selecting the energy entries in the program output.

2.1 Solution of the uncoupled vibrational problem

The main method of solving the radial equation used by Duo is the so-called sinc DVR (discrete variable representation); this method (or closely related ones) has been independently applied to the one-dimensional Schrödinger equation by various authors Lund and Riley (1984); Guardiola and Ros (1982a, b); Colbert and Miller (1992).

In this method the coordinate is truncated to an interval and discretized in a grid of uniformly spaced points (where ) with grid step . The Schrödinger Eq. (1) is then transformed to an ordinary matrix eigenvalue problem


where is the matrix representing the kinetic energy and is given in the sinc method by Colbert and Miller (1992); Tannor (2007)


and while the vector contains the values of at the grid points. The resulting real symmetric matrix is then diagonalized (see section 2.8 for details). The sinc DVR method usually provides very fast (faster than polynomial) convergence of the calculated energies and wave functions with respect to the number of the grid points, . Figure 1 a) shows the convergence for three energy levels of a Morse potential, showing a rate of convergence approximately exponential with respect to the number of grid points.

Figure 1: Illustrative examples of the fast convergence for energies (plot a) and matrix elements (plot b) using the sinc DVR method; in both cases the rate of convergence is approximately exponential with respect to the number of grid points. Results are for a Morse potential which approximately models the ground electronic state of the CO molecule,  Å,  cm,  cm with atomic masses for carbon and oxygen. A uniformly spaced grid was used, keeping fixed  Å,  Å. In plot a) we show absolute errors for the , and energy levels; in plot b) we show relative errors of matrix elements of the type ; the flattening of the error for large numbers of grid points is due to the numerical error present in floating point calculations (see text).

Duo obtains all integrals over vibrational coordinates by summation over the grid points:


The rectangle rule is simple and exponentially accurate for integration over infinite range of functions which decay fast (exponentially or faster) and which do not have singularities in the complex plane close to the real axis Trefethen and Weideman (2014). We illustrate in fig. 1 b) the quick convergence of matrix elements of the type for a Morse potential; analytical formulae for matrix elements of this kind are available from the literature Gallas (1980); Rong et al. (356) and were used to obtain exact reference values. In plot 1 b) it is apparent that the accuracy of matrix elements does not improve beyond a certain value; for example, the matrix elements always has less than about 10 significant digits no matter how many points are used. This behaviour is completely normal and expected when performing floating-point calculations with a fixed precision; Duo uses double precision numbers with a machine epsilon Higham (1993) and the expected relative error due to the finite precision in the sum given by eq. (25) is given by, indicating with the value of the sum performed with infinite accuracy and with the value obtained with finite accuracy:


The expression above implies that whenever the matrix element of a function comes out very small with respect to the value of significant digits will be lost; there are techniques such as Kahan compensated summation Higham (1993) which reduce the error above by a factor but these have not been implemented at this time.

A prime example of this situation if given by the line intensities of very high vibrational overtones; in a recent study Medvedev et al. (2015) observed that matrix elements of the type for the CO molecule when computed with double precision floating-point arithmetic decrease approximately exponentially (as expected on the basis of theoretical models and as confirmed by quadruple precision calculations) for , when they reach the value of about  D. This situation is fully expected on the basis of the considerations above but it should never constitute a problem in practice.

Apart from the sinc DVR, Duo implements finite-difference (FD) methods for solving the uncoupled vibrational problem, where the kinetic energy operator in Eq. (23) can be approximated using, for example, a 5-point central FD5 formulae:


and furthermore with . Note that the expression above gives incorrect results for the first two and last two grid points, but this does not matter as long as the grid truncation parameters and are chosen so that near the borders of the grid.

The formulae (27) lead to a symmetric pentadiagonal banded matrix, which can in principle be diagonalized more efficiently than a generic dense matrix. However, the convergence of the eigenvalues is much slower, with error decreasing as instead of .

2.2 Levels lying close to dissociation

A general requirement for convergence is that both the inner and the outer grid truncation values and should be chosen such that and are both much larger than . A problem arises when one is trying to converge states very close to the dissociation limit, as such loosely bound states can extend to very large values of and therefore require an excessive number of points when a uniformly spaced grid is used; this is illustrated in fig. 2.

Figure 2: Illustrative example of the effect of the outer grid truncation parameter on energy levels close to dissociation. Data are relative to a Morse potential with  cm,  cm,  Å,  Da. This potential supports 15 bound states ( to ) and we consider in this example the three highest-energy ones, with energies  cm,  cm,  cm. In all calculation we fixed  Å and the grid step to 0.05 Å. The dotted vertical lines are the outer turning points for the three states, i.e. the points such that ; the error in the computed energy levels is expected to decrease exponentially when . The plot shows that to converge the last energy level a very large is required, which in turn leads to a large number of grid points when they are uniformly spaced. Specifically, to converge to  cm it is sufficient to choose  Å, leading to 120 points; for we need  Å and 180 points; for we need  Å and 1500 points.

Excited states of alkali diatoms such as Li Coxon and Melville (2006), Na Qi et al. (2007) or K Falke et al. (2006) constitute an important class of systems for which large are needed; such systems are prime choices for studies of ultracold atoms and molecules Tiemann and Knöckel (2013) and often require grids extending up to several tens Coxon and Melville (2006); Qi et al. (2007) or even hundreds Falke et al. (2006) of Angstroms.

In such cases it may be beneficial to use a non uniform grid; Duo implements the adaptive analytical mapping approach of Meshkov et al Meshkov et al. (2008) and offers several mapping choices, which are described in the manual. However, at this time support for non uniform grids should be considered experimental and they cannot be used in combination with the sinc DVR method but only with the less efficient 5-point finite-difference one. Indicatively we recommend considering non-uniform grids only when is required to be larger than 50 Å.

2.3 States beyond the dissociation limit

Potential curves with local maxima higher than the dissociation limit of the potential for may support shape resonances, i.e. metastable states in which the two atoms are trapped for a finite time in the potential well but eventually dissociate. Such states are also known as quasibound or tunnelling predissociation states. For the rotational potential will practically always introduce such a maximum, and the corresponding quasibound levels are known as orbiting resonances or rotationally predissociating states, see fig. 3 a) for an example.

Figure 3: Example of a quasibound, orbiting resonance state. Plot a): Morse potential for and (same parameters as in fig. 2); the potential for has a local maximum higher than the dissociation limit and supports one quasibound state with energy 11 931.1 cm. Plot b): eigenvalues for as a function of the outer grid truncation . The quasibound state manifests itself as a series of avoided crossings.

Several techniques have been developed to deal with quasibound states, most notably in the context of diatomic molecules by Le Roy and co-workers Le Roy and Bernstein (5114); Le Roy and Liu (1978); Connor and Smith (1981); Pryce (1994); Riss and Meyer (1993); Čížek and Horáček (1996); Sidky and Ben-Itzhak (1999); Huang and Le Roy (2003, 2007). At the moment Duo does not provide any explicit functionality to treat quasibound states, although we plan to rectify this deficiency in future versions.

Nevertheless, long-lived quasibound states (i.e., narrow resonance) can be identified using the present version of Duo by using the so-called stabilization method Hazi and Taylor (1970); Simons (1981); Levebvre (1985); Garcia-Sucre and Levebvre (1986); Mandelshtam et al. (1993); Martín (1999). In one version of this approach energy levels are computed for increasing values of the outer grid truncation and then plotted as function of ; quasibound states manifests themselves by being relatively stable with respect to increase of and undergo a series of avoided crossings, see fig. 3 b) for an example. From an analysis of these curves it is also possible to compute the lifetime of the quasibound state Mandelshtam et al. (1993).

2.4 Printing the wave functions

Both the vibrational basis functions , see eq. (11), and the final (, electronically coupled or both) rovibronic wave functions coefficients , see eq. (21), can be written to a file for further analysis, e.g. for plotting purposes or for the computation of factors; see the manual for details.

2.5 Convergence of rotationally excited states

In our approach calculations are performed using a basis expansion in terms of the wave functions. As a guideline it was found by numerical experimentation that in order to obtain converged results for rotationally excited states up to one has to use a vibrational basis set of size only slightly larger than and that a reasonable minimum value for the size of the vibrational basis set is given by . For example, to converge rotationally excited levels up to it should be sufficient to use a vibrational basis set of size 40.

2.6 Additional terms in the Hamiltonian

Duo supports the inclusion of a number of terms additional to the non-relativistic Hamiltonian (2) caused by spin-orbit , spin-rotational , spin-spin and -doubling interactions Lefebvre-Brion and Field (2004); Kato (1993); Brown et al. (1987); Davis et al. (1988); Brown and Merer (1979); Richards et al. (1981):

  1. The Breit-Pauli spin-orbit operator Marian (2001); Fedorov et al. (2003); Veseth (1970); Richards et al. (1981) has non-zero matrix elements between electronic states obeying the following coupling rules Lefebvre-Brion and Field (2004): ; ; ; if and the matrix elements is zero (this last rule implies that singlet-to-singlet matrix elements are zero); electronic states may have non-zero matrix elements with states but matrix elements are zero; finally, in case of homonuclear diatomics, only and matrix elements are non-zero.

    The diagonal SO matrix elements determine the spin-orbit splitting of a multiplet , where and . Both diagonal and off-diagonal matrix elements of the spin-orbit Hamiltonian can be obtained as functions of using quantum chemistry programs.

  2. The nonzero diagonal and off-diagonal matrix elements of operator are given by


    where is a dimensionless function of .

  3. The diagonal matrix elements of the operator are taken in the phenomenological form


    Both and functions can be obtained either ab initio or empirically.

  4. The lambda-doubling (LD) couplings for a state in the -representation (Eq. (11)) are of the following three types Brown and Merer (1979):


    where and are related to the conventional terms as given by Brown and Merer (1979):

  5. It is now well-established that, at least for states, the small shifts to energy levels due to non-adiabatic interactions with remote states (as opposed to near-degenerate ones) can be accurately modelled by modifying the vibrational and rotational energy operators in the Hamiltonian Pachucki and Komasa (2008); Watson (2004); Le Roy (1999); Herman and Ogilvie (1998); Bunker and Moss (1977); Herman and Asgharian (1968); specifically, the vibrational energy operator in Eq. (5) is replaced by


    while the rotational kinetic energy operator in Eq. (7) should be replaced by


    The functions and are sometimes referred to as Born-Oppenheimer breakdown (BOB) curves Le Roy and Huang (2002) and can also be interpreted as introducing position-dependent vibrational and rotational masses; they are sometimes expressed in terms of the dimensionless -factor functions and by and . The rotational function can be determined experimentally by analysis of the Zeeman splitting of energy levels due to an external magnetic field Ogilvie et al. (2000).

2.7 Representation of the couplings

Duo assumes that the coupling matrix elements and the transition dipole moments are given in the representation of the basis functions (11) corresponding to Hund’s case (a). In this representation the component is diagonal and has a signed value (see Eq.(13)) and therefore it will be referred to as the -representation. It can be shown that by choosing appropriate phase factors for the electronic wave functions all coupling matrix elements in the -representation can be made real; note that in this representation the electronic wave functions are complex numbers, as they contain a factor of the kind , where is the angle corresponding to rotation around the axis Barr et al. (1975). On the other hand quantum chemistry programs such as Molpro Werner et al. (2012) normally work with real wave functions and consequently compute matrix elements in this representation, which we call Cartesian as the electronic wave functions are ultimately expressed in terms of atom-centred Cartesian components Marian (2001) , , , etc.

Duo can accept input in either the Cartesian or the -representation. For the Cartesian-representation Duo will then transform these inputs to the -representation as follows:


where and denote Cartesian components that correspond to the and Abelian point group symmetries, respectively. are the elements of the unitary transformation from the Cartesian to the -representation. The obvious way to reconstruct this transformation is to diagonalize the Cartesian representation of the matrix. Thus the transformed matrix elements in the -representation are given by


or, in tensorial form where correspond to a () electronic state with .

In principle all Cartesian matrix elements must be provided to perform the transformation in Eq. (36). However, by means of the coupling rules all non-zero matrix elements can be related to only one, non-zero reference matrix element. For example, the matrix element between and is zero because it corresponds to a simultaneous change of and by . This property together with the help of Eq. (35) allows one to use the non-zero spin-orbit matrix elements as a reference and to reconstruct all other non-zero Cartesian component by


as required for Eq. (36).

Off-diagonal matrix elements of the various operators included into our model, i.e. the various couplings between electronic states, are subject in actual calculations to arbitrary changes of sign due to the sign indeterminacy of the electronic wave functions computed at different geometries. Often the phases of each ab initio coupling have to be post-processed in order to provide a consistent, smooth function of . It is important that the relative phases between different elements preserved. This issue is illustrated graphically by Patrascu et al. (2014), where different ab initio coupling curves of AlO obtained with Molpro were presented. Transition dipole moment functions, discussed in the next section, also may exhibit phase changes  Tennyson (2014), which should be corrected using the same phase convention used for other matrix elements (Patrascu et al., 2014).

2.8 Computational considerations

Duo uses the matrix diagonalization routines DSYEV or, optionally, DSYEVR from the LAPACK library (Anderson et al., 1999). The subroutine DSYEVR uses the multiple relatively robust representations algorithm and is expected to be faster than DSYEV, which is based on the QR algorithm Demmel et al. (2008); Van Zee and Van De Geijnand G. Quintana-Ortí (2014); however, the current version of DSYEVR is poorly parallelized and therefore not recommended for parallel environments.

The dimension of the final rovibrational Hamiltonian matrix depends on the number of vibrational functions selected, the number of electronic states present, the spin multiplicities of the electronic states and the quantum number. For example, for electronic states, vibrational functions are retained for each of them and denoting with the average spin multiplicity, the size of the Hamiltonian matrix is approximately given by and the size of the parity matrix to be diagonalized is half of this value. The size of each block of the Hamiltonian reaches dimensions of the order of a thousand only for rather complicated cases (e.g., , and ) and consequently the time taken to compute the energy levels for a given is usually only a small fraction of a second.

3 Line intensities and line lists

The Einstein coefficient (in ) for a transition is computed as


where are the electronically averaged body-fixed components of the electric dipole moment (in Debye) in the irreducible representation


and the index is defined by Eq. (22). The vibrationally averaged transition dipole moments are computed using the vibrational wave functions .

The absorption line intensity is then given by


where is the partition function defined as


is the nuclear statistical weight factor, is the second radiation constant, is the term value, and is the temperature in K. For heteronuclear molecules is a total number of combinations of nuclear spins as given by where and are the corresponding nuclear spins. For a homonuclear molecules, these combinations are distributed among the four symmetries , , , , where is the parity of the molecule with respect to and is the property of the total rovibronic wave function to be symmetric/asymmetric upon upon inversion Kato (1993). In the representation of point group symmetry, this corresponds to , , , and . Thus, for the case two different values are necessary and these depend on whether the nuclei are fermions ( is half-integer) or bosons ( integer) Herzberg (1950); Šimečková et al. (2006):


For example, carbon C has and therefore for the C molecule are 1 for , and 0 for , states, respectively.

The computed Einstein coefficients can be used to compute radiative lifetimes of individual states and cooling functions in a straightforward manner Tennyson et al. (2016b).

3.1 Line list format

A line list is defined as a catalogue of transition frequencies and intensities Tennyson et al. (2013). In the basic ExoMol format Tennyson et al. (2013), adopted by Duo, a line list consists of two files: ‘States’ and ‘Transitions’; an example for the molecule AlO is given in Tables 2 and 3. The ‘States’ file contains energy term values supplemented by the running number , total degeneracy , rotational quantum number (all obligatory fields) as well as quantum numbers , , parity (), , and the electronic state label (e.g. X2Sigma+). The ‘Transitions’ file contains three obligatory columns, the upper and lower state indexes and which are running numbers from the ‘State’ file, and the Einstein coefficient . For the convenience we also provide the wavenumbers as the column 4. The line list in the ExoMol format can be used to simulate absorption or emission spectra for any temperature in a general way. Note that ExoMol format has recently been significantly extended Tennyson et al. (2016c) but structure of the States and Transitions file has been retained.

1 0.000000 12 0.5 + e X2SIGMA+ 0 0 0.5 0.5
2 965.435497 12 0.5 + e X2SIGMA+ 1 0 0.5 0.5
3 1916.845371 12 0.5 + e X2SIGMA+ 2 0 0.5 0.5
4 2854.206196 12 0.5 + e X2SIGMA+ 3 0 0.5 0.5
5 3777.503929 12 0.5 + e X2SIGMA+ 4 0 0.5 0.5
6 4686.660386 12 0.5 + e X2SIGMA+ 5 0 0.5 0.5
7 5346.116382 12 0.5 + e A2PI 0 1 0.5 0.5
8 5581.906844 12 0.5 + e X2SIGMA+ 6 0 0.5 0.5
9 6066.934830 12 0.5 + e A2PI 1 1 0.5 0.5
10 6463.039443 12 0.5 + e X2SIGMA+ 7 0 0.5 0.5
11 6778.997803 12 0.5 + e A2PI 2 1 0.5 0.5
12 7329.427637 12 0.5 + e X2SIGMA+ 8 0 0.5 0.5
13 7483.145675 12 0.5 + e A2PI 3 1 0.5 0.5
14 8159.170405 12 0.5 + e A2PI 4 1 0.5 0.5
15 8201.467744 12 0.5 + e X2SIGMA+ 9 0 0.5 0.5
16 8857.266385 12 0.5 + e A2PI 5 1 0.5 0.5
17 9029.150380 12 0.5 + e X2SIGMA+ 10 0 0.5 0.5
18 9535.195842 12 0.5 + e A2PI 6 1 0.5 0.5
19 9854.882567 12 0.5 + e X2SIGMA+ 11 0 0.5 0.5
20 10204.019475 12 0.5 + e A2PI 7 1 0.5 0.5
21 10667.668381 12 0.5 + e X2SIGMA+ 12 0 0.5 0.5
22 10864.560220 12 0.5 + e A2PI 8 1 0.5 0.5
23 11464.897083 12 0.5 + e X2SIGMA+ 13 0 0.5 0.5
24 11519.212123 12 0.5 + e A2PI 9 1 0.5 0.5
25 12156.974798 12 0.5 + e A2PI 10 1 0.5 0.5
26 12257.694655 12 0.5 + e X2SIGMA+ 14 0 0.5 0.5
27 12793.671660 12 0.5 + e A2PI 11 1 0.5 0.5
28 13030.412255 12 0.5 + e X2SIGMA+ 15 0 0.5 0.5
29 13421.583651 12 0.5 + e A2PI 12 1 0.5 0.5
30 13790.933964 12 0.5 + e X2SIGMA+ 16 0 0.5 0.5

: State counting number.

: State energy in cm.

: State degeneracy.

: Total angular momentum.

: Total parity.

: Rotationless parity.

State: Electronic state label.

: State vibrational quantum number.

: Absolute value of (projection of the electronic angular momentum).

: Absolute value of (projection of the electronic spin).

: Absolute value of (projection of the total angular momentum).

Table 2: Extract from the output ‘State’ file produced by Duo for the AlO molecule Patrascu et al. (2015).
173 1 4.2062E-06
174 1 1.3462E-02
175 1 1.3856E-02
176 1 9.7421E-03
177 1 1.2511E-06
178 1 1.1544E-02
179 1 6.7561E+02
180 1 4.1345E+00
181 1 2.4676E+03
182 1 3.5289E+01
183 1 4.6960E+03
184 1 1.9771E+02
185 1 6.1540E+03
186 1 4.8062E+03
187 1 1.9401E+03

: Upper state counting number.
: Lower state counting number.
: Einstein-A coefficient in s.
: Transition wavenumber in cm (optional).

Table 3: Sample extracts from the output ‘Transition’ files produced by Duo for the AlO molecule Patrascu et al. (2015).

4 Inverse problem

The inverse problem consists in finding the potential energy and coupling curves which best reproduce a given set of energy levels, , or frequencies (i.e., differences between energy levels), typically extracted from experiment. In the following we will call this optimization process empirical refinement.

4.1 Implementation

The refinement problem can be formulated as a non-linear least-squares problem where one seeks to minimize the objective function Dennis and Schnabel (1996):


where are the calculated energies or frequencies and implicitly depend on the parameters defining the potential and coupling curves. The are weighting factor assigned to each value and may be chosen as where is the experimental uncertainty of . The input weights are automatically renormalized by Duo so that .

Duo uses the non-linear conjugate gradient method for the optimization; in particular, the linearized least-square problem is solved by default using the LAPACK subroutine DGELSS, although the alternative built-in subroutine LINUR is also available. For each curve appearing in Duo it is possible to specify if any given parameter should be refined (fitted) or should be kept fixed to the value given in the input file. The first derivatives with respect to the fitting parameters required for the non-linear least squares are computed using finite differences with a step size taken as 0.1% of the initial values or 0.001 if is initially zero.

4.2 Constrained minimization

In order to avoid unphysical behaviour and also to avoid problems when the amount of experimental data provided is insufficient for determining all the parameters, the shapes of the curves can be contrained to be as close as possible to some reference curves provided in the input (typically ab initio ones)  Tikhonov and Arsenin (1979); Yurchenko et al. (2003, 2011); Szidarovszky and Császár (2014). This is done by including into the fitting objective function not only differences of the computed energy levels but also differences between the refined curves and the reference ones as follows:


where refers to the -th curve, counts over the grid points, are the corresponding weight factors of the individual points normalized to one and are further weight factors defining the relative importance of the corresponding curve. The weights in Eq. (45) are normalized as follows: