Exact and efficient calculation of Lagrange multipliers in constrained biological polymers: Proteins and nucleic acids as example cases

In order to accelerate molecular dynamics simulations, it is very common to impose holonomic constraints on their hardest degrees of freedom. In this way, the time step used to integrate the equations of motion can be increased, thus allowing, in principle, to reach longer total simulation times. The imposition of such constraints results in an aditional set of Nc equations (the equations of constraint) and unknowns (their associated Lagrange multipliers), that must be solved in one way or another at each time step of the dynamics. In this work it is shown that, due to the essentially linear structure of typical biological polymers, such as nucleic acids or proteins, the algebraic equations that need to be solved involve a matrix which is banded if the constraints are indexed in a clever way. This allows to obtain the Lagrange multipliers through a non-iterative procedure, which can be considered exact up to machine precision, and which takes O(Nc) operations, instead of the usual O(Nc3) for generic molecular systems. We develop the formalism, and describe the appropriate indexing for a number of model molecules and also for alkanes, proteins and DNA. Finally, we provide a numerical example of the technique in a series of polyalanine peptides of different lengths using the AMBER molecular dynamics package.Comment: 29 pages, 10 figures, 1 tabl

The canonical equilibrium of constrained molecular models

In order to increase the efficiency of the computer simulation of biological molecules, it is very common to impose holonomic constraints on the fastest degrees of freedom; normally bond lengths, but also possibly bond angles. Since the maximum time step required for the stability of the dynamics is proportional to the shortest period associated with the motions of the system, constraining the fastest vibrations allows to increase it and, assuming that the added numerical cost is not too high, also increase the overall efficiency of the simulation. However, as any other element that affects the physical model, the imposition of constraints must be assessed from the point of view of accuracy: both the dynamics and the equilibrium statistical mechanics are model-dependent, and they will be changed if constraints are used. In this review, we investigate the accuracy of constrained models at the level of the equilibrium statistical mechanics distributions produced by the different dynamics. We carefully derive the canonical equilibrium distributions of both the constrained and unconstrained dynamics, comparing the two of them by means of a “stiff” approximation to the latter. We do so both in the case of flexible and hard constraints, i.e., when the value of the constrained coordinates depends on the conformation and when it is a constant number. We obtain the different correcting terms associated with the kinetic energy mass-metric tensor determinants, but also with the details of the potential energy in the vicinity of the constrained subspace (encoded in its first and second derivatives). This allows us to directly compare, at the conformational level, how the imposition of constraints changes the thermal equilibrium of molecular systems with respect to the unconstrained case. We also provide an extensive review of the relevant literature, and we show that all models previously reported can be considered special cases of the most general treatments presented in this work. Finally, we numerically analyze a simple methanol molecule in order to illustrate the theoretical concepts in a practical case.This work has been supported by the research projects E24/3 (DGA, Spain), FIS2009- 13364-C02-01 (MICINN, Spain), 200980I064 (CSIC, Spain) and ARAID and Ibercaja grant for young researchers (Spain). P. G.-R. is supported by a JAE PREDOC grant (CSIC, Spain)

How accurate does Newton have to be?

We analyze the convergence of quasi-Newton methods in exact and finite precision arithmetic. In particular, we derive an upper bound for the stagnation level and we show that any sufficiently exact quasi-Newton method will converge quadratically until stagnation. In the absence of sufficient accuracy, we are likely to retain rapid linear convergence. We confirm our analysis by computing square roots and solving bond constraint equations in the context of molecular dynamics. We briefly discuss implications for parallel solvers.Comment: 12 pages, 2 figures, preprint accepted by PPAM 2022, expected to appear in LNCS vol. 13826 during 202

Implementación e integración en GROMACS de un algoritmo eficiente y preciso para imponer ligaduras en simulaciones de dinámica molecular

Las simulaciones de dinámica molecular describen la evolución en el tiempo de un sistema de partículas. Herramientas computacionales para realizar dichas simulaciones son fundamentales para los investigadores de campos tan diversos como la medicina (búsqueda de tratamientos para enfermedades como el Alzheimer, la fibrosis quística o el cáncer; diseño computacional de medicamentos), la química (diseño de catalizadores) o la ingeniería de materiales. GROMACS es un extendido y versátil paquete de dinámica molecular escrito en el lenguaje de programación C. Originalmente desarrollado en la Universidad de Groningen (Países Bajos), actualmente se mantiene por desarrolladores de universidades y centros de investigación de todo el mundo, y cuenta con miles de usuarios. Para realizar simulaciones eficientes, una opción común es imponer ligaduras, es decir, restricciones a la longitud de los enlaces atómicos o al valor de los distintos ángulos entre átomos. Esto permite incrementar el paso temporal (time step) de las simulaciones y con ello lograr realizar simulaciones más duraderas, con lo que se incrementa el poder predictivo y el realismo de la simulación. Los dos algoritmos más utilizados para la implementación de ligaduras, los denominados SHAKE y LINCS (LINear Constraint Solver), presentan limitaciones a la hora de imponer ligaduras en los ángulos. Además, están basados en aproximaciones, lo cual afecta negativamente a su exactitud, eficiencia y estabilidad numérica. El presente proyecto fin de carrera se basa en la reciente propuesta teórica del Dr. Pablo García Risueño, la cual propone mejorar la imposición de ligaduras en las simulaciones de dinámica molecular mediante el uso de cálculos analíticos. Por lo tanto, el objetivo de este proyecto es implementar el algoritmo denominado ILVES-S, propuesta alternativa al algoritmo SHAKE, e integrarlo en el paquete de dinámica molecular GROMACS

A review of High Performance Computing foundations for scientists

The increase of existing computational capabilities has made simulation emerge as a third discipline of Science, lying midway between experimental and purely theoretical branches [1, 2]. Simulation enables the evaluation of quantities which otherwise would not be accessible, helps to improve experiments and provides new insights on systems which are analysed [3-6]. Knowing the fundamentals of computation can be very useful for scientists, for it can help them to improve the performance of their theoretical models and simulations. This review includes some technical essentials that can be useful to this end, and it is devised as a complement for researchers whose education is focused on scientific issues and not on technological respects. In this document we attempt to discuss the fundamentals of High Performance Computing (HPC) [7] in a way which is easy to understand without much previous background. We sketch the way standard computers and supercomputers work, as well as discuss distributed computing and discuss essential aspects to take into account when running scientific calculations in computers.Comment: 33 page

Linearly scaling direct method for accurately inverting sparse banded matrices

In many problems in Computational Physics and Chemistry, one finds a special kind of sparse matrices, termed "banded matrices". These matrices, which are defined as having non-zero entries only within a given distance from the main diagonal, need often to be inverted in order to solve the associated linear system of equations. In this work, we introduce a new O(n) algorithm for solving such a system, being n X n the size of the matrix. We produce the analytical recursive expressions that allow to directly obtain the solution, as well as the pseudocode for its computer implementation. Moreover, we review the different options for possibly parallelizing the method, we describe the extension to deal with matrices that are banded plus a small number of non-zero entries outside the band, and we use the same ideas to produce a method for obtaining the full inverse matrix. Finally, we show that the New Algorithm is competitive, both in accuracy and in numerical efficiency, when compared to a standard method based in Gaussian elimination. We do this using sets of large random banded matrices, as well as the ones that appear when one tries to solve the 1D Poisson equation by finite differences.Comment: 24 pages, 5 figures, submitted to J. Comp. Phy

Improving simulation of biological molecules: refining mathematical, physical and computational tools

[ES]: El objetivo de las investigaciones correspondientes a este documento de tesis doctoral es la mejora de los métodos para la simulación de moléculas, especialmente moléculas biológicas, como las proteínas. Las capacidades de los ordenadores han aumentado sobremanera en las últimas décadas, y actualmente son una herramienta poderosa para el cálculo de magnitudes físicas y químicas. Debido a esto, un gran número de investigadores se dedica al desarrollo de métodos computacionales precisos y eficientes, a su uso en simulaciones con fines científicos, a una mezcla de ambos o a otros campos dentro de la ciencia computacional. El trabajo de investigación resumido en esta tesis tiene como finalidad la mejora de la eficiencia y exactitud de los métodos de simulación de moléculas biológicas.[EN]: The goal of the research summarized in this Ph. D. dissertation is the improvement of the computational methods for simulation of molecules, specially biological molecules such as proteins. The performance of computers has strongly increased during the last decades, and at present they are a powerful tool for the calculation of physical and chemical quantities. Because of this, there exist a large number of researchers who work either as developers of reliable and efficient computational methods to make calculations on physical and chemical systems, or as users of existing algorithms and programs for scientific purposes. The work presented in this dissertation is expected to improve the efficiency and accuracy of simulations of biological molecules.Peer reviewe