Single molecules of the solute guaiacol (2-methoxy phenol), Fig. 1, and ethanol (EtOH), were initially built in Avogadro22, pre-minimised using MMFF9423 and further optimised in Gaussian09 (v. 09, revision A.02, Gaussian, Inc., Wallingford, USA)24. The opt = tight convergence criteria was used as implemented in the software at the HF/6–31G* level of theory as recommended for use with the Amber family of force fields25. Atomic partial charges were derived through the restrained electrostatic potential (RESP) procedure26. Molecular coordinates representing a series of EtOH-water mixtures were built in PACKMOL27 with molar ratios of EtOH (X EtOH ) ranging from 0 to 1.0, using 0.1 increments, both in the absence and presence of the solute guaiacol. In the Supplementary Tables S1 and S2 more details are presented with regards to the mixture design, the starting box coordinates and the coordinates obtained after 10 ns of simulation.

All-atom MD simulations were performed using the Amber software (v. 10, USCF, San Francisco, USA)28, 29 together with the Amber12SB and the compatible GAFF30 force fields, which were developed for organic molecules. The SPC/E31 water model, which includes the average effect of polarisation in polar liquids, was used because the model is known to correctly reproduce the density, radial distribution function and self-diffusion constant of water32. Moreover, the model has been shown to be adequate for simulations at the air-water interface33, 34. In order to test the validity of the adapted parameters, a series of reference simulations of guaiacol, the pure solvents, and of the water-EtOH mixtures in the absence of guaiacol were performed. Simulations at the melting temperature of guaiacol (T = 301.15 K) resulted in enthalpy of vaporisation of 67.6 kJ mol−1, which is in good agreement with the experimental value of 60.4 kJ mol−1 35. The corresponding value for EtOH at 298.15 K was calculated to 46.1 kJ mol−1, which is also in agreement with the experimental value of 42.3 kJ mol−1 35. The densities of the liquid-phases of guaiacol and ethanol were calculated to 1.15 g cm−3 (301.15 K) and 0.80 g cm−3 (298.15 K), which is in general agreement with the experimental data36 of 1.11 and 0.785 g cm−3, respectively.

Calculations of water-EtOH excess volumes for mixtures in the absence of guaiacol reproduced the experimental finding that solutions with a molar fraction of X EtOH = 0.3–04 (59–69 vol% of EtOH) give rise to the greatest deviation from ideal behaviour, see Supplementary Table S3 and Figure S1. As in earlier simulations reported by others37, however, the differences in excess volume were higher in the experiment than in the simulations. The concentration of EtOH in each mixture is also presented as volume-% (vol-%). The approximate value of vol-% of ethanol in each mixture was obtained from a calculation of the molecular volume of ethanol, as explained in Table S3 and added to the text to aid in the interpretation of data.

Likewise, estimations of the static dielectric constant for each water-EtOH mixture through calculations of the fluctuation of the total dipole moment (<M2>-<M>2)38, 39 reproduced the decrease of the dielectric constant as observed in dielectric relaxation experiments40. The absolute value of the obtained dielectric constant, however, was lower than in the experiments, see Supplementary Figure S2. Calculations of the Stokes-Einstein self-diffusion rates of water and ethanol in each mixture, both in absence and presence of guaiacol, revealed no influence of the presence of the solute in the mixtures, see Supplementary Table S4.

Bulk-phase simulations

Each system was initially energy minimised to remove unfavourable van der Waals contacts using a total of 10000 steps. These were divided into 5000 steepest descent steps followed by 5000 steps of conjugate gradient minimisation. After the energy minimisation, velocities for all atoms were randomly assigned from a Maxwell-Boltzmann distribution of velocities at T = 298.15 K, and the systems were simulated for 100 ps under NVT conditions (i.e. constant number of particles, volume and temperature). An additional simulation of 500 ps was then conducted at a condition of constant number of particles, a pressure of 1 bar and a temperature of T = 298.15 K. The production phase data were collected for each system from 10 ns of simulation at NPT conditions.

During all simulations, the temperature was held constant using Langevin dynamics with a collision frequency set to 1.0 ps−1. The NPT simulation steps were conducted using an isotropic pressure scaling with the Berendsen barostat and a pressure relaxation constant τ p of 2.0 ps. All bonds to hydrogen were constrained using the SHAKE algorithm41 that allowed a time step set to 0.002 ps. Periodic boundary conditions were applied in all directions using a 12-Å non-bonded interaction cut-off. Long-range electrostatics were treated using the particle mesh Ewald (PME) summation method, and long-range van der Waals interactions were corrected using a continuum model correction of both energy and pressure. Data was collected every 2 ps.

Liquid-air interface simulations

Boxes representing molecular coordinates of the series of water-EtOH-guaiacol mixtures that were obtained after 10 ns of simulations were inserted and centred as a liquid-slab in new boxes. Each of these boxes had the same x and y dimensions as before, whereas the z-axis was elongated by half its dimension at its positive and negative sides. The mixtures were then subsequently simulated at conditions of constant volume for further 50 ns typically collecting data every 10 ps.

Analysis

Spatial distribution functions (SDFs) were calculated using the Grid command implemented in the AmberTools suite of programs (v. 15, USCF, San Francisco, USA). The positions of selected atoms of water and EtOH around guaiacol were binned from root-mean-square coordinate fit frames over all guaiacol atoms at intervals of 10 ps intervals for 50 ns liquid-air interface simulations, and at intervals of 2 ps for 10 ns bulk-phase simulations. A grid spacing of (0.5)3 Å3 and a box size of 1003 Å3 were used in all calculations.

Graphical representations of the occupancy of selected solvent atoms around the solute guaiacol were shown using either 3 or 8 times the contour value expected for the bulk density of the pure solvents. The values used as reference bulk densities (Density Bulk ) of water and ethanol in all mixtures were calculated to 20.9 (0.03346 atoms/Å3 × 0.125 Å3 × 5000 frames) and 6.5 (0.01046 × 0.125 Å3 × 5000 frames), respectively.

One-dimensional densities were calculated with the Density Profile Tool42 as implemented in VMD (version 1.9.1, University of Illinois at Urbana-Champaign, USA)43.

All z-axis atomic densities were calculated as the mean ± the standard deviation over the total simulation time of 50 ns with the centre of the box along the z-axis set to z = 0, taking into account each density value along the z-axis by using absolute values of the z-axis coordinates. For reasons of clarity, the final density profile has been symmetrized and shows both negative and positive values of z.

When error estimates are reported, values are presented as the mean ± standard deviation from five separate 30 ns blocks of data covering the total simulation time of 50 ns (0–30 ns, 5–35 ns, 10–40 ns, 15–45 ns, and 20–50 ns).