lammps和gromacs聚合物建模的方法

Kangtor

一. 本文更新时间：第一次 2020.7.17  10.46                                          
第二次 2020.7.20   19.25   经过较为复杂的体系（聚合物）测试，证明了topo的准确性。同时在原操作基础上补充了一些步骤。感谢LacrimosaL大佬
第三次  2021.8.12   21.14  忙了一年，主要是做聚合物相关领域（不是纯聚合物体系），总结一些经验。目前一直在用gromacs
第四次  2022.2.22  （无数个2，） 本次详细地更新一下聚合物界面的构造过程（上次写的也仓促了）


二. lammps方法简介
1. 具有图形窗口界面的软件        +      packmol        +    VMD  +  data文件检查
2. msi2lmp +  data文件检查
3.moltemplate + data文件检查
注：在第一次更新时间的文章中，本人只学习了第一种方法。另外两种之后会更新，或者哪位大佬有空补充一下


方法详解。
方法1   
各软件的作用：具有图形窗口界面的软件的作用-------生成.pdb 文件       packmol-----搭建所需的结构    VMD----data文件生成
具有图形窗口界面的软件----MS，gaussian等          packmol学习方法-----见参考来源5    tcl了解----见参考来源7   VMD学习----见参考来源8 9
本人感觉此方法需要tcl的知识，所以有必要学习了解tcl。本方法主要操作难点在有相关tcl文本命令的输入，即“VMD----data文件生成” 这一步
tcl输入可采用以下两种方法：（1）使用tcl脚本文件    （2）在VMD TKConsle（命令台）下输入
第一种输入方法详解：
采用的例子来源于参考1。本人感觉参考1中的例子不错（乘机夸夸大佬）
步骤：将water.pdb文件和water.tcl文件放在VMD的目录下，然后在TKConsle中输入source water.tcl ，然后端起你的小茶杯喝口水等待几秒，data文件就生成了
问题：新手不会写脚本文件 ，这真的是"女的看了流泪，男的看了默哀" （原谅我的黑色幽默）
第二种输入方法详解：
直接在VMD TKConsle输入。虽然比较麻烦，但是相对于第一种来说，本人觉得学习较快。
步骤：先在VMD中加载构型water.pdb文件。如果没有立方体的那个框框，在TKConsle中输入pbc box。 如下：


在这一步中应当设置好原子类型，可以看看最下面的测试

接下来，通过atomselect命令来编辑原子属性参数如质量，电荷等。 我将通过一个于本贴主题无关的操作实例来演示atomselect命令如何使用。首先，加载构型是一个水分子, 在Graphics representation中设置如下Drawing method-----VDW  coloring method-----name。所以水分子是这个样子。    然后，通过atomselect命令来改变不同种类的原子半径，从而改变图形构型（只改变氧原子即可）。 在tkconsole中的命令如下：

  图形构造变成了 。这样就完成了。可以通过下图中的命令查看半径的信息。

那如何设置氧的电荷是不是很明白了，只要将radius改为charge 即可


再接下来开始解决第三张图中 Angles: 0  Dihedrals: 0  Impropers: 0的问题，同样还是以一个水分子为例。在未开始之前，先看看一个水分子的结构的分析 ，可以看出原先的pdb文件中含有键的信息。而要推出角度，二面角等信息，是需要键的信息的，例如 。
本人在阅读各位大佬的帖子后发现，需要在原先pbd文件的基础上删除掉所有键的信息，再重新建立bond topology。但是实际操作后发现不需要（所以在这里打个问号）？
按照各位大佬们的帖子，输入的命令可为：虚线后面为解释
topo clearbonds ---------deletes all bonds
mol bondsrecalc top--------will guess bonding from atomic radii
topo guessangles-------- will derive angledefinitions from the bond topology
topo guessdihedrals-------will derive dihedraldefinitions from the bond topology
topo guessimpropers-------- will derive improperdefinitions from bonds (three bonds to one atomand less than 5% out of plane).
topo retypebonds-------resets all bond types，下同
topo retypeangles
topo retypedihedrals
到这一步算是完成了， 可以通过mol reanalyze top从新查看结构信息。则一个水分子的结构信息变为

最后输出data文件，可通过此命令   例如：topo writelammpsdata mysystem.data

完成后，一定要检查data文件信息


较为复杂的体系（聚合物链）的测试过程结果:
选择作为测试案例的文章为 Quantum Chemistry Based Force Field for Simulations of Poly(vinylidene fluoride) Oleksiy G. Byutner and Grant D. Smith*
选择的聚合物体系为聚偏二氯乙烯（pvdf）   合成单体：vdf
首先，在MS构建了单体单元，编辑了相应的“name”，并在label显示出来
     
原子的name属性是和原子种类应该是相同的



接下来，测试过程和结果
   此图为命令的输入
  上图为clearbonds之后的显示，下图命令输入完后得到的结果。（可以实际证明 键的拓扑结构 的结果是直接显示在图形窗口的。这样，就可以直观地观察到guessbonds的正确性）

通过与论文中的结果比较：（以下所有左图为论文内容，右图为测试结果）
原子种类：      注：图中Ch3 和 Cf3为分子链末端的碳原子。如果事先设置这俩种原子类型会使后续过程变麻烦

键的类型：   
键角：   

二面角    

可以从以上数据看出topo的结果还不错
但是，
   从图中可以看出分子id出现了问题。所以后续检查非常重要。


三. gromacs聚合物建模方法


本人一直研究的体系是线性聚合物链。以下给出我常常使用的方法以及一些研究心得。


gromac模拟最不好得到的文件就是itp文件。对于聚合物来说，十分的麻烦，自己写是不切实际的。


利用pdb2gmx去产生itp文件。一般来说，只需要利用r2b和rtp文件就够了，个人不喜欢利用hdb文件来自动补氢，因为我最初的文件的原子名是按照我自己习惯设定的。


r2b文件是设定残基名转换关系的。它将原文件中的残基名改成了rtp中设定的残基名。


rtp文件是利用pdb文件来构建拓扑信息的。这个文件中的构建聚合物拓扑文件的内容是要自己写的。最简单的例子是需要写三种残基（两个末端残基和中间重复的残基。）中各个原子的信息以及键接关系。每个残基中至少需要两个字段【atom】和【bonds】。对于【atom】字段，主要是记录：第一列---原子名（要和pdb文件中的一样，顺序也一样），第二列----力场的原子类型（例如opls-aa ff的opls_186），第三列是原子电荷，第四列是序数。【bonds】字段是需要将键接关系写出来的，要注意利用  +  或  -  来指定与下一个或上一个残基的键接关系。
在这里说一些我个人的关于聚合物原子电荷的设定的感受。我用的基本上就是opls-aa ff，这个力场的原子电荷有 最初开发的；1.25*CM5;1.14* CM1A或1.14*CM1A-LBCC 三种。我个人对聚合物用的都是最初开发的；研究小分子采用CM1A或CM1A-LBCC。


再说下对于聚合物和某些物质（如纳米粒子）的界面构造的问题。
一般有两种构造方式。
一种是 slab模型+聚合物+真空
第二种是真空+slab模型+聚合物+slab模型+真空

感觉各有利弊。第一种遇到聚合物和溶剂的混合体系会出现小分子往上飞的情况。研究聚合物的动力学是需要高温加速链段运动，但是小分子就会被蒸发，这种情况就利用第二种。第二种的话盒子会建的很大，计算比较慢，但也能接受。



Gromacs的几个建模的命令还是挺好使的，比如editconf，solvate，insert-molecules等。熟练使用这几个命令基本上都能够构造大多数模型了。（做科学问题最难受的就是缺程序，尤其是完美契合自己的程序，希望自己以后能够写出一个相关的计算/分析程序来。）

我个人认为聚合物-界面模型的构造最麻烦的点在于：将聚合物分子链插入到固定体积的盒子中。因此，我建议了一种 Gromacs基多步构造/松弛方案  供大家选择。

以下通过聚合物A和slabB的 “真空+slab模型+聚合物+slab模型+真空” 界面模型来进行讲解：

（1）首先，搭建一个“B+真空+B” 界面模型。 真空区域的体积采用均相A模型的密度来进行构建
（2）利用insert-molecules命令将一部分A插入到“B+真空+B” 界面模型中。 可以选择多种A的构象去插入
（3）将“B+部分A+B” 界面模型在一定温度下进行NVT run。此时，会发现聚合物链吸附到了板上，留出了一部分真空区。
（4）重复步骤 2和3，直到获得一个完整的结构，即，“B+A+B”模型。
（5）通过一个合理的平衡方法对体系进行平衡




         
参考来源：
1.使用VMD中的topo模块构建 .data, .psf, .top文件 http://bbs.keinsci.com/forum.php ... =4753&fromuid=19927 (出处: 计算化学公社)
2.聚合物模拟的力场参数如何选取和获得的问题 http://bbs.keinsci.com/forum.php ... =6522&fromuid=19927 (出处: 计算化学公社)
3.LAMMPS学习资源整理 http://bbs.keinsci.com/forum.php ... id=73&fromuid=19927 (出处: 计算化学公社)
4.msi2lmp工具的使用问题 http://bbs.keinsci.com/forum.php ... 16032&fromuid=19927 (出处: 计算化学公社)
5.分子动力学初始结构构建程序Packmol的使用 http://bbs.keinsci.com/forum.php ... 12549&fromuid=19927 (出处: 计算化学公社)
6.如何得到模拟固体的原子坐标，类型，电荷以及键的信息？ http://bbs.keinsci.com/forum.php ... =9328&fromuid=19927 (出处: 计算化学公社)
7.Tcl迅速入门教程 http://bbs.keinsci.com/forum.php ... d=157&fromuid=19927 (出处: 计算化学公社)8.VMD里原子选择语句的语法和例 http://bbs.keinsci.com/forum.php ... 14267&fromuid=19927 (出处: 计算化学公社)
9.谁有VMD的使用教程，能否分享哈，谢谢 http://bbs.keinsci.com/forum.php ... =1403&fromuid=19927 (出处: 计算化学公社)

fineren

不错，收藏了，楼主不要中断更新

sobereva

贴图方式不对，别人看不到。不要直接把图片往窗口里复制，注意看置顶的新社员必读贴了解怎么正确上传和插入图片，请重新编辑帖子

Kangtor

sobereva 发表于 2020-7-18 10:18
贴图方式不对，别人看不到。不要直接把图片往窗口里复制，注意看置顶的新社员必读贴了解怎么正确上传和插入 ...

好的，社长

Lacrimosa

对于方法1补充一下，水的体系中只有一种bond/angle type用topo guessbond/angles没有问题，但是对于较复杂的体系就需要一个个添加成键关系。举个例子，用于CLAYFF中的硅羟基Si-O-Htopotools中bondtype应该是根据原子类型进行判断的，而粘土矿物中原子类型比较多，并且有一部分是没有成键作用的，所以直接guess bond会出现很多多余的键。所以才会有以下特例
set ilist [[atomselect top "type oh"] get index] ；选择所有羟基氧（原子类型为oh）的index
foreach i $ilist {    ；对ilist中的每个index进行循环
set j [atomselect top "type ho and within 2 of index $i"];  选择距离index为i的羟基氧原子 0.2nm范围内的羟基氢原子（type ho）
set k [atomselect top "{not type ho oh} and within 2 of index $i"] ;选择距离index为i的羟基氧原子 0.2nm范围内的非羟基氢原子（type ho）
topo addbond [$j get index] $i -bondtype 2 ; 连接H-O键，bondtype设为2
topo addangle [$j get index] $i [$k get index] 2 ;连接M-O-H键键角，angletype设为2，与bondtype不同此处不需要加-angletype直接输入2即可
}

Kangtor

Lacrimosa 发表于 2020-7-18 12:31
对于方法1补充一下，水的体系中只有一种bond/angle type用topo guessbond/angles没有问题，但是对于较复杂 ...

我用乙烷尝试了上述的循环命令，得到了这样的结果
>Main< (VMD) 27 % set a [[atomselect top "element C"] get index]
1 4
>Main< (VMD) 28 % foreach b $a {
set j [atomselect top "element H and within 1.6 of index $b"]
set k [atomselect top "{not element H} and within 1.6 of index $b"]
topo addbond [$j get index] $b -bondtype 2
topo addangle [$j get index] $a [$k get index] 2
}
atomsel : setbonds: Need one bondlist for each selected atom
然后不知道怎么输入命令了

Lacrimosa

Kangtor 发表于 2020-7-19 17:59
我用乙烷尝试了上述的循环命令，得到了这样的结果
>Main< (VMD) 27 % set a [[atomselect top "element  ...

因为一个碳连接了三个氢，在addbond中[$j get index]表示的是三个氢组成的list，所以会报错。显然这样做更合适一些
set a [[atomselect top "element H"] get index]
foreach b $a {
set j [atomselect top "element C and exwithin 1.6 of index $b"]
topo addbond [$j get index] $b -bondtype 2 ;
}
topo addbond list [[atomselect top "element C"] get index] -bondtype 1;

Lacrimosa

另外topotools直接猜出来的拓扑文件好像也挺准确的，直接运行以下命令就好了
topo clearbonds
topo clearangles
topo cleardihedrals
topo clearimpropers

topo guessbonds
topo guessangles
topo guessdihedrals

topo retypebonds
topo retypeangles
topo retypedihedrals

对于比较复杂的有机物，最好先根据成键关系判断好原子类型，然后再用guessbonds

Kangtor

Lacrimosa 发表于 2020-7-20 10:45
因为一个碳连接了三个氢，在addbond中[$j get index]表示的是三个氢组成的list，所以会报错。显然这样做 ...

topo addbond list [[atomselect top "element C"] get index] -bondtype 1;
感觉这一句是在建立碳碳键，但是这一句我一直在报错

Lacrimosa

Kangtor 发表于 2020-7-20 17:41
topo addbond list [[atomselect top "element C"] get index] -bondtype 1;
感觉这一句是在建立碳碳键 ...

是建立C-C键的，这里这样写只适用于乙烷，因为乙烷只有两个碳。

zxs1127

感谢，收货很大

kazusa

十分感谢，很有帮助

tjuptz

补充下，topotools的原始教程有较好的生成lammps data及拓扑处理实例。
另外，对于单个结构比较复杂，比如交联聚合物或者生物大分子或者很复杂的有机物、无机物，其实用topotools里guessbonds或者vmd加载分子时autobonds自己根据距离标准判断的成键关系，不一定靠谱，肉眼又不好发现，需要谨慎一点。

Heritis

请问一下，我在TK Console中输入“pbc box”，回车后提示“ERROR) Suspicious pbc side length (a=0.000000 b=0.000000 c=0.000000). Have you forgotten to set the pbc parameters?”是什么原因呢？谢谢。

Lacrimosa

Heritis 发表于 2021-3-15 12:44
请问一下，我在TK Console中输入“pbc box”，回车后提示“ERROR) Suspicious pbc side length (a=0.000000 ...

说明你的结构里缺少周期性边界参数