La dinámica molecular ( MD ) es un método de simulación por computadora para analizar los movimientos físicos de átomos y moléculas . Se permite que los átomos y las moléculas interactúen durante un período de tiempo fijo, lo que proporciona una visión de la "evolución" dinámica del sistema. En la versión más común, las trayectorias de los átomos y las moléculas se determinan resolviendo numéricamente las ecuaciones de movimiento de Newton para un sistema de partículas en interacción, donde las fuerzas entre las partículas y sus energías potenciales a menudo se calculan utilizando potenciales interatómicos o campos de fuerza mecánicos moleculares . El método se aplica principalmente en física química , ciencia de los materiales y biofísica .
Dado que los sistemas moleculares suelen estar compuestos por una gran cantidad de partículas, es imposible determinar las propiedades de sistemas tan complejos de forma analítica; la simulación MD evita este problema mediante el uso de métodos numéricos . Sin embargo, las simulaciones MD largas están mal condicionadas matemáticamente y generan errores acumulativos en la integración numérica que pueden minimizarse con una selección adecuada de algoritmos y parámetros, pero no eliminarse.
En el caso de los sistemas que obedecen a la hipótesis ergódica , la evolución de una simulación de dinámica molecular puede utilizarse para determinar las propiedades termodinámicas macroscópicas del sistema: los promedios temporales de un sistema ergódico corresponden a los promedios de conjunto microcanónicos . La MD también se ha denominado " mecánica estadística por números" y " la visión de Laplace de la mecánica newtoniana " de predecir el futuro animando las fuerzas de la naturaleza [1] y permitiendo comprender el movimiento molecular a escala atómica.
La MD se desarrolló originalmente a principios de la década de 1950, a raíz de éxitos anteriores con simulaciones de Monte Carlo (que se remontan al siglo XVIII, en el problema de la aguja de Buffon , por ejemplo), pero fue popularizada para la mecánica estadística en el Laboratorio Nacional de Los Álamos por Marshall Rosenbluth y Nicholas Metropolis en lo que hoy se conoce como el algoritmo Metropolis-Hastings . El interés en la evolución temporal de los sistemas de N cuerpos data de mucho antes del siglo XVII, comenzando con Isaac Newton , y continuó hasta el siglo siguiente en gran medida con un enfoque en la mecánica celeste y cuestiones como la estabilidad del sistema solar . Muchos de los métodos numéricos utilizados hoy en día se desarrollaron durante este período de tiempo, que es anterior al uso de computadoras; por ejemplo, el algoritmo de integración más común utilizado hoy en día, el algoritmo de integración de Verlet , fue utilizado ya en 1791 por Jean Baptiste Joseph Delambre . Los cálculos numéricos con estos algoritmos pueden considerarse MD hechos "a mano".
Ya en 1941 se llevó a cabo la integración de las ecuaciones de movimiento de muchos cuerpos con computadoras analógicas . Algunas emprendieron el laborioso trabajo de modelar el movimiento atómico mediante la construcción de modelos físicos, por ejemplo, utilizando esferas macroscópicas. El objetivo era organizarlas de tal manera que replicaran la estructura de un líquido y usar esto para examinar su comportamiento. JD Bernal describe este proceso en 1962, escribiendo: [2]
... Tomé varias pelotas de goma y las pegué con varillas de distintas longitudes, que iban desde 7 a 10 centímetros. Intenté hacerlo de la forma más informal posible, trabajando en mi propia oficina, siendo interrumpido cada cinco minutos más o menos y sin recordar lo que había hecho antes de la interrupción.
Tras el descubrimiento de las partículas microscópicas y el desarrollo de las computadoras, el interés se expandió más allá del campo de pruebas de los sistemas gravitacionales hacia las propiedades estadísticas de la materia. En un intento por comprender el origen de la irreversibilidad , Enrico Fermi propuso en 1953, y publicó en 1955, [3] el uso de la computadora MANIAC I , también en el Laboratorio Nacional de Los Álamos , para resolver la evolución temporal de las ecuaciones de movimiento para un sistema de muchos cuerpos sujeto a varias elecciones de leyes de fuerza. Hoy, este trabajo seminal se conoce como el problema de Fermi–Pasta–Ulam–Tsingou . La evolución temporal de la energía a partir del trabajo original se muestra en la figura de la derecha.
En 1957, Berni Alder y Thomas Wainwright utilizaron una computadora IBM 704 para simular colisiones perfectamente elásticas entre esferas duras . [4] En 1960, en quizás la primera simulación realista de la materia, JB Gibson et al . simularon el daño por radiación del cobre sólido utilizando un tipo de interacción repulsiva de Born-Mayer junto con una fuerza de superficie cohesiva . [5] En 1964, Aneesur Rahman publicó simulaciones de argón líquido que usaban un potencial de Lennard-Jones ; los cálculos de las propiedades del sistema, como el coeficiente de autodifusión , se compararon bien con los datos experimentales. [6] Hoy en día, el potencial de Lennard-Jones sigue siendo uno de los potenciales intermoleculares más utilizados . [7] [8] Se utiliza para describir sustancias simples (también conocido como Lennard-Jonesio [9] [10] [11] ) para estudios conceptuales y de modelos y como un bloque de construcción en muchos campos de fuerza de sustancias reales. [12] [13]
Utilizado por primera vez en física teórica , el método de dinámica molecular ganó popularidad en la ciencia de los materiales poco después, y desde la década de 1970 también se ha utilizado comúnmente en bioquímica y biofísica . La MD se utiliza con frecuencia para refinar las estructuras tridimensionales de proteínas y otras macromoléculas basadas en restricciones experimentales de cristalografía de rayos X o espectroscopia de RMN . En física, la MD se utiliza para examinar la dinámica de fenómenos a nivel atómico que no se pueden observar directamente, como el crecimiento de películas delgadas y la subimplantación de iones, y para examinar las propiedades físicas de dispositivos nanotecnológicos que no se han creado o aún no se pueden crear. En biofísica y biología estructural , el método se aplica con frecuencia para estudiar los movimientos de macromoléculas como proteínas y ácidos nucleicos , lo que puede ser útil para interpretar los resultados de ciertos experimentos biofísicos y para modelar interacciones con otras moléculas, como en el acoplamiento de ligandos . En principio, MD se puede utilizar para la predicción ab initio de la estructura de proteínas simulando el plegamiento de la cadena polipeptídica a partir de una bobina aleatoria .
Los resultados de las simulaciones de MD se pueden probar mediante la comparación con experimentos que miden la dinámica molecular, de los cuales un método popular es la espectroscopia de RMN. Las predicciones de estructura derivadas de MD se pueden probar a través de experimentos de toda la comunidad en la Evaluación Crítica de la Predicción de la Estructura de Proteínas ( CASP ), aunque el método ha tenido históricamente un éxito limitado en esta área. Michael Levitt , quien compartió el Premio Nobel en parte por la aplicación de MD a las proteínas, escribió en 1999 que los participantes de CASP generalmente no usaban el método debido a "... una vergüenza central de la mecánica molecular , a saber, que la minimización de la energía o la dinámica molecular generalmente conduce a un modelo que es menos parecido a la estructura experimental". [14] Las mejoras en los recursos computacionales que permiten más y más largas trayectorias de MD, combinadas con mejoras modernas en la calidad de los parámetros del campo de fuerza , han producido algunas mejoras tanto en la predicción de la estructura como en el refinamiento del modelo de homología , sin llegar al punto de utilidad práctica en estas áreas; muchos identifican los parámetros del campo de fuerza como un área clave para un mayor desarrollo. [15] [16] [17]
Se ha informado sobre la simulación MD para el desarrollo de farmacóforos y el diseño de fármacos . [18] Por ejemplo, Pinto et al . implementaron simulaciones MD de complejos Bcl-xL para calcular las posiciones promedio de los aminoácidos críticos involucrados en la unión del ligando. [19] Carlson et al . implementaron simulaciones de dinámica molecular para identificar compuestos que complementan un receptor mientras causan una interrupción mínima en la conformación y flexibilidad del sitio activo. Se superpusieron instantáneas de la proteína a intervalos de tiempo constantes durante la simulación para identificar regiones de unión conservadas (conservadas en al menos tres de once marcos) para el desarrollo de farmacóforos. Spyrakis et al . se basaron en un flujo de trabajo de simulaciones MD, huellas dactilares para ligandos y proteínas (FLAP) y análisis discriminante lineal (LDA) para identificar las mejores conformaciones ligando-proteína para actuar como plantillas de farmacóforos basándose en el análisis ROC retrospectivo de los farmacóforos resultantes. En un intento de mejorar el modelado del descubrimiento de fármacos basado en la estructura, en vista de la necesidad de muchos compuestos modelados, Hatmal et al . propusieron una combinación de simulación MD y análisis de contactos intermoleculares ligando-receptor para discernir los contactos intermoleculares críticos (interacciones de unión) de los redundantes en un único complejo ligando-proteína. Los contactos críticos pueden luego convertirse en modelos de farmacóforos que pueden usarse para el cribado virtual. [20]
Un factor importante son los enlaces de hidrógeno intramoleculares , [21] que no se incluyen explícitamente en los campos de fuerza modernos, sino que se describen como interacciones de Coulomb de cargas puntuales atómicas . [ cita requerida ] Esta es una aproximación burda porque los enlaces de hidrógeno tienen una naturaleza parcialmente mecánica cuántica y química . Además, las interacciones electrostáticas generalmente se calculan utilizando la constante dieléctrica del vacío , aunque la solución acuosa circundante tiene una constante dieléctrica mucho más alta. Por lo tanto, el uso de la constante dieléctrica macroscópica a distancias interatómicas cortas es cuestionable. Finalmente, las interacciones de van der Waals en MD generalmente se describen mediante potenciales de Lennard-Jones [22] [23] basados en la teoría de Fritz London que solo es aplicable en el vacío. [ cita requerida ] Sin embargo, todos los tipos de fuerzas de van der Waals son en última instancia de origen electrostático y, por lo tanto, dependen de las propiedades dieléctricas del entorno . [24] La medición directa de las fuerzas de atracción entre diferentes materiales (como la constante de Hamaker ) muestra que "la interacción entre hidrocarburos a través del agua es aproximadamente el 10% de la que se produce a través del vacío". [24] La dependencia del entorno de las fuerzas de van der Waals se descuida en las simulaciones estándar, pero se puede incluir desarrollando campos de fuerza polarizables.
El diseño de una simulación de dinámica molecular debe tener en cuenta la potencia computacional disponible. El tamaño de la simulación ( n = número de partículas), el intervalo de tiempo y la duración total del tiempo deben seleccionarse de modo que el cálculo pueda finalizar en un período de tiempo razonable. Sin embargo, las simulaciones deben ser lo suficientemente largas como para ser relevantes para las escalas de tiempo de los procesos naturales que se están estudiando. Para sacar conclusiones estadísticamente válidas de las simulaciones, el lapso de tiempo simulado debe coincidir con la cinética del proceso natural. De lo contrario, es análogo a sacar conclusiones sobre cómo camina un humano cuando solo se observa menos de un paso. La mayoría de las publicaciones científicas sobre la dinámica de las proteínas y el ADN [25] [26] utilizan datos de simulaciones que abarcan desde nanosegundos (10 −9 s) hasta microsegundos (10 −6 s). Para obtener estas simulaciones, se necesitan varios días de CPU a años de CPU. Los algoritmos paralelos permiten distribuir la carga entre las CPU ; un ejemplo es el algoritmo de descomposición espacial o de fuerza. [27]
Durante una simulación clásica de dinámica molecular, la tarea que requiere más recursos de la CPU es la evaluación del potencial como función de las coordenadas internas de las partículas. Dentro de esa evaluación de energía, la más costosa es la parte no enlazada o no covalente. En la notación O grande , las simulaciones de dinámica molecular comunes escalan según si se deben tener en cuenta explícitamente todas las interacciones electrostáticas y de van der Waals por pares . Este costo computacional se puede reducir empleando métodos electrostáticos como la suma de Ewald de malla de partículas ( ), la malla partícula-partícula-partícula ( P 3 M ) o buenos métodos de corte esférico ( ). [ cita requerida ]
Otro factor que afecta el tiempo total de CPU necesario para una simulación es el tamaño del paso de tiempo de integración. Este es el tiempo transcurrido entre evaluaciones del potencial. El paso de tiempo debe elegirse lo suficientemente pequeño para evitar errores de discretización (es decir, menor que el período relacionado con la frecuencia vibracional más rápida en el sistema). Los pasos de tiempo típicos para la MD clásica son del orden de 1 femtosegundo (10 −15 s). Este valor se puede extender mediante el uso de algoritmos como el algoritmo de restricción SHAKE , que fija las vibraciones de los átomos más rápidos (por ejemplo, hidrógenos) en su lugar. También se han desarrollado métodos de escala de tiempo múltiple, que permiten tiempos extendidos entre actualizaciones de fuerzas de largo alcance más lentas. [28] [29] [30]
Para simular moléculas en un disolvente , se debe elegir entre un disolvente explícito e implícito . Las partículas de disolvente explícito (como los modelos de agua TIP3P , SPC/E y SPC-f ) deben calcularse de forma costosa mediante el campo de fuerza, mientras que los disolventes implícitos utilizan un enfoque de campo medio . El uso de un disolvente explícito es costoso desde el punto de vista computacional, ya que requiere la inclusión de aproximadamente diez veces más partículas en la simulación. Pero la granularidad y la viscosidad del disolvente explícito son esenciales para reproducir ciertas propiedades de las moléculas de soluto. Esto es especialmente importante para reproducir la cinética química .
En todo tipo de simulaciones de dinámica molecular, el tamaño del cuadro de simulación debe ser lo suficientemente grande como para evitar artefactos de condiciones de contorno . Las condiciones de contorno se tratan a menudo eligiendo valores fijos en los bordes (que pueden causar artefactos) o empleando condiciones de contorno periódicas en las que un lado de la simulación vuelve al lado opuesto, imitando una fase en masa (que también puede causar artefactos).
En el conjunto microcanónico , el sistema está aislado de los cambios en moles (N), volumen (V) y energía (E). Corresponde a un proceso adiabático sin intercambio de calor. Una trayectoria de dinámica molecular microcanónica puede verse como un intercambio de energía potencial y cinética, con conservación de la energía total. Para un sistema de N partículas con coordenadas y velocidades , el siguiente par de ecuaciones diferenciales de primer orden puede escribirse en la notación de Newton como
La función de energía potencial del sistema es una función de las coordenadas de las partículas . En física, se la denomina simplemente potencial o campo de fuerza en química. La primera ecuación proviene de las leyes de movimiento de Newton ; la fuerza que actúa sobre cada partícula del sistema se puede calcular como el gradiente negativo de .
Para cada paso de tiempo, la posición y velocidad de cada partícula se pueden integrar con un método de integración simpléctica como la integración de Verlet . La evolución temporal de y se denomina trayectoria. Dadas las posiciones iniciales (por ejemplo, a partir del conocimiento teórico) y las velocidades (por ejemplo, gaussianas aleatorias ), podemos calcular todas las posiciones y velocidades futuras (o pasadas).
Una fuente frecuente de confusión es el significado de la temperatura en MD. Comúnmente tenemos experiencia con temperaturas macroscópicas, que involucran una gran cantidad de partículas, pero la temperatura es una cantidad estadística. Si hay una cantidad suficientemente grande de átomos, la temperatura estadística se puede estimar a partir de la temperatura instantánea , que se obtiene al igualar la energía cinética del sistema a nk B T /2, donde n es el número de grados de libertad del sistema.
Un fenómeno relacionado con la temperatura surge debido al pequeño número de átomos que se utilizan en las simulaciones MD. Por ejemplo, considere simular el crecimiento de una película de cobre comenzando con un sustrato que contiene 500 átomos y una energía de deposición de 100 eV . En el mundo real, los 100 eV del átomo depositado se transportarían rápidamente y se compartirían entre una gran cantidad de átomos ( o más) sin un gran cambio en la temperatura. Sin embargo, cuando solo hay 500 átomos, el sustrato se vaporiza casi inmediatamente por la deposición. Algo similar sucede en las simulaciones biofísicas. La temperatura del sistema en NVE aumenta naturalmente cuando las macromoléculas como las proteínas experimentan cambios conformacionales exotérmicos y unión.
En el conjunto canónico , se conservan la cantidad de sustancia (N), el volumen (V) y la temperatura (T). A veces también se denomina dinámica molecular a temperatura constante (CTMD). En NVT, la energía de los procesos endotérmicos y exotérmicos se intercambia con un termostato .
Hay una variedad de algoritmos de termostato disponibles para agregar y quitar energía de los límites de una simulación MD de una manera más o menos realista, aproximándose al conjunto canónico . Los métodos populares para controlar la temperatura incluyen el reescalamiento de velocidad, el termostato Nosé-Hoover , las cadenas Nosé-Hoover, el termostato Berendsen , el termostato Andersen y la dinámica de Langevin . El termostato Berendsen podría introducir el efecto del cubo de hielo volador , que conduce a traslaciones y rotaciones no físicas del sistema simulado.
No es fácil obtener una distribución de conjuntos canónicos de conformaciones y velocidades utilizando estos algoritmos. La forma en que esto depende del tamaño del sistema, la elección del termostato, los parámetros del termostato, el paso de tiempo y el integrador es el tema de muchos artículos en el campo.
En el conjunto isotérmico-isobárico se conservan la cantidad de sustancia (N), la presión (P) y la temperatura (T). Además de un termostato, se necesita un barostato . Se corresponde más estrechamente con las condiciones de laboratorio con un matraz abierto a temperatura y presión ambiente.
En la simulación de membranas biológicas , el control de la presión isotrópica no es adecuado. En el caso de las bicapas lipídicas , el control de la presión se produce con un área de membrana constante (NPAT) o una tensión superficial constante "gamma" (NPγT).
El método de intercambio de réplicas es un conjunto generalizado. Fue creado originalmente para abordar la dinámica lenta de los sistemas de espín desordenado. También se lo denomina templado paralelo. La formulación de MD de intercambio de réplicas (REMD) [31] intenta superar el problema de los mínimos múltiples intercambiando la temperatura de réplicas no interactuantes del sistema que funcionan a varias temperaturas.
Una simulación de dinámica molecular requiere la definición de una función potencial o una descripción de los términos mediante los cuales interactuarán las partículas en la simulación. En química y biología, esto se conoce generalmente como campo de fuerza y, en física de materiales, como potencial interatómico . Los potenciales se pueden definir en muchos niveles de precisión física; los que se utilizan con más frecuencia en química se basan en la mecánica molecular y representan un tratamiento de la mecánica clásica de las interacciones entre partículas que pueden reproducir cambios estructurales y conformacionales , pero que, por lo general, no pueden reproducir reacciones químicas .
La reducción de una descripción completamente cuántica a un potencial clásico implica dos aproximaciones principales. La primera es la aproximación de Born-Oppenheimer , que establece que la dinámica de los electrones es tan rápida que se puede considerar que reaccionan instantáneamente al movimiento de sus núcleos. En consecuencia, se pueden tratar por separado. La segunda trata los núcleos, que son mucho más pesados que los electrones, como partículas puntuales que siguen la dinámica newtoniana clásica. En la dinámica molecular clásica, el efecto de los electrones se aproxima como una superficie de energía potencial, que generalmente representa el estado fundamental.
Cuando se necesitan niveles más finos de detalle, se utilizan potenciales basados en la mecánica cuántica ; algunos métodos intentan crear potenciales híbridos clásicos/cuánticos donde la mayor parte del sistema se trata de manera clásica pero una pequeña región se trata como un sistema cuántico, generalmente experimentando una transformación química.
Los potenciales empíricos utilizados en química se denominan frecuentemente campos de fuerza , mientras que los utilizados en física de materiales se denominan potenciales interatómicos .
La mayoría de los campos de fuerza en química son empíricos y consisten en una suma de fuerzas enlazadas asociadas con enlaces químicos , ángulos de enlace y diedros de enlace , y fuerzas no enlazadas asociadas con fuerzas de van der Waals y carga electrostática . [32] Los potenciales empíricos representan efectos mecánico-cuánticos de manera limitada a través de aproximaciones funcionales ad hoc. Estos potenciales contienen parámetros libres como carga atómica , parámetros de van der Waals que reflejan estimaciones de radio atómico y longitud de enlace de equilibrio , ángulo y diedro; estos se obtienen mediante el ajuste contra cálculos electrónicos detallados (simulaciones químicas cuánticas) o propiedades físicas experimentales como constantes elásticas , parámetros de red y mediciones espectroscópicas .
Debido a la naturaleza no local de las interacciones no enlazadas, implican al menos interacciones débiles entre todas las partículas del sistema. Su cálculo es normalmente el cuello de botella en la velocidad de las simulaciones MD. Para reducir el costo computacional, los campos de fuerza emplean aproximaciones numéricas como radios de corte desplazados, algoritmos de campo de reacción , suma de Ewald de malla de partículas o el más nuevo método partícula-partícula-partícula-malla ( P3M ).
Los campos de fuerza de la química emplean comúnmente disposiciones de enlace preestablecidas (una excepción es la dinámica ab initio ), y por lo tanto no pueden modelar el proceso de ruptura de enlaces químicos y reacciones explícitamente. Por otro lado, muchos de los potenciales utilizados en física, como los basados en el formalismo de orden de enlace, pueden describir varias coordinaciones diferentes de un sistema y ruptura de enlaces. [33] [34] Ejemplos de tales potenciales incluyen el potencial de Brenner [35] para hidrocarburos y sus desarrollos posteriores para los sistemas C-Si-H [36] y COH [37] . El potencial ReaxFF [38] puede considerarse un híbrido completamente reactivo entre los potenciales de orden de enlace y los campos de fuerza de la química.
Las funciones potenciales que representan la energía no enlazada se formulan como una suma de las interacciones entre las partículas del sistema. La opción más simple, empleada en muchos campos de fuerza populares , es el "potencial de par", en el que la energía potencial total se puede calcular a partir de la suma de las contribuciones de energía entre pares de átomos. Por lo tanto, estos campos de fuerza también se denominan "campos de fuerza aditivos". Un ejemplo de un potencial de par de este tipo es el potencial de Lennard-Jones no enlazado (también denominado potencial 6-12), utilizado para calcular las fuerzas de van der Waals.
Otro ejemplo es el modelo de Born (iónico) de la red iónica. El primer término en la siguiente ecuación es la ley de Coulomb para un par de iones, el segundo término es la repulsión de corto alcance explicada por el principio de exclusión de Pauli y el término final es el término de interacción de dispersión. Por lo general, una simulación solo incluye el término dipolar, aunque a veces también se incluye el término cuadrupolar. [39] [40] Cuando n l = 6, este potencial también se llama potencial de Coulomb-Buckingham .
En los potenciales de muchos cuerpos , la energía potencial incluye los efectos de tres o más partículas que interactúan entre sí. [41] En simulaciones con potenciales por pares, también existen interacciones globales en el sistema, pero ocurren solo a través de términos por pares. En los potenciales de muchos cuerpos, la energía potencial no se puede encontrar mediante una suma sobre pares de átomos, ya que estas interacciones se calculan explícitamente como una combinación de términos de orden superior. En la visión estadística, la dependencia entre las variables no se puede expresar en general utilizando solo productos por pares de los grados de libertad. Por ejemplo, el potencial de Tersoff , [42] que se usó originalmente para simular carbono , silicio y germanio , y desde entonces se ha usado para una amplia gama de otros materiales, implica una suma sobre grupos de tres átomos, siendo los ángulos entre los átomos un factor importante en el potencial. Otros ejemplos son el método de átomos incrustados (EAM), [43] el EDIP, [41] y los potenciales de aproximación de segundo momento de enlace fuerte (TBSMA), [44] donde la densidad electrónica de los estados en la región de un átomo se calcula a partir de una suma de contribuciones de los átomos circundantes, y la contribución de energía potencial es entonces una función de esta suma.
Los potenciales semiempíricos hacen uso de la representación matricial de la mecánica cuántica. Sin embargo, los valores de los elementos de la matriz se encuentran a través de fórmulas empíricas que estiman el grado de superposición de orbitales atómicos específicos. Luego, la matriz se diagonaliza para determinar la ocupación de los diferentes orbitales atómicos y se utilizan nuevamente fórmulas empíricas para determinar las contribuciones energéticas de los orbitales.
Existe una amplia variedad de potenciales semiempíricos, denominados potenciales de enlace fuerte , que varían según los átomos que se modelan.
La mayoría de los campos de fuerza clásicos incluyen implícitamente el efecto de la polarizabilidad , por ejemplo, aumentando las cargas parciales obtenidas a partir de cálculos químicos cuánticos. Estas cargas parciales son estacionarias con respecto a la masa del átomo. Pero las simulaciones de dinámica molecular pueden modelar explícitamente la polarizabilidad con la introducción de dipolos inducidos a través de diferentes métodos, como partículas de Drude o cargas fluctuantes. Esto permite una redistribución dinámica de la carga entre átomos que responde al entorno químico local.
Durante muchos años, las simulaciones MD polarizables se han promocionado como la próxima generación. Para líquidos homogéneos como el agua, se ha logrado una mayor precisión mediante la inclusión de la polarizabilidad. [45] [46] [47] También se han logrado algunos resultados prometedores para las proteínas. [48] [49] Sin embargo, todavía no se sabe cuál es la mejor manera de aproximar la polarizabilidad en una simulación. [ cita requerida ] El punto se vuelve más importante cuando una partícula experimenta diferentes entornos durante su trayectoria de simulación, por ejemplo, la translocación de un fármaco a través de una membrana celular. [50]
En la dinámica molecular clásica, una superficie de energía potencial (normalmente el estado fundamental) se representa en el campo de fuerza. Esto es una consecuencia de la aproximación de Born-Oppenheimer . En estados excitados, reacciones químicas o cuando se necesita una representación más precisa, el comportamiento electrónico se puede obtener a partir de los primeros principios utilizando un método mecánico cuántico, como la teoría del funcional de la densidad . Esto se denomina dinámica molecular ab initio (AIMD). Debido al coste de tratar los grados de libertad electrónicos, la carga computacional de estas simulaciones es mucho mayor que la de la dinámica molecular clásica. Por este motivo, la AIMD suele limitarse a sistemas más pequeños y tiempos más cortos.
Los métodos mecánicos y químicos cuánticos ab initio pueden utilizarse para calcular la energía potencial de un sistema sobre la marcha, según sea necesario para las conformaciones en una trayectoria. Este cálculo se realiza normalmente en la vecindad cercana de la coordenada de reacción . Aunque se pueden utilizar varias aproximaciones, estas se basan en consideraciones teóricas, no en ajustes empíricos. Los cálculos ab initio producen una gran cantidad de información que no está disponible con los métodos empíricos, como la densidad de estados electrónicos u otras propiedades electrónicas. Una ventaja significativa de utilizar métodos ab initio es la capacidad de estudiar reacciones que implican la ruptura o formación de enlaces covalentes, que corresponden a múltiples estados electrónicos. Además, los métodos ab initio también permiten recuperar efectos más allá de la aproximación de Born-Oppenheimer utilizando enfoques como la dinámica clásica-cuántica mixta .
Los métodos de mecánica cuántica (QM) son muy potentes, pero son costosos desde el punto de vista computacional, mientras que los métodos de mecánica clásica o molecular (MM) son rápidos pero tienen varias limitaciones (requieren una parametrización extensa; las estimaciones de energía obtenidas no son muy precisas; no se pueden usar para simular reacciones en las que se rompen o forman enlaces covalentes; y tienen una capacidad limitada para proporcionar detalles precisos sobre el entorno químico). Ha surgido una nueva clase de método que combina los puntos positivos de los cálculos de mecánica cuántica (precisión) y de mecánica molecular (velocidad). Estos métodos se denominan métodos mixtos o híbridos de mecánica cuántica y de mecánica molecular (QM/MM híbridos). [51]
La ventaja más importante del método híbrido QM/MM es la velocidad. El costo de hacer dinámica molecular clásica (MM) en el caso más sencillo escala O(n 2 ), donde n es el número de átomos en el sistema. Esto se debe principalmente al término de interacciones electrostáticas (cada partícula interactúa con todas las demás partículas). Sin embargo, el uso del radio de corte, actualizaciones periódicas de listas de pares y, más recientemente, las variaciones del método de Ewald de malla de partículas (PME) han reducido esto a entre O(n) y O(n 2 ). En otras palabras, si se simula un sistema con el doble de átomos, se necesitaría entre dos y cuatro veces más potencia de cálculo. Por otro lado, los cálculos ab initio más simples generalmente escalan O(n 3 ) o peor (se han sugerido cálculos restringidos de Hartree-Fock para escalar ~O(n 2.7 )). Para superar el límite, una pequeña parte del sistema se trata de forma cuántica-mecánica (normalmente el sitio activo de una enzima) y el sistema restante se trata de forma clásica.
En implementaciones más sofisticadas, existen métodos QM/MM para tratar tanto los núcleos ligeros susceptibles a los efectos cuánticos (como los hidrógenos) como los estados electrónicos. Esto permite generar funciones de onda del hidrógeno (similares a las funciones de onda electrónicas). Esta metodología ha sido útil para investigar fenómenos como el efecto túnel del hidrógeno. Un ejemplo en el que los métodos QM/MM han proporcionado nuevos descubrimientos es el cálculo de la transferencia de hidruros en la enzima alcohol deshidrogenasa hepática . En este caso, el efecto túnel cuántico es importante para el hidrógeno, ya que determina la velocidad de reacción. [52]
En el otro extremo de la escala de detalle se encuentran los modelos de grano grueso y de red. En lugar de representar explícitamente cada átomo del sistema, se utilizan "pseudoátomos" para representar grupos de átomos. Las simulaciones MD en sistemas muy grandes pueden requerir recursos informáticos tan grandes que no pueden estudiarse fácilmente con los métodos tradicionales de todos los átomos. De manera similar, las simulaciones de procesos en escalas de tiempo largas (más allá de aproximadamente 1 microsegundo) son prohibitivamente caras, porque requieren muchos pasos de tiempo. En estos casos, a veces se puede abordar el problema utilizando representaciones reducidas, que también se denominan modelos de grano grueso . [53]
Ejemplos de métodos de granulado grueso (CG) son la dinámica molecular discontinua (CG-DMD) [54] [55] y los modelos Go. [56] El granulado grueso se realiza a veces tomando pseudoátomos más grandes. Estas aproximaciones de átomos unidos se han utilizado en simulaciones MD de membranas biológicas. La implementación de dicho enfoque en sistemas donde las propiedades eléctricas son de interés puede ser un desafío debido a la dificultad de utilizar una distribución de carga adecuada en los pseudoátomos. [57] Las colas alifáticas de los lípidos están representadas por unos pocos pseudoátomos al reunir de 2 a 4 grupos metileno en cada pseudoátomo.
La parametrización de estos modelos de grano grueso debe realizarse de manera empírica, haciendo coincidir el comportamiento del modelo con datos experimentales apropiados o simulaciones de todos los átomos. Lo ideal sería que estos parámetros tuvieran en cuenta tanto las contribuciones entálpicas como entrópicas a la energía libre de manera implícita. [58] Cuando el granulado grueso se realiza en niveles superiores, la precisión de la descripción dinámica puede ser menos confiable. Pero los modelos de grano grueso se han utilizado con éxito para examinar una amplia gama de cuestiones en biología estructural, organización de cristales líquidos y vidrios poliméricos.
Ejemplos de aplicaciones del granulado grueso:
La forma más simple de granulado grueso es el átomo unido (a veces llamado átomo extendido ) y se utilizó en la mayoría de las primeras simulaciones MD de proteínas, lípidos y ácidos nucleicos. Por ejemplo, en lugar de tratar explícitamente los cuatro átomos de un grupo metilo CH3 (o los tres átomos del grupo metileno CH2 ) , se representa todo el grupo con un pseudoátomo. Por supuesto, debe estar parametrizado correctamente para que sus interacciones de van der Waals con otros grupos tengan la dependencia de la distancia adecuada. Se aplican consideraciones similares a los enlaces, ángulos y torsiones en los que participa el pseudoátomo. En este tipo de representación de átomos unidos, normalmente se eliminan todos los átomos de hidrógeno explícitos excepto aquellos que tienen la capacidad de participar en enlaces de hidrógeno ( hidrógenos polares ). Un ejemplo de esto es el campo de fuerza CHARMM 19.
Los hidrógenos polares suelen conservarse en el modelo, porque el tratamiento adecuado de los enlaces de hidrógeno requiere una descripción razonablemente precisa de la direccionalidad y las interacciones electrostáticas entre los grupos donadores y aceptores. Un grupo hidroxilo, por ejemplo, puede ser tanto un donador de enlaces de hidrógeno como un aceptor de enlaces de hidrógeno, y sería imposible tratar esto con un pseudoátomo de OH. Aproximadamente la mitad de los átomos de una proteína o un ácido nucleico son hidrógenos no polares, por lo que el uso de átomos unidos puede proporcionar un ahorro sustancial en tiempo de computación.
Los campos de fuerza de aprendizaje automático (MLFF) representan un enfoque para modelar interacciones interatómicas en simulaciones de dinámica molecular. [59] Los MLFF pueden lograr una precisión cercana a la de los métodos ab initio . Una vez entrenados, los MLFF son mucho más rápidos que los cálculos mecánicos cuánticos directos. Los MLFF abordan las limitaciones de los campos de fuerza tradicionales al aprender superficies de energía potencial complejas directamente a partir de datos mecánicos cuánticos de alto nivel. Varios paquetes de software ahora admiten MLFF, incluidos VASP [60] y bibliotecas de código abierto como DeePMD-kit [61] [62] y SchNetPack. [63] [64]
En muchas simulaciones de un sistema soluto-disolvente, el enfoque principal se centra en el comportamiento del soluto con poco interés en el comportamiento del disolvente, particularmente en aquellas moléculas de disolvente que residen en regiones alejadas de la molécula de soluto. [65] Los disolventes pueden influir en el comportamiento dinámico de los solutos a través de colisiones aleatorias e imponiendo un arrastre de fricción en el movimiento del soluto a través del disolvente. El uso de condiciones de contorno periódicas no rectangulares, límites estocásticos y capas de disolvente pueden ayudar a reducir la cantidad de moléculas de disolvente necesarias y permitir que una mayor proporción del tiempo de cálculo se dedique a simular el soluto. También es posible incorporar los efectos de un disolvente sin necesidad de que haya moléculas de disolvente explícitas presentes. Un ejemplo de este enfoque es utilizar una fuerza media potencial (PMF) que describe cómo cambia la energía libre a medida que varía una coordenada particular. El cambio de energía libre descrito por PMF contiene los efectos promediados del disolvente.
Sin incorporar los efectos de solventes, las simulaciones de macromoléculas (como las proteínas) pueden producir un comportamiento poco realista e incluso las moléculas pequeñas pueden adoptar conformaciones más compactas debido a fuerzas de van der Waals favorables e interacciones electrostáticas que se verían atenuadas en presencia de un solvente. [66]
Una interacción de largo alcance es una interacción en la que la interacción espacial disminuye no más rápido que donde es la dimensionalidad del sistema. Algunos ejemplos incluyen interacciones carga-carga entre iones e interacciones dipolo-dipolo entre moléculas. Modelar estas fuerzas presenta un gran desafío ya que son significativas en una distancia que puede ser mayor que la mitad de la longitud de la caja con simulaciones de muchos miles de partículas. Aunque una solución sería aumentar significativamente el tamaño de la longitud de la caja, este enfoque de fuerza bruta es menos que ideal ya que la simulación se volvería computacionalmente muy costosa. Truncar esféricamente el potencial también está fuera de cuestión ya que se puede observar un comportamiento poco realista cuando la distancia está cerca de la distancia de corte. [67]
Las simulaciones de dinámica molecular dirigida (SMD), o simulaciones de sonda de fuerza, aplican fuerzas a una proteína para manipular su estructura arrastrándola a lo largo de los grados de libertad deseados. Estos experimentos se pueden utilizar para revelar cambios estructurales en una proteína a nivel atómico. La SMD se utiliza a menudo para simular eventos como el desdoblamiento o estiramiento mecánico. [68]
Existen dos protocolos típicos de SMD: uno en el que la velocidad de tracción se mantiene constante y otro en el que la fuerza aplicada es constante. Normalmente, una parte del sistema estudiado (por ejemplo, un átomo en una proteína) está restringida por un potencial armónico. Luego, se aplican fuerzas a átomos específicos a una velocidad constante o a una fuerza constante. El muestreo de paraguas se utiliza para mover el sistema a lo largo de la coordenada de reacción deseada variando, por ejemplo, las fuerzas, distancias y ángulos manipulados en la simulación. A través del muestreo de paraguas, se muestrean adecuadamente todas las configuraciones del sistema, tanto de alta energía como de baja energía. Luego, el cambio de energía libre de cada configuración se puede calcular como el potencial de fuerza media . [69] Un método popular para calcular PMF es a través del método de análisis de histograma ponderado (WHAM), que analiza una serie de simulaciones de muestreo de paraguas. [70] [71]
Muchas aplicaciones importantes de SMD se encuentran en el campo del descubrimiento de fármacos y las ciencias biomoleculares. Por ejemplo, SMD se utilizó para investigar la estabilidad de las protofibrillas de Alzheimer, [72] para estudiar la interacción proteína-ligando en la quinasa dependiente de ciclina 5 [73] e incluso para mostrar el efecto del campo eléctrico en el complejo de trombina (proteína) y aptámero (nucleótido) [74] entre muchos otros estudios interesantes.
La dinámica molecular se utiliza en muchos campos de la ciencia.
Los siguientes ejemplos biofísicos ilustran esfuerzos notables para producir simulaciones de sistemas de tamaño muy grande (un virus completo) o tiempos de simulación muy largos (hasta 1,112 milisegundos):
Otra aplicación importante del método MD se beneficia de su capacidad de caracterización tridimensional y análisis de la evolución microestructural a escala atómica.
El modelado molecular en GPU es la técnica de utilizar una unidad de procesamiento gráfico (GPU) para simulaciones moleculares. [87]
En 2007, NVIDIA presentó tarjetas de video que podían usarse no solo para mostrar gráficos sino también para cálculos científicos. Estas tarjetas incluyen muchas unidades aritméticas (a partir de 2016 [actualizar], hasta 3.584 en Tesla P100) trabajando en paralelo. Mucho antes de este evento, la potencia computacional de las tarjetas de video se usaba exclusivamente para acelerar los cálculos gráficos. Lo novedoso es que NVIDIA hizo posible desarrollar programas paralelos en una interfaz de programación de aplicaciones (API) de alto nivel llamada CUDA . Esta tecnología simplificó sustancialmente la programación al permitir escribir programas en C / C++ . Más recientemente, OpenCL permite la aceleración de GPU multiplataforma .