# Ion pairing in aqueous tetramethylammonium-acetate solutions by neutron scattering and molecular dynamics simulations

Ngoc Lan Le Nguyen, Ondrej Tichacek, Pavel Jungwirth, Hector Martinez-Seara, Philip Mason, Elise Duboué-Dijon

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

## ► To cite this version:

Ngoc Lan Le Nguyen, Ondrej Tichacek, Pavel Jungwirth, Hector Martinez-Seara, Philip Mason, et al.. Ion pairing in aqueous tetramethylammonium-acetate solutions by neutron scattering and molecular dynamics simulations. Physical Chemistry Chemical Physics, 2024, 27, pp.2553-2562. 10.26434/chemrxiv-2024-b804r. [https://doi.org/10.26434/chemrxiv-2024-b804r](https://doi.org/10.26434/chemrxiv-2024-b804r) [hal-04778485v2](https://hal.science/hal-04778485v2)

HAL Id: hal-04778485

https://hal.science/hal-04778485v2

Submitted on 7 Feb 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Cite this: DOI: 00.0000/xxxxxxxxxx

# Ion pairing in aqueous tetramethylammonium-acetate solutions by neutron scattering and molecular dynamics simulations†

Ngoc Lan Le Nguyen,[a] Ondrej Tichacek,[a] Pavel Jungwirth,[a] Hector Martinez-Seara,[a] Philip E. Mason,[*][a] Elise Duboué-Dijon,[*][b]

Received Date

Accepted Date

DOI: 00.0000/xxxxxxxxxx

Tetramethylammonium (TMA) is a ubiquitous cationic motif in biochemistry, found in the charged choline headgroup of membrane phospholipids and in tri-methylated lysine residues, which modulate histone-DNA interactions and impact epigenetic mechanisms. TMA interactions with anionic species, particularly carboxylate groups of amino acid residues and extracellular sugars, are of substantial biological relevance, as these interactions mediate a wide range of cellular processes. This study investigates the molecular interactions between TMA and acetate, representing carboxylate-containing groups, using neutron scattering experiments complemented by force field and ab initio molecular dynamics (MD) simulations. Neutron diffraction with isotopic substitution reveals specific ion pairing signatures between TMA and acetate, with simulations providing a detailed interpretation of the ion pairing structures. Force fields, notably CHARMM36 with the electronic continuum correction (ECC) (by a factor of 0.85) and AMBER99SB, capture essential pairing characteristics, but only revPBE-based ab initio MD simulations accurately model specific experimental features such as the low Q peak intensity in reciprocal space. Our study delivers a refined molecular model of TMA-carboxylate interactions, guiding the selection of force fields for complex biological systems where such interactions are of a significant importance.

## 1 Introduction

The cationic tetramethylammonium (TMA) motif is pervasive in biochemistry, particularly as the charged choline group. Choline-containing lipids are the most abundant phospholipids in biological membranes.¹,² Beyond its biological importance, choline-based ionic liquids have recently become of interest as eco-friendly and biocompatible solvents.³⁻⁵ Tetramethylammonium is also a key motif in other biological species such as betaine and lecithin, all of which play important roles in choline metabolism.⁶⁻¹⁵ The R-N(CH₃)₃⁺ motif is also found in trimethylated lysine. Lysine methylation is a translational modification that occurs on lysine residues within a wide range of proteins. Depending on the degree of methylation, the mono-, di-, and tri-methylated lysine structures are formed with different biochemical and functional consequences.¹⁶,¹⁷ Methylation patterns

modulate the interactions between histone proteins and DNA and constitute the core mechanisms of epigenetics.¹⁸⁻²³

<figure><img src="image_3.png" /><figcaption>Fig. 1 Chemical structures of choline as the head group of phosphatidylcholine and trimethyllysine, with their cationic tetramethylammonium (TMA) motif; hyaluronic acid dimer with the carboxylate groups (gray arrows) and the acetate (Ac) molecular anionic species.</figcaption></figure>

Given the ubiquity of the TMA motif—especially with choline groups decorating the surface of virtually all membranes as headgroup of phosphatidylcholine lipids²⁴—it is important to under-

<footnote>a Institute of Organic Chemistry and Biochemistry of the Czech Academy of Sciences, Flemingovo nám. 542, 160 00 Praha 6, Czech Republic, Email: philip.mason@uochb.cas.cz</footnote>

<footnote>b Université Paris Cité, CNRS, Laboratoire de Biochimie Théorique, 13 rue Pierre et Marie Curie, 75005, Paris, France, E-mail: elise.duboue-dijon@cnrs.fr</footnote>

<footnote>† Electronic Supplementary Information (ESI) available: Further details of the tested force fields. Additional neutron scattering and FFMD analysis. See DOI: 10.1039/cXCP00000x/</footnote>
stand the molecular details of its interaction with other biological species.

This study thus focuses on the interactions between TMA and anionic species, specifically the carboxylate groups present, for example, in anionic amino acid side chains (Asp and Glu) and anionic sugars. Most cell surfaces are decorated with glycosylations rich in carboxylated molecular groups, either attached to plasma proteins through terminal sialic acid groups of their N-glycosylations or present in the extracellular space surrounding plasma membranes, as in hyaluronic acid.²⁵,²⁶ Therefore, the interaction of the cationic TMA motifs of phosphatidylcholine lipids with these carboxylate groups in anionic sugars and amino acids is of direct biological relevance.

From the molecular structure point of view, TMA is an intriguing cation, with very low charge density. TMA is often denoted as the simplest "hydrophobic ion", being closer to hydrophobic species in certain respects than simple ions. For instance, choline binding sites are often rich in hydrophobic and aromatic residues rather than negatively charged groups.²⁷⁻³¹ A recurring question in the literature is then whether TMA behaves more like an ion or a hydrophobic species, e.g., in terms of a hydration structure.³²⁻³⁵ It is thus still unclear whether the geometry of the ion pair between TMA and an acetate molecular anion would resemble that of a simple atomic ion pair with the carboxylate group, or whether TMA would rather interact with the acetate leaning towards the hydrophobic CH₃ moiety.

Here, we employ a combination of neutron scattering experiments, namely neutron diffraction with isotopic substitution, with force field and *ab initio* molecular dynamics (MD) simulations to probe at the molecular level the interactions between two model ions—TMA as a model for the tri-methylated Lys or the choline headgroup—and acetate as a proxy for carboxylate groups, as shown in Figure 1. As a key result, we obtain an experimental signature of the ion pair and refine the molecular interpretation of the experimental signal with the help of force field and *ab initio* molecular dynamics simulations, providing a comprehensive molecular picture of the ion pair structure and stability.

## 2 Methods

### 2.1 Neutron scattering measurements

Two sets of neutron scattering measurements were performed on aqueous TMA-acetate (TMA-Ac) solutions, one on the 7C2 diffractometer at the Laboratoire Louis Brillouin (LLB) in Saclay and the second on the D4C diffractometer at the Institut Laue Langevin (ILL) in Grenoble (France). The former experiment acquired enough data to obtain information on the H<sub>TMA</sub>H<sub>Ac</sub> and H<sub>Ac</sub>H<sub>W</sub> correlations, where H<sub>TMA</sub>, H<sub>Ac</sub>, and H<sub>W</sub> representing hydrogen atoms on TMA, acetate, and water. The latter experiment involved fewer samples and could only be used to examine the H<sub>Ac</sub>H<sub>W</sub> correlations. Our dual experiment approach is of value for two reasons. First, neutron diffraction experiments are typically so expensive to perform that they are rarely repeated. Therefore, doubling the measurements yields a precious insight into the reproducibility of the technique over different samples run on different diffractometers in sample cans of different geometry and

neutrons of different energy. Second, the D4C results contain a small amount of lower Q data (in the reciprocal space), which turns out to be critical to assessing the relative performance of different force fields used in molecular dynamics simulations of the same systems.

All the samples measured involved 2 m TMACl (here, we use the notation 2 m as a shorthand for 2 moles of solute per 55.5 moles of either H₂O or D₂O). The isotopic compositions of the 4 samples measured at D4C were *h*₁₂TMA and *h*₃Ac in H₂O, and *h*₁₂TMA and *d*₃Ac in H₂O with the same samples repeated in D₂O. The 8 solutions measured at 7C2 involved the same solutions as measured at D4C, with an additional set of samples containing *d*₁₂TMA.

The general method used to prepare all the samples was identical to ensure the greatest isotopic parity between the solutions. Two TMACl solutions (one *d*₁₂TMACl, one *h*₁₂TMACl) were run down a Dowex 1X8 chloride ion exchange column (all purchased from Sigma Aldrich) and an aliquot of each of the resulting solutions was calibrated via pH titration with a common HCl solution. Four other equimolar aliquots were taken from the two solutions and were titrated versus identical concentration solutions of *h*₃HOAc and *d*₃HOAc to yield four solutions containing equimolar concentrations of TMA-Ac (*d*₁₂TMA-*d*₃Ac, *d*₁₂TMA-*h*₃Ac, *h*₁₂TMA-*d*₃Ac, and *h*₁₂TMA-*h*₃Ac). The water was first evaporated and then added to give four chemically identical solutions at a 2 m concentration in H₂O. An identical volume of each sample was prepared (~ 2.6 mL). After the acquisition of the neutron scattering data from the solution, H₂O was then removed at 1 atm by a custom-made closed-loop apparatus via evaporation. Typically >90% of H₂O was removed before ~ 6 mL of D₂O was added and the cycle repeated. This cycle was repeated 3 more times before the correct amount of D₂O was added gravimetrically to reach the desired concentration of the D₂O solutions.

For the 7C2 experiments neutron scattering patterns were acquired for all 8 solutions in a thin-walled vanadium cell (inner diameter of 5.94 mm and wall thickness of 0.14 mm). The neutron wavelength was 0.7203 Å. For the D4C experiments, neutron scattering patterns were acquired on all four solutions in a null scattering Titanium/Zirconium cell (inner diameter of 5 mm and wall thickness of 0.75 mm). The neutron wavelength was 0.5003 Å. The total neutron scattering patterns of these solutions were recorded for about 2 hours for the D₂O solutions and 4 hours for the H₂O solutions. The data was corrected for multiple scattering and absorption and normalized versus a standard vanadium rod.

The total scattering pattern of these solutions (Figures 2a-b) is a complex sum of structural correlations. When the neutron diffraction patterns of two isotopically different but chemically identical solutions (assumed to be structurally isomorphic) are subtracted from each other, all of the correlations cancel out except those of the substituted nuclei. This is the core of the neutron diffraction technique with isotopic substitution (NDIS).³⁶⁻³⁹

For instance, subtracting the neutron scattering patterns of two D₂O solutions, one containing D<sub>TMA</sub>D<sub>Ac</sub> and the other H<sub>TMA</sub>D<sub>Ac</sub> yields the so-called the first order difference ΔS<sub>H<sub>TMA</sub></sub><sup>X<sub>D<sub>Ac</sub></sub></sup>(Q) that can
<figure><img src="image_4.png" /><figcaption>Fig. 2 Total diffraction patterns for the solutions of the isotopic compositions of the 8 samples measured at 7C2 were (a) $h_{12}$TMA–$h_3$Ac ($H_{TMA}H_{Ac}$, green), $h_{12}$TMA–$d_3$Ac ($H_{TMA}D_{Ac}$, magenta), $d_{12}$TMA–$h_3$Ac ($D_{TMA}H_{Ac}$, violet), and $d_{12}$TMA–$d_3$Ac ($D_{TMA}D_{Ac}$, orange) in $H_2O$. (b) with the same samples repeated in $D_2O$. (c) First order difference $\Delta S_{H_{TMA}}^{X_{D_{Ac}}}(Q)$ (light blue) obtained by subtracting the neutron scattering patterns of two $H_2O$ solutions, one containing $D_{TMA}D_{Ac}$ and the other $H_{TMA}D_{Ac}$. Subtracting the signals from containing $H_{Ac}$ rather than $D_{Ac}$ yields the first order difference $\Delta S_{H_{TMA}}^{X_{H_{Ac}}}(Q)$ (red). (d) First order differences of $\Delta S_{H_{TMA}}^{X_{D_{Ac}}}(Q)$ (light blue) and $\Delta S_{H_{TMA}}^{X_{H_{Ac}}}(Q)$ (red) in $D_2O$ solutions. (e) Second order difference $6.8 \cdot (S_{H_{TMA}H_{Ac}}(Q) - 1)$, obtained through the difference of the two first order differences shown in (c) (red) and in (d) (blue). (f) Second order differences after the background subtraction and terminating data up to $6 \AA^{-1}$.</figcaption></figure>

be expressed as a function of only 6 structural correlations (Equation 1). Similarly, subtracting the signals from solutions containing $H_{Ac}$ rather than $D_{Ac}$ yields the first order difference $\Delta S_{H_{non}}^{X_{H_{Ac}}}(Q)$ (see Equation 2). These two first order differences (See Figures 2c-d) report on the correlations between non-exchangeable H on TMA and every other atom (X) in the system in units of mbarns (mb). Using standard literature methods,⁴⁰ the prefactors in the equations were computed from the atomic concentration and neutron scattering lengths of the various elements in the system. The differences between Equations 1 and 2 are highlighted in bold.

$$ 
\begin{aligned}
\Delta S_{H_{TMA}}^{X_{D_{Ac}}}(Q) &= S[d_{12}\text{TMA--d}_3\text{Ac}](Q) - S[h_{12}\text{TMA--d}_3\text{Ac}](Q) \\
&= 80.4 \cdot S_{H_{TMA}D_W}(Q) + 37.5 \cdot S_{H_{TMA}O}(Q) \\
&\quad + 8.6 \cdot S_{H_{TMA}C}(Q) + 2.0 \cdot S_{H_{TMA}N}(Q) \\
&\quad + 4.4 \cdot S_{H_{TMA}H_{Ac}}(Q) + 3.8 \cdot S_{H_{TMA}H_{TMA}}(Q) - 136.8
\end{aligned}
\quad (1) $$

$$ 
\begin{aligned}
\Delta S_{H_{TMA}}^{X_{H_{Ac}}}(Q) &= S[d_{12}\text{TMA--h}_3\text{Ac}](Q) - S[h_{12}\text{TMA--h}_3\text{Ac}](Q) \\
&= 80.4 \cdot S_{H_{TMA}D_W}(Q) + 37.5 \cdot S_{H_{TMA}O}(Q) \\
&\quad + 8.6 \cdot S_{H_{TMA}C}(Q) + 2.0 \cdot S_{H_{TMA}N}(Q) \\
&\quad - 2.4 \cdot S_{H_{TMA}H_{Ac}}(Q) + 3.8 \cdot S_{H_{TMA}H_{TMA}}(Q) - 130.0
\end{aligned}
\quad (2) $$

The data can be further simplified by performing a second order difference $\Delta\Delta S_{H_{TMA}H_{Ac}}(Q) = \Delta S_{H_{TMA}}^{X_{D_{Ac}}}(Q) - \Delta S_{H_{TMA}}^{X_{H_{Ac}}}(Q)$ (see Equation 3). This second order difference then reports on a single correlation between the first nuclei that was substituted (H/D on the TMA) and the second nuclei (H/D on acetate), see Figure 2e.

$$ 
\begin{aligned}
\Delta\Delta S_{H_{TMA}H_{Ac}}(Q) &= S[d_{12}\text{TMA--d}_3\text{Ac}](Q) - S[h_{12}\text{TMA--d}_3\text{Ac}](Q) \\
&\quad - (S[d_{12}\text{TMA--h}_3\text{Ac}](Q) - S[h_{12}\text{TMA--h}_3\text{Ac}](Q)) \\
&= 6.8 \cdot (S_{H_{TMA}H_{Ac}}(Q) - 1)
\end{aligned}
\quad (3) $$
Following a similar strategy, but subtracting the signals obtained in H₂O and D₂O solutions, we also extracted the HₐcHw pair correlation, ΔΔS<sub>HₐcHw</sub>(Q) (see details in Equations S1-3 in the Supporting Information (SI)), which reports on acetate hydration and, indirectly, on ion pairing.

The difficulty of these experiments scales approximately with
the size of the second order difference versus the total scattering
of the solution (~400mb). While the correlation between H<sub>Ac</sub>
and H<sub>w</sub> is relatively straightforward to measure (~31mb, see
Equation S3 in SI), the H<sub>TMA</sub> to H<sub>Ac</sub> contrast is small (~7mb).
Nevertheless, we can obtain the H<sub>TMA</sub>H<sub>Ac</sub> correlation from both
H<sub>2</sub>O and D<sub>2</sub>O samples, which provides us with an internal consis-
tency check as the same function is measured twice.

Due to strong inelastic scattering from ¹H and the Placzek ef-
fect,³⁷ hydrogen-rich samples exhibit significant background, es-
pecially in light water H₂O, see Figure 2e. This effect is sig-
nificantly reduced in heavy water D₂O. The inelastic scattering
scales with ¹H atomic concentration, so the first-order differences
should display the same Placzek background, which disappears
in the second-order difference. After the background subtraction
(See Figure 2f) and terminating the data using a window function
up to 6 Å⁻¹ for the two second order differences, the constant
value in the last figure confirms the absence of background and
indicates consistent chemical composition across the solutions.
Data truncation at 6 Å⁻¹ aims to reduce ringing in the Fourier
transform due to experimental noise. Its impact was evaluated
by truncating the MD data for CHARMM36-ECC85 at the same
Q-value, before transforming it back into R-space, which had no
significant effect on the Fourier transform (Figure S2).

## 2.2 Force field molecular dynamics (FFMD) simulations

All force field-based MD simulations were prepared and carried
out using the GROMACS 2021.4 software.⁴¹ The system of inter-
est is a 2 m TMA–Ac aqueous solution, that comprises 50 TMA
molecules, 50 acetate molecules, and 1388 water molecules in a
~3.7 nm cubic unit cell with periodic boundary conditions.

In order to identify the force fields giving the description of
ion pairing in best agreement with experiment, we tested sev-
eral combinations of force fields for TMA, acetate, and water
molecules. For the molecular ions, we tested both the broadly
employed CHARMM36⁴² and AMBER99SB⁴³,⁴⁴ force fields, as
well as scaled-charge variants using the Electronic Continuum
Correction (ECC) approach⁴⁵⁻⁴⁸ with different charge-scaling
factors (0.75 and 0.85), which have been suggested in litera-
ture to improve the description of ion pairing in solution. As
none of these force fields fully captured the experimental sig-
nal, we also tested variants of CHARMM36 with modified CO-
Ac or CH-TMA dipole moments in the hope to favor different
ion pairing geometries. Simulations were performed with the
corresponding TIP3P water model.⁴⁹,⁵⁰ The main text reports
results obtained with the force fields CHARMM36—standard
(CHARMM36-FULL) and ECC with charge scaling factors of
0.85 (CHARMM36-ECC85) and 0.75 (CHARMM36-ECC75)—and
AMBER99SB—standard (AMBER99SB-FULL) and ECC with the
charge scaling factor of 0.85 (AMBER99SB-ECC85). The Am-

ber force field with 0.75 (AMBER99SB-ECC85) scaling and the
dipole changes (CHARMM36-Low CH dipole, CHARMM36-Low
CO dipole, and CHARMM36-High CO dipole) are presented in
the SI. The list of all tested combinations and the details of all
force fields are provided in Tables S1-3.

The system was prepared with GROMACS tools before mini-
mization and equilibrations in the NPT ensemble (0.5 ns). The
simulations were propagated over 200 ns with a 2 fs time step in
the constant temperature (310 K) and constant pressure (1 bar)
ensemble (NPT), employing the Nosé-Hoover thermostat⁵¹⁻⁵³
(τ<sub>T</sub> = 1 ps) and Parrinello-Rahman barostat⁵⁴ (τ<sub>P</sub> = 5 ps), respec-
tively. The trajectory was saved every 0.2 ps for subsequent analy-
ses. All the bonds containing hydrogen were constrained with the
LINCS algorithm⁵⁵ and the Particle Mesh Ewald (PME) method⁵⁶
was used to calculate the long-range electrostatic interactions. A
1.2 nm cut-off distance was used for short-range electrostatic and
van der Waals interactions. Both potentials were shifted to zero
at the cutoff distance. For van der Waals interactions, the shift
starts at 1.0 nm.

## 2.3 Ab initio molecular dynamics (AIMD) simulations

To further characterize the TMA–Ac ion pairing structure, we
complemented our experiments and force field-based simulations
with Born–Oppenheimer *ab initio* molecular dynamics simula-
tions (AIMD). Due to the high computational cost of such sim-
ulations, the simulated system then included a single ion pair
with a TMA cation and an acetate anion, solvated by 128 water
molecules with periodic boundary conditions.

Force field molecular dynamics simulations of the same system
were employed to pre-equilibrate the system and generate initial
configurations for the AIMD simulations in different ion pair ge-
ometries. After minimization, equilibration in the NPT ensemble
was done for 0.5 ns to determine the average size of the simu-
lation cell, followed by a 20 ns simulation in the NPT ensemble.
We extracted from there 20 representative initial configurations
of the ion pair.

Starting from these pre-equilibrated structures, AIMD simula-
tions were conducted with the CP2K program package⁵⁷ version
9.1, using the generalized gradient approximation revPBE DFT
functional⁵⁸⁻⁶⁰ incorporating the D3 dispersion correction.⁶¹⁻⁶³
The potential energy of core electrons was substituted with GTH
pseudopotentials,⁶⁴,⁶⁵ while valence electrons were described us-
ing the triple-ζ basis set TZV2P with polarization functions.⁶⁶
The auxiliary plane wave basis set in the Gaussian Plane Wave
Method (GPW) method⁶⁷ was used with a cutoff of 400 Ry. The
equilibration of all systems was performed using AIMD Langevin
dynamics for a minimum of 10 ps with a damping constant of
γ = 0.02 ps⁻¹,⁶⁸ with the temperature maintained at 310 K us-
ing a global CVSR thermostat⁶⁹ with a time constant of 1 ps.
A cubic simulation cell with dimensions of 16.8 Å was used.
To enhance sampling, twenty 100 ps parallel AIMD simulations
were conducted in the NVT ensemble, resulting in an overall
AIMD simulation time of 2000 ps used for structural analyses.
All simulation data can be found at https://zenodo.org (DOI:
https://doi.org/10.5281/zenodo.14050921).
## 2.4 MD simulations analyses

Standard analyses, e.g., calculation of radial distribution functions or analysis of contacts, were all performed using GROMACS tools. Contacts between TMA and acetate were defined based on the distance between the nitrogen atom of the TMA cation (N<sub>TMA</sub>) and the carboxyl carbon of acetate (C<sub>Ac</sub>). Based on the radial distribution function of the N<sub>TMA</sub> and C<sub>Ac</sub> atoms, we determined the typical N<sub>TMA</sub>-C<sub>Ac</sub> distances used to distinguish direct contact and solvent-shared ion pairs. Accordingly, the number of contacts and then its distribution for each force field simulation were calculated. Moreover, to quantify the orientations of acetate relative to TMA ion, we computed the distribution of angle θ, which is the angle between TMA nitrogen, acetate carboxyl carbon, and acetate methyl carbon atoms.

# 3 Results and Discussion

## 3.1 TMA-acetate ion pair formation from NDIS experiments

The structure of TMA-Ac 2 m in aqueous solutions was probed using neutron scattering at the 7C2 diffractometer, employing the technique of neutron diffraction with isotopic substitution (NDIS). The double difference signal ΔΔS<sub>H<sub>TMA</sub>H<sub>Ac</sub></sub>(Q), known as the second order difference (see Equation 3 and Figure 3a), reports on the single correlation between TMA and acetate non-exchangeable hydrogens, which makes it a very sensitive probe of the formation of TMA-Ac ion pairs and their structure. The Fourier-transformation of the experimental signal to real-space ΔΔG<sub>H<sub>TMA</sub>H<sub>Ac</sub></sub>(r) is the single radial distribution function between H<sub>TMA</sub> and H<sub>Ac</sub>:

$$ \Delta\Delta G_{H_{TMA}H_{Ac}}(r) = 6.8 \cdot (g_{H_{TMA}H_{Ac}}(r) - 1) \quad (4) $$

<figure><img src="image_5.png" /><figcaption>Fig. 3 Experimental second order difference with background subtraction (at 7C2) (a) ΔΔS<sub>H<sub>TMA</sub>H<sub>Ac</sub></sub>(Q) and (b) ΔΔG<sub>H<sub>TMA</sub>H<sub>Ac</sub></sub>(r) of H<sub>2</sub>O (red) and D<sub>2</sub>O (blue). These two functions should be identical and represent an internal consistency check for the method and data.</figcaption></figure>

This correlation was obtained twice from completely independent samples either in D<sub>2</sub>O or in H<sub>2</sub>O (see Figure 3), which provides a stringent reproducibility test and allows to assess the robustness of individual signal features. As expected from the high incoherent scattering from the H<sub>2</sub>O sample, the noise is higher on the H<sub>2</sub>O than D<sub>2</sub>O signal. Comparison of the two signals in Figure 3a shows that while the location of the initial raise in the reciprocal space (Q-space) signal and the height and location of the first peak are very robust, subsequent oscillations are not. In real space, the corresponding robust features are the position, height, and width of the peak around 6 Å. Features at low distances below 2 Å are mostly noise and Fourier-transform artifacts.

## 3.2 Insights from molecular simulations

The neutron scattering experiments described above provide a simple and robust signal that reports exclusively on a single correlation between the two ions, and thus probes very sensitively ion pairing in the solution. However, molecular interpretation of this signal is not straightforward. To gain additional insights, we complemented these experiments with molecular simulations, in line with several prior works where this combination provided valuable insights into the structure of electrolytes and ion pairing in solution.<sup>71–80</sup>

Molecular dynamics simulations of TMA-Ac at the experimental concentration have been performed using various force fields, based either on the CHARMM36 or AMBER99SB force fields. In addition to the original full charge version of these force fields, we derived variants following the Electronic Continuum Correction (ECC) strategy of using scaled charges to effectively take into account electronic polarization (see Methods).<sup>45–48</sup> This approach was previously demonstrated to significantly improve the description of ion pairing in electrolyte solutions.<sup>71–80</sup> Two different scaling factors (0.75 or 0.85) were tested, as previous work suggested that the original 0.75 scaling factor may in some cases underestimate the strength of ion pairing.<sup>81</sup>

For each force field we evaluated the radial distribution function ΔΔG<sub>H<sub>TMA</sub>H<sub>Ac</sub></sub>(r). A Fourier transform then provides ΔΔS<sub>H<sub>TMA</sub>H<sub>Ac</sub></sub>(Q), which can be directly compared with the experimental double difference signal (Figures 4a-d). Comparison between experimental and computed signal is more robust in Q-space than in real space (R-space) since the limited experimental Q-range causes artifacts in the Fourier transformed signal. Overall, among the used force fields, the scaled charge CHARMM36 with a 0.85 scaling factor (CHARMM36-ECC85) and the full charge AMBER99SB force fields provide the best fit of the experimental data, especially in the low Q region with the location of the initial sharp rise of the signal, which is a very robust experimental feature. The computed density maps of H<sub>TMA</sub> atoms around acetate (Figure 4) are very similar for all force fields employed and show that in these simulations TMA interacts with acetate from above and below the -COO<sup>-</sup> plane, rather than along the direction of the carboxylate group. However, none of the present force fields could properly capture the height of the first peak at Q = 1.2 Å<sup>-1</sup> (Figures 4a-b), nor, in real space, the features below 5 Å. Neither did additional tests of force field variants with
<figure><img src="image_6.png" /></figure>

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

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

<figure><img src="image_9.png" /><figcaption>(e)</figcaption></figure>

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

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

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

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

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

CHARMM36 FULL

CHARMM36 ECC85

CHARMM36 ECC75

AMBER99SB FULL

AMBER99SB ECC85

<caption>Fig. 4 Comparison of the second order difference in Q-space $\Delta\Delta S_{H_{TMA}H_{Ac}}(Q)$ of (a) CHARMM36-FULL (green), CHARMM36-ECC85 (blue), CHARMM36-ECC75 (magenta), and experiment (black); (b) AMBER99SB-FULL (red), AMBER99SB-ECC85 (orange), and experiment (black). Comparison of the second order difference in real space (R-space) $\Delta\Delta G_{H_{TMA}H_{Ac}}(r)$ of (c) CHARMM36-FULL (green), CHARMM36-ECC85 (blue), CHARMM36-ECC75 (magenta), and experiment (black); (d) AMBER99SB-FULL (red), AMBER99SB-ECC85 (orange), and experiment (black). The colors used to represent the different force field sectors are consistent throughout this article. The experiment was performed at 7C2. (e) The density maps, prepared with VMD,⁷⁰ show the density of $H_{TMA}$ around acetate corresponding to the used force fields.</caption>

decreased or increased dipole moments on both CH-TMA and CO-Ac bonds (see Figure S3).

While none of the force fields properly capture all the experimental features they, however, recover a range of different ion pairing behaviors with varying ion pair structures. For instance, simulations performed with CHARMM36 (both the standard version and scaled charges) present a higher proportion of solvent-shared ion pairs and less direct contacts between the two ions, while AMBER99SB-ECC85 is the model with the highest peak in $g(r)$ corresponding to direct contact ion pairs (see Figure 5a). The two models that best fit the experimental data (CHARMM36-ECC85 and AMBER99SB-FULL) exhibit a very similar propensity to form contact ion pairs. The internal structures of the ion pairs are also very similar to each other, with the dominant interaction

mode being the acetate lying flat against the faces of TMA. Nevertheless, we see a broad distribution of relative orientations of the two ions from the carboxylic group pointing towards the TMA to an interaction through the acetate hydrophobic tail (see Figure 5b).

To understand how the real ion pairing in solution correlates with the picture obtained from MD simulations, which never fully capture the experimental signal, we can use another combination of neutron signals, that offers complementary information about ion pairing in the solution. Namely, while we so far focused on the $\Delta\Delta S_{H_{TMA}H_{Ac}}(Q)$ which directly probes the ion-ion correlation, the second order difference $\Delta\Delta S_{H_{Ac}H_{W}}(Q)$—obtained from a different combination of the raw Q-space signals in $D_{2}O$ and $H_{2}O$ (see Methods)—also reports on ion pairing, though more indirectly,
see Equation 3.

$\Delta\Delta S_{H_A c_{H_W}}(Q)$ reports again on a single pair-correlation in the solution, this time between the acetate non-exchangeable H and water exchangeable H. This correlation is not directly relevant for ion pairing, but since the formation of TMA–Ac ion pairs leads to water depletion in the vicinity of acetate, the formation of ion pairs impacts the $H_{Ac}H_W$ pair correlation. More precisely, the amount of ion pairing is reflected in the experimental $\Delta\Delta S_{H_A c_{H_W}}(Q)$ mostly in the height of the low-Q plateau (Figure 6). This plateau and, overall, the $\Delta\Delta S_{H_A c_{H_W}}(Q)$, is very well described by the two force fields that so far gave the best description of the TMA–Ac solution (i.e., AMBER99SB-FULL and CHARMM36-ECC85). This suggests that these two force fields properly capture the amount of TMA–Ac ion pairing in the solution, and that other tested force fields over- or under-estimate it.

<figure><img src="image_15.png" /><figcaption>Fig. 5 (a) Radial distribution functions from the TMA nitrogen atom to the acetate carboxyl carbon atom, $g(r_{N_{TMA}-C_{Ac}})$ from the selected force fields and (b) Distribution of angle $\theta$ quantifying the orientation of the acetate relative to the TMA ion. The snapshots illustrating the different orientations were prepared with VMD.⁷⁰</figcaption></figure>

What can explain, then, the distinguished feature in the $H_{TMA}H_{Ac}$ signal that could not be captured by our force fields? If the overall amount of association between the ions is correctly described in our simulations, then this would suggest that it is the internal geometry of the ion pair (relative orientation of the two ions) that is not fully captured by simulations. However, characterizing how the orientation of the ions should change to give a better fit to experimental data is very challenging.

To investigate in more detail the structure of contact ion pairs,

we thus performed multiple *ab initio* MD simulations, at the revPBE DFT level of theory. The large computational cost of such simulations limits the size of the simulated system to a single TMA–Ac ion pair and the length of the simulations to 100 ps. However, starting from 20 different geometries of the ion pairs, these simulations offer enough sampling to yield the $H_{TMA}H_{Ac}$ pair correlation function, and obtain, after Fourier-transform, the double difference in Q-space, $\Delta\Delta S_{H_{TMA}H_{Ac}}(Q)$ (Figure 7a).

<figure><img src="image_16.png" /><figcaption>Fig. 6 Second order difference $\Delta\Delta S_{H_A c_{H\_water}}(Q)$ of experiments at 7C2 (black) and D4C (gray) versus FFMD simulations with different force fields: (a) CHARMM36-FULL (green), CHARMM36-ECC85 (blue), and CHARMM36-ECC75 (magenta); (b) AMBER99SB-FULL (red) and AMBER99SB-ECC85 (orange).</figcaption></figure>

In contrast with all the empirical force fields we tested, the signal computed from our AIMD simulation provides an excellent fit of the experimental signal. Remarkably, it captures perfectly the height of the peak at $Q = 1.2 \AA^{-1}$, which remained out-of-reach of empirical force fields. This success suggests that the internal geometry of the ion pair—i.e., the relative orientation of the two ions—is properly described by the AIMD simulations. If we now compare (Figure 7b) the distribution of angle $\theta$ that quantifies the orientation of the acetate relative to the TMA ion, it seems distinct in AIMD simulations (despite the higher numerical noise due to limited sampling) from all the empirical force field simulations, with a preference for larger angles. This means that the carboxylate group, instead of lying flat against TMA as in force field simulations, would have a larger tendency to point towards the TMA cation. In other words, in AIMD, TMA behaves less "hy-
drophobic" in ion pairing than what is depicted in force field sim-
ulations.

<figure><img src="image_17.png" /><figcaption>Fig. 7 (a) Comparison of the second order difference $\Delta\Delta S_{H_{TMA}H_{Ac}}(Q)$ of AIMD simulations (pink), FFMD CHARMM36-ECC85 (blue) and experiment at 7C2 (black). (b) Distribution of angle $\theta$ quantifying the orientation of the acetate relative to the TMA ion in AIMD simulations versus the experiment.</figcaption></figure>

4 Conclusions

Neutron diffraction experiments with double isotopic substitution allowed us to characterize ion pairing in TMA–Ac solutions. The obtained signal is robust and reflects a single correlation between $H_{TMA}$ and $H_{Ac}$. Comparison of the experimental signal with that computed from force field and *ab initio* molecular dynamics simulations enabled us to obtain a molecular interpretation of the experimental signal, accurately characterizing the amount of ion pairing and the ion pair molecular structure.

Among the tested force fields, CHARMM36 with ECC charge scaling with a 0.85 scaling factor and the standard AMBER99SB full charges perform the best, in line with the fact that charge scaling tends not to be a decisive factor for the pairing of a low charge density ions like TMA. The two force fields give a very similar picture of weak ion pairing with the acetate lying flat against the TMA surface.

However, none of the investigated force fields—despite trying
several force field variants with modified charges and CO or CH
dipole moments—captures all the features of the experimental
signal, in particular, the intensity of the peak at low Q. In contrast,
AIMD simulations based on the revPBE DFT functional provide an
excellent fit to experiment, properly capturing also this feature.
Analysis of the AIMD ion pair structure shows that the acetate

carboxylate group points slightly more towards the TMA cation
in the ion pair compared to the force field-based simulations.

The combination of experiments with various flavors of molec-
ular simulations thus allowed us to reach an in depth molecular
understanding of the amount of ion pairing and on the internal
geometry of the ion pair formed between two biologically rel-
evant moieties, TMA and acetate. This study thus helps refine
our molecular understanding of biological processes that involve
these interactions. Additionally, the present study provides guide-
lines for the choice of adequate force fields for simulations of
complex biological systems where such interactions play impor-
tant roles.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was performed using neutron beam time allocated on
the 7C2 diffractometer at Laboratoire Louis Brillouin (LLB, Saclay,
France; proposal 602) and on the D4C beamline at the Institut
Laue Langevin (ILL, Grenoble, France; proposal 8-03-807). We
especially thank Henry Fischer for helping run the experiments
on D4C. N.L.L.N and O.T. acknowledge the Faculty of Mathemat-
ics and Physics of the Charles University (Prague, Czech Repub-
lic) where they are enrolled as PhD students, and the Interna-
tional Max Planck Research School (IMPRS) for "Quantum Dy-
namics and Control" and "Many-Particle Systems in Structured
Environments" (Dresden, Germany) for support. E.D.D. acknowl-
edges support from the "Initiative d'Excellence" program from the
French State (Grant "DYNAMO", ANR11-LABX-0011). P.J. ac-
knowledges support from the European Research Council via an
ERC Advanced Grant no. 101095957.

References

1 S. Tayebati, G. Marucci, C. Santinelli, M. Buccioni and F. Amenta, *Curr. Med. Chem.*, 2015, **22**, 4328–4340.

2 P. L. Yeagle, in *The Membranes of Cells (Third Edition), Chapter 12 - Lipid-Protein Interactions in Membranes*, ed. P. L. Yeagle, Academic Press, Boston, 2016, pp. 291–334.

3 L. Gontrani, *Biophys. Rev.*, 2018, **10**, 873–880.

4 B. L. Gadilohar and G. S. Shankarling, *J. Mol. Liq.*, 2017, **227**, 234–261.

5 X. Li, N. Ma, L. Zhang, G. Ling and P. Zhang, *Int J Pharm.*, 2021, **612**, 121366.

6 D. Dobrijević, K. Pastor, N. Nastić, F. Özogul, J. Krulj, B. Kokić, E. Bartkiene, J. M. Rocha and J. Kojić, *Mol.*, 2023, **28**, 4824.

7 M. K. Arumugam, M. C. Paal, T. M. Donohue, M. Ganesan, N. A. Osna and K. K. Kharbanda, *Biology*, 2021, **10**, 456.

8 J. Du, L. Shen, Z. Tan, P. Zhang, X. Zhao, Y. Xu, M. Gan, Q. Yang, J. Ma, A. Jiang, G. Tang, Y. Jiang, L. Jin, M. Li, L. Bai, X. Li, J. Wang, S. Zhang, and L. Zhu, *Nutrients*, 2018, **10**, 131–147.

9 P. M. Ueland, P. I. Holm and S. Hustad, *Clin. Chem. Lab. Med.*, 2005, **43**, 1069–1075.

10 P. M. Ueland, *J. Inherit. Metab. Dis.*, 2011, **34**, 3–15.
11 C. Wang, C. Ma, L. Gong, S. Dai and Y. Li, *Eur. J. Pharmacol.*, 2021, **912**, 174604.

12 C. Robert, L. Couëdelo, C. Vaysse and M.-C. Michalski, *Biochim.*, 2020, **169**, 121-132.

13 M.-J. LeBlanc, S. Brunet, G. Bouchard, T. Lamireau, I. M. Yousef, V. Gavino, E. Lévy and B. Tuchweber, *J. Nutr. Biochem.*, 2003, **14**, 40-48.

14 D. J. Canty and S. H. Zeisel, *Nutr. Rev.*, 1994, **52**, 327-339.

15 S. M. Ezzat, M. A. Salem, N. M. E. Mahdy and M. M. Mahfouz, in *Antioxidants Effects in Health, Chapter 4.9 - Lecithin*, ed. S. M. Nabavi and A. S. Silva, Elsevier, 2022, pp. 375-386.

16 M. Luo, *Chem. Rev.*, 2018, **118**, 6656-6705.

17 A. Kupai, R. M. Vaughan, B. M. Dickson and S. B. Rothbart, *Front. Cell Dev. Biol.*, 2020, **8**, 1-13.

18 A. Akter, S. Takahashi, W. Deng, D. J. Shea, E. Itabashi, M. Shimizu, N. Miyaji, K. Osabe, N. Nishida, Y. Suzuki, C. A. Helliwell, M. Seki, W. J. Peacock, E. S. Dennis and R. Fujimoto, *DNA Res.*, 2019, **26**, 433-443.

19 M. Foroozani, M. P. Vandal and A. P. Smith, *Planta*, 2021, **253**, 1-17.

20 F. Turck, F. Roudier, S. Farrona, M.-L. Martin-Magniette, E. Guillaume, N. Buisine, S. Gagnot, R. A. Martienssen, G. Coupland and V. Colot, *PLoS Genet.*, 2007, **3**, 855-866.

21 B. E. Collins, C. B. Greer, B. C. Coleman and J. D. Sweatt, *Epigenetics Chromatin*, 2019, **12**, 433-443.

22 N. R. Rose and R. J. Klose, *Biochim. Biophys. Acta. Gene. Regul. Mech.*, 2014, **1839**, 1362-1372.

23 A. Agredo and A. L. Kasinski, *Front. genet.*, 2023, **14**, 1243395.

24 I. Pasquali-Ronchetti, D. Quaglino, G. Mori, B. Bacchelli and P. Ghosh, *J. Struct. Biol.*, 1997, **120**, 1-10.

25 M. Suniti, H. V. C., M. R. R. and G. Shibnath, *Front. immunol.*, 2015, **6**, 201.

26 S. Garantziotis and R. C. Savani, *Matrix Biol.*, 2019, **78-79**, 1-10.

27 Z. Zhang, W. Li, C. Frolet, R. Bao, A.-M. di Guilmi, T. Vermet and Y. Chen, *Acta Crystallogr. Sect. F Struct. Biol. Cryst. Commun.*, 2009, **65**, 757-761.

28 B. Maestro, I. Velasco, I. Castillejo, M. Arévalo-Rodríguez, Ángel Cebolla and J. M. Sanz, *J. Chromatogr. A*, 2008, **1208**, 189-196.

29 J. Gutiérrez-Fernández, M. Saleh, M. Alcorlo, A. Gómez-Mejía, D. Pantoja-Uceda, M. A. Treviño, F. Voß, M. R. Abdullah, S. Galán-Bartual, J. Seinen, P. A. Sánchez-Murcia, F. Gago, M. Bruix, S. Hammerschmidt and J. A. Hermoso, *Sci. Rep.*, 2016, **6**, 38094.

30 C. Yuan and C. Kent, *J. Biol. Chem.*, 2004, **279**, 17801-17809.

31 B. Maestro, H. Zamora-Carreras, M. Ángeles Jiménez and J. M. Sanz, *Int. J. Biol. Macromol.*, 2021, **190**, 679-692.

32 P. E. Mason, T. Martinek, B. Fábián, M. Vazdar, P. Jungwirth, O. Tichacek, E. Duboué-Dijon and H. Martinez-Seara, *Phys. Chem. Chem. Phys.*, 2024, **26**, 3208-3218.

33 J. Turner, A. Soper and J. Finney, *Mol. Phys.*, 1990, **70**, 679-700.

34 N. G. Polydorou, J. D. Wicks and J. Z. Turner, *J. Chem. Phys.*, 1997, **107**, 197-204.

35 E. W. Lang, S. Bradl, W. Fink, H. Radkowitsch and D. Girlich, *J. Phys.: Condens. Matter*, 1990, **2**, SA195-SA200.

36 P. H. K. de Jong and G. W. Neilson, *J. Phys. Condens. Matter*, 1996, **8**, 9275.

37 P. E. Mason, G. W. Neilson, C. E. Dempsey and J. W. Brady, *J. Am. Chem. Soc.*, 2006, **128**, 15136-15144.

38 P. E. Mason, S. Ansell and G. W. Neilson, *J. Phys. Condens. Matter*, 2006, **18**, 8437.

39 G. W. Neilson, P. E. Mason, S. Ramos and D. Sullivan, *Phil. Trans. R. Soc. A.*, 2001, **359**, 1575-1591.

40 J. E. Enderby, D. E. Williams and J. Randall, *Proc. R. Soc. Lond. A Math. Phys. Sci.*, 1975, **345**, 107-117.

41 M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, *SoftwareX*, 2015, **1-2**, 19-25.

42 J. B. Klauda, R. M. Venable, J. A. Freites, J. W. O'Connor, D. J. Tobias, C. Mondragon-Ramirez, I. Vorobyov, A. D. J. MacKerell and R. W. Pastor, *J. Phys. Chem. B*, 2010, **114**, 7830-7843.

43 J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, *J. Comput. Chem.*, 2004, **25**, 1157-1174.

44 V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg and C. Simmerling, *Proteins*, 2006, **65**, 712-725.

45 I. V. Leontyev and A. A. Stuchebrukhov, *J. Chem. Phys.*, 2009, **130**, 085102.

46 B. J. Kirby and P. Jungwirth, *J. Phys. Chem. Lett.*, 2019, **10**, 7531-7536.

47 E. Duboué-Dijon, M. Javanainen, P. Delcroix, P. Jungwirth and H. Martinez-Seara, *J. Chem. Phys.*, 2020, **153**, 050901.

48 R. Nencini, C. Tempra, D. Biriukov, M. Riopedre-Fernandez, V. Cruces Chamorro, J. Polák, P. E. Mason, D. Ondo, J. Heyda, O. H. S. Ollila, P. Jungwirth, M. Javanainen and H. Martinez-Seara, *J. Chem. Theory Comput.*, 2024, **20**, 7546-7559.

49 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, *J. Chem. Phys.*, 1983, **79**, 926-935.

50 A. D. J. MacKerell, D. Bashford, M. Bellott, R. L. J. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin and M. Karplus, *J. Phys. Chem. B*, 1998, **102**, 3586-3616.

51 S. Nosé, *J. Chem. Phys.*, 1984, **81**, 511-519.

52 S. Nosé, *Mol. Phys.*, 1984, **52**, 255-268.

53 W. G. Hoover, *Phys. Rev. A*, 1985, **31**, 1695-1697.

54 M. Parrinello and A. Rahman, *J. Appl. Phys.*, 1981, **52**, 7182-7190.

55 B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, *J. Comput. Chem.*, 1997, **18**, 1463-1472.

56 U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, *J. Chem. Phys.*, 1995, **103**, 8577-8593.

57 T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiff-
mann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, *J. Chem. Phys.*, 2020, **152**, 194103.

58 J. P. Perdew, K. Burke and M. Ernzerhof, *Phys. Rev. Lett.*, 1996, **77**, 3865–3868.

59 Y. Zhang and W. Yang, *Phys. Rev. Lett.*, 1998, **80**, 890–890.

60 J. P. Perdew, K. Burke and M. Ernzerhof, *Phys. Rev. Lett.*, 1998, **80**, 891–891.

61 S. Grimme, *J. Comput. Chem.*, 2006, **27**, 1787–1799.

62 S. Grimme, J. Antony, S. Ehrlich and H. Krieg, *J. Chem. Phys.*, 2010, **132**, 154104.

63 S. Grimme, S. Ehrlich and L. Goerigk, *J. Comput. Chem.*, 2011, **32**, 1456–1465.

64 S. Goedecker, M. Teter and J. Hutter, *Phys. Rev. B*, 1996, **54**, 1703–1710.

65 M. Krack, *Theor. Chem. Acc.*, 2005, **114**, 145–152.

66 J. VandeVondele and J. Hutter, *J. Chem. Phys.*, 2007, **127**, 114105.

67 G. Lippert, J. Hutter and M. Parrinello, *Mol. Phys.*, 1997, **92**, 477–488.

68 W. van Gunsteren and H. Berendsen, *Mol. Phys.*, 1982, **45**, 637–647.

69 G. Bussi, D. Donadio and M. Parrinello, *J. Chem. Phys.*, 2007, **126**, 014101.

70 W. Humphrey, A. Dalke and K. Schulten, *J. Mol. Graph.*, 1996, **14**, 33–38.

71 P. E. Mason, E. Wernersson and P. Jungwirth, *J. Phys. Chem. B*, 2012, **116**, 8145–8153.

72 M. Vazdar, P. Jungwirth and P. E. Mason, *J. Phys. Chem. B*, 2013, **117**, 1844–1848.

73 E. Pluhařová, P. E. Mason and P. Jungwirth, *J. Phys. Chem. A*, 2013, **117**, 11766–11773.

74 M. Kohagen, E. Pluhařová, P. E. Mason and P. Jungwirth, *J. Phys. Chem. Lett.*, 2015, **6**, 1563–1567.

75 M. Kohagen, P. E. Mason and P. Jungwirth, *J. Phys. Chem. B*, 2016, **120**, 1454–1460.

76 E. Duboué-Dijon, P. E. Mason, H. E. Fischer and P. Jungwirth, *J. Chem. Phys.*, 2017, **146**, 185102.

77 T. Martinek, E. Duboué-Dijon, T. Štěpán, P. E. Mason, K. Baxová, H. E. Fischer, B. Schmidt, E. Pluhařová and P. Jungwirth, *J. Chem. Phys.*, 2018, **148**, 222813.

78 E. Duboué-Dijon, P. E. Mason, H. E. Fischer and P. Jungwirth, *J. Phys. Chem. B*, 2018, **122**, 3296–3306.

79 P. E. Mason, L. Tavagnacco, M.-L. Saboungi, T. Hansen, H. E. Fischer, G. W. Neilson, T. Ichiye and J. W. Brady, *J. Phys. Chem. B*, 2019, **123**, 10807–10813.

80 D. Biriukov, H.-W. Wang, N. Rampal, C. Tempra, P. Kula, J. C. Neuefeind, A. G. Stack and M. Předota, *J. Chem. Phys.*, 2022, **156**, 194505.

81 D. Mendes de Oliveira, S. R. Zukowski, V. Palivec, J. Hénin, H. Martinez-Seara, D. Ben-Amotz, P. Jungwirth and E. Duboué-Dijon, *Phys. Chem. Chem. Phys.*, 2020, **22**, 24014–24027.