How PGOPHER works or, how to write your own version of PGOPHER

To see how PGOPHER works, we will work through the calculations involved in one fairly simple example, but one that is complicated enough to show the full potential complexities. The example we pick here is the A²Σ⁺ ← X²Π transition in the OH radical. This is a strong transition in the UV, and is often visible in flames. We simulate here the origin band, v'=0 - v"=0. The required steps break down as follows:

We start with the excited state, as this is simpler than the ground state in OH. The essential pattern of the calculation is shown with just three constants - the state origin (32402 cm^-1), the rotational constant, B=16.96 cm^-1 and the spin rotation constant, γ=0.2 cm^-1. The most recent set of constants has many more constants, but these three will give the basic pattern of energy levels. The first step in calculating energy levels is to choose a basis set; for linear molecules PGOPHER uses a Hund's case (a) basis, in which the basis states have well defined values of Λ (here 0), electron spin S=½ and its projection Σ=±½ on the molecular axis. A Σ state would actually be considered classic Hund's case (b), but a Hund's case (a) basis will give identical results, provided the basis used is complete, which is always the case for PGOPHER (unless the OmegaSelect option is not All) . Given the choice of basis set, the next step is to work out the matrix elements. The formulae PGOPHER uses can be found by right clicking on a state and selecting "Matrix Elements", which gives the following results in this case:

Matrix Elements in A v=0

<J,Lambda=0,Omega=-0.5|H|J,Lambda=0,Omega=-0.5> = 
+  Origin*1
+  B*((4*J*(J+1)+1)/4)
+  gamma*(-1/2)

<J,Lambda=0,Omega=0.5|H|J,Lambda=0,Omega=0.5> = 
+  Origin*1
+  B*((4*J*(J+1)+1)/4)
+  gamma*(-1/2)

<J,Lambda=0,Omega=-0.5|H|J,Lambda=0,Omega=0.5> =
+  B*((-sqrt(4*J*(J+1)+1))/2)
+  gamma*(sqrt(4*J*(J+1)+1)/4)

<J,Lambda=0,Omega=0.5|H|J,Lambda=0,Omega=-0.5> = 
+  B*((-sqrt(4*J*(J+1)+1))/2)
+  gamma*(sqrt(4*J*(J+1)+1)/4)

Any matrix elements not listed should be taken as zero. PGOPHER uses a basic symbolic algebra system to work these out for printing, so the results look a little cumbersome.  A little rearranging gives the following more compact form:

The first is the diagonal matrix element, and contains the state origin and the basic BJ(J+1) rotational dependence.  The second mixes the two states in the basis, and it is the presence of this term that means that the state is best described as case (b), rather than case (a).

The first step in calculating the energy levels using these matrix elements is to take symmetry adapted combinations.  Then basis states used above do not have well defined parity, but taking plus and minus combinations does give  states with specific parity:

|J ± Ω=½> = 1/sqrt(2){|J,Lambda=0,Omega=+0.5> ± |J,Lambda=0,Omega=-0.5>}

The 1/sqrt(2) ensures the resulting combination is normalised - all the basis states are set up to be normalised and  orthogonal. There are two different ways of specifying the symmetry for linear molecules. The parity is the symmetry of the wavefunction  with respect to inversion of (space fixed) co-ordinates; a properly constructed wavefunction will stay the same (+) or simply change sign (-).  The parity of, say, the + combination above will alternate with J, so e and f notation is more normal; for the A state of OH the + combination has e parity always, but + for even J and - for odd J. Provided JAdjustSym is set at the molecule level, PGOPHER will display both of these, as above, and normally accept either on input. If false, e and f are not displayed. The matrix elements given above can be transformed to matrix elements between the symmetry adapted basis functions. PGOPHER will also display matrix elements between these symmetry adapted linear combinations (for linear molecules only at the moment). For the A state we obtain:

Symmetrized Matrix Elements in A v=0
<J,Lambda=0,Omega=0.5 + |H| + J,Lambda=0,Omega=0.5> = 
+  Origin*1
+  B*((4*J*(J+1)-2*sqrt(4*J*(J+1)+1)+1)/4)
+  gamma*((sqrt(4*J*(J+1)+1)-2)/4)
<J,Lambda=0,Omega=0.5 - |H| - J,Lambda=0,Omega=0.5> = 
+  Origin*1
+  B*((4*J*(J+1)+2*sqrt(4*J*(J+1)+1)+1)/4)
+  gamma*((-sqrt(4*J*(J+1)+1)-2)/4)

which with a little rearranging gives the clearer form:

<J,Ω=±½|H|J,Ω=±½> = Origin + BJ(J+1) + B/4 - ½γ ± (B-½γ)Sqrt(J(J+1)+1/4)

Note that there are no matrix elements between the + and - combinations, as required by symmetry (which is of course why the combinations are used). We can compare the results of this formula with values worked out by PGOPHER - try "View, States", and ensure that the A state is selected in the box in the top line. If, say J = 4.5 and e parity is selected the display will look something like this:

Hamiltonian Matrix for  A, J = 4.5, Symmetry = e (+)
                 |A v=0 4.5 + 0.5>
<A v=0 4.5 + 0.5| 32741.6

(This output is actually taken from the log window with the print level on the right set to "Matrices"). Note the general form of the state  specification with the key given in the states window: |Name J +- Omega>. Name is the manifold and state name (here A and v=0), J is as above,  +/- indicates whether the + or - combination is taken, and Omega is actually |Ω| here. (Λ is omitted in the label as it is always taken as positive in the symmetrised wavefunctions, and is thus implied by the state name.)

In this case there is only one basis symmetrised basis state for any given J and symmetry, so the wavefunctions  will simply be equal to one of the basis states and the energy will be equal to one of the matrix elements above.  The states window will show the resulting energy levels and wavefunctions; for J = 4.5 e parity we have:

Eigenvalues for A, J = 4.5, Symmetry = e (+)
                                  |A v=0 4.5 + 0.5>
  1  A v=0 4.5  4 F1e    32741.600         1.000000

This should be read as the state name on the left has energy 32741.6 cm^-1 and the overall wavefunction is equal to the (symmetrised)  basis sate given. The state name key is given in the states window as Name J N Fn p. The "Name" is the manifold and state name as before,  but now an additional quantum number, N, is given (here 4). (The vector J is the vector sum of N and S, so in this case the quantum number J = N ± ½.) Also given is the F number, here F1 or F2, and alternative way of specifying the  spin orientation, and p indicates the parity (here e). For linear molecules, settings at the molecule level can be used to select  the quantum numbers displayed - see "ShowJ", "ShowOmega", "ShowN", "ShowFNumber" and "Showef". By default these are all set on apart from  "ShowOmega".

The ground state - X²Π

The ground state is more complicated, as the orbital angular momentum, Λ is now not zero, so a full matrix diagonalisation is required. To obtain the full spitting pattern we must also consider lambda doubling. The required constants are the spin-orbit coupling constant, A = -139.21 cm^-1, the rotational constant, B = 18.55 cm^-1, two lambda doubling constants, p = 0.235 cm^-1 and q = -0.0391 cm^-1. The Hund's case (a) basis now includes two possible values of Λ (=±1) in addition to the electron spin S=½ and its projection Σ=±½ as for the excited states. Given two possible values for each of Λ and Σ there are four basis states. The matrix elements are:

Matrix Elements in X v=0

<J,Lambda=-1,Omega=-1.5|H|J,Lambda=-1,Omega=-1.5> = 
+  B*((4*J*(J+1)-3)/4)
+  A*(1/2)

<J,Lambda=-1,Omega=-1.5|H|J,Lambda=-1,Omega=-0.5> = 
   B*((-sqrt(4*J*(J+1)-3))/2)

<J,Lambda=-1,Omega=-0.5|H|J,Lambda=-1,Omega=-1.5> = 
   B*((-sqrt(4*J*(J+1)-3))/2)

<J,Lambda=-1,Omega=-0.5|H|J,Lambda=-1,Omega=-0.5> = 
+  B*((4*J*(J+1)+5)/4)
+  A*(-1/2)

<J,Lambda=1,Omega=0.5|H|J,Lambda=1,Omega=0.5> = 
+  B*((4*J*(J+1)+5)/4)
+  A*(-1/2)

<J,Lambda=1,Omega=0.5|H|J,Lambda=1,Omega=1.5> = 
   B*((-sqrt(4*J*(J+1)-3))/2)

<J,Lambda=1,Omega=1.5|H|J,Lambda=1,Omega=0.5> = 
   B*((-sqrt(4*J*(J+1)-3))/2)

<J,Lambda=1,Omega=1.5|H|J,Lambda=1,Omega=1.5> = 
+  B*((4*J*(J+1)-3)/4)
+  A*(1/2)

<J,Lambda=-1,Omega=-0.5|H|J,Lambda=1,Omega=0.5> = 
+  p*((-sqrt(4*J*(J+1)+1))/4)
+  q*((-sqrt(4*J*(J+1)+1))/2)

<J,Lambda=1,Omega=0.5|H|J,Lambda=-1,Omega=-0.5> = 
+  p*((-sqrt(4*J*(J+1)+1))/4)
+  q*((-sqrt(4*J*(J+1)+1))/2)

<J,Lambda=-1,Omega=-1.5|H|J,Lambda=1,Omega=0.5> = 
   q*(sqrt((4*J*(J+1)-3)*(4*J*(J+1)+1))/8)

<J,Lambda=-1,Omega=-0.5|H|J,Lambda=1,Omega=1.5> = 
   q*(sqrt((4*J*(J+1)-3)*(4*J*(J+1)+1))/8)

<J,Lambda=1,Omega=0.5|H|J,Lambda=-1,Omega=-1.5> = 
   q*(sqrt((4*J*(J+1)-3)*(4*J*(J+1)+1))/8)

<J,Lambda=1,Omega=1.5|H|J,Lambda=-1,Omega=-0.5> = 
   q*(sqrt((4*J*(J+1)-3)*(4*J*(J+1)+1))/8)

Note that Σ = Ω - Λ, so is not specified separately in the basis states. As before,  we can simplify the matrices by taking symmetry adapted combinations.

Symmetrized Matrix Elements in X v=0

<J,Lambda=1,Omega=1.5 + |H| + J,Lambda=1,Omega=1.5> = 
+  B*((4*J*(J+1)-3)/4)
+  A*(1/2)

<J,Lambda=1,Omega=1.5 - |H| - J,Lambda=1,Omega=1.5> = 
+  B*((4*J*(J+1)-3)/4)
+  A*(1/2)

<J,Lambda=1,Omega=1.5 + |H| + J,Lambda=1,Omega=0.5> = 
+  B*((-sqrt(4*J*(J+1)-3))/2)
+  q*(sqrt((4*J*(J+1)-3)*(4*J*(J+1)+1))/8)

<J,Lambda=1,Omega=1.5 - |H| - J,Lambda=1,Omega=0.5> = 
+  B*((-sqrt(4*J*(J+1)-3))/2)
+  q*((-sqrt((4*J*(J+1)-3)*(4*J*(J+1)+1)))/8)

<J,Lambda=1,Omega=0.5 + |H| + J,Lambda=1,Omega=1.5> = 
+  B*((-sqrt(4*J*(J+1)-3))/2)
+  q*(sqrt((4*J*(J+1)-3)*(4*J*(J+1)+1))/8)

<J,Lambda=1,Omega=0.5 - |H| - J,Lambda=1,Omega=1.5> = 
+  B*((-sqrt(4*J*(J+1)-3))/2)
+  q*((-sqrt((4*J*(J+1)-3)*(4*J*(J+1)+1)))/8)

<J,Lambda=1,Omega=0.5 + |H| + J,Lambda=1,Omega=0.5> = 
+  B*((4*J*(J+1)+5)/4)
+  A*(-1/2)
+  p*((-sqrt(4*J*(J+1)+1))/4)
+  q*((-sqrt(4*J*(J+1)+1))/2)

<J,Lambda=1,Omega=0.5 - |H| - J,Lambda=1,Omega=0.5> = 
+  B*((4*J*(J+1)+5)/4)
+  A*(-1/2)
+  p*(sqrt(4*J*(J+1)+1)/4)
+  q*(sqrt(4*J*(J+1)+1)/2)

which can be rearranged to give the more compact form:

<J,Ω=3/2,±|H|J,Ω=3/2,±> = A/2 + BJ(J+1) - 3B/4
<J,Ω=1/2,±|H|J,Ω=1/2,±> = -A/2 + BJ(J+1) + 5B/4 - ±½(p+2q)(J+1/2)
<J,Ω=3/2,±|H|J,Ω=1/2,±> = -B(J(J+1)-¾)^½ ± ½q{(J(J+1)-¾)(J(J+1)+¼)}^½

The symmetry adapted combinations mean that the four basis states split into two pairs of states for any given J, one for each parity. A sample Hamiltonian matrix for J = 4.5, f (-) parity shows a non zero matrix elements mixing the Ω=1/2 and Ω=3/2 states:

Hamiltonian Matrix for  X, J = 4.5, Symmetry = f (-)
                 |X v=0 4.5 - 1.5>|X v=0 4.5 - 0.5>
<X v=0 4.5 - 1.5| 375.595          -90.397194212542
<X v=0 4.5 - 0.5| -90.397194212542 552.297

The resulting energy levels are thus mixtures of these two states:

Eigenvalues for X, J = 4.5, Symmetry = f (-)
                                  |X v=0 4.5 - 1.5>|X v=0 4.5 - 0.5>
  1  X v=0 4.5  4 F1f      337.544        -0.921674        -0.387965
  2  X v=0 4.5  5 F2f      590.348        -0.387965         0.921674

This should be read as the lowest state, labelled 1, has energy 318.994 and is mainly Ω=3/2 - the coefficient in the wavefunction is -0.92. The fact that it is not 1 reflects partial Hund's case (b) character; for higher values of J the matrix elements mixing the two values of Ω increase, and the magnitudes of the coefficients become more nearly equal. In the case (b) limit at high J the coefficients are ±1/2^½. At the lowest J=½ the Ω=3/2 state is not present, so the matrix is 1x1.

Transitions

Given the upper and lower energy levels and wavefunctions the transition frequencies and intensities can now be calculated. The transition frequency is straightforward, just the difference between the energies calculated above. The intensity is more involved. The first step is the transition dipole moments between the basis states; PGOPHER will print these if you right click on the transition moment (<v=0|T(1)|v=0>) and select "Matrix Elements". The first set is:

Reduced Matrix Elements of <v=0|T(1)|v=0>
<J-1,Lambda=0,Omega=-0.5|T(1)|J,Lambda=-1,Omega=-1.5> = sqrt(((2*J+1)*(2*J+3))/(2*J))/2
<J-1,Lambda=0,Omega=-0.5|T(1)|J,Lambda=1,Omega=0.5> = sqrt(((2*J-1)*(2*J+1))/(2*J))/2
<J-1,Lambda=0,Omega=0.5|T(1)|J,Lambda=-1,Omega=-0.5> = sqrt(((2*J-1)*(2*J+1))/(2*J))/2
<J-1,Lambda=0,Omega=0.5|T(1)|J,Lambda=1,Omega=1.5> = sqrt(((2*J+1)*(2*J+3))/(2*J))/2
<J,Lambda=0,Omega=-0.5|T(1)|J,Lambda=-1,Omega=-1.5> = ((-2*J-1)*sqrt(((2*J-1)*(2*J+3))/(2*J*(2*J+1)*(J+1))))/2
<J,Lambda=0,Omega=-0.5|T(1)|J,Lambda=1,Omega=0.5> = (2*J+1)^2/(2*sqrt(2*J*(2*J+1)*(J+1)))
<J,Lambda=0,Omega=0.5|T(1)|J,Lambda=-1,Omega=-0.5> = ((-2*J-1)*(2*J+1))/(2*sqrt(2*J*(2*J+1)*(J+1)))
<J,Lambda=0,Omega=0.5|T(1)|J,Lambda=1,Omega=1.5> = ((2*J+1)*sqrt(((2*J-1)*(2*J+3))/(2*J*(2*J+1)*(J+1))))/2
<J+1,Lambda=0,Omega=-0.5|T(1)|J,Lambda=-1,Omega=-1.5> = sqrt(((2*J-1)*(2*J+1))/(2*(J+1)))/2
<J+1,Lambda=0,Omega=-0.5|T(1)|J,Lambda=1,Omega=0.5> = sqrt(((2*J+1)*(2*J+3))/(2*(J+1)))/2
<J+1,Lambda=0,Omega=0.5|T(1)|J,Lambda=-1,Omega=-0.5> = sqrt(((2*J+1)*(2*J+3))/(2*(J+1)))/2
<J+1,Lambda=0,Omega=0.5|T(1)|J,Lambda=1,Omega=1.5> = sqrt(((2*J-1)*(2*J+1))/(2*(J+1)))/2

These are Reduced transition matrix elements as the M dependence has been taken out; see Transition Moments and Line Strengths for more details. In this case the formula used is:

<J', Ω'||T(1,±)||J,Ω> = (-1)^J'-Ω'[(2J'+1)(2J+1)]^½(J' 1 J -Ω' Ω'-Ω, Ω)

Note the dependence on the change in J and Ω. The line strengths are also printed, which are simply the squares of the reduced matrix elements above. In the absence of state mixing these are Honl-London factors, but in the current case there is significant mixing so they are not useful. (Any calculation of intensity must be done with the unsquared matrix elements.)

The transition matrix elements between the symmetry adapted states are also given:

Symmetrized Reduced Matrix Elements of <v=0|T(1)|v=0>
<J-1,Lambda=0,Omega=0.5 + |T(1)| + J,Lambda=1,Omega=1.5> = sqrt(((2*J+1)*(2*J+3))/(2*J))/2
<J-1,Lambda=0,Omega=0.5 - |T(1)| - J,Lambda=1,Omega=1.5> = sqrt(((2*J+1)*(2*J+3))/(2*J))/2
<J-1,Lambda=0,Omega=0.5 + |T(1)| + J,Lambda=1,Omega=0.5> = sqrt(((2*J-1)*(2*J+1))/(2*J))/2
<J-1,Lambda=0,Omega=0.5 - |T(1)| - J,Lambda=1,Omega=0.5> = (-sqrt(((2*J-1)*(2*J+1))/(2*J)))/2
<J,Lambda=0,Omega=0.5 + |T(1)| - J,Lambda=1,Omega=1.5> = ((2*J+1)*sqrt(((2*J-1)*(2*J+3))/(2*J*(2*J+1)*(J+1))))/2
<J,Lambda=0,Omega=0.5 - |T(1)| + J,Lambda=1,Omega=1.5> = ((2*J+1)*sqrt(((2*J-1)*(2*J+3))/(2*J*(2*J+1)*(J+1))))/2
<J,Lambda=0,Omega=0.5 + |T(1)| - J,Lambda=1,Omega=0.5> = (2*J+1)^2/(2*sqrt(2*J*(2*J+1)*(J+1)))
<J,Lambda=0,Omega=0.5 - |T(1)| + J,Lambda=1,Omega=0.5> = ((-2*J-1)*(2*J+1))/(2*sqrt(2*J*(2*J+1)*(J+1)))
<J+1,Lambda=0,Omega=0.5 + |T(1)| + J,Lambda=1,Omega=1.5> = sqrt(((2*J-1)*(2*J+1))/(2*(J+1)))/2
<J+1,Lambda=0,Omega=0.5 - |T(1)| - J,Lambda=1,Omega=1.5> = sqrt(((2*J-1)*(2*J+1))/(2*(J+1)))/2
<J+1,Lambda=0,Omega=0.5 + |T(1)| + J,Lambda=1,Omega=0.5> = sqrt(((2*J+1)*(2*J+3))/(2*(J+1)))/2
<J+1,Lambda=0,Omega=0.5 - |T(1)| - J,Lambda=1,Omega=0.5> = (-sqrt(((2*J+1)*(2*J+3))/(2*(J+1))))/2

As before,the + and - in the bra and ket indicate the sign of the combination. The selection rules are + to - if J does not change, and + to + and - to - if J changes. PGOPHER also gives the line strengths, which are simply the squares of the matrix elements. As before the

The numerical values of the above can be displayed with "View, Transitions" and selecting the required upper and lower state J and parity, or looking in the log window with log window with the print level on the right set to "Matrices". The matrix elements between the two J= 4.5 states shown above are:

Rank 1 Transition Matrix for  <X, J = 4.5, Symmetry = f (-)| |A, J = 4.5, Symmetry = e (+)>
                 |A v=0 4.5 + 0.5>
<X v=0 4.5 - 1.5| 1.2712834523275
<X v=0 4.5 - 0.5| 1.2974982402692

These are between the symmetry adapted basis states, but as we saw above the X state energy levels are mixtures of these basis states, so the above matrix need to be transformed with these. The above matrix becomes:

Transformed Rank 1 Transition Matrix for  <X, J = 4.5, Symmetry = f (-)| |A, J = 4.5, Symmetry = e (+)>

                        |A v=0 4.5  4 F1e>
      <X v=0 4.5  4 F1f| -1.67509286382788 
      <X v=0 4.5  5 F2f|  0.702657240207629

Note the switch it notation - the first matrix was between basis states (labelled with J, phase and Ω) and the second is between the final energy levels, labelled with J, N and F1/F2. The squares of the numbers above are n ow proportional to transition intensity; to calculate absolute values the population and some units conversion must be done as described in Intensity Formulae. For the states above the following linelist entries are generated:

Mo M'   J' S'#'M"   J" S"#"    Position  I               Eupper        Elower S        A        Width
OH A  4.5  e 1 X  4.5  f 1    32404.0563 .0124568    32741.6000      337.5437 2.805936 8982545  0.000000 : qQ1(4.5) : A v=0 4.5  4 F1e - X v=0 4.5  4 F1f
OH A  4.5  e 1 X  4.5  f 2    32151.2517 .0006520    32741.6000      590.3483 .4937272 1543847  0.000000 : pQ12(4.5) : A v=0 4.5  4 F1e - X v=0 4.5  5 F2f

Note #' labels the upper state eigenvalue number, here always 1 as a 1x1 matrix was used for the A state. #" is the lower state eigenvalue number, here 1 or 2 to indicate the which of the two states in the matrix was used.