EjemploGromacs+4.6.4+GPUs

Ejemplo de una simulación en Gromacs 4.6.4 para GPUs

Usando las siguientes instrucciones se puede llevar a cabo la simulación de un sistema usando Gromacs 4.6.4

Más ejemplos se pueden hallar en: http://www.gromacs.org/Documentation/Tutorials

Se ejemplifica el mecanismo usando la proteína BAX (PDB 1F16) pero se puede partir de otro archivo en formato PDB.

Para obtener el archivo PDB desde el cluster es posible descargarlo de internet con el siguiente comando (sustituyendo el código pdb, para Bax es 1F16):


 * wget http://www.rcsb.org/pdb/files/1F16.pdb **

El comando es una sola línea, por si no se ve completa.

De manera un poco informal, a los comandos utilizados siguen algunos comentarios y advertencias que pueden ser útiles.

Pasos a seguir

1.- Construcción de la topología y la estructura a partir del archivo PDB. 2.- Edición de la caja del sistema. 3.- Generar las moléculas de solvente. 4.- Pre-procesar la equilibración del sistema. 5.- Equilibrar la carga del sistema añadiendo iones. 6.- Pre-procesar la equilibración del sistema con iones añadidos. 7.- Correr la equilibración del sistema. 8.- Pre-procesar una simulación con restricción en las posiciones. 9.- Correr una simulación con restricción en las posiciones 10.- Pre-procesar la simulación para el tiempo deseado. 11.- Correr la simulación.

Archivos necesarios:

1F16.pdb Archivo PDB de la proteína BAX em.mdp Archivo de configuración de la equilibración pr.mdp Archivo de configuración de la simulación restringida full.mdp Archivo de configuración de la simulación

El único proceso que es necesario correr en paralelo usando varios CPUs es el último (correr la simulación). El paso 9 (correr una simulación con restricción en las posiciones) también puede llegar a requerir un tiempo importante de cómputo y también puede ser ejecutado en paralelo.

Para usar más de un nodo en una ejecución es necesario tener acceso por **ssh **sin dar password, al final de este documento se incluyen instrucciones para conseguir esto.

En resumen, los comandos a ejecutar son los siguientes:

pdb2gmx -ignh -f 1F16.pdb -p bax.top -o bax.gro

editconf -f bax.gro -o baxedit.gro -d 1.0

genbox -cp baxedit.gro -cs -p bax.top -o before_em.gro

grompp -v -f em.mdp -c before_em.gro -o em.tpr -p bax.top

genion -np 3 -s em.tpr -p bax.top -o baxion.gro

grompp -v -f em.mdp -c baxion.gro -o em.tpr -p bax.top

mdrun -v -s em.tpr -o em.trr -c after_emion.gro -g emlog.log

grompp -f pr.mdp -o pr.tpr -c after_emion.gro -r after_emion.gro -p bax.top

mdrun -v -s pr.tpr -e pr.edr -o pr.trr -c after_pr.gro -g 1> out_pr.txt 2> err_pr.txt &

grompp -v -f full.mdp -o full.tpr -c after_pr.gro -p bax.top

nohup mdrun -v -gpu_id 0 -s full.tpr -e full.edr -o full.trr -c after_full.gro 1>outfull.txt 2>errfull.txt &

1.- Construcción de la topología y la estructura a partir del archivo PDB.

Conectarse a una máquina con gromacs 4.6.4

Ejecutar:

pdb2gmx -ignh -f 1F16.pdb -p bax.top -o bax.gro

Pide seleccionar un campo de fuerza, hemos estado usando el OPLS-AA, que suele ser bueno para simular proteínas.

se generan los archivos bax.gro bax.top posre.itp

La opción **-ignh ** resuelve un problema que parece ser causado por un exceso de átomos de hidrógeno, ignorándolos, al parecer las librerías de topología tienen definidos estados de protonación que pueden conflictuar con el pdb.

2.- Edición de la caja del sistema.

La proteína se encuentra en una caja ajustada a su alrededor. Antes de rodearla de agua conviene editar esta caja de manera que haya un margen entre los extremos de la proteína y los de la caja. En gromacs se usan nanómetros para las distancias, 1 nm = 10 Angstroms.

editconf -f bax.gro -o baxedit.gro -d 1.0

se genera el archivo

baxedit.gro

3.- Generar las moléculas de solvente.

genbox -cp baxedit.gro -cs -p bax.top -o before_em.gro

se actualiza el archivo bax.top se respalda el archivo se genera el archivo before_em.gro
 * 1) bax.top.1#

La topología bax.top nos servirá para todo lo siguiente, el archivo before_em.gro puede ser abierto en VMD para ver a la proteína con las moléculas de agua que se añadieron.

4.- Pre-procesar la equilibración del sistema.

Hace falta el archivo de configuración em.mdp, para este ejemplo utilizo los mismos que vienen en el tutorial de gromacs (dentro de un nodo con gromacs: /usr/share/gromacs/tutor/ speptide/ ). En el archivo em.mdp fue necesario hacer una modificación, en la línea nsteps = 1000 el valor original era de 100.

En el comando se incluyen el archivo de opciones (-f em.mdp), la estructura de referencia (-c before_em.gro), el nombre del archivo a generar (-o em.tpr) y la topología de referencia (-p bax.top). El comando completo es:

grompp -v -f em.mdp -c before_em.gro -o em.tpr -p bax.top

se generan los archivos em.tpr mdout.mdp

Los archivos tpr son usados como entrada para la simulación.

Dentro de la salida de esta ejecución de g_grompp, hay una nota que nos dice que nuestro sistema no tiene carga neutra:


 * <span style="font-family: 'courier new',monospace;">NOTE 1 [file bax.top, line 27833]: **
 * <span style="font-family: 'courier new',monospace;">System has non-zero total charge: -3.000001e+00 **

Es recomendable realizar la simulación con carga neutra, eso lo conseguiremos en el siguiente paso.

5.- Equilibrar la carga del sistema añadiendo iones.

Usando el archivo em.tpr y la topología bax.top vamos a sustituir moléculas de solvente para añadir iones, de manera que se equilibre la carga del sistema.

El comando indica que se añadirán 3 iones positivos (-np 3), que serán de tipo Na, el binario de opciones de la simulación (-s em.tpr), la topología de referencia (-p bax.top) y el archivo de estructura que se va a generar incluyendo los iones (-o baxion.gro). El comando completo es:

<span style="background-color: #ffffff; color: #222222; font-family: arial,sans-serif; font-size: 12.666666984558105px;">genion -np 3 -s em.tpr -p bax.top -o baxion.gro

Nos pide seleccionar un grupo para sustituir elementos de ese grupo por los iones, seleccionamos:

Group 12 ( SOL) has 21585 elements

se generan los archivos

baxion.gro bax.top genion.log

Si vemos las últimas líneas de bax.top se nota la sustitución realizada.

tail bax.top [ system ] PROTEIN (APOPTOSIS REGULATOR BAX, MEMBRANE; 3 ISOFORM ALPHA) in water
 * Name

[ molecules ] Protein_A 1 SOL 7192 Na+ 3 Cl- 0
 * Compound #mols

También podemos buscar los iones dentro del archivo baxion.gro, el cual puede ser visualizado con VMD.


 * <span style="font-family: 'courier new',monospace;">[root@localhost bax2]# grep Na baxion.gro **
 * <span style="font-family: 'courier new',monospace;">7385Na+ Na24554 1.650 3.582 1.394 **
 * <span style="font-family: 'courier new',monospace;">7386Na+ Na24555 2.093 5.437 0.208 **
 * <span style="font-family: 'courier new',monospace;">7387Na+ Na24556 6.192 0.964 1.985 **

6.- Pre-procesar la equilibración del sistema con iones añadidos.

Ahora se debe generar un nuevo archivo tpr para simular el sistema con iones.

<span style="background-color: #ffffff; color: #222222; font-family: arial,sans-serif; font-size: 12.666666984558105px;">grompp -v -f em.mdp -c baxion.gro -o em.tpr -p bax.top

Se actualiza el archivo em.tpr se respalda el archivo se actualiza el archivo mdout.mdp
 * 1) em.tpr.1#

7.- Correr la equilibración del sistema.

En este paso se corre una dinámica molecular que trata de acomodar a los átomos de manera que no se traslapen y que la mayor fuerza en el sistema no sea mayor que Fmax.

<span style="background-color: #ffffff; color: #222222; font-family: arial,sans-serif; font-size: 12.666666984558105px;">mdrun -v -s em.tpr -o em.trr -c after_emion.gro -g emlog.log

Al final de la ejecución se debe ver algo como lo siguiente.


 * <span style="font-family: 'courier new',monospace;">Steepest Descents converged to Fmax < 2000 in 218 steps **
 * <span style="font-family: 'courier new',monospace;">Potential Energy = -3.7150081e+05 **
 * <span style="font-family: 'courier new',monospace;">Maximum force = 1.4393900e+03 on atom 2169 **
 * <span style="font-family: 'courier new',monospace;">Norm of force = 7.4190405e+03 **

Como se ve, el objetivo se alcanza en 218 pasos, de ahi la necesidad del cambio en el parámetro nsteps mencionado en el paso 4.

Se generan los archivos after_emion.gro emlog.log em.trr ener.edr

8.- Pre-procesar una simulación con restricción en las posiciones.

Para este paso usamos el archivo pr.mdp del tutorial, que incluye una opción para restringir todos los enlaces: Esto permite que los átomos se acomoden de manera que no haya energías artificialmente elevadas.
 * <span style="font-family: 'courier new',monospace;">constraints = all-bonds **

El comando es:

<span style="background-color: #ffffff; color: #222222; font-family: arial,sans-serif; font-size: 12.666666984558105px;">grompp -f pr.mdp -o pr.tpr -c after_emion.gro -r after_emion.gro -p bax.top

se generan los archivos mdout.mdp pr.tpr

9.- Correr una simulación con restricción en las posiciones

<span style="background-color: #ffffff; color: #222222; font-family: arial,sans-serif; font-size: 12.666666984558105px;">mdrun -v -s pr.tpr -e pr.edr -o pr.trr -c after_pr.gro -g 1> out_pr.txt 2> err_pr.txt &

El archivo que encuentro más útil es err_pr.txt, se puede ver cómo va avanzando con el comando

tail -f err_pr.txt

El ampersand final (&) sirve para que se ejecute en segundo plano, de manera que los regresa al prompt y ya pueden cerrar su ventana de terminal, el proceso está corriendo. Este procesamiento se tarda algunos minutos, pero la simulación completa (pasos siguientes) se tarda generalmente varias horas. 10.- Pre-procesar la simulación para el tiempo deseado.

<span style="background-color: #ffffff; color: #222222; font-family: arial,sans-serif; font-size: 12.666666984558105px;">grompp -v -f full.mdp -o full.tpr -c after_pr.gro -p bax.top

Hay que cambiar full.mdp para ajustar el tiempo que se desea simular.


 * <span style="font-family: 'courier new',monospace;">nsteps = 1000000 ; total 2000 ps = 2 ns. **

11.- Correr la simulación.

<span style="background-color: #ffffff; color: #222222; font-family: arial,sans-serif; font-size: 12.666666984558105px;">nohup mdrun -v -gpu_id 0 -s full.tpr -e full.edr -o full.trr -c after_full.gro 1>outfull.txt 2>errfull.txt &

Todo igual que en la simulación con rescricciones. La simulación de 2ns con Bax en agua tardó un par de horas. El archivo trr es el más grande, full.trr pesa 2.2 GB

Eso es todo! Sobre el trr ya se puede visualizar la simulación, aunque es lento por ser tan pesado el archivo. También con el archivo trr podemos calcular el RMSD, SAS y hacer otros análisis.

Para continuar una simulación, digamos para añadir otro nanosegundo a partir de donde terminó una simulación previa, la forma del archivo mdp puede ser la misma, solo modificando el tiempo deseado. También se pueden modificar la temperatura y otros parámetros. Para tener en cuenta la simulación precedente, se debe considerar el archivo de estado state.cpt.

Modificación en el mdp para 1ns (otro_ns.mpd):


 * <span style="font-family: 'courier new',monospace;">nsteps = 500000 ; total 1000 ps = 1 ns. **

Aunque no es indispensable, se puede asignar un tiempo de inicio (por default inicia en t=0), por ejemplo, para empezar en 2000 ps :

tinit = 2000

Para preprocesar la simulación:


 * grompp -v -f otro_ns.mdp -o otro_ns.tpr -c after_full.gro -p bax.top** -t state.cpt

Se genera el archivo otro_ns.tpr, y se puede correr de igual manera que la simulación precedente:

nohup **<span style="font-family: 'courier new',monospace;">/usr/lib64/openmpi/1.2.7-gcc/bin/mpirun -np 16 --host node05,node06 g_mdrun_mpi -v -s otro_ns.tpr -e otro_ns.edr -o otro_ns.trr -c after_otro_ns.gro -g otro_ns.log >& otro_ns.job 1> out_otro_ns.txt 2> err_otro_ns.txt & **

Para obtener el rms y el radio de giro se necesitan el archivo de la trayectoria (trr) y el de entrada para la simulación (tpr), también se puede asignar un nombre al archivo de salida:

g_rms -f full.trr -s full.tpr -o full_rms.xvg

g_gyrate -f full.trr -s full.tpr -o full_gyrate.xvg

Para obtener una visualización en la cual la proteína se vea completa y en el centro de la caja, conviene primero pasar a una representación más compacta (el ejemplo conserva un solo decimal de precisión y solamente uno de cada 200 frames):

g_trjconv -f full.trr -o full.xtc -s full.tpr -pbc nojump -ndec 1 -skip 200

Corregir moléculas "discontinuas", centrar permite elegir un grupo, por ejemplo los C-alfa.

g_trjconv -f full.xtc -o fullw.xtc -s full.tpr -pbc whole -center -boxcenter rect

Colocar todas las moléculas en la caja.

g_trjconv -f fullw.xtc -o fullwm.xtc -s full.tpr -pbc mol -ur compact