计算之道的博客

在完成分子动力学模拟之后，需要对模拟轨迹进行分析，计算径向分布函数（Radial Distribution Function, RDF）即是一种常见的分析。大多分子动力学模拟软件提供了计算RDF的功能，但往往不够灵活，也不能适应不同的轨迹格式。是一个用于处理分子动力学轨迹的Python包，它能够接受多种轨迹格式，将不同格式的轨迹作为统一的轨迹对象处理；同时，它也提供了许多性质计算的功能，其中就包括RDF；最后，通过使用Python代码，使得性质计算灵活了许多。

cpptraj输出的第一行是标题行。从第二行开始，每一行包含每个氨基酸对应的二级结构，包含上述8种二级结构，其中，0指的是Coil。#echo $tCount $len_I #len_I为蛋白质氨基酸总数，需要提前获取。#vmd识别TEBHGIC7种，这里将loop区域定为Turn(T)和Coil(C)轨迹的话可以用cpptraj来计算（内置dssp），单纯PDB的话可以直接用VMD来计算，

在对M86-S1体系进行显式溶剂模拟时的原子数为67644，而隐式溶剂模拟时的原子数为6296。原本以为隐式溶剂MD模拟要比显式的快得多，但是查看速度发现前者的速度为261 ns/day，而后者的速度为135 ns/day。总之就是，对于小体系的话，隐式比显式快，而对于大体系的话，显式模型可能更胜一筹。原因是，速度快慢取决于非键相互作用的计算，显式溶剂模拟有非键截断距离，而隐式溶剂模拟需要计算所有原子的非键相互作用。附：prod.in文件。

命令：

越大的体系，一开始水盒子得加的越厚，因为有空隙会收缩，需要保证收缩后仍大于cutoff值+2A，否则一个周期内的蛋白会和另一个周期内的自己作用，而发生形变。所以水盒子要厚，一般我默认盒子12A起，density也要多步，至少4次5000步，最后看密度波动，波动值小于1%了就可以走下一步了。NPT时，密度平衡会使得盒子的周期性边界收缩。NVT时水的聚集导致形成“真空气泡”。

参考：http://jerkwin.github.io/2015/12/08/%E4%BD%BF%E7%94%A8AmberTools+ACPYPE+Gaussian%E5%88%9B%E5%BB%BA%E5%B0%8F%E5%88%86%E5%AD%90GAFF%E5%8A%9B%E5%9C%BA%E7%9A%84%E6%8B%93%E6%89%91%E6%96%87%E4%BB%B6/Gaussian 09C.01及后续版本恢复了误删的代码并且加上了gesp的代码，所以以上关键字全部可以使用。

（Xi(tj)-xi）就是t时刻某个原子的位置减去初始时刻它的位置（也是位置偏移量），然后取所有时刻（可以理解为50w个2fs）的偏移量的平方和，然后对时间T取平均，然后开方，就是这个原子在时间T内，相对于初始时刻的RMSF。δi就是某一帧的第i个原子的位置减去参考构象中它的位置（位置偏移量），然后取所有原子的偏移量的平方和，然后对原子数N取平均，然后开方，就是这一帧结构相对于参考构象的RMSD。在轨迹分析中，最经常用，最简单，也最有用的就是这两巨头，二者都是对位移的平方和再求平方根，最后求得均值。

整体名就是complex，是自己定义的复合物的名字；残基编号在此复合物中必须是独一无二的（且从第一个残基开始依次编号，如第一个残基的编号为2，则即使第2个残基的编号为888，导入tleap之后默认残基编号也是3）；原子名就是要bond的原子的名字，如CA。bbb = sequence{WAT} # usage: sequence ，LIST中可以用逗号或者用空格。bond 整体名.残基编号.原子名 整体名.残基编号.原子名。1，使用tleap创建一个短肽。

解决：猜想会不会是openmpi环境的问题（搞MSDFT的时候将1.6.3版本的openmpi换成了1.4.4），于是将openmpi的版本由1.4.4换成了1.6.3，成功运行。很久没用过amber跑动力学了，在做6ef8动力学的时候想尝试一次，但是没有办法用sander.mpi，有报错。

（2）默认使用最新版本，目前的版本是2.2.0（使用的CgenFF版本是4.0）,可以选择使用1.0版本（使用的CgenFF版本是3.0.1），输出结果与charmm-gui结果是一样的。（1）使用方法是先用gv保存分子的mol2文件，然后将mol2文件中的所有Ar替换为ar，在Molecule Name一行改为小分子的名字，然后导入分子，就可以直接生成。2、下载与生成小分子相匹配的CGenFF力场，如：使用的CgenFF版本是2.2.0，则需要使用CgenFF的4.0版本。1、下载最新版本的蛋白质力场。

用了rigidbonds，就把与氢原子相连的化学键的键长固定在平衡值附近了，因为timestep的选取是和体系的振动频率有关，而氢的振动频率最快，因此，加上rigidbond就降低了振动频率，所以才能用2fs，如果不用rigidbond，则2fs的步长太大，用0.5fs比较合适。——>>>修改：实战中发现体系平行一段时间之后，会有部分小分子远离原来的位置，所>以在体系升温之后，就应该做bond。####上面的是设置系统压力的类型，下面的是设置控压方法，压力应该有，控压方法不一定有（只有平衡时才用控压）。

gmx_mpi genion -s ions.tpr -o pro2_solv_ions.gro -p topol.top -pname NA -nname CL -nn 0 #没添加。3、定义单位盒子并填充溶剂。

>show representation #其中representation可以为：cartoon, ribbon, dots, spheres, surface和mesh。Pymol> log_open script-file-name.pml #记录一个文本文档，该文件的后缀名应为.pml。一、Snapshot（用不同颜色显示蛋白不同二级结构，用不同颜色显示不同小分子，显示连键）Pymol> png file-name # 图片被保存在PyMOL安装默认的文件夹中。

十、可用两个库文件调用力场参数，然后合并（如果主要想用parm数据库，或者把parm数据库中的need revision删除，然后先调用gaff，再调用parm。这样可以防止因为不可抗力导致的程序停止之后，仍能导出与能识别的最后一帧匹配的，带速度的rst。#对于没有原子跑出盒子边界的体系。注意：整体名就是complex，是自己定义的复合物的名字，原子序号在此复合物中必须是独一无二的，原子名就是要bond的原子的名字，如CA。六、将跑出水盒子的镜像回来，用到cpptraj的autoimage命令。

TK：键角力常数（The harmonic force constants for the angle "ITT"-"JTT"- "KTT" in units of kcal/mol/(rad**2) (radians are the traditional unit for angle parameters in force fields).）ITT-JTT-KTT：成键角原子 （The atom symbols for the atoms making an angle.）PN：二面角势垒的周期性。

Windows下安装：下载安装文件压缩包，解压。jmol的运行需要系统安装了java运行环境（JRE）。在此条件下，win用jmol.bat，linux用jmol.sh, mac用jmol.mac。shell中运行可以直接在jmol.sh后加结构文件名，win环境中可直接将结构文件往jmol界面拖即可加载该结构文件。其中，0.02是轨道等值面的数值。30是轨道序数，也可以用homo和lumo关键词。translucent是半透明选项。介绍：rcsb数据库上也是直接应用的jmol。