Propensity, free energy contributions and conformation of primary n-alcohols at a water surface

Abstract

Atmospheric aerosols contain organic molecules that serve as cloud condensation nucleation sites and affect the climate. Several experimental and simulation studies have been dedicated to investigate their surface propensity, but the mechanisms that drive them to the water surface are still not fully understood. In this molecular dynamics (MD) simulation study, primary alcohols are considered as a model system representing polar organic molecules. We find that the surface affinity of n-alcohols increases linearly with the length of the hydrophobic tail. By decomposing the adsorption free energy into enthalpy and entropy contributions, we find that the transition from bulk to surface is entropically driven, compatible with the fact that the hydrophobic effect of small solutes is of entropic origin. The enthalpy of surface adsorption is nearly invariant among different n-alcohols because the loss of solvent–alcohol interactions is balanced by a gain in solvent–solvent interactions. Structural analysis shows that, at the surface, the linear alcohols prefer an orientation with the hydrophobic tail pointing out from the surface, whereas the hydroxyl group remains buried in the water. This general behaviour is likely transferable to other small molecules with similar structures but other functional groups that are present in the atmosphere. Therefore, the present study is a step forward toward a general description of organic molecules in aerosols.

---

1 Introduction

Water covers a majority of our planet. As wind and waves whip the water surfaces, droplets and mist are created and carried up into the troposphere – the layer of the atmosphere closest to the surface of Earth. The tropospheric water has an important role in the climate on Earth and affects all atmospheric chemistry. Water originating from the oceans, lakes or seas carry organic material into the atmosphere,^1–6 that affect the chemical processes at the surface of the water droplets. Many important chemical reactions take place at these water–air interfaces,^7,8 and their composition therefore indirectly influences the climate.⁹ Apart from effects on atmospheric chemistry, organic molecules affect the climate by serving as nucleation sites for cloud condensation. As a consequence of its importance, organic molecules on water aerosols and their effect on atmospheric chemistry has been studied for decades.^1–6

Surface-sensitive X-ray photo-electron spectroscopy (XPS) has been widely used to study water surfaces.^10–17 Many of these studies have investigated the surface propensity of organic molecules, such as alcohols,^13,14,18 carboxylic acids,¹⁷ alkyle amines¹⁷ and esters.¹⁶

As a complement to experiments, molecular dynamics (MD) simulations have been used to study organic molecules at water interphases,^{13,14,16–21} providing an atomistic interpretation of the experimental signals. For instance, Walz et al.¹³ used MD simulations to reveal structural arrangements of surface-enriched pentanol that could give rise to experimental XPS spectra. However, MD simulations may suffer from inaccurate force fields or finite-size effects, and idealized simulations systems may not fully reflect more complex experimental conditions. Therefore, combining experimental observations with simulations is a fruitful approach to overcome limitations of each of the methods alone.

Here, using n-alcohols as a test case, we study how the length of a carbon chain modulates the thermodynamics of surface adsorption of organic molecules. We computed free energies, enthalpies and entropies of surface adsorption with MD simulations, and we compared our results with previous XPS experiments of the free energy of adsorption, ΔG_ads. By decomposing the enthalpy contributions further into solute–solvent and solvent–solvent components, we provide additional insight into the driving forces of surface propensity of primary linear alcohols.

In a previous study by Houriez et al.,²² carboxylate ions of increasing length were investigated at the water/vacuum interface by MD simulations using a polarizable model. Their results suggest that the hydrocarbon chain and ionic headgroup solvate independently, and that only ions with more than 4 carbons retain propensity in large water droplets. Like carboxylate ions, primary alcohols contain a polar headgroup and a hydrocarbon tail. The difference is that the polarity of the headgroup is not as pronounced, and the consequences of that is here to be addressed.

2 Methods

MD simulations were used to compute the adsorption free energy ΔG_ads of n-alcohols to a planar water–vapor interface. Two sets of seven different alcohol/water systems were simulated: (i) alcohols molecules in bulk water with no restraints on the alcohol molecule, and (ii) umbrella sampling simulations with the alcohol molecules restrained at different distances from a water surface to obtain ΔG_ads.

2.1 General simulation parameters

The GROMACS 2016.4²³ software was used for all simulations. Classical non-polarizable force fields have previously been successful in modeling experimentally observed behavior of organic molecules on water surfaces in multiple earlier studies.^13,16,19 We thus used the well-documented OPLS-AA force field²⁴ for the alcohols and the SPC/E water model²⁵ for our simulations. The combination of OPLS-AA and SPC/E water was previously found to perform excellently in reproducing experimental hydration enthalpies and entropies, as well as showing good performance for the calculation of solvation-free energies.²⁶

In both bulk and umbrella simulations, Newton's second equation of motion was solved using a leap-frog integration scheme with a time step of 2 fs. Periodic boundary conditions were applied in all directions. All bonds with hydrogen atoms within the alcohol molecules were converted to constraints and held at their equilibrium distance using the P-LINCS algorithm.²⁷ Water molecules were kept rigid with the SETTLE²⁸ algorithm. Temperature was implicitly controlled by a stochastic dynamics integrator with τ_t = 1.0 ps. Electrostatic interactions were computed with the particle-mesh Ewald method.^29,30 Lennard-Jones and direct-space Coulomb interactions were truncated at 1.2 nm. Dispersion correction of the energy due to long-range van der Waals interactions was applied.

2.2 Bulk simulations

Each alcohol molecule was centered inside a box with the shape of a rhombic dodecahedron with a minimum distance of 1 nm between the solute and the box edges and subsequently solvated by 250–300 water molecules using the gmx solvate module. The structure was minimized by steepest descent followed by an NPT MD simulation of 1 ns. These simulations were used to obtain basic structural parameters of the alcohol molecules in water, as listed in Table 1. The solvent accessible surface area and volume were computed with the gmx sasa module and the end-to-end distances with gmx distance.

Table 1 Structural properties of alcohols in water computed from bulk simulations

Molecule	Carbons [number]	End-to-end [nm]	Area [nm^2]	Volume [nm^3]
Methanol	1	0.19	1.63	0.186
Ethanol	2	0.29	1.98	0.243
Propanol	3	0.35	2.26	0.294
Butanol	4	0.49	2.59	0.349
Pentanol	5	0.60	2.90	0.402
Hexanol	6	0.73	3.21	0.455
Heptanol	7	0.81	3.51	0.508

---

2.3 Umbrella simulations

The simulation with a water slab was set up as follows. A cubic simulation box of 4 × 4 × 4 nm³ was filled with SPC/E water molecules²⁵ and equilibrated by means of an NPT simulations. The box extended in z-direction to 12 nm while keeping the water slab at the box center along the z direction. The adsorption free energies ΔG_ads were calculated with umbrella sampling and potential of mean force (PMF) calculations.³¹ The alcohol molecules, described in Table 1, were inserted into the center of the box, i.e. in the center of the water slab. Steepest descent energy minimization and a short NVT equilibration run of 100 ps was performed to relax the system to create starting structures for the umbrella simulations. The alcohol molecules were position restrained in the center during equilibration.

To avoid evaporation of water molecules, a spherical flat-bottomed quadratic position restraint acting on each water molecule was applied with a radius r_fb = 1.35 nm and force constant k_fb = 500 kJ mol⁻¹ nm⁻², using the following general potential. Here, H is the Heaviside step function and R_i is a reference position, from the initial, equilibrated structure.

The reaction coordinate r, was defined as the center of mass (COM) distance between of the water slab and the respective alcohol molecule in z-direction, i.e. normal to the water surface. A pictorial view of the reaction coordinate is shown in Fig. 1.

Fig. 1 Pictorial view of methanol at various positions along the reaction coordinate r.

Finally, for each molecule, 50 different positions from r = 0 to r = 5 nm, with 0.1 nm separation, were chosen as reference positions for the umbrella sampling. The position of the alcohol molecule was restrained along r by an harmonic umbrella potential with a force constant k = 1000 kJ mol⁻¹ nm⁻² using the umbrella pull code of GROMACS. Each windows was simulated 100 ns, while excluding the first 2 nanoseconds from analysis due to equilibration. From the series of umbrella simulations, the PMF along the reaction coordinate, with bulk free energy set to zero, was computed with the weighted histogram analysis method.³² Statistical errors were estimated by bootstrapping complete histograms as implemented in the gmx wham module of GROMACS 2016.4.³³

2.4 Free energy decomposition

The free energy of adsorption (ΔG_ads), was obtained from the PMF minimum, located slightly below the Gibbs dividing surface. It can be decomposed into contributions from enthalpy (ΔH_ads) and entropy (−TΔS_ads), as follows:

ΔG_ads = ΔH_ads − TΔS_ads (2)

ΔH_ads profiles were computed from the simulations by averaging the total potential energy within each umbrella window and henceforth defining ΔH_ads to zero in bulk water. Then, −TΔS_ads was computed from eqn (2). The ΔH_ads profile was further analyzed by decomposing it into contributions from solvent–solvent (ΔH_{solvent–solvent}), solvent–solute (ΔH_{solvent–solute}) and solute–solute (ΔH_{solute–solute}):

ΔH_ads = ΔH_{solvent–solvent} + ΔH_{solvent–solute} + ΔH_{solute–solute} (3)

3 Results and discussion

In the following, the free energy profiles (or potentials of mean force, PMFs) of seven primary alcohols of increasing length are presented. The contributions from enthalpy and entropy to the PMFs are analyzed, and the enthalpy is further decomposed into contributions from solute–solvent and solvent–solvent interactions. Finally, the preferred orientations of the alcohol molecules are presented.

3.1 Free energy of adsorption

The PMF along a reaction coordinate corresponding to the COM distance between a single alcohol molecule and the water slab was calculated from umbrella sampling simulations (Fig. 2). Evidently, the PMF of each alcohol molecule exhibits a pronounced free energy minimum in the vicinity of the surface, where the depth ΔG_ads of the miniumum corresponds to the (negative) adsorption free energy. As shown in Fig. 3, ΔG_ads decreases linearly with increasing chain length of the alcohol by 3.4 kJ mol⁻¹ per carbon atom and the location of the minimum is at 2.0 nm for methanol and gradually shifting toward 2.1 nm for heptanol. The contribution of a CH₂ group to the hydration free energy of an n-alkane chain was previously found to increase by 0.63 kJ mol⁻¹.³⁴ The desolvation of the alkane chain part of the alcohol can therefore only partially explain the adsorption, and other advantageous mechanisms must be involved.

Fig. 2 Potentials of mean force (PMFs) along the reaction coordinate r defined as the COM distance between water slab and alcohol molecule. The grey vertical line indicates the Gibbs dividing surface.

Fig. 3 Free energy of adsorption ΔG_ads as a function of number of carbons in alcohol molecule calculated from data in Fig. 2 and previous experimental data.^14,35,36 Both data sets fitted to linear functions.

The linear decrease in Fig. 3 is in qualitative agreement with previous experiments,^14,35,36 that show a decrease by 2.3 kJ mol⁻¹ per carbon atom. This means that the force field slightly overestimates the surface affinity of longer alcohols. In a separate series of control simulations, where free energy of hydration was calculated (see ESI,† Fig. S1), this discrepancy between experimental and simulated results was not seen, indicating that known drawbacks of the modelled SPC/E water surface (e.g. a surface tension underestimated by approximately 10 m N m⁻¹ compared to experimental values³⁷) are influencing the results.

3.1.1 Enthalpy penalty vs. entropy gain. To explore the thermodynamic driving forces of the observed surface propensities, the respective contributions from enthalpy and entropy were obtained from eqn (2), shown in Fig. 4. At the PMF minima at r ≈ 2 nm, the balance between enthalpy and entropy of the different solutes is depicted as stacked bars.

Fig. 4 Contributions to from enthalpy and entropy for the different solutes. Top: Along the reaction coordinate r. Bottom: Contributions to ΔG_ads, i.e. at the PMF minima.

The pattern is similar for all molecules; moving the alcohols from bulk water to the surface leads to a large gain in entropy (decrease of the −TΔS term) and, for most alcohols, only to a small increase of enthalpy ΔH. As the reaction coordinate increase above the Gibbs dividing surface, the enthalpy penalty becomes predominant, causing the PMF to increase above its bulk level. The enthalpy contributions to the PMF minimum are positive and seem independent of chain length above 3 carbons. Methanol, being a small molecule, is an exception and shows both negative enthalpy and entropy contribution at its minimum. The entropy contribution is substantially increasing with the number of carbons. This demonstrates that for longer alcohols (number of carbons ≥2), the increasing surface affinity is entropy-driven.

3.1.2 Enthalpy in the interfacial region. To rationalize why the enthalpy contributes much less to the surface affinity compared to the entropy, we further decomposed ΔH into contributions from solvent–solvent (ΔH_{solvent–solvent}), solute–solute (ΔH_{solute–solute}) and solute–solvent interactions (ΔH_{solute–solvent}). To this end, contributions from short-range Coulomb and Lennard-Jones interaction were averaged, according to the decomposition in eqn (3). The decomposition at the PMF minimum is shown in Fig. 5, where another balance act becomes evident. A negative contribution from solvent–solvent interactions on one hand counteracts the positive contribution from solute–solvent interactions on the other (Fig. 5, red and blue bars). Hence, during surface solvation, the large loss of solute–solvent interactions is nearly compensated by a gain in solvent–solvent interactions. Both sides increase in magnitude with increasing number of carbons in the solute, therefore the total enthalpy term remains constant. Notably, the change in intramolecular non-bonded energy (ΔH_{solute–solute}) was found to be negligible upon transition from bulk to surface (Fig. 5, yellow bars).

Fig. 5 Contributions to adsorption enthalpy (ΔH_ads) from solvent–solvent (blue), solvent–solute (red) and solute–solute (yellow, barely visible) for the different solutes.

3.1.3 Solute orientation and conformation. In order to investigate the orientation of the solutes, especially close to the surface, we calculated an orientation parameter, cos(θ), where θ is the angle between the end-to-end vector of each alcohol, pointing towards the hydroxyl group, and the normal vector of the water surface (the z-direction of the computational box). In Fig. 6, orientation of some of the solutes is presented as a distribution landscape of this orientation parameter. Evidently, the orientation is random (cos(θ) evenly distributed) when the solutes are in the bulk (r < 1.5 nm). As the solutes approach the surface (1.5 nm < r< 3.0 nm), their orientation parameter is more likely to be close to −1, i.e. with the hydroxyl group pointing toward the surface and the tail stretched away from the slab into the vacuum. With increasing length of hydrocarbon tail, the effect is even more pronounced. For heptanol, it is especially clear in the direct vicinity of Gibbs dividing surface (r = 2.5 nm). The distance from the Gibbs dividing surface at which the molecule first becomes affected, is increasing from approximately 0.8 nm for methanol to approximately 1.3 nm for heptanol, indicating that longer alcohols are immediated oriented as their hydroxyl group is able to form contact with the water surface. Snapshots from umbrella simulations near the PMF miniumum at r = 2 nm are shown in Fig. 7, to provide a pictorial view of the oriented solutes at the surface.

Fig. 6 Population landscape (color bar shows likelihood) of orientational parameter cos(θ), where (θ is defined as the angle between the end-to-end vector of each solute, pointing towards the hydroxyl group, and the normal vector of the water surface (the z-direction of the computational box) from umbrella simulations of methanol, propanol, pentanol and heptanol at different COM distance r between solute and water slab. The Gibbs dividing surface of the water slab is found at r = 2.5 nm in all systems. Here, cos(θ) = −1 means that the end-to-end vector is aligned with the z-axis (slab normal) and the hydroxyl group is pointing towards the water.

Fig. 7 Snapshots of methanol, propanol, pentanol and heptanol at the surface.

Furthermore, the rotational conformation around carbon–carbon bonds in the alcohols with at least 3 carbon atoms were investigated, some of the results are presented in Fig. 8. The rotation around the carbon–carbon bond closest to the hydroxyl group and characterized by the dihedral angle C3–C2–C1–O (See Fig. 8) in all solutes showed gauche as well as trans conformation. Moreover, their distributions were affected by the location of the solute with respect to the surface, where trans conformation became increasingly populated near the surface, less pronounced for the shortest chain having this dihedral angle, i.e. propanol. Hence, at the surface, especially the larger alcohols preferred a more linear, elongated arrangement in the vicinity of their hydroxyl groups, thereby allowing them to remain fully solvated in water while avoiding unfavorable hydrophobic contacts between the hydrocarbon tail and water. This can also be seen in Fig. 9, showing that alcohols with 3 or more carbons increase their length as they approach the surface and adopt a more collapsed structure above it. The smaller alcohols, methanol and ethanol also stretch out upon surface approach, although more modestly. Ethanol collapses in vacuum as the larger alcohols, but methanol keeps its stretched conformation above the surface.

Fig. 8 Top: Distribution of the rotation about C1–C2 torsion angle in propanol, pentanol and heptanol at different reaction coordinate r indicated by the numbers in white corresponding to the rotation around the carbon–carbon bond closest to the hydroxyl group. Bottom: Structure of hexanol with atom names and the investigated torsion angle.

Fig. 9 The alcohol average end to end distance normalized to its length in bulk water at different reaction coordinates r. The Gibbs dividing surface is indicated by a vertical line.

Molecules with multiple carbon–carbon torsion angles (pentanol and heptanol) showed no dependence on location with respect to surface on their mid- and outer C–C bonds, where the trans conformation was the most populated at all reaction coordinates. A previous study on molecular conformation of alkane chains in water³⁴ showed an absence of solvation effects of torsional angles of alkane chains up to 20 carbons, and our results suggest that a similar behavior is seen in the hydrocarbon tails of short primary alcohols, but with the addition that the conformation close to the hydroxyl group is affected by the presence of a water surface.

4 Conclusions

In coherence with previous experimental results, this simulation study shows that the length of the carbon tail in organic molecules strongly affects the surface propensity for smaller linear alcohols with one to seven carbons. Both simulations and experiments show a linear relation between ΔG_ads and the length of the hydrocarbon chain. While previous simulation results for carboxylate ions showed propensity only when the number of carbons exceeded 3, even the shortest primary alcohol (methanol) is surface active, as shown in previous experiments^36,38 and corroborated by the simulations in this study. By disentangling contributions from enthalpy and entropy it was shown that increasing surface affinity (more negative ΔG_ads) with increasing tail length is driven by a large gain in entropy, whereas the change in enthalpy is positive and fairly constant for medium sized chains. The enthalpy term at the PMF minima turned out to be almost a zero sum game between a gain in solvent–solvent interactions (due to the restored water hydrogen-bond network) and a loss of solvent–solute interactions. Together with structural observations, the driving force of the increasing surface propensity is of entropic nature, nevertheless depending on the enthalpy balance. Climate modelling is a very complex task,³⁹ and uncertainties in the models affects our ability to predict the future climate. Improvement of the our understanding of the atmospheric chemistry is important when building reliable climate models in order to predict climate changes. The general behaviour described here – the free energy of adsorption decreases linearly with the number of carbons in the chain – is most likely transferable to molecules of similar size and structure but other functional groups. Hence, the results presented here suggest the possibility to create simple empirical models that can be used in more general descriptions of how different organic molecules are enriched at the surface of aerosols.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at HPC2N and UPPMAX partially funded by the Swedish Research Council through projects 2019-8-268 and 2019-3-680 and a Small Storage project 2020-16-22. C. C. acknowledges the Swedish Research Council (2018-00740) and the Helmholtz Association through the Center for Free-Electron Laser Science at DESY.

Notes and references

1. R. Duce, V. Mohnen, P. Zimmerman, D. Grosjean, W. Cautreels, R. Chatfield, R. Jaenicke, J. Ogren, E. Pellizzari and G. Wallace, Rev. Geophys., 1983, 21, 921–952.
2. J. Hoff, D. Mackay, R. Gillham and W. Shiu, Environ. Sci. Technol., 1993, 27, 2174–2180.
3. K. Hunter and P. Liss, Mar. Chem., 1977, 5, 361–379.
4. T. Novakov and J. Penner, Nature, 1993, 365, 823–826.
5. H. Tervahattu, J. Juhanoja and K. Kupiainen, J. Geophys. Res.: Atmos., 2002, 107, 4319.
6. H. Tervahattu, J. Juhanoja, V. Vaida, A. Tuck, J. Niemi, K. Kupiainen, M. Kulmala and H. Vehkamaki, J. Geophys. Res.: Atmos., 2005, 110, D06207.
7. Y. Jung and R. A. Marcus, J. Am. Chem. Soc., 2007, 129, 5492–5502.
8. O. Acevedo and K. Armacost, J. Am. Chem. Soc., 2010, 132, 1966–1975.
9. V. Vaida, J. Chem. Phys., 2011, 135, 020901.
10. N. Ottosson, E. Wernersson, J. Söderström, W. Pokapanich, S. Kaufmann, S. Svensson, I. Persson, G. Öhrwall and O. Björneholm, Phys. Chem. Chem. Phys., 2011, 13, 12261.
11. N. L. Prisle, N. Ottosson, G. Öhrwall, J. Söderström, M. Dal Maso and O. Björneholm, Atmos. Chem. Phys., 2012, 12, 12227–12242.
12. G. Öhrwall, N. L. Prisle, N. Ottosson, J. Werner, V. Ekholm, M.-M. Walz and O. Björneholm, J. Phys. Chem. B, 2015, 119, 4033–4040.
13. M.-M. Walz, C. Caleman, J. Werner, V. Ekholm, D. Lundberg, N. L. Prisle, G. Öhrwall and O. Björneholm, Phys. Chem. Chem. Phys., 2015, 17, 14036–14044.
14. M.-M. Walz, J. Werner, V. Ekholm, N. L. Prisle, G. Öhrwall and O. Björneholm, Phys. Chem. Chem. Phys., 2016, 18, 6648–6656.
15. J. Werner, M. Dalirian, M. M. Walz, V. Ekholm, U. Wideqvist, S. J. Lowe, G. Öhrwall, I. Persson, I. Riipinen and O. Björneholm, Environ. Sci. Technol., 2016, 50, 7434–7442.
16. V. Ekholm, C. Caleman, N. B. Prytz, M.-M. Walz, J. Werner, G. Öhrwall, J.-E. Rubensson and O. Björneholm, Phys. Chem. Chem. Phys., 2018, 20, 27185.
17. J. Werner, I. Persson, O. Björneholm, D. Kawecki, C.-M. Saak, M.-M. Walz, V. Ekholm, I. Unger, C. Valtl, C. Caleman, G. Öhrwall and N. L. Prisle, Phys. Chem. Chem. Phys., 2018, 20, 23281–23293.
18. S. Johansson, J. Lovrić, X. Kong, E. Thomson, P. Papagiannakopoulos, S. Briquez, C. Toubin and J. Pettersson, Phys. Chem. Chem. Phys., 2018, 21, 1141–1151.
19. V. Ekholm, M. Vazdar, P. Mason, E. Bialik, M.-M. Walz, G. Ohrwall, J. Werner, J.-E. Rubensson, P. Jungwirth and O. Björneholm, J. Chem. Phys., 2018, 148, 144508.
20. J. S. Hub, C. Caleman and D. van der Spoel, Phys. Chem. Chem. Phys., 2012, 14, 9537–9545.
21. R. Zubillaga, A. Labastida, B. Cruz, J. CarlosMartínez, E. Sánchez and J. Alejandre, J. Chem. Theory Comput., 2013, 9, 1611.
22. C. Houriez, M. Meot-Ner (Mautner) and M. Masella, J. Phys. Chem. B, 2015, 119, 12094–12107.
23. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
24. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
25. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1952, 1987(91), 6269–6271.
26. B. Hess and N. F. A. van der Vegt, J. Phys. Chem. B, 2006, 110, 17616–17626.
27. B. Hess, J. Chem. Theory Comput., 2008, 4, 116–122.
28. S. Miyamoto and P. Kollman, J. Comput. Chem., 1992, 13, 952–962.
29. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
30. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
31. G. Torrie and J. Valleau, J. Comput. Phys., 1977, 23, 187–199.
32. S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman and J. M. Rosenberg, J. Comput. Chem., 1992, 13, 1011–1021.
33. J. S. Hub, B. L. de Groot and D. van der Spoel, J. Chem. Theory Comput., 2010, 6, 3713–3720.
34. A. Ferguson, P. Debenedetti and A. Panagiotopoulus, J. Phys. Chem. B, 2009, 113, 6405–6414.
35. R. Marinho, M.-M. Walz, V. Ekholm, G. Öhrvall, O. Björneholm and A. de Brito, J. Phy. Chem. B, 2017, 121, 7916–7923.
36. H. Chen, W. Gan, R. Lu, Y. Guo and H.-F. Wang, J. Phys. Chem. B, 2005, 109, 8064–8075.
37. F. Chen and P. E. Smith, J. Chem. Phys., 2007, 126, 221101.
38. J. Sung, K. Park and D. Kim, J. Phys. Chem. B, 2005, 109, 18507–18514.
39. B. J. Soden, W. D. Collins and D. R. Feldman, Science, 2018, 361, 326–327.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1cp02818a

---