Application of computer simulation free-energy methods to compute the free energy of micellization as a function of micelle composition. 2. Implementation

In this paper, the implementation of the CS-FE/MT model introduced in article 1 is discussed, and computer simulations are performed to evaluate the feasibility of the new theoretical approach. As discussed in article 1, making predictions of surfactant/solubilizate aqueous solution behavior using the CS-FE/MT model requires evaluation of DeltaDeltaG for multiple surfactant-to-solubilizate or surfactant-to-cosurfactant transformations. The central goal of this article is to evaluate the quantitative accuracy of the alchemical computer simulation method used in the CS-FE/MT modeling approach to predict DeltaDeltaG for a single surfactant-to-solubilizate or for a single surfactant-to-cosurfactant transformation. A hybrid single/dual topology approach was used to morph the ionic surfactant sodium dodecyl sulfate (SDS) into the ionic solubilizate ibuprofen (IBU), and a dual topology approach was used to morph the nonionic surfactant octyl glucoside (OG) into the nonionic solubilizate p-aminobenzoate (PAB). In addition, a single topology approach was used to morph the nonionic surfactant n-decyl dimethylphosphine oxide (C10PO) into the nonionic cosurfactant n-decyl methyl sulfoxide (C10SO), the nonionic surfactant octylsulfinyl ethanol (C8SE) into the nonionic cosurfactant decylsulfinyl ethanol (C10SE), and the nonionic surfactant n-decyl methyl sulfoxide (C10SO) into the nonionic cosurfactant n-octyl methyl sulfoxide (C8SO). Each DeltaDeltaG value was computed by using thermodynamic integration to determine the difference in free energy associated with (i) transforming a surfactant molecule of type A into a cosurfactant/solubilizate molecule of type B in a micellar environment (referred to as DeltaG2) and (ii) transforming a surfactant molecule of type A into a cosurfactant/solubilizate molecule of type B in aqueous solution (referred to as DeltaG1). CS-FE/MT model predictions of DeltaDeltaG for each alchemical transformation were made at a number of simulation conditions, including (i) different equilibration times at each value of the coupling parameter lambda, (ii) different data-gathering times at each lambda value, and (iii) simulation at a different number of lambda values. For the three surfactant-to-cosurfactant transformations considered here, the DeltaDeltaG values predicted by the CS-FE/MT model were compared with DeltaDeltaG values predicted by an accurate molecular thermodynamic (MT) model developed by fitting to experimental CMC data. Even after performing lengthy equilibration and data gathering at each lambda value, physically unrealistic values of DeltaDeltaG were predicted by the CS-FE/MT model for the transformations of SDS into IBU and of OG into PAB. However, more physically realistic DeltaDeltaG values were predicted for the transformation of C10PO into C10SO, and reasonable free-energy predictions were obtained for the transformations of C8SE into C10SE and C10SO into C8SO. Each of the surfactant-to-cosurfactant transformations considered here involved less extensive structural changes than the surfactant-to-solubilizate transformations. As computer power increases and as improvements are made to alchemical free-energy methods, it may become possible to apply the CS-FE/MT model to make accurate predictions of the free-energy changes associated with forming multicomponent surfactant and solubilizate micelles in aqueous solution where the chemical structures of the surfactants, cosurfactants, and solubilizates differ significantly.     

Comparison of two simulation methods to compute solvation free energies and partition coefficients

The thermodynamic integration (TI) and expanded ensemble (EE) methods are used here to calculate the hydration free energy in water, the solvation free energy in 1‐octanol, and the octanol‐water partition coefficient for a six compounds of varying functionality using the optimized potentials for liquid simulations (OPLS) all‐atom (AA) force field parameters and atomic charges. Both methods use the molecular dynamics algorithm as a primary component of the simulation protocol, and both have found wide applications in fields such as the calculation of activity coefficients, phase behavior, and partition coefficients. Both methods result in solvation free energies and 1‐octanol/water partition coefficients with average absolute deviations (AAD) from experimental data to within 4 kJ/mol and 0.5 log units, respectively. Here, we find that in simulations the OPLS‐AA force field parameters (with fixed charges) can reproduce solvation free energies of solutes in 1‐octanol with AAD of about half that for the solute hydration free energies using a extended simple point charge (SPC/E) model of water. The computational efficiency of the two simulation methods are compared based on the time (in nanoseconds) required to obtain similar standard deviations in the solvation free energies and 1‐octanol/water partition coefficients. By this analysis, the EE method is found to be a factor of nine more efficient than the TI algorithm. For both methods, solvation free energy calculations in 1‐octanol consume roughly an order of magnitude more CPU hours than the hydration free energy calculations. © 2012 Wiley Periodicals, Inc.     

Calculation of the hydration free energy difference between pyridine and its methyl-substituted derivatives by computer simulation methods

The difference of the hydration free energy of pyridine and its methyl- and symmetrically dimethyl-substituted derivatives has been calculated by the method of free energy perturbation. To check the precision of the results obtained, we have repeated the calculations using thermodynamic integration over different paths. Besides the hydration free energy, the difference in the energy and entropy of hydration between pyridine and monomethyl- and dimethylpyridines has also been determined. The obtained results clearly show that the hydration free energy of the pyridine derivatives becomes more negative with each additional methyl group. However, the accuracy of the calculation does not allow us to draw any conclusion about the dependence of the hydration free energy on the location of the methyl group. The analysis of the Coulomb and Lennard-Jones contributions to the hydration free energy differences has shown the dominance of the latter term. The comparison of the hydration energy and free energy values has shown that there is a strong compensation effect between the energetic and entropic terms of the free energy. The hydration energy of the solute becomes considerably more negative with each additional methyl group due to the dispersion attraction between the methyl group and the surrounding water molecules. The introduction of a methyl group results in an approximately 30 J/(mol K) decrease of the entropy of hydration, and hence, at 300 K, the entropic contribution to the hydration free energy increases by about 9 kJ/mol. Due to their opposite signs, the entropic and energetic contributions largely cancel each other, resulting in approximately an order of magnitude smaller value for the free energy.     

On the use of computer simulation to determine the excess free energy in fluid mixtures

Haile, J.M., 1986. On the use of computer simulation to determine the excess free energy in fluid mixtures. Fluid Phase Equilibria, 26: 103–127This paper first reviews the use of Kirkwood's coupling parameter for determining residual chemical potentials, activity coefficients, and Henry's constants from Monte Carlo and molecular dynamics computer simulations. A new version of the method is then developed for obtaining the excess Gibbs free energy from isothermal-isobaric simulations. New expressions are also given for the excess volume and excess entropy.The revised method is demonstrated by computing excess free energies for 14 mixtures of repulsive soft spheres. Isothermal-isobaric molecular dynamics was used to generate the necessary simulation data. Although the excess free energies for these particular mixtures are small in magnitude (|G^E/NkT| < 0.1), the simulation method generally gives G^E within 5% of the values calculated by thermodynamic perturbation theory.     

Thermodynamic cycle integration by computer simulation as a tool for obtaining free energy differences in molecular chemistry

Summary A new and promising development in the field of computer simulation of molecular systems is the socalled thermodynamic cycle integration technique, which combines well-known results from statistical thermodynamics with powerful computer simulation methods. The basic formulas, the development and the applications in the areas of drug design, protein engineering and conformational analysis of this elegant technique are discussed.     

An approximate but efficient method to calculate free energy trends by computer simulation: Application to dihydrofolate reductase-inhibitor complexes

Summary Derivatives of free energy differences have been calculated by molecular dynamics techniques. The systems under study were ternary complexes of Trimethoprim (TMP) with dihydrofolate reductases of E. coli and chicken liver, containing the cofactor NADPH. Derivatives are taken with respect to modification of TMP, with emphasis on altering the 3-, 4- and 5-substituents of the phenyl ring. A linear approximation allows the encompassing of a whole set of modifications in a single simulation, as opposed to a full perturbation calculation, which requires a separate simulation for each modification. In the case considered here, the proposed technique requires a factor of 1000 less computing effort than a full free energy perturbation calculation. For the linear approximation to yield a significant result, one has to find ways of choosing the perturbation evolution, such that the initial trend mirrors the full calculation. The generation of new atoms requires a careful treatment of the singular terms in the non-bonded interaction. The result can be represented by maps of the changed molecule, which indicate whether complex formation is favoured under movement of partial charges and change in atom polarizabilities. Comparison with experimental measurements of inhibition constants reveals fair agreement in the range of values covered. However, detailed comparison fails to show a significant correlation. Possible reasons for the most pronounced deviations are given.     

Direct calculation of the crystal-melt interfacial free energy via molecular dynamics computer simulation

We review our recent work on the direct calculation of the interfacial free energy, gamma, of the crystal-melt interface via molecular dynamics computer simulation for a number of model systems. The value of gamma as a function of crystal orientation is determined using a thermodynamic integration technique employing moving cleaving walls [Phys. Rev. Lett. 2000, 85, 4751]. The calculation is sufficiently precise to resolve the small anisotropy in gamma, which is crucial in determining the kinetics and morphology of dendritic growth. We report values of gamma for the hard-sphere and Lennard-Jones systems, as well as recent results on the series of inverse-power potentials. For the inverse sixth-, seventh-, and eighth-power systems, we determine gamma for both fcc and bcc crystal structures. For these systems, the bcc-melt gamma is lower than that for fcc crystals by about 25%, consistent with recent experiments and computer simulations on fcc-forming systems that show preferential formation of bcc nuclei in the initial stages of crystallization. In addition, we show that our results give a molecular interpretation of Turnbull's rule, which is the empirical relationship between gamma and the enthalpy of fusion.     

Revisiting the Frenkel-Ladd method to compute the free energy of solids: the Einstein molecule approach

In this paper a new method to evaluate the free energy of solids is proposed. The method can be regarded as a variant of the method proposed by Frenkel and Ladd [J. Chem. Phys. 81, 3188 (1984)]. The main equations of the method can be derived in a simple way. The method can be easily implemented within a Monte Carlo program. We have applied the method to determine the free energy of hard spheres in the solid phase for several system sizes. The obtained free energies agree within the numerical uncertainty with those obtained by Polson et al. [J. Chem. Phys. 112, 5339 (2000)]. The fluid-solid equilibria has been determined for several system sizes and compared to the values published previously by Wilding and Bruce [Phys. Rev. Lett. 85, 5138 (2000)] using the phase switch methodology. It is shown that both the free energies and the coexistence pressures present a strong size dependence and that the results obtained from free energy calculations agree with those obtained using the phase switch method, which constitutes a cross-check of both methodologies. From the results of this work we estimate the coexistence pressure of the fluid-solid transition of hard spheres in the thermodynamic limit to be p*=11.54(4), which is slightly lower than the classical value of Hoover and Ree (p*=11.70) [J. Chem. Phys. 49, 3609 (1968)]. Taking into account the strong size dependence of the free energy of the solid phase, we propose to introduce finite size corrections, which allow us to estimate approximately the free energy of the solid phase in the thermodynamic limit from the known value of the free energy of the solid phase with N molecules. We have also determined the free energy of a Lennard-Jones solid by using both the methodology of this work and the finite size correction. It is shown how a relatively good estimate of the free energy of the system in the thermodynamic limit is obtained even from the free energy of a relatively small system.     

Molecular simulation of an aerosol OT reverse micelle: 1. The shape and structure of a micelle

The simulation of an Aerosol OT micelle in the apolar environment is performed via the molecular dynamics method in the approximation of a coarse-grain model. The mean size and shape of a micelle, as well as its molecular structure, are determined as functions of the water-surfactant ratio and aggregation number. Geometric parameters of aggregates are estimated through calculation of the inertia tensors of its internal portion under the assumption of an ellipsoidal shape of a micelle. Radial profiles of the partial density and pair correlation functions are obtained, which are used to calculate coordination numbers for water molecules, counterions, and surface-active ions. The most probable arrangement of water molecules and surfactant anions are found on the basis of orientation distribution functions.     

Calculation of absolute ligand binding free energy to a ribosome-targeting protein as a function of solvent model

A comparative analysis is provided of the effect of different solvent models on the calculation of a potential of mean force (PMF) for determining the absolute binding affinity of the small molecule inhibitor pteroic acid bound to ricin toxin A-chain (RTA). Solvent models include the distance-dependent dielectric constant, several different generalized Born (GB) approximations, and a hybrid explicit/GB-based implicit solvent model. We found that the simpler approximation of dielectric screening and a GB model, with Born radii fitted to a switching-window dielectric-boundary surface Poisson solvent model, severely overpredicted the binding affinity as compared to the experimental value, estimated to range from -4.4 to -6.0 kcal/mol. In contrast, GB models that are parametrized to fit the Lee-Richards molecular surface performed much better, predicting binding free energy within 1-3 kcal/mol of experimental estimates. However, the predicted free-energy profiles of these GB models displayed alternative binding modes not observed in the crystal structure. Finally, the most rigorous and computationally costly approach in this work, which used a hybrid explicit/implicit solvent model, correctly determined a binding funnel in the PMF near the crystallographic bound state and predicted an absolute binding affinity that was 2 kcal/mol more favorable than the estimated experimental binding affinity.     

The thickness of a liquid layer on the free surface of ice as obtained from computer simulation

Molecular dynamic simulations were performed for ice I(h) with a free surface by using four water models, SPC/E, TIP4P, TIP4P/Ice, and TIP4P/2005. The behavior of the basal plane, the primary prismatic plane, and of the secondary prismatic plane when exposed to vacuum was analyzed. We observe the formation of a thin liquid layer at the ice surface at temperatures below the melting point for all models and the three planes considered. For a given plane it was found that the thickness of a liquid layer was similar for different water models, when the comparison is made at the same undercooling with respect to the melting point of the model. The liquid layer thickness is found to increase with temperature. For a fixed temperature it was found that the thickness of the liquid layer decreases in the following order: the basal plane, the primary prismatic plane, and the secondary prismatic plane. For the TIP4P/Ice model, a model reproducing the experimental value of the melting temperature of ice, the first clear indication of the formation of a liquid layer, appears at about -100 degrees C for the basal plane, at about -80 degrees C for the primary prismatic plane, and at about -70 degrees C for the secondary prismatic plane.     

A comparison of methods to compute the potential of mean force.

Most processes occurring in a system are determined by the relative free energy between two or more states because the free energy is a measure of the probability of finding the system in a given state. When the two states of interest are connected by a pathway, usually called reaction coordinate, along which the free-energy profile is determined, this profile or potential of mean force (PMF) will also yield the relative free energy of the two states. Twelve different methods to compute a PMF are reviewed and compared, with regard to their precision, for a system consisting of a pair of methane molecules in aqueous solution. We analyze all combinations of the type of sampling (unbiased, umbrella-biased or constraint-biased), how to compute free energies (from density of states or force averaging) and the type of coordinate system (internal or Cartesian) used for the PMF degree of freedom. The method of choice is constraint-bias simulation combined with force averaging for either an internal or a Cartesian PMF degree of freedom.     

Application of ESCA to polymer chemistry. XVI. Electron mean free paths as a function of kinetic energy in polymeric films determined by means of ESCA

Electron mean free paths as a function of kinetic energy have been measured by the substrate overlayer technique for in situ-polymerized films of poly(p-xylylene) and the monochloro- and monobromo-substituted derivatives. The results are compared with previous estimates of mean free paths available in the literature for organic materials. Comparison is also drawn with corresponding experimental data for typical metals and semiconductors, and it is shown that organic polymers fit into a consistent picture which may be rationalized on the basis of existing theory. For electrons of kinetic energy ～969 eV, ～1170 eV, 1202 eV, and 1403 eV, mean free paths of ～14 Å, ～22 Å, ～23 Å, and ～29 Å, respectively, are obtained for the poly(p-xylylene) polymer films studied in this work.     

Application of the multiensemble sampling to the hydration free energy

We demonstrate the feasibility of using multiensemble sampling method (MESM) to determine the free energy difference between two far states for which the configurational distributions do not overlap at all. The MESM is a recently developed non‐Boltzmann sampling technique. The free energy of charging a sodium ion in water is accurately calculated in a single simulation, introducing nine intermediate ionic states. This is due to the ability of the method to explore the relevant parts of configuration space equally for every state, and this ability comes from the universality of weighting function W and the simplicity in adjusting its parameters. Detailed procedures of adjusting the parameters are presented. The comparison with a free energy perturbation method (FEPM) shows that the MESM is more reliable and efficient. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1004–1009, 2001     

Inhomogeneous polymer solutions: Theory of the interfacial free energy and the composition profile near the consolute point

The Flory–Huggins formulation of the combinatorial entropy, supplemented with residual free energy, is applied locally to obtain the interfacial free energy and the concentration profile of polymer in the interface between two demixed polymer solution phases. Two choices were investigated for the residual free energy: a “regular solution” formulation and an empirical formulation of Koningsveld for polystyrene in cyclohexane. Asymptotic, analytical solutions of the equations near the critical solution point and solutions obtained by numerical calculations are given as a function of temperature for several molecular weights. At temperatures farther below the critical temperature the equations have no solutions. The reason for this is not entirely clear. The local formulation of the free energy used here is an improved version of a previous one, which gave wrong results for asymmetric systems (polymer in a low molecular weight solvent). This newer version is consistent with our theory of critical opalescence and gives a relation between the interface “thickness” and the correlation range of the concentration fluctuations. The calculated correlation ranges were in good accord with those found experimentally by Debye, Chu, and Woerman. That the newer version of our equations for an interface gives no acceptable solutions at lower temperatures could be caused by a “collapse” of a diffuse to a sharp interface as suggested by Nose.     

Rapid prediction of solvation free energy. 3. Application to the SAMPL2 challenge

The SAMPL2 hydration free energy blind prediction challenge consisted of a data set of 41 molecules divided into three subsets: explanatory, obscure and investigatory, where experimental hydration free energies were given for the explanatory, withheld for the obscure, and not known for the investigatory molecules. We employed two solvation models for this challenge, a linear interaction energy (LIE) model based on explicit-water molecular dynamics simulations, and the first-shell hydration (FiSH) continuum model previously calibrated to mimic LIE data. On the 23 compounds from the obscure (blind) dataset, the prospectively submitted LIE and FiSH models provided predictions highly correlated with experimental hydration free energy data, with mean-unsigned-errors of 1.69 and 1.71 kcal/mol, respectively. We investigated several parameters that may affect the performance of these models, namely, the solute flexibility for the LIE explicit-solvent model, the solute partial charging method, and the incorporation of the difference in intramolecular energy between gas and solution phases for both models. We extended this analysis to the various chemical classes that can be formed within the SAMPL2 dataset. Our results strengthen previous findings on the excellent accuracy and transferability of the LIE explicit-solvent approach to predict transfer free energies across a wide spectrum of functional classes. Further, the current results on the SAMPL2 test dataset provide additional support for the FiSH continuum model as a fast yet accurate alternative to the LIE explicit-solvent model. Overall, both the LIE explicit-solvent model and the FiSH continuum solvation model show considerable improvement on the SAMPL2 data set over our previous continuum electrostatics-dispersion solvation model used in the SAMPL1 blind challenge.     

Use of computer simulation to determine the triplet distribution function in dense fluids

In this paper a new method is described for determining the triplet distribution function g₃ from computer simulations of dense fluids. The method involves expressing g₃ in terms of hyperspherical coordinates (?,θ,φ) and then expanding g₃ in spherical harmonics of the angular variables θ and φ. We first test the convergence of the spherical harmonic expansion by applying the method to the superposition approximation for g₃. We then show how the expansion coefficients may be obtained from an analysis of computer simulation data. Lastly, we report new results for g₃ obtained by applying the method to data from a molecular dynamics simultion of the Lennard-Jones fluid.     

Application of semiempirical quantum chemical methods as a scoring function in docking

A crucial point in docking simulations is the scoring function used for estimation of the target-ligand interaction energy. The usual practice is to employ fast but simplified empirical scoring functions. Rigorous quantum chemical methods are too slow to screen virtual combinatorial libraries consisting of thousands of molecules, but they can be used in the final step of the simulations for assessing the results obtained. At this stage quantum chemical calculations can be performed only for the 10–100 top binders predicted by simplified scoring functions, and only using linear-scaling semiempirical quantum chemical methods such as MOZYME. The possibilities and potentialities of the quantum chemical methods for estimation of the binding affinities in docking simulations are a largely unexplored area, so the main goal of this study is a detailed evaluation of the potential and limitations of the MOZYME methodology for estimation of the target-ligand binding energies and its comparison with available experimental data.Proceedings of the 11th International Congress of Quantum Chemistry satellite meeting in honor of Jean-Louis Rivail