Surface Excess Free Energy: An Elaboration with Particular Insight for Use as a Predictor of Solvophilicity in Molecular Simulation

This manuscript details the method to determine the surface excess from readily derivable ensemble properties, namely the pressure tensor, via computational molecular dynamics. It will then expand upon the theoretical and practical uses of quantities in Gibbs-Duhem like relationships for the surface excess and molecular concentration at the interface. Furthermore, it details several limitations of computational molecular dynamics, mainly to determine force field parameters natively and also to determine criteria for switching the bond order at certain temperatures. The goal in predicting surface presence is in inter-relating the relative surface excess free energies of each species with respect to the total system relative to the free energy of hydration of that system. *Corresponding author: Mongelli GF, Department of Chemical Engineering, Case Western Reserve University10900 Euclid Ave, AW Smith 116, Cleveland, OH 44120, USA, Tel: 12163682000; E-mail: Gfm12@case.edu Received December 14, 2016; Accepted December 21, 2016; Published December 31, 2016 Citation: Mongelli GF (2016) Surface Excess Free Energy: An Elaboration with Particular Insight for Use as a Predictor of Solvophilicity in Molecular Simulation. J Material Sci Eng 6: 310. doi: 10.4172/2169-0022.1000310 Copyright: © 2016 Mongelli GF. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Introduction
While the fundamental thermodynamic equations of state have been well documented, particularly relations for entropy and free energy of hydration for various systems, such relationships have not been explored with the use computational molecular dynamics. This is especially the case for parameters of unique interest to polymeric systems within mixed solvents. Such systems are a growing and multibillion dollar global annual industry [1].
Two methods will be discussed to determine the free energy of hydration or solvation, and their relative complexity will be compared. The equations of relevance for calculation of the surface excess free energy will be explored in detail within this manuscript. Which equations are best suited for experimental or theoretical determination will be discussed. The surface excess free energy will be defined. How to compute this important quantity with computational molecular dynamics software will be elucidated. Additionally, some potential explanations for which aspects of molecular structure lead to reduced surface tension will be discussed.

Potentially Simpler Free Energy of Hydration Computation
An alternative method will employ the pull code within molecular dynamics simulations to obtain the free energy landscape for the surfactant molecule with respect to both the vacuum gas state, the bulk solubilized, and the interface. The pulling constant can be altered to perturb the polymer or target molecule of interest to various z-space values, where the interfaces to the various states are at z-axis parallels. In such a manner, the intermolecular forces present in different z-slabs can be determined [2]. As molecular structuring around the interface has been known to be different from the bulk, the density profile can be mapped onto the free energy of hydration in cases where the pulling constant is such that it allows the polymer or target molecule to swing throughout all z-space. Such density profiles can be proportional probability of a polymer or target molecule taking a certain position z in with respect to the interface is given by the following relation, and therefore indicates the relative free energies of the system: where C is a constant and Prob[z] comes from umbrella sampled systems. The free energy landscape, G(z), is thus obtained by inverting this equation, The values of Prob(z) are obtained as histograms constructed from a molecular dynamics trajectory. With the energy present at k BT , umbrella sampling may be used to force the system into relatively unstable states. The free energy can then be related to the un-perturbed systems via Li-Makarov analysis [3]. This is required at low alcohol contents in co-solvated systems, where the intermolecular forces and free energies of solvation are of energies with magnitude greater than kB*T. To obtain better statistics, umbrella sampling will be used to bias the surfactant to certain positions with respect to the liquid surface; this bias is removed in the free energy calculation using the self-consistent histogram method [4,5]. Several simulations may need to be performed to start the system in the high energy state, i.e., solvated for insoluble polymers and in a surface state for soluble polymers. Soluble polymers might have even higher energy configurations when placed more than the cut-off radius away from the solvent interface in the vacuum region.
By establishing the differences between the relative free energy relationships, it becomes very computationally cheap to predict the solubility and surface activity properties of new molecules and new solvents. A single simulation of a polymer studied in this work with a new solvent allows for the prediction of the surface parameter of all of the molecules within such work in the studied solvents. Correspondingly, a simulation of a new solute in a previously studied solvent allows for the prediction of the surface parameter in each of the studied solvents. This is determinable because each solvent has a one-to-one mapping of free energy of hydration relative to surface parameter.
After mapping out the free energy changes between solvents, and between polymers, we can do a single simulation and predict the surface parameter in all other contents for a particular solvents and indeed different solvents as well. The argument for exploring ΔG[z] beyond the cutoff radius in the 100% alcohol case allows us to estimate when polymers will go into the bulk solution that have less repulsion from the solvents than the least repulsive solute in a series. Additionally, we will only get all of the ΔG[z] information we need when the bulk state parameter goes to one and not non-zero --as stated earlier. This would imply that the surface parameter would need to go to zero in the 100% alcohol limit, which we have not seen for many systems. That being said, we need the free energy increases as a function of z beyond those estimated by these materials in the low alcohol content limit to predict when more soluble materials will go in and at what critical free energy values.
The free energy of hydration may be more simply determined by treatment of its time averaged quantity as a probability distribution in NVT simulations. The following functions are written for the isothermal-isobaric ensemble, but translate readily into the NVT/ Canonical ensemble by replacing Gibb's Free energy, G[z], with Helmholtz Free Energy, A[z]. Then: Inverting to free energy space: Therefore, relative free energy relationships can be determined between NVT systems. It is not immediately clear how to compare the relative free energy of solvation between isothermic-isobaric simulations and NVS/Canonical simulations. In considering entropy changes from the simulation results, derivations are readily possible from general statistical mechanics theory. In the NPT and NVT systems, for bulk non-interface systems without a vacuum: In the case of NPT systems, the Helmholtz Free Energy is the negative of the Legendre transform with respect to entropy.
This entropy change implies a change in the number of configurations available to the system. The number of configurations available to the system can be compared in ideal gas states, interface active, bulk solvated states. The simple relation required to go from entropy, which can be established from (5) or (6) via an ensemble potential, is simply Boltzmann's Entropy formula:

Surface Excess
One quantity of interest which is available from these results via a derivation and many computational molecular dynamicists have not computed is the surface excess and associated Gibb's isotherms of hydrophobic and hydrophilic molecules. Such a surface excess [6] tells the area at the interface per molecule, α, as described by: N is Avogadro's number. Γ is the surface excess in units of moles/ m 2 . This quantity offers key information regarding how the surface tension, surface parameter, and chemical potentials of the solvent and the polymer vary together. Maximizing surface tension changes with minimum chemical potential changes represents an improved surfactantability, since the same surface tension can be achieved with reduced polymer concentrations. That is what should be observed from computing the other information within this section.
The following equations describe surface phenomena in the Gibb's surface formalism [7]: This tells how the surface tension will be impacted as a function of the surface excess quantities of various molecules and their respective chemical potentials. It is useful in estimation of whether molecules will increase or decrease the surface tension.
The above relation [8] indicates how the surface excess is calculated as a function of the number of molecules of a specific type per unit surface area. It is simple to study more complex, multi-component systems by increasing the number of types of molecules detracting from n i TOTAL in the numerator. It is most useful for determining the surface excess via a single computational molecular dynamics simulation. Although the author is not aware of any codes presently which do this computation natively in one step, the number of molecules of each species within ca. one nanometer -or the cutoff length of non-bonded computations --of the interface should be determined, and the relative number of molecules will allow for the computation of the surface excess via the above described method.
The above relation is most useful for the determination of the surface excess from experimental methods. The alcohol or polymeric content is varied and the surface tension is measured. The surface tension decrease will be linear up to a breakdown point, at which additional concentration increases do not result in surface tension decreases. For immiscible materials, this means that the added molecules are no longer being added to the surface, but instead take bulk states in the preferred phase, and do not impact the surface tension. This holds true for polymers. For miscible materials, the free energy of solvation changes for each species as a function of the relative mass contents of the mixture. Therefore, adding more content of one species or another will change the free energies of solvation for either species and the energies of moving any individual species from the bulk to the interface. This relation also is useful for determining the surface excess from comparing the surface tensions of multiple simulations.
If the surface excess of one species is known in a binary mixture, then the surface excess of the second member of the mixture may be determined from an experiment or simulation in a different content -assuming that the surface tension decrease has not occurred or is negligible.  The full NPT/isothermal-isobaric ensemble theory is written above, including terms for surface tension and area change effects. Of increased interest in this theory is the NVT/Canonical ensemble derivation, which utilizes the canonical partition function, Z, as defined earlier in an exponential of the Helmholtz free energy relative to the inverse temperature k B T, rather than free energy.
For a system with a constant number of molecules and a constant interfacial area, as the alcohol content and alcohol structure varies, the area of the van der Walls surface varies as does the chemical potential in the vicinity of the interface. The surface excess refers to the energy per unit area that the molecule contributes to the interfacial chemical potential. The surface excess then tells how many molecules are required at an interface to alter the surface tension to a particular value and how many are required to be replaced to maintain a constant chemical potential. The solubility and intermolecular forces dictate whether such replacements actually occur.
The surface excess and its relationship to predicting the surface state presence of a molecule has an analogous relationship with fugacity. The surface excess tells of the surface presence of molecules spatial position in the system relative to the free energy of hydration whereas the fugacity predicts the vapor pressure and phase change fraction relative to the enthalpy of vaporization.
It is expected that there exists a relation similar and analogous to that of chemical potential, surface excess and surface tension which interrelates the volumetric hydrogen bond density and the free energy of hydration. Furthermore, a program should be written to will address the spatial correlation of the hydrogen bonds to classify them as within the cutoff radius of the interface and beyond the cutoff radius of the interface. These results can be related to the surface partitioning results discussed above, and will be used to assess the explanation of polyether solubility trends [5].

Possible Relation between Molecular Parameterizations and Surface Tension Reduction Capacity
The presence of a particular molecule in a mixture can only affect the surface tension if said solute is insoluble in the mixed solvent. Therefore, the foaming capability of insoluble polyethers is lost when adding alcohol co-solvents due to the loss of low surface tensionassociated molecules from occupying the interface. The loss of foaming capability is not from alcohol molecules displacing polymer molecules from the interfacial surface caused by their own solvophobic repulsions. Since the surface tension of a mixture is dependent on the surface tensions of the individual components and what fraction of the surface they occupy, the addition of alcohols have low surface tension would decrease the overall mixture surface tension.
What is left to determine is which types of forces within the polymer or target molecule system result in the surface tension reduction relative to their non-polymer-or non-target-molecule-containing co-solvents. This phenomenon is likely the result of a needle or any other surface interacting body pushing against the dihedral interaction energies along the polymer backbone. Once a certain critical energy state is reached, the Lennard-Jones and Coulombic interactions combined with the physical external force of the body overcome the dihedral energy barrier. This is a sort of rotational energy problem. If you want to spin a water wheel so you can put grain through a mill or generate electricity but the wheel has some trash stuck in it and water isn't strong enough to fully rotate the wheel. If you give it a good push, then you will get the trash up to the top of the wheel and then the wheel will rotate down until it gets stuck again. In this case, the trash is sterically hindered groups on molecules, which actually increases the rotation of the wheel after it hits the top of the wheel, and the water flow is k BT .

Conclusions
With new rigorous computational methods, such as those discussed in this work, it becomes easier to determine the free energies of solvation. We can receive results faster or more cheaply than with experiments in a large number case in the present force field. Faster, in the sense that doing many simulations scales fiscally much more readily than does doing many experiments. At a point beyond that, it may be determined that some presently non-parameterized atoms have such force field embodiments and whose associated to-be-determinedproperties are desired. Though it is much more desired that molecules associated with certainly real force field parameterizations are studies via simulation when they cannot be synthesized above and before purely theoretically molecules which result from perturbed force field parameterizations of known molecules. Again, in the limit of more computational power than scholars can come up with short-term coding applications for -which may be the case given the acceleration of Moore's Law-and that complex coding and theory mathematics to run on supercomputers is limited by human derivations, perhaps the latter is desirable.
A possible explanation for the surface tension reducing capability of a molecule is suggested. One potential method to test for this explanation is to dial up the dihedral interaction constants for a surface active molecule and see if it will change the ensemble average of the pressure tensor and associated surface tension. Perhaps a more ambitious goal, then, is to determine which classes of structures could be surface active and maintain foaming capability from molecular simulation.

Future Work
Once polyethers, polysilicones and their polyalkane analogue materials can be simulated each of the quantities of interest can be determined via molecular simulation and from experimental studies. Any discrepancies in their derivations from these two methods may lead to changes in computational molecular dynamics theory for polymeric materials or, rather, to changes in the molecular parameterization constants for these materials within the Optimized Potentials for Liquid Simulations-All Atoms (OPLS-AA) standard.