EjemploGromacs4

Ejemplo de una simulación en Gromacs 4.0.4-1

Usando las siguientes instrucciones se puede llevar a cabo la simulación de un sistema usando Gromacs 4.0.4-1.

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.

Actualmente en el cluster Gromacs 4.0.4-1 solo está instalado en los nodos 05 y 06, para ejecutar los comandos es necesario estar conectado a uno de ellos. Por ejemplo, para conectarse al node05:


 * ssh node05 **

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

Los archivos de configuración mdp se modificarán a partir de los que vienen en el tutorial de gromacs. Está instalado en el nodo 5 en:


 * /usr/share/gromacs/tutor/speptide/**

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:

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

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

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

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

g_genion -np 3 -pname Na+ -nname Cl- -s em.tpr -p bax.top -o baxion.gro

g_grompp -v -f em.mdp -c baxion.gro -o em.tpr -p bax.top -maxwarn 1

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

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

/usr/lib64/openmpi/1.2.7-gcc/bin/mpirun -np 8 g_mdrun_mpi -v -s pr.tpr -e pr.edr -o pr.trr -c after_pr.gro -g pr.log >& pr.job **1> out.txt 2> err.out **&

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

nohup **/usr/lib64/openmpi/1.2.7-gcc/bin/mpirun -np 16 --host node05,node06 g_mdrun_mpi -v -s full.tpr -e full.edr -o full.trr -c after_full.gro -g full.log >& full.job 1>outfull.txt 2>errfull.txt & **

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

Conectarse a un nodo con gromacs 4.0.4

Ejecutar:


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

Pide seleccionar un campo de fuerza, hemos estado usando el 5, 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.


 * g_editconf -f bax.gro -o baxedit.gro -d 1.0 **

se genera el archivo

baxedit.gro

3.- Generar las moléculas de solvente.

g_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:


 * g_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:


 * NOTE 1 [file bax.top, line 27833]: **
 *  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.

Los valores por default Na y Cl no funcionan en OPLS, hay que usar Na+ y Cl-

El comando indica que se añadirán 3 iones positivos (-np 3), que serán de tipo Na+ (-pname Na+), no se añaden iones negativos, pero como referencia se pone el nombre que se usaría para iones de Cl- en ese caso (-nname Cl-), 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:


 * g_genion -np 3 -pname Na+ -nname Cl- -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.


 * [root@localhost bax2]# grep Na baxion.gro **
 *  7385Na+ Na24554 1.650 3.582 1.394 **
 *  7386Na+ Na24555 2.093 5.437 0.208 **
 *  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="font-family: courier new,monospace;">g_grompp -v -f em.mdp -c baxion.gro -o em.tpr -p bax.top -maxwarn 1 **


 * <span style="font-family: courier new,monospace;">Warning: atom name 38645 in bax.top and baxion.gro does not match (NA - Na) **
 * <span style="font-family: courier new,monospace;"> Warning: atom name 38646 in bax.top and baxion.gro does not match (NA - Na) **
 * <span style="font-family: courier new,monospace;"> Warning: atom name 38647 in bax.top and baxion.gro does not match (NA - Na) **


 * <span style="font-family: courier new,monospace;"> WARNING 1 [file bax.top, line 27835]: **
 * <span style="font-family: courier new,monospace;"> 3 non-matching atom names **
 * <span style="font-family: courier new,monospace;"> atom names from bax.top will be used **
 * <span style="font-family: courier new,monospace;"> atom names from baxion.gro will be ignored **

No aparece NA ni en bax.top ni en baxion.gro... pero si en em.tpr, no le he podido quitar ese Warning pero al parecer no afecta. Sin embargo, para generar la salida a pesar del warning es necesario incluir la opción -maxwarn 1.

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="font-family: courier new,monospace;">g_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 **

Hay que cambiar pr.mdp para añadir Na+


 * <span style="font-family: courier new,monospace;">tc-grps = Protein SOL Na+ **
 * <span style="font-family: courier new,monospace;"> tau_t = 0.1 0.1 0.1 **
 * <span style="font-family: courier new,monospace;"> ref_t = 300 300 300 **
 * <span style="font-family: courier new,monospace;"> ; Energy monitoring **
 * <span style="font-family: courier new,monospace;"> energygrps = Protein SOL Na+ **

El comando es:


 * g_grompp -v -f pr.mdp -o pr.tpr -c after_emion.gro -r after_emion.gro -p bax.top**

Las notas que se generan sería interesante comentarlas un día de estos.


 * <span style="font-family: courier new,monospace;">NOTE 1 [file pr.mdp, line unknown]: **
 * <span style="font-family: courier new,monospace;"> The Berendsen thermostat does not generate the correct kinetic energy **
 * <span style="font-family: courier new,monospace;"> distribution. You might want to consider using the V-rescale thermostat. **


 * <span style="font-family: courier new,monospace;">NOTE 2 [file pr.mdp, line unknown]: **
 * <span style="font-family: courier new,monospace;"> You are using a plain Coulomb cut-off, which might produce artifacts. **
 * <span style="font-family: courier new,monospace;"> You might want to consider using PME electrostatics **.

se generan los archivos mdout.mdp pr.tpr

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

Este proceso puede llevar algo de tiempo, conviene correrlo en varios procesadores.

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

En 8 procesadores
 * <span style="font-family: courier new,monospace;">nohup /usr/lib64/openmpi/1.2.7-gcc/ ****<span style="font-family: courier new,monospace;">bin/mpirun -np 8 g_mdrun_mpi -v -s pr.tpr -e pr.edr -o pr.trr -c after_pr.gro -g prlog.log >& pr.job 1> out.txt 2> err.out ** &

En el archivo err.out aparece una nota interesante, y también tiene medidas del desempeño del cluster para este cálculo. Entre proteína y agua hay 38647 átomos y se consiguen 21 ns/día en 16 procesadores.

NOTE: 5 % of the run time was spent communicating energies, you might want to use the -nosum option of mdrun

Parallel run - timing based on wallclock.

NODE (s) Real (s) (%) Time: 41.000 41.000 100.0 (Mnbf/s) (GFlops) (ns/day) (hour/ns) Performance: 981.781 29.782 21.077 1.139

Algunos comentarios sobre la instrucción utilizada. Esta versión de Gromacs usa openMPI, que es una versión distinta al MPI que se instaló en el cluster originalmente. Para usar mpirun hay que llamarlo con su ruta completa:


 * <span style="font-family: courier new,monospace;">/usr/lib64/openmpi/1.2.7-gcc/ ****<span style="font-family: courier new,monospace;">bin/mpirun **

Luego sigue el número de procesadores, que usando completos los dos nodos son 16, y especificamos los nodos donde queremos que corra, en este caso el 5 y el 6.

-np 16 --host node05,node06

A continuación se llama a mdrun, que en esta versión se le antepone g_ y la versión en paralelo termina en _mpi. Lo demás es lo que se usa en Gromacs desde versiones anteriores. El binario con las opciones de simulación -s pr.tpr, el archivo donde se guardará el monitoreo de la energía -e pr.edr, el archivo donde se guardará la trayectoria pr.trr, el archivo para la estructura alcanzada al final de la simulación -c after_pr.gro, y finalmente los archivos de bitácora que nos permiten ver cómo va la simulación. El que encuentro más útil es err.out, se puede ver cómo va avanzando con el comando

tail -f err.out

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. El comando completo es:


 * <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 pr.tpr -e pr.edr -o pr.trr -c after_pr.gro -g pr.log >& pr.job 1> out.txt 2> err.out &

El **nohup** al principio ayuda a que sea posible cerrar la terminal sin que el proceso se detenga.

10.- Pre-procesar la simulación para el tiempo deseado.


 * g_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 y añadir Na+


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


 * <span style="font-family: courier new,monospace;">tc-grps = Protein SOL Na+ **
 * <span style="font-family: courier new,monospace;"> tau_t = 0.1 0.1 0.1 **
 * <span style="font-family: courier new,monospace;"> ref_t = 300 300 300 **
 * <span style="font-family: courier new,monospace;"> ; Energy monitoring **
 * <span style="font-family: courier new,monospace;"> energygrps = Protein SOL Na+ **

11.- Correr la simulación.


 * <span style="font-family: courier new,monospace;">nohup /usr/lib64/openmpi/1.2.7-gcc/bin/mpirun -np 16 --host node05,node06 g_mdrun_mpi -v -s full.tpr -e full.edr -o full.trr -c after_full.gro -g full.log >& full.job 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 deben considerar los archivos de trayectoria (ttr) y de energía (edr).

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:


 * g_grompp -v -f otro_ns.mdp -o otro_ns.tpr -c after_full.gro -p bax.top** -t full.trr -e full.edr

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

Instrucciones para acceder a otros nodos sin necesidad de indicar la clave de usuario cada vez que se conecta con ssh.

1- Conectarse al cluster 2- Ejecutar el comando


 * <span style="font-family: courier new,monospace;">ssh-keygen -b 4096 -t rsa **

pide un nombre de archivo y luego una frase, dar ENTER simplemente para aceptar valores por default.

Eso genera **<span style="font-family: courier new,monospace;">$HOME/.ssh/id_rsa **<span style="font-family: arial,helvetica,sans-serif;">y **<span style="font-family: courier new,monospace;"> $HOME/.ssh/id_rsa.pub ** que son las claves privada y pública para el sistema RSA.

3- Añadir la clave generada al archivo authorized_keys y asegurar que los permisos de este archivo sean 600.

Esto se consigue con los siguientes tres comandos


 * <span style="font-family: courier new,monospace;">cd .ssh/ **


 * <span style="font-family: courier new,monospace;"> cat id_rsa.pub >> authorized_keys **


 * <span style="font-family: courier new,monospace;"> chmod 600 authorized_keys **

Eso es todo. El proceso suele incluir algunos pasos más, pero como el /home es el mismo en todos los nodos basta con estos pasos.

Como prueba, se puede intentar una conexión a uno de los nodos, por ejemplo:


 * <span style="font-family: courier new,monospace;">ssh node05 **

Se debe conectar sin pedir password.

Referencia: []