计算化学公社
标题: 使用Gromacs模拟碳纳米管的一个简单例子 [打印本页]
作者Author: sobereva 时间: 2014-12-14 22:54
标题: 使用Gromacs模拟碳纳米管的一个简单例子
注：本文的做法如今已经完全过时，如今强烈建议用Sobtop（http://sobereva.com/soft/Sobtop）建立纳米球/管/板的拓扑文件，比x2top灵活得多，参考网页里的例子。另外，北京科音分子动力学与GROMACS培训班（http://www.keinsci.com/workshop/KGMX_content.html）里有模拟纳米管、纳米板的完整讲解和详细的例子。
使用Gromacs模拟碳纳米管的一个简单例子
A simple example of using Gromacs to simulate carbon nanotubes
文/Sobereva @北京科音 2014-Dec-15
很早以前写过一个帖子《amber与gromacs读入碳纳米管的方法》（http://sobereva.com/39）介绍了怎么用gromacs处理碳纳米管，不过貌似还是有很多人不会。这里就提供一个在真空中和在水中模拟碳纳米管的简单且完整的例子。这个小教程本来是给一个国际友人写的，所以都是英文。本文用的是Gromacs 4.6.5，文中的CNT10指的是一个很普通的碳纳米管，cnt10.pdb是其结构文件，这可以用Nanotube Modeler很容易地产生。
文中涉及到的文件在此：
(, 下载次数 Times of downloads: 404)
(, 下载次数 Times of downloads: 328)
(, 下载次数 Times of downloads: 413)
(, 下载次数 Times of downloads: 372)
Create a new file named atomname2type.n2t in "share/gromacs/top/gromos54a7.ff" folder in gromacs directory, and then copy below content into it
C CR1 0.0 12.011 3 C 0.142 C 0.142 C 0.142
C CR1 0.0 12.011 2 C 0.142 C 0.142
C CR1 0.0 12.011 1 C 0.142
That means when we use x2top to generate topology file under G54A7 forcefield, if a carbon atom has one or two or three neighbours with a distance about 0.142nm (normal C-C bond length in CNT), then this carbon will be recognized as CR1 atom, which is the atomtype corresponding to the carbon in aromatic CH-group of G54A7. This atomtype is suitable for representing vdW interaction between CNT carbons and environmental atoms.
Run below command to generate the CNT10.top:
g_x2top -f cnt10.pdb -o CNT10.top -ff select -nopbc -name CNT10 -kb 400000 -kt 600 -kp 150
Select G54A7 from the forcefield list, then CNT10.top will be yielded in current folder, and the molecular name is CNT10 (specified by "-name"). -nopbc have to be specified, otherwise x2top can't work normally. -kb, -kt and -kp are used to define the force constant of bond, angle and dihedral terms.
In the CNT10.top, the [ bonds ] field looks like below
1 3 1 1.420000e-01 4.000000e+05 1.420000e-01 4.000000e+05
2 3 1 1.420000e-01 4.000000e+05 1.420000e-01 4.000000e+05
2 4 1 1.420000e-01 4.000000e+05 1.420000e-01 4.000000e+05
2 27 1 1.420000e-01 4.000000e+05 1.420000e-01 4.000000e+05
3 218 1 1.420000e-01 4.000000e+05 1.420000e-01 4.000000e+05
...
The 4th column is the equilibrium length determined based on the input geometry (cnt10.pdb), the 5th column is the force constant set by -kb. The last two columns are redundant, you can simply ignore or even directly delete them. The content of [ angles ] and [ dihedrals ] fields are similar to [ bonds ].
The bond, angle and dihedral forcefield parameters we set above are not important, they are mainly used to keep the nearly rigid structure of CNT during simulation, so the force constant can be simply set to a large value; however, too large values will cause high-frequency vibrations and make the simulation unstable. If you want to make the CNT more flexible (rigid), you can decrease (increase) the dihedral force contant. Currently, AFAIK, no forcefield contains bond, angle and dihedral parameters specifically optimized for CNT modelling.
Then there are two cases, you can follow either one
1 MD simulation in vaccum
grompp -f md_vacuum.mdp -c cnt10.pdb -p CNT10.top -o md_vacuum.tpr
mdrun -v -deffnm md_vacuum
2 MD simulation in solvated box
editconf -bt triclinic -f cnt10.pdb -o cnt10_box.gro -d 2
genbox -cp cnt10_box.gro -cs spc216.gro -o cnt10_wat.gro -p CNT10.top
Add below sentence at the head of CNT10.top
#include "gromos54a7.ff/spce.itp"
Start simulation:
grompp -f md.mdp -c cnt10_wat.gro -p CNT10.top -o md.tpr
mdrun -v -deffnm md
模拟和石墨烯、碳球的过程和此文没有任何区别，把文中的cnt10.pdb替换为相应的结构文件即可。石墨烯结构也可以由Nanotube Modeler产生。Nanotube Modeler自带的富勒烯库里有大量碳球结构，各种原子数的各种富勒烯异构体也可以用CaGe产生，CaGe的基本使用方法见《生成富勒烯的螺旋算法简介以及使CaGe中的编号与Fowler-Manolopoulos编号相符的方法》（http://sobereva.com/104）
作者Author: ruanyang 时间: 2014-12-15 07:33
强烈顶贴！Sob老师就是强！
作者Author: shanshan 时间: 2015-3-5 09:18
Sob老师，问一个比较傻的问题：这里如果用
g_x2top -f cnt10.pdb -o CNT10.top -ff select -nopbc -name CNT10 -kb 400000 -kt 600 -kp 150
atomname2type.n2t可以让距离0.142nm的碳看成是CR1的碳， 但是怎么得到键，角的信息呢？
因为我在用的时候老是显示： can not find forcefield for atom C with 0 bonds
非常感谢！！！
作者Author: sobereva 时间: 2015-3-5 18:50
你先严格、完整地重复一下帖子中的例子，需要的文件都提供了。如果不行，换成文中的gmx版本（有的版本在判断时可能有bug）
只要能根据n2t正确判断出原子，在top里bond、angle、dihedral项就都会有了
作者Author: shanshan 时间: 2015-3-6 03:24
非常感谢~~~试了无数遍终于发现是版本问题。
作者Author: lingogo 时间: 2015-4-9 10:04
你好，我想请问一下，如果初始的pdb是一个蛋白质和纳米管的复合物，然后再加溶剂进行动力学，那么生成top的过程是怎么样的？
作者Author: sobereva 时间: 2015-4-9 13:57
单独构造一个纳米管的itp文件。然后用一般方法构建蛋白质的top文件，把纳米管的include进去并把纳米管的结构和蛋白质结构合并。然后加水。
作者Author: xiaowu759 时间: 2015-4-9 18:00
谢谢分享！
作者Author: lingogo 时间: 2015-4-14 10:34
sob老师,再请问一下,碳纳米管的itp文件如何构造呢?
作者Author: sobereva 时间: 2015-4-18 05:16
按此贴的方法用x2top就得到了
作者Author: lingogo 时间: 2015-4-18 10:14
sob老师，不好意思我比较笨，我用了这个方法，最后只得到了CNT10.top，itp文件需要用什么命令生成或者转化吗？
作者Author: ruanyang 时间: 2015-4-18 14:25
你把CNT10.top改名为CNT10.itp,然后自己写一个top文件将这个itp文件include进去就可以了！建议看看手册第5章的内容
作者Author: lingogo 时间: 2015-4-19 08:55
非常感谢
作者Author: qiaobeiming 时间: 2017-10-29 20:35
看来这个软件用途广，学习
作者Author: wh.shen 时间: 2019-6-14 23:16
本帖最后由 wh.shen 于 2019-6-14 23:38 编辑
sob老师，按照您的方法在gromos54a7.ff力场下可以做碳纳米管的模拟。我要做的是在amber03力场下丝蛋白加入纳米管，这个时候我在主top文件中include CNT.itp（CNT.itp是在gromos54a7.ff力场下生成的），但是在生成运行文件时出现截图中错误。
当我尝试在amber03力场下生成CNT.itp时，（同样将CNT.itpinclude到主top中）错误变为找不到CR1类型原子。
请教老师！@sobereva
作者Author: sobereva 时间: 2019-6-14 23:56
top文件里字段顺序不对。forcefield.itp必须是第一个引入的
作者Author: wh.shen 时间: 2019-6-15 10:02
sob老师，我尝试了改字段位置还是没有搞定，麻烦您帮我看下，fws.top是主top，CNT.itp是纳米管的top@sobereva
结构文件传不上来。。。结构文件https://gitee.com/silkcnt/silkcnt/blob/master/silkcnt.gro我放到链接里了
作者Author: sobereva 时间: 2019-6-16 00:49
如置顶的新社员必读贴和论坛首页的公告栏所示，上传超过500KB的文本型文件上传前必须先压缩再上传，以节约论坛空间、加快下载速度。已经把你的违规超标的附件删了，以后注意。
你的fws.top里根本就没引入CNT.itp。如果把后者引入前者，显然得把后者对forcefield.itp的引用删掉，否则[defaults]明显会出现两次
作者Author: lao7 时间: 2019-7-16 23:26
请问，我执行这一步的时候：“g_x2top -f cnt10.pdb -o CNT10.top -ff select -nopbc -name CNT10 -kb 400000 -kt 600 -kp 150”为什么显示“bash:g_x2top: command not found...” gromacs版本是2018.4. 上面的atomname2type.n2t文件夹已经建立。执行GMX命令显示软件安装正常。我是把sob提供的四个附件文件拷贝到home/user下，和gromacs安装文件放在一个文件夹下。
好些年学的，忘完了！谢谢！
作者Author: sobereva 时间: 2019-7-17 02:17
这个教程是对于5.0版以前的，不适合当前版本
作者Author: 杨洁 时间: 2019-9-2 16:29
本帖最后由 杨洁 于 2019-9-2 16:31 编辑
sob老师，您好,刚接触GMX。按您的方法进行操作（我的版本是2019.3），在md_vacuum里添加了cutoff-scheme = group,改到和4.6的默认一致
grompp -f md_vacuum.mdp -c cnt10.pdb -p CNT10.top -o md_vacuum.tpr
没有得到tar结果，请问这是什么原因呢?
作者Author: sobereva 时间: 2019-9-3 00:29
我不知道你说的tar是什么
加上-maxwarn 10，忽略最多10个warning就完了
作者Author: 杨洁 时间: 2019-9-3 11:38
谢谢，好了！不好意思,tpr打错了。
作者Author: 杨洁 时间: 2019-9-3 17:09
sob老师，还有一个问题。如果我想模拟管壁带电的碳纳米管，需要对力场的选择和n2t文件进行改动吗？
作者Author: sobereva 时间: 2019-9-4 01:45
不一定非得改。但如果是那种因为带电都能令结构扭曲得比较厉害的情况（可以用量化程序优化一下看看结构），那得做特殊考虑
作者Author: adong 时间: 2019-9-23 12:00
sob老师，请问如果我想做有机溶剂在金属表面的动力学模拟，如何添加有机溶剂？是跟添加水一样吗？有机溶剂的itp文件是自己构建吗？
作者Author: sobereva 时间: 2019-9-25 03:21
如果对初始要求高，先对有机溶剂做NPT跑出一个平衡后的溶剂盒子，用gmx solvate对真空区加溶剂
要求不高的话用packmol就行
是，用这里提到的工具。其中acpype我个人最推荐
几种生成有机分子GROMACS拓扑文件的工具
http://sobereva.com/266（http://bbs.keinsci.com/thread-428-1-1.html）
作者Author: adong 时间: 2019-9-25 11:29
好的，谢谢老师
作者Author: liuyjhx 时间: 2019-9-28 22:45
Program: gmx x2top, version 2018.4
Source file: src/gromacs/gmxpreprocess/x2top.cpp (line 546)
Fatal error:
No or incorrect atomname2type.n2t file found (looking for gromos54a7.ff)
老师请问我把atomname2type.n2t放在了gromos54a7.ff下面，但是为什么运行gmx x2top -f cnt10.pdb -o CNT10.top -ff select -nopbc -name CNT10 -kb 400000 -kt 600 -kp 150
后出现这样的错误？应该怎么修改？谢谢
作者Author: sobereva 时间: 2019-9-29 02:39
不在现场说不清楚
确认文件内容符合格式规矩，gmx安装过程完全正确，包括已经source了GMXRC文件
作者Author: adong 时间: 2019-10-6 21:23
老师，请问像LiPF6这样的化合物如何添加到混合溶剂里面？是分别插入Li+离子和PF6-离子吗？PF6-离子的拓扑信息如何添加呢？
作者Author: sobereva 时间: 2019-10-7 07:09
是。packmol插入
先想办法构造出PF6-的拓扑信息，然后在主top文件中include
作者Author: 黑色幽默pxj 时间: 2020-6-28 10:45
请问在gmx里面使用interface力场也这么操作嘛？我需要做一个蛋白和二氧化硅的复合物，文献中说二氧化硅是用的interface力场。
作者Author: DwyaneWan 时间: 2020-12-22 15:11
sob 老师如果要生成表面带羧基的CNT的话，在n2t文件中应该怎么去填写C.H.O的次序呢？
作者Author: sobereva 时间: 2020-12-23 03:31
n2t里不用管次序的事，只需要考虑内容
作者Author: DwyaneWan 时间: 2020-12-27 04:17
sob老师我有几个问题，
1.如果按照你的例子用gromos力场生成的羧基化的碳纳米管的top，原子电荷算下来是负的感觉不对劲。能否用RESP电荷呢？
2.对于碳纳米管，原子数为500个，如果我用xtb优化那么事先接枝上的羧基就会远离碳纳米管，那我还需要优化嘛？
3.对于500个原子的羧基碳纳米管，我用orca算单点，我核数给到16，内存给到64g还是会爆内存错误。能否用xtb算出的单点能去计算resp电荷呢，用于之后的分子动力学模拟呢？我用xtb算的单点，拟合的resp电荷除了和羧基相连的c电荷数高一些，碳纳米管的碳原子电荷都在零附近波动。所以想请问老师这样可行吗？
作者Author: sobereva 时间: 2020-12-27 06:38
1 我的文中并没有说原子电荷的事，不知道你怎么得到的
羧基整体的吸电子能量并不多强，如果纳米管部分为轻微的正值并无异常
用Multiwfn在恰当级别下算RESP电荷更好。如果算不动，索性让整个纳米管部分原子电荷都为0。
2 可能是初始结构不理想。另外，xtb支持不同方法，说清楚用的是哪种。
亦可尝试用半经验方法优化。
3 GFN-xTB的波函数质量太低劣，算RESP电荷的质量无法接受
作者Author: DwyaneWan 时间: 2020-12-28 00:44
老师你在博文中的碳纳米管全部都是C，所以你你在n2t文件里所有的的碳原子的电荷都为0。 我搜到一篇马普所做的关于羧基化的碳纳米管（https://doi.org/10.1039/C8FD00011E）。他得到RESP电荷的实验步骤the geometry was optimized and restraint electrostatic potential (RESP) charges 67 were calculated using R.E.D. (version III.52) with Gaussian 09 (RESP-A1: HF/6-31G* Connolly surface algo., 2 stage RESP fit qwt=.0005/.001). 由于我的碳纳米管有400多个原子，如果要用精度比较好的基组去算实在算不动，所以我用GFN-xTB算点拟合了的RESP电荷和他的作比较。
他的羧基C的charge为0.72，羧基上的=O的charge为-0.55，羧基上(-OH)的O和H的charge为-0.6和0.44。
我的拟合出的结果是羧基上的C的charge为0.68到0.88（因为我的羧基是randomly分布在碳管表面，他的羧基只在碳纳米管两段），羧基上的=O的charge为-0.71到-0.57，羧基（-OH）的O和H上的charge分别为-0.75到-0.55以及0.52到0.57。
如果我按照他的电荷用在我的拓扑上会使得我的碳纳米管的整体charge小于0，但是如果GFN-xTB算出的单点拟合出来的RESP的话，按老师的意思质量无法接受。我也跟类似的论文中做氧化石墨烯的含氧官能团的RESP电荷做了对比，发现我拟合出的电荷整体趋势上跟别人的是类似的。
如果我拟合出的RESP电荷质量不行的化，不知道还有别的方法呢？还是说就用acpype生成的电荷使用呢？
作者Author: sobereva 时间: 2020-12-28 06:59
acpype直接给的AM1-BCC电荷明显更不妥，那都是对一批有机小分子训练的校正参数，而且本身对静电势重现性就比RESP明显要差，此文第7节已经充分体现了：http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml
直接把羧基接上一个较小（或较短）的纳米管用Multiwfn算RESP电荷。然后把羧基以及相连那一些原子的电荷挪到实际体系上，其它纳米管原子电荷当成0。
作者Author: DwyaneWan 时间: 2020-12-28 12:56
好的，感谢老师。
作者Author: jun.ren 时间: 2021-1-15 17:13
sob老师，我执行第一串命令时出错：No or incorrect atomname2type.n2t file found (looking for gromos54a7.ff)，这是怎么回事？谢谢先！
作者Author: sobereva 时间: 2021-1-16 01:36
你若选择GROMOS54A7力场，必须自行创建和编辑gromos54a7.ff目录下的atomname2type.n2t之后x2top才能用
作者Author: 吧唧爱吃糖 时间: 2022-9-17 16:44
本帖最后由 吧唧爱吃糖 于 2022-9-17 17:28 编辑
老师您好，按照您的方法
1.在n2t文件上原子类型处写的是不属于oplsaa力场包含的原子类型CX，生成top文件之后在top文件中include了修改过的ffnonbonded-cnt.itp文件，文件中的C原子类型用的是top文件中的CX。之后grompp用得就是这个生成的top文件。之后报错 Atomtype CX not found应该如何解决？直接在atomtypes.atp中添加吗？
2.我在x2top命令中选择的是oplsaa力场，如果我想在n2t文件中选择力场中包含的原子类型，oplsaa力场中的atomtypes.atp中包含了很多碳的原子类型，我应该如何选择？（文献中未提到）麻烦sob老师了，若得回复十分感谢
作者Author: sobereva 时间: 2022-9-18 04:13
如果没有特殊要求，别基于OPLS-AA，原子类型定义得乱糟糟的
如今我都建议用sobtop（http://sobereva.com/soft/Sobtop）创建这类体系的拓扑文件，比x2top灵活强大得多
作者Author: 吧唧爱吃糖 时间: 2022-9-20 10:00
本帖最后由 吧唧爱吃糖 于 2022-9-20 10:02 编辑
老师您好，我是gromacs的新手，我直接进行nvt模拟的时候会出现这样的错误然后就停止了。
我想可能是结构不够好，所以我进行了能量最小化，但是会出现模型中本应在纳米管内的水分子全部都移动到纳米管外的情况。
我想问老师
1.nvt出现这样的错误是否是由于结构不够优化？如果直接进行nvt应该怎么解决报错？
2.为什么能量最小化会出现这样的问题，应该怎么解决。
附上两次的mdp文件。十分感谢sob老师。
作者Author: sobereva 时间: 2022-9-21 22:42
没结构文件不好说
肉眼看上去，垂直于管的方向盒子尺寸太小了。当前不该用EnerPres
作者Author: 吧唧爱吃糖 时间: 2022-9-22 22:06
本帖最后由 吧唧爱吃糖 于 2022-9-22 22:11 编辑
老师，这个是我的结构文件，我在min.mdp文件中按照您的知道把EnerPres项注释掉了，但是问题还是没有解决。而且输出力原子编号出现了47091，但是我Gro文件里没有这么多原子啊。。。
在SWNT.mdp文件中注释掉EnerPres后问题也没有解决
作者Author: sobereva 时间: 2022-9-23 05:18
除非文件搞错了，否则不可能提示47091号原子
注意拓扑文件和结构文件的对应关系
而且用这么长的管很莫名其妙
作者Author: RES 时间: 2022-10-14 15:33
sob老师，我第一次尝试做分子模拟，对于整个模拟过程还是很模糊，想跟您请教，我输入第一个命令，为什么显示的是没有g_x2top这个命令，请问这是为什么，是我Gromacs装错了吗？但是我输入gxm可以正常显示版本。
作者Author: sobereva 时间: 2022-10-14 20:07
g_x2top是5.0以前的版本才用的，早已过时
这种体系目前都建议用sobtop建立拓扑文件，看http://sobereva.com/soft/Sobtop里面的例子
作者Author: RES 时间: 2022-10-14 21:12
谢谢sob老师
作者Author: YPL 时间: 2022-10-14 22:11
我想请问一下，建立好碳纳米管后，里面塞入水分子。做能量最小化的时候，碳纳米管分裂成四块，这个是什么问题呢？
作者Author: sobereva 时间: 2022-10-14 22:16
搞清楚PBC是什么自然就明白了
不想分成四部分的话，一开始的坐标不要出现负值，否则会被挪到盒子另一头
作者Author: YPL 时间: 2022-10-14 22:18
好的，谢谢sob老师
作者Author: YPL 时间: 2022-10-14 22:41
再请教一下,使用Nanotube Modeler软件建立的碳纳米管默认坐标位置就带负值，请问这个怎么更改呢？
作者Author: sobereva 时间: 2022-10-15 00:31
gmx editconf结合-translate可以平移体系
作者Author: YPL 时间: 2022-10-15 13:14
好的，谢谢sob老师
作者Author: RES 时间: 2022-10-15 15:53
sob老师，我是想要进行一个简单的模拟，是在真空中，对于无氢的开口碳管,想要用分子动力学模拟查看氢原子或者离子更易于跟那个位置的C成键，老师您看我这种思路适合用分子动力学来模拟吗？
作者Author: sobereva 时间: 2022-10-16 04:06
至少不适合GROMACS，GROMACS不支持反应力场，没法描述动力学动力学过程中的成键
作者Author: RES 时间: 2022-10-16 21:03
谢谢sob老师
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |