layout: post title: 使用GROMACS计算径向分布函数 categories:
使用GROMACS计算径向分布函数和配位数比较简单, 只要了解一下gmx rdf程序的使用即可. 但即便是简单的工作, 有时候也需要我们稍微深入思考一下, 如何才能做得更快更好, 如何才能尽量减少重复运行所需的时间. 只有这样, 我们才能不断突破自己的心理舒适区, 提高自己. 否则的话, 每次都暴虎冯河, 自己既做得累, 也易生厌. 毕竟, 人, 是要做些有技巧性, 艺术性的东西的.
就以我最近要处理的问题来说明吧. 我跑了一个纤维素二聚体的模拟, 这个分子共有45个原子,
<figure><script>var Mol1=new ChemDoodle.TransformCanvas3D('Mol-1',642,396);Mol1.specs.shapes_color='#fff';Mol1.specs.backgroundColor='black';Mol1.specs.set3DRepresentation('Ball and Stick');Mol1.specs.projectionPerspective_3D=false;Mol1.specs.compass_display=true; /*//Mol1.specs.atoms_resolution_3D=15; //Mol1.specs.bonds_resolution_3D=15; //Mol1.specs.crystals_unitCellLineWidth=1.5;*/ Mol1.nextFrame=function(delta){var matrix=[];ChemDoodle.lib.mat4.identity(matrix);var change=delta*Math.PI/15000;ChemDoodle.lib.mat4.rotate(matrix,change,[1,0,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,1,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,0,1]);ChemDoodle.lib.mat4.multiply(this.rotationMatrix, matrix)}; Mol1.startAnimation=ChemDoodle._AnimatorCanvas.prototype.startAnimation;Mol1.stopAnimation=ChemDoodle._AnimatorCanvas.prototype.stopAnimation;Mol1.isRunning=ChemDoodle._AnimatorCanvas.prototype.isRunning;Mol1.dblclick=ChemDoodle.RotatorCanvas.prototype.dblclick;Mol1.timeout=5;Mol1.handle=null; var Fmol='45\ncellulose\nH 15.370 11.250 9.500\nO 15.240 11.230 10.460\nC 14.510 12.410 10.800\nH 14.350 13.020 9.920\nO 15.360 13.260 11.670\nC 14.710 14.550 12.130\nH 14.490 15.120 11.220\nC 15.760 15.220 12.970\nH 15.350 16.280 13.230\nH 15.900 14.680 13.850\nO 17.060 15.280 12.410\nH 17.490 15.940 12.980\nC 13.470 14.270 12.950\nH 13.770 13.750 13.840\nC 12.480 13.310 12.160\nH 12.160 13.990 11.410\nO 11.300 12.920 12.890\nH 10.880 12.190 12.470\nC 13.240 12.020 11.560\nH 13.530 11.330 12.350\nO 12.190 11.510 10.650\nH 12.400 10.580 10.430\nO 12.870 15.540 13.250\nC 12.800 16.080 14.610\nH 13.770 15.880 15.040\nO 11.760 15.260 15.250\nC 11.330 15.660 16.560\nH 12.080 15.500 17.340\nC 10.140 14.750 17.020\nH 9.690 14.870 18.000\nH 9.400 14.800 16.220\nO 10.700 13.430 17.100\nH 11.230 13.220 16.330\nC 10.910 17.140 16.610\nH 10.110 17.290 15.910\nO 10.410 17.480 18.040\nH 10.200 18.460 18.030\nC 12.080 18.100 16.120\nH 12.960 18.060 16.720\nO 11.580 19.430 16.050\nH 12.240 20.030 16.210\nC 12.480 17.620 14.660\nH 11.540 17.700 14.020\nO 13.510 18.380 14.100\nH 13.430 18.550 13.120\n'; Mol1.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol1.startAnimation();Mol1.stopAnimation();function setProj1(yesPers){Mol1.specs.projectionPerspective_3D=yesPers;Mol1.setupScene();Mol1.repaint()}function setModel1(model){Mol1.specs.set3DRepresentation(model);Mol1.setupScene();Mol1.repaint()}function setSpeed1(){Mol1.timeout=500-document.getElementById('spd1').value;Mol1.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol1.startAnimation()}</script><br><span class='meta'>视图: <input type='radio' name='group2' onclick='setProj1(true)'>投影 <input type='radio' name='group2' onclick='setProj1(false)' checked=''>正交 速度: <input type='range' id='spd1' min='1' max='500' onchange='setSpeed1()'/><br>模型: <input type='radio' name='model' onclick='setModel1('Ball and Stick')' checked=''>球棍 <input type='radio' name='model' onclick='setModel1('van der Waals Spheres')'>范德华球 <input type='radio' name='model' onclick='setModel1('Stick')'>棍状 <input type='radio' name='model' onclick='setModel1('Wireframe')'>线框 <input type='radio' name='model' onclick='setModel1('Line')'>线型 <input type='checkbox' onclick='Mol1.specs.atoms_displayLabels_3D=this.checked;Mol1.repaint()'>名称<br>左键: 转动 滚轮: 缩放 双击: 自动旋转开关 Alt+左键: 移动</span><br><figurecaption>Fig.1</figurecaption></figure>跑完后, 要计算二聚体中每个原子与溶剂水中的氧原子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脚本
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #AA22FF; font-weight: bold">for</span><span style="color: #666666">((</span><span style="color: #B8860B">i</span><span style="color: #666666">=</span>1; i<<span style="color: #666666">=</span>45; i++<span style="color: #666666">))</span>; <span style="color: #AA22FF; font-weight: bold">do</span> <span style="color: #AA22FF">echo</span> -e <span style="color: #B8860B">$i</span><span style="color: #BB4444">'</span><span style="color: #BB6622; font-weight: bold">\n</span>46<span style="color: #BB4444">'</span> | gmx rdf -f -n <span style="color: #AA22FF; font-weight: bold">done</span> </pre></div>当然, 如果你bash版本支持, 你也可以使用
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #AA22FF; font-weight: bold">for</span> i in <span style="color: #666666">{</span>1..45<span style="color: #666666">}</span>; <span style="color: #AA22FF; font-weight: bold">do</span> <span style="color: #AA22FF">echo</span> -e <span style="color: #B8860B">$i</span><span style="color: #BB4444">'</span><span style="color: #BB6622; font-weight: bold">\n</span>46<span style="color: #BB4444">'</span> | gmx rdf -f -n <span style="color: #AA22FF; font-weight: bold">done</span> </pre></div>如果你还知道parallel, 那就更好了, 你可以并行执行RDF计算, 这会大大节省时间
所以, 学习一点简单的bash语句, 你的生活可以更轻松. 学习的过程并不舒服, 可投入一段时间的不舒服, 换取以后更多时间的舒服, 你要如何选择呢?
上面的示例虽然用了gmx 5.x版本, 但做法还是基于4.x版本的. gmx 5.x版本的rdf实际支持一次计算多个组的RDF, 只要连续给出组编号就可以了. 这样得到的输出文件中, 含有多个RDF.