Hydrogen bonding and effects on the accommodation of methylamine at the air water interface

Hydrogen bonding and orientation effects on the accommodation of methylamine at the air-water interface Methylamine is an abundant amine compound detected in the atmosphere which can affect the nature of atmospheric aerosol surfaces, changing their chemical and optical properties. Molecular dynamics simulation results show that methylamine accommodation on water is close to unity with the hydrophilic head group solvated in the interfacial environment and the methyl group pointing into the air phase. A detailed analysis of the hydrogen bond network indicates stronger hydrogen bonds between water and the primary amine group at the interface, suggesting that atmospheric trace gases will likely react with the methyl group instead of the solvated amine site. These findings suggest new chemical pathways for methylamine acting on atmospheric aerosols in which the methyl group is the site of orientation specific c hemistry i nvolving i ts c onversion i nto a c arbonyl s ite providing hydrophilic groups for uptake of additional water. This conversion may explain the tendency of aged organic aerosols to form cloud condensation nuclei. At the same time, formation of NH2 radical and formaldehyde is suggested to be a new source for NH2 radicals at aerosol surfaces, other than by reaction of absorbed NH3. The results have general implications for the chemistry of other amphiphilic organics, amines in particular, at the surface of atmospherically relevant aerosols. I. INTRODUCTION Amines are one of the most common nitrogen-based organic compounds present in the atmosphere, appearing as a result of emissions from industrial processes, agricultural practices, animal husbandry, fires, and other anthropogenic and natural sources.1 Amines are semi-volatile derivatives of ammonia, in which one hydrogen (or more) atom has been replaced by an aryl or alkyl group, resulting in a molecular structure with a hydrophilic amine group attached to one (or more) hydrophobic ones. This substitution results in an amphiphilic character which makes these compounds easily detectable on the surface of liquid droplets.2–4 The overall impact of amines on our atmosphere is tremendous. Due to their capacity at acid-neutralization, amines are capable of modifying the acid/base nature of aerosol and cloud droplet interfaces5 and, consequentially, the formation of new aerosols at both the local and the large scale.6–8 A deeper understanding of the behavior associated with these compounds at the gas/liquid interface and of their accommodation at the surface of atmospheric particulates as a function of atmospheric conditions is needed to quantitatively address the impact of anthropogenic and natural emissions on the dynamical and chemical evolution of our atmosphere.9 After ammonia, methylamine (hereafter referred to as MA) is another abundant amine detected within the 1 atmosphere.10 The partitioning of MA into cloud droplets involves several steps: gas phase interaction with the water interface of the cloud droplet, transport across the vapor-liquid water interface, hydrolysis in the aqueous phase, aqueous phase diffusion, and chemical reaction inside the cloud droplet. The chemical species that reside at the boundary of atmospheric particulate play an important role in determining the fate and the chemical properties of atmospheric aerosols. For example, MA seems to promote the growth of secondary organic aerosols,11,12 SOAs, especially by interacting with glyoxal and methylglyoxal.6,13 Current chemical transport models for atmospheric modeling underestimate the formation of SOAs, a process which has an important impact on our atmosphere since SOAs can be converted to brown carbon particulates, i.e., aerosols capable of strongly absorbing sunlight in the near-UV and visible range.6,14 The reliability of weather, climate, and air-quality predictions is still generally affected by the lack of kinetic and thermodynamic data for heterogeneous mixtures and chemistry of atmospheric amines (MA especially).1,3,5,6 Apart from its atmospheric importance, MA is an interesting solute from a more general and physicochemical perspective since it combines one hydrophobic group (the methyl) connected to a more hydrophilic (the primary amine) one. Information about MA uptake, accommodation, and solvation environment would certainly generate fundamental knowledge that can be applied to other atmospheric amines, as well as other compounds which possess both hydrophobic and hydrophilic groups. In this work, we report novel results concerning the uptake, accommodation, and solvation environment of MA at the vapor/liquid-water interface from molecular dynamics (MD) simulations. We determine the mass accommodation coefficients for MA at the surface of liquid water together with a possible interpretation of the observed quantities in terms of the MA free energy of hydration and the molecular details of its solvation environment. The mass accommodation coefficients reported here can be immediately used as input in regional air-pollution models, addressing the impact of MA chemistry on air-quality, while the molecular insights on MA uptake can be used to draw a more general picture of amines and amphiphilic trace gas uptake in atmosphere and on more general interfaces. II. METHODOLOGY AND COMPUTATIONAL DETAILS Herein, we investigate the uptake and bulk/surface solvation conditions of MA at the surface of liquid water aerosols through the use of classical molecular dynamics simulations. MD simulations are based on an approximate description of the intramolecular and intermolecular potentials, such as the van der Waals (vdW) and Coulombic potentials. These approximations permit proper description of the dynamics at interfacial environments with feasible computational time. The set of physically relevant parameters required to describe both the intramolecular and the intermolecular interactions, as well as their functional forms, comprises the force field.15 As in previous studies on the solvation of organic molecules at air/condense water interfaces,16–21 in this work the intramolecular and vdW interaction parameters for MA were adopted from the Generalized Amber Force Field, GAFF.22 The vdW interactions were described through Lennard-Jones potentials with interaction sites located on all the atoms of MA. The electrostatic nonbonding interaction is described by a Coulombic interaction potential between individual partial charges located on each atomic site. Following common practices for GAFF, partial charges were calculated by fitting the electrostatic potential obtained from an ab-initio optimized molecular structure of MA at HF/6-31g* level using the Restrained Electrostatic Potential (RESP) method with a Merz-Singh-Kollman scheme.23 Gaussian0924 and Antechamber25 were used to perform both the ab-initio optimization and the charge fitting procedure. Table SI in the supplementary material reports the LennardJones parameters for MA,84 while Figure 1 reveals the charges obtained following the above procedure for MA. The TIP3P,26 TIP4P/2005,27 and TIP5P-Ew28 water models were selected to describe the liquid phase. TIP3P is a three-site water model with positive charges located on the hydrogen atoms and a negatively charged oxygen. TIP3P is a reference model for GAFF and is used for the simulation of organic molecules16–21 also of biophysical interest.29 In contrast, the TIP4P/2005 model shifts the negative charge from the oxygen atom onto a virtual site along the molecular bisector; this expedient resulted in a remarkably accurate description of the water phase diagram, and consequently correctly reproduces the liquid-vapor coexistence of water at the interface.30 Lastly, TIP5P-Ew employs two virtual sites FIG. 1. Charges for nitrogen (blue), carbon (light green), and hydrogens (light grey) in MA molecule, obtained by RESP procedure. corresponding to the lone pair of water molecule, for a total of 5 interaction sites (explicitly, these sites are oxygen and hydrogen atoms and lone pair sites). The presence of the lone pair sites is crucial to the resulting formation of the tetrahedral water arrangement, especially in solid phase.31 It is for that reason that the TIP5P-Ew model is widely used in simulations of ice,31,32 as it reproduces the exact melting temperature33,34 and the water diffusivity within the supercool regime.35 The structural and interaction parameters for these three models are reported in Table SII of the supplementary material.84 We adopted this selection of water models because they span a wide range of possible parameterizations (from a three-site water model to a five-site model) and applications (from simulations of interesting biophysical problems to those within environmental and atmospheric chemistry). The MD simulations conducted herein use nonpolarizable force fields, and therefore do not account for possible effects due to induced dipoles. Polarizability has been proven to be important for the description of air/liquid water interfaces, especially in presence of highly polarizable species (e.g., ions).36 However, MA is a poorly polarizable molecule and polarization is implicitly incorporated within all water models used in this work, as their permanent dipoles are enlarged compared to the gas phase one, taking into account the induced dipole in the condense phase.26–28 For these reasons, we performed MD simulations without polarizable force fields, where their use would have resulted in computationally expensive runs, limited simulation length and, thus, reduced sampling. The choice of force field must be properly validated through comparison of MD-derived predictions with either experimental measurements or high-level ab-initio quantum mechanics calculations for appropriately selective properties of a system. The free energy of hydration is a common benchmark in evaluating force fields because experimental values can be easily compared with those obtained by free energy methods in MD.16,37,38 Experimental free energy of hydration for organic molecules in liquid bulk water can be inferred from the Henry’s constant, K H , which is defined as the ratio between the molar concentration of the solute in the condense phase and its partial pressure in the vapor phase.39,40 K H is related to the free energy of solvation by39,40 ∆G0 = −RT ln 0 ) ( p0 K H cl0 , (1) where ∆G0 is the change in free energy associated with the transfer of a single solute molecule from the gas to the liquid bulk at standard conditions (i.e., gas pressure, p0 = 1 atm, and aqueous concentration cl0 = 1 mol/l, temperature of T = 298.15 K). Under standard conditions, the Henry’s 0 constant is referred to as K H . For comparison within this work, standard reference values of the Henry’s constant are taken from a compilation of such values by Sander.41 The free energy of hydration can be calculated in MD at the infinite dilution limit, i.e., computing the free energy cost of moving one MA molecule from the gas phase toward the bulk water. The MD-computed free energy of hydration, ∆G∗, differs from that calculated at standard conditions, ∆G0; these values are related through38,41,42 ∆G∗ = ∆G0 + RT ln ( p0 ) = ∆G0 − 7.9 kJ/mol. cl0 RT (2) Equation (2) was used to compare the MD-computed free energy and the experimental values, thus validating the choice of our force field. In order to avoid any possible confusion in reference and comparison, hereafter we shall report all free energy profiles as those obtained at infinite dilution conditions. We performed 200 MD simulations of 1 ns length each, impinging MA on the surface of liquid water with different velocities and angles, in order to calculate the mass accommodation coefficients. We determine four possible events for a molecule arriving at the surface, which are also schematically reported in Figure 2: scattering, adsorption, absorption, and desorption. The definitions above were inspired by similar studies of Veiceli et al.43 and Julin et al.,44 however, our criteria do present some notable differences, which are discussed in detail in the supplementary material.84 There are two possible definitions of the mass accommodation coefficients. The surface mass accommodation, α s , accounts for the fraction of sticking events, while α b is the fraction of events that resulted in an absorption of the impinging molecule within the liquid bulk. A detailed description of the definition of these two coefficients is also given in the supplementary material.84 For each water model considered in this work (TIP3P, TIP4P/2005, and TIP5P-Ew), we built a cubic box of 5 nm FIG. 2. Schematic view of the possible events associated to the impinging of MA on the air/liquid water interface. The possible events are scattering ((i) and (ii)), adsorption (iii), desorption (iv), and absorption (v). Further details are provided in the supplementary material.84 side length containing 4097 water molecules. Afterwards, we equilibrated these water boxes by performing 2 ns constant pressure and temperature (NPT) simulation at 1 bar and 298.15 K. At the end of the equilibration, the simulation boxes were elongated to 16 nm along the Z-direction, resulting in water slabs with two equilibrated air-liquid water interfaces exposed to the vacuum, and periodic boundary conditions in all directions. Starting from the equilibrated water slab interfaces, we calculated the free energy cost, ∆G∗, associated with moving a single MA molecule from the gas phase through the bulk water by combining MD simulation with umbrella sampling method. The procedure shares similarities to those utilized in previous studies.16,38,45,46 Firstly, we generated 160 umbrella configurations by pulling a MA molecule across the water slab, starting from a distance of 1.5 nm from the water-vacuum interface (see Figure 3) by applying a harmonic constraint with force constant of 1000 kJ mol−1 nm−2 and a pull rate of 0.005 nm/ps along the Z-direction over a 1.6 ns NVT run. Afterwards, for each 0.05 nm window, we performed 6 ns NVT simulations biasing the MA center of mass position with a harmonic constraint of 1000 kJ mol−1 along the z-direction. In all runs, the center of mass motion of the liquid water slab was removed every time step using a cosine weighted scheme.47 For each umbrella window, we built a population histogram recording the MA position over the last 5 ns of the NVT runs which was used to extract the free energy profile (Figure 4) by weighted histogram analysis method (WHAM).48,49 Bootstrapping analysis over a total of 100 cycles, implemented in the g-wham46 utility of GROMACS 4.6.7, was used to evaluate the statistical error in the reported free energy of hydration. The mass accommodation coefficient was calculated over the course of 200 NVT simulations at T = 298.15 K. These simulations were initialized with a MA molecule randomly placed at 2 nm distance from the air/water interface and assigned the three components of its center of mass velocity according to a Maxwell-Boltzmann probability distribution at the designated temperature. The total simulation time was 1 ns for each trajectory. The MA molecule was disentangled from the thermostat during these simulations; this was done in order to not interfere with the accommodation process and FIG. 3. Schematic view for the umbrella sampling simulation. One MA molecule was placed at a distance of 1.5 nm from the interface. Afterwards MA was pulled from the vacuum through the water bulk collecting the configurations for the umbrella procedure. The horizontal black lines indicate the simulation box sides: since the system is periodic in the XY in-plane direction, our simulation box mimics an infinite 2D water slab. FIG. 4. Free energy profile for MA solvated in TIP3P (red line), TIP4P/2005 (green line), and TIP5P-Ew (blue line) water slab. The water density profile (black), obtained collecting the Z position of all the water oxygen during the simulation time, is reported in arbitrary units and Z = 0 corresponds to the center of the water slab. The vertical dashed red line corresponds to the interfacial environment defined at the region between 10% and 90% of water density. ∆G ∗ and ∆G sl are the free energy of hydration and of transferring MA from the interface to the bulk, respectively. ∆G s the free energy difference between the interfacial minima and the reference gas value. to prevent velocity rescaling of MA by thermostatting.44 The average position over the last 20 ps of the 1 ns trajectory was used to determine whether the MA was adsorbed or absorbed. Moreover, the water slab system in our simulation box was subject to 3D standard periodic boundary conditions. These periodic boundary conditions permit observation of multiple events during the 1 ns run with, e.g., MA being desorbed from one interface and subsequently adsorbed by the other. We scanned the 200 MD trajectories, recording the number and type of events for the MA within our simulations. All MD simulations were performed using GROMACS 4.6.7,50 and using a time step of 2 fs. Temperature was controlled by v-rescale thermostat51 of coupling time of 0.1 ps. For the NPT runs the pressure was kept at 1 bar by a Berendsen pressure barostat52 with a coupling constant of 2 fs. The realspace Coulomb and the vdW interactions were cut-off at a distance of 1.0 nm in the case of water slab simulations with TIP3P and at 0.9 nm for water slabs using the TIP4P/2005 and TIP5P-Ew models, following the cut-off suggestions in the original references for these water models.26–28 The particle mesh Ewald method53 with a relative tolerance of 10−5, fourth order cubic interpolation, and a Fourier spacing parameter of 0.16 were used to evaluate the long range part of the Coulomb interactions. Consistently with GAFF practice,22 Lorentz-Berthelot combination rules were used to describe vdW interactions between solute and solvent. MA bonds were constrained using LINCS54 while SETTLE55 was used for constraining the water molecules. III. RESULTS AND DISCUSSION A. Free energy results Figure 4 shows the free energy profile associated with moving one MA molecule from the vacuum (used as reference zero-point value for the free energy) through the bulk. The figure also reports the water density profile (in arbitrary units) along the Z-direction (perpendicular to the interface) while collecting the position of each water’s oxygen atom during each step of the simulation time: the plateau in the water density profile is associated with the bulk region of the water slab while the steep portions indicate the interfacial regions. As the MA molecule approaches the interface, the free energy decreases reaching a constant value in the bulk. Now, returning from the bulk to the gas phase, the free energy profile again shows a minimum at the interface, and finally returning to zero in the gas phase. Figure 4 shows the hydration profile for MA in TIP3P (red curve), TIP4P/2005 (green), and TIP5P-Ew (blue) water slabs. These profiles are roughly symmetric with respect to the center of the water slab, which indicates that the sampling period was of sufficient length. The difference in free energy between the gas phase and the bulk defines the free energy of hydration, ∆G∗. ∆G sl is the free energy cost for the transition of an MA molecule from the interface to the bulk, while ∆G s is the free energy difference between the interfacial minima and the gas reference value. The values of ∆G s , ∆G sl , and ∆G∗ reported in Table I indicate that the maximum discrepancy in the free energy of hydration is between TIP4P/2005 and TIP5P-Ew water slabs, which is about 5.8 kJ/mol. However, the statistical uncertainty associated with ∆G∗ as reported by the single standard deviation is ∼2 kJ/mol; thus, all MD-computed ∆G∗ values agree fairly well with the experimental values as reported by the horizontal violet line within Figure 4. Moreover, inspection of Figure 4 suggests a higher presence of MA at the interface regions due to a minimum in the free energy that is present using all the three water models. Nevertheless, there are differences in the intensity of this interfacial minimum. The free energy profile for TIP4P/2005 experiences both the largest surface minima ∆G s = −26.0 kJ/mol and the largest ∆G sl = 5.2 kJ/mol, indicating that the interfacial enrichment of MA should be greater for a TIP4P/2005 water slab than one built with another water model. Free energy profiles, like the ones in Figure 4, are useful to determine the physicochemical behavior of a solute during the solvation process while additionally providing a validation test for the force field choice. Indeed, the comparison between the MD-computed and the experimental ∆G∗ is widely used to demonstrate the reliability of the interaction parameter used to describe physical and chemical systems in classical MD, especially for accurate modeling of gas uptake.16,17,37,38,43,56 It is of interest to note that, even if for all the three water models we have a reasonable description of the free energy of hydration, there are some differences in the intensity of the surface minima and for the free energy cost of transferring TABLE I. Free energy difference from gas to interface, ∆G s , free energy of hydration, ∆G ∗, and free energy of transition from water surface to bulk liquid, ∆G sl for MA in TIP3P, TIP4P/2005, and TIP5P-Ew water slab. Water model ∆G s (kJ/mol) ∆G sl (kJ/mol) ∆G ∗ (kJ/mol) TIP3P TIP4P/2005 TIP5P-Ew −20.2 −26.0 −18.1 2.5 5.2 3.0 −17.5 −20.8 −15.0 TABLE III. Residence times for desorption events observed during the collision of MA with TIP5P-Ew water slab. one MA from the interface to the gas; these differences suggest that ∆G∗ could not be (in general) the only target for force field benchmark in MD simulations of gas uptake on liquid water. In Secs. III B–III E we investigate the behavior of MA on TIP4P/2005 and TIP5P-Ew liquid water slab. These two models represent one of the state of the art models for simulations of liquid water and ice using classical non-polarizable MD.31,32,57,58 Moreover, the lower and higher values of free energy of hydration reported in Table I correspond for MA solvated in these water models. For all these reasons, these two models can provide us a good sample helping in drawing conclusions that can be as general as possible. Residence time, τ res (ps) 24 25 93 100 162 200 218 229 241 266 275 510 555 658 866 B. Mass accommodation coefficient The values for the surface and bulk accommodation, α s and α b , for MA colliding with a TIP4P/2005 or a TIP5PEw vapor/liquid water interface are reported in Table II. For both water models, the values for the surface mass accommodation, α s , agree fairly well and they are near unity: 0.93 and 0.97 for TIP4P/2005 and TIP5P-Ew, respectively. Interestingly, the bulk accommodation coefficients α b are not as uniform as for the surface. Indeed, in the case of a TIP4P/2005 water slab no desorption events were observed, which implies that all the adsorbed MA will eventually enter in the bulk, leading to an α b = α s = 0.93 (based on the definition of the two accommodation coefficients given in the supplementary material84). The impinging of MA toward the TIP5P-Ew water slab has permitted the recording of an α b = 0.85 with 15 desorption events and residence times on the liquid water interface reported in Table III. We note that in this case the residence time can range broadly from 24 ps to 866 ps; this finding comes in spite of the fact that we employed both a larger number of simulations and longer timed simulations compared to similar studies of gas molecule uptakes present in the literature.42–44 The statistics obtained for the residence time are insufficient to properly sketch a population distribution, and are therefore insufficient to draw any reliable conclusions concerning the average time of surface residence. Nevertheless, these results provide a range of values that can be used to better constrain nonreactive uptake coefficient in chemical transport models.59–61 An interpretation of the accommodation coefficients in Table II and of the desorption residence time in Table III can be made according to the free energy profile reported in Figure 4. The observation that no desorption events are recorded from impinging MA on TIP4P/2005 water slab correlates well with the deeper surface minima ∆G s detected for this water slab and all the adsorbed molecules will eventually enter the bulk, leading to α s = α b . This is a reasonable statement since the barrier to return to the gas phase ∆G s is much larger than that for the bulk ∆G sl . On the other hand, the ∆G s being lower (in module) for TIP5P-Ew, suggests some possibility of MA to escape the interface to the gas phase, leading to a smaller α b . However, the difference among these coefficients is small, and suggests caution concerning any quantitative conclusions being drawn. C. Surface preference The free energy pattern is, indeed, reflected in the probability distribution for the nitrogen and carbon positions of MA over the course of 200 ns NVT water slab simulation in Figure 5. The nitrogen and carbon atoms are likely to be preferentially located at the interfaces even though there is also the possibility for MA to explore the bulk region. The strongest free energy minima at interfaces for MA in TIP4P/2005 water slab show up in the density plot with higher peaks for the interfacial positions of MA, Figures 5(a) and 5(b). In the TIP5P-Ew water slab case, MA can more easily sample the bulk since the barrier between interface and bulk, ∆G sl is only few 3 kJ/mol, which can be overcome by thermal motion. Furthermore, Figure 5 shows that, even with 200 ns trajectory there is still a slight left-right asymmetry in the height of the density peaks for nitrogen and carbon in the water slab: this supports our choice of using free energy methods in order TABLE II. Number of scattering, desorption, adsorption and absorption events and the corresponding values for surface mass accommodation (α s ) and bulk accommodation (α b ) calculated over the course of 200 runs with MA colliding with TIP4P/2005 and TIP5P-Ew water slabs. The mass accommodation are reported with the 95% binomial confidential level. Water model Scattering Desorption Adsorption Absorption αs αb TIP4P/2005 TIP5P-Ew 15 7 0 15 136 92 64 108 0.93 ± 0.03 0.97 ± 0.02 0.93 ± 0.03 0.85 ± 0.10 with the presence in the MA molecule of a hydrophobic aliphatic core and of a hydrophilic primary amine group in the molecule. D. Orientation of MA at the interface FIG. 5. Probability distribution plots for the Z-position of the nitrogen (blue) and carbon atom (green) during the course of 200 ns NVT simulation on a TIP4P/2005, panel (a), and TIP5P-Ew, panel (b), water slab. The black density profile reports, in arbitrary units, the position of the water oxygen while the vertical dashed lines the interfacial environment defined at the region between 10% and 90% of water density. The solvation environment of MA at the interface can be further elucidated by looking at its orientation at the interface. Figure 6(a) shows the bivariate distribution for the angle between the MA molecular axis, i.e., the vector connecting the nitrogen to the carbon atom, and the interface normal as function of the Z-position of the nitrogen atom. Figure 6(a) indicates a strong peak in the probability distribution corresponding to MA sticking out perpendicularly from the interface, even if tilted angles up to 30◦ are as well explored. This implies that the methyl group sticks out to the gas phase once the primary amine group is solvated at the interfacial region. The fact that MA is also exploring tilted configurations means that, even if hydrophobic, the aliphatic group can still interact with the liquid water at interface. Figure 6(b) shows a snapshot of the trajectory with MA, indeed, tilted at the interface. The inlet in 6(b) also reports the regions explored by the water oxygen and hydrogen around MA at the interfacial region in the closest solvation shell. The inlet suggests that the nitrogen atom of MA molecules acts as a hydrogen bond donor and acceptor with interfacial water but the first solvation shell barely interacts with the aliphatic group. Only after including the next solvation shells is it possible to detect the presence of water around the aliphatic group. Thus, the interactions between the aliphatic group and the water are sufficient to slightly tilt its axis. E. Hydrogen bonding network surrounding MA to calculate the free energy of hydration instead of direct sampling techniques. Remarkably, the positions of the density peaks in Figure 5 show that the nitrogen atom remains closer to the liquid bulk than the carbon atom. This is consistent Figure 7 reports the h-bond correlation function, calculated according the procedure described in Ref. 62: if a water oxygen and the nitrogen of the primary amine group are bounded by an hydrogen bond at time t = 0, the FIG. 6. Bivariate and spatial distributions for MA in TIP5P-Ew liquid water slab. (a) Bivariate distribution for the tilting angle between the MA molecular axis and the interface normal cos(θ) and the position of the nitrogen atom, Z. The vertical solid lines indicate the boundaries of the interfacial region. (b) Simulation snapshot for methylamine at the surface of liquid water. The inlet shows the average position of hydrogen (white) and oxygen (red) atoms during the MD trajectory, showing the methyl group (light green) exposed to the air. FIG. 7. Hydrogen bond correlation function calculated as the life-time probability of an hydrogen bond. The top panels ((a) and (b)) show the correlation function for a h-bond between water oxygen and the amine nitrogen (N-acceptor) while the lower panels ((c) and (d)) for the h-bond between amine hydrogen and water oxygen (N-donor). Black circle corresponds to simulation of MA fully solvated in bulk water, red squares in water slab, and green diamond with MA constrained at the interface. Left panels ((a) and (c)) refers to TIP4P/2005 water and right panels ((b) and (d)) to TIP5P-Ew water. correlation function indicates the probability that a hydrogen bond will survive at a later time t. The lifetime of the h-bonds is a weighted measure of the importance of the energetic and entropic factors determining the average properties of the system. For the bulk, the availability of donors and acceptors is larger than at the interface. The inhomogeneity of the interfacial environment also affects the lifetime of the h-bonds. We have considered two specific cases: (i) the nitrogen atom in the amine group is an acceptor for water hydrogen, i.e., HOH–NH2 and (ii) the nitrogen of the primary amine group is hydrogen donor for the water oxygen (N-donor), i.e., HNH–OH2. A hydrogen bond is detected when the distance between the donor and the acceptor atoms is smaller than 0.35 nm and the angle between the hydrogen, the donor and the acceptor is smaller than 30◦.62,63 In Figure 7 we considered the strength of the hydrogen bond for MA in three different solvation environments: bulk, slab, and interface. For the h-bond in the bulk, we solvated MA in a fully 3D periodic bulk water box, while for the slab the MA was free to sample both the bulk and the interfacial environment. For the interfacial environment we performed a slab simulation constraining the Z-position of MA at the center of the interfacial regions using a mild harmonic constraint of force constant 700 kJ mol−1 nm−2. In this way, MA was able to sample all the interfacial environment, defined as the 10% and 90% of the bulk water density (see Figure 5), without leaving it. Figure 7 shows the results for MA in TIP4P/2005 and TIP5P-Ew water: for both water models, the h-bond lifetime is longer in the bulk (black circle curve) than at the interfacial environment (green diamond), while the slab case (red squares) shows an intermediate behavior. The fact that the h-bond lifetimes in the slab and the interface are similar is expected since MA is preferentially located in the interfacial environment as shown in Figure 5, indicating that the mild harmonic potential used to constrain the molecule at the interface does not significantly affect the h-bond network around MA. The observed difference between bulk and interface suggests that the asymmetry in the solvation environment at the interface significantly affects the h-bond lifetime. The longer h-bond lifetime at the interface could suggest that the water and MA interaction is energetically more favorable at interface, which correlates well with the interfacial minima observed in Figure 4. In order to clarify the possible link and to rationalize differences among used water models, we performed an electronic structure optimization for MA and one water molecule at MP2 level, running a steep energy minimization on the MP2 optimized structure using the classical force fields. The results are reported in Figure 8. In all these cases, the water molecule donates one hydrogen bond to the MA but small differences are observable depending on the water model used to describe the liquid phase. The nitrogen-water hydrogen (H–N) intermolecular distance is TABLE IV. Probabilities of different h-bond configurations in bulk and at the interface. XA indicates the number of h-bonds with the nitrogen of MA as acceptor while XD the number of h-bonds with nitrogen as donor. FIG. 8. Comparison among optimized structures of the MA and water dimer obtained at MP2/auc-cc-pvTz level, and steep energy minimization using TIP4P/2005 and TIP5P-Ew water models. The figure reports the nitrogenwater hydrogen, H–N, distance and the root mean square displacement (RMSD) of the atomic position fitted to the MP2 optimized structure. shorter when modeling water with TIP4P/2005 compared to either the MP2 level or using the TIP5P-Ew model. This effect seems to be the result of different partial charges on the water hydrogen for the two models, which are also reported in Table SII of the supplementary material.84 In TIP4P/2005 the hydrogen charge is +0.5564 e, which is larger than that in TIP5P-Ew (+0.241e), favoring stronger Coulombic interaction with the negative nitrogen of MA (see Figure 1). All water models considered in this work (TIP3P, TIP4P/2005 and TIP5P-Ew) have similar Lennard-Jones parameters for the van der Walls interactions (see Table SII in the supplementary material84) and, therefore, the difference in the HOH–NH2 interaction will be mainly determined by the partial charge assigned to the hydrogen atoms. Interestingly, the partial charge on the hydrogen atoms in TIP3P has an intermediate value (+0.4170e) between those of either TIP4P/2005 or TIP5P-Ew, which correlates with the order of the depths in the surface free energy minima in Figure 4. Thus, the deeper minima at the interface observed in TIP4P/2005 water slab seem to correspond to a more favorable interfacial HOH–NH2 interaction, resulting in a longer lifetime for the h-bond at the interface as shown in Figure 7(a). The differences in the partial charge distribution between the two water models affect the hydrogen bond network around MA. Indeed, the HOH–NH2 bond has a longer lifetime with TIP4P/2005 model than with TIP5P-Ew especially at the interface (Figures 7(a) and 7(b)) but the HNH–OH2 lifetime is longer in the case of TIP5P-Ew water than TIP4P/2005 (Figures 7(c) and 7(d)). In TIP4P/2005 the negative charge (−1.1128 e) is placed on a negative virtual side along the molecule bisector, weakening the electrostatic interaction between the hydrogen of the primary amine group and the neutral oxygen of the water molecule. On the other hand, in the TIP5P-Ew model the negative charge is assigned to two virtual sites (−0.214 e) mimicking the lone pairs of the water molecule, favoring a more balanced lifetime for the nitrogen donor and acceptor h-bonds. 0D 1D 2D 3D 0A 1A 2A 3A TIP4P/2005 Bulk 0.010 0.131 0.047 0.001 0.018 0.323 0.136 0.004 0.009 0.198 0.110 0.004 0.000 0.005 0.002 0.000 0A 1A 2A 3A TIP5P-Ew Bulk 0.009 0.034 0.012 0.001 0.050 0.207 0.091 0.006 0.059 0.314 0.196 0.016 0.001 0.004 0.001 0.000 0A 1A 2A 3A TIP4P/2005 Interface 0.006 0.324 0.042 0.000 0.005 0.419 0.086 0.000 0.001 0.090 0.026 0.001 0.000 0.001 0.001 0.000 0A 1A 2A 3A TIP5P-Ew Interface 0.037 0.118 0.013 0.000 0.105 0.365 0.069 0.001 0.043 0.187 0.059 0.002 0.000 0.000 0.000 0.000 All these features are reflected in the number of h-bonds between water and the primary amine group reported in Table IV. In this table, we reported the probability of finding h-bonds from water toward the nitrogen of the primary amine group, XA, and from nitrogen to water oxygen, XD. In the TIP4P/2005 bulk water, the most likely configuration is the primary amine group accepting and donating one hydrogen bond (1A-1D entry with probability 0.323) even if 1A-2D configurations are still possible (probability 0.198). On the contrary, the lone pairs in TIP5P-Ew model favor the configurations in which the primary amine group donating two h-bonds (1A-2D, probability 0.314) in the bulk. At the interface, the configurations with nitrogen donating only one h-bond become the most populated in both TIP4P/2005 and TIP5P-Ew water slab. Figure 9 reports some typical configurations with different number of water hydrogen bond donors and acceptors toward the primary amine group. It is interesting to connect this figure with the optimized geometries for MA solvated in small water clusters (up to 5 water molecules) obtained by recent electronic structure calculations of Sha-Sha et al.64 According to Sha-Sha et al. the most stable geometries were tetramer and pentamer water clusters with MA located outside the water ring and the nitrogen atom both donating and accepting one hydrogen bond. Since the 1D-1A configuration is the most common configuration at the interface, this observation provides a further molecular explanation for the surface preference of MA. The solvation environment around MA can be further investigated by performing a bivariate distribution for two angular variables describing the orientation of water around the solute, following the approach suggested in Ref. 65. For FIG. 9. Snapshots of 8 typical h-bonding configurations taken from 50 ns bulk simulation of MA in TIP5P-Ew water. XA indicates the number of h-bonds with nitrogen acting as acceptor while XD the number of h-bonds with nitrogen acting as donor. each water we defined a local frame of reference with the Z-axis point along the water dipole and the Y-axis along the vector connecting the water hydrogen. The X-axis is perpendicular to the water plane, such that the X, Y, Z axes form a right-hand frame of reference. θ µ is the angle between the water dipole and the vector connecting water oxygen and the nitrogen atom of the MA, ⃗r O w N while φ is the angle formed by the projection of ⃗r O w N on XY-plane and the X-axis. These angular variables are used to describe the water arrangement around the solute and to obtain insights about the hydrogen bond network between solute and water.65–68 As an initial reference, Figure 10 shows the bivariate distribution of an H2O around a water molecule over the course of 50 ns NPT simulation, recording the (θ µ , φ) angular position for the vector connecting the water oxygen and the oxygen of H2O within a distance of 0.35 nm from the water’s oxygen. The pattern shows two peaks, one about (θ µ , φ) (50◦, 90◦) and another at (θ µ , φ) (130◦, 0◦). The former dense area corresponds to water molecules donating h-bonds to the H2O while the second area to water accepting hbonds from H2O. Figures 10(b) and 10(d) show, in 3D, the regions explored by the oxygen of the selected H2O in the local frame of water molecules. Figures 10(a) and 10(c) show equally intense peaks for the two water molecules at (θ µ , φ) = (50◦, 90◦) while the lobes at φ = 0◦ look remarkably different, with a more spread populated area in TIP4P/2005 case. This asymmetry in the h-bond pattern is also visible in Figures 10(b) and 10(d) with a more equally distributed population of H2O around the water hydrogen and the lone pairs of TIP5P-Ew. This result shares some similarities with the findings of Remsing et al.,66 which show how the lone pairs in the TIP5P-water model families produce a symmetric h-bond pattern around hydrophobic cavity, resembling the results of accurate DFT-based MD in bulk water. Similarly but around a hydrophilic H2O, we observe an asymmetric h-bond pattern for TIP4P/2005 and a symmetric pattern for TIP5P-Ew in Figure 10(b). FIG. 10. Bivariate distribution of the (θ µ, φ) variables for a H2O around TIP4P/2005 (panel (a)) and TIP5P-Ew (panel (c)) bulk water. On the right side, the 3D location of the selected H2O oxygen around TIP4P/2005 water, panel (b), and TIP5P-Ew, panel (d). In panels (b) and (d), the brown and blue spheres qualitatively represent the position of the virtual site and lone pairs in TIP4P/2005 and TIP5P-Ew, respectively. FIG. 11. Bivariate distribution of the (θ µ, φ) variables for the nitrogen atom of the primary amine group around TIP4P/2005 (panel (a)) and TIP5P-Ew (panel (c)) in the bulk water simulation. On the right side, the 3D location of the nitrogen atom oxygen around TIP4P/2005 water, panel (b), and TIP5P-Ew, panel (d). In panels (b) and (d), the brown and blue spheres qualitatively represent the position of the virtual site and lone pairs in TIP4P/2005 and TIP5P-Ew, respectively. Both Figures 11 and 12 show the same type of analysis with the MA located in the bulk and at the interface, respectively, based on the angular position of ⃗r OW N . For both figures, we selected water molecules that were within 0.35 nm from the nitrogen atom, since this distance corresponds to the first solvation shell radius. In both solvation environments, the most populated region corresponds to (θ µ , φ) ∼ (50◦, 90◦), which corresponds to water molecules donating hydrogen bonds to the nitrogen atom of the primary amine group, i.e., HOH–NH2. These highly populated areas are located at the exact same position in TIP4P/2005 and TIP5P-Ew at about θ µ ∼ 50◦, consistently with the same hydrogen structural arrangement in the topology of the two water models. Looking at Table IV, we may also infer how the HOH–NH2 interaction mainly consists of the amine group accepting one h-bond from the water oxygen (1A lines in Table IV for both bulk and interface runs). On the other hand, there is a significant difference in the less populated maximum at φ ∼ 0◦, which would correspond to HNH–OH2 interaction. Indeed, the lobe associated to this interaction looks smeared out with TIP4P/2005 water and centered around θ µ ∼ 100◦, while in TIP5P-Ew is centered at θ mu ∼ 130◦. This pattern resembles the lack of lone pairs in the TIP4P/2005: the hydrogen of the primary amine group can point toward the water oxygen without any preferential orientation in TIP4P/2005 water while the tetrahedral charge arrangement in TIP5P-Ew induces a preferential direction for the electrostatic interaction with the positive charged hydrogen of MA. Figures 11 and 12 point toward a weaker HNH–OH2 interaction which may explain the MA preference for the interfacial environments seen in Figure 5. Interestingly, the same feature has been observed in accurate (yet extremely expensive and short simulation time) Born Oppenheimer molecular dynamics of the NH2 radical on a small liquid droplet,68 also showing that NH2 is preferentially located at the interface. Moreover, we observed the primary amine group bounded to a small number of water molecules, especially at the interface. This small coordination number is also consistent with recently proposed electronic structure calculations for hydrated MA in small water clusters64 and on similar studies on surface-bound state for ammonia.69 Nevertheless, the solute-solvent interaction appears to not be the only factor favoring the surface preference. Carignano et al.70 have shown that, in the case of ammonia, there is also an energetic cost to rearrange the solvent molecules around the solute. The solvent-solvent energy contribution could be equally important, especially in the case of MA as the hydrophobic cavity around the methyl group can induce an energy penalty for rearrangement of the solvent molecules while MA is in the bulk, favoring the migration of MA to the interface. This picture is consistent with the methyl group protruding from the interface while the primary amine group remains solvated, as shown in Figure 6. In conclusion, both solute-solvent and solvent-solvent energetic contributions are relevant for the preferential behavior of MA toward the interface. Nevertheless, this work reveals a need for new methods for better describing the solvent FIG. 12. Bivariate distribution of the (θ µ, φ) variables for the nitrogen atom of the primary amine group around TIP4P/2005 (panel (a)) and TIP5P-Ew (panel (c)), constraining the nitrogen position at the interfacial region. On the right side, the 3D location of the nitrogen atom oxygen around TIP4P/2005 water, panel (b), and TIP5P-Ew, panel (d). In panels (b) and (d), the brown and blue spheres qualitatively represent the position of the virtual site and lone pairs in TIP4P/2005 and TIP5P-Ew, respectively. around solute molecules than that presented here using the (θ µ , φ) bivariate distribution. This set of angular variables was designed to study simple solutes (water,65,66 NH2 radical,68 Cl−,67 etc.) in which the water distribution around the solute can be simply described using only two angular coordinates. The presence of both a hydrophilic group attached to a hydrophobic group breaks the symmetry, therefore requiring additional (or different) variables to properly describe the solvation around the MA. The hydrogen bond analysis proposed here suggests that the h-bond pattern could be somewhat dependent on the choice of force field, even if either force field combination used here (i.e., one using TIP4P/2005 and the other TIP5P-Ew to describe the liquid phase) predicts a free energy of hydration for MA close to experimental values. The free energy of hydration is a popular and often used benchmarking target for force fields, but our results show that the details regarding the solvation structure around a solute could vary with the model especially at the interface, even if the free energy of hydration is correctly reproduced. For this reason, proper care should be considered. Nevertheless, the results suggest that regardless of the force field combination used, the mass accommodation coefficients suggest that MA accommodates on the surface once it collides with it. To the best of our knowledge, experimental values for both α s and α b associated with MA on liquid droplets are not available, but our results are consistent with other near unity mass accommodation coefficients for other organic compounds.44 F. Chemistry implications of MA at the interface The fate of any MA that is emitted into the atmosphere from various anthropogenic sources has been suggested to be determined by either gas-phase oxidation or heterogeneous reactions with aerosols.71,72 The partitioning of MA into cloud droplets involves several steps: gas phase interaction with the water interface of cloud droplets, transport across the vapor-liquid water interface, hydrolysis in the aqueous phase, aqueous phase diffusion, and chemical reaction inside the cloud droplet. Methylamine belongs to a class of organics referred as amphiphilic organics, which are known to be partitioned at the interface creating a hydrophobic film on an aqueous aerosol73,74 and determining many chemical properties of aerosols. Recent experimental and theoretical research75–77 has shown that atmospheric reactions of organic molecules can be enhanced on a surface compared to the gas phase. As will any surface active molecule, amphiphilic organics will reduce the surface tension of water, affecting the droplet growth and trace-gas uptake, finally determining the reactivity towards oxidative gases and the ability of the aerosol to absorb or scatter radiation.78,79 Understanding the orientation of methylamine on the surface of water gives important insight into how methylamines at the gas-water interface may react with gasphase trace gases in the atmosphere. The molecular dynamics simulations studies from this work show that a property of methylamine is that it contains two regions of directly opposing polarities: a hydrophobic tail and hydrophilic head group. Results from this work clearly state that the MA uptake on the water surface is close to unitary and MA accommodates at the interface with the hydrophilic head group, i.e., NH2 group, solvated by superficial water molecules, and the hydrophobic tail, the CH3 group, pointing into the air phase (as shown in Figure 6). Specifically, the primary amine group acts as hydrogen bond donor and acceptor with interfacial water, but the first solvation shell barely interacts with the methyl group. The hydrogen bond analysis reported in Figure 7 clearly shows that the hydrogen bonds between water and the primary amine group are longer living at the interface than in the bulk, making unlikely the reaction of trace gases with the solvated amine group. This leaves the methyl group prone to react with ozone and OH radicals in the gas phase. The results regarding the accommodation and orientation of MA at the interface of liquid water point to unique chemistry pathways involving MA on the surface of liquid water aerosols. The gas phase reaction of OH with methylamine has been suggested to proceed via hydrogen abstraction from either the CH or a NH bond.80 The two abstraction paths have been modeled and branching ratios of the channels of H abstraction from the CH3 and NH2 groups show that in the gas phase the NH site accounts for about 20% of the overall reaction. However, this work shows that in the case FIG. 13. Suggested chemistry pathway for methylamine at the surface of liquid water aerosols. Through reactions with trace atmospheric gases the methyl group is converted in a carbonyl site or gives rise to the formation of solvated NH2 radical and formaldehyde. of methylamine at the cloud water-interface, the CH3 group is exposed to attack by gas-phase molecules since the NH2 is solvated and strongly hydrogen bonded by superficial waters. Indeed, OH gas phase radicals need to be taken, diffuse, and break the hydrogen bond network around NH2 in order to interact with it. On the other hand and as it has been already suggested by Vaida et al.,81 the CH3 group is immediately exposed to the attack of gas OH radical followed by the subsequent reaction of the newly formed radical centers with atmospheric O2, NO, and NO2, which leads to the formation of additional new carbonyl sites, as shown in Figure 13. Once formed, these new carbonyl sites functionalize the surface and provide hydrophilic sites for uptake of additional water. This chemistry can explain the tendency of aged organic aerosols to form CCNs and nucleate clouds.14 Moreover, these new carbonyl sites can accommodate additional methylamine molecules via the facile water-catalyzed addition reaction described in the recent work of Louie et al.,82 promoting aerosol growth. Given that just a single water molecule is sufficient to make this reaction cycle energetically feasible and facile,82 this growth mechanism can occur under the conditions experienced by atmospheric aerosol. Nevertheless, Figure 13 points also to the possibility of another channel with the formation of NH2 radical and formaldehyde, which are also known to be preferentially located at the surface of liquid water.68,77 It is known from studies83 of substituted alkoxy radicals that C–X bonds are weaker than C–H ones because of the lower activation energy barrier for bond cleavage, suggesting that the release of NH2 radical is also possible. Interestingly, this proposed mechanism for the chemistry of MA at the air/liquid water interface also suggests a new source for NH2 radicals at aerosol surfaces, other than by reaction of absorbed NH3. This provides important support for the recent work of Martin-Costa et al.77 and Zhong et al.68 on the interaction of NH2 radical and volatile organic compounds at the surface of water droplet. Future work is planned to examine the influence of both temperature and humidity on the position/orientation of such by-products at aerosol droplets. IV. CONCLUSIONS MA is an amine compound present in the atmosphere, which has been recently indicated as a catalyzer in the formation of SOAs and cloud condensation nuclei. Methylamine is amphiphilic and, thus, it can partition at the interface between the gas and liquid phase, affecting water and gases uptake on atmospheric aerosols, the heterogeneous reactivity of oxidative gases, and the ability of aerosol to absorb or scatter radiation. Novel results about the accommodation and orientation of MA on the air/liquid water interface from molecular dynamics simulations provide new insight into the reactivity of MA on water interfaces. An almost unity surface mass accommodation, α s , was determined irrespectively by the force field choice. Important new findings from this work show that the primary amine group of MA results in being a better hydrogen bond acceptor than donor toward water, especially at the interface: this and the energy cost required to rearrange water around the hydrophobic methyl group in the bulk could be at the origin of MA preference for the interface. The preferential location of MA at the interface with the methyl group protruding into the gas environment was, herein, verified by molecular dynamics simulations. These findings may have important atmospheric implications: the methyl group exposed to the gas phase can enhance the formation of new carbonyl sites by the attack of gas-phase OH followed by the subsequent reaction with atmospheric O2 and NO. These new carbonyl sites functionalize the surface and provide more hydrophilic sites for uptake of additional water and trace gases, which could explain the tendency of aged organic aerosols to be good cloud condensation nuclei. Moreover, a new channel for the reaction is the formation of NH2 radical and formaldehyde on the interface of liquid water aerosols. This channel becomes a new source of NH2 on aerosol interfaces, in addition to that resulting from NH3. The results about uptake and accommodation of MA can be used to draw similar and more general conclusions for the chemistry of other amphiphilic organics, amines in particular, at the surface of atmospherically relevant aerosols. ACKNOWLEDGMENTS The authors would like to thank Professor Sergey Nizkorodov for his thoughtful reading and feedback on this manuscript, which greatly improved it. 1X. Ge, A. S. Wexler, and S. L. Clegg, “Atmospheric amines. Part I. A review,” Atmos. Environ. 45, 524–546 (2011). 2K. Sellegri, B. Umann, M. Hanke, and F. Arnold, “Deployment of a ground- based cims apparatus for the detection of organic gases in the boreal forest during the quest campaign,” Atmos. Chem. Phys. 5, 357–372 (2005). 3X. Ge, A. S. Wexler, and S. L. Clegg, “Atmospheric amines. Part II. Thermodynamic properties and gas/particle partitioning,” Atmos. Environ. 45, 561–577 (2011). 4C. Qiu and R. Zhang, “Multiphase chemistry of atmospheric amines,” Phys. Chem. Chem. Phys. 15, 5738–5752 (2013). 5D. Lee and A. S. Wexler, “Atmospheric amines. Part III. Photochemistry and toxicity,” Atmos. Environ. 71, 95–103 (2013). 6A. Laskin, J. Laskin, and S. A. Nizkorodov, “Chemistry of atmospheric brown carbon,” Chem. Rev. 115, 4335–4382 (2015). 7V. Loukonen, T. Kurten, I. K. Ortega, H. Vehkamaki, A. Padua, K. Sellegri, and M. Kulmala, “Enhancing effect of dimethylamine in sulfuric acid nucleation in the presence of water-a computational study,” Atmos. Chem. Phys. 10, 4961–4974 (2010). 8M. Hallquist, J. C. Wenger, U. Baltensperger, Y. Rudich, D. Simpson, M. Claeys, J. Dommen, N. M. Donahue, C. George, A. H. Goldstein, J. F. Hamilton, H. Herrmann, T. Hoffmann, Y. Linuma, M. Jang, M. E. Jenkin, J. L. Jimenez, A. Kiendler-Scharr, W. Maenhaut, G. McFiggans, Th. F. Mentel, A. Monod, A. S. H. Prevot, J. H. Seinfeld, J. D. Surratt, R. Szmigielski, and J. Wildt, “The formation, properties and impact of secondary organic aerosol: Current and emerging issues,” Atmos. Chem. Phys. 9, 5155–5236 (2009). 9C. Leng, J. D. Kish, J. E. Roberts, I. Dwebi, N. Chon, and Y. Liu, “Temperature-dependent Henry’s law constants of atmospheric amines,” J. Phys. Chem. A 119, 8884–8891 (2015). 10G. W. Schade and P. J. Crutzen, “Emission of aliphatic amines from animal husbandry and their reactions: Potential source of N2O and HCN,” J. Atmos. Chem. 22, 319–346 (1995). 11H. Chen, M. E. Varner, R. B. Gerber, and B. J. Finlayson-Pitts, “Reactions of methanesulfonic acid with amines and ammonia as a source of new particles in air,” J. Phys. Chem. B 120, 1526 (2016). 12M. L. Dawson, M. E. Varner, V. Perraud, M. J. Ezell, B. R. Gerber, and B. J. Finlayson-Pitts, “Simplified mechanism for new particle formation from methanesulfonic acid, amines, and water via experiments and ab initio calculations,” Proc. Natl. Acad. Sci. U. S. A. 109, 18719–18724 (2012). 13J. Kua, H. E. Krizner, and D. O. D. Haan, “Thermodynamics and kinetics of imidazole formation from glyoxal, methylamine, and formaldehyde: A computational study,” J. Phys. Chem. A 115, 1667–1675 (2011). 14C. George, M. Ammann, B. D’Anna, D. J. Donaldson, and S. A. Nizkorodov, “Heterogeneous photochemistry in the atmosphere,” Chem. Rev. 115, 4218–4258 (2015). 15D. Frenkel and B. Smit, Understanding Molecular Dynamics (Academic Press, 2002). 16R. Vacha, P. Jungwirth, J. Chen, and K. T. Valsaraj, “Adsorption of polycyclic aromatic hydrocarbons at the air–water interface: Molecular dynamics simulations and experimental atmospheric observations,” Phys. Chem. Chem. Phys. 8, 4461 (2006). 17R. Vacha, L. Cwiklik, J. Rezac, P. Hobza, P. Jungwirth, K. Valsaraj, S. Bahr, and V. Kempter, “Adsorption of aromatic hydrocarbons and ozone at environmental aqueous surfaces,” J. Phys. Chem. A 112, 4942 (2008). 18D. Heger, D. Nachtigallova, F. Surman, J. Krausko, B. Magyarova, M. Brumovsky, M. Rubes, I. Gladich, and P. Klan, “Self-organization of 1-methylnaphthalene on the surface of artificial snow grains: A combined experimental-computational approach,” J. Phys. Chem. A 115, 11412 (2011). 19T. P. Liyana-Arachchi, K. T. Valsaraj, and F. R. Hung, “Molecular simulation study of the adsorption of naphthalene and ozone on atmospheric air/ice interfaces,” J. Phys. Chem. A 115, 9226 (2011). 20T. P. Liyana-Arachchi, K. T. Valsaraj, and F. R. Hung, “Adsorption of naphthalene and ozone on atmospheric air/ice interfaces coated with surfactants: A molecular simulation study,” J. Phys. Chem. A 116, 2519 (2012). 21R. Kania, J. K. Malongwe, D. Nachtigallova, J. Krausko, I. Gladich, M. Roeselova, D. Heger, and P. Klan, “Spectroscopic properties of benzene at the air–ice interface: A combined experimental–computational approach,” J. Phys. Chem. A 118, 7535–7547 (2014). 22J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, “Development and testing of a general amber force field,” J. Comput. Chem. 25, 1157 (2004). 23C. I. Bayly, P. Cieplak, W. D. Cornell, and P. A. Kollman, “A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: The RESP model,” J. Phys. Chem. 97, 10269 (1993). 24M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox,  09, Revision E.01, Gaussian Inc., Wallingford, CT, 2009. 25J. Wang, W. Wang, P. A. Kollman, and D. A. Case, “Automatic atom type and bond type perception in molecular mechanical calculations,” J. Mol. Graphics Modell. 25, 247–260 (2006). 26W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, “Comparison of simple potential functions for simulating liquid water,” J. Chem. Phys. 79, 926–935 (1983). 27J. L. F. Abascal and C. Vega, “A general purpose model for the condensed phases of water: TIP4P/2005,” J. Chem. Phys. 123, 234505 (2005). 28S. W. A. Rick, “Reoptimization of the five-site water potential (TIP5P) for use with Ewald sums,” J. Chem. Phys. 120, 6085 (2004). 29I. Gladich, A. Rodriguez, R. P. H. Enriquez, F. Guida, F. Berti, and A. Laio, “Designing high-affinity peptides for organic molecules by explicit solvent molecular dynamics,” J. Phys. Chem. B 119, 12963–12969 (2015). 30J. L. F. Abascal and C. Vega, “Dipole-quadrupole force ratios determine the ability of potential models to describe the phase diagram of water,” Phys. Rev. Lett. 98, 237801 (2007). 31I. Gladich and M. Roeselova, “Comparison of selected polarizable and nonpolarizable water models in molecular dynamics simulations of ice Ih,” Phys. Chem. Chem. Phys. 14, 11371–11385 (2012). 32I. Gladich, A. Oswald, N. Bowens, S. Naatz, P. Rowe, M. Roeselova, and S. Neshyba, “Mechanism of anisotropic surface self-diffusivity at the prismatic ice-vapor interface,” Phys. Chem. Chem. Phys. 17, 22947–22958 (2015). 33C. Vega, E. Sanz, and J. L. F. Abascal, “The melting temperature of the most common models of water,” J. Chem. Phys. 122, 114507 (2005). 34R. G. Fernandez, J. L. F. Abascal, and C. Vega, “The melting point of ice Ih for common water models calculated from direct coexistence of the solidliquid interface,” J. Chem. Phys. 124, 144506 (2006). 35G. C. Picasso, D. Malaspina, M. A. Carignano, and I. Szleifer, “Cooperative dynamic and diffusion behavior above and below the dynamical crossover of supercooled water,” J. Chem. Phys. 139, 044509 (2013). 36D. J. Tobias, A. C. Stern, M. D. Baer, Y. Levin, and C. J. Mundy, “Simulation and theory of ions at atmospherically relevant aqueous liquid-air interfaces,” Annu. Rev. Phys. Chem. 64, 339–359 (2013). 37I. Gladich, P. Shepson, I. Szleifer, and M. Carignano, “Halide and sodium ion parameters for modeling aqueous solutions in TIP5P-Ew water,” Chem. Phys. Lett. 489, 113–117 (2010). 38I. Gladich, A. Habartova, and M. Roeselova, “Adsorption, mobility, and self-association of naphthalene and 1-methylnaphthalene at the water-vapor interface,” J. Phys. Chem. A 118, 1052–1066 (2014). 39D. A. McQuarrie and J. S. Simon, Physical Chemistry: A Molecular Approach (University Science Books, Herndon, VA, 1997). 40P. Atkins and J. DePaula, Physical Chemistry (Oxford University Press, UK, 2010). 41R. Sander, “Compilation of Henry’s law constants (version 4.0) for water as solvent,” Atmos. Chem. Phys. 15, 4399–4981 (2015). 42M. Roeselova, J. Vieceli, L. X. Dang, B. C. Garrett, and D. J. Tobias, “Hydroxyl radical at the airwater interface,” J. Am. Chem. Soc. 126, 16308–16309 (2004). 43J. Vieceli, M. Roeselova, N. Potter, L. X. Dang, B. C. Garrett, and D. J. Tobias, “Molecular dynamics simulations of atmospheric oxidants at the airwater interface: Solvation and accommodation of OH and O3,” J. Phys. Chem. B 109, 15876–15892 (2005). 44J. Julin, P. M. Winkler, N. M. Donahue, P. E. Wagner, and I. Riipinen, “Nearunity mass accommodation coefficient of organic molecules of varying structure,” Environ. Sci. Technol. 48, 12083–12089 (2014). 45J. A. Lemkul and D. R. Bevan, “Assessing the stability of Alzheimer’s amyloid protofibrils using molecular dynamics,” J. Phys. Chem. B 114, 1652–1660 (2010). 46J. Hub, B. L. de Groot, and D. van Der Spoel, “g-wham—A free weighted histogram analysis implementation including robust error and autocorrelation estimates,” J. Chem. Theory Comput. 6, 3713 (2010). 47O. Engin, A. Villa, M. Sayar, and B. Hess, “Driving forces for adsorption of amphiphilic peptides to the airwater interface,” J. Phys. Chem. B 114, 11093 (2010). 48S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, and P. A. Kollman, “The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method,” J. Comput. Chem. 13, 1011–1021 (1992). 49S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, and P. A. Kollman, “Multidimensional free-energy calculations using the weighted histogram analysis method,” J. Comput. Chem. 16, 1339–1350 (1995). 50B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, “GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation,” J. Chem. Theory Comput. 4, 435–447 (2008). 51G. Bussi, D. Donadio, and M. Parrinello, “Canonical sampling through velocity rescaling,” J. Chem. Phys. 126, 014101 (2007). 52H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola, and J. R. Haak, “Molecular-dynamics with coupling to an external bath,” J. Chem. Phys. 81, 3684–3690 (1984). 53U. Essmann, L. Perera, M. L. Berkowitz, T. D. H. Lee, and L. G. Pedersen, “A smooth particle mesh Ewald method,” J. Chem. Phys. 103, 8577–8593 (1995). 54B. Hess, H. Bekker, H. J. C. Berendsen, and J. Fraaije, “LINCS: A linear constraint solver for molecular simulations,” J. Comput. Chem. 18, 1463–1472 (1997). 55S. Miyamoto and P. A. Kollman, “Settle: An analytical version of the shake and rattle algorithm for rigid water models,” J. Comput. Chem. 13, 952–962 (1992). 56D. J. Donaldson and K. T. Valsaraj, “Adsorption and reaction of trace gasphase organic compounds on atmospheric water film surfaces: A critical review,” Environ. Sci. Technol. 44, 865 (2010). 57C. Vega, J. L. F. Abascal, M. M. Conde, and J. L. Aragones, “What ice can teach us about water interactions: A critical comparison of the performance of different water models,” Faraday Discuss. 141, 251 (2009). 58C. Vega and J. L. F. Abascal, “Simulating water with rigid non-polarizable models: A general perspective,” Phys. Chem. Chem. Phys. 13, 19663–19688 (2011). 59M. Shiraiwa, A. Zuend, A. Bertram, and J. H. Seinfeld, “Gas-particle partitioning of atmospheric aerosols: Interplay of physical state, non-ideal mixing and morphology,” Phys. Chem. Chem. Phys. 15, 11441–11453 (2013). 60R. Kumar, M. C. Barth, S. Madronich, M. Naja, G. R. Carmichael, G. G. Pfister, C. Knote, G. P. Brasseur, N. Ojha, and T. Sarangi, “Effects of dust aerosols on tropospheric chemistry during a typical pre-monsoon season dust storm in northern India,” Atmos. Chem. Phys. 14, 6813–6834 (2014). 61B. Ervens, “Modeling the processing of aerosol and trace gases in clouds and fogs,” Chem. Rev. 115, 4157–4198 (2015). 62A. Luzar and D. Chandler, “Effect of environment on hydrogen bond dynamics in liquid water,” Phys. Rev. Lett. 76, 928–931 (1996). 63O. Markovitch and N. Agmon, “Reversible geminate recombination of hydrogen-bonded water molecule pair,” J. Chem. Phys. 129, 084505 (2008). 64S.-S. Lv, Y.-R. Liu, T. Huang, Y.-J. Feng, S. Jiang, and W. Huang, “Stability of hydrated methylamine: Structural characteristics and H2N · · · H–O hydrogen bonds,” J. Phys. Chem. A 119, 3770–3779 (2015). 65P. Jedlovszky, M. Predota, and I. Nezbeda, “Hydration of apolar solutes of varying size: A systematic study,” Mol. Phys. 104, 2465–2476 (2006). 66R. C. Remsing, M. D. Baer, G. K. Schenter, C. J. Mundy, and J. D. Weeks, “The role of broken symmetry in solvation of a spherical cavity in classical and quantum water models,” J. Phys. Chem. Lett. 5, 2767–2774 (2014). 67M. D. Baer, D. J. Tobias, and C. J. Mundy, “Investigation of interfacial and bulk dissociation of HBr, HCl, and HNO3 using density functional theory-based molecular dynamics simulations,” J. Phys. Chem. C 118, 29412–29420 (2014). 68J. Zhong, Y. Zhao, L. Li, H. Li, J. S. Francisco, and X. C. Zeng, “Interaction of the NH2 radical with the surface of a water droplet,” J. Am. Chem. Soc. 137, 12070–12078 (2015). 69D. J. Donaldson, “Adsorption of atmospheric gases at the airwater interface. I. NH3,” J. Phys. Chem. A 103, 62–70 (1999). 70M. A. Carignano, M. M. Jacob, and E. E. Avila, “On the uptake of ammonia by the water-vapor interface,” J. Phys. Chem. A 112, 3676–3679 (2008). 71C. J. Nielsen, H. Herrmann, and C. Weller, “Atmospheric chemistry and environmental impact of the use of amines in carbon capture and storage (CCS),” Chem. Soc. Rev. 41, 6684–6704 (2012). 72S. M. Murphy, A. Sorooshian, J. H. Kroll, N. L. Ng, P. Chhabra, C. Tong, J. Surratt, E. Knipping, R. C. Flagan, and J. H. Seinfeld, “Secondary aerosol formation from atmospheric reactions of aliphatic amines,” Atmos. Chem. Phys. 7, 2313–2337 (2007). 73G. B. Ellison, A. F. Tuck, and V. Vaida, “Atmospheric processing of organic aerosols,” J. Geophys. Res.: Atmos. 104, 11633–11641, doi:10.1029/1999JD900073 (1999). 74D. J. Donaldson and V. Vaida, “The influence of organic films at the airaqueous boundary on atmospheric processes,” Chem. Rev. 106, 1445–1461 (2006). 75Y. Paz, S. Trakhtenberg, and R. Naaman, “Reaction between O(3P) and organized organic thin films,” J. Phys. Chem. 98, 13517–13523 (1994). 76A. K. Bertram, A. V. Ivanov, M. Hunter, L. T. Molina, and M. J. Molina, “The reaction probability of OH on organic surfaces of tropospheric interest,” J. Phys. Chem. A 105, 9415–9421 (2001). 77M. T. C. Martins-Costa, J. M. Anglada, J. S. Francisco, and M. F. RuizLopez, “Reactivity of volatile organic compounds at the surface of a water droplet,” J. Am. Chem. Soc. 134, 11821–11827 (2012). 78B. J. Finlayson-Pitts and J. N. Pitts, Jr., Chemistry of the Upper and Lower Atmosphere: Theory, Experiments, and Applications (Academic Press, 1999). 79J. H. Seinfeld and S. N. Pandis, Atmospheric Chemistry and Physics: From Air Pollution to Climate Change (John Wiley & Sons, 2012). 80A. Galano and J. R. Alvarez-Idaboy, “Branching ratios of aliphatic amines + OH gas-phase reactions: A variational transition-state theory study,” J. Chem. Theory Comput. 4, 322–327 (2008). 81V. Vaida, A. Tuck, and G. Ellison, “Optical and chemical properties of atmospheric organic aerosols,” Phys. Chem. Earth 25, 195–198 (2000). 82M. K. Louie, J. S. Francisco, M. Verdicchio, S. J. Klippenstein, and A. Sinha, “Dimethylamine addition to formaldehyde catalyzed by a single water molecule: A facile route for atmospheric carbinolamine formation and potential promoter of aerosol growth,” J. Phys. Chem. A 120, 1358 (2016). 83A. C. Davis and J. S. Francisco, “Reactivity trends within alkoxy radical reactions responsible for chain branching,” J. Am. Chem. Soc. 133, 18208–18219 (2011). 84See supplementary material at http://dx.doi.org/10.1063/1.4950951 for description of mass accommodation coefficient and water model parameters.