计算化学公社

标题: TRAPPE_small力场包 [打印本页]

作者
Author: 喵星大佬 时间: 2023-4-5 16:03
标题: TRAPPE_small力场包
本帖最后由喵星大佬于 2024-1-12 18:25 编辑

==========2023-04-08更新===========
由于有用户反映该力场常与水模型合用以模拟水溶液的情况，因此更新了力场中的设置并添加了水模型的相关支持

对于默认使用SETTLE的硫化氢和环氧乙烷分子，添加了使用SHAKE约束的参数

如果需要模拟这两种分子在水溶液中的情况并对水使用SETTLE约束，可以在mdp文件中设置

define = -DNOSETTLES

复制代码

选项激活相关设置以保证可以正确对水使用SETTLE

另一项更新是在二氧化碳中使用了虚拟质量点。
由于二氧化碳为直线分子，如果直接使用SHAKE/RATTLE算法对3个键长进行约束将，180度的键角会导致约束步的Jacobian矩阵奇异而无法施加约束。因此，为了保证直线约束的稳定性，重新调整了约束方式。采用新参数后，可以正确的对二氧化碳分子施加LINCS约束
该参数已通过二氧化碳水溶液NPT模拟中的RDF进行了验证

此外，由于SHAKE约束方法无法在能量极小化过程中使用，因此建议如果需要做能量极小化，可以使用布朗动力学进行近似能量极小化

============分割线==============

TRAPPE力场是明尼苏达大学Siepmann课题组开发的一系列力场，主要用于包括相平衡在内的一系列模拟
其中的Small版本包含了7种简单分子，为硫化氢，氮气，氧气，环氧乙烷，氨气，二氧化碳和氦气
其官网为： http://trappe.oit.umn.edu/#small

本人根据官网和原文中提供的数据(不包含氦气，不管原文还是官网都没见到参数，但是页面上又说有)，将其转换为Gromacs的格式力场包，可进行相关模拟

力场包中包含了所有6种分子的itp文件和pdb文件

6种气体单独存在的情况下已经过测试(测试条件为常温常压，盒子大小约30 nm，内含1000个各种分子)无误，注意其建议范德华截断半径为1.4 nm，且须结合PME计算长程静电作用

由于其均为刚性模型，考虑到约束的运行效率，对硫化氢和环氧乙烷设置了SETTLE约束，因此不能同时加入这两种分子，当然由于本来这两种分子之间就会发生反应，不应该共存，所以实际并没有影响。但是注意，此时不能和任意刚性水模型合用，除非修改水模型的约束，因为Gromacs中的水模型默认用SETTLE的参数，但是一个模拟中只能有一个分子使用SETTLE

对其他的默认为LINCS约束，对氨气和二氧化碳由于涉及非对角元，LINCS算法的近似矩阵求逆不成立，只能使用SHAKE约束迭代求解验证中CO2较同样原子数的H2S慢了相当多(约20-30%)

在forcefield.doc中对各个参数的来源作了说明

欢迎使用

[attach]66605[/attach]

作者
Author: ByTheWay 时间: 2023-4-5 16:20
谢谢老师，帮大忙了

作者
Author: KRAS 时间: 2023-4-5 20:46

谢谢老师

作者
Author: 喵星大佬 时间: 2023-4-8 00:44
更新了版本，欢迎大家测试和使用

作者
Author: 喵星大佬 时间: 2023-4-21 13:49
之前CO2的键长设置错误，现已修正

作者
Author: xiaoruoqiu 时间: 2023-11-2 06:20
请问这些力场参数能否移植到诸如oplsaa、amber、charmm力场里？谢谢

作者
Author: 喵星大佬 时间: 2023-11-2 15:12
应该不可以吧

如果需要其他分子可以先看看trappe的网站有没有

作者
Author: 喵星大佬 时间: 2023-11-18 16:45
本帖最后由喵星大佬于 2023-11-18 16:47 编辑

有同学反应在刚结束的gromacs培训班文件包中trappe_small的力场文件有误，经检查为4月5日的初版，各位同学如需使用请在本帖中下载最新版，在新版中已对错误进行修正并更改约束形式以使动力学稳定性增加，修改内容见本帖一楼

感谢卢老师分享并同时提醒更新一下文件，谢谢

作者
Author: KRAS 时间: 2023-11-20 22:59

喵星大佬发表于 2023-11-18 16:45
有同学反应在刚结束的gromacs培训班文件包中trappe_small的力场文件有误，经检查为4月5日的初版，各位同学 ...

感谢楼主大大在力场文件中留了邮箱

作者
Author: lmch 时间: 2024-1-10 21:54
请问
1. CO2分子能否使用LINCS约束?
2. 拓扑文件中的虚拟位点参数是如何确定的？
3. 能够推广使用到其他形式类似的co2力场参数吗？

作者
Author: 喵星大佬 时间: 2024-1-11 03:00

lmch 发表于 2024-1-10 21:54
请问
1. CO2分子能否使用LINCS约束?
2. 拓扑文件中的虚拟位点参数是如何确定的？

可以使用lincs约束
仔细看文件，实际上这个CO2是由两个质量22的实际原子构成的，距离通过使其转动惯量与原力场定义的CO2相等来设置的，也就是说实际上有运动方程的只是一个双原子分子，而且这两个假原子不直接受到任何力
3个在力场设置的原子的位置是虚拟作用位点，其通过两个假原子的坐标线性组合来确定位置，其受力会通过线性组合给到两个假原子上(计算上就是链式求导)。
这种方式可以扩展到任何刚体线性分子上，因为实际上这个刚体3原子分子只有5个自由度(有3个约束条件和两个简并自由度)，而这和刚体双原子分子是没有区别的(本质上就是3个质心位置坐标和两个欧拉角)。因此这个方法完全可以扩展，保证所有原始原子(相互作用位点)的力矩，分子总质量，分子转动惯量与原始分子一样即可，所以原则上任线性刚体分子都可以转换成无相互作用双原子分子+虚拟相互作用位点的形式

作者
Author: lmch 时间: 2024-1-11 15:45

喵星大佬发表于 2024-1-11 03:00
可以使用lincs约束
仔细看文件，实际上这个CO2是由两个质量22的实际原子构成的，距离通过使其转动惯量与 ...

感谢您的解答
还有个问题哈，您在文件包中的这些设置是否是正确的？

ffnonbonded.itp中
CC 6 0.0000 0.0000 A 0.280000000 0.224478000 ;C as virtual site in CO2

元素质量也是0吗？

对其他的默认为LINCS约束，对氨气和二氧化碳由于涉及非对角元，LINCS算法的近似矩阵求逆不成立，只能使用SHAKE约束迭代求解验证中CO2较同样原子数的H2S慢了相当多(约20-30%)

这里的意思是指？

作者
Author: 喵星大佬 时间: 2024-1-11 21:59
本帖最后由喵星大佬于 2024-1-12 00:21 编辑

lmch 发表于 2024-1-11 15:45
感谢您的解答
还有个问题哈，您在文件包中的这些设置是否是正确的？

当然是0，真正在原子位置的是虚拟点，只是用来提供真实相互作用的力矩。真正用来描述质心运动和分子转动的是那两个22质量啥相互作用都没有的原子。lincs约束的也是这个玩意的运动

对于CO2那句已经作废了，当时用SHAKE算法约束三点共线，后来发现在求解约束方程的时候矩阵不是绝对对角占优且不保证正定，此时SHAKE方法采用的高斯-赛德尔迭代不成立，产生的结果都是异常的，结构会变得很奇怪，因此改用了后来的方法

作者
Author: lmch 时间: 2024-1-12 17:29

喵星大佬发表于 2024-1-11 21:59
当然是0，真正在原子位置的是虚拟点，只是用来提供真实相互作用的力矩。真正用来描述质心运动和分子转动 ...

OC 8 16.00 0.0000 A 0.305000000 0.656806000 ;O as virtual site in CO2

那请问为何这里不为0呢？

作者
Author: 喵星大佬 时间: 2024-1-12 18:07
标题: g
本帖最后由喵星大佬于 2024-1-12 18:13 编辑

lmch 发表于 2024-1-12 17:29
那请问为何这里不为0呢？

可能当时写错了。。。。。我查了下我这里的是对的，我重新传一下

不过无所谓其实，因为在CO2的itp里面重新写了一遍质量，是0，所以在[ atomtypes ]里面写的质量会被覆盖掉

作者
Author: 3163582842 时间: 2024-3-13 20:56
请问大佬，如何写packmol能识别的trappe-small氮气分子的.pdb文件呢？

作者
Author: 喵星大佬 时间: 2024-3-14 02:19

3163582842 发表于 2024-3-13 20:56
请问大佬，如何写packmol能识别的trappe-small氮气分子的.pdb文件呢？

文件包里不是已经提供了各种分子的pdb格式文件么

作者
Author: 3163582842 时间: 2024-3-14 17:25

喵星大佬发表于 2024-3-14 02:19
文件包里不是已经提供了各种分子的pdb格式文件么

大佬，这个pdb格式，packmol可以读取吗？后面想用packmol+moltemplate建模

作者
Author: 喵星大佬 时间: 2024-3-14 17:53

3163582842 发表于 2024-3-14 17:25
大佬，这个pdb格式，packmol可以读取吗？后面想用packmol+moltemplate建模

标准的pdb格式为啥读取不了

作者
Author: 3163582842 时间: 2024-3-14 21:52

喵星大佬发表于 2024-3-14 17:53
标准的pdb格式为啥读取不了

好的感谢大佬我试试去

作者
Author: raku 时间: 2024-11-28 21:36
专业专业

作者
Author: raku 时间: 2024-12-4 22:18
请教下大佬，在N2的itp中是否需要设置两个原子之间的[ exclusions ]，盼复，谢谢啦！

作者
Author: 喵星大佬 时间: 2024-12-5 14:56

raku 发表于 2024-12-4 22:18
请教下大佬，在N2的itp中是否需要设置两个原子之间的[ exclusions ]，盼复，谢谢啦！

函数类型1的constraint本来就排除了，而且我应该写了exclusions用来排除和虚拟点的

作者
Author: zhishidishu 时间: 2024-12-25 10:31
请问下大佬，我想要模拟液氨的密度，设置盒子大小是4.5nm^3，采用npt控压1MPa，首先采用GAFF力场，计算液氨密度为850kg/m3，实验密度应为（600kg/m3），采用LINCS，计算速度大概1000ns/day，更改采用trappe-small力场，更改算法为SHAKE，液氨密度为580kg/m3接近实验值，但速度下降到大概50ns/day。
1.更改算法后这样速度正常吗？
2.如果正常，更改哪些设置可以使速度变快？
3.如果在液氨中加入有机分子，能否对其他分子设置lincs算法，单独对NH3采用SHAKE算法？

作者
Author: 喵星大佬 时间: 2024-12-25 12:25
shake要迭代肯定慢，没法混用不同的约束方法，但可以适当降低shake算法的误差容限

作者
Author: zhishidishu 时间: 2024-12-25 15:11

喵星大佬发表于 2024-12-25 12:25
shake要迭代肯定慢，没法混用不同的约束方法，但可以适当降低shake算法的误差容限

感谢老师的回复，我还想再问下
1.手册上说比LINCS稍微慢一点，我这个慢了20倍，是不是太慢了？
2.还有老师您说的是更改shake-tol吗，我尝试从0.0001增大到了0.01，速度也仅从50ns/day增加到70ns/day，还有其他加速的方法吗？

作者
Author: 喵星大佬 时间: 2024-12-25 15:23

zhishidishu 发表于 2024-12-25 15:11
感谢老师的回复，我还想再问下
1.手册上说比LINCS稍微慢一点，我这个慢了20倍，是不是太慢了？
2.还有 ...

稍慢指的是只约束孤立键长，你这个要约束的是6个耦合的自由度，计算量本来就很大。

你可以用类似CO2的方法把氨变成3个等效质点，保证3个转动自由度的转动惯量和总质量与原分子一致，所有原子的位置放虚拟点，用SATTLE算法约束这个等边三角形

作者
Author: zhishidishu 时间: 2024-12-25 17:41

喵星大佬发表于 2024-12-25 15:23
稍慢指的是只约束孤立键长，你这个要约束的是6个耦合的自由度，计算量本来就很大。

你可以用类似CO2的 ...

好的十分感谢老师的回复，刚接触gromacs对虚拟点相关知识不太了解，手册上看的也是一头雾水，网上有没有相关计算过程，老师能否推荐学习一下，再次感谢老师！

作者
Author: thor 时间: 2025-2-3 20:12

喵星大佬发表于 2024-1-11 21:59
当然是0，真正在原子位置的是虚拟点，只是用来提供真实相互作用的力矩。真正用来描述质心运动和分子转动 ...

这个结果异常是指什么呢，我看到很多文章（气体分离）模拟CO2用的TraPPE力场，没有提到增加额外的虚拟质量点。

作者
Author: raku 时间: 2025-3-20 19:38
请教老师，在使用您制作的itp和pdb模拟NH3时，
总是提示SHAKE is not supported with domain decomposition and constraints that cross
domain boundaries, use LINCS
-------------------------------
mdp设置如下：
constraints             =  none / all-bonds
constraint-algorithm    = SHAKE
continuation          = no
Shake-SOR             = yes
shake-tol             = 0.0001
-------md 命令如下-----------------
gmx mdrun -v -s be.tpr -deffnm eq -ntmpi 1 -nt 1 ；使用-ntmpi 1 和 -nt 1防止域分解
---------------------------------------
不知道这样设置还有什么不对的地方导致提示报错提示。请帮指导，谢谢！

作者
Author: 喵星大佬 时间: 2025-3-20 22:44

thor 发表于 2025-2-3 20:12
这个结果异常是指什么呢，我看到很多文章（气体分离）模拟CO2用的TraPPE力场，没有提到增加额外的虚拟质 ...

那个只是为了能顺利模拟进行的修改，跟原版是等价的，实际上并没有可以约束三点共线并兼容热浴的方法，所以不这么改应该是不能跑起来的，模拟完了作图不展示虚拟点也是很正常的

作者
Author: 喵星大佬 时间: 2025-3-20 22:44

raku 发表于 2025-3-20 19:38
请教老师，在使用您制作的itp和pdb模拟NH3时，
总是提示SHAKE is not supported with domain decompositio ...

我并没有碰到过这种问题，你是不是和其他分子一起用了

作者
Author: raku 时间: 2025-3-21 13:09

喵星大佬发表于 2025-3-20 22:44
我并没有碰到过这种问题，你是不是和其他分子一起用了

没有的，纯NH3就这样提示

作者
Author: raku 时间: 2025-3-21 13:16

喵星大佬发表于 2025-3-20 22:44
我并没有碰到过这种问题，你是不是和其他分子一起用了

找到问题了：Apple M1芯片不能这样。使用Windows电脑或者Linux就好了

欢迎光临计算化学公社 (http://bbs.keinsci.com/)