The table 3 lists the computed and experimental anharmonic vibrational frequencies, expressed in cm -1 , for the H 2 35 Cl + and H 2 37 Cl + isotopologues. The differences between the experimental and the computed values are 0.5 cm -1 for the bending mode (ν 2 ) and about 5 cm -1 for the two stretching modes (ν 1 and ν 3 ). This shows a good agreement between our calculated values and the experimental ones. Furthermore, these values show a clear improvement with respect to those computed by Botschwina with his unscaled potential , in particular, for the bending mode. Indeed, Botschwina developed two potentials: (i) the first one (unscaled) is obtained after the fitting of the ab-initio calculations; (ii) the second one is derived from the first one, but some parameters are scaled to reproduce the experimental fundamental transitions. Therefore, we compare our full ab-initio potential with the first Botschwina potential, and our results clearly show the quality of our potential. Table gives the vibrational frequencies for the other isotopologues of this ion, namely HD 35 Cl + , HD 37 Cl + , D 2 35 Cl + , D 2 37 Cl + . The equilibrium molecular structure was therefore determined by taking the expectation values of each coordinate (⟨R 1 ⟩, ⟨R 2 ⟩ and ⟨Θ⟩) for the vibrational ground-state. The values resulting from these calculations are reported in the table 5. These values obtained for the isotopologues, H 2 35 Cl + , HD 35 Cl + and H 2 37 Cl + agree very well with the experimental values (the distances and the angle are extracted from the experimental rotational constants) . Indeed, for the three isotopologues, the angle values are very close to the experimental ones and the errors are less than 0.1 • . With respect to the Cl-H distance of H 2 35 Cl + obtained from ref , the difference with our calculation is about 0.1 Å, which is surprisingly large. However, we are confident with our values. Indeed, our equilibrium values agree very well with Botschwina ones and also from the experimental ones . Furthermore, with respect to a more recent study of Araki et al. , the experimental Cl-H and Cl-D distances (associated to the zero point average structure ) are in good agreement with our calculated values since the largest difference with respect to the experimental value is only 0.005 Å with respect to ref . Furthermore, our calculations nicely reproduce the shortening of the Cl-D bond length with respect to the Cl-H one for the HD 35 Cl + ion . In the Table , we also give the expectation values for other isotopologues, D 2 35 Cl + , HD 37 Cl + and D 2 37 Cl + . In addition, we were able to compute the three rotational constants associated with the vibrational ground state. However, this has to be done with care. Indeed, since we did not perform a full ro-vibrational calculation, the rotational Eckart conditions must be applied to minimize the Coriolis coupling. The Cartesian equilibrium geometry (obtained from R 1 = R 2 = 2.46729 bohr and θ = 1.64588 rad) is used as the reference geometry to apply the Eckart conditions. They have been implemented in Tnum as described in ref . The rotational constants have been computed with the Fortran codes ElVibRot coupled with Tnum using the following procedure: (i) Calculate the average values (⟨µ α,β ⟩, with α and β being x, y or z) of the contravariant rotational components of the metric tensor. This is done with the numerical conditions (basis set and grid parameters) used to compute the energy levels. (ii) Diagonalize the previous 3 × 3 matrix to get the rotational constants as its eigenvalues. The comparison between the computed and the experimental values given in Table shows a nice agreement. Indeed, the largest error is about 1. GHz. Furthermore, the rotational constant shifts due to the isotopic substitution (from H 2 35 Cl + to H 2 37 Cl + ) are well reproduced. In this computational work, we have generated the 3D-potential energy surface of the ground state of the chloronium ion using a state-of-the-art electronic structure level, namely the F12-explicitly correlated coupled cluster method with a F12-triple zeta atomic basis set. This fitted polynomial expansion was incorporated into a home-made variational Fortran code to compute the energy levels, the zero-point average coordinate values and the rotational constants of several isotopologues, H