Thermal motion from monte carlo simulations of aqueous ionic solutions

The statistical pattern recognition procedure of Marchese and Beveridge for the analysis of Monte Carlo simulations of aqueous ionic solutions [J. Am. Chem. Soc. 106 , 3713 (1984)] has been extended to include average hydrogen positions as well as oxygen positions. In addition, thermal ellipsoids have been calculated for each atom and displayed graphically. Application of this procedure to the analysis of a dilute aqueous solution of Zn⁺⁺ reveals an octahedral arrangement of water molecules within the ion's first hydration shell. The thermal ellipsoids show that most of the water motion is manifest in the hydrogen atoms.     

On the efficiency of exchange in parallel tempering monte carlo simulations

We introduce the concept of effective fraction, defined as the expected probability that a configuration from the lowest index replica successfully reaches the highest index replica during a replica exchange Monte Carlo simulation. We then argue that the effective fraction represents an adequate measure of the quality of the sampling technique, as far as swapping is concerned. Under the hypothesis that the correlation between successive exchanges is negligible, we propose a technique for the computation of the effective fraction, a technique that relies solely on the values of the acceptance probabilities obtained at the end of the simulation. The effective fraction is then utilized for the study of the efficiency of a popular swapping scheme in the context of parallel tempering in the canonical ensemble. For large dimensional oscillators, we show that the swapping probability that minimizes the computational effort is 38.74%. By studying the parallel tempering swapping efficiency for a 13-atom Lennard-Jones cluster, we argue that the value of 38.74% remains roughly the optimal probability for most systems with continuous distributions that are likely to be encountered in practice.     

Calculation of the entropy and free energy of peptides by molecular dynamics simulations using the hypothetical scanning molecular dynamics method

Hypothetical scanning (HS) is a method for calculating the absolute entropy S and free energy F from a sample generated by any simulation technique. With this approach each sample configuration is reconstructed with the help of transition probabilities (TPs) and their product leads to the configuration's probability, hence to the entropy. Recently a new way for calculating the TPs by Monte Carlo (MC) simulations has been suggested, where all system interactions are taken into account. Therefore, this method--called HSMC--is in principle exact where the only approximation is due to insufficient sampling. HSMC has been applied very successfully to liquid argon, TIP3P water, self-avoiding walks on a lattice, and peptides. Because molecular dynamics (MD) is considered to be significantly more efficient than MC for a compact polymer chain, in this paper HSMC is extended to MD simulations as applied to peptides. Like before, we study decaglycine in vacuum but for the first time also a peptide with side chains, (Val)(2)(Gly)(6)(Val)(2). The transition from MC to MD requires implementing essential changes in the reconstruction process of HSMD. Results are calculated for three microstates, helix, extended, and hairpin. HSMD leads to very stable differences in entropy TDeltaS between these microstates with small errors of 0.1-0.2 kcal/mol (T=100 K) for a wide range of calculation parameters with extremely high efficiency. Various aspects of HSMD and plans for future work are discussed.     

Calculation of free energy differences for water from computer simulations

Free energy differences for water at different temperatures have been calculated from Monte Carlo computer simulations using both ratio overlap and umbrella sampling methods. The problems of calculating precise values from these methods are discussed. Three models for water (ST2, TIPS2 and PE) have been used in these calculations which have been compared with experimental estimates for the configurational free energy. The importance of being able to predict accurate mean potential energies as well as accurate energy distributions is emphasised.     

Calculation of weight function for reactive differential cross sections by monte carlo methods

A weight function which depends only on the center-of-mass (c.m.) angle and velocity is used to relate the c.m. differential cross section to a measured reactive scattering contour map. The weight function is calculated by Monte Carlo methods. Its calculation includes the viewing factor correction.     

Thermodynamic characterization of fluids confined in heterogeneous pores by monte carlo simulations in the grand canonical and the isobaric-isothermal ensembles

Materials presenting nanoscale porosity are able to condense gases in their structure. This "capillary condensation" phenomenon has been studied for more than one century. Theoretical models help to understand experimental results but fail in explaining all experimental features. Most of the time, the difficulties in making quantitative or even qualitative predictions are due to the geometric complexity of the porous materials, such as large pore size distribution, chemical heterogeneities, or pore interconnections. Numerical calculations (lattice gas models or molecular simulations) are of considerable interest to calculate the adsorption properties of a fluid confined in a porous model with characteristic sizes up to several tens of nanometers. For instance, the grand canonical Monte Carlo method allows one to compute the average amount of fluid adsorbed in the porous model as a function of the temperature and the chemical potential of the fluid. However, the grand potential, necessary for a complete characterization of the system, is not a direct output of the algorithm. It is shown in this paper that the use of the isobaric-isothermal (NPT) ensemble allows one to circumvent this problem; that is, it is possible to get in one single Monte Carlo run the absolute grand potential for any given thermodynamic state of the fluid. A simplified thermodynamic integration scheme is then used to evaluate the grand potential over the whole isotherm branch passing through this initially given point. Since the usual NPT technique is a priori limited to homogeneous pores, it is proposed, for the first time, to generalize this procedure to a pore presenting a chemical heterogeneity along its axis. The new method gives the same results as the previous for homogeneous pores and allows new predictions for chemically heterogeneous pores. Comparison with the full integration scheme shows that the proposed direct calculation is faster since it avoids multiple Monte Carlo runs and more precise because it avoids the possible cumulative errors of the integration procedure.     

Fluorescence correlation spectroscopy simulations of photophysical phenomena and molecular interactions: a molecular dynamics/monte carlo approach

Fluorescence correlation spectroscopy (FCS) is being applied increasingly to study diffusion and interactions of fluorescently labeled macromolecules in complex biological systems. Fluctuations in detected fluorescence, deltaF(t), are expressed as time-correlation functions, G(tau), and photon-count histograms, P(k;DeltaT). Here, we developed a generalized simulation approach to compute G(tau) and P(k;DeltaT) for complex systems with arbitrary geometry, photophysics, diffusion, and macromolecular interactions. G(tau) and P(k;DeltaT) were computed from deltaF(t) generated by a Brownian dynamics simulation of single-molecule trajectories followed by a Monte Carlo simulation of fluorophore excitation and detection statistics. Simulations were validated by comparing analytical and simulated G(tau) and P(k;DeltaT) for diffusion of noninteracting fluorophores in a three-dimensional Gaussian excitation and detection volume. Inclusion of photobleaching and triplet-state relaxation produced significant changes in G(tau) and P(k;DeltaT). Simulations of macromolecular interactions and complex diffusion were done, including transient fluorophore binding to an immobile matrix, cross-correlation analysis of interacting fluorophores, and anomalous sub- and superdiffusion. The computational method developed here is generally applicable for simulating FCS measurements on systems complicated by fluorophore interactions or molecular crowding, and experimental protocols for which G(tau) and P(k;DeltaT) cannot be computed analytically.     

Adsorption mechanism of carbon dioxide in faujasites: grand canonical monte carlo simulations and microcalorimetry measurements

Molecular simulations have been coupled with adsorption microcalorimetry measurements in order to understand more deeply the interactions between carbon dioxide and various types of faujasite surfaces. The modeling studies, based on newly derived interatomic potentials for describing the interactions within the whole system, provide isotherms and evolutions of the differential enthalpy of adsorption as a function of coverage for DAY, NaY, and NaLSX which are in very good accordance with those obtained experimentally. The microscopic mechanism of CO2 adsorption was carefully analyzed, with different behaviors proposed, depending on the energetic characteristics of each faujasite surface, which are consistent with the trends observed for the differential enthalpies of adsorption.     

The salt dependence of the preferential interaction coefficient in dna solutions as determined by grand canonical monte carlo simulations

Grand canonical Monte Carlo simulations have been used to calculate the preferential interaction coefficient for aqueous solutions of sodium DNA. A cell model representation was used in these simulations, which were performed over a range of DNA phosphate (1–31 mM) and NaCl (0.2–111 mM) concentrations. Results are in close agreement with the values of the preferential interaction coefficient calculated using the Poisson-Boltzmann cell model, and with published values of experimental Donnan coefficients. The high degree of non-ideality indicated by the calculated preferential interaction coefficients is interpreted in terms of the classical Donnan equilibrium and the concept of thermodynamic association.     

Absolute hydration free energy scale for alkali and halide ions established from simulations with a polarizable force field

A polarizable potential function for the hydration of alkali and halide ions is developed on the basis of the recent SWM4-DP water model [Lamoureux, G.; MacKerell, A. D., Jr.; Roux, B. J. Chem. Phys. 2003, 119, 5185]. Induced polarization is incorporated using classical Drude oscillators that are treated as auxiliary dynamical degrees of freedom. The ions are represented as polarizable Lennard-Jones centers, whose parameters are optimized to reproduce the binding energies of gas-phase monohydrates and the hydration free energies in the bulk liquid. Systematic exploration of the parameters shows that the monohydrate binding energies can be consistent with a unique hydration free energy scale if the computed hydration free energies incorporate the contribution from the air/water interfacial electrostatic potential (-540 mV for SWM4-DP). The final model, which can satisfyingly reproduce both gas and bulk-phase properties, corresponds to an absolute scale in which the intrinsic hydration free energy of the proton is -247 kcal/mol.     

Calculation of the entropy and free energy by the hypothetical scanning Monte Carlo method: application to peptides

A new approach, the hypothetical scanning Monte Carlo (HSMC), for calculating the absolute entropy, S, and free energy, F, has been introduced recently and applied first to fluids (argon and water) and later to peptides. In this paper the method is further developed for peptide chains in vacuum. S is calculated from a given MC sample by reconstructing each sample conformation i step-by-step, i.e., calculating transition probabilities (TPs) for the dihedral and bond angles and fixing the related atoms at their positions. At step k of the process the chain's coordinates that have already been determined are kept fixed (the "frozen past") and TP(k) is obtained from a MC simulation of the "future" part of the chain whose TPs as yet have not been determined; when the process is completed the contribution of conformation i to the entropy is, S(i) approximately -ln Pi(k) TP(k). In a recent paper we studied polyglycine chains, modeled by the AMBER force field with constant bond lengths and bond angles (the rigid model). Decaglycine [(Gly)(10)] was studied in the helical, extended, and hairpin microstates, while (Gly)(16) was treated only in the first two microstates. In this paper the samples are increased and restudied, (Gly)(16) is also investigated in the hairpin microstate, and for (Gly)(10) approximations are tested where only part of the future is considered for calculating the TPs. We calculate upper and lower bounds for F and demonstrate that like for fluids, F can be obtained from multiple reconstructions of a single conformation. We also test a more realistic model of (Gly)(10) where the bond angles are allowed to move (the flexible model). Very accurate results for S and F are obtained which are compared to results obtained by the quasiharmonic approximation and the local states method. Thus, differences in entropy and free energy between the three microstates are obtained within errors of 0.1-0.3 kcal/mol. The HSMC method can be applied to a macromolecule with any degree of flexibility, ranging from local fluctuations to a random coil. The present results demonstrate that the difference in stability, DeltaF(mn)=F(m)-F(n) between significantly different microstates m and n, can be obtained from two simulations only without the need to resort to thermodynamic integration. Our long-term goal is to extend this method to any peptide and apply it to a peptide immersed in a box with explicit water.     

Statics and dynamics of dense polymer systems studied by monte carlo simulation

Monte Carlo simulations of coarse–grained models of macromolecules offer a unique tool to study the interplay between coil conformations, thermodynamic properties, and chain configurational relaxation and diffusion. Two examples are discussed where the chain conformation strongly differs from a gaussian coil: (i) collapsed chains in a bad solvent, where anomalous diffusion occurs in the Rouse limit and the relaxation time increases at least with the third power of chain length. (ii) Expulsion of a chain from a semidilute polymer brush. The initially stretched chain contracts to a gaussian coil and the center of mass moves outward with constant velocity until it reaches the region of the “last blob” where crossover to diffusive behavior occurs.     

Thermodynamic study of alkali sulphates: Calculation of enthalpy,entropy and Gibbs free energy

Joined with measures of specific heat at low temperature , a calorimetric study of alkali sulphates from ordinary temperature to 1500 K has allowed the calculation of enthalpy, entropy and Gibbs free energy of these salts. The value of entropy at melting point is compared with the value assessed by an acoustical method.     

Calculation of ligand binding free energies from molecular dynamics simulations

A recently developed method for predicting binding affinities in ligand–receptor complexes, based on interaction energy averaging and conformational sampling by molecular dynamics simulation, is presented. Polar and nonpolar contributions to the binding free energy are approximated by a linear scaling of the corresponding terms in the average intermolecular interaction energy for the bound and free states of the ligand. While the method originally assumed the validity of electrostatic linear response, we show that incorporation of systematic deviations from linear response derived from free energy perturbation calculations enhances the accuracy of the approach. The method is applied to complexes of wild-type and mutant human dihydrofolate reductases with 2,4-diaminopteridine and 2,4-diaminoquinazoline inhibitors. It is shown that a binding energy accuracy of about 1 kcal/mol is attainable even for multiply ionized compounds, such as methotrexate, for which electrostatic interactions energies are very large. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 69: 77–88, 1998     

Calculation of phase coexistence properties and surface tensions of n-alkanes with grand-canonical transition-matrix monte carlo simulation and finite-size scaling

Grand-canonical transition-matrix Monte Carlo is combined with configurational-bias and expanded ensemble Monte Carlo techniques to obtain saturated densities and vapor pressures of select n-alkanes. Surface tension values for butane, hexane, and octane are also computed via the finite-size scaling method of Binder. The exponential-6 model of Errington and Panagiotopoulos is used to describe the molecular interactions. The effect of the number of configurational-bias trial conformations on the efficiency of phase equilibra calculations is studied. We find that a broad range of trial conformation numbers give reasonable performance, with the optimal value increasing with decreasing temperature for a fixed chain length. Phase coexistence properties are in good agreement with literature values and are obtained with very reasonable computing resources. Similar to other recently developed n-alkane force fields, the exponential-6 model overestimates the surface tension relative to experimental values. Statistical uncertainties for coexistence properties obtained with the current approach are relatively small compared to existing methods.     

Calculation of solvation free energy using RISM theory for peptide in salt solution

We developed a robust, highly efficient algorithm for solving the full reference interaction site model (RISM) equations for salt solutions near a solute molecule with many atomic sites. It was obtained as an extension of our previously reported algorithm for pure water near the solute molecule. The algorithm is a judicious hybrid of the Newton–Raphson and Picard methods. The most striking advantage is that the Jacobian matrix is just part of the input data and need not be recalculated at all. To illustrate the algorithm, we solved the full RISM equations for a dipeptide (NH₂(SINGLE BOND)CHCH₃(SINGLE BOND)CONH(SINGLE BOND)CHCH₃(SINGLE BOND)COOH) in a 1 M NaCl solution. The extended simple point charge (SPC/E) model was employed for water molecules. Two different conformations of the dipeptide were considered. It was assumed for each conformation that the dipeptide was present either as an un-ionized form or as a zwitterion. The structure of the salt solution near the dipeptide and salt effects on the solvation free energy were also discussed. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1724–1735, 1998     

Nucleation free energy of pore formation in an amphiphilic bilayer studied by molecular dynamics simulations

The formation of a pore in a membrane requires a considerable rearrangement of the amphiphilic molecules about to form the bilayer edge surrounding the pore, and hence is accompanied by a steep increase of the free energy. Recent rupture and conductance experiments suggest that this reshuffling process is also responsible for a small energy barrier that stabilizes "prepores" with diameters of less than 1 nm, rendering both the opening and closing of pores an activated process. We use the potential of mean constraint force method to study this free energy profile, as a function of pore radius, in a coarse grained bilayer model. The calculations show that the free energy rises by (15-20) kT during pore opening, making it an extremely rare nucleation event. Although we do not observe a barrier to pore closure, the results do make the existence of such a barrier plausible. For larger pores we find a smooth transition to Litster's model, from which a line tension coefficient of about 3.7 x 10(-11) J m(-1) is deduced.     

A monte carlo model of network formation by step-growth: Star-shaped structures developed at an isolated crosslink site

A Monte Carlo lattice model of step growth network formation from systems of divalent and multivalent molecules has been used to investigate the manner in which the latter function as nodes when diluted considerably be the former. Star molecules develop in the simulation, giving nodes that have pendent tails and loops. We show how the number and size of tails and loops is rather insensitive to the valency, f, of the nodes, and examine the distribution functions for the numbers of tails and loops present in the systems at the end of the reaction. Topological considerations created by a high number of valencies that might be expected to enhance molecule size are frustrated by the presence of adjacent segments of chains from developing the full structural potential.