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:

# 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


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

     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?

  • 1
    $\begingroup$ Just a few hints (and I don't know if any of that is actually helpful to the problem at hand): Specify #P to get more output, there might be more stuff hidden between the defaults. Use checkpoint files (%Chk). Split the calculation into one optimisation, and a frequency calculation (you can use %OldChk). I have made the experience that specifying keywords twice may cause problems for g16, so combine the option stacks. I try look into this deeper, but i can't promise $\endgroup$ – Martin - マーチン Mar 14 '19 at 17:31
  • $\begingroup$ I tried running a separate file, but that seems to mess up the optimized anharmonic frequencies without fixing the missing value for $\alpha_e$. #P doesn't seem to give me any additional relevant information. I did go ahead and remove the freq=ReadIsotopes and just specified directly in the atom specifications H(Iso=2). $\endgroup$ – Dustin Wheeler Mar 14 '19 at 21:55


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.


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.

| improve this answer | |
  • $\begingroup$ Thanks for the clarification, this was a concept with which I was unfamiliar. Nice to eliminate one red herring. Any insight into the rot-vib coupling constant, $\alpha_e$? $\endgroup$ – Dustin Wheeler Mar 26 '19 at 14:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.