The free energy of a reaction coordinate at multiple constraints: a concise formulation

The free energy as a function of the reaction coordinate (rc) is the key quantity for the computation of equilibrium and kinetic quantities. When it is considered as the potential of mean force, the problem is the calculation of the mean force for given values of the rc. We reinvestigate the PMCF (potential of mean constraint force) method which applies a constraint to the rc to compute the mean force as the mean negative constraint force and a metric tensor correction. The latter allows for the constraint imposed to the rc and possible artefacts due to multiple constraints of other variables which for practical reasons are often used in numerical simulations. Two main results are obtained that are of theoretical and practical interest. First, the correction term is given a very concise and simple shape which facilitates its interpretation and evaluation. Secondly, a theorem describes various rcs and possible combinations with constraints that can be used without introducing any correction to the constraint force. The results facilitate the computation of free energy by molecular dynamics simulations.     

Free energy from molecular dynamics with multiple constraints

In molecular dynamics simulations of reacting systems, the key step to determining the equilibrium constant and the reaction rate is the calculation of the free energy as a function of the reaction coordinate. Intuitively the derivative of the free energy is equal to the average force needed to constrain the reaction coordinate to a constant value, but the metric tensor effect of the constraint on the sampled phase space distribution complicates this relation. The appropriately corrected expression for the potential of mean constraint force method (PMCF) for systems in which only the reaction coordinate is constrained was published recently. Here we will consider the general case of a system with multiple constraints. This situation arises when both the reaction coordinate and the ‘hard’ coordinates are constrained, and also in systems with several reaction coordinates. The obvious advantage of this method over the established thermodynamic integration and free energy perturbation methods is that it avoids the cumbersome introduction of a full set of generalized coordinates complementing the constrained coordinates. Simulations of n-butane and n-pentane in vacuum illustrate the method.     

Reactions in viscous media: potential and free energy surfaces in solvent–solute coordinates

We propose a universal definition of reaction-specific solvent coordinate, described in terms of molecular-dynamic trajectories of the solvent and suitable for a variety of reactions in viscous solvents. We then use molecular dynamics simulations to obtain two-dimensional potential energy and free energy surfaces of model liquid phase isomerisation reactions in solvent–solute coordinates and discuss factors determining their topography.     

Vaidya spacetime in the diagonal coordinates

We have analyzed the transformation from initial coordinates (v, r) of the Vaidya metric with light coordinate v to the most physical diagonal coordinates (t, r). An exact solution has been obtained for the corresponding metric tensor in the case of a linear dependence of the mass function of the Vaidya metric on light coordinate v. In the diagonal coordinates, a narrow region (with a width proportional to the mass growth rate of a black hole) has been detected near the visibility horizon of the Vaidya accreting black hole, in which the metric differs qualitatively from the Schwarzschild metric and cannot be represented as a small perturbation. It has been shown that, in this case, a single set of diagonal coordinates (t, r) is insufficient to cover the entire range of initial coordinates (v, r) outside the visibility horizon; at least three sets of diagonal coordinates are required, the domains of which are separated by singular surfaces on which the metric components have singularities (either g ₀₀ = 0 or g ₀₀ = ∞). The energy–momentum tensor diverges on these surfaces; however, the tidal forces turn out to be finite, which follows from an analysis of the deviation equations for geodesics. Therefore, these singular surfaces are exclusively coordinate singularities that can be referred to as false fire-walls because there are no physical singularities on them. We have also considered the transformation from the initial coordinates to other diagonal coordinates (η, y), in which the solution is obtained in explicit form, and there is no energy–momentum tensor divergence.     

The increase of entropy upon release of a constraint in classical systems

Computer simulations have focused renewed attention on holonomic constraints imposed on selected coordinates for simplified effective calculations, and also for maintaining and parameterizing equilibrium states. Novel theoretical means for calculating thermodynamic potentials have recently been developed. Here, spontaneous transitions caused by the release of a constraint are treated on the same basis. To describe spontaneous changes, a parameter-like coordinate is converted into a dynamic generalized coordinate. The analysis yields, depending on the boundary conditions employed during the transition, a decrease in the free energy or, at constant energy, an increase in the entropy. This finding is in agreement with observations in computer simulation, and compatible with general theorems on the effect of interventions. The applicability of constraints to experimental setups is discussed using examples.     

Constrained reaction coordinate dynamics for systems with constraints

In the context of molecular dynamics simulations of rare events, the application of constraints on a suitable reaction coordinate has often been found useful for sampling of the free energy barrier. The efficiency of these calculations is hampered by geometrical difficulties, related to the metric factor and inertial forces. Some years ago Mulders et al. [1996, J. chem. Phys., 104, 48691 suggested a way to simplify the approach. Their idea was demonstrated shortly afterwards by Sprik and Ciccotti [1998, J. chem. Phys., 109, 77371. The present paper extends these results to vector reaction coordinate and molecular systems modelled with holonomic constraints.     

Sampling along reaction coordinates with the Wang-Landau method

The multiple range random walk algorithm recently proposed by Wang and Landau [2001, Phys. Rev. Lett., 86, 2050] is adapted to the computation of free energy profiles for molecular systems along reaction coordinates. More generally, we show how to extract partial averages in various statistical ensembles without invoking simulations with constraints, biasing potentials or unknown parameters. The method is illustrated on a model 10-dimensional potential energy surface, for which analytical results are obtained. It is then applied to the potential of mean force associated with the dihedral angle of the butane molecule in the gas phase and in carbon tetrachloride solvent. Finally, isomerization in a small rocksalt cluster, (NaF)₄, is investigated in the microcanonical ensemble, and the results are compared to those of parallel tempering Monte Carlo.     

The effect of a gravitational induction force on the stellar orbits in galaxies

This paper explores the effect of a gravitational induction force, similar to gravito-electromagnetism but of several orders of magnitude larger, on a spiral galaxy. This would provide, because of the huge amount of mass currents associated with spirals, a rather large gravitational induction force. Standard gravitational dynamics is modified by adding an anti-symmetric force tensor analogous to electromagnetism. By choosing a special free space solution to the Einstein field equations, it was found that the resulting force field would drop off by an r ⁻¹ rather than an r ⁻² rule. Thus this new induction force would grow to dominance over the large distances associated with galaxies. Any r ⁻¹ force will produce the flat rotation curves that are observed in disk galaxies as well as explain the Tully-Fisher relationship. It would also explain why the inner core of the spirals remain Newtonian. Quite surprisingly, this simple hypothesis also seems to offer explanations for many other observed phenomena as well.     

More efficient Brownian dynamics algorithms

We examine the relative efficiencies of three- algorithms for performing Brownian Dynamics simulations without many-body hydrodynamics. We compare the conventional Brownian Dynamics algorithm of Ermak (CBD), Smart Monte Carlo (SMC) which incorporates Boltzmann sampling into essentially a CBD procedure, and the Stochastic Runge Kutta (SRK) method. We show, using the repulsive potential φ(r) = ε(σ/r) ⁿ , where n = 36 and 72, that the SRK algorithm gives the most accurate short-time dynamics for the mean-square displacements. The SRK algorithm static and dynamical properties converge better with a reducing time step to the exact values, than those generated by the CBD algorithm; giving efficiency gains typically of a factor of 3–4. Both CBD and SMC have the incorrect sign for the first correction term to the mean square displacement in a time step, whereas the SRK algorithm gives essentially the exact solution to order Δt ², where Δt is the simulation time step. In fact, these correction terms are almost equal and opposite in sign. Expressions for these terms were derived in terms of the average interaction energy per particle. The force, shear and bulk stress autocorrelation functions were calculated. The average energy per particle and time correlation functions at short time have values in excess of the exact values, while the corresponding quantities for SRK are below this. This difference in behaviour can be traced back to the extent of compliance of the particle trajectories with the exact expansion of the Smoluchowski equation. The accuracy, at a given value of the time step, of the stochastic algorithms can significantly depend on the form of the interaction potential between particles. It is also demonstrated that the long time limits of various correlation functions are fairly insensitive to a particular scheme (SRK or CBD) used in the simulations. All the correlation functions have a stretched exponential region at intermediate to long times, and the values of the exponents on density and force law steepness have been determined.     

Soft-Mode Scenarios of Shear Localization: Atomic-Level Landscapes

Using an elementary case study of stress-induced structural instability in a crystal we discuss how collective response to large-strain deformation at the atomic level can give insights into the mechanistic nature of shear localization. Stability criteria and reaction pathway sampling are used to probe saddle-point configurations, both local spatial coordinates and activation energy barriers, to optimally extract the information available from molecular dynamics simulations.     

Splitting probabilities as a test of reaction coordinate choice in single-molecule experiments

To explain the observed dynamics in equilibrium single-molecule measurements of biomolecules, the experimental observable is often chosen as a putative reaction coordinate along which kinetic behavior is presumed to be governed by diffusive dynamics. Here, we invoke the splitting probability as a test of the suitability of such a proposed reaction coordinate. Comparison of the observed splitting probability with that computed from the kinetic model provides a simple test to reject poor reaction coordinates. We demonstrate this test for a force spectroscopy measurement of a DNA hairpin.     

Deuteron form factors and the tensor force

Results of a non-relativistic calculation of deuteron form factors are presented for separable potentials with and without tensor force. The tensor term in triplet state is added in such a way as to keep the values of deuteron binding energy,a _t andr _0t unaltered, so that the difference in the form factors can be regarded as the effect of tensor force only. The calculation has been performed for two different shapes of separable potentials and for three differentD-state probabilities to study their comparative effect.     

Transitions between secondary structures in isolated polyalanines

Monte Carlo simulations of gas-phase polyalanine peptides have been carried out with the Amber ff96 force field. A low-temperature structural transition takes place between the α-helix stable conformation and β-sheet structures, followed by the unfolding phase change. The transition state ensembles connecting the helix and sheet conformations are investigated by sampling the energy landscape along specific geometric order parameters as putative reaction coordinates, namely the electric dipole μ, the end-to-end distance d, and the gyration radius R_g. By performing series of shooting trajectories, the committor probabilities and their distributions are obtained, revealing that only the electric dipole provides a satisfactory transition coordinate for the α↔β interconversion. The nucleus at the transition is found to have a high helical content.     

粗粒化自由能模拟研究ABC输出蛋白质构象态转变

本文应用了粗粒化分子动力学方法模拟计算了ATP结合盒式输出蛋白沿其构象转变途径反应坐标的平均力势,这个反应坐标被定义为内门和外门质量中心距离之差. 计算得到的平均力势能很好地描述不同的向内构象态、向外构象态和阻塞构象状态,以及它们之间的转变. 粗粒化分子动力学自由能模拟显示,在向内构象态到向外构象态转变过程中,内门在外门打开之前先行关闭;反之,在向外构象态到向内构象态转变过程中,外门在内门大开之前先行关闭. 因此,在向内构象态和向外构象态两种转变过程中,都经过了阻塞构象状态. 模拟结果揭示了ATP结合盒式输出蛋白的传输单向性,这种特性在生命体系功能实现中具有十分重要的意义. 这些结果与先前晶体结构实验［Proc. Natl. Acad. Sci. USA 104,19005 (2007)］发现有根本的不同,这些实验结果显示了内外门同时打开的不合理结果. 本文通过计算模拟阐明了ABC输出蛋白构象态变化的分子机理.     

Complete Determination of Nitrogen Quadrupole and Hyperfine Tensors in an Oxovanadium Complex by Simultaneous Fitting of Multifrequency ESEEM Powder Spectra

One-dimensional C- and X-band as well as two-dimensional X-band ESEEM experiments were performed on the complex oxobis(2-methylquinolin-8-olato) vanadium(IV) in frozen solution. A¹⁴N ESEEM simulation strategy based on initial first- and second-order perturbation analysis of peak positions in orientationally selected ESEEM spectra is presented. The constraint parameters extracted enable one to reduce the number of free fitting parameters for each nitrogen from 10 to 4. These are the α, β resp. the φ, θ Euler angles of the NQI and the HFI tensor defined in the coordinate system of the axialgtensor. The local symmetry of the complex allows one to reduce the number of free parameters to two angles only. Subsequently, a grid search in the remaining Euler space produced the starting parameters for the final fit of the¹⁴N hyperfine and quadrupole tensors. The anisotropic nitrogen hyperfine interaction tensor was found to be strongly nonaxial (0.06, 0.51, −0.57) MHz with the components significantly smaller than the isotropic hyperfine constant −6.18 MHz. In contrast, the quadrupole tensor withK= 0.58 MHz is close to axial (η = 0.13). These tensors share the principal axis normal to the ligand plane (as imposed by the local symmetry). The axes in the ligand plane are, however, rotated 50° with respect to each other. The orientation of the quadrupole tensor axes correlate within 10° with the orientation of the ligand plane following from the X-ray structure.     

A simulation study of microwave field effects on a 3D orthorhombic lattice of rotating dipoles: short-range potential energy variation

Variation of the short-range potential energy of interaction of nearest dipoles in a three-dimensional (3D) orthorhombic lattice exposed to microwave electric fields is studied by means of the Langevin dynamics simulations. The global increase of the mean potential energy is typical for all the frequencies and intensities at lower temperatures, whereas separate potential energy peaks or peak chains are observed at intermediate temperatures. A simple statistical model proposed to account for the temperature dependence of the field intensity for potential energy peaks suggests the concerted collective rotation of the dipoles. The temperature dependence of the peak frequency is explained using a combination of the one-dimensional Kramers and the resonant activation theories applied to the field-driven collective rotation, with the nearly degenerate angular coordinates of the dipoles being used as a single effective coordinate.     

Efficient schemes to compute diffusive barrier crossing rates

The formulation of the classical barrier-crossing problem is reviewed in the context of numerical simulations, with the focus on barrier crossing problems where the reaction coordinate depends in a non-trivial way on the Cartesian coordinates of many particles. Often it is convenient to measure the barrier height using constrained dynamics. Such a calculation requires a knowledge of the Jacobian for the coordinate transformation between Cartesian and generalized (‘reaction’) coordinates, and it is shown that the calculation of this Jacobian can be simplified. The conventional expression for the crossing rate is found to become computationally inefficient when the barrier crossing is diffusive. An alternative formulation of the barrier-crossing rate is given that leads to much better statistical accuracy in the computed crossing rates.     

Relationship between crystalline order and melting mechanisms of solids

The mechanism for homogeneous nucleation of the liquid phase in Lennard-Jones solids is studied by combining the Landau free energy approach with some of the methodology developed to characterise transition path ensembles. The second-order bond orientational order parameter, Q ₆ which indexes the overall degree of crystalline order, is shown to provide a dynamically significant collective coordinate describing the melting process. Trajectories generated from configurations sampled in the vicinity of the maximum in the Landau free energy curve, F(Q ₆), are shown to have equal likelihood of teminating in either the solid or liquid-like free energy minima. It is also demonstrated that Q ₆ is necessary but not sufficient as a dynamical coordinate to describe melting and it is necessary to explore possiblities for additional coordinates which are critical for initiating melting. Our sudy suggests that the additional coordinates for describing the melting process would be some type of localised defect, much smaller in spatial extent than the size of the critical nucleus predicted by classical nucleation theory.