Modeling the Vibrational Circular Dichroism
Spectroscopy of Phenylcyclohexanediol Solvated in
Dimethyl Sulfoxide Using Polarizable Molecular
Dynamics

Emeline Boyer, Lou D'haese, Jessica Bowles, Sascha Jähnigen, Katia Le Barbu-Debus, Rodolphe Vuilleumier, Anne Zehnacker, Florent Calvo, Carine Clavaguéra

► To cite this version:

Emeline Boyer, Lou D'haese, Jessica Bowles, Sascha Jähnigen, Katia Le Barbu-Debus, et al. Modeling the Vibrational Circular Dichroism Spectroscopy of Phenylcyclohexanediol Solvated in Dimethyl Sulfoxide Using Polarizable Molecular Dynamics. Journal of Physical Chemistry B, 2026, 10.1021/acs.jpcb.5c05306. [hal-05451847](https://HAL Id: hal-05451847)

HAL Id: hal-05451847

https://universite-paris-saclay.hal.science/hal-05451847v1

Submitted on 10 Jan 2026

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Modeling the Vibrational Circular Dichroism
Spectroscopy of Phenylcyclohexanediol
Solvated in Dimethyl Sulfoxide Using
Polarizable Molecular Dynamics

Emeline Boyer,† Lou C. G. D'haese,†,⊥ Jessica Bowles,† Sascha Jähnigen,‡ Katia Le Barbu-Debus,¶ Rodolphe Vuilleumier,§ Anne Zehnacker,¶ Florent Calvo,*,|| and Carine Clavaguéra*,†

†Institut de Chimie Physique, Université Paris Saclay, CNRS UMR8000, 91 405 Orsay,
France

‡Freie Universität Berlin, Fachbereich Biologie, Chemie, Pharmazie, 14195 Berlin,
Germany

¶Institut des Sciences Moléculaires d'Orsay, Université Paris Saclay, CNRS UMR8214, 91
405 Orsay, France

§CPCV Laboratory, Département de Chimie, Ecole Normale Supérieure, PSL University,
Sorbonne Université, CNRS, 75005 Paris, France

||LIPhy, Université Grenoble Alpes, CNRS, 38 000 Grenoble, France

⊥Namur Institute of Structured Matter (NISM), Université de Namur, Laboratoire de Chimie Théorique (LCT), Namur, 5000, Belgique

E-mail: florent.calvo@univ-grenoble-alpes.fr; carine.clavaguera@universite-paris-saclay.fr
Abstract

Modeling vibrational circular dichroism (VCD) spectra in solution remains a com-
plex task due to the intricate interactions involved. Recent advances in computational
chemistry, particularly the implementation of polarizable force fields, have greatly im-
proved the description of environmental effects in classical molecular dynamics simula-
tions. In this study, we use the AMOEBA polarizable model to calculate VCD spectra
through the expression of the electric and magnetic dipole moments developed recently
[Bowles J. et al. ChemPhysChem 2024, 25, e202300982] in combination with the use
of induced dipole models. This methodology enables the reliable treatment of flexible
molecules and has been applied to the 1-Phenyl-1,2-Cyclohexanediol (PC) molecule in
dimethyl sulfoxide (DMSO). The computed infrared and VCD spectra provide detailed
insights into molecular conformations and solute-solvent interactions. An in-depth in-
vestigation was conducted into the possible effects of system size, concentration, and
spectral convergence. Comparison with experimental vibrational frequencies in deuter-
ated solvent further supports the interpretation of spectral features, although residual
differences in the spectral shifts and broadenings prove to be challenging for theory,
especially in the case of VCD.

**Introduction**

Vibrational circular dichroism (VCD) is a chiroptical technique that allows determining ab-
solute configurations of chiral molecules by discerning absorption differences in right- and
left- circularly polarized light within the infrared (IR) range. This subtle technique unveils
molecular conformational variations and is sensitive to both intra- and intermolecular in-
teractions, especially those that involve hydrogen bonding.¹,² Even more than IR spectra,
interpreting differential VCD signals is challenging and often requires a thorough comparison
with theoretical models. In particular, accounting for environmental effects due to the pres-
ence of a solvent has been historically recognized as a major source of additional difficulties,
as is the case of cooperative effects arising in molecular crystals. However, recent advances
in computational modeling have paved the way to address these challenges.³,⁴ The most
common methods employed are static harmonic calculations at the density functional theory
level (DFT) incorporating continuum solvation models, as well as the "clusters-in-a-liquid"
model, which treats the solvent and solute as interconnected clusters in cases where strong
intermolecular interactions, such as hydrogen bonds, are present.²,⁵⁻⁷

Alternatively, molecular dynamics (MD) approaches are powerful tools for sampling the
conformational space of flexible molecular systems at finite temperature. In this context, IR
and VCD spectra can be determined without relying on the harmonic approximation but
from the Fourier transform of appropriate time correlation functions.⁸ Specifically, the IR
spectrum is calculated from the autocorrelation function of the electric dipole moment, while
the VCD spectrum is obtained from the cross-correlation function of the electric and magnetic
dipole moments.⁹ Within this framework, first principles MD simulations – particularly
when enhanced with nuclear velocity perturbation theory (NVPT) to obtain the required
magnetic dipole moments¹⁰⁻¹² – have shown robust performance in the interpretation of VCD
spectra, with applications in both liquid and solid phases.²,³,¹²,¹³ However, this method has
a high computational cost and is constrained by system size, time scales, as well as sampling
efficiency.

Another approach relies on MD simulations with classical force fields,⁸,¹⁴⁻¹⁶ although
care has to be taken to ensure sufficient accuracy for vibrational spectroscopy. In recent
years, advanced force fields have been developed to better account for important physical
effects, such as refined electrostatics and polarization forces. Properly parameterized polar-
izable force fields have been shown to effectively handle large systems with long sampling
times while preserving chemical accuracy.¹⁷⁻²¹ A key ingredient for VCD spectra is the mag-
netic dipole moment. The proposed methodology based on a polarizable force field extends
the classical treatment of magnetic dipole moments by incorporating several contributions
from the induced electric dipole moment and its time derivative,²² drawing inspiration from
previous work utilizing NVPT in VCD spectroscopy.¹,²³,²⁴

In this context, highly flexible molecules in polar solvents offer attractive testbeds for the-
oretical VCD spectroscopy. The case of 1-Phenyl-1,2-cyclohexanediol (PC, C₁₂H₁₆O₂), which
is permanently chiral, is particular as this molecule is often used as a unit for chiral station-
ary phases.²⁵ As a result, it has been the subject of a number of experimental studies²⁵⁻²⁷
for example, for the synthesis of homochiral ligands with applications in pharmaceuticals
and cosmetics. A detailed theoretical investigation on the spectral characteristics of weak
interactions between hydroxyl groups in vicinal diols in the gas phase was recently carried
out by Hartwig and Suhm.²⁸ These authors emphasized the impact of substitution and ring
effects on the vibrational frequencies. The role of the phenyl substituent was determined to
be essential in studies conducted under supersonic expansion conditions for the investigation
of the structure of PC.²⁹ Understanding these weak interactions is key for elucidating the
behavior of such systems in different environments, including solvent effects that can either
enhance or disrupt the hydrogen bonding network.³⁰,³¹

In this work, we use the recently proposed methodology²² for calculating vibrational
spectra that combines classical molecular dynamics simulations with the AMOEBA polar-
izable force field.³² Specifically, we implement this approach to investigate PC in dimethyl
sulfoxide (DMSO) solvent, leading to the determination of both IR and VCD spectra. The
modeling of these spectra sheds valuable light onto the molecular properties and in particular
the complex solute-solvent interactions.

The paper is organized as follows: Section 2 describes the computational strategy em-
ployed, presenting the AMOEBA force field, the methodology for extracting VCD spectra
from time correlation functions, and the details of the molecular dynamics simulations. Sec-
tion 3 is dedicated to the analysis of the results, beginning with a study of the structural
parameters of the PC solute and the organization of the DMSO molecules around it. The
IR and VCD spectra are then analyzed by examining the influence of the DMSO concentra-
tion and the need for a careful selection of electric and magnetic dipole moments to reveal
the PC signatures is critically discussed. Finally, some of the spectral bands are qualitatively assigned, in an attempt to compare with experimental data involving PC in the fully deuterated DMSO solvent.

# Computational strategy

## The AMOEBA force field

The AMOEBA (Atomic Multipole Optimized Energetics for Biomolecular Simulation) polarizable force field³² was used here to describe the molecular systems of interest through MD simulations. Originally developed for organic molecules and biomolecules in condensed phases,³³⁻³⁵ it has been shown to perform very well for vibrational spectra alike.¹⁷,²²,³⁶ It reproduces satisfactorily the relative conformational energies compared to other force fields, providing the required accuracy for structural and thermodynamic properties. This good agreement originates from the use of permanent electrostatic multipolar moments on each atom up to the quadrupole, as well as induced dipole moments $\vec{\mu}_i^{\text{ind}}$ reacting self-consistently to the presence of local electric fields on each polarizable site.

In practice, the induced dipole moment $\vec{\mu}_i^{\text{ind}}$ on atomic site i is related to the local electric field $\vec{E}_i$ acting on this site as:

$$ \vec{\mu}_i^{\text{ind}} = \alpha_i \vec{E}_i, \qquad (1) $$

with $\alpha_i$ being the atomic polarizability, assumed to be isotropic in AMOEBA. The electric field $\vec{E}_i$ has a contribution from the induced dipole moments on other sites, and the entire set of induced dipole moments must be solved self-consistently, which in AMOEBA is achieved iteratively.

For the systems presently studied, the major part of AMOEBA parameters are available from the literature, except for the multipoles of PC. Following the standard AMOEBA procedure,³⁷ a set of atomic multipoles was generated from the MP2/6-311G(d,p) wavefunction
using the distributed multipole analysis³⁸,³⁹ and then refined by electrostatic potential fitting at the MP2/aug-cc-pVTZ level. Furthermore, a comparative study of the relative energies and structures of the low-energy conformations of the isolated molecule confirms the good agreement between quantum chemistry results and those of the force field (see Table S1 and Figure S1 in SI). The set of AMOEBA parameters for the PC molecule are provided as Supplementary Information.

## Vibrational circular dichroism and infrared spectroscopy

IR spectroscopy measures the response of the electric dipole moment to nuclear motion. In practice, IR intensities are obtained from the Fourier transform of time autocorrelation functions of the electric dipole moment $\vec{\mu}$. Integration by parts straightforwardly shows that the IR spectral line shape can also be obtained from the time autocorrelation function of the time derivative of $\vec{\mu}$, also called the velocity form of the electric dipole moment:

$$ \dot{\vec{\mu}}(t) = \frac{d\vec{\mu}(t)}{dt}, \qquad (2) $$

thereby giving the absorption intensity line shape of the IR spectrum:

$$ A_{\dot{\mu}}(\omega) = \frac{2\pi\beta}{3Vcn(\omega)} \int_{-\infty}^{+\infty} \langle \dot{\vec{\mu}}(0) \cdot \dot{\vec{\mu}}(t) \rangle e^{-i\omega t} dt, \qquad (3) $$

in which $\omega$ denotes the frequency, $V$ the volume of the sample, $n(\omega)$ the refraction index of the medium at temperature $T = \frac{1}{k_B\beta}$ with $k_B$ the Boltzmann constant, and $c$ the velocity of light in vacuum. The form expressed in Eq. (3), based on the hypervirial theorem, avoids the origin dependence that arises from using the position form of the electric dipole moment alone.

The total electric dipole moment $\vec{\mu}$ contains the partial charges $q_i$ distributed at the
respective positions $\vec{r}_i$, and the induced dipoles $\vec{\mu}_i^{\text{ind}}$ at each polarizable site $i$:

$$ \vec{\mu}(t) = \sum_{i}^{N} [q_i \vec{r}_i(t) + \vec{\mu}_i^{\text{ind}}(t)] \quad (4) $$

where $N$ designates the number of particles carrying charges and induced dipoles in the simulation box.

VCD spectral intensities are obtained from the Fourier transform of the cross-correlation function between the electric dipole moment with the magnetic dipole moment $\vec{m}$. The resulting expression for the VCD intensity reads$^{8,14,15,40}$:

$$ \Delta A_{\dot{\vec{\mu}}}(\omega) = \frac{8\pi\beta\omega}{3Vcn(\omega)} \int_{-\infty}^{+\infty} \langle \dot{\vec{\mu}}(0) \cdot \vec{m}(t) \rangle e^{-i\omega t} dt. \quad (5) $$

The VCD intensity $\Delta A$ characterizes the rotation of the electric current around the axes of propagation of the circularly polarized light, which is sensitive to chiral information. $\vec{m}$ itself originates from the circular motion of the electric charges and dipole moments around a common reference point. A classical approximation is used for the present purely atomistic model provided by the AMOEBA force field, in which $\vec{m}$ is defined as the sum of the moving point charges, with additional contributions from the fixed and induced electric dipole moments in their motion, as well as an extra contribution from the time-varying electric dipole moment moving around fixed positions:$^{22}$

$$ \vec{m}(t) = \sum_{i=1}^{N} \left[ \frac{1}{2c} q_i \vec{r}_i(t) \times \vec{v}_i(t) + \frac{1}{c} \vec{\mu}_i^{\text{ind}}(t) \times \vec{v}_i(t) + \frac{1}{2c} \vec{r}_i(t) \times \dot{\vec{\mu}}_i^{\text{ind}}(t) \right], \quad (6) $$

where $\vec{v}_i(t)$ denotes the velocity vector of atom $i$ at time $t$.

## Simulation details

Molecular dynamics trajectories based on the AMOEBA force field were propagated with a modified version of the Tinker software package,$^{41}$ for the isolated PC molecule, its com-
plex with one solvent molecule (1:1 PCDMSO), for PC embedded into 50 (PC50DMSO) or 140 (PD140DMSO) solvent molecules, and for a box of 50 pure DMSO solvent molecules (50DMSO). The different boxes were constructed using the xyzedit tool available in Tinker by combining the two coordinate files, one from the PC molecule and one from the DMSO equilibrated box. A local optimization of the final box was performed to relax the system before starting the MD simulations. To build the PC50DMSO system, the PC molecule was inserted into the 50DMSO box after removing 4 DMSO molecules.

Initially, a 1 ns (10 ns in the PC50DMSO case) trajectory was generated for each system described in Table 1 within the canonical (NVT) ensemble at 298 K, employing the Nosé-Hoover thermostat with a 1 ps coupling time to equilibrate the initial configurations. Phase space configurations were selected every 10 ps of the last 400 ps equilibrated trajectory. The purpose of this canonical trajectory is twofold: one is to analyze the structure of the system at equilibrium (see Section 3.1), and the second is to collect a set of configurations along the trajectory that can be used as starting points for calculating and averaging IR and VCD spectra. Both IR and VCD intensities are absolute intensities.

Subsequently, each starting point underwent an additional but shorter NVT trajectory for equilibration, followed by a microcanonical (NVE) trajectory without the thermostat and from which the electric and magnetic dipole moments were recorded every time step as a function of time. Each microcanonical trajectory results in one IR spectrum and one VCD spectrum. Both NVT and NVE simulations employed a time step of 0.5 fs and were propagated for 50 ps. Such a time step had been previously validated as sufficiently short and suitable for the vibrational range of 0–4000 cm⁻¹.¹⁷,²²,⁴² The total number of independent trajectories was set at 40 for all systems, corresponding to a total simulation time of 2 ns per spectrum, with occasionally longer simulations of 80 trajectories for convergence evaluation purposes (see details in Table 1). Besides normal DMSO solvent, additional simulations with fully deuterated DMSO were carried out with the purpose of comparing the spectra with experimental measurements performed with this solvent, the resulting system being denoted
as PC50DMSOd (vide infra).

The electric and magnetic dipole moments recorded along the NVE trajectories were processed using the ChirPy program.¹²,⁴³ The issue of the origin dependence of the magnetic dipole moment has been solved by considering the nearest periodic image between molecular units, following the approach introduced by Jähnigen and coworkers.⁴⁴ The resulting IR and VCD spectra were averaged over the entire set of individual spectra, providing reference quantities from which global convergence can be statistically assessed using appropriate error measures.

<table><caption>**Table 1:** Computational details of the simulated systems, all except 50DMSO containing a single PC molecule in the simulation box. PC50DMSOd indicates that the DMSO molecules are all fully deuterated. Simulation times are given for the first NVT trajectory from which initial configurations are generated (ns), for the NVT trajectory used to equilibrate these starting points (ps), and for the subsequent NVE trajectory (ps).</caption><thead><tr><th>System</th><th>Box length (Å)</th><th>Number of atoms</th><th>Density (g/cm³)</th><th>Number of DMSO(d) molecules</th><th>Simulation time (ns, ps, ps)</th><th>Number of trajectories</th></tr></thead><tbody><tr><td>Isolated PC</td><td>(no box)</td><td>30</td><td>gaseous phase</td><td>0</td><td>1, 50, 50</td><td>40</td></tr><tr><td>1:1 PCDMSO</td><td>(no box)</td><td>40</td><td>gaseous phase</td><td>1</td><td>1, 50, 50</td><td>40</td></tr><tr><td>PC50DMSO(d)</td><td>18</td><td>490</td><td>1.57 (1.65)</td><td>46</td><td>1 (10), 50, 50</td><td>40, 80</td></tr><tr><td>PC140DMSO</td><td>24</td><td>1430</td><td>1.34</td><td>140</td><td>1, 50, 50</td><td>40, 80</td></tr><tr><td>50DMSO</td><td>18</td><td>500</td><td>1.11</td><td>50</td><td>1, 50, 50</td><td>40</td></tr></tbody></table>

The 1:1 PCDMSO system was further studied at the DFT level. Following previous works,²,²⁹ the conformational search was performed in the MacroModel program.⁴⁵ The lowest energy conformations obtained in a 10 kcal.mol⁻¹ window were subsequently optimized at the B3LYP/6-311++G(d,p) level associated with an implicit polarizable continuum model to account for the DMSO solvent environment. The procedure resulted in 6 structures within the 1.1 kcal.mol⁻¹ energy range (see Table S2 and Figure S2 in the Supplementary Material). Harmonic IR and VCD spectra were computed at the same level of theory, and the resulting spectra were convoluted using appropriate Boltzmann weights (vide infra).
Experimental methods

The IR absorption and VCD spectra were measured using a Fourier transform infrared
(FTIR) spectrometer Vertex 70 equipped with a VCD module PMA 50 (Bruker). The signal
was measured by a MCT IR detector cooled with liquid nitrogen and equipped with a BaF₂
window. The spectral resolution used for both absorption and VCD spectra was 4 cm⁻¹.
The IR radiation was first spectrally filtered by a 2000 cm⁻¹ low-pass filter then polarized
with a linear polarizer. The linearly-polarized light was modulated by a 50 kHz ZnSe photo-
elastic modulator (Hinds). The signal of the MCT detector was demodulated and amplified
by means of a lock-in amplifier (Stanford Research Systems SR 830).

The spectra of DMSO-d6 solutions were measured in a 56 µm cell with an acquisition
time of 3 hours using concentrations ranging between 0.1 and 1 M. The largest concentration
used corresponds to the saturation limit at room temperature and to an optical density of
~0.7 for the larger-intensity band of the fingerprint region at 1446 cm⁻¹. The experimen-
tal conditions were controlled by optimizing the position of the cell to minimize artefacts
due to birefringence and by checking the mirror-image relation between the VCD spectra
of the two enantiomers of pinene (neat liquid). DMSO-d6 was purchased by Aldrich and
used without further purification. Both enantiomers of PC were obtained from Santa Cruz
Biotechnology (99% purity; enantiopurity >98%). The spectra presented below correspond
to a concentration of 0.55 M, which provided the best signal to noise ratio.

**Results and discussion**

Structural properties

The minima of PC in gas phase were explored via MD simulations. Among 1000 snap-
shots optimized at the AMOEBA level and reoptimized at the DFT level, only five distinct
conformers were identified. Their structures and relative energies are reported in SI. The
accuracy of the AMOEBA force field for determining the PC conformational landscape was assessed by comparing the energy ranking obtained by static calculations at both force field and B3LYP/6-311++G(d,p) levels (see Table S1). The latter level reproduced well the experimental results in terms of most stable structure and frequencies.²⁹ All the stable calculated structures listed in Table S1 have the cyclohexane ring in a chair geometry and the phenyl group in equatorial position, with the O₂H₂ group in equatorial position, and the O₁H₁ group in axial position. The two methods give very similar relative energies. All of the listed structures are stabilized by a sterically constrained OH···O intramolecular interaction but differ by the direction of the hydrogen bond (O₁H₁···O₂ or O₂H₂···O₁ interaction) and a rotation of the accepting OH. The most stable form, to which the experimentally observed structure is assigned, shows an O₁H₁···O₂ interaction accompanied by a weak O₂H₂···π (see Figure S1).

From the 10 ns NVT simulation of PC solvated in a box of 50 DMSO molecules (referred to as PC50DMSO in Table 1) at 298 K, the basic structural features can be inferred from looking at the solute itself, its first solvation shell, or structural data, primarily radial distribution functions (RDFs), $g(r)$, such that $g(r)dr$ measures the probability density of finding specific pairs at a distance between $r - dr/2$ and $r + dr/2$ (see Figure S3). Concerning the PC molecule itself, the cyclohexane ring is found to adopt a chair conformation with the phenyl group being positioned equatorially, thereby minimizing steric constraints. The two hydroxyl groups are in the cis configuration, allowing for the possible formation of an intramolecular hydrogen bond that stabilizes the structure between the O₁ and H₂ atoms, separated by 2.5 Å from each other, in average.

Likewise, the arrangement between the solute and the immediate solvent molecules can also be described from typical configurations visited along the NVT trajectory. We define the hydrogen bonds by labeling as O₁H₁ and O₂H₂ the two hydroxyl groups of the PC solute that are each hydrogen bonded with a different DMSO molecule. In Figure 1, these intermolecular hydrogen bonds are characterized by the distances $d_1$ and $d_2$ and the angles
$\Theta_1 = \widehat{O_1H_1O_{DMSO}}$ and $\Theta_2 = \widehat{O_2H_2O_{DMSO}}$, with $O_{DMSO}$ denoting the accepting oxygen atom of the DMSO molecule.

<figure><img src="image_2.png" /><figcaption>Figure 1: The two hydrogen bonding interactions between the PC solute molecule and one DMSO solvent molecule: (a) bonding to $O_1H_1$; (b) bonding to $O_2H_2$.</figcaption></figure>

Examining the radial distribution function of the oxygen atoms $O_1$ and $O_2$ from PC with the oxygen atoms $O_{DMSO}$ of all DMSO molecules (see Figure S3) reveals distinct peaks at 2.95 Å and 2.85 Å, respectively, as well as around 5.35 Å for $O_1$. However, the latter peak is absent from the RDF associated with $O_2$. The peak at the shortest distance of 2.95 Å or 2.85 Å confirms the interaction of DMSO molecules with the hydroxyl groups in the first solvation shell of the PC solute. However, the integrated RDF $N(r) = \int_0^r 4\pi s^2 g(s) ds$, also shown in Figure S3, indicates that only a fraction of DMSO molecule is found at this distance, implying some fluxionality in the solute-solvent interaction and the occasional breaking of the associated hydrogen bond. The peak at 5.35 Å for $O_1$ is related to the intramolecular hydrogen bond. This further suggests an oscillation between the two -OH terminations, with a stronger preference for the $O_2H_2$ group, presumably due to steric effects and partial retention of the intramolecular interaction.
<figure><img src="image_3.png" /><figcaption>Figure 2: Combined distribution functions between the angles Θ₁ and Θ₂ and the distances d₁ and d₂, (a) for OH₁; and (b) for OH₂. Red boxes highlight the hydrogen bond regions (see the text for details).</figcaption></figure>

A far more complete picture of the hydrogen bonding pattern between PC and DMSO solvent molecules is provided by the combined distribution function (CDF) which measures the propensity for solvent molecules to be hydrogen-bonded as oxygen donors to the two anchor sites O₁H₁ and O₂H₂ of PC, as a double function of the OₖHₖ distance (k = 1 or 2) with the associated angle Θₖ of interest. Figure 2 shows such joint distributions, obtained from considering the DMSO molecules identified as being within 10 Å of the PC solute. Below 6 Å, the CDF shows two areas of high density around Θ₁ = 100–140° and d₁ = 4.5–5.2 Å for O₁H₁ [Figure 2(a)], and Θ₂ = 145–180° and d₁ = 1.7–2.0 Å for O₂H₂ [Figure 2(b)]. These values allow setting the following criteria for hydrogen bonding, namely that it should fulfill conditions on the angle to be at least 130° and the O···H distance to be at most 3.5 Å. The DMSO molecules maintain the hydrogen bond with O₂H₂ during a significant part of the trajectory but also involve O₁H₁ as a secondary interaction. To better interpret this result, the residence time associated with the hydrogen bond between O₂H₂ and the DMSO molecules was evaluated along the trajectory, from the time series of 50 fs and their statistics. The residence time refers to the continuous time a molecule spends within a predefined area of 5 Å around the PC molecule and the corresponding results are provided in the SI as
Figure S4. A mean residence time of 107 ps was determined, after discarding short events of less than 2 ps, with occasional events persisting up to 200 ps, thereby confirming the highly fluxional nature of the DMSO-O₂H₂ interactions.

## Extracting the spectroscopic contribution of the solute from the signal

Following the methodology described in Section 2, classical MD simulations were performed for the PC molecule in a DMSO box in order to collect the time-resolved dipole moments of the system during the simulation and to calculate the IR and VCD spectra. The aim is to understand the influence of the solvent on the spectroscopic signals, which implies being able to decouple the signatures of the PC solute from those of the DMSO solvent molecules. This issue is particularly acute when the spectroscopic response of the solvent falls onto the same range as the one of the solute. One natural way to address it, which has been used in the past,²²,²³,⁴⁰ consists of selecting only a subset of solvent molecules around the solute in the evaluation of dipole moments.

In Figure 3, the IR and VCD spectra of the PC50DMSO system are shown using different conditions in the selection process. Figure 3(a) corresponds to the spectra obtained for the entire simulation box, i.e. without any selection. The IR spectrum is found to be quite comparable to the corresponding spectra determined for the pure solvent box, without a PC molecule (50DMSO system in Table 1), shown as Supplementary Information in Figure S5. For both systems, four main spectral regions of interest can be identified as 900–1100 cm⁻¹, ~1500 cm⁻¹, ~3100 cm⁻¹ and 3450–3700 cm⁻¹. These bands appear consistent with the experimental frequencies of liquid DMSO⁴⁶ (see Table S3), confirming the reliability of the present computational methodology to describe the DMSO solvent. The VCD spectra show greater differences between PC50DMSO and 50DMSO, as expected for a solvent composed of achiral molecules. The PC50DMSO spectrum displays notable activity around 1500 cm⁻¹, and even more pronounced activity between 3000 and 3700 cm⁻¹. The deviation from the
VCD spectrum obtained for the pure, achiral solvent is another indication that VCD spectra are more difficult to converge than IR spectra, as well known from earlier work.²² This was anticipated as the average of individual VCD spectra over MD trajectories, and the average itself in the right hand side of Eq. (5), involve quantities with positive and negative signs in similar numbers. The comparison between the spectra obtained for 50DMSO and PC50DMSO shows that, at this level of dilution, no clear signal specific to the PC solute can be inferred from the simulations if all dipole moments are to be included when calculating the spectra.

To smoothly remove the deleterious contribution from the solvent molecules when evaluating the IR and VCD spectra from the appropriate dipole moments, we scale the contribution of each solvent molecule by an exponentially decaying function of its distance to the PC solute. This radial cutoff function, evaluated from the center of mass of the PC molecule, is implemented in the ChirPy program as post processing of the simulation.⁴⁷ In practice, the decaying exponent is chosen so as to select appropriate solvation shells around the solute.

For the present simulations, the spectra were determined for a vanishing cutoff radius (thus accounting solely for the PC molecular response), for a cutoff radius of 5 Å (corresponding to DMSO molecules in the first solvation shell), or without any cutoff (entire box), as represented in Figure S6. Here it should be stressed again that varying the cutoff radius only affects the post-processing of the MD simulations, the same trajectories being used but with different contributions of the various molecules to the dipole moments.

The resulting spectra are represented in Figure 3. On the IR spectra computed with a vanishing or finite cutoff, peaks at 1100–1400 cm⁻¹, ~1600 cm⁻¹, and 3000–3250 cm⁻¹ are predominant. By removing the contribution from most of the solvent to the dipole moments, the intensity of the signal has decreased by a factor of almost 40, magnifying smaller peaks that were barely distinguishable on the spectrum obtained for the entire simulation box.

Despite their similarities, the spectra obtained for PC alone and with a 5 Å cutoff tend to differ as well, only the latter exhibiting an intense narrow peak near 3000 cm⁻¹, suggesting
that this is an active mode of the DMSO solvent.

As far as the VCD spectra are concerned, the differences with and without cutoff are also very clear. However such is not the case between the spectra obtained with finite and vanishing cutoff radii that are very similar to each other. Both spectra show a significant number of peaks, which were hidden by the solvent bands in the spectrum calculated without cutoff. This should help to highlight more PC vibrational modes, for example in the 1400-2000 cm⁻¹ range. Based on the spectra analyzed thus far, the interpretation of the VCD signal at this stage needs more investigation. Unless otherwise stated, we will now consider spectra calculated with a cutoff of 5 Å.

To better understand the influence of the various parameters on the calculated spectra, the effect of the box size was varied by inserting the PC molecule into a much larger box of DMSO, preserving the total volumetric density to be around 1.3–1.5 g.cm⁻³ but reducing effectively the PC molar density by a factor of 42%. The resulting PC140DMSO simulations were processed similarly as with the PC50DMSO system, and the IR and VCD spectra obtained using a 5 Å cutoff radius are represented in Figure 4. Evidently, both systems display similar IR profiles, but the peaks in the PC50DMSO spectrum are slightly better resolved. Moreover, larger bands and slight blueshifts can be noticed, particularly at 1000, 1500, and 3100 cm⁻¹ in the PC140DMSO spectrum, with respect to the smaller system.

Despite such overall similarities, the VCD spectra display quite stronger differences. The band around 1500 cm⁻¹ is better resolved and split into two bands for PC140DMSO, whereas the bands around 3500 cm⁻¹ present a sign inversion between the two systems. Such differences could naturally be caused by the apparently higher concentration in PC in the smaller box, however we cannot exclude convergence issues in the spectra themselves at this stage.
<figure><img src="image_4.png" /><figcaption>Figure 3: IR and VCD spectra of PC50DMSO obtained from MD trajectories at 298 K, determined by varying the contribution of solvent molecules around the PC solute depending on their distance. (a) No cutoff (entire box considered); (b) cutoff radius of 5 Å; (c) vanishing cutoff (no DMSO).</figcaption></figure>

## Statistical convergence of the spectra

As suggested by Eqs. (3) and (5) themselves, the IR and VCD spectra result from statistical accumulation of time series over a sufficient number of initial conditions of the MD trajectories. The global convergence of the IR and VCD spectra can be quantified by defining appropriate error measures when using a subset of the available data, accumulating the
<figure><img src="image_5.png" /><figcaption>Figure 4: IR and VCD spectra of (a) PC50DMSO and (b) PC140DMSO, both obtained from 40 MD trajectories at 298 K and processed with a cutoff radius of 5 Å when evaluating the contribution of solvent molecules to the dipole moments. For both spectra, the average spectra are shown with thick lines, standard deviations being superimposed as light backgrounds.</figcaption></figure>

time correlation functions as described previously²² but assuming in the absence of higher knowledge that the entire set of available data provides the reference spectra to be compared with. The standard deviation around the average IR spectra, superimposed in Figure 4, is generally proportional to the IR intensity. However, the corresponding deviation for the VCD spectra is much higher, suggesting much slower convergence of this quantity.

Following Ref. 22, convergence relative to these reference spectra is quantified as the number of MD trajectories included in the sample increases towards the entire set, by comparing the reference spectrum with partial spectra inferred from including only a random subset of the entire sample. The error $E_{IR}$ between two IR spectra is calculated by the quadratic difference, normalized by the square sum of individual spectra. VCD spectra, which exhibit positive and negative signals, require a different measure, provided here by
the SimVCD method,⁴⁸ and using more specifically 1-SimVCD in such a way that identical spectra give a vanishing error $E_{VCD}$. The two errors for the IR and VCD spectra are thus formally defined as

$$E_{\text{IR}} = \frac{A_{12}}{A_{11} + A_{22}}, \qquad (4.3)$$

and

$$E_{\text{VCD}} = 1 - \frac{Z_{12}}{Z_{11} + Z_{22} - |Z_{12}|} = 1 - \text{Sim VCD}, \qquad (4.4)$$

where $A_{ij}$ is the quadratic difference between the IR spectra i and j, $A_{ii}$ their quadratic norm, $Z_{ij}$ the overlap integral between the VCD spectra i and j (also for $i = j$):

$$A_{ij} = \int |I_i(\omega) - I_j(\omega)|^2 d\omega, \qquad (4.5)$$

$$A_{ii} = \int |I_i(\omega)|^2 d\omega, \qquad (4.6)$$

$$Z_{ij} = \int \text{VCD}_i(\omega) \cdot \text{VCD}_j(\omega) d\omega, \qquad (4.7)$$

$$Z_{ii} = \int |\text{VCD}_i(\omega)|^2 d\omega, \qquad (4.8)$$

and with $I_i(\omega)$ being the IR intensity at frequency $\omega$ and $\text{VCD}_i(\omega)$ the VCD intensity.

Given a number of time series provided by individual MD trajectories, a Monte Carlo procedure is employed to generate partial spectra from a random subset and compare them to the reference spectrum, resulting in a specific error $E_{\text{IR}}$ or $E_{\text{VCD}}$. By repeating this process 10,000 independent times, average errors can be calculated as a function of the number of individual MD trajectories.²²

Figure 5 shows the convergence of the IR and VCD spectra for the PC50DMSO and PC140DMSO systems as a function of the total simulation time covered by the subset of trajectories out of 40 and 80 in the entire sets, respectively. The PC50DMSO system shows a faster decrease of the IR error compared to PC140DMSO (Figure 5.a), thus indicating that larger simulation boxes also require longer simulation times to converge the IR spectra.
The VCD spectra generally display a slower convergence compared to the IR spectra for both system sizes. The associated error also decreases faster for the smaller system, confirming the trend noticed just above for the IR signal. In all cases, a total simulation time of 2 ns (40 trajectories) seems to be sufficient to converge the IR and VCD spectra for these systems.

<figure><img src="image_6.png" /><figcaption>Figure 5: Normalized error of the average in the IR and VCD spectra of PC50DMSO and PC140DMSO systems for 40 and 80 MD trajectories of 50 ps each, performed at 298 K. (a) and (b) IR spectra; (c) and (d) VCD spectra. For all errors, the maximum number of trajectories is indicated between brackets.</figcaption></figure>

Qualitative attribution and comparison with experiment

Various vibrational modes can be tentatively assigned to specific molecular groups, exploiting mainly the IR spectra but also, to some useful extent, the VCD spectra. Figure 6 compares
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 6: IR and VCD spectra in the 500–2500 cm⁻¹ range, obtained from 40 MD trajectories at 298 K. (a), (b) and (c) IR spectra; (e), (f) and (g) VCD spectra of PC50DMSO (5 Å cutoff radius), isolated PC and 1:1 PCDMSO (no cutoff), and pure 50DMSO box (no cutoff), respectively. Colored vertical bands are included as visual guides, with the yellow and purple zones corresponding to the DMSO and PC vibrational modes, respectively.</figcaption></figure>
the spectra obtained for PC50DMSO (5 Å cutoff), isolated PC, 1:1 PCDMSO and 50DMSO, within the range of 500–2500 cm⁻¹, giving insight into the distinct roles of molecular interactions and solvent effects on the spectral features of the PC molecule. The spectra of pure DMSO have been obtained without a cutoff.

As usual, the IR bands potentially carry less information but are more robust and thus discussed first. In the PC50DMSO system, the most intense peak lies near 1500 cm⁻¹ and matches a peak in the 50DMSO spectrum. This signal is less intense in the 1:1 PCDMSO spectrum, as it is likely hidden by the PC intrinsic signals, and no clear peak appears near this wavenumber for the isolated solute molecule. This band is characteristic of the umbrella mode of the CH₃ groups in DMSO.

The peaks between 1550 and 1750 cm⁻¹ in the IR spectrum are present in PC50DMSO, 1:1 PCDMSO and isolated PC spectra with different intensities, but they are absent in 50DMSO. They are better defined for the isolated system, and are thus characteristic of the PC molecule. These peaks can be associated with the OH bending modes but also with the stretching mode of the aromatic C=C groups of PC. The peaks around 1300 cm⁻¹ are present in PC50DMSO, 1:1 PCDMSO and isolated PC spectra and are characteristics of the CH₂ bending modes of PC. Additionally, beyond 2250 cm⁻¹, all spectra with a PC molecule show a bump that is absent from the 50DMSO solvent spectrum. This feature is characteristic of the stretching mode of the CH₂ groups of PC. Following the same reasoning, the peaks around 1000 cm⁻¹ can be attributed to the contribution of the S=O stretching mode of the DMSO molecule, as well as the C-C-C and C-O-H deformation modes of the PC molecule.

The VCD spectrum exhibits relatively similar profiles for the systems that contain a PC molecule, although the intensity of the signal generally decreases with increasing DMSO concentration. The 50DMSO spectrum shows a very weak signal in this frequency range, with the exception of the band group around 1500 cm⁻¹, which exhibits a slightly more intense signal. Near 1000 cm⁻¹, a sign inversion is found between the spectra of PC50DMSO, 1:1 PCDMSO, and that of isolated PC. This suggests that DMSO influences the VCD signal of
the system, as the sign depends on its presence: at higher DMSO concentration, the sign tends to become positive at $1000 \text{ cm}^{-1}$.

<figure><img src="image_13.png" /><figcaption>Figure 7: IR (left) and VCD (right) spectra calculated for PC50DMSOd and 1:1 PCDMSO in an implicit solvent at the DFT level (see computational details), and experimentally measured for PC solvated by fully deuterated DMSO. The IR and VCD spectra of PC50DMSOd were calculated using the Welch method with a window length of 800 fs.⁴³ The PC50DMSOd and DFT spectra were scaled by 0.91 and 0.98, respectively. For better readability, the VCD spectra are vertically shifted, and the dashed horizontal lines indicate the zero baseline.</figcaption></figure>

Experimental spectra were recorded in the 1100–1700 $\text{cm}^{-1}$ range in a deuterated DMSO solvent, and additional simulations were thus performed for this system, using 50 solvent molecules in the PC50DMSOd system listed in Table 1. The calculated and measured spectra are compared with each other in Figure 7. Remarkably, the convergence of the IR and VCD spectra for the deuterated solvent occurs significantly faster than that of the normal solvent, only about 15 trajectories appearing sufficient (Figure S7). For the present classical systems, the structural details for the deuterated and non deuterated solvents are very similar under the same simulation conditions, as shown from the comparison of the radial distribution functions in Figure S3.

The simulated spectrum is broader than the experimental one (see Table S4 for temperature fluctuations during the NVE simulations), which could probably be corrected by using a
lower temperature, but this would be detrimental to the convergence of the simulation owing
to different isomer statistics. Moreover, theoretical spectra were scaled to correct their sys-
tematic shift from the experimental data. Despite the broader signal, one can notice a good
correspondence between the main features observed in the IR experimental and computed
spectra, which confirms the general reliability of our modeling. Furthermore, static spectra
of the 1:1 PCDMSO system obtained from our DFT calculations including implicit solvent
effects are also reported in Figure 7. It is worth noting that the IR DFT spectrum accurately
reproduces the main features of the experimental spectrum.

The combination of theoretical and experimental data allows the main bands to be as-
signed. The transitions at ~1200 cm⁻¹ correspond to superimposed CH₂ bending modes of
PC. The most intense features are calculated at 1450-1550 cm⁻¹; they correspond to the
experimental features where the coupled OH and CH bending modes of PC are expected.
The bands in the range 1600 to 1700 cm⁻¹ are found to be more intense in the MD spectrum
than in DFT and correspond to the experimental bands assigned to in-plane aromatic CH
bending modes.

Concerning the VCD spectra, the signal pattern is only partly reproduced: the bands are
predominantly negative in the 1150-1300 cm⁻¹ region and positive in the 1300-1500 cm⁻¹
region. The DFT spectrum shows a higher level of resolution, enabling most of the exper-
imental features to be captured. The similarity between the MD spectra in DMSO and
DMSOd confirms the molecular origin of the signal as being mainly due to the PC molecule.

Conclusions and perspectives

In this work, we challenged and further validated a previously developed method for obtaining
IR and VCD spectra using classical molecular dynamics simulations with a polarizable atomic
multipole force field. This approach was applied to a flexible chiral molecule in solution,
specifically 1-phenylcyclohexane-1,2-diol in dimethyl sulfoxide, under different simulation
conditions.

We first analyzed the structural properties obtained from NVT trajectories of PC in DMSO. The radial distribution function analysis revealed a preferential clustering of DMSO molecules around the $O_2H_2$ group of PC. This was also supported by the combined distribution function of the angle between the PC hydroxyl group and the oxygen atom of DMSO, indicating that a single DMSO molecule predominantly, but not exclusively, forms a strong hydrogen bond with the $O_2H_2$ group.

The calculated IR spectrum of DMSO in solution showed correct agreement with experimental data after applying an appropriate scaling factor. Comparing the IR spectra of PC in a DMSO box to that of pure DMSO suggests that solvent molecules hide the PC signal. To isolate the PC contribution, a cutoff function was applied, encapsulating at 5 Å the first solvation shell, thus providing a reliable description of the vibrational spectra. Further comparisons with simulated gas phase PC including static DFT calculations of the 1:1 complex under implicit DMSO solvent, allowed for a more quantitative interpretation of the IR spectrum of PC in DMSO. However, interpreting in depth the experimental VCD spectrum was found to remain more challenging and would definitely require deeper investigations.

Larger systems with more DMSO molecules were also explored in order to shed light onto the possible size and concentration effects on the IR and VCD spectra. These systems were found to require more trajectories for accurate spectral calculations and it was also confirmed²² that VCD spectra demand a higher cumulated simulation time to reach convergence compared to IR spectra. Such difficulties are partly associated with the greater range of variation of VCD spectra, with their positive and negative peaks but also, more intrinsically, to their differential nature with respect to IR spectra that are already not straightforwardly converged.

Despite the noted difficulties, the present methodology has the virtue of being able to treat strongly anharmonic systems such as molecular liquids. It would then be valuable to use it for even condensed phase systems for which structural characterization is known to
be challenging. Chiral ionic liquids (ILs) and deep eutectic solvents (DES), among which pyrrolidinium-based ILs⁴⁹ and menthol-based DES⁵⁰ that have already received some attention from the computational side, could particularly benefit from the presently used tools. These green solvents are known for their exceptional properties but pose significant challenges in the understanding of their physical chemistry properties.

## Acknowledgement

This work was supported by the ANR DICHROSOL, Grant No. ANR-23-CE29-0006 of the French Agence Nationale de la Recherche. This collaborative work was also supported by GDR 2035 SolvATE and GDR 3533 EMIE, funded by CNRS. This work was performed using computational resources from the “Mésocentre” computing center of Université Paris-Saclay, CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Ile-de-France. This work was performed using HPC resources from the GENCI (CINES/IDRIS, Grant No. A0140806830).

## Supporting Information Available

Supporting Information presents the following material: relative energies and structures of the low-energy conformations of PC in a gas phase, relative energies and structures of the 1:1 PCDMSO conformations, radial distribution functions, variation of distances for hydrogen bonded molecules to PC, IR and VCD spectra for pure liquid DMSO, assignment of IR frequencies for bulk DMSO, representation of the cutoff region to scale dipole moments, convergence of the IR and VCD spectra of the PC50DMSOd system, average temperature and standard deviation for NVE trajectories. The coordinates of the initial simulation boxes are provided as individual files.
References

(1) Nafie, L. A. Infrared and Raman vibrational optical activity: theoretical and experi-
mental aspects. *Annu. Rev. Phys. Chem* **1997**, *48*, 357–386.

(2) Le Barbu-Debus, K.; Bowles, J.; Jähnigen, S.; Clavaguéra, C.; Calvo, F.; Vuilleumier, R.; Zehnacker, A. Assessing cluster models of solvation for the description of vibrational circular dichroism spectra: synergy between static and dynamic approaches. *Phys. Chem. Chem. Phys.* **2020**, *22*, 26047.

(3) Bloino, J.; Jähnigen, S.; Merten, C. After 50 Years of Vibrational Circular Dichroism Spectroscopy: Challenges and Opportunities of Increasingly Accurate and Complex Experiments and Computations. *J. Phys. Chem. Lett.* **2024**, *15*, 8813–8828.

(4) Giovannini, T.; Olszòwka, M.; Cappelli, C. Effective fully polarizable QM/MM approach to model vibrational circular dichroism spectra of systems in aqueous solution. *J. Chem. Theory Comput.* **2016**, *12*, 5483–5492.

(5) Weirich, L.; Merten, C. Solvation and self-aggregation of chiral alcohols: how hydrogen bonding affects their VCD spectral signatures. *Phys. Chem. Chem. Phys* **2019**, *21*, 13494–13503.

(6) Perez-Mellor, A.; Le Barbu-Debus, K; Zehnacker, A. Solid-state synthesis of cyclo LD-diphenylalanine: A chiral phase built from achiral subunits. *Chirality* **2020**, *32*, 693–703.

(7) Perera, A. S.; Thomas, J.; Poopari, M. R.; Xu, Y. The clusters-in-a-liquid approach for solvation: new insights from the conformer specific gas phase spectroscopy and vibrational optical activity spectroscopy. *Front. Chem.* **2016**, *4*, 9.

(8) Calvo, F. Atomistic Modeling of IR Action Spectra Under Circularly Polarized Electromagnetic Fields: Toward Action VCD Spectra. *Chirality* **2015**, *27*, 253.
(9) Frömbgen, T.; Drysch, K.; Tassaing, T.; Buffeteau, T.; Hollóczki, O.; Kirchner, B. Induced Chirality and Vibrational Optical Activity in an Ionic-Liquid Anion. *Angew. Chem. Int. Ed.* **2025**, e202502885.

(10) Scherrer, A.; Agostini, F.; Sebastiani, D.; Gross, E. K. U.; Vuilleumier, R. Nuclear velocity perturbation theory for vibrational circular dichroism: An approach based on the exact factorization of the electron-nuclear wave function. *J. Chem. Phys.* **2015**, *143*, 074106.

(11) Jähnigen, S.; Scherrer, A.; Vuilleumier, R.; Sebastiani, D. Chiral crystal packing induces enhancement of vibrational circular dichroism. *Angew. Chem. Int. Ed.* **2018**, *57*, 13344–13348.

(12) Jähnigen, S. Vibrational Circular Dichroism Spectroscopy of Chiral Molecular Crystals: Insights from Theory. *Angew. Chem. Int. Ed.* **2023**, *62*, e202303595.

(13) Jähnigen, S.; Le Barbu-Debus, K.; Guillot, R.; Vuilleumier, R.; Zehnacker, A. How Crystal Symmetry Dictates Non-Local Vibrational Circular Dichroism in the Solid State. *Angew. Chem. Int. Ed.* **2023**, *62*, e202215599.

(14) Yang, S.; Cho, M. Direct calculations of vibrational absorption and circular dichroism spectra of alanine dipeptide analog in water: Quantum mechanical/molecular mechanical molecular dynamics simulations. *J. Chem. Phys.* **2009**, *131*, 135102.

(15) Abbate, S.; Longhi, G.; Kwon, K.; Moscowitz, A. The use of cross-correlation functions in the analysis of circular dichroism spectra. *J. Chem. Phys.* **1998**, *108*, 50–62.

(16) Horníček, J.; Kaprálová, P.; Bouř, P. Simulations of vibrational spectra from classical trajectories: Calibration with ab initio force fields. *The Journal of Chemical Physics* **2007**, *127*, 084502.
(17) Semrouni, D.; Sharma, A.; Dognon, J.-P.; Ohanessian, G.; Clavaguéra, C. Finite Temperature Infrared Spectra from Polarizable Molecular Dynamics Simulations. *J. Chem. Theory Comput.* **2014**, *10*, 3190–3199.

(18) Schultheis, V.; Bernhard, R. R.; Schropp; Tavan, P. A Polarizable Force Field for Computing the Infrared Spectra of the Polypeptide Backbone. *J. Phys. Chem. B* **2008**, *112*, 12217–12230.

(19) Liu, J.; Miller, W. H.; Fanourgakis, G. S.; Xantheas, S. S.; Imoto, S.; Saito, S. Insights in quantum dynamical effects in the infrared spectroscopy of liquid water from a semiclassical study with an ab initio-based flexible and polarizable force field. *J. Chem. Phys.* **2011**, *135*, 244503.

(20) Koner, D.; Salehi, S. M.; Mondal, P.; Meuwly, M. Non-conventional force fields for applications in spectroscopy and chemical reaction dynamics. *J. Chem. Phys.* **2020**, *153*, 010901.

(21) Jansen, T. L. C. Linear Absorption and Two-Dimensional Infrared Spectra of N-Methylacetamide in Chloroform Revisited: Polarizability and Multipole Effects. *J. Phys. Chem. B* **2014**, *118*, 8162–8169.

(22) Bowles, J.; Jähnigen, S.; Agostini, F.; Vuilleumier, R.; Zehnacker, A.; Calvo, F.; Clavaguéra, C. Vibrational circular dichroism spectroscopy with a classical polarizable force field: alanine in the gas and condensed phases. *ChemPhysChem* **2024**, *25*, e202300982.

(23) Scherrer, A.; Vuilleumier, R.; Sebastiani, D. Nuclear Velocity Perturbation Theory of Vibrational Circular Dichroism. *J. Chem. Theory Comput.* **2013**, *9*, 5305–5312.

(24) Nafie, L. A.; Freedman, T. B. Vibronic coupling theory of infrared vibrational transitions. *J. Chem. Phys.* **1983**, *78*, 7108–7116.
(25) Hirose, K.; Yongzhu, J.; Nakamura, T.; Nishioka, R.; Ueshige, T.; Tobe, Y. Preparation and evaluation of a chiral stationary phase covalently bound with chiral pseudo-18-crown-6 ether having 1-phenyl-1,2-cyclohexanediol as a chiral unit. *J. Chromatogr. A* **2005**, *1078*, 35–41.

(26) Tobe, Y.; Iketani, H.; Tsuchiya, Y.; Konishi, M.; Naemura, K. Homochiral ligands derived from cis-1-phenylcyclohexane-1,2-diol and cis-2-azido-2-phenylcyclohexanol. *Tetrahedron: Asymmetry* **1997**, *8*, 3735–3744.

(27) Naemura, K.; Ueno, K.; Takeuchi, S.; Tobe, Y.; Kaneda, T.; Sakata, Y. Azophenolic acerands having chiral 1-phenyl-cis-1,2-cyclohexanediol units: a correlation between enantiorecognitive coloration and host-guest complementarity. *J. Am. Chem. Soc.* **1993**, *115*, 8475–8476.

(28) Hartwig, B.; Suhm, M. A. Subtle hydrogen bonds: benchmarking with OH stretching fundamentals of vicinal diols in the gas phase. *Phys. Chem. Chem. Phys.* **2021**, *23*, 21623–21640.

(29) Dupont, J.; Hartwig, B.; Le Barbu-Debus, K.; Lepere, V.; Guillot, R.; Suhm, M. A.; Zehnacker, A. Homochiral vs. heterochiral preference in chiral self-recognition of cyclic diols. *Phys. Chem. Chem. Phys.* **2024**, *26*, 10610–10621.

(30) Weirich, L.; Magalhães de Oliveira, J.; Merten, C. How many solvent molecules are required to solvate chiral 1,2-diols with hydrogen bonding solvents? A VCD spectroscopic study. *Physical Chemistry Chemical Physics* **2020**, *22*, 1525–1533.

(31) Le Barbu-Debus, K.; Zehnacker, A. Competition between inter and intramolecular hydrogen bond evidenced by vibrational circular dichroism spectroscopy: The case of (1S,2R)-(-)-cis-1-amino-2-indanol. *Chirality* **2021**, *33*, 858–874.

(32) Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.;
Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A. et al. Current Status of the AMOEBA Polarizable Force Field. *J. Phys. Chem. B* **2010**, *114*, 2549–25.

(33) Rasmussen, T. D.; Ren, P.; Jay W. Ponder, F. J. Force field modeling of conformational energies: Importance of multipole moments and intramolecular polarization. *Spectroc. Acta A* **2007**, *107*, 1390–1395.

(34) Ren, P.; Wu, C.; Ponder, J. W. Polarizable Atomic Multipole-Based Molecular Mechanics for Organic Molecules. *J. Chem. Theory Comput.* **2011**, *7*, 3143–3161.

(35) Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P. Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins. *J. Chem. Theory Comput.* **2013**, *9*, 4046–4063.

(36) Thaunay, F.; Jana, C.; Clavaguéra, C.; Ohanessian, G. Strategy for Modeling the Infrared Spectra of Ion-Containing Water Drops. *J. Phys. Chem. A* **2018**, *122*, 832–842.

(37) Walker, B.; Liu, C.; Wait, E.; Ren, P. Automation of AMOEBA polarizable force field for small molecules: Poltype 2. *J. Comput. Chem.* **2022**, *43*, 1530–1542.

(38) Stone, A. J. Distributed multipole analysis, or how to describe a molecular charge distribution. *Chem. Phys. Lett.* **1981**, *83*, 233–239.

(39) Stone, A.; Alderton, M. Distributed multipole analysis. *Mol. Phys.* **1985**, *56*, 1047–1064.

(40) Jähnigen, S.; Sebastiani, D.; Vuilleumier, R. The important role of non-covalent interactions for the vibrational circular dichroism of lactic acid in aqueous solution. *Phys. Chem. Chem. Phys.* **2021**, *23*, 17232–17241.

(41) Rackers, J. A.; Wang, Z.; Lu, C.; Laury, M. L.; Lagardère, L.; Schnieders, M. J.; Piquemal, J.-P.; Ren, P.; Ponder, J. W. Tinker 8: Software Tools for Molecular Design. *J. Chem. Theory Comput.* **2018**, *14*, 5273–5289.
(42) Bowles, J.; Jähnigen, S.; Vuilleumier, R.; Calvo, F.; Clavaguéra, C.; Agostini, F. Influence of the environment on the infrared spectrum of alanine: An effective mode analysis. *J. Chem. Phys.* **2023**, *158*, 094305.

(43) Jähnigen, S. ChirPy - A Python Package for Chirality, Dynamics, and Molecular Vibrations (version 0.23.2). 2022; https://doi.org/10.5281/zenodo.4775330 (accessed October 1 2023).

(44) Jähnigen, S.; Zehnacker, A.; Vuilleumier, R. Computation of Solid-State Vibrational Circular Dichroism in the Periodic Gauge. *J. Phys. Chem. Lett.* **2021**, *12*, 7213-7220.

(45) Mohamadi, F.; Richards, N. G. J.; Guida, W. C.; Liskamp, R.; Lipton, M.; Caufield, C.; Chang, G.; Hendrickson, T.; Still, W. C. Macromodel, an integrated software system for modeling organic and bioorganic molecules using molecular mechanics. *J. Comput. Chem.* **1990**, *11*, 440-467.

(46) Śmiechowski, M. The influence of intermolecular correlations on the infrared spectrum of liquid dimethyl sulfoxide. *Spectrochimica Acta, Part A* **2021**, *260*, 119869.

(47) Scherrer, A.; Vuilleumier, R.; Sebastiani, D. Vibrational circular dichroism from ab initio molecular dynamics and nuclear velocity perturbation theory in the liquid phase. *J. Chem. Phys.* **2016**, *145*, 084101.

(48) Shen, J.; Zhu, C.; Reiling, S.; Vaz, R. A novel computational method for comparing vibrational circular dichroism spectra. *Spectroc. Acta A* **2010**, *76*, 418-422.

(49) Torabifard, H.; Reed, L.; Berry, M. T.; Hein, J. E.; Menke, E.; Cisneros, G. A. Computational and experimental characterization of a pyrrolidinium-based ionic liquid for electrolyte applications. *J. Phys. Chem. A* **2021**, *125*, 2402-2412.

(50) Ali, M. A.; Rahman, M. S.; Roy, R.; Gambill, P.; Raynie, D. E.; Halim, M. A. Structure
Elucidation of Menthol-Based Deep Eutectic Solvent using Experimental and Computational Techniques. *J. Phys. Chem. A* **2021**, *125*, 2402–2412.
TOC Graphic

<figure><img src="image_14.png" /></figure>