Repositorio Institucional de Documentos

Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es Proyecto Fin de Carrera Implementación e integración en GROMACS de un algoritmo eficiente y preciso para imponer ligaduras en simulaciones de dinámica molecular Autora: María Astón Serrano Gracia Directores: Pablo García Risueño Institut für Physik und IRIS Adlershof Humboldt Universität zu Berlin Jesús Alastruey Benedé Dpto. Informática e Ingeniería de Sistemas Universidad de Zaragoza Ingeniería informática Escuela de Ingeniería y Arquitectura Universidad de Zaragoza Curso 2012/2013 Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia Implementación e integración en GROMACS de un algoritmo eficiente y preciso para imponer ligaduras en simulaciones de Dinámica Molecular RESUMEN Las simulaciones de dinámica molecular describen la evolución en el tiempo de un sistema de partículas. Herramientas computacionales para realizar dichas simulaciones son fundamentales para los investigadores de campos tan diversos como la medicina (búsqueda de tratamientos para enfermedades como el Alzheimer, la fibrosis quística o el cáncer; diseño computacional de medicamentos), la química (diseño de catalizadores) o la ingeniería de materiales. GROMACS es un extendido y versátil paquete de dinámica molecular escrito en el lenguaje de programación C. Originalmente desarrollado en la Universidad de Groningen (Países Bajos), actualmente se mantiene por desarrolladores de universidades y centros de investigación de todo el mundo, y cuenta con miles de usuarios. Para realizar simulaciones eficientes, una opción común es imponer ligaduras, es decir, restricciones a la longitud de los enlaces atómicos o al valor de los distintos ángulos entre átomos. Esto permite incrementar el paso temporal (time step) de las simulaciones y con ello lograr realizar simulaciones más duraderas, con lo que se incrementa el poder predictivo y el realismo de la simulación. Los dos algoritmos más utilizados para la implementación de ligaduras, los denominados SHAKE y LINCS (LINear Constraint Solver), presentan limitaciones a la hora de imponer ligaduras en los ángulos. Además, están basados en aproximaciones, lo cual afecta negativamente a su exactitud, eficiencia y estabilidad numérica. El presente proyecto fin de carrera se basa en la reciente propuesta teórica del Dr. Pablo García Risueño, la cual propone mejorar la imposición de ligaduras en las simulaciones de dinámica molecular mediante el uso de cálculos analíticos. Por lo tanto, el objetivo de este proyecto es implementar el algoritmo denominado ILVES-S, propuesta alternativa al algoritmo SHAKE, e integrarlo en el paquete de dinámica molecular GROMACS. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia AGRADECIMIENTOS Con este proyecto se cierra una importante etapa y muchas son las personas que han pasado por ella, quisiera que todas se sintieran agradecidas pero algunas se merecen una especial mención. Primero, quiero agradecer el apoyo y la inspiración que he recibido de las personas con las que he trabajado en este proyecto y con las que tanto he aprendido: el Dr. Pablo García Risueño, el Dr. Jesús Alastruey Benedé y el Dr. Carl Christian Kjelgaard Mikkelsen. Gracias a todos mis amigos por los buenos momentos inolvidables. Víctor, gracias por estar ahí siempre y ser mi gran apoyo. Finalmente, gracias a mi familia por todo el cariño que recibo de todos ellos pero muy especialmente gracias a mis padres por su lucha y sacrificio continuo para que mis hermanos y yo tengamos un futuro mejor. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 1 Índice general 1. INTRODUCCIÓN .................................................................................................................................................. 3 1.1 CONTEXTO DEL PROYECTO Y MOTIVACIÓN .................................................................................................................. 3 1.2 OBJETIVOS ........................................................................................................................................................... 4 1.3 PLANIFICACIÓN Y TRABAJO REALIZADO ....................................................................................................................... 4 1.4 ORGANIZACIÓN DEL DOCUMENTO ............................................................................................................................ 6 2. DINÁMICA MOLECULAR ...................................................................................................................................... 7 2.1 SIMULACIONES DE DINÁMICA MOLECULAR .................................................................................................................. 7 2.2 PROPUESTA TEÓRICA: ILVES-S .............................................................................................................................. 13 3. GROMACS ......................................................................................................................................................... 15 3.1 INTRODUCCIÓN ................................................................................................................................................... 15 3.2 ANÁLISIS DEL ALGORITMO SHAKE EN GROMACS. ................................................................................................... 16 4. IMPLEMENTACIÓN DEL ALGORITMO ILVES-S .................................................................................................... 21 4.1 INTRODUCCIÓN ................................................................................................................................................... 21 4.2 INICIALIZACIÓN ................................................................................................................................................... 22 4.3 PROCESO ITERATIVO ............................................................................................................................................ 27 5. RESULTADOS .................................................................................................................................................... 29 5.1 ENTORNO EXPERIMENTAL ..................................................................................................................................... 29 5.2 VERIFICACIÓN DE LOS RESULTADOS ......................................................................................................................... 30 5.3 ANÁLISIS DE LOS RESULTADOS ................................................................................................................................ 30 6. CONCLUSIONES ................................................................................................................................................ 35 6.1 CONCLUSIONES ................................................................................................................................................... 35 6.2 TRABAJO FUTURO ................................................................................................................................................ 35 6.3 VALORACIÓN PERSONAL ....................................................................................................................................... 36 ANEXO I: CAMPOS DE FUERZA .............................................................................................................................. 37 ANEXO II: INSTALACIÓN Y EJECUCIÓN DE GROMACS ............................................................................................ 41 ANEXO III: ANÁLISIS DE LAS MATRICES DISPERSAS ............................................................................................... 55 ANEXO IV: TRABAJO DESARROLLADO EN EL ICHEC (IRLANDA) .............................................................................. 59 REFERENCIAS ........................................................................................................................................................ 63 Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 2 Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 3 Capítulo 1 INTRODUCCIÓN 1.1 Contexto del proyecto y motivación La dinámica molecular es una técnica de simulación por ordenador en la que átomos y moléculas interactúan durante un período de tiempo, dando lugar a una descripción visual y detallada del movimiento de las partículas. Estas simulaciones se basan en aproximaciones de la física conocida y se utilizan con frecuencia en el estudio de muchas moléculas bioquímicas como proteínas, lípidos y ácidos nucleicos. La dinámica molecular permite a los científicos escudriñar en el movimiento de los átomos de un modo que no es viable en experimentos de laboratorio: es posible tomar instantáneas de las estructuras permitiendo el acceso a todas las escalas de tiempo y movimiento con resolución atómica. Debido a la complejidad de los sistemas biológicos, estas simulaciones conllevan un alto consumo de recursos, tanto de memoria como de capacidad de procesamiento. Es por tanto habitual la utilización de ordenadores de gran rendimiento y con grandes capacidades para la computación paralela o distribuida. Gracias a los avances tecnológicos esta es una ciencia en constante cambio y sobre la que se buscan mejoras en el rendimiento y la precisión de los resultados de las simulaciones. Una técnica muy común que permite realizar simulaciones más eficientes es aumentar el paso temporal (time step) cuyos valores típicos son del orden de un femtosegundo (10−15 s). Este valor puede ser ampliado si se usan algoritmos como SHAKE [1], que permiten imponer ligaduras para fijar las vibraciones de los átomos más rápidos. El presente proyecto pretende conseguir mejorar las simulaciones de dinámica molecular mediante el estudio del algoritmo SHAKE y la posterior implementación de un algoritmo alternativo, denominado ILVES-S, más eficiente y estable. Este proyecto surge como fruto de la colaboración entre el físico Doctor Pablo García Risueño (Humboldt Universität zu Berlin, Alemania) y el grupo de Arquitectura de Computadores de la Universidad de Zaragoza. El trabajo desarrollado parte de la investigación realizada por el Dr. Pablo García Risueño. Además, en el desarrollo del proyecto colabora el investigador Dr. Carl Christian Kjelgaard Mikkelsen (Universidad de Umea, Suecia), matemático Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 10 Fuerza de ligadura sobre el átomo en el instante Multiplicador de Lagrange asociado a la ligadura Ligadura La existencia de las ligaduras convierte un sistema de ecuaciones diferenciales con incógnitas en un sistema de ecuaciones diferenciales con incógnitas. Estas ecuaciones se resuelven en dos fases: en primer lugar se calculan las nuevas posiciones de los átomos sin tener en cuenta las ligaduras impuestas al sistema y en segundo lugar se calculan las fuerzas de ligadura y se corrigen las posiciones. Por lo tanto, una vez resueltas las ecuaciones clásicas del movimiento de Newton, las ecuaciones por las ligaduras son las nuevas ecuaciones que se deben resolver y, los multiplicadores de Lagrange, son las nuevas incógnitas cuyos valores han de ser encontrados para resolver el sistema. Las ligaduras en los enlaces o en los ángulos de enlace pueden ser expresadas del siguiente modo: Ecuación 6. Ecuación que define una ligadura donde es una constante (distancia entre los átomos i y j) y representa la ligadura que implica a los átomos y . Un algoritmo clásico para el tratamiento de las ligaduras es el algoritmo SHAKE. Algoritmo SHAKE El algoritmo SHAKE [1] corrige el conjunto de coordenadas calculadas sin aplicar las ligaduras, , por un conjunto de coordenadas que cumplen con la lista de ligaduras, usando las coordenadas como referencia. En el algoritmo SHAKE se calcula una aproximación de los multiplicadores de Lagrange , resolviendo la siguiente ecuación para cada una de las ligaduras (desde hasta ): Ecuación 7. Sistema de ecuaciones para obtener los multiplicadores de Lagrange Todos los términos excepto son conocidos y: Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 11 donde y representan la delta de Kronecker, función matemática de dos variables, que vale 1 si son iguales, y 0 si son diferentes: Ecuación 8. Delta de Kronecker Se tiene un sistema de ecuaciones cuadráticas en que se puede resolver de forma iterativa. En cada iteración se resuelve, por lo tanto, el siguiente sistema lineal: Ecuación 9. Sistema lineal para obtener de forma iterativa . La resolución de este sistema lineal se lleva a cabo en SHAKE de forma aproximada, resolviendo en primer lugar la primera ecuación, después la segunda y así sucesivamente. La resolución de la ecuación k-ésima garantiza la satisfacción de la ligadura késima una vez refrescadas las posiciones; pero asimismo, rompe parcialmente la satisfacción de las ligaduras con índices menores a . Se obtiene primero , considerando que son iguales a cero en la primera iteración. En las sucesivas iteraciones se usan los obtenidos en la iteración anterior. De esta forma, una vez obtenidas las aproximaciones de los multiplicadores de Lagrange, se sustituyen en la ECUACIÓN 10 para conseguir las nuevas coordenadas corregidas. La posición del átomo en el instante se obtiene a partir de la posición del siguiente modo: Ecuación 10. Ecuación para obtener las posiciones corregidas del los átomos Si las nuevas posiciones corregidas no cumplen las restricciones dentro de una tolerancia relativa, shake_tol, se realiza una nueva iteración para obtener nuevos multiplicadores de Lagrange. El objetivo de este proyecto es mejorar este algoritmo resolviendo el proceso de obtención de las aproximaciones de los multiplicadores de Lagrange de forma más precisa. Es decir, se usan técnicas de resolución de sistemas lineales con matrices banda para resolver el sistema lineal . Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 12 Esquema del algoritmo completo de dinámica molecular A continuación se muestra un esquema general del algoritmo de dinámica molecular: ALGORITMO GLOBAL DE DINÁMICA MOLECULAR 1. Input: Condiciones iniciales Energía potencial del sistema, Posiciones de todos los átomos del sistema, Velocidades de todos los átomos del sistema repetir los pasos 2, 3, 4, 5, 6 para el número de pasos (time steps) requeridos 2. Calcular fuerzas Para cada átomo se calcula sin tener en cuenta las restricciones 3. Actualizar configuración: algoritmo de integración numérica 3.1 Calcular velocidades 3.2 Calcular posiciones 4. Algoritmo Shake repetir 4.1 y 4.2 hasta que todas cumplan las restricciones dentro de una tolerancia 4.1 Resolver el sistema lineal ( y dependen de las posiciones actualizadas ) para obtener los multiplicadores de Lagrange 4.2 Calcular posiciones corregidas 5. Corregir velocidades de acuerdo a las nuevas posiciones 6. Output Escribir las coordenadas, energía… en los ficheros que recogen la trayectoria Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 13 2.2 Propuesta teórica: ILVES-S El Dr. Pablo García Risueño propone un algoritmo alternativo al algoritmo SHAKE, denominado ILVES-S [3]. La motivación para el desarrollo de este nuevo algoritmo es que SHAKE puede ser fuente de inestabilidades y tener problemas de convergencia. Estos problemas son especialmente importantes si se imponen ligaduras a los ángulos entre enlaces (bond angles). El proceso para resolver el sistema lineal puede ser más preciso si se usan técnicas de resolución de matrices banda, ya que la matriz A (véase ECUACIÓN 9) es una matriz dispersa para las moléculas biológicas y puede convertirse en una matriz banda reordenando las ligaduras apropiadamente. Además, el formalismo ofrece la posibilidad de imponer ligaduras en ángulos diedros, algo que actualmente no está contemplado con el algoritmo SHAKE. Por lo tanto, se esperan alcanzar mejoras en la eficiencia, exactitud y estabilidad de las simulaciones de dinámica molecular usando el algoritmo ILVES-S para resolver de forma analítica el sistema lineal de la ECUACIÓN 9. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 14 Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 15 Capítulo 3 GROMACS En este capítulo se introduce el paquete software de dinámica molecular GROMACS y se expone el análisis realizado del algoritmo SHAKE. Como ejemplo para el análisis se presenta la glicina, el aminoácido más pequeño. 3.1 Introducción GROMACS (GROningen MAchine for Chemical Simulations) [2] es un versátil paquete de dinámica molecular. Es utilizado, entre otras cosas, para simular las ecuaciones newtonianas (F = ma) de movimiento para sistemas con cientos de millones de partículas. GROMACS fue originalmente desarrollado en la Universidad de Groningen (Países Bajos) pero en la actualidad se mantiene por desarrolladores de diversas universidades y centros de investigación (Suecia, Alemania, EEUU, Australia) [14]. Las rutinas de dinámica molecular están escritas en el lenguaje de programación C y basadas en el programa GROMOS, desarrollado por el mismo grupo de la Universidad de Groningen. GROMACS es uno de los paquetes de software más rápido y popular disponibles y se puede ejecutar tanto en CPUs como en GPUs. Es de código abierto y libre, publicado bajo una Licencia Pública General de GNU (General Public License). GROMACS está incluido en SPEC CPU 2006 [18], un conjunto de programas intensivos en cálculo que se utilizan para medir y comparar el rendimiento sistemas informáticos. Además es ampliamente utilizado en los denominados proyectos de computación distribuida como Folding@home [15], de la Universidad de Stanford, o EvoGrid [16]. También es utilizado como una de las aplicaciones que forman parte del ScalaLife Competence Center [17]. Proyectos como el Folding@home contribuyen en la investigación de enfermedades como el alzhéimer, la fibrosis quística o el cáncer. Para ello, se realizan simulaciones del plegamiento de proteínas1. La función biológica de una proteína depende de su correcto plegamiento. Si una proteína no se pliega correctamente no será funcional y, por lo tanto, no será capaz de cumplir su función biológica. Las simulaciones permiten determinar las causas de un mal plegamiento. 1 Proceso por el que una proteína alcanza su estructura tridimensional final. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 16 En menor medida, estos proyectos también son capaces de determinar la estructura final de una proteína o cómo éstas interactúan con otras moléculas, lo que tiene aplicaciones útiles en el diseño computacional de medicamentos. 3.2 Análisis del algoritmo SHAKE en GROMACS. La implementación del algoritmo ILVES-S se desarrolla en el paquete software GROMACS debido a su extendido uso y la importante repercusión que puede tener. Para ello es necesario realizar un análisis de una simulación de dinámica molecular primero, y del código del actual algoritmo que implementa la imposición de ligaduras (SHAKE) después. La ejecución detallada de una simulación de dinámica molecular puede encontrarse en el Anexo II. Para realizar este análisis es necesario recompilar GROMACS y modificar su fichero de configuración (configure) para cambiar algunos flags de compilación. Se elimina el flag -O3 y se añaden los flags -O0 y -g. Esto permite depurar el código usando el depurador GNU DDD [8]. Además también es necesario eliminar el flag -fomit-frame-pointer puesto que es incompatible con el flag -g. A continuación se detalla el análisis realizado del código de GROMACS que implementa el algoritmo SHAKE. Este análisis es importante para entender como se ha llevado a cabo la implementación del algoritmo que se va a sustituir y también para familiarizarse con el estilo de codificación y las variables de GROMACS. Datos de entrada Dado que en GROMACS se implementan dos algoritmos de tratamiento de ligaduras (SHAKE y LINCS), es posible elegir uno de ellos mediante un parámetro en el fichero de entrada donde se establecen los parámetros de la simulación: el fichero .mdp (Ver Anexo II): constraint_algorithm = SHAKE ó LINCS Otros parámetros interesantes que también se incluyen en este fichero son: integrator Algoritmo de integración numérica nsteps Número máximo de pasos temporales (time steps) de la simulación dt Valor del paso temporal en picosegundos shake_tol Tolerancia para el algoritmo SHAKE (por defecto su valor es 0,0001) constraints Tipo de ligaduras que se implementan Este último parámetro es muy importante ya que indica que tipo de ligaduras se desean imponer. Existen cinco posibles valores para este parámetro: none No se imponen ligaduras hbonds Se imponen ligaduras sólo a los enlaces que involucran átomos de Hidrógeno all-bonds Se imponen ligaduras a todos los enlaces Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 17 h-angles Se imponen ligaduras a todos los enlaces y además a los ángulos que involucran átomos de Hidrógeno all-angles Se imponen ligaduras a todos los enlaces y a todos los ángulos La glicina Para el análisis de la implementación del algoritmo SHAKE en GROMACS se realizan simulaciones de la glicina (el aminoácido más sencillo) en su estado zwitterion. Es un residuo con 10 átomos y 9 enlaces. La FIGURA 3 muestra la composición y una representación del compuesto escogido. Figura 3. Representación de la Glicina en su estado zwitterion La composición de la glicina es la siguiente:  Átomo azul: Nitrógeno  Átomos blancos: Hidrógeno  Átomos negros: Carbono  Átomos rojos: Oxígeno La FIGURA 4 muestra la glicina tal y como se almacena en GROMACS: los números de átomos y enlaces se utilizan en el acceso a los vectores que contienen su descripción (masas, posiciones, velocidad, tipo…) y el tipo de enlace determina la ligadura asociada al mismo. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 18 Figura 4. Representación de la Glicina en GROMACS Estructuras de datos y variables En esta sección se presentan las estructuras de datos que se usan en GROMCAS para la implementación del algoritmo SHAKE. Para describirlas y mostrar sus valores (si son constantes a lo largo de toda la simulación) se usa como ejemplo la glicina y se imponen ligaduras a todos sus enlaces, es decir, el parámetro constraints, descrito al comienzo de este capítulo es igual a all-bonds. Primero se enumeran las variables y después se muestra una imagen descriptiva de las mismas:  int ncons Constante que representa el número de ligaduras impuestas al sistema (number of constraints). En el ejemplo de la glicina su valor es 9 (una ligadura por cada enlace).  int natoms Número total de átomos (incluyendo los átomos del disolvente, agua por ejemplo). Para el ejemplo que nos ocupa, su valor es 718.  real invdt Inversa de dt. dt es un parámetro de entrada en el fichero *.mdp definido como el paso temporal (time step) en picosegundos. Para el ejemplo, dt tiene el valor de 0.002 ps, por lo tanto, invdt es 500. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 19  real shake_tol Valor usado para calcular la tolerancia relativa bajo la cual todas las ligaduras han de satisfacerse. Para este ejemplo se toma su valor por defecto: 0,0001.  t_iatom *iatoms (t_iatom = real) Vector que describe las ligaduras mediante grupos de tres elementos. El primero, describe el tipo de ligadura y el segundo y tercer elementos son los dos átomos involucrados en la ligadura. El tamaño del vector es, por lo tanto, tres veces el número de ligaduras (3*ncons). Si k es el número de la ligadura de tipo t que involucra a los átomos i y j entonces: iatoms[3*k] = t; iatoms[3*k+1] = i; iatoms[3*k+2] = j  real invmass[] Vector que contiene la inversa de la masa de los átomos. En el caso de la glicina el tamaño del vector es 10 y sus valores pueden verse en la FIGURA 5.  real lagr Vector que contiene los multiplicadores de Lagrange. Su tamaño coincide con el número de ligaduras. En el caso de la glicina, su tamaño es 9  real dist2[] Vector que contiene el valor de las ligaduras al cuadrado, es decir, (véase ECUACIÓN 6). El tamaño del vector es igual al número de ligaduras; para la glicina es 9.  real x[] Vector que contiene las posiciones (x,y,z) de los átomos al inicio del paso temporal, es decir, contiene los valores . El tamaño del vector es 3 veces el número de átomos. En el caso de la glicina son 10 átomos, por lo tanto el tamaño del vector es 30.  real xp[] Vector que contiene las posiciones (x,y,z) de los átomos al final del paso temporal. Antes de ejecutar el algoritmo SHAKE contiene las posiciones antes de ser corregidas, es decir, y después de ejecutar SHAKE contiene las posiciones corregidas, es decir, . Al igual que el vector x[], su tamaño es 30.  real xij[] Vector que contiene la distancia (x,y,z) entre los átomos ligados en el paso temporal anterior, es decir cuando existe una ligadura k que relaciona los átomos i y j. El tamaño de este vector para el ejemplo de la glicina es 3 veces el número de ligaduras del sistema: 3*9 = 27. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 26 ii -> k(i,j) jj -> k’(l,m) weight 0 (0,1) 0 (0,1) 0 (0,1) 1 (0,3) 0 (0,1) 2 (0,2) 0 (0,1) 3 (0,4) 1 (0,3) 0 (0,1) 1 (0,3) 1 (0,3) 1 (0,3) 2 (0,2) 1 (0,3) 3 (0,4) 2 (0,2) 0 (0,1) 2 (0,2) 1 (0,3) 2 (0,2) 2 (0,2) 2 (0,2) 3 (0,4) 3 (0,4) 0 (0,1) 3 (0,4) 1 (0,3) 3 (0,4) 2 (0,2) 3 (0,4) 3 (0,4) 3 (0,4) 4 (4,6) 3 (0,4) 5 (4,5) 3 (0,4) 6 (4,7) 4 (4,6) 3 (0,4) 4 (4,6) 4 (4,6) 4 (4,6) 5 (4,5) 4 (4,6) 6 (4,7) 5 (4,5) 3 (0,4) 5 (4,5) 4 (4,6) 5 (4,5) 5 (4,5) 5 (4,5) 6 (4,7) 6 (4,7) 3 (0,4) 6 (4,7) 4 (4,6) 6 (4,7) 5 (4,5) 6 (4,7) 6 (4,7) 6 (4,7) 7 (7,9) 6 (4,7) 8 (7,8) 7 (7,9) 6 (4,7) 7 (7,9) 7 (7,9) 7 (7,9) 8 (7,8) 8 (7,8) 6 (4,7) 8 (7,8) 7 (7,9) 8 (7,8) 8 (7,8) Figura 12. Elementos no nulos de la matriz A para el ejemplo de la glicina. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 27 4.3 Proceso iterativo Generación de la matriz A Para resolver el sistema lineal se hace uso de una librería externa denominada LAPACK (Linear Algebra PACKage) [9] que proporciona rutinas para, entre otras cosas, resolver sistemas de ecuaciones lineales. Una vez generada la tabla presentada anteriormente, es muy sencillo construir la matriz A: simplemente hay que añadir a cada entrada el término . El mayor problema en este paso es almacenar la matriz en el formato de entrada para las funciones LAPACK. En la siguiente figura se puede observar a la derecha la matriz A para la glicina, y a la izquierda, la misma matriz en formato LAPACK: Figura 13. A la derecha, matriz A para la glicina. A la izquierda, la misma matriz en formato LAPACK. Generación del vector b El vector se obtiene a partir de los vectores dist2[] y xp[] del siguiente modo: Ecuación 13. Vector b del sistema lineal Resolución del sistema lineal La resolución del sistema lineal se lleva a cabo mediante la llamada a dos funciones de la librería externa LAPACK. En concreto, para precisión doble se usan las funciones dgbtrf() y dgbtrs(). La primera rutina calcula la factorización LU de la matriz A, mientras que la Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 28 segunda resuelve el sistema de ecuaciones utilizando la factorización LU realizada por la primera función. Corrección de posiciones Por último, la corrección de posiciones se realiza aplicando la siguiente ecuación a cada átomo implicado en las ligaduras: Ecuación 14. Corrección de posiciones Condición de parada Una vez aplicadas todas las operaciones hay que comprobar si se aceptan las nuevas posiciones, es decir si las ligaduras impuestas al sistema se satisfacen dentro de una tolerancia dada. Si no es así, será necesario realizar una nueva iteración. Para ello se hace uso de la misma tolerancia usada por el algoritmo SHAKE: shake_tol, que por defecto vale 0.0001. A partir de esta variable, que se puede cambiar en el fichero de entrada (.mdp), se define para cada ligadura: Ecuación 15. Condición de parada para el algoritmo ILVES-S Si es mayor que 1 para alguna de las ligaduras, entonces es necesario realizar una nueva iteración. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 29 Capítulo 5 RESULTADOS Este capítulo presenta el entorno experimental (máquina en la que se ejecutan las simulaciones, molécula simulada) y los principales resultados obtenidos. 5.1 Entorno experimental Maquina de trabajo La máquina de trabajo en la que se ejecutan los experimentos tiene las siguientes características:  Sistema operativo: CentOS release 6.3  Kernel: Linux 2.6.32-220.2.1.el6.x86_64  Procesador: Intel(R) Core(TM) i7-2600 CPU @ 3.40GHz (microarquitectura Sandy Bridge): 4 núcleos físicos con hyperthreading (8 núcleos lógicos).  3 niveles de cache integrados en el chip: o L1I + L1D: 32 kB + 32 kB (por núcleo) o L2: 256 kB (por núcleo) o L3: 8 MB (compartida)  Compilador: gcc versión 4.4.6 20120305 (Red Hat 4.4.6-4) (GCC) Simulaciones de dinámica molecular La verificación del correcto funcionamiento del algoritmo se ha realizado mediante sencillos experimentos. Primero con la simulación de la glicina, el aminoácido más pequeño, y después con cadenas de tres aminoácidos. Para evaluar el resultado de la implementación de la primera versión del algoritmo ILVES-S y cuantificar su rendimiento y precisión, se realizan simulaciones de dinámica molecular de la arginina, uno de los aminoácidos esenciales que forma parte de las proteínas. Esta molécula tiene 27 átomos y 26 enlaces que dan lugar a 26 ligaduras. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 30 5.2 Verificación de los resultados Para realizar la verificación del correcto funcionamiento de ILVES-S, además de realizar sencillos test en los que se comparan los multiplicadores de Lagrange obtenidos con el algoritmo SHAKE y los que se obtienen con ILVES-S, es importante comprobar que las posiciones corregidas cumplen las restricciones dentro de una cierta tolerancia. Así, el error relativo se calcula del siguiente modo: Ecuación 16. Error absoluto cometido en el cálculo de las nuevas posiciones Además, el código implementado e integrado en GROMACS y el código GROMACS original se han modificado para volcar en un fichero los datos necesarios que permiten: (a) comparar el algoritmo original de GROMACS (SHAKE) y el nuevo algoritmo implementado (ILVES-S), y (b) comprobar el correcto funcionamiento de este último. También se han programado scripts que obtienen de forma automatizada todas estas medidas con diferentes entradas y configuraciones del algoritmo. 5.3 Análisis de los resultados Precisión Para comparar las precisiones de ambos algoritmos hemos utilizado el error relativo máximo y el error relativo medio. En cada paso temporal de simulación, se ha registrado el error relativo máximo y el error relativo medio de todas las ligaduras (valores absolutos, ver ECUACIÓN 16). Posteriormente se han promediado los errores registrados en cada paso de simulación. La FIGURA 14 muestra el error relativo máximo y el error relativo medio de cada algoritmo para distintos valores de tolerancia relativa, shake_tol, especificada en los parámetros de simulación. En primer lugar, se observa que ILVES-S resuelve correctamente la imposición de ligaduras en simulaciones de dinámica molecular. Además, para todos los valores de tolerancia con los que se ha experimentado, ILVES es más preciso que SHAKE, tanto en error relativo máximo como en error relativo medio. Reseñar que para una tolerancia de 1e-4 (el valor de shake_tol por defecto en GROMACS), el error cometido por ILVES-S es del orden de tres órdenes de magnitud menor que el de SHAKE. Por lo tanto, las simulaciones con ILVES-S son más fieles a la realidad física que modelan. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 31 Figura 14. Error relativo máximo y error relativo medio para distintos valores de tolerancia relativa (shake_tol). Observar que la escala del eje y es logarítmica. Las variaciones de precisión observadas en la anterior gráfica pueden explicarse analizando el número medio de iteraciones ejecutado por cada algoritmo. En cada paso temporal de simulación, se ha registrado el número de iteraciones ejecutadas por ambos algoritmos de cálculo (SHAKE e ILVES-S). Posteriormente se han promediado los valores registrados en cada paso de simulación. En la FIGURA 15 se puede observar el número medio de iteraciones de cada algoritmo para distintos valores de tolerancia relativa shake_tol. Para el caso de tolerancias menores o iguales a 1e-6, ILVES-S es capaz de cumplir las ligaduras ejecutando una media de en torno a dos iteraciones, mientras que las tolerancias mayores (1e-7 y 1e-8) requieren la ejecución de tres iteraciones. El comportamiento de SHAKE es distinto, ya que el número de iteraciones necesario para cumplir las ligaduras crece al aumentar la precisión requerida. Para una tolerancia de 1e-6, SHAKE ejecuta una media de casi 11 iteraciones frente a las 2 que ejecuta ILVES-S. Es de esperar que en moléculas con muchas ligaduras esto se traduzca en un menor tiempo de ejecución por parte de ILVES-S. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 32 Figura 15. Número medio de iteraciones efectuadas por los algoritmos de cálculo para distintos valores de tolerancia relativa (shake_tol). Tiempo Durante cada simulación, GROMACS registra el tiempo (medido en ciclos de CPU) que dedica a distintas tareas. Una de ellas es el cumplimiento de las ligaduras. La FIGURA 16 muestra el número de ciclos de CPU dedicados a los cálculos relacionados con las ligaduras. Para todas las tolerancias, SHAKE es más rápido, con diferencias que oscilan entre el 4.8% y el 7.3%. Estas diferencias son debidas a varios factores:  La implementación de SHAKE en GROMACS está muy optimizada, mientras que la versión evaluada de ILVES-S es una primera implementación centrada en analizar su potencial y no en optimizar su rendimiento.  La versión implementada del algoritmo ILVES-S utiliza funciones de la librería LAPACK para resolver el sistema lineal de ecuaciones. Actualmente, el equipo de investigación (en concreto, el matemático Carl Christian Kjelgaard Mikkelsen) está trabajando en una función para resolver sistemas de ecuaciones lineales banda, lo que permitirá reducir de forma significativa el tiempo de ejecución de ILVES-S. Así pues, se espera que próximamente ILVES-S resuelva las ligaduras más rápidamente que SHAKE. Además, a medio plazo se pretende implementar una versión paralela de ILVES-S, lo cual permitirá mejorar todavía más sus prestaciones (la implementación de SHAKE en GROMACS es secuencial). Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 33 Por otra parte, puede observarse que el tiempo de ejecución no escala con el número de iteraciones. Esto puede ser debido al reducido número de ligaduras. Próximamente, esperamos realizar experimentos con moléculas más grandes y con más ligaduras para validar esta hipótesis. Figura 16. Número de ciclos de CPU (en miles de millones) invertidos en cálculos relacionados con las ligaduras. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 34 Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 35 Capítulo 6 CONCLUSIONES Este último capítulo recoge las conclusiones del proyecto fin de carrera, el trabajo futuro y la valoración personal. 6.1 Conclusiones La dinámica molecular es una técnica usada en todo el mundo en numerosas áreas relacionadas con la física, la química y la biología, donde es necesario comprender, estudiar y analizar el comportamiento de los átomos y las moléculas. Es un campo multidisciplinario que requiere el conocimiento de expertos en diversas áreas como las matemáticas, la física, la química y las ciencias computacionales. Ejemplo de ello es el proyecto del que forma parte este proyecto fin de carrera y en el que participa el grupo de Arquitectura de Computadores de la Universidad de Zaragoza, el Dr. Pablo García Risueño (Universidad Humboldt de Berlín, Alemania) experto en física teórica y el Dr. Carl Christian Kjelgaard Mikkelsen (Universidad de Umea, Suecia), matemático especializado en cálculos analíticos y resolución de grandes sistemas de ecuaciones lineales. Se ha integrado en GROMACS la primera versión de un algoritmo, ILVES-S, que permite imponer ligaduras en simulaciones de dinámica molecular con ligaduras. Para una tolerancia dada, ILVES-S es más preciso que el algoritmo original, SHAKE. La solución que se obtiene ahora se aproxima más a la realidad física. En cuanto al tiempo de ejecución, ILVES-S es ligeramente más lento que SHAKE (entre un 4.8% y un 7.3% según las distintas tolerancias). Esta diferencia es debida en gran parte al uso de rutinas de resolución de sistemas lineales LAPACK que no explotan eficientemente el hecho de que son sistemas banda. 6.2 Trabajo futuro Actualmente hay dos líneas abiertas de trabajo en el algoritmo ILVES-S: - Sustituir las funciones de resolución de sistemas de ecuaciones lineales de la librería LAPACK por funciones desarrolladas específicamente resolver los sistema de ecuaciones que surgen Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 42  GROMACS v.4.5.5 (Doble precisión): ./configure --prefix=/home/maston/gromacs --disable-float --enableshared make make install Los flags incluidos en la configuración de las instalaciones son: --prefix para cambiar el directorio destino de instalación. --enable-float para compilar con precisión simple. --disable-float para compilar con doble precisión. -–enable-shared librería dinámica [24]. 2. Ejecución En esta sección se detalla en qué consisten los distintos pasos para llevar a cabo, en GROMACS, una simulación de dinámica molecular partiendo del denominado fichero PDB (Protein Data Bank [25]). Este fichero recoge la estructura de la molécula que se desea simular. Como ejemplo de ejecución se usa una lisozima, una encima presente en los huevos de gallina que tiene 1079 átomos y cuyo código PDB es 1AKI. Su estructura, recogida en el fichero .pdb, es la que se observa en la FIGURA 18. Figura 18. Estructura de la lysozima Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 43 PASO 1. Preparar la topología En primer lugar, partiendo del fichero .pdb, se generan algunos ficheros necesarios para la simulación con la herramienta pdb2gmx: pdb2gmx -f 1aki.pdb -o 1aki _processed.gro -water spce Se selecciona el force field número 14, que corresponde a OPLS-AA/L all-atom force field (2001 aminoacid dihedrals). Los ficheros obtenidos son:  Topología de la molécula (fichero topol.top): contiene toda la información necesaria para definir la molécula en una simulación. Incluye información no relacionada con los enlaces (tipos de átomos y sus cargas) así como parámetros relacionados con los enlaces (enlaces covalentes, ángulos de enlace y ángulos diedros).  Fichero con la definición de las ligaduras (fichero posre.itp)  Estructura post-procesada (fichero 1aki_processed.gro) Además, se obtienen algunos datos importantes en la salida estándar: ... Read 'LYSOZYME', 1079 atoms Analyzing pdb file Splitting chemical chains based on TER records or chain id changing. There are 1 chains and 1 blocks of water and 207 residues with 1079 atoms chain #res #atoms 1 'A' 129 1001 2 ' ' 78 78 (only water) … Now there are 129 residues with 1960 atoms Making bonds... Number of bonds was 1984, now 1984 Generating angles, dihedrals and pairs... Before cleaning: 5142 pairs Before cleaning: 5187 dihedrals Keeping all generated dihedrals Making cmap torsions...There are 5187 dihedrals, 426 impropers, 3547 angles, 5106 pairs, 1984 bonds and 0 virtual sites Total mass 14313.197 a.m.u. Total charge 8.000 e … Processing chain 2 (78 atoms, 78 residues) … Now there are 78 residues with 234 atoms Making bonds... Number of bonds was 156, now 156 Generating angles, dihedrals and pairs... Making cmap torsions...There are 0 dihedrals, 0 impropers, 78 angles, 0 pairs, 156 bonds and 0 virtual sites Total mass 1405.201 a.m.u. Total charge 0.000 e Including chain 1 in system: 1960 atoms 129 residues Including chain 2 in system: 234 atoms 78 residues Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 44 Now there are 2194 atoms and 207 residues Total mass in system 15718.398 a.m.u. Total charge in system 8.000 e … PASO 2. Definición del recipiente. A continuación se define el recipiente que contendrá la lisozima (en este caso un sencillo recipiente con forma cúbica), usando para ello la herramienta editconf: editconf -f 1aki_processed.gro -o 1aki_nexbox.gro -c -d 1.0 -bt cubic Se obtiene un nuevo fichero *.gro, en este caso 1aki_nexbox.gro, con la nueva estructura y la siguiente información: ... Read 1960 atoms Volume: 123.376 nm^3, corresponds to roughly 55500 electrons No velocities found system size : 3.817 4.234 3.454 (nm) diameter : 5.010 (nm) center : 2.781 2.488 0.017 (nm) box vectors : 5.906 6.845 3.052 (nm) box angles : 90.00 90.00 90.00 (degrees) box volume : 123.38 (nm^3) shift : 0.724 1.017 3.488 (nm) new center : 3.505 3.505 3.505 (nm) new box vectors : 7.010 7.010 7.010 (nm) new box angles : 90.00 90.00 90.00 (degrees) new box volume : 344.48 (nm^3) ... PASO 3. Definición del disolvente En este paso es necesario “añadir” al recipiente creado en el paso anterior, el disolvente deseado, en este caso agua. Para ello se usa la herramienta genbox: genbox -cp 1aki_nexbox.gro -cs spc216.gro -o 1aki_solv.gro -p topol.top Se obtiene otra vez un nuevo fichero con la estructura de la molécula, 1aki_solv.gro, y se actualiza el fichero que contiene la topología, topol.top. Como salida obtenemos la siguiente información relevante: ... Reading solvent configuration "216H2O,WATJP01,SPC216,SPC-MODEL,300K,BOX(M)=1.86206NM,WFVG,MAR. 1984" solvent configuration contains 648 atoms in 216 residues Initialising van der waals distances... Will generate new solvent configuration of 4x4x4 boxes Generating configuration Sorting configuration Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 45 Found 1 molecule type: SOL ( 3 atoms): 13824 residues Calculating Overlap... box_margin = 0.315 Removed 2736 atoms that were outside the box Neighborsearching with a cut-off of 0.48 Table routines are used for coulomb: FALSE Table routines are used for vdw: FALSE Cut-off's: NS: 0.48 Coulomb: 0.48 LJ: 0.48 System total charge: 0.000 Grid: 16 x 16 x 16 cells Successfully made neighbourlist nri = 61846, nrj = 1388188 Checking Protein-Solvent overlap: tested 36470 pairs, removed 1968 atoms. Checking Solvent-Solvent overlap: tested 390293 pairs, removed 4272 atoms. Added 10832 molecules Generated solvent containing 32496 atoms in 10832 residues Writing generated configuration to 1AKI_solv.gro LYSOZYME Output configuration contains 34456 atoms in 10961 residues Volume : 344.484 (nm^3) Density : 1173.99 (g/l) Number of SOL molecules: 10832 ... PASO 4. Añadir iones En este paso se añaden iones al sistema, para ello se usa la herramienta genion pero antes es necesario ensamblar (grompp) los parámetros especificados en un nuevo fichero de entrada (ions.mdp) con la topología y las coordenadas del sistema para obtener un fichero que contiene todos los parámetros para todos los átomos del sistema (ions.tpr). Para este ejemplo, el fichero de entrada usado es: ; ions.mdp - used as input into grompp to generate ions.tpr ; Parameters describing what to do, when to stop and what to save Integrator = steep ; Algorithm (steep = steepest descent minimization) emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm emstep = 0.01 ; Energy step size nsteps = 50000 ; Maximum number of (minimization) steps to perform ; Parameters describing how to find the neighbors of each atom and how to calcúlate ; the interactions nstlist = 1 ; Frequency to update the neighbor list and long range forces ns_type = grid ; Method to determine neighbor list (simple, grid) rlist = 1.0 ; Cut-off for making neighbor list (short range forces) coulombtype = PME ; Treatment of long range electrostatic interactions rcoulomb = 1.0 ; Short-range electrostatic cut-off rvdw = 1.0 ; Short-range Van der Waals cut-off pbc = xyz ; Periodic Boundary Conditions (yes/no) Ensamblado: grompp -f ions.mdp -c 1aki_solv.gro -p topol.top -o ions.tpr Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 46 La salida obtenida es: … Analysing residue names: There are: 129 Protein residues There are: 10824 Water residues There are: 8 Ion residues Analysing Protein... Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups... Number of degrees of freedom in T-Coupling group rest is 70845.00 Largest charge group radii for Van der Waals: 0.163, 0.157 nm Largest charge group radii for Coulomb: 0.164, 0.163 nm Calculating fourier grid dimensions for X Y Z Using a fourier grid of 60x60x60, spacing 0.117 0.117 0.117 Estimate for the relative computational load of the PME mesh part: 0.30 … Añadir iones: genion_d -s ions.tpr -o 1aki_solv_ions.gro -p topol.top -pname NA -nname CL –nn 8 En este caso se selecciona la opción 13 (Group 13 (SOL) has 32496 elements). Información relevante en la salida estándar: … Number of (3-atomic) solvent molecules: 10832 Processing topology Replacing 8 solute molecules in topology file (topol.top) by 0 NA and 8 CL ions. Back Off! I just backed up topol.top to ./#topol.top.2# Replacing solvent molecule 1270 (atom 5770) with CL Replacing solvent molecule 8207 (atom 26581) with CL Replacing solvent molecule 5287 (atom 17821) with CL Replacing solvent molecule 9164 (atom 29452) with CL Replacing solvent molecule 3607 (atom 12781) with CL Replacing solvent molecule 1735 (atom 7165) with CL Replacing solvent molecule 4182 (atom 14506) with CL Replacing solvent molecule 9597 (atom 30751) with CL … PASO 5: Minimización energética Ahora es necesario asegurarse de que el sistema no tiene la geometría inapropiada para ello se “relaja” la estructura mediante el proceso denominado minimización energética (Energy Minimization, EM). El proceso es similar al anterior: primero se realiza un ensamblado de varios ficheros de entrada en un único fichero (1aki_em.tpr) y después se realiza el proceso de minimización energética. El fichero de entrada minim.mdp es: Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 47 ; minim.mdp - used as input into grompp to generate 1aki_em.tpr ; Parameters describing what to do, when to stop and what to save integrator = steep ; Algorithm (steep = steepest descent minimization) emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm emstep = 0.01 ; Energy step size nsteps = 50000 ; Maximum number of (minimization) steps to perform ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions nstlist = 1 ; Frequency to update the neighbor list and long range forces ns_type = grid ; Method to determine neighbor list (simple, grid) rlist = 1.0 ; Cut-off for making neighbor list (short range forces) coulombtype = PME ; Treatment of long range electrostatic interactions rcoulomb = 1.0 ; Short-range electrostatic cut-off rvdw = 1.0 ; Short-range Van der Waals cut-off pbc = xyz ; Periodic Boundary Conditions (yes/no) Ensamblado: grompp -f minim.mdp -c 1aki_solv_ions.gro -p topol.top -o 1aki_em.tpr Minimización energética: mdrun_d -v -deffnm 1aki_em La salida estándar es: … Steepest Descents converged to Fmax < 1000 in 979 steps Potential Energy = -6.0590838e+05 Maximum force = 9.3957776e+02 on atom 736 Norm of force = 1.9575663e+01 … PASO 6: Equilibrado del disolvente y los iones que rodean a la proteína En este paso está dividido en dos fases y se equilibra la temperatura y la presión del sistema. Ambos procesos se realizan de forma similar a los anteriores: ensamblado + ejecución. Por lo tanto se necesitan dos nuevos ficheros de entrada nvt.mdp para la primera fase y npt.mdp para la segunda fase. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 48 title = OPLS Lysozyme NVT equilibration define = -DPOSRES ; position restrain the protein ; Run parameters integrator = md ; leap-frog integrator nsteps = 50000 ; 2 * 50000 = 100 ps dt = 0.002 ; 2 fs ; Output control nstxout = 100 ; save coordinates every 0.2 ps nstvout = 100 ; save velocities every 0.2 ps nstenergy = 100 ; save energies every 0.2 ps nstlog = 100 ; update log file every 0.2 ps ; Bond parameters continuation = no ; first dynamics run constraint_algorithm = lincs ; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs rlist = 1.0 ; short-range neighborlist cutoff (in nm) rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm) rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.16 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = Protein Non-Protein ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 300 300 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = no ; no pressure coupling in NVT ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr = EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp = 300 ; temperature for Maxwell distribution gen_seed = -1 ; generate a random seed Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 49 title = OPLS Lysozyme NPT equilibration define = -DPOSRES ; position restrain the protein ; Run parameters integrator = md ; leap-frog integrator nsteps = 50000 ; 2 * 50000 = 100 ps dt = 0.002 ; 2 fs ; Output control nstxout = 100 ; save coordinates every 0.2 ps nstvout = 100 ; save velocities every 0.2 ps nstenergy = 100 ; save energies every 0.2 ps nstlog = 100 ; update log file every 0.2 ps ; Bond parameters continuation = yes ; Restarting after NVT constraint_algorithm = lincs ; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs rlist = 1.0 ; short-range neighborlist cutoff (in nm) rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm) rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.16 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = Protein Non-Protein ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 300 300 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1 refcoord_scaling = com ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr = EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = no ; Velocity generation is off Ensamblado y equilibrio de la temperatura: grompp_d -f nvt.mdp -c 1aki_em.gro -p topol.top -o 1aki_nvt.tpr mdrun_d -deffnm 1aki_nvt Ensamblado y equilibrio de la presión: grompp_d -f npt.mdp -c 1aki_nvt.gro -t 1aki_nvt.cpt -p topol.top -o 1aki_npt.tpr mdrun_d -deffnm 1aki_npt Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 50 PASO 5. Simulación de dinámica molecular El último paso es la simulación de dinámica molecular y, de nuevo, es necesario realizar un ensamblado previo utilizando un fichero *.mdp (md.mdp): title = OPLS Lysozyme MD ; Run parameters integrator = md ; leap-frog integrator nsteps = 500000 ; 2 * 500000 = 1000 ps, 1 ns dt = 0.002 ; 2 fs ; Output control nstxout = 1000 ; save coordinates every 2 ps nstvout = 1000 ; save velocities every 2 ps nstxtcout = 1000 ; xtc compressed trajectory output every 2 ps nstenergy = 1000 ; save energies every 2 ps nstlog = 1000 ; update log file every 2 ps ; Bond parameters continuation = yes ; Restarting after NPT constraint_algorithm = lincs ; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs rlist = 1.0 ; short-range neighborlist cutoff (in nm) rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm) rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.16 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = Protein Non-Protein ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 300 300 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr = EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = no ; Velocity generation is off Ensamblado: grompp -f md.mdp -c 1aki_npt.gro -t 1aki_npt.cpt -p topol.top -o 1aki_md.tpr Simulación: mdrun -deffnm 1aki_md Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 51 La ejecución de mdrun genera cuatro ficheros de salida:  Fichero con todas las variables de estado (en intervalos temporales) del sistema simulado (fichero 1aki_md.cpt)  Información acerca de la energía del sistema (fichero 1aki_md.edr)  Estructura post-procesada (fichero 1aki_md.gro)  Log de la ejecución (fichero 1aki_md.log)  Trayectoria de la molécula: coordenadas, velocidades y fuerzas sobre el sistema (fichero 1aki_md.trr)  Versión comprimida de la trayectoria que contiene solo información sobre las coordenadas, tiempo y el recipiente (fichero 1aki_md.xtc) A continuación se muestra en forma de esquema un resumen de los distintos pasos de ejecución necesarios para llevar a cabo una simulación de dinámica molecular. Los recuadros en gris indican la herramienta que ha de utilizarse y junto a ellos se describen los ficheros de entrada y los de salida. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 58 Los índices de las filas y columnas de la matriz A son los índices de las ligaduras, es decir, son matrices cuadradas. Además, de la ecuación 9 se deduce que las matrices A representan con un valor distinto de 0 dos ligaduras que tienen algún átomo en común y con un 0, dos ligaduras que no tienen átomos comunes. Por lo tanto, debido a la estructura característica de las moléculas, las matrices A son matrices dispersas que, como se ha mencionado anteriormente, se pueden convertir en matrices banda reordenando las ligaduras. Un espacio en blanco en la matriz representa un valor distinto de cero. Por ejemplo, si nos fijamos en el elemento , es un valor distinto de cero porque las ligaduras y tienen en común el átomo . Sin embargo, si nos fijamos en el elemento , su valor es cero porque las ligaduras y no tienen ningún átomo en común. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 59 ANEXO IV Trabajo desarrollado en el ICHEC (Irlanda) Los conocimientos de dinámica molecular adquiridos durante el desarrollo de este proyecto fin de carrera, me han ayudado a conseguir una beca competitiva convocada por PRACE (Partnership for Advanced Computing in Europe) [19], para realizar una estancia durante el verano en el Ireland's High-Performance Computing Centre (ICHEC) [20] en Dublín, Irlanda. El resultado del proyecto desarrollado denominado “Profiling and optimization of th the hybrid Molecular Dynamics code, DL_POLY, on heterogeneous clusters”, puede consultarse a partir de la siguiente página. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 60 Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 61 Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 62 © Copyright PRACE AISBL, 2013, all rights reserved. Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 63 REFERENCIAS [1] Ryckaert, J. P., Ciccotti, G., Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints; molecular dynamics of n-alkanes. J.Comp. Phys. 23:327–341, 1977. [2] http://www.gromacs.org Fecha última consulta: 17-agosto-2013 [3] García-Risueño, P. Constraint implementation based on analytical calculations: a posible way to improve widely used solvers. Technical Report. [4] Hess, B.; Bekker, H.; Berendsen, H. J. C. and Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J.Comp. Chem. 18: 1463-1472, 1997 [5] http://www.fftw.org Fecha última consulta: 17-agosto-2013 [6] http://www.cgl.ucsf.edu/chimera/ Fecha última consulta: 17-agosto-2013 [7] http://avogadro.openmolecules.net/wiki/Main_Page Fecha última consulta: 17-agosto-2013 [8] http://www.gnu.org/software/ddd/ Fecha última consulta: 17-agosto-2013 [9] http://www.netlib.org/lapack/ Fecha última consulta: 17-agosto-2013 [10] http://www.doxygen.org Fecha última consulta: 19-agosto-2013 [11] http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/index.html Fecha última consulta: 19-agosto-2013 [12] Hockney, R. W., Goel, S. P., Eastwood, J. Quiet High Resolution Computer Models of a Plasma. J. Comp. Phys. 14:148–158, 1974. [13] Swope, W. C., Andersen, H. C., Berens, P. H., Wilson, K. R. A computer-simulation method for the calculation of equilibrium-constants for the formation of physical clusters of molecules: Application to small water clusters. J. Chem. Phys. 76:637–649, 1982. [14] http://www.gromacs.org/About_Gromacs/People Fecha última consulta: 1-Agosto-2013 [15] http://folding.stanford.edu/English/HomePage Fecha última consulta: 1-Agosto-2013 [16] http://www.evogrid.org/index.php/Main_Page Fecha última consulta: 1-Agosto-2013 Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 64 [17] http://www.scalalife.eu/ Fecha última consulta: 1-Agosto-2013 [18] http://www.spec.org/cpu2006/Docs/435.gromacs.html Fecha última consulta: 28-agosto-2013 [19] http://www.prace-ri.eu/ Blog del programa de verano (Summer of HPC): https://summerofhpc.prace-ri.eu/ Fecha última consulta: 28-agosto-2013 [20] http://www.ichec.ie/ Fecha última consulta: 28-agosto-2013 [21] Lennard-Jones, J. E. Cohesion. Proceedings of the Physical Society 1931, 43, 461-482 [22] http://www.gromacs.org/Documentation/Installation_Instructions_4.5 Fecha última consulta: 16-agosto-2013 [23] http://es.wikipedia.org/wiki/Interfaz_de_Paso_de_Mensajes Fecha última consulta: 16-agosto-2013 [24] http://es.wikipedia.org/wiki/Enlace_din%C3%A1mico Fecha última consulta: 16-agosto-2013 [25] http://www.rcsb.org/pdb/home/home.do Fecha última consulta: 21-agosto-2013 Proyecto Fin de Carrera – Ingeniería Informática Mª Astón Serrano Gracia 65