Largely Reduced Grid Densities in a Vibrational Self-Consistent Field Treatment Do Not Significantly Impact the Resulting Wavenumbers

Especially for larger molecules relevant to life sciences, vibrational self-consistent field (VSCF) calculations can become unmanageably demanding even when only first and second order potential coupling terms are considered. This paper investigates to what extent the grid density of the VSCF’s underlying potential energy surface can be reduced without sacrificing accuracy of the resulting wavenumbers. Including single-mode and pair contributions, a reduction to eight points per mode did not introduce a significant deviation but improved the computational efficiency by a factor of four. A mean unsigned deviation of 1.3% from the experiment could be maintained for the fifteen molecules under investigation and the approach was found to be applicable to rigid, semi-rigid and soft vibrational problems likewise. Deprotonated phosphoserine, stabilized by two intramolecular hydrogen bonds, was investigated as an exemplary application.


Introduction
The continuously increasing availability of computational resources as well as the development of accurate and efficient quantum chemical approaches have made computational vibrational spectroscopy an indispensable field complementing experimental techniques. Nowadays, almost every quantum chemical software package enables the analysis of second-order properties of the energy (e.g., IR and Raman absorptions, IR intensities, electric dipole polarizabilities, nuclear magnetic resonance chemical shifts, spin-spin coupling constants). Arguably, the prediction of IR and Raman absorptions is among the most important applications considering the prevalence of these technique. While most fundamental IR absorptions can be assigned quite satisfyingly on an empirical basis or via normal coordinate analysis, especially the near-infrared region is cluttered by overtones and combination excitations that are cumbersome to assign.
With regard to computational approaches, the harmonic oscillator approximation (HOA) is the most fundamental technique for obtaining vibrational spectroscopic data. However, due to the rigorous assumption introduced, in that the bond potential exhibits a harmonic shape, significant deviation from experiment is frequently observed. The most simple solution accounting for the lack of anharmonicities is to introduce empirical scaling factors that are multiplied with the harmonically approximated absorptions. Scott and Radom derived scaling factors for a vast number of ab initio and semi-empirical methods and a large number of basis sets [1,2]. It has to be stressed, however, that an empirical scaling may not be applied on a system-independent basis, even though the scaling factors have been derived for a rather large set of molecules. Noteworthy, a harmonic bond potential cannot be assumed a proper basis for considerations towards excitations involving more than one quantum of energy.
Especially due to these fundamental deficiencies, further corrective techniques have been developed that account for anharmonicities in an explicit manner. The two most prominent approaches are the vibrational self-consistent field (VSCF) method [3][4][5][6][7][8] and the vibrational second order perturbative ansatz (VPT2) [9,10]. VPT2 is a technique that relies on the computation of higher-order derivatives of the energy. Gaussian [11] is probably the most prominent commercial software package incorporating a VPT2 algorithm which relies on third and semi-diagonal fourth derivatives of the energy with respect to the nuclear coordinates. This method is applied regularly [12,13] but it may only be employed safely to vibrational problems where the harmonic part of the potential is dominant. Moreover, the computationally demanding generation of third and fourth derivatives limits the approach to rather small systems. VSCF on the other hand utilizes a separability ansatz which enables the anharmonicity of each mode to be accounted for by screening the respective potential in a point-wise manner. Mode interactions may be evaluated by computing d dimensional grids and the maximum possible dimensionality is determined by the number of vibrational degrees of freedom (N ) featured by the molecule under investigation. In a conventional VSCF calculation, each mode's potential is characterized by 16 grid points [14]. In the past, scientific papers have been reported that employed grids of higher (r = 32) [15] and lower resolutions (r = 8) [16][17][18]. Roy et al. [19] recently presented data based on a variety of grid densities ranging from 8 to 16 points. Since they reported overall deviations summed over a set of reference molecules, one could conclude that this rather diverse choice of grid resolutions does not impose a significant error on the wavenumbers. However, to date there still seem to exist ambiguities with regard to a generally applicable grid density that is computationally feasible. Hence, and considering the massive computational effort involved in high resolution VSCF calculations, it seemed promising to conduct a systematic study addressing these issues. 15 reference molecules are investigated in detail and an application to the lowest energy conformer of phosphoserine is presented. At this point it seems worth mentioning that while many quantum chemical software packages such as GAMESS [20], NWChem [21] and MOLPRO [22] incorporate VSCF implementations, GAMESS is the program employed in this work due to the fact that it is freely available and compatible with the most prominent computer operating systems.

Methods
In this section, a brief overview describing the procedures underlying a VSCF evaluation is given. Particularities of the involved techniques are presented while reference to the original literature is given for further details.

Energy Minimization
A prerequisite for every computational spectroscopic analysis is a proper structural ensemble. For the majority of cases, the user would want to obtain absorption data corresponding to an equilibrium structure and thus, an energetically favorable reference geometry is required. This implies that at least a local energy minimum (or the global energy minimum) is obtained and this is reflected by a nearly zero gradient of the energy g(q)∆q with respect to the nuclear coordinates q. While the program's default criteria suffice for energetical and structural considerations in many cases, a VSCF calculation requires a thoroughly minimized geometry. This condition is realized by setting up rigorous cutoff tolerances for g(q)∆q. We considered a molecular geometry properly minimized when each gradient contribution is smaller than 0.000001 E h · Bohr −1 . Every calculation reported herein was performed at the Møller-Plesset level of theory, accounting for electron correlation effects via a perturbative ansatz [23][24][25]. The frozen-core approximation [26] was employed, explicitly correlating all but the core electrons. Dunning's correlation-consistent polarized basis set of triple-ζ quality has been utilized since it can be considered a routinely employed set of functions that is known to deliver results of acceptable accuracy [27]. The appropriate point group was imposed onto the molecular geometry where applicable in order to make use of the efficient symmetry optimized SCF algorithm in GAMESS.

Second Derivatives of the Energy
A potential V (q) may be expressed as a Taylor series and for the case that the reference geometry resembles an energy minimum and assuming a harmonic potential shape, all but the first and third terms vanish. This simplified approach enables the evaluation of the elements of H r,s by analytical or numerical analysis: A transformation into mass-weighted coordinates and diagonalization of the Hessian matrix H r,s yields the eigenvectors corresponding to the normal modes. At this point, computational spectroscopic data at the HOA level is available. It seems worth mentioning that the further description of extended vibrational analysis is based on mass-weighted normal coordinates or other non-redundant coordinates for the sake of clarity, whereas the expansion of the potential energy surfaces was realized in rectilinear coordinates (Section 3.1). VSCF calculations require a Hessian corresponding to the equilibrium geometry, as will be discussed in the next paragraph. The Hessian matrices in this work have been obtained semi-numerically involving two displacements of ±0.01 Bohr about the reference geometry for each atom.

The VSCF Routine
Even though VSCF accounts for anharmonicities explicitly, it still requires H r,s as a reference state. This is owed to the fact, that VSCF assumes a molecule's vibrational wave function Ψ(Q 1,..,N ) to be separable into single mode wave functions ψ (n) i : While this product ansatz can be implemented rather efficiently, it inherently limits the accuracy of the VSCF technique since it represents decoupled vibrational modes.
There are, besides the methods mentioned in Section 2.4, techniques available which are able to circumvent this limitation [28,29]. Introducing the variational principle [30,31], the single mode VSCF equation [6,32] is formulated as: where V (n) i (Q i ) is mode Q i 's effective potential. Equation (4) neglects parts of the full Watson Hamiltonian for non-linear molecules [33] which accounts for N vibrational degrees of freedom: The second and third part of the Hamiltonian represent the generalized inverse of the effective moment of inertia µ αβ and the vibrational angular momentum operatorπ α . Given that the full Watson Hamiltonian is cumbersome to implement computationally, most implementations resort to neglecting the last two terms ofĤ. The error introduced is small especially for fundamental excitations as has been shown by Bowman et al. [34].
The VSCF equations are solved in an iterative manner until self-consistency is achieved. Illustratively, VSCF may be seen as an approach where each vibration is influenced by the mean field of the surrounding modes. However, especially when certain mode interactions exhibit unique properties, this mean field ansatz does no longer properly describe the vibrational problem sufficiently well. A number of post-VSCF techniques accounting for this deficiency have been presented and they are discussed in the subsequent section.
Since the VSCF technique involves a grid-based screening of bond potentials and their interactions, a discussion of the so-called hierarchical expansion is justified. For a system exhibiting N normal modes, the VSCF expansion is denoted as: A complete interaction scheme is realized by imposing N -dimensionality onto Equation (6) but obviously, such a treatment scales unfavorably with the studied system size. Noteworthy, most applications of VSCF theory truncate the expansion largely by including only [35][36][37][38]. This simplification is valid for a large number of applications and the inclusion of higher-order interaction terms is hardly justifiable for other than the smallest molecules or for cases, where convergence is only achieved when including such terms. Programs such as MULTIMODE [34,39], Molpro [22] and MIDAScpp enable coupling orders of d > 3 yielding IR bands of exceptional accuracy [7,40]. In this article, we confine our discussion to contributions since this study aims at larger molecules where an incorporation of higher coupling terms would be prohibitive. The number of grid points due (N p ) is calculated via Equation (7) .
Assuming r = 16 and taking glycine as an example, the pair-wise approximation would require 71,040 points to be computed while a VSCF calculation involving also three-mode interactions already gives rise to 8,361,344 energy evaluations, a nearly 120 fold increase in computational burden. Herein, we will examine the impact of the grid density on the quality of the absorption data by conducting VSCF calculations involving between 6 and 16 grid points per mode. For the displacements underlying the PES scan, a symmetric grid range of [−4ω −0.5 i , +4ω −0.5 i ], with ω i being the harmonic frequency of the i-th normal mode, was chosen. The grid points within these boundaries have been set-up in an equidistant manner which implies that for even grid densities (i.e., 6, 8, 10, 12, 14 and 16 points), equilibrium is located between the two innermost grid points. For odd grid densities (i.e., 7, 9, 11, 13 and 15 points), equilibrium is described by one distinct grid point. For every resolution r < 16, the grid points are interpolated to the original resolution of r = 16 by a polynomial fit.

Extensions of the VSCF Approach
As mentioned earlier, the VSCF assumption does not hold for many applications. Therefore, a perturbative approach has been suggested, accounting for the error introduced by the mean-field VSCF ansatz [15,41,42]. A prerequisite is that the difference between the "true" energy and the mean-field VSCF energy is small. An expansion known from electronic structure theory [23] is then formulated for the energy and wave function of a state n that is truncated at second order: with H = H 0 n + λ∆V and an insertion of E n and Ψ n in HΨ n = E n Ψ n , the energy contribution at second order may be formulated as: Herein, Ψ 0 n and Ψ 0 m are the unperturbed vibrational product wave functions and E 0 n and E 0 m are the unperturbed VSCF energies. This technique, which is commonly abbreviated as second order perturbation theory augmented (PT2)-VSCF or correlation-corrected (CC)-VSCF, has proven useful for many applications. However, since obtaining correlation-corrected VSCF energies requires significant computational resources, the basic PT2-VSCF approach is limited to small systems. Due to the inherent need for an efficient solution of Equation (10), Gerber's work group developed a solution based on the assumption of orthogonal vibrational single mode wave functions (Equation (3)) which leads to the annihilation of diagonal elements in Equation (10) [43]. Tests involving Gly n peptides showed that the runtime improvement is greater for larger molecules: for monomeric Gly, a speedup factor close to 6 was observed while the correlation-corrected wavenumbers of tetraglycine have been evaluated more than 16 times faster as with the conventional PT2-VSCF ansatz. This acceleration technique, and also the fact that this accelerated PT2-VSCF technique is readily implemented in GAMESS [20], are probably the main reasons for PT2-VSCF being a routinely employed VSCF correction. It has to be stressed, however, that the evaluation of the VSCF equations with its further corrections may not be confused with the preceding and highly demanding evaluation of the potential energy grid (vide supra).
Problems may arise when degenerate vibrational states are present. PT2-VSCF can fail here due to a close to zero denominator in Equation (4) which can lead to a largely overestimated perturbative correction. Hence, the degenerate PT2-VSCF method has been developed [44] but to date it has been implemented in GAMESS exclusively for degeneracies arising from fundamental excitations. A more generally available simple solution available in the GAMESS program code is that contributions involving denominators falling below a critical value are excluded from the treatment. Most applications reporting PT2-VSCF derived data rely on this simplification and therefore, we will also confine our calculations to this simplified technique.
Due to the popularity of the PT2-VSCF method and the fact that results of good quality at manageable computational cost are available also for larger molecules, all data presented in this paper is corrected exclusively with this perturbative ansatz.

Results and Discussion
The fifteen molecules under investigation gave rise to 176 distinct vibrational degrees of freedom, considering all fundamental stretching, deformation and torsional vibrations. PT2-VSCF can yield questionable or sometimes even divergent results for very low-lying and floppy torsions, which is due to the fact that the PES expansion is conventionally carried out in Cartesian coordinates [18,19,54].
The error induced through an anharmonic correction of such a vibration can exceed the boundaries of accuracy known from the HOA [54]. Hence, the usual procedure is to either treat such vibrations harmonically or to describe the underlying displacements in internal coordinates. A substitution for physically more meaningful internal coordinates was proposed by Njegic and Gordon [54] and could be shown to yield good results for formamide and thioformamide [55] as well as for H 2 O 2 but by introducing a new expansion technique for the kinetic energy operator [56]. Nonetheless, setting up internal coordinates that properly describe the displacements underlying a VSCF treatment is by no means a trivial task. The GAMESS code is able to identify each normal mode's contributions to particular internal coordinates [57], but the user still has to input a balanced description of each vibrational degree of freedom which can become unmanageably difficult for larger systems. Hence, and especially when large molecules are investigated, the majority of users of VSCF theory resort to a PES expansion in Cartesian coordinates and the contributions of critical torsions are omitted when self-consistency is not achieved.
For glycine (C 2 H 5 NO 2 ), three normal modes (i.e., the N-C α -C carb -O torsion, the NH 2 group torsion and δ C carb O 2 ,oop ) had to be excluded from the VSCF treatment since they are known to lead to divergence during a perturbation theory corrected VSCF evaluation [19]. Similarly, for methanol (CH 3 OH) the CO axis torsion has been omitted from the PT2-VSCF treatment due to an inadequate description of this particularly floppy torsional mode within the VSCF framework [58]. Dimethylether also exhibits two floppy torsions involving the C-O axes that have been excluded likewise. For ethane, one normal mode near 290 cm −1 [59][60][61][62] was omitted due to its floppy character but the other six missing modes did arise from degenerate states in ν CH 3 ,as , δ CH 3 ,as and ρ CH 3 . Experimentally, these degeneracies cannot be distinguished and since the PT2-VSCF derived values did not exhibit significant numerical discrepancies, each pair of degeneracies is presented as a single mean value for the sake of visibility.

Performance of the PT2-VSCF Approach
The computationally obtained absorption data are compared to experimental data in Table 1. The column headers indicate the employed number of grid points during the VSCF evaluations. As a measure of quality, the mean absolute percentage error µ for each molecule and each grid density is calculated according to Equation (11): The obtained values for µ are generally smaller than 2%, which is in good agreement with the recent work by Roy et al. [19] who concluded that MP2/cc-pVTZ based PT2-VSCF evaluations exhibit a mean unsigned error of under 2%. MP2 is a topical ab initio method that accounts for electron correlation effects to a certain extent and it seems as if for most molecules, this method indeed delivers satisfactory results. Acetonitrile is an exception due to its CN triple bond. However, the particularities of this molecule are discussed elsewhere [63] and it was found that state-of-the-art quantum chemical methods [64][65][66][67][68][69] are required for a proper description of ν CN . Importantly, the more or less ubiquitous stretching motions arising from a methyl group are not described in a reliable manner and this has been recently ascribed to the nature of the MP2 method [19]. Conversely, it was shown that when higher order coupling terms (i.e., N i<j<k V triples i,j,k (Q i , Q j , Q k )) are included in a DPT2-VSCF [44] treatment, MP2 indeed seems to be a viable ab initio method of choice [63]. VSCF data relies on the HOA and hence makes use of the rigid rotor approximation [70] which inherently is co-responsible for discrepancies between experiment and theory. Considering that still a large number of approximations are involved even in state-of-the-art VSCF calculations and their underlying ab initio energy evaluations, it must be concluded that computationally derived spectroscopic data are prone to a certain degree of fortuitous error compensation. However, both the data presented by Roy et al. [19] and our results indicate that triple-ζ based MP2 calculations in a VSCF treatment rather reliably deliver data with no more than 2% of unsigned error.

Performance at Reduced Grid Densities
As a VSCF evaluation with 16 grid points per mode is computationally highly demanding, a reduction of the grid density is desirable. From the data in Table 1 it becomes evident that for most cases, a reduction to 10 grid points yields identical results as the corresponding calculation with 16 grid points, indicating that the computational demand in a VSCF treatment of contributions can be reduced safely by more than 50%. It has to be kept in mind, however, that for every set of reduced densities, the grid is interpolated to the original resolution and the displacement boundaries are not affected by the reduction. Hence, as long as the points on the potential energy surface are chosen properly, the accuracy of the VSCF results is not affected. Noteworthy, the default routine in GAMESS chooses the displacement boundaries in dependence to the underlying normal mode's wavenumber but the user is nonetheless able to alter these settings if desired. Table 2. Obtained values for µ in %, their arithmetic means and the corresponding t-values according to Equation (12). n.a. means "not applicable". Reviewing the values for µ in Table 2 permits the conclusion that major errors are starting to appear at grid densities <8 points. The GAMESS manual states that the VSCF code is "thought to give accuracy to 50 cm −1 for the larger fundamentals" when MP2 with a triple-ζ basis set is employed and hence, the boundaries of accuracy are wider than the error inflicted by a reduction of the grid density to 8 points. In order to provide statistical evidence as to what extent a reduction of the grid density is viable, paired t-tests [114] of the obtained values for µ have been conducted with the following equation:

Number of Grid Points
In Equation (12), µ ref is the reference value at a grid density of 16 points. µ calc is the mean absolute percentage deviation at the examined grid density. χ is the number of investigated molecules and the t-value is calculated with s d as the corrected sample standard deviation of the respective grid density. For a two-tailed problem, critical values for t can be specified according to Table 3. Except for the data obtained at grid densities of 6 and 7 points, no significant deviations from the reference data are observed, indicating that even a reduction to 8 grid points is a very viable option. The computational effort involved in a pair-approximated VSCF evaluation (Equation (7)) may thus be reduced by a factor of nearly 4. It has to be stressed, however, that the obtained t-value for r = 8 (i.e., 0.09) may not be confused as to indicate the objectively best agreement with the reference values since the statistical analysis is based on unsigned error values. Moreover, the anharmonic correction seems to be slightly over-estimated at r = 8, yielding more red-shifted absorptions, as can be learned from Table 1. This over-correction is responsible for a slightly improved agreement with experiment for most of the molecules, especially since it is known that MP2 tendentially yields blue-shifted IR absorptions [19]. Overall, it may be concluded that VSCF treatments involving at least 8 grid points reproduce the experimental values within the same boundaries of accuracy as the reference values with a mean absolute percentage deviation from experiment of ∼1.3% (Table 2).
In their recent paper, Roy et al., reported that soft and semi-rigid molecules have to be treated at higher grid densities than rigid molecules and they mention ethylene oxide as a particular example where the VSCF equations converge at no less than 12 grid points. Our results contradict this statement since our VSCF calculations of ethylene oxide did indeed converge for grid densities ≥6 points and a density of 8 points yielded the same quality of results as observed for the other molecules. Given the fact that the identical code, the same ab initio method and the same basis set were employed, this outcome is rather puzzling. This disagreement may be ascribed to the use of symmetry during our VSCF calculations and the underlying reference geometries having been minimized to the most stringent criterion available. Overall, none of our VSCF evaluations exhibited convergence problems apart from the few excluded torsional degrees of freedom that have been mentioned earlier.

Application to Deprotonated Phosphoserine
Phosphorylation of proteins is considered a major signal transduction mechanism, mainly occurring at the OH terminus of the amino acids serine, threonine and tyrosine. As a major constituent of biofluids, phosphoserine was subject of a gas-chromatographic investigation of human urine samples [115]. In autopsied Alzheimer's disease brain tissue, L-phosphoserine was found in elevated concentrations [116] and only very recently, plasma phosphoserine levels were found to be upregulated in sepsis patients [117]. Deprotonated phosphoserine, abbreviated as [pSer-H] − , was investigated via infrared multiple photon dissociation (IRMPD) spectroscopy and hence, experimental data of a few major gas phase IR absorptions of the molecule is available [118]. The authors identified the lowest energy conformer of [pSer-H] − as the one exhibiting hydrogen bonds between the carboxylic OH and the phosphate O as well as between the phosphate OH and the amino N. This conformer was chosen as a benchmark for the 8 point VSCF analysis due to the fact that it represents rigid and weak structural motifs likewise. Figure 1 shows the minimum geometry obtained at the MP2/cc-pVTZ level of theory which served as a reference state for the VSCF evaluation. Table 4 lists the PT2-VSCF computed absorptions of [pSer-H] − , the respective intensities and the available experimental data [118]. Due to the cage-like structure, many normal modes couple among each other and hence, especially lower vibrations cannot be ascribed to distinct normal modes. In such cases, the main vibrational contributions are given in Table 4. For illustrative purposes, the calculated IR spectrum of [pSer-H] − is shown in Figure 2. The intensities have been computed using harmonic dipole derivatives at the MP2/cc-pVTZ level of theory and band broadening was introduced via a Lorentzian function using a band width at half height of 20 cm −1 .   Table 4. Considering the 11 available experimental values, µ was computed for the corresponding VSCF data as 1.38%. Since we did calculate the IR absorptions of [pSer-H] − exclusively using 8 grid points per mode, no definitive conclusion can be drawn whether higher grid densities would enable an improved accuracy for this larger molecule. A re-evaluation of the potential energy surface at grid densities up to 16 points would be unmanageably expensive but nonetheless, the data presented fits well into the error boundaries as observed for the 15 molecules discussed earlier.
The band assignment in Table 4 clearly shows that empirical considerations are only applicable to the large fundamentals. For a proper description of low lying vibrations that are determined by a number of torsions, computational approaches are an important option aiding the spectroscopist. In their experimental work, Scuderi et al., have employed computational techniques during their band assignment of [pSer-H] − as well, but they resorted to scaled HOA absorption data at the B3LYP/6-311+G(d,p) level of theory.

Conclusions
VSCF theory has grown an important field in computational spectroscopy and this is owed to the increasingly available computational resources and to efficient and easily applicable algorithms. With this systematic study it has been demonstrated that a largely reduced grid density in a VSCF evaluation does not only spare considerable resources but also does not significantly affect the resulting absorption data. It was found that the convergence of VSCF equations is not impaired by such a reduced-effort technique even when highly problematic densities of <8 grid points are employed. While further investigations are required with regard to the role of symmetry in a VSCF treatment, it was found that reduced grid densities may be safely applied to a wide set of molecular entities. Application to [pSer-H] − showed that a 8 point VSCF calculation yields accurate data for a molecule exhibiting weak interactions. The routine presented herein should, in conjunction with other recent developments [35], prove valuable for larger molecules relevant to life sciences where a conventional VSCF treatment with 16 grid points per mode is not feasible.

Conflicts of Interest
The authors declare no conflict of interest.