The theoretical study of methylamine isotopologues was started from the results of a previous work devoted to the main isotopologue CH₃NH₂ (Senent 2018). In this earlier paper, the structural parameters of the minimum energy structure and a three-dimensional ab initio potential energy surface (3D-PES) were computed using explicitly correlated coupled cluster theory with single and double substitutions augmented by a perturbative treatment of triple excitations (CCSD(T)-F12b) (Adlet et al., 2007) (Knizia et al., 2009) using the MOLPRO package default options (Werner et al., 2012). The procedure was applied in connection with the AVTZ-F12 basis set, which contains the Dunning’s type aug-cc-pVTZ atomic orbitals (AVTZ) (Kendall et al., 1992) and the corresponding functions for the density fitting and the resolutions of the identity. These previous computed data are mass independent properties that can be used for the different isotopic species.

To determine the core-valence electron correlation effects on the rotational constants, the structure was optimized using CCSD(T) (coupled-cluster theory with single and double substitutions, augmented by a perturbative treatment of triple excitations) (Hampel et al., 1992) and the cc-pCVTZ basis set (CVTZ) (Woon and Dunning Jr 1995).

The full-dimensional anharmonic force field and the vibrational corrections of the potential energy surface are mass dependent properties that must be computed for each isotopologue. For this reason, new electronic structure calculations have been performed in the present work. That properties were determined using second order Möller-Plesset theory (MP2) (Møller & Plesset 1934) implemented in GAUSSIAN (Frisch et al., 2016). Anharmonic force fields allow obtain spectroscopic properties using second order perturbation theory (VPT2) (Barone 2005) (Bloino et al., 2012). The vibrationally corrected surfaces were employed to construct Hamiltonians for the isotopologues. The energy levels corresponding to the large amplitude vibrations and to the HNH bending mode were computed using a variational procedure implemented in ENEDIM (Senent 1998a; 1998b, 2001).

The main isotopologue, as well as ¹³CH₃NH₂ and CH₃ ¹⁵NH₂, can be classified in the G₁₂ molecular symmetry group (MSG) (Ohashi & Hougen 1987) and in the C_s point group. However, the H →D substitution carries out changes in the symmetry properties. CH₃NDH must be classified in the C₁ point group and in the G₆ MSG, due to the absence of the symmetry plane. In CDH₂NH_2, the D atom can replace the in-plane H atom (C_s-CDH₂NH₂) or one out-of plane H atom (C₁-CDH₂NH₂). If VPT2 is applied and a unique minimum is considered, the molecule is assumed to be semi-rigid and all the vibrations are described as small displacements around the equilibrium. Two different point groups C₁ and C_s are used. However, if the internal rotation is taken into account, C_s-CDH₂NH₂ and C₁-CDH₂NH₂ represent different minima of the same potential energy surface and they can be inter-converted. Then, both are classified in the same G₄ MSG.

The G₁₂ MSG contains six irreducible representations, four non-degenerate, A₁, A₂, B₁ and B₂, and two double-degenerate E₁ and E₂. The G₆ MSG contains three irreducible representations, two non-degenerate, A₁, and A₂, and one double-degenerate E. The G₄ MSG contains four non-degenerate irreducible representations, A₁, A₂, B₁, and B₂.

In the earlier paper (Senent 2018), the CCSD(T)-F12/AVTZ structural parameters of the methylamine equilibrium geometry, are detailed. The structure is shown in Figure 1, that helps to understand the atom labelling and the isotopic substitutions.

For all the isotopologues, the vibrational ground state rotational constants shown in Table 1, were computed from the CCSD(T)-F12 equilibrium rotational constants using the following equation: B 0 =  B e   ( CCSD ( T ) - F 12 / AVTZ - F 12 )   +   △ B e ( CCSD ( T ) / CVTZ ) core + △ B vib ( MP 2 / AVTZ ) (1)

Here, ΔB_e ^core collects the core-valence electron correlation effects on the equilibrium structure and ΔB^vib represents the vibrational contribution derived from the second order perturbation theory (VPT2) α _r ⁱ vibration-rotation interaction parameters. These last were determined using the MP2/AVTZ cubic force fields and vibrational second order perturbation theory. ΔB_e ^core was determined from the CCSD(T)/CVTZ parameters B_e (CV) and B_e(V), calculated correlating both core and valence electrons (CV) or just the valence electrons (V) in the post-SCF process. Then: Δ B e core   =  B e ( CV )   -  B e ( V ) (2)

This approximation has been corroborated in previous studies of other non-rigid molecules providing really accurate parameters, whose deviations with respect available experimental data, represent few MHz (Boussesi et al., 2016) (Dalbouha et al., 2016, 2021). In Table 1, the computed rotational constants of CH₃NH₂, ¹³CH₃NH_2, and CH₃NDH are compared with available experimental parameters (Ilyushin et al., 2005) (Motiyenko et al., 2016) (Ohashi et al., 1991). The MP2/AVTZ quartic centrifugal distortion constants corresponding to the asymmetrically reduced Hamiltonian, are shown in Table 2 where they are compared with previous experimental data (Ilyushin et al., 2005) (Motiyenko et al., 2016) (Ohashi et al., 1991). Disagreements between experimental and computed data can be correlated with the level of ab initio calculations used to compute the anharmonic force field. In addition, in methyl amine, the interaction between the internal and global rotation causes deviations. Isotopic shifts are more reliable.

The anharmonic fundamental frequencies shown in Table 3, were computed using VPT2 theory (Barone 2005) (Bloino et al., 2012) implemented in Gaussian (Frisch et al., 2016) and the MP2/AVTZ force fields. The modes are ordered following the criteria used for the main isotopologue that helps to make visible the isotopic shifts. Although VPT2 does not represent the proper treatment for the study of the vibrations responsible for the non-rigidity, it provides a good description of the mid- and near-infrared regions and a useful first description of the far-infrared region. In addition, it allows predict possible band displacements due to Fermi resonances. VPT2 theory ignores the inter-conversion of minima and treats the molecule as a semi-rigid species with a single minimum. If the existence of a single minimum is assumed, the resulting VPT2 properties are different for C_s-CDH₂NH₂ than for C₁-CDH₂NH₂.

The frequencies corresponding to the main isotopologue are compared with experimental data measured in the gas phase (Ohashi et al., 1989) (Kreglewski & Wlodarczak 1992) (Gulaczyk et al., 2017) (Hirakawa et al., 1972) [58]. Previous results are available for CH₃ ¹⁵NH₂ (Hirakawa et al., 1972). Deviation for several modes are significant, whereas the isotopic shits computed at the MP2 level of theory are reliable.

In Table 3, emphasized in bold, are the fundamental frequencies for which resonances can be relevant. Displacements due to the Fermi interactions were found to be relevant for the ν₃ fundamental (CH₃ st), that interacts with two overtones (2ν₆ and 2ν₁₂). The NH₂ bending fundamental is predicted to interact strongly with the NH₂ wagging overtone. Since both amine vibrations behave as inseparable modes, the variational model used for exploring the far infrared region, includes explicitly the bending coordinate.

As was assumed in the previous paper devoted to the main isotopologue (Senent 2018), the low-lying vibrational energy levels corresponding to the two large amplitude motions, the methyl torsion (θ) and the amine NH₂ wagging (α) can be determined by solving variationally a three-dimensional Hamiltonian where a third coordinate, the HNH bending angle (β), is considered to be an independent variable. The Hamiltonian obeys the formula: H ( β , α , θ ) = − ∑ i = 1 3 ∑ j = 1 3 ( ∂ ∂ q i ) B q i q j ( β , α , θ ) ( ∂ ∂ q j ) + V e f f ( β , α , θ ) (3)

This Hamiltonian was defined by taking into consideration the predictions of the test of resonances described in the previous section and in the previous paper (Senent 2018). Significant interactions between the NH₂ bending and wagging vibrational modes were predicted. This fact suggests the prerequisite of a 3D-model. In Eq. 3, B q i q j and V^eff represent the kinetic energy parameters and the effective potential defined as the sum of three contributions: V eff ( β , α , θ )    = V ( β , α , θ )   +  V ′ ( β , α , θ ) + V ZPVE ( β , α , θ ) (4)

Here, V( β , α , θ   ) is the mass independent ab initio three-dimensional potential energy surface; V’( β , α , θ   ) and V^ZPVE( β , α , θ   ) represent the Podolsky pseudopotential and the zero point vibrational energy correction (Dalbouha et al., 2021). The two last contributions must be computed for all the isotopologues because they are mass dependent properties. β, α, and θ, the HNH bending, the NH2 wagging and the torsional coordinates, are defined using curvilinear internal coordinates: β = HNH - HNH e α = 180 .0 - γ θ = ( H 5 C 4 N 1 X + H6C4N1X + H7C4N1X  - 2 Π ) / 3 (5)

HNH^e is the value of the HNH bending angle corresponding to the equilibrium geometry; γ represents the angle between the C-N bond and the HNH plane (see Figure 1); X denotes a ghost atom lying in the HNH plane perpendicular to the HNH angle bisector. The set of internal coordinates were chosen taking into consideration the procedure for the determination of the 3D-PES which demands a partial optimization of the geometry. Three internal coordinates, NHN, γ and H5C4N1X distinguish the selected conformations whereas twelve “dependent coordinates” are allowed to be relaxed in all the structures.

The ab initio three-dimensional potential energy surface, V( β , α , θ   ), was computed for the study of the main isotopologue (Senent 2018). It was constructed using the CCSD(T)-F12/AVTZ energies of 131 geometries defined for selected values of the independent coordinates that were fitted to the following series: V ( β , α , θ ) =   ∑ K , L , M A K M L C C β K ⁡ cos ⁡ L α ⁡ cos ⁡ 6 ⁡ M θ + A K M L S S β K ⁡ sin ⁡ L α ⁡ sin ⁡ 3 ( 2 M + 1 ) θ (6)

This analytical expression transforms as the totally symmetric representation of the G₁₂ MSG. Formally identical expressions can be employed for V’, V^ZPVE, V^eff( β , α , θ   ) and the diagonal kinetic energy parameters B_qiqi of the main isotopologue, ¹³CH₃NH₂ and CH₃ ¹⁵NH₂. However, since the H →D substitution carries out symmetry changes, the effective potential V^eff( β , α , θ   ) and the diagonal kinetic parameters must be expressed using less-symmetric analytical expressions. For CH₃NDH (G₆): V e f f ( β , α , θ ) =   ∑ K , L , M A K M L C C β K ⁡ cos ⁡ L α ⁡ cos ⁡ 3 ⁡ M θ + A K M L S S β K ⁡ sin ⁡ L α ⁡ sin ⁡ 3 ⁡ M θ (7) and for CDH₂NH₂ (G₄) V e f f ( β , α , θ ) =   ∑ K , L , M A K M L C C β K ⁡ cos ⁡ L α ⁡ cos ⁡ 2 ⁡ M θ + A K M L S S β K ⁡ sin ⁡ L α ⁡ sin ( 2 M + 1 ) θ (8)

To construct the effective potential using Eq. 4, two mass-dependent properties V′ and V^ZPVE must be computed for all the isotopologues and for all the geometries. The V’ pseudopotential is very small. However, V^ZPVE has important effects on the levels. It was determined within the harmonic approximation at the MP2/AVTZ level of theory. To obtain the mass-dependent properties of the low-symmetry varieties, more than 131 geometries and more than 131 sets of harmonic frequencies need to be computed. For example, in the case of CDH₂NH₂, 131x3 geometries are required because the three hydrogen atoms of the methyl group are not identical.

The ground vibrational state potential energy surface contains six equivalent minima corresponding to a single conformer. The contours of Figures 2, 3 represents layers of the 3D-surface of the main isotopologue containing the minimum energy structure. Figure 2 corresponds to V^eff (α, β; θ = 270°) and Figure 3 to V^eff (α, θ; β = 106°). Figures emphasize the coupling between coordinates.

The kinetic energy parameters were also computed for all the selected geometries and for all the isotopologues. The number of selected geometries required for their computation in the deuterated forms was 171 and 393 for CH₃NDH and CDH₂NH₂, respectively. For all the symmetries, the diagonal terms B_ββ, B_αα, and B_θθ transform as the totally symmetric representation A₁. However, the symmetry properties of the off-diagonal elements vary with the MSG:

B_αθ transforms as B₁ (G₁₂, G₄) and A₁(G₆)

B_αβ transforms as B₂ (G₁₂, G₄) and A₂(G₆)

B_θβ transforms as A₂ (G₁₂, G₄, G₆)

The non-zero coefficients A₀₀₀(B_qiqz) of the kinetic energy expressions are shown in Table 4. For the main isotopologue, they are compared with previous data (Ohashi et al., 1988, 1992), although in works based in experiments, these coefficients are considered to be constants. The potential energy barriers, V^tor and V^inv were estimated using the effective potentials. For the main isotopologue, they are in reasonable good agreement with previous data (Ohashi et al., 1988, 1992) (Kreglewski 1993). Isotopic shifts of all the potential parameters are only important for the deuterated forms.

Symmetry adapted series were employed as trial functions for the variational calculations. Products of harmonic oscillator solutions X_K (for the bending coordinate) and double Fourier series (for the wagging and torsional coordinates) were employed. Table 5 shows the symmetry eigenvectors. The convergence of the low energy levels requires long basis sets leading to Hamiltonian matrices of 18,755 x18,755 elements. In the case of the G₁₂ species, the matrices factorize by symmetry into eight blocks which dimensions are 1815 (A₁, B₂), 1,518 (B₁), 1,507 (A₂), and 3,025 (E_1x, E_1y, E_2x and E_2y). For the G₆ species, the corresponding submatrix dimensions were 3,333 (A₁), 3,322 (A₂), and 6,050 (E), whereas for the G₄ species, the dimensions were 4,840 (A₁, B₂), 4,543 (B₁), and 4,532 (A₂).

The resulting energy levels are shown in Table 6 and they are classified using symmetry and the v₇, v₉ and v₁₅ quantum numbers. For the main isotopologue, the energies are compared with those of Kreglewski (1989) obtained from experimental parameters. The computed levels denote a slight improvement with respect to the work of Senent (2018), after using longer expansions for the kinetic energy parameters. The aim was to increase precision considering that isotopic shifts are relatively small. We observed that the vibrational energies are very sensitive to the kinetic contributions. It can be pointed out that their computations in the deuterated forms is not straightforward.

Each energy level splits into six components corresponding to the six minima of the potential energy surface. Their distributions are represented in Figure 4. In the G₁₂ species, the levels split into two non-degenerate and two double-degenerated sublevels. The components of the ground vibrational state were computed to lie at 0.000 (A₁), 0.163 (B₂), 0.325 (E₁), and 0.407 (E₂) cm⁻¹. Very small shifts are found for ¹³CH₃-NH₂, whereas for CH₃-¹⁵NH₂, the subcomponents are close in energy (0.000 (A₁), 0.153 (B₂), 0.322 (E₁), and 0.398 (E₂)). The non-degenerated components B₁ and A₂ of the ʋ₁₅ fundamental (0 0 1) were obtained to lie at 265.572 and 266.117 cm⁻¹ in the main isotopologue and at 264.441 and 264.995 cm−1 in ¹³CH₃-NH₂, and at 264.428 and 264.936 cm⁻¹ in CH₃-¹⁵NH_2. For ʋ_9, the corresponding components of the (0 1 0) level were obtained to lie at 771.083 and 775.011 in the main isotopologue and at 776.413 and 777.617 cm−1 in ¹³CH₃-NH₂, and at 763.413 and 769.602 cm−1 in CH₃-¹⁵NH_2. It may be concluded that the effects of isotopic substitutions on the heavy atoms are less relevant for the torsional excitation than for inversion excitations.

As was expected, isotopic effects on the low-lying energies are more noticeable for the deuterated species. For CH₃-NDH, the nondegenerate components of the ʋ₉ and ʋ₁₅ fundamentals have been computed to be 236.260 and 236.470 cm⁻¹, and to be 715.178 and 718.848 cm⁻¹. The gaps among subcomponents of the ground vibrational state are smaller than in the hydrogenated species. The isotopic substitution in one methyl group hydrogen breaks ten the degeneracy of the CDH₂NH₂ levels. The ground vibrational state splits into two A₁, two B₂, one B₁ and one A₂ components lying in the 0.000–1.672 cm⁻¹ range.