Relationship between Solvation Thermodynamics from IST and DFT Perspectives.

The potential of mean force between two large parallel hydrophobic oblate ellipsoidal plates in liquid water is determined by molecular dynamics. Each ellipsoid displaces approximately 40 water molecules and has major and minor axes of 3.1 and 9.3 A, respectively, has a surface area of 650 A*, and interacts repulsively with the solvent water molecules. The potential of mean force is calculated from thermodynamic perturbation theory for a series of decreasing plate separations, using constant-pressure molecular dynamics. As the plates are moved together, they are first separated by three water layers and then by two, but for shorter distances, a dewetting transition occurs, and one water layer is never observed despite the fact that one can fit. As the plates are brought together, there is a corresponding weak oscillation in the potential of mean force corresponding to the removal of each water layer until the dewetting transition takes place, and for closer separations, the surrounding water molecules induce a constant average attractive force of 25 (kJ/mol)/A between the plates. This hydrophobic attraction is largely entropic in character, and the potential of mean force is found to be proportional to the area of the water-vacuum surface in this dewetting regime. The constant of proportionality is found to be smaller than the gas-liquid surface tension of the water model used. There is a very strong short-range driving force toward contact pairing. 1.

Water molecules at the binding interface of biomolecular complexes or water molecules displaced from hydrophobic cavities have lately been recognized as important modulators of binding affinity. One approach to computing the contribution of these water molecules to solvation thermodynamics is inhomogeneous fluid solvation theory (IFST). Over the past few years this approach has been applied to interfacial water molecules, both individual and in clusters. Our implementation of IFST resulted in the computational package Solvation Thermodynamics of Ordered Water (STOW). This chapter gives an overview of the theory and its applications and describes how to calculate the thermodynamic contributions of ordered water molecules using STOW.

We propose the thermodynamic integration along a spatial reaction coordinate using the molecular dynamics simulation combined with the three-dimensional reference interaction site model theory. This method provides a free energy calculation in solution along the reaction coordinate defined by the Cartesian coordinates of the solute atoms. The proposed method is based on the blue moon algorithm which can, in principle, handle any reaction coordinate as far as it is defined by the solute atom positions. In this article, we apply the present method to the complex formation process of the crown ether 18-Crown-6 (18C6) with the potassium ion in an aqueous solution. The separation between the geometric centers of these two molecules is taken to be the reaction coordinate for this system. The potential of mean force (PMF) becomes the maximum at the separation between the molecular centers being ∼4 Å, which can be identified as the free energy barrier in the process of the molecular recognition. In a separation further than the free energy barrier, the PMF is slightly reduced to exhibit a plateau. In the region closer than the free energy barrier, approach of the potassium ion to the center of 18C6 also decreases the PMF. When the potassium ion is accommodated at the center of 18C6, the free energy is lower by -5.7 ± 0.7 kcal/mol than that at the above mentioned plateau or converged state. By comparing the results with those from the free energy calculation along the coupling parameters obtained in our previous paper [T. Miyata, Y. Ikuta, and F. Hirata, J. Chem. Phys. 133, 044114 (2010)], it is found that the effective interaction in water between 18C6 and the potassium ion vanishes beyond the molecular-center-separation of 10 Å. Furthermore, the conformation of 18C6 is found to be significantly changed depending upon the 18C6-K(+) distance. A proper conformational sampling and an accurate solvent treatment are crucial for realizing the accurate PMF, and we believe that the proposed method is useful to evaluate the PMF in a solution. A discussion upon the PMF in terms of the three-dimensional distribution function for the solvent is also presented.

The swift progression and expansion of machine learning (ML) have not gone unnoticed within the realm of statistical mechanics. In particular, ML techniques have attracted attention by the classical density-functional theory (DFT) community, as they enable automatic discovery of free-energy functionals to determine the equilibrium-density profile of a many-particle system. Within classical DFT, the external potential accounts for the interaction of the many-particle system with an external field, thus, affecting the density distribution. In this context, we introduce a statistical-learning framework to infer the external potential exerted on a classical many-particle system. We combine a Bayesian inference approach with the classical DFT apparatus to reconstruct the external potential, yielding a probabilistic description of the external-potential functional form with inherent uncertainty quantification. Our framework is exemplified with a grand-canonical one-dimensional classical particle ensemble with excluded volume interactions in a confined geometry. The required training dataset is generated using a Monte Carlo (MC) simulation where the external potential is applied to the grand-canonical ensemble. The resulting particle coordinates from the MC simulation are fed into the learning framework to uncover the external potential. This eventually allows us to characterize the equilibrium density profile of the system by using the tools of DFT. Our approach benchmarks the inferred density against the exact one calculated through the DFT formulation with the true external potential. The proposed Bayesian procedure accurately infers the external potential and the density profile. We also highlight the external-potential uncertainty quantification conditioned on the amount of available simulated data. The seemingly simple case study introduced in this work might serve as a prototype for studying a wide variety of applications, including adsorption, wetting, and capillarity, to name a few.

A new three-parameter (valency, ion size, and charge asymmetry) model, the asymmetric restricted primitive model (ARPM) of ionic liquids, has recently been proposed. Given that ionic liquids generally are composed of monovalent species, the ARPM effectively reduces to a two-parameter model. Monte Carlo (MC) simulations have demonstrated that the ARPM is able to reproduce key properties of room temperature ionic liquids (RTILs) in bulk and at charged surfaces. The relatively modest complexity of the model raises the possibility, which is explored here, that a classical density functional theory (DFT) could resolve its properties. This is relevant because it might generate great improvements in terms of both numerical efficiency and understanding in the continued research of RTILs and their applications. In this report, a DFT for rod-like molecules is proposed as an approximate theoretical tool for an ARPM fluid. Borrowing data on the ion pair fraction from a single bulk simulation, the ARPM is modelled as a mixture of dissociated ions and connected ion pairs. We have specifically studied an ARPM where the hard-sphere diameter is 5 Å, with the charge located 1 Å from the hard-sphere centre. We focus on fluid structure and electrochemical behaviour of this ARPM fluid, into which a model electrode is immersed. The latter is modelled as a perfect conductor, and surface polarization is handled by the method of image charges. Approximate methods, which were developed in an earlier study, to take image interactions into account, are also incorporated in the DFT. We make direct numerical comparisons between DFT predictions and corresponding simulation data. The DFT theory is implemented both in the normal mean field form with respect to the electrostatic interactions and in a correlated form based on hole formation by both steric repulsions and ion-ion Coulomb interactions. The results clearly show that ion-ion correlations play a very important role in the screening of the charged surfaces by our ARPM ionic liquid. We have studied electrostatic potentials and ion density profiles as well the differential capacitance. The mean-field DFT fails to reproduce these properties, but the inclusion of ion-ion correlation by a simple approximate treatment yields quite reasonable agreement with the corresponding simulation results. An interesting finding is that there appears to be a surface phase transition at relatively low surface charge which is readily explored by DFT, but seen also in the MC simulations at somewhat higher asymmetry.

The ongoing scientific interest in the properties and structure of electric double layers (EDLs) stems from their pivotal role in (super)capacitive energy storage, energy harvesting, and water treatment technologies. Classical density functional theory (DFT) is a promising framework for the study of the in- and out-of-plane structural properties of double layers. Supported by molecular dynamics simulations, we demonstrate the adequate performance of DFT for analyzing charge layering in the EDL perpendicular to the electrodes. We discuss charge storage and capacitance of the EDL and the impact of screening due to dielectric solvents. We further calculate, for the first time, the in-plane structure of the EDL within the framework of DFT. While our out-of-plane results already hint at structural in-plane transitions inside the EDL, which have been observed recently in simulations and experiments, our DFT approach performs poorly in predicting in-plane structure in comparison to simulations. However, our findings isolate fundamental issues in the theoretical description of the EDL within the primitive model and point towards limitations in the performance of DFT in describing the out-of-plane structure of the EDL at high concentrations and potentials.

This article reviews the recent progress on predicting the adsorption properties of metal-organic framework by using classical density functional theory and focused on the application of the classical density functional theory to the high-throughput screening, which is accelerated by fast Fourier Transform. Comparing to the conventional molecular simulations, the advantage of the accelerated classical density functional theory is the calculation speed, especially for simple small molecule systems, which makes the high-throughput screening on MOF materials feasible. However, it appears that there is a lack of efficient method to deal with the complicated molecules. How to construct a reasonable free energy functional of complicated fluid is the main challenge to state of art classical density functional theory. In a word, the improvement of CDFT theory and the combination of CDFT and molecular simulation are the two main ways for CDFT to predict gas adsorption in MOF.

Canfosfamide (TLK286, TELCYTA) is a prodrug that upon activation by glutathione transferase P1-1 (GST P1-1) yields an anticancer alkylating agent and a glutathione derivative. The rationale underlying the use of TLK286 in chemotherapy is that tumor cells overexpressing GST P1-1 will be locally exposed to the released alkylating agent with limited collateral toxicity to the surrounding normal tissues. TLK286 has demonstrated clinical effects in phase II and III clinical trials for the treatment of malignancies, such as ovarian cancer, nonsmall cell lung cancer, and breast cancer, as a single agent and in combination with other chemotherapeutic agents. In spite of these promising results, the detailed mechanism of GST P1-1 activation of the prodrug has not been elucidated. Here, we propose a mechanism for the TLK286 activation by GST P1-1 on the basis of density functional theory (DFT) and on potential of mean force (PMF) calculations. A catalytic water molecule is instrumental to the activation by forming a network of intermolecular interactions between the active-site Tyr7 hydroxyl and the sulfone and COO(-) groups of TLK286. The results obtained are consistent with the available experimental kinetic data and provide an atomistic understanding of the TLK286 activation mechanism.

The dissociative electron transfer reaction CH3Cl+e−→CH3•+Cl− in aqueous solution is studied by using a QM/MM method. In this work the quantum subsystem (a methylchloride molecule plus an electron) is described using density functional theory while the solvent (300 water molecules) is described using the TIP3P classical potential. By means of molecular dynamics simulations and the thermodynamic integration technique we obtained the potential of mean force (PMF) for the carbon–chlorine bond dissociation of the neutral and radical anion species. Combining these two free energy curves we found a quadratic dependence of the activation free energy on the reaction free energy in agreement with Marcus’ relationship, originally developed for electron transfer processes not involving bond breaking. We also investigated dynamical aspects by means of 60 dissociative trajectories started with the addition of an extra electron to different configurations of a methylchloride molecule in solution. The PMF shows the existence of a very flat region, in which the system is trapped during some finite time if the quantum subsystem quickly losses its excess kinetic energy transferring it to the solvent molecules. One of the most important factors determining the effectiveness of this energy transfer seems to be the existence of close contacts (hydrogen bonds) between the solute and the solvent.

We study by molecular dynamics simulations the driving force for the hydrophobic interaction between graphene sheets of different sizes down to the atomic scale. Similar to the prediction by Lum, Chandler, and Weeks for hard-sphere solvation [J. Phys. Chem. B 1999, 103, 4570-4577], we find the driving force to be length-scale dependent, despite the fact that our model systems do not exhibit dewetting. For small hydrophobic solutes, the association is purely entropic, while enthalpy favors dissociation. The latter is demonstrated to arise from the enhancement of hydrogen bonding between the water molecules around small hydrophobes. On the other hand, the attraction between large graphene sheets is dominated by enthalpy which mainly originates from direct solute-solute interactions. The crossover length is found to be inside the range of 0.3-1.5 nm(2) of the surface area of the hydrophobe that is eliminated in the association process. In the large-scale regime, different thermodynamic properties are scalable with this change of surface area. In particular, upon dimerization, a total and a water-induced stabilization of approximately 65 and 12 kJ/mol/nm(2) are obtained, respectively, and on average around one hydrogen bond is gained per 1 nm(2) of graphene sheet association. Furthermore, the potential of mean force between the sheets is also scalable except for interplate distances smaller than 0.64 nm which corresponds to the region around the barrier for removing the last layer of water. It turns out that, as the surface area increases, the relative height of the barrier for association decreases and the range of attraction increases. It is also shown that, around small hydrophobic solutes, the lifetime of the hydrogen bonds is longer than in the bulk, while around large hydrophobes it is the same. Nevertheless, the rearrangement of the hydrogen-bond network for both length-scale regimes is slower than in bulk water.

Changes in free energy provide valuable information for molecular recognition, including both ligand-receptor binding thermodynamics and kinetics. Umbrella sampling (US), a widely used free energy calculation method, has long been used to explore the dissociation process of ligand-receptor systems and compute binding free energy. In existing publications, the binding free energy computed from the potential of mean force (PMF) with US simulation mostly yielded "ball park" values with experimental data. However, the computed PMF values are highly influenced by factors such as initial conformations and/or trajectories provided, the reaction coordinate, and sampling of conformational space in each US window. These critical factors have rarely been carefully studied. Here we used US to study the guest aspirin and 1-butanol dissociation processes of β-cyclodextrin (β-CD) and an inhibitor SB2 dissociation from a p38α mitogen-activated protein kinase (MAPK) complex. For β-CD, we used three different β-CD conformations to generate the dissociation path with US windows. For p38α, we generated the dissociation pathway by using accelerated molecular dynamics followed by conformational relaxing with short conventional MD, steered MD, and manual pulling. We found that, even for small β-CD complexes, different β-CD conformations altered the height of the PMF, but the pattern of PMF was not affected if the MD sampling in each US window was well-converged. Because changing the macrocyclic ring conformation needs to rotate dihedral angles in the ring, a bound ligand largely restrains the motion of cyclodextrin. Therefore, once a guest is in the binding site, cyclodextrin cannot freely change its initial conformation, resulting in different absolute heights of the PMF, which cannot be overcome by running excessively long MD simulations for each US window. Moreover, if the US simulations were not converged, the important barrier and minimum were missed. For ligand-protein systems, our studies also suggest that the dissociation trajectories modeled by an enhanced sampling method must maintain a natural molecular movement to avoid biased PMF plots when using US simulations.

Enhanced removal efficiency of heavy metal ions from wastewater through functionalized silicon carbide membrane: A theoretical study

The interactions of amino acid side-chains confer diverse energetic contributions and physical properties to a protein's stability and function. Various computational tools estimate the effect of changing a given amino acid on the protein's stability based on parametrized (free) energy functions. When parametrized for the prediction of protein stability in water, such energy functions can lead to suboptimal results for other solvents, such as ionic liquids (IL), aqueous ionic liquids (aIL), or salt solutions. However, to our knowledge, no comprehensive data are available describing the energetic effects of aIL on intramolecular protein interactions. Here, we present the most comprehensive set of potential of mean force (PMF) profiles of pairwise protein-residue interactions to date, covering 50 relevant interactions in water, the two biotechnologically relevant aIL [BMIM/Cl] and [BMIM/TfO], and [Na/Cl]. These results are based on a cumulated simulation time of >1 ms. aIL and salt ions can weaken, but also strengthen, specific residue interactions by more than 3 kcal mol-1, depending on the residue pair, residue-residue configuration, participating ions, and concentration, necessitating considering such interactions specifically. These changes originate from a complex interplay of competitive or cooperative noncovalent ion-residue interactions, changes in solvent structural dynamics, or unspecific charge screening effects and occur at the contact distance but also at larger, solvent-separated distances. This data provide explanations at the atomistic and energetic levels for complex IL effects on protein stability and should help improve the prediction accuracies of computational tools that estimate protein stability based on (free) energy functions.

In classical density functional theory (DFT), the part of the Helmholtz free energy functional arising from attractive inter-particle interactions is often treated in a mean-field or van der Waals approximation. On the face of it, this is a somewhat crude treatment as the resulting functional generates the simple random phase approximation (RPA) for the bulk fluid pair direct correlation function. We explain why using standard mean-field DFT to describe inhomogeneous fluid structure and thermodynamics is more accurate than one might expect based on this observation. By considering the pair correlation function g(x) and structure factor S(k) of a one-dimensional model fluid, for which exact results are available, we show that the mean-field DFT, employed within the test-particle procedure, yields results much superior to those from the RPA closure of the bulk Ornstein-Zernike equation. We argue that one should not judge the quality of a DFT based solely on the approximation it generates for the bulk pair direct correlation function.

Classical density functional theory (DFT) is the primary method for investigations of inhomogeneous fluids in external fields. It requires the excess Helmholtz free energy functional as input to an Euler-Lagrange equation for the one-body density. A variant of this methodology, the force-DFT, uses instead the Yvon-Born-Green equation to generate density profiles. It is known that the latter are consistent with the virial route to the thermodynamics, while DFT is consistent with the compressibility route. In this work we will show an alternative DFT scheme using the Lovett-Mou-Buff-Wertheim (LMBW) equation to obtain density profiles, that are shown to be also consistent with the compressibility route. However, force-DFT and LMBW DFT can both be implemented using a closure relation on the level of the two-body correlation functions. This is proven to be an advantageous feature, opening the possibility of an optimization scheme in which the structural inconsistency between different routes to the density profile is minimized. ("Structural inconsistency" is a generalization of the notion of thermodynamic inconsistency, familiar from bulk integral equation studies). Numerical results are given for the density profiles of two-dimensional systems of hard-core Yukawa particles with a repulsive or an attractive tail, in planar geometry.