利用gromacs软件绘制自由能形貌图

一、准备gromacs能够识别的文件：

模型的拓扑文件和坐标文件，在amber中，使用的后缀分别是*.prmtop和*.crd

1.利用python的parm模块（要是装有amber，就有这个模块）进行top和crd文件格式转换：

编写如在的格式转换py脚本convert.py：

#!/usr/bin/python3
import parmed as pmd
import sys

# convert GROMACS to AMBER :
def conver_gromacs_to_amber(top,crd):
	gmx_top=pmd.load_file(top,xyz=crd)
	gmx_top.save("myfile.prmtop",format="amber")
	gmx_top.save("myfile.inpcrd",format="rst7")


# convert amber to gromacs:
def convert_amber_to_gromacs(top,crd):
	amber=pmd.load_file(top,crd)
	amber.save("myfile.top")
	amber.save("myfile.gro")

if __name__ == "__main__":
	
	if sys.argv[1]=="Y":
		amber_crd_file=sys.argv[2]
		amber_top_file=sys.argv[3]
		convert_amber_to_gromacs(amber_top_file,amber_crd_file)
	else:
		#such as "pmaawaterFE20mer2.top"
		gmx_top_file=sys.argv[1]
		#such as "pmaawaterFE20mer2.gro"
		gmx_crd_file=sys.argv[2]
		convert_amber_to_gromacs(gmx_top_file,gmx_crd_file)

输入如下的命令，将amber的top和crd文件转换为gromacs格式：

python3 convert.py Y xx.inpcrd xx.prmtp

2.运行vmd软件将轨迹文件*nc转化为gromac支持的格式trr

二、准备绘制FEL的数据:

1.利用gromacs计算轨迹的RMSD和Rg数据：

我这里使用的坐标文件是转换脚本得到的xx.gro，轨迹文件是VMD得到的xx.trr

#去除整体的平动和转动：
#*.tpr是坐标文件，*.xtc是轨迹文件
gmx trjconv –s pro.tpr –f pro.xtc –o pro.xtc –fit rot+trans

#计算CA的RMSD值，输出数据为rmsd.xvg:
gmx rms –s pro.tpr –f pro.xtc –o rmsd.xvg 

#计算CA的Rg值，输出数据为gyrate.xvg:
gmx gyrate –s pro.tpr –f pro.xtc –o gyrate.xvg

这部分可以用cpptraj来做，*.xvg文件用gedit文本编辑器可以打开

2.将RMSD数据和Rg数据整合为一个文件out.xvg:

编写PY脚本combine.py：

import pandas as pd

#read data:
rg=pd.read_csv(r'tyr_rg.xvg',sep=",")
rms=pd.read_csv(r'tyr_rmsd.xvg',sep=",")

#extract specified data:
frame=rg["frame"]
rmsd=rms["rmsd"]
Rg=rg["rg"]

#整合数据：
with open ("out.xvg",mode="w+",encoding="utf-8") as f:
      for i in range(len(frame)):
            f.write(str(frame[i])+","+str(rmsd[i])+","+str(Rg[i])+"\n")

在运行脚本前，需要将输入文件*.xvg格式改为逗号分隔的格式，同时添加行名

得到用逗号分隔的输出数据文件out.xvg

out.xvg数据列之间用7个空格分割，第一行信息用#号注释

3.gromacs的sham绘制自由能形貌图： 

gmx sham –tsham 310 –nlevels 100 –f out.xvg –g gibbs.log –ls gibbs.xpm –g gibbs.log \ 
–lsh  enthalpy.xpm -lss entropy.xpm
#或是：
gmx sham -f out.xvg

这里得到的由gromacs得到的四张xpm图片是灰度图 ，利用脚本xpm2png.py将图片转为彩色的png

脚本xpm2png.py的下载：xpm2png.py下载链接

输入命令，利用xpm2png.py，将灰度图转换成彩色图：

python xpm2png.py -f Gibbs.xpm -o Gibbs.png -ip yes -show yes

不加-ip yes 输出的还是灰度图，但-show可加可不加，得到的Gibbs.png

查看得到的gibbs.png：

​

三、FEL图片细节的调整

1.对图像标签的调整： 

因为数据不是主成分，而是RMSD和Rg，就在Gibbs.xpm里修改x-label和y-label ：

​

2.对图像色板的调整：

若是要调整色板和图像距离，在xpm2png.py的代码段plt.colorbar()上加上pad=0.12，这里的<float>是一个浮点数，表示右侧色板和中央能量形貌图的距离。

若是要调整色板为水平状态，在plt.colorbar()里加上orientation="horizontal"，将色板转为水平状态。

想要划分出色板渐变的颜色，在plt.colorbar()里加上boundaries=[i\10 for i in range(1,61)]，这里的list是色板上的刻度值。想要为色板渐变的颜色添加划线，同样在plt.colorbar()里加上参数drawedges=True

一般最好不要修改这个脚本。

查看修改的效果: 

  ​

图中只有一个能量井，说明这个轨迹里出现一个能量最低的稳定构象