Calculating Free Energies Using a Scaled-Force Molecular Dynamics Algorithm

We propose and test a family of methods to calculate the free energy along a generalized coordinate, , based on computing the force acting on this coordinate. First, we derive a formula that connects the free energy in unconstrained simulations with the force of constraint that can be readily calculated numerically. Then, we consider two methods, which improve the efficiency of the free energy calculation by yielding uniform or nearly uniform sampling of . Both rely on modifying the force acting on . In one method, this force is replaced by a force with zero mean and is advanced quasistatically. In the second method, the force is augmented adaptively by a biasing force. We provide formulas for calculating the free energy of the unmodified system from the forces acting in these modified, non-Hamiltonian systems. Using conformational transitions in 1,2-dichloroethane as a test case, we show that both methods perform very well.     

Classical Molecular Dynamics Simulation of Kappa Squared Factor in Resonance Energy Transfer for Linear Dipole Models

Molecular dynamics simulations of linear models interacting through a dipolar Kihara intermolecular potential are presented. Molecular orientation correlations are used to calculate the orientational factor kappa squared in the resonance energy transfer (RET) as a function of the intermolecular separation. The distance, R 0 (2/3), at which the simulated systems show an isotropic behavior is calculated and an analysis of the dependence of R 0 (2/3) on microscopic properties (molecular aspect ratio and dipole moment) as well on thermodynamics (temperature and density) is presented. An explanation of the use of metallic cations as probes in RET is given and some relations of our models with biological molecules are pointed out.     

Modeling Confined Fluids: An NhPT Molecular Dynamics Method

We present an NhPT MD method developed for systematic investigation of both the structural and dynamical properties of confined fluids without resorting to chemical potential or explicit reservoir. This method allows confined fluids to expand or contract transversely and the same number of fluid molecules to be simulated throughout all surface separations. Its first implementation using confined Lennard-Jones fluid yields step-like changes in surface density, layered configurations, in-plane ordering, and oscillatory perpendicular pressures and transverse diffusivities that are consistent with previous studies. Additionally, a pseudo-Poisson's ratio and transverse isothermal compressibility were calculated. Like other properties, they oscillate at smaller surface separations and approach constant values when the surface separation becomes sufficiently large. The limiting value of the pseudo-Poisson's ratio is interestingly equivalent to that of incompressible continua.     

Numerical Calculation of the Bridge Function for Soft-Sphere Supercooled Fluids via Molecular Dynamics Simulations

Abstract

Isokinetic molecular dynamics simulations have been performed for 13,500 soft-spheres interacting through the inverse-power potential, ε([sgrave]/r)ⁿ , near and below the freezing temperature. The bridge function for the integral equation of the theory of liquids is extracted from the pair distribution function (PDF) obtained by the computer simulations for n = 6 and 12. The result is compared with that of approximate theories, i.e., the Rogers-Young (RY) approximation and a modified hypernetted-chain approximation for supercooled soft-sphere fluids (MHNCS approximation). Below the freezing temperature, the bridge function obtained by the computer simulation begins to oscillate around zero at intermediate distances where the second peak of the PDF appears. Such oscillatory behavior of the bridge function is well reproduced by the MHNCS approximation which includes correlations given by the leading elementary diagram, in remarkable contrast to that of the RY approximation. The present result suggests that the split second peak of the PDF for highly supercooled liquids is essentially dominated by the intermediate-distance-range correlation of the leading elementary diagram.     

Molecular Dynamics Simulations for Metallic Nanosystems

Nanotechnology is a crucial field for future scientific development where many different disciplines meet. Computational modelization of nanometer-sized structures is a key issue in this development because (i) it allows a considerable saving of resources and costly experimental setups intended to fabricate nanometric test devices and (ii) nowadays the study of nanometric sized systems is feasible with thoroughly designed computational codes and relatively low cost computational resources. This article describes how molecular dynamics simulations, in combination with potentials obtained in the framework of the embedded atom method, are able to describe the properties of two systems of interest for the development of future nanoelectronic devices: metallic nanowires and metallic nanofilms. Our results show that nanowire stretching results in a series of well-defined geometric structures (shells) and that thin films experiment a crystallographic phase transition for a decreasing number of layers. In both cases, good agreement with experiments is found.     

On the Accuracy Of Some Common Molecular Dynamics Algorithms

Abstract

Attention is drawn to the fact that some of the algorithms used in the simulation of molecular dynamics are less accurate than is commonly believed. In particular, we show that many of the “Verlet-equivalent” integration schemes are not equivalent to the Verlet algorithm, and consequently are not necessarily third order schemes which exhibit exact time-reversal symmetry. Of this class of algorithms, only Beeman's technique is found to generate the optimal positions and velocities for a third order technique. It is also pointed out that the method of constraints introduces errors of O(τ³) into the calculated position, and hence limits the accuracy of simulations that employ this method to second order.     

Ewald Summation in the Molecular Dynamics Simulation of Large Ionic Systems

Abstract

The accuracy and efficiency of the direct Ewald summation are discussed in terms of the size of a Molecular Dynamics (MD) ionic system and the ranges of the r-space and q-space summations. The dependence of the convergence parameter α on the size of the system and on the choice of cut-off radius for the short-range potential is given. The possibility of neglecting the q-space term for large ionic systems is discussed in terms of the accuracy and efficiency of the simulation.     

生物大分子的自由能计算方法

在本文中,我们总结了计算生物物理中广泛使用的经典的自由能计算方法,主要可以归纳为三大类:第一类是基于过渡态理论,利用路径系综得到自由能的方法;第二类方法采用一些经验势来补偿自由能表面达到均匀采样的目标;最后一种是将整个转换径分段计算,然后拼接完成得到完整自由能的方法.通过详细介绍自由能计算的重要性和实施细节,希望能够让...     

Estimation of Free Energy Systematic Errors in Molecular Simulations of Globular Proteins Surrounded by Finite Water Clusters. One Center Multipole Expansion of Reaction Field Differences

Free energy calculated in simulations on the atomic level (Monte Carlo or Molecular Dynamics) has a systematic error, if the water shell surrounding a globular protein is finite. The error (“cluster error”) is equal to a difference of free energies obtained in simulations with an infinite and finite water shell. In this work a continuum dielectric model was used to estimate the “cluster error”. A multipole expansion of the estimate was performed for a water shell with a spherical outer boundary. The expansion has very simple form. Each term is a product of two functions, one of them depending only on the charge's conformation, and the other one only on dielectric properties of the system. There are two practical uses of the expansion. First, it may be used to estimate the “cluster error” in a simulation already made; second, it may be used to plan a simulation in such a way that the “cluster error” is minimal. Numerical values of the largest terms in the multipole expansion corresponding to a typical system in simulations of globular proteins are given.     

Multiple Time Step Brownian Dynamics for Long Time Simulation of Biomolecules

We report a multiple time step algorithm applied to an atomistic Brownian dynamics simulation for simulating the long time scale dynamics of biomolecules. The algorithm was based on the original multiple time step method; a short time step was used to keep faster motions in local equilibrium. When applied to a 28-mer # # ! folded peptide, the simulation gave stable trajectories and the computation time was reduced by a factor of 160 compared to a conventional molecular dynamics simulation using explicit water molecules. We applied it for the folding simulation of a 13-mer ! -helical peptide, giving a successful folding simulation. These results indicate that the Brownian dynamics with the multiple time step algorithm is useful for studies of biomolecular motions by long time simulation.     

Molecular Dynamics and Free Energy Calculations of the B and Z Forms of C8-Arylguanine Modified Oligonucleotides

Abstract

Arylhydrazines found in the mushroom Agaricus bisporus have been shown to be carcinogenic. Upon metabolic activation, arylhydrazines are transformed into aryl radicals, forming 8-arylpurines, which may play a role in arylhydrazine carcinogenesis. These adducts are poorly read and inhibit chain extension but do alter the conformational preferences of oligonucleotides. We have shown that C8-phenylguanine modification of d(CGCGCG*CGCG) (G*= 8-phenylguanine) stabilizes it in the Z-DNA conformation (B/Z-DNA=1:1, 200 mM NaCl, pH 7.4). Here we have conducted molecular dynamics and free energy calculations to determine the sources(s) of these conformational affects and to predict the affect of the related C8- tolyl and C8-hydroxymethylphenyl guanine adducts on B/Z-DNA equilibrium. Force field parameters for the modified guanines were first developed using Guassian98 employing the B3LYP method and the standard 6–31G* basis set and fit to the Cornell 94 force field with RESP. Molecular dynamics simulations and free energy calculations, using the suite of programs contained in Amber 6 and 7 with the Cornell 94 force field, were used to determine the structural and thermodynamic properties of the DNA. The principal factors that drive conformation are stacking of the aryl group over the 5′-cytosine in the phenyl and tolyl modified oligonucleotides while hydrogen bonding opposes stacking in the hydroxymethylphenyl derivative. The phenyl and tolyl-modified DNA's favored the Z-DNA form as did the hydroxymethylphenyl derivative when hydrogen bonding was not present. The B-DNA conformation was preferred by the unmodified oligonucleotide and by the hydroxymethylphenyl-modified oligonucleotide when hydrogen bonding was considered. Z-DNA stability was not found to directly correlated with carcinogenicity and additional biological factors, such as recognition and repair, may also need to be considered in addition to Z-DNA formation.     

Helix Unfolding/Refolding Characterizes the Functional Dynamics of Staphylococcus aureus Clp Protease

The ATP-dependent Clp protease (ClpP) plays an essential role not only in the control of protein quality but also in the regulation of bacterial pathogen virulence, making it an attractive target for antibacterial treatment. We have previously determined the crystal structures of Staphylococcus aureus ClpP (SaClpP) in two different states, extended and compressed. To investigate the dynamic switching of ClpP between these states, we performed a series of molecular dynamics simulations. During the structural transition, the long and straight helix E in the extended SaClpP monomer underwent an unfolding/refolding process, resulting in a kinked helix very similar to that in the compressed monomer. As a stable intermediate in the molecular dynamics simulation, the compact state was suggested and subsequently identified in x-ray crystallographic experiment. Our combined studies also determined that Ala¹⁴⁰ acted as a “hinge” during the transition between the extended and compressed states, and Glu¹³⁷ was essential for stabilizing the compressed state. Overall, this study provides molecular insights into the dynamics and mechanism of the functional conformation changes of SaClpP. Given the highly conserved sequences of ClpP proteins among different species, these findings potentially reflect a switching mechanism for the dynamic process shared in the whole ClpP family in general and thus aid in better understand the principles of Clp protease assembly and function.     

Molecular Dynamics Simulation of Hexamine and Suberic Acid

In order to perform a molecular dynamics (MD) simulation of the incommensurate crystalline structure hexamethylenetetramine suberate (C 6 H 12 N 4 )(HOOC-(CH 2 ) 6 -COOH), we present in a first step the separate simulations of the crystalline structure of each of the two pure components, hexamethylenetetramine (HMT) and suberic acid. The domain decomposition parallel MD program ddgmq is used for this purpose. A second-generation consistent force field (CFF91) is employed to describe the interactions between atoms. Starting from experimental crystal structures, both pure components were heated from low to high temperatures. Our MD results show that the HMT system can be well represented by CFF91. In the case of suberic acid the layered structure of the crystal is largely preserved although deviations in the unit cell lengths from the experimental values are ～10%. Rather than attempt a complex re-parametrisation of CFF91 we chose to impose a fixed compensating external pressure tensor to correct for the deficiencies of the chosen force field. After optimising this compensating external pressure tensor at one temperature we find that experimental lattice constants and angles can be well reproduced over a range of temperatures.     

Molecular Dynamics Simulation of Collective Motions in Binary Liquids

Collective dynamic properties of different kind of binary liquid mixtures have been investigated by molecular dynamics simulation. The study includes both the longitudinal and the transverse current spectra in simple liquid alloys, 1:1 molten salts and liquid binary mixtures of neutral particles with an ionic-like structure. These systems were chosen as representative of binary liquids with different static structures in order to analyse the effects of structural ordering on the mechanisms of dynamic collective properties. The effect of the mass asymmetry between the two species in the mixture has been also discussed from the results for two different mass ratios for each kind of structure. Two length scales have been considered. On the one hand, the hydrodynamic scale (low wave numbers), where the modes for the partial currents of the two species are characterised by very close frequencies. On the other hand, the molecular scale (higher wave numbers), where the characteristic frequencies for the two species show noticeable differences. Vibrational concentration current modes (optic modes) have been found in neutral mixtures though their influence is rather weak, being the collective dynamic properties of this kind of systems dominated by the mass current modes (acoustic modes). On the contrary, in mixtures of charged particles such as molten salts the contribution of the concentration (charge) currents to the collective dynamics is important and optic modes can be characterised by a well-defined frequency for a wide range of wave numbers. It has been observed that heavy particles have a more relevant role on the mass current correlations whereas light particles play a dominant role on the concentration current correlations. The overall results for the three kinds of liquid mixtures analysed in this paper show that both the longitudinal and transverse current spectra are little dependent on the static structure of the system whereas marked differences are revealed when the particles in the system are either neutral or carry an electric charge.     

Molecular Dynamics Simulation of Model Lipid Membranes: Structural Effects of Impurities

Abstract

The gel to fluid phase transition or ordered to disordered phase transition observed in biological membranes are simulated by using constant energy Molecular Dynamics. The surface part of the membrane is modelled as a two-dimensional matrix formed by the head groups of the phospholipid molecules. Head molecules which are modelled as three spheres fused with three force centers, interact with each other via van der Waals and Coulomb type interactions. The -so called- impurity or foreign molecule embedded in the surface represents the protein type molecule which is present in biological membranes and control its activity. It is modelled as a pentagon having one force centers in each corner. It also interacts with the surface molecules again via van der Waals and Coulomb type interactions. The surface density is kept constant in the simulations of the systems with or without impurity. Structural and orientational changes due to impurity were observed and proved by monitoring two-dimensional order parameter. It has been shown that melting of the surface or breakage of the ordering of the surface molecules becomes easier and ordered to disordered phase transition temperature was lowered by 100 K if the impurity is present.     

Molecular Dynamics Simulations of a Nucleosome and Free DNA

Abstract

All atom molecular dynamics simulations (10ns) of a nucleosome and of its 146 basepairs of DNA free in solution have been conducted. DNA helical parameters (Roll, Tilt, Twist, Shift, Slide, Rise) were extracted from each trajectory to compare the conformation, effective force constants, persistence length measures, and fluctuations of nucleosomal DNA to free DNA. The conformation of DNA in the nucleosome, as determined by helical parameters, is found to be largely within the range of thermally accessible values obtained for free DNA. DNA is found to be less flexible on the nucleosome than when free in solution, however such measures are length scale dependent. A method for disassembling and reconstructing the conformation and dynamics of the nucleosome using Fourier analysis is presented. Long length variations in the conformation of nucleosomal DNA are identified other than those associated with helix repeat. These variations are required to create a proposed tetrasome conformation or to qualitatively reconstruct the 1.75 turns of the nucleosome's superhelix. Reconstruction of free DNA using selected long wavelength variations in conformation can produce either a left-handed or a right-handed superhelix. The long wavelength variations suggest 146 basepairs is a natural length of DNA to wrap around the histone core.     

An Iterative Variable-timestep Algorithm for Molecular Dynamics Simulations

A method for performing variable-timestep molecular dynamics integration is described, in which an iterative algorithm is used to select the largest timestep consistent with the desired simulation accuracy. Accuracy in this context is defined in terms of energy conservation, rather than trajectory correctness. Specifically, a timestep-independent measure of the rate of "diffusion" of the total energy is used. This variable timestep approach is compared to fixed-timestep integration for three different hydrocarbon systems (polyethylene, liquid benzene and ethylene), which are modeled with a reactive bond-order potential. These systems represent both equilibrium and highly non-equilibrium systems at temperatures ranging from 298 to 2500 K. The variable-timestep method is found to be approximately twice as computationally efficient as fixed-timestep integration for the non-equilibrium sputtering of polyethylene, and the two methods were competitive for the equilibrium systems. The algorithm requires the specification of two parameters controlling the rates of timestep growth and decay, but it is found that one set of values is appropriate for all three systems studied, and there is reason to believe that the parameters are transferable to other systems. The algorithm was developed specifically for simulations involving disparate timescales, such as are encountered with the reactive bond-order model used here, but it should also prove beneficial for a wide range of molecular dynamics applications.     

A Molecular Dynamics Investigation of Compressed Aqueous Alkanoate Monolayers

Molecular dynamics are conducted on a dodecanoic acid monolayer/aqueous surface. Surface pressure is controlled by imposing constant-volume conditions for series of lengths of the square slab constituting the MD cell. The response of the alkanoate chains to the pressure is followed by examining various computed quantities that monitor their conformational order. These include atom-pair radial distribution functions, chain torsional angles, energies, atomic densities perpendicular to the interface, diffusivities and atomic plots. These quantities lead to chain separations which in the range 4-5 Å implying order when the alkanoate chains have a mean area of 0.18 nm 2.     

Theoretical Investigation on the Binding Specificity of Sialyldisaccharides with Hemagglutinins of Influenza A Virus by Molecular Dynamics Simulations

Recognition of cell-surface sialyldisaccharides by influenza A hemagglutinin (HA) triggers the infection process of influenza. The changes in glycosidic torsional linkage and the receptor conformations may alter the binding specificity of HAs to the sialylglycans. In this study, 10-ns molecular dynamics simulations were carried out to examine the structural and dynamic behavior of the HAs bound with sialyldisaccharides Neu5Acα(2–3)Gal (N23G) and Neu5Acα(2–6)Gal (N26G). The analysis of the glycosidic torsional angles and the pair interaction energy between the receptor and the interacting residues of the binding site reveal that N23G has two binding modes for H1 and H5 and a single binding mode for H3 and H9. For N26G, H1 and H3 has two binding modes, and H5 and H9 has a single binding mode. The direct and water-mediated hydrogen bonding interactions between the receptors and HAs play dominant roles in the structural stabilization of the complexes. It is concluded from pair interaction energy and Molecular Mechanic-Poisson-Boltzmann Surface Area calculations that N26G is a better receptor for H1 when compared with N23G. N23G is a better receptor for H5 when compared with N26G. However, H3 and H9 can recognize N23G and N26G in equal binding specificity due to the marginal energy difference (≈2.5 kcal/mol). The order of binding specificity of N23G is H3 > H5 > H9 > H1 and N26G is H1 > H3 > H5 > H9, respectively. The proposed conformational models will be helpful in designing inhibitors for influenza virus.     

Molecular Dynamics Simulation of Limiting Conductance for Li + Ion in Supercritical Water using Polarizable Models

We report results of molecular dynamics simulations of the limiting conductance of Li + ion in ambient water and in supercritical water using polarizable models for water and Li + . The limiting conductances of Li + in ambient water calculated from mean square displacement (MSD) using four points transferable intermolecular potential model (TIP4P), extended simple point charge model (SPC/E), and revised polarizable model 1 (RPOL1) are larger than the experimental value. The behavior of the limiting conductance of Li + in supercritical water using the RPOL models results in good agreement with experimental results for the limiting conductance of LiCl. The agreement of the RPOL1 model with the experimental results is much better than the RPOL2 model in the higher-density regime, whereas that of the RPOL2 model is much better than the RPOL1 model in the lower-density regime. Using the RPOL models (in contrast to the SPC/E model), the number of hydration water molecules around Li + is the dominating contributor to the limiting conductance in the higher-density regime. In agreement with the SPC/E model, the interaction strength between Li + and the hydration water molecules is a non-factor in the lower-density region since the potential energy per hydration water molecule decreases with decreasing water density at the lowest water densities.