Bibliographic Review on Molecular Dynamics: Modeling, Algorithms, and Applications.

Escuela Técnica Superior de Ingenierías Informática y de Telecomunicación y Facultad de Ciencias Grado en Ingeniería Informática y Matemáticas trabajo de fin de grado Revisión bibliográfica sobre dinámica molecular. Modelización, algoritmos y aplicaciones Presentado por: Joaquín Arcila Pérez Curso académico 2024-2025 Revisión bibliográfica sobre dinámica molecular. Modelización, algoritmos y aplicaciones Joaquín Arcila Pérez Joaquín Arcila Pérez Revisión bibliográfica sobre dinámica molecular. Modelización, algoritmos y aplicaciones. Trabajo de fin de Grado. Curso académico 2024-2025. Responsable de tutorización Juan Calvo Yagüe Departamento de Matemática Aplicada Lázaro René Izquierdo Fábregas Departamento de Matemática Aplicada Grado en Ingeniería Informática y Matemáticas Escuela Técnica Superior de Ingenierías Informática y de Telecomunicación y Facultad de Ciencias Universidad de Granada Declaración de originalidad D. Joaquín Arcila Pérez Declaro explícitamente que el trabajo presentado como Trabajo de Fin de Grado (TFG), correspondiente al curso académico 2024-2025, es original, entendido esto en el sentido de que no he utilizado para la elaboración del trabajo fuentes sin citarlas debidamente. En Granada a 14 de julio de 2025 Fdo: Joaquín Arcila Pérez A mis padres, por su esfuerzo, dedicación y la confianza que siempre han depositado en mí. A mis hermanos, por ser un ejemplo a seguir y por haberme guiado a lo largo del camino. A mi familia y a mis amigos, por estar siempre ahí y por animarme cada día a alcanzar mis metas. A todos mis profesores, especialmente a mis tutores, Juan y René, y a mi tío Alberto, por su paciencia y el interés mostrado en mi proceso de aprendizaje. Índice general Summary VII Resumen IX Introducción XI Planificación y Presupuesto XIII 1. Contexto histórico 1 2. Fundamentos de Dinámica Molecular 3 2.1. Principios de la mecánica clásica ........................... 3 2.1.1. Leyes de Newton ................................ 3 2.1.2. Principio de conservación de la energía ................... 4 2.1.3. Introducción a las ecuaciones de movimiento de Newton ........ 5 3. Modelización en Dinámica Molecular 7 3.1. Modelos de interacción ................................. 7 3.1.1. Campos de fuerza ............................... 7 3.2. Representación de sistemas moleculares ....................... 11 3.2.1. Modelos atomísticos vs coarse-graining ................... 11 3.2.2. Cajas de simulación .............................. 12 3.2.3. Condiciones de contorno ........................... 12 3.2.4. Métodos de truncamiento ........................... 13 4. Algoritmos en Dinámica Molecular 15 4.1. Integración de ecuaciones de movimiento ...................... 15 4.1.1. Algoritmo de Verlet .............................. 19 4.1.2. Velocity-Verlet .................................. 20 4.1.3. Leap-Frog .................................... 21 4.2. Control de temperatura y presión .......................... 23 4.2.1. Termostatos ................................... 24 4.2.2. Barostatos .................................... 27 4.3. Algoritmos de optimización .............................. 29 4.3.1. Minimización de energía ........................... 30 4.3.2. Métodos de Monte Carlo ........................... 33 4.4. Simulaciones paralelas y aceleración ......................... 35 5. Implementación de una Simulación de Dinámica Molecular 37 5.1. Inicialización del sistema ................................ 39 5.2. Cálculo de fuerzas ................................... 41 5.3. Integración de las ecuaciones de movimiento .................... 43 5.4. Programas de dinámica molecular .......................... 44 v Introducción En el Capítulo 2se introducen los fundamentos físicos que sustentan la dinámica molecular, analizándose las leyes y ecuaciones de Newton, así como los conceptos clave relacionados con la energía, la temperatura y las condiciones iniciales del sistema. El Capítulo 3trata sobre la modelización en dinámica molecular. Se profundiza en los distintos potenciales de interacción y condiciones de contorno, y se evidencia la importancia de la representación del sistema en MD. El Capítulo 4está dedicado al estudio y análisis de los principales algoritmos empleados en dinámica molecular. Se presentan los métodos numéricos utilizados para integrar las ecuaciones de movimiento, con especial atención al algoritmo de Verlet y sus variantes. Asimismo, se abordan técnicas para el control de las condiciones termodinámicas del sistema mediante termostatos (como los de Berendsen y Nosé–Hoover) y barostatos (como Andersen y Parrinello–Rahman). También se exploran algoritmos de optimización, cuyo objetivo es hallar configuraciones de mínima energía. Finalmente, se discuten estrategias para mejorar el rendimiento computacional mediante el uso de simulaciones paralelas y técnicas de aceleración por hardware, como el empleo de unidades de procesamiento gráfico (GPU) y arquitecturas de alto rendimiento (HPC). En el Capítulo 5se presenta una posible implementación de un programa de dinámica molecular, exponiendo un pseudo-algoritmo. Este capítulo permite ilustrar de forma concreta los aspectos técnicos tratados en los anteriores. En el Capítulo 6se exploran algunas de las principales aplicaciones de la dinámica molecular, con especial atención a su uso en biología computacional, simulación de proteínas y diseño de fármacos. Se destaca cómo la MD ha contribuido significativamente a avances científicos y tecnológicos. Finalmente, el Capítulo 7aborda las tendencias actuales y futuras del campo. Se discute la integración de la inteligencia artificial y el aprendizaje automático en la mejora de potenciales y análisis de resultados, así como los enfoques híbridos que combinan dinámica molecular clásica con métodos cuánticos y redes neuronales. A lo largo del documento, se verá reflejado el carácter interdisciplinar de la dinámica molecular, combinando física, química, matemáticas e informática con el fin de afrontar problemas complejos en ciencia e ingeniería. xii Planificación y Presupuesto La planificación de este proyecto se ha organizado dividiendo el trabajo en diferentes tareas, las cuales se han distribuido a lo largo de las semanas de duración del proyecto. Para mejorar la claridad, dicha planificación se ha representado mediante un diagrama de Gantt, permitiendo así visualizar las dependencias entre tareas y asegurar una gestión eficiente del tiempo. Figura 2.: Diagrama de Gantt del proyecto En el diagrama se puede observar cómo se han planificado y distribuido las distintas tareas del proyecto a lo largo de 22 semanas. En primer lugar, se lleva a cabo un estudio general de los objetivos del proyecto y, a continuación, se realiza una revisión bibliográfica relacionada con dichos objetivos. Luego, se estudian las bases teóricas de la dinámica molecular, profundizando en la modelización, modelos de interacción y representación de sistemas moleculares. Estas tareas se extienden, aproximadamente, hasta la semana 9, y sirven como base para los siguientes bloques del proyecto. Seguidamente, se introduce la etapa en la que se invierte más tiempo: el estudio de los algoritmos empleados en dinámica molecular. Esta fase se desarrolla durante 4semanas, diviendo su contenido en el análisis de los métodos de integración utilizados, los mecanismos de control de temperatura y presión, los algoritmos de optimización y una revisión de las simulaciones paralelas y el uso de aceleración por hardware. Entre las semanas 14 y17 se diseña un algoritmo para realizar una simulación de dinámica molecular, incluyendo el desarrollo del pseudo-código correspondiente y un análisis de las unidades físicas utilizadas en dichas simulaciones. Finalmente, se lleva a cabo un estudio de las aplicaciones de la dinámica molecular, así como de los desafíos actuales y futuras posibles líneas de investigación. Tras dicho estudio, xiii Planificación y Presupuesto se completan las secciones restantes de la memoria y se realiza una revisión general del proyecto. Esta planificación organizada ayuda a que el proyecto avance de forma clara y ordenada, permitiendo ir incorporando conocimientos de forma progresiva. En cuanto al presupuesto del proyecto, se tienen en cuenta los siguientes factores: Coste de personal: El proyecto ha sido desarrollado por una única persona, asumiendo el rol de ingeniero informático junior. Considerando un salario bruto medio de 2.250€ mensuales y una duración estimada del proyecto de 5meses, el coste total de personal asciende a 11.250€. Coste de Hardware y Software: El único coste a considerar en este caso es el del equipo empleado para desarrollar el proyecto, un MSI GF63 Thin 10SCXR, cuyo precio aproximado es de 800€. Por tanto, el coste total del proyecto es de, aproximadamente, 12.050€. xiv 1. Contexto histórico El desarrollo de la dinámica molecular ha estado ligado al avance de la computación y la física teórica. Las primeras simulaciones se dieron en el año 1957, cuando Alder y Wainwright simularon un gas de esferas duras, demostrando la existencia de transiciones de fase mediante métodos computacionales [ 3 ]. Seguidamente, en el año 1964, Rahman realizó la primera simulación con un potencial de Lennard-Jones, modelando el comportamiento del argón líquido [4]. En los años 70 se introdujeron los primeros métodos de integración numérica, como el algoritmo de Verlet [ 5 ], que permitió cálculos más precisos y estables. Con ello, la dinámica molecular se empezó a aplicar en sólidos, líquidos y sistemas biológicos en los años 19761979. Luego, en el año 1985, Car y Parrinello desarrollaron la dinámica molecular ab initio, integrando métodos cuánticos en la simulación clásica [ 6 ]. Dicho método es útil si se busca comparar directamente los resultados de la simulación con mediciones experimentales en materiales específicos, sin embargo, tiene un gran costo computacional. Por ello, para lograr un equilibrio entre precisión y eficiencia computacional, nos centraremos en la dinámica molecular clásica, adecuada para el análisis de fenómenos generales y la comparación de diferentes teorías, siempre que el modelo utilizado represente correctamente los principios físicos esenciales del sistema en estudio. En la década de 1990 y principios de los 2000, la MD se convirtió en una herramienta fundamental para la simulación de proteínas y ADN, impulsando avances en biomedicina y el desarrollo de fármacos. A partir de 2010, la incorporación de la computación de alto rendimiento (HPC) y los aceleradores GPU ha revolucionado el campo, permitiendo la simulación de sistemas con millones de átomos y reduciendo drásticamente los tiempos de cálculo. Mientras que una simulación de 10 nanosegundos solía requerir aproximadamente una semana, ahora puede completarse en tan solo 12 horas, lo que ha posibilitado simulaciones más precisas y extendidas hasta escalas del orden de los microsegundos. Debido al gran impacto que ha tenido la dinámica molecular, siendo clave en proyectos galardonados con el Premio Nobel de Química, como el de Karplus, Levitt y Warshel en 2013 “for the development of multiscale models for complex chemical systems”, hoy en día se están integrando la inteligencia artificial y el aprendizaje automático con el objetivo de mejorar la precisión de los potenciales de interacción, además de explorar nuevos enfoques híbridos, combinando MD con métodos cuánticos y redes neuronales. 1 2. Fundamentos de Dinámica Molecular Para comprender el funcionamiento de la dinámica molecular y los algoritmos presentados en este trabajo, es fundamental conocer sus principios básicos. En este capítulo, se introducen los conceptos esenciales que sustentan la dinámica molecular. 2.1. Principios de la mecánica clásica La dinámica molecular se fundamenta en la mecánica clásica, una rama de la física desarrollada principalmente por Isaac Newton en el siglo XVII, que estudia el movimiento de los cuerpos bajo la acción de fuerzas. Su formulación se apoya en las leyes de Newton [ 7 ], así como en los principios de conservación de la energía y del momento, y en las ecuaciones de movimiento que rigen la evolución de un sistema en el tiempo. 2.1.1. Leyes de Newton En MD, cada átomo o molécula de un sistema se modela como una partícula clásica que se mueve según las leyes de Newton, también conocidas como leyes del movimiento de Newton. Por tanto, las tres leyes del movimiento de Newton son los fundamentos sobre los cuales se basa la dinámica molecular: Primera Ley de Newton La Primera Ley de Newton, conocida como el Principio de Inercia, establece que un cuerpo mantiene su estado de reposo o de movimiento rectilíneo uniforme a menos que una fuerza externa actúe sobre él. En [ 8 ], Newton lo enuncia como “Corpus omne perseverare in statu suo quiescendi vel movendi uniformiter in directum, nisi quatenus illud a viribus impressis cogitur statum suum mutare”, lo que implica que un objeto no cambiará su estado de movimiento a menos que una fuerza lo obligue a hacerlo. Aplicado al contexto de la dinámica molecular, esto significa que, en ausencia de interacciones con otras partículas, los átomos seguirían trayectorias rectilíneas con velocidad constante dentro de la simulación. Segunda Ley de Newton Dado que en la dinámica molecular clásica la masa de cada partícula se considera un parámetro constante, la Segunda Ley de Newton, también conocida como el Principio Fundamental de la Dinámica, establece que la aceleración de un cuerpo es proporcional a la fuerza neta aplicada e inversamente proporcional a su masa. En palabras de Newton: “Mutationem motus proportionalem esse vi motrici impressæ, & fieri secundum lineam rectam qua vis illa imprimitur.” [ 8 ]. Por tanto, la relación fundamental que rige el movimiento es:  F=md v dt=md2 r dt2=m a, (1) 3 2. Fundamentos de Dinámica Molecular donde:  Frepresenta la fuerza aplicada sobre el cuerpo, mes la masa del cuerpo,  ves la velocidad del cuerpo,  res la posición del cuerpo,  aes la aceleración resultante del cuerpo. Tercera Ley de Newton La Tercera Ley de Newton, o Principio de Acción y Reacción, establece que si un cuerpo ejerce una fuerza sobre otro, entonces el segundo cuerpo ejerce una fuerza de igual magnitud, pero dirección opuesta, sobre el primero: “Actioni contrariam semper & æqualem esse reactionem: sive corporum duorum actiones in se mutuo semper esse æquales & in partes contrarias dirigi.” [8] Es decir, sea  Fij la fuerza ejercida por un cuerpo i sobre un cuerpo j , y sea  Fji la fuerza ejercida por el cuerpo jsobre el cuerpo i, entonces  Fij =− Fji. (2) 2.1.2. Principio de conservación de la energía El principio de conservación de la energía establece que, en un sistema cerrado y aislado, la energía no puede crearse ni destruirse, solo transformarse entre sus distintas formas [ 9 ]. En el contexto de la MD, este principio es fundamental para comprender cómo se comportan los sistemas a lo largo del tiempo. En particular, este principio se manifiesta claramente en simulaciones que se realizan bajo las condiciones del colectivo 1 microcanónico (NVE) [ 10 ]. En este tipo de simulaciones, el número de partículas ( N ), el volumen ( V ) y la energía total ( E ) permanecen constantes durante toda la simulación, ya que el sistema está completamente aislado del entorno. Sin embargo, en ocasiones resulta necesario simular condiciones más cercanas a la realidad, lo que requiere introducir mecanismos de control externos, que regulen la temperatura o la presión, por ejemplo. Con este fin, se recurre a los siguientes acomplamientos externos, que permiten la entrada o salida de energía del sistema: Termostato: regula la temperatura ( T ), permitiendo realizar simulaciones en el colectivo canónico (NVT), donde se mantiene constante la temperatura media del sistema [ 11 , 12 ]. Barostato: regula la presión ( P ). Combinado con un termostato, permite simulaciones en el colectivo isóbaro-isotermo (NPT), en el que se conservan tanto la presión media como la temperatura [13,14]. 1 En mecánica estadística una colectividad representa todas las posibles configuraciones microscópicas de un sistema bajo determinadas condiciones externas. 4 2.1. Principios de la mecánica clásica En estos casos, como se ha comentado, el sistema puede intercambiar energía con el exterior, por lo que la energía interna total deja de ser constante. Por tanto, el balance energético general se describe mediante la expresión: UT=Ui+W+Q, donde: UTes la energía interna total del sistema, Uies la energía interna inicial del sistema, Wes el trabajo realizado por o sobre el sistema, Qes el calor añadido o eliminado del sistema. En la Tabla 2.1se muestra un resumen de los colectivos estadísticos más empleados en dinámica molecular, junto con las variables que permanecen constantes, los mecanismos de control necesarios, y algunos ejemplos típicos de aplicación. Colectivo Variables constantes Control externo Ejemplos de uso NVE N,V,ENinguno Análisis energético puro, validación de integradores, estudios teóricos sin influencia externa [15] NVT N,V,TTermostato Procesos biológicos a temperatura constante, simulaciones de proteínas, análisis estructural en equilibrio térmico [11] NPT N,P,TTermostato y barostato Estudio de fases (cristalización, fusión), compresión, materiales a presión ambiente, simulaciones biomoleculares con entorno acuoso [14] Tabla 2.1.: Colectividades estadísticas comunes en dinámica molecular. Por último, es importante señalar que la elección del colectivo depende del tipo de proceso que se desea simular, de las propiedades físicas que se quieren medir y de la disponibilidad de datos experimentales con los que contrastar los resultados. En la práctica, los colectivos NVT y NPT son los más utilizados, ya que permiten replicar con mayor realismo las condiciones experimentales habituales. 2.1.3. Introducción a las ecuaciones de movimiento de Newton Las ecuaciones de movimiento constituyen el núcleo matemático de la dinámica molecular, ya que, a partir de ellas, se obtiene la evolución de las posiciones y velocidades de las partículas a lo largo del tiempo bajo la acción de fuerzas. Dichas ecuaciones derivan directamente de la Segunda Ley de Newton (1). En la práctica, debido a la complejidad que supone resolver estas ecuaciones de forma analítica para sistemas que involucran un gran número de partículas, se recurre a métodos 5 2. Fundamentos de Dinámica Molecular numéricos, que permiten alcanzar soluciones aproximadas mediante la discretización temporal. Estos métodos, fundamentales para el desarrollo de simulaciones computacionales en dinámica molecular, se abordan en detalle en el Capítulo 4. 6 3. Modelización en Dinámica Molecular La modelización constituye una etapa fundamental en cualquier simulación de dinámica molecular, ya que determina el grado de fidelidad con el que el sistema físico real será representado computacionalmente. En esta sección se abordan los principales componentes que conforman la modelización: los modelos de interacción, la representación de los sistemas moleculares y las condiciones de contorno. 3.1. Modelos de interacción Los sistemas están compuestos por partículas (átomos, iones o moléculas), que interaccionan entre sí mediante funciones matemáticas llamadas potenciales de interacción, las cuales están diseñadas para aproximar fuerzas físicas entre las partículas. Estos potenciales describen cómo varía la energía potencial en función de la distancia entre las partículas. Las fuerzas se obtienen a partir de dichos potenciales mediante la siguiente expresión matemática:  Fij =−∇V(rij), (3) donde:  Fij representa la fuerza ejercida por un cuerpo isobre un cuerpo j, V(rij)es el valor del potencial de interacción según la distancia entre los cuerpos iyj. 3.1.1. Campos de fuerza Un campo de fuerza agrupa los distintos potenciales de interacción presentes en un sistema con el objetivo de describir todas las fuerzas internas y externas del mismo. Estas fuerzas se dividen, en general, en dos categorías: interacciones de corto alcance e interacciones de largo alcance. Lennard-Jones Las fuerzas de Van der Waals [ 16 ] son aquellas que se producen entre átomos y moléculas. Tienen carácter atractivo y repulsivo. Las móleculas y átomos se atraen hasta cierta distancia, pero si se acercan demasiado, se repelen. El potencial más comúnmente usado para modelar dichas interacciones es el potencial de Lennard-Jones. Este potencial capta dos contribuciones principales: una repulsiva a distancias muy cortas y otra atractiva a distancias intermedias. Su forma típica es: VLJ(r) = 4ϵσ r12 −σ r6, (4) donde: ϵrepresenta la profundidad del pozo de potencial, 7 3. Modelización en Dinámica Molecular pueda entrar en el radio de corte, la lista debe ser reconstruida para mantener la precisión del modelo. Figura 7.: Distintas etapas de la formación de la lista de Verlet. Se puede observar el radio de corte (círculo sólido) y el radio de vecindad (círculo discontinuo). La lista ha de ser reconstruida antes de que las partículas negras (inicialmente fuera de la lista) entren en el radio de corte [2]. 14 4. Algoritmos en Dinámica Molecular En dinámica molecular, como ya se comentó en el Capítulo 2, el comportamiento de las partículas que componen un sistema se determina mediante la resolución de un conjunto de ecuaciones diferenciales ordinarias que derivan de (1) . Estas ecuaciones describen la evolución de las posiciones y velocidades de dichas partículas a lo largo del tiempo. Sin embargo, ante la imposibilidad de resolverlas de forma analítica debido a la complejidad de los sistemas estudiados y el gran número de grados de libertad implicados, se recurre a algoritmos numéricos que permiten calcular la evolución de las partículas paso a paso, obteniendo así una simulación del sistema. 4.1. Integración de ecuaciones de movimiento Tomando como punto de partida las ecuaciones (1) y (3) , la base matemática de toda simulación de dinámica molecular viene dada por la integración de las siguientes ecuaciones diferenciales: mid2 ri dt2= Fi=−∇ riV( r1, . . .,  rN). La elección de un método de integración adecuado resulta fundamental, ya que debe garantizar la conservación de la energía, estabilidad numérica a largo plazo y un coste computacional razonable. Los métodos numéricos que se presentarán discretizan el tiempo en intervalos de tamaño ∆t , permitiendo así obtener una evolución aproximada de las trayectorias. Por ello, para garantizar la estabilidad numérica a largo plazo, es necesario que el método sea capaz de controlar los errores de redondeo y truncamiento a lo largo de la simulación. Si el tamaño del paso temporal ∆t es demasiado grande, el error acumulado puede aumentar, provocando inestabilidades, como explosiones de energía o trayectorias físicamente irreales. Por tanto, la elección del parámetro ∆t es crucial, ya que debe ser lo suficientemente pequeño como para capturar las dinámicas rápidas del sistema, como las vibraciones intramoleculares, pero lo bastante grande como para optimizar el tiempo de cálculo, pues a menor valor de ∆t, mayor coste computacional. En general, se asume que las aceleraciones (y, por tanto, las fuerzas) se mantienen aproximadamente constantes durante cada paso temporal ∆t , lo que permite que los integradores numéricos utilicen la fuerza calculada al inicio del paso para predecir la evolución del sistema [ 24 ]. De esta forma, cuanto más ligeros sean los átomos, más rápidas serán sus oscilaciones, y menor deberá ser el paso ∆tpara evitar errores significativos. En particular, el límite superior para ∆t está determinado por las frecuencias más altas del sistema, normalmente asociadas a vibraciones intramoleculares rápidas. Una práctica común 15 4. Algoritmos en Dinámica Molecular consiste en seleccionar ∆t al menos 10 veces menor que el período de oscilación más corto del sistema [ 24 ]. Por ejemplo, para una vibración de período τ≈10−14 s, se recomienda tomar ∆t≤10−15 s =1 femtosegundo (fs). Superar este umbral puede provocar un crecimiento exponencial de los errores numéricos. A pesar de ello, existe una cierta tensión entre la necesidad de utilizar pasos de tiempo pequeños para asegurar la precisión numérica y el deseo de alcanzar escalas temporales largas, típicas de muchos procesos microscópicos. Una simulación de 1 microsegundo con ∆t=1 fs requiere del orden de 106 ciclos de integración, lo que implica un coste computacional considerable, ya que en cada paso es necesario recalcular todas las fuerzas. Por ello, se busca elegir el mayor valor de ∆t posible que mantenga la estabilidad numérica y la conservación de la energía [24]. Además, aunque se utilicen representaciones de alta precisión (por ejemplo, doble precisión en punto flotante), los errores de redondeo se acumulan gradualmente a lo largo de la trayectoria, especialmente si se ejecutan simulaciones con millones de pasos. Por ello, la elección adecuada del integrador y del tamaño del paso de integración es fundamental para garantizar la estabilidad y precisión global de la simulación a largo plazo. Propiedades geométricas: métodos simplécticos La dinámica molecular se enmarca dentro del contexto de los sistemas hamiltonianos [ 25 ], los cuales describen la evolución temporal de un sistema mediante ecuaciones que dependen de las coordenadas y los momentos generalizados. Estas ecuaciones derivan de un Hamiltoniano H(  r, p), que representa la energía total del sistema: H(  r, p) = N ∑ i=i ∥ pi∥2 2mi+V(  r1, . . ., rN), y cuya evolución temporal se determina por:            d ri dt =∂H ∂ pi= pi mi= vi, d pi dt =−∂H ∂  ri= Fi, donde: H(  r, p)es la energía total del sistema,  rrepresenta las coordenadas generalizadas2,  p=m· vson los momentos conjugados. 2 En este trabajo, por simplicidad, se considera únicamente el movimiento traslacional de las partículas, es decir,  ri corresponde a las coordenadas espaciales del centro de masa de la partícula i . Sin embargo, en dinámica molecular también pueden intervenir grados de libertad rotacionales y vibracionales, especialmente relevantes en moléculas poliatómicas o sistemas rígidos. 16 4.1. Integración de ecuaciones de movimiento Una propiedad fundamental de estas ecuaciones es su naturaleza simpléctica, lo que implica que las trayectorias del sistema conservan el volumen en el espacio de fases, según el teorema de Liouville [ 25 ]. Es decir, aunque las configuraciones del sistema pueden evolucionar deformando la forma de un volumen en el espacio de fases, el volumen total ocupado permanece constante [ 24 ]. Esta característica tiene importantes consecuencias en el ámbito numérico. Así, un algoritmo de integración de las ecuaciones de Hamilton debe cumplir: Aunque se produzcan pequeñas oscilaciones en la energía o el momento, los errores no se acumulen sistemáticamente, permitiendo conservar estos invariantes a largo plazo. Garantizar la estabilidad estructural del sistema en simulaciones prolongadas. Los integradores simplécticos están diseñados para respetar estas propiedades geométricas, imitando el comportamiento del flujo hamiltoniano exacto, haciéndolos especialmente adecuados para simular sistemas conservativos y estudiar propiedades de equilibrio. De hecho, se ha demostrado que estos métodos pueden interpretarse como soluciones exactas de un sistema hamiltoniano modificado, ligeramente perturbado, cuya forma depende del paso de integración utilizado [ 26 ]. Esta característica explica su capacidad para conservar invariantes del sistema, como la energía o el volumen de fase, a lo largo de millones de pasos. Para ilustrar las diferencias entre integradores simplécticos y no simplécticos, a continuación se comparan los resultados obtenidos al simular un oscilador armónico utilizando los métodos de Euler [ 27 ] (ver Figura 8) y Verlet (ver Figura 9). En todos los casos, se considera una partícula sometida a una fuerza lineal del tipo  F=−k r , donde k es la constante de elasticidad. Figura 8.: Simulación con el método de Euler (no simpléctico). (a) Evolución temporal de la posición; (b) Trayectoria en el espacio de fases. Se observa una trayectoria espiral divergente, que indica una ganancia artificial de energía. 17 4. Algoritmos en Dinámica Molecular Figura 9.: Simulación con el método de Verlet (simpléctico). (a) Evolución temporal de la posición; (b) Trayectoria en el espacio de fases. La energía se conserva a largo plazo y el sistema describe una órbita cerrada. Sin embargo, a pesar de sus ventajas, los integradores simplécticos no son incondicionalmente estables, pues, como se ha comentado anteriormente, existe un límite superior para el tamaño del paso temporal ∆t más allá del cual las simulaciones pueden volverse inestables o generar trayectorias inconsistentes con la física del sistema. Esto puede observarse claramente en la Figura 10, al analizar el comportamiento del algoritmo de Verlet frente a diferentes valores de ∆t . Para pasos suficientemente pequeños, el método conserva su estabilidad y mantiene las propiedades geométricas del sistema. Sin embargo, al superar el umbral crítico, el integrador pierde estabilidad y deja de reflejar adecuadamente la dinámica del sistema. Figura 10.: Evolución temporal de la posición con el algoritmo de Verlet para dos tamaños de paso temporal. (a) ∆t=0.02 , donde el algoritmo es estable; (b) ∆t=2 , donde es inestable. Por tanto, aunque los integradores simplécticos tienen excelentes propiedades de conservación de invariantes geométricos, es fundamental elegir adecuadamente el paso de tiempo 18 4.1. Integración de ecuaciones de movimiento para garantizar su estabilidad numérica. Por otro lado, esta propiedad está estrechamente relacionada con la reversibilidad temporal, pues un integrador simpléctico reversible es capaz de reproducir la trayectoria inversa exacta al invertir los momentos (en ausencia de errores de redondeo). Este comportamiento refleja fielmente la simetría del flujo hamiltoniano, y contribuye a la fiabilidad del método [26]. A continuación se presentan algunos de los algoritmos simplécticos reversibles más utilizados en simulaciones de MD, debido a su equilibrio entre eficiencia, precisión y simplicidad. 4.1.1. Algoritmo de Verlet El algoritmo de Verlet constituye la base de muchos de los métodos más utilizados en dinámica molecular, destacando por su simplicidad y por sus excelentes propiedades de conservación de la energía [ 15 ]. Su formulación matemática se apoya en un desarrollo en serie de Taylor, mediante la cual se aproxima la posición de una partícula en los instantes t+∆t y t−∆t , a partir de su valor en t . Esta simetría temporal lo convierte en un método reversible en el tiempo [26]. En primer lugar, se desarrolla ri(t+∆t)en serie de Taylor alrededor de t:  ri(t+∆t) = ∞ ∑ n=0  ri(n)(t) n!(∆t)n= ri(t) + vi(t)∆t+1 2 ai(t)∆t2+1 6˙  ai(t)∆t3+O(∆t4). (5) Análogamente, se obtiene el desarrollo hacia atrás:  ri(t−∆t) = ∞ ∑ n=0  ri(n)(t) n!(−∆t)n= ri(t)− vi(t)∆t+1 2 ai(t)∆t2−1 6˙  ai(t)∆t3+O(∆t4). (6) Finalmente, combinando (5) y (6) , se cancelan los términos impares y se obtiene la expresión del algoritmo de Verlet: s  ri(t+∆t) = 2 ri(t)− ri(t−∆t) + ∆t2 miFi(t). (7) La aproximación obtenida presenta un error local de orden O(∆t4) , luego, al tratarse de una ecuación diferencial de segundo orden, se tiene un error global de orden O(∆t2) . Por tanto, el algoritmo de Verlet se convierte en una herramienta muy precisa. Además, el algoritmo de Verlet presenta una excelente estabilidad numérica. Al ser un integrador simpléctico, no conserva exactamente la energía total en cada paso, pero las fluctuaciones que introduce permanecen acotadas y no se acumulan de forma sistemática. Esta propiedad resulta especialmente importante en simulaciones realizadas bajo condiciones del conjunto microcanónico (NVE), donde la conservación de la energía es un criterio fundamental [28]. 19 4. Algoritmos en Dinámica Molecular Una de las limitaciones del algoritmo de Verlet es que no actualiza las velocidades de las partículas de forma explícita. No obstante, si se desea, estas se pueden estimar a partir de las posiciones en distintos pasos de tiempo, aprovechando la simetría del método. Para ello, restando las ecuaciones (5)y(6), se obtiene:  ri(t+∆t)− ri(t−∆t) = 2 vi(t)∆t+O(∆t3). Despejando la velocidad:  vi(t) =  ri(t+∆t)− ri(t−∆t) 2∆t+O(∆t2). (8) Proporcionando así una estimación de las velocidades con un error local de orden O(∆t2) , lo cual resulta suficientemente preciso para la mayoría de aplicaciones prácticas, como el cálculo de energía cinética o temperatura en simulaciones bajo condiciones del colectivo microcanónico. 4.1.2. Velocity-Verlet El algoritmo Velocity-Verlet es una de las variantes más populares del método de Verlet, ya que calcula de forma explícita, además de la posición, la velocidad de la partícula en cada paso de tiempo, lo que resulta muy útil para calcular propiedades como la energía cinética o la temperatura [15,29]. Al igual que el algoritmo de Verlet clásico, este algoritmo también se basa en un desarrollo en serie de Taylor de las posiciones y velocidades, convirtiéndose así en un método reversible en el tiempo [26]. El desarrollo en serie de Taylor de la velocidad hacia adelante nos da:  vi(t+∆t) =  vi(t) + ai(t)∆t+O(∆t2). Sin embargo, dicha expresión solo utiliza la aceleración actual, ignorando la que la partícula alcanzará durante el desplazamiento. Por tanto, en busca de obtener una mayor precisión y simetría temporal, se calcula la velocidad en dos pasos: una mitad antes de actualizar la posición, y otra mitad después de calcular la nueva aceleración [30]. 1. Primero se actualiza la velocidad a mitad de paso:  vit+∆t 2= vi(t) + 1 2 ai(t)∆t+O(∆t2). (9) 2. Seguidamente, tomando (9), se actualiza la posición completa:  ri(t+∆t) = ri(t) + vi(t)∆t+1 2 ai(t)∆t2+O(∆t3) = ri(t) + vit+∆t 2∆t+O(∆t3). 20 4.1. Integración de ecuaciones de movimiento 3. Tras actualizar la posición, se recalcula la aceleración en el nuevo tiempo t+∆t a partir de (1):  ai(t+∆t) = 1 mi  Fi(t+∆t). 4. Finalmente, se completa la otra mitad de la velocidad:  vi(t+∆t) =  vit+∆t 2+1 2 ai(t+∆t)∆t+O(∆t3). Obteniendo así las expresiones del algoritmo Velocity-Verlet:  vit+∆t 2= vi(t) + 1 2 ai(t)∆t,  ri(t+∆t) = ri(t) + vit+∆t 2∆t,  vi(t+∆t) =  vit+∆t 2+1 2 ai(t+∆t)∆t. Se obtiene un error local de orden O(∆t3) tanto en el cálculo de la posición como de velocidad, luego, al tratarse de ecuaciones diferenciales de primer orden, presenta un error global de orden O(∆t2)en ambos cálculos. Además, al igual que el algoritmo de Verlet, el algoritmo Velocity-Verlet muestra una excelente estabilidad numérica para sistemas conservativos. Por tanto, el algoritmo Velocity-Verlet representa un perfecto equilibrio entre precisión numérica, fidelidad física y eficiencia computacional, haciéndolo así uno de los algoritmos más empleados en simulaciones de dinámica molecular, destacando sobre todo en simulaciones biomoleculares [30]. 4.1.3. Leap-Frog El método Leap-Frog, o “salto de rana”, debe su nombre al hecho de que las posiciones y velocidades se actualizan en pasos de tiempo intercalados. Aunque esta característica introduce cierta dificultad a la hora de calcular magnitudes como la energía cinética o la temperatura en tiempos enteros, su simpleza y eficiencia lo convierten en una opción adecuada para muchas simulaciones, especialmente en combinación con algunos controladores de temperatura [15,29]. Para obtener la expresión de la velocidad, se parte del desarrollo en serie de Taylor centrada en t para los instantes t±∆t 2 , conviertiéndose así en un método reversible en el tiempo [ 26 ]:  vit+∆t 2= vi(t) + 1 2 ai(t)∆t+1 8˙  ai(t)∆t2+O(∆t3),  vit−∆t 2= vi(t)−1 2 ai(t)∆t+1 8˙  ai(t)∆t2+O(∆t3). 21 4. Algoritmos en Dinámica Molecular Restando ambas expresiones y despejando, se obtiene:  vit+∆t 2= vit−∆t 2+ ai(t)∆t. (10) Una vez conocida la velocidad en el instante t+∆t 2, se actualiza la posición con:  ri(t+∆t) = ri(t) + vit+∆t 2∆t. (11) Las ecuaciones (10) y (11) constituyen el núcleo del algoritmo Leap-Frog, con un error local de orden O(∆t3)tanto en posición como en velocidad. En caso de que se desee obtener la velocidad en un paso de tiempo completo t , se puede aproximar mediante el promedio de las velocidades:  vi(t) = 1 2 vit−∆t 2+ vit+∆t 2. Al igual que el algoritmo Velocity-Verlet, Leap-Frog presenta un error local de orden O(∆t3) tanto en el cálculo de la posición como de velocidad, luego, al tratarse de ecuaciones diferenciales de primer orden, presenta un error global de orden O(∆t2)en ambos cálculos. Por otro lado, el algoritmo Leap-Frog, al ser simpléctico, también presenta una excelente estabilidad numérica para sistemas conservativos. Finalmente, una ventaja destacable del algoritmo Leap-Frog es que proporciona una forma eficiente y estable de integrar las velocidades en pasos intermedios, lo cual puede resultar útil en combinación con ciertos termostatos que requieren acceso frecuente a las velocidades del sistema [28]. Conclusión Los algoritmos de Verlet y sus variantes son ampliamente utilizados en MD debido a: Su derivación simple (por desarrollo de Taylor). Su reversibilidad temporal. Su estructura simpléctica, que, para un valor adecuado de ∆t , garantiza estabilidad y conservación del volumen en el espacio de fases. Su bajo coste computacional, ya que solo requieren una evaluación de fuerzas por paso. En la Tabla 4.1se puede ver un resumen de las propiedades de los algoritmos introducidos. 22 4.2. Control de temperatura y presión Algoritmo numérico Magnitud Error local de truncamiento Error global de truncamiento Reversible Simpléctico Verlet  r∆t4∆t2Sí Sí  v∆t2∆t2 Velocity-Verlet  r∆t3∆t2Sí Sí  v∆t3∆t2 Leap-Frog  r∆t3∆t2Sí Sí  v∆t3∆t2 Tabla 4.1.: Propiedades del algoritmo de Verlet y de sus variantes. 4.2. Control de temperatura y presión Las simulaciones de dinámica molecular no siempre se realizan bajo las condiciones del colectivo microcanónico (NVE), es decir, en condiciones de energía constante. En la práctica, con el fin de representar de forma realista las condiciones experimentales, suele ser necesario trabajar bajo condiciones del colectivo canónico (NVT), manteniendo constante la temperatura, o del isóbaro-isotermo (NPT), manteniendo constante tanto la presión como la temperatura. Para ello, se introducen mecanismos de control externos como termostatos y barostatos, que permiten acoplar el sistema a un baño térmico o barométrico, respectivamente [10,15]. Desde un punto de vista físico-matemático, estos mecanimos modifican las ecuaciones de movimiento del sistema, ya sea mediante la introducción de términos adicionales o mediante un reescalado dinámico de las variables, en busca de establecer las condiciones deseadas. La incorporación de termostatos y barostatos es necesaria porque, sin mecanismos de control, la temperatura y la presión del sistema pueden alejarse considerablemente de los valores deseados, ya sea por fluctuaciones estadísticas o por el propio proceso de inicialización. Por ejemplo, aplicando el teorema de equipartición de la energía para un sistema tridimensional con Npartículas se tiene que: 1 N N ∑ i=1 1 2mi∥ vi∥2=3 2kBT⇒T=1 3NkB N ∑ i=1 mi∥ vi∥2(K), (12) donde: mies la masa de la partícula i,  vies la velocidad de la partícula i, kBes la constante de Boltzmann (≈1.38065 ×10−23 J/K), Tes la temperatura instantánea del sistema. Por tanto, sin termostato, la temperatura del sistema puede derivar progresivamente debido a errores numéricos acumulados, especialmente en simulaciones largas. Esto puede provocar resultados físicamente incorrectos o inconsistentes con condiciones experimentales. Por ejemplo, si se desea modelar una proteína a 300 K pero no se regula la temperatura, la simulación puede desviarse de este valor, afectando a la estructura y la dinámica del sistema [12,31]. 23 4. Algoritmos en Dinámica Molecular artificiales o solapamientos entre átomos que pueden haber sido introducidos durante la construcción inicial del sistema [29]. Por otro lado, los métodos de Monte Carlo ofrecen una alternativa interesante a la dinámica molecular clásica. Mientras que la MD resuelve ecuaciones diferenciales para obtener la evolución temporal del sistema, el enfoque de Monte Carlo introduce aleatoriedad y se centra en el muestreo de configuraciones según una determinada distribución de probabilidad. Este enfoque resulta especialmente útil en estudios termodinámicos o en simulaciones donde la evolución temporal explícita no es relevante [15]. 4.3.1. Minimización de energía La etapa de minimización de energía, como se ha mencionado anteriormente, constituye un paso previo habitual en muchas simulaciones de dinámica molecular. Su objetivo es encontrar una configuración estable del sistema, es decir, una disposición de las partículas que minimice la energía potencial total. Esta configuración corresponde, en general, a un mínimo local de la superficie de energía potencial, ya que el mínimo global suele ser inaccesible computacionalmente para sistemas con muchos grados de libertad [15,29]. Además de garantizar la estabilidad mecánica local, esta etapa es crucial para evitar que fuerzas no físicas generen aceleraciones extremas o inestabilidades numéricas al comienzo de la simulación [ 15 ]. En particular, iniciar la integración desde una configuración alejada del equilibrio puede dar lugar a explosiones numéricas o a errores de integración acumulados. En ciertos contextos, como el estudio del plegamiento de proteínas o el diseño de materiales, se recurre a técnicas de optimización global, como las metaheurísticas, que permiten escapar de mínimos locales y explorar regiones más amplias del paisaje energético [15,36]. Desde un punto de vista matemático, esto equivale a resolver el siguiente problema de optimización: min  r1,..., rN V(  r1, . . ., rN), donde V(  r1, . . ., rN) representa la energía potencial total del sistema en función de las posiciones de los átomos. Para alcanzar dicho mínimo local, muchos de los métodos más utilizados se basan en el gradiente del potencial, ya que los mínimos locales satisfacen las siguientes condiciones de primer y segundo orden: ∇ riV=0, HessV(  ri)>0, (13) donde la primera condición indica que la fuerza neta sobre cada partícula es nula, y la segunda que la matriz Hessiana es definida positiva en ese punto [37]. Recordemos que la matriz Hessiana de Vviene dada por: HessV(  r) =         ∂2V ∂r2 11 ∂2V ∂r11∂r12 ··· ∂2V ∂r11∂rN3 ∂2V ∂r12∂r11 ∂2V ∂r2 12 ··· ∂2V ∂r12∂rN3 . . .. . ..... . . ∂2V ∂rN3∂r11 ∂2V ∂rN3∂r12 ··· ∂2V ∂r2 N3         , (14) 30 4.3. Algoritmos de optimización donde  r= (r11,r12,r13, . . .,rN1,rN2,rN3)T∈R3N , para un sistema de N partículas, cada una con posición tridimensional  ri= (ri1,ri2,ri3) . Lo habitual en el contexto de dinámica molecular es que esta matriz sea simétrica, ya que, en general, se trabaja con funciones suaves [26]. En la práctica, se considera que el sistema ha alcanzado un equilibrio mecánico local cuando la fuerza neta sobre cada átomo es inferior a un umbral prefijado. A continuación se presentan algunos de los métodos de optimización más utilizados en simulaciones de dinámica molecular: Descenso por gradiente Este método, estudiado en la asignatura Aprendizaje Automático, consiste en actualizar iterativamente las posiciones de las partículas en la dirección del gradiente negativo del potencial hasta que la fuerza neta sobre cada átomo quede por debajo del umbral prefijado [15,37]: Algoritmo 1Descenso por gradiente  r← r0 while ∥∇V(  r)∥>εdo  r← r−α∇V(  r) end while return r donde  r0 representa la configuración inicial del sistema, α es el paso de aprendizaje o tasa de actualización, y εel umbral prefijado. Este método es especialmente útil cuando la configuración inicial se encuentra lejos del mínimo local. Sin embargo, el método puede volverse ineficiente al aproximarse al mínimo, ya que si el valor de α es demasiado grande, se producirán oscilaciones, mientras que si el valor de αes muy pequeño, la convergencia será muy lenta (ver Figura 11). Figura 11.: Comparación del descenso por gradiente según la tasa de aprendizaje [38]. 31 4. Algoritmos en Dinámica Molecular Gradiente conjugado El método del gradiente conjugado mejora la eficiencia del descenso por gradiente al generar, en cada iteración, una dirección de búsqueda conjugada respecto a la matriz Hessiana del sistema [37]. Esto permite una convergencia más rápida, especialmente cerca del mínimo. Desde un punto de vista matemático, desarrollando V(  r) alrededor de la configuración inicial r0: V(  r)≈V(  r0) + ∇V(  r0)T(  r− r0) + 1 2(  r− r0)THessV(  r0)(  r− r0), donde HessV(  r0)es la matriz Hessiana (14) evaluada en r0. Ahora, definiendo  x= r− r0 y observando que el término constante V(  r0) no afecta al cálculo del mínimo local (véase (13)), pues su derivada es nula, se escribe: V(  r)≈1 2 xTHessV(  r0) x+∇V(  r0)T x=1 2 xTA x− xTb, tomando A=HessV(  r0) y b=−∇V(  r0) . Por tanto, el método del gradiente conjugado busca minimizar la función cuadrática: f( x) = 1 2 xTA x− xTb. El algoritmo del gradiente conjugado viene dado por: Algoritmo 2Gradiente conjugado  x0←A r0− b  p0← − x0 k←0 while ∥ xk∥>εdo  αk← xT k xk  pT kA pk  rk+1← rk+ αk pk  xk+1← xk+ αkA pk  βk← xT k+1 xk+1  xT k xk  pk+1← − xk+1+ βk+1 pk k←k+1 end while return rk En simulaciones reales se suele emplear una combinación de ambos métodos, comenzando con descenso por gradiente para escapar rápidamente de regiones de alta energía, y continuando con gradiente conjugado cerca del mínimo local [15]. 32 4.3. Algoritmos de optimización Metaheurísticas y optimización global Si bien el descenso por gradiente y el gradiente conjugado son ampliamente utilizados por su eficacia y simplicidad, por lo que son preferibles en simulaciones rutinarias, tienen limitaciones importantes en escenarios donde la función objetivo posee múltiples mínimos locales, como ocurre en biomoléculas plegadas. Para estos casos, han ganado relevancia las metaheurísticas, técnicas de optimización global que no dependen únicamente de información local del gradiente, sino que están diseñadas para explorar globalmente el espacio de búsqueda, permitiendo escapar de mínimos locales a través de mecanismos estocásticos o evolutivos. Entre las más utilizadas se encuentran los algoritmos genéticos y evolutivos [ 39 ], estudiados en la asignatura Metaheurísticas, el método de Basin Hopping [ 36 ] y el descenso estocástico del gradiente (Stochastic Gradient Descent, SGD) [ 40 ], estudiado en la asignatura Aprendizaje Automático. Estas técnicas, aunque no garantizan encontrar el mínimo global, aumentan significativamente la probabilidad de alcanzarlo. Si bien exigen un mayor coste computacional en comparación con los métodos deterministas basados en gradientes, su uso está justificado en situaciones donde es prioritario encontrar mínimos globales, como en el estudio del plegamiento de proteínas o el diseño de materiales. En resumen, la elección del método de optimización depende de la topología de la superficie de energía potencial, del conocimiento previo del sistema y de los recursos computacionales disponibles. 4.3.2. Métodos de Monte Carlo El enfoque de Monte Carlo (MC), diseñado principalmente para el colectivo canónico (NVT), ofrece una alternativa eficaz y conceptualmente diferente a la dinámica molecular tradicional. Mientras que la dinámica molecular integra las ecuaciones de movimiento para obtener la evolución temporal de las partículas, los métodos de MC generan configuraciones aleatorias del sistema y deciden si se aceptan o no según un criterio probabilístico. El objetivo es muestrear configuraciones de acuerdo con una distribución de probabilidad del tipo Boltzmann: P(  r) = 1 Zexp−V(  r) kBT, (15) donde Zes la constante de normalización o función de partición: Z=Zexp−V(  r) kBTd r. (16) A partir de esta expresión, las configuraciones con menor energía tienen mayor probabilidad de ser aceptadas. Sin embargo, también se permite aceptar configuraciones con mayor energía para evitar quedar atrapados en mínimos locales. Los métodos de Monte Carlo tienen la ventaja de ser más fáciles de implementar y suelen ser más eficientes cuando se trata de muestrear grandes espacios de configuración, especialmente en sistemas con restricciones geométricas, como sólidos, redes cristalinas o polímeros [41]. 33 4. Algoritmos en Dinámica Molecular El algoritmo más conocido en este contexto es el algoritmo de Metropolis [ 42 ], desarrollado en 1953. Este método permite muestrear eficientemente el espacio de configuraciones de un sistema a temperatura constante, generando un conjunto de estados que siguen la distribución de Boltzmann (15) . Su ventaja fundamental es que no requiere calcular la función de partición (16) , una cantidad computacionalmente inalcanzable en sistemas con muchos grados de libertad, ya que implica una integral en un espacio de alta dimensionalidad [ 15 , 43 ]. En su lugar, el algoritmo trabaja con el cociente de probabilidades entre estados consecutivos: P(  r′) P(  r)= 1 Zexp−V(  r′) kBT 1 Zexp−V(  r) kBT=expV(  r)−V(  r′) kBT=exp−∆V kBT=4exp−∆E kBT. Aunque este cociente no representa una probabilidad en sentido estricto, el algoritmo lo emplea como probabilidad de aceptación cuando ∆E>0 . En tal caso, se define p= exp−∆E kBT∈(0,1) , y se compara con un número aleatorio r∈[0,1] , lo que permite al sistema aceptar configuraciones energéticamente desfavorables con cierta probabilidad, facilitando así la salida de mínimos locales. En cambio, si ∆E<0 , se acepta directamente la nueva configuración, asegurando que el sistema puede evolucionar hacia estados de menor energía. Por tanto, su implementación viene dada por: Algoritmo 3Algoritmo de Metropolis Inicializar configuración ral azar for paso = 1hasta Ndo Elegir una partícula ial azar Calcular r′a partir de r, modificando la posición de la partícula i Calcular ∆E=V(  r′)−V(  r) if ∆E<0then Aceptar la nueva configuración: r← r′ else Calcular p=exp−∆E kBT Generar número aleatorio r∈[0,1] if r<pthen Aceptar: r← r′ else Rechazar: mantener configuración actual end if end if end for 4 En el contexto del algoritmo de Metropolis aplicado a simulaciones Monte Carlo bajo el colectivo canónico (NVT), los términos cinéticos no se consideran explícitamente, ya que la energía cinética no influye en la distribución de probabilidad configuracional. Por ello, se asume que el incremento de energía ∆E equivale al cambio en la energía potencial ∆V[15]. 34 4.4. Simulaciones paralelas y aceleración Este método garantiza que, tras un número suficiente de pasos, las configuraciones visitadas siguen la distribución de Boltzmann, lo cual es útil para el cálculo de propiedades termodinámicas en equilibrio. Por tanto, aunque el método no permite estudiar trayectorias o dinámicas temporales, es especialmente útil para explorar estados de equilibrio, detectar configuraciones estables y estimar propiedades macroscópicas del sistema. En resumen, mientras que la minimización de energía busca una configuración estable inicial resolviendo un problema de optimización determinista, los métodos de Monte Carlo exploran distintas configuraciones del sistema mediante una estrategia probabilística. 4.4. Simulaciones paralelas y aceleración Las simulaciones de dinámica molecular requieren una elevada carga computacional debido a la complejidad de los cálculos implicados. A medida que el número de partículas crece o se desea simular procesos en escalas temporales más largas, el número de operaciones necesarias se incrementa de forma exponencial. En particular, el cálculo de fuerzas entre partículas y la integración de las ecuaciones de movimiento representan los principales cuellos de botella computacionales, siendo estos cálculos de complejidad O(N2) en el caso más general [15,29]. Para poder simular sistemas de interés biológico o materiales realistas con millones de átomos durante tiempos físicamente relevantes, es fundamental paralelizar el proceso de simulación. Este objetivo se logra mediante diferentes estrategias de paralelización, que permiten distribuir el trabajo entre múltiples procesadores o unidades de cómputo: Descomposición espacial Una de las técnicas más utilizadas es la descomposición espacial, donde el espacio simulado se divide en subdominios (cajas) y cada procesador es responsable de calcular las fuerzas e integrar las trayectorias de las partículas dentro de su región. Para mantener la coherencia global, los procesadores deben intercambiar información sobre partículas cercanas a los límites de sus dominios, denominadas partículas fantasma, lo cual requiere comunicación eficiente entre nodos [21]. Replicación de datos En sistemas de menor tamaño o con baja comunicación entre procesos, se puede optar por replicar los datos en todos los procesadores, de manera que cada uno tenga una copia completa del sistema. Aunque esto incrementa el uso de memoria, permite minimizar la latencia de comunicación, a cambio de realizar cálculos redundantes [15]. Paralelización por fuerza o por átomo Otra opción es paralelizar por tipo de cálculo, asignando por ejemplo a cada núcleo una fracción de las interacciones a calcular (paralelización por fuerzas) o el seguimiento de un conjunto de partículas específicas (paralelización por átomos). Aunque menos escalables, 35 4. Algoritmos en Dinámica Molecular estas estrategias pueden ser útiles en arquitecturas específicas o en algoritmos especializados [26]. Uso de arquitecturas HPC y GPUs El uso de arquitecturas de computación de alto rendimiento (High-Performance Computing, HPC) es fundamental para llevar a cabo simulaciones de dinámica molecular a gran escala. Sistemas como supercomputadores o clústeres permiten dividir el trabajo entre múltiples núcleos o nodos mediante bibliotecas paralelas como OpenMP [ 44 ] o MPI (Message Passing Interface) [45], acelerando así la ejecución de las simulaciones. Por otro lado, las unidades de procesamiento gráfico (GPUs) han supuesto un avance significativo en este ámbito. A diferencia de las CPUs tradicionales, que están optimizadas para ejecutar unos pocos hilos de manera muy eficiente, las GPUs están diseñadas para manejar miles de hilos de manera simultánea. Esta arquitectura masivamente paralela es ideal para realizar operaciones vectoriales repetitivas, como el cálculo de fuerzas entre partículas, que constituye el núcleo computacional de la dinámica molecular [46]. Gracias a estas ventajas, muchos programas populares de simulación, como AMBER, GROMACS [ 47 ] o NAMD, han incorporado versiones optimizadas para su ejecución en GPU, permitiendo alcanzar aceleraciones de hasta un orden de magnitud (es decir, 10 veces más rápidas) con respecto a implementaciones basadas únicamente en CPU. En conclusión, la paralelización y el uso de arquitecturas aceleradas son elementos fundamentales en la dinámica molecular actual. No solo permiten simular sistemas más grandes o durante más tiempo, sino que permiten llevar a cabo estudios que serían computacionalmente inviables con recursos secuenciales. 36 5. Implementación de una Simulación de Dinámica Molecular Una vez establecidos los conceptos teóricos fundamentales de la dinámica molecular en los capítulos anteriores, una de las mejores formas de comprender una simulación es plantear cómo se estructuraría un programa simple bajo las siguientes condiciones: Condiciones termodinámicas y parámetros característicos de la simulación La simulación presentada en esta sección se lleva a cabo bajo las condiciones del colectivo microcanónico (NVE). En este conjunto estadístico se mantiene constante el número de partículas N, el volumen Vy la energía total E(ver Capítulo 2). Durante la evolución temporal, se calculan diversas propiedades termodinámicas del sistema, como la temperatura instantánea y la energía total por partícula. Unidades reducidas En simulaciones de dinámica molecular es habitual emplear un sistema de unidades reducidas con el fin de simplificar los cálculos, mejorar la estabilidad numérica y evitar errores derivados de trabajar con constantes físicas muy pequeñas (por ejemplo, kB∼10−23 J/K). Esta estrategia también permite escribir código más limpio y eficiente, eliminando la necesidad de introducir factores dimensionales en cada paso del algoritmo [15]. La idea fundamental consiste en escoger magnitudes características del sistema como escalas de referencia. En particular, cuando se emplea el potencial de Lennard-Jones, se fijan como unidades base las siguientes constantes: σ: longitud característica (diámetro efectivo de las partículas), se mide en metros [m]. ε: energía de interacción (profundidad del pozo de potencial), se mide en julios [J]. m: masa de las partículas, se mide en kilogramos [kg]. kB: constante de Boltzmann. Esto permite expresar todas las demás magnitudes relevantes en forma adimensional, es decir, sin unidades físicas explícitas. A partir de estas definiciones base, se pueden reescalar el resto de variables del sistema, tal y como se resume en la Tabla 5.1. Por ejemplo, supongamos que en una simulación se emplean los valores σ=ε=m= kB=1 (elección habitual en unidades reducidas), y se obtiene una temperatura reducida T∗=1.2 . Entonces, empleando los valores físicos reales del argón, se obtiene la siguiente temperatura instantánea: T=T∗·ε kB≈1.2 ·1.65 ×10−21 1.38 ×10−23 ·J J/K ≈1.98 ×10−21 1.38 ×10−23 K≈143.5K. 37 5. Implementación de una Simulación de Dinámica Molecular De igual forma, si se desea que, en valores físicos reales, el paso de tiempo empleado sea de ∆t=10−15 s, entonces el paso de tiempo tomado en la simulación vendrá dado por: t=t∗·σrm ε=t∗·3.405 ×10−10s6.63 ×10−26 1.65 ×10−21 =10−15 ⇒t∗≈4.636 ×10−4. Magnitud Magnitud reducida Reescalado a unidades físicas t t∗=t σ√m/εt=t∗·σrm ε(s) T T∗=kBT εT=T∗·ε kB(K)  v v∗= v √ε/m v= v∗·rε m(m/s) P P∗=P·σ3 εP=P∗·ε σ3(Pa)  F F∗=  F·σ ε  F= F∗·ε σ(N) Tabla 5.1.: Conversión entre magnitudes físicas y reducidas Aunque las variables numéricas utilizadas en simulaciones con unidades reducidas son adimensionales y no poseen unidades físicas explícitas, sí representan de forma coherente las proporciones y relaciones entre las propiedades del sistema. Esto permite interpretar los resultados de forma cualitativa, comparar distintos sistemas moleculares y, si se desea, convertir los valores obtenidos nuevamente a unidades físicas mediante las expresiones de reescalado correspondientes [15]. Estructura del programa En la Tabla 5.2se muestra la secuencia típica de etapas que componen una simulación de dinámica molecular. Por tanto, el pseudocódigo desarrollado tendrá la siguiente estructura: 1. Lectura de los parámetros que especifican las condiciones iniciales de la simulación (número de partículas, temperatura inicial, paso de tiempo, etc.). 2. Inicialización del sistema: incluye la preparación, el calentamiento y el equilibrado del sistema. 3. Cálculo de las fuerzas que actúan sobre cada partícula. 4. Integración de las ecuaciones de movimiento de Newton. 5. Cálculo y visualización de propiedades termodinámicas del sistema. 38 5.1. Inicialización del sistema Etapa Descripción Preparación del sistema Definición de las coordenadas iniciales y asignación de velocidades, habitualmente mediante una distribución de MaxwellBoltzmann. Calentamiento Escalado de las velocidades a la temperatura deseada. Equilibrado El sistema se lleva a una situación de equilibrio a partir de su configuración inicial. Producción Generación de las trayectorias del sistema, a partir de las cuales se calculan propiedades físicas y termodinámicas. Tabla 5.2.: Etapas del desarrollo de una simulación de dinámica molecular. [48] A partir de esta estructura, se obtiene el siguiente pseudocódigo, que implementa una simulación de dinámica molecular para un sistema atómico tridimensional sencillo, compuesto por Npartículas: Algoritmo 4Esquema general de una simulación de dinámica molecular Leer parámetros ▷Paso 1 Inicializar sistema ▷Paso 2 t←0 while t<tmax do Calcular fuerzas ▷Paso 3 Integrar ecuaciones de movimiento ▷Paso 4 Calcular y mostrar propiedades ▷Paso 5 t←t+∆t end while Las subrutinas inicializar sistema, calcular fuerzas e integrar ecuaciones de movimiento serán descritas en los algoritmos 5,6y7, respectivamente. 5.1. Inicialización del sistema En esta fase se llevan a cabo tres pasos esenciales: la preparación, el calentamiento y el equilibrado del sistema. En primer lugar, se definen las coordenadas iniciales de las partículas, que suelen colocarse sobre una red regular (como una red cúbica) para evitar solapamientos, y se asignan velocidades iniciales extraídas de una distribución aleatoria. Luego, se realiza un reescalado de las velocidades con el objetivo de que la energía cinética total del sistema corresponda a la temperatura deseada, lo que se conoce como calentamiento del sistema. Finalmente, se elimina la posible velocidad del centro de masas y se ajustan las condiciones para alcanzar un estado de equilibrio termodinámico. Estas tareas pueden implementarse de forma conjunta dentro de una única rutina de inicialización (ver Algoritmo 5). A continuación se detallan los cálculos realizados en ella: Preparación del sistema: Las posiciones iniciales  ri de las partículas se colocan sobre una red regular tridimensional (habitualmente cúbica), mediante una función lattice_pos(i) . Este paso evita solapamientos entre partículas que podrían producir fuerzas repulsivas extremadamente grandes en los primeros pasos de simulación [15]. 39 5. Implementación de una Simulación de Dinámica Molecular LAMMPS: Código versátil orientado a la simulación de materiales y polímeros, pero también capaz de modelar sistemas biológicos. Soporta campos de fuerza clásicos y avanzados, lo que le permite simular fenómenos complejos como la formación y ruptura de enlaces o efectos electrónicos. Además, admite simulaciones a nivel atómico y mesoscópico [53]. OpenMM: Librería orientada a la programación de simulaciones en Python con soporte nativo para GPU. Permite construir simulaciones flexibles, usar campos de fuerza estándar y desarrollar algoritmos personalizados de integración o energía [54]. 46 6. Aplicaciones de la Dinámica Molecular La dinámica molecular se ha convertido en una herramienta fundamental en numerosos campos científicos y tecnológicos. Su aplicabilidad se extiende a múltiples disciplinas, desde la biología estructural hasta la ciencia de materiales, lo que la convierte en una técnica con gran relevancia en la ciencia computacional moderna [15,55]. Antes de entrar en las distintas áreas de aplicación, es importante distinguir entre dos grandes tipos de simulaciones: aquellas centradas en la evolución temporal del sistema (como las que estudian transporte, difusión o reacciones dinámicas), y aquellas cuyo objetivo es obtener propiedades de equilibrio, lo que se denomina un estudio termodinámico 5 . En este segundo caso, no interesa cómo cambia el sistema con el tiempo, sino cómo se comporta en promedio, en condiciones de equilibrio. Sin embargo, incluso en estos casos se emplean algoritmos de dinámica molecular, como el integrador de Verlet, con el fin de obtener trayectorias representativas del colectivo estadístico considerado. 6.1. Bioquímica y biología estructural En esta sección se presentan las principales macromoléculas biológicas y se analiza cómo la dinámica molecular permite estudiar sus propiedades, interacciones y funciones, proporcionando una aproximación computacional que complementa las técnicas experimentales clásicas de la biología estructural. Para facilitar la comprensión, se parte de una breve descripción de los componentes clave, antes de abordar sus aplicaciones concretas en simulación computacional. 6.1.1. Moléculas biológicas: conceptos básicos En el contexto biológico, los protagonistas de las simulaciones suelen ser macromoléculas como las siguientes: Proteínas Las proteínas son cadenas de aminoácidos que se pliegan en formas tridimensionales específicas. Este plegamiento determina su función biológica, que puede ir desde actuar como enzimas que catalizan reacciones hasta transportar oxígeno o servir como receptores en membranas celulares. Estudiar su dinámica es clave para entender cómo se pliegan, cómo interactúan con otras moléculas, o cómo ciertas mutaciones afectan a su comportamiento [55,56]. 5 Un estudio termodinámico se centra en calcular propiedades macroscópicas de equilibrio del sistema, sin considerar su evolución temporal. 47 6. Aplicaciones de la Dinámica Molecular Ácidos nucleicos (ADN y ARN) Los ácidos nucleicos son las moléculas encargadas de almacenar y transmitir la información genética. El ADN es conocido por su estructura de doble hélice, pero su conformación puede variar en función del entorno o de su interacción con proteínas. El ARN, además de su papel en la síntesis de proteínas, puede adoptar estructuras tridimensionales complejas con funciones catalíticas o regulatorias [57]. Membranas celulares Las membranas celulares son estructuras flexibles que rodean y delimitan las células, separando su interior del entorno exterior. Están compuestas principalmente por una doble capa de moléculas similares a grasas, llamadas lípidos, que actúan como una barrera semipermeable. Esta barrera no solo protege la célula, sino que también regula qué sustancias pueden entrar o salir. Mediante simulaciones de dinámica molecular, es posible modelar el comportamiento de estas membranas a nivel atómico y estudiar procesos esenciales como la difusión de pequeñas moléculas, el transporte controlado de sustancias, o cómo ciertas proteínas se insertan y funcionan dentro de la membrana [58]. Complejos biomoleculares Los complejos biomoleculares son conjuntos de varias moléculas (proteínas, ADN, lípidos) que se ensamblan de forma funcional. Por ejemplo, un ribosoma o un canal iónico. Simular su comportamiento permite comprender mecanismos dinámicos que a menudo no son observables directamente mediante técnicas experimentales [57]. 6.1.2. Aplicaciones en simulación biomolecular Gracias a la MD, es posible investigar numerosos procesos relevantes en biología y bioquímica, que tienen un papel central en el diseño de fármacos. A continuación, se presentan algunas de las aplicaciones más relevantes de la dinámica molecular en este contexto: Plegamiento de proteínas El plegamiento es el proceso mediante el cual una cadena lineal de aminoácidos adopta su conformación tridimensional funcional. Las simulaciones se centran en la evolución temporal del sistema, permitiendo explorar este fenómeno paso a paso, identificar intermedios estructurales e incluso estudiar plegamientos erróneos asociados a enfermedades neurodegenerativas como el Alzheimer o el Parkinson [56]. Dinámica estructural Las proteínas y los ácidos nucleicos no son estructuras rígidas, sino que sufren cambios térmicos y reordenamientos que pueden activar o inhibir sus funciones. La dinámica molecular se centra en la evolución temporal del sistema, permitiendo estudiar estos cambios en función del tiempo y en distintos contextos fisiológicos. 48 6.1. Bioquímica y biología estructural Interacción ligando-receptor Uno de los usos más relevantes de la dinámica molecular en farmacología es el estudio del acoplamiento molecular, que describe cómo una molécula pequeña (ligando) se une a una proteína diana (receptor). Este análisis se centra en la evolución temporal del sistema, permitiendo predecir afinidades de unión, identificar sitios activos y guiar el diseño racional de fármacos [59]. Estos mecanismos fundamentales son aprovechados conjuntamente en estrategias de descubrimiento de fármacos: Diseño de fármacos y descubrimiento molecular Uno de los usos más avanzados de la dinámica molecular en biomedicina es el diseño racional de fármacos. Esta estrategia, conocida como drug discovery in silico, consiste en predecir mediante simulaciones qué compuestos podrían unirse eficazmente a una proteína implicada en una enfermedad, la cual actúa como diana biológica. A diferencia de los métodos tradicionales de cribado experimental, que requieren analizar físicamente miles de moléculas, la simulación permite filtrar virtualmente grandes bibliotecas de ligandos, evaluando su interacción con dianas biológicas en función de propiedades estructurales, energéticas y dinámicas [56,59]. Una técnica ampliamente utilizada en el diseño de fármacos es el docking, que consiste en predecir computacionalmente cómo se une un ligando a un receptor. El objetivo es determinar la orientación y posición más probable del ligando dentro del sitio activo del receptor, evaluando su ajuste espacial y la energía asociada a la interacción. Sin embargo, este procedimiento suele asumir modelos rígidos, en los que tanto el ligando como el receptor permanecen estáticos durante el proceso de acoplamiento. Esta simplificación limita la capacidad del docking para capturar la flexibilidad estructural de las moléculas y la complejidad real del entorno molecular. Al combinar docking con simulaciones de dinámica molecular, es posible refinar las predicciones iniciales y estudiar cómo evoluciona la interacción a lo largo del tiempo. Este enfoque combinado permite detectar inestabilidades en la unión, explorar conformaciones alternativas del complejo y validar si la interacción predicha se mantiene estable en un entorno solvado y dinámico [ 57 ]. Esto es especialmente útil para reducir falsos positivos, es decir, aquellos casos en los que una predicción inicial sugiere una buena afinidad entre un ligando y una proteína, pero que en realidad no resulta estable cuando se simula con más realismo. Una vez identificadas moléculas con actividad prometedora frente a una diana biológica, estas se consideran compuestos líderes, los cuales representan puntos de partida clave en el proceso de desarrollo de fármacos. Las simulaciones de dinámica molecular permiten analizar en detalle cómo interactúa un compuesto líder con su diana, observando los contactos atómicos, los ajustes estructurales y la estabilidad de la unión. A partir de esta información, es posible diseñar modificaciones estructurales (por ejemplo, añadir o sustituir grupos químicos específicos) para mejorar su eficacia, la selectividad (es decir, que actúe únicamente sobre la diana deseada) y otras 49 6. Aplicaciones de la Dinámica Molecular propiedades farmacológicas como la solubilidad o la estabilidad en el organismo. Este proceso, conocido como optimización del compuesto líder, es esencial para avanzar desde un candidato preliminar hasta un fármaco potencial. Finalmente, esta estrategia ayuda a seleccionar candidatos farmacológicos, es decir, moléculas que, tras un filtrado computacional riguroso, presentan mayor probabilidad de éxito en etapas posteriores como los ensayos preclínicos o clínicos. Este enfoque ha sido aplicado con éxito al desarrollo de terapias frente a enfermedades como el cáncer, el VIH, trastornos neurodegenerativos o infecciones virales como el SARS-CoV-2[60]. En los últimos años, la combinación de MD con modelos de inteligencia artificial ha impulsado aún más el diseño computacional de fármacos. Modelos basados en redes neuronales, como ANI [ 61 ] o DeepMD [ 62 ], permiten representar funciones de energía aprendidas a partir de cálculos cuánticos, alcanzando una precisión similar pero con un coste computacional mucho menor [ 63 ]. Estas herramientas permiten generar nuevos compuestos, predecir afinidades de unión o identificar patrones relevantes directamente a partir de grandes bases de datos moleculares. 6.2. Materiales y nanotecnología La dinámica molecular ha adquirido un papel fundamental en el campo de la ciencia de materiales, al permitir simular el comportamiento de materiales a nivel atómico. En lugar de describir los materiales mediante propiedades medias o macroscópicas, la MD permite observar directamente cómo se mueven los átomos, cómo se rompen enlaces o cómo se forman defectos internos bajo diferentes condiciones. Esto es especialmente útil en contextos donde los métodos tradicionales (como modelos continuos o ensayos experimentales) no permiten acceder a escalas nanométricas o tiempos extremadamente cortos. Esta técnica proporciona una herramienta eficaz para investigar la estructura interna de sólidos cristalinos, materiales amorfos, polímeros o aleaciones metálicas [26]. Simulación de propiedades de materiales Mediante simulaciones de dinámica molecular es posible determinar muchas de las propiedades físicas clave de un material. En este contexto, el objetivo principal no es estudiar la evolución temporal del sistema, sino realizar un análisis termodinámico que permita calcular magnitudes a partir del comportamiento promedio del sistema en equilibrio. Por ejemplo, al simular una red cristalina, como la del cobre o el silicio, es posible aplicar una deformación virtual y observar cómo responden los átomos que la componen. A partir de esta respuesta, se pueden calcular magnitudes relevantes como el módulo de elasticidad, que mide la rigidez del material y cuantifica cuánto se deforma ante una tensión aplicada; la energía de cohesión, que representa la energía necesaria para separar los átomos del sólido y refleja la estabilidad interna de la estructura; y el límite elástico, que señala el punto a partir del cual el material deja de comportarse de forma elástica y comienza a sufrir deformaciones permanentes. Estas propiedades permiten caracterizar el comportamiento mecánico del material a nivel atómico, algo fundamental para el diseño y la optimización de nuevos compuestos. Además, se pueden estudiar defectos como dislocaciones (fallos lineales en la estructura), vacantes (átomos que faltan) o intersticiales (átomos extra entre posiciones regulares), y ana50 6.2. Materiales y nanotecnología lizar cómo estos afectan a la resistencia mecánica, la difusión de átomos o la conductividad térmica [64]. Otra aplicación importante es el estudio de procesos de fractura y fatiga. Simulando un material con grietas o sometido a tensiones cíclicas, se puede observar cómo se propagan fallos internos, lo cual es muy útil para diseñar materiales más resistentes o duraderos. Modelado de superficies e interfaces Muchas aplicaciones industriales y tecnológicas implican materiales con superficies o interfaces (por ejemplo, un recubrimiento sobre una pieza metálica). La MD permite simular qué ocurre en esas regiones de frontera, es decir, cómo se adsorben moléculas sobre una superficie, cómo se produce el crecimiento de capas de material (crecimiento epitaxial) o cómo se comportan los materiales en contacto cuando hay diferencias de estructura o composición [65]. Este tipo de modelado es especialmente importante en áreas como la eléctronica, donde se usan capas finas de materiales con propiedades específicas; en la nanofabricación, donde se necesita entender procesos de deposición y litografía a escala atómica; o en fenómenos como la corrosión y la adhesión, en los que se simula cómo interactúan los materiales con el entorno o con otros sólidos. Nanotecnología y materiales avanzados La nanotecnología se basa en diseñar y manipular materiales a escala nanométrica (1nanómetro =10−9 metros). En este contexto, la dinámica molecular permite modelar estructuras novedosas llamadas nanomateriales, que pueden tener propiedades muy diferentes a sus contrapartes macroscópicas. Algunos ejemplos importantes son: Nanotubos de carbono (CNTs): estructuras cilíndricas formadas por átomos de carbono, extremadamente resistentes y ligeras. Con MD se puede estudiar su flexibilidad, resistencia mecánica o capacidad de conducción térmica. Fullerenos y nanocápsulas: esferas de carbono que pueden encapsular otras moléculas (como fármacos), y cuya estabilidad y reactividad se pueden analizar por simulación. Grafeno y materiales 2D: láminas de un átomo de grosor con propiedades electrónicas y mecánicas excepcionales. La MD permite estudiar cómo vibran sus átomos (fonones), cómo se comporta térmicamente o cómo se deforma bajo tensión [66]. Además, la dinámica molecular permite modelar fenómenos como el autoensamblaje de nanopartículas, es decir, cómo ciertas moléculas se organizan espontáneamente formando estructuras útiles. Estos fenómenos son clave en el diseño de nanodispositivos aplicados a la medicina personalizada. La interacción entre nanopartículas y sistemas biológicos (por ejemplo, nanopartículas diseñadas para transportar fármacos) también ha sido ampliamente estudiada mediante dinámica molecular, permitiendo predecir aspectos clave como la toxicidad, la absorción celular 51 6. Aplicaciones de la Dinámica Molecular o la compatibilidad de los materiales utilizados en nanomedicina. Gracias a su capacidad para incorporar entornos realistas (como agua, membranas o iones), la dinámica molecular permite analizar estas interacciones con un nivel de detalle que complementa y profundiza la información obtenida experimentalmente. [67]. En resumen, la dinámica molecular aporta una herramienta versátil para explorar y diseñar materiales desde su estructura atómica, anticipando su comportamiento real incluso antes de fabricarlos físicamente. Esta capacidad resulta especialmente útil en contextos de investigación aplicada, innovación tecnológica y desarrollo de nuevos materiales con propiedades a medida. 6.3. Química y física de fluidos La dinámica molecular también se emplea para estudiar el comportamiento de fluidos a escala molecular, es decir, considerando directamente la interacción entre las moléculas individuales que componen el líquido o gas. En este tipo de simulaciones, el objetivo principal es caracterizar el sistema en equilibrio, por lo que se centran en el estudio termodinámico del sistema, permitiendo calcular propiedades macroscópicas como la viscosidad (resistencia al flujo), el coeficiente de difusión (cómo se dispersan las moléculas en el medio), la tensión superficial (energía que se requiere para aumentar la superficie de un líquido) o la conductividad térmica (capacidad del fluido para transferir calor) [15]. Esta forma de simular resulta especialmente útil cuando los modelos tradicionales, que tratan los fluidos como medios continuos, no son suficientes para describir correctamente el sistema. Esto ocurre, por ejemplo, cuando se trabaja a escalas muy pequeñas (nanométricas), donde las propiedades del fluido dependen directamente del comportamiento individual de las moléculas. En estos casos, tener en cuenta la naturaleza discreta de la materia y las interacciones entre partículas es esencial para obtener resultados realistas [68]. Un ejemplo claro es el caso de los líquidos complejos, las soluciones iónicas (como la sal disuelta en agua), los fluidos confinados en canales muy estrechos o las mezclas con muchos tipos de moléculas. En todos ellos, las interacciones específicas entre moléculas pueden dar lugar a efectos que los modelos continuos no predicen, como patrones de organización local, separación de fases o transporte anómalo. Por otro lado, la dinámica molecular también permite estudiar reacciones químicas en medios líquidos utilizando enfoques híbridos, como el método QM/MM (Quantum Mechanics/Molecular Mechanics). En este tipo de simulaciones, una parte del sistema (normalmente donde tiene lugar la reacción) se modela con métodos de química cuántica, mientras que el entorno restante se simula con dinámica clásica (ver Subsección 7.2.1). Este tipo de técnica ha sido utilizado con éxito en el estudio de mecanismos enzimáticos, reacciones ácido-base en disolución o procesos de transferencia electrónica [33,56]. 6.4. Estudios de energía y catalizadores La dinámica molecular también se utiliza ampliamente en el estudio de tecnologías relacionadas con la producción, almacenamiento y conversión de energía. En particular, se ha aplicado 52 6.4. Estudios de energía y catalizadores al modelado de dispositivos como baterías de ion-litio, supercondensadores y celdas de combustible [ 29 ]. Estos sistemas dependen críticamente de cómo se mueven los iones (átomos cargados) dentro del material, de cómo interactúan con los electrodos y de los fenómenos que ocurren en las interfaces entre diferentes fases (por ejemplo, entre un sólido y un líquido). Mediante simulaciones de dinámica molecular es posible observar estos procesos a nivel atómico, proporcionando información sobre la velocidad a la que se difunden los iones, cómo afectan las impurezas o defectos estructurales al transporte iónico, o cómo varía la estructura del material bajo diferentes condiciones de carga o temperatura. Esto es fundamental para diseñar baterías más eficientes y duraderas, optimizar la estabilidad térmica de los materiales activos o mejorar la conductividad de los electrolitos [64]. Otro campo muy importante en el que se aplica la dinámica molecular es la catálisis, es decir, el estudio de materiales que aceleran las reacciones químicas sin consumirse en el proceso. En catálisis heterogénea (donde el catalizador y los reactivos están en fases distintas, como un sólido y un gas), la MD permite estudiar cómo se adsorben los reactivos sobre la superficie del catalizador, cómo cambian de posición y reaccionan, y cómo se liberan los productos. Estas simulaciones ayudan a entender la estructura y dinámica de los llamados “sitios activos”, que son las regiones del material donde ocurren las reacciones [69]. En el caso de la catálisis homogénea (donde el catalizador está disuelto en el mismo medio que los reactivos), la dinámica molecular también es útil para analizar el entorno solvatado 6 del catalizador y cómo influye en su reactividad. Combinada con métodos cuánticos (como el enfoque QM/MM), se pueden simular directamente las transformaciones químicas que ocurren durante la catálisis [33,56]. 6 Proceso por el cual las moléculas del disolvente (por ejemplo, agua) rodean a una especie química disuelta, formando una capa de solvatación. Estas interacciones afectan a la estructura, estabilidad y reactividad de dicha especie. 53 7. Desafíos y futuras direcciones A pesar de su consolidación como herramienta fundamental en la investigación computacional, la dinámica molecular continúa enfrentando importantes desafíos, tanto desde el punto de vista metodológico como computacional. Aunque en las últimas décadas se han logrado avances significativos en técnicas numéricas robustas y en la capacidad computacional disponible, sigue teniendo limitaciones que restringen su aplicabilidad en escenarios más complejos o a escalas mayores. Estas limitaciones, sumadas al crecimiento exponencial de los recursos de hardware y a los recientes avances en algoritmos inteligentes, han abierto nuevas líneas de desarrollo con el objetivo de ampliar la aplicabilidad y mejorar la precisión de esta técnica. En este capítulo se presentan dichas restricciones, se identifican tendencias emergentes en el campo y se exponen diversas oportunidades de mejora que podrían marcar la evolución futura de la MD. 7.1. Limitaciones actuales A continuación, se describen algunas de las principales limitaciones actuales, tanto desde el punto de vista computacional como metodológico: 7.1.1. Escalabilidad Una de las principales limitaciones de la dinámica molecular clásica es su limitada escalabilidad computacional. Si bien los métodos de paralelización y el uso de arquitecturas GPU han permitido avances significativos en rendimiento, las simulaciones a gran escala siguen estando condicionadas por el elevado coste computacional asociado al cálculo de fuerzas y a la integración temporal. Esto se vuelve especialmente restrictivo cuando se desea simular sistemas que contienen millones de átomos o que requieren tiempos de simulación del orden de micro o milisegundos [15,46]. Este cuello de botella impide que muchos fenómenos de interés, como el plegamiento de proteínas, puedan ser simulados con suficiente resolución. Además, incluso en arquitecturas de alto rendimiento, la comunicación entre procesos puede convertirse en un factor limitante, especialmente cuando se emplean esquemas de descomposición espacial con frecuentes intercambios de información entre nodos. 7.1.2. Precisión de los modelos de fuerza Otra limitación crítica proviene de la precisión de los modelos de interacción utilizados, los conocidos campos de fuerza. Los campos de fuerza clásicos, como AMBER, CHARMM o OPLS, se basan en parametrizaciones empíricas que, si bien son eficientes, pueden carecer de precisión cuando se aplican a entornos no estándar, como interfaces complejas, metales o moléculas sintéticas, ya que en estos casos pueden darse tipos de enlace o interacciones que no están correctamente representados en el modelo original [ 29 ]. Esta aproximación limita 55 8. Conclusiones sus posibles dianas terapéuticas. Esta capacidad resulta clave en la búsqueda de tratamientos para enfermedades complejas como el cáncer, el VIH o los trastornos neurodegenerativos, cuya cura representaría, sin duda, un hito histórico en la medicina. Además, la dinámica molecular también se aplica en campos como la ciencia de materiales, la nanotecnología o la biología estructural. Finalmente, se ponen de manifiesto las limitaciones actuales de la dinámica molecular, entre las que destacan la limitada escalabilidad computacional, la falta de precisión de los modelos de campos de fuerza clásicos para describir fenómenos de naturaleza cuántica y los elevados tiempos de simulación requeridos por ciertos experimentos. Frente a estas barreras, se presentan líneas de investigación abiertas y prometedoras, como el uso de inteligencia artificial, la integración de métodos multiescala o el aprovechamiento de arquitecturas de alto rendimiento, que reflejan el gran potencial de la dinámica molecular para seguir evolucionando. 62 A. Código Para apoyar la explicación de la simplecticidad y de la influencia del valor del paso de integración ∆ten la estabilidad del integrador, se han implementado dos scripts: 1import numpy as np 2import matplotlib . pyplot as plt 3 4# Par á metros del sistema 5dt = 0.02 # paso de tiempo 6tend = 20 # tiempo de finalizaci ón de la simulaci ón 7n = int(tend / dt) # n ú meros de iteraciones en la simulaci ón 8x0 , v0 = 0.0, 1.0 # posici ón y velocidades iniciales 9k, m = 1.0 , 1.0 # constante de boltzmann y masa del cuerpo 10 11 # Vector de tiempos 12 tvec = np. linspace (0 , tend , n) 13 14 # Algoritmo de Verlet 15 x_verlet = np . zeros (n) # vector de posici ón 16 v_verlet = np . zeros (n) # vector de velocidad 17 x_verlet [0] = x0 # posici ón inicial 18 v_verlet [0] = v0 # velocidad inicial 19 x_prev = x0 - v0 * dt + 0.5 * (-k * x0 / m) * dt **2 # x(-dt) 20 21 # Integraci ón con Verlet ( posici ón) 22 for iin range (0, n - 1): 23 force = -k * x_verlet [i] 24 x_next = 2 * x_verlet [i] - x_prev + (dt **2 / m) * force 25 x_prev = x_verlet [i] 26 x_verlet [i + 1] = x_next 27 28 # C á lculo de velocidades con Verlet 29 for iin range (1, n - 1): 30 v_verlet [i] = ( x_verlet [i + 1] - x_verlet [i - 1]) / (2 * dt) 31 32 # Recorte para espacio de fases , ya que v[n] no se calcula 33 x_verlet_cut = x_verlet [: -1] 34 v_verlet_cut = v_verlet [: -1] 35 36 # Algoritmo de Euler 37 x_euler = np. zeros (n) # vector de posici ón 38 v_euler = np. zeros (n) # vector de velocidad 39 x_euler [0] = x0 # posici ón inicial 40 v_euler [0] = v0 # velocidad inicial 41 42 # C á lculo de posiciones y velocidades 43 for iin range (n - 1): 63 A. Código 44 x_euler [i + 1] = x_euler [i] + v_euler [i] * dt 45 v_euler [i + 1] = v_euler [i] - k * x_euler [i] * dt 46 47 # Soluci ón exacta 48 x_exact = v0 * np. sin( tvec) + x0 * np. cos( tvec) 49 50 # Grá fica conjunta del algoritmo de Verlet 51 fig , (ax1 , ax2) = plt . subplots (1 , 2, figsize =(12 , 5)) 52 53 # Posici ón 54 ax1 . plot (tvec , x_verlet , label = 'Verlet', linewidth=2) 55 ax1 . plot (tvec , x_exact , '--', label = 'Exacta', linewidth=2) 56 ax1.set_title(f" Algoritmo de Verlet ( dt = { dt })") 57 ax1.set_xlabel(" Tiempo ") 58 ax1.set_ylabel(" Posici ón") 59 ax1 . legend () 60 ax1 . grid (True) 61 62 # Espacio de fases 63 ax2 . plot ( x_verlet_cut , v_verlet_cut , color='tab :blue ') 64 ax2.set_xlabel(" Posici ón") 65 ax2.set_ylabel(" Momento ") 66 ax2.set_title(" Verlet - Espacio de fases ") 67 ax2. axis (" equal ") 68 ax2 . grid (True) 69 70 plt . tight_layout (rect =[0 , 0.03 , 1, 0.95]) 71 plt .show () 72 73 # Grá fica conjunta del algoritmo de Euler 74 fig , (ax1 , ax2 ) = plt . subplots (1, 2, figsize =(12 , 5)) 75 76 # Posici ón 77 ax1 . plot (tvec , x_euler , label ='Euler', linewidth=2) 78 ax1 . plot (tvec , x_exact , '--', label = 'Exacta', linewidth=2) 79 ax1.set_title(f" Algoritmo de Euler ( dt = { dt }) ") 80 ax1.set_xlabel(" Tiempo ") 81 ax1.set_ylabel(" Posici ón") 82 ax1 . legend () 83 ax1 . grid (True ) 84 85 # Espacio de fases 86 ax2 . plot ( x_euler , v_euler , color ='tab:red') 87 ax2.set_xlabel(" Posici ón") 88 ax2.set_ylabel(" Momento ") 89 ax2.set_title(" Euler - Espacio de fases " ) 90 ax2 . axis (" equal ") 91 ax2 . grid (True ) 92 93 plt . tight_layout ( rect =[0 , 0.03 , 1, 0.95]) 94 plt .show () 64 Código A.1: Script para reflejar la diferencia entre los algoritmos simplécticos y no simplécticos. Para ello, se ha simulado un oscilador armónico con el método de Euler (no simpléctico) y con el de Verlet (simpléctico). En ambas simulaciones se muestra la evolución temporal de la posición y la trayectoria en el espacio de fases obtenidas. 1import numpy as np 2import matplotlib . pyplot as plt 3 4def verlet_simulation (ax , dt , t_end =20 , k=1.0 , m=1.0 , x0 =0.0 , v0 =1.0): 5n = int( t_end / dt ) # nú meros de iteraciones en la simulación 6x = np . zeros (n) # vector de posici ón 7v = np . zeros (n) # vector de velocidad 8tvec = np. linspace (0 , t_end , n) # vector de tiempos 9 10 # Condiciones iniciales 11 x[0] = x0 12 x_prev = x0 - v0 * dt + 0.5 * (-k * x0 / m) * dt **2 # x(-dt) 13 14 # Integraci ón con Verlet ( posici ón) 15 for iin range (0 , n - 1) : 16 force = -k * x[i] 17 x_next = 2 * x[i] - x_prev + (dt **2 / m) * force 18 x_prev = x[i] 19 x[i + 1] = x_next 20 21 # Soluci ón exacta 22 x_exact = x0 * np. cos (tvec) + v0 * np. sin ( tvec ) 23 24 # Gr á fico 25 ax . plot (tvec , x , label = " Verlet ") 26 ax . plot ( tvec , x_exact , '--', label = " Exacta ") 27 ax.set_xlabel(" Tiempo ") 28 ax.set_ylabel(" Posici ón ") 29 ax.set_title(f"dt = { dt}") 30 ax . legend () 31 ax . grid ( True ) 32 33 # Crear una figura con 2 subplots 34 fig , (ax1 , ax2 ) = plt . subplots (2 , 1, figsize =(10 , 8) ) 35 36 # Ejecutar para dt = 0.02 y dt = 2 37 verlet_simulation (ax1 , dt =0.02) 38 verlet_simulation (ax2 , dt =2) 39 40 plt . tight_layout ( rect =[0 , 0.03 , 1, 0.95]) 41 plt .show () 65 A. Código Código A.2: Script para comprobar la estabilidad del algoritmo de Verlet para distintos valores del paso de integración ∆t . Para ello, se ha simulado un oscilador armónico y se ha mostrado la evolución temporal obtenida de la posición para ∆t=0.02 y para ∆t=2. El código mostrado se puede encontrar en el siguiente repositorio: https://github.com/Joarpe02/TFG_DinamicaMolecular.git 66 Glosario Adhesión Capacidad de dos materiales diferentes para mantenerse unidos en una interfaz común, debido a interacciones físicas (como fuerzas de Van der Waals) o químicas (como enlaces covalentes o iónicos). Aminoácido Molécula orgánica que contiene un grupo amino ( –NH2 ), un grupo carboxilo ( –COOH ) y una cadena lateral específica (grupo R). Los aminoácidos son las unidades básicas que forman las proteínas mediante enlaces peptídicos. Átomo Es la unidad básica de la materia. Está compuesto por un núcleo, que contiene protones (con carga positiva) y neutrones (sin carga), y por electrones (con carga negativa) que giran alrededor del núcleo. Baño térmico Entorno idealizado que intercambia energía térmica con el sistema simulado, manteniendo su temperatura constante. Baño barométrico Entorno idealizado que permite el intercambio de volumen con el sistema, manteniendo la presión constante. Barostato Mecanismo que regula la presión del sistema en una simulación de dinámica molecular, ajustando el volumen de la celda simulada para mantener la presión deseada. Campo de fuerza (force field) Modelo matemático que describe las fuerzas internas y externas actuando sobre un sistema molecular. Incluye términos para enlaces, ángulos, torsiones, y fuerzas de Van der Waals y electrostáticas. Catálisis Estudio de materiales que aceleran las reacciones químicas sin consumirse en el proceso. Catalizador Sustancia que aumenta la velocidad de una reacción química sin consumirse en el proceso. Centro de masas Punto que representa el movimiento del sistema como un todo. Por ejemplo, si no hay fuerza neta, entonces el centro de masas tiene un movimiento rectilíneo y uniforme. Coarse-graining Técnica de simplificación utilizada en simulaciones en la que se agrupan varios átomos en una única partícula efectiva. Colectivo En mecánica estadística, un colectivo (o ensamble) representa el conjunto de todas las configuraciones microscópicas posibles de un sistema que cumplen ciertas condiciones macroscópicas impuestas externamente, como número de partículas, volumen, energía o temperatura. Colectivo canónico (NVT) Colectivo en el que se mantienen constantes el número de partículas ( N ), el volumen ( V ) y la temperatura ( T ) del sistema. Este tipo de simulación requiere el uso de un termostato. 67 Glosario Colectivo isóbaro-isotermo (NPT) Colectivo en el que se mantienen constantes el número de partículas ( N ), la presión ( P ) y la temperatura ( T ) del sistema. Este tipo de simulación requiere el uso tanto de un termostato como de un barostato. Colectivo microcanónico (NVE) Colectivo en el que se mantienen constantes el número de partículas (N), el volumen (V) y la energía total (E) del sistema. Complejos biomoleculares Conjuntos de macromoléculas biológicas, como proteínas, ácidos nucleicos o lípidos, que interactúan de manera específica y estable para llevar a cabo funciones biológicas concretas. Compuesto líder Molécula con actividad biológica prometedora que actúa sobre una diana terapéutica concreta y que sirve como punto de partida para el desarrollo y optimización de nuevos fármacos. Configuración estable Estado del sistema en el que las partículas se encuentran en un mínimo local de la superficie de energía potencial. Corrosión Proceso químico o electroquímico mediante el cual un material, normalmente metálico, se degrada debido a su interacción con el entorno. CPU (Central Processing Unit) Unidad central de procesamiento de un ordenador. Es el componente encargado de ejecutar instrucciones secuenciales y coordinar el funcionamiento del resto del sistema. Deformación anisotrópica Cambio en la forma y/o volumen del sistema en el que las distintas direcciones espaciales pueden escalarse de manera diferente. Deformación isotrópica Cambio en el volumen de un sistema en el que todas las dimensiones espaciales se escalan por igual, preservando la forma del sistema. Deposición Técnica utilizada en ciencia de materiales y nanotecnología para depositar capas delgadas de material sobre una superficie. Diana biológica Molécula del organismo (generalmente una proteína, como un receptor o una enzima) cuya modulación por parte de un compuesto químico produce un efecto terapéutico deseado. Docking Técnica computacional utilizada para predecir la orientación y afinidad de unión entre una molécula pequeña (ligando) y una macromolécula (receptor o enzima). Enlace peptídico Enlace covalente que une el grupo carboxilo ( –COOH ) de un aminoácido con el grupo amino (–NH2) de otro, liberando una molécula de agua. Entorno solvatado Sistema en el que una o más moléculas (solutos) están rodeadas por moléculas de disolvente, comúnmente agua. Enzima Proteína especializada que actúa como catalizador biológico, acelerando reacciones químicas específicas sin consumirse en el proceso. Ergodicidad Propiedad estadística según la cual el promedio temporal de una magnitud a lo largo de una trayectoria del sistema coincide con su promedio en el colectivo estadístico. 68 Glosario Espacio de fases Espacio matemático multidimensional en el que cada punto representa un estado completo del sistema, especificado por las coordenadas y los momentos de todas las partículas. Fase de equilibrado Etapa de la simulación en la que se lleva al sistema a una situación de equilibrio a partir de su configuración inicial. Fase de producción Etapa de la simulación en la que se generan las trayectorias del sistema, a partir de las cuales se calculan propiedades físicas y termodinámicas. Fluido Sustancia que puede fluir y adaptarse a la forma del recipiente que la contiene. Incluye líquidos, gases y plasmas. Los fluidos carecen de forma fija y pueden sufrir deformaciones continuas bajo la acción de una fuerza, por pequeña que sea. GPU (Graphics Processing Unit) Unidad de procesamiento gráfico especializada en operaciones de cálculo en paralelo. HPC (High-Performance Computing) Conjunto de técnicas y recursos computacionales que permiten resolver problemas de elevada complejidad mediante el uso de arquitecturas paralelas, clústeres de ordenadores o supercomputadores. Ligando Molécula que se une de manera específica y reversible a una macromolécula, como una proteína o un receptor. Lípido Molécula orgánica que desempeña funciones estructurales y de almacenamiento energético. Litografía Técnica de microfabricación utilizada para transferir patrones definidos sobre una superficie mediante el uso de radiación (habitualmente luz ultravioleta). Membranas celulares Estructuras lipídicas que delimitan las células, regulando el transporte de sustancias y la comunicación entre compartimentos. Molécula Agrupación de dos o más átomos unidos mediante enlaces químicos. Metaheurística Estrategia de optimización general que guía y controla algoritmos de búsqueda para resolver problemas complejos donde los métodos deterministas resultan ineficaces. Las metaheurísticas no garantizan encontrar el óptimo global, pero son útiles para explorar y explotar adecuadamente el espacio de búsqueda. Paisaje energético Representación conceptual de la energía potencial de un sistema en función de sus configuraciones microscópicas. Cada punto del paisaje corresponde a un estado del sistema, y la topología del mismo refleja la estabilidad relativa de las configuraciones, así como las trayectorias posibles de transición entre ellas. Paso de integración Unidad de tiempo que separa dos estados consecutivos en una simulación numérica. Período Tiempo que tarda una partícula en completar un ciclo completo de oscilación. En simulaciones de dinámica molecular, se usa como referencia para seleccionar un paso de integración adecuado, asegurando una representación precisa del movimiento atómico. 69 Glosario PBC (Periodic Boundary Conditions) Condiciones de contorno periódicas utilizadas en simulaciones para evitar efectos de borde artificiales. Consisten en replicar la celda de simulación en todas las direcciones del espacio, de modo que cuando una partícula sale por un lado, entra por el lado opuesto. Esto permite simular sistemas infinitos a partir de un número finito de partículas. Potencial de interacción Función matemática que describe la energía potencial entre pares de partículas en función de su distancia o configuración relativa. Ejemplos comunes son el potencial de Lennard-Jones y el potencial de Coulomb. Proteína Molécula biológica formada por cadenas de aminoácidos unidas mediante enlaces peptídicos. Su estructura tridimensional determina su función. Radio de corte Distancia máxima a la que se consideran las interacciones entre partículas en una simulación. Más allá de este radio, las fuerzas (como las de Lennard-Jones o electrostáticas) se despreciarán para reducir el coste computacional, asumiendo que su contribución es insignificante. Radio de vecindad Distancia superior al radio de corte que se utiliza para construir listas de vecinos en simulaciones. Permite anticipar qué partículas podrían entrar en la zona de interacción en los siguientes pasos de la simulación, evitando recalcular interacciones en cada paso y mejorando la eficiencia computacional. Receptor Proteína que reconoce y se une de manera específica a un ligando, desencadenando una respuesta bioquímica o celular. Reversibilidad temporal Propiedad de las ecuaciones de movimiento de la mecánica clásica por la cual, si se invierte el sentido del tiempo y los momentos de las partículas, el sistema sigue una trayectoria compatible con las leyes físicas originales. Esta simetría implica que el sistema puede evolucionar hacia atrás en el tiempo siguiendo las mismas ecuaciones que rigen su evolución hacia adelante. Simplecticidad Propiedad matemática de ciertos integradores que preservan la estructura geométrica del sistema hamiltoniano, evitando la acumulación sistemática de errores. En particular, conservan el volumen en el espacio de fases. Termostato Mecanismo que regula la temperatura del sistema en una simulación de dinámica molecular, asegurando que permanezca próxima a un valor objetivo. Tiempo de relajación Tiempo que tarda un sistema en volver a un estado de equilibrio o cuasi-equilibrio tras una perturbación. Trade-off Compromiso inherente en el diseño de modelos o algoritmos, donde la mejora de una característica (como la precisión) suele implicar una pérdida en otra (como el coste computacional). 70 Bibliografía [1] Dirk P. Kroese, Thomas Taimre, and Zdravko I. Botev. Handbook of Monte Carlo Methods. Wiley Series in Probability and Statistics. John Wiley & Sons, 2011. [2] Norbert Attig, Kurt Binder, Helmut Grubmüller, and Kurt Kremer. Computational Soft Matter: From Synthetic Polymers to Proteins, volume 23 of NIC Series. John von Neumann Institute for Computing, Jülich, Germany, 2004. [3] B. J. Alder and T. E. Wainwright. Phase transition for a hard sphere system. The Journal of Chemical Physics,27(5):1208–1209,1957. [4] Aneesur Rahman. Correlations in the motion of atoms in liquid argon. Physical Review, 136(2A):A405–A411,1964. [5] Loup Verlet. Computer experiments on classical fluids. i. thermodynamical properties of lennardjones molecules. Physical Review,159(1):98–103,1967. [6] Roberto Car and Michele Parrinello. Unified approach for molecular dynamics and densityfunctional theory. Physical Review Letters,55(22):2471–2474,1985. [7] Isaac Newton. The Principia: Mathematical Principles of Natural Philosophy. University of California Press, Berkeley, CA, 1999. Original work published in 1687. [8] Isaac Newton. Philosophiæ Naturalis Principia Mathematica. Royal Society, London, 1687. [9] Paul A. Tipler and Gene Mosca. Physics for Scientists and Engineers. W. H. Freeman and Company, New York, 6th edition, 2008. Chapter 8: Conservation of Energy. [10] M.P. Allen and D.J. Tildesley. Computer Simulation of Liquids. Oxford University Press, Oxford, 1989. [11] Shuichi Nosé. A unified formulation of the constant temperature molecular dynamics methods. The Journal of Chemical Physics,81(1):511–519,1984. [12] William G. Hoover. Canonical dynamics: Equilibrium phase-space distributions. Physical Review A,31(3):1695–1697,1985. [13] Glenn J. Martyna, Douglas J. Tobias, and Michael L. Klein. Constant pressure molecular dynamics algorithms. The Journal of Chemical Physics,101(5):4177–4189,1994. [14] Michele Parrinello and Aneesur Rahman. Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied Physics,52(12):7182–7190,1981. [15] Daan Frenkel and Berend Smit. Understanding Molecular Simulation: From Algorithms to Applications. Academic Press, 2nd edition, 2001. [16] Jacob N. Israelachvili. Intermolecular and surface forces. Academic Press,3,2011. [17] Tom Darden, Darrin York, and Lee Pedersen. Particle mesh ewald: An n·log(n) method for ewald sums in large systems. The Journal of Chemical Physics,98(12):10089–10092,1993. [18] William D. Cornell, Piotr Cieplak, Christopher I. Bayly, Kenneth M. Merz, David M. Ferguson, David C. Spellmeyer, Thomas Fox, James W. Caldwell, and Peter A. Kollman. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. Journal of the American Chemical Society,117(19):5179–5197,1995. [19] Alexander D. MacKerell Jr., D. Bashford, M. Bellott, R. L. Dunbrack Jr., J. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. 71