Application of the constrained fluid lambda-integration path to the calculation of high temperature Au(110) surface free energies

Recently a method termed constrained fluid lambda-integration was proposed for calculating the free energy difference between bulk solid and liquid reference states via the construction of a reversible thermodynamic integration path; coupling the two states in question. The present work shows how the application of the constrained fluid lambda-integration concept to solid/liquid slab simulation cells makes possible a generally applicable computer simulation methodology for calculating the free energy of any surface and/or surface defect structure, including surfaces requiring variations in surface atom or density number, such as the (1 x 5) Au(100) or (1 x 2) missing row Au(110) reconstructed surfaces or excess adatom/vacancy/step populated surfaces. We evaluate the methodology by calculating the free energy of various disordered high temperature Au(110) embedded atom method surfaces constrained to differing excess surface atom numbers [including those corresponding to the (1 x 2) missing row reconstructed surface] and obtained the interesting result that at 1000 K (as distinct from lower temperatures) the free energy difference between these surfaces is reduced to zero; a result which is consistent with an expected order-disorder phase transition for the Au(110) surface at such high temperatures.     

On the calculation of surface energies

A simple formalism based on response theory is used to evaluate expressions for the surface energy of a uniform dielectric or metallic medium and the change in surface energy of a dielectric solvent when a dilute electrolyte is dissolved in it. The latter result coincides with the result of Onsager and Samaras.     

On the calculation of absolute free energies from molecular-dynamics or Monte Carlo data

We propose a method for calculating absolute free energies from Monte Carlo or molecular-dynamics data. The method is based on the identity that expresses the partition function Q as a Boltzmann average: 1Q=w(p,x)exp[betaH(p,x)], where w(p,x) is an arbitrary weight function such that its integral over the phase space is equal to 1. In practice, to minimize statistical errors the weight function is chosen such that the regions of the phase space where sampling statistics are poor are excluded from the average. The "ideal" weight function would be the equilibrium phase-space density exp[-betaH(p,x)]Q itself. We consider two methods for constructing the weight function based on different estimates of the equilibrium phase-space density from simulation data. In the first method, it is chosen to be a Gaussian function, whose parameters are obtained from the covariance matrix of the atomic coordinates. In the second, a clustering algorithm is used to attempt partitioning the data into clusters corresponding to different basins of attraction visited by the system. The weight function is then constructed as a superposition of Gaussians calculated for each cluster separately. We show that these strategies can be used to improve upon previous methods of estimating absolute entropies from covariance matrices.     

On the calculation of Auger energies

“Frozen-orbital” expressions within an average-energy formalism are presented for the calculation of Auger transition energies based on the orbitals determined for the initial state. Trial calculations at the minimal-basis level are reported for the K-LL transition in Ne, Na and Mg. We also discuss an ambiguity in the use of certain models for the calculation of both semi-empirical and ab initio Auguer energies. For illustration we give some numerical data on the atoms Ne, Na, Mg and S.     

Potential dependence of step and kink energies on Au(100) electrodes in sulfuric acid

Using electrochemical STM we studied monolayer high Au islands on Au(100) electrodes in sulfuric acid as a function of the electrode potential. We made use of theoretical and experimental methods recently developed for UHV experiments on metal islands. It is demonstrated that these models are likewise applicable to islands on metal electrodes in a liquid environment. From a quantitative analysis of the equilibrium island shape and of the island shape fluctuations we determined the step free energy (line tension) as a function of orientation and the kink energy, and the dependence of these quantities on the electrode potential. In a first approach to a theoretical understanding the electrostatic contributions to the line tension are considered. It is concluded that these contributions should add significantly to the observed variation with the potential. This fails however to provide essential features of the experimental result.     

Direct calculation of the binding free energies of FKBP ligands

Direct calculations of the absolute free energies of binding for eight ligands to FKBP protein were performed using the Fujitsu BioServer massively parallel computer. Using the latest version of the general assisted model building with energy refinement (AMBER) force field for ligand model parameters and the Bennett acceptance ratio for computing free-energy differences, we obtained an excellent linear fit between the calculated and experimental binding free energies. The rms error from a linear fit is 0.4 kcal/mol for eight ligand complexes. In comparison with a previous study of the binding energies of these same eight ligand complexes, these results suggest that the use of improved model parameters can lead to more predictive binding estimates, and that these estimates can be obtained with significantly less computer time than previously thought. These findings make such direct methods more attractive for use in rational drug design.     

A water-swap reaction coordinate for the calculation of absolute protein-ligand binding free energies

The accurate prediction of absolute protein-ligand binding free energies is one of the grand challenge problems of computational science. Binding free energy measures the strength of binding between a ligand and a protein, and an algorithm that would allow its accurate prediction would be a powerful tool for rational drug design. Here we present the development of a new method that allows for the absolute binding free energy of a protein-ligand complex to be calculated from first principles, using a single simulation. Our method involves the use of a novel reaction coordinate that swaps a ligand bound to a protein with an equivalent volume of bulk water. This water-swap reaction coordinate is built using an identity constraint, which identifies a cluster of water molecules from bulk water that occupies the same volume as the ligand in the protein active site. A dual topology algorithm is then used to swap the ligand from the active site with the identified water cluster from bulk water. The free energy is then calculated using replica exchange thermodynamic integration. This returns the free energy change of simultaneously transferring the ligand to bulk water, as an equivalent volume of bulk water is transferred back to the protein active site. This, directly, is the absolute binding free energy. It should be noted that while this reaction coordinate models the binding process directly, an accurate force field and sufficient sampling are still required to allow for the binding free energy to be predicted correctly. In this paper we present the details and development of this method, and demonstrate how the potential of mean force along the water-swap coordinate can be improved by calibrating the soft-core Coulomb and Lennard-Jones parameters used for the dual topology calculation. The optimal parameters were applied to calculations of protein-ligand binding free energies of a neuraminidase inhibitor (oseltamivir), with these results compared to experiment. These results demonstrate that the water-swap coordinate provides a viable and potentially powerful new route for the prediction of protein-ligand binding free energies.     

Density functional calculations of surface free energies

We propose a general method of thermodynamic integration to find the free energy of a surface, where our integration parameter is taken to be the strain on the unit cell of the system (which in the example presented in this paper is simply the extension of the unit cell along the normal to the surface), and the integration is performed over the thermal average stress from a molecular dynamics run. In order to open up a vacuum gap in a continuous and reversible manner, an additional control interaction has been introduced. We also use temperature integration to find a linear relation for the temperature dependence of the free surface energy. These methods have been applied to the titanium dioxide (110) surface, using first principles density functional theory. A proof of principle calculation for zero temperature shows excellent agreement between the integral calculation and the difference in energy calculated by the DFT program. Calculations that have been performed at 295 and 1000 K give excellent agreement between the two integration methods.     

NO adsorption on reconstructed and unreconstructed Pt(100) surface at 300 K TDS studies

The consequence of filling NO and N₂ thermodesorption states and the adsorption rate are found to depend on the initial surface structure. The initial sticking coefficients for (1×1) and (hex) structures are 1 and 0.35, respectively. If _NO < 0.3 ML, the dissociation probability of NO is shown to be higher when adsorption occurs on the unreconstructed surface.     

Rosenbluth-sampled nonequilibrium work method for calculation of free energies in molecular simulation

We present methods that introduce concepts from Rosenbluth sampling [M. N. Rosenbluth and A. W. Rosenbluth, J. Chem. Phys. 23, 356 (1955)] into the Jarzynski nonequilibrium work (NEW) free-energy calculation technique [C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)]. The proposed hybrid modifies the way steps are taken in the NEW process. With it, each step is selected from a range of alternatives, with bias given to steps that contribute the least work. The definition of the work average is modified to account for the bias. We introduce two variants of this method, lambda-bias sampling and configuration-bias sampling, respectively; a combined lambda- and configuration-bias method is also considered. By reducing the likelihood that large nonequilibrated work values enter the ensemble average, the Rosenbluth sampling aids in remedying problems of inaccuracy of the calculation. We demonstrate the performance of the proposed methods through a model system of N independent harmonic oscillators. This model captures the difficulties involved in calculating free energies in real systems while retaining many tractable features that are helpful to the study. We examine four variants of this model that differ qualitatively in the nature of their phase-space overlap. Results indicate that the lambda-bias sampling method is most useful for systems with entropic sampling barriers, while the configuration-bias methods are best for systems with energetic sampling barriers. The Rosenbluth-sampling schemes yield much more accurate results than the unbiased nonequilibrium work method. Typically the accuracy can be improved by about an order of magnitude for a given amount of sampling; this improvement translates into two or more orders of magnitude less sampling required to obtain a given level of accuracy, owing to the generally slow convergence of the NEW calculation when the inaccuracy is large.     

Rapid calculation of partition functions and free energies of fluids

The partition function (Q) is a central quantity in statistical mechanics. All the thermodynamic properties can be derived from it. Here we show how the partition function of fluids can be calculated directly from simulations; this allows us to obtain the Helmholtz free energy (F) via F = -k(B)T ln Q. In our approach, we divide the density of states, assigning half of the configurations found in a simulation to a high-energy partition and half to a low-energy partition. By recursively dividing the low-energy partition into halves, we map out the complete density of states for a continuous system. The result allows free energy to be calculated directly as a function of temperature. We illustrate our method in the context of the free energy of water.     

An efficient method for the calculation of quantum mechanics/molecular mechanics free energies

The combination of quantum mechanics (QM) with molecular mechanics (MM) offers a route to improved accuracy in the study of biological systems, and there is now significant research effort being spent to develop QM/MM methods that can be applied to the calculation of relative free energies. Currently, the computational expense of the QM part of the calculation means that there is no single method that achieves both efficiency and rigor; either the QM/MM free energy method is rigorous and computationally expensive, or the method introduces efficiency-led assumptions that can lead to errors in the result, or a lack of generality of application. In this paper we demonstrate a combined approach to form a single, efficient, and, in principle, exact QM/MM free energy method. We demonstrate the application of this method by using it to explore the difference in hydration of water and methane. We demonstrate that it is possible to calculate highly converged QM/MM relative free energies at the MP2/aug-cc-pVDZ/OPLS level within just two days of computation, using commodity processors, and show how the method allows consistent, high-quality sampling of complex solvent configurational change, both when perturbing hydrophilic water into hydrophobic methane, and also when moving from a MM Hamiltonian to a QM/MM Hamiltonian. The results demonstrate the validity and power of this methodology, and raise important questions regarding the compatibility of MM and QM/MM forcefields, and offer a potential route to improved compatibility.     

SAMPL4, a blind challenge for computational solvation free energies: the compounds considered

For the fifth time I have provided a set of solvation energies (1 M gas to 1 M aqueous) for a SAMPL challenge. In this set there are 23 blind compounds and 30 supplementary compounds of related structure to one of the blind sets, but for which the solvation energy is readily available. The best current values of each compound are presented along with complete documentation of the experimental origins of the solvation energies. The calculations needed to go from reported data to solvation energies are presented, with particular attention to aspects which are new to this set. For some compounds the vapor pressures (VP) were reported for the liquid compound, which is solid at room temperature. To correct from VP_{subcooled liquid} to VP_sublimation requires ΔS_fusion, which is only known for mannitol. Estimated values were used for the others, all but one of which were benzene derivatives and expected to have very similar values. The final compound for which ΔS_fusion was estimated was menthol, which melts at 42 °C so that modest errors in ΔS_fusion will have little effect. It was also necessary to look into the effects of including estimated values of ΔCp on this correction. The approximate sizes of the effects of inclusion of ΔCp in the correction from VP_{subcooled liquid} to VP_sublimation were estimated and it was noted that inclusion of ΔCp invariably makes ΔG_S more positive. To extend the set of compounds for which the solvation energy could be calculated we explored the use of boiling point (b.p.) data from Reaxys/Beilstein as a substitute for studies of the VP as a function of temperature. B.p. data are not always reliable so it was necessary to develop a criterion for rejecting outliers. For two compounds (chlorinated guaiacols) it became clear that inclusion represented overreach; for each there were only two independent pressure, temperature points, which is too little for a trustworthy extrapolation. For a number of compounds the extrapolation from lowest temperature at which the VP was reported to 25 °C was long (sometimes over 100°) so that it was necessary to consider whether ΔCp might have significant effects. The problem is that there are no experimental values and possible intramolecular hydrogen bonds make estimation uncertain in some cases. The approximate sizes of the effects of ΔCp were estimated, and it was noted that inclusion of ΔCp in the extrapolation of VP down to room temperature invariably makes ΔGs more negative.     

NO adsorption on unreconstructed and reconstructed Pt(100) surface at 300 K: HREELS studies

Dissociative adsorption at 300 K is shown to proceed only on unreconstructed Pt(100)–(1×1) surface. When heated, NO (ad) layers behave in a similar way on both surfaces. Dissociation products are N₂ and O(ad). Adsorption on the (hex) surface is supported to occur in islands. The island size in the saturation layer appears to be 40 Å.     

Structural analysis of the reconstructed Si(001)-C surface

The atomic structure of reconstructed Si(001)c(4 x 4)-C surface has been studied by coaxial impact collision ion scattering spectroscopy. When the 100L of ethylene (C(2)H(4)) molecules have been exposed on Si(001)-(2 x 1) surface at 700 degrees C, it is found that C atoms cause the ordering of missing Si dimer defects and occupy the fourth layer of Si(001) directly below the bridge site. Our results provide the support for the previous model in which a missing dimer structure is accompanied by C incorporation into the subsurface.     

The surface structure of flame-annealed Au(100) in aqueous solution: An STM study

The initial structure of flame-annealed Au(100) surfaces has been studied in air and in 0.1 M H₂SO₄ by scanning tunnelling microscopy (STM). It is shown that before, during and after contact with the electrolyte, at potentials sufficiently negative to prevent specific adsorption of anions, the flame-annealed Au(100) surface is reconstructed into exactly the same “hex” form as a surface which has been prepared by annealing in ultrahigh vacuum (UHV). However, the quality of the reconstructed surface depends sensitively on the sample preparation and on the experimental conditions of the flame-annealing procedure. The influence of the cooling procedure after flame annealing on the initial surface structure of the Au(100) electrode is demonstrated and briefly discussed in the light of results published previously.     

Au adatoms in self-assembly of benzenethiol on the Au(111) surface

Self-assembly of benzenethiol at low coverage on Au(111) was studied using low-temperature scanning tunneling microscopy. Phenylthiolate species (PhS), formed by thermal dehydrogenation of the parent PhSH molecule, was found to self-assemble into surface-bonded complexes with gold adatoms. Each complex involves two PhS species and one gold adatom. The PhS species form either cis- or trans-geometry relative to each other. At a higher coverage, the complexes coalesce, most likely due to the formation of weak C-H...S hydrogen bonds facilitated by the spatial arrangement of the PhS groups. Our findings thus establish that the self-assembly of arenethiols on the Au(111) surface is driven by gold adatom chemistry, which has recently been found to be the key ingredient in the self-assembly of alkanethiols on gold.     

On the calculation of multiplet energies by the hartree-fock-slater method

It is shown that a consistent application of the p ^1/3 approximation of the Hartree-Fock-Slater method requires the use of one specific procedure, the sum method, for the calculation of the energy E _s ¹ of singlet excited states of closed shell molecules. Further, E _s ¹ is found to be in reasonable agreement with experiment for a number of molecules, contrary to the energy E _s ² obtained according to another method discussed in the literature. The calculation of other multiplet splittings than singlet-triplet in the Hartree-Fock-Slater method is also considered.     

On the importance of surface structure in electrochemistry: the well-defined Au (111) single crystal electrode surface

The electrochemical nucleation of cadmium onto vitreous carbon and onto tin oxide electrodes in studied using potentiostatic methods. The influence of supporting electrolyte is investigated and in every case new and interesting features arise. The deposition of cadmium on tin oxide is of particular interest and evidently occurs via an initial process of 2D nucleation followed, at higher overpotentials, by 3D nucleation.