【求助&计算资料整理】晶体-有机小分子双层异质结data文件生成

PLwang

【摘要】我是lammps的新手，最近正在使用Lammps做一个“多个有机小分子在晶体表面的吸附和分散性模拟”。我希望使用moltemplate将多个lt文件整合，生成完整的lammps输入文件。但是，在生成lt文件时，我无法准确地通过各类软件向其中提供合理的力场和电荷参数，具体的问题见下帖子。我会详细的列出我的建模和计算过程，一方面希望各位老师对我的计算过程进行检查，指出一些常见错误；另一方面，我也希望我的建模和计算过程（主要是生成data文件的过程，特别是一些我作为初学者走过的弯路），可以对后续的学习lammps的新人有帮助。对于后一点，我无法承诺自己的计算方法完全正确，但我会尽量提供我能获取到的原始资料，对其进行整理，希望对新手入门有一定帮助。为了方便各位阅读，我在每一段前均有标准【整理分享】或【问题求助】，前者主要是方便后来的计算新手入门，里面会写一些常见的方法、我对方法的理解和评价；后者主要是我目前存在问题的地方，希望能得到大家的解惑。非常感谢！
一、【整理分享】通过materials studio和msi2lmp生成data文件
具体方法：

• 使用MS导出car文件：对于已经建好的模型，选择Module-Forcite-Calculation-Energy选项卡内，首先选择Forcefiled，对于力场的选择（分享一篇关于力场对比的总结 - 分子模拟 (Molecular Modeling) - 计算化学公社）。选择力场后，点击More，取消所有的Calculate automatically，点击Calculate。随后，一定要检查体系内的分子是否被分配了合适的力场，MS直接分配的力场往往是不完全正确的。对于不正确的力场，请手动调整，并重新Assign。确认力场无误之后，对于电荷，使用Use current即可。最后，导出*car和*mdf文件；
• 使用msi2lmp生成data文件：生成的car和mdf文件应存放在同一个目录下。使用CMD，cd到对应的文件夹，输入：
  

1. msi2lmp.exe FILENAME -i -class 1 -frc cvff > data.FILENAME

复制代码

提醒：

• MS导出car和mdf文件，只需要选择car即可，两个文件会一同生成；
• 第二步使用CMD，不能使用Power Shell。后者经常会错误识别<（该提醒同样适用于后续）；
• msi2lmp支持的力场比较有限，往往只能使用MS和msi2lmp两者都支持的力场才能完成该步骤；
• 对于不完整的力场文件，需要后期手动添加其他数据。
  

在我看来，MS的forcite是一个对新手用户友好、但往往没有那么高效的工具。对于我来说，我的体系是NiO晶体和有机小分子（含有N C S P H等元素），就无法使用该方法直接生成比较合理的data文件，例如，MS无法识别二价Ni；msi2lmp支持的有限的力场，也无法满足我的使用需求。
参考：Molecular Dynamics Simulation——MS建模后转data文件以及势参数处理的小技巧 - 知乎



二、【整理分享】通过packmol随机生成多个有机小分子的data文件
具体方法：

• MS/Gaussian导出有机小分子的pdb文件；
• 按照自己的实际需求，编写合适的inp文件，具体方法请参考：分子动力学初始结构构建程序Packmol的使用 - 思想家公社的门口：量子化学·分子模拟·二次元
• 将二者放在同一目录，输入代码：
  

1. packmol < FILENAME.inp

复制代码

提醒：

• 编写inp文件时，一定要注意细节。例如，该文件中的，"0."的小数点一定不能忽略，因为这标识你的上输入是一个浮点数，否则系统会报错
  

1. tolerance 2.0
   
2. add_box_sides 1.2
   
3. output methanol_box.pdb
   
4. structure methanol.pdb
   
5.   number 1000
   
6.   inside cube 0. 0. 0. 50.
   
7. end structure

复制代码

   2. 用packmol建模，相当于只是把原子放在了盒子里，原子是不带任何信息的（比如此种原子的电荷、周期性排列方式等），packmol只给出含有坐标信息的结构文件，其它都不体现；


三、【整理分享】使用VMD计算、生成具有连接关系等信息的data文件

具体方法：

• 使用MS建模，结果导出PDB文件。在VMD中打开；
• 打开VMD-Extensions-TK console；
• 输入以下代码，以下代码会生成键、角、二面角关系并写入data文件；
  

1. topo clearbonds ---------deletes all bonds
   
2. mol bondsrecalc top--------will guess bonding from atomic radii
   
3. topo guessangles-------- will derive angledefinitions from the bond topology
   
4. topo guessdihedrals-------will derive dihedraldefinitions from the bond topology
   
5. topo guessimpropers-------- will derive improperdefinitions from bonds (three bonds to one atomand less than 5% out of plane).
   
6. topo retypebonds-------resets all bond types，下同
   
7. topo retypeangles
   
8. topo retypedihedrals

复制代码

提醒：

• VMD生成的data文件，缺少力场信息，需手动补充；
• topotools需要下载，并将其存储在VMD对应的目录下TopoTools Plugin, Version 1.8。
  

参考：lammps和gromacs聚合物建模的方法 - 分子模拟 (Molecular Modeling) - 计算化学公社；使用VMD中的topo模块构建 .data , .psf , .top文件 - 分子模拟 (Molecular Modeling) - 计算化学公社

四、【整理分享】使用amber生成有机分子的data文件
具体方法：

• Ambertools和acpype的安装 - 知乎。ambertools比较麻烦的是，其支持linux环境。对于完全没有linux操作技术、也没有linux系统的新手来说，可以尝试这个网站。该网站可以为你提供一个虽然不是很好用、但免费的linux环境，新手可以在上面任意尝试而不用担心会引发服务器的问题。精选项目课程_IT热门课程_蓝桥云课课程 - 蓝桥云课。
• 准备好分子的pdb结构，输入：
  

1. antechmber -i FILENAME.pdb -fi pdb -o FILENAME.data -fo data -c bcc -at gaff

复制代码

提醒：

• amber是一个适用于有机体系的力场，关于如何选择amber力场，请参阅AMBER程序支持的力场一览 - 思想家公社的门口：量子化学·分子模拟·二次元；
• amber力场不能对packmol产生的复杂体系生产力场。如果你体系中含有多个分子，其会报错。
  

五、【整理分享】使用ltemplify.py将data文件转化为lt文件，方便后续moltemoplate进行处理
具体方法：

• moltemplate/doc/doc_ltemplify.md at master · jewettaij/moltemplate · GitHub。请阅读该文件的README文档。
• 参考方法，输入以下代码：
  

1. ltemplify.py -name FILENAME FILENAME.data > FILENAME.lt

复制代码

六、【整理分享】使用moltemplate生成复杂体系的data文件

• 一方面，可以参阅互联网上已有的教程，例如LAMMPS复杂体系分子模拟packmol、moltemplate建模方法 - 知乎；LAMMPS之Moltemplate建模教程 - 哔哩哔哩；
• 另一方面，可以参阅moltemplate的官方教程。8.6.4. Moltemplate Tutorial — LAMMPS documentation
  

七、【问题求助】关于体系data文件的生成
需求：我需要建立起一个异质结，第一层为NiO晶体，第二层为20 A的有机小分子，要求有机小分子随机分布在其中。随后，我需要对体系进行分子动力学模拟，在400K下初始化，随后降温到298K，最后在298K下平衡10 ns。
我的总体方案：采用各种方案生成NiO单胞和小分子层的data文件，随后使用ltemplify.py将其转化为lt文件，最后使用moltemplate将其组合。
存在的问题：
1. 我希望有机小分子层使用amber-gaff力场；我希望有机小分子层使用packmol进行随机分散。但问题是，amber不支持对于多个分子的力场构建。其只允许我输入一个分子，并为其生成data文件，否则就会报错，我该如何对多个分子进行生成力场？
   对于该问题，我尝试了使用amber生成一个分子力场文件，随后使用moltemplate进行旋转和移动，但是该方法无法保证分子间不冲突、也无法保证其能够准确平移到我的限制范围内（move移动质心，无法保证整个分子都在合理的区域）。
2. NiO晶胞的结构较为简单，我希望使用interface力场。但是由于MS和msi2lmp对力场支持受限，我无法使用msi2lmp。我通过VMD得到了其data文件，该文件包含了连接关系，但是没有对连接关系的力场参数（键参数、角参数等）进行描述。我该如何做，可以将力场参数自动添加至data文件？
我的lt文件：

1. import "sam.lt"
   
2. import "NIO.lt"
   
3. 
   
4. write_once("Data Boundary") {
   
5.    0.0 35.0 xlo xhi
   
6.    0.0 38.0 ylo yhi
   
7.    0.0 52.50 zlo zhi
   
8. }
   
9. 
   
10. # --- SAM Layer ---
    
11. # 创建SAMs小分子随机分布的中间层
    
12. sam = new sam
    
13. sam.move(0.0,0.0,12.60)
    
14. nio = new NIO [8].move(4.20,0.00,0.00)
    
15.               [9].move(0.00,4.20,0.00)
    
16.               [3].move(0.00,0.00,4.20)
    
17. 
    
18. 
    
19. # --- Assemble the Complete System ---
    
20. # system {
    
21. #     subunit nio = new NIO_layer
    
22. #     subunit sam = new SAM_mol
    
23. # }

复制代码

sobereva

游客，本帖隐藏的内容需要积分高于 25 才可浏览，您当前积分为 0

sobereva

用acpype产生拓扑文件是早已过时的做法，用sobtop产生拓扑文件（http://sobereva.com/soft/Sobtop）灵活、强大、普适、方便得多，还不需要安装python环境、不需要装臃肿的Ambertools，在Windows下直接双击就能运行。这使得acpype产生AMBER/GAFF力场的拓扑文件方面已经没有使用价值了。

PLwang

sobereva 发表于 2025-1-22 06:19
用acpype产生拓扑文件是早已过时的做法，用sobtop产生拓扑文件（http://sobereva.com/soft/Sobtop）灵活、 ...

sob老师，我其实前期已经注意到您开发的这个非常高效的工具了。但是我注意到，该工具目前主要是应用在gromacs领域，没有直接输出data或其他lammps文件的选项。我也注意到您直接提到的，gromacs可以实现lammps的绝大多数功能，甚至比lammps做的更好。但是我们还是准备继续使用lammps完成这个模拟。目前我的打算是，准备好两个lt文件，随后通过move和rot，得到一个差不多的结构（我们表面的分子数量比较少，所以人工调整也是完全ok的）。这个工作结束之后，我准备参加一下gromacs的培训班，详细学一下gromacs。非常非常感谢

sobereva

PLwang 发表于 2025-1-22 15:59
sob老师，我其实前期已经注意到您开发的这个非常高效的工具了。但是我注意到，该工具目前主要是应用在gro ...

用以下程序都可以把gromacs的拓扑文件转成lammps的
GRO2LAM：https://github.com/hernanchavezthielemann/GRO2LAM
InterMol：https://github.com/shirtsgroup/InterMol

PLwang

补充一下最近遇到的计算问题和解决方案
Bond atom missing\Atom missing
对于该问题，其是由于某个原子在某一步内得到了非常高的速度，从而超过了正常键的范围\直接超出了正常的边界范围。
该问题网络上有一些解决方案，比如：1. 调整边界大小；2. 提高邻居搜索范围；3.提高力的范围；
但是我尝试了上述的多种方法，均没有解决问题。其往往仅能延后问题报错的步数，但无法解决该问题.

【给其他人（主要是新手）的建议】在遇到上述问题时，请一定要优先检查自己的建模是否正确，力场参数是否正确设置。我在建模中出现的问题主要有以下两个：
1. MS FORCITE模块没有自动计算出准确的结果。请千万不要过于相信MS的模块自动分配的力场类型。哪怕是非常简单的模型也必须人工进行检查，请确保所有的元素都被准确的分配了力场类型。
解决方案：手动更正力场类型。可以直接使用select功能选中错误分配的原子，并重新手动分配力场类型。
2. 请检查是否正确设置了角参数。我发现我使用msi2lmp生成的data文件中，部分角参数是缺失的，其角度直接被设置成了默认的0度。比如，A-B-C三个原子，其角度原应是180度，但如果被错误的设置为0度之后，其会导致AC两个原子距离过近。特别是当AB和BC键长相同时，其会导致两个原子直接重合。最终导致在初始时，两个原子获得非常高的受力，最终直接断裂对应的化学键；

ddddnight

请问根据mol2文件或者pdb文件是不是可以转换成lt文件啊，或者通过packmol或者ms的ac盒子建模后怎么转换成lt文件啊