二氧化硅力场文件

Lacrimosa

参照charmm27力场的形式以及 10.1021/cm500365c 此文的参数写了一个Silica的力场文件供gmx使用，放入gmx主目录下的share/gromacs/top文件夹中即可通过x2top生成相应的top文件。在编写过程中发现在写入非键作用参数（ffnonbonded.itp）的时候对照了一下 10.1021/cm500365c文中charmm27.ff和gmx手册中的LJ作用函数，发现二者形式略有差异，原文中的参数sigma需要乘2^(-1/6)才符合gmx中定义的sigma。而INTERFACEFF中给出的charmm27_interface_v1_5.prm中却直接将此参数除以了2，应该是出错了。
这次对之前的文件稍微进行了修改，删除了多余的文件，简化了冗长的原子类型，参数没有变。
链接：https://pan.baidu.com/s/1v72ZMjgcASvC4Gfas-qBkw
提取码：b9c0

2023.3.14日更新，本文所述方法较为繁琐且容易出错，当前已不建议使用。


2024.5.18日更新, 在最下方新增了Sobtop输入文件，该输入文件暂时不适用于解离体系，并且文件中不含二面角参数，如有需要请参考原文自行添加。

67楼提供了二氧化硅表面接枝其他官能团（如：-RCH3，-RCOOH，-RNH2等）体系的方案







2020.10.6日更新，对ffnonbonded.itp中的σ参数进行了修改，并新增了charmm36-mar2019-Silica.ff和amber14sb_OL15-Silica.ff, 这两个文件放入gmx可以直接通过x2top生成top文件，无需额外修改，这二者的差别在于ffbonded.itp中Si-O的参数r0, 详情见原文。
在下面对使用方法做一详细说明：
压缩包内有以下文件：
1.atomtypes.atp:用于查阅原子类型以及电荷
2.ffnonbonded.itp:内部包含原子类型 原子序数 质量 范德华参数
3.ffbonded.itp:键伸缩项、键角项参数
4.silicaiff.n2t:用于生成top文件
5.silica.top:top的示例文件
1.获取二氧化硅的pdb文件
获取二氧化硅的pdb方法很多，可以从晶体数据库（例如COD、AMCSD等）下载cif，然后用M$打开，对晶胞进行切面扩胞，用氢原子把表面悬挂的氧补足。具体做法：
Build-Surface-Cleave Surface可以选择切的面、位置和厚度，根据自己的意向切好后，全选原子，点击上方的Adjust Hydrogen按钮，就是一个H，把表面的悬挂氧用氢补上。Build-Symmetry-Supercell进行扩胞，扩展到你需要的大小，最后Build-Crystal-Build Vacuum Slab选择晶面的位置、上方的真空区之类的参数，设定好后就可以导出成pdb文件了。


2.生成top文件
使用的时候需要配合charmm27力场,首先将ffbonded.itp和ffnonbonded.itp中的参数复制到charmm27ff内相应文件的相应位置。将n2t文件放入charmm27.ff文件夹内。使用x2top（gmx x2top -f name.pdb -o name.top -ff select）生成silica的文件。需要注意的是，如果silica中有离子需要先把离子删除（vmd输入：::TopoTools::selections2mol [atomselect top "not element Na"]即可）之后再使用x2top指令。

如果出现以下报错：
Can not find forcefield for atom O-xx with x bond （x代表数字）
这说明该原子距离其他原子过远，解决方法是：
测量该原子与其他原子的键长，参照n2t文件的写法把键长相应的修改好
例如：Si      SC4  1.10    28.0860  4  O   0.16         O  0.16     O 0.16     O 0.16  ;for SiO2 in bulk or SiOH ；如果Si-O距离为0.14则把0.16改为0.14即可，单位为nm
生成top文件之后，把[moleculetype]到[system]之前的部分复制到新的文件中，命名为silica.itp。top文件的格式参照silica.top即可。
对于带有钠离子的二氧化硅，按前述方法删除钠离子生成top文件之后，再把钠离子全部放到silica的后面。用VMD可以快速实现。
::TopoTools::selections2mol [atomselect top "not element Na"]
::TopoTools::selections2mol [atomselect top "element Na"]
::TopoTools::mergemols [list 1 2] ;1和2分别对应的是silica和钠离子的ID
输入：[atomselect top "element Na"] num查看钠离子的数量，打开charmm27中的ions.itp找到与钠离子对应的[moleculetype],将其名称和数量填在top文件[molecules]中silica的下方。
例：

[ molecules ]                                ；这里各种分子上下的顺序需要与pdb文件一致，否则会出现报错
; Compound        #mols     
silica             1                             ；silica的[moleculetype]+数量
NA                 96                          ；钠离子的
SOL               2980                       ；水分子的

tjuptz

您好，可否把atomname2type.n2t分享一下，谢谢

Lacrimosa

tjuptz 发表于 2019-10-21 21:47
您好，可否把atomname2type.n2t分享一下，谢谢

Si      silica_1si  1.10    28.0860  4  O   0.16         O  0.16     O 0.16     O 0.16  ;SiO2 or SiOH
;Si      silica_1si  1.10    28.0860  3  O   0.16         O  0.16     O 0.16             ;SiO2 or SiOH
;Si      silica_1si  1.10    28.0860  2  O   0.16         O  0.16                        ; SiO2 i or SiOH
Si      silica_2si  0.425   28.0860  4  O   0.16         O  0.16     O 0.16     H 0.16  ;3O-Si-H
O       silica_4o  -0.900   15.9994  1  Si  0.16                                                                     ; SiO-
O       silica_1o  -0.675   15.9994  2  Si  0.16         H   0.11                                             ;Si-OH
O       silica_2o  -0.55    15.9994  2  Si  0.16         Si  0.16                                             ;SiO2
H       silica_1h   0.40    1.0080   1  O   0.11                                                                     ;Si-OH
;H       silica_1h   0.40    1.0080   1  O   0.10                                                                     ;Si-OH
H       silica_2h   0.40    1.0080   1  Si  0.16                                                                   ;Si-H
                    

tjuptz

再请教，关于体相原子电荷设置为0的合理性，楼主是否得到结论了呢？

Lacrimosa

tjuptz 发表于 2019-10-22 15:19
再请教，关于体相原子电荷设置为0的合理性，楼主是否得到结论了呢？

关于这点我询问了一下作者，对方的意思大概是他们主要想研究表面和有机物的相互作用，所以为了便捷就做了那样的处理，电荷参数合理的话定性研究似乎是没什么大问题的。就算不把体相中的原子电荷设为0，在freeze住体相之后，也就不考虑体相中的相互作用了吧？那样的话效果是不是一样呢？

tjuptz

Lacrimosa 发表于 2019-10-22 17:42
关于这点我询问了一下作者，对方的意思大概是他们主要想研究表面和有机物的相互作用，所以为了便捷就做了 ...

但对算静电和范德华能量还是有影响的，如果要算能量就不能设为零吧

Lacrimosa

tjuptz 发表于 2019-10-22 19:15
但对算静电和范德华能量还是有影响的，如果要算能量就不能设为零吧

嗯，定量的话可能还是会有影响，但是体相中的原子距离表面比较远，产生的非键相互作用不会那么强，定性研究应该足够准确了。

tjuptz

Lacrimosa 发表于 2019-10-23 11:55
嗯，定量的话可能还是会有影响，但是体相中的原子距离表面比较远，产生的非键相互作用不会那么强，定性研 ...

有道理

tjuptz

Lacrimosa 发表于 2019-10-22 14:13
Si      silica_1si  1.10    28.0860  4  O   0.16         O  0.16     O 0.16     O 0.16  ;SiO2 or S ...

请问，这个n2t中 ; 的注释是冗余的规则吗？为何写冗余的规则

Lacrimosa

tjuptz 发表于 2019-10-24 22:47
请问，这个n2t中 ; 的注释是冗余的规则吗？为何写冗余的规则

直接写在后面的是方便自己看懂每个原子类型代表什么，有两行注释掉了是因为当时输入的结构不完整（出现了只连接三个氧或两个氧的硅原子）我为了处理方便就那样写了，其实删掉比较好

tjuptz

Lacrimosa 发表于 2019-10-24 23:40
直接写在后面的是方便自己看懂每个原子类型代表什么，有两行注释掉了是因为当时输入的结构不完整（出现了 ...

谢谢楼主

tjuptz

楼主，我看了下charmm27和charmm36的ff，似乎力场形式是一致的，那么您提供的该力场也可以与charmm36兼容咯？
另，我看charmm36里自带了石英和硅酸盐的参数，也很方便了

Lacrimosa

tjuptz 发表于 2019-10-26 17:16
楼主，我看了下charmm27和charmm36的ff，似乎力场形式是一致的，那么您提供的该力场也可以与charmm36兼容咯 ...

非键作用函数一样的话一般都是可以通用的，而且根据原文，这个力场参数在拟合的时候就考虑了charmm的函数形式。charmm36力场好像没有n2t文件，生成top的时候不知道用pdb2gmx是否可行

tjuptz

Lacrimosa 发表于 2019-10-28 20:35
非键作用函数一样的话一般都是可以通用的，而且根据原文，这个力场参数在拟合的时候就考虑了charmm的函数 ...

用 x2top 自己写n2t就行啦

黑色幽默pxj

请问楼主，蛋白质和二氧化硅的复合体能用这个文件么？怎么跟gmx里的charmm27使用呢