Finite size effects and optimization of the calculation of the surface tension in surfactant mixtures at liquid/vapour interfaces

The surface tension of monolayers with mixtures of anionic and nonionic surfactant at the liquid/vapour interface is studied. Previous works have observed that calculations of the surface tension of simple ﬂuids show artiﬁcial oscillations for small interfacial areas, indicating that the surface tension data ﬂuctuate due to the ﬁnite size effects and periodic boundary conditions. In the case of simulations of monolayers composed of surfactant mixtures, the surface tension not only oscillates for small areas but can also give non-physical data, such as negative values. Analysis of the monolayers with different surfactant mixtures, ionic (DTAB, CTAB, SDS) and nonionic (SB3-12), was done for density proﬁles, parameters of order and pair correlation functions for small and large box areas and all of them present similar behaviour. The ﬂuctuations and the non-physical values of the surface tension are corrected when boxes with large interfacial areas are considered. The results indicate that in order to obtain reliable values of the surface tension, in computer simulations, it is important to choose not only the correct force ﬁeld but also the appropriate size of the simulation box.


Introduction
Surfactant molecules have been extensively studied in different interfacial problems not only for their scientific interest but also for their numerous industrial applications.Therefore, several experimental techniques have been used to investigate their properties [1][2][3][4][5][6], in particular the surface tension [7][8][9][10][11][12].Although, most of the studies are conducted for one type of surfactant, a lot of commercial products consist of a mixture of molecules which present richer properties.For example, anionic surfactants are usually mixed with non-anionic ones in products such as shampoos, dish washing liquids, washing powders, whereas mixtures of cationic with non-ionic surfactants are good for disinfectants, cleaners and sanitizers.On the other hand, the capability of surfactants to be adsorbed at liquid interfaces is an important property and it is widely used for many technological processes, such as detergents, foam, and emulsion stabilizers.In fact, adsorption of surfactant solutions at liquid interfaces has been investigated in many papers using mixtures of non-ionic with anionic or cationic surfactants [13][14][15].
In particular, synergistic effects of mixtures of anionic and cationic surfactants are relevant in several industrial applications.In experiments, the synergy can be determined by measuring the surface tension as a function of concentration at different surfactant ratios, which is time consuming.Therefore, the combination of computer methods and experiments will reduce the experimental times and provide microscopic insight into the mechanism of the synergistic effect.
Since computer simulations are relevant to understanding the interfacial problems, they have appeared as an alternative to the study of such complex systems.From the computer simulations it is possible to obtain more data about thermodynamic properties, such as the surface tension, that are not easy to collect from actual experiments.However, to achieve this goal, it is imperative to have good computer methods to gain reliable results, i.e., not only the force field is important but also the conditions to carry out the simulations.For instance, the calculation of the surface tension can depend on the number of molecules and the cutoff values used in the simulations [16][17][18][19].Therefore, some papers have investigated the effects of truncated potentials [16], the inclusion of long-range corrections [20], the finite size and the periodic boundary conditions in Lennard-Jones fluids at the liquid/vapour interface.Moreover, some authors have shown that the surface tension data present an oscillatory behaviour [20] and they concluded that reasonable values of the surface tension can be obtained with large interfacial areas [20] and cutoff values of about 2.5 nm [21].
In the present work we extend previous studies and calculate the surface tension of systems composed of monolayers with one type and mixtures of surfactant molecules at the liquid/vapour interface.We show that choosing incorrect simulations boxes can give erroneous surface tension results, including values that are not physically acceptable as negative values.The errors can be corrected if boxes with large interfacial areas are taken.Simulations were carried out by using different ionic and nonionic surfactants.

Computational model
Different surfactant molecules were simulated, the sodium dodecyl sulfate (SDS), the dodecyltrimethylammonium bromide (DTAB), the cetyltrimethylammonium bromide (CTAB) and the zwitterionic lauryl sulfobetaine (SB-3), called Betaine in figure 1.The first one is anionic, the next two are cationic and the last one is nonionic surfactant, respectively (see figure 1).In all cases the atoms of the surfactant headgroups were explicitly simulated whereas the CH  groups in the hydrocarbon tails were modelled using the united atom model (see figure 1).Counterions, Na + and Br − , were also introduced into the water phase to neutralize the SDS and DTAB (and CTAB) systems.For water, the TIP4P/ model was used [22] and all the unlike interactions were handled with the Lorentz-Berthelot rules [23].
Initially, simulations of monolayers with one type of surfactant were carried out to validate the force field parameters of all the surfactants.For these simulations four monolayers were prepared, each one with 45 SDS, 45 DTAB, 45 CTAB and 45 Betaine molecules, which were placed at the liquid/vapour interface in a rectangular box with dimensions   =   = 5 nm and   = 15 nm.Then, the surface tensions were calculated and compared with the experimental data.In some cases, the Lennard-Jones (LJ) parameters of the surfactants were reparametrized to obtain a better agreement between the simulations and the experiments (see the next section).With the LJ parameters obtained from previous calculations other series of simulations were executed for monolayers prepared with mixtures of CTAB/SDS, CTAB/DTAB and CTAB/Betaine surfactants at different compositions.
All simulations were carried out using the Gromacs software [24] in the canonical (NVT) ensemble at constant temperature,  = 298 K, using the Nosé-Hoover thermostat [25] with a relaxation time constant of   = 1 ps.Electrostatic interactions were handled with the particle mesh Ewald method [26] and the short range interactions were cutoff at 2.5 nm as suggested in reference [21].Bond lengths were constrained using the Lincs algorithm [27].Then, simulations were carried out up to 50 ns after 10 ns of equilibration with a time step of d = 0.002 ps, and collecting data, over the last 20 ns, were taken for data analysis.Equilibration of the systems was monitored with the configurational energy and the calculated structural properties (see the results section), which did not change significantly along the simulation time.

Monolayers with a single type of surfactant
As stated above, initial simulations of monolayers with a single type of surfactant were performed to reparametrize the surfactant LJ parameters and obtain reasonable values of the surface tensions.Reparametrization was carried out by the three step systematic parametrization procedure, 3SSPP [28], where all the  and  LJ parameters were scaled from their original values until the surface tension measurements were in good agreement with the experiments, typically within 5% error.The procedure was conducted for the CTAB and SDS monolayers whereas for the Betaine monolayer the parameters were taken from reference [29].In the molecular dynamics, the surface tension was calculated using the expression, where   is the length of the simulation box and   are the components of the pressure tensor,   = (1/) (  −   ) with   the kinetic energy,  the volume and    the virial expression [24], r ij is the position vector between atom  and atom , and F ij is the force vector between those atoms.Equation (3.1) is appropriate to calculate the surface tension with "planar" interfaces, i.e., without pronounced curvatures.In our case we have two interfaces, water/air at one side of the box and water/surfactant/air at the other.Then, from the gromacs software the surface tension is obtained,  =  water/air +  water/surfactant/air . (3.3) By knowing the  water/air , the  water/surfactant/air can be estimated.Therefore, simulations for the water/air interface were carried out and a value of 69 mN/m was found, in good agreement with the experimental data (72 mN/m).The new surfactant Lennard-Jones parameters and the surface tensions for the monolayers with single types of surfactants are given in tables 1-4.

Monolayers with surfactant mixtures
With the current surfactant LJ parameters, three monolayers were prepared with mixtures of CTAB/SDS, CTAB/DTAB and CTAB/Betaine at distinct compositions, 75/25, 50/50 and 25/75.Initial simulations were conducted with the same box dimensions as described in the previous section, i.e.,   =   = 5 nm and   = 15 nm using an area per molecule of 26 Å 2 /molecule for the CTAB/SDS and 40 Å 2 /molecule for the CTAB/DTAB and CTAB/Betaine mixtures.The area per molecule is defined as the interfacial area divided by the number of surfactants (  ×   / Num.surfactants).With these simulation boxes the surface tension was calculated, although at these conditions the results lead to large errors.Even negative values were found for the asymmetric compositions, which have no physical meaning (figure 2a).Exploring possible sources of errors, from the snapshots of the last configuration of the monolayers at small interfacial areas it is observed that they are distorted, i.e., they are not planar (top of figure 3).Therefore, equation (3.1) cannot be used, or in other words its application yields the not-physical results.Interesting are the surfactant molecules in the bulk water phase forming micelle-like structures for some surfactant compositions (bottom of figure 3).We analyze the size effects on the surface tension data by increasing the interfacial area of the box in different simulations.Figure 2 shows the values of the surface tension as function of the box length, by keeping the area per molecule constant.The simulations were carried out for all the mixtures, for different compositions, and in all cases similar trends were observed, i.e., the surface tension fluctuates for small interfacial areas and it goes to constant values at large areas, for box lengths above   =   = 120 Å.It is worth mentioning that in these simulations, the number of surfactant molecules changes (to keep the area per molecule fixed), and also changes the number of water molecules to have a sufficient bulk water phase, at least 3 nm layer thick in the -direction.We also imposed   = 3  and a cutoff radius of 2.5 nm, to have reliable surface tension data as suggested in previous works [21].Simulations were run up to 80 ns and similar behaviour was observed, the monolayers lost their planar structure for small box lengths (  and   ).At large areas, they are flat enough to consider that the equation (3.1) can be used correctly.Moreover, for large box lengths, the components of the pressure tensor,    and   have similar values as expected.It is important to mention that in the simulations, the pressure tensor components show large fluctuations, although the average value remains constant (figure 4).
In order to better understand the behaviour of the monolayers at the interface, in the next sections we study only the CTAB/SDS system more in detail, similar tendencies were observed for all the

Density profiles for the monolayers of surfactant mixtures
Mass density profiles were calculated to determine the structure of the monolayers at the interface, defined as where  () is the mass of the calculated species in the volume   ×   × Δ.Since we are only interested in the monolayers at the interface, the analysis of the CTAB/SDS mixture at the 25/75 and 75/25 compositions, was done without considering the micelles formed in the water phase at large areas (see bottom of figure 3).Charge distribution profiles can also give us information about the distribution of the monolayers at the water/air interface.These profiles are calculated as, where () is the total charge in the volume   ×   × Δ.The studies of the profiles were also carried out for the monolayers simulated with small and large interfacial areas.In figure 8 it is possible to observe that the charge density profiles are nearly zero, far from the interfaces (blue lines), for the systems simulated with a large area (  =   = 120 Å).For the monolayers simulated with the small area (  =   = 60 Å), an excess of charge in the middle of the systems (red lines) is noted, most likely caused by the surfactants headgroups.This last result suggests that not all surfactants are located at the interface, i.e., planar monolayers are not formed.

Order of the surfactant tails and tilt angle
Additional information about the structure of the monolayers with small and large box areas can be obtained from an order parameter which tells us how the surfactants are arranged within them.In experiments, the ordering of the surfactant hydrocarbon tails is characterized by the so-called deuterium order parameter,   , related to the average inclination of the C-H bond, in the CH  groups, with respect to the surface normal.However, in computer simulations the order parameter is calculated, in the united atom model, by the following equation [30], 13605-8 Charge Density (e nm with, where ,  = , ,  and   is the angle between the -th molecular axis and the normal to the interface.
In this work we calculated the   order parameter.It is worth mentioning that   = −0.5 correspond to complete order parallel to the interface whereas   = 1.0 is complete order in the direction normal to the interface.However, here it is more convenient to use the average of the order parameter, ⟨  ⟩, to study the order in the tails [38].Then, the analysis of the order parameter was carried out over the carbons that are at the water/air interface.In table 5 the ⟨  ⟩ values for the first four carbons are given for both surfactants, CTAB and SDS, for a small and for a large interfacial area.It is noted that the tails are more ordered for large areas than for small ones, suggesting that the surfactants accommodate better in the monolayer at the large interfacial area.
The structure of the monolayers can also be analyzed by the tilt angle of the surfactant chains.The angle is measured as,   =   /  , where   is the average projection of the chains along the normal to the interface and   is the total length of the tails (from the last to the first carbon atom).The results are  shown in table 6.It is observed that the tilt angle in general is a bit lower for the mixtures in the small area than in the large one, regardless of the surfactant, i.e., the tails are straighter with respect to the interface at the large box area allowing better arrangement of the surfactants in the monolayers, in agreement with the order parameter ⟨  ⟩ results.The arrangement between the different surfactants inside the monolayer was also studied in terms of pair distribution functions, ().The () were calculated between the headgroups of each surfactant and the results are given in figure 9.For all the CTAB/SDS compositions, the peaks in the () are higher for the systems with the large area (length of   =   = 120 Å) than for the small area (  =   = 60 Å).These results suggest that there are more surfactants close to each other in the monolayers simulated with large areas than with small ones and consequently, the monolayers are more "compact" in the first case.
One last analysis to compare the monolayers for small and large interfacial areas can be given in terms of the configuration energies, calculated over the monolayers only, and the results are shown in table 5. Energies are negatively greater for the monolayers with large areas than for the small areas, indicating 13605-10 that the systems in large areas are more stable.The values of energies are divided by the total number of atoms in the monolayer.

Conclusions
The mixtures of surfactants are of great relevance to investigate in several industrial and technological applications (foaming, detergency, etc.), where the surface tension is used as a parameter to determine their synergistic effects [39,40].In the present paper, we carried out molecular dynamics of monolayers with ionic and non-ionic surfactants to study the surface tension at the water/air interface.We pointed out the importance of choosing right simulation boxes to obtain reliable surface tension values of surfactant monolayers.Unsuitable conditions, such as calculations with small interfacial areas, can lead to artificial effects in the surface tension values.For instance, fluctuations in the surface tension data are observed due to the periodic boundary conditions and the finite box size.In the case of simulations of surfactant mixtures with small areas, the surface tension calculations can give non-physical values such as negative ones.It was found that reasonable values of the surface tension are obtained when the simulations are conducted with large box areas (box lengths of above 120 Å), with cutoffs of 2.5 nm and using a sufficient number of water molecules to form bulk phases at least 30 Å thick.
Several interfacial problems, e.g., the description of the solid-liquid interface requires the knowledge of the surface tension.Then, the wetting properties of such complex interfaces can be studied.Therefore, we expect that our results may be useful and stimulate theoretical developments along lines of research discussed in [41][42][43].

Figure 1 .
Figure 1.(Colour online) Molecular structure for each surfactant used in the simulations within the united atom model.

Figure 2 .
Figure 2. (Colour online) Surface tension data for the monolayers with mixtures of a) CTAB/SDS, b) CTAB/DTAB and CTAB/Betaine at different surfactant compositions as function of the simulation box lengths.

Figure 3 .
Figure 3. (Colour online) Snapshots of the last configurations of the monolayers with mixtures of CTAB/SDS surfactants at different compositions, at different interfacial areas of the simulation box ( =   ×   ).Top figures,   =   = 60 Å.Bottom figures   =   = 120 Å. Pink color for the CTAB and blue for the SDS surfactants.Counterions are represented by dots and water is removed for better visualization of the monolayers.

Figure 4 .
Figure 4. (Colour online) Components of the pressure tensor as function of the simulation time for the monolayer with the CTAB/SDS mixture, a) P   , b) P  , c) P  .

Figure 5 .
Figure 5. (Colour online) density profiles for the CTAB/SDS monolayer at the 25/75 composition.Top: small box area.Bottom: large box area.

Figure 6 .
Figure 6.(Colour online) Mass density profiles for the CTAB/SDS monolayer at the 50/50 composition.Top: small box area.Bottom: large box area.

Figure 7 .
Figure 7. (Colour online) Mass density profiles for the CTAB/SDS monolayer at the 75/25 composition.Top: small box area.Bottom: large box area.

Table 4 .
Surface tension data of all the surfactants, experiments, calculated in this work.

Table 5 .
Configurational energy and the average order parameter, ⟨  ⟩, for the CTAB/SDS mixture at