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

Pathogenesis of Staphylococcus epidermidis in prosthetic joint infections : can identification of virulence genes differentiate between infecting and commensal strains?

Background: Staphylococcus epidermidis is a commensal of human skin flora and a frequent causative micro-organism in prosthetic joint infections (PJIs). To date, no single marker has been identified to distinguish infecting strains from commensal S. epidermidis populations. Aim: To find possible genetic markers to distinguish between the two populations. Methods: Fifty S. epidermidis strains from patients with PJIs were analysed, 50 from the skin of healthy individuals (commensal strains) and 17 from the surgical field of patients undergoing primary arthroplasty. In these three groups the antimicrobial susceptibility profile, sequence type, biofilm formation, and virulence factors were studied. Strains from the surgical field have not been compared previously with strains from the other two groups. Findings: S. epidermidis strains from PJI patients were significantly more antibiotic resistant than commensal strains and surgical field strains. A wide variety of sequence types was found in commensal and surgical field strains. The predominant sequence type was ST2 and it was only present in PJI strains (44%). Differences in biofilm production did not differ between populations. Virulence genes sdrF and bhp, the complete ica operon, and the insertion sequence IS256 were significantly predominant in PJI strains. By contrast, embp and hld genes and the arginine catabolic mobile element (ACME) were more prevalent in commensal strains. Surgical field strains could be a valid control group to discriminate between infecting and commensal strains. Conclusion: A combination of characteristic features can differentiate between infecting and commensal S. epidermidis strains in PJI, whereas a single marker cannot

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