Articles | Volume 23, issue 10
Research article
01 Jun 2023
Research article |  | 01 Jun 2023

Collision-sticking rates of acid–base clusters in the gas phase determined from atomistic simulation and a novel analytical interacting hard-sphere model

Huan Yang, Ivo Neefjes, Valtteri Tikkanen, Jakub Kubečka, Theo Kurtén, Hanna Vehkamäki, and Bernhard Reischl

Kinetics of collision-sticking processes between vapor molecules and clusters of low-volatility compounds govern the initial steps of atmospheric new particle formation. Conventional non-interacting hard-sphere models underestimate the collision rate by neglecting long-range attractive forces, and the commonly adopted assumption that every collision leads to the formation of a stable cluster (unit mass accommodation coefficient) is questionable for small clusters, especially at elevated temperatures. Here, we present a generally applicable analytical interacting hard-sphere model for evaluating collision rates between molecules and clusters, accounting for long-range attractive forces. In the model, the collision cross section is calculated based on an effective molecule–cluster potential, derived using Hamaker's approach. Applied to collisions of sulfuric acid or dimethylamine with neutral bisulfate–dimethylammonium clusters composed of 1–32 dimers, our new model predicts collision rates 2–3 times higher than the non-interacting model for small clusters, while decaying asymptotically to the non-interacting limit as cluster size increases, in excellent agreement with a collision-rate-theory atomistic molecular dynamics simulation approach. Additionally, we calculated sticking rates and mass accommodation coefficients (MACs) using atomistic molecular dynamics collision simulations. For sulfuric acid, a MAC ≈1 is observed for collisions with all cluster sizes at temperatures between 200 and 400 K. For dimethylamine, we find that MACs decrease with increasing temperature and decreasing cluster size. At low temperatures, the MAC ≈1 assumption is generally valid, but at elevated temperatures MACs can drop below 0.2 for small clusters.

1 Introduction

In the first steps of atmospheric new particle formation (NPF), condensed-phase clusters (Vehkamäki and Riipinen2012) form from super-saturated vapors (Guo et al.2014; Kirkby et al.2011; Kulmala et al.2000, 2013; Zhang et al.2012). Precise quantification of the cluster formation rates is important to accurately assess global aerosol budgets and their impact on climate (Hallquist et al.2009; Ramanathan et al.2001; Seinfeld et al.2016) and health (Pope and Dockery2006; Pöschl2005). Cluster formation is conventionally approximated as a microscopic kinetic process where the size evolution of the new-born clusters is driven by competing events of monomer acquisition and loss. This results in a set of pseudo-first-order kinetic equations (Schenter et al.1999) (written here for a one-component system for simplicity):

(1) d N i d t = α i - 1 , 1 k i - 1 , 1 N i - 1 N 1 - α i , 1 k i , 1 N i N 1 - γ i N i + γ i + 1 N i + 1 , i = 1 , 2 , ,

where Ni is the number concentration of clusters composed of i monomers (atoms or molecules), ki,1 the monomer–cluster collision rate coefficient, γi the evaporation rate coefficient, and αi,1 the mass accommodation coefficient which is a measure of the average probability for the monomer and cluster to “stick together” after a collision. The collision rate coefficient and mass accommodation coefficient can be combined into the sticking rate coefficient si,1=αi,1ki,1. In principle, any cluster formation process can be accurately described by Eq. (1) and appropriate extensions (e.g., considering cluster dissociation or coagulation) if the associated rate coefficients are correctly determined. This work focuses particularly on the collision and sticking rate coefficients.

Precise determination of the collision and sticking rate coefficients is non-trivial, as they are sensitive to multiple variables including cluster size, structure, composition, system temperature, and vapor concentration. Conventional models, like the commonly used classical nucleation theory (CNT; Becker and Döring1935; Farkas1927), resort to crude assumptions to achieve generality. This inevitably sacrifices precision and leads to predictions agreeing with experiments only within a narrow range of saturation ratios and temperatures. In most conventional approaches, sticking rate coefficients are approximated by non-interacting hard-sphere collision rate coefficients derived from kinetic gas theory. By neglecting intermolecular interactions, kinetic gas theory may not predict accurate collision rate coefficients, as intermolecular long-range attractive forces can enhance collisions, especially for polar or polarizable molecules common in atmospheric aerosol formation. Moreover, using collision rate coefficients to approximate sticking rate coefficients may be questionable, especially for small clusters, as the formation of a stable product cluster after a collision may be unsuccessful in a significant number of events; i.e., the mass accommodation coefficient can be lower than unity. The non-interacting hard-sphere approximation thus potentially leads to significant errors in modeling the initial steps of cluster formation.

Molecular dynamics (MD) simulations allow us to study the time evolution of collision systems, where all intra- and intermolecular interactions are described through force fields. Recently, we have developed several atomistic molecular dynamics simulation frameworks (Yang et al.2018; Halonen et al.2019; Neefjes et al.2022) to determine collision probabilities from sampling trajectories of atoms or molecules, which predict enhanced collision rate coefficients in good agreement with experiments (Lehtipalo et al.2016; Stolzenburg et al.2020). However, as the adequate sampling of collision probabilities requires simulating a large number of binary collisions, these approaches are computationally expensive. More efficient ways need to be proposed to systematically calculate rate coefficients, in particular for collisions involving molecules and clusters of different sizes.

Rate coefficients can also be obtained through analytical models that consider the interactions between the collision partners. These models usually require a fraction of the computational cost of MD simulations but generally still rely on significant approximations. In the central field approach, the collision partners are approximated by point particles interacting through an isotropic attractive potential. Along this line, Gioumousis and Stevenson (1958) developed a model to calculate the secondary ion–molecule reaction rates in mass spectrometers. Su and Bowers (1973) derived a theory to describe ion–permanent dipole interactions and combined it with the central field approach to model ion–polar molecule collisions. Several other authors (Moran and Hamill1963; Su et al.1978; Clary1985) have also theoretically studied collision rate coefficients between ions and molecules by applying the central field approach.

We have previously shown that the central field approach yields very similar results to those from sampling collision trajectories using MD, provided that the attractive interaction used in the central field model is fitted to the tail of the potential of mean force (PMF) between the collision partners (Neefjes et al.2022). The largest collision partners considered in that study were, however, dimers. For collisions involving larger clusters, it becomes increasingly difficult to obtain accurate PMFs. Additionally, in the central field model, a collision is only considered to have occurred when the center of mass distance between the collision partners is infinitesimal. Such a collision criterion may underestimate cluster formation rates, as in realistic systems, cluster formation can already occur at larger center-to-center distances, where the extended structures of the collision partner begin to overlap.

In this study, we present two extensions to the central field approach to efficiently predict collision rate coefficients between monomers and clusters of arbitrary sizes.

  1. We consider an interacting hard-sphere model, in which the collision partners are still treated as point particles, interacting through an effective long-ranged isotropic attractive potential, but a collision is considered to have occurred if the center of mass distance between the collision partners is less than the sum of their hard-sphere radii to account for their extended structures.

  2. We integrate the monomer–monomer attractive potential following the approach of Hamaker (Hamaker1937) to obtain an effective monomer–cluster attractive potential for which a well-defined collision cross section (CCS) can be determined and the monomer–cluster collision rate coefficients can be calculated. Parameters of the monomer–monomer attractive potential are taken from the PMF between the collision partners, calculated from well-tempered metadynamics simulations using atomistic force fields. To validate our analytical approach, we compared it to both atomistic MD collision simulations and MD simulations between two point particles interacting through the same effective potential.

Furthermore, we provide arguments from MD collision simulations for when a unit mass accommodation coefficient is valid and when this approximation fails. Sticking probabilities and sticking rate coefficients inferred from MD collision simulations can be used in Eq. (1) to model acid–base-induced atmospheric NPF and compare it to experimental results.

As test systems, we considered collisions between the monomers of sulfuric acid (H2SO4) and dimethylamine ((CH3)2NH) and clusters of dimethylammonium–bisulfate dimers ([HSO4- ⋅ (CH3)2NH2+]n, n=1,2,4,8,16,32). Sulfuric acid is known to be a crucial component of continental NPF, while dimethylamine likely facilitates the formation of clusters through acid–base interactions (Kurtén et al.2008; Sipilä et al.2010). This paper is structured as follows: in Sect. 2.1, we present our analytical interacting hard-sphere collision rate model and show how the CCSs can be calculated for monomer–monomer and monomer–cluster collisions using a simple van der Waals potential and an effective potential derived from Hamaker's approach. In Sect. 2.2, we describe the well-tempered metadynamics calculations and the atomistic and point particle MD simulations. In Sect. 3, we validate the analytical model and present and discuss simulation results of the collision and sticking probabilities and rate coefficients for sulfuric acid–dimethylamine clusters. Section 4 summarizes and concludes the paper with an emphasis on its implications to atmospheric NPF modeling.

2 Theoretical and computational methods

2.1 Analytical methods

2.1.1 Central field approach

In the central field model, both collision partners are represented by point particles. We consider the initial condition where one particle remains at rest while the other particle approaches from infinitely far away with some initial relative velocity v0. The perpendicular distance between the initial velocity of the moving particle and the center of the particle at rest defines the impact parameter b. The geometry of an example system is illustrated in Fig. 1. This typical two-dimensional two-body motion can be first reduced to a two-dimensional single-body orbital motion in a central field and then to a one-dimensional single-body motion along the particle center-to-center direction (Landau and Lifshitz1976), i.e.,

(2) 1 2 μ v 0 2 = U ( r ) + 1 2 μ v 2 = U eff ( r ) + 1 2 μ v r 2 ,

where Ueff(r)=U(r)+μvt2/2=U(r)+μv02b2/2r2 is an effective potential resulting from reducing the two-dimensional orbital motion to the one-dimensional motion along the center-to-center direction, U(r) the attractive potential between the point particles, r the center-to-center distance, v the instantaneous relative velocity, vr the component of the instantaneous relative velocity along the center-to-center direction, vt the component of the instantaneous relative velocity perpendicular to the center-to-center direction, and μ the reduced mass.

Figure 1Central field model. The collision partners are described by point particles. In the frame of reference of the target particle, the collision geometry is fully described by the incoming particle's initial velocity v0 and the impact parameter b.


Reducing the two-dimensional single-body orbital motion to a one-dimensional motion along the center-to-center direction introduces a centrifugal energy barrier (the repulsive part of the effective potential) between the point particles. In the central field model, a “collision” is defined by zero distance between the point particles. Crossing the centrifugal energy barrier in the center-to-center direction is, therefore, a necessary and sufficient condition for a collision. This condition is equivalent to the kinetic energy in the center-to-center direction 1/2μvr2 being positive for all values of r>0; i.e., for a collision to happen, the following condition must hold for arbitrary r>0:

(3) U eff ( r ) - 1 2 μ v 0 2 0 .

After substituting in the expression of the effective potential and rearranging the above equation, we obtain

(4) b 2 r 2 1 - 2 U ( r ) μ v 0 2 ω v ( r ) ,

where the function ωv(r) is defined for convenience. As Eq. (4) must hold for arbitrary r>0, it is apparent that the minimum value of function ωv(r) in interval r>0 sets an upper boundary for the impact parameter b, below which a collision will happen and above which a collision is not possible. This upper boundary defines the critical impact parameter bc,CF in the central field model:

(5) b c , CF 2 = ω v ( R m ) ,

where ωv(Rm) is the minimum of function ωv(r) with Rm being the corresponding distance where the minimum is achieved. Physically, Rm is the closest distance that the collision partners can approach each other when the centrifugal barrier cannot be crossed.

Depending on the attractive potential U(r), Eq. (5) might, however, not provide a satisfactory definition for the critical impact parameter. Figure 2a shows the shape of ωv(r) for general attractive potentials of the form U(r)=-A(r/r0)a. For this general form, we can distinguish three cases.

  • If a<-2, ωv(r) is a convex function exhibiting a single minimum at r=Rm=-r0[A(a+2)/μv02)]-1/a>0 and ωv0,a<-2(Rm)=-ar02[A(a+2)/μv02)]-2/a/(a+2).

  • If a=-2, ωv(r) monotonically increases with a minimum at r=Rm=0 and ωv,a=-2(Rm)=2Ar02/μv02.

  • If -2<a<0, ωv(r) is a monotonically increasing function with a minimum at r=Rm=0 and ωv,-2<a<0(Rm)=0.

These cases indicate that (1) for -2<a<0, such as for the Coulombic potential, bc,CF2=ωv,-2<a<0(Rm)=0 and the CCS is equally zero; the centrifugal barrier can never be crossed, and a collision is impossible regardless of the initial velocity and impact parameter. (2) For a<-2 and large initial velocities v0, bc,CF2=ωv,a<-2(Rm)=-ar02[A(a+2)/μv02)]-2/a/(a+2), and the CCS likewise becomes vanishingly small. These issues arise because the collision partners are approximated as point particles in the central field approach; a collision is only considered to have occurred when the distance between the particles is zero. In realistic systems, the collision partners have an extended structure. As a result, the point particle approximation is expected to give a lower collision rate than when the extent of the system is considered.

Figure 2Shape of ωv(r). (a) The general shape of the function ωv(r) as a function of r≥0 depending on the interaction exponential a. (b) The general shape of ωv(r) for a<-2, together with two possible values for the critical impact parameter bc. In the case Ri+Rj<Rm, the critical impact parameter is determined solely by the interaction potential. However, if Ri+RjRm, the critical impact parameter is determined by the interaction potential and sum of the hard-sphere radii.


2.1.2 Interacting hard-sphere model

In the interacting hard-sphere model, the aforementioned situations are accounted for by using a revised collision criterion; i.e., if the distance between the collision partners is at any point smaller than the sum of their hard-sphere radii, Ri+Rj, a collision is considered to have occurred. It is apparent that when Ri+Rj<Rm, this revised criterion will give an identical result as that in the central field model because in this case “reaching a distance smaller than Ri+Rj” is equivalent to “crossing the centrifugal barrier” for the point particles. The new criterion will make a difference only when Ri+RjRm, where even if the centrifugal barrier is not crossed, there is still a chance for the point particles to reach a distance smaller than Ri+Rj (see Fig. 2b). Using this revised criterion, we can define two cases for the critical impact parameter bc:

(6) b c 2 = ω v ( R m ) , if R m > R i + R j , ω v ( R i + R j ) , if R m R i + R j .

The concept of determining the critical impact parameter by examining the aforementioned two cases has been discussed by other authors (Ouyang et al.2012; Fuchs and Sutugin1965), but explicit expressions using the hard-sphere radii, like Eq. (6), have not been presented in the existing literature to the best of our knowledge. Once the appropriate critical impact parameter is determined, the dynamical collision cross section (CCS) is given by the usual expression,

(7) Ω ( v 0 ) = 2 π 0 b c ( v 0 ) b d b = π b c 2 ( v 0 ) ,

and the collision rate coefficient for interacting hard spheres is obtained as

(8) k IHS = 0 Ω ( v 0 ) v 0 f ( v 0 ) d v 0 ,

where f(v0) is the Maxwell–Boltzmann distribution. In the case of non-interacting hard spheres, the CCS is independent of the relative velocity, and the collision rate coefficient can be directly calculated by substituting Ω(v)=π(Ri+Rj)2 into the above equation:

(9) k HS = 8 k B T π μ 1 / 2 π ( R i + R j ) 2 ,

where kB is the Boltzmann constant and T is the temperature.

To determine the critical impact parameter bc in the interacting hard-sphere model, we first need to find the value of Rm for a given interaction potential U(r), which will be done in the following sections.

Solution for the monomer–monomer potential

Vapors relevant to atmospheric cluster formation are mostly composed of polar and polarizable molecules interacting through a van der Waals potential (Leite et al.2012) which scales with the separation distance as

(10) U mm ( r ) = - 4 ϵ r σ - 6 ,

where ϵ and σ are the characteristic energy and length of the attractive potential and the subscript “mm” denotes the monomer–monomer interaction. Substituting Eq. (10) into Eq. (4), the positive minimum of ωv(r,Umm) is found to be

(11) ω v ( R m , U mm ) = 3 σ 2 2 ϵ μ v 0 2 1 / 3 ,

located at

(12) R m ( U mm ) = σ 16 ϵ μ v 0 2 1 / 6 .

For a given pair of monomers, bc can be readily evaluated from Eqs. (6) and (12) with respect to any given relative speed v0. We obtained values for ϵ and σ from the potential of mean force (PMF) along the center of mass distance between the collision partners calculated from well-tempered metadynamics calculations using an empirical atomistic force field, presented in Sect. 2.2 (results are included in the Supplement).

Solution for the monomer–cluster potential

Following the approach of Hamaker (Hamaker1937), an effective monomer–cluster potential Umc(r) can be obtained by integrating the monomer–monomer potential over the volume of a spherical cluster, in which uniform monomer number density is assumed:


where ρc is the monomer number density, Rc the radius of the cluster, dV=ρ2sin (ϕ)dθdϕdρ the volume element, Vc the total volume of the cluster, and nc=ρcVc the total number of monomers in the cluster. In Eq. (13), the effective monomer–cluster potential is essentially a sum over all individual contributions of the monomer–monomer potentials, which has two implications. First, the monomer–monomer potentials must be pairwise so that their individual contributions are additive. Hence, this approach is not suitable for non-pairwise interactions. Second, the approach can be conveniently extended to multi-component clusters containing several types of monomers. By using Eq. (13) in a similar manner, the corresponding multi-component monomer–cluster potential is -4inc,iϵiσi6/r2-Rc3, where nc,i is the number of monomers of type i in the cluster and ϵi and σi are the interaction parameters between a monomer of type i in the cluster and the free monomer colliding with the cluster. Note that further integrating Umc(r) over the volume of another cluster recovers the well-known result of Hamaker for a van der Waals attractive potential between two spherical objects, which can be invoked to treat cluster–cluster collision rates. Integrated cluster–cluster potentials have previously been used to model collision rates, but the models usually resort to fitting against experimental data (Fuchs and Sutugin1965; Sceats1989). Replacing U(r) in Eq. (4) with Umc(r), we find ωv(r,Umc) to be a convex function with a positive minimum (note that we require r>Rc). By taking the derivative of ωv(r,Umc) with respect to r, the minimum is found to be located at r=Rm(Umc) and satisfies the following:

(14) i = 0 4 a i R m 2 i = 0 ,

with a0=Rc2(Rc6-lc6), a1=-2lc6-4Rc6, a2=6Rc4, a3=-4Rc2, and a4=1, where lc[8ncϵσ6/(μv02)]1/6 is a cluster size- and potential-dependent characteristic length. Equation (14) is a fourth-order (quartic) equation after writing it in terms of Rm2, and it can be shown that only one real root exists for Rm>Rc (see Supplement). The solution is readily available:

(15) R m 2 ( U mc ) = R c 2 + M + - M 2 - q 4 M ,

where q=-2lc6, M=(N+Δ0/N)/3/2, and N=(Δ1+Δ12-4Δ03)/23 with Δ0=-36Rc2lc6 and Δ1=108lc12. With Eqs. (6) and (15), we can predict bc for any monomer and cluster pair, as long as the monomer–monomer interaction parameters ϵ and σ are provided.

2.2 Computational methods

Computer simulations were used to

  1. compute the PMF between monomer pairs to obtain the ϵ and σ parameter values required for the analytical collision rate coefficient model;

  2. validate the analytical model by comparison with atomistic and point particle molecular dynamics (MD) collision simulations;

  3. calculate predictions for the sticking probability, sticking rate coefficient, and mass accommodation coefficient by analyzing the atomistic MD collision simulations over a range of relative velocities and impact parameters.

2.2.1 Systems and force field

We considered the atmospherically relevant sulfuric acid (H2SO4) and dimethylamine ((CH3)2NH) molecules, as well as clusters consisting of 1, 2, 4, 8, 16, or 32 bisulfate–dimethylammonium dimers ([HSO4- ⋅ (CH3)2NH2+]). Clusters composed of n [HSO4- ⋅ (CH3)2NH2+] dimers were obtained by sintering and equilibrating two smaller clusters with n/2 dimers.

In the PMF and atomistic MD simulations, these systems are described by a force field fitted according to the OPLS (optimized potentials for liquid simulations) all-atom procedure (Jorgensen et al.1996). In the OPLS force field, the intramolecular interactions consist of harmonic bond potentials between covalently bonded atoms, harmonic angle potentials between atoms separated by two covalent bonds, and dihedral angle potentials between atoms separated by three covalent bonds,

(16) U intra OPLS = i = 1 N bonds k i b 2 r i - r i 0 2 + j = 1 N angles k j θ 2 θ j - θ j 0 2 + k = 1 N dihedrals n = 1 4 V n 2 1 + cos n ϕ k - ϕ n k ,

where kib, ri, and ri0 are the force constant, instantaneous, and equilibrium length of bond i; kjθ, θj, and θj0 are the force constant, instantaneous, and equilibrium value of angle j; and Vn, ϕnk, and ϕk are the Fourier coefficients, phase angles, and instantaneous value of the dihedral angle k.

The intermolecular interactions, as well as intramolecular interactions between atoms separated by more than three covalent bonds, are described by Lennard-Jones potentials between atoms i and j separated by a distance rij, with distance and energy parameters σij and ϵij, as well as Coulomb interactions between the atoms' partial charges qi and qj:

(17) U inter = i = 1 N 1 j = 1 N 2 4 ϵ i j σ i j r i j 12 - σ i j r i j 6 + i = 1 N 1 j = 1 N 2 1 4 π ϵ 0 q i q j r i j ,

where ϵ0 is the vacuum permittivity.

The OPLS force field parameters used in this study were obtained from Loukonen et al. (2010). We note that in the original OPLS force field, Lennard-Jones and Coulomb interactions between atoms separated by three covalent bonds (“1–4 interactions”) are scaled by a factor 0.5. Loukonen et al. (2010) set this scaling factor to zero when fitting the force field parameters. For consistency, we have also set these interactions to zero in our simulations.

We have validated these force field parameters for MD collision simulations in our previous work (Halonen et al.2019) by comparing predicted molecular structures, binding energies, and vibrational spectra of sulfuric acid molecules with ab initio results.

2.2.2 Potential of mean force calculations

To determine the thermally averaged interaction potential between the collision partners, we calculated the PMF as a function of their center of mass distance from well-tempered meta-dynamics simulations (Barducci et al.2008), using the PLUMED plug-in (Tribello et al.2014) for LAMMPS (Plimpton1995). Calculations were performed for five combinations of monomer pairs at 300 K: (CH3)2NH(CH3)2NH, H2SO4H2SO4, H2SO4(CH3)2NH, H2SO4–[HSO4- ⋅ (CH3)2NH2+]1, and (CH3)2NH–[HSO4- ⋅ (CH3)2NH2+]1. We used the velocity Verlet integrator with a time step of 1 fs, where the Lennard-Jones interactions were cut off at 14 Å, while electrostatic interactions were evaluated with a cut-off at 40 Å. A total of 40 random walkers were employed to deposit Gaussians with a width of 0.1 Å and initial height of 2kBT (kBT/10 for the (CH3)2NH(CH3)2NH system) every 500 time steps along the collective variable; a harmonic wall was used to restrict the collective variable to values below 35 Å. A bias factor of 20 (5 for the (CH3)2NH(CH3)2NH system) was chosen, and a stochastic velocity rescaling thermostat with a time constant of 0.1 ps was used to maintain a constant temperature. Note that we used the PMF at 300 K to approximate cases at 200 and 400 K in deriving corresponding collision rate coefficients. Although there are slight differences in the PMFs at these temperatures, their influence on the analytical collision rate coefficients was found to be negligible. Results are summarized in Fig. S1 and Table S1 of the Supplement. The attractive interaction between (CH3)2NH and (CH3)2NH was too small to be distinguished from thermal noise in our calculations. It was therefore directly set to zero and is not shown in Fig. S1. The monomer–monomer interactions were assumed to take the form of a 12–6 Lennard-Jones potential. The depth of the potential well (or binding energy) ϵ was directly taken as the minimum energy of the PMF curve, and σ was determined from σ=rϵ/26, where rϵ is the corresponding distance for the minimum energy point on the PMF curve. The binding energy for the H2SO4H2SO4 system was found to be −0.29 eV, which is in excellent agreement with previous ab initio results (Temelso et al.2012) of −0.3 eV obtained from Boltzmann averaging over four energy minimum structure dimers at 298.15 K and is in good agreement with recent calculations (Elm et al.2016; Myllys et al.2017) of −0.23 and −0.26 eV at a higher level of theory. As the outcome of the collisions is already decided at a significant separation between the collision partners for the systems considered here, the role of proton transfer after a successful collision and other quantum effects can be ignored, and the classical force field approach is expected to capture the most important features of the collision.

We note that, though in the current paper the monomer–monomer interacting parameters were obtained from the PMF calculated by atomistic simulations, they could in principle be taken directly from literature values if available.

2.2.3 Atomistic molecular dynamics collision simulations

Atomistic simulations were performed with the LAMMPS code (Plimpton1995) for collisions of either H2SO4 or (CH3)2NH with [HSO4- ⋅ (CH3)2NH2+]n, where n=1,2,4,8,16, or 32 at temperatures of 200, 300, and 400 K. Equations of motion were integrated using a velocity Verlet algorithm with a time step of 1 fs, where Lennard-Jones and electrostatic interactions were evaluated with a cut-off at 280 Å. The collision partners are originally placed 550 Å apart along the x axis, well beyond the cut-off of the intermolecular potentials. At these positions, the collision partners were equilibrated for 50 ps using separate Langevin thermostats for each collision partner with a damping factor of 0.1 ps. The atomic velocities were drawn from the Maxwell–Boltzmann distribution at the target temperature with the center of mass motion of each collision partner and angular momentum of the combined system removed. After equilibration, the distance between the now orientationally randomized collision partners was decreased to 150 Å along the x axis, bringing them within range of the long-range intermolecular potentials. Both collision partners were then given a velocity along the x direction of vx=±v0/2 towards each other, where v0 is the initial relative velocity. For each system, the range of relative velocities started at 50 m s−1 and increased in steps of 50 m s−1. The highest relative velocity was determined so that at least 99 % of the Maxwell–Boltzmann distribution was sampled. We sampled impact parameters starting from 0 up to 36 Å, in steps of 1 Å along the z axis. At initial relative velocities of 50 m s−1, there is still some probability for a collision even beyond an impact parameter of 36 Å. These probabilities, however, do not contribute significantly to the collision rate coefficient.

For the collision systems with clusters containing 1, 2, 4, 8, 16, and 32 [HSO4- ⋅ (CH3)2NH2+] dimers, we performed N=1000, 1000, 800, 400, 200, and 100 independent trajectory simulations, respectively, for each combination of relative velocity and impact parameter. These result in statistically significant estimations of the collision probability. The resulting trajectories were analyzed to determine the collision probability: if the collision partners’ center of mass distance was smaller than Ri+Rj (determined once again from the mass and bulk density by assuming spherical shapes) in at least one time frame, then a collision event was identified. This definition is consistent with that in the analytical model, and it leads to the collision probability Pc(b,v0)=Nc/N, where Nc is the numbers of identified collision events. The definition of the sticking probability is similar to that of the collision probability, i.e., Ps(b,v0)=Ns/N, where the criterion for determining Ns is discussed in detail in Sect. 3.3.1.

2.2.4 Point particle molecular dynamics collision simulations

In the point particle molecular dynamics (PP MD) simulations, collision partners were treated as point masses interacting through effective pair potentials. The purpose of performing such simulations was to verify the critical impact parameters derived from the analytical interacting hard-sphere model. The collision simulation setup in this case is similar to that in the atomistic MD simulation. First, the point particles were placed 200 Å apart on the x axis, with a certain impact parameter along the y axis, and given an initial relative velocity along the x direction. Second, the equations of motion of the point particles were solved with the velocity Verlet algorithm. A dynamic time step interval was used to ensure that the displacement of the point particles at each time step cannot exceed 0.1 Å. A collision simulation was stopped when either of the following two situations happened: (1) the distance between the point particle pair is less that the sum of their hard-sphere radii (given in the analytical interacting hard-sphere model), or (2) the distance between the point particle pair is larger than the starting distance. Situation 1 corresponds to a collision event, while situation 2 indicates that a collision will not happen in that simulation. The critical impact parameter for a given initial relative speed was found by gradually increasing the impact parameter (starting from the summation of hard-sphere radii given in the interacting hard-sphere model) until the transition from situation 1 to situation 2 was observed. This procedure was used to simulate two example cases where point particle pairs interact through either a Coulomb or a Lennard-Jones potential. Simulated critical impact parameters as a function of initial relative speeds and their comparisons with those predicted from the analytical interacting hard-sphere model are discussed in Sect. 3.1.

2.2.5 Collision rate, sticking rate, and mass accommodation coefficients

The MD collision rate coefficient kMD was calculated by numerically integrating the collision probability as a function of impact parameter and relative initial velocity obtained from the MD trajectory simulations over the impact parameter and relative speed distribution at the target temperature:

(18) k MD = 2 π 0 d v 0 0 v 0 f ( v 0 ) b P c ( b , v 0 ) d b .

Note that Eq. (18) reduces to Eq. (8) by letting Pc(b,v0)=1 for b<bc and Pc(b,v0)=0 for b>bc. Similarly, the sticking rate coefficient is calculated through

(19) s MD = 2 π 0 d v 0 0 v 0 f ( v 0 ) b P s ( b , v 0 ) d b .

The mass accommodation coefficient is defined as

(20) α MD = s MD k MD ,

and it characterizes the average probability of sticking after collisions at the specified temperature.

3 Results and discussion

3.1 Validation of the analytical interacting hard-sphere model

To validate our analytical interacting hard-sphere model, we first compared the analytical critical impact parameter bc predicted by Eq. (6) to molecular dynamics simulations of two point particles interacting through the same effective pair potential (i.e., PP MD). Two cases were considered: (1) collisions between the (CH3)2NH2+ ion and HSO4- ion under the attractive Coulomb potential U=-e2/(4πϵ0r), where e is the elementary charge and ϵ0 the vacuum permittivity (Fig. 3a), and (2) collisions between the neutral sulfuric acid H2SO4 monomer and the neutral acid–base [HSO4- ⋅ (CH3)2NH2+]1 dimer under the van der Waals potential U=-4ϵ(σ/r)6, with parameters obtained from the potential of mean force (PMF) along the center of mass (COM) distance between the molecules (Fig. 3b). In the PP MD simulations, bc was found for a specific relative speed by gradually increasing b until the distance r between the particles is never smaller than the sum of the hard-sphere radii. In the analytical model, Rm=0 for the attractive Coulomb potential and Rm=σ16ϵ/(μv02)1/6 for the van der Waals potential.

Excellent agreement is observed in both cases (note that the dependent variable bc is set as the x axis in Fig. 3). For the Coulomb potential, Ri+Rj is always larger than Rm, and bc=ωv(Ri+Rj)1/2 in the analytical model. The Coulomb potential case discussed here is only used to verify Eq. (6) but should not be implemented further in Eqs. (7) and (8) to calculate collision rate coefficients, as the Coulomb potential is non-negligible at distances comparable to the mean free path of the colliding ions at atmospheric pressures and hence violates the assumption of free molecular regime. For the van der Waals potential, Rm>Ri+Rj at low relative speeds, but at relative velocities over  800 m s−1 the situation is reversed and Rm<Ri+Rj. For velocities above  800 m s−1, the analytical bc transitions from ωv(Rm)1/2 to ωv(Ri+Rj)1/2. It is important to note that for the H2SO4–[HSO4- ⋅ (CH3)2NH2+]1 system, relative speeds over  800 m s−1 account for less than 1 % of the Boltzmann relative speed distribution. In this case, bc=ωv(Rm)1/2 is a good approximation. However, as Rmσ(ϵ/μv02)1/6, smaller values of ϵ and σ or larger values of ω and v0 can cause this transition to occur at lower relative speeds and increase the importance of the bc=ωv(Ri+Rj)1/2 case in the calculation of collision rate coefficients. In the Supplement, we show examples of collisions of carbon dioxide and water molecules, where considering ωv(Ri+Rj) is important to obtain reasonable collision cross sections. These systems will be studied in further detail in a future publication.

Figure 3Critical impact parameter. Critical impact parameters bc from the analytical model and from point particle molecular dynamics simulations (PP MD) for (a) Coulomb attractive potential for the (CH3)2NH2+ and HSO4- ion pair and (b) van der Waals attractive potential for the H2SO4 and [HSO4- ⋅ (CH3)2NH2+]1 neutral pair. As the relative initial velocity v0 increases, the system in panel (b) transitions from a regime where the critical impact parameter bc=ωv(Rm) to bc=ωv(Ri+Rj). Example orbits for v0=2000 m s−1 and v0=150 m s−1 are shown in panels (c) and (d), respectively.


In Fig. 4, we compare the critical impact parameters obtained from the analytical model (Eq. 6) to the collision probability Pc(b,v0) obtained from atomistic collision MD simulations, sampling the relevant range of impact parameters and relative velocities. We examine both “monomer–monomer” (H2SO4 and [HSO4- ⋅ (CH3)2NH2+]1) and “monomer–cluster” (H2SO4 and [HSO4- ⋅ (CH3)2NH2+]16) collisions. For the former, the effective interaction in the analytical model is obtained from the PMF as a function of the monomer–monomer COM distance (see Fig. S1 and Table S1 in the Supplement). For the latter, the effective interaction is obtained from Hamaker's approach (Eq. 13). Note that the atomistic force field describing the molecules and ions in the collision MD simulations is the same as in the PMF calculation, from which we derive the effective interaction parameters used to calculate the analytical critical impact parameters.

In general, Pc(b,v0) exhibits a clear boundary between values of b and v0 where collisions are highly probable (white region) and values where the probability is essentially 0 (black region). The location of this boundary depends on the intermolecular forces between the collision partners. An exception to this clear boundary can be seen at low relative speeds. Here, Pc(b,v0)<1 even for small values of b, while at the same time, Pc(b,v0) remains non-zero even for large values of b. This issue arises from the asymmetric structure of the molecules: due to the low relative speed, small but periodic repulsive interactions at certain collision angles can result in the molecules rebounding from each other before the orbit crosses the sum of the hard-sphere radii (Halonen et al.2019). These low relative speeds, however, only account for a very small percentage of the Boltzmann speed distribution.

For both systems, the analytical bc values obtained from ωv(Rm) are in excellent agreement with the collision probabilities obtained from atomistic collision MD simulations over the relevant range of relative velocities at 300 K. The transition to Rm>Ri+Rj occurs for relative velocities around 800 ms−1 for H2SO4 and [HSO4- ⋅ (CH3)2NH2+]1 and around 2300 m s−1 for H2SO4 and [HSO4- ⋅ (CH3)2NH2+]16. Once again, this indicates that the bc=ωv(Ri+Rj)1/2 case is not relevant for these systems at the studied temperature of 300 K. Most importantly, the excellent agreement between analytical bc values and collision probabilities from atomistic collision MD simulations for the “monomer–cluster” collision demonstrates the validity of Hamaker's approach in modeling the effective monomer–cluster interactions of the studied systems.

Figure 4Collision probability and critical impact parameter. Collision probability Pc(b,v0) from atomistic collision MD simulation and analytical critical impact parameters bc (solid lines) for the collision of (a) H2SO4 and [HSO4- ⋅ (CH3)2NH2+]1 and (b) H2SO4 and [HSO4- ⋅ (CH3)2NH2+]16 at 300 K. For the analytical monomer–cluster critical impact parameter in panel (b), the effective monomer–cluster interaction potential was obtained using Eq. (13). The dotted line represents the sum of the hard-sphere radii of the collision partners, obtained by assuming bulk density and spherical shape.


3.2 Collision rate enhancement factors predicted by the analytical interacting hard-sphere model

Having validated the analytical interacting hard-sphere model by comparison with MD simulations, we now use the model to predict enhancement factors η over the non-interacting hard-sphere model. We define η=kIHS/kHS as the ratio between the collision rate coefficient obtained from the interacting and non-interacting hard-sphere model. Figure 5 shows η for H2SO4–[HSO4- ⋅ (CH3)2NH2+]n collisions and (CH3)2NH–[HSO4- ⋅ (CH3)2NH2+]n collisions at 300 K. Both monomer types have similar η: while the H2SO4 monomer has stronger intermolecular interactions, the (CH3)2NH monomer has larger thermal mean speeds. For both cases, the monomer–cluster collisions gradually approach non-interacting hard-sphere behavior as the size of the cluster increases, i.e., η→1 as n→∞. At smaller size ranges the non-interacting hard-sphere model can underestimate the collision rate coefficient by a factor of 2 to 3, potentially introducing non-negligible systematic errors in atmospheric NPF models. Hence, the correction due to long-range attractive forces may need to be included in NPF models.

Figure 5Collision rate enhancement factor. Monomer–cluster collision rate enhancement factors gradually decay to 1 with increasing cluster radii, indicating a transition to hard-sphere-like behavior. The main plot shows the data with logarithmic radius axis and the inset with linear radius axis.


3.3 Analysis of molecular dynamics collision simulations

3.3.1 Timescale analysis and issues related to classical force field models

In this work, we use classical, atomistic molecular dynamics simulations to determine parameters used in the analytical model and to validate the analytical model by comparing them to collision MD simulations. It is therefore important to recall the limitations of classical force fields regarding the physico-chemical processes occurring in gas-phase collisions and the timescales on which these different processes occur.

Whether or not a collision will occur is typically determined at a relatively large intermolecular distance, where the collision partners can accurately be modeled using classical force field descriptions in a molecular dynamics simulation. Modeling the probability of the collision partners forming a stable product using MD simulations is, however, significantly more complex, as the “sticking” probabilities are influenced by the possibility of bond formation between the collision partners, which is not captured by non-reactive classical force fields. The acid–base clusters considered in this study can form hydrogen bonds upon collision, which is captured by the classical force field employed in the atomistic simulation; however, proton transfer between acids and bases cannot happen.

After a collision, the newly formed cluster will generally have some excess energy due to these bond formations. In the atmosphere, this energy will be removed from the cluster through collisions with surrounding gas molecules. In the collision MD simulations, however, no additional carrier gas is modeled. As such, the excess energy after a collision cannot be released.

Lastly, it is likely that the cluster dissociation process is unphysically enhanced due to the fact that in the classical atomistic model employed, each vibrational mode possesses kBT/2 energy on average, while in a quantum-mechanical description of the same system, some high-frequency degrees of freedom would remain “frozen” at atmospherically relevant temperatures.

The mentioned limitations of MD simulations in modeling the sticking probability all enhance the dissociation of the formed cluster. As such, the sticking probability obtained from MD simulations can be seen as a limiting value. If the collision partners are able to stick together for a certain time after collision in the MD simulations, then it is likely that the same collision partners would be able to stick together for the same amount of time in the atmosphere. MD simulations are, however, limited in the timescale that they can model due to the computational cost. We, therefore, discuss some characteristic timescales in atmospheric cluster formation after collision.

We first differentiate the two types of monomer–cluster interactions influencing the cluster formation process: (1) the vibration and diffusion of a condensable vapor monomer on the cluster surface immediately after a collision, which, if successful, drives the cluster formation while storing excess energy due to bond formation (hydrogen bonds in this study), and (2) collisions between background carrier gas and the cluster, dissipating the excess energy and equilibrating/thermalizing the nascent cluster. We consider the following scenario: after a condensable vapor–cluster collision, the vapor monomer experiences a few rounds of unsteady vibrations on the cluster surface before forming stable hydrogen bonds (i.e., sticking) or rebounding from the cluster surface. In the former case, a stable cluster with excess thermal energy, distributed over its degrees of freedom, is produced. At this point, the captured monomer can only escape from the cluster through thermal fluctuation (i.e., evaporation), as the monomer has become thermally indistinguishable from other molecules in the cluster. At typical atmospheric conditions, the timescale ts for the formation of stable bonds only requires a few vibrations of the impinging monomer on the cluster surface. Hence, ts is commonly smaller than the mean free time tc for the carrier gas–cluster collision and the time teq required to fully thermalize the cluster. Moreover, ts, tc, and teq can all be assumed to be orders of magnitude smaller than the time tev for a monomer to evaporate (see the schematic diagram in Fig. 6).

For example, for H2SO4–[HSO4- ⋅ (CH3)2NH2+] and (CH3)2NH–[HSO4- ⋅ (CH3)2NH2+] surrounded by air at standard atmospheric pressure and 300 K, ts2Lc/v¯10 ps, estimated based on its molecular neighbor separation distance Lc and the mean thermal speed v¯. Based on the molecular collision frequency, tc∼102 ps. As thermalization requires tens to hundreds of collisions (Yang et al.2022), teq∼103 ps. For these systems, by far the shortest evaporation timescale tev∼106 ps is observed for the evaporation of a (CH3)2NH from the (CH3)2NH-[HSO4- ⋅ (CH3)2NH2+] cluster. The timescales of all other possible evaporation processes are significantly longer, tev1017-1022 ps (Ortega et al.2012). Therefore, for the systems considered here, even the shortest tev is significantly larger than ts, tc, and teq.

The timescale analysis strongly indicates that a stable cluster is very likely to form if the impinging condensable vapor monomer can survive a very short window of unsteady vibration right after a collision.

Figure 6Events of interest in atmospheric cluster growth. The schematic shows events of interest in the acid–base cluster growth process and the corresponding timescales at atmospherically relevant pressures and temperatures.


3.3.2 Collision and sticking criteria

To calculate collision and sticking probabilities in atomistic molecular dynamics simulations we need simple, but robust, criteria for when two collision partners have collided and when they are sticking together, as the large number of individual trajectories makes a manual inspection unfeasible. Collision and sticking criteria can be based on geometric and/or energetic properties of the system (Yang et al.2018; Halonen et al.2019; Neefjes et al.2022). Here, we considered a criterion based on the center of mass distance between the collision partners, as this does not require knowledge of the detailed bonding between monomers and clusters of different sizes and can be directly related to the analytical interacting hard-sphere model.

Figure 7Classification of MD collision trajectories. Trajectories can be classified as “flyby” (no collision), “rebounding collision” (no formation of a metastable complex), collisions leading to the formation of a metastable complex that will re-dissociate on a timescale <25 ps, and “sticking collisions” leading to a cluster that remains stable for longer than 25 ps. The dashed black line indicates the sum of hard-sphere radii of the collision partners.


In analyzing the MD collision simulations, we broadly encounter four different kinds of trajectories: flyby, rebounding collision, collision with metastable product, and collision with stable product. An example of each of these types of trajectories is shown in Fig. 7. A flyby is characterized by a clearly parabolic curve when plotting the center of mass distance between the collision partners as a function of time. The curves of all collision trajectories show a bend towards a faster decrease in center of mass distance over time due to attractive interactions between the collision partners. In rebounding collisions, the trajectory sharply reverses direction after reaching its lowest intermolecular distance. In the two other collision types, a cluster is formed and the trajectory oscillates around some intermolecular distance. The collision partners can either dissociate after some time (metastable) or stay together longer than the assumed timescale for cluster equilibration (stable). As the impact parameter decreases, it becomes more difficult to distinguish between a flyby and a collision.

A distance criterion that is too strict runs the risk of miscounting collisions with relatively large center of mass distances arising from certain relative orientations as flybys. Conversely, a distance criterion that is too loose results in certain close flybys being counted as rebounding collisions. Generally, a good distance criterion would lie just above the largest amplitudes of the oscillating sticking trajectories in Fig. 7 and would need to be found for each combination of T, v0, and b, which is not possible due to the large number of combinations that are analyzed. We found the sum of hard-sphere radii of the studied collision partners presented a surprisingly good distance criterion for the studied systems over a wide range of impact parameters and relative velocities. A comparison between the hard-sphere radii distance criterion and a manually optimized distance criterion yielded differences in collision probabilities of at most 1 %.

Based on our analyses, a collision is considered to have occurred if the center of mass distance between the collision partners is less than the sum of their hard-sphere radii in at least one time frame of the trajectory. Based on our timescale analyses in Sect. 3.3.1, a formed cluster is considered “sticking” if the center of mass distance between the collision partners is less than the sum of hard-sphere radii at any time between 25 ps after collision and the end of the simulation.

3.4 Collision and sticking rate coefficients from molecular dynamics simulations

Based on the above analysis, we monitored the center of mass distance between the impinging monomer and cluster during the atomistic MD simulations. Collision probabilities and sticking probabilities 25 ps or 50 ps after collision are shown in Fig. 8 for two different cases: (A) collisions between H2SO4 and [HSO4- ⋅ (CH3)2NH2+]1 and (B) collisions between (CH3)2NH and [HSO4- ⋅ (CH3)2NH2+]1. The figure also shows the 1000 underlying trajectories at the relative velocity and impact parameter closest to the mode of the Maxwell–Boltzmann distribution and the sum of hard-sphere radii, respectively, for both systems. There are several implications of the results in Fig. 8. First, the collision and sticking probabilities are almost identical for the system of H2SO4 and [HSO4- ⋅ (CH3)2NH2+]1, suggesting that most collisions lead to the formation of stable product clusters. For the system of (CH3)2NH and [HSO4- ⋅ (CH3)2NH2+]1, the sticking probability is visibly lower than the collision probability, indicating that a significant portion of formed clusters dissociate in the first 25 ps after collision. This can largely be attributed to the fact that the hydrogen bond between the (CH3)2NH monomer and [HSO4- ⋅ (CH3)2NH2+]n cluster is weaker than that between the H2SO4 monomer and [HSO4- ⋅ (CH3)2NH2+]n cluster. Second, even though some dissociations can still be observed in the underlying trajectories of the system of (CH3)2NH and [HSO4- ⋅ (CH3)2NH2+]1 later than 25 ps after collision, the sticking probability at 25 ps is almost identical to that at 50 ps for both systems, indicating that increasing ts beyond 25 ps has little effect on the sticking probability; 25 ps is sufficient for the formation of stable product clusters, which is consistent with our characteristic timescale analysis.

Figure 8Collision and sticking probability. Collision probability Pc and sticking probabilities Ps after 25 and 50 ps for (a) the H2SO4 and [HSO4- ⋅ (CH3)2NH2+]1 pair and (b) the (CH3)2NH and [HSO4- ⋅ (CH3)2NH2+]1 pair. The 1000 collision MD trajectories obtained for the most probable relative velocity and an impact parameter equal to the sum of hard-sphere radii are shown for both systems to illustrate the underlying dynamics of the collision. Dotted black lines represent the sum of collision pairs’ radii obtained by assuming the bulk density and spherical shape.


The sticking rate coefficients and rate enhancement factors for different clusters sizes and different temperatures are shown in Fig. 9a for H2SO4 and [HSO4- ⋅ (CH3)2NH2+]n collisions and Fig. 9b for (CH3)2NH and [HSO4- ⋅ (CH3)2NH2+]n collisions, where circles and squares correspond to results obtained at 25 and 50 ps, respectively. The sticking rate enhancement factor is an ad hoc quantity defined as the ratio between the sticking rate coefficient and the hard-sphere collision rate coefficient. For comparison, the collision rate coefficients from both the analytical model (solid line) and MD simulations (stars) are also shown in Fig. 9. For collisions between H2SO4 and [HSO4- ⋅ (CH3)2NH2+]n (Fig. 9a), the collision rate coefficients and the sticking rate coefficients obtained at the two different times are all identical, as the hydrogen bonds are strong enough to “capture” all impinging monomers colliding with the cluster. For collisions between (CH3)2NH and [HSO4- ⋅ (CH3)2NH2+]n (Fig. 9b), several interesting trends can be observed. At a low system temperature of 200 K, the collision rate coefficients are almost identical to the sticking rate coefficients except for a few data points at small cluster sizes. At this temperature, the binding energy of the hydrogen bond is large enough compared to the low kinetic energy of the impinging monomer to prevent rebounds. However, the deviation between the collision and sticking rate coefficients becomes visible with increasing system temperature, which leads to higher kinetic energies of the impinging monomer, resulting in more rebounding. At a system temperature of 300 K, this deviation is expected to be solely a consequence of the unsuccessful vibrational coupling events (i.e., rebounding), whereas at T=400 K, the characteristic evaporation time tev of the nascent product cluster has decreased to a value comparable to ts. At 400 K, the deviation between collision and sticking rate coefficients is caused by the coupled effects of the unsuccessful vibrational coupling events and the evaporation of the nascent product cluster. Consequently, visible differences between sticking rate coefficients obtained at 25 and 50 ps appear at 400 K. Moreover, the deviation between sticking and collision rate coefficients decreases with increasing cluster size: larger sizes favor the formation of stable clusters. This is not surprising as larger clusters offer more potential sites for hydrogen bonds and more molecular neighbors to enhance mutual attractions. They also serve as more efficient sinks for dissipating the latent heat released during hydrogen bond formation upon collision (Yang et al.2019, 2022).

Figure 9Rate coefficient and enhancement factor vs. cluster size (n) at various temperatures. Collision rate coefficient (CR), sticking rate coefficient (SR), hard-sphere collision rate coefficient (HSCR), collision rate enhancement factor (CEF = CR/HSCR), and sticking rate enhancement factor (SEF = SR/HSCR) for collisions of (a) H2SO4 and [HSO4- ⋅ (CH3)2NH2+]n and (b) (CH3)2NH and [HSO4- ⋅ (CH3)2NH2+]n at T=200, 300, and 400 K. The collision rate results are from both analytical predictions (Theo.) and molecular dynamics simulations (M.D.). The sticking rate results from molecular dynamics simulations were obtained with two different sticking timescale criteria, 25 and 50 ps.


It is important to remember the limitations of MD simulations using a classical force field description and the timescale analysis performed in Sect. 3.3.1 when predicting sticking rate coefficients. Overall, we expect our simulation approach to overestimate the probability of dissociation. Based on our simulations, it is therefore very likely that the H2SO4–[HSO4- ⋅ (CH3)2NH2+] system indeed has a mass accommodation coefficient near unity. For the (CH3)2NH–[HSO4- ⋅ (CH3)2NH2+] system, our simulations suggest a mass accommodation coefficient significantly less than unity at room temperature and above. Based on the weaker interactions between (CH3)2NH and [HSO4- ⋅ (CH3)2NH2+] this is plausible, but the limitations of our simulation approach do not permit us to make a definitive statement. Ideally, collision would be studied using ab initio or at least reactive MD simulations in the presence of carrier gas, and for much longer timescales, to properly model all post-collision processes. Such calculations are, however, beyond the scope of current computational possibilities.

4 Conclusions

Precise rate coefficients are essential in modeling the gas-phase collisions which are the initial step in atmospheric new particle formation (NPF). The non-interacting hard-sphere model commonly adopted for neutral molecules and clusters may significantly underestimate collision rate coefficients due to neglecting long-range attractive forces. On the other hand, the assumption of a unit mass accommodation coefficient can overestimate the number of stable clusters formed, in particular for small clusters at elevated temperatures. It is important to note, however, that these two errors will not necessarily cancel each other out.

Simulation approaches based on atomistic modeling of the collision partners have recently been developed and shown to give reasonable estimates for collision rates and enhancement factors for collisions of individual molecules or ions. To systematically calculate rate coefficients for collisions between molecules and clusters of atmospherically relevant acid–base systems, we have developed a new analytical interacting hard-sphere model. The molecule–cluster interactions are obtained from Hamaker's approach by integrating monomer–monomer interactions over the volume of the cluster. Here, the underlying monomer–monomer interacting parameters were obtained from fitting Lennard-Jones potential to the monomer–monomer potential of mean force calculated from atomistic simulations, but we note that the monomer–monomer interacting parameters could also be obtained by other methods or taken directly from literature values if available. The critical impact parameter is determined using the sum of hard-sphere radii or the minimum distance between the point-like collision partners which cannot be crossed due to the centrifugal barrier in the central field model. The accuracy of the analytical model was validated against molecular dynamics collision simulations using the full atomistic model of the collision partners. The analytical model has an accuracy comparable to atomistic molecular dynamics simulations but can be applied efficiently for the systematic calculation of molecule–cluster collision rates required in atmospheric NPF models. The same approach can also be applied to calculate effective cluster–cluster interactions.

For collisions of sulfuric acid (H2SO4) or dimethylamine ((CH3)2NH) molecules with clusters consisting of 1-32 bisulfate–dimethylammonium ([HSO4- ⋅ (CH3)2NH2+]) dimers, we find that for cluster radii smaller than 2 nm, the non-interacting hard-sphere model underestimates the collision rate coefficient by a factor of 2 to 3 due to the neglect of the long-range attractive forces. The enhancement factor drops below 1.2 for cluster radii larger than 25 nm. This deviation from the non-interacting hard-sphere model can be large enough to introduce systematic errors in atmospheric NPF models. For all the systems considered in this study, the critical impact parameter was well described using the standard minimum distance criterion of the central field model.

The analytical interacting hard-sphere model cannot give the probability of sticking or rebounding after collisions. To assess the fraction of collisions leading to a stable product cluster, we analyzed the large data set of collisions obtained during the validation of the analytical interacting hard-sphere model, taking into account the limitations of the classical force fields employed, as well as the timescales of equilibration processes in the atmosphere, not simulated explicitly in the binary collision setup.

Our analysis shows that if a monomer can stay on the cluster surface for more than a characteristic time (typically tens of picoseconds) after a collision, then a stable product cluster is formed. With this criterion, the sticking probability, sticking rate coefficient, and mass accommodation coefficient can be calculated by analyzing a set of molecular dynamics collision trajectories where impact parameters and relative speeds are properly sampled.

The mass accommodation coefficient decreases with the increasing temperature and decreasing cluster size. For collisions between H2SO4 and [HSO4- ⋅ (CH3)2NH2+]n, the unit mass accommodation coefficient assumption is generally valid in an atmospherically relevant temperature range of 200–400 K, regardless of the cluster size. However, the assumption may not hold if the temperature is further increased. For collisions between (CH3)2NH and [HSO4- ⋅ (CH3)2NH2+]n, the unit mass accommodation coefficient assumption is approximately valid at 200 K but becomes significantly worse at elevated temperatures; e.g., at 400 K it ranges from 0.1–0.6 for clusters with 1–32 monomers.

Future research directions include the investigation of systems where the transition between the two regimes defining the critical impact parameter occurs at lower and more relevant relative initial velocities, as well as the application of the analytical interacting hard-sphere model to cluster–cluster collisions. The model presented in this paper is expected to provide a useful tool for the atmospheric NPF community due to its relative simplicity, demonstrated good accuracy, and widespread applicability.

Data availability

The data generated in this study are fully represented in the figures and tables shown in the paper and the Supplement. Input files for simulations are available from the authors upon reasonable request.


The supplement related to this article is available online at:

Author contributions

HY conceived the research and derived the theoretical framework of the model. IN and HY performed the molecular dynamics collision simulations. IN provided theoretical background for the central field approach. BR and VT performed the metadynamics simulations. VT and JK helped to prepare the clusters used in the simulation. TK and HV helped to plan the project. HY, IN, and BR analyzed the simulation data and wrote the manuscript with comments from all other authors.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Computational resources were provided by the CSC–IT Center for Science Ltd., Finland. The authors wish to thank the Finnish Computing Competence Infrastructure (FCCI) for supporting this project with computational and data storage resources.

Financial support

This research has been supported by the European Research Council (project no. 692891 DAMOCLES), the Academy of Finland Flagship Program (grant no. 337549 ACCC) and Centers of Excellence Program (grant no. 346368 VILMA).

Open-access funding was provided by the Helsinki University Library.

Review statement

This paper was edited by Hinrich Grothe and reviewed by two anonymous referees.


Barducci, A., Bussi, G., and Parrinello, M.: Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method, Phys. Rev. Lett., 100, 020603,, 2008. a

Becker, R. and Döring, W.: Kinetische Behandlung der Keimbildung in übersättigten Dämpfen, Ann. Phys., 416, 719–752,, 1935. a

Clary, D.: Calculations of rate constants for ion-molecule reactions using a combined capture and centrifugal sudden approximation, Mol. Phys., 54, 605–618,, 1985. a

Elm, J., Jen, C. N., Kurtén, T., and Vehkamäki, H.: Strong Hydrogen Bonded Molecular Interactions between Atmospheric Diamines and Sulfuric Acid, J. Phys. Chem. A, 120, 3693–3700,, 2016. a

Farkas, L.: Keimbildungsgeschwindigkeit in übersättigten Dämpfen, Z. Phys. Chem., 125U, 236–242,, 1927. a

Fuchs, N. and Sutugin, A.: Coagulation rate of highly dispersed aerosols, J. Colloid Sci., 20, 492–500,, 1965. a, b

Gioumousis, G. and Stevenson, D. P.: Reactions of Gaseous Molecule Ions with Gaseous Molecules. V. Theory, J. Chem. Phys., 29, 294–299,, 1958. a

Guo, S., Hu, M., Zamora, M. L., Peng, J., Shang, D., Zheng, J., Du, Z., Wu, Z., Shao, M., Zeng, L., Molina, M. J., and Zhang, R.: Elucidating severe urban haze formation in China, P. Natl. Acad. Sci. USA, 111, 17373–17378,, 2014. a

Hallquist, M., Wenger, J. C., Baltensperger, U., Rudich, Y., Simpson, D., Claeys, M., Dommen, J., Donahue, N. M., George, C., Goldstein, A. H., Hamilton, J. F., Herrmann, H., Hoffmann, T., Iinuma, Y., Jang, M., Jenkin, M. E., Jimenez, J. L., Kiendler-Scharr, A., Maenhaut, W., McFiggans, G., Mentel, Th. F., Monod, A., Prévôt, A. S. H., Seinfeld, J. H., Surratt, J. D., Szmigielski, R., and Wildt, J.: The formation, properties and impact of secondary organic aerosol: current and emerging issues, Atmos. Chem. Phys., 9, 5155–5236,, 2009. a

Halonen, R., Zapadinsky, E., Kurtén, T., Vehkamäki, H., and Reischl, B.: Rate enhancement in collisions of sulfuric acid molecules due to long-range intermolecular forces, Atmos. Chem. Phys., 19, 13355–13366,, 2019. a, b, c, d

Hamaker, H.: The London—van der Waals attraction between spherical particles, Physica, 4, 1058–1072,, 1937. a, b

Jorgensen, W. L., Maxwell, D. S., and Tirado-Rives, J.: Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids, J. Am. Chem. Soc., 118, 11225–11236,, 1996. a

Kirkby, J., Curtius, J., Almeida, J., Dunne, E., Duplissy, J., Ehrhart, S., Franchin, A., Gagné, S., Ickes, L., Kürten, A., Kupc, A., Metzger, A., Riccobono, F., Rondo, L., Schobesberger, S., Tsagkogeorgas, G., Wimmer, D., Amorim, A., Bianchi, F., Breitenlechner, M., David, A., Dommen, J., Downard, A., Ehn, M., Flagan, R. C., Haider, S., Hansel, A., Hauser, D., Jud, W., Junninen, H., Kreissl, F., Kvashin, A., Laaksonen, A., Lehtipalo, K., Lima, J., Lovejoy, E. R., Makhmutov, V., Mathot, S., Mikkilä, J., Minginette, P., Mogo, S., Nieminen, T., Onnela, A., Pereira, P., Petäjä, T., Schnitzhofer, R., Seinfeld, J. H., Mikko Sipilä, Y. S., Stratmann, F., Tomé, A., Vanhanen, J., Viisanen, Y., Vrtala, A., Wagner, P. E., Walther, H., Weingartner, E., Wex, H., Winkler, P. M., Carslaw, K. S., Worsnop, D. R., Baltensperger, U., and Kulmala, M.: Role of sulphuric acid, ammonia and galactic cosmic rays in atmospheric aerosol nucleation, Nature, 476, 429–433,, 2011. a

Kulmala, M., Pirjola, L., and Mäkelä, J. M.: Elucidating severe urban haze formation in China, P. Natl. Acad. Sci. USA, 404, 66–69,, 2000. a

Kulmala, M., Kontkanen, J., Junninen, H., Lehtipalo, K., Manninen, H. E., Nieminen, T., Petäjä, T., Sipilä, M., Schobesberger, S., Rantala, P., Franchin, A., Jokinen, T., Järvinen, E., Äijälä, M., Kangasluoma, J., Hakala, J., Aalto, P. P., Paasonen, P., Mikkilä, J., Vanhanen, J., Aalto, J., Hakola, H., Makkonen, U., Ruuskanen, T., Mauldin, R. L., Duplissy, J., Vehkamäki, H., Bäck, J., Kortelainen, A., Riipinen, I., Kurtén, T., Johnston, M. V., Smith, J. N., Ehn, M., Mentel, T. F., Lehtinen, K. E. J., Laaksonen, A., Kerminen, V.-M., and Worsnop, D. R.: Direct Observations of Atmospheric Aerosol Nucleation, Science, 339, 943–946,, 2013. a

Kurtén, T., Loukonen, V., Vehkamäki, H., and Kulmala, M.: Amines are likely to enhance neutral and ion-induced sulfuric acid-water nucleation in the atmosphere more effectively than ammonia, Atmos. Chem. Phys., 8, 4095–4103,, 2008. a

Landau, L. D. and Lifshitz, E. M.: Mechanics, vol. 1, Butterworth-Heinemann,, 1976. a

Lehtipalo, K., Rondo, L., Kontkanen, J., Schobesberger, S., Jokinen, T., Sarnela, N., Kürten, A., Ehrhart, S., Franchin, A., Nieminen, T., Riccobono, F., Sipilä, M., Yli-Juuti, T., Duplissy, J., Adamov, A., Ahlm, L., Almeida, J., Amorim, A., Bianchi, F., Breitenlechner, M., Dommen, J., Downard, A. J., Dunne, E. M., Flagan, R. C., Guida, R., Hakala, J., Hansel, A., Jud, W., Kangasluoma, J., Kerminen, V.-M., Keskinen, H., Kim, J., Kirkby, J., Kupc, A., Kupiainen-Määttä, O., Laaksonen, A., Lawler, M. J., Leiminger, M., Mathot, S., Olenius, T., Ortega, I. K., Onnela, A., Petäjä, T., Praplan, A., Rissanen, M. P., Ruuskanen, T., Santos, F. D., Schallhart, S., Schnitzhofer, R., Simon, M., Smith, J. N., Tröstl, J., Tsagkogeorgas, G., Tomé, A., Vaattovaara, P., Vehkamäki, H., Vrtala, A. E., Wagner, P. E., Williamson, C., Wimmer, D., Winkler, P. M., Virtanen, A., Donahue, N. M., Carslaw, K. S., Baltensperger, U., Riipinen, I., Curtius, J., Worsnop, D. R., and Kulmala, M.: The effect of acid–base clustering and ions on the growth of atmospheric nano-particles, Nat. Commun., 7, 11594,, 2016. a

Leite, F. L., Bueno, C. C., Da Róz, A. L., Ziemath, E. C., and Oliveira, O. N.: Theoretical Models for Surface Forces and Adhesion and Their Measurement Using Atomic Force Microscopy, Int. J. Mol. Sci., 13, 12773–12856,, 2012. a

Loukonen, V., Kurtén, T., Ortega, I. K., Vehkamäki, H., Pádua, A. A. H., Sellegri, K., and Kulmala, M.: Enhancing effect of dimethylamine in sulfuric acid nucleation in the presence of water – a computational study, Atmos. Chem. Phys., 10, 4961–4974,, 2010. a, b

Moran, T. F. and Hamill, W. H.: Cross Sections of Ion–Permanent-Dipole Reactions by Mass Spectrometry, J. Chem. Phys., 39, 1413–1422,, 1963. a

Myllys, N., Olenius, T., Kurtén, T., Vehkamäki, H., Riipinen, I., and Elm, J.: Effect of Bisulfate, Ammonia, and Ammonium on the Clustering of Organic Acids and Sulfuric Acid, J. Phys. Chem. A, 121, 4812–4824,, 2017. a

Neefjes, I., Halonen, R., Vehkamäki, H., and Reischl, B.: Modeling approaches for atmospheric ion–dipole collisions: all-atom trajectory simulations and central field methods, Atmos. Chem. Phys., 22, 11155–11172,, 2022. a, b, c

Ortega, I. K., Kupiainen, O., Kurtén, T., Olenius, T., Wilkman, O., McGrath, M. J., Loukonen, V., and Vehkamäki, H.: From quantum chemical formation free energies to evaporation rates, Atmos. Chem. Phys., 12, 225–235,, 2012. a

Ouyang, H., Gopalakrishnan, R., and Hogan, C. J.: Nanoparticle collisions in the gas phase in the presence of singular contact potentials, J. Chem. Phys., 137, 064316,, 2012. a

Plimpton, S.: Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 117, 1–19,, 1995. a, b

Pope III, C. A. and Dockery, D. W.: Aerosols, Climate, and the Hydrological Cycle, J. Air Waste Manage., 56, 709–742,, 2006. a

Pöschl, U.: Atmospheric Aerosols: Composition, Transformation, Climate and Health Effects, Angew. Chem. Int. Edit., 44, 7520–7540,, 2005. a

Ramanathan, V., Crutzen, P. J., Kiehl, J. T., and Rosenfeld, D.: Aerosols, Climate, and the Hydrological Cycle, Science, 294, 2119–2124,, 2001. a

Sceats, M. G.: Brownian coagulation in aerosols—the role of long range forces, J. Colloid Interface Sci., 129, 105–112,, 1989. a

Schenter, G. K., Kathmann, S. M., and Garrett, B. C.: Dynamical Nucleation Theory: A New Molecular Approach to Vapor-Liquid Nucleation, Phys. Rev. Lett., 82, 3484–3487,, 1999. a

Seinfeld, J. H., Bretherton, C., Carslaw, K. S., Coe, H., DeMott, P. J., Dunlea, E. J., Feingold, G., Ghan, S., Guenther, A. B., Kahn, R., Kraucunas, I., Kreidenweis, S. M., Molina, M. J., Nenes, A., Penner, J. E., Prather, K. A., Ramanathan, V., Ramaswamy, V., Rasch, P. J., Ravishankara, A. R., Rosenfeld, D., Stephens, G., and Wood, R.: Improving our fundamental understanding of the role of aerosol-cloud interactions in the climate system, P. Natl. Acad. Sci. USA, 113, 5781–5790,, 2016. a

Sipilä, M., Berndt, T., Petäjä, T., Brus, D., Vanhanen, J., Stratmann, F., Patokoski, J., Mauldin, R. L., Hyvärinen, A.-P., Lihavainen, H., and Kulmala, M.: The Role of Sulfuric Acid in Atmospheric Nucleation, Science, 327, 1243–1246,, 2010. a

Stolzenburg, D., Simon, M., Ranjithkumar, A., Kürten, A., Lehtipalo, K., Gordon, H., Ehrhart, S., Finkenzeller, H., Pichelstorfer, L., Nieminen, T., He, X.-C., Brilke, S., Xiao, M., Amorim, A., Baalbaki, R., Baccarini, A., Beck, L., Bräkling, S., Caudillo Murillo, L., Chen, D., Chu, B., Dada, L., Dias, A., Dommen, J., Duplissy, J., El Haddad, I., Fischer, L., Gonzalez Carracedo, L., Heinritzi, M., Kim, C., Koenig, T. K., Kong, W., Lamkaddam, H., Lee, C. P., Leiminger, M., Li, Z., Makhmutov, V., Manninen, H. E., Marie, G., Marten, R., Müller, T., Nie, W., Partoll, E., Petäjä, T., Pfeifer, J., Philippov, M., Rissanen, M. P., Rörup, B., Schobesberger, S., Schuchmann, S., Shen, J., Sipilä, M., Steiner, G., Stozhkov, Y., Tauber, C., Tham, Y. J., Tomé, A., Vazquez-Pufleau, M., Wagner, A. C., Wang, M., Wang, Y., Weber, S. K., Wimmer, D., Wlasits, P. J., Wu, Y., Ye, Q., Zauner-Wieczorek, M., Baltensperger, U., Carslaw, K. S., Curtius, J., Donahue, N. M., Flagan, R. C., Hansel, A., Kulmala, M., Lelieveld, J., Volkamer, R., Kirkby, J., and Winkler, P. M.: Enhanced growth rate of atmospheric particles from sulfuric acid, Atmos. Chem. Phys., 20, 7359–7372,, 2020.  a

Su, T. and Bowers, M. T.: Theory of ion‐polar molecule collisions. Comparison with experimental charge transfer reactions of rare gas ions to geometric isomers of difluorobenzene and dichloroethylene, J. Chem. Phys., 58, 3027–3037,, 1973. a

Su, T., Su, E. C., and Bowers, M. T.: Ion–polar molecule collisions. Conservation of angular momentum in the average dipole orientation theory. The AADO theory, J. Chem. Phys., 69, 2243–2250,, 1978. a

Temelso, B., Phan, T. N., and Shields, G. C.: Computational Study of the Hydration of Sulfuric Acid Dimers: Implications for Acid Dissociation and Aerosol Formation, J. Phys. Chem. A, 116, 9745–9758,, 2012. a

Tribello, G. A., Bonomi, M., Branduardi, D., Camilloni, C., and Bussi, G.: PLUMED 2: New feathers for an old bird, Comput. Phys. Commun., 185, 604–613,, 2014. a

Vehkamäki, H. and Riipinen, I.: Thermodynamics and kinetics of atmospheric aerosol particle formation and growth, Chem. Soc. Rev., 41, 5160–5173,, 2012. a

Yang, H., Goudeli, E., and Hogan, C. J.: Condensation and dissociation rates for gas phase metal clusters from molecular dynamics trajectory calculations, J. Chem. Phys., 148, 164304,, 2018. a, b

Yang, H., Drossinos, Y., and Hogan, C. J.: Excess thermal energy and latent heat in nanocluster collisional growth, J. Chem. Phys., 151, 224304,, 2019. a

Yang, H., Song, G., and Hogan, C. J.: A molecular dynamics study of collisional heat transfer to nanoclusters in the gas phase, J. Aerosol Sci., 159, 105891,, 2022. a, b

Zhang, R., Khalizov, A., Wang, L., Hu, M., and Xu, W.: Nucleation and Growth of Nanoparticles in the Atmosphere, Chem. Rev., 112, 1957–2011,, 2012. a

Short summary
We present a new analytical model for collision rates between molecules and clusters of arbitrary sizes, accounting for long-range interactions. The model is verified against atomistic simulations of typical acid–base clusters participating in atmospheric new particle formation (NPF). Compared to non-interacting models, accounting for long-range interactions leads to 2–3 times higher collision rates for small clusters, indicating the necessity of including such interactions in NPF modeling.
Final-revised paper