Duo: a general program for calculating spectra of diatomic molecules ^{1}^{1}1This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect
Abstract
Duo is a general, userfriendly 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 closedshell diatomics) but also for the general case of an arbitrary number and type of couplings between electronic states (typical for openshell diatomics and excited states). Possible couplings include spinorbit, angular momenta, spinrotational and spinspin. Corrections due to nonadiabatic effects can be accounted for by introducing the relevant couplings using socalled BornOppenheimer breakdown curves.
Duo requires userspecified 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.
keywords:
diatomics, spectroscopy, onedimensional Schrödinger equation, excited electronic states, intramolecular perturbation, coupledchannel radial equations, transition probabilities, intensitiesProgram 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, spinorbit 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 BornOppenheimer or adiabatic approximation Cederbaum (2004) the rotationalvibrational (rovibrational) energy levels of a diatomic molecule with nuclei and and in a electronic state are given by the solution of the onedimensional Schrödinger equation:
(1) 
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 onedimensional Schrödinger equation is a wellstudied 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 BalintKurti (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” CooleyNumerov 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 semiclassical RydbergKleinRees (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, gridbased 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 coworkers 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 coupledstate 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, userfriendly 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 spinorbit, electronicrotational, spinrotational and spinspin couplings. Duo also has auxiliary capabilities such as interpolating and extrapolating curves and calculating lists of line positions and line intensities (socalled line lists). Duo is currently being used as part of the ExoMol project (Tennyson and Yurchenko, 2012), whose aim is to generate hightemperature 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 centreofmass motion and having introduced a bodyfixed set of Cartesian axes with origin at the centre of nuclear mass and with the axis along the internuclear direction the nonrelativistic Hamiltonian of a diatomic molecule can be written as Islampour and Miralinaghi (2015); Sutcliffe (2007); Kato (1993); Pack and Hirschfelder (1968); Bunker (1968):
(2) 
where the meaning of the various terms is as follows. is the electronic Hamiltonian and is given by
(3) 
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 masspolarisation term given by
(4) 
where is the total nuclear mass; is the vibrational kinetic energy operator and is given by
(5) 
where is the reduced mass of the molecule. is the rotational Hamiltonian and can be expressed in terms of the bodyfixed rotational angular momentum (AM) operator as
(6) 
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 laboratoryfixed and the bodyfixed 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
(7)  
The approach used by Duo to solve the total rovibronic Schrödinger equation with the Hamiltonian (2) follows closely the standard coupledsurface BornOppenheimer treatment Cederbaum (2004); Kutzelnigg (1997); Hutson and Howard (1980). It is assumed that one has preliminary solved the electronic motion problem with clamped nuclei
(8) 
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):
(9)  
(10) 
where is the symmetry operator corresponding to a reflection through the bodyfixed 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 () onedimensional 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
(11) 
where is a symmetrictop 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 nonvanishing matrix elements of the angular momentum operators in the rotational Hamiltonian (7) are given by the standard rigidrotor expressions Schwenke (2015):
(12)  
(13)  
(14) 
while matrix elements of the spin operators between electronic wave functions (omitting the ‘state’ label for simplicity) are given by
(15)  
(16)  
(17) 
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 nonzero only for . The term containing is called Suncoupling and is nonzero for . The term containing is called Luncoupling and is nonzero for . Finally, the term containing is called spinelectronic and is nonzero 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 semiempirically, 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 masspolarisation Hamiltonian using the electronic wave functions gives rise Herman and Asgharian (1968); Kutzelnigg (1997) to the socalled BornOppenheimer diagonal correction (also called adiabatic correction), which can be added to the BornOppenheimer 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 spinorbit coupling (see section 2.6 for a list of possible additional terms to the Hamiltonian). The vibrational matrix elements
(19) 
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 welldefined parity; parity (even or odd) is defined with respect to inversion of all laboratoryfixed coordinates Kato (1993); Barr et al. (1975); Røeggen (1971); Pack and Hirschfelder (1968) and is equivalent to the reflection operation through the moleculefixed plane, . The parity properties of the basis functions of Eq. (11) are given by Kato Kato (1993)
(20) 
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)
(21) 
where the are expansion coefficients and here is a shorthand index for the basis set labels ‘state’, , , , , , and :
(22) 
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 spinelectronicrotational contribution. In situations where some couplings are absent some approximate quantum numbers can become exact; for example, in the absence of spinorbit and spinspin 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 rovibrational 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).
2.1 Solution of the uncoupled vibrational problem
The main method of solving the radial equation used by Duo is the socalled sinc DVR (discrete variable representation); this method (or closely related ones) has been independently applied to the onedimensional 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
(23) 
where is the matrix representing the kinetic energy and is given in the sinc method by Colbert and Miller (1992); Tannor (2007)
(24) 
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.
Duo obtains all integrals over vibrational coordinates by summation over the grid points:
(25) 
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 floatingpoint 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:
(26) 
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 floatingpoint 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 finitedifference (FD) methods for solving the uncoupled vibrational problem, where the kinetic energy operator in Eq. (23) can be approximated using, for example, a 5point central FD5 formulae:
(27) 
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.
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 5point finitedifference one. Indicatively we recommend considering nonuniform 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.
Several techniques have been developed to deal with quasibound states, most notably in the context of diatomic molecules by Le Roy and coworkers 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 BenItzhak (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, longlived quasibound states (i.e., narrow resonance) can be identified using the present version of Duo by using the socalled stabilization method Hazi and Taylor (1970); Simons (1981); Levebvre (1985); GarciaSucre 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
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 nonrelativistic Hamiltonian (2) caused by spinorbit , spinrotational , spinspin and doubling interactions LefebvreBrion and Field (2004); Kato (1993); Brown et al. (1987); Davis et al. (1988); Brown and Merer (1979); Richards et al. (1981):

The BreitPauli spinorbit operator Marian (2001); Fedorov et al. (2003); Veseth (1970); Richards et al. (1981) has nonzero matrix elements between electronic states obeying the following coupling rules LefebvreBrion and Field (2004): ; ; ; if and the matrix elements is zero (this last rule implies that singlettosinglet matrix elements are zero); electronic states may have nonzero matrix elements with states but matrix elements are zero; finally, in case of homonuclear diatomics, only and matrix elements are nonzero.
The diagonal SO matrix elements determine the spinorbit splitting of a multiplet , where and . Both diagonal and offdiagonal matrix elements of the spinorbit Hamiltonian can be obtained as functions of using quantum chemistry programs.

The nonzero diagonal and offdiagonal matrix elements of operator are given by
(28) (29) where is a dimensionless function of .

The diagonal matrix elements of the operator are taken in the phenomenological form
(30) Both and functions can be obtained either ab initio or empirically.

It is now wellestablished that, at least for states, the small shifts to energy levels due to nonadiabatic interactions with remote states (as opposed to neardegenerate 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
(33) while the rotational kinetic energy operator in Eq. (7) should be replaced by
(34) The functions and are sometimes referred to as BornOppenheimer breakdown (BOB) curves Le Roy and Huang (2002) and can also be interpreted as introducing positiondependent 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 atomcentred Cartesian components Marian (2001) , , , etc.
Duo can accept input in either the Cartesian or the representation. For the Cartesianrepresentation Duo will then transform these inputs to the representation as follows:
(35) 
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
(36) 
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 nonzero matrix elements can be related to only one, nonzero 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 nonzero spinorbit matrix elements as a reference and to reconstruct all other nonzero Cartesian component by
(37) 
as required for Eq. (36).
Offdiagonal 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 postprocessed 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. QuintanaOrtí (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
(39)  
where are the electronically averaged bodyfixed components of the electric dipole moment (in Debye) in the irreducible representation
(40) 
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
(41) 
where is the partition function defined as
(42) 
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 halfinteger) or bosons ( integer) Herzberg (1950); Šimečková et al. (2006):
(43) 
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.
State  

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).
173  1  4.2062E06  
174  1  1.3462E02  
175  1  1.3856E02  
176  1  9.7421E03  
177  1  1.2511E06  
178  1  1.1544E02  
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.
: EinsteinA coefficient in s.
:
Transition wavenumber in cm (optional).
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 nonlinear leastsquares problem where one seeks to minimize the objective function Dennis and Schnabel (1996):
(44) 
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 nonlinear conjugate gradient method for the optimization; in particular, the linearized leastsquare problem is solved by default using the LAPACK subroutine DGELSS, although the alternative builtin 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 nonlinear 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:
(45) 
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: