Ab Initio Study of the Large Amplitude Motions of Various Monosubstituted Isotopologues of Methylamine (CH3-NH2)

CCSD(T)-F12 theory is applied to determine electronic ground state spectroscopic parameters of various isotopologues of methylamine (CH3-NH2) containing cosmological abundant elements, such as D, 13C and 15N. Special attention is given to the far infrared region. The studied isotopologues can be classified in the G12, G6 and G4 molecular symmetry groups. The rotational and centrifugal distortion constants and the anharmonic fundamentals are determined using second order perturbation theory. Fermi displacements of the vibrational bands are predicted. The low vibrational energy levels corresponding to the large amplitude motions are determine variationally using a flexible three-dimensional model depending on the NH2 bending and wagging and the CH3 torsional coordinates. The model has been defined assuming that, in the amine group, the bending and the wagging modes interact strongly. The vibrational levels split into six components corresponding to the six minima of the potential energy surface. The accuracy of the kinetic energy parameters has an important effect on the energies. Strong interactions among the large amplitude motions are observed. Isotopic effects are relevant for the deuterated species.


INTRODUCTION
Methylamine (CH 3 -NH 2 ) plays important roles in the gas phase chemistry in the terrestrial and extraterrestrial atmospheres. The presence in the Earth's atmosphere has both natural and anthropogenic causes (Ge et al., 2011). In air quality studies, it is considered to be a Volatile Organic Compound (VOC) that can be a precursor of secondary organic aerosols (SOA) in the presence of glyoxal (De Haan et al., 2009). In 1974, it was detected in the interstellar medium and it is contemplated as a relatively abundant species (Kaifu et al., 1974) (Fourikis et al., 1974). Recent studies consider it a precursor of glycine and a building block of life (Ohoshi et al., 2019). Recently, methylamine has been detected in the quasar PKS 1830-211 (Muller et al., 2011) and together with other simple N-bearing species, it has been observed in the hot cores NGC 6334I MM1-3 (Bøgelund et al., 2019). Fourikis et al. (1977) have reported the probable detection of deuterated methylamine (CH 3 NHD) in Sgr B2.
The aim of the present work is the theoretical study of probably detectable methylamine isotopologues. Monosubstituted isotopologues were detected for many astrophysical molecules such as dimethyl-ether and methyl-formate as it is described in the references provided by the papers of Fernández et al. (2019) and Gámez et al. (2019). In a recent study of the methylamine main isotopologue, highly correlated ab initio methods were employed to simulate the far infrared spectra (Senent 2018). The low-lying vibrational energy levels in and their tunneling splitting components were computed, providing relevant information for rotational spectrum assignments, which are mandatory for the detection using radio-astronomy. Very accurate results were obtained by comparing with previous experimental data. A detailed review of previous theoretical and experimental works can be found in Senent (2018).
Previous studies devoted to n-methyl amines describe theoretical techniques and symmetry concepts useful for the present work   (Smeyers et al., , 1998 (Senent 2018)On the basis of previous ab initio results (Senent 2018), performed using explicitly correlated coupled cluster theory, CCSD(T)-F12 (Adlet et al., 2007) (Knizia et al., 2009), in this new paper, we attend to several monosubstituted isotopologues containing abundant cosmological elements. Although, to our knowledge, a unique isotopologue CH 3 NHD has been probably detected (Fourikis et al., 1977), other species are considered to be detectable species. Four isotopic species, 13 CH 3 NH 2 , CH 3 15 NH 2 , CH 3 NHD, and CH 2 DNH 2 , are studied and compared with the main isotopologue for predicting theoretically isotopic shifts. Recently, interstellar amines and their fragments have been studied using quantum-chemical computations   ).
An earliest CCSD(T)-F12 three-dimensional potential energy surface is revisited in the present work (Senent 2018) because it is mass independent. It is employed for constructing mass dependent effective potential energy surfaces for the different isotopologues. The surfaces present six minima separated by relatively low potential energy barriers. If the minimum interconversion is taken into consideration, the most abundant isotopologue can be classified in the G 12 molecular symmetry group . The isotopic substitutions carry out changes in the symmetry. Details concerning the followed procedure can be found in our previous paper devoted to the acetone isotopologues (Dalbouha et al., 2021). The effective surfaces allow to construct Hamiltonians depending on three interacting coordinates, two interacting large amplitude motions, the NH 2 wagging and the CH 3 torsion, and the HNH bending. Then, both the bending and wagging of the amine group are treated together. The final levels are computed variationally.

Electronic Structure Calculations
The theoretical study of methylamine isotopologues was started from the results of a previous work devoted to the main isotopologue CH 3 NH 2 (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 Symmetry of the Isotopologues
The main isotopologue, as well as 13 CH 3 NH 2 and CH 3 15 NH 2 , can be classified in the G 12 molecular symmetry group (MSG)  and in the C s point group. However, the H →D substitution carries out changes in the symmetry properties. CH 3 NDH must be classified in the C 1 point group and in the G 6 MSG, due to the absence of the symmetry plane. In CDH 2 NH 2, the D atom can replace the in-plane H atom (C s -CDH 2 NH 2 ) or one out-of plane H atom (C 1 -CDH 2 NH 2 ). 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 1 and C s are used. However, if the internal rotation is taken into account, C s -CDH 2 NH 2 and C 1 -CDH 2 NH 2 represent different minima of the same potential energy surface and they can be inter-converted. Then, both are classified in the same G 4 MSG.
The G 12 MSG contains six irreducible representations, four non-degenerate, A 1 , A 2 , B 1 and B 2 , and two double-degenerate E 1 and E 2 . The G 6 MSG contains three irreducible representations, two non-degenerate, A 1 , and A 2 , and one double-degenerate E. The G 4 MSG contains four non-degenerate irreducible representations, A 1 , A 2 , B 1 , and B 2 .

Rovibrational Parameters
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: 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 i 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: 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(Dalbouha et al., , 2021. In Table 1, the computed rotational constants of CH 3 NH 2 , 13 CH 3 NH 2, and CH 3 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 2 NH 2 than for C 1 -CDH 2 NH 2 .
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 3 15 NH 2 (Hirakawa et al., 1972). Deviation for several modes are significant, whereas the isotopic shits computed at the MP2 level of theory are reliable. Frontiers in Chemistry | www.frontiersin.org September 2021 | Volume 9 | Article 751203 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 ν 3 fundamental (CH 3 st), that interacts with two overtones (2ν 6 and 2ν 12 ). The NH 2 bending fundamental is predicted to interact strongly with the NH 2 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.

The far Infrared Spectrum
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 2 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: 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 2 bending and wagging vibrational modes were predicted. This fact suggests the prerequisite of a 3D-model. In Eq. 3, B qiqj and V eff represent the kinetic energy parameters and the effective potential defined as the sum of three contributions: Here, V(β, α, θ ) is the mass independent ab initio threedimensional 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 − γ θ (H5C4N1X + H6C4N1X + H7C4N1X − 2Π)/3 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: This analytical expression transforms as the totally symmetric representation of the G 12 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, 13 CH 3 NH 2 and CH 3 15 NH 2 . 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 3 NDH (G 6 ): FIGURE 2 | V eff (α, β; θ 270°) two dimensional potential energy surface (in cm −1 ).
Frontiers in Chemistry | www.frontiersin.org September 2021 | Volume 9 | Article 751203 6 and for CDH 2 NH 2 (G 4 ) To construct the effective potential using Eq. 4, two massdependent 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 2 NH 2 , 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. 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 3 NDH and CDH 2 NH 2 , respectively. For all the symmetries, the diagonal terms B ββ , B αα , and B θθ transform as the totally symmetric representation A 1 . However, the symmetry properties of the off-diagonal elements vary with the MSG: B αθ transforms as B 1 (G 12 , G 4 ) and A 1 (G 6 ) B αβ transforms as B 2 (G 12 , G 4 ) and A 2 (G 6 ) B θβ transforms as A 2 (G 12 , G 4 , G 6 ) The non-zero coefficients A 000 (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(Ohashi et al., , 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(Ohashi et al., , 1992 (Kreglewski 1993). Isotopic shifts of all the potential parameters are only important for the deuterated forms.
Frontiers in Chemistry | www.frontiersin.org September 2021 | Volume 9 | Article 751203  Table 6 and they are classified using symmetry and the v 7 , v 9 and v 15 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 12 species, the levels split into two non-degenerate and two doubledegenerated sublevels. The components of the ground vibrational state were computed to lie at 0.000 (A 1 ), 0.163 (B 2 ), 0.325 (E 1 ), and 0.407 (E 2 ) cm −1 . Very small shifts are found for 13 CH 3 -NH 2 , whereas for CH 3 -15 NH 2 , the subcomponents are close in energy (0.000 (A 1 ), 0.153 (B 2 ), 0.322 (E 1 ), and 0.398 (E 2 )). The non-degenerated components B 1 and A 2 of the ʋ 15 fundamental (0 0 1) were obtained to lie at 265.572 and 266.117 cm −1 in the main isotopologue and at 264.441 and 264.995 cm−1 in 13 CH 3 -NH 2 , and at 264.428 and 264.936 cm −1 in CH 3 -15 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 13 CH 3 -NH 2 , and at 763.413 and 769.602 cm−1 in CH 3 -15 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 3 -NDH, the nondegenerate components of the ʋ 9 and ʋ 15 fundamentals have been computed to be 236.260 and 236.470 cm −1 , and to be 715.178 and 718.848 cm −1 . 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 2 NH 2 levels. The ground vibrational state splits into two A 1 , two B 2 , one B 1 and one A 2 components lying in the 0.000-1.672 cm −1 range.

CONCLUSION
This work describes the shifts of spectroscopic parameters and the symmetry changes due to the isotopic substitutions for various probably detectable methylamine isotopologues, 13 CH 3 NH 2 , CH 3 15 NH 2 , CH 3 NHD, and CDH 2 ND 2 . A variational procedure and VPT2 theory are employed for describing rovibrational properties with a special attention to the far infrared region. For all the isotopologues, the levels up to 1,500 cm −1 over the ground vibrational state are determine variationally and classified using the G 12 , G 6 and G 4 MSG properties. For the main isotopologue, the ground vibrational state splits into six components computed to lie at 0.000 (A 1 ), 0.163 (B 2 ), 0.325 (E 1 ), and 0.407 (E 2 ) cm −1 . Very small differences are found for 13 CH 3 -NH 2 , whereas for CH 3 -15 NH 2 , the computed subcomponents are close in energy (0.000 (A 1 ), 0.153 (B 2 ), 0.322 (E 1 ), and 0.398 (E 2 )). Isotopic shifts are relevant for the deuterated forms, whereas the effects of substitution of heavy atoms are less relevant for the torsional excitation than for inversion excitations. Small variations of the kinetic energy parameters carry out substantial displacements of the levels. It can be pointed out that their computations in the deuterated forms is not straightforward.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

ACKNOWLEDGMENTS
The author acknowledges the CTI (CSIC) and CESGA for computing facilities.