- 2017-06-22 21:54:18
使用GROMACS计算径向分布函数和配位数比较简单, 只要了解一下gmx rdf程序的使用即可. 但即便是简单的工作, 有时候也需要我们稍微深入思考一下, 如何才能做得更快更好, 如何才能尽量减少重复运行所需的时间. 只有这样, 我们才能不断突破自己的心理舒适区, 提高自己. 否则的话, 每次都暴虎冯河, 自己既做得累, 也易生厌. 毕竟, 人, 是要做些有技巧性, 艺术性的东西的.
就以我最近要处理的问题来说明吧. 我跑了一个纤维素二聚体的模拟, 这个分子共有45个原子,
跑完后, 要计算二聚体中每个原子与溶剂水中的氧原子OW, 氢原子HW的RDF, 从而确定每个原子的第一溶剂化层厚度.
要算RDF么, 首先要做出索引文件, 把二聚体的每个原子, OW, HW放到不同的组中. 简单. 执行
gmx make_ndx -f conf.gro
再使用下面的指令
a OW: 选择水中的氧原子
a HW: 选择水中的氢原子
splitat 1: 将组1拆分, 每个原子一个组. 组1是other也就是二聚体.
q: 保存索引文件, 退出程序.
这样我们就得到了索引文件, 可用于指定要计算的RDF了.
打开得到的索引文件看看, 各个组的排列顺序如下
[ System ]
[ Other ]
[ ROH ]
[ 4GB ]
[ 0GB ]
[ Water ]
[ SOL ]
[ non-Water ]
[ OW ]
[ HW ]
[ Other_HO1_1 ]
[ Other_O1_2 ]
[ Other_C1_3 ]
...
[ Other_H2_43 ]
[ Other_O2_44 ]
[ Other_H2O_45 ]
使用这个索引文件指定要计算的RDF不方便, 因为二聚体的原子编号和索引组编号不一致, 我们把它们调整到一致(注意, 索引组编号是从0开始的), 调整后, 索引文件顺序如下
[ System ]
[ Other_HO1_1 ]
[ Other_O1_2 ]
[ Other_C1_3 ]
...
[ Other_H2_43 ]
[ Other_O2_44 ]
[ Other_H2O_45 ]
[ OW ]
[ HW ]
[ Other ]
[ ROH ]
[ 4GB ]
[ 0GB ]
[ Water ]
[ SOL ]
[ non-Water ]
这样在计算RDF时, 只要指定1 46就可以算出二聚体的1号原子与OW的RDF了. 类似操作45次, 就能得到45个原子与OW的RDF了. 再类似操作45次, 得到45个与HW的RDF. 完成.
如果你每次都这么做上90次, 还要保证不输错数字, 那还是很能锻炼耐心的. 我却没有这么大的耐心, 所以就只能偷个懒了, 使用bash脚本
for((i=1; i<=45; i++)); do
echo -e $i'\n46' | gmx rdf -f -n
done
当然, 如果你bash版本支持, 你也可以使用
for i in {1..45}; do
echo -e $i'\n46' | gmx rdf -f -n
done
如果你还知道parallel, 那就更好了, 你可以并行执行RDF计算, 这会大大节省时间
paralle echo -e {}'\\n46' | gmx rdf -f -n ::: {1..45}
所以, 学习一点简单的bash语句, 你的生活可以更轻松. 学习的过程并不舒服, 可投入一段时间的不舒服, 换取以后更多时间的舒服, 你要如何选择呢?
补充
上面的示例虽然用了gmx 5.x版本, 但做法还是基于4.x版本的. gmx 5.x版本的rdf实际支持一次计算多个组的RDF, 只要连续给出组编号就可以了. 这样得到的输出文件中, 含有多个RDF.
