<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 FIMDA
Script que realiza el análisis de dinámica para una trayectoria con 2 puentes di sulfuro.
El presente script genera los archivos resultantes del análisis de una trayectoria haciendo uso
Step1:  Ruta de la trayectoria
Escribir después de la diagonal, la ruta de la trayectoria seleccionada con el rmsd más bajo.
Step2:  Convirtiendo la trayectoria DCD -> XTC
Los siguientes comandos convierten la trayectoria DCD contenida en la carpeta seleccionada a formato de XTC
Crear la nueva ruta para enviar las trayectorias convertidas
Step3:  Realizando la conversión de la trayectoria
Step4:  Cargando la nueva trayectoria en VMD para su revisión
Step5:  Calculando el RMSD con Gromacs 5
El siguiente script obtiene el RMSD de la trayectoria haciendo uso de Gromacs 5
Creando la carpeta de RMSD
Step6:  Calculando el RMSD con la opción 3 'C-Alpha'
Select group for least squares fit
Group     3 (        C-alpha)
Select a group
Step7:  Creando archivo rmsd.dat para su visualización en XMGRACE
Se genera el archivo de salida rmsd.dat, éste se deberá visualizar con Xmgrace para guardarlo en formato PNG.
Step8:  Creando el archivo rmsd_residue.dat para visualizar con xmgrace
Se crea el archivo rmsd_residue.dat formateado para su visualización en Xmgrace, en donde se deberá guardar como imagen PNG.
Step9:  Creando archivo rmsd.dat para su visualización en Matplotlib
Se genera el gráfico de salida para matplotlib
Step10:  Creando archivo rmsd_residue.dat para su visualización en Matplotlib
Se genera el gráfico de salida para matplotlib
Step11:  RMSF
Se crea una carpeta RMSF para guardar los archivos generados.
Step12:  Calculando el RMSF con la opción 3 'C-Alpha'
Step13:  Creando archivo rmsf.dat para su visualización en XMGRACE
Se genera el archivo de salida rmsf.dat, éste se deberá visualizar con Xmgrace para guardarlo en formato PNG.
Step14:  Creando archivo rmsf.dat para su visualización en Matplotlib
Se genera el gráfico de salida para matplotlib
Step15:  B-factors
Generando archivo para visualizarlo con XMGRACE
Step16:  Generando archivo para visualizar con Matplotlib
Step17:  Secondary Structure
Se crea la carpeta para cálculo de la estructura
Step18:  Calculando la estructura secundaria
Se necesita contar con el programa dssp en la ruta /usr/local/bin, el cual se enlaza con Gromacs 5
Step19:  R-GYRATE
Se crea una carpeta rgiro para guardar los archivos generados.
Step20:  Calculando el r-gyro con la opción (3) - C-alpha
Se calcula para los carbonos alfa.
Step21:  Generando el archivo rgyro.dat para su análisis con XMGRACE
Step22:  Ploteando el archivo gyrate.xvg con matplotlib
Step23:  RMSD Helix Alfa
Para realizar este análisis se debe cargar el pdb original de la proteina que se encuentra en la carpeta 01_BUILD.
Cargarlo con VMD y dirigirse al Menú EXTENSIONS -> ANALYSIS -> SEQUENCE VIEWER, en la cual se tomará el rango de átomos del campo Struct (H), el cual se proporcionará de la forma "resid X1 to X2" donde X1 es primer átomo de la helix y X2 el último átomo de la helix.
Step24:  Entrada de datos
Para la entrada se deberá dar con la opción "resid X to X".
Step25:  SASA
Creando la estructura de carpeta para el cálculo
Step26:  Ejecutando el análisis de SASA con Gromacs4
Step27:  Creando el archivo sasa_residuo.dat para salida con XMGRACE
Step28:  Cargando archivo residue-sas.xvg para su visualización en Matplotlib
Se genera el gráfico de salida para matplotlib
Step29:  Creando el archivo sasa.dat para salida con XMGRACE
Step30:  Cargando archivo solven-accessible-surface.xvg para graficar con Matplotlib
Step31:  MATRIZ DE RMSD
Step32:  Cargar el scrit rmsd_matrix con vmd en la nueva trayectoria
Arrancar VMD, dirigirse al manú Extensions -> Tk Console, copiar y ejecutar la siguiente secuencia de comandos
Step33:  Graficando el archivo de salida
Step34:  Matriz de distancia mínima
Step35:  Calculando la matriz de distancia mínima
Seleccionar el backbone (opción 4)
Step36:  Generando los archivos  para visualizarlos
Step37:  Free Energy
Para el cálculo de la energía libre se requiere el valor mínimo y máximo del RMSD y del radio de gyro, así como el valor de la temperatura a la cual se realizó la simulación. Estos datos son de entrada para el script del cálculo del mismo.
Step38:  Calculando el rmsd y el r-gyro para obtener el mínimo y máximo de cada uno de ellos.
Step39:  Escribiendo script a /tmp para utilizar en el cálculo
Step40:  Realizando los cálculos para la Free Energy
Step41:  Ploteando con GNUplot
Step42:  PCA
Step43:  Una vez calculada la matriz el eigenvalues y eigenvectors sirven de entrada para generar el pca.
El siguiente comando representa el movimiento del primer y segundo eigenvector.
Step44:  Ploteando el archivo con gnuplot
Step45:  Análisis de puentes di sulfuro
Este aplica para 2 puente, para lo cual se utiliza el software HTMD.
Step46:  Creando la ruta
Ruta para el análisis de los datos.
Step47:  Cargando de los puentes di sulfuro
Para este análisis se deberá revisar el archivo psf_charmm.tcl de la carpeta 01_BUILD, en el cual se tiene la definición de los puentes como la siguiente
Step48:  Obteniendo los índices de los puentes
Step49:  Ordenando los puentes  de la forma ['N', 'CA', 'CB', 'SG', 'SG', 'CB', 'CA', 'N']
Step50:  Creando los archivos tcl para el cálculo del RMSD de los puentes
Se crean los archivos de salida en formato tcl.
Step51:  Ejecutando los archivos rmsd en tcl con vmd
Ejecutando los archivos en VMD
Step52:  Generando los gráficos RMSD en matplotlib
Step53:  FREE ENERGY  DIHEDRAL INTRAMOLECULAR
Se calculan las distancias de los ángulos diedros para el cálculo de la free energy intramolecular
Step54:  Creación de los archivos tcl para el cálculo de los ángulos diedros
Step55:  Ejecutando los archivos de los ángulos diedros tcl generados con VMD
Step56:  Calculando la Free Energy Intramolecular para el Puente 1
Step57:  Ploteando con GNUPLOT el Puente 1
Step58:  Calculando la Free Energy intramolecular para el Puente 2
Step59:  Ploteando con GNUPLOT el Puente 2
Step60:  Free Energy Intermolecular
Step61:  Ploteando la Free Energy Intermolecular puentes DB1 y DB2
Step62:  Calcular los histogramas de los diedros
Step63:  Ángulos de Enlace de los puentes Intermolecular
Step64:  Copiando el archivo de generación de FES
Step65:  Generando los archivos Tcl para el cálculo de los ángulos.
Step66:  Ejecutando los archivos tcl generados con VMD
Step67:  Calculando la Free Energy de los Bonds de los puentes
Step68:  Ploteando la Free Energy de los ángulos con gnuplot
Step69:  Calculando los histogramas de los bonds
Step70:  Generación de clusters
Crear la nueva ruta para calcular los clusters
Step71:  Calculando los clusters con la opción (1= Protein)
Step72:  Cargando los clusters para su visualización en VMD
Se cargan los clusters en VMD y se guardan sus coordenadas para cada uno de ellos haciendo uso de VMD
Step73:  colorByRMSF
Creando la carpeta para salida de datos
Step74:  Copiando el archivo a la carpeta de datos
Step75:  Calculando el RMSF para el análisis de la proteína con la opción (1) Protein
Step76:  Cargar el scrit colorByRMSF.vmd en VMD
Arrancar VMD, dirigirse al menú Extensions -> Tk Console, copiar y ejecutar la siguiente secuencia de comandos en el cual pondremos los valores del Mínimo_RMSF y Máximo_RMSF calculado en la celda anterior
Step77:  Graficando B-Factors con Chimera
Step78:  Adecuando archivo bfac.pdb para obtener la columna de B-factors
Step79:  Revisando la estructura del archivo generado.
Revisar que los campos se encuentren completamente alineados en la estructura de los campos.
Guardar y salir.
Step80:  Generando el archivo de Bfactors para todos los átomos FALTA ADECUAR PARA SACAR EL MAYOR POR RESIDUO
Step81:  Cargando el archivo pdb con Chimera para realizar la coloración de Bfactors
Step82:  Instrucciones para generar la imagen de B-factors
 ESTABLECER EL MODO DE VISUALIZACIÓN
1. Seleccionar del menú principal Presets -> Interactive 2 (all atoms).
2. Seleccionar del menú principal Actions -> Surface -> Show.
3. Ajustar el tamaño de la ventana principal.
4. Ajustar el tamaño y posición de la figura haciendo uso de la tecla CTRL+ Button wheel mouse.
 COLOREAR LOS B-FACTORS 
Seleccionar  Tools -> Depiction -> Render by Attribute.
Nos desplegará una ventana Render/Select by Attribute.
  1. Del campo Attribute seleccionar bfactor.
  2. En el histograma que se muestra, seleccionar la barra blanca y cambiar el color de blanco a amarillo en el campo color.
  3. Pulsar el botón Apply para visualizar los cambios de coloración.
  4. Pulsar OK para finalizar.
 FONDO BLANCO
Para aplicar el fondo blanco
Step83:  Graficando SASA
Step84:  Coloreando el SASA
Arrancar VMD.
Ventana vmdICE
Dirigirse al menú Extensions -> Analysis -> vmdICE, se presentará una ventana y se deberán cambiar los valores de los siguientes campos
Step85:  Restaurando configuración default de VMD
Step86:  Graficando el RGYRO
Step87:  Coloreando el RGYRO
Arrancar VMD, dirigirse al manú Extensions -> Tk Console, copiar y ejecutar la siguiente secuencia de comandos
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import pylab as pl
import matplotlib.patches as mpatches
import matplotlib.ticker as ticker
import os
import shutil
from IPython.display import Image
from matplotlib.ticker import FormatStrFormatter
Explanation: FIMDA
Script que realiza el análisis de dinámica para una trayectoria con 2 puentes di sulfuro.
El presente script genera los archivos resultantes del análisis de una trayectoria haciendo uso:
Gromacs 5 y Gromacs 4
Xmgrace
VMD
Chimera
catdcd
trjconv.
Se deberá contar con  la trayectoria que contenga el rmsd más bajo para realizar el análisis sobre ella.
Librerías a utilizar
End of explanation
ruta=os.getcwd()
c=input('Nombre de la trayectoria para realizar el análisis... Ejemplo: run001....')
if os.path.isdir(c):
    indir = '/'+c
    print (indir)
    ruta_old_traj=ruta+indir
    print (ruta)
    print (ruta_old_traj)
else:
    print ('La carpetac'+c+' no existe...')
#
ruta_scripts=ruta+'/scripts_fimda'
print (ruta_scripts)
if os.path.exists(ruta_scripts): 
    print ('Ruta identificada para búsqueda de scripst adicionales  ===>',ruta_scripts)
else:
        print ('La carpeta de scripst adicionales no existe, copiar en '+ruta_scripts+' ..!!!')
Explanation: Ruta de la trayectoria
Escribir después de la diagonal, la ruta de la trayectoria seleccionada con el rmsd más bajo.
End of explanation
#Verificando que exista la nueva carpeta para la conversión de trayectorias
#nuevaruta = ruta+'/'+indir+'_XTC'
nuevaruta = ruta+indir+'_Dinamica'
print ( nuevaruta )
if not os.path.exists(nuevaruta): 
    os.makedirs(nuevaruta)
    print ('Se ha creado la ruta ===>',nuevaruta)
else:
        print ("La ruta "+nuevaruta+" existe..!!!")
    
Explanation: Convirtiendo la trayectoria DCD -> XTC
Los siguientes comandos convierten la trayectoria DCD contenida en la carpeta seleccionada a formato de XTC
Crear la nueva ruta para enviar las trayectorias convertidas
End of explanation
print ('Obtenemos los archivos a convertir')
#Buscamos el archivo DCD, PDB y PSF para realizar las operaciones
for filename in os.listdir(ruta_old_traj):
    if filename.endswith('.dcd'):
        dcd_file=filename
    if filename.endswith('.psf'):
        psf_file=filename
    if filename.endswith('.pdb'):
        pdb_file=filename
print ('pdb file =>', pdb_file)
print ('psf file =>', psf_file)
print ('dcd file =>', dcd_file)
print ( 'Nos vemos a ....', ruta_old_traj )
os.chdir( ruta_old_traj )
print ('\nEjecutando CATDCD para convertir la trayectoria....')
output_catdcd=!catdcd -otype trr -o output.trr $dcd_file 
print (output_catdcd.n)
print ('\nEjecutando TRJCONV para convertir la trayectoria....')
output_trjconv=!trjconv -f output.trr -o output.xtc -timestep 20
#print (output_trjconv.n)
print ('\nBorrando archivos temporales de conversión...')
output_rm=!rm output.trr
print ('\nMoviendo los archivos de salida al directorio '+nuevaruta)
source_file=ruta_old_traj+'/output.xtc'
dest_file=nuevaruta+'/output.xtc'
shutil.move(source_file,dest_file)
print ('\Copiando el archivo ionized.pdb a '+nuevaruta)
source_file=ruta_old_traj+'/ionized.pdb'
dest_file=nuevaruta+'/ionized.pdb'
shutil.copy(source_file,dest_file)
print ('\nCopiando el archivo ionized.psf a '+nuevaruta)
source_file=ruta_old_traj+'/ionized.psf'
dest_file=nuevaruta+'/ionized.psf'
shutil.copy(source_file,dest_file)
print('\nTrayectoria convertida, regresando a '+ruta)
os.chdir( ruta )
Explanation: Realizando la conversión de la trayectoria
End of explanation
print ('Visualizando la nueva trayectoria')
file_psf=nuevaruta+'/'+psf_file
traj = nuevaruta+'/output.xtc'
!vmd $file_psf $traj
Explanation: Cargando la nueva trayectoria en VMD para su revisión
End of explanation
### Creando el directorio para el análisis del RMSD
#Verificando que exista la nueva carpeta para la conversión de trayectorias
#nuevaruta = ruta+'/'+indir+'_XTC'
ruta_rmsd = nuevaruta+'/rmsd'
print ( ruta_rmsd )
if not os.path.exists(ruta_rmsd): 
    os.makedirs(ruta_rmsd)
    print ('Se ha creado la ruta ===>',ruta_rmsd)
else:
        print ("La ruta "+ruta_rmsd+" existe..!!!")
    
print ( 'Nos vamos a ....', ruta_rmsd )
os.chdir( ruta_rmsd )
Explanation: Calculando el RMSD con Gromacs 5
El siguiente script obtiene el RMSD de la trayectoria haciendo uso de Gromacs 5
Creando la carpeta de RMSD
End of explanation
print ('Ejecutando el análisis de rmsd...')
!echo 3 3 | g_rms -f ../output.xtc -s ../ionized.pdb -a avgrp.xvg
Explanation: Calculando el RMSD con la opción 3 'C-Alpha'
Select group for least squares fit
Group     3 (        C-alpha)
Select a group: 3
Selected 3: 'C-alpha'
Select group for RMSD calculation
Group     3 (        C-alpha)
Select a group: 3
Selected 3: 'C-alpha'
End of explanation
#Inicializando vector
rmsd=[]
try:
     archivo = open( 'rmsd.xvg' )
except IOError:
     print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in archivo.readlines():
    fila = linea.strip()
    sl = fila.split()
    cadena=sl[0]
    if (not '#' in cadena) and (not '@' in cadena):
        num=float(sl[0])
        #num2=float(sl[1])
        num=num/1000
        rmsd.append(repr(num)+'\t'+sl[1]+'\n')
        i=i+1
#Escribiendo el archivo RMSD
f = open('rmsd.dat', 'w')
#f.write('@    title "RMSD" \n')
f.write('@    xaxis  label " Time (ns)" \n')
f.write('@    xaxis  label char size 1.480000\n')
f.write('@    xaxis  bar linewidth 3.0\n')
f.write('@    xaxis  ticklabel char size 1.480000\n')
f.write('@    yaxis  label " RMSD (nm)" \n')
f.write('@    yaxis  label char size 1.480000\n')
f.write('@    yaxis  bar linewidth 3.0\n')
f.write('@    yaxis  ticklabel char size 1.480000\n')
f.write('@    s0 line linewidth 1.5\n')
f.write('@TYPE xy \n')
#f.write('@ subtitle "C-alpha after lsq fit to C-alpha" \n')
f.write("".join(rmsd))
f.close()
#Cargando el archivo para visualizar en xmgrace
!xmgrace rmsd.dat
#Cargando la imagen generada en xmgrace
Image(filename='rmsd.png')
Explanation: Creando archivo rmsd.dat para su visualización en XMGRACE
Se genera el archivo de salida rmsd.dat, éste se deberá visualizar con Xmgrace para guardarlo en formato PNG.
End of explanation
#Inicializando vector
rmsd_residue=[]
try:
     archivo_rmsd = open( 'aver.xvg' )
except IOError:
     print ('No se pudo abrir el archivo o no existe·..')
i=1
for linea in archivo_rmsd.readlines():
    fila = linea.strip()
    sl = fila.split()
    cadena=sl[0]
    if (not '#' in cadena) and (not '@' in cadena):
        num=int(sl[0])
        print ('Residuo =>',num+1)
        rmsd_residue.append(repr(num+1)+'\t'+sl[1]+'\n')
        i=i+1
#Escribiendo el archivo RMSD_RESIDUE
f = open('rmsd_residue.dat', 'w')
#f.write('@    title "C-alpha" \n')
f.write('@    xaxis  label "Residue" \n')
f.write('@    xaxis  label char size 1.480000\n')
f.write('@    xaxis  bar linewidth 3.0\n')
f.write('@    xaxis  ticklabel char size 1.480000\n')
f.write('@    yaxis  label " RMSD (nm)" \n')
f.write('@    yaxis  label char size 1.480000\n')
f.write('@    yaxis  bar linewidth 3.0\n')
f.write('@    yaxis  ticklabel char size 1.480000\n')
f.write('@    s0 line linewidth 2.5\n')
f.write('@    s0 symbol 1\n')
f.write('@    s0 symbol size 1.000000\n')
f.write('@    s0 symbol color 1\n')
f.write('@    s0 symbol pattern 1\n')
f.write('@    s0 symbol fill color 2\n')
f.write('@    s0 symbol fill pattern 1\n')
f.write('@    s0 symbol linewidth 1.0\n')
f.write('@TYPE xy \n')
f.write("".join(rmsd_residue))
f.close()
          
!xmgrace rmsd_residue.dat
#Cargando la imagen generada en xmgrace
Image(filename='rmsd_residue.png')
Explanation: Creando el archivo rmsd_residue.dat para visualizar con xmgrace
Se crea el archivo rmsd_residue.dat formateado para su visualización en Xmgrace, en donde se deberá guardar como imagen PNG.
End of explanation
data_rmsd=np.loadtxt('rmsd.xvg',comments=['#', '@'])
#Engrosar marco    
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
  ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
pl.plot(data_rmsd[:,0]/1000, data_rmsd[:,1], linewidth = 2, markeredgewidth=3, color='black')
pl.xlabel("Time (ns)", fontsize = 40)
pl.ylabel('RMSD (nm)', fontsize = 40)
#pl.suptitle('RMSD', fontsize=50)
#pl.title('C-alpha after lsq fit to C-alpha', fontsize=30)
pl.xticks(fontsize=30) 
pl.yticks(fontsize=30) 
Explanation: Creando archivo rmsd.dat para su visualización en Matplotlib
Se genera el gráfico de salida para matplotlib
End of explanation
data_rmsd_res=np.loadtxt('aver.xvg',comments=['#', '@'])
    
#Engrosar marco    
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
  ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
pl.plot(data_rmsd_res[:,0]+1, data_rmsd_res[:,1],  '-o', color='black', markersize=25,
        markerfacecolor='red',markeredgecolor='black',markeredgewidth=3, linewidth = 4, )
pl.xlabel("Residue", fontsize = 40)
pl.ylabel('RMSD (nm)', fontsize = 40)
#pl.title('C-alpha', fontsize=40)
pl.xticks(fontsize=30) 
pl.yticks(fontsize=30) 
pl.xlim(0, len(data_rmsd_res[:,1]))
Explanation: Creando archivo rmsd_residue.dat para su visualización en Matplotlib
Se genera el gráfico de salida para matplotlib
End of explanation
### Creando el directorio para el análisis del RMSF
#Verificando que exista la nueva carpeta para la conversión de trayectorias
ruta_rmsf = nuevaruta+'/rmsf'
print ( ruta_rmsf )
if not os.path.exists(ruta_rmsf): 
    os.makedirs(ruta_rmsf)
    print ('Se ha creado la ruta ===>',ruta_rmsf)
else:
        print ("La ruta "+ruta_rmsf+" existe..!!!")
print ( 'Nos vamos a ....', ruta_rmsf )
os.chdir( ruta_rmsf )
Explanation: RMSF
Se crea una carpeta RMSF para guardar los archivos generados.
End of explanation
print ('Ejecutando el análisis de rmsf...')
!echo 3 | g_rmsf -f ../output.xtc -s ../ionized.pdb -oq bfac.pdb -o rmsf.xvg -res
Explanation: Calculando el RMSF con la opción 3 'C-Alpha'
End of explanation
#Inicializando vector
rmsf=[]
rmsf_x=[]
rmsf_y=[]
try:
     file_rmsf = open( 'rmsf.xvg' )
except IOError:
     print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in file_rmsf.readlines():
    fila = linea.strip()
    sl = fila.split()
    cadena=sl[0]
    if (not '#' in cadena) and (not '@' in cadena):
        print ('Residue =>',cadena)
        rmsf.append(sl[0]+'\t'+sl[1]+'\n')
        rmsf_x.append(int(sl[0]))
        rmsf_y.append(float(sl[1]))
        i=i+1
file_rmsf.close()
#Escribiendo el archivo RMSD
f = open('rmsf.dat', 'w')
#f.write('@    title "RMSF fluctuation" \n')
f.write('@    xaxis  label " Residue" \n')
f.write('@    xaxis  label char size 1.480000\n')
f.write('@    xaxis  bar linewidth 3.0\n')
f.write('@    xaxis  ticklabel char size 1.480000\n')
f.write('@    yaxis  label "RMSF (nm)" \n')
f.write('@    yaxis  label char size 1.480000\n')
f.write('@    yaxis  bar linewidth 3.0\n')
f.write('@    yaxis  ticklabel char size 1.480000\n')
f.write('@    s0 line linewidth 2.5\n')
f.write('@    s0 symbol 1\n')
f.write('@    s0 symbol size 1.000000\n')
f.write('@    s0 symbol color 1\n')
f.write('@    s0 symbol pattern 1\n')
f.write('@    s0 symbol fill color 2\n')
f.write('@    s0 symbol fill pattern 1\n')
f.write('@    s0 symbol linewidth 1.0\n')
f.write('@TYPE xy \n')
f.write("".join(rmsf))
f.close()
!xmgrace rmsf.dat
#Cargando la imagen generada en xmgrace
Image(filename='rmsf.png')
Explanation: Creando archivo rmsf.dat para su visualización en XMGRACE
Se genera el archivo de salida rmsf.dat, éste se deberá visualizar con Xmgrace para guardarlo en formato PNG.
End of explanation
data_rmsf=np.loadtxt('rmsf.xvg',comments=['#', '@'])
    
#Engrosar marco    
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
  ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
pl.plot(data_rmsf[:,0], data_rmsf[:,1],  '-o', color='black', markersize=25,
        markerfacecolor='red',markeredgecolor='black',markeredgewidth=3, linewidth = 4, )
pl.xlabel("Residue", fontsize = 40)
pl.ylabel('RMSF (nm)', fontsize = 40)
#pl.title('RMSF Fluctuation', fontsize=40)
pl.xticks(fontsize=30) 
pl.yticks(fontsize=30) 
pl.xlim(0, len(data_rmsf[:,1]))
Explanation: Creando archivo rmsf.dat para su visualización en Matplotlib
Se genera el gráfico de salida para matplotlib
End of explanation
#Inicializando vector
bfactors=[]
try:
     file_bfactor = open( 'bfac.pdb' )
except IOError:
     print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in file_bfactor.readlines():
    fila = linea.strip()
    sl = fila.split()
    if (sl[0]=='ATOM'):
        #print (sl[0])
        idresidue=fila[23:26]
        bfactor=fila[60:66]
        print (idresidue + '\t'+bfactor)
        bfactors.append(idresidue+'\t'+bfactor+'\n')
        #i=i+1
#Escribiendo el archivo BFACTOR.dat
f = open('bfactor.dat', 'w')
#f.write('@    title "B-factors" \n')
foo = 'baz "\\"'
f.write('@    xaxis  label " Residue" \n')
f.write('@    xaxis  label char size 1.480000\n')
f.write('@    xaxis  bar linewidth 3.0\n')
f.write('@    xaxis  ticklabel char size 1.480000\n')
f.write('@    yaxis  label "B-factors (' +"\\"+'cE'+"\\"+'C)"\n')
f.write('@    yaxis  label char size 1.480000\n')
f.write('@    yaxis  bar linewidth 3.0\n')
f.write('@    yaxis  ticklabel char size 1.480000\n')
f.write('@    s0 line linewidth 2.5\n')
f.write('@    s0 symbol 1\n')
f.write('@    s0 symbol size 1.000000\n')
f.write('@    s0 symbol color 1\n')
f.write('@    s0 symbol pattern 1\n')
f.write('@    s0 symbol fill color 2\n')
f.write('@    s0 symbol fill pattern 1\n')
f.write('@    s0 symbol linewidth 1.0\n')
f.write('@TYPE xy \n')
f.write("".join(bfactors))
f.close()
!xmgrace bfactor.dat
#Cargando la imagen generada en xmgrace
Image(filename='bfactor.png')
Explanation: B-factors
Generando archivo para visualizarlo con XMGRACE
End of explanation
#Inicializando vector
bfactors=[]
try:
     file_bfactor = open( 'bfac.pdb' )
except IOError:
     print ('No se pudo abrir el archivo o no existe·..')
i=0
print ('Residuo' + '\t'+'bfactor')
for linea in file_bfactor.readlines():
    fila = linea.strip()
    sl = fila.split()
    if (sl[0]=='ATOM'):
        #print (sl[0])
        idresidue=fila[23:26]
        bfactor=fila[60:66]
        print (idresidue + '\t'+bfactor)
        bfactors.append(idresidue+'\t'+bfactor+'\n')
        #i=i+1
#Escribiendo el archivo BFACTOR.dat
f = open('bfactor.dat', 'w')
f.write("".join(bfactors))
f.close()
data_bfactor=np.loadtxt('bfactor.dat',comments=['#', '@'])
#Engrosar marco    
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
  ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
#ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
pl.plot(data_bfactor[:,0], data_bfactor[:,1],  '-o', color='black', markersize=25,
        markerfacecolor='red',markeredgecolor='black',markeredgewidth=3, linewidth = 4, )
pl.xlabel('Residue', fontsize = 40)
pl.ylabel('B-factors ('+ r'$\AA$'+')' , fontsize = 40)
#pl.title('B-Factors', fontsize=40)
pl.xticks(fontsize=30) 
pl.yticks(fontsize=30) 
pl.xlim(0, len(data_bfactor[:,1]))
Explanation: Generando archivo para visualizar con Matplotlib
End of explanation
### Creando el directorio para el análisis del RMSF
#Verificando que exista la nueva carpeta para la conversión de trayectorias
ruta_ss = nuevaruta+'/estructura'
print ( ruta_ss )
if not os.path.exists(ruta_ss): 
    os.makedirs(ruta_ss)
    print ('Se ha creado la ruta ===>',ruta_ss)
else:
        print ("La ruta "+ruta_ss+" existe..!!!")
        
print ( 'Nos vamos a ....', ruta_ss )
os.chdir( ruta_ss )
Explanation: Secondary Structure
Se crea la carpeta para cálculo de la estructura
End of explanation
print ('Ejecutando el análisis de esctructura secundaria...') 
!echo 5 | do_dssp -f ../output.xtc -s ../ionized.pdb -o sec_est.xpm -tu ns
print ('\n Convirtiendo el archivo a ps...')
!xpm2ps -f sec_est.xpm -by 6 -bx .1 -o est_sec.eps
print('\nConvirtiendo a png...')
!convert -density 600 est_sec.eps -resize 1024x1024 est_sec.png
print ('Cargando el archivo...')
Image(filename='est_sec.png', width=1024)
Explanation: Calculando la estructura secundaria
Se necesita contar con el programa dssp en la ruta /usr/local/bin, el cual se enlaza con Gromacs 5
End of explanation
### Creando el directorio para el análisis del r-gyro
#Verificando que exista la nueva carpeta para la conversión de trayectorias
ruta_rgyro = nuevaruta+'/rgyro'
print ( ruta_rgyro )
if not os.path.exists(ruta_rgyro): 
    os.makedirs(ruta_rgyro)
    print ('Se ha creado la ruta ===>',ruta_rgyro)
else:
        print ("La ruta "+ruta_rgyro+" existe..!!!")
print ( 'Nos vamos a ....', ruta_rgyro)
os.chdir( ruta_rgyro )
Explanation: R-GYRATE
Se crea una carpeta rgiro para guardar los archivos generados.
End of explanation
print ('Ejecutando el análisis de rgyro...')
!echo 3 | g_gyrate -f ../output.xtc -s ../ionized.pdb -o gyrate.xvg
Explanation: Calculando el r-gyro con la opción (3) - C-alpha
Se calcula para los carbonos alfa.
End of explanation
#Inicializando vector
rgyro=[]
try:
     file_rmsf = open( 'gyrate.xvg' )
except IOError:
     print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in file_rmsf.readlines():
    fila = linea.strip()
    sl = fila.split()
    cadena=sl[0]
    if (not '#' in cadena) and (not '@' in cadena):
        num=float(sl[0])
        #num2=float(sl[1])
        num=num/1000
        rgyro.append(repr(num)+'\t'+sl[1]+'\n')
        i=i+1
#Escribiendo el archivo RGYRO.DAT
f = open('rgyro.dat', 'w')
#f.write('@    title "Radius of gyration" \n')
f.write('@    xaxis  label " Time (ns)" \n')
f.write('@    xaxis  label char size 1.480000\n')
f.write('@    xaxis  bar linewidth 3.0\n')
f.write('@    xaxis  ticklabel char size 1.480000\n')
f.write('@    yaxis  label "Rg (nm)" \n')
f.write('@    yaxis  label char size 1.480000\n')
f.write('@    yaxis  bar linewidth 3.0\n')
f.write('@    yaxis  ticklabel char size 1.480000\n')
f.write('@    s0 line linewidth 2.5\n')
f.write('@TYPE xy \n')
f.write("".join(rgyro))
f.close()
!xmgrace rgyro.dat
#Cargando la imagen generada en xmgrace
Image(filename='rgyro.png')
Explanation: Generando el archivo rgyro.dat para su análisis con XMGRACE
End of explanation
data_rgyro=np.loadtxt('gyrate.xvg',comments=['#', '@'])
    
#Engrosar marco    
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
  ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
pl.plot(data_rgyro[:,0]/1000, data_rgyro[:,1], linewidth = 2, color='black')
pl.xlabel("Time (ns)", fontsize = 40)
pl.ylabel('Rg (nm)', fontsize = 40)
#pl.suptitle('Radius of gyration', fontsize=50)
pl.xticks(fontsize=30) 
pl.yticks(fontsize=30) 
Explanation: Ploteando el archivo gyrate.xvg con matplotlib
End of explanation
### Creando el directorio para el análisis del RMSF
#Verificando que exista la nueva carpeta para la conversión de trayectorias
ruta_helix = nuevaruta+'/rmsd_helix'
print ( ruta_helix )
if not os.path.exists(ruta_helix): 
    os.makedirs(ruta_helix)
    print ('Se ha creado la ruta ===>',ruta_helix)
else:
        print ("La ruta "+ruta_helix+" existe..!!!")
print ( 'Nos vamos a ....', ruta_helix)
os.chdir( ruta_helix )
Explanation: RMSD Helix Alfa
Para realizar este análisis se debe cargar el pdb original de la proteina que se encuentra en la carpeta 01_BUILD.
Cargarlo con VMD y dirigirse al Menú EXTENSIONS -> ANALYSIS -> SEQUENCE VIEWER, en la cual se tomará el rango de átomos del campo Struct (H), el cual se proporcionará de la forma "resid X1 to X2" donde X1 es primer átomo de la helix y X2 el último átomo de la helix.
End of explanation
num=input('Número de hélices con las que cuenta la proteína:')
print (num)
if (int(num)==1):
    indices_ha1=input('Proporciona el rango de índices de la Hélice 1:')
    print (indices_ha1)
    r_helix_1=1
    r_helix_2=0
    r_helix_3=0
    r_helix_4=0
if (int(num)==2):
    indices_ha1=input('Proporciona el rango de índices de la Hélice 1:')
    print (indices_ha1)
    indices_ha2=input('Proporciona el rango de índices de la Hélice 2:')
    print (indices_ha2)
    r_helix_1=1
    r_helix_2=1
    r_helix_3=0
    r_helix_4=0
if (int(num)==3):
    indices_ha1=input('Proporciona el rango de índices de la Hélice 1:')
    print (indices_ha1)
    indices_ha2=input('Proporciona el rango de índices de la Hélice 2:')
    print (indices_ha2)
    indices_ha3=input('Proporciona el rango de índices de la Hélice 3:')
    print (indices_ha3)
    r_helix_1=1
    r_helix_2=1
    r_helix_3=1
    r_helix_4=0
if (int(num)==4):
    indices_ha1=input('Proporciona el rango de índices de la Hélice 1:')
    print (indices_ha1)
    indices_ha2=input('Proporciona el rango de índices de la Hélice 2:')
    print (indices_ha2)
    indices_ha3=input('Proporciona el rango de índices de la Hélice 3:')
    print (indices_ha3)
    indices_ha4=input('Proporciona el rango de índices de la Hélice 4:')
    print (indices_ha4)
    r_helix_1=1
    r_helix_2=1
    r_helix_3=1
    r_helix_4=1
#Script para vmd de la Hélice Alfa 2
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
if (r_helix_1==1):
    f = open('ha1.tcl', 'w')
    print(f)
    f.write('set psfFile '+ psf+' \n')
    f.write('set dcdFile '+ dcd+' \n')
    f.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f.write('set outfile ' +'[open ' +'rmsd_ha1.dat'+' w]\n')
    f.write('set nf [molinfo top get numframes]\n')
    f.write('\n#RMSD calculation loop\n')
    f.write('set f1 [atomselect top "'+indices_ha1+' " frame 0]\n')
    f.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f.write('    set sel [atomselect top "'+indices_ha1+' " frame $i]\n')
    f.write('    $sel move [measure fit $sel $f1]\n')
    f.write('    set time [expr {$i +1}]\n')
    f.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f.write('    puts $outfile "$time $time"\n')
    f.write('}\n')
    f.write('close $outfile')
    f.close()
if (r_helix_2==1):
    f = open('ha2.tcl', 'w')
    print(f)
    f.write('set psfFile '+ psf+' \n')
    f.write('set dcdFile '+ dcd+' \n')
    f.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f.write('set outfile ' +'[open ' +'rmsd_ha2.dat'+' w]\n')
    f.write('set nf [molinfo top get numframes]\n')
    f.write('\n#RMSD calculation loop\n')
    f.write('set f1 [atomselect top "'+indices_ha2+' " frame 0]\n')
    f.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f.write('    set sel [atomselect top "'+indices_ha2+' " frame $i]\n')
    f.write('    $sel move [measure fit $sel $f1]\n')
    f.write('    set time [expr {$i +1}]\n')
    f.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f.write('    puts $outfile "$time $time"\n')
    f.write('}\n')
    f.write('close $outfile')
    f.close()
if (r_helix_3==1):
    f = open('ha3.tcl', 'w')
    print(f)
    f.write('set psfFile '+ psf+' \n')
    f.write('set dcdFile '+ dcd+' \n')
    f.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f.write('set outfile ' +'[open ' +'rmsd_ha3.dat'+' w]\n')
    f.write('set nf [molinfo top get numframes]\n')
    f.write('\n#RMSD calculation loop\n')
    f.write('set f1 [atomselect top "'+indices_ha3+' " frame 0]\n')
    f.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f.write('    set sel [atomselect top "'+indices_ha3+' " frame $i]\n')
    f.write('    $sel move [measure fit $sel $f1]\n')
    f.write('    set time [expr {$i +1}]\n')
    f.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f.write('    puts $outfile "$time $time"\n')
    f.write('}\n')
    f.write('close $outfile')
    f.close()
if (r_helix_4==1):
    f = open('ha4.tcl', 'w')
    print(f)
    f.write('set psfFile '+ psf+' \n')
    f.write('set dcdFile '+ dcd+' \n')
    f.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f.write('set outfile ' +'[open ' +'rmsd_ha4.dat'+' w]\n')
    f.write('set nf [molinfo top get numframes]\n')
    f.write('\n#RMSD calculation loop\n')
    f.write('set f1 [atomselect top "'+indices_ha4+' " frame 0]\n')
    f.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f.write('    set sel [atomselect top "'+indices_ha4+' " frame $i]\n')
    f.write('    $sel move [measure fit $sel $f1]\n')
    f.write('    set time [expr {$i +1}]\n')
    f.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f.write('    puts $outfile "$time $time"\n')
    f.write('}\n')
    f.write('close $outfile')
    f.close()
if (r_helix_1==1):
    #Calculando con VMD hélice 1
    !vmd -dispdev text < ha1.tcl
if (r_helix_2==1):
    #Calculando con VMD hélice 2
    !vmd -dispdev text < ha2.tcl
if (r_helix_3==1):
    #Calculando con VMD hélice 3
    !vmd -dispdev text < ha3.tcl
if (r_helix_4==1):
    #Calculando con VMD hélice 4
    !vmd -dispdev text < ha4.tcl
if (int(num)==1):
    #Graficando
    data_ha1=np.loadtxt('rmsd_ha1.dat',comments=['#', '@'])
    #Engrosar marco    
    fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
    ax = fig.add_subplot(111)
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(4)
    #Formateando los valores de los ejes
    ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    #pl.plot(data_ha1[:,0], data_ha1[:,1], linewidth = 3)
    pl.plot(data_ha1[:,1]*0.02, data_ha1[:,0]/10, linewidth = 3, color='black')
    pl.xlabel("Time (ns)", fontsize = 40)
    pl.ylabel('RMSD (nm)', fontsize = 40)
    #pl.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
    #pl.title('RMSD Helix Alfa', fontsize=50)
    pl.xticks(fontsize=30) 
    pl.yticks(fontsize=30) 
    
if (int(num)==2):
    #Graficando
    data_ha1=np.loadtxt('rmsd_ha1.dat',comments=['#', '@'])
    data_ha2=np.loadtxt('rmsd_ha2.dat',comments=['#', '@'])
    #Engrosar marco    
    fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
    ax = fig.add_subplot(111)
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(4)
    #Formateando los valores de los ejes
    ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    #pl.plot(data_ha1[:,0], data_ha1[:,1], linewidth = 3)
    pl.plot(data_ha1[:,1]*0.02, data_ha1[:,0]/10, linewidth = 3, color='black')
    pl.plot(data_ha2[:,1]*0.02, data_ha2[:,0]/10, linewidth = 3, color='red')
    pl.xlabel("Time (ns)", fontsize = 40)
    pl.ylabel('RMSD (nm)', fontsize = 40)
    #pl.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
    #pl.title('RMSD Helix Alfa', fontsize=50)
    pl.xticks(fontsize=30) 
    pl.yticks(fontsize=30) 
if (int(num)==3):
    #Graficando
    data_ha1=np.loadtxt('rmsd_ha1.dat',comments=['#', '@'])
    data_ha2=np.loadtxt('rmsd_ha2.dat',comments=['#', '@'])
    data_ha3=np.loadtxt('rmsd_ha3.dat',comments=['#', '@'])
    #Engrosar marco    
    fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
    ax = fig.add_subplot(111)
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(4)
    #Formateando los valores de los ejes
    ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    #pl.plot(data_ha1[:,0], data_ha1[:,1], linewidth = 3)
    pl.plot(data_ha1[:,1]*0.02, data_ha1[:,0]/10, linewidth = 3, color='black')
    pl.plot(data_ha2[:,1]*0.02, data_ha2[:,0]/10, linewidth = 3, color='red')
    pl.plot(data_ha3[:,1]*0.02, data_ha3[:,0]/10, linewidth = 3, color='green')
    pl.xlabel("Time (ns)", fontsize = 40)
    pl.ylabel('RMSD (nm)', fontsize = 40)
    #pl.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
    #pl.title('RMSD Helix Alfa', fontsize=50)
    pl.xticks(fontsize=30) 
    pl.yticks(fontsize=30)
if (int(num)==4):
    #Graficando
    data_ha1=np.loadtxt('rmsd_ha1.dat',comments=['#', '@'])
    data_ha2=np.loadtxt('rmsd_ha2.dat',comments=['#', '@'])
    data_ha3=np.loadtxt('rmsd_ha3.dat',comments=['#', '@'])
    data_ha4=np.loadtxt('rmsd_ha4.dat',comments=['#', '@'])
    #Engrosar marco    
    fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
    ax = fig.add_subplot(111)
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(4)
    #Formateando los valores de los ejes
    ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
    #pl.plot(data_ha1[:,0], data_ha1[:,1], linewidth = 3)
    pl.plot(data_ha1[:,1]*0.02, data_ha1[:,0]/10, linewidth = 3, color='black')
    pl.plot(data_ha2[:,1]*0.02, data_ha2[:,0]/10, linewidth = 3, color='red')
    pl.plot(data_ha3[:,1]*0.02, data_ha3[:,0]/10, linewidth = 3, color='green')
    pl.plot(data_ha4[:,1]*0.02, data_ha4[:,0]/10, linewidth = 3, color='blue')
    pl.xlabel("Time (ns)", fontsize = 40)
    pl.ylabel('RMSD (A)', fontsize = 40)
    #pl.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
    #pl.title('RMSD Helix Alfa', fontsize=50)
    pl.xticks(fontsize=30) 
    pl.yticks(fontsize=30)
Explanation: Entrada de datos
Para la entrada se deberá dar con la opción "resid X to X".
End of explanation
### Creando el directorio para el análisis del SASA
### NOTA: se calcula con gromacs4  ya que arroja bien los resultados comparado con gromacs5
ruta_sasa = nuevaruta+'/sasa'
print ( ruta_sasa )
if not os.path.exists(ruta_sasa): 
    os.makedirs(ruta_sasa)
    print ('Se ha creado la ruta ===>',ruta_sasa)
else:
        print ("La ruta "+ruta_sasa+" existe..!!!")
print ( 'Nos vamos a ....', ruta_sasa )
os.chdir( ruta_sasa )
Explanation: SASA
Creando la estructura de carpeta para el cálculo
End of explanation
print ('Ejecutando el análisis de sasa con Gromacs 4 utilizando la opción 1 (protein)...')
!echo 1 1 | /opt/gromacs4/bin/g_sas -f ../output.xtc -s ../ionized.pdb -o solven-accessible-surface.xvg -oa atomic-sas.xvg -or residue-sas.xvg
Explanation: Ejecutando el análisis de SASA con Gromacs4
End of explanation
#Inicializando vector
sasa_residuo=[]
try:
     residue_sas = open( 'residue-sas.xvg' )
except IOError:
     print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in residue_sas.readlines():
    fila = linea.strip()
    sl = fila.split()
    cadena=sl[0]
    if (not '#' in cadena) and (not '@' in cadena):
        print ('Residue =>',cadena)
        sasa_residuo.append(sl[0]+'\t'+sl[1]+'\n')
        i=i+1
#Escribiendo el archivo RMSD
f = open('sasa-residuo.dat', 'w')
#f.write('@    title "Area per residue over the trajectory" \n')
f.write('@    xaxis  label " Residue " \n')
f.write('@    xaxis  label char size 1.480000\n')
f.write('@    xaxis  bar linewidth 3.0\n')
f.write('@    xaxis  ticklabel char size 1.480000\n')
f.write('@    yaxis  label "Area (nm' +"\\"+'S2'+"\\N"+')"\n')
f.write('@    yaxis  label char size 1.480000\n')
f.write('@    yaxis  bar linewidth 3.0\n')
f.write('@    yaxis  ticklabel char size 1.480000\n')
f.write('@    s0 line linewidth 2.5\n')
f.write('@    s0 symbol 1\n')
f.write('@    s0 symbol size 1.000000\n')
f.write('@    s0 symbol color 1\n')
f.write('@    s0 symbol pattern 1\n')
f.write('@    s0 symbol fill color 2\n')
f.write('@    s0 symbol fill pattern 1\n')
f.write('@    s0 symbol linewidth 1.0\n')
f.write('@TYPE xy \n')
f.write("".join(sasa_residuo))
f.close()
!xmgrace sasa-residuo.dat
#Cargando la imagen generada en xmgrace
Image(filename='sasa-residuo.png')
Explanation: Creando el archivo sasa_residuo.dat para salida con XMGRACE
End of explanation
data_sasa_residue=np.loadtxt('residue-sas.xvg',comments=['#', '@'])
    
#Engrosar marco    
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f'))
pl.plot(data_sasa_residue[:,0], data_sasa_residue[:,1],  '-o', color='black', markersize=25,
        markerfacecolor='red',markeredgecolor='black',markeredgewidth=3, linewidth = 4, )
pl.xlabel("Residue", fontsize = 30)
#pl.ylabel('Area (nm2)', fontsize = 30)
pl.ylabel('Area ( nm'+ r'$\ ^2$'+')' , fontsize = 40)
#pl.title('Area per residue over the trajectory', fontsize=40)
pl.xticks(fontsize=30) 
pl.yticks(fontsize=30) 
pl.xlim(0, len(data_sasa_residue[:,1]))
Explanation: Cargando archivo residue-sas.xvg para su visualización en Matplotlib
Se genera el gráfico de salida para matplotlib
End of explanation
#Inicializando vector
sasa=[]
try:
     sasafile = open( 'solven-accessible-surface.xvg' )
except IOError:
     print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in sasafile.readlines():
    fila = linea.strip()
    sl = fila.split()
    cadena=sl[0]
    if (not '#' in cadena) and (not '@' in cadena):
        #print (cadena)
        num=float(sl[0])
        num=num/1000
        sasa.append(repr(num)+'\t'+sl[1]+'\t'+sl[2]+'\t'+sl[3]+'\n')
        i=i+1
cel2=float(sl[2])
print(cel2)
#Escribiendo el archivo RMSD
f = open('sasa.dat', 'w')
#f.write('@    title "Solven Accessible Surface" \n')
f.write('@    xaxis  label " Time (ns) " \n')
f.write('@    xaxis  label char size 1.480000\n')
f.write('@    xaxis  bar linewidth 3.0\n')
f.write('@    xaxis  ticklabel char size 1.480000\n')
f.write('@    yaxis  label "Area (nm' +"\\"+'S2'+"\\N"+')"\n')
f.write('@    yaxis  label char size 1.480000\n')
f.write('@    yaxis  bar linewidth 3.0\n')
f.write('@    yaxis  ticklabel char size 1.480000\n')
#f.write('@    s0 legend "Hydrophobic"\n')
#if (cel2>0):
    #f.write('@    s1 legend "Hydrophilic"\n')
f.write('@TYPE xy \n')
f.write("".join(sasa))
f.close()
!xmgrace sasa.dat
#Cargando la imagen generada en xmgrace
Image(filename='sasa.png')
Explanation: Creando el archivo sasa.dat para salida con XMGRACE
End of explanation
data_sasa=np.loadtxt('solven-accessible-surface.xvg',comments=['#', '@'])
    
#Engrosar marco    
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
#ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f'))
pl.xlabel("Time (ns)", fontsize = 40)
pl.ylabel('Area ( nm'+ r'$\ ^2$'+')' , fontsize = 40)
#pl.title('Solvent Accessible Surface', fontsize=50)
pl.xticks(fontsize=30) 
pl.yticks(fontsize=30) 
dato=data_sasa[:,2]
dato2=dato[0]
if (dato2>0):
    pl.plot(data_sasa[:,0]/1000, data_sasa[:,1], linewidth = 2, color='black')
    pl.plot(data_sasa[:,0]/1000, data_sasa[:,2], linewidth = 2, color='red')
else:
    pl.plot(data_sasa[:,0]/1000, data_sasa[:,1], linewidth = 2, color='black')
Explanation: Cargando archivo solven-accessible-surface.xvg para graficar con Matplotlib
End of explanation
### Creando el directorio para el análisis del SASA
### NOTA: se calcula con gromacs4  ya que arroja bien los resultados comparado con gromacs5
ruta_m_rmsd = nuevaruta+'/matriz'
print ( ruta_m_rmsd )
if not os.path.exists(ruta_m_rmsd): 
    os.makedirs(ruta_m_rmsd)
    print ('Se ha creado la ruta ===>',ruta_m_rmsd)
else:
        print ("La ruta "+ruta_m_rmsd+" existe..!!!")
print ( 'Nos vamos a ....', ruta_m_rmsd )
os.chdir( ruta_m_rmsd )
print ('\nCopiando el archivo rmsd_matrix.tcl a '+ruta_m_rmsd)
source_file=ruta_scripts+'/rmsd_matriz/rmsd_matrix.tcl'
dest_file=ruta_m_rmsd+'/rmsd_matrix.tcl'
shutil.copy(source_file,dest_file)
#print ( 'Nos vemos a ....', ruta_old_traj )
#os.chdir( ruta_old_traj )
file_dcd=ruta_old_traj+'/'+dcd_file
file_psf=ruta_old_traj+'/'+psf_file
print (file_dcd)
print ('\nEjecutando CATDCD para obtener 100 frames de la trayectoria original....')
output_catdcd=!catdcd -o 100.dcd -stride 50 $file_dcd
print (output_catdcd.n)
Explanation: MATRIZ DE RMSD
End of explanation
#Arrancando VMD para cargar el script rmsd_matrix.tcl
!vmd 100.dcd $file_psf
ruta_matriz=os.getcwd()
if os.path.isfile('salida.dat'):
    print ('El archivo salida.dat existe')
else:
    print ('El archivo salida.dat no existe.. ejecutar desde MATRIZ DE RMSD...')
Explanation: Cargar el scrit rmsd_matrix con vmd en la nueva trayectoria
Arrancar VMD, dirigirse al manú Extensions -> Tk Console, copiar y ejecutar la siguiente secuencia de comandos:
tcl
source rmsd_matrix.tcl
rmsd_matrix -mol top -seltext "name CA" -frames all -o salida.dat
exit
End of explanation
#Creando el gráfico
data_matriz=np.loadtxt('salida.dat',comments=['#', '@'])
print(data_matriz.shape)
pl.figure(figsize=(20, 12), dpi=100)
imgplot = pl.imshow(data_matriz, origin='lower', cmap=pl.cm.Greens, interpolation='nearest')
#imgplot = pl.imshow(data_matriz, origin='lower', cmap=pl.cm.coolwarm, interpolation='nearest')
pl.xlabel("Time (ns)", fontsize = 30)
pl.ylabel('Time (ns)', fontsize = 30)
#pl.suptitle('RMSD', fontsize=50)
#pl.title('C-Alpha RMSD matrix', fontsize=40)
pl.xticks(fontsize=20) 
pl.yticks(fontsize=20) 
pl.xlim(0, 100)
pl.ylim(0, 100)
pl.colorbar()
Explanation: Graficando el archivo de salida
End of explanation
### Creando el directorio para el análisis del RMSF
#Verificando que exista la nueva carpeta para la conversión de trayectorias
ruta_matriz_dm = nuevaruta+'/matriz_dm'
print ( ruta_matriz_dm )
if not os.path.exists(ruta_matriz_dm): 
    os.makedirs(ruta_matriz_dm)
    print ('Se ha creado la ruta ===>',ruta_matriz_dm)
else:
        print ("La ruta "+ruta_matriz_dm+" existe..!!!")
        
print ( 'Nos vamos a ....', ruta_matriz_dm )
os.chdir( ruta_matriz_dm )
Explanation: Matriz de distancia mínima
End of explanation
!echo 4 | g_mdmat -f ../output.xtc -s ../ionized.pdb -mean average -frames frames -dt 10000
Explanation: Calculando la matriz de distancia mínima
Seleccionar el backbone (opción 4)
End of explanation
!xpm2ps -f frames.xpm -o frames.eps
!xpm2ps -f average.xpm -o average.eps
print('\nConvirtiendo a png...')
!convert -density 600 frames.eps -resize 1024x1024 frames.png
!convert -density 600 average.eps -resize 1024x1024 average.png
print ('Cargando el archivo average...')
Image(filename='average.png', width=800)
Explanation: Generando los archivos  para visualizarlos
End of explanation
### Creando el directorio para el análisis de la libre energía
ruta_f_energy = nuevaruta+'/free_energy'
print ( ruta_f_energy )
if not os.path.exists(ruta_f_energy): 
    os.makedirs(ruta_f_energy)
    print ('Se ha creado la ruta ===>',ruta_f_energy)
else:
        print ("La ruta "+ruta_f_energy+" existe..!!!")
print ( 'Nos vamos a ....', ruta_f_energy )
os.chdir( ruta_f_energy )
#Solicita la temperatura
t=input('Temperatura a la cual se realizó la simulación:')
temperatura=int(t)
print ('Temperatura=>',temperatura)
Explanation: Free Energy
Para el cálculo de la energía libre se requiere el valor mínimo y máximo del RMSD y del radio de gyro, así como el valor de la temperatura a la cual se realizó la simulación. Estos datos son de entrada para el script del cálculo del mismo.
End of explanation
print ('Ejecutando el análisis de rmsd...')
!echo 3 3 | g_rms -f ../output.xtc -s ../ionized.pdb -a avgrp.xvg
print ('Ejecutando el análisis de rgyro...')
!echo 3 | g_gyrate -f ../output.xtc -s ../ionized.pdb -o gyrate.xvg
Explanation: Calculando el rmsd y el r-gyro para obtener el mínimo y máximo de cada uno de ellos.
End of explanation
print ('\nCopiando el archivo generateFES.py a '+ruta_f_energy)
source_file=ruta_scripts+'/free_energy/generateFES.py'
dest_file=ruta_f_energy+'/generateFES.py'
shutil.copy(source_file,dest_file)
#Cambiando permisos de ejecución
!chmod +x generateFES.py
Explanation: Escribiendo script a /tmp para utilizar en el cálculo
End of explanation
#Cargando valores del RMSD
data_rmsd=np.loadtxt('rmsd.xvg',comments=['#', '@'])
#Cargnaod valores del R-GYRO
data_rgyro=np.loadtxt('gyrate.xvg',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del rmsd
min_rmsd=np.amin(data_rmsd[:,1])
max_rmsd=np.amax(data_rmsd[:,1])
print ('Minimo RMSD=>',min_rmsd)
print ('Máximo RMSD=>',max_rmsd)
#Obteniendo los valores máximo y mínimo del r-gyro
min_rgyro=np.amin(data_rgyro[:,1])
max_rgyro=np.amax(data_rgyro[:,1])
print ('Minimo RGYRO=>',min_rgyro)
print ('Máximo RGYRO=>',max_rgyro)
#Creando los archivos de entrada para el script
np.savetxt('rmsd.dat',data_rmsd[:,1], fmt='%1.7f')
np.savetxt('rgyro.dat',data_rgyro[:,1], fmt='%1.7f')
!paste rgyro.dat rmsd.dat  > fes.dat
#Ejecutando el script de FES
!python generateFES.py fes.dat $min_rgyro $max_rgyro $min_rmsd $max_rmsd  200 200 $temperatura FEES.dat
#Cargando el archivo generado para plotear con matplotlib
data_fes=np.loadtxt('FEES.dat',comments=['#', '@'])
Explanation: Realizando los cálculos para la Free Energy
End of explanation
# This loads the magics for gnuplot
%load_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "free_energy.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "Rg (nm)
set ylabel "RMSD (nm)"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "FEES.dat" with pm3d
Explanation: Ploteando con GNUplot
End of explanation
### Creando el directorio para el análisis del PCA
ruta_pca = nuevaruta+'/pca'
print ( ruta_pca )
if not os.path.exists(ruta_pca): 
    os.makedirs(ruta_pca)
    print ('Se ha creado la ruta ===>',ruta_pca)
else:
        print ("La ruta "+ruta_pca+" existe..!!!")
print ( 'Nos vamos a ....', ruta_pca )
os.chdir( ruta_pca )
#Calculando matriz de covarianza
!echo 1 1 | g_covar -s ../ionized.pdb -f ../output.xtc -o eigenvalues.xvg -v eigenvectors.trr -xpma covar.xpm
Explanation: PCA
End of explanation
!echo 1 1 | g_anaeig -s ../ionized.pdb -f ../output.xtc -v eigenvectors.trr -eig eigenvalues.xvg -first 1 -last 2 -2d 2dproj_1_2.xvg
#pcaX, pcaY=np.loadtxt('2dproj_1_2.xvg',comments=['#', '@'], unpack=True)
data_pca=np.loadtxt('2dproj_1_2.xvg',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del pca
min_pcaX=np.amin(data_pca[:,0])
max_pcaX=np.amax(data_pca[:,0])
print ('Minimo PCA_X=>',min_pcaX)
print ('Máximo PCA_X=>',max_pcaX)
min_pcaY=np.amin(data_pca[:,1])
max_pcaY=np.amax(data_pca[:,1])
print ('Minimo PCA_Y=>',min_pcaY)
print ('Máximo PCA_Y=>',max_pcaY)
#Creando los archivos de entrada para el script
np.savetxt('PCA.dat',data_pca, fmt='%1.5f')
#Copiando el script generateFES de la carpeta Free_energy
print ('\nCopiando el archivo generateFES.py a '+ruta_pca+ ' desde '+ ruta_f_energy)
source_file=ruta_f_energy+'/generateFES.py'
dest_file=ruta_pca+'/generateFES.py'
shutil.copy(source_file,dest_file)
#Ejecutando el script de FES
!python generateFES.py PCA.dat $min_pcaX $max_pcaX $min_pcaY $max_pcaY  200 200 $temperatura FEES_PCA.dat
Explanation: Una vez calculada la matriz el eigenvalues y eigenvectors sirven de entrada para generar el pca.
El siguiente comando representa el movimiento del primer y segundo eigenvector.
End of explanation
#Volver a cargar el kernel de gnuplot para limpiar su buffer
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "pca.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "projection on eigenvector 1 (nm)"
set ylabel "projection on eigenvector 2 (nm)"
set title " "
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "FEES_PCA.dat" with pm3d
Explanation: Ploteando el archivo con gnuplot
End of explanation
from htmd import *
Explanation: Análisis de puentes di sulfuro
Este aplica para 2 puente, para lo cual se utiliza el software HTMD.
End of explanation
### Creando el directorio para el análisis de los RMSD de los puentes
ruta_rmsd_diedros = nuevaruta+'/rmsd_diedros'
print ( ruta_rmsd_diedros )
if not os.path.exists(ruta_rmsd_diedros): 
    os.makedirs(ruta_rmsd_diedros)
    print ('Se ha creado la ruta ===>',ruta_rmsd_diedros)
else:
        print ("La ruta "+ruta_rmsd_diedros+" existe..!!!")
print ( 'Nos vamos a ....', ruta_rmsd_diedros)
os.chdir( ruta_rmsd_diedros )
Explanation: Creando la ruta
Ruta para el análisis de los datos.
End of explanation
# Cargando la molécula
mol = Molecule('../ionized.pdb')
# Solicitando los datos de entrada
px1l=input('Índice del DB1 izquierdo:')
px1r=input('Índice del DB1 derecho:')
px2l=input('Índice del DB2 izquierdo:')
px2r=input('Índice del DB2 derecho:')
revisa1=1
revisa2=1
Explanation: Cargando de los puentes di sulfuro
Para este análisis se deberá revisar el archivo psf_charmm.tcl de la carpeta 01_BUILD, en el cual se tiene la definición de los puentes como la siguiente:
    patch DISU A:4   A:22
    patch DISU A:8   A:18
El número del puente se determinará de acuerdo al orden en que se encuentran definidos en este archivo, por ejemplo, la nota anterior:
    DB1 4-22
    DB2 8-18
La entrada de datos será por los índices del lado izquierdo y derecho respectivamente, con los cuales se creará la estructura completa de cada uno de ellos tomando los valores de los indices para su respectivo análisis.
End of explanation
if (revisa1>0):
    #Obteniendo lado izquierdo del DB1
    x1l_name=mol.get('name','resname CYS and noh and resid '+px1l)
    x1l_index=mol.get('index','resname CYS and noh and resid '+px1l)
    x1l_resid=mol.get('resid','resname CYS and noh and resid '+px1l)
    #Obteniendo lado derecho del DB1
    x1r_name=mol.get('name','resname CYS and noh and resid '+px1r)
    x1r_index=mol.get('index','resname CYS and noh and resid '+px1r)
    x1r_resid=mol.get('resid','resname CYS and noh and resid '+px1r)
if (revisa2>0):
    #Obteniendo el lado izquierdo del DB2
    x2l_name=mol.get('name','resname CYS and noh and resid '+px2l)
    x2l_index=mol.get('index','resname CYS and noh and resid '+px2l)
    x2l_resid=mol.get('resid','resname CYS and noh and resid '+px2l)
    #Obteniendo el lado derecho del DB2
    x2r_name=mol.get('name','resname CYS and noh and resid '+px2r)
    x2r_index=mol.get('index','resname CYS and noh and resid '+px2r)
    x2r_resid=mol.get('resid','resname CYS and noh and resid '+px2r)
    
#Obteniendo la lista de índices de los puentes
print ('Generando la lista de los índices para enviarlos')
db1x1l=[]
db1x2l=[]
db1x3m=[]
db1x2r=[]
db1x1r=[]
db1l_name_l=[]
db1l_index_l=[]
db1r_name_l=[]
db1r_index_l=[]
db2l_name_l=[]
db2l_index_l=[]
db2r_name_l=[]
db2r_index_l=[]
db3l_name_l=[]
db3l_index_l=[]
db3r_name_l=[]
db3r_index_l=[]
if (revisa1>0):
    #Obteniendo los índices del DB1
    for i in range(len(x1l_name)):
        if (x1l_name[i]=='N' or x1l_name[i]=='CA' or x1l_name[i]=='CB' or x1l_name[i]=='SG'):
            db1l_name_l.append(str(x1l_name[i]))
            db1l_index_l.append(str(x1l_index[i]))
    for i in range(len(x1r_name)):
        if (x1r_name[i]=='N' or x1r_name[i]=='CA' or x1r_name[i]=='CB' or x1r_name[i]=='SG'):
            db1r_name_l.append(str(x1r_name[i]))
            db1r_index_l.append(str(x1r_index[i])) 
    print ('DB1 X1L =>',db1l_name_l)
    print (db1l_index_l)
    print ('DB1 X1R =>',db1r_name_l)
    print (db1r_index_l)
if (revisa2>0):
    #Obteniendo los índices del DB2
    for i in range(len(x2l_name)):
        if (x2l_name[i]=='N' or x2l_name[i]=='CA' or x2l_name[i]=='CB' or x2l_name[i]=='SG'):
            db2l_name_l.append(str(x2l_name[i]))
            db2l_index_l.append(str(x2l_index[i]))
    for i in range(len(x2r_name)):
        if (x2r_name[i]=='N' or x2r_name[i]=='CA' or x2r_name[i]=='CB' or x2r_name[i]=='SG'):
            db2r_name_l.append(str(x2r_name[i]))
            db2r_index_l.append(str(x2r_index[i])) 
    print ('DB2 X1L =>',db2l_name_l)
    print (db2l_index_l)
    print ('DB2 X1R =>',db2r_name_l)
    print (db2r_index_l)
Explanation: Obteniendo los índices de los puentes
End of explanation
#Generando el DB1 completo ordenado
filas=8
col=2
DB1_i=[]
DB1_N=[]
DB2_i=[]
DB2_N=[]
DB3_i=[]
DB3_N=[]
for i in range(0,filas):
    DB1_N.append([' '])
    DB1_i.append(['0'])
    DB2_N.append([' '])
    DB2_i.append(['0'])
    DB3_N.append([' '])
    DB3_i.append(['0'])
if (revisa1>0):
    #Cargando índices para el puente 1
    for i in range(len(db1l_name_l)):
        if db1l_name_l[i]=='N':
            DB1_N[0] = db1l_name_l[i]
            DB1_i[0]='index '+db1l_index_l[i]
        if db1l_name_l[i]=='CA':
            DB1_N[1] = db1l_name_l[i]
            DB1_i[1]='index '+db1l_index_l[i]
        if db1l_name_l[i]=='CB':
            DB1_N[2] = db1l_name_l[i]
            DB1_i[2]='index '+db1l_index_l[i]
        if db1l_name_l[i]=='SG':
            DB1_N[3] = db1l_name_l[i]
            DB1_i[3]='index '+db1l_index_l[i]
        
    for i in range(len(db1r_name_l)):
        if db1r_name_l[i]=='SG':
            DB1_N[4] = db1r_name_l[i]
            DB1_i[4]='index '+db1r_index_l[i]
        if db1r_name_l[i]=='CB':
            DB1_N[5] = db1r_name_l[i]
            DB1_i[5]='index '+db1r_index_l[i]
        if db1r_name_l[i]=='CA':
            DB1_N[6] = db1r_name_l[i]
            DB1_i[6]='index '+db1r_index_l[i]
        if db1r_name_l[i]=='N':
            DB1_N[7] = db1r_name_l[i]
            DB1_i[7]='index '+db1r_index_l[i]
    
    print ('Puente DB1 = resid '+px1l+':'+px1r)
    print ('Names DB1=>',DB1_i)
    print ('Index DB1=>',DB1_N)
    print ('\n')
if (revisa2>0):
    #Cargando índices para el puente 2
    for i in range(len(db2l_name_l)):
        if db2l_name_l[i]=='N':
            DB2_N[0] = db2l_name_l[i]
            DB2_i[0]='index '+db2l_index_l[i]
        if db2l_name_l[i]=='CA':
            DB2_N[1] = db2l_name_l[i]
            DB2_i[1]='index '+db2l_index_l[i]
        if db2l_name_l[i]=='CB':
            DB2_N[2] = db2l_name_l[i]
            DB2_i[2]='index '+db2l_index_l[i]
        if db2l_name_l[i]=='SG':
            DB2_N[3] = db2l_name_l[i]
            DB2_i[3]='index '+db2l_index_l[i]
        
    for i in range(len(db2r_name_l)):
        if db2r_name_l[i]=='SG':
            DB2_N[4] = db2r_name_l[i]
            DB2_i[4]='index '+db2r_index_l[i]
        if db2r_name_l[i]=='CB':
            DB2_N[5] = db2r_name_l[i]
            DB2_i[5]='index '+db2r_index_l[i]
        if db2r_name_l[i]=='CA':
            DB2_N[6] = db2r_name_l[i]
            DB2_i[6]='index '+db2r_index_l[i]
        if db2r_name_l[i]=='N':
            DB2_N[7] = db2r_name_l[i]
            DB2_i[7]='index '+db2r_index_l[i]
            
    print ('Puente DB2 = resid '+px2l+':'+px2r)
    print ('Names DB2=>',DB2_i)
    print ('Index DB2=>',DB2_N)
    print ('\n')
Explanation: Ordenando los puentes  de la forma ['N', 'CA', 'CB', 'SG', 'SG', 'CB', 'CA', 'N']
End of explanation
if (revisa1>0):
  
    #Creando script para DB1_x1l
    psf=ruta_old_traj+'/'+psf_file
    dcd=ruta_old_traj+'/'+dcd_file
    print(psf)
    f1 = open('DB1_x1l.tcl', 'w')
    print(f1)
    f1.write('set psfFile '+ psf+' \n')
    f1.write('set dcdFile '+ dcd+' \n')
    f1.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f1.write('set outfile ' +'[open ' +'db1_x1l.dat'+' w]\n')
    f1.write('set nf [molinfo top get numframes]\n')
    f1.write('\n#RMSD calculation loop\n')
    f1.write('set f1 [atomselect top "'+DB1_i[0]+' or '+DB1_i[1]+' or '+DB1_i[2]+' or '+DB1_i[3]+' " frame 0]\n')
    f1.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f1.write('    set sel [atomselect top "'+DB1_i[0]+' or '+DB1_i[1]+' or '+DB1_i[2]+' or '+DB1_i[3]+' " frame $i]\n')
    f1.write('    $sel move [measure fit $sel $f1]\n')
    f1.write('    set time [expr {$i +1}]\n')
    f1.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f1.write('    puts $outfile "     $time"\n')
    f1.write('}\n')
    f1.write('close $outfile')
    f1.close()
    
    #Creando script para DB1_x2l
    psf=ruta_old_traj+'/'+psf_file
    dcd=ruta_old_traj+'/'+dcd_file
    print(psf)
    f2 = open('DB1_x2l.tcl', 'w')
    print(f2)
    f2.write('set psfFile '+ psf+' \n')
    f2.write('set dcdFile '+ dcd+' \n')
    f2.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f2.write('set outfile ' +'[open ' +'db1_x2l.dat'+' w]\n')
    f2.write('set nf [molinfo top get numframes]\n')
    f2.write('\n#RMSD calculation loop\n')
    f2.write('set f1 [atomselect top "'+DB1_i[1]+' or '+DB1_i[2]+' or '+DB1_i[3]+' or '+DB1_i[4]+' " frame 0]\n')
    f2.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f2.write('    set sel [atomselect top "'+DB1_i[1]+' or '+DB1_i[2]+' or '+DB1_i[3]+' or '+DB1_i[4]+' " frame $i]\n')
    f2.write('    $sel move [measure fit $sel $f1]\n')
    f2.write('    set time [expr {$i +1}]\n')
    f2.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f2.write('    puts $outfile "     $time"\n')
    f2.write('}\n')
    f2.write('close $outfile')
    f2.close()
    
    #Creando script para DB1_x3m
    psf=ruta_old_traj+'/'+psf_file
    dcd=ruta_old_traj+'/'+dcd_file
    print(psf)
    f3 = open('DB1_x3m.tcl', 'w')
    print(f3)
    f3.write('set psfFile '+ psf+' \n')
    f3.write('set dcdFile '+ dcd+' \n')
    f3.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f3.write('set outfile ' +'[open ' +'db1_x3m.dat'+' w]\n')
    f3.write('set nf [molinfo top get numframes]\n')
    f3.write('\n#RMSD calculation loop\n')
    f3.write('set f1 [atomselect top "'+DB1_i[2]+' or '+DB1_i[3]+' or '+DB1_i[4]+' or '+DB1_i[5]+' " frame 0]\n')
    f3.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f3.write('    set sel [atomselect top "'+DB1_i[2]+' or '+DB1_i[3]+' or '+DB1_i[4]+' or '+DB1_i[5]+' " frame $i]\n')
    f3.write('    $sel move [measure fit $sel $f1]\n')
    f3.write('    set time [expr {$i +1}]\n')
    f3.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f3.write('    puts $outfile "     $time"\n')
    f3.write('}\n')
    f3.write('close $outfile')
    f3.close()
    
    #Creando script para DB1_x2r
    psf=ruta_old_traj+'/'+psf_file
    dcd=ruta_old_traj+'/'+dcd_file
    print(psf)
    f4 = open('DB1_x2r.tcl', 'w')
    print(f4)
    f4.write('set psfFile '+ psf+' \n')
    f4.write('set dcdFile '+ dcd+' \n')
    f4.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f4.write('set outfile ' +'[open ' +'db1_x2r.dat'+' w]\n')
    f4.write('set nf [molinfo top get numframes]\n')
    f4.write('\n#RMSD calculation loop\n')
    f4.write('set f1 [atomselect top "'+DB1_i[3]+' or '+DB1_i[4]+' or '+DB1_i[5]+' or '+DB1_i[6]+' " frame 0]\n')
    f4.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f4.write('    set sel [atomselect top "'+DB1_i[3]+' or '+DB1_i[4]+' or '+DB1_i[5]+' or '+DB1_i[6]+' " frame $i]\n')
    f4.write('    $sel move [measure fit $sel $f1]\n')
    f4.write('    set time [expr {$i +1}]\n')
    f4.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f4.write('    puts $outfile "     $time"\n')
    f4.write('}\n')
    f4.write('close $outfile')
    f4.close()
    
    #Creando script para DB1_x1r
    psf=ruta_old_traj+'/'+psf_file
    dcd=ruta_old_traj+'/'+dcd_file
    print(psf)
    f5 = open('DB1_x1r.tcl', 'w')
    print(f5)
    f5.write('set psfFile '+ psf+' \n')
    f5.write('set dcdFile '+ dcd+' \n')
    f5.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f5.write('set outfile ' +'[open ' +'db1_x1r.dat'+' w]\n')
    f5.write('set nf [molinfo top get numframes]\n')
    f5.write('\n#RMSD calculation loop\n')
    f5.write('set f1 [atomselect top "'+DB1_i[4]+' or '+DB1_i[5]+' or '+DB1_i[6]+' or '+DB1_i[7]+' " frame 0]\n')
    f5.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f5.write('    set sel [atomselect top "'+DB1_i[4]+' or '+DB1_i[5]+' or '+DB1_i[6]+' or '+DB1_i[7]+' " frame $i]\n')
    f5.write('    $sel move [measure fit $sel $f1]\n')
    f5.write('    set time [expr {$i +1}]\n')
    f5.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f5.write('    puts $outfile "     $time"\n')
    f5.write('}\n')
    f5.write('close $outfile')
    f5.close()
if (revisa2>0):
    ##########################################################################################
    ## Creando los archivos para DB2
    #######################################################################################
    #Creando script para DB2_x1l
    psf=ruta_old_traj+'/'+psf_file
    dcd=ruta_old_traj+'/'+dcd_file
    print(psf)
    f6 = open('DB2_x1l.tcl', 'w')
    print(f6)
    f6.write('set psfFile '+ psf+' \n')
    f6.write('set dcdFile '+ dcd+' \n')
    f6.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f6.write('set outfile ' +'[open ' +'db2_x1l.dat'+' w]\n')
    f6.write('set nf [molinfo top get numframes]\n')
    f6.write('\n#RMSD calculation loop\n')
    f6.write('set f1 [atomselect top "'+DB2_i[0]+' or '+DB2_i[1]+' or '+DB2_i[2]+' or '+DB2_i[3]+' " frame 0]\n')
    f6.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f6.write('    set sel [atomselect top "'+DB2_i[0]+' or '+DB2_i[1]+' or '+DB2_i[2]+' or '+DB2_i[3]+' " frame $i]\n')
    f6.write('    $sel move [measure fit $sel $f1]\n')
    f6.write('    set time [expr {$i +1}]\n')
    f6.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f6.write('    puts $outfile "     $time"\n')
    f6.write('}\n')
    f6.write('close $outfile')
    f6.close()
    
    #Creando script para DB1_x2l
    psf=ruta_old_traj+'/'+psf_file
    dcd=ruta_old_traj+'/'+dcd_file
    print(psf)
    f7 = open('DB2_x2l.tcl', 'w')
    print(f7)
    f7.write('set psfFile '+ psf+' \n')
    f7.write('set dcdFile '+ dcd+' \n')
    f7.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f7.write('set outfile ' +'[open ' +'db2_x2l.dat'+' w]\n')
    f7.write('set nf [molinfo top get numframes]\n')
    f7.write('\n#RMSD calculation loop\n')
    f7.write('set f1 [atomselect top "'+DB2_i[1]+' or '+DB2_i[2]+' or '+DB2_i[3]+' or '+DB2_i[4]+' " frame 0]\n')
    f7.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f7.write('    set sel [atomselect top "'+DB2_i[1]+' or '+DB2_i[2]+' or '+DB2_i[3]+' or '+DB2_i[4]+' " frame $i]\n')
    f7.write('    $sel move [measure fit $sel $f1]\n')
    f7.write('    set time [expr {$i +1}]\n')
    f7.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f7.write('    puts $outfile "     $time"\n')
    f7.write('}\n')
    f7.write('close $outfile')
    f7.close()
    
    #Creando script para DB1_x3m
    psf=ruta_old_traj+'/'+psf_file
    dcd=ruta_old_traj+'/'+dcd_file
    print(psf)
    f8 = open('DB2_x3m.tcl', 'w')
    print(f8)
    f8.write('set psfFile '+ psf+' \n')
    f8.write('set dcdFile '+ dcd+' \n')
    f8.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f8.write('set outfile ' +'[open ' +'db2_x3m.dat'+' w]\n')
    f8.write('set nf [molinfo top get numframes]\n')
    f8.write('\n#RMSD calculation loop\n')
    f8.write('set f1 [atomselect top "'+DB2_i[2]+' or '+DB2_i[3]+' or '+DB2_i[4]+' or '+DB2_i[5]+' " frame 0]\n')
    f8.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f8.write('    set sel [atomselect top "'+DB2_i[2]+' or '+DB2_i[3]+' or '+DB2_i[4]+' or '+DB2_i[5]+' " frame $i]\n')
    f8.write('    $sel move [measure fit $sel $f1]\n')
    f8.write('    set time [expr {$i +1}]\n')
    f8.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f8.write('    puts $outfile "     $time"\n')
    f8.write('}\n')
    f8.write('close $outfile')
    f8.close()
    #Creando script para DB1_x2r
    psf=ruta_old_traj+'/'+psf_file
    dcd=ruta_old_traj+'/'+dcd_file
    print(psf)
    f9 = open('DB2_x2r.tcl', 'w')
    print(f9)
    f9.write('set psfFile '+ psf+' \n')
    f9.write('set dcdFile '+ dcd+' \n')
    f9.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f9.write('set outfile ' +'[open ' +'db2_x2r.dat'+' w]\n')
    f9.write('set nf [molinfo top get numframes]\n')
    f9.write('\n#RMSD calculation loop\n')
    f9.write('set f1 [atomselect top "'+DB2_i[3]+' or '+DB2_i[4]+' or '+DB2_i[5]+' or '+DB2_i[6]+' " frame 0]\n')
    f9.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f9.write('    set sel [atomselect top "'+DB2_i[3]+' or '+DB2_i[4]+' or '+DB2_i[5]+' or '+DB2_i[6]+' " frame $i]\n')
    f9.write('    $sel move [measure fit $sel $f1]\n')
    f9.write('    set time [expr {$i +1}]\n')
    f9.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f9.write('    puts $outfile "     $time"\n')
    f9.write('}\n')
    f9.write('close $outfile')
    f9.close()
    
    #Creando script para DB1_x1r
    psf=ruta_old_traj+'/'+psf_file
    dcd=ruta_old_traj+'/'+dcd_file
    print(psf)
    f10 = open('DB2_x1r.tcl', 'w')
    print(f10)
    f10.write('set psfFile '+ psf+' \n')
    f10.write('set dcdFile '+ dcd+' \n')
    f10.write('\nmol load psf $psfFile dcd $dcdFile\n')
    f10.write('set outfile ' +'[open ' +'db2_x1r.dat'+' w]\n')
    f10.write('set nf [molinfo top get numframes]\n')
    f10.write('\n#RMSD calculation loop\n')
    f10.write('set f1 [atomselect top "'+DB2_i[4]+' or '+DB2_i[5]+' or '+DB2_i[6]+' or '+DB2_i[7]+' " frame 0]\n')
    f10.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    f10.write('    set sel [atomselect top "'+DB2_i[4]+' or '+DB2_i[5]+' or '+DB2_i[6]+' or '+DB2_i[7]+' " frame $i]\n')
    f10.write('    $sel move [measure fit $sel $f1]\n')
    f10.write('    set time [expr {$i +1}]\n')
    f10.write('    puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
    f10.write('    puts $outfile "     $time"\n')
    f10.write('}\n')
    f10.write('close $outfile')
    f10.close()
Explanation: Creando los archivos tcl para el cálculo del RMSD de los puentes
Se crean los archivos de salida en formato tcl.
End of explanation
if (revisa1>0):
    #Calculando con VMD rmsd DB1 X1L
    !vmd -dispdev text < DB1_x1l.tcl
    #Calculando con VMD DB1 X2L
    !vmd -dispdev text < DB1_x2l.tcl
    #Calculando con VMD DB1 X3M
    !vmd -dispdev text < DB1_x3m.tcl
    #Calculando con VMD DB1 X2R
    !vmd -dispdev text < DB1_x2r.tcl
    #Calculando con VMD DB1 X1R
    !vmd -dispdev text < DB1_x1r.tcl
if (revisa2>0):
    #Calculando con VMD rmsd DB2 X1L
    !vmd -dispdev text < DB2_x1l.tcl
    #Calculando con VMD DB2 X2L
    !vmd -dispdev text < DB2_x2l.tcl
    #Calculando con VMD DB2 X3M
    !vmd -dispdev text < DB2_x3m.tcl
    #Calculando con VMD DB2 X2R
    !vmd -dispdev text < DB2_x2r.tcl
    #Calculando con VMD DB2 X1R
    !vmd -dispdev text < DB2_x1r.tcl
Explanation: Ejecutando los archivos rmsd en tcl con vmd
Ejecutando los archivos en VMD
End of explanation
escale_y=[]
fig = pl.figure(figsize=(25,8))
fig.subplots_adjust(hspace=.4, wspace=0.3)
#Formateando los valores de los ejes
#Engrosando marcos
ax = fig.add_subplot(2,5,1)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,2)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,3)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,4)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,5)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,6)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,7)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,8)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,9)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,10)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
if (revisa1>0):
    #Datos de DB1
    data_db1_x1l=np.loadtxt('db1_x1l.dat',comments=['#', '@'])
    data_db1_x2l=np.loadtxt('db1_x2l.dat',comments=['#', '@'])
    data_db1_x3m=np.loadtxt('db1_x3m.dat',comments=['#', '@'])
    data_db1_x2r=np.loadtxt('db1_x2r.dat',comments=['#', '@'])
    data_db1_x1r=np.loadtxt('db1_x1r.dat',comments=['#', '@'])
    sub1 = fig.add_subplot(251) # instead of plt.subplot(2, 2, 1)
   
    #sub1.set_title('DB1_X1L') 
    sub1.set_xlabel('Time (ns)')
    sub1.set_ylabel('RMSD (nm)')
    sub1.plot(data_db1_x1l[:,1]*0.02, data_db1_x1l[:,0]/10, color='black', linewidth = 1, label='DB1_X1L')
    x1,x2,y1,y2=sub1.axis()
    escale_y.append(y2)
    sub2 = fig.add_subplot(252)
    #sub2.set_title('DB1_X2L')
    sub2.set_xlabel('Time (ns)')
    sub2.set_ylabel('RMSD (nm)')
    sub2.plot(data_db1_x2l[:,1]*0.02, data_db1_x2l[:,0]/10, color='black', linewidth = 1, label='DB1_X2L')
    x1,x2,y1,y2=sub2.axis()
    escale_y.append(y2)
    sub3 = fig.add_subplot(253)
    #sub3.set_title('DB1_X3M')
    sub3.set_xlabel('Time (ns)')
    sub3.set_ylabel('RMSD (nm)')
    sub3.plot(data_db1_x3m[:,1]*0.02, data_db1_x3m[:,0]/10, color='black', linewidth = 1, label='DB1_X3M')
    x1,x2,y1,y2=sub3.axis()
    escale_y.append(y2)
    sub4 = fig.add_subplot(254)
    #sub4.set_title('DB1_X2R')
    sub4.set_xlabel('Time (ns)')
    sub4.set_ylabel('RMSD (nm)')
    sub4.plot(data_db1_x2r[:,1]*0.02, data_db1_x2r[:,0]/10, color='black', linewidth = 1, label='DB1_X2R')
    x1,x2,y1,y2=sub4.axis()
    escale_y.append(y2)
    sub5 = fig.add_subplot(255)
    #sub5.set_title('DB1_X1R')
    sub5.set_xlabel('Time (ns)')
    sub5.set_ylabel('RMSD (nm)')
    sub5.plot(data_db1_x1r[:,1]*0.02, data_db1_x1r[:,0]/10, color='black', linewidth = 1, label='DB1_X1R')
    x1,x2,y1,y2=sub5.axis()
    escale_y.append(y2)
if (revisa2>0):
    #DAtos de DB2
    data_db2_x1l=np.loadtxt('db2_x1l.dat',comments=['#', '@'])
    data_db2_x2l=np.loadtxt('db2_x2l.dat',comments=['#', '@'])
    data_db2_x3m=np.loadtxt('db2_x3m.dat',comments=['#', '@'])
    data_db2_x2r=np.loadtxt('db2_x2r.dat',comments=['#', '@'])
    data_db2_x1r=np.loadtxt('db2_x1r.dat',comments=['#', '@'])
    #Ploteando DB2
    sub6 = fig.add_subplot(256)
    #sub6.set_title('DB2_X1L')
    sub6.set_xlabel('Time (ns)')
    sub6.set_ylabel('RMSD (nm)')
    sub6.plot(data_db2_x1l[:,1]*0.02, data_db2_x1l[:,0]/10, color='red', linewidth = 1, label='DB2_X1L')
    x1,x2,y1,y2=sub6.axis()
    escale_y.append(y2)
    sub7 = fig.add_subplot(257)
    #sub7.set_title('DB2_X2L')
    sub7.set_xlabel('Time (ns)')
    sub7.set_ylabel('RMSD (nm)')
    sub7.plot(data_db2_x2l[:,1]*0.02, data_db2_x2l[:,0]/10, color='red', linewidth = 1, label='DB2_X2L')
    x1,x2,y1,y2=sub7.axis()
    escale_y.append(y2)
    sub8 = fig.add_subplot(258)
    #sub8.set_title('DB2_X3M')
    sub8.set_xlabel('Time (ns)')
    sub8.set_ylabel('RMSD (nm)')
    sub8.plot(data_db2_x3m[:,1]*0.02, data_db2_x3m[:,0]/10, color='red', linewidth = 1, label='DB2_X3M')
    x1,x2,y1,y2=sub8.axis()
    escale_y.append(y2)
    sub9 = fig.add_subplot(259)
    #sub9.set_title('DB2_X2R')
    sub9.set_xlabel('Time (ns)')
    sub9.set_ylabel('RMSD (nm)')
    sub9.plot(data_db2_x2r[:,1]*0.02, data_db2_x2r[:,0]/10, color='red', linewidth = 1, label='DB2_X2R')
    x1,x2,y1,y2=sub9.axis()
    escale_y.append(y2)
    sub10 = fig.add_subplot(2,5,10)
    #sub10.set_title('DB2_X1R')
    sub10.set_xlabel('Time (ns)')
    sub10.set_ylabel('RMSD (nm)')
    sub10.plot(data_db2_x1r[:,1]*0.02, data_db2_x1r[:,0]/10, color='red', linewidth = 1, label='DB2_X1R')
    x1,x2,y1,y2=sub10.axis()
    escale_y.append(y2)
    
#escale_y
escale_y.sort(reverse=True)
escale_y
##Cambiando los ejes de las y
sub1.axis((x1,x2,y1,escale_y[0]))
sub2.axis((x1,x2,y1,escale_y[0]))
sub3.axis((x1,x2,y1,escale_y[0]))
sub4.axis((x1,x2,y1,escale_y[0]))
sub5.axis((x1,x2,y1,escale_y[0]))
sub6.axis((x1,x2,y1,escale_y[0]))
sub7.axis((x1,x2,y1,escale_y[0]))
sub8.axis((x1,x2,y1,escale_y[0]))
sub9.axis((x1,x2,y1,escale_y[0]))
sub10.axis((x1,x2,y1,escale_y[0]))
Explanation: Generando los gráficos RMSD en matplotlib
End of explanation
### Creando el directorio para el análisis de las distancias de enlace de los puentes
ruta_diedros = nuevaruta+'/diedros_intra'
print ( ruta_diedros )
if not os.path.exists(ruta_diedros): 
    os.makedirs(ruta_diedros)
    print ('Se ha creado la ruta ===>',ruta_diedros)
else:
        print ("La ruta "+ruta_diedros+" existe..!!!")
print ( 'Nos vamos a ....', ruta_diedros)
os.chdir( ruta_diedros )
Explanation: FREE ENERGY  DIHEDRAL INTRAMOLECULAR
Se calculan las distancias de los ángulos diedros para el cálculo de la free energy intramolecular
End of explanation
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
if (revisa1>0):
  
    #Creando script para DB1_x1l
    d1 = open('dihed_DB1_x1l.tcl', 'w')
    print(d1)
    d1.write('set psfFile '+ psf+' \n')
    d1.write('set dcdFile '+ dcd+' \n')
    d1.write('\nmol load psf $psfFile dcd $dcdFile\n')
    d1.write('set outfile ' +'[open ' +'dihed_db1_x1l.dat'+' w]\n')
    d1.write('set nf [molinfo top get numframes]\n')
    d1.write(' \n')
    d1.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[0]+'"] get index]\n')
    d1.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[1]+'"] get index]\n')
    d1.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[2]+'"] get index]\n')
    d1.write('set selatoms4 [[atomselect top "protein and chain A and '+DB1_i[3]+'"] get index]\n')
    d1.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
    d1.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    d1.write('    set x [measure dihed $dihed frame $i]\n')
    d1.write('    set time [expr {$i +1}]\n')
    d1.write('    puts $outfile "$time     $x"\n')
    d1.write('}\n')
    d1.close()
    
    #Creando script para DB1_x2l
    d2 = open('dihed_DB1_x2l.tcl', 'w')
    print(d2)
    d2.write('set psfFile '+ psf+' \n')
    d2.write('set dcdFile '+ dcd+' \n')
    d2.write('\nmol load psf $psfFile dcd $dcdFile\n')
    d2.write('set outfile ' +'[open ' +'dihed_db1_x2l.dat'+' w]\n')
    d2.write('set nf [molinfo top get numframes]\n')
    d2.write(' \n')
    d2.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[1]+'"] get index]\n')
    d2.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[2]+'"] get index]\n')
    d2.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[3]+'"] get index]\n')
    d2.write('set selatoms4 [[atomselect top "protein and chain A and '+DB1_i[4]+'"] get index]\n')
    d2.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
    d2.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    d2.write('    set x [measure dihed $dihed frame $i]\n')
    d2.write('    set time [expr {$i +1}]\n')
    d2.write('    puts $outfile "$time     $x"\n')
    d2.write('}\n')
    d2.close()
       
    #Creando script para DB1_x3m
    d3 = open('dihed_DB1_x3m.tcl', 'w')
    print(d3)
    d3.write('set psfFile '+ psf+' \n')
    d3.write('set dcdFile '+ dcd+' \n')
    d3.write('\nmol load psf $psfFile dcd $dcdFile\n')
    d3.write('set outfile ' +'[open ' +'dihed_db1_x3m.dat'+' w]\n')
    d3.write('set nf [molinfo top get numframes]\n')
    d3.write(' \n')
    d3.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[2]+'"] get index]\n')
    d3.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[3]+'"] get index]\n')
    d3.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[4]+'"] get index]\n')
    d3.write('set selatoms4 [[atomselect top "protein and chain A and '+DB1_i[5]+'"] get index]\n')
    d3.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
    d3.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    d3.write('    set x [measure dihed $dihed frame $i]\n')
    d3.write('    set time [expr {$i +1}]\n')
    d3.write('    puts $outfile "$time     $x"\n')
    d3.write('}\n')
    d3.close()
    
    #Creando script para DB1_x2r
    d4 = open('dihed_DB1_x2r.tcl', 'w')
    print(d4)
    d4.write('set psfFile '+ psf+' \n')
    d4.write('set dcdFile '+ dcd+' \n')
    d4.write('\nmol load psf $psfFile dcd $dcdFile\n')
    d4.write('set outfile ' +'[open ' +'dihed_db1_x2r.dat'+' w]\n')
    d4.write('set nf [molinfo top get numframes]\n')
    d4.write(' \n')
    d4.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[3]+'"] get index]\n')
    d4.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[4]+'"] get index]\n')
    d4.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[5]+'"] get index]\n')
    d4.write('set selatoms4 [[atomselect top "protein and chain A and '+DB1_i[6]+'"] get index]\n')
    d4.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
    d4.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    d4.write('    set x [measure dihed $dihed frame $i]\n')
    d4.write('    set time [expr {$i +1}]\n')
    d4.write('    puts $outfile "$time     $x"\n')
    d4.write('}\n')
    d4.close()
    
    #Creando script para DB1_x1r
    d5 = open('dihed_DB1_x1r.tcl', 'w')
    print(d5)
    d5.write('set psfFile '+ psf+' \n')
    d5.write('set dcdFile '+ dcd+' \n')
    d5.write('\nmol load psf $psfFile dcd $dcdFile\n')
    d5.write('set outfile ' +'[open ' +'dihed_db1_x1r.dat'+' w]\n')
    d5.write('set nf [molinfo top get numframes]\n')
    d5.write(' \n')
    d5.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[4]+'"] get index]\n')
    d5.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[5]+'"] get index]\n')
    d5.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[6]+'"] get index]\n')
    d5.write('set selatoms4 [[atomselect top "protein and chain A and '+DB1_i[7]+'"] get index]\n')
    d5.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
    d5.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    d5.write('    set x [measure dihed $dihed frame $i]\n')
    d5.write('    set time [expr {$i +1}]\n')
    d5.write('    puts $outfile "$time     $x"\n')
    d5.write('}\n')
    d5.close()
    
if (revisa2>0):
    #####################################################################
    ########## Puente 2
    ##########################################3
    #Creando script para DB2_x1l
    d6 = open('dihed_DB2_x1l.tcl', 'w')
    print(d6)
    d6.write('set psfFile '+ psf+' \n')
    d6.write('set dcdFile '+ dcd+' \n')
    d6.write('\nmol load psf $psfFile dcd $dcdFile\n')
    d6.write('set outfile ' +'[open ' +'dihed_db2_x1l.dat'+' w]\n')
    d6.write('set nf [molinfo top get numframes]\n')
    d6.write(' \n')
    d6.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[0]+'"] get index]\n')
    d6.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[1]+'"] get index]\n')
    d6.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[2]+'"] get index]\n')
    d6.write('set selatoms4 [[atomselect top "protein and chain A and '+DB2_i[3]+'"] get index]\n')
    d6.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
    d6.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    d6.write('    set x [measure dihed $dihed frame $i]\n')
    d6.write('    set time [expr {$i +1}]\n')
    d6.write('    puts $outfile "$time     $x"\n')
    d6.write('}\n')
    d6.close()
    
    #Creando script para DB2_x2l
    d7 = open('dihed_DB2_x2l.tcl', 'w')
    print(d7)
    d7.write('set psfFile '+ psf+' \n')
    d7.write('set dcdFile '+ dcd+' \n')
    d7.write('\nmol load psf $psfFile dcd $dcdFile\n')
    d7.write('set outfile ' +'[open ' +'dihed_db2_x2l.dat'+' w]\n')
    d7.write('set nf [molinfo top get numframes]\n')
    d7.write(' \n')
    d7.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[1]+'"] get index]\n')
    d7.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[2]+'"] get index]\n')
    d7.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[3]+'"] get index]\n')
    d7.write('set selatoms4 [[atomselect top "protein and chain A and '+DB2_i[4]+'"] get index]\n')
    d7.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
    d7.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    d7.write('    set x [measure dihed $dihed frame $i]\n')
    d7.write('    set time [expr {$i +1}]\n')
    d7.write('    puts $outfile "$time     $x"\n')
    d7.write('}\n')
    d7.close()
       
    #Creando script para DB2_x3m
    d8 = open('dihed_DB2_x3m.tcl', 'w')
    print(d8)
    d8.write('set psfFile '+ psf+' \n')
    d8.write('set dcdFile '+ dcd+' \n')
    d8.write('\nmol load psf $psfFile dcd $dcdFile\n')
    d8.write('set outfile ' +'[open ' +'dihed_db2_x3m.dat'+' w]\n')
    d8.write('set nf [molinfo top get numframes]\n')
    d8.write(' \n')
    d8.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[2]+'"] get index]\n')
    d8.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[3]+'"] get index]\n')
    d8.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[4]+'"] get index]\n')
    d8.write('set selatoms4 [[atomselect top "protein and chain A and '+DB2_i[5]+'"] get index]\n')
    d8.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
    d8.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    d8.write('    set x [measure dihed $dihed frame $i]\n')
    d8.write('    set time [expr {$i +1}]\n')
    d8.write('    puts $outfile "$time     $x"\n')
    d8.write('}\n')
    d8.close()
    
    #Creando script para DB2_x2r
    d9 = open('dihed_DB2_x2r.tcl', 'w')
    print(d9)
    d9.write('set psfFile '+ psf+' \n')
    d9.write('set dcdFile '+ dcd+' \n')
    d9.write('\nmol load psf $psfFile dcd $dcdFile\n')
    d9.write('set outfile ' +'[open ' +'dihed_db2_x2r.dat'+' w]\n')
    d9.write('set nf [molinfo top get numframes]\n')
    d9.write(' \n')
    d9.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[3]+'"] get index]\n')
    d9.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[4]+'"] get index]\n')
    d9.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[5]+'"] get index]\n')
    d9.write('set selatoms4 [[atomselect top "protein and chain A and '+DB2_i[6]+'"] get index]\n')
    d9.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
    d9.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    d9.write('    set x [measure dihed $dihed frame $i]\n')
    d9.write('    set time [expr {$i +1}]\n')
    d9.write('    puts $outfile "$time     $x"\n')
    d9.write('}\n')
    d9.close()
    
    #Creando script para DB2_x1r
    d10 = open('dihed_DB2_x1r.tcl', 'w')
    print(d10)
    d10.write('set psfFile '+ psf+' \n')
    d10.write('set dcdFile '+ dcd+' \n')
    d10.write('\nmol load psf $psfFile dcd $dcdFile\n')
    d10.write('set outfile ' +'[open ' +'dihed_db2_x1r.dat'+' w]\n')
    d10.write('set nf [molinfo top get numframes]\n')
    d10.write(' \n')
    d10.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[4]+'"] get index]\n')
    d10.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[5]+'"] get index]\n')
    d10.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[6]+'"] get index]\n')
    d10.write('set selatoms4 [[atomselect top "protein and chain A and '+DB2_i[7]+'"] get index]\n')
    d10.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
    d10.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    d10.write('    set x [measure dihed $dihed frame $i]\n')
    d10.write('    set time [expr {$i +1}]\n')
    d10.write('    puts $outfile "$time     $x"\n')
    d10.write('}\n')
    d10.close()
Explanation: Creación de los archivos tcl para el cálculo de los ángulos diedros
End of explanation
if (revisa1>0):
    #Calculando con VMD rmsd DB1 X1L
    !vmd -dispdev text < dihed_DB1_x1l.tcl
    #Calculando con VMD DB1 X2L
    !vmd -dispdev text < dihed_DB1_x2l.tcl
    #Calculando con VMD DB1 X3M
    !vmd -dispdev text < dihed_DB1_x3m.tcl
    #Calculando con VMD DB1 X2R
    !vmd -dispdev text < dihed_DB1_x2r.tcl
    #Calculando con VMD DB1 X1R
    !vmd -dispdev text < dihed_DB1_x1r.tcl
if (revisa2>0):
    #Calculando con VMD rmsd DB2 X1L
    !vmd -dispdev text < dihed_DB2_x1l.tcl
    #Calculando con VMD DB2 X2L
    !vmd -dispdev text < dihed_DB2_x2l.tcl
    #Calculando con VMD DB2 X3M
    !vmd -dispdev text < dihed_DB2_x3m.tcl
    #Calculando con VMD DB2 X2R
    !vmd -dispdev text < dihed_DB2_x2r.tcl
    #Calculando con VMD DB2 X1R
    !vmd -dispdev text < dihed_DB2_x1r.tcl
print ('\nCopiando el archivo generateFES.py a '+ruta_diedros)
source_file=ruta_f_energy+'/generateFES.py'
dest_file=ruta_diedros+'/generateFES.py'
shutil.copy(source_file,dest_file)
#Cambiando permisos de ejecución
!chmod +x generateFES.py
Explanation: Ejecutando los archivos de los ángulos diedros tcl generados con VMD
End of explanation
if (revisa1>0):
    #Cargando valores del DB1_X1L
    data_db1_x1l=np.loadtxt('dihed_db1_x1l.dat',comments=['#', '@'])
    #Cargando valores del DB1_X1R
    data_db1_x1r=np.loadtxt('dihed_db1_x1r.dat',comments=['#', '@'])
    
    #Obteniendo los valores máximo y mínimo del DB1_X1L
    min_x1l=np.amin(data_db1_x1l[:,1])
    max_x1l=np.amax(data_db1_x1l[:,1])
    print ('Minimo DB1_X1L=>',min_x1l)
    print ('Máximo DB1_X1L=>',max_x1l)
    #Obteniendo los valores máximo y mínimo del DB1_X1R
    min_x1r=np.amin(data_db1_x1r[:,1])
    max_x1r=np.amax(data_db1_x1r[:,1])
    print ('Minimo DB1_X1R=>',min_x1r)
    print ('Máximo DB1_X1R=>',max_x1r)
    
    #Creando los archivos de entrada para el script
    np.savetxt('db1_x1l.dat',data_db1_x1l[:,1], fmt='%1.14f')
    np.savetxt('db1_x1r.dat',data_db1_x1r[:,1], fmt='%1.14f')
    !paste db1_x1l.dat db1_x1r.dat  > DB1_x1_lr.dat
    
    #Ejecutando el script de FES
    !python generateFES.py DB1_x1_lr.dat $min_x1l $max_x1l $min_x1r $max_x1r  200 200 $temperatura XL1_XR1.dat
    
    ###################################################################
    #Cargando valores del DB1_X2l
    data_db1_x2l=np.loadtxt('dihed_db1_x2l.dat',comments=['#', '@'])
    #Cargando valores del DB1_X1R
    data_db1_x2r=np.loadtxt('dihed_db1_x2r.dat',comments=['#', '@'])
    
    #Obteniendo los valores máximo y mínimo del DB1_X1L
    min_x2l=np.amin(data_db1_x2l[:,1])
    max_x2l=np.amax(data_db1_x2l[:,1])
    print ('Minimo DB1_X2L=>',min_x2l)
    print ('Máximo DB1_X2L=>',max_x2l)
    #Obteniendo los valores máximo y mínimo del DB1_X1R
    min_x2r=np.amin(data_db1_x2r[:,1])
    max_x2r=np.amax(data_db1_x2r[:,1])
    print ('Minimo DB1_X2R=>',min_x2r)
    print ('Máximo DB1_X2R=>',max_x2r)
    
    #Creando los archivos de entrada para el script
    np.savetxt('db1_x2l.dat',data_db1_x2l[:,1], fmt='%1.14f')
    np.savetxt('db1_x2r.dat',data_db1_x2r[:,1], fmt='%1.14f')
    !paste db1_x2l.dat db1_x2r.dat  > DB1_x2_lr.dat
    
    #Ejecutando el script de FES
    !python generateFES.py DB1_x2_lr.dat $min_x2l $max_x2l $min_x2r $max_x2r  200 200 $temperatura XL2_XR2.dat
    
    ######################################################################################
    #Generando los archivos para X3M
    data_db1_x3m=np.loadtxt('dihed_db1_x3m.dat',comments=['#', '@'])
    
    #Obteniendo los valores máximo y mínimo del DB1_X1L
    min_x3m=np.amin(data_db1_x3m[:,1])
    max_x3m=np.amax(data_db1_x3m[:,1])
    
    print ('Minimo DB1_X3M=>',min_x3m)
    print ('Máximo DB1_X3M=>',max_x3m)
    print ('Minimo DB1_X1L=>',min_x1l)
    print ('Máximo DB1_X1L=>',max_x1l)
    print ('Minimo DB1_X2L=>',min_x2l)
    print ('Máximo DB1_X2L=>',max_x2l)
    print ('Minimo DB1_X1R=>',min_x1r)
    print ('Máximo DB1_X1R=>',max_x1r)
    print ('Minimo DB1_X2R=>',min_x2r)
    print ('Máximo DB1_X2R=>',max_x2r)
    
    #Creando los archivos de entrada para el script
    np.savetxt('db1_x3m.dat',data_db1_x3m[:,1], fmt='%1.14f')
    !paste db1_x3m.dat db1_x1l.dat  > DB1_x3m_x1l.dat
    !paste db1_x3m.dat db1_x2l.dat  > DB1_x3m_x2l.dat
    !paste db1_x3m.dat db1_x1r.dat  > DB1_x3m_x1r.dat
    !paste db1_x3m.dat db1_x2r.dat  > DB1_x3m_x2r.dat
    
    #Ejecutando el script de FES
    !python generateFES.py DB1_x3m_x1l.dat $min_x3m $max_x3m $min_x1l $max_x1l  200 200 $temperatura XM3_XL1.dat
    !python generateFES.py DB1_x3m_x2l.dat $min_x3m $max_x3m $min_x2l $max_x2l  200 200 $temperatura XM3_XL2.dat
    !python generateFES.py DB1_x3m_x1r.dat $min_x3m $max_x3m $min_x1r $max_x1r  200 200 $temperatura XM3_XR1.dat
    !python generateFES.py DB1_x3m_x2r.dat $min_x3m $max_x3m $min_x2r $max_x2r  200 200 $temperatura XM3_XR2.dat
Explanation: Calculando la Free Energy Intramolecular para el Puente 1
End of explanation
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db1_xl1_vs_xr1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^L_1}"
set ylabel "{/=30 X@^R_1}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XL1_XR1.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db1_xl2_vs_xr2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^L_2}"
set ylabel "{/=30 X@^R_2}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XL2_XR2.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db1_xm3_vs_xl1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^L_1}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XM3_XL1.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db1_xm3_vs_xl2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^L_2}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XM3_XL2.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db1_xm3_vs_xr2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^R_2}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XM3_XR2.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db1_xm3_vs_xr1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^R_1}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XM3_XR1.dat" with pm3d
Explanation: Ploteando con GNUPLOT el Puente 1
End of explanation
if (revisa2>0):
    #Cargando valores del DB2_X1L
    data_db2_x1l=np.loadtxt('dihed_db2_x1l.dat',comments=['#', '@'])
    #Cargando valores del DB1_X1R
    data_db2_x1r=np.loadtxt('dihed_db2_x1r.dat',comments=['#', '@'])
    
    #Obteniendo los valores máximo y mínimo del DB2_X1L
    min_db2_x1l=np.amin(data_db2_x1l[:,1])
    max_db2_x1l=np.amax(data_db2_x1l[:,1])
    print ('Minimo DB2_X1L=>',min_db2_x1l)
    print ('Máximo DB2_X1L=>',max_db2_x1l)
    
    #Obteniendo los valores máximo y mínimo del DB2_X1R
    min_db2_x1r=np.amin(data_db2_x1r[:,1])
    max_db2_x1r=np.amax(data_db2_x1r[:,1])
    print ('Minimo DB2_X1R=>',min_db2_x1r)
    print ('Máximo DB2_X1R=>',max_db2_x1r)
    
    #Creando los archivos de entrada para el script
    np.savetxt('db2_x1l.dat',data_db2_x1l[:,1], fmt='%1.14f')
    np.savetxt('db2_x1r.dat',data_db2_x1r[:,1], fmt='%1.14f')
    !paste db2_x1l.dat db2_x1r.dat  > DB2_x1_lr.dat
    
    #Ejecutando el script de FES
    !python generateFES.py DB2_x1_lr.dat $min_db2_x1l $max_db2_x1l $min_db2_x1r $max_db2_x1r  200 200 $temperatura DB2_XL1_XR1.dat
    
    ###################################################################
    #Cargando valores del DB2_X2l
    data_db2_x2l=np.loadtxt('dihed_db2_x2l.dat',comments=['#', '@'])
    #Cargando valores del DB2_X1R
    data_db2_x2r=np.loadtxt('dihed_db2_x2r.dat',comments=['#', '@'])
    
    #Obteniendo los valores máximo y mínimo del DB2_X1L
    min_db2_x2l=np.amin(data_db2_x2l[:,1])
    max_db2_x2l=np.amax(data_db2_x2l[:,1])
    print ('Minimo DB2_X2L=>',min_db2_x2l)
    print ('Máximo DB2_X2L=>',max_db2_x2l)
    
    #Obteniendo los valores máximo y mínimo del DB2_X1R
    min_db2_x2r=np.amin(data_db2_x2r[:,1])
    max_db2_x2r=np.amax(data_db2_x2r[:,1])
    print ('Minimo DB2_X2R=>',min_db2_x2r)
    print ('Máximo DB2_X2R=>',max_db2_x2r)
    
    #Creando los archivos de entrada para el script
    np.savetxt('db2_x2l.dat',data_db2_x2l[:,1], fmt='%1.14f')
    np.savetxt('db2_x2r.dat',data_db2_x2r[:,1], fmt='%1.14f')
    !paste db2_x2l.dat db2_x2r.dat  > DB2_x2_lr.dat
    
    #Ejecutando el script de FES
    !python generateFES.py DB2_x2_lr.dat $min_db2_x2l $max_db2_x2l $min_db2_x2r $max_db2_x2r  200 200 $temperatura DB2_XL2_XR2.dat
    
    ######################################################################################
    #Cargando valores del DB2_X3M
    data_db2_x3m=np.loadtxt('dihed_db2_x3m.dat',comments=['#', '@'])
    
    #Obteniendo los valores máximo y mínimo del DB2_X3M
    min_db2_x3m=np.amin(data_db2_x3m[:,1])
    max_db2_x3m=np.amax(data_db2_x3m[:,1])
    
    print ('Minimo DB2_X3M=>',min_db2_x3m)
    print ('Máximo DB2_X3M=>',max_db2_x3m)
    
    print ('Minimo DB2_X1R=>',min_db2_x1r)
    print ('Máximo DB2_X1R=>',max_db2_x1r)
    
    print ('Minimo DB2_X2R=>',min_db2_x2r)
    print ('Máximo DB2_X2R=>',max_db2_x2r)
    
    print ('Minimo DB2_X1L=>',min_db2_x1l)
    print ('Máximo DB2_X1L=>',max_db2_x1l)
    
    print ('Minimo DB2_X2L=>',min_db2_x2l)
    print ('Máximo DB2_X2L=>',max_db2_x2l)
    
    #Creando los archivos de entrada para el script
    np.savetxt('db2_x3m.dat',data_db2_x3m[:,1], fmt='%1.14f')
    !paste db2_x3m.dat db2_x1r.dat  > DB2_x3m_x1r.dat
    !paste db2_x3m.dat db2_x2r.dat  > DB2_x3m_x2r.dat
    !paste db2_x3m.dat db2_x1l.dat  > DB2_x3m_x1l.dat
    !paste db2_x3m.dat db2_x2l.dat  > DB2_x3m_x2l.dat
    
    #Ejecutando el script de FES
    !python generateFES.py DB2_x3m_x1r.dat $min_db2_x3m $max_db2_x3m $min_db2_x1r $max_db2_x1r  200 200 $temperatura DB2_XM3_XR1.dat
    !python generateFES.py DB2_x3m_x2r.dat $min_db2_x3m $max_db2_x3m $min_db2_x2r $max_db2_x2r  200 200 $temperatura DB2_XM3_XR2.dat
    !python generateFES.py DB2_x3m_x1l.dat $min_db2_x3m $max_db2_x3m $min_db2_x1l $max_db2_x1l  200 200 $temperatura DB2_XM3_XL1.dat
    !python generateFES.py DB2_x3m_x2l.dat $min_db2_x3m $max_db2_x3m $min_db2_x2l $max_db2_x2l  200 200 $temperatura DB2_XM3_XL2.dat
Explanation: Calculando la Free Energy intramolecular para el Puente 2
End of explanation
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db2_xl1_vs_xr1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^L_1}"
set ylabel "{/=30 X@^R_1}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XL1_XR1.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db2_xl2_vs_xr2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^L_2}"
set ylabel "{/=30 X@^R_2}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XL2_XR2.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db2_xm3_vs_xl1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^L_1}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XM3_XL1.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db2_xm3_vs_xl2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^L_2}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XM3_XL2.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db2_xm3_vs_xr2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^R_2}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XM3_XR2.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db2_xm3_vs_xr1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^R_1}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XM3_XR1.dat" with pm3d
Explanation: Ploteando con GNUPLOT el Puente 2
End of explanation
############################################
#### Intermolecular DB1- DB2 - X1L
############################################
#Creando el DB1-DB2-X1L
!paste db1_x1l.dat db2_x1l.dat  > DB1_DB2_x1l.dat
print('Minimo DB1-X1L=>',min_x1l)
print('Máximo DB1-X1L=>',max_x1l)
print('Minimo DB2-X1L=>',min_db2_x1l)
print('Máximo DB2-X1L=>',max_db2_x1l)   
#Ejecutando el script de FES
!python generateFES.py DB1_DB2_x1l.dat $min_x1l $max_x1l $min_db2_x1l $max_db2_x1l  200 200 $temperatura DB1_DB2_X1L.dat
    
#########################################
#### Intermolecular DB1- DB2 - X2L
############################################
#Creando el DB1-DB2-X2L
!paste db1_x2l.dat db2_x2l.dat  > DB1_DB2_x2l.dat
print('Minimo DB1-X2L=>',min_x2l)
print('Máximo DB1-X2L=>',max_x2l)
print('Minimo DB2-X2L=>',min_db2_x2l)
print('Máximo DB2-X2L=>',max_db2_x2l)
#Ejecutando el script de FES
!python generateFES.py DB1_DB2_x2l.dat $min_x2l $max_x2l $min_db2_x2l $max_db2_x2l  200 200 $temperatura DB1_DB2_X2L.dat
############################################
#### Intermolecular DB1- DB2 - X3M
############################################
#Creando el DB1-DB2-X3M
!paste db1_x3m.dat db2_x3m.dat  > DB1_DB2_x3m.dat
print('Minimo DB1-X3M=>',min_x3m)
print('Máximo DB1-X3M=>',max_x3m)
print('Minimo DB2-X3M=>',min_db2_x3m)
print('Máximo DB2-X3M=>',max_db2_x3m)
#Ejecutando el script de FES
!python generateFES.py DB1_DB2_x3m.dat $min_x3m $max_x3m $min_db2_x3m $max_db2_x3m  200 200 $temperatura DB1_DB2_X3M.dat
############################################
#### Intermolecular DB1- DB2 - X2R
############################################
#Creando el DB1-DB2-X2R
!paste db1_x2r.dat db2_x2r.dat  > DB1_DB2_x2r.dat
print('Minimo DB1-X2R=>',min_x2r)
print('Máximo DB1-X2R=>',max_x2r)
print('Minimo DB2-X2R=>',min_db2_x2r)
print('Máximo DB2-X2R=>',max_db2_x2r)
#Ejecutando el script de FES
!python generateFES.py DB1_DB2_x2r.dat $min_x2r $max_x2r $min_db2_x2r $max_db2_x2r  200 200 $temperatura DB1_DB2_X2R.dat
############################################
#### Intermolecular DB1- DB2 - X1R
############################################
#Creando el DB1-DB2-X1R
!paste db1_x1r.dat db2_x1r.dat  > DB1_DB2_x1r.dat
print('Minimo DB1-X1R=>',min_x1r)
print('Máximo DB1-X1R=>',max_x1r)
print('Minimo DB2-X1R=>',min_db2_x1r)
print('Máximo DB2-X1R=>',max_db2_x1r)
#Ejecutando el script de FES
!python generateFES.py DB1_DB2_x1r.dat $min_x1r $max_x1r $min_db2_x1r $max_db2_x1r  200 200 $temperatura DB1_DB2_X1R.dat
Explanation: Free Energy Intermolecular
End of explanation
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "DB1_DB2_X1L.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 DB1 X@^L_1}"
set ylabel "{/=30 DB2 X@^L_1}"
set title "Free Energy Surface Intermolecular DB1-DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB1_DB2_X1L.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "DB1_DB2_X2L.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 DB1 X@^L_2}"
set ylabel "{/=30 DB2 X@^L_2}"
set title "Free Energy Surface Intermolecular DB1-DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB1_DB2_X2L.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "DB1_DB2_X3M.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 DB1 X@^M_3}"
set ylabel "{/=30 DB2 X@^M_3}"
set title "Free Energy Surface Intermolecular DB1-DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB1_DB2_X3M.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "DB1_DB2_X2R.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set xyplane 0
set pm3d interpolate 0,0
set xlabel "{/=30 DB1 X@^R_2}"
set ylabel "{/=30 DB2 X@^R_2}"
set title "Free Energy Surface Intermolecular DB1-DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB1_DB2_X2R.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "DB1_DB2_X1R.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 DB1 X@^R_1}"
set ylabel "{/=30 DB2 X@^R_1}"
set title "Free Energy Surface Intermolecular DB1-DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
set cbrange[8:10]
splot "DB1_DB2_X1R.dat" with pm3d
Explanation: Ploteando la Free Energy Intermolecular puentes DB1 y DB2
End of explanation
hist_escale_y=[]
fig = pl.figure(figsize=(25,8))
fig.subplots_adjust(hspace=.4, wspace=.3)
#subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)
#left  = 0.125  # the left side of the subplots of the figure
#right = 0.9    # the right side of the subplots of the figure
#bottom = 0.1   # the bottom of the subplots of the figure
#top = 0.9      # the top of the subplots of the figure
#wspace = 0.2   # the amount of width reserved for blank space between subplots
#hspace = 0.2   # the amount of height reserved for white space between subplots
#Formateando los valores de los ejes
#Engrosando marcos
ax = fig.add_subplot(2,5,1)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
    ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
ax = fig.add_subplot(2,5,2)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
    ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
ax = fig.add_subplot(2,5,3)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
    ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
ax = fig.add_subplot(2,5,4)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
    ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
ax = fig.add_subplot(2,5,5)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(3)
    ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
#Cargando valores del DB1
data_h_db1_x1l=np.loadtxt('db1_x1l.dat',comments=['#', '@'])
data_h_db1_x2l=np.loadtxt('db1_x2l.dat',comments=['#', '@'])
data_h_db1_x3m=np.loadtxt('db1_x3m.dat',comments=['#', '@'])
data_h_db1_x2r=np.loadtxt('db1_x2r.dat',comments=['#', '@'])
data_h_db1_x1r=np.loadtxt('db1_x1r.dat',comments=['#', '@'])
#Cargando valores del DB2
data_h_db2_x1l=np.loadtxt('db2_x1l.dat',comments=['#', '@'])
data_h_db2_x2l=np.loadtxt('db2_x2l.dat',comments=['#', '@'])
data_h_db2_x3m=np.loadtxt('db2_x3m.dat',comments=['#', '@'])
data_h_db2_x2r=np.loadtxt('db2_x2r.dat',comments=['#', '@'])
data_h_db2_x1r=np.loadtxt('db2_x1r.dat',comments=['#', '@'])
sub1 = fig.add_subplot(251) # instead of plt.subplot(2, 2, 1)
sub1.set_xlabel('Angle (Degree) ', fontsize=10)
sub1.set_ylabel('P(Angle)')
n1, bins1, rectangles1 = sub1.hist(data_h_db1_x1l,100, normed=True, color='black',histtype='step', linewidth=3)
n2, bins2, rectangles2 = sub1.hist(data_h_db2_x1l,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=sub1.axis()
hist_escale_y.append(y2)
sub2 = fig.add_subplot(252) # instead of plt.subplot(2, 2, 1)
sub2.set_xlabel('Angle (Degree) ', fontsize=10)
sub2.set_ylabel('P(Angle)')
n1, bins1, rectangles1 = sub2.hist(data_h_db1_x2l,100, normed=True, color='black',histtype='step', linewidth=3)
n2, bins2, rectangles2 = sub2.hist(data_h_db2_x2l,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=sub2.axis()
hist_escale_y.append(y2)
sub3 = fig.add_subplot(253) # instead of plt.subplot(2, 2, 1)
sub3.set_xlabel('Angle (Degree) ', fontsize=10)
sub3.set_ylabel('P(Angle)')
n1, bins1, rectangles1 = sub3.hist(data_h_db1_x3m,100, normed=True, color='black',histtype='step', linewidth=3)
n2, bins2, rectangles2 = sub3.hist(data_h_db2_x3m,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=sub3.axis()
hist_escale_y.append(y2)
sub4 = fig.add_subplot(254) # instead of plt.subplot(2, 2, 1)
sub4.set_xlabel('Angle (Degree) ', fontsize=10)
sub4.set_ylabel('P(Angle)')
n1, bins1, rectangles1 = sub4.hist(data_h_db1_x2r,100, normed=True, color='black',histtype='step', linewidth=3)
n2, bins2, rectangles2 = sub4.hist(data_h_db2_x2r,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=sub4.axis()
hist_escale_y.append(y2)
sub5 = fig.add_subplot(255) # instead of plt.subplot(2, 2, 1)
sub5.set_xlabel('Angle (Degree) ', fontsize=10)
sub5.set_ylabel('P(Angle)')
n1, bins1, rectangles1 = sub5.hist(data_h_db1_x1r,100, normed=True, color='black',histtype='step', linewidth=3)
n2, bins2, rectangles2 = sub5.hist(data_h_db2_x1r,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=sub5.axis()
hist_escale_y.append(y2)
#escale_y
hist_escale_y.sort(reverse=True)
hist_escale_y
##Cambiando los ejes de las y
sub1.axis((x1,x2,y1,hist_escale_y[0]))
sub2.axis((x1,x2,y1,hist_escale_y[0]))
sub3.axis((x1,x2,y1,hist_escale_y[0]))
sub4.axis((x1,x2,y1,hist_escale_y[0]))
sub5.axis((x1,x2,y1,hist_escale_y[0]))
Explanation: Calcular los histogramas de los diedros
End of explanation
### Creando el directorio para el análisis de las distancias de enlace de los puentes INTERMOLECULAR
ruta_bonds_puentes = nuevaruta+'/bonds_puentes'
print ( ruta_bonds_puentes )
if not os.path.exists(ruta_bonds_puentes): 
    os.makedirs(ruta_bonds_puentes)
    print ('Se ha creado la ruta ===>',ruta_bonds_puentes)
else:
        print ("La ruta "+ruta_bonds_puentes+" existe..!!!")
print ( 'Nos vamos a ....', ruta_bonds_puentes)
os.chdir( ruta_bonds_puentes )
Explanation: Ángulos de Enlace de los puentes Intermolecular
End of explanation
print ('\nCopiando el archivo generateFES.py a '+ruta_bonds_puentes)
source_file=ruta_scripts+'/free_energy/generateFES.py'
dest_file=ruta_bonds_puentes+'/generateFES.py'
shutil.copy(source_file,dest_file)
#Cambiando permisos de ejecución
!chmod +x generateFES.py
Explanation: Copiando el archivo de generación de FES
End of explanation
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print ('Puente DB1=>',DB1_N)
print ('Puente DB1=>',DB1_i)
print ('Puente DB2=>',DB2_N)
print ('Puente DB2=>',DB2_i)
puente=2
if (int(puente)==2):
  
    #Creando script para Bond X1 Left
    b1 = open('bond_DB1_left.tcl', 'w')
    print(b1)
    b1.write('set psfFile '+ psf+' \n')
    b1.write('set dcdFile '+ dcd+' \n')
    b1.write('\nmol load psf $psfFile dcd $dcdFile\n')
    b1.write('set outfile ' +'[open ' +'bond_db1_left.dat'+' w]\n')
    b1.write('set nf [molinfo top get numframes]\n')
    b1.write(' \n')
    b1.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[1]+'"] get index]\n')
    b1.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[2]+'"] get index]\n')
    b1.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[3]+'"] get index]\n')
    b1.write('set angle [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] ]\n')
    b1.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    b1.write('    set x [measure angle $angle frame $i]\n')
    b1.write('    set time [expr {$i +1}]\n')
    b1.write('    puts $outfile "$time     $x"\n')
    b1.write('}\n')
    b1.close()
    
    #Creando script para Bond X1 Right
    b2 = open('bond_DB1_right.tcl', 'w')
    print(b2)
    b2.write('set psfFile '+ psf+' \n')
    b2.write('set dcdFile '+ dcd+' \n')
    b2.write('\nmol load psf $psfFile dcd $dcdFile\n')
    b2.write('set outfile ' +'[open ' +'bond_db1_right.dat'+' w]\n')
    b2.write('set nf [molinfo top get numframes]\n')
    b2.write(' \n')
    b2.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[4]+'"] get index]\n')
    b2.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[5]+'"] get index]\n')
    b2.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[6]+'"] get index]\n')
    b2.write('set angle [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] ]\n')
    b2.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    b2.write('    set x [measure angle $angle frame $i]\n')
    b2.write('    set time [expr {$i +1}]\n')
    b2.write('    puts $outfile "$time     $x"\n')
    b2.write('}\n')
    b2.close()
    
    #Creando script para Bond DB2 X1 Left
    b3 = open('bond_DB2_left.tcl', 'w')
    print(b3)
    b3.write('set psfFile '+ psf+' \n')
    b3.write('set dcdFile '+ dcd+' \n')
    b3.write('\nmol load psf $psfFile dcd $dcdFile\n')
    b3.write('set outfile ' +'[open ' +'bond_db2_left.dat'+' w]\n')
    b3.write('set nf [molinfo top get numframes]\n')
    b3.write(' \n')
    b3.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[1]+'"] get index]\n')
    b3.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[2]+'"] get index]\n')
    b3.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[3]+'"] get index]\n')
    b3.write('set angle [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] ]\n')
    b3.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    b3.write('    set x [measure angle $angle frame $i]\n')
    b3.write('    set time [expr {$i +1}]\n')
    b3.write('    puts $outfile "$time     $x"\n')
    b3.write('}\n')
    b3.close()
    
    #Creando script para Bond DB2 X1 Right
    b4 = open('bond_DB2_right.tcl', 'w')
    print(b4)
    b4.write('set psfFile '+ psf+' \n')
    b4.write('set dcdFile '+ dcd+' \n')
    b4.write('\nmol load psf $psfFile dcd $dcdFile\n')
    b4.write('set outfile ' +'[open ' +'bond_db2_right.dat'+' w]\n')
    b4.write('set nf [molinfo top get numframes]\n')
    b4.write(' \n')
    b4.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[4]+'"] get index]\n')
    b4.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[5]+'"] get index]\n')
    b4.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[6]+'"] get index]\n')
    b4.write('set angle [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] ]\n')
    b4.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
    b4.write('    set x [measure angle $angle frame $i]\n')
    b4.write('    set time [expr {$i +1}]\n')
    b4.write('    puts $outfile "$time     $x"\n')
    b4.write('}\n')
    b4.close()
Explanation: Generando los archivos Tcl para el cálculo de los ángulos.
End of explanation
#Calculando con VMD bond DB1 Left
!vmd -dispdev text < bond_DB1_left.tcl
#Calculando con VMD bond DB1 Right
!vmd -dispdev text < bond_DB1_right.tcl
#Calculando con VMD bond DB2 Left
!vmd -dispdev text < bond_DB2_left.tcl
#Calculando con VMD bond DB2 Right
!vmd -dispdev text < bond_DB2_right.tcl
Explanation: Ejecutando los archivos tcl generados con VMD
End of explanation
#Cargando valores del DB1
data_bond_db1_left=np.loadtxt('bond_db1_left.dat',comments=['#', '@'])
#Cargando valores del DB1_X1R
data_bond_db1_right=np.loadtxt('bond_db1_right.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del DB1 Left
min_bond1_left=np.amin(data_bond_db1_left[:,1])
max_bond1_left=np.amax(data_bond_db1_left[:,1])
print ('Minimo DB1_Left=>',min_bond1_left)
print ('Máximo DB1_Left=>',max_bond1_left)
#Obteniendo los valores máximo y mínimo del DB1 Right
min_bond1_right=np.amin(data_bond_db1_right[:,1])
max_bond1_right=np.amax(data_bond_db1_right[:,1])
print ('Minimo DB1_Right=>',min_bond1_right)
print ('Máximo DB1_Right=>',max_bond1_right)
#Creando los archivos de entrada para el script
np.savetxt('bond_DB1_left.dat',data_bond_db1_left[:,1], fmt='%1.14f')
np.savetxt('bond_DB1_right.dat',data_bond_db1_right[:,1], fmt='%1.14f')
!paste bond_DB1_left.dat bond_DB1_right.dat  > angles_DB1.dat
#Ejecutando el script de FES
!python generateFES.py angles_DB1.dat $min_bond1_left $max_bond1_left $min_bond1_right $max_bond1_right  200 200 $temperatura Angles_DB1.dat
###################################################################3
#Cargando valores del DB2
data_bond_db2_left=np.loadtxt('bond_db2_left.dat',comments=['#', '@'])
#Cargando valores del DB1_X1R
data_bond_db2_right=np.loadtxt('bond_db2_right.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del DB2 Left
min_bond2_left=np.amin(data_bond_db2_left[:,1])
max_bond2_left=np.amax(data_bond_db2_left[:,1])
print ('Minimo DB2_Left=>',min_bond2_left)
print ('Máximo DB2_Left=>',max_bond2_left)
#Obteniendo los valores máximo y mínimo del DB2 Right
min_bond2_right=np.amin(data_bond_db2_right[:,1])
max_bond2_right=np.amax(data_bond_db2_right[:,1])
print ('Minimo DB2_Right=>',min_bond2_right)
print ('Máximo DB2_Right=>',max_bond2_right)
#Creando los archivos de entrada para el script
np.savetxt('bond_DB2_left.dat',data_bond_db2_left[:,1], fmt='%1.14f')
np.savetxt('bond_DB2_right.dat',data_bond_db2_right[:,1], fmt='%1.14f')
!paste bond_DB2_left.dat bond_DB2_right.dat  > angles_DB2.dat
#Ejecutando el script de FES
!python generateFES.py angles_DB2.dat $min_bond2_left $max_bond2_left $min_bond2_right $max_bond2_right  200 200 $temperatura Angles_DB2.dat
Explanation: Calculando la Free Energy de los Bonds de los puentes
End of explanation
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db1_a1_a2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 C@^1_{/Symbol a}}-{/=30 C@^1_{/Symbol b}}-{/=30 S@^1_{/Symbol g}}"
set ylabel "{/=30 C@^2_{/Symbol a}}-{/=30 C@^2_{/Symbol b}}-{/=30 S@^2_{/Symbol g}}"
set title "Free Energy Surface Angles DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "Angles_DB1.dat" with pm3d
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo  transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
%%gnuplot
set output "db2_a1_a2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
                      1 '#000fff',\
                      2 '#0090ff',\
                      3 '#0fffee',\
                      4 '#90ff70',\
                      5 '#ffee00',\
                      6 '#ff7000',\
                      7 '#ee0000',\
                      8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 C@^1_{/Symbol a}}-{/=30 C@^1_{/Symbol b}}-{/=30 S@^1_{/Symbol g}}"
set ylabel "{/=30 C@^2_{/Symbol a}}-{/=30 C@^2_{/Symbol b}}-{/=30 S@^2_{/Symbol g}}"
set title "Free Energy Surface Angles DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "Angles_DB2.dat" with pm3d
Explanation: Ploteando la Free Energy de los ángulos con gnuplot
End of explanation
bonds_escale_y=[]
#Cargando valores del DB1
data_h_db1_left=np.loadtxt('bond_DB1_left.dat',comments=['#', '@'])
data_h_db1_right=np.loadtxt('bond_DB1_right.dat',comments=['#', '@'])
#Cargando valores del DB2
data_h_db2_left=np.loadtxt('bond_DB2_left.dat',comments=['#', '@'])
data_h_db2_right=np.loadtxt('bond_DB2_right.dat',comments=['#', '@'])
#Engrosar marco    
figb=pl.figure(figsize=(12, 10), dpi=100, linewidth=3.0)
figb.subplots_adjust(hspace=.5)
ax = figb.add_subplot(221)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(4)
ax = figb.add_subplot(222)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(4)
ax = figb.add_subplot(223)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(4)
ax = figb.add_subplot(224)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
bond1 = figb.add_subplot(221) # instead of plt.subplot(2, 2, 1)
#bond1.set_title('CA1 - CB1 - SY1') 
bond1.set_xlabel('Angle (Degree)')
bond1.set_ylabel('P (Angle)')
n, bins, rectangles = bond1.hist(data_h_db1_left,100, normed=True, color='black',histtype='step',linewidth=3)
x1,x2,y1,y2=bond1.axis()
bonds_escale_y.append(y2)
bond2 = figb.add_subplot(222) # instead of plt.subplot(2, 2, 1)
#bond2.set_title('CA2 - CB2 - SY2') 
bond2.set_xlabel('Angle (Degree)')
bond2.set_ylabel('P (Angle)')
n, bins, rectangles = bond2.hist(data_h_db1_right,100, normed=True, color='black',histtype='step', linewidth=3)
x1,x2,y1,y2=bond2.axis()
bonds_escale_y.append(y2)
bond3 = figb.add_subplot(223) # instead of plt.subplot(2, 2, 1)
#bond3.set_title('CA1 - CB1 - SY1') 
bond3.set_xlabel('Angle (Degree)')
bond3.set_ylabel('P (Angle)')
n, bins, rectangles = bond3.hist(data_h_db2_left,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=bond3.axis()
bonds_escale_y.append(y2)
bond4 = figb.add_subplot(224) # instead of plt.subplot(2, 2, 1)
#bond4.set_title('CA2 - CB2 - SY2') 
bond4.set_xlabel('Angle (Degree)')
bond4.set_ylabel('P (Angle)')
n, bins, rectangles = bond4.hist(data_h_db2_right,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=bond4.axis()
bonds_escale_y.append(y2)
#escale_y
bonds_escale_y.sort(reverse=True)
bonds_escale_y
##Cambiando los ejes de las y
sub1.axis((x1,x2,y1,bonds_escale_y[0]))
sub2.axis((x1,x2,y1,bonds_escale_y[0]))
sub3.axis((x1,x2,y1,bonds_escale_y[0]))
sub4.axis((x1,x2,y1,bonds_escale_y[0]))
Explanation: Calculando los histogramas de los bonds
End of explanation
### Creando el directorio para el análisis de los puentes
ruta_clusters = nuevaruta+'/clusters'
print ( ruta_clusters )
if not os.path.exists(ruta_clusters): 
    os.makedirs(ruta_clusters)
    print ('Se ha creado la ruta ===>',ruta_clusters)
else:
        print ("La ruta "+ruta_clusters+" existe..!!!")
print ( 'Nos vamos a ....', ruta_clusters)
os.chdir( ruta_clusters )
    
Explanation: Generación de clusters
Crear la nueva ruta para calcular los clusters
End of explanation
 !echo 1 1 | g_cluster -f ../output.xtc -s ../ionized.pdb  -method gromos -cl out.pdb -g out.log -cutoff 0.2
Explanation: Calculando los clusters con la opción (1= Protein)
End of explanation
!vmd out.pdb
Explanation: Cargando los clusters para su visualización en VMD
Se cargan los clusters en VMD y se guardan sus coordenadas para cada uno de ellos haciendo uso de VMD
End of explanation
### Creando el directorio para el análisis de colorByRMSF
ruta_colorByRMSF = nuevaruta+'/colorByRMSF'
print ( ruta_colorByRMSF )
if not os.path.exists(ruta_colorByRMSF): 
    os.makedirs(ruta_colorByRMSF)
    print ('Se ha creado la ruta ===>',ruta_colorByRMSF)
else:
        print ("La ruta "+ruta_colorByRMSF+" existe..!!!")
print ( 'Nos vamos a ....', ruta_colorByRMSF)
os.chdir( ruta_colorByRMSF )
    
Explanation: colorByRMSF
Creando la carpeta para salida de datos
End of explanation
print ('\nCopiando el archivo colorByRMSF.vmd  a '+ruta_colorByRMSF)
source_file=ruta_scripts+'/colorByRMSF/colorByRMSF.vmd'
dest_file=ruta_colorByRMSF+'/colorByRMSF.vmd'
shutil.copy(source_file,dest_file)
Explanation: Copiando el archivo a la carpeta de datos
End of explanation
print ('Ejecutando el análisis de rmsf...')
!echo 1 | g_rmsf -f ../output.xtc -s ../ionized.pdb -oq bfac.pdb -o rmsf.xvg
#Calculando el mínimo y máximo del rmsf
#Cargando valores del RMSF
data_rmsf_gcolor=np.loadtxt('rmsf.xvg',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del RMSF
min_rmsf_gcolor=np.amin(data_rmsf_gcolor[:,1])
max_rmsf_gcolor=np.amax(data_rmsf_gcolor[:,1])
print ('Minimo_RMSF=>',min_rmsf_gcolor)
print ('Máximo_RMSF=>',max_rmsf_gcolor)
Explanation: Calculando el RMSF para el análisis de la proteína con la opción (1) Protein
End of explanation
# Cargando el pdb con VMD
!vmd ../ionized.pdb
Explanation: Cargar el scrit colorByRMSF.vmd en VMD
Arrancar VMD, dirigirse al menú Extensions -> Tk Console, copiar y ejecutar la siguiente secuencia de comandos en el cual pondremos los valores del Mínimo_RMSF y Máximo_RMSF calculado en la celda anterior:
tcl
source colorByRMSF.vmd
colorByRMSF top rmsf.xvg Mínimo_RMSF Máximo_RMSF
ESCALA DE COLOR
Dirigirse al menú Extensions -> Visualization -> Color Scale Bar y cambiar los valores de los siguientes campos:
1. Colocar el valor calculado de Mínimo_RMSF en el campo Mínimum scale value
2. Colocar el valor calculado de Máximo_RMSF en el campo Maximum scale value.
3. Seleccionar el color Black en el campo Color of labels.
CAMBIAR EL COLOR DE FONDO
Dirigirse al menú Graphics -> Colors , y realizar las siguientes selecciones:
1. Categories seleccionar Display
2. Names  seleccionar Background
3. Colors seleccionar 8 White
REMOVER EJE X,Y,Z
Dirigirse al menú Display -> Axes -> Off, con el cual eliminaremos el eje de X,Y,Z.
End of explanation
print ( 'Nos vamos a ....', ruta_colorByRMSF )
os.chdir( ruta_colorByRMSF )
Explanation: Graficando B-Factors con Chimera
End of explanation
#Inicializando vector
rmsf=[]
rmsf_x=[]
rmsf_y=[]
try:
    file_Bfactor = open( 'bfac.pdb' )
    new_bfactor=open('bfac_new.pdb','w')
    
except IOError:
     print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in file_Bfactor.readlines():
    fila = linea.strip()
    sl = fila.split()
    cadena=sl[0]
    if (cadena=='ATOM'):
        if (len(sl)==12):
            new_bfactor.write(linea)
        else:
            x=linea[0:60]
            tempFactor=linea[60:66]
            #print (x)
            #print(tempFactor)
            y=fila[67:]
            #print (y)
            enviar=x+' '+tempFactor+y
            new_bfactor.write(enviar+'\n')
            #print(enviar)
            
    else:
        #print (linea)
        new_bfactor.write(linea)
new_bfactor.close()
Explanation: Adecuando archivo bfac.pdb para obtener la columna de B-factors
End of explanation
!gedit bfac_new.pdb
Explanation: Revisando la estructura del archivo generado.
Revisar que los campos se encuentren completamente alineados en la estructura de los campos.
Guardar y salir.
End of explanation
#Inicializando vector
bfactors_color=[]
try:
     file_bfactor_color = open( 'bfac_new.pdb' )
except IOError:
     print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in file_bfactor_color.readlines():
    fila = linea.strip()
    sl = fila.split()
    if (sl[0]=='ATOM'):
        #print (sl[0])
        idresidue=fila[23:26]
        bfactor=fila[60:66]
        #print (idresidue + '\t'+bfactor)
        bfactors_color.append(idresidue+'\t'+bfactor+'\n')
        #i=i+1
#Escribiendo el archivo BFACTOR.dat
f = open('protein_bfactor.dat', 'w')
#f.write('@    title "B-factors" \n')
f.write('@    xaxis  label " Residue" \n')
f.write('@    xaxis  label char size 1.480000\n')
f.write('@    xaxis  bar linewidth 5.0\n')
f.write('@    xaxis  ticklabel char size 1.480000\n')
f.write('@    yaxis  label "B-factors (' +"\\"+'cE'+"\\"+'C)"\n')
f.write('@    yaxis  label char size 1.480000\n')
f.write('@    yaxis  bar linewidth 5.0\n')
f.write('@    yaxis  ticklabel char size 1.480000\n')
f.write('@    s0 line linewidth 7\n')
f.write('@    s0 symbol 1\n')
f.write('@    s0 symbol size 1.000000\n')
f.write('@    s0 symbol color 1\n')
f.write('@    s0 symbol pattern 1\n')
f.write('@    s0 symbol fill color 2\n')
f.write('@    s0 symbol fill pattern 1\n')
f.write('@    s0 symbol linewidth 1.0\n')
f.write('@TYPE xy \n')
f.write("".join(bfactors_color))
f.close()
        
!xmgrace protein_bfactor.dat
#Cargando la imagen generada en xmgrace
Image(filename='protein_bfactor.png')
#Calculando el mínimo y máximo del rmsf
#Cargando valores del RMSF
data_bfactor_color=np.loadtxt('protein_bfactor.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del RMSF
min_bfactor_color=np.amin(data_bfactor_color[:,1])
max_bfactor_color=np.amax(data_bfactor_color[:,1])
print ('Minimo_B-Factor=>',min_bfactor_color)
print ('Máximo_B-Factor=>',max_bfactor_color)
Explanation: Generando el archivo de Bfactors para todos los átomos FALTA ADECUAR PARA SACAR EL MAYOR POR RESIDUO
End of explanation
!chimera bfac_new.pdb
Explanation: Cargando el archivo pdb con Chimera para realizar la coloración de Bfactors
End of explanation
##Cargando la imagen generada
print ('Cargando el archivo...')
Image(filename='image.png') 
Explanation: Instrucciones para generar la imagen de B-factors
 ESTABLECER EL MODO DE VISUALIZACIÓN
1. Seleccionar del menú principal Presets -> Interactive 2 (all atoms).
2. Seleccionar del menú principal Actions -> Surface -> Show.
3. Ajustar el tamaño de la ventana principal.
4. Ajustar el tamaño y posición de la figura haciendo uso de la tecla CTRL+ Button wheel mouse.
 COLOREAR LOS B-FACTORS 
Seleccionar  Tools -> Depiction -> Render by Attribute.
Nos desplegará una ventana Render/Select by Attribute.
  1. Del campo Attribute seleccionar bfactor.
  2. En el histograma que se muestra, seleccionar la barra blanca y cambiar el color de blanco a amarillo en el campo color.
  3. Pulsar el botón Apply para visualizar los cambios de coloración.
  4. Pulsar OK para finalizar.
 FONDO BLANCO
Para aplicar el fondo blanco:
1.  Seleccionar del menú principal Presets->Publication_1.
 SALVAR LA POSICIÓN DE LA IMAGEN
Una vez que se ha obtenido la imagen coloreada, ajustar la visualización rotando la imagen, con la finalidad de dejar los espacios adecuados para la inclusión de las etiquetas y la barra de color.
Para salvar la posición final de la imagen:
  1. Seleccionar Favorites -> Command Line.
  2. En la línea de comando teclear savepos p1
Si por alguna razón movemos la posición, para restaurarla hacer lo siguiente:
  1. Seleccionar Favorites -> Command Line.
  2. En la línea de comando teclear reset p1 
 TITULO Y BARRA DE COLOR 
Seleccionar del menú principal Tools -> Utilities -> Color Key. El cual desplegará la ventana 2D Labels/Color Key.
Para desplegar la barra de color:
Seleccionar la pestaña Color Key.
Cambiar el color blanco por amarillo pulsando en el botón correspondiente.
Cambiar la palabra min por el valor mínimo calculado del bfactor.
Cambiar la palabra max por el valor máximo calculado del bfactor.
Dar click con el mouse en la parte inferior de la imagen en donde se desea visualizar la escala. Arrastrar el mouse para definir el largo y ancho de la escala.
Para desplegar el título de la barra:
Seleccionar la pestaña  Labels.
Dar click en la parte superior de barra de color para incrustar el título.
Escribir el título de la barra como B-Factors(Å).
Para ajustar el tamaño de letra, en el campo Font size cambiar el valor adecuado.
Para desplegar el título de la imagen:
Seleccionar la pestaña  Labels.
Dar click en la parte superior de la imagenr para incrustar el título.
Escribir el título con el nombre correspondiente.
Ajustar el tamaño de letra, en el campo Font size cambiar el valor adecuado.
Para el título en negrita, en el campo Font style seleccionar bold.
Notas:
1. Si desea cambiar una etiqueta de posición, deberá estar en la pestaña Labels, mantener pulsado el botón izquierdo del mouse sobre la etiqueta y moverla a la posición deseada.
2. Si desea eliminar una etiqueta, deberá seleccionarla en el campo de Labels y desmarcar la opción Show.
 SALVAR LA IMAGEN 
Seleccionar del menú principal File -> Save Image. El cual desplegará la ventana Save image, en el cual en el campo File name dar el nombre de image.png.
SALVAR LA SESIÓN DE QUIMERA
Seleccionar del menú principal File -> Save Session as. El cual desplegará la ventana Choose Session Save File, en el cual en el campo File name colocar el nombre con la extensión .py.
End of explanation
### Creando el directorio para el análisis del SASA en el directorio de VMD
print ('Nos vamos a ', ruta)
os.chdir( ruta )
output_find=!find   /usr/local -maxdepth 2 -type d -name vmd
print (output_find)
ruta_vmd=output_find[0]
print (ruta_vmd)
ruta_vmd_sasa = ruta_vmd+'/plugins/noarch/tcl/iceVMD1.0'
print ( ruta_vmd_sasa )
if not os.path.exists(ruta_vmd_sasa): 
    os.makedirs(ruta_vmd_sasa)
    print ('Se ha creado la ruta ===>',ruta_vmd_sasa)
else:
        print ("La ruta "+ruta_vmd_sasa+" existe..!!!")
print ( 'Nos vamos a ....', ruta_vmd_sasa )
os.chdir( ruta_vmd_sasa )
#Copiando los archivos generados a la carpeta plugins de VMD
print ('\nCopiando los archivos generados a '+ruta_vmd_sasa)
source_file=ruta_scripts+'/iceVMD1.0/colorplot.tcl'
dest_file=ruta_vmd_sasa+'/colorplot.tcl'
shutil.copy(source_file,dest_file)
source_file=ruta_scripts+'/iceVMD1.0/multiplot.tcl'
dest_file=ruta_vmd_sasa+'/multiplot.tcl'
shutil.copy(source_file,dest_file)
source_file=ruta_scripts+'/iceVMD1.0/pkgIndex.tcl'
dest_file=ruta_vmd_sasa+'/pkgIndex.tcl'
shutil.copy(source_file,dest_file)
source_file=ruta_scripts+'/iceVMD1.0/vmdICE.tcl'
dest_file=ruta_vmd_sasa+'/vmdICE.tcl'
shutil.copy(source_file,dest_file)
print('\nArchivos copiados.. Regresando a... '+nuevaruta)
os.chdir( nuevaruta )
### Creando el directorio para la graficación del sasa
ruta_sasaColor = nuevaruta+'/sasaColor'
print ( ruta_sasaColor )
if not os.path.exists(ruta_sasaColor): 
    os.makedirs(ruta_sasaColor)
    print ('Se ha creado la ruta ===>',ruta_sasaColor)
else:
        print ("La ruta "+ruta_sasaColor+" existe..!!!")
    
print ( 'Nos vamos a ....', ruta_sasaColor )
os.chdir( ruta_sasaColor )
print ('\nCopiando el archivo de configuracion a '+ruta_sasaColor)
source_file=ruta_scripts+'/iceVMD1.0/vmdrc'
dest_file=ruta_sasaColor+'/.vmdrc'
shutil.copy(source_file,dest_file)
Explanation: Graficando SASA
End of explanation
!vmd ../ionized.psf ../output.xtc
Explanation: Coloreando el SASA
Arrancar VMD.
Ventana vmdICE
Dirigirse al menú Extensions -> Analysis -> vmdICE, se presentará una ventana y se deberán cambiar los valores de los siguientes campos:
1. To: Colocar el rango máximo de frames de la trayectoria.
2. Selection for Calculation: agregar a chain A and protein.
3. Pulsar en el botón SASA Single Atom y esperar a que termine el cálculo.
CAMBIAR EL COLOR DE FONDO
Dirigirse al menú Graphics -> Colors , y realizar las siguientes selecciones:
1. Categories seleccionar Display
2. Names  seleccionar Background
3. Colors seleccionar 8 White
CAMBIAR RESOLUCIÓN DE ESFERAS
Dirigirse al menú Graphics - Representations, y en el campo Sphere Resolution cambiamos al valor de 50.
ROTAR LA IMAGEN PARA PRESENTAR UNA MEJOR VISTA Y GUARDARLA.
End of explanation
#Borrando los archivos del vmd
!rm -r $ruta_vmd_sasa
Explanation: Restaurando configuración default de VMD
End of explanation
### Creando el directorio para la graficación del rgyro
ruta_gyroColor = nuevaruta+'/color_rgyro'
print ( ruta_gyroColor )
if not os.path.exists(ruta_gyroColor): 
    os.makedirs(ruta_gyroColor)
    print ('Se ha creado la ruta ===>',ruta_gyroColor)
else:
        print ("La ruta "+ruta_gyroColor+" existe..!!!")
    
print ( 'Nos vamos a ....', ruta_gyroColor )
os.chdir( ruta_gyroColor )
print ('\nCopiando el script colorRgyro.tcl a '+ruta_gyroColor)
source_file=ruta_scripts+'/colorRgyro/colorRgyro.tcl'
dest_file=ruta_gyroColor+'/colorRgyro.tcl'
shutil.copy(source_file,dest_file)
Explanation: Graficando el RGYRO
End of explanation
!vmd ../ionized.psf ../output.xtc
Explanation: Coloreando el RGYRO
Arrancar VMD, dirigirse al manú Extensions -> Tk Console, copiar y ejecutar la siguiente secuencia de comandos:
tcl
source colorRgyro.tcl
CAMBIAR EL COLOR DE FONDO
Dirigirse al menú Graphics -> Colors , y realizar las siguientes selecciones:
1. Categories seleccionar Display
2. Names  seleccionar Background
3. Colors seleccionar 8 White
ROTAR LA IMAGEN PARA PRESENTAR UNA MEJOR VISTA Y GUARDARLA.
End of explanation
<END_TASK>