# Hydration of biologically relevant tetramethylammonium cation by neutron scattering and molecular dynamics

Phil Mason, Tomáš Martinek, Balázs Fábián, Mario Vazdar, Pavel Jungwirth, Ondrej Tichacek, Elise Duboué-Dijon, Hector Martinez-Seara

► To cite this version:

Phil Mason, Tomáš Martinek, Balázs Fábián, Mario Vazdar, Pavel Jungwirth, et al.. Hydration of biologically relevant tetramethylammonium cation by neutron scattering and molecular dynamics. Physical Chemistry Chemical Physics, In press, 10.1039/D3CP05449G. hal-04381431

HAL Id: hal-04381431

https://hal.science/hal-04381431v1

Submitted on 9 Jan 2024

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>
Cite this: DOI: 00.0000/xxxxxxxxxx

# Hydration of biologically relevant tetramethylammonium cation by neutron scattering and molecular dynamics†

Philip Mason*<sup>a</sup>, Tomas Martinek<sup>a</sup>, Balázs Fábián<sup>a‡</sup>, Mario Vazdar<sup>b</sup>, Pavel Jungwirth<sup>a</sup>, Ondrej Tichacek<sup>a</sup>, Elise Duboué-Dijon<sup>c</sup>, and Hector Martinez-Seara*<sup>a</sup>

Received Date

Accepted Date

DOI: 00.0000/xxxxxxxxxx

Neutron scattering and molecular dynamics studies were performed on a concentrated aqueous tetramethylammonium (TMA) chloride solution to gain insight into the hydration shell structure of TMA, which is relevant for understanding its behavior in biological contexts of, e.g., properties of phospholipid membrane headgroups or interactions between DNA and histones. Specifically, neutron diffraction with isotopic substitution experiments were performed on TMA and water hydrogens to extract the specific correlation between hydrogens in TMA ($H_{TMA}$) and hydrogens in water ($H_{W}$). Classical molecular dynamics simulations were performed to help interpret the experimental neutron scattering data. Comparison of the hydration structure and simulated neutron signals obtained with various force field flavors (e.g. overall charge, charge distribution, polarity of the CH bonds and geometry) allowed us to gain insight into how sensitive the TMA hydration structure is to such changes and how much the neutron signal can capture them. We show that certain aspects of the hydration, such as the correlation of the hydrogen on TMA to hydrogen on water, showed little dependence on the force field. In contrast, other correlations, such as the ion-ion interactions, showed more marked changes. Strikingly, the neutron scattering signal cannot discriminate between different hydration patterns. Finally, *ab initio* molecular dynamics was used to examine the three-dimensional hydration structure and thus to benchmark force field simulations. Overall, while neutron scattering has been previously successfully used to improve force fields, in the particular case of TMA we show that it has only limited value to fully determine the hydration structure, with other techniques such as *ab initio* MD being of a significant help.

## 1 Introduction

Tetramethylammonium (TMA) represents an important and ubiquitous motif in biological systems. It is found in the cellular membranes of almost all cells as the phosphorylcholine group in many phospholipids. Also, successive methylation of the amino acid lysine eventually results in a TMA functional group. This methylation is vital in histone-DNA binding and the epigenetic expression of DNA.¹ Tetraalkylammonium salts are also powerful denaturants² and are widely used as phase transfer catalysts³ and are also seen frequently in ionic liquids.⁴

TMA is one of the simplest and most spherical representatives of the so-called “hydrophobic ions”, as compared to more complex variants, such as tetraphenylphosphonium or tetrabutylammonium. TMA has short hydrophobic chains, so an intriguing question arises. How much does its hydration properties differ from those of a hypothetical perfectly spherical large ion?⁵⁻⁷ Small spherical alkali cations have a fairly generic hydration structure, with the first hydration shell water oxygen atom pointing towards the cation and the hydrogen atoms away. This orientation is the strongest for lithium, getting weaker upon moving down the periodic table. How does the solvation shell structure change by the time we get to a large cation like TMA? Also, given the hydrophobic character of the methyl groups, does hydrophobicity-driven cation-cation aggregation take place in aqueous TMA solutions?⁸

Neutron scattering experiments have provided insight into TMA hydration structure in TMA chloride or bromide solutions.⁵,⁶,⁹⁻¹¹ In aqueous solutions, the neutron scattering signal is dominated by water-water correlations. The structural correlations of interest, such as those between TMA and water nuclei, are thus hidden among the dominant water-water correlations.

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

<footnote><sup>b</sup> Department of Mathematics, Informatics, and Cybernetics, University of Chemistry and Technology Prague, Technická 5, 16628 Prague, Czech Republic</footnote>

<footnote><sup>c</sup> Université Paris Cité, CNRS, Laboratoire de Biochimie Théorique, 13 rue Pierre et Marie Curie, 75005, Paris, France</footnote>

<footnote><sup>‡</sup> Present address: Department of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany</footnote>

<footnote><sup>†</sup> Electronic Supplementary Information (ESI) available: ESI contains additional neutron scattering and AIMD analysis. See DOI: 10.1039/cXCP00000x/</footnote>
A way to deal with this problem is to employ neutron diffraction with isotopic substitution (NDIS).¹²⁻¹⁵ This method takes advantage of the fact that different isotopes have different neutron scattering properties and relies on the assumption that the mass of a nucleus does not affect the structure of the solution. This assumption of neglecting nuclear quantum effects has proven fairly robust even for materials with the largest isotope effects. Even H₂O and D₂O vary in number density by less than half a percent despite substituting two-thirds of the nuclei in the system. Still, hydrogen and deuterium are excellent nuclei to use in NDIS. Both are easily experimentally available and have among the largest neutron scattering contrasts for any element. The technique requires two identical solutions to be prepared, which differ only in isotopic concentration of one nucleus. As the structure of these solutions is assumed to be identical, the subtraction of their diffraction patterns cancels out for all components unrelated to the substituted nucleus. With its 12 identical hydrogen nuclei and the large contrast between H and D isotopes, the TMA ion is an ideal candidate to be studied by neutron scattering.⁵,⁹,¹⁰ However, even if NDIS makes isolating a few structural correlations of interest possible, not all pairwise distributions are equally easy to understand. The intuitive way of examining TMA hydration is to look at the spatial distribution of water oxygen atoms around the central TMA nitrogen atom. Unfortunately, NDIS cannot provide this specific structure factor due to lack of suitable oxygen isotopes. The relatively easily experimentally accessible HₜMA-Hw structure factor corresponds neither to the center of the TMA nor that of the water molecule, which makes intuitive interpretations of the data significantly more complicated.

Despite these limitations, early NDIS experiments⁵,⁹,¹⁰,¹⁶ were interpreted as evidence of apolar hydration around TMA, with an edge-on orientation of water molecules. More recently, Monte-Carlo-based empirical potential structure refinement (ESPR) simulations were used to assist the interpretation of the neutron scattering signal, confirming the apolar character of TMA hydration but suggesting that water arranges tetrahedrally around TMA.¹¹ In the past years, molecular dynamics simulations have proven to be highly valuable to help further the interpretation of neutron scattering signals,¹⁷,¹⁸ which is also the strategy adopted in the present work. In particular, we performed NDIS experiments on a 2 m TMACl solution, allowing us to extract a single structural correlation between hydrogens of the cation HₜMA and those of water Hw. We employed DFT-based ab initio molecular dynamics (AIMD) and force field molecular dynamics (FFMD) simulations to assist in interpreting the experimental signal. Simulations with various force fields were performed to test the sensitivity of the measured signal to various changes in the hydration structure, which turned out to be the key to understanding the exact amount of information contained in the neutron scattering signal and avoiding over-interpretations. In addition, we obtained new insights into optimizing force fields to better represent the hydration shells around solutes.

## 2 Methods

### 2.1 Neutron scattering measurements

NDIS measurements were performed 23 °C using the D4C diffractometer at the nuclear reactor at the Institut Laue-Langevin in Grenoble, France.¹⁹ All the samples were loaded into the same cylindrical null scattering titanium-zirconium cell, loaded in an identical geometry in the diffractometer. The sample diameter was 5.0 mm, the wall thickness 0.75 mm and the beam height 24 mm. The neutron wavelength was 0.4985 Å. Four chemically identical solutions of 2 m TMACl in water were prepared, which differed only in H/D substitution on the TMA (h12-TMA and d12-TMA) and H/D substitution on water. The four diffraction patterns (Fig. 1a,b) were recorded for about 2 hrs for each D₂O solution and for 4 hrs for each H₂O solution²⁰ The results were then corrected for multiple scattering and absorption and normalized against a standard vanadium scatterer.²¹

Taking the difference between the diffraction patterns associated with solutions that differ only by the H/D substitution on TMA (both in D₂O and H₂O solutions) yields the first-order differences $\Delta S_{H_{non}}^{X_{D2O}}(Q)$ and $\Delta S_{H_{non}}^{X_{H2O}}(Q)$ (Fig. 1c), that report on the correlation between non-exchangeable H on TMA and every other atom (X) in the system. They are respectively defined as (in units of mbarns):

$$ \Delta S_{H_{non}}^{X_{D2O}}(Q) = 90.2 \cdot S_{H_{TMA}H_{w}}(Q) + 39.3 \cdot S_{H_{TMA}O}(Q) \\ + 6.5 \cdot S_{H_{TMA}C}(Q) + 2.3 \cdot S_{H_{TMA}N}(Q) \\ + 2.3 \cdot S_{H_{TMA}Cl}(Q) + 4.3 \cdot S_{H_{TMA}H_{TMA}}(Q) - 144.8, \qquad (1) $$

$$ \Delta S_{H_{non}}^{X_{H2O}}(Q) = -50.6 \cdot S_{H_{TMA}H_{w}}(Q) + 39.3 \cdot S_{H_{TMA}O}(Q) \\ + 6.5 \cdot S_{H_{TMA}C}(Q) + 2.3 \cdot S_{H_{TMA}N}(Q) \\ + 2.3 \cdot S_{H_{TMA}Cl}(Q) + 4.3 \cdot S_{H_{TMA}H_{TMA}}(Q) - 4.0. \qquad (2) $$

These prefactors were calculated from the atomic concentration and neutron scattering lengths of the various elements in the system by standard literature methods.²²

The difference between eq. 1 and eq. 2 yields the second order difference $\Delta\Delta S_{H_{non}}(Q)$ (Fig. 1d), which reports on the single correlation between the TMA non-exchangeable H atoms and the water H atoms.

$$ \Delta\Delta S_{H_{non}}(Q) = \Delta S_{H_{non}}^{X_{D2O}}(Q) - \Delta S_{H_{non}}^{X_{H2O}}(Q) = 140.8 \cdot (S_{H_{TMA}H_{w}}(Q) - 1) \qquad (3) $$

This function provides a useful internal consistency check for the accuracy of the solutions and the multiple scattering and absorption corrections performed on the data. Due to the large inelastic scattering of ¹H and the Placzek effect,²³ hydrogen-containing samples always present a dominant background. The higher the atomic concentration of ¹H, the larger the effect, primarily visible in light water samples such as in Fig. 1a. The effect
<figure><img src="image_2.png" /><figcaption>Fig. 1 a) Total diffraction patterns for H₂O solutions of d12-TMACl(black) and h12-TMACl(grey). b) Total diffraction patterns for D₂O solutions of d12-TMACl(black) and h12-TMACl(grey). c) First order differences $\Delta S_{H_{non}}^{X_{H_2O}}(Q)$ (red line, obtained by the difference of the two diffraction patterns shown in a) and $\Delta S_{H_{non}}^{X_{D_2O}}(Q)$ (blue line, obtained by the difference of the two diffraction patterns shown in b). d) Second order difference $140.8 \cdot (S_{H_{TMA,H_2O}}(Q) - 1)$, obtained through the difference of the two first-order differences shown in c).</figcaption></figure>

is greatly diminished for the heavy water samples, as shown in Fig. 1b. The amount of inelastic scattering is largely determined by the number of ¹H nuclei per unit volume, so the first-order differences (Fig. 1c) should have exactly the same Placzek background, which should vanish completely in the second order difference (Fig. 1d). The largely constant value of the second-order difference in Fig. 1d indicates the absence of a visible background and proves that the four solutions were prepared with the same chemical composition.

## 2.2 Force field molecular dynamics (FFMD) simulations

We performed classical FFMD simulations of seven systems comprising 50 TMA cations or, for comparison, 50 neutral neopentane molecules matching the experimental concentration of 2 m. Force field details for each system are available in Table 1. Systems containing TMA were neutralized with chloride counterions of the CHARMMCl_o type²⁴ (−1 charge) or CL_2s type 13 (−0.75 scaled charge). Each system contained 1388 TIP3P water molecules with Lennard-Jones parameters also on the hydrogens²⁵. All systems were assembled using GROMACS 2021.2 and 2021.5 tools.²⁶

For the TMA moiety, several parametrizations used in the head group of lipids are available in the literature, which differ mainly in their partial charges on the three atom types of the TMA group. First, we employed the CHARMM36 parameters with the full +1 charge (denoted CHARMM). Next, we used a scaled variant, prosECCo, where charges were scaled by 0.75. This physi-

cally based Electronic Continuum Correction (ECC) approach accounts for electronic polarization in a mean-field way by scaling charges²⁷. ECC has been shown to improve significantly the description of ion pairing behavior in solution.²⁸,²⁹ Based on this scaled charge model with an overall charge of +0.75, we designed three additional variants with different charge distributions but the same overall +0.75 charge: 1) low CH dipole, where partial charges are reduced; 2) central-N, where the whole charge is concentrated on the nitrogen; 3) surface-H, where all charge is redistributed on the hydrogens (see Table 1). All these TMA models have the same Lennard-Jones parameters, as prescribed in CHARMM36.

For comparison, we also performed simulations with the neutral neopentane with the same geometry as TMA using the CHARMM36 force field. Note that this compound is highly hydrophobic and aggregates quickly when mixed in water. Finally, we modeled a single charged sphere similar to the coarse-grained model of TMA, center-bead. In this coarse-grained-like simulations, water still has an atomistic resolution while the entire tetramethylammonium group is replaced by a single “extended atom”. The size of this so-called bead is set to match the one of TMA in the first peak of its radial distribution with the surrounding water.³⁰

All systems were simulated with GROMACS 2021,²⁶ using simulation parameters as provided by CHARMM-GUI.³² The LINCS and Settle algorithms were used to constrain the geometry of TMA hydrogens³³ and water molecules, respectively, allowing the use of 2 fs time step in the simulations.³⁴ The coordinates were saved every 10 ps for each 2 µs simulation. For neopentane that aggregates due to its hydrophobic nature, we have simulated 200 ns biased simulations placing a lower wall between neopentanes at 0.75 nm distance (KAPPA=20000.0 EXP=2 EPS=1 OFFSET=0) using the lowest collective variable³⁵ to gather enough water orientation statistics around neopentane by excluding neopentane-neopentane contacts. Simulations were performed in the isothermal-isobaric (NpT) ensemble. Temperature was controlled using a Nosé-Hoover thermostat³⁶ with a time constant of 1 ps and a reference value of 310 K, and a constant pressure of 1 bar was maintained by an isotropically coupled Parrinello-Rahman³⁷ barostat with a time constant of 5 ps. Van der Waals interactions were treated using a cutoff of 1.2 nm with a force-switch at 1.0 nm using Verlet cutoff scheme³⁸ for neighbors. Long-range Coulomb interactions were accounted for using particle mesh Ewald (PME) with a 1.2 nm cutoff as implemented in GROMACS.²⁶ All simulation data can be found at zenodo.org (DOI: 10.5281/zenodo.10406618).

## 2.3 Ab Initio molecular dynamics (AIMD) simulations

As a benchmark for the TMA solvation structure, we complemented the neutron diffraction experiments and the force field-based simulations with Born-Oppenheimer ab initio molecular dynamics simulations (AIMD) of a single TMA cation with 64 water molecules under periodic boundary conditions. The large computational cost of AIMD simulations precludes using larger systems, such as those in FFMD. The present system was not neu-
<table><caption>Table 1 TMA and neopentane simulation models used in this study. In parenthesis, the used CHARMM36 atoms types³¹ are listed.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Geom.</th><th colspan="3">Atom (atom type) partial charge</th><th rowspan="2">Overall charge</th></tr><tr><th>N (NTL)</th><th>C (CTL5)</th><th>H (HL)</th></tr></thead><tbody><tr><td>CHARMM</td><td>TMA</td><td>-0.60</td><td>-0.35</td><td>0.25</td><td>+1.00</td></tr><tr><td>prosECCo</td><td>TMA</td><td>-0.61</td><td>-0.35</td><td>0.23</td><td>+0.75</td></tr><tr><td>low CH dipole</td><td>TMA</td><td>-0.05</td><td>-0.10</td><td>0.10</td><td>+0.75</td></tr><tr><td>center-N</td><td>TMA</td><td>+0.75</td><td>0.00</td><td>0.00</td><td>+0.75</td></tr><tr><td>surface-H</td><td>TMA</td><td>0.00</td><td>0.00</td><td>0.0625</td><td>+0.75</td></tr><tr><td colspan="6">center-bead (σ = 0.550 nm, ε = 0.83680 kJ/mol)</td></tr><tr><td>center-bead</td><td>sphere</td><td></td><td>+0.75</td><td></td><td>+0.75</td></tr><tr><td colspan="6">center-bead (C (CT), C (CT3), H (HA3))</td></tr><tr><td>neopentane</td><td>TMA</td><td>0.00</td><td>-0.30</td><td>0.10</td><td>0.0</td></tr></tbody></table>

tralized by any counterion. Force field molecular dynamics was used to preequilibrate the system and to prepare the initial configurations for the subsequent AIMD simulation, as follows, with a constant pressure simulation used to estimate the average size of the cubic simulation cell of 12.552 Å. A constant volume simulation was then used to prepare 10 initial configurations separated by 2 ns and equilibrated for 1 ns using FFMD. From these 10 equilibrated structures, AIMD simulations were performed using the generalized gradient approximation revPBE DFT functional³⁹⁻⁴¹ with the D3 dispersion correction⁴²⁻⁴⁴. Core electrons were replaced by GTH pseudopotentials,⁴⁵,⁴⁶, and the triple-ζ basis set TZV2P with polarization functions was used for valence electrons.⁴⁷ A cutoff of 400 Ry was used for the auxiliary plane wave basis set in the GPW method.⁴⁸ The system was first equilibrated using AIMD Langevin dynamics for at least 16 ps with a damping constant of γ = 0.02 ps⁻¹⁴⁹. During the production simulation, the temperature was set to 300 K via a global CVSR thermostat with a time constant of 1 ps.⁵⁰ To enhance sampling, ten 50 ps parallel AIMD simulations totaling 500 ps of trajectory were used for further structural analysis. All AIMD simulations were performed at constant volume, using the CP2K program package, version 7.1.⁵¹⁻⁵³ All simulation data can be found at zenodo.org (DOI: 10.5281/zenodo.10406618).

## 2.4 Density maps

Data from FFMD and AIMD simulations were processed with the help of an in-house developed software for unbiased alignment and density analysis. The analysis included atoms found within a 10 Å of the nitrogen atom in TMA, covering all TMA molecules in the system and all simulation frames. The TMA “neighborhoods” were aligned in two steps. First, a representative conformation was found for each simulation, and then the neighborhoods were aligned to that conformation using the positions of carbon atoms. In the case of center-bead approximation of TMA, the neighboring oxygen/chlorine atoms were used for the alignment instead. The alignment was performed similarly as in Ref. 54, that is, using a permutation-based unconstrained alignment, such that the reference atoms used for the alignment (TMA carbons or oxygens and chlorines within its first solvation shell) do not need to be labeled and sorted. This way, the thermal noise is uniformly distributed in the density maps of all systems.

## 3 Results and Discussion

Neutron scattering patterns were obtained for four identical TMACl solutions that differed only by their isotopic composition on the nonexchangeable H of TMA and the exchangeable water H (see Methods section). While the diffraction patterns of the solution are dominated by the signal coming from the water signal, taking differences between pairs of solutions leads to the cancellation of the signal part that does not depend on the isotopic composition. Hence, first-order differences (see Methods) report the TMA(H/D) correlation to every other atom in the system. While rich in information, this signal still contains all the intramolecular correlation peaks, which hides the information about hydration (see ESI). Hence, we proceeded to obtain the double difference signal $\Delta\Delta S_{H_{non}}(Q)$ (see Methods and Fig. 2), which reports on the single correlation between $H_{TMA}$ and $H_{W}$. It thus directly probes the TMA hydration structure and is much easier to interpret than the total diffraction patterns. The obtained $\Delta\Delta S_{H_{non}}(Q)$ exhibits neatly resolved features below 10 Å which characterize TMA hydration. Interpretation of the signal tends to be more intuitive in direct space. Hence, we computed the inverse Fourier-transform $\Delta\Delta G_{H_{non}}(r) = \mathcal{F}^{-1}[\Delta\Delta S_{H_{non}}(Q)]$, which is directly related to the single pair-correlation function $g_{H_{TMA}H_{W}}$:

$$ \Delta\Delta G_{H_{non}}(r) = \Delta G_{H_{non}}^{X_{D_2O}}(Q) - \Delta G_{H_{non}}^{X_{H_2O}}(Q) = 140.8 \cdot (g_{H_{TMA}H_{W}}(r) - 1) \quad (4) $$

$\Delta\Delta G_{H_{non}}(r)$ shown in Fig. 2 presents a characteristic shoulder around 3–4 Å, and a peak at 6 Å. While direct molecular interpretation of these features is not straightforward, the same pair correlation function can be easily computed from molecular dynamics simulations with two objectives in mind. First, we want to assist in interpreting the neutron diffraction signal and obtain a molecular-level picture of the TMA hydration. Second, we want to investigate how sensitive TMA hydration is to variations in the intermolecular interactions, and whether the resulting patterns in hydration structure would be distinguishable in the neutron signal. To this aim, we performed FFMD simulations of a TMACl solution, at the same 2 m concentration as used in the neutron scattering experiment, using different force field variants (see Methods, Table 1), characterizing in each case the hydration structure and the associated neutron signal.
<figure><img src="image_3.png" /><figcaption>Fig. 2 Upper the function $\Delta\Delta S_{H_{non}}(Q)$ (grey) and lower the direct Fourier transform of this function $\Delta\Delta G_{H_{non}}(r)$ (grey). Shown in red is the function $\Delta\Delta S_{H_{non}}(Q)$ used for the rest of this paper, terminating the data using a windows function up to $\sim 10 \AA^{-1}$. The lower red function is the real space version of the upper red function.</figcaption></figure>

## 3.1 What is the effect of overall charge on TMA solvation?

TMA is a large cation with a relatively low charge density. Hence, a question arises on how different its hydration structure is from that of neutral solutes, such as neopentane (beyond the absence of counterions and possible aggregation, as discussed below). In addition, within standard force fields, TMA, as intuitively expected, is assigned a global charge of +1. However, electronic polarization, which further screens interactions between ions, is missing in FFMD simulations using non-polarizable force fields, which may lead to artefacts such as excessive ion pairing.²⁷,²⁸ A mean-field strategy to implicitly account for electronic polarization in FFMD simulations is the Electronic Continuum Correction (ECC) approach, which is mathematically equivalent to scaling partial charges by a factor $1/\sqrt{\epsilon_\infty} \approx 0.75$, where $\epsilon_\infty$ is the high-frequency dielectric constant of water. Hence, we designed an

ECC version of the TMA force field (denoted as prosECCo, see Table 1), with an overall charge of +0.75.²⁸,⁵⁵ Here, we compared the structure of the solution simulated using this scaled-charge force field with that obtained using the standard full charge CHARMM36 force field, as well as to that of a solution containing electrically neutral neopentane molecules for comparison. Note that neopentane is insoluble and starts precipitating fast in the simulations. Consequently, we used a biased simulation to avoid their aggregation.

For each simulation, we obtained radial distribution functions from the central nitrogen atom to the surrounding O, H<sub>w</sub>, and Cl<sup>-</sup>, as well as density maps of chloride, water oxygen, and water hydrogen around TMA (see Fig. 3). Fig. 3 clearly shows that changing the charge from +1 to +0.75 has only a minor effect on the hydration structure, resulting in nearly identical radial distribution functions. If we view the TMA ion as a tetrahedron with four faces, six edges, and four corners, the water oxygen atoms (as well as the chloride counterion) are located more at the faces, less at the edges, and not at all at the corners. The amount of TMA<sup>+</sup>-Cl<sup>-</sup> ion pairing is also very similar for the two systems, exhibiting only a small excess of solvent-shared ion pairs in the full charge force field compared to that with scaled charges (see second peak at 6-8 Å in g<sub>NCI</sub>(r)). This is contrary to what was previously observed for small monovalent cations such as Li<sup>+</sup>⁵⁶, and with divalent ions such as Ca<sup>2+</sup>⁵⁷, where charge scaling changed qualitatively the number of ion pairs. In contrast, due to the low charge density of TMA, only minor polarization effects are observed for ion hydration and pairing.

If we now compare TMA hydration with that of neutral neopentane, the main difference is due to a size effect with the hydration layer shifted slightly further away from the central atom. The orientation of the water molecules is also different – in neopentane, the water OH bonds point very slightly towards the central atom; in contrast, in TMA, the water dipole orients towards the ion, and the OH points slightly outwards. Hence, while much less strongly oriented than around smaller ions such as lithium, the hydration of TMA still appears in these simulations significantly different from that of a hydrophobic solute. Note that the observations discussed are robust, despite neopentane statistics being worse due to less sampling due to their strong clustering propensity.

## 3.2 Effect of charge distribution on TMA solvation

For a given overall charge, we investigated how the charge distribution within the TMA impacted its hydration structure and its propensity to form TMA–Cl ion pairs. For water as strong hydrogen bonding moiety, the partial charge on the central oxygen is about –0.8 with +0.4 on the hydrogen, while for TMA, the partial charges in most force fields are about –0.3 on the carbon and +0.2 on the hydrogen. Starting from the scaled charge prosECCo model with an overall TMA charge of +0.75, we compared four different charge distributions – the prosECCo force field (with results presented above), a low CH dipole force field variant where the charges on the methyl C and H atoms are respectively –0.1 and +0.1, a center-N force field where all the charge is placed on the central nitrogen, and a surface-H force
<figure><img src="image_4.png" /><figcaption>Fig. 3 Comparison charge effects in first solvation shell of TMA geometries. The first column shows charge +1.0 TMA (CHARMM), the second column +0.75 TMA (prosECCo), and the third column (neopentane). The first row shows chloride ions and the lower oxygen and hydrogen atoms. In all density maps, the density of chloride ion (green) is six times bulk density, oxygen (red) is three times bulk density, and hydrogen (white) is two times bulk density. The lower row shows the RDFs from the central atom for each model with oxygen (red), hydrogen (blue), and chloride (green). The fourth column compares the N-Cl⁻ and N/C-O RDFs for the three models.</figcaption></figure>

field where the charge is equally distributed over the surface hy-
drogens of the TMA only. Note that while the charge distribution
is rather different between the center-N and surface-H models,
they both result in a very low polarity of the C–H bonds and zero
or very small charge on hydrogens. In each of the four cases,
we characterized the hydration structure and ion pairing via the
radial distribution functions g$_{NH}$(r), g$_{NO}$(r) and g$_{NCI}$(r) and the
density maps of H$_{W}$, O and Cl around TMA (Fig. 4).

Our calculations show that changes in the charge distribution
within TMA (see Fig. 4 have a much larger effect on hydration
structure than the reduction of the overall charge of the ion from
+1.0 to +0.75 (see Fig. 3. Namely, lowering the charge on
the surface hydrogens qualitatively changes the hydration pat-
tern. The density maps (Fig. 4) clearly show that the water oxy-
gens (as well as the chloride counterions, which follow the same
trends) then move to the center of the faces of the TMA tetra-
hedron (i.e., away from the H atoms), forming bridges over the
tetrahedron edges. This is manifested in the radial distribution

functions g$_{NH}$(r), g$_{NO}$(r) and g$_{NCI}$(r) as a subtle increase in the
bimodality of the first peak for g$_{NO}$(r) and g$_{NCI}$(r) in center-N.
At the same time, changing the charge distribution does not sig-
nificantly change the orientation of the hydration water OH bonds
(see Fig. 4).

**3.3 Effect of molecular geometry on TMA solvation**

Coarse-grained models of TMA reduce the molecular geometry to
a single spherical bead. Is this loss of molecular structure impor-
tant or, in other words, how much does TMA behave as a sim-
ple charged sphere? To examine this issue we performed a set
of three simulations with the all-atom prosECCo force field, the
center-N force field variant where the charge is localized on the
central N atom, and the coarse-grained force field where all the
charge is at the center of a single spherical bead. Again, the den-
sity maps for O, H<sub>W</sub> and Cl were calculated around TMA, as well
as the radial distribution functions g<sub>NH</sub>(r), g<sub>NO</sub>(r) and g<sub>NCI</sub>(r), see
Fig. 5.
<figure><img src="image_5.png" /><figcaption>Fig. 4 For all density maps, the density of chloride ion (green) is 6x bulk density, oxygen (red) is 3x bulk density and hydrogen (white) is 2x bulk density. All the TMA cations in this figure have a +0.75 charge. Left column, prosECCo force field, center left, the polarity of the CH bond is comparable to that of neopentane, with the remaining charge on the central nitrogen. Center right, all of the charge of the ion is on the central nitrogen, and far right, all the charge is spread evenly on the hydrogen atoms. Lower are shown the RDFs from the central atom of the TMA to oxygen (red), hydrogen (blue), and chloride (green). Right are shown the same RDFs but grouped for each of the four FF. Upper for the chloride ion and lower for the oxygen atom.</figcaption></figure>

Interestingly, the all-atom prosECCo model is more similar in terms of the NO and N-Hw radial distribution functions to the coarse-grained charged sphere than to the center-N model (Fig. 5). Notably, the radial distribution functions of the center-N model are much more bimodal than the two others due to the combination of structural arrangements around the TMA tetrahedral features of faces, edges, and corners. This means that much of the TMA hydration structure reflects the constraints imposed by the charged sphere of a given size on the water H-bond network. At the same time, the atomistic TMA model is strikingly different from a simple charged sphere when looking at its interaction with the chloride counterion. In the former, the counterions adopt a tetrahedral geometry at the center of the faces similar to the oxygen atoms, while the latter tends to form linear Cl-TMA-Cl structures (Fig. 5).

## 3.4 Comparison of neutron scattering data to MD simulations

The double difference signal, $\Delta\Delta G_{H_{non}}(r)$, obtained from neutron scattering experiments after Fourier transform is composed of a single radial distribution function of $g_{H_{TMA}Hw}(r)$, which can be directly compared to the same radial distribution function computed from FFMD simulations with different force fields (Fig. 6). All investigated force fields (CHARMM: full charge, prosECCo:

scaled, and center-N: +0.75 on N) capture the location of the main peak at about 6 Å, even if the CHARMM force field seems to provide a somewhat worse fit than the two others. Small differences are visible in the low-$r$ range, but the comparison does not allow us to decide on the best force field. The steep rise of the experimental signal is slightly shifted for all force fields, and the shoulders around 4-5 Å while all slightly different, are never the same as in the experimental signal.

As neutron scattering data cannot differentiate between the different 3D hydration arrangements found by FFMD, we performed DFT-based AIMD simulations and used them as a benchmark. Since the obtained FFMD density maps differ qualitatively from each other, we aimed to employ the AIMD results to determine which force field is more accurate. The very high computing cost of such simulations (see Methods for more details) necessarily limits both the size of the simulated system (a single TMA ion in a small box of 64 water molecules) and the length of the simulations (500 ps). Due to the limited box size, the calculation of radial distribution functions is thus limited to the small $r$ range (Fig. 2 in ESI). We also obtained 3D density maps for water oxygen atoms around TMA. Despite limited statistics, these plots tend to show that the highest density of water oxygen atoms is located at the center of the faces of the TMA tetrahedron, with bridges across the sides. The hydration density maps from AIMD (see
<figure><img src="image_6.png" /><figcaption>Fig. 5 For all density maps, density of chloride ion (green) is 6x bulk density, oxygen (red) is 3x bulk density, and hydrogen (white) is 2x bulk density. All the ions in this figure have +0.75 charge. Left column, prosECCo force field, center all of the charge of the ion is on the central nitrogen, and right, as center but all the tetrahedral structure of the ion has been replaced by a single VDW sphere. Lower are shown the RDFs from the central atom of the TMA to oxygen (red), hydrogen (blue), and chloride (green). Right are shown the same RDFs but group for each three FFs. Upper for the chloride ion and lower for the oxygen atom.</figcaption></figure>

Fig. 7) display both similarities and differences with respect to the FFMD density maps. A feature common to all density maps is that the oxygen clouds found over the faces and edges are closer to TMA than those of the $H_W$ clouds. There is, however, a difference in how these clouds are arranged for the prosECCo and low CH dipole force fields. In each case, the $H_W$ clouds are similar in shape but vary such that the low CH dipole force field has a greater tendency to spill over the edges of the TMA tetrahedron. For the low CH dipole force field the O clouds match the orientation of the $H_W$ cloud, while for the prosECCo force field, an opposite pattern is observed with the triangles of the $H_W$ and O cloud being anticorrelated (Figure 5 and 7). Hydration of the AIMD TMA is more similar to that of low CH dipole than that of prosECCo, where prosECCo and CHARMM models have similar hydrations. Interestingly, at lower atomic densities, AIMD density maps show overlapping $H_W$ and $O_W$ density clouds at the corners of TMA, which is not replicated in any force field-based FFMD simulation. Taken all together, these results suggest that the orientation of water molecules in the TMA hydration shell is better

captured by low CH dipole than by the other force field variants tested here.

To further investigate the sensitivity of the neutron scattering signal to different aspects of TMA hydration, we compared $g_{H_{TMA}H_W}(r)$ (which is directly related to the experimentally measured quantity) among all the different variants of the TMA force fields. While we previously showed (Fig. 3) that different charge distributions lead to strikingly different hydration patterns around TMA and different ion-pairing behavior, resulting in different patterns visible at the $H_W$ density map around TMA (Fig. 8), these differences do not significantly modify the $g_{H_{TMA}H_W}(r)$ (see Fig. 8, where the $g_{H_{TMA}H_W}(r)$ computed with different force fields are compared). Neutron scattering experiments examining the $g_{H_{TMA}H_W}(r)$ correlation are thus unable to distinguish the different hydration patterns and ion pairing propensities that we have shown to exist for these different force fields. It is rather unexpected that these different three-dimensional water orientations around the different TMA force fields give such similar RDFs for $g_{H_{TMA}H_W}(r)$. The origin of the similarity of the RDFs de-
<figure><img src="image_7.png" /><figcaption>Fig. 6 The function $\Delta\Delta G_{H_{non}}(r)$ is shown in black, with the same function calculated from FFMD simulations in red for the CHARMM (upper), prosECCo (middle) and Center-N force fields (lower).</figcaption></figure>

spite the significant differences in the density maps is as follows. The hydrogen density clouds have symmetric ordering around the center of the TMA molecule. However, the substituted nuclei are not at this symmetry center, and there is a correlation between each of the clouds of H<sub>w</sub> density and each of the 12 substituted hydrogens, which makes the function g<sub>H<sub>TMAH<sub>w</sub></sub></sub>(r) rather broad.

<figure><img src="image_8.png" /><figcaption>Fig. 7 O<sub>w</sub> and H<sub>w</sub> density maps around TMA in various simulations. Upper O<sub>w</sub> density is shown in red at 3.3x bulk density (approximated for AIMD simulations) and lower 2.8x bulk density for O<sub>w</sub> (red) and 2.1x bulk density for H<sub>w</sub> (white). The orientation of the O<sub>w</sub> around TMA shows an inverse relationship between prosECCo and low CH dipole FFs (highlighted by blue triangles). The H<sub>w</sub> densities are very similar for prosECCo, low CH dipole, and AIMD (highlighted by green triangles). The O<sub>w</sub> densities for AIMD and low CH dipole are similar, however, none of the FFMDs replicate the H<sub>w</sub> and O<sub>w</sub> clouds off the corner of the TMA tetrahedron (highlighted in blue) at any density contour level.</figcaption></figure>

The g<sub>H<sub>TMAH<sub>w</sub></sub></sub>(r) component of the structural data, despite being a large fraction of the total scattering, is thus not very informative of the relevant ion hydration structure. Even if we had examined the correlations between the central nitrogen of the TMA (N<sub>TMA</sub>) and H<sub>w</sub>, the corresponding RDFs would still be all remarkably similar to each other for all of these force fields (see Fig. 6), and would have produced a far smaller NDIS signal. Only, if the correlation between the central nitrogen and the oxygen in water could be measured, this would enable differentiating between these force fields. Unfortunately, isolating this component of the total scattering signal with isotopic substitution is not possible, as no suitable oxygen isotopes exist for such an experiment.

## 4 Conclusions

Neutron scattering experiments with double isotopic substitution on both H<sub>TMA</sub> and H<sub>water</sub> were performed allowing us to single out the correlation between H<sub>TMA</sub> and H<sub>water</sub>. Thanks to the very high contrast of the H/D substitution and the large number of H atoms in TMA, the signal is well above the noise level. It should thus allow for a detailed characterization of the hydration structure. However, we show here that molecular interpretation of the experimental signal proves to be very challenging as all experimentally measurable RDFs (g<sub>H<sub>TMAH<sub>w</sub></sub></sub>(r) and g<sub>NH</sub>(r)) are not very sensitive to different hydration patterns caused by changes in the employed force fields. We thus can neither fully infer the most probable hydration structure from such a comparison nor validate the preferential choice of the tested force field variants. Nevertheless, the simulation results provide important insights into TMA hydration and its sensitivity to force field parameters. Shifting from CHARMM (+1.0) to prosECCo (+0.75) charges had a relatively small effect on the ion hydration. Changing the po-
<figure><img src="image_9.png" /><figcaption>Fig. 8 Density maps of water hydrogen atom around TMA, at twice bulk density. Shown to the right of each density map is the corresponding $H_{TMA}-H_W$ radial distribution function, compared to the same function for the prosECCo TMA force field to help comparison.</figcaption></figure>

larity of the C–H bond from –0.35 C and +0.23 H (prosECCo) to –0.1 and +0.1 (low CH dipole) largely inverts the hydration structure of oxygens around the TMA. AIMD simulations show a hydration structure very similar to low CH dipole at the faces and edges of the TMA tetrahedron. The AIMD hydration structure seen at the corners of the TMA tetrahedron is, however, not replicated by any FFMD. Thus, we suggest that the C–H bond polarity in standard CHARMM and its variants is too high to capture the TMA hydration properly. Removal of the tetrahedral structure by employing a center-bead force field for TMA and converting it to a single large bead has a relatively minor effect on the radial hydration structure of the ion. Still, it strongly changes the form of the counterion interaction. Namely, our findings imply that longer-range ion–ion ordered structures may not be accurately replicated as the structure of TMA is simplified to the level of a single bead within a coarse-grained force field. All these results suggest that caution should be taken when simulating moieties with TMA groups where hydration may play a relevant role, such as for common phospholipids and methylated lysines extensively present in biological systems.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

This work was performed using neutron beam time allocated at the Institut Laue Langevin (proposal 8–3–844). The authors thank Henry Fischer for his help and support. H.M.-S. acknowledges the support of the Czech Science Foundation (project 19-19561S). O.T. acknowledges the Faculty of Mathematics and Physics of the Charles University (Prague, Czech Republic) where he is enrolled as a PhD student, and the International Max Planck Research School for “Many-Particle Systems in Structured Environments” (Dresden, Germany) for support. M. V. acknowledges support by the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254), Project OPEN-28-18. EDD acknowledges support from the “Initiative d’Excellence” program from the French State (Grant “DY-NAMO”, ANR11-LABX-0011). P.J. acknowledges support from the European Research Council via an ERC Advanced Grant no. 101095957.

## Notes and references

1. P. Cheung and P. Lau, *Molecular Endocrinology*, 2005, **19**, 563–573.

2. P. H. von Hippel and K.-Y. Wong, *Journal of Biological Chemistry*, 1965, **240**, 3909–3923.
3 A. Brändström, *Advances in Physical Organic Chemistry*, 1977, 267–330.

4 A. Mustain, B. S. Gupta, M. Taha and M.-J. Lee, *New Journal of Chemistry*, 2023, **47**, 12304–12313.

5 J. Turner, A. K. Soper and J. L. Finney, *Molecular Physics*, 1990, **70**, 679–700.

6 N. G. Polydorou, J. D. Wicks and J. Z. Turner, *The Journal of Chemical Physics*, 1997, **107**, 197–204.

7 E. W. Lang, S. Bradl, W. Fink, H. Radkowitsch and D. Girlich, *Journal of Physics: Condensed Matter*, 1990, **2**, SA195–SA200.

8 D. Bhowmik, N. Malikova, G. Mériguet, O. Bernard, J. Teixeira and P. Turq, *Phys. Chem. Chem. Phys.*, 2014, **16**, 13447–13457.

9 J. Turner, A. Soper and J. Finney, *Molecular Physics*, 1992, **77**, 411–429.

10 J. Z. Turner, A. K. Soper and J. L. Finney, *The Journal of Chemical Physics*, 1995, **102**, 5438–5443.

11 E. J. Nilsson, V. Alfredsson, D. T. Bowron and K. J. Edler, *Physical Chemistry Chemical Physics*, 2016, **18**, 11193–11201.

12 P. H. K. D. Jong and G. W. Neilson, *J. Phys.: Condens. Matter*, 1996, **8**, 9275–9279.

13 P. E. Mason, G. W. Neilson, C. E. Dempsey and J. W. Brady, *Journal of the American Chemical Society*, 2006, **128**, 15136–15144.

14 P. E. Mason, S. Ansell and G. W. Neilson, *Journal of Physics Condensed Matter*, 2006, **18**, 8437–8447.

15 G. W. Neilson, P. E. Mason, S. Ramos and D. Sullivan, *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 2001, **359**, 1575–1591.

16 J. Turner and A. K. Soper, *The Journal of Chemical Physics*, 1994, **101**, 6116–6125.

17 E. Pluhařová, H. E. Fischer, P. E. Mason and P. Jungwirth, *Molecular Physics*, 2014, **112**, 1230–1240.

18 T. Martinek, E. Duboué-Dijon, Štěpán Timr, P. E. Mason, K. Baxová, H. E. Fischer, B. Schmidt, E. Pluhařová and P. Jungwirth, *Journal of Chemical Physics*, 2018, **148**, 222813.

19 H. E. Fischer, G. J. Cuello, P. Palleau, D. Feltin, A. C. Barnes, Y. S. Badyal and J. M. Simonson, *Applied Physics A: Materials Science and Processing*, 2002, **74**, s160–s162.

20 MASON Philip, Henry E. Fischer, JUNGWIRTH Pavel and TIMR Stepan, *Towards a fuller understanding of protein lipid interactions*, 2015.

21 G. J. Herdmants and G. W. Neilsont, *Journal of Physics: Condensed Matter*, 1996, **4**, 627–638.

22 J. E. Enderby, D. E. Williams and J. Randall, *Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences*, 1975, **345**, 107–117.

23 P. E. Mason, G. W. Neilson, C. E. Dempsey and J. W. Brady, *Journal of the American Chemical Society*, 2006, **128**, 15136–15144.

24 J. Chandrasekhar, D. C. Spellmeyer and W. L. Jorgensen, *Journal of the American Chemical Society*, 1984, **106**, 903–910.

25 A. D. MacKerell, D. Bashford, M. Bellott, R. L. 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, *The Journal of Physical Chemistry B*, 1998, **102**, 3586–3616.

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

27 I. Leontyev and A. Stuchebrukhov, *Physical Chemistry Chemical Physics*, 2011, **13**, 2613.

28 E. Duboué-Dijon, M. Javanainen, P. Delcroix, P. Jungwirth and H. Martinez-Seara, *The Journal of Chemical Physics*, 2020, **153**, 050901.

29 V. Kostal, P. E. Mason, H. Martinez-Seara and P. Jungwirth, *The Journal of Physical Chemistry Letters*, 2023, 4403–4408.

30 S. V. Bennun, M. I. Hoopes, C. Xing and R. Faller, *Chemistry and Physics of Lipids*, 2009, **159**, 59–66.

31 J. B. Klauda, R. M. Venable, J. A. Freites, J. W. O’Connor, D. J. Tobias, C. Mondragon-Ramirez, I. Vorobyov, A. D. MacKerell and R. W. Pastor, *The Journal of Physical Chemistry B*, 2010, **114**, 7830–7843.

32 J. Lee, X. Cheng, J. M. Swails, M. S. Yeom, P. K. Eastman, J. A. Lemkul, S. Wei, J. Buckner, J. C. Jeong, Y. Qi, S. Jo, V. S. Pande, D. A. Case, C. L. Brooks, A. D. MacKerell, J. B. Klauda and W. Im, *Journal of Chemical Theory and Computation*, 2016, **12**, 405–413.

33 B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, *Journal of Computational Chemistry*, 1997, **18**, 1463–1472.

34 S. Miyamoto and P. A. Kollman, *Journal of Computational Chemistry*, 1992, **13**, 952–962.

35 G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, *Computer Physics Communications*, 2014, **185**, 604–613.

36 S. Nosé, *Molecular Physics*, 1984, **52**, 255–268.

37 M. Parrinello and A. Rahman, *Journal of Applied Physics*, 1981, **52**, 7182–7190.

38 S. Páll and B. Hess, *Computer Physics Communications*, 2013, **184**, 2641–2650.

39 J. P. Perdew, K. Burke and M. Ernzerhof, *Physical Review Letters*, 1996, **77**, 3865–3868.

40 Y. Zhang and W. Yang, *Physical Review Letters*, 1998, **80**, 890–890.

41 J. P. Perdew, K. Burke and M. Ernzerhof, *Physical Review Letters*, 1998, **80**, 891–891.

42 S. Grimme, *Journal of Computational Chemistry*, 2006, **27**, 1787–1799.

43 S. Grimme, J. Antony, S. Ehrlich and H. Krieg, *Journal of Chemical Physics*, 2010, **132**, 154104.

44 S. Grimme, S. Ehrlich and L. Goerigk, *Journal of Computational Chemistry*, 2011, **32**, 1456–1465.

45 S. Goedecker, M. Teter and J. Hutter, *Physical Review B*, 1996, **54**, 1703–1710.

46 M. Krack, *Theoretical Chemistry Accounts*, 2005, **114**, 145–152.
47 J. VandeVondele and J. Hutter, *Journal of Chemical Physics*, 2007, **127**, 114105.

48 G. Lippert, J. Hutter and M. Parrinello, *Molecular Physics*, 1997, **92**, 477–488.

49 W. F. V. Gunsteren and H. J. C. Berendsen, *Molecular Physics*, 1982, **45**, 637–647.

50 G. Bussi, D. Donadio and M. Parrinello, *Journal of Chemical Physics*, 2007, **126**, 014101.

51 J. Hutter, M. Iannuzzi, F. Schiffmann and J. Vandevondele, *Wiley Interdisciplinary Reviews: Computational Molecular Science*, 2014, **4**, 15–25.

52 T. D. Kühne, M. Iannuzzi, M. D. Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belle-

flamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, *Journal of Chemical Physics*, 2020, **152**, 194103.

53 J. Vandevondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, *Computer Physics Communications*, 2005, **167**, 103–128.

54 M. T. H. Nguyen, O. Tichacek, H. Martinez-Seara, P. E. Mason and P. Jungwirth, *The Journal of Physical Chemistry B*, 2021, **125**, 3153–3162.

55 R. Nencini, C. Tempra, D. Biriukov, J. Polák, D. Ondo, J. Heyda, S. O. Ollila, M. Javanainen and H. Martinez-Seara, *Biophysical Journal*, 2022, **121**, 157a.

56 E. Pluhařová, P. E. Mason and P. Jungwirth, *The Journal of Physical Chemistry A*, 2013, **117**, 11766–11773.

57 M. Kohagen, P. E. Mason and P. Jungwirth, *Journal of Physical Chemistry B*, 2014, **118**, 7902–7909.