SHAKE,rattle, and roll: Efficient constraint algorithms for linked rigid bodies

We present an iterative constraint algorithm, QSHAKE, for use with semirigid molecules in molecular dynamics simulations. The algorithm is based on “SHAKE-ing” bond constraints between rigid bodies, whose equations of motion are solved in the quaternion framework. The algorithm is derived and its performance compared with SHAKE for liquid octane. QSHAKE is significantly more efficient whenever SHAKE requires triangles (or tetrahedra) of constraints to maintain molecular shape. Efficiencies are gained because QSHAKE reduces the number of holonomic constraints that must be solved iteratively and requires fewer iterations to obtain convergence. The gains in efficiency are most noticeable when a high degree of precision is imposed on the constraint criteria. QSHAKE is also stable at larger time steps than SHAKE, thus allowing for even more efficient exploration of phase space for semirigid molecules. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 102–111, 1998     

MILCH SHAKE: An efficient method for constraint dynamics applied to alkanes

We describe a method to impose constraints in a molecular dynamics simulation. A technique developed to solve the special case of a linear topology (MILC SHAKE) is hybridized with the SHAKE algorithm. The methodology, which we term MILC‐hybridized SHAKE (or MILCH SHAKE), applies to more complex topologies. Here we consider the important case of all atom models of alkanes. Exploiting the mass difference between carbon and hydrogen we show that for higher alkanes MILCH SHAKE can be an order of magnitude faster than SHAKE. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

A method for constrained energy minimization of macromolecules

Two algorithms for the local energy minimization of the structure of macromolecules in the presence of constraints are proposed. They are a combination of the method of steepest descents and the method of conjugate gradients with the procedure SHAKE, by which distance constraints can be satisfied. The two algorithms are tested by applying them to a small protein, the bovine pancreatic trypsin inhibitor (BPTI), and compared with the penalty function method for conserving constraints. The efficiency of the proposed methods depends on the level of interdependence of the constraints. For bond-length constraints, the use of SHAKE is superior to the penalty function method. However, when bond-angle constraints are included, SHAKE is more efficient only if the curvature of the penalty function is considerably greater than that of the potential function being minimized. The results indicate that with bond-length constraints the minimization behavior is similar to that without constraints. However, the simultaneous application of bond-length and bond-angle constraints appears to confine the molecule to a very limited part of configuration space, very different from the part covered by an unconstrained minimization. This conclusion calls into question energy minimizations of protein systems in which only the dihedral angles are allowed to vary.     

Algorithms for constrained molecular dynamics

In molecular dynamics simulations, the fastest components of the potential field impose severe restrictions on the stability and hence the speed of computational methods. One possibility for treating this problem is to replace the fastest components with algebraic length constraints. In this article the resulting systems of mixed differential and algebraic equations are studied. Commonly used discretization schemes for constrained Hamiltonian systems are discussed. The form of the nonlinear equations is examined in detail and used to give convergence results for the traditional nonlinear solution technique SHAKE iteration and for a modification based on successive overrelaxation (SOR). A simple adaptive algorithm for finding the optimal relaxation parameter is presented. Alternative direct methods using sparse matrix techniques are discussed. Numerical results are given for the new techniques, which have been implemented in the molecular modeling software package CHARMM and show as much as twofold improvement over SHAKE iteration. © 1995 John Wiley & Sons, Inc.     

LINCS: A linear constraint solver for molecular simulations

In this article, we present a new LINear Constraint Solver (LINCS) for molecular simulations with bond constraints. The algorithm is inherently stable, as the constraints themselves are reset instead of derivatives of the constraints, thereby eliminating drift. Although the derivation of the algorithm is presented in terms of matrices, no matrix matrix multiplications are needed and only the nonzero matrix elements have to be stored, making the method useful for very large molecules. At the same accuracy, the LINCS algorithm is three to four times faster than the SHAKE algorithm. Parallelization of the algorithm is straightforward. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1463–1472, 1997     

Stochastic generator of chemical structure. 4. Building polymeric systems with specified properties

A novel computational technique to generate close‐to‐equilibrium crosslinked polymeric systems is proposed. Compared to the current state‐of‐the‐art equilibration methods, the new technique appears to be faster by several orders of magnitude. The main advantage of the technique is that one can circumvent the bottlenecks in configuration space that inhibit relaxation in molecular dynamics or Monte Carlo simulations. The problem of polymer equilibration described by continuous equations in molecular dynamics is reduced to a discrete representation where solutions are approximated by simple algorithms. In the current study, a series of coarse‐grained, united‐atom, and fully atomistic crosslinked networks has been generated. Network statistics and topology, X‐ray scattering intensities, and elastic properties are tested vs. experimental results and similar models generated using molecular dynamics and Monte Carlo simulations. The results demonstrate the efficiency of this new method for generating large realistic polymeric systems up to 1.4 M atoms. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 580–590, 2001     

Collective motion artifacts arising in long‐duration molecular dynamics simulations

We tested a variety of molecular dynamics simulation strategies in long‐duration (up to several nanoseconds) constant‐temperature simulations of liquid water under periodic boundary conditions. Such long durations are necessary to achieve adequate conformational sampling in simulations of membrane assemblies and other large biomolecular systems. Under a variety of circumstances, serious artifacts arise in the form of spurious collective behavior that becomes obvious only after the simulation has gone at least several hundred picoseconds. The potential energy of the system drops and the system changes from a liquid to an icy or glassy state. The underlying cause is accumulated center‐of‐mass motion of the system, coupled with velocity rescaling associated with constant‐temperature control. The velocity rescaling in the constant‐temperature algorithm reduces the thermal velocity as the net center‐of‐mass velocity grows, effectively causing the kinetic energy of the system to drain from thermal motions into coordinated motions. We found that the incidence and magnitude of the underlying artifactual motion leading to the spurious transition is mediated by: choice of method for computing electrostatic interactions; choice of ensemble; size of the simulation cell; SHAKE tolerance; frequency of nonbonded pairlist updating; and closeness of coupling to the temperature bath. The appearance of the spurious transition can be avoided by periodically subtracting net center‐of‐mass motion during the dynamics, or by improving the accuracy of the simulation by means of tightening SHAKE tolerance and updating nonbonded pairlists every timestep. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 121–131, 2000     

Constraint dynamics algorithm for simulation of semiflexible macromolecules

Semiflexible models are often used to study macromolecules containing stable structural elements. Based on rigid body dynamics, we developed a rigid fragment constraint dynamics algorithm for the simulation of semiflexible macromolecules. Stable structural elements are treated as rigid fragments. Rigid fragment constraints, defined as combinations of distance constraints and position constraints, are introduced to limit internal molecular motion to the required mode. The constraint forces are solved separately for each rigid fragment constraint and iteratively until all constraint conditions are satisfied within a given tolerance at each time step, as is done for the bond length constraint in the SHAKE algorithm. The orientation of a rigid fragment is represented by the quaternion parameters, and both translation and rotation are solved by the leap-frog formulation. We tested the algorithm with molecular dynamics simulations of a series of peptides and a small protein. The computation cost for the constraints is roughly proportional to the size of the molecule. In the microcanonical ensemble simulation of polyvalines, the total energy was conserved satisfactorily with time steps as large as 20 fs. A helix folding simulation of a synthetic peptide was carried out to show the efficiency of the algorithm in a conformational search. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1555–1566, 1998     

Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models

An analytical algorithm, called SETTLE, for resetting the positions and velocities to satisfy the holonomic constraints on the rigid water model is presented. This method is still based on the Cartesian coordinate system and can be used in place of SHAKE and RATTLE. We implemented this algorithm in the SPASMS package of molecular mechanics and dynamics. Several series of molecular dynamics simulations were carried out to examine the performance of the new algorithm in comparison with the original RATTLE method. It was found that SETTLE is of higher accuracy and is faster than RATTLE with reasonable tolerances by three to nine times on a scalar machine. Furthermore, the performance improvement ranged from factors of 26 to 98 on a vector machine since the method presented is not iterative. © 1992 by John Wiley & Sons, Inc.     

Free energy calculations using flexible-constrained, hard-constrained and non-constrained molecular dynamics simulations.

A comparison of different treatments of bond-stretching interactions in molecular dynamics simulation is presented. Relative free energies from simulations using rigid bonds maintained with the SHAKE algorithm, using partially rigid bonds maintained with a recently introduced flexible constraints algorithm, and using fully flexible bonds are compared in a multi-configurational thermodynamic integration calculation of changing liquid water into liquid methanol. The formula for the free energy change due to a changing flexible constraint in a flexible constraint simulation is derived. To allow for a more direct comparison between these three methods, three different pairs of models for water and methanol were used: a flexible model (simulated without constraints and with flexible constraints), a rigid model (simulated with standard hard constraints), and an alternative flexible model (simulated with flexible constraints and standard hard constraints) in which the ideal or constrained bond lengths correspond to the average bond lengths obtained from a short simulation of the unconstrained flexible model. The particular treatment of the bonds induces differences of up to 2 % in the liquid densities, whereas (excess) free energy differences of up to 5.7 (4.3) kJ mol(-1) are observed. These values are smaller than the differences observed between the three different pairs of methanol/water models: up to 5 % in density and up to 8.5 kJ mol(-1) in (excess) free energy.     

Some issues on the calculation of interfacial properties by molecular simulation

Some of the pitfalls that may befall molecular simulations of interfaces are discussed. They are all related to the calculation of the pressure tensor profiles, which are needed in order to compute surface tensions. We focus on three controversial points: (1) the calculation of the pressure tensor profiles for polyatomic systems, in particular, when the SHAKE algorithm is employed, (2) the addition of long-range corrections to compensate the truncation of the potential, and (3) the importance of including adequate error intervals with the results. Most of the conclusions are general, but some specifically apply to multiple site molecular-dynamics simulations.     

Use of Cartesian Rotation Vectors in Brownian Dynamics Algorithms: Theory and Simulation Results

Summary: We have shown that the components of Cartesian rotation vectors can be used successfully as generalized coordinates describing angular orientation in Brownian dynamics simulations of non‐spherical nanoparticles. For this particular choice of generalized coordinates, we rigorously derived the conformation‐space diffusion equations from kinetic theory for both free nanoparticles and nanoparticles interconnected by springs or holonomic constraints into polymer chains. The equivalent stochastic differential equations were used as a foundation for the Brownian dynamics algorithms. These new algorithms contain singularities only for points in the conformation‐space where both the probability density and its first coordinate derivative equal zero (weak singularities). In addition, the coordinate values after a single Brownian dynamics time step are throughout the conformation‐space equal to the old coordinate values plus the respective increments. For some parts of the conformation‐space these features represent a major improvement compared to the situation when Eulerian angles describe rotational dynamics. The presented simulation results of the equilibrium probability density for free nanoparticles are in perfect agreement with the results from kinetic theory.

Simulation of p^(eq)(Φ) for free nanoparticles.     

Rigid-body dynamics in the isothermal-isobaric ensemble: a test on the accuracy and computational efficiency

We have developed a time-reversible rigid-body (rRB) molecular dynamics algorithm in the isothermal-isobaric (NPT) ensemble. The algorithm is an extension of rigid-body dynamics [Matubayasi and Nakahara, J Chem Phys 1999, 110, 3291] to the NPT ensemble on the basis of non-Hamiltonian statistical mechanics [Martyna, G. J. et al., J Chem Phys 1994, 101, 4177]. A series of MD simulations of water as well as fully hydrated lipid bilayer systems have been undertaken to investigate the accuracy and efficiency of the algorithm. The rRB algorithm was shown to be superior to the state-of-the-art constraint-dynamics algorithm SHAKE/RATTLE/ROLL, with respect to computational efficiency. However, it was revealed that both algorithms produced accurate trajectories of molecules in the NPT as well as NVT ensembles, as long as a reasonably short time step was used. A couple of multiple time-step (MTS) integration schemes were also examined. The advantage of the rRB algorithm for computational efficiency increased when the MD simulation was carried out using MTS on parallel processing computer systems; total computer time for MTS-MD of a lipid bilayer using 64 processors was reduced by about 40% using rRB instead of SHAKE/RATTLE/ROLL.     

Single Alternating Group Explicit（SAGE） Method for Electrochemical Finite Difference Digital Simulation

The four different schemes of Group Explicit Method (GEM):GER,GEL,SAGE and DAGE have been claimed to be unstable when employed for electrochemical digital simulations with large model diffusion coefficient DM.However,in this investigation,in spite of the conditional stability of GER and GEL,the SAGE scheme,which is a combination of GEL and GER,was found to be unconditionally stable when used for simulations of electrochemical reaction-diffusion and had a performance comparable with or even better than the fast Quasi Explicit Finite Difference Method(FQEFD) in some aspects.Corresponding differential equations of SAGE scheme for digital simulations of various electrochemical mechanisms with both uniform and exponentially expanded space units were established.The effectiveness of the SAGE method was further demonstrated by the simulations of an EC and a catalytic mechanism with very large homogeneous rate constants.     

Nucleotide‐Independent Copper(II)‐Based Distance Measurements in DNA by Pulsed ESR Spectroscopy

A site‐specific Cu²⁺ binding motif within a DNA duplex for distance measurements by ESR spectroscopy is reported. This motif utilizes a commercially available 2,2′‐dipicolylamine (DPA) phosphormadite easily incorporated into any DNA oligonucleotide during initial DNA synthesis. The method only requires the simple post‐synthetic addition of Cu²⁺ without the need for further chemical modification. Notably, the label is positioned within the DNA duplex, as opposed to outside the helical perimeter, for an accurate measurement of duplex distance. A distance of 2.7 nm was measured on a doubly Cu²⁺‐labeled DNA sequence, which is in exact agreement with the expected distance from both DNA modeling and molecular dynamic simulations. This result suggests that with this labeling strategy the ESR measured distance directly reports on backbone DNA distance, without the need for further modeling. Furthermore, the labeling strategy is structure‐ and nucleotide‐independent.     

Modeling of Current Distribution on Smooth and Columnar Platinum Structures

Studying the growth and stability of anisotropic or isotropic disordered surfaces in electrodeposition is of importance in catalytic electrochemistry. In some cases, the metallic nature of the electrode defines the topography and roughness, which are also controlled by the experimental time and applied external potential. Because of the experimental restrictions in conventional electrochemical techniques and ex situ electron microscopies, a theoretical model of the surface geometry could aid in understanding the electrodeposition process and current distributions. In spite of applying a complex theory such as dynamic scaling method or perturbation theories, the resolution of mixed mass‐/charge‐transfer equations (tertiary distribution) for the electrodeposition process would give reliable information. One of the main problems with this type of distribution is the mathematics when solving the spatial n‐dimensional differential equations. Use of a primary current distribution is proposed here to simplify the differential equations; however it limits wide application of the first assumption. Distributions of concentration profile, current density, and electrode potential are presented here as a function of the distance normal to the surface for the cases of smooth and rough platinum growth. In the particular case of columnar surfaces, cycloid curves are used to model the electrode, from which the concentration profile is presented in a parameterized form after solving a first‐type curvilinear integral. The concentration contour results in a combination of a trigonometric inverse function and a linear distribution leading to a negative concavity curve. The calculation of the current density and electrode potential contours also show trigonometric shapes exhibiting forbidden imaginary values only at the minimal values of the trochoid curve.     

Multiple‐Well,multiple‐path unimolecular reaction systems. II. 2‐methylhexyl free radicals

Vibrationally excited 2‐methylhexyl radicals formed by shock wave activation or by chemical activation can isomerize by multiple pathways to form any of six stable isomers, can fragment by multiple C? H and C? C bond fission pathways, and can be collisionally stabilized. Master equation simulations of chemical activation and of shock wave activation are used to explore the generic behavior of this complicated coupled system. Selecting the argon pressure in chemical activation systems that produce the 2‐methyl‐1‐hexyl radical isomer ( 1 ) can control the yield of specific isomers. Shock heating of 1 also shows a highly regular sequence of isomer formation. This regular behavior is because the first isomerization steps are faster than subsequent steps. Other radical isomers, such as 2‐methyl‐3‐hexyl ( 3 ), do not show such regular behavior, because the first isomerization step is slower than subsequent steps. Incubation and unimolecular rate‐constant fall‐off are observed in the shock wave simulations. The unimolecular rate‐constant fall‐off for the coupled system produces low‐pressure limiting rate constants proportional to [M]ⁿ, where n can be greater than unity. The fact that n can be greater than unity is a natural feature of multichannel coupled unimolecular reaction systems, but detection of the effect in experiments may be very demanding. © 2001 John Wiley & Sons, Inc. Int J Chem Kinet 33: 246–261, 2001     

A 45-ns molecular dynamics simulation of hemoglobin in water by vectorizing and parallelizing COSMOS90 on the earth simulator: dynamics of tertiary and quaternary structures

Molecular dynamics (MD) simulations of human adult hemoglobin (HbA) were carried out for 45 ns in water with all degrees of freedom including bond stretching and without any artificial constraints. To perform such large-scale simulations, one of the authors (M.S.) accelerated his own software COSMOS90 on the Earth Simulator by vectorization and parallelization. The dynamical features of HbA were investigated by evaluating root-mean-square deviations from the initial X-ray structure (an oxy T-state hemoglobin with PDB code: 1GZX) and root-mean-square fluctuations around the average structure from the simulation trajectories. The four subunits (alpha(1), alpha(2), beta(1), and beta(2)) of HbA maintained structures close to their respective X-ray structures during the simulations even though no constraints were applied to HbA in the simulations. Dimers alpha(1)beta(1) and alpha(2)beta(2) also maintained structures close to their respective X-ray structures while they moved relative to each other like two stacks of dumbbells. The distance between the two dimers (alpha(1)beta(1) and alpha(2)beta(2)) increased by 2 A (7.4%) in the initial 15 ns and stably fluctuated at the distance with the standard deviation 0.2 A. The relative orientation of the two dimers fluctuated between the initial X-ray angle -100 degrees and about -105 degrees with intervals of a few tens of nanoseconds.