Xmgrace学习笔记
已有 3466 次阅读
|系统分类:
整理: 李卫星对于用惯Linux的人, 在Linux下画图有点不方便, gnuplot虽然挺不错, 但是不能直接在图上进行操作, 这是它的不足吧. 为此, 简单介绍一下xmgrace.xmgrace可以直接用命令来安装：sudo apt-get install xmgrace也可以下载软件包自己编译. 官网.有一个整理好的压缩包, 里面包含了xmgrace–5.1.25的源代码以及一些资料.grace和xmgrace差不多, 具体区别自己看网上.几个网址：, 个人觉得这个是不错的下载上面的grace–5.1.25后, Grace Tutorial教程实例文件在grace-5.1.25/doc下面（Tutorial.pdf也在）. 还有一些做好的图例在grace-5.1.25/examples, 从那里可以看到xmgrace也是可以的. 下面是随便挑的几张, 图片的的字体都是用xmgrace加上去的.本文档主要是学习Grace Tutorials的总结, 如有不对, 请指正.可以在工作目录下打开终端, 输入xmgrace file.如果数据文件里面是多列数据, 即X Y1 Y2…, 可以输入xmgrace -nxy file. 加上-nxy就可以是多条曲线了.如果要对坐标轴取对数, 可以添加选项 -log x|y|xy, 如：xmgrace -log x -nxy file. （好像要放前面, 放后面会出错！）也可以在终端只输入xmgrace, 不带文件名, 打开空白的xmgrace再自己导入. 点击菜单Data | Import | ASCII在files下拉框选择文件, 点一下OK就可以了.注意, 好像只有扩展名为.dat的才显示, 而.xvg格式的不显示, 需要自己改一下扩展名. 还有就是, 如果文件路径有中文的话显示会出现乱码, 但仍然可以选择, 只要你确定哪个是你的文件.也可以自己用表格创建绘图数据. 在终端输入xmgrace, 打开空白的xmgrace. 点击菜单Main | Edit/Data_sets...窗口, 在上面的空白框中, 一直右击按住出现菜单移动（右击还不要放开）, 选creat new选in spreadsheet（这里还有另外三个, 如编写公式输入, 你可以自己摸索或看手册）打开表格, 就可以输入x y点数值了.画好了图片就保存, 点击File | print_setup, 在device选格式, 我常选png, 对不起的是好像没有tiff格式可选. 输入想要保存的名字, 点aceept退出窗口, 按一下快捷键ctrl+p就保存了. 其实也可以用File | print, 但最好记住快捷键. 看看目录内是否有图片文件.xmgrace不像origin那样有什么工程文件, 它的设置不会自动保存, 所以你上一次对数据的操作是不会保存的. 这很麻烦. 所以如果我们以后还要对文件进行操作（如坐标轴, 曲线粗细不够的等问题）, 很麻烦, 其他的设置要重新做一遍. 为避免这些麻烦, 就得保存设置. 点击菜单File | save或File | save as. File | save直接把设置加入到你的数据文件开头, 自己操作完可以去看看. File | save as就是另存一份, 原数据文件开头添加设置. 如果开头保存有设置的话, 直接右击用xmgrace查看, 就可以看见你原来设置的图片. 你修改后, 重新File | save. 这样也有一些好处, 如果多个文件都是一样设置, 做好一个文件, 把设置保存到数据文件开头, 然后把它们复制到其他文件, 就搞定了, 然后再用xmgrace输出图片.你可以去grace-5.1.25/exemple打开各个文件学习一下各种图是如何进行设置的.第一点说过如何使用-nxy选项绘制多列数据, 但如果只想用第1和第3列数据画图, 咋办? 可以使用xmgrace -block file.dat -bxy 1:3. 其中-bxy就是选1和3两列. 也可以在-nxy打开多条曲线后, 点击图上的曲线, 出现set appearance窗口, 在select set下, 选中不要的数据列, 长按右击移动选hide就可以了（数据列多的话, 就麻烦一些）;说到这里, 说一下对文件的数据处理, 看示例：xmgrace -nosafe -nxy box.xvg -pexec &s0.y=(s0.y*s1.y)/64& -pexec &kill s1& -pexec &autoscale&我的文件box.xvg中有三列数据, 时间帧、盒子x大小、盒子y大小, 计算膜表面的APL(area prea lipid)是用x*y/64. 用xmgrace怎么计算呢? 像上面那样用-pexec选项输入参数命令就可以了. 输出的图像就是x轴是时间, y轴是经过处理的得到的APL. 更具体的请看手册教程.还是再说一些多列数据的情况. data | Import/ASCII选好dat文件后在中间那里看到load as了吧, 点一下选block data, 后点ok跳出来框框, x from column自己选y from column, 选好就可以. 如果都选1, 会是什么图像, 猜一下, xy都是一样的值, 当然是45度斜率的直线y=x. 那里可以选偏差条, 就是xydy, 自己摸索了.数据选择操作说的差不多, 说一下图片的外观设置（也就是具体的坐标轴文字大小、间距、范围, 图片标题, 线条粗、细颜色等）. 我们可以双击相应的地方（和origin类似, 自己体会）也可以菜单里选择plot | plot appearance.说一下一些基本的操作：Legends: 图例, 就是多条曲线时小框用不同颜色标记出来frame: 画出的图的框框tick label: 就是坐标的下面的坐标间隔标度1 2 3 4 5 6tick mark: 就是坐标轴的突起的刻度, 如 1 头上对应轴有个突起leg.box: 在标题窗口里面. 可以调整图例的位置, 里面的location设置就可以了axis placement: 用处在 y 轴 坐标在–1.51.5 范围时, 坐标轴若在y=–1.5 处 就去zreo axis点一下就可以了.图像的四个角可以拉大拉小的, 如果不合适自己调.如果有多条曲线, 出现了图例与曲线重叠, 可以试试这样操作：拖动图例的方法, Ctrl+L单击, 就可以使箭头变手形, 拖动图例了. 如还不行, 双击再试试看可以了吗. 也可以在appearance更改.记住set appearance要在窗口里面的select set位置选相应的s0 s1...., main标签下的Legend就是图例, 你可以给每条曲线标记不同名字以便区分.还有需要说一下的画布左边的按钮, 就说一个吧, AS表示恢复（图片设置的外观不变的, 好）. 所以你按钮按错了, 点它就行了.如果用不同坐标刻度值时, 也就是跳跃很大, 如μs, ms, s, 在Axes（双击坐标轴就出来了）的窗口的tick properties框里面选Format, 里面好几种格式, 选Compute(K, M, G....)或Engnieering, 这样就可以一个坐标轴多个单位了. 不用显示那么长的数字串.有时候要一个画布放多张图, 打开菜单Edit | Arrange_graphs对话框, 填写需要几rows几columns, 点应用, 就出来四个小框了（假设2*2格式）. 具体外观, 可以在刚才的框框里继续调节, 点击各个小框点击, 像前面的第2点那样选好文件进来就可以. 具体其他调节参看前面的6步骤操作.也可以使用命令输入, 加上指令就可以了：-pexec &arrange (2, 2, .1, .1, .1, ON, ON, ON)&前面的两个2代表2*2排列, 加上-graph更好, 代表的是后面接的文件放在哪个位置, 从0开始数：xmgrace -graph 2 10_rdf.xvg -graph 1 11.dat -graph 0 rdf01.dat -graph 3 rdf.dat -pexec &arrange(2,2,.1,.1,.1,ON,ON,ON )&说实话, xmgrace做内嵌图像, 我也没懂, 希望会的人能介绍一下. 参考方法：打开菜单Edit | overlay_graphs, 点中在overlay graph的G0, 然后右击, 选creat new, 出现了另一个空数据集. 然后overlay graph和on to, 反正个数据集就可以, 后点同意, 就两个图的, 当然看起来是一个, 自己点一下四周点移动, 就可以看见两个的, 只是重叠而已. 然后一个一个添加数据了. 刚才说了内嵌图形还记得吗, 这里也可以把其中一个图拉小移动作为一个内嵌图, 再修改一下需要显示的坐标范围. 具体靠自己去实践了.有时候我们需要同一图像两条不同y轴的曲线, 如例子：坐标x表示时间, 左边y轴代表压力, 右边y轴代表是温度. 注意调整坐标刻度, 在Axes窗口tick label和tick mark的draw on, 一个选正常normal side, 一个选oppsite side相反, 就可以了. 对边的刻度就不要显示了.选项在data | tramsformation.线性拟合打开图像后点击菜单DATA | TRANSFORMATION | REGRESSION选择SET选择LINEAR FIT按下ACCEPTn出来一个窗口显示拟合的直线方程表达式了. 如果需要斜率, 就自己记下方程. Save the information about the slope and intercept that appear in a blue console as slope.dat.二次方程点击DATA | TRANSFORMATION | INTERPOLATION/SPLINE选择SET | METHOD | CUBIC SPLINE, START设为1, STOP设为10, LENGTH设为500或1000, 然后点击`ACCEPT.非线性拟合输入方程形式, 还要输入参数初始值.如果像gnuplot文件那样修改设置的话, 我们直接修改文本就可以了. 可以参考grace-5.1.25/exemple下的示例文件去学习, 吃透了就是xmgrace的高手了. 下面是其中部分参数的意义.@ title && 标题@ title font 0 标题字体格式@ title size 1.500000 标题字体大小@ xaxis bar linewidth 3.0 坐标轴粗细@ xaxis label “E （KJ/mol）” x坐标轴表示的物理量@ xaxis tick major 100 单位刻度100@ yaxis label && y坐标轴表示的物理量@ yaxis tick major 2 每隔2画一个标度, 也就是精度为1@ s5 line linewidth 3.0 第6条曲线的的粗细@ xaxis tick major 轴上的标度@ xaxis ticklabel char size 2.000000 轴上的标度 显示大小@ s3 hidden false 是否隐藏曲线@ s3 legend “Interface” 第四条曲线的图例的名称标示, 就是表示哪条曲线@ s3 symbol 5 曲线的点的符号, 如圆圈, 星号, 倒三角@ s3 line type 5 曲线变成长虚线, 1是直线@ page size 792, 612 白色画布的大小颜色代码的意义@map color 0 to (255, 255, 255), “white”@map color 1 to (0, 0, 0), “black”@map color 2 to (255, 0, 0), “red”@map color 3 to (0, 255, 0), “green”@map color 4 to (0, 0, 255), “blue”@map color 5 to (255, 255, 0), “yellow”@map color 6 to (188, 143, 143), “brown”@map color 7 to (220, 220, 220), “grey”@map color 8 to (148, 0, 211), “violet”@map color 9 to (0, 255, 255), “cyan”@map color 10 to (255, 0, 255), “magenta”@map color 11 to (255, 165, 0), “orange”@map color 12 to (114, 33, 188), “indigo”@map color 13 to (103, 7, 72), “maroon”@map color 14 to (64, 224, 208), “turquoise”@map color 15 to (0, 139, 0), “green4”上面的这些在.xvg文件开头写出来给我们参考. 让我们明白这些数字代表什么颜色. 还有下面的字体格式部分字体代码的意义：@map font 8 to “Courier”, “Courier”@map font 10 to “Courier-Bold”, “Courier-Bold”@map font 11 to “Courier-BoldOblique”, “Courier-BoldOblique”@map font 9 to “Courier-Oblique”, “Courier-Oblique”@map font 14 to “Courier-Regular”, “Courier-Regular”@map font 15 to “Dingbats-Regular”, “Dingbats-Regular”@map font 4 to “Helvetica”, “Helvetica”@map font 6 to “Helvetica-Bold”, “Helvetica-Bold”@map font 7 to “Helvetica-BoldOblique”, “Helvetica-BoldOblique”@map font 5 to “Helvetica-Oblique”, “Helvetica-Oblique”@map font 20 to “NimbusMonoL-Bold”, “NimbusMonoL-Bold”@map font 21 to “NimbusMonoL-BoldOblique”, “NimbusMonoL-BoldOblique”@map font 22 to “NimbusMonoL-Regular”, “NimbusMonoL-Regular”@map font 23 to “NimbusMonoL-RegularOblique”, “NimbusMonoL-RegularOblique”@map font 24 to “NimbusRomanNo9L-Medium”, “NimbusRomanNo9L-Medium”@map font 25 to “NimbusRomanNo9L-MediumItalic”, “NimbusRomanNo9L-MediumItalic”@map font 26 to “NimbusRomanNo9L-Regular”, “NimbusRomanNo9L-Regular”@map font 27 to “NimbusRomanNo9L-RegularItalic”, “NimbusRomanNo9L-RegularItalic”@map font 28 to “NimbusSansL-Bold”, “NimbusSansL-Bold”@map font 29 to “NimbusSansL-BoldCondensed”, “NimbusSansL-BoldCondensed”@map font 30 to “NimbusSansL-BoldCondensedItalic”, “NimbusSansL-BoldCondensedItalic”@map font 31 to “NimbusSansL-BoldItalic”, “NimbusSansL-BoldItalic”@map font 32 to “NimbusSansL-Regular”, “NimbusSansL-Regular”@map font 33 to “NimbusSansL-RegularCondensed”, “NimbusSansL-RegularCondensed”@map font 34 to “NimbusSansL-RegularCondensedItalic”, “NimbusSansL-RegularCondensedItalic”@map font 35 to “NimbusSansL-RegularItalic”, “NimbusSansL-RegularItalic”@map font 36 to “StandardSymbolsL-Regular”, “StandardSymbolsL-Regular”@map font 12 to “Symbol”, “Symbol”@map font 38 to “Symbol-Regular”, “Symbol-Regular”@map font 2 to “Times-Bold”, “Times-Bold”@map font 3 to “Times-BoldItalic”, “Times-BoldItalic”@map font 1 to “Times-Italic”, “Times-Italic”@map font 0 to “Times-Roman”, “Times-Roman”@map font 43 to “URWBookmanL-DemiBold”, “URWBookmanL-DemiBold”@map font 44 to “URWBookmanL-DemiBoldItalic”, “URWBookmanL-DemiBoldItalic”@map font 45 to “URWBookmanL-Light”, “URWBookmanL-Light”@map font 46 to “URWBookmanL-LightItalic”, “URWBookmanL-LightItalic”@map font 47 to “URWChanceryL-MediumItalic”, “URWChanceryL-MediumItalic”@map font 48 to “URWGothicL-Book”, “URWGothicL-Book”@map font 49 to “URWGothicL-BookOblique”, “URWGothicL-BookOblique”@map font 50 to “URWGothicL-Demi”, “URWGothicL-Demi”@map font 51 to “URWGothicL-DemiOblique”, “URWGothicL-DemiOblique”@map font 52 to “URWPalladioL-Bold”, “URWPalladioL-Bold”@map font 53 to “URWPalladioL-BoldItalic”, “URWPalladioL-BoldItalic”@map font 54 to “URWPalladioL-Italic”, “URWPalladioL-Italic”@map font 55 to “URWPalladioL-Roman”, “URWPalladioL-Roman”@map font 56 to “Utopia-Bold”, “Utopia-Bold”@map font 57 to “Utopia-BoldItalic”, “Utopia-BoldItalic”@map font 58 to “Utopia-Italic”, “Utopia-Italic”@map font 59 to “Utopia-Regular”, “Utopia-Regular”@map font 13 to “ZapfDingbats”, “ZapfDingbats”最后还是说一句, 有些方法不懂话, 看看grace-5.1.25/examples下面的例子, 右击用grace打开, 看看人家的图的设置, 有你想要的效果吗, 仔细琢磨一下人家怎么设置的. 这是最好的学习方法, 比看教程好很多, 当然这是个人认为.◆本文地址: , 转载请注明◆
转载本文请联系原作者获取授权，同时请注明本文来自李继存科学网博客。链接地址：
上一篇：下一篇：
当前推荐数：1
评论 ( 个评论)
扫一扫，分享此博文
作者的精选博文
作者的其他最新博文
热门博文导读
Powered by
Copyright &Date: June 7, 2015
Views: 1,943
浅谈PCA与g_covar+g_anaeig+ddtpd+sigmaplot做自由能面图的方法
文/Sobereva&& 2010-Nov-11
其实这个问题两年前就写过帖子讨论过，但由于近来又有人问，也觉得在作图的时候有些问题值得补充说明一下，故撰此贴。
N个原子的柔性大体系运动轨迹需要3N维笛卡尔坐标描述，然而这样高维度的数据很难直观分析，尤其是随着计算能力越来越高，模拟的体系也越来越大、模拟的时间越来越长、内容越来越复杂，这就需要发展更多的分析方法去挖掘出海量信息中感兴趣的信息，而使数据“噪音”被过滤，比如从复杂轨迹中了解蛋白的折叠/去折叠主要路径、酶重要域的大尺度运动行为、体系构象分布倾向、配体/残基突变/质子化等效应对体系构象的影响等等。一种解决途径是将坐标维度约化得尽可能低，同时尽可能使原始高维空间下的信息最大程度保留，这样就方便人直观考察体系的运动过程。统计学上的Principle component analysis(PCA)正适合解决这个问题，一般是将变量选为3N个笛卡尔坐标（也有其它选择方式，诸如二面角PCA等），要处理的数据集就是各帧结构。首先由MD轨迹构建协方差矩阵，然后计算其本征值和本征向量。这3N个本征向量也就是由原始笛卡尔坐标组合而成的新坐标，它们彼此间从协方差的意义上说是不相关的（也就是把轨迹转换到这些新坐标上后，协方差均为0）。本征向量对应的本征值越大说明体系的运动行为越多地体现在这个本征向量坐标上，若最大的几个本征值的和除以协方差矩阵的迹等于0.85，说明在这几个维度上就可以描述体系85%的运动信息。这样有高本征值的本征向量也称主成分(Principle component, PC)。
实际上一般不将大分子全部N个原子都考虑，比如研究蛋白骨架行为就考虑alpha-C即可，研究某一个区域的运动就选择这个区域的原子即可，不仅计算会省时、省内存，而且如果把过多无关原子也考虑进去，若它们在模拟中运动也比较明显，则感兴趣的原子的运动反倒被掩盖了。另外，做PCA前必须将整个轨迹向指定帧的结构做align以消除平动和转动的影响，前文所谓“运动”皆是指体系的内运动，否则分子内运动行为会被整体运动掩盖住。（实际上即便做了align，整体运动的成分也并没彻底消除，因而出现了基于内坐标的PCA方法，即dihedral PCA）。限于本文篇幅和目的不再对PCA进行更具体讨论，可以参考Principal Components Analysis: A Review of its Application on Molecular Dynamics Data。下载地址：
2 自由能面图(Free energy landscape)
一般所谓自由能面图就是通过一张图来描述大分子各种构象自由能的大小。这就需要知道各种构象{x}出现的概率密度P(x)，然后通过Boltzmann关系计算G(x)：
G(x)=-kT*Ln(P(x))+常数
常数项不用管，它和配分函数有关，是极难计算的。我们关心的只是各种构象的相对自由能而不是绝对自由能，常数项与此无关，所以做自由能面图时常数项可以取任意值。
在平面上最能清晰展现出来的是二维坐标的图，通常将自由能面图的两个坐标轴分别取为PC1、PC2（最大两个本征向量），若选取的原子在模拟中运动范围大，它们一般已经能描述体系70%以上的运动行为。如果PC3的本征值也不小，可以再加一张图以PC1和PC3为轴。（对某些问题也可以用其它描述体系特征的量为轴，比如转动半径）
选定了坐标轴后需要将原始轨迹根据相应的本征向量投影到这两个坐标轴上。如果做散点图，图上就会显现疏密不均的小点，点越密集处说明构象出现概率越大，相应的构象自由能越低。接下来将图上的散点分布转换成概率分布P(x)，再由前式获得G(x)，做它的等值线图或者填色图就得到了自由能面图。接下来将给出实例，通过gromacs的g_covar和g_anaeig子程序、笔者开发的ddtpd以及sigmaplot绘图程序绘制一个10ns模拟中酶蛋白骨架的自由能面图。
3 使用g_covar
g_covar用来生成轨迹的协方差矩阵并计算本征向量和本征值。运行方法（方括号代表可选的项）：
g_covar -f 输入的轨迹文件 -s 输入的结构文件 -n [输入的index文件] -o 输出的本征值文件 -v 输出的本征向量文件 -av 输出的平均结构文件 -l 输出的日志文件 -ascii [输出的协方差矩阵数据文件] -xpm [图形描述的N阶协方差矩阵] -xpma [图形描述的3N阶协方差矩阵]
例如g_covar -s a-sol-prnowat.gro -f allnowat.xtc -o eigenvalues.xvg -v eigenvectors.trr -xpma covapic.xpm
其中allnowat.xtc是一个10000帧的10ns的轨迹，a-sol-prnowat.gro是与它对应的某一帧结构文件。运行后程序会问将哪些范围原子做align，选C-alpha，说明每帧的C-alpha将会朝着输入的结构文件做align；然后会问对哪些范围原子做PCA，选C-alpha（也可以与刚才选的不同）。
由于alpha-C共有232个，所以协方差矩阵是3*232=696阶的方阵，故eigenvalues.xvg含有696个值，本征值顺序是从大到小。由于这个体系处于平衡态，骨架整体运动不明显，所以由eigenvalues.xvg信息绘制的下面的“本征值大小vs本征值序号”的图上看，虽然随序号增大本征值立刻降低，但是收敛不快，超过5号后降低平缓，前十个本征向量总共只能解释体系58.7%的运动信息，PC1和PC2一起只能解释32.2%的信息，说明对此体系全部alpha-C原子做PCA意义不大。由于本文目的主要是说明如何做自由能面图，所以这无所谓。
696个本征向量分别记录在eigenvectors.trr中的每一帧中，其原子坐标代表原3N个笛卡尔坐标是如何组合成本征向量的，这个轨迹虽然也可以用可视化软件观看，但没什么意义。实际上得到的轨迹有698帧，是因为第一帧保存的是align时的参考结构，第二帧是轨迹平均结构。
covapic.xpm是图形化描述的协方差矩阵，在Linux下可以用比如GNOME之眼打开；若在windows下打开，可以先用gromacs自带的xpm2ps转换成eps文件，再用诸如再用诸如photoshop、GhostView查看（也可以用adobe distiller等软件转成pdf）。此文件中3N*3N协方差矩阵已经被转化成N*N协方差矩阵，即坐标轴刻度对应原子序号。越红的区域说明矩阵元数值大，白的区域对应矩阵元数值为0，越蓝代表越负。矩阵对角线就是方差，必大于零，故是白色或偏红色。图中坐标为(28,28)的区域较红，说明26~29号alpha-C在模拟中运动比较明显，实际上这几个残基正处于蛋白质柔性末端，直接考察轨迹也能得出相同结论。非对角元区域可正可负，越正/越负说明相应两个原子运动方式越趋于正相关/负相关。图中虽有深蓝色区域，但从色彩刻度上看其值其实并不很负。alpha-C原子间运动多呈负相关和残基侧链间静电作用导致排斥/吸引而使相关原子产生反相运动有关。若是研究同一残基的原子运动的协方差矩阵，由于某原子运动会拉动周围原子同向运动，所以对角线附近也会偏红。
4 使用g_anaeig
g_anaeig的主要功能之一是将轨迹投影到选取的本征向量上，这里要投影到PC1和PC2上，执行：
g_anaeig -f allnowat.xtc -s a-sol-prnowat.gro -v eigenvectors.trr -first 1 -last 2 -2d 2d.xvg
首先allnowat.xtc各帧会被align到eigenvectors.trr里记录的参考结构上，然后投影到-first和-last所选的本征向量，命令中1、2就是指前两个本征向量，即PC1和PC2。-2d 2d.xvg代表将每帧结构在PC1和PC2上的投影值输出到2d.xvg。
如果装了xmgrace，运行xmgr 2d.xvg，调整显示方式得到下面的散点图：
每个点代表一帧结构，越密集的区域说明这个区域对应的分子构象能量越低。在密集区域中随便取一个点，若这个点的编号是i，就说明轨迹中第i帧的构象是较稳定、能量较低的，反之稀疏区域的点对应的帧是不稳定构象。在xmgrace中获得点的编号首先要将感兴趣的区域放大，然后在设定符号/图例的页面中将显示的符号设为Index，数据点就以其编号显示了。利用PCA可以做簇分析，比如图中左侧一团密集的点对应的构象就可以归为同一个簇，即它们构象相似。下图是投影到PC10和PC11的结果。
比较可见，投影到PC10、11后数据点都凑在一起成球状，看不出任何特征，构象的差异引起的能量差异得不到充分显现，也不可能做簇分析。这就是为什么轨迹要投影到PC1、2的原因，PC1、2最能描述体系运动模式，或者说PC1、2能将体系的自由能面形貌最大程度地展现，使其细节特征能够充分地暴露出来。
使用g_anaeig时如果用“-3d 文件名”，会将轨迹投影到从-first到-last的3个主成分上，输出的是.gro文件，原子坐标的X/Y/Z值分别代表在三个PC上的投影，原子编号就是帧号。如果用“-proj 文件名”，可以同时将轨迹投影到从-first到-last的任意多个主成分上。比如投影到前六个PC，用xmgr对得到的.xvg文件作图得到：
图中表现轨迹在PC1至PC6上投影值随时间的变化。前两个PC波动比较大，说明其数据方差大，这也说明为何在协方差矩阵中其对应的本征值大（本征值=方差）。后面几个PC的波动依次渐缓，尤其是PC5和PC6，完全就像是体系在平衡态的RMSD曲线，体系的运动信息用它们根本描述不出来。
5 使用ddtpd做自由能面图
5.1 ddtpd简介
ddtpd全称Converting dot distribution to probability distribution，是笔者开发的将上述-2d关键词生成的点分布.xvg文件转换成概率或自由能分布的小程序。原理是根据数据最大、最小值设定空间范围，然后根据用户输入的两方向格点数将空间划分为一个个微小的格子，根据 P(x,y)=此网格内数据点数量/总数据点数量/网格面积 来计算不同位置的概率密度，进而转化为自由能。ddtpd v1.3还支持高斯展宽功能，用于解决数据点较少时图像难看的问题。
ddtpd v1.3的下载地址：。
压缩包里的test.xvg就是下文要用的xvg文件。
5.2 作图步骤
启动ddtpd后，依次输入
2d.xvg& //文件名
100& //将X轴划分的格子数
100& //将Y轴划分的格子数
2& //选择输出方式。2代表输出-LnP
y& //令P=0的点也输出。此时数据被输出到当前目录的result.txt下。由于没有数据点分布的区域P=0，没法求对数，对这样的区域ddtpd认为其P是一个很小的常数，自由能直接设为-Ln0.000001。
y& //把-LnP最负的值，此处为-0.821262设为自由能的零点。我们感兴趣的只是相对值，所以可以加减任意常数，这样所有值加上0.821262之后数据最小值就是0，设定色彩刻度会比较方便。此时数据将输出到当前目录的result2.txt下，前两列是PC1和PC2的坐标，第三列是-LnP的值。
sigmaplot比较易用，适合做填色图描绘自由能面。为了让色彩刻度中色彩变化明显的范围覆盖-LnP数据主要分布范围，以使图上数值大小不同的区域能够被清晰地区分开，需要对ddtpd输出的数据做一下调整。注意程序最后在屏幕上输出“Now maximum result is&&& 2.772589”（注：这个最大值是指除P=0以外区域的值）以及“14.636773 means P=0 in this minival area”。打开result2.txt，将P=0的值替换为最大值加上1.0~1.5得到的整数或半整数（凑整数/半整数只是为了方便而已），此例中就是将14.636773替换为4。由于数据较多，不同文本编辑器替换速度由于算法不同差异明显，建议用Ultraedit。
然后将result2.txt放到sigmaplot里做图，得到下图（色彩刻度为默认）：
注意这类图的单位是kT，k是Boltzmann常数(J/K)，T是当前模拟时的温度(K)，所以kT=T*1.^-23J。如果想转化到kcal/mol，就除以4186再乘以阿伏伽德罗常数。比如对于图中的4位置，假设模拟是在310K进行的，就是4*310*1.^-23/*10^23=2.4633kcal/mol.
5.3 格点数目的选择对图像的影响
上图整体效果不错，越蓝区域说明自由能越低，图中显示构象主要分布在左、右上和右下三部分，自由能依次升高。不过图也略有缺点，也就是诸如右上角低能量区域有很多红色或黄色的窟窿，蓝色区域没那么连贯，显得细碎。这与格点数目选取有关，格点数太多不仅图像细碎，做图还慢；格点数太少虽然会连贯，但是缺乏细节，图像边缘时常呈现棱角。一般来说，.xvg文件中数据点越多可以用越多的格点，图像边缘细节会更丰富，且不使图像细碎；当数据点不够多的时候用较大数目格点只有坏处，得到的图甚至就像散点图。下图是此xvg文件50*50格点和200*200格点的图，后者凌乱很不好看，而前者又有些模糊，由于前面100*100做的图也略有零碎，所以75*75格点比较适合做这个xvg文件的填色图。对不同体系应当反复摸索最适合的格点设定，一般不超过60*60至120*120范围。另外，调整色彩刻度上下限也能有效改善图像效果，这里就不谈了，请自行摸索。
5.4 P=0的点对图像的影响
为何要让ddtpd输出P=0的点需要说明一下，尽管这些点看似没有意义。sigmaplot会对输入的数据点的空隙进行差值使图像平滑。如果P=0的区域都不输出，这么大范围的数据空缺区域sigmaplot都会通过插值来弥补。但由于空缺区域往往过大，插值出的结果经常很糟糕。下图是75*75格点，不含P=0的点时sigmaplot给出的结果，和前面的图相对比，会发现在数据点集中分布的区域周围多出了一些的东西（如图的左下角），尤其是在图的中心偏左下的位置多出现一块深蓝区域，如果将它解释为自由能很低的区域就明显错了，这都是因为sigmaplot插值所致。
5.5 对数据点进行高斯展宽
如果用g_anaeig的时候加上-skip 10，就代表每10帧输出一次投影值，前例有10000个点，此时就只有1000个点了，。这样少的数据点很难做出满意的图，哪怕只用50*50格点，图像仍然零碎，尤其是右上角的一片区域，如下图所示。如果进一步降低格点数，图像就不好看了。这时需要利用ddtpd的高斯展宽功能来改善。
所谓高斯展宽，就是把每个点改成用高斯函数来表达（经过归一化），简单来说，比如一个格子内没有数据点，而附近的某个格子内有1个数据点，那么这个格子也能分享到零点几个数据点。这样就不会造成数据点分布密集区域中由于恰好有一些格子没有数据点（格点密度高的时候这种情况尤甚），而使图像出现一堆窟窿或者显得零碎。下图是经过高斯展宽后的图，明显数值连续性比上图好很多。为了效果更好，图中色彩刻度设为了-1~4.5。
此图做法是启动ddtpd，依次输入：
2dskip10.xvg& //文件名
50& //将X轴划分的格子数
50& //将Y轴划分的格子数
4& //输出将数据点经过高斯展宽后的-LnP
1& //展宽系数，默认是1，其值越大展宽效果越明显，但会造成细节区域越难以显现，可反复尝试。假设输入的是k，就代表展开成的高斯函数的半高宽是k*sqrt(dx^2+dy^2)，其中dx和dy分别是X/Y方向的网格宽度。
y& //输出P=0的点
y& //调整G的零点
注意屏幕上提示数据的最大值为5.11，20.305551对应P=0区域的值。对于利用高斯展宽得到的结果，作图前将P=0区域的值替换成最大值就行了（不用再加1.0~1.5），此例即把result2.txt里的20.305551都替换为5.11，然后放到sigmaplot里作图，调调色彩刻度即得到上图。}