Regularized Bennett and Zwanzig free energy estimators.

We consider the problem of free energy estimation from the general perspective of regularization and Bayes estimation theory. We try to take advantage of an assumed a priori knowledge of the free energy. We reformulate the original Bennett acceptance ratio method, in this perspective, devise a numerical algorithm to solve it, and give a closed formula to estimate the confidence in the prior. Finally, we test the derived estimators by applying them to a toy model.

Mass-zero constrained molecular dynamics for electrode charges in simulations of electrochemical systems.

Classical molecular dynamics simulations have recently become a standard tool for the study of electrochemical systems. State-of-the-art approaches represent the electrodes as perfect conductors, modeling their responses to the charge distribution of electrolytes via the so-called fluctuating charge model. These fluctuating charges are additional degrees of freedom that, in a Born-Oppenheimer spirit, adapt instantaneously to changes in the environment to keep each electrode at a constant potential. Here, we show that this model can be treated in the framework of constrained molecular dynamics, leading to a symplectic and time-reversible algorithm for the evolution of all the degrees of freedom of the system. The computational cost and the accuracy of the new method are similar to current alternative implementations of the model. The advantage lies in the accuracy and long term stability guaranteed by the formal properties of the algorithm and in the possibility to systematically introduce additional kinematic conditions of arbitrary number and form. We illustrate the performance of the constrained dynamics approach by enforcing the electroneutrality of the electrodes in a simple capacitor consisting of two graphite electrodes separated by a slab of liquid water.

Adaptive resolution molecular dynamics simulation through coupling to an internal particle reservoir.

For simulation studies of (macro) molecular liquids it would be of significant interest to be able to adjust or increase the level of resolution within one region of space, while allowing for the free exchange of molecules between open regions of different resolution or representation. We generalize the adaptive resolution idea and suggest an interpretation in terms of an effective generalized grand canonical approach. The method is applied to liquid water at ambient conditions.

Path integral based calculations of symmetrized time correlation functions. I.

In this paper, we examine how and when quantum evolution can be approximated in terms of (generalized) classical dynamics in calculations of correlation functions, with a focus on the symmetrized time correlation function introduced by Schofield. To that end, this function is expressed as a path integral in complex time and written in terms of sum and difference path variables. Taylor series expansion of the path integral's exponent to first and second order in the difference variables leads to two original developments. The first order expansion is used to obtain a simple, path integral based, derivation of the so-called Schofield's quantum correction factor. The second order result is employed to show how quantum mechanical delocalization manifests itself in the approximation of the correlation function and hinders, even in the semiclassical limit, the interpretation of the propagators in terms of sets of guiding classical trajectories dressed with appropriate weights.

Path integral based calculations of symmetrized time correlation functions. II.

Schofield's form of quantum time correlation functions is used as the starting point to derive a computable expression for these quantities. The time composition property of the propagators in complex time is exploited to approximate Schofield's function in terms of a sequence of short time classical propagations interspersed with path integrals that, combined, represent the thermal density of the system. The approximation amounts to linearization of the real time propagators and it becomes exact with increasing number of propagation legs. Within this scheme, the correlation function is interpreted as an expectation value over a probability density defined on the thermal and real path space and calculated by a Monte Carlo algorithm. The performance of the algorithm is tested on a set of benchmark problems. Although the numerical effort required is considerable, we show that the algorithm converges systematically to the exact answer with increasing number of iterations and that it is stable for times longer than those accessible via a brute force, path integral based, calculation of the correlation function. Scaling of the algorithm with dimensionality is also examined and, when the method is combined with commonly used filtering schemes, found to be comparable to that of alternative semiclassical methods.

Transient hydrodynamical behavior by dynamical nonequilibrium molecular dynamics: the formation of convective cells.

We present a method based on dynamical nonequilibrium molecular dynamics (D-NEMD) that allows one to produce rigorous ensemble averages for the transient regimes. We illustrate the method by describing the formation of convective cells within a two-dimensional fluid system of soft disks in which a gravity field and a thermal gradient are present. We analyze two different physical settings, with the thermal gradient orthogonal or parallel to the gravity field. In both settings, we follow the formation of the convective flows from the initial time, when the perturbation is turned on, to the steady state. In the first setting (orthogonal fields) we investigate several different cases, varying the initial stationary ensemble and the perturbing field. We find that the final steady-state convective cell is independent of the specific sequence of perturbation fields, which only affects the transient behavior. In all cases, we find that the convective roll is formed through a sequence of damped oscillations of the local fields (density, temperature, and velocity), superimposed to an overall relaxation toward the local steady-state values. Then, we show how D-NEMD can be applied to the Rayleigh-Bénard (RB) setting (parallel fields). In these conditions, the convective flow only establishes above a threshold, without a preferred verse of rotation. We analyze only the response to the ignition of the gravity field in a stationary system under the action of a vertical thermal gradient. Also in this case we characterize the transient response by following the evolution of the density, temperature, and velocity fields until the steady-state RB convective cell is formed. The observed transients are similar to those observed in the case of orthogonal fields. However, the final steady states are quite different. Finally, we briefly discuss the conditions for the general applicability of the D-NEMD method.

Ab initio simulation of carbon clustering on an Ni(111) surface: a model of the poisoning of nickel-based catalysts.

In this work, we studied the poisoning of a nickel surface due to carbon. Performing ab initio simulations, within the framework of density functional theory, we computed the surface energy of the nickel (111) surface as a function of carbon coverage. On the basis of these results, we can assert that the stable state of the nickel/carbon surface is either a clean nickel surface or a fully carbon-covered nickel surface, which has a graphitic configuration. The relative stability of the two states depends on the temperature and partial pressure of the carbon gas. At fixed nominal carbon coverage, the most stable configurations are those forming carbon clusters. However, the nickel sites hosting these clusters change from hexagonal close packed/face centered cubic (hcp/fcc) sites to on-top sites when going toward larger clusters. This indicates that poisoning due to graphitic patches occurs on on-top sites.

Molecular dynamics simulation of sucrose- and trehalose-coated carboxy-myoglobin.

We performed a room temperature molecular dynamics (MD) simulation on a system containing 1 carboxy-myoglobin (MbCO) molecule in a sucrose-water matrix of identical composition (89% [sucrose/(sucrose + water)] w/w) as for a previous trehalose-water-MbCO simulation (Cottone et al., Biophys J 2001;80:931-938). Results show that, as for trehalose, the amplitude of protein atomic mean-square fluctuations, on the nanosecond timescale, is reduced with respect to aqueous solutions also in sucrose. A detailed comparison as a function of residue number evidences mobility differences along the protein backbone, which can be related to a different efficacy in bioprotection. Different heme pocket structures are observed in the 2 systems. The joint distribution of the magnitude of the electric field at the CO oxygen atom and of the angle between the field and the CO unit vector shows a secondary maximum in sucrose, absent in trehalose. This can explain the CO stretching band profile (A substates distribution) differences evidenced by infrared spectroscopy in sucrose- and trehalose-coated MbCO (Giuffrida et al., J Phys Chem B 2004;108:15415-15421), and in particular the appearance of a further substate in sucrose. Analysis of hydrogen bonds at the protein-solvent interface shows that the fraction of water molecules shared between the protein and the sugar is lower in sucrose than in trehalose, in spite of a larger number of water molecules bound to the protein in the former system, thus indicating a lower protein-matrix coupling, as recently observed by Fourier transform infrared (FTIR) experiments (Giuffrida et al., J Phys Chem B 2004;108:15415-15421).

Soret coefficient for liquid argon-krypton mixtures via equilibrium and nonequilibrium molecular dynamics: a comparison with experiments.

The Soret and the other transport coefficients, characterizing the heat and mass transport in binary mixtures, have been obtained by equilibrium and nonequilibrium molecular dynamics (EMD and NEMD, respectively). Two state points of the argon-krypton mixture are considered, for which experimental values of the Soret coefficient are available. To attempt a comparison between simulations and experiments the common enthalpy-diffusion-free expression for the heat flux has been chosen. The comparison of the simulations with the experiments shows a remarkable agreement, for all the several utilized EMD and NEMD techniques (dynamical and stationary). The techniques, used over 0.3 micros of total simulation time span, are slow convergent but have comparable performances.

Mapping the hydropathy of amino acids based on their local solvation structure.

In spite of its relevant biological role, no general consensus exists on the quantitative characterization of amino acid's hydropathy. In particular, many hydrophobicity scales exist, often producing quite different rankings for the amino acids. To make progress toward a systematic classification, we analyze amino acids' hydropathy based on the orientation of water molecules at a given distance from them as computed from molecular dynamics simulations. In contrast with what is usually done, we argue that assigning a single number is not enough to characterize the properties of an amino acid, in particular when both hydrophobic and hydrophilic regions are present in a residue. Instead we show that appropriately defined conditional probability densities can be used to map the hydrophilic and hydrophobic groups on the amino acids with greater detail than possible with other available methods. Three indicators are then defined based on the features of these probabilities to quantify the specific hydrophobicity and hydrophilicity of each amino acid. The characterization that we propose can be used to understand some of the ambiguities in the ranking of amino acids in the current scales. The quantitative indicators can also be used in combination with standard bioinformatics tools to predict the location of transmembrane regions of proteins. The method is sensitive to the specific environment of the amino acids and can be applied to unnatural and modified amino acids, as well as to other small organic molecules.

Molecular dynamics simulation of carboxy-myoglobin embedded in a trehalose-water matrix.

We report on a molecular dynamics (MD) simulation of carboxy-myoglobin (MbCO) embedded in a water-trehalose system. The mean square fluctuations of protein atoms, calculated at different temperatures in the 100-300 K range, are compared with those from a previous MD simulation on an H2O-solvated MbCO and with experimental data from Mössbauer spectroscopy and incoherent elastic neutron scattering on trehalose-coated MbCO. The results show that, for almost all the atomic classes, the amplitude of the nonharmonic motions stemming from the interconversion among the protein's conformational substates is reduced with respect to the H2O-solvated system, and their onset is shifted toward higher temperature. Moreover, our simulation shows that, at 300 K, the heme performs confined diffusive motions as a whole, leaving the underlying harmonic vibrations unaltered.

Hydration of beta-cyclodextrin: a molecular dynamics simulation study.

We study by molecular dynamics simulations the hydration of beta-cyclodextrin. Our simulations show that within these barrel-shaped molecules hydrophobicity dominates, while at the top and bottom sides of the barrel interactions with water are mostly hydrophilic in nature. These results agree with crystallographic data at 120 K and, in particular, with the spontaneous hydration process of a cyclodextrin crystal in wet atmosphere. The predicted structure of the hydration shells is discussed and compared with previous molecular mechanics calculations which report an overall hydrophobic behavior. Moreover, the temperature dependence of the hydration process is discussed.

Structural and dynamic properties of the homodimeric hemoglobin from Scapharca inaequivalvis Thr-72-->Ile mutant: molecular dynamics simulation, low temperature visible absorption spectroscopy, and resonance Raman spectroscopy studies.

Molecular dynamics simulations, low temperature visible absorption spectroscopy, and resonance Raman spectroscopy have been performed on a mutant of the Scapharca inaequivalvis homodimeric hemoglobin, where residue threonine 72, at the subunit interface, has been substituted by isoleucine. Molecular dynamics simulation indicates that in the Thr-72-->Ile mutant several residues that have been shown to play a role in ligand binding fluctuate around orientations and distances similar to those observed in the x-ray structure of the CO derivative of the native hemoglobin, although the overall structure remains in the T state. Visible absorption spectroscopy data indicate that in the deoxy form the Soret band is less asymmetric in the mutant than in the native protein, suggesting a more planar heme structure; moreover, these data suggest a similar heme-solvent interaction in both the liganded and unliganded states of the mutant protein, at variance with that observed in the native protein. The "conformation sensitive" band III of the deoxy mutant protein is shifted to lower energy by >100 cm-1 with respect to the native one, about one-half of that observed in the low temperature photoproducts of both proteins, indicating a less polar or more hydrophobic heme environment. Resonance Raman spectroscopy data show a slight shift of the iron-proximal histidine stretching mode of the deoxy mutant toward lower frequency with respect to the native protein, which can be interpreted in terms of either a change in packing of the phenyl ring of Phe-97, as also observed from the simulation, or a loss of water in the heme pocket. In line with this latter interpretation, the number of water molecules that dynamically enters the intersubunit interface, as calculated by the molecular dynamics simulation, is lower in the mutant than in the native protein. The 10-ns photoproduct for the carbonmonoxy mutant derivative has a higher iron-proximal histidine stretching frequency than does the native protein. This suggests a subnanosecond relaxation that is slowed in the mutant, consistent with a stabilization of the R structure. Taken together, the molecular dynamics and the spectroscopic data indicate that the higher oxygen affinity displayed by the Thr-72-->Ile mutant is mainly due to a local perturbation in the dimer interface that propagates to the heme region, perturbing the polarity of the heme environment and propionate interactions. These changes are consistent with a destabilization of the T state and a stabilization of the R state in the mutant relative to the native protein.