计算化学公社

标题: 利用Gromacs计算有机物混合物沸点的一些疑问 [打印本页]

作者
Author: williamwyx 时间: 2023-9-15 22:03
标题: 利用Gromacs计算有机物混合物沸点的一些疑问
最近在尝试利用GROMACS对一些有机混合物沸点进行计算，主要流程是对体系进行能量最小化、nvt模拟然后npt模拟读密度，根据不同温度下重复进行npt模拟后得到的密度来判断沸点，有一些疑问如下：
（1）这个方法是否合理，还是需要每一次修改温度后完整再从能量最小化开始？
（2）npt模拟是不是密度的平均值几乎不变了就可以认为是平衡了？
（3）对于一些沸点低的纯有机物，发现计算的密度值只有10左右，有的时候温度也控制不住，但这个温度下有机物应该明显是液体的，对于沸点高的有机物模拟结果就是正常的，请问这是什么原因？结构都是通过LigParGen直接生成的，选择的OPLSAA力场

作者
Author: sobereva 时间: 2023-9-16 01:02
1 em做一次就够了
更高温度的MD可以基于之前更低温度、发现没达到沸点的MD最终状态接着跑

2 我没见过基于考察密度来算沸点的文章。超过沸点后达到平衡不太可能，因为达到沸点后体系会一直气化、盒子尺寸会不断迅速增大。用真空-液态两相的NVT模型，考察什么温度下液态分子会充分气化并充满整个盒子，是更可行的。

3 可能当前力场对当前物质的沸点描述得本来就不准，或者计算有什么硬伤。即实验上没超过沸点，但此温度已经超过当前计算模型下的沸点了。

作者
Author: williamwyx 时间: 2023-9-16 10:47

sobereva 发表于 2023-9-16 01:02
1 em做一次就够了
更高温度的MD可以基于之前更低温度、发现没达到沸点的MD最终状态接着跑

谢谢sob老师的回复，还有一些进一步的疑问：
（1）之前选择密度一方面是因为觉得气化后体系应该会剧烈波动，密度值会变小，另外一方面是因为想用npt来控制压力，如果使用nvt模拟的话是不是就是先构建液相盒子，再把盒子某个方向扩大一倍就可以实现真空-液态两相了？请问用nvt怎么控制压力呢，害怕沸点对于的压力不是一开始想要的压力了。
（2）用nvt计算沸点是不是就是每次从em.gro更换温度直接跑nvt就可以？请问用哪个指标来看分子是否气化充满盒子呢？
（3）就是不考虑沸点的话，同一个力场参数和模拟流程对小分子模拟的不准是不是说明力场文件有问题？因为我试了一下碳数多一点的有机物，npt给出的密度计算结果误差是完全可以接受的。

作者
Author: sobereva 时间: 2023-9-18 08:58

williamwyx 发表于 2023-9-16 10:47
谢谢sob老师的回复，还有一些进一步的疑问：
（1）之前选择密度一方面是因为觉得气化后体系应该会剧烈波 ...

1 气化之后靠压浴很难理想地控压
NVT本来也不需要你控制压力，又不是计算特定压力下的属性，只需要根据轨迹判断是否已经汽化了就完了

2 是。统计不同截面密度分布曲线，或者直接肉眼看轨迹

3 大概率是

作者
Author: williamwyx 时间: 2023-9-18 11:18

sobereva 发表于 2023-9-18 08:58
1 气化之后靠压浴很难理想地控压
NVT本来也不需要你控制压力，又不是计算特定压力下的属性，只需要根据 ...

谢谢sob老师的解答
（1）我想的是nvt的话是不是需要设置一个比较合适的尺寸，因为其实是想模拟常压下的沸腾，计算1bar下的沸点。肉眼看轨迹是不是就是把gro文件放到vmd里就可以，看分子是不是充满了盒子？请问这个统计密度分布可以直接用gmx的代码输出嘛？
（2）这两天模拟发现好像npt模拟在切换温度时因为要选择bren控制压力，好像最后平均的压力也很难控制在1bar左右，平均值稳定在8-9bar左右，然后密度和提交都稳定了，这样是不是意味着压力其实被控制在了一个更高的值所以导致没有汽化，这种有办法精准的控制压力吗？

作者
Author: sobereva 时间: 2023-9-18 11:23

williamwyx 发表于 2023-9-18 11:18
谢谢sob老师的解答
（1）我想的是nvt的话是不是需要设置一个比较合适的尺寸，因为其实是想模拟常压下的 ...

1 是。看整个轨迹更好，有助于判断气化过程
gmx density也可以对特定时间范围内的轨迹计算平均密度分布曲线用于判断

2 没有叫bren的东西

作者
Author: williamwyx 时间: 2023-9-18 13:07

sobereva 发表于 2023-9-18 11:23
1 是。看整个轨迹更好，有助于判断气化过程
gmx density也可以对特定时间范围内的轨迹计算平均密度分布 ...

谢谢sob老师
就是npt调整压力现在用的是berendsen控压，也是因为这个过程应该是不稳定的，然后发现压力好像很难压到1bar，是不是还是应该平稳后在转一下PR控压，观察一下体积会不会快速增加

作者
Author: sobereva 时间: 2023-9-18 18:42

williamwyx 发表于 2023-9-18 13:07
谢谢sob老师
就是npt调整压力现在用的是berendsen控压，也是因为这个过程应该是不稳定的，然后发现压力 ...

气体的可压缩系数远大于凝聚相，一般的压浴很难控压到稳定，所以不让你用NPT来跑

作者
Author: williamwyx 时间: 2023-9-20 22:49

sobereva 发表于 2023-9-18 18:42
气体的可压缩系数远大于凝聚相，一般的压浴很难控压到稳定，所以不让你用NPT来跑

sob老师好，我这两天用packmol和vmd构建了液相-真空两部分的盒子，但发现能量最小化后，1ns左右的nvt模拟在不同温度下调用gmx density得到的相同的结果，观察gro文件好像也是相同的，请问这是说明我应该增加模拟时长嘛，还是参数设置的有问题？然后好像gro文件显示盒子里有分子，但density对于的洁面显示的是0，这是不是说明计算的不对？

作者
Author: sobereva 时间: 2023-9-21 01:36

williamwyx 发表于 2023-9-20 22:49
sob老师好，我这两天用packmol和vmd构建了液相-真空两部分的盒子，但发现能量最小化后，1ns左右的nvt模拟 ...

自己看一下轨迹自然就清楚了

作者
Author: williamwyx 时间: 2023-9-23 19:41

sobereva 发表于 2023-9-21 01:36
自己看一下轨迹自然就清楚了

sob老师好，之前发现是读取density分布的文件用错了，现在解决了，但有了两个新的疑惑：
（1）目前是把一开始的分子聚集在盒子中间，两侧z方向上下扩大盒子大小产生真空区，然后进行nvt模拟，但发现在低温（例如室温下）分子依然会上下扩散，请问这个一般是什么原因，因为我原来可能认为维持液相的话分子会聚集在中间不动
（2）如果扩散不对的话是不是就需要修改模拟参数？我现在用的是gmx软件自带的opls力场参数，如果修改的话应该怎么修改呢？
谢谢sob老师

作者
Author: sobereva 时间: 2023-9-24 03:46

williamwyx 发表于 2023-9-23 19:41
sob老师好，之前发现是读取density分布的文件用错了，现在解决了，但有了两个新的疑惑：
（1）目前是把 ...

正常情况。到了沸点以上分子会迅速大量逸出到真空区，而即便在沸点以下，MD过程中照样会时不时蒸发出个别分子

作者
Author: williamwyx 时间: 2023-9-24 12:58
本帖最后由 williamwyx 于 2023-9-24 13:00 编辑

sobereva 发表于 2023-9-24 03:46
正常情况。到了沸点以上分子会迅速大量逸出到真空区，而即便在沸点以下，MD过程中照样会时不时蒸发出个别 ...

sob老师好，我模拟nvt模拟发现在室温和高温下分子都会扩散的很开，还有些团簇，这是什么原因呢，图片是我设置的初始和中间帧的盒子分子分布，是不是盒子xy方向太窄导致分子穿盒子和对面的产生相互作用？

作者
Author: sobereva 时间: 2023-9-24 13:28

williamwyx 发表于 2023-9-24 12:58
sob老师好，我模拟nvt模拟发现在室温和高温下分子都会扩散的很开，还有些团簇，这是什么原因呢，图 ...

跟窄不窄没关系
大概率是LJ参数或者原子电荷不合理，或者mdp、拓扑文件有硬伤等

作者
Author: williamwyx 时间: 2023-9-24 20:05

sobereva 发表于 2023-9-24 13:28
跟窄不窄没关系
大概率是LJ参数或者原子电荷不合理，或者mdp、拓扑文件有硬伤等

谢谢sob老师，我目前用的是ligpargen直接产生的结构和力场文件，应该这些原子类型不需要调用原始gromacs里的力场参数，唯一做的操作可能就是手动修改了一下原子电荷使其和为0，不知道是不是这个的原因？然后mdp文件就使用的正常的em和nvt的mdp文件，用的v-res控制温度

作者
Author: williamwyx 时间: 2023-9-25 16:29
sob老师好，请问ligpargen直接生成的itp文件中哪部分体现了LJ参数呢？因为我现在体系只有一种纯物质，应该拓扑文件、原子电荷都是合理的，可能问题就出在mdp模拟参数上，下面是我使用的mdp文件，请sob老师帮我看看有啥问题，谢谢。

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