A many-body dissipative particle dynamics parametrisation scheme to study behaviour at air–water interfaces †

In this article, we present a general parametrisation scheme for many-body dissipative particle dynamics (MDPD). The scheme is based on matching model components to experimental surface tensions and chemical potentials. This allows us to obtain the correct surface and mixing behaviours of complex, multicomponent systems. The methodology is tested by modelling the behaviour of nonionic polyoxyethylene alkyl ether surfactants at an air/water interface. In particular, the influence of the number of ethylene oxide units in the surfactant head group is investigated. We find good agreement with many experimentally obtained parameters, such as minimum surface area per molecule; and a decrease in the surface tension with increasing surfactant surface density. Moreover, we observe an orientational transition, from surfactants lying directly on the water surface at low surface coverage, to surfactants lying parallel or tilted with respect to the surface normal at high surface coverage. The parametrisation scheme is also extended to cover the zwitterionic surfactant lauryldimethylamine oxide (LDAO), where we provide good predictions for the surface tension at maximum surface coverage. Here, if we exceed this coverage, we are able to demonstrate the spontaneous production of micelles from the surface surfactant layer.


Introduction
The mesoscopic simulation method of dissipative particle dynamics (DPD 1 ) has found application in modelling systems where more computationally expensive methods struggle.In many soft matter systems, techniques such as molecular dynamics often have limited applicability due to the phenomena of interest occurring on time and length scales that are inaccessible using contemporary computers.Consequently, DPD has found uses in the study of block co-polymers, 2-5 lipid bilayers, [6][7][8] nanoparticles 9 thermotropic 3,10-15 and lyotropic liquid crystalline phases, 14,[16][17][18][19][20][21][22][23][24] which often involve large scale molecular organisation and slow dynamics.Molecules in DPD are coarsegrained, replacing multiple atomic sites with a single bead.This approach achieves computational efficiency via a reduction in the complexity of the system, and the use of a smooth, softrepulsive force acting between beads.
One of the main drawbacks of standard DPD is that it can not be used to simulate vapour-liquid co-existence.This is due to the repulsive nature of the interaction force between bead pairs, meaning that DPD can only be used to simulate fluid behaviour in confined spaces, and therefore cannot be used to study systems with free surfaces.In order to simulate vapour-liquid interfaces in standard DPD, researchers have often had to resort to the creation of fictitious 'vapour beads'. 25To address this limitation, many-body dissipative particle dynamics (MDPD) was developed. 26MDPD is a variation of standard DPD, in which the conservative interaction force between beads is altered to allow for liquid and vapour formation in the same simulation box.This is achieved with the addition of an attractive term in the conservative force.Since its inception, MDPD has found application in areas of research such as droplets, [27][28][29] flow over surfaces, 30 bubbles, 31 and liquid bridges. 32n additional drawback of the DPD method is that no standard parametrisation scheme exists for assigning the strength of interactions between different bead species (akin to the concept of a force field in molecular dynamics).A number of excellent DPD parametrisation schemes have been developed for various molecules, 16,[33][34][35][36][37][38][39][40][41][42][43][44] however, they are usually developed for application to a small selection of molecules (e.g.alkyl ethoxylate surfactants 19,43,44 ), rather than attempting to be more general.A lack of a standard approach makes it difficult to assign off-the-shelf parameters so that DPD can be immediately applied to any particular real system.Furthermore, there is limited existing work developing models for MDPD use, although attempts have been made to connect the interaction parameters in MDPD to real systems. 45n this work, we develop a new MDPD parametrisation scheme, with a particular focus on predicting surface tensions.Moreover, while MDPD has been applied to a variety of research areas, there is limited existing literature using MDPD to study surfactants at an interface, which will be a key focus of this study.Our MDPD model is developed by a top-down coarse-graining scheme, matching to experimental surface tension values for simple systems, in order to find the correct interaction parameters between beads.We also make use of activity coefficients at infinite dilution, in order to find the correct interactions between different types of bead species.We then test this approach by applying it to complex systems containing a mixture of bead species.In particular, we apply our MDPD model to nonionic polyoxyethylene alkyl ether surfactants at an air/water interface.These surfactants are widely used in consumer products, such as laundry detergents and shampoos, and therefore their behaviour is of particular interest.The properties of polyoxyethylene alkyl ethers are strongly dependent on the number of ethylene oxide units in the molecule, and therefore we investigate the impact of varying the size of the head group in our simulations.
5][56][57] However, MDPD provides a desirable approach to studying these systems and can provide valuable additional insights into what is often complex behaviour.For example, surfactants at an interface can be very difficult to study experimentally.There is often disagreement for surface coverage at the CMC, as determined using neutron reflection, compared with those obtained by surface tension isotherms 58,59 (where the surface tension is related to the surface excess by the Gibbs equation).
Reasons for these discrepancies are often attributed to the fact that, in order to use the surface tension method, a number of assumptions have to be made about the composition of the surface layer.This is often difficult due to the presence of impurities. 59Molecular dynamics simulations are also typically conducted with a relatively small interfacial area, and can lead to incorrect surface tension calculations depending on the force field used. 55n this article, we begin by outlining the many-body dissipative particle dynamics method, before explaining our approach for creating a new MDPD parametrisation.We investigate the performance of the model by applying it to a variety of simple systems.Following this, we test the parametrisation scheme by studying the effect that the number of ethylene oxide units in polyoxyethylene alkyl ether surfactants has on surfactants at an air/water interface.Finally, we also test a second surfactant, to highlight the transferability of the parametrisation to different types of molecules.

Many-body dissipative particle dynamics
In standard DPD, the motion of non-bonded beads is described by a combination of three forces: the conservative force, F C ij , dissipative force, F D ij and random force, F R ij .The conservative force takes the form where a ij is the tunable parameter that takes positive values and is used to capture the correct behaviour of the interaction between particles i and j. r ij is the distance separating the particles and r C is the interaction cut-off, commonly taken as unity.Since 0 o a ij , this force is purely repulsive in nature.This leads to an inability to simulate coexisting liquid and vapour phases.
The remaining two non-bonding forces are a pairwise friction-based dissipative (or drag) force F D ij and a random pairwise force F R ij which, when coupled together, act as the simulation thermostat and maintain hydrodynamics.The dissipative force F D ij and random force F R ij are given by where o D and o R are weight functions, g is a friction coefficient and s is the noise amplitude.v ij is the velocity between beads i and j, Dt is the time step and z ij (t) is a randomly fluctuating Gaussian variable with zero mean and unit variance.It was shown by Espan ˜ol and Warren 60 that the relationship between the weight functions and amplitudes must be where k B is the Boltzmann constant and T is the temperature, in order to ensure that the fluctuation-dissipation theorem 60  These potentials take the harmonic form where l 0 is an equilibrium bond length and C and D are constants defining the strength of the bond and the rigidity, which we set to values C = 150 and D = 5, which are comparable to choices made in previous works. 14,43y ijk is a bond angle between beads i, j and k, and y 0 is an equilibrium bond which we set to y 0 = 1801.In this work, we set equilibrium bond length l 0 = 0.5, for reasons which will be discussed later in the paper.Therefore, the total force f i acting bead i can be written as where F ij are the forces acting on bead i by bead j.
MDPD uses the same random, dissipative, and bonding forces as the standard DPD method.The difference between DPD and MDPD comes in the form of the conservative force, which is modified to include a density-dependent component.The new conservative force takes the form where r are local densities for each particle, and are calculated as For MDPD, the cutoff r C is usually taken as r C = 1 and the many-body cutoff r d is taken as 0.75. 27The conservative force parameter a ij takes a negative value, acting as an inter-particle attraction.The B parameter takes a value 0 o B, meaning that the density dependent part of the conservative force is repulsive.The attractive pairwise term and repulsive many-body term together create a potential that has both attractive and repulsive regions, which can reproduce a van der Waals loop.Hence, MDPD has the capacity to simulate vapour liquid interfaces. 27The no-go theorem 61 imposes global B and r d parameters for all particle types, and therefore, in this work, we chose to take the common choices of B = 25 and r d = 0.75. 27,62,63Tuning the a ij parameter has an impact on both the surface tension and the mixing energies of each particle.Therefore in order to develop our parametrisation scheme, we first need to obtain relationships between a ij , the chemical potential and the surface tension.All simulations described are performed using the DL_MESO software package. 64Calculating a ij : a parametrisation scheme For a system which consists of S different bead species, this requires the calculation of S self-interaction values for the a ii parameter, and S(S À 1)/2 cross interaction values a ij .We approach the calculation of the self-interaction and crossinteraction values differently, as discussed in this section.

Self-interactions
The a ii self-interaction parameters are calculated by finding the a ii value which reproduces (in MDPD) the experimentally known surface tension for a particular fluid.We wish to create a general approach, such that we have a method of finding an a ii value for any possible liquid we may wish to model, provided we have the experimental surface tension.We aim to model water such that 1 bead represents three water molecules, and that long chain molecules are modelled as a series of bonded beads.In this work we focus on molecules which can be represented as a linear chain of beads (i.e. with no branching).
In order to find a general relationship between a ii and the surface tension, we perform a parameter sweep in the range 25 r a ii r 52, with an interval of Da ii = 3.These simulations are conducted in a box with dimensions 10r C Â 10r C Â 100r C using n = 10 000 simulation beads, which are initialised with a random initial configuration.We also investigate the influence of bonding on pure systems, by bonding N beads together in a linear chain (N À 1 bonds per chain).We vary N in the range 1 r N r 8, keeping the overall total number of beads n approximately constant.Simulations are performed in the NVT (constant volume and temperature) ensemble, the time step used is Dt = 0.01, and the mass of all beads m = 1.
After running the simulation from an initially random placement of molecules, all simulated cases with different a ii values eventually form separated regions of liquid and vapour in the simulation box, allowing surface tension to be calculated.An example of the formation of liquid-vapour coexistence from an initial configuration is shown in Fig. 1, where the interface forms in the x-y plane (note that throughout this work we use VMD 65 software for creating visualisations).We run these simulations for I = 1.5 Â 10 6 time-steps and, on average, the beads separate into a bulk fluid and vapour after I E 1 Â 10 5 , at which point we begin collecting data for calculating the surface tension.The surface tension g can be calculated using the following definition [66][67][68] g ¼ where p N (z) = p zz (z) is the local normal pressure and p T (z) = p xx (z) = p zz (z) is the local tangential pressure.Calculation of g via this method requires the calculation of the local values of pressure tensor components.Alternatively, for a system with two interfaces, it can be shown 68 that the surface tension is related to the macroscopic averages where P N and P T are macroscopic normal and tangential components.Therefore, we use eqn (12) to calculate g in this work.The pressure in DL_MESO is calculated using the virial theorem, and the individual components are calculated as where a = X, Y, Z. Pressure components P N = P Z and P T = (P X + P Y )/2, and V(t) is box volume.From our parameter sweep, it is determined that surface tension g can be fitted by an equation of the following form where b is a function of interaction parameter a ii (for more details of this fitting see ESI † Section S1).We find that b can be reasonably fit by a linear expression of a ii : b = ma ii + c where m = À0.4262 and c = À7.272.However, a second-order polynomial is required for more accurate surface tension prediction at high and low values of a ii .This is particularly important for high surface tension fluids such as water.Therefore g is predicted using a second-order expression for b, resulting in the following relationship between g and a ii : where g is in units of k B T/r C 2 .A similar relationship between a ii and the bead number density r can be obtained (for more details of this fitting see ESI † Section S2), taking the form Traditionally, conversion from DPD units into real units is performed by matching the DPD mass density, to the room temperature experimentally known value, thereby defining a value of the length scale r C in real units.Therefore, it is difficult to match both the surface tension and the mass density of every particle species in the system, to those known from experiment.However, for surfactant systems, water is the universally most important solvent, and we, therefore, chose to experimentally match to both the surface tension and density for water at room temperature, in order to obtain a value for r C .Using a value of 72.0 mN m À1 for the surface tension of water at room temperature, k B T = 4.11 Â 10 À21 J, N = 1, and a coarse graining of 3 water molecules per bead (mass of bead m = 54 g mol À1 ), allows us to solve eqn (15) and ( 16) to obtain values r C = 8.53 Å and a WW = À50.683.For all other bead types, we continue to use r C = 8.53 Å for the purpose of converting into real units, and instead obtain values of a ii for other bead types by using eqn (15) only.

Cross interactions
A different approach is required to find values for the mixing interaction parameters of a ij when i a j, since it is unlikely that we will always have an experimental surface tension for a mixture of every bead pairing required.In order to find a ij we calculate activity coefficients at infinite dilution (g N ).We calculate ln(g N A ) using the relation where m É A is the chemical potential of pure component A, and m Ã A is the chemical potential of A in the mixture at infinite dilution (where the units of chemical potential m are energy per particle).For computational purposes, the chemical potential of component A in a mixture can be represented using the sum where m ideal A is the chemical potential of an ideal solution and m excess A is the excess chemical potential.The ideal component can be calculated using m ideal A = k B T ln r A L 3 where L is the thermal de Broglie wavelength and r is the density.The excess potential energy m excess is calculated via the Widom insertion method as where the change in potential energy DU due to the insertion of a single DPD bead is measured and an ensemble average can be used to calculate the excess chemical potential per particle.Therefore (setting where DU AA is the energy change that results from insertion of particle A into a bath of particle A and DU AB is the energy change from inserting particle A into a bath of particle B. In order to calculate ln(g N A ) as a function of a AA , a BB and a AB , we perform a parameter sweep in the range 25 r a ij r 52 for all three interaction parameters.In contrast to the previous simulations for calculating a ii , these simulations are performed using the NPT ensemble in a cubic simulation box.For these simulations, there is no interface and the domain is entirely bulk fluid.Performing NPT simulations allows us to simulate bulk fluid at the correct density for different interaction parameter values.The simulation box is initialised with dimensions 10r C Â 10r C Â 10r C using n = 10 000 simulation beads.We equilibrate for I = 1.5 Â 10 6 time-steps using a time step of Dt = 0.001.The domain equilibrates extremely rapidly (I o 3000) to a constant average density before Widom insertions are performed.
Following density equilibration, Widom insertion allows us to determine that ln(g N A ), as a function of the interaction parameters, can be represented by an equation of the form where x i are fitted constants: x 1 = 0.008331, x 2 = À0.001588,x 3 = À0.01282,x 4 = À0.2158,x 5 = À0.5868,x 6 = 0.3641 and x 7 = À7.684(see ESI † Section S3 for details).Therefore, if one has the computationally (or experimentally) obtained value for the activity coefficient g N A for a particular bead species, and also the values of a AA and a BB from using eqn (15), then eqn (21) can be used to calculate the required a AB parameter for any bead pair.Therefore, the combination of eqn (15) along with eqn (21), allows us to calculate the required a ij parameters for all bead interactions in the simulation.
This journal is © The Royal Society of Chemistry 2023 3.3 a ij parameters for C 12 E j surfactants Since we aim to be able to apply the parametrisation to predict properties of C 12 E j surfactants at an air/water interface, in this section we calculate the required a ij parameters for these molecules.Polyoxyethylene alkyl ethers H(CH 2 ) i (OCH 2 CH 2 ) j OH (often abbreviated as C i E j ) are extremely common nonionic surfactants.We coarse grain the surfactant tail such that one alkyl bead (C3) represents three carbon atoms and one head bead (EO) represents a single ethylene oxide unit [OCH 2 CH 2 ].The final tail bead of an alkyl chain (C3T) is slightly different to the C3 beads, as it also includes an additional hydrogen atom.Likewise, the end bead (EOT) of the ethylene glycol portion is treated differently, in order to represent the [OCH 2 CH 2 OH] unit which terminates the molecule.This coarse graining procedure is shown in Fig. 2.
We coarse-grain water beads (W) such that one bead represents three water molecules.
Assuming that the C-C-C angles are tetrahedral (approximately 109.51) and that the C-C experimental bond is 1.543 Å, 69 the separation between three carbon atoms in the alkyl chain (or two adjacent beads) is expected to be around 3.8 Å.The choice of bond length l 0 = 0.5 (see Section 2) results in a [CH 2 CH 2 CH 2 ]-[CH 2 CH 2 CH 2 ] bead separation of l 0 = 4.27 Å, which is close to the experimental value.The bond length l 0 = 0.5 is also used between beads in the polyethylene glycol portion of the surfactant.Experimentally the distance between successive ethylene oxide units is expected to be around 3.7 Å, 70 making this a suitable choice.We note that we also investigated the effect of varying l 0 on eqn (15), finding that the bond length actually has very little impact on the surface tension (see ESI † Section S5 for more information).
The values for the activity coefficients at infinite dilution are calculated via COSMO-RS [71][72][73] calculations, and the values are presented in Table 1.COSMO-RS uses quantum chemical calculations in order to model the distribution of charge on molecules, and when combined with statistical thermodynamics can be used to calculate the activity coefficient of a component. 74nce the self-interaction values for each bead species have been determined, the values in Table 1 and eqn (21) are used to determine cross interaction values.
In order to find the self-interaction for the EOT bead, a EOT,EOT , we match to the experimental surface tension g = 45.65 mN m À1 75 for diethylene glycol (note that all experimental surface tensions are taken at room temperature which we define as 25 1C).We approximate the coarse-grained representation of diethylene glycol as two bonded EOT beads (N = 2), using eqn ( 14) and ( 15) to determine a EOT,EOT = À39.951.
For the tail beads, we match to experimental surface tension g = 25.3 mN m À1 for dodecane, [76][77][78] using four beads (N = 4) to represent this molecule in the coarse-grained representation.We choose to set the self-interactions between all tail beads to be equivalent (i.e. a C3,C3 = a C3T,C3T = a C3,C3T ) for ease of matching to experimental surface tension data.Therefore, we determine a self-interaction a C3,C3 = À29.796.
For single ethylene oxide beads, we chose to match to the surface tension of tetraethylene glycol g = 45.13 mN m À1 , 79 represented as EOT-EO-EO-EOT.This molecule is more complex than eqn ( 14) allows for, since the initial parametrisation was developed for bonding between only like bead types.This was applicable for the alkyl chain, as we do not expect there to be a significant difference in the strength of the interaction between the C3 and C3T beads (note that this is supported by the fact that the activity coefficient at infinite dilution of C3T in C3 is found to be nearly zero ln(g N A ) = 0.0004, indicating that these beads can be treated as equivalent).However, it is expected that there is a great deal of difference between EOT and EO beads.We expect that the hydroxyl groups (OH) play a large part in influencing the surface tension, as interactions between these groups display strong hydrogen bonding.This is highlighted by the relatively large surface tension of ethylene glycol (g E 48 mN m À1 75 ), relative to longer polyethylene glycols.
Therefore in order to calculate a value for a EO,EO using experimental data for tetraethylene glycol, we extend the parametrisation to be able to represent molecules bonded with more than one bead type.We perform an additional parameter sweep to represent molecules of form D-(C) l -D where the number l of middle beads C is varied.We perform the parameter sweep for variables a CC , a DD and a CD in the range 28 r a ij r 44.The simulations are conducted in the same way as previously described in Section 3.1.When a CD % o 33 the molecules are observed to crystallise into neat periodic layers, and so these cases are excluded from our general formula fitting, and we fit to simulation cases with 34 r a CD only.
We determine that the surface tension can be fit using an equation of the form where N = l + 2 and the variables fitted are given in Table 2 (see ESI † Section S1 for fit details).Using a combination of eqn ( 21) and ( 22), along with the experimental surface tension for tetraethylene glycol and the value of the activity coefficient at infinite dilution in Table 1, allows us to find values a EO,EO = À37.182and a EO,EOT = À38.847.A table of the final calculated a ij parameters for every bead pair is given in Table 3.Note that for each bead pair we must choose one bead to be the solute (particle A in eqn (21)) and one to be  the solvent (particle B in eqn ( 21)) when calculating the cross interaction.We generally choose to define the solvent as the bead which is more abundant (when compared with the abundance of the other bead) in each bead pairing.For example, for a EO,EOT the EOT bead is selected as the solute and EO bead is the solvent.However, since we are varying the number of ethylene oxide groups in the hydrocarbon chain, the relative abundance of head to tail beads varies for our different simulations in this study.For consistency, we chose to use the head groups as the solute and tail beads as the solvent when calculating the a EO,C3 and a EOT,C3T interactions, therefore using the same interaction parameters across all C i E j surfactants.

Simulation set-up
Before performing simulations of C 12 E j surfactants at an interface, we test the parameters in Table 3 on simpler cases, such as the surface tension of pure fluids or binary mixtures.These simulations are conducted in a box with dimensions 22r C Â 22r C Â 100r C using n = 100 000 simulation beads.We note that despite the fact that this surface area is larger than that used to generate eqn (15), the surface tension calculated is independent of the size of the interface (see ESI † Section S4 for details).This is consistent with what is concluded by other researchers. 80he simulations are initialised in a similar way to that described in Section 3.1, that is with a random initial configuration.Molecules form coexisting bulk and vapour phases with two interfaces after around I E 1.5 Â 10 5 time-steps, and in total we run I E 1 Â 10 6 time-steps for collecting data.An exception to this box size is for when we perform a single simulation of a mixture of dodecane and water, which is expected to exhibit phase separation.In this case, we shorten the length of the simulation box and, we use dimensions 22r C Â 22r C Â 34.03r C using n W = 50 000 and 12500 dodecane molecules (n = 100 000 simulation beads in total).This box size is chosen such that the water and dodecane fill the simulation box to avoid the formation of dodecane/gas interfaces.
For simulations of surfactants at an air-water interface, surfactants are placed at a pre-equilibrated air/water interface, where the bulk water is generated in the same way as previously described for pure systems (dimensions 22r C Â 22r C Â 100r C using n = 100 000 simulation beads).Following the generation of the interface, the surfactants are arranged at both surfaces in a grid-like configuration, as shown in Fig. 3.
The initial configuration places the head groups partially submerged in bulk water, while the tails protrude outwards into the vapour phase, perpendicular to the surface of the water.We then run the simulation for I = 1.5 Â 10 6 time steps, collecting stress tensor values every time step in order to calculate the surface tension.We collect positional data every 1000 time steps, which is used to calculate various properties such as the density profile and orientation of molecules relative to the surface.

Complex molecules and binary mixtures
First, the surface tension is calculated for a variety of molecules as pure systems (i.e.single component molecule fluids), in order to confirm that complex molecules can be simulated.Table 4 shows calculated surface tension values for dodecane (C 12 ), tetraethylene    glycol (E 4 ) and three polyoxyethylene alkyl ether surfactants, alongside a comparison with experimental values.Overall, we observe good agreement with the available experimental data.Additionally, the parametrisation scheme is also tested on its ability to predict the surface tension of aqueous mixtures of polyethylene glycol.The behaviour of these mixtures is of importance to confirm since the ethylene oxide head group in C i E j surfactants is weakly hydrophilic.Therefore, it can be difficult to predict the correct interaction parameters for these bead types.Here, simulations are set up and initialised in the same manner as the pure simulation case.Fig. 4 shows the surface tension calculated for aqueous mixtures of diethylene glycol, triethylene glycol, and tetraethylene glycol.The surface tensions predicted show very good agreement for triethylene glycol with the experimental data across the composition range and reasonably good agreement for the other two compounds, indicating the correct interaction between ethylene glycol beads and water.The predicted surface tension is most accurate at high and low concentrations, with the mid-range of concentrations showing the most deviation from the experimental measurements.This is to be expected due to the approach of the parametrisation, where self-interactions dominate at high and low concentrations.Additionally, our cross-interactions are calculated using beads at infinite dilution.Therefore, these values may not be fully applicable in the mid-concentration range, and are expected to perform best at lower concentrations.Theoretically, this means that the bulk concentration (the middle of the fluid region) may be different to the overall concentration which is plotted in Fig. 4.However, we find that the thickness of the layer of excess polyethylene glycol molecules near the surface is extremely thin, and therefore the bulk concentration is fairly comparable to the overall concentration.Similarly, in experimental measurements, the surface tension is reported in relation to the bulk concentration because the fraction of surfactant/solute molecules at the interface is minuscule compared to the number of molecules in solution.The accumulation of molecules near the interface also implies that the MDPD surface tension results may be different when the bulk region is extended in the z direction.Therefore, we investigated the impact of box length.Once again, as the excess layer is so thin, we find that there is minimal impact on the results (see ESI † Section S6 for more information).
We find that there is a slight preference for solute molecules to gather near the surface, particularly at lower concentrations.An example of this for diethylene glycol is shown in Fig. 5.
Finally, we also test a mixture of dodecane and water.Fig. 6 shows the expected phase separation which results.We observe there are a small number of water beads in the dodecane bulk, while no dodecane beads reside in the water bulk.This is to be expected given that, while the solubility of each species in the other is very low, the solubility of water in dodecane (65 ppm 83 ) is significantly larger than that of dodecane in water (3.7 ppb 84 ).We find a ratio of E1 water bead for every 930 dodecane beads in the dodecane region, resulting in a mass ratio of E0.0014 g water per g dodecane, meaning that the concentration of water in dodecane in simulation is somewhat larger than that of experiment (E20 times larger than experiment).The cross interactions between water and alkane beads are calculated by performing insertions using single C3 beads.We believe that the discrepancy between simulated and experimental solubility may be, in part, due to the fact that we use single (rather than bonded) beads in the parametrisation method.This is supported by the values we calculate for the solubility of water in other alkanes.We performed additional Fig. 4 The surface tension calculated using MDPD for diethylene glycol, triethylene glycol, and tetraethylene glycol aqueous solutions.Calculated values are compared with experimental data. 81,82g. 5 Density profiles for various bead types in aqueous mixtures of diethylene glycol.Profiles are plotted as a function of distance from the centre of the water slab and are symmetric about z = 0. simulations, similar to the one described here for dodecane and water, for nonane and hexane.For reference, the experimental values of water solubility are 79 ppm 83 and 94 ppm 85 in nonane and hexane respectively.We calculate the water solubility to be 0.00093 (E12 times larger than experiment) in nonane and 0.00070 (E7 times larger) in hexane.This shows that increasing the length of the molecule widens the gap between the simulated and experimental results.Therefore, it may be possible to improve the calculated solubility by basing the a ij values off the activity coefficients for longer molecules, e.g.dodecane.
6 Application to polyoxyethylene alkyl ethers (C 12 E j ) In this section, the parametrisation is tested on C 12 E j surfactants at the air/water interface for different surfactants with varying j.

Density profiles
The surfactant behaviour at the interface is shown for two different surface area coverages in Fig. 7.We observed that when the interface has a low density of surfactants (high surface area per surfactant), the surfactant tails tend to remain in contact with the surface of the water.When the number of surfactants at the surface is increased (low surface area per surfactant) the alkyl tails re-orientate such that they are directed towards the gas phase, perpendicular to the liquid-gas interface.This is reflected in the density profiles as shown in Fig. 8, where there is a much greater overlap of the head and tail profiles when the surface coverage is low.When the surface density is increased we also see an enhanced separation between the water and surfactant profiles (particularly for tails).This orientational behaviour at different surface densities is consistent with that observed in similar molecular dynamics studies. 55In all cases, the profiles for tail groups overlap to some extent with those of the head groups.This implies that there is no complete segregation between surfactant head and tails, as is often assumed in theoretical models. 86Theoretical models typically require some knowledge of the behaviour (and thickness) of the adsorption layer, and incorrect assumptions in models can lead to differences between theoretical predictions and experimental results.

Surface tension
The surface tension is calculated using eqn (12) and plotted as a function of surface coverage for various surfactants in Fig. 9.We observe a smooth decrease in the surface tension with increasing surface coverage.As the number of surfactants at the surface increases, we expect to reach a point where maximum packing is exceeded.The surface tension decreases until the surface reaches what we conclude to be its maximum packing.When the number of ethylene oxide units in the surfactant molecule is large (6 Z j), we find that there is a plateau in the surface tension, and after this plateau, the surface tension begins to increase with increasing surfactant surface density.In contrast, at lower numbers of ethylene oxide units ( j = 2 and j = 4), when we attempt to initialise simulations with surface densities greater than a maximum density cutoff point (o0.0388mol Å À2 ) the simulations become unstable and the surface breaks apart.In these cases, no plateau is observed.
Generally, as the number of ethylene oxide units increases, the surface area at maximum packing also decreases, which is experimentally the expected behaviour.However, the overall  surface tension drop caused by the surfactants at the interface is lower than that calculated via experiment.The lowest surface tension found for a variety of surfactants is compared with experiment in Table 5.For surfactants with j = 2 and j = 4 we present in the table the surface tension calculated at surface density 0.0388 mol Å À2 , however, it is important to note that this is not necessarily comparable with surfactants with 4 o j since there is evidence that the surface is unstable at this high packing.
In principle, for large enough simulations, the critical micelle concentration could be established from the levelling of the surface tension results.However, due to the relatively small amount of bulk water, this isn't a practical approach to determining the CMC.The CMC calculated would differ greatly from that of an experimentally determined value, as in our simulations the surface layer-to-bulk ratio is much larger than that of an experimental vessel.It can be calculated using the results shown in Fig. 9, that this would result in a CMC value which considerably exceeds the reported values for C 12 E j surfactants. 94Furthermore, even a single molecule in the bulk water would have a concentration of E0.009 wt%, exceeding the expected CMC value of around 0.006 wt%.

Surfactant layer thickness
Due to the spread of partial density profiles, it is nontrivial to define a layer boundary (and thickness) of the surfactant layer at the surface.However, in order to best compare with existing experimental data, 48,52 we obtain the layer thickness by fitting a Gaussian distribution to the density profiles.This profile takes the form where s is the full width at 1/e of the maximum height, and z is the position of the center of the peak.The thickness can then be characterised using the fitted s parameter.
In each simulation, a value is calculated for the thickness of the ethylene oxide group (s h ), the alkyl chain (s t ), and the overall thickness (s s ).This is plotted as a function of surface coverage in Fig. 10 (plotted up to the maximum surface coverage Fig. 9 Surface tension as a function of surface coverage for a variety of C 12 E j surfactants.Note that dashed lines are intended to be a guide to the eye.chieved before the surface tension begins to increase).It is shown that the thickness of each component increases largely linearly with the number of molecules at the surface.The thickness of the alkyl chain sublayer is relatively consistent between different C 12 E j surfactants, which is to be expected since the tail length is constant for all surfactants.The thickness of the head group portion increases with additional head group beads, this then influences the overall thickness of the surfactant layer.
It is of note that, as a result of the overlap of the head and tail profiles demonstrated in Fig. 8, the total thickness is less than the sum of the head and tail thicknesses.Table 6 compares a selection of calculated layer thicknesses with experimental values, showing good agreement between simulation and experiment.The roughness correction provided to the experimental data originates from the fact that the experimentally measured thickness is contributed to by both the intrinsic thickness of the surfactant monolayer (in the direction normal to the surface) and a surface roughness.The surface roughness is thought to primarily originate from capillary waves; thus, the contribution level can be estimated from the surface tension.Therefore, the removal of this roughness contribution allows one to more accurately determine the actual thickness of the layer.This is discussed in more detail in Lu et al. 484 Surface coverage Fig. 9 highlights the differences in the equilibrium surface coverage between the different surfactants tested.We chose to estimate our minimum area per molecule for j 4 4 surfactants as the minimum area per molecule before surface tension begins to increase.The calculated values, along with comparable experimental values, are presented in Table 7.The minimum surface area per surfactant increases as the number of ethylene oxide units increases, in agreement with experimental trends.
For surfactants with j r 4, the surface becomes unstable at high surfactant density, although in these cases we observe that the surface can be overloaded more easily.At high density, the surface tension continues to decrease, however the surface becomes heavily distorted so that it can no longer be considered a planar interface.Instead, the surface can become extremely convexed, such that the effective surface area of the surfactant layer is increased.This causes the surface tension to appear to be artificially low, due to the increase in area.Furthermore, eqn (12)  assumes that the interface is planar to the z-axis, and therefore isn't valid for calculating the surface tension in these cases.Experimentally, we expect the surface to have a surface area per molecule of around 33 Å 2 and 44 Å 2 for C 12 E 2 and C 12 E 4 surfactants respectively. 52These experimental values are in reasonable agreement with the surface coverage value when the smooth relationship between surface tension and surface coverage begins to break down in Fig. 9.We also observe that the surface becomes unstable at this point from the thickness of the head group layer shown in Fig. 10.

Tilt angle
In order to analyse the orientation of surfactant molecules at the interface, we calculate the average director vector of the surfactant head and tail components, relative to the surface normal (z-direction).The head vector ( v h ) is defined as the vector which connects the hydroxyl group (EOT) and the EO atom that is connected to the alkyl tail.Similarly, the tail vector ( v t ) is defined as that connecting the first carbon bead (C3) and the final bead in the alkyl tail (C3T).We observe that the tilt angle rotates from its initial configuration (where y h = y t = 0), Fig. 10 The thickness of the total surfactant layer, as well as the thickness of the head and tail groups individually.This journal is © The Society of Chemistry 2023 and therefore we begin collecting data for this calculation after I = 5 Â 10 5 time-steps when the tilt angle has been observed to reach a stable value.Fig. 11 shows the angle between the average vector and the surface normal, where generally y h and y t decrease as the number of surfactants at the interface increases.This is expected from the observed behaviour (discussed in Section 6.1), where, at low surface densities surfactants lay directly on the water surface.It is interesting to note that the tail angle y t is almost always larger than the head group angle y h , indicating that the surfactant heads are more perpendicular to the interface than the tails are.This somewhat counter-intuitive observation has also been reported in MD studies. 55For surfactants C 12 E 2 and C 12 E 4 , the surfactant molecules become almost perfectly perpendicular to the surface at high surface density.However, for other surfactants, the tilt angle generally decreases to a minimum, before increasing again with greater surfactant density.
Experimentally, the tilt angle can be difficult to accurately determine, particularly at low surface coverages.However, at the CMC surface coverage, the alkyl chain tilt angle for C 12 E 12 is reported to be around y t = 451, 48 which is in agreement with the plateau at higher surface coverages for C 12 E 12 in this study.The plateau to larger tilt angles that we find are also consistent with MD studies, where values of y h = 56.81 and y t = 541 for C 12 E 6 surfactants are reported at CMC packing. 54Application to other surfactants One might expect that when the expected maximum packing at the surface is exceeded, molecules are removed from the surface into the bulk.This would increase the minimum area per surfactant and stabilise the surface.The formation of micelles within the bulk should also be expected from the ratio of surfactant molecules to water molecules in the simulation box.For example the expected CMC for C 12 E 8 is around 0.08 mM, 94 corresponding to around 10 surfactant molecules in our simulations.While this phenomenon is not observed for C 12 E j at the air/water interface, we observe this behaviour when different molecules are tested with other head groups.This surface removal is followed by the formation of micelles in the bulk water.We theorise that the lack of surfactant removal from the surface for C 12 E j surfactants is related to the (relatively) weak hydrophilicity of the ethylene oxide head group.
Consider zwitterionic surfactant lauryldimethylamine oxide (LDAO), an amine oxide based surfactant with a dodecyl alkyl tail.This molecule is coarse grained as shown in Fig. 12, where the [C 2 H 6 NO] head group is represented by a single bead (AO).Here we also demonstrate an approach to finding a value for the selfinteraction a ii when the surface tension of a pure component is unavailable from existing literature or is unobtainable from experiments.The self-interaction parameter for the AO head group is calculated by matching to experimental data for trimethylamine N-oxide (TMAO), however, anhydrous TMAO exists as a solid at room temperature.Therefore we must calculate a value for a AO,AO by matching to the surface tension for aqueous TMAO, where at a concentration of 15 wt% TMAO solutions have an experimental surface tension of 69.3mN m À1 . 97,98The values for ln(g N A ) are once again calculated using COSMO-RS, where ln(g N A ) = À10.394for AO beads at infinite dilution in water, and ln(g N A ) = 7.0388 for AO beads in baths of hydrocarbon tails at infinite dilution.The amine oxide head group, therefore, displays much more hydrophilic behaviour than the ethylene oxide head groups.
In order to find a AO,AO we perform a parameter sweep over a AO,AO in the range À34 r a AO,AO r À26 at intervals of Da AO,AO = 1, and for each value the corresponding value of a AO,W is calculated using eqn (21).Simulations are performed at a concentration of 15 wt% TMAO.This allows us to determine values a AO,AO = À30.461and a AO,W = À48.505.Finally, the interactions with the hydrocarbon tail are calculated as a AO,C3 = a AO,C3T = À21.747.
The calculated surface tension as a function of surface coverage is given in Fig. 13.There is a plateau in the surface tension observed at higher concentrations (initial surface densities Z0.0328), at which point the surface removes surfactants into the bulk.Fig. 14 shows the formation of micelles at high surfactant densities.
The surface tension at the CMC is experimentally reported as between 25-35 mN m À1 , 99-101 therefore our plateau to Fig. 11 The tilt angle of the head y h and tail y t groups (relative to the surface) of various C 12 E j surfactants at different surface coverages.around E26 mN m À1 at higher concentrations is consistent with the available experimental data.There is a brief drop (to E19 mN m À1 ) in the surface tension before this plateau at high concentrations.This 'over drop' in surface tension is interpreted to be a result of no (or limited) surfactant removal from the interface.Table 8 shows the number of surfactants at the interface following the removal of surfactants into the bulk water.We observe that surfactants are removed from the surface to produce a nearly constant surface density of E0.026 mol Å À1 .This results in a consistent value calculated for the surface tension.For the data points where the surface tension drops below 20 mN m À1 (where the number of molecules placed at the surface are 1024 and 1089 molecules/interface), the density is significantly higher due to a lack of surfactant removal to the bulk.Therefore the surface needs to have an initialised surface density, which is significantly larger than its desired equilibrium density, in order to see surfactant removal from the surface on the time scale of the simulations.

Conclusions
In this work, we have presented a parametrisation scheme for many-body dissipative particle dynamics simulations, which can be used to calculate the surface tension for complex systems.Our application of the parametrisation to C 12 E j molecules at an air/water interface allows us to easily study their behaviour at different surface coverages.We obtain good agreement with experimental results, in particular, the trend for surface tension drop and CMC surface coverage for increasing numbers of ethylene oxide groups j.These systems can be difficult to study both experimentally and by using existing molecular dynamics techniques, and therefore MDPD provides a valuable tool for studying surfactant behaviour at the air/water interface.
We also find that the parametrisation scheme introduced here provides good predictions for quantities such as the thickness of the surfactant layer and the angle of the surfactants relative to the surface, gaining insight into how the molecules pack at the surface.
We extended our parametrisation scheme to cover lauryldimethylamine oxide surfactants at the air/water interface.In this case, the surface spontaneously removes surfactants into the bulk water when the surface is initialised with a large surface density, in order to optimise the number of surfactants at the surface.This behaviour is difficult to observe in traditional molecular dynamics simulations, where it is far more common to observe an over-drop in the surface tension associated with a distortion of the interface.
It is useful to make some general comments on the applicability of this scheme to other molecules.In this work, we chose to use a bond length of l 0 = 0.5r C , which was chosen because it would produce a reasonable bead-bead bond length for a number of different systems.When applying the parametrisation to other molecules, one should take care that the coarse graining chosen produces a bond length which is relatively close to the one which the parametrisation was developed for.Eqn (15) and (22) in particular may differ if a significantly different bond length is chosen.Furthermore, bonded beads   are not only subject to bonded interactions (eqn ( 6) and ( 7)) but also non-bonded interactions with each other.These nonbonded interactions between bonded beads may result in an alteration of the bond length in the simulation to a value other than exactly l 0 .Similarly, this parametrisation was developed for linear molecules only but could be extended relatively simply for branched molecules through alteration of eqn (22).Moreover, in this work, we have not included any separate electrostatic interactions, and therefore the parametrisation has not been tested on ionic surfactants or other systems with charged molecules.The extension of this model to charged systems is currently under investigation.Finally, we note that we have not attempted to use the model to study high surfactant concentrations beyond the (post-CMC) concentration regime when the surface is fully packed and micelles form.However, the emergence of micelles from overpacked LDAO surfaces suggests the possibility of using this parametrisation scheme for studying micellar aggregation numbers and micelle shape.
is satisfied.In this work, we set k B T = 1 and set the values of constants to be s = 3 and g = 4.5.Bonding forces are calculated from potentials U B ij and U A ij .U B ij chemically bonds beads to form long chain molecules, while U A ij introduces rigidity within the chain or between the beads.

Fig. 1
Fig. 1 Formation of coexisting liquid and vapour phases from an initially random configuration.

Fig. 3
Fig. 3 Initial placement of surfactants at the surface: example given for C 12 E 4 surfactant with 625 surfactants per surface (56 Å 2 per molecule) where beads are coloured according to type: water (blue), head groups (green) and tail groups (purple).

Fig. 6
Fig. 6 Phase separation in the simulation of water (blue) and dodecane (purple).

Fig. 7 C
Fig. 7 C 12 E 8 surfactants at the air/water interface where beads are coloured according to type: water (blue), head groups (green) and tail groups (purple).Plots at surface area coverages 61 Å 2 (a) and 352 Å 2 (b) per molecule.

Fig. 8
Fig. 8 Density profiles for various bead types in C 12 E 8 surfactants at an air/water interface with average surface coverages 61 Å 2 (a) and 352 Å 2 (b) per molecule.Profiles are plotted as a function of distance from the centre of the water slab and are symmetric about z = 0 (profile of one surface shown).

Fig. 13
Fig. 13 Surface tension as a function of surface coverage for LDAO.

Fig. 14
Fig.14The emergence of micelles into the bulk when LDAO surfactants are placed at an air/water interface.Example given is for a system initialised with 1369 surfactants per surface (26.5 Å 2 per molecule) where beads are coloured according to type: water (blue), head groups (yellow) and tail groups (purple).

Table 1
Activity coefficient at infinite dilution lng of bead i in bead j

Table 2
Fitted parameters for calculating surface tension using eqn(22)

Table 3
Calculated a ij parameters for beads making up surfactant molecules of the form C i E j

Table 4
Surface tension calculated for simulations of pure molecules, compared with experimental values This journal is © The Royal Society of Chemistry 2023

Table 5
Lowest surface tension achieved for C 12 E j surfactants

Table 6
A comparison of a selection of simulated and experimental values at the CMC for the ethylene oxide group thickness s h and the thickness of the alkyl chain s t .Experimental values are corrected for 'capillary wave roughness' producing s Ã h and s Ã t , which are expected to be closer to a realistic value for the thickness.Simulation thicknesses are obtained at surface coverage values: 32.3 Å 2 (C 12 E 2 ), 45.0 Å 2 (C 12 E 4 ), and 72.8 Å 2 (C 12 E 12 )

Table 7
Minimum surface area per molecule for C 12 E j surfactants at maximum packing

Table 8
The number of surfactants per interface (and corresponding density) when the surfactants are initialised at the surface, compared with when molecules have been removed into the bulk following equilibration