This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matri... more This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matrix to square-reduced form and the approximation of all its eigenvalues using the implicit version of Van Loan's method. The transformation of the Hamiltonian matrix to a square-reduced form transforms a Hamiltonian eigenvalue problem of order 2n to a Hessenberg eigenvalue problem of order n. The eigenvalues of the Hamiltonian matrix are the square roots of those of the Hessenberg matrix. Symplectic scaling and norm scaling are provided, which, in some cases, improve the accuracy of the computed eigenvalues. We demonstrate the performance of the subroutines for several examples and show how they can be used to solve some control-theoretic problems.
This article describes the design and implementation of a variable timestep method for simulating... more This article describes the design and implementation of a variable timestep method for simulating time-reversible constrained dynamical systems. Based on the Adaptive Verlet method of Huang and Leimkuhler 14], and the SHAKE 23] and RATTLE 1] discretizations, the new method (VRATTLE) deenes a mapping of the constraint manifold which preserves the reversible structure. It achieves this through the solution of a single additional scalar nonlinear equation at each timestep, together with the equations of constraint. As a nontrivial application, we simulate the dynamics of an elastic (inextensible, unshearable) rod undergoing large deformations and collisions with the sides of a bounding box. Numerical experiments indicate that adapting the stepsize using VRATTLE can smooth the numerical energy and improve the overall eeciency of the simulation.
ABSTRACT This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamilton... more ABSTRACT This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matrix to square-reduced form and the approximation of all its eigenvalues using the implicit version of Van Loan's method. The transformation of the Hamiltonian matrix to a square-reduced form transforms a Hamiltonian eigenvalue problem of order 2n to a Hessenberg eigenvalue problem of order n. The eigenvalues of the Hamiltonian matrix are the square roots of those of the Hessenberg matrix. Symplectic scaling and norm scaling are provided, which, in some cases, improve the accuracy of the computed eigenvalues. We demonstrate the performance of the subroutines for several examples and show how they can be used to solve some control-theoretic problems. Key words: eigenvalues, (square-reduced) Hamiltonian matrix, skew-Hamiltonian matrix, algebraic Riccati equation AMS(MOS) subject classifications: 65F15, 93B40 Some of this work was complete while this author was with the Fakultat fur Mathemati...
Siam Journal on Scientific Computing, Nov 16, 1997
. This article describes the design and implementation of a variable timestep method for simulati... more . This article describes the design and implementation of a variable timestep method for simulating time-reversible constrained dynamical systems. Based on the Adaptive Verlet method of Huang and Leimkuhler [13], and the SHAKE [20] and RATTLE [1] discretizations, the new method ...
Page 1. CHAPTER 7 Extending the Time Scale in Atomically Detailed Simulations Alfredo E. Ca´rdena... more Page 1. CHAPTER 7 Extending the Time Scale in Atomically Detailed Simulations Alfredo E. Ca´rdenasa and Eric Barthb aDepartment of Chemistry, University of South Florida, Tampa, Florida bDepartment of Mathematics and ...
This article describes Fortran 77 subroutines for computing eigenvalues and invariant subspaces o... more This article describes Fortran 77 subroutines for computing eigenvalues and invariant subspaces of Hamiltonian and skew-Hamiltonian matrices. The implemented algorithms are based on orthogonal symplectic decompositions, implying numerical backward stability as well as symmetry preservation for the computed eigenvalues. These algorithms are supplemented with balancing and block algorithms, which can lead to considerable accuracy and performance improvements. As a by-product, an efficient implementation for computing symplectic QR decompositions is provided. We demonstrate the usefulness of the subroutines for several, practically relevant examples.
Besides preserving the energy, the ow of a conservative multibody system possesses important geom... more Besides preserving the energy, the ow of a conservative multibody system possesses important geometric (symplectic) invariants. Symplectic discretization schemes that mimic the corresponding feature of the true ow have been shown to be e ective alternatives to standard methods for many conservative problems. For systems of rigid bodies, the development of such schemes can be complicated or costly to implement, depending on the choice of problem formulation. In this article, we demonstrate that a special formulation of the multibody system (based on a particle representation) together with a symplectic discretization for constrained problems borrowed from molecular dynamics o ers an e cient alternative to standard approaches. Numerical experiments illustrating this approach are described. 1991 Mathematics Subject Classi cation. 65L05, 70-08, 70F20. Key words and phrases. symplectic discretizations, constrained systems, di erential-algebraic equations, multibody systems, molecular dynamics, modelling of mechanisms, modelling of manipulators.
Innovative algorithms have been developed during the past decade for simulating Newtonian physics... more Innovative algorithms have been developed during the past decade for simulating Newtonian physics for macromolecules. A major goal is alleviation of the severe requirement that the integration timestep be small enough to resolve the fastest components of the motion and thus guarantee numerical stability. This timestep problem is challenging if strictly faster methods with the same all-atom resolution at small timesteps are sought. Mathematical techniques that have worked well in other multiple-timescale contexts-where the fast motions are rapidly decaying or largely decoupled from others-have not been as successful for biomolecules, where vibrational coupling is strong.
In molecular dynamics simulations, the fastest components of the potential eld impose severe rest... more In molecular dynamics simulations, the fastest components of the potential eld 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 paper, the resulting systems of mixed di erential 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 modi cation based on Successive OverRelaxation (SOR). A simple adaptive algorithm for nding the optimal relaxation parameter is presented. Alternative direct methods using sparse matrix techniques are discussed. Numerical results are given for the new techniques, implemented in the molecular modeling software package CHARMM, showing as much as twofold improvement over SHAKE iteration. matrix methods 1. Introduction. In molecular dynamics, the length of timestep for numerically integrating the equations of motion is dictated by the contributions to the force vector which maintain pairs of atoms near some equilibrium distance. The imposition of algebraic constraints that x these lengths removes the associated rapid vibrational modes, enabling the use of longer timesteps without substantially altering important physical characteristics of the motion 1]. Although we treat only length constraints in the present work, constrained techniques are also of interest for conformational search and conformational free energy simulations 2]. In 3] the SHAKE iteration was described for solving the nonlinear equations at each timestep of a constrained version of the Verlet discretization, and a similar scheme was proposed in 4] for use with the RATTLE discretization.
We present an efficient new method termed LN for propagating biomolecular dynamics according to t... more We present an efficient new method termed LN for propagating biomolecular dynamics according to the Langevin equation that arose fortuitously upon analysis of the range of harmonic validity of our normal-mode scheme LIN. LN combines force linearization with force splitting techniques and disposes of LIN's computationally intensive minimization ͑anharmonic correction͒ component. Unlike the competitive multiple-timestepping ͑MTS͒ schemes today-formulated to be symplectic and time-reversible-LN merges the slow and fast forces via extrapolation rather than ''impulses;'' the Langevin heat bath prevents systematic energy drifts. This combination succeeds in achieving more significant speedups than these MTS methods which are limited by resonance artifacts to an outer timestep less than some integer multiple of half the period of the fastest motion ͑around 4-5 fs for biomolecules͒. We show that LN achieves very good agreement with small-timestep solutions of the Langevin equation in terms of thermodynamics ͑energy means and variances͒, geometry, and dynamics ͑spectral densities͒ for two proteins in vacuum and a large water system. Significantly, the frequency of updating the slow forces extends to 48 fs or more, resulting in speedup factors exceeding 10. The implementation of LN in any program that employs force-splitting computations is straightforward, with only partial second-derivative information required, as well as sparse Hessian/vector multiplication routines. The linearization part of LN could even be replaced by direct evaluation of the fast components. The application of LN to biomolecular dynamics is well suited for configurational sampling, thermodynamic, and structural questions.
Force splitting or multiple timestep ͑MTS͒ methods are effective techniques that accelerate biomo... more Force splitting or multiple timestep ͑MTS͒ methods are effective techniques that accelerate biomolecular dynamics simulations by updating the fast and slow forces at different frequencies.
We present a general molecular-dynamics simulation scheme, based on the Nosé thermostat, for samp... more We present a general molecular-dynamics simulation scheme, based on the Nosé thermostat, for sampling according to arbitrary phase space distributions. We formulate numerical methods based on both Nosé-Hoover and Nosé-Poincaré thermostats for two specific classes of distributions; namely, those that are functions of the system Hamiltonian and those for which position and momentum are statistically independent. As an example, we propose a generalized variable temperature distribution that is designed to accelerate sampling in molecular systems.
This article considers the design and implementation of variable-timestep methods for simulating ... more This article considers the design and implementation of variable-timestep methods for simulating holonomically constrained mechanical systems. Symplectic variable stepsizes are brie y discussed, we then consider time-reparameterization techniques employing a time-reversible (symmetric) integration method to solve the equations of motion. We give several numerical examples, including a simulation of an elastic (inextensible, unshearable) rod undergoing large deformations and collisions with the sides of a bounding box. Numerical experiments indicate that adaptive stepping can signi cantly smooth the numerical energy and improve the overall e ciency of the simulation.
This article describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian mat... more This article describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matrix to square-reduced form and the approximation of all its eigenvalues using the implicit version of Van Loan's method. The transformation of the Hamiltonian matrix to a square-reduced form transforms a Hamiltonian eigenvalue problem of order 2n to a Hessenberg eigenvalue problem of order n. The eigenvalues of the Hamiltonian matrix are the square roots of those of the Hessenberg matrix. Symplectic scaling and norm scaling are provided, which, in some cases, improve the accuracy of the computed eigenvalues. We demonstrate the performance of the subroutines for several examples and show how they can be used to solve some control-theoretic problems.
This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matri... more This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matrix to square-reduced form and the approximation of all its eigenvalues using the implicit version of Van Loan's method. The transformation of the Hamiltonian matrix to a square-...
This article describes the design and implementation of a variable timestep method for simulating... more This article describes the design and implementation of a variable timestep method for simulating time-reversible constrained dynamical systems. Based on the Adaptive Verlet method of Huang and Leimkuhler 13], and the SHAKE 20] and RATTLE 1] discretizations, the new method VRATTLE de nes a mapping of the constraint manifold which preserves the reversible structure. It achieves this through the solution of a single scalar nonlinear equation at each timestep, together with the equations of constraint. As a nontrivial application, we simulate the dynamics of an elastic (inextensible, unshearable) rod undergoing large deformations and collisions with the sides of a bounding box. Numerical experiments indicate that adapting the stepsize using VRATTLE can smooth the numerical energy and improve the overall e ciency of the simulation.
This article describes a collection of model problems for aiding numerical analysts, code develop... more This article describes a collection of model problems for aiding numerical analysts, code developers and others in the design of computational methods for molecular dynamics (MD) simulation. Common types of calculations and desirable features of algorithms are surveyed, and these are used to guide selection of representative models. By including essential features of certain classes of molecular systems, but otherwise limiting the physical and quantitative details, it is hoped that the test set can help to facilitate cross-disciplinary algorithm and code development e orts.
This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matri... more This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matrix to square-reduced form and the approximation of all its eigenvalues using the implicit version of Van Loan's method. The transformation of the Hamiltonian matrix to a square-reduced form transforms a Hamiltonian eigenvalue problem of order 2n to a Hessenberg eigenvalue problem of order n. The eigenvalues of the Hamiltonian matrix are the square roots of those of the Hessenberg matrix. Symplectic scaling and norm scaling are provided, which, in some cases, improve the accuracy of the computed eigenvalues. We demonstrate the performance of the subroutines for several examples and show how they can be used to solve some control-theoretic problems.
This article describes the design and implementation of a variable timestep method for simulating... more This article describes the design and implementation of a variable timestep method for simulating time-reversible constrained dynamical systems. Based on the Adaptive Verlet method of Huang and Leimkuhler 14], and the SHAKE 23] and RATTLE 1] discretizations, the new method (VRATTLE) deenes a mapping of the constraint manifold which preserves the reversible structure. It achieves this through the solution of a single additional scalar nonlinear equation at each timestep, together with the equations of constraint. As a nontrivial application, we simulate the dynamics of an elastic (inextensible, unshearable) rod undergoing large deformations and collisions with the sides of a bounding box. Numerical experiments indicate that adapting the stepsize using VRATTLE can smooth the numerical energy and improve the overall eeciency of the simulation.
ABSTRACT This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamilton... more ABSTRACT This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matrix to square-reduced form and the approximation of all its eigenvalues using the implicit version of Van Loan's method. The transformation of the Hamiltonian matrix to a square-reduced form transforms a Hamiltonian eigenvalue problem of order 2n to a Hessenberg eigenvalue problem of order n. The eigenvalues of the Hamiltonian matrix are the square roots of those of the Hessenberg matrix. Symplectic scaling and norm scaling are provided, which, in some cases, improve the accuracy of the computed eigenvalues. We demonstrate the performance of the subroutines for several examples and show how they can be used to solve some control-theoretic problems. Key words: eigenvalues, (square-reduced) Hamiltonian matrix, skew-Hamiltonian matrix, algebraic Riccati equation AMS(MOS) subject classifications: 65F15, 93B40 Some of this work was complete while this author was with the Fakultat fur Mathemati...
Siam Journal on Scientific Computing, Nov 16, 1997
. This article describes the design and implementation of a variable timestep method for simulati... more . This article describes the design and implementation of a variable timestep method for simulating time-reversible constrained dynamical systems. Based on the Adaptive Verlet method of Huang and Leimkuhler [13], and the SHAKE [20] and RATTLE [1] discretizations, the new method ...
Page 1. CHAPTER 7 Extending the Time Scale in Atomically Detailed Simulations Alfredo E. Ca´rdena... more Page 1. CHAPTER 7 Extending the Time Scale in Atomically Detailed Simulations Alfredo E. Ca´rdenasa and Eric Barthb aDepartment of Chemistry, University of South Florida, Tampa, Florida bDepartment of Mathematics and ...
This article describes Fortran 77 subroutines for computing eigenvalues and invariant subspaces o... more This article describes Fortran 77 subroutines for computing eigenvalues and invariant subspaces of Hamiltonian and skew-Hamiltonian matrices. The implemented algorithms are based on orthogonal symplectic decompositions, implying numerical backward stability as well as symmetry preservation for the computed eigenvalues. These algorithms are supplemented with balancing and block algorithms, which can lead to considerable accuracy and performance improvements. As a by-product, an efficient implementation for computing symplectic QR decompositions is provided. We demonstrate the usefulness of the subroutines for several, practically relevant examples.
Besides preserving the energy, the ow of a conservative multibody system possesses important geom... more Besides preserving the energy, the ow of a conservative multibody system possesses important geometric (symplectic) invariants. Symplectic discretization schemes that mimic the corresponding feature of the true ow have been shown to be e ective alternatives to standard methods for many conservative problems. For systems of rigid bodies, the development of such schemes can be complicated or costly to implement, depending on the choice of problem formulation. In this article, we demonstrate that a special formulation of the multibody system (based on a particle representation) together with a symplectic discretization for constrained problems borrowed from molecular dynamics o ers an e cient alternative to standard approaches. Numerical experiments illustrating this approach are described. 1991 Mathematics Subject Classi cation. 65L05, 70-08, 70F20. Key words and phrases. symplectic discretizations, constrained systems, di erential-algebraic equations, multibody systems, molecular dynamics, modelling of mechanisms, modelling of manipulators.
Innovative algorithms have been developed during the past decade for simulating Newtonian physics... more Innovative algorithms have been developed during the past decade for simulating Newtonian physics for macromolecules. A major goal is alleviation of the severe requirement that the integration timestep be small enough to resolve the fastest components of the motion and thus guarantee numerical stability. This timestep problem is challenging if strictly faster methods with the same all-atom resolution at small timesteps are sought. Mathematical techniques that have worked well in other multiple-timescale contexts-where the fast motions are rapidly decaying or largely decoupled from others-have not been as successful for biomolecules, where vibrational coupling is strong.
In molecular dynamics simulations, the fastest components of the potential eld impose severe rest... more In molecular dynamics simulations, the fastest components of the potential eld 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 paper, the resulting systems of mixed di erential 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 modi cation based on Successive OverRelaxation (SOR). A simple adaptive algorithm for nding the optimal relaxation parameter is presented. Alternative direct methods using sparse matrix techniques are discussed. Numerical results are given for the new techniques, implemented in the molecular modeling software package CHARMM, showing as much as twofold improvement over SHAKE iteration. matrix methods 1. Introduction. In molecular dynamics, the length of timestep for numerically integrating the equations of motion is dictated by the contributions to the force vector which maintain pairs of atoms near some equilibrium distance. The imposition of algebraic constraints that x these lengths removes the associated rapid vibrational modes, enabling the use of longer timesteps without substantially altering important physical characteristics of the motion 1]. Although we treat only length constraints in the present work, constrained techniques are also of interest for conformational search and conformational free energy simulations 2]. In 3] the SHAKE iteration was described for solving the nonlinear equations at each timestep of a constrained version of the Verlet discretization, and a similar scheme was proposed in 4] for use with the RATTLE discretization.
We present an efficient new method termed LN for propagating biomolecular dynamics according to t... more We present an efficient new method termed LN for propagating biomolecular dynamics according to the Langevin equation that arose fortuitously upon analysis of the range of harmonic validity of our normal-mode scheme LIN. LN combines force linearization with force splitting techniques and disposes of LIN's computationally intensive minimization ͑anharmonic correction͒ component. Unlike the competitive multiple-timestepping ͑MTS͒ schemes today-formulated to be symplectic and time-reversible-LN merges the slow and fast forces via extrapolation rather than ''impulses;'' the Langevin heat bath prevents systematic energy drifts. This combination succeeds in achieving more significant speedups than these MTS methods which are limited by resonance artifacts to an outer timestep less than some integer multiple of half the period of the fastest motion ͑around 4-5 fs for biomolecules͒. We show that LN achieves very good agreement with small-timestep solutions of the Langevin equation in terms of thermodynamics ͑energy means and variances͒, geometry, and dynamics ͑spectral densities͒ for two proteins in vacuum and a large water system. Significantly, the frequency of updating the slow forces extends to 48 fs or more, resulting in speedup factors exceeding 10. The implementation of LN in any program that employs force-splitting computations is straightforward, with only partial second-derivative information required, as well as sparse Hessian/vector multiplication routines. The linearization part of LN could even be replaced by direct evaluation of the fast components. The application of LN to biomolecular dynamics is well suited for configurational sampling, thermodynamic, and structural questions.
Force splitting or multiple timestep ͑MTS͒ methods are effective techniques that accelerate biomo... more Force splitting or multiple timestep ͑MTS͒ methods are effective techniques that accelerate biomolecular dynamics simulations by updating the fast and slow forces at different frequencies.
We present a general molecular-dynamics simulation scheme, based on the Nosé thermostat, for samp... more We present a general molecular-dynamics simulation scheme, based on the Nosé thermostat, for sampling according to arbitrary phase space distributions. We formulate numerical methods based on both Nosé-Hoover and Nosé-Poincaré thermostats for two specific classes of distributions; namely, those that are functions of the system Hamiltonian and those for which position and momentum are statistically independent. As an example, we propose a generalized variable temperature distribution that is designed to accelerate sampling in molecular systems.
This article considers the design and implementation of variable-timestep methods for simulating ... more This article considers the design and implementation of variable-timestep methods for simulating holonomically constrained mechanical systems. Symplectic variable stepsizes are brie y discussed, we then consider time-reparameterization techniques employing a time-reversible (symmetric) integration method to solve the equations of motion. We give several numerical examples, including a simulation of an elastic (inextensible, unshearable) rod undergoing large deformations and collisions with the sides of a bounding box. Numerical experiments indicate that adaptive stepping can signi cantly smooth the numerical energy and improve the overall e ciency of the simulation.
This article describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian mat... more This article describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matrix to square-reduced form and the approximation of all its eigenvalues using the implicit version of Van Loan's method. The transformation of the Hamiltonian matrix to a square-reduced form transforms a Hamiltonian eigenvalue problem of order 2n to a Hessenberg eigenvalue problem of order n. The eigenvalues of the Hamiltonian matrix are the square roots of those of the Hessenberg matrix. Symplectic scaling and norm scaling are provided, which, in some cases, improve the accuracy of the computed eigenvalues. We demonstrate the performance of the subroutines for several examples and show how they can be used to solve some control-theoretic problems.
This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matri... more This paper describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matrix to square-reduced form and the approximation of all its eigenvalues using the implicit version of Van Loan's method. The transformation of the Hamiltonian matrix to a square-...
This article describes the design and implementation of a variable timestep method for simulating... more This article describes the design and implementation of a variable timestep method for simulating time-reversible constrained dynamical systems. Based on the Adaptive Verlet method of Huang and Leimkuhler 13], and the SHAKE 20] and RATTLE 1] discretizations, the new method VRATTLE de nes a mapping of the constraint manifold which preserves the reversible structure. It achieves this through the solution of a single scalar nonlinear equation at each timestep, together with the equations of constraint. As a nontrivial application, we simulate the dynamics of an elastic (inextensible, unshearable) rod undergoing large deformations and collisions with the sides of a bounding box. Numerical experiments indicate that adapting the stepsize using VRATTLE can smooth the numerical energy and improve the overall e ciency of the simulation.
This article describes a collection of model problems for aiding numerical analysts, code develop... more This article describes a collection of model problems for aiding numerical analysts, code developers and others in the design of computational methods for molecular dynamics (MD) simulation. Common types of calculations and desirable features of algorithms are surveyed, and these are used to guide selection of representative models. By including essential features of certain classes of molecular systems, but otherwise limiting the physical and quantitative details, it is hoped that the test set can help to facilitate cross-disciplinary algorithm and code development e orts.
Uploads
Papers by Eric Barth