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 (出处: 计算化学公社)

youknowdcf

Olive.Lv 发表于 2023-12-25 17:34
我用vmd给packmol建的模型附电荷，C原子我附的是6，因为用vmd打开其他优化好的文件，查C的电荷显示的是6， ...

正常原始石墨烯的C都是不带电的LJ粒子

Olive.Lv

我用vmd给packmol建的模型附电荷，C原子我附的是6，因为用vmd打开其他优化好的文件，查C的电荷显示的是6，但是我今天在一篇文献中看到说石墨烯的电荷设为0，请问这是为什么？

Olive.Lv

Lacrimosa 发表于 2023-12-15 16:34
你输入topo就可以看到手册呀，里面明明写的是topo retypedihedrals

是的，是我没看仔细，后来改过来了

Lacrimosa

Olive.Lv 发表于 2023-12-15 14:07
老师你好，我输入topo retypedihedral，然后出现报错Unknown topotools command 'retypedihedral'，如下 ...

你输入topo就可以看到手册呀，里面明明写的是topo retypedihedrals

Olive.Lv

Lacrimosa 发表于 2020-7-20 10:50
另外topotools直接猜出来的拓扑文件好像也挺准确的，直接运行以下命令就好了
topo clearbonds
topo clear ...

老师你好，我输入topo retypedihedral，然后出现报错Unknown topotools command 'retypedihedral'，如下图

neocc

怎么向cvff力场的聚合物的盒子里面添加spce或者tip4p模型水分子做溶剂化。


minimize最小化，跑md都是水分子报错。

报错有：
最小化过程中能量在30多步暴增到inf
md过程出现bonds atoms xxx xxx missing ...

sigusigu

tjuptz 发表于 2021-2-16 17:19
补充下，topotools的原始教程有较好的生成lammps data及拓扑处理实例。
另外，对于单个结构比较复杂，比如 ...

请问topotools的原始教程在哪里看呢？

Hollytree

，在使用nve系综下是没有报错的，但是使用nvt或npt时就会报错

Hollytree

尝试做了最简单的一个模型，但是这个模型和我之前在MS里面构建的是不一样的，之前是水分子全在slab模型的上面，不知道为何在生成data文件的以后就变成下面这种图形了

Hollytree

Kangtor

Hollytree 发表于 2021-11-30 22:48
可以详细的说一下聚合物+slab模型+真空构建部分吗？在做这方面的工作，但是一直报错，不知道是哪里出了问题 ...

具体报错的原因是？

Hollytree

可以详细的说一下聚合物+slab模型+真空构建部分吗？在做这方面的工作，但是一直报错，不知道是哪里出了问题。

mgqqlwq

xhzhangwsu 发表于 2021-9-26 13:55
A.U.>200的聚合物的话，可以试试PolyParGen：http://polypargen.com/。发表的文章： https://doi.org/10. ...

谢谢您的回复！但是PolyParGen有个问题就是输出文件中的原子顺序会完全打乱，我有多条链的时候就需要重新恢复顺序，但是我没有明白输出文件和输入文件中原子顺序的关系。

xhzhangwsu

mgqqlwq 发表于 2021-8-13 12:04
楼主gromacs部分我看的不是很懂，可否跟您请教一下。
我也是想生成聚合物的opls力场的拓扑文件，比如我的 ...

A.U.>200的聚合物的话，可以试试PolyParGen：http://polypargen.com/。发表的文章： https://doi.org/10.2477/jccjie.2018-0034