计算化学公社

 找回密码 Forget password
 注册 Register
Views: 38237|回复 Reply: 41

[Lammps] lammps和gromacs聚合物建模的方法

  [复制链接 Copy URL]

108

帖子

0

威望

2716

eV
积分
2824

Level 5 (御坂)

发表于 Post on 2020-7-17 12:10:20 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 Kangtor 于 2022-2-22 18:11 编辑

一. 本文更新时间第一次 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。 如下:
贴1.PNG 贴2.PNG 贴3.PNG

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


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

下载1.png   图形构造变成了 贴5.PNG 。这样就完成了。可以通过下图中的命令查看半径的信息。 下载2.png

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


再接下来开始解决第三张图 Angles: 0  Dihedrals: 0  Impropers: 0的问题,同样还是以一个水分子为例。在未开始之前,先看看一个水分子的结构的分析 下载3.png ,可以看出原先的pdb文件中含有键的信息。而要推出角度,二面角等信息,是需要键的信息的,例如 下载4.png
本人在阅读各位大佬的帖子后发现,需要在原先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从新查看结构信息。则一个水分子的结构信息变为 下载5.png

最后输出data文件,可通过此命令 下载6.png   例如: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显示出来
测试11.PNG      
原子的name属性是和原子种类应该是相同的
测试13.PNG 测试2.PNG


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

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

键的类型: 测试15.PNG    测试7.PNG
键角:    测试8.PNG

二面角 测试16.PNG     测试9.PNG

可以从以上数据看出topo的结果还不错
但是
测试12.PNG 测试10.PNG    从图中可以看出分子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 (出处: 计算化学公社)













贴6.PNG
测试1.PNG

评分 Rate

参与人数
Participants 16
eV +70 收起 理由
Reason
Gzh_NJ + 5 谢谢
风清行 + 4 赞!
woller + 3 谢谢
panernie + 5 精品内容
zuoyongkang1 + 3 赞!
遗世之光 + 3 赞!
zero! + 3 谢谢
wwj + 5 好物!
学!!! + 3 谢谢分享
peachRL + 4 谢谢分享!
agent99 + 5 GJ!
bobosiji + 5 满分鼓励 期待接下来的更新
ABetaCarw + 5 满分鼓励 期待接下来的更新
tiandikuoyuan + 5 谢谢分享
sobereva + 10
Lacrimosa + 2 希望能持续分享经验

查看全部评分 View all ratings

116

帖子

0

威望

1608

eV
积分
1724

Level 5 (御坂)

发表于 Post on 2020-7-18 10:17:12 | 显示全部楼层 Show all
不错,收藏了,楼主不要中断更新

4万

帖子

99

威望

4万

eV
积分
89888

管理员

公社社长+计算化学玩家

发表于 Post on 2020-7-18 10:18:50 | 显示全部楼层 Show all
贴图方式不对,别人看不到。不要直接把图片往窗口里复制,注意看置顶的新社员必读贴了解怎么正确上传和插入图片,请重新编辑帖子
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

108

帖子

0

威望

2716

eV
积分
2824

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2020-7-18 11:20:27 | 显示全部楼层 Show all
sobereva 发表于 2020-7-18 10:18
贴图方式不对,别人看不到。不要直接把图片往窗口里复制,注意看置顶的新社员必读贴了解怎么正确上传和插入 ...

好的,社长

240

帖子

0

威望

2070

eV
积分
2310

Level 5 (御坂)

Nerv

发表于 Post on 2020-7-18 12:31:15 | 显示全部楼层 Show all
本帖最后由 Lacrimosa 于 2020-7-20 11:20 编辑

对于方法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即可
}
God's in his heaven,all is right with the world

108

帖子

0

威望

2716

eV
积分
2824

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2020-7-19 17:59:09 | 显示全部楼层 Show all
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
然后不知道怎么输入命令了

240

帖子

0

威望

2070

eV
积分
2310

Level 5 (御坂)

Nerv

发表于 Post on 2020-7-20 10:45:40 | 显示全部楼层 Show all
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;

God's in his heaven,all is right with the world

240

帖子

0

威望

2070

eV
积分
2310

Level 5 (御坂)

Nerv

发表于 Post on 2020-7-20 10:50:29 | 显示全部楼层 Show all
本帖最后由 Lacrimosa 于 2020-7-20 11:18 编辑

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

topo guessbonds
topo guessangles
topo guessdihedrals

topo retypebonds
topo retypeangles
topo retypedihedrals

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

God's in his heaven,all is right with the world

108

帖子

0

威望

2716

eV
积分
2824

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2020-7-20 17:41:47 | 显示全部楼层 Show all
Lacrimosa 发表于 2020-7-20 10:45
因为一个碳连接了三个氢,在addbond中[$j get index]表示的是三个氢组成的list,所以会报错。显然这样做 ...

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

240

帖子

0

威望

2070

eV
积分
2310

Level 5 (御坂)

Nerv

发表于 Post on 2020-7-21 17:48:57 | 显示全部楼层 Show all
Kangtor 发表于 2020-7-20 17:41
topo addbond list [[atomselect top "element C"] get index] -bondtype 1;
感觉这一句是在建立碳碳键 ...

是建立C-C键的,这里这样写只适用于乙烷,因为乙烷只有两个碳。
God's in his heaven,all is right with the world

13

帖子

0

威望

143

eV
积分
156

Level 3 能力者

发表于 Post on 2020-8-19 11:55:56 | 显示全部楼层 Show all
本帖最后由 zxs1127 于 2020-11-4 21:03 编辑

感谢,收货很大

7

帖子

0

威望

173

eV
积分
180

Level 3 能力者

发表于 Post on 2021-2-1 11:09:38 | 显示全部楼层 Show all
十分感谢,很有帮助

489

帖子

1

威望

3218

eV
积分
3727

Level 5 (御坂)

发表于 Post on 2021-2-16 17:19:15 | 显示全部楼层 Show all
本帖最后由 tjuptz 于 2021-2-16 17:23 编辑

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

11

帖子

0

威望

107

eV
积分
118

Level 2 能力者

发表于 Post on 2021-3-15 12:44:39 | 显示全部楼层 Show all
请问一下,我在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?”是什么原因呢?谢谢。

240

帖子

0

威望

2070

eV
积分
2310

Level 5 (御坂)

Nerv

发表于 Post on 2021-3-15 13:08:28 | 显示全部楼层 Show all
Heritis 发表于 2021-3-15 12:44
请问一下,我在TK Console中输入“pbc box”,回车后提示“ERROR) Suspicious pbc side length (a=0.000000 ...

说明你的结构里缺少周期性边界参数
God's in his heaven,all is right with the world

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-2 22:42 , Processed in 0.262460 second(s), 25 queries .

快速回复 返回顶部 返回列表 Return to list