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, dando una visión de la "evolución" dinámica del sistema. En la versión más común, las trayectorias de átomos y moléculas se determinan resolviendo numéricamente las ecuaciones de movimiento de Newton para un sistema de partículas que interactúan, 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 materiales y biofísica .
Debido a que los sistemas moleculares normalmente constan de una gran cantidad de partículas, es imposible determinar analíticamente las propiedades de sistemas tan complejos ; La simulación MD evita este problema mediante el uso de métodos numéricos . Sin embargo, las simulaciones MD largas están matemáticamente mal condicionadas , generando errores acumulativos en la integración numérica que pueden minimizarse con la selección adecuada de algoritmos y parámetros, pero no eliminarse.
Para sistemas que obedecen a la hipótesis ergódica , la evolución de una simulación de dinámica molecular se puede utilizar para determinar las propiedades termodinámicas macroscópicas del sistema: los promedios temporales de un sistema ergódico corresponden a promedios de conjuntos microcanónicos . MD también ha sido denominada " mecánica estadística por números" y " 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.
MD se desarrolló originalmente a principios de la década de 1950, tras los éxitos anteriores con las simulaciones de Monte Carlo , que a su vez se remontan al siglo XVIII, en el problema de la aguja de Buffon, por ejemplo, pero Rosenbluth y Metropolis lo popularizaron para la mecánica estadística en el Laboratorio Nacional de Los Álamos. en lo que hoy se conoce como algoritmo de Metropolis-Hastings . El interés por la evolución temporal de los sistemas de N cuerpos data mucho antes del siglo XVII, comenzando con Newton, y continuó hasta el siglo siguiente centrándose en gran medida en la mecánica celeste y cuestiones como la estabilidad del sistema solar. Muchos de los métodos numéricos que se utilizan hoy en día se desarrollaron durante este período, anterior al uso de las 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 "a mano".
Ya en 1941 la integración de las ecuaciones de movimiento de muchos cuerpos se llevó a cabo con ordenadores analógicos. Algunos 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 disponerlos de tal manera que replicaran la estructura de un líquido y utilizarlo para examinar su comportamiento. JD Bernal dijo, en 1962: [2]
... Tomé varias pelotas de goma y las uní con varillas de una selección de diferentes longitudes que iban desde 2,75 a 4 pulgadas. En primer lugar, traté de hacer esto de la manera más informal posible, trabajando en mi propia oficina, siendo interrumpido aproximadamente cada cinco minutos 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, Fermi propuso en 1953, y publicó en 1955, [3] el uso de MANIAC I , también en el Laboratorio Nacional de Los Álamos , para resolver la evolución temporal de las ecuaciones de movimiento para muchos sistema corporal sujeto a varias elecciones de leyes de fuerza; hoy, este trabajo fundamental se conoce como el problema de Fermi-Pasta-Ulam-Tsingou . La evolución temporal de la energía del trabajo original se muestra en la figura de la derecha.
En 1957, Alder y Wainwright [4] 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, Gibson et al. Daño por radiación simulado de cobre sólido mediante el uso de una interacción repulsiva del tipo Born-Mayer junto con una fuerza superficial cohesiva. [5] En 1964, Rahman [6] publicó simulaciones de argón líquido que utilizaban 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 elemento básico en muchos campos de fuerza de sustancias reales. [12] [13]
Utilizado por primera vez en física teórica , el método MD ganó popularidad en la ciencia de materiales poco después y, desde la década de 1970, también es común en bioquímica y biofísica . La MD se utiliza con frecuencia para refinar estructuras tridimensionales de proteínas y otras macromoléculas basándose en limitaciones experimentales de la cristalografía de rayos X o la 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 subplantación de iones, y también 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 , que pueden ser útiles 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, la MD se puede utilizar para la predicción ab initio de la estructura de la proteína simulando el plegamiento de la cadena polipeptídica a partir de una espiral aleatoria .
Los resultados de las simulaciones de MD se pueden probar comparándolos con experimentos que miden la dinámica molecular, de los cuales un método popular es la espectroscopia de RMN. Las predicciones de estructuras derivadas de MD se pueden probar mediante experimentos en toda la comunidad en Evaluación crítica de predicción de estructuras de proteínas ( CASP ), aunque históricamente el método ha tenido un éxito limitado en esta área. Michael Levitt , que compartió el Premio Nobel en parte por la aplicación de la MD a las proteínas, escribió en 1999 que los participantes del CASP normalmente no utilizaban el método debido a "... una vergüenza central de la mecánica molecular, a saber, que la minimización de energía o la dinámica molecular en general conduce a un modelo que se parece menos a la estructura experimental". [14] Las mejoras en los recursos computacionales que permiten trayectorias MD más largas y más largas, 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 alcanzar el 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 de 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] Por otro lado, Carlson et al. Implementaron simulación de dinámica molecular para identificar compuestos que complementan el receptor y al mismo tiempo causan una alteración mínima de 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 fotogramas) para el desarrollo del farmacóforo. Spyrakis et al. se basó en un flujo de trabajo de simulaciones MD, huellas dactilares para ligandos y proteínas (FLAP) y análisis discriminante lineal 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 por mejorar el modelado de descubrimiento de fármacos basado en estructuras, frente a 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. Luego, los contactos críticos se pueden convertir en modelos farmacóforos que se pueden utilizar para la detección virtual. [20]
Un factor importante son los enlaces de hidrógeno intramoleculares , [21] que no están incluidos explícitamente en los campos de fuerza modernos, sino que se describen como interacciones de Coulomb de cargas puntuales atómicas. [ cita necesaria ] Esta es una aproximación burda porque los enlaces de hidrógeno tienen una naturaleza química y mecánica parcialmente cuántica . Además, las interacciones electrostáticas suelen calcularse utilizando la constante dieléctrica del vacío , aunque la solución acuosa circundante tiene una constante dieléctrica mucho mayor. 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 basados en la teoría de Fritz London que solo es aplicable en el vacío. [ cita necesaria ] 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 medio ambiente . [22] La medición directa de las fuerzas de atracción entre diferentes materiales (como la constante de Hamaker ) muestra que "la interacción entre los hidrocarburos a través del agua es aproximadamente el 10% de la interacción a través del vacío". [22] La dependencia del entorno de las fuerzas de Van der Waals se ignora 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. Se deben seleccionar el tamaño de la simulación ( n = número de partículas), el paso de tiempo y la duración total para que el cálculo pueda finalizar dentro de un período de tiempo razonable. Sin embargo, las simulaciones deberían ser lo suficientemente largas como para ser relevantes para las escalas de tiempo de los procesos naturales que se estudian. 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 ser humano cuando sólo mira menos de un paso. La mayoría de las publicaciones científicas sobre la dinámica de las proteínas y el ADN [23] [24] utilizan datos de simulaciones que abarcan nanosegundos (10 −9 s) hasta microsegundos (10 −6 s). Para obtener estas simulaciones, se necesitan desde varios días de CPU hasta 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. [25]
Durante una simulación MD clásica, la tarea que consume más CPU es la evaluación del potencial en función de las coordenadas internas de las partículas. Dentro de esa evaluación energética, la más cara es la parte no enlazada o no covalente. En la notación Big O , las simulaciones de dinámica molecular comunes se escalan si todas las interacciones electrostáticas y de van der Waals por pares deben tenerse en cuenta explícitamente. Este costo computacional se puede reducir empleando métodos electrostáticos como la suma de Ewald de malla de partículas ( ), malla de partículas-partícula-partícula ( P3M ) o buenos métodos de corte esférico ( ). [ cita necesaria ]
Otro factor que afecta el tiempo total de CPU que necesita una simulación es el tamaño del paso de tiempo de integración. Este es el período de tiempo entre evaluaciones del potencial. El paso de tiempo debe elegirse lo suficientemente pequeño para evitar errores de discretización (es decir, más pequeño que el período relacionado con la frecuencia vibratoria más rápida del sistema). Los pasos de tiempo típicos para la MD clásica son del orden de 1 femtosegundo (10 −15 s). Este valor puede ampliarse 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últiples métodos de escala de tiempo, que permiten tiempos prolongados entre actualizaciones de fuerzas más lentas de largo alcance. [26] [27] [28]
Para simular moléculas en un disolvente, se debe elegir entre disolvente explícito e implícito . Las partículas de disolvente explícitas (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 y requiere la inclusión de aproximadamente diez veces más partículas en la simulación. Pero la granularidad y viscosidad del disolvente explícito es esencial para reproducir ciertas propiedades de las moléculas del 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 para evitar artefactos de condiciones límite . Las condiciones de contorno a menudo se tratan eligiendo valores fijos en los bordes (lo que puede causar artefactos), o empleando condiciones de contorno periódicas en las que un lado de la simulación regresa al lado opuesto, imitando una fase masiva (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, conservándose la energía total. Para un sistema de N partículas con coordenadas y velocidades , el siguiente par de ecuaciones diferenciales de primer orden se pueden escribir en la notación de Newton como
La función de energía potencial del sistema es función de las coordenadas de las partículas . Se le conoce simplemente como potencial en física o campo de fuerzas en química. La primera ecuación proviene de las leyes del 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 integrador simpléctico como la integración de Verlet . La evolución temporal de y se llama 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 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 encuentra igualando 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 a la pequeña cantidad 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 serían rápidamente transportados y compartidos entre un gran número de átomos ( o más) sin grandes cambios de temperatura. Sin embargo, cuando sólo hay 500 átomos, el sustrato se vaporiza casi inmediatamente por la deposición. Algo similar ocurre en las simulaciones biofísicas. La temperatura del sistema en NVE aumenta naturalmente cuando macromoléculas como las proteínas experimentan cambios conformacionales y unión exotérmica.
En el conjunto canónico se conservan la cantidad de sustancia (N), el volumen (V) y la temperatura (T). A veces también se le llama 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 disponibles una variedad de algoritmos de termostato para agregar y eliminar energía de los límites de una simulación MD de una manera más o menos realista, acercándose al conjunto canónico . Los métodos populares para controlar la temperatura incluyen el cambio de escala de velocidad, el termostato Nosé-Hoover , las cadenas Nosé-Hoover, el termostato Berendsen , el termostato Andersen y la dinámica Langevin . El termostato Berendsen podría introducir el efecto del cubito de hielo volador , que provoca traslaciones y rotaciones no físicas del sistema simulado.
No es trivial obtener una distribución de conjunto canónico de conformaciones y velocidades utilizando estos algoritmos. Cómo 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 este 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 del termostato, se necesita un barostato. Corresponde más estrechamente a 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 apropiado. Para las bicapas lipídicas, el control de la presión se produce bajo un área de membrana constante (NPAT) o 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 desordenados. También se le llama templado paralelo. La formulación MD de intercambio de réplicas (REMD) [29] intenta superar el problema de mínimos múltiples intercambiando la temperatura de réplicas que no interactúan del sistema que se ejecuta 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 suele denominarse campo de fuerzas y en física de materiales potencial interatómico . Los potenciales pueden definirse en muchos niveles de precisión física; los más comúnmente utilizados en química se basan en la mecánica molecular y representan un tratamiento mecánico clásico de las interacciones entre partículas que pueden reproducir cambios estructurales y conformacionales , pero generalmente no pueden reproducir reacciones químicas .
La reducción de una descripción totalmente 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, pueden ser tratados por separado. El segundo 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 de detalle más finos, 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 forma clásica pero una pequeña región se trata como un sistema cuántico, que generalmente sufre 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 . [30] Los potenciales empíricos representan efectos mecánico-cuánticos de forma limitada a través de aproximaciones funcionales ad hoc. Estos potenciales contienen parámetros libres como la carga atómica , los parámetros de van der Waals que reflejan estimaciones del radio atómico y la longitud , el ángulo y el diédrico del enlace de equilibrio; estos se obtienen ajustándolos con 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 suele ser 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 la nueva malla partícula-partícula-partícula ( P3M ).
Los campos de fuerza químicos comúnmente emplean disposiciones de enlace preestablecidas (una excepción es la dinámica ab initio ) y, por lo tanto, no pueden modelar explícitamente el proceso de ruptura de enlaces químicos y reacciones. Por otro lado, muchos de los potenciales utilizados en física, como los basados en el formalismo del orden de los enlaces, pueden describir varias coordinaciones diferentes de un sistema y la ruptura de enlaces. [31] [32] Ejemplos de tales potenciales incluyen el potencial de Brenner [33] para hidrocarburos y sus desarrollos posteriores para los sistemas C-Si-H [34] y COH [35] . El potencial ReaxFF [36] puede considerarse un híbrido totalmente reactivo entre potenciales de orden de enlace y campos de fuerza químicos.
Las funciones potenciales que representan la energía no enlazada se formulan como una suma de interacciones entre las partículas del sistema. La opción más sencilla, 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 tal par de potenciales 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 Born (iónico) de red iónica. El primer término de 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 último término es el término de interacción de dispersión. Normalmente, una simulación sólo incluye el término dipolar, aunque a veces también se incluye el término cuadrupolar. [37] [38] 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í. [39] 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 la suma de pares de átomos, ya que estas interacciones se calculan explícitamente como una combinación de términos de orden superior. Desde el punto de vista estadístico, la dependencia entre las variables en general no se puede expresar utilizando únicamente productos por pares de los grados de libertad. Por ejemplo, el potencial de Tersoff , [40] que se utilizó originalmente para simular carbono , silicio y germanio , y desde entonces se ha utilizado para una amplia gama de otros materiales, implica una suma de grupos de tres átomos, con los ángulos entre los Los átomos son un factor importante en el potencial. Otros ejemplos son el método del átomo integrado (EAM), [41] el EDIP, [39] y los potenciales de aproximación del segundo momento de enlace estrecho (TBSMA), [42] donde la densidad electrónica de los estados en la región de un átomo es 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 utilizan la representación matricial de la mecánica cuántica. Sin embargo, los valores de los elementos de la matriz se encuentran mediante fórmulas empíricas que estiman el grado de superposición de orbitales atómicos específicos. Luego se diagonaliza la matriz para determinar la ocupación de los diferentes orbitales atómicos, y se utilizan nuevamente fórmulas empíricas para determinar las contribuciones de energía de los orbitales.
Existe una amplia variedad de potenciales semiempíricos, denominados potenciales de unión estrecha , 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 mediante diferentes métodos, como partículas Drude o cargas fluctuantes. Esto permite una redistribución dinámica de 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 polarizabilidad. [43] [44] [45] También se han logrado algunos resultados prometedores para las proteínas. [46] [47] Sin embargo, todavía no está claro cómo aproximar mejor la polarizabilidad en una simulación. [ cita necesaria ] 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. [48]
En la dinámica molecular clásica, una superficie de energía potencial (normalmente el estado fundamental) está representada 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 primeros principios utilizando un método de la mecánica cuántica, como la teoría funcional de la densidad . Esto se llama Dinámica Molecular Ab Initio (AIMD). Debido al costo 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, AIMD suele limitarse a sistemas más pequeños y tiempos más cortos.
Se pueden utilizar métodos químicos y mecánicos cuánticos ab initio 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 generalmente se realiza en las inmediaciones de la coordenada de reacción . Aunque se pueden utilizar varias aproximaciones, éstas 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 mediante 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 mixta cuántica-clásica .
Los métodos QM (mecánicos cuánticos) son muy poderosos. Sin embargo, son costosos desde el punto de vista computacional, mientras que los métodos MM (mecánica clásica o molecular) son rápidos pero adolecen de varias limitaciones (requieren una parametrización extensa; las estimaciones de energía obtenidas no son muy precisas; no se pueden utilizar para simular reacciones en las que se rompen/forman enlaces covalentes). ; y tienen capacidades limitadas 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 QM (precisión) y MM (velocidad). Estos métodos se denominan métodos mixtos o híbridos de mecánica cuántica y mecánica molecular (híbridos QM/MM). [49]
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, las actualizaciones periódicas de la lista de pares y, más recientemente, las variaciones del método de malla de partículas de Ewald (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 ha sugerido que los cálculos restringidos de Hartree-Fock escalan ~O(n 2,7 )). Para superar el límite, una pequeña parte del sistema se trata de forma mecánica cuántica (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 núcleos ligeros susceptibles a efectos cuánticos (como los hidrógenos) como estados electrónicos. Esto permite generar funciones de onda de hidrógeno (similares a las funciones de onda electrónicas). Esta metodología ha resultado útil en la investigación de fenómenos como los túneles de hidrógeno. Un ejemplo en el que los métodos QM/MM han proporcionado nuevos descubrimientos es el cálculo de la transferencia de hidruro en la enzima alcohol deshidrogenasa hepática . En este caso, el túnel cuántico es importante para el hidrógeno, ya que determina la velocidad de reacción. [50]
En el otro extremo de la escala de detalle se encuentran los modelos de grano grueso y de celosía. 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 mediante métodos tradicionales de todos los átomos. De manera similar, las simulaciones de procesos en escalas de tiempo largas (más allá de 1 microsegundo aproximadamente) son prohibitivamente costosas 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 . [51]
Ejemplos de métodos de grano grueso (CG) son la dinámica molecular discontinua (CG-DMD) [52] [53] y los modelos Go. [54] El grano 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 este enfoque en sistemas donde las propiedades eléctricas son de interés puede resultar desafiante debido a la dificultad de utilizar una distribución de carga adecuada en los pseudoátomos. [55] Las colas alifáticas de los lípidos están representadas por unos pocos pseudoátomos reuniendo de 2 a 4 grupos metileno en cada pseudoátomo.
La parametrización de estos modelos de grano muy grueso debe realizarse empíricamente, haciendo coincidir el comportamiento del modelo con datos experimentales apropiados o simulaciones de todos los átomos. Idealmente, estos parámetros deberían tener en cuenta las contribuciones tanto entálpicas como entrópicas a la energía libre de forma implícita. [56] Cuando el análisis de grano grueso se realiza en niveles más altos, la precisión de la descripción dinámica puede ser menos confiable. Pero se han utilizado con éxito modelos de grano muy grueso para examinar una amplia gama de cuestiones en biología estructural, organización de cristales líquidos y vidrios poliméricos.
Ejemplos de aplicaciones de grano grueso:
La forma más simple de grano 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 CH 3 (o los tres átomos del grupo metileno CH 2 ), se representa el grupo completo con un pseudoátomo. Por supuesto, debe parametrizarse adecuadamente para que sus interacciones de Van der Waals con otros grupos tengan la dependencia adecuada de la distancia. 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 generalmente se retienen 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 donador y aceptor. Un grupo hidroxilo, por ejemplo, puede ser tanto un donante como un aceptor de enlaces de hidrógeno, y sería imposible tratarlo con un pseudoátomo de OH. Aproximadamente la mitad de los átomos de una proteína o ácido nucleico son hidrógenos no polares, por lo que el uso de átomos unidos puede suponer un ahorro sustancial de tiempo en el ordenador.
En muchas simulaciones de un sistema soluto-solvente, el enfoque principal está en el comportamiento del soluto con poco interés en el comportamiento del solvente, particularmente en aquellas moléculas de solvente que residen en regiones alejadas de la molécula del soluto. [57] Los disolventes pueden influir en el comportamiento dinámico de los solutos mediante colisiones aleatorias e imponiendo un arrastre de fricción sobre 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 solvente pueden ayudar a reducir la cantidad de moléculas de solvente requeridas y permitir que se dedique una mayor proporción del tiempo de cálculo a simular el soluto. También es posible incorporar los efectos de un disolvente sin necesidad de que esté presente ninguna molécula de disolvente explícita. 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 solvente.
Sin incorporar los efectos de las simulaciones de macromoléculas (como las proteínas) con disolventes, se puede producir un comportamiento poco realista e incluso las moléculas pequeñas pueden adoptar conformaciones más compactas debido a las fuerzas favorables de Van der Waals y a las interacciones electrostáticas que se amortiguarían en presencia de un disolvente. [58]
Una interacción de largo alcance es una interacción en la que la interacción espacial no disminuye más rápido que la dimensionalidad del sistema. Los ejemplos incluyen interacciones carga-carga entre iones e interacciones dipolo-dipolo entre moléculas. Modelar estas fuerzas presenta todo un 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 no es ideal ya que la simulación resultaría computacionalmente muy costosa. Tampoco es posible truncar el potencial de forma esférica, ya que se puede observar un comportamiento poco realista cuando la distancia es cercana a la distancia de corte. [59]
Las simulaciones de dinámica molecular dirigida (SMD), o simulaciones de sondas de fuerza, aplican fuerzas a una proteína para manipular su estructura tirando de ella 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. SMD se utiliza a menudo para simular eventos como el despliegue o estiramiento mecánico. [60]
Hay 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, parte del sistema estudiado (p. ej., un átomo de una proteína) está restringido por un potencial armónico. Luego se aplican fuerzas a átomos específicos a una velocidad constante o a una fuerza constante. El muestreo general 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 general, se muestrean adecuadamente todas las configuraciones del sistema, tanto de alta como de baja energía. Luego, el cambio de energía libre de cada configuración se puede calcular como el potencial de la fuerza media . [61] 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 generales. [62] [63]
Muchas aplicaciones importantes de SMD se encuentran en el campo del descubrimiento de fármacos y las ciencias biomoleculares. Por ejemplo, se utilizó SMD para investigar la estabilidad de las protofibrillas de Alzheimer, [64] para estudiar la interacción del ligando de proteína en la quinasa 5 dependiente de ciclina [65] e incluso para mostrar el efecto del campo eléctrico sobre la trombina (proteína) y el aptámero (nucleótido). complejo [66] 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 muy gran tamaño (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 de gráficos (GPU) para simulaciones moleculares. [79]
En 2007, NVIDIA introdujo tarjetas de vídeo que podían usarse no sólo para mostrar gráficos sino también para cálculos científicos. Estas tarjetas incluyen muchas unidades aritméticas (a partir de 2016 [update], hasta 3584 en el Tesla P100) que funcionan en paralelo. Mucho antes de este evento, el poder computacional de las tarjetas de video se usaba exclusivamente para acelerar los cálculos gráficos. Lo nuevo 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 que los programas se escribieran en C / C++ . Más recientemente, OpenCL permite la aceleración de GPU multiplataforma .