Ejercicio para reproducir la curva de la figura 2 del JCP "Constant force approach ..."

Contenido de la carpeta y descripción:
run - escript para generar los archivos necesarios y ejecutar los diferentes programas del gromacs
EM.mdp - archivo de entrada al grompp que establece los parámetros de la simulación (preproceso de archivos de entrada de gromacs). En este caso para realizar una rápida minimización de energía (para evitar solapes de la configuración inicial)
T_0.mdp - archivo de entrada al grompp para realizar la corrida formal con la mínima temperatura. 
table_l100.xvg - tabla de energía potencial con pendiente alpha=100 y lambda = 1.5 (no la usamos pero por si acaso)
table_l400.xvg - tabla de energía potencial con pendiente alpha=400 y lambda = 1.5
topol.top - topología del sistema (tipos de partícula, número de moléculas, etc.)
Tr.gro - posición de una partícula tipo Tr (definida en topol.top) que se usa cómo referencia.
temperaturas.ods - Archivo auxiliar tipo planilla de cálculo donde vemos las temperaturas a correr para reproducir la curva 2 del artículo.
vdwradii.dat - Aquí simplemente se agrega un radio para la especie Tr que se utiliza para la generación de la configuración inicial únicamente.
2013-06-JCP.pdf - Artículo "Constant force approach ..."

Lectura relacionada en el manual de gromacs ...
genbox Ejecutable que genera la caja
grompp Preprocesador de archivos de entrada
mdrun Ejecutable de gromacs (construye las trayectorias de la dinámica, es decir, la simulación propiamente dicha)
editconf Editor de la configuración (la usamos para agrandar la caja y dejar vacío a los lados)


archivo run paso a paso.
Este archivo es el que vamos a utilizar para correr las diferentes instrucciones a gromacs.
# -> significa que es una línea comentada (al principio está todo comentado).

------------------------------
Bloque 1 (generar una configuración inicial 
Descomentar la línea 
#gmx solvate -cs Tr.gro -maxsol 1200 -box 10 10 20 -o conf.gro
y correr run para obtener conf.gro.
Luego correr 
#genbox -cp Tr.gro -ci Tr.gro -nmol 1199 -box 10 10 20 -o conf.gro
guardar el archivo run (desde ahora descomentar significará quitar "#" de la línea y guardar)
correr ./run en una terminal sobre el directorio en cuestión
si no funciona hacer 
chmod +x run (para hacerlo ejecutable)

al final del proceso debe decir algo cómo 

Output configuration contains 1200 atoms in 1200 residues
Volume                 :        2000 (nm^3)
Density                :     11.9668 (g/l)
Number of SOL molecules:      0 

y luego una línea con alguna frace célebre (a mi me tocó) gcq#149: "I Ripped the Cord Right Out Of the Phone" (Capt. Beefheart)
hay otras más interesantes

se habrá generado conf.gro  (un archivo que contiene la ubicación de las moléculas)
Aquí el programa genbox "solvata" a una molécula de Tr con 1199 moléculas del solvente Tr (está pensado para solvatar una proteina con agua, por ej.).

Comentamos la línea que acabamos de usar y descomentamos la siguiente: grompp -f EM.mdp
Corremos ./run (guardarlo antes!)
al final obtenemos ...
Number of degrees of freedom in T-Coupling group System is 3597.00 (lo que es 3*1200 -3, por lo que es correcto)
se genera el archivo topol.tpr (contiene toda la información para que mdrun pueda correr).

Corremos las siguiente línea: nice -n +19 mdrun -v -table table_400.xvg
Aquí hace una minimización de la energía. nice -n +19 significa que correra la instrucción mdrun con la prioridad más baja (19). -v es que diga lo que está pasando.

cp confout.gro conf.gro -> para que la configuración que jale el próximo grompp sea la minimizada (la anterior la tiramos).

gmx editconf -f conf.gro -box 10 10 40 -> agranda la caja y genera out.gro
cp out.gro conf.gro -> me quedo con la configuración cambiada.
------------------------------

------------------------------
Bloque 2.
Comentamos todo y descomentamos lo siguiente:

j=57.7304
for i in {1..11..1}
do
cp T_0.mdp T_$i.mdp
j=`echo "$j*1.0235394508" | bc`
sed -i 's/ref-t                    = 57.7304/ref-t                    = '$j'/g' T_$i.mdp
sed -i 's/gen-temp                 = 57.7304/gen-temp                 = '$j'/g' T_$i.mdp
done

corremos ./run.
Esto nada tiene que ver con gromacs. Aquí corremos un escript para generar 11 copias de T_0.mdp (con su nombre cambiado de acuerdo a i) y que luego editamos para cambiarle la temperatura 57.7304 a su valor correspondiente j. Una vez que se corre el escript es conveniente editar T_11.mdp para verificar que la temperatura de este archivo corresponda a la temperatura más alta que queremos correr (ver el archivo temperaturas.ods para la lista de todas las temperaturas). Hay que notar que la temperatura de gromacs no coincide con la temperatura adimensional (ver el manual en la parte de unidades). 
------------------------------


------------------------------
Bloque 3.
Comentamos todo y descomentamos lo siguiente:

for i in {0..11..1}
do
grompp -f T_$i.mdp -c conf.gro -o triangle_$i.tpr
nice -n +19 mdrun -s triangle_$i.tpr -c conf_$i.gro -e ener$i.edr -o traj$i.trr -table table_400 -v
done

Este es otro escript donde se realizarán, de forma serial, 12 simulaciones. 
Aquí grompp trabaja con los diferentes T_$i.mdp que creamos (idénticos con excepción de la temperatura) y genera los correspondientes triangle_$i.tpr.
Estos archivos serán la entrada a mdrun para generar configuraciones finales conf_$i.gro, trayectorias traj$i.trr, y archivos de energía ener$i.edr. 
Las corridas aquí son cortas (10000 pasos). Esto se estipula en la línea "nsteps = 10000   ;40000000" del archivo T_0.mdp (a la derecha de los puntos y comas son comentarios)


------------------------------
Bloque 4.
Comentamos todo y descomentamos lo siguiente:

for i in {0..11..1}
do
grompp -f T_$i.mdp -c conf_$i.gro -o triangle_$i.tpr
done

Aquí se preprocesan los diferentes T_$i.mdp y conf_$i.gro. Estos últimos, generados en el bloque 3, son las configuraciones equilibradas a las diferentes temperaturas.
Por último se corren todas en paralelo permitiendo intercambios de las configuraciones con la siguiente línea.

nice -n +19 mpirun -np 12 gmx mdrun_mpi -s triangle_.tpr -multi 12 -replex 500 -table table_400.xvg &

luego le das enter, y top. Verás que están corriendo 12 procesos. Debes notar que la corrida será muy corta y no vas a poder reproducir con esto los datos de la figura 2. Lo que se debe hacer es, luego de generados los conf_$i.gro, cambiar el archivo T_0.mdp (y todos los T_$i.mdp mediante el Bloque 2) para darle más pasos. Esto se haría antes de correr el Bloque 4. Nosotros utilizamos los 40000000 que vienen luego del ";" en la líena que define a nsteps. Para una prueba puedes correr 10 veces menos y obtendrás algo parecido. 
Yo acabo de seguir los pasos para estar seguro de que aquí funciona. Por último, para analizar los datos utilizarás las líneas 
gmx g_energy -f ener0.edr (evolución de la energía, presión, tensión superficial, etc.) donde en el ejemplo se analiza la temperatura más baja y gmx g_density -f traj0.xtc -sl 12 (perfil de densidad, -sl número de puntos del perfil)

  























