Predicting Vibrational Spectroscopy for Flexible Molecules and Molecules with Non‐Idle Environments

Flexible molecules and non‐idle environments can present severe problems in the prediction of vibrational spectra. This work focuses on infrared and vibrational circular dichroism for which the latter is extremely sensitive to such complicated situations. Here two possible ways to perform spectroscopy calculations for such situations are investigated. The first approach is based on static quantum chemical calculations employing cluster‐weighting. The second approach is rooted in ab initio molecular dynamics simulations using the time correlations approach. For the present example ((R)‐butan‐2‐ol), excellent spectra from simulations are obtained for gas and bulk, when the former is averaged over trajectories with all possible starting conformers and these are scaled to match the experimental spectrum. The cluster‐weighted approach is inferior to the simulations but still reaches very good results at much less computational cost. A simplified, but computationally less expensive simulation approach is considered by approximating the electric and magnetic moments, which are input quantities in the correlation function, by classical equations and different populations analysis for the required partial charges. The results are inferior to the full simulations but still give satisfying results at the advantage of being much faster to calculate.


Introduction
The calculation of bulk phase spectra [1] often requires the accurate treatment of the environment. This is especially the case if the environment behaves non-idle towards the solute or observed molecules. Consequently, two questions arise, first what is considered an "accurate treatment" and second what counts as a "non-idle environment"? As non-idle environment we simply refer to such an environment that directly or indirectly changes the properties (here the spectrum) of the solute or the molecule under investigation. Therefore, it needs to be considered in the calculations on a certain level of method and model -often on DOI: 10.1002/adts.202000223 the same level as the observed molecule or only on a slightly reduced one. Examples for such an environmental influence can be electronic or geometrical effects due to non-negligible intermolecular interactions between solute and media as for example hydrogen bonds. [2] However, the accurate treatment of these influences is not always trivial. For instance a large size of the observed molecule (mostly biomolecules) may cause the environment to be not only around the molecule but contained within it. Already in isolated systems, complicated cases occur where large amplitude motion [3] plays a role which becomes even more complicated when an active environment has to be considered.
Coming back to the first point of the accurate treatment, there are several ways to deal with such an intrusive environment. In this work we will consider two approaches, the first being based on static quantum chemical calculations of the observed molecules including explicitly some of the media molecules or fragments. We call the whole assembly a cluster. [4,5] The individual properties of several clusters can be calculated and averaged to arrive at the desired observable property. Besides the method of quantum cluster equilibrium (QCE) [4][5][6] or cluster weighting (CW) [7] there are many variants of cluster models. The first work, which utilized clusters obtained from snapshots of trajectories, was carried out by the Hermansson et al. [8] and little later by the Huber group. [9] Pliego and Riveros [10] call their approach a cluster-continuum model, while Xu and coworkers [11] have developed the cluster-in-a-solvent method. A cluster can contain explicit (micro-solvation) or implicit (continuum) solvent molecules [12] and the best results are obtained by the combination of both, [13] for an excellent review see ref. [11]. A crucial step in the application of all these different models is the construction of the clusters since the description of the often very involved bulk phase has to be reduced to a small set of representative structures. For this reason, several construction schemes for clusters have been proposed. The clusters can be generated based on approaches of neural network or machine learning, they can be intuitively guessed or extracted from molecular dynamics simulations. [8,9] Alternatively they can be obtained by conformer rotamer ensemble sampling as recently presented by Grimme. [14] When subsequently the Boltzmann weighting (BW) is applied, clusters of equal size and composition are usually selected. However, in the case of QCE, clusters of different size and composition [15] can serve as input for the mational changes [56,57] and environmental effects induced by solute-solvent or solute-solute interactions which thus have to be included in the calculations to match the experiment. [58,59] Since the calculation of VCD spectra based on AIMD simulations relies on the Fourier transformed cross-correlation function of the electric dipole moment and the magnetic dipole moment, [37,60,61] these properties have to be calculated. Traditionally, the electric dipole moments are accessible via Wannier centers [62,63] and more recently by integrating the electron density over a molecular cell, which is determined by a radical Voronoi tessellation. [64] The magnetic dipole moments can be obtained by classical approaches [26,38] or by quantum mechanical approaches. [53,54] To quantify the compliance, the overlap between the experimental and calculated spectra can be calculated, for example, with the SimVCD measure. [65] A value of +1 represents a perfect agreement between experiment and calculation, and a value of −1 indicates that both spectra are mirrored to each other. Thus, −1 implies, in principle, opposite enantiomers. Starting from 0.2, the value allows a safe (>90%) determination of the absolute configuration.
The article is structured as follows: The cluster weighting methodology, cluster generation, cluster weighting input details, and the quantum chemical methodology used are described in detail in Section 2. We then give details about the AIMD simulations and the trajectory analysis (Section 3). In Section 4, we present the results. The following Section 4.1 gives a detailed analysis of the AIMD spectra for the gas phase and the bulk. In Section 4.2, we discuss functional and basis set including CW input parameter dependencies on the spectra obtained from CW in order to compare in Section 4.3 AIMD with CW spectrum. An analysis of the moments in Section 4.4 follows. Finally, Section 5 gives a general discussion and conclusions.

The Cluster-Weighting Method
The QCE approach [4,6] begins its description of gaseous and liquid phases by a so-called cluster gas, which is in thermodynamic equilibrium and consists of non-point particles. [4] The chemical equation for the cluster equilibrium of a pure substance is with C 1 and C ℘ being the monomer and a cluster of size i(℘), respectively. In 2011 we extended the QCE to binary systems (bQCE) [15] : A cluster C ℘ is composed of i(℘) and j(℘) monomers of the components 1 and 2. The primary goal of the model is to obtain the equilibrium cluster populations N ℘ that minimize the free energy A of a system with volume V, temperature T and total number of monomers N tot = N tot 1 + N tot 2 , where N tot 1 and N tot 2 is the number of monomers of the components 1 and 2. If the cluster www.advancedsciencenews.com www.advtheorysimul.com populations are known, the systems' total partition function Q tot can be calculated as follows where {N ℘ } denotes the full set of total cluster populations N ℘ and q tot ℘ is a single cluster partition function given by The rotational, vibrational and translational partition functions q rot ℘ (T), q vib ℘ (T) and q trans ℘ (V, T) are calculated from standard equations for the rigid rotator, harmonic oscillator and particle in a box, respectively. [66] Based on the non-point-shaped particles, the volume is corrected by the translational partition functions, since the total volume is not accessible for the particles. Therefore, the exclusion volume V ex , given by is subtracted from the phase volume V. Thereby v ℘ denotes the volume of cluster C ℘ and b xv is an empirical parameter which corrects the cluster volumes. The electronic cluster partition function q elec ℘ (V, T) is calculated from the referenced electronic energy elec ℘ and is supplemented by a term that takes into account interactions between clusters described by a mean-field-type attractive energy. The electronic partition function reads as follows where k B is the Boltzmann constant and a mf is the mentioned mean-field parameter, thus scaling the interactions between clusters. Thus, a mf is the second empirical parameter in the theory.
The systems' total partition function Q tot is a function of the cluster populations {N ℘ } obtained by minimizing the free energy This function should be minimized under the conditions of a fixed total number of particles and a fixed total mass. Two polynomial equations can be derived from both constraints, namely the population polynomial and the volume polynomial Both polynomials result in a non-linear equation system, which can be solved numerically. The two empirical parameters a mf and b xv are optimized on the experimental reference (usually one density and the boiling point). The obtained partition functions and populations can be used in further calculations. Applications of the QCE method involve many neat and binary systems, including water, [4,67,68] aqueous solutions [15,69] and even ionic liquids. [70] In addition, we calculated activity coefficients [71] and ionicity [17,18,72] with binary QCE and we applied QCE or cluster weighting CW to vibrational circular dichroism [7] .

Cluster Generation
Since the description of the complex bulk phase is based on a limited set of clusters, the choice of the input structures for the cluster-weighting is an essential step [7] . First, it was shown that for the calculation of VCD spectra of conformationally flexible molecules more than one conformer has to be considered. [55] Furthermore, the CW calculation allows the simultaneous inclusion of different oligomers and for these, constitutional isomers can be selected. As constitutional isomers we consider such clusters with different motifs, that is, for a trimer a ring and a chain can be constructed, for larger clusters even more structure motifs appear. All clusters considered are taken from our previous work. [7] We created the global minimum structure for the monomer up to the decamer by applying a structure optimization based on genetic algorithms on the classical force field level of theory. For this purpose we combined the software Ogolem, [73,74] with the Amber 2016 molecular dynamics package [75] and the generalized Amber force field (GAFF). [76]

Quantum Chemical Calculations
Since the goal of the work is to compare CW to AIMD, we calculated the clusters quantum chemically on different theoretical levels. All clusters were optimized using the BP86, [77,78] the BLYP, [77,78] and the B3LYP [77,78] functional. For the generalized gradient approximation (GGA) functionals we used both the def2-SVP and the def2-TZVP basis set, [79] for the B3LYP functional only the def2-TZVP basis set. For each calculation dispersion corrections according to the D3 dispersion scheme [80,81] were used. The neglect of these corrections was tested and some examples show improved results, but these results were found for the conformers of rotaxanes. [82] Nevertheless, we expect in principle that the dispersion correction will play an important role already for a dimer. All clusters were treated with implicit solvation, which was considered by the solvation model Cosmo [83] using the dielectric constant of (R)-butan-2-ol at 293.15 K corresponding to 17.26. [84] In general, we have chosen 10 −4 a.u. as convergence criteria for the maximum norm of the cartesian gradient and 10 −8 a.u. for the energy of the SCF. We use a fine multigrid (m4) [85] for the numerical integration of the functional. Some tests for the sensitivity of the electronic energy with respect to the mentioned parameters are given in Table 1. Note that in some problematic cases we increased the multi-grid to m5 for some larger clusters, all information is given in the supporting information. For all calculations, we list the total energy E h in Hartree and (ΔE) interaction energy in kJ mol −1 with one digit after the decimal. All calculations were performed with the BP86 functional and the SVP as well as the TZVP basis set.
The results of the BLYP functional and some results of the BP86 functional were compared with unscaled and scaled AIMD calculations. For details see methodology of AIMD calculations.
All quantum chemical optimizations as well as the calculations of the individual cluster VCD spectra were performed using the Turbomole 7.4 software package [44,[86][87][88] and the broadening of the VCD spectra was modeled by a Lorentzian band shape employing a full width at half maximum of 8 cm -1 .
Since both the Boltzmann and the cluster weighting are very sensitive to the energy gained, we also tested various parameters of the quantum chemical calculations for the dimer. The corresponding spectra are shown in the supporting information (Figures S4 and S5, Supporting Information) and do not show a pronounced sensitivity. Table 1 shows the absolute interaction energies for the dimer and monomer and the interaction energies for the dimer. We observe a good agreement between the different energies, so we are sure that our parameters are adequately chosen. The interaction energies show deviations of less than 1 kJ mol −1 when comparing the different criteria.

Cluster Weighting Calculations
The QCE equilibrium populations were evaluated with Peacemaker, [4,5,89,90] a publicly available software package for neat and binary QCE calculations. [90] All calculations were performed for a fixed pressure of 101.325 kPa and in a temperature range of 200-500 K. The cluster volumes were defined as van der Waals volumes with radii from Bondi's compilation. [91] As experimental input the boiling point T b = 372.0 K [84] and the density at 293.15 K 293.15 = 0.806 g cm −3 [84] of (R)-butan-2-ol were employed. The empirical parameters a mf and b xv were optimized between values of 0.0 and 3.0 to reproduce the boiling point and density at 293.15 K as accurately as possible.

AIMD Simulations
Several systems were simulated from AIMD, nine single molecules (based on all possible monomer conformations, referenced as gas) and 16 molecules (coined bulk) of (R)-butan-2-ol. While the single molecule trajectories were created in this work, the 16 molecule trajectory we took from ref. [26]. The starting conformations can be understood according to the dihedral angles defined in Figure 1 with gas1 ≡ ( 1 = 300; 2 = 300), gas2 ≡ ( 1 = 180; 2 = 180), gas3 ≡ ( 1 = 300; 2 = 180), gas4 ≡ ( 1 = 180; 2 = 60), gas5 ≡ ( 1 = 60; 2 = 60), gas6 ≡ ( 1 = 300; 2 = 60), gas7 ≡ ( 1 = 180; 2 = 300), gas8 ≡ ( 1 = 60; 2 = 300), gas9 ≡ ( 1 = 60; 2 = 180) and angles in degrees. Of course these rearrange partly in the equilibration runs and partly stay close to their local minima. All AIMD simulations are based on density functional theory using the QUICKSTEP module [92] within the CP2K [93] program package. The hybrid Gaussian and plane waves (GPW) approach was used to calculate the energies and forces on the atoms. As basis set we choose the molecularly optimized double-basis set (MOLOPT-DZVP-SR-GTH). [94] Note that this generally contracted basis set contains diffuse primitives and is obtained by minimizing a linear combination of the total energy and the condition number of the overlap matrix. It is used in combinations with pseudopotentials and the hybrid approach and is therefore not comparable with standard quantum chemical basis sets. Therefore, for a given accuracy in the total energy, considerably fewer basis functions are required than in the usual split valence scheme. This can make the MOLOPT-DZVP-SR-GTH comparable to the def2-TZVP basis set, which is important for the comparison between CW and AIMD spectra. For the functional we choose here and in the previous work [26] the generalized gradient approximations (GGA) functional BLYP [77,78] and the corresponding BLYP Goedecker-Teter-Hutter pseudopotentials for core electrons, [95][96][97] since hybrid functionals are outside the scope of standard AIMD simulations and BLYP showed good performance for condensed phase systems in several comparative studies. [98,99] The choice of this functional in the AIMD simulations is the reason why we have also repeated the static quantum chemical calculations of all clusters with the same BLYP functional. A 280 Ry density CUTOFF criterion with the finest grid level was used together with multi-grids number 5 (NGRID 5 and REL_CUTOFF 40) using the smoothing for the electron density (NN10_SMOOTH) and its derivative (NN10). [92] Dispersion interactions were accounted for by applying the DFT-D3 type of a pair potential van der Waals density functional. [80,81] For the SCF calculation, the default value (10 −5 ) was used as target accuracy for the SCF convergence. The DIIS minimizer [93] was employed to achieve a faster orbital transformation by direct inversion in the iterative subspace. The maximum number of SCF iterations to be performed for one iteration was set to 200.
The following equilibration and production runs were carried out applying periodic boundary conditions to avoid boundary effects. The canonical (NVT) ensemble was applied using the Nosé-Hoover-chain thermostats [100,101] with a time constant of the thermostat chain of 50 fs (first equilibrium run), 10 fs (second equilibrium run), and 50 fs (production run). The first equilibration runs were performed over 5 ps (10 000 steps with a time step of 0.5 fs) and the temperature was increased to 600 K to achieve faster equilibration. The second equilibration runs were performed over 5 ps (10 000 steps with a time step of 0.5 fs) defining the keyword REGION MASSIVE (thermostating each individual atom) and lowering the temperature to 400 K to relax the system to the correct temperature. The subsequent production runs were performed over 30 ps (60 000 steps à 0.5 fs) for the bulk and gas2 to gas9 systems and over 50 ps (100 000 steps) for the gas1 system at a temperature of 400 K. Electric and magnetic dipole moments were analyzed for the gas1 trajectory. For this trajectory also different populations analysis were carried out to obtain partial charges, namely Blöchl, [102] Mulliken, [103] and Löwdin [104] as implemented in the CP2K package. [93] Please note here that the temperature has been set to 400 K to improve the sampling. In addition, several AIMD studies [99,105] have shown that the temperature is 20-50 K too low compared to experiment. The box length for the bulk system was set to 13.45 Å, which corresponds to a density of 0.81 g cm −3 at 293.15 K. For the gas phase simulation, the box length was set to 12.00 Å.

Analysis
Our open source program TRAVIS (Trajectory Analyzer and Visualizer) was used to investigate the time-dependent structural evolution of the simulated systems and to calculate IR and VCD spectra from the AIMD trajectories. [106,107] The results are visualized with the Xmgrace software [108] or, alternatively with the program gnuplot (version 5.0) [109] to generate 3D data plots.

Calculation of Vibrational Spectra
The AIMD calculation of different spectroscopic properties is based on the Green and Kubo, [21,22] approach. To transfer the correlation functions into a frequency spectrum, a Fourier transformation according to the Wiener-Khintchine theorem is applied: [110,111] Using the time derivative, which is advantageous for numerical reasons, the IR intensities A( ) are then given by For the comparison with experimental data the corresponding wavenumber-dependent representation is Abbate et al. [60] developed the time-correlation ansatz for the calculation of VCD spectra from molecular dynamics trajectories, which is based on a cross-correlation of the electric and magnetic dipole moment: The corresponding wavenumber-dependent expression is Please note that the VCD intensity has a slightly different prefactor with respect to the speed of light c compared to the IR intensity due to the unit of the magnetic moment (which is A m 2 ). Nevertheless, also the VCD intensity has the unit m 2 mol −1 . The magnetic dipole moment is calculated by the electric current density, for more details, see ref. [26] The radical Voronoi tessellation is used to calculate the electric dipole moments. [64] All Voronoi cells C i of all atoms belonging to the same molecule are united to obtain molecular cells M k . The recognition of the molecules is available in TRAVIS. [106,107] Molecular dipole moments are obtained by integrating the electron density (r), which is accessible for each time step in an AIMD simulation, as a function of the spatial coordinate r over the previously defined molecular cells: Adv. Theory Simul. 2021, 4, 2000223 Table 2. Overlap of from SimVCD tool [65] between experimental gas and bulk phase VCD spectra of (R)-butan-2-ol with calculated VCD spectra from AIMD simulations. The calculated frequencies were scaled with factors between 0.95 and 1.05, which are given in the first row. Bold marked overlaps show the best agreement with the experimental references. Please note, that gas is the mean of the nine trajectories gas1 to gas9. Conformational plots for the simulations gas1 to gas9 can be found in Figure 1.

Comparison of Ab Initio Molecular Dynamics Simulations and Experiment
In general, ab initio (and traditional) molecular dynamics simulations take anharmonicity into account, since molecular dynamics does not use a harmonic potential approximation as in frequency analysis in standard static quantum chemical calculations. However, the motion of the nuclei in standard simulations is treated classically. Therefore, nuclear quantum effects are not considered. The electronic structure is based on the Born-Oppenheimer approximation and the electronic structure is based on generalized gradient approximation (GGA) DFT, which has been shown to introduce another error to AIMD (see e.g., the work of Klein [99] and the discussion of the self-interaction error therein). The influence of the electronic structure on the spectrum (BLYP vs BLYP-D3) was also demonstrated. [98] In addition, the step-wise integration of the equation of motion has an influence over the time step size, which was demonstrated, for example, by Bourȃnd co-workers. [61] In quantum chemistry, besides the treatment of anharmonicity, such problems are sometimes solved by a simple approach using frequency scaling. [43,112] To maximize the overlap with the experiment, we have also determined scaling factors for the calculated frequencies in the gas and in the bulk phase. The results are listed in Table 2.
Our scaling factors for the gas (1.02) and for the bulk (1.03) lead to overlaps of gas = 0.536 and bulk = 0.456, indicating a quite satisfactory agreement with the experiment. The gas phase spectrum was obtained from nine simulations with different initial conditions using all possible conformers (see Figure 1 and computational details) as the arithmetic mean between the nine trajectories (see Figure 2). Both the C-O bond and the central C-C bond exhibit three distinct torsional minima on the potential energy surface, resulting in nine conformations of the molecule, all of which seem to be sufficiently visited in the AIMD simulations, as can be seen from the conformation plots previously created in ref. [113], see Figure 2.
The excellent agreement between the averaged gas phase simulations and experiment suggests to follow this approach, where the simulations are started from all possible monomer conformations in order to correctly sample the conformational space. It is also clear that only averaging leads to the complete exploration of the conformational space, compare Figure 1 with 2. This is in line with the philosophy of the more recent research topic . VCD spectra for gas (left) and bulk (right) in the range between 900 and 1400 cm −1 . The experimental gas spectrum is a diluted solution of (R)-butan-2-ol in carbon disulfide, both experimental gas and bulk are from ref. [56] Top: Comparison between experiment and unscaled AIMD data. Middle: Experiment and scaled AIMD data. Bottom: Scaled and unscaled AIMD data. Note that the gas indicates the mean of the nine trajectories from gas1 and gas9. Blue: unscaled; red: scaled; black: experiment. For corresponding overlaps see Table 2.
of enhanced sampling. [114][115][116][117][118][119] For the averaged simulations, it is obvious that the conformation with 1 = 180 • and 2 = 300 and 180 • are the most frequently visited conformations. This is fully consistent with our previously published results from applying cluster weighting to the monomer conformations in static quantum chemical calculations. [7] Interestingly, due to intermolecular interactions, the dominance of the conformers shifts in the bulk to those with 1 = 60 and 180 • and 2 = 60 • .
In Figure 3, we show the gas and bulk AIMD and experimental VCD spectra for the scaled and unscaled results. As discussed earlier, [26] due to the theoretical standard, a noticeable red-shift of the calculated data occurs (Figure 3 top and bottom panels). However, we observe an excellent agreement between experiment and calculations, when applying our scaling factors. In particular, AIMD agrees exceptionally well with the −/+ VCD pattern in the gas phase at 1145 and 1242 cm −1 . For the bulk AIMD, the spectrum agrees well with a positive peak at around 1109 cm −1 followed by a strong negative peak at 1152 cm −1 and another weak positive peak at 1290 cm −1 of the experiment. However, below 1100 cm −1 the agreement between experiment and simulation is worse. The calculated overlap between bulk simulation and experiment in this range is only 0.112. One of the influences for the lack of the correct description of the bandwidth could be the length of the simulations. In Figure S2, Supporting Informa- tion, we confirm the change of the bands with simulation time. We also show that our trajectory length appears sufficiently converged for bands above 1100 cm −1 . Furthermore, the intensity ratios are also very well reproduced in the bulk. This is also reflected in the high overlaps of 0.536 (gas) and 0.456 (bulk) ( Table 2). That the overlap numbers should be treated with caution and can only be interpreted in combination with the spectrum is shown by the small overlaps between the scaled and the unscaled spectrum (see Table S1, Supporting Information, bulk value), which are the same spectra shifted only by the scaling factor. Overall, given the short simulation time (30 ps) and system size (16 molecules) the agreement between AIMD simulation and experiment is excellent. Such a locality behavior of the electric dipole moment has been investigated previously. [120] It seems that all important conformations are sampled (see Figures 1 and 2), especially in the averaged gas phase.
Finally, we show the gas and bulk IR spectra in Figure 4. Again, we observe that scaling leads to a better agreement between bulk and experiment. Since AIMD covers anharmonicity effects the discrepancy between experiment and simulations can be attributed to the incorrect description of GGA functionals, the small basis set or the choice of time step.
The full scaling range from 0.95 to 1.05 for the gas and the bulk AIMD is given in Table S2, Supporting Information. Furthermore, the comparison between VCD and IR is possible in Table 3. Overlap of from SimVCD tool [65] between experimental bulk phase VCD spectrum and calculated VCD spectra employing CW for different sized cluster sets on different theoretical levels.  Figure S1, Supporting Information, where we have plotted the spectra together.

Comparison of Different Functionals in Cluster Weighting
We now turn to the results of the cluster weighting, that is, static quantum chemical calculations based on the QCE approach, and in this section we compare the performance of different functionals. Previously, [7] we build a cluster set, which we have coined S4 gm . The superscript "gm" stands for global minimum and the four stands only for the fourth. It is important, that this cluster set contains all global minimum oligomers from one to ten. Furthermore it is important to understand, that all clusters included in the calculation are given a weight by the QCE calculation. The individual spectrum of each cluster is then weighted to obtain a total spectrum. For these clusters, we tested the agreement with the experiment step by step, that is, by increasing the oligomers contained.
In Table 3, the cluster label in the cluster column indicates the largest cluster included for the overlap calculations, all smaller clusters are always included.
The clusters were calculated using the D3 dispersion correction and implicit solvation with Cosmo, see Section 2.3. It is disappointing that a strict systematic trend is missing, although increasing the cluster size seems to lead to slightly better results. Compared to the unscaled bulk AIMD, the overlaps are much better, with the BLYP/SVP combination performing best. The early convergence of the SimVCD value for some level of theory suggests that inclusion of clusters up to the decamer is unnecessary. Indeed, visual inspection of the AIMD bulk trajectories shows that mainly bent chain trimers and tetramers are formed, a pentamer and a hexamer occur only rarely. Table 3 contains two more columns, in which we artificially overbind the interaction energy of the clusters by a factor of two. The reason is simply that the worse basis set SVP exceeds the TZVP basis set, see column two and three as well as four and five in Table 3, and one difference is Figure 5. Cluster weighting bulk VCD spectra for different level of theory in the range between 900 and 1400 cm −1 . The experimental spectrum is from ref. [56]. Blue: unscaled; red: scaled; black: experiment. For corresponding overlaps see Table 4.
the overbinding of SVP compared to TZVP. As a result both functionals BLYP and BP86 perform better than the original TZVP data, but only in the case of BP86 an extreme improvement (even compared to SVP) can be observed. These trends are also not systematic, which is easy to understand. In QCE calculations a difficult equilibrium of the input data leads to significantly different populations (see Supporting Information), making error tracking and the understanding error compensation more difficult than in a case where only electronic energies of one molecule are compared.
The corresponding spectra of the S4 gm set for different level of theory are exhibited in Figure 5 blue curves.
Similar to the AIMD simulations we also scale the spectra, which is even more serious, because here we counteract the missing description of the anharmonicity besides incompleteness of the basis set, level of theory, and remaining errors. The results are presented in Table 4 and in Figure 5, red curves.
All spectra based on GGA functionals have to be shifted to higher wavenumbers, while the B3LYP data have to be scaled by a factor lower than one. The best performance can be found for the B3LYP/TZVP closely followed by the BP86/TZVP results. This is in accordance with the findings found in the literature for gas phase calculations. [121] Besides this the artificially doubled energy BP86/TZVP ⋆ slightly exceeds the good B3LYP/TZVP performance, but this is not the case for the BLYP/TZVP ⋆ , which shows that there is an important influence besides the energy. The overall very good performance of BP86/TZVP ⋆ and BP86/TZVP is interesting and could be due to the good frequency calculation Table 4. Overlap of from SimVCD tool [65] between bulk experimental and calculated VCD spectra employing CW for different levels of theory. For all CW spectra the complete S4 gm set has been employed. ⋆ indicates that the referenced binding energies have been doubled in the CW calculation. Bold marked overlaps show the best agreement with the experiment. Figure 6. Cluster weighting bulk IR spectra for different levels of theory in the range between 900 and 1400 cm −1 . The experimental spectrum is taken from ref. [56]. Blue: unscaled; red: scaled; black: experiment.
rather than in overbinding, that is not observed for these levels of theory compared to the other methods, see Table S5, Supporting Information. The remaining data from different level of theory are also above an overlap with experiment of 0.2. Therefore all cluster weighting calculations are sufficient and serve well to determine the absolute configuration in bulk phase systems. As with the AIMD calculations the performance of the CW below 1100 cm −1 is not as good as when considering the whole range. To give an estimate of Table S4, Supporting Information, we list the overlaps below 1100 cm −1 . Interestingly, BP86/TZVP outperforms the AIMD simulation and B3LYP/TZVP is in the same range. For completeness, Figure 6 displays the corresponding IR spectra. Finally, we show in Figure 7 the received CW populations based on different functionals and basis set.
The variety of curves is deflating, but also not surprising considering that the variety already exists only for the interaction energies (Table S5, Supporting Information). On closer examination, however, some systematic trends can be derived from this. All methods occupy only the monomer in the gas phase with the exception of BLYP/SVP, which also occupy a majority of monomer population, but still have a certain amount of dimer, tetramer, and pentamer in the gas phase and BP86/TZVP ⋆ which occupies the dimer and the trimer. In the liquid phase the clusters above the hexamer (at 300 K) seem to be insignificantagain with the BLYP/SVP and BP86/TZVP ⋆ exception. This strongly suggests that higher clusters than the hexamer are dispensable for calculations of hydrogen bonded systems. Unfortunately, the most well-performing functional/basis set combinations, where certain clusters play a role, do not agree. While B3LYP/TZVP populates the pentamer followed by the hexamer and slightly followed by the trimer and tetramer, the BP86/TZVP combination favors the monomer, trimer, and dimer. These results are certainly due to the slightly different structures obtained by geometry optimizations (Table S3, Supporting Information) together with the extreme differences in interaction energies, but they also suggest what was observed before. [68,122] Namely, that structural motifs contained in clusters are more important than a certain cluster size, that is, a particular www.advancedsciencenews.com www.advtheorysimul.com oligomer. The artificially increased energy combination BP86/TZVP ⋆ apparently contributes to the curiosity of this graph in that here the decamer is highly and the heptamer and trimer are slightly populated. But this is also understandable if one considers the interaction energies (Table S5, Supporting Information). It reflects that higher energies contribute to populate certain clusters, see BLYP/SVP and BP86/TZVP ⋆ . The outlier BP86/SVP behaves differently because of very low frequencies.
Frequencies are also an input quantity in the vibrational partition function for QCE calculations, but Grimme showed that these low frequencies must enter the rotational partition function instead of the vibrational because they are prevented from rotating. [123] Therefore, it is not justified to believe that a certain cluster or the choice of a large size of clusters makes the spectrum good. It remains to be seen how much other isomeric clusters, that is, other local minima for a particular oligomer, or clusters with more conformations change the picture. However, since CW is strongly dominated by energies an improvement by isomeric clusters with higher-energies (at least much higher energies) seems unlikely.
In summary, we can give the following advice when performing static quantum chemical calculations with cluster weighting: Given that it is affordable in terms of computer time, it is best to choose the hybrid functional B3LYP/TZVP combination, followed by the BP86/TZVP combination for bulk calculations. If computer time is a problem, BLYP/SVP calculation may be an alternative. This relatively good performance of BLYP/SVP for the bulk strongly suggests that large error compensations are involved.

Comparison of Ab Initio Molecular Dynamics Simulations with Cluster Weighting
We now compare the CW and the AIMD VCD spectra to see how much the former can model the latter, see Figure 8. To reduce the influence of different parameters, we mainly compare the AIMD bulk simulation (scaled and unscaled) with the BLYP/SVP and BLYP/TZVP, but we also add the B3LYP/TZVP data, since it achieved the best overlap values with the experiment.
All CW data compare well with the AIMD in the strong negative peak at around 1152 cm −1 . While the BLYP with smaller SVP basis set shows at least a hint of the positive peak at 1109 cm −1 , this hint is missing in the BLYP/TZVP combination, which already in the previous subsection 4.2 led to a bad comparison to the experiment. An exception to the poor performance of BLYP/TZVP with respect to AIMD is the weak peak at 1290 cm −1 , which is almost reached by it. All three static calculations strongly disagree with AIMD in the lower wavenumber region, showing an opposite (positive) peak to the negative one at around 975 cm −1 . While B3LYP/TZVP-CW also lacks these negative peaks, it shows on the other hand a weak positive peak at around 965 cm −1 which is missing in the AIMD. It remains to be seen whether AIMD, which is based on other functionals changes the spectrum extensively and for the better, but for the time being the CW approach is better based on the BP86 or B3LYP functional, as has been shown several times in literature for Boltzmann weighting in the gas phase. [55] Figure 8. Cluster weighting (cyan) and AIMD (red) bulk VCD spectra for different levels of theory in the range between 900 and 1400 cm −1 . [56]

Magnetic Moments from Ab Initio Molecular Dynamics Simulations
The calculation of the VCD spectra from the correlation approach requires the knowledge of the moments along the trajectories. In Figure 9, we show a short snapshot of the moments' time development (upper panels) together with their distribution (middle panels) and the x-, y-, and z-component distributions (lower panels) for the gas1 trajectory.
Since the Thomas-Kirchner ansatz [26] to obtain the magnetic moments is quite computationally intensive, we also test here Cho's long established approach for classical or QM/MM trajectories. Instead of using QM/MM or classical simulations, we apply Cho's ansatz to AIMD trajectories, where the electric dipole moment ⃗(t) is obtained via: with q i (t) being the partial charges at time t and ⃗ r i (t) the position vector of the molecular atoms. Similarly, the magnetic dipole moment ⃗ M(t) is given by: with c being the speed of light and ⃗ v i being the according velocity.
In Figure 9, we show on the left side the corresponding electric dipole moments based on various population analysis marked in different colors. In the upper picture we see that all population methods follow the same trends for , with the largest dispersion  [26] Red, blue and orange are the classical approach developed by Cho and co-workers [37,38] on the basis of Blöchl, [102] Mulliken, [103] and Löwdin [104] partial charges. from the main trends coming from the Löwdin population analysis. This agreement is not surprising since we choose the same trajectory and electronic structure method on which the analysis is based. For the Voronoi approach we observe a dipole moment in the range of 1 and 2.5 D with a mean of 1.71 D, Blöchl lies between 1 and 2.25 D with a mean of 1.53 D, Mulliken occupies the range between 0.5 and 2 D with the average of 1.25 D, and Löwdin provides the lowest dipole moment distribution between zero and 1.5 D and an average of 0.81 D. A comparison of these values with the experiment of 1.66 D [124] clearly shows the superiority of the Voronoi together with the Blöchl approach.
The situation is different with the magnetic moment. The agreement between Voronoi and the three other approaches is no longer in phase, because the calculation of the moment in the case of Voronoi is based on a completely different method. Nevertheless the range covered seems to be comparable and the moments based on the population analysis obviously fluctuate. Therefore, we now focus on the calculated spectra in Figure 10.
Again, the IR spectrum based on electric dipole moments shows a large agreement between the different approaches for the peak locations. Only the intensities differ more between both the Mulliken and the Löwdin compared to the Voronoi approach. Blöchl, [102] and Voronoi compare very well in both location and in- Figure 10. IR and VCD AIMD gas spectra. Black from the Voronoi approach. [26] Red, blue and orange are the classical approach established by Cho and co-workers [37,38] based on Blöchl, [102] Mulliken, [103] and Löwdin [104] partial charges. tensity. The difference in the intensity for different electric dipole moment methods is nothing new. It was observed earlier, when we compared our Voronoi and the Wannier function approach. [64] This intensity problem is naturally more pronounced in the VCD case. While our Thomas-Kirchner [26] approach has the experimental +/− feature between 1020 and 1150 cm −1 , Löwdin fails completely, Mulliken underestimates the extremity in the second peak and Blöchl performs well. The second positive peak around 1200 cm −1 is only slightly reproduced by Mulliken, the other two approaches of population analysis neglect to show this peak. In summary, we can say that these tests are very promising, since the computer time and data handling are much less complex. Furthermore, the outcome especially the case of the Blöchl charges [102] is very satisfactory.

Discussion and Conclusion
Using the example of (R)-butan-2-ol in the gas and bulk phase, we investigate the problems that arise when predicting the vibrational spectrum of a flexible molecule in a non-idle environment. Two approaches are investigated in detail. The methods used are ab initio molecular dynamics simulations and cluster weighting based on static quantum chemical calculations.
From AIMD simulations we calculated the gas and the bulk phase spectra, with excellent agreement when we applied a scaling factor of 1.03 (bulk) and 1.02 (gas) for the selected www.advancedsciencenews.com www.advtheorysimul.com BLYP/MOLOPT-DZVP-SR-GTH basis set combination. (R)butan-2-ol shows several conformations connected by barriers low enough to be influential in experiment but too high to be visited by only one initial conformation in the ab initio simulation. To scan all the important conformational regions, we simulated from the different initial conditions of all monomer conformers. A subsequent conformation analysis shows, that taking the average of all trajectories covers the whole conformation range apparently with the right proportion resulting in an excellent spectrum. Therefore, this ansatz, where multiple trajectories are simulated starting from all monomer conformers and then averaged over them, is a valid approach for an accurate spectrum when considering a flexible species in gas phase or a single solute molecule is solvated in a different solvent than the solute molecule. It should be noted that AIMD is computationally very demanding.
The second approach used in this work, is much less computationally demanding. Therefore, we could compare different functionals (BP86, BLYP, B3LYP) and basis sets (SVP and TZVP) and confirm for the bulk, what is valid and already known for the gas phase infrared and vibrational circular dichroism spectra calculation. The best results are obtained for the BP86/TZVP and B3LYP/TZVP functional and basis set combination. This superiority of these two level of theory may be due more to the good frequencies (input for cluster weighting) than the energies and geometries. We found a rough converging trend with increasing cluster size indicating that for hydrogen bonded systems clusters up to tetramers or hexamers might be sufficient. Again the scaling leads to very good agreement with the experiment. In total, all methods resulted in an agreement with the experiment above 0.2, which allows a safe assignment of the absolute configuration.
In the last part we analyzed the moments of one gas phase trajectory based on our previously developed method and on different population analyses. The Voronoi electric dipole moment is well matched to the experiment followed by that of the Blöchl charges. There is a perfect agreement between Voronoi and population analysis even in the time development and the resulting IR spectra agree very well in peak location. The magnetic moments from the population analysis and from our Voronoi approach do not fluctuate in phase, but still show a similar distributions. Interestingly, the Blöchl vibrational circular dichroism spectrum reproduces the Voronoi characteristics very well with much less computer time and makes it a valuable approach for the calculation of such spectra.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.