HCl rota-vibrational lab: ab initio simulation of rotation constants

Question

I'm trying to set up the simulation portion of the classic undergraduate lab studying the rotavibrational spectrum of HCl and DCl. In terms of simulation (using Gaussian 16), the equilibrium geometry and vibrational energy calculations are trivial. Including anharmonicity is also simple through the use of the freq=anharm keyword. I'm running into trouble with the inclusion of vibrational-rotational parameters, which the program is supposed to support, but I can't find appropriate output information. I'm either missing the output data because I don't know what I'm looking for, or something is wrong in the calculation. Specifically, I'd like the students to get ab initio values for $D_e$, $\widetilde{\nu}_0$, $B_e$, and $\alpha_e$ in the following equation ($m$ is a substitution for $-J''$ or $J''+1$ for the $P$ and $R$ branches, respectively): $$\widetilde{\nu}(m) = \widetilde{\nu}_0 + (2 B_e - 2 \alpha_e) m - \alpha_e m^2 - 4 D_e m^3$$ I'm currently able to find $D_e$ and $\widetilde{\nu}_0$.

The only clue I can find in the calculation is the following lines:

WARNING: Anharmonic treatment of linear tops is experimental. 
         Moreover, an hybrid treatment is used to simulate spectra:
         - Energy: equations including degenerate modes are used.
         - Intensity: summation done on N' modes, considering only one mode
           per couple of degenerate modes. No variational correction done. 

    ==================================================
                  Coriolis Couplings
    ==================================================
[No Coriolis Coupling for diatomic molecules

Later, in the entry for Vibro-Rot alpha Matrix (where I'd expect to see entries for $\alpha_e$), I get:

Vibro-Rot alpha Matrix (in cm^-1)
---------------------------------
            A(z)        B(x)        C(y)
Q(    1)    -0.00000        NaN        NaN

I'm assuming the message for Coriolis Couplings is connected to the zero-valued alpha matrix. Is there a way to pull out these parameters in a Gaussian simulation? If not, is there another ab initio package that will perform these calculations? I've tried (without luck) to look through the NWChem documentation.

The input might be useful to reproduce the calculation. Unfortunately, the output is larger than SE.com will allow. Let me know if you require additional lines from it (or a full copy).

Input file:

%nproc=4
%mem=400MB
# MP2 def2tzvp
# opt freq=(anharm,vibrot,noraman) scf=tight freq=ReadIsotopes

Vibrational analysis for D35Cl

0   1
 Cl 0.0 0.0 0.0
 H  0.0 0.0 1.27

298.15  1.0
35
2

Update:

After more reading through the Gaussian references, I've found the rotational constants (Ba(x), Ca(y)) under the following section:

     ==================================================
          Vibrational Energies at Anharmonic Level
     ==================================================

 Units: Vibrational energies and rotational constants in cm^-1.
 NOTE: Transition energies are given with respect to the ground state.

 Reference Data
 --------------
                              E(harm)  E(anharm)     Ba(x)      Ca(y)
 Equilibrium Geometry                             10.643785  10.643785
 Ground State                 1525.714  1515.023        NaN        NaN

 Fundamental Bands
 -----------------
     Mode(n)      Status      E(harm)  E(anharm)     Ba(x)      Ca(y)
        1(1)      active      3051.429  2951.452        NaN        NaN

 Overtones
 ---------
     Mode(n)                  E(harm)  E(anharm)     Ba(x)      Ca(y)
        1(2)                  6102.858  5802.928        NaN        NaN

 WARNING: Anharmonic transition moments for symmetric and linear tops
          are not yet fully implemented.

However, I still can't find an entry for $\alpha_e$, the coupling constant between vibrational and rotational states. Is Gaussian capable of putting out this constant? If so, am I looking for it under the wrong name?

Answer 1

Disclaimer:

I never used Gaussian, so I am just guessing about the program.

Coriolis coupling:

Coriolis coupling couples bending modes with rotation and there is no bending in diatomic molecules.

There is no Coriolis coupling in Gaussian for HCl, because there can't be one. The relevant output could also be NaN, it is simply not defined.

Alpha "matrix":

For diatomic molecules you have two degenerate rotations and one stretching mode. Hence your alpha "matrix" is actually well represented by a scalar and used as such in your equations. For this reason one would expect only one alpha value.

That there is written -0.000 and not 0.000 is an indication for a truncated negative number. In addition other values are properly set to NaN. I think that your one reported alpha value is simply close to zero. You can try a verbose output to verify this.

Solution:

If you want to get a non zero alpha I would recommend a molecule with a shallower potential hyper surface, i.e. a lower force constant. Or you could use this example to discuss floating point arithmetics and truncation with your students.