gromacs读入碳纳米管方法.docVIP

与gromacs读入碳纳米管的方法 有人问我leap怎么读入碳纳米管，这里说一下，顺便也说一下gromacs读入的方法。只说操作方法，不涉及选用力场等问题。 碳纳米管结构比较大，而结构很有规律，一般仅有C原子。若不考虑边界问题，可以将全部电荷当成0，这种近似并不会对结果真实性有多少影响。对于碳纳米管这样的特征，适合当成一般的小分子来处理，但不宜用antechamber而且也没必要用。 先用Nanotube modeler生成一段碳纳米管，含240个C原子，导出结构为tube.pdb。但是会发现，原子名都是相同的，而leap要求结构中相同残基中每个原子都有独立的名字，所以需要先改名。比较方便的方法是用ultraedit，打开后开启列模式，选定从第1个原子到最后一个原子的15和16列(即原子名C后面的两列)，选Column-Insert Number，OK，即把原子名改为了C1至C240，保存为tubename.pdb。 这里假定C原子使用gaff力场的c2原子类型对应的力场参数，依次运行xleap -f leaprc.gaffa=loadpdb /sob/tubename.pdbbondbydistance a 默认是2埃内成键，故碳纳米管的碳原子都会正常成键，可用edit a 检查。然后把碳纳米管内全部240个原子都设成c2原子类型，可以用批处理脚本。由于leap只能一条一条运行批处理文件中的内容而不支持shell脚本的循环，需要将循环转化成普通命令脚本再在leap运行。写脚本：for ((i=1;i=240;i=i+1))doecho set a.1.$i type c2 leapdodone在shell下运行，得到leapdo，然后在leap里运行source leapdo即可check a 进行参数检查，应该OK。如果想改力场参数或缺少某些参数，可以edit gaff往里添加和修改。如果要加溶剂和普通情况无异，这里略过。最后saveamberparm a tube.top tube.inpcrd =============gromacs读入碳纳米管还是用上面的tube.pdb为例子。小分子在gromacs中常用prodrg来处理，但是理由和antechamber一样并不适合，而且有更方便的办法。gromacs提供了一个构建拓扑文件的工具x2top，专适用于构建结构规律性很强的体系。x2top中有一些bug和“规矩”，而且在不同gromacs版本中bug和“规矩”还不一样，这里使用gromacs4.0.4的x2top。 首先需要写n2t文件，比如使用ffgmx2力场，就在力场文件所在文件夹里写一个ffgmx2.n2t，内容如下：C??? CX??? 0.0??? 12.011 3??? C 0.14?? C 0.14 C 0.14C??? CX??? 0.0??? 12.011 2??? C 0.14?? C 0.14C??? CX??? 0.0??? 12.011 1??? C 0.14第一行说明如果体系中任何一个C原子，与周围3个C原子的距离都在0.14nm左右（判断标准大概为正负10%），就把它当作CX原子类型（原子类型名字自定，可以不属于力场包含的原子类型），电荷为0.0，相对原子质量为12.011。第二行和第三行与第一行意义类似，即代表与两个C原子和与一个C原子成键的C原子被当成什么类型、多少电荷、多少原子质量，此例中将碳纳米管所有原子都当成CX。 运行x2top -f tube.pdb -o tube.top -ff gmx2这样就得到tube.top，里面包含这个碳纳米管全部键结项。-ff代表用什么力场，此处即ffgmx2，如果-ff select则是出现列表自行选择。这里如果使用ffgmx力场即便正确写了ffgmx.n2t也可能无法正确执行，即不能正确根据距离判断成键，此时应换用别的力场来执行x2top。如果运行后程序卡住不动，可尝试加-nopbc解决，这是一个bug。默认情况下，x2top会自动在.top里相应键结项加上力场参数，这种方式自动加入的力场参数并不是根据力场中相应原子类型间的键结参数加入的。平衡距离/键角/二面角就是tube.pdb中的相应项的距离/键角/二面角，键的力常数/键角力常数/二面角力常数默认是400000/400/5，可以分别用-kb、-kt、-kp设定。如果不想让其自动加入，可以加上-noparam。 如果n2t中把C原子转换为ffgmx2力场中已有的类型，比如CB，那么此时这个.top搭配tube.pdb已经可以直接用于模拟了。此例将C设为力场中没有的原子类型CX是为了假设需要加入专门的力场参数的情况，此时需要定义新的原子类型