layout: post title: 使用GROMACS进行团簇分析 categories:
GROMACS自带了一个团簇分析工具cluster, 但这个工具主要用于对蛋白的构象进行分类, 支持的通用距离文件为xpm格式, 基本没法使用由其他程序生成的距离矩阵, 除非修改源码.
在合金材料的模拟中, 有时需要对合金进行团簇分析. 自己写代码有点麻烦, 我们可以组合GROMACS自带的两个工具, mdmat和cluster来实现这个功能.
具体来说, mdmat可以获得轨迹中所有残基对之间的距离, 并将其写为xpm格式. 根据得到的xpm距离矩阵, 就可以调用cluster进行团簇分析. 显然, 这种方法需要先使用mdmat输出xpm文件, 再由cluster读入xpm文件. 当体系中原子很多时, xpm文件会变得相当大, 写出和读入都要耗费不少时间. 因此, 这种方法虽然不需要改动任何代码, 实施起来简单, 但处理效率就差强人意了. 如果你想要提高处理速度, 可以修改GROMACS的源代码, 将mdmat和cluster两个函数进行整合, 让距离矩阵的计算和分析都在程序内部完成. 根据测试, 对10k原子的体系, 在我的Intel I5-3550/8G台式机上, 不改代码的情况下, 处理一帧需要40秒左右, 修改代码后大约需要10秒, 提速还可以. 如果再对代码进行下优化, 应该能达到3秒一帧的速度, 基本满足需要了.
简单示例一下命令吧. 假设有了合金模拟的轨迹, 先将其转换为GROMACS格式的轨迹conf.gro, 然后使用make_ndx创建index.ndx文件, 再写个简单的topol.top文件和grompp.mdp文件. 这样输入文件就准备好了. 接下来执行下面的命令就可以得到团簇分析结果了.
如果你有很多帧轨迹, 那就写个脚本来连续调用后面两个步骤就好了. 当然, 在进行团簇分析的时候还有很多选项可用的, 具体的可以查阅GROMACS的文档.
有了团簇的分析结果, 可以写个脚本将每个团簇提取出来, 进行更详细的分析. 这些就不再细说了.