Surfactants are amphiphilic molecules with multiple uses and industrial applications as detergents, wetting agents, emulsifiers, and so forth. They can be divided into three main categories: nonionic, ionic, and zwitterionic. The development of a universal computational framework able to predict key properties such as their critical micelle concentration (cmc) and the size of the micelles they form and to ultimately extract phase diagrams for their aqueous solutions, possibly in the presence of salts and oils, using their chemical constitution as input, would provide valuable information for the design and the production of these materials. In this work, we focus on ionic surfactants and investigate a possible route toward the development of such a framework based on coarse-grained simulations using the MARTINI forcefield in two versions: its implicit solvent version, called Dry MARTINI, and its explicit solvent version, called Wet MARTINI. The surfactants considered in our efforts are the anionic sodium dodecyl sulfate (SDS) and the three cationic cetyl, dodecyl, and octyl trimethyl ammonium bromide (CTAB, DTAB, and OTAB, respectively). First, we choose their mapping onto coarse-grained MARTINI beads. Next, we estimate their cmc’s, their peak aggregation numbers, N agg, and in the case of SDS, its small angle neutron scattering pattern at low concentrations, applying the Dry MARTINI forcefield. With a single modification to the Lennard-Jones energy parameter between hydrophobic beads and invoking Ewald summation with a physically meaningful dielectric constant for electrostatic interactions, our estimates are in very good agreement with experimental results. Furthermore, we predict the phase behavior of SDS/water and CTAB/water binary solutions using Wet MARTINI and find semiquantitative agreement with experimental phase diagrams. We conclude that the MARTINI forcefield, with careful treatment of electrostatic interactions and appropriate modification of parameters for some key functional groups, can be a powerful ally in the quest for a universal computational framework for the design of new surfactants with improved properties.

A complete thermodynamic analysis of mixtures consisting of molecules with complex chemical constitution can be rather demanding. The Kirkwood–Buff theory of solutions allows the estimation of thermodynamic properties, which cannot be directly extracted from atomistic simulations, such as the Gibbs energy of mixing (Δmix G). In this work, we perform molecular dynamics simulations of n-hexane–ethanol binary mixtures in the liquid state under two temperature–pressure conditions and at various mole fractions. On the basis of the recently published methodology of Galata Galata Fluid Phase Equilib.20184702537, we first calculate the Kirkwood–Buff integrals in the isothermal–isobaric (NpT) ensemble, identifying how system size affects their estimation. We then extract the activity coefficients, excess Gibbs energy, excess enthalpy, and excess entropy for the n-hexane–ethanol binary mixtures we simulate. We employ two approaches for quantifying composition fluctuations: one based on counting molecular centers of mass and a second one based on counting molecular segments. Results from the two approaches are practically indistinguishable. We compare our results against predictions of vapor–liquid equilibria obtained in a previous simulation work using the same force field, as well as with experimental data, and find very good agreement. In addition, we develop a simple methodology to identify the hydrogen bonds between ethanol molecules and analyze their effects on mixing properties.

We introduce a physics-based model for calculating partition coefficients of solutes between water and alkanes, using a combination of a semi-empirical method for COSMO charge density calculation and statistical sampling of internal hydrogen bonds (IHBs). We validate the model on the experimental partition data (∼3500 molecules) of small organics, drug-like molecules, and statistical assessment of modeling of proteins and ligand drugs. The model combines two novel algorithms: a bond-correction method for improving the calculation of COSMO charge density from AM1 calculations and a sampling method to deal with IHBs. From a comparison of simulated and experimental partition coefficients, we find a root-mean-square deviation of roughly one log 10 unit. From IHB analysis, we know that IHBs can be present in two states: open (in water) and closed (in apolar solvent). The difference can lead to a shift of as much as two log 10 units per IHB; not taking this effect into account can lead to substantial errors. The method takes a few minutes of calculation time on a single core, per molecule. Although this is still much slower than quantitative structure–activity relationship, it is much faster than molecular simulations and can be readily incorporated into any screening method.