计算化学公社

标题: MSD计算离子电导率的问题 [打印本页]
作者
Author: jlyjlysjd    时间: 2024-6-12 14:30
标题: MSD计算离子电导率的问题
在使用gromacs计算离子电导率时遇到了一些问题：
根据Electrical Conductivity based on Diffusion Coefficients (aqion.de) 里的Nernst-Einstein公式可得：
(, 下载次数 Times of downloads: 67)
我想计算1M MgCl2水溶液的Mg2+的离子电导率，通过gmx msd输出的文件可以得到Di=0.4184*10^-5 ，
带入进去后确实可以计算得出一个数值，但是计算出来感觉比较大为62.83，类比实验获得的1M ZnCl2的电导率约为25，感觉不太正常。
于是我进行公式的单位推导，但推导一半后，无法继续推导，如下所示：
(, 下载次数 Times of downloads: 66)
请各位看看问题是在哪里呢？这样的计算方法是对的吗？


感谢各位！
作者
Author: 二分音符    时间: 2024-6-12 14:42
本帖最后由 二分音符 于 2024-7-4 12:44 编辑

对于浓溶液，NE方程不再适用，应该用电流自相关计算电导率。

这是因为，浓溶液中存在束缚阴阳离子对，他们的扩散不会贡献电导率，只有自由离子的扩散会贡献电导率。所以MSD计算出来的结果会高估电导率。

详细原理可以看这篇文章
作者
Author: jlyjlysjd    时间: 2024-6-13 10:34
本帖最后由 jlyjlysjd 于 2024-6-13 10:35 编辑
二分音符 发表于 2024-6-12 14:42
对于浓溶液，NE方程不再适用，应该用速度自相关计算电导率。

这是因为，浓溶液中存在束缚阴阳离子对，他 ...

好的，谢谢。我才上手gromacs不久，那请问是不是使用gmx velacc以及analyze分析平衡后的prd.trr文件，得到扩散系数呢
作者
Author: jlyjlysjd    时间: 2024-6-16 23:40
二分音符 发表于 2024-6-12 14:42
对于浓溶液，NE方程不再适用，应该用速度自相关计算电导率。

这是因为，浓溶液中存在束缚阴阳离子对，他 ...

您好老师，就是我现在使用速度自相关了，但是有一个问题就是结合文章和您所说的，那怎么确定自由离子的数量呢，我是计算1M MgCl2水溶液，动力学下来确实有六合水镁离子，还有[MgClH2O]+离子，那自由离子是不是可以理解为没有与阴离子结合的六合水镁离子呢，gromacs里面能不能统计出来不同水合离子的数量能用index索引呢
作者
Author: 二分音符    时间: 2024-6-18 12:32
本帖最后由 二分音符 于 2024-7-4 12:44 编辑
jlyjlysjd 发表于 2024-6-16 23:40
您好老师，就是我现在使用速度自相关了，但是有一个问题就是结合文章和您所说的，那怎么确定自由离子的数 ...

电流自相关可以直接算出准确的电导率，无需再经过msd和自由离子比例，除非你想研究相关性质。

如果算出来依然低估，那么说明力场不够好，这个时候一般缩放电荷，比如所有离子的电荷＊0.8。这一步需要做一系列值直到电导率和实验一致。需要注意的是，缩放电荷后可能其他性质不准了，这是经典力场的局限，很难兼顾。如果侧重结构准，就用之前的电荷，如果侧重动力学，就用缩放的电荷。

上段都是tricks，足够发文章，但存在显而易见的问题。追求更准的描述需要更高级的模拟工具，根据自己的资源和水平来了。
作者
Author: jlyjlysjd    时间: 2024-6-18 14:47
二分音符 发表于 2024-6-18 12:32
速度自相关可以直接算出准确的电导率，无需再经过msd和自由离子比例，除非你想研究相关性质。

如果算 ...

请问老师，通过速度自相关函数怎么算出电导率呢
作者
Author: 二分音符    时间: 2024-6-18 17:59
jlyjlysjd 发表于 2024-6-18 14:47
请问老师，通过速度自相关函数怎么算出电导率呢

抱歉之前说错了，要用gmx current算电流自相关，后面公式原文里有
作者
Author: jlyjlysjd    时间: 2024-6-18 18:16
二分音符 发表于 2024-6-18 17:59
抱歉之前说错了，要用gmx current算电流自相关，后面公式原文里有

十分感谢老师，我后面也查到了使用gmx current 。使用以后结果输出为：
Mg离子的电导率为：
Einstein-Helfand fit to the MSD of the translational dipole moment yields:

sigma=11.3085
translational fraction of M^2: 4.2404
Dielectric constant using EH: 65.5655

11.3085 mS/cm 结果就十分正常了。
谢谢老师!
作者
Author: jlyjlysjd    时间: 2024-7-3 19:31
本帖最后由 jlyjlysjd 于 2024-7-3 22:22 编辑
二分音符 发表于 2024-6-18 12:32
速度自相关可以直接算出准确的电导率，无需再经过msd和自由离子比例，除非你想研究相关性质。

如果算 ...

老师您好，我最近又产生了一些问题：我使用gmx current 输出了不同浓度的ZnCl2溶液的Zn2+的电导率，但是我发现并没有产生这样的规律——随着浓度提高电导率提高。输出的电导率时大时小，没有规律性。
所以老师，请问：
1、电导率是看整个系统的还是可以看Zn离子的呢？
2、gmx current 这个后面的参数需不需要调整呢，比如-bfit(default:100)和-efit(default:400) 我尝试调整了这两个的值，发现sigma会有很大的变动，应该怎么确定呢，另外我有看了一些有关电荷缩放的文献，我有一个问题就是：
3、是不是在itp文件里面将每一个离子的电荷都缩小0.8倍就可以了呢，其他的参数是不是不用调节？
谢谢老师！

作者
Author: slxc920113    时间: 2024-7-4 10:01
二分音符 发表于 2024-6-12 14:42
对于浓溶液，NE方程不再适用，应该用速度自相关计算电导率。

这是因为，浓溶液中存在束缚阴阳离子对，他 ...

但是对于电化学反应来说，其实要的是离子的传质，是否是阴阳离子对其实不重要。不同的溶剂化结构只是影响去溶剂化过程的难易而已。这个时候离子扩散的通量反而是有真实的物理意义的。
作者
Author: 二分音符    时间: 2024-7-4 13:13
jlyjlysjd 发表于 2024-7-3 19:31
老师您好，我最近又产生了一些问题：我使用gmx current 输出了不同浓度的ZnCl2溶液的Zn2+的电导率，但是 ...

1. 可以看系统和每个离子的贡献
2. 最好不要用gmx直接输出的值，而是-o的xvg里取平衡段做拟合，参考文献里的公式
3. 是的

另外需要注意，这个值确实会波动很大，所以要多跑几个independent cases做平均，才能确保算出了准确的值。
作者
Author: 二分音符    时间: 2024-7-4 13:15
slxc920113 发表于 2024-7-4 10:01
但是对于电化学反应来说，其实要的是离子的传质，是否是阴阳离子对其实不重要。不同的溶剂化结构只是影响 ...

扩散和电导是两个独立的概念，都有对应的物理意义。只是浓溶液里不能简单地把扩散系数换算成电导率。
作者
Author: jlyjlysjd    时间: 2024-7-4 16:42
本帖最后由 jlyjlysjd 于 2024-7-4 16:54 编辑
二分音符 发表于 2024-7-4 13:13
1. 可以看系统和每个离子的贡献
2. 最好不要用gmx直接输出的值，而是-o的xvg里取平衡段做拟合，参考文献 ...

非常感谢老师耐心的回答！
我还有一个问题，就是输出文件中有一个dsp.xvg文件，打开后看到这个文件就是文献中的公式：
(, 下载次数 Times of downloads: 66)
(, 下载次数 Times of downloads: 65)
而文献中说可以从均方偶极子位移随时间的斜率中获得电导率，所以正如老师您说的只需要拟合就行，但是我使用gmx current -s prd.tpr -f prd.trr  -dt 1 ，输出了全部的产生相时间，获得了一张均方偶极子位移和时间的图：
(, 下载次数 Times of downloads: 59)
但是没有表现出线性，所以请问老师，应该去拟合哪一段的呢。



作者
Author: 二分音符    时间: 2024-7-4 22:44
jlyjlysjd 发表于 2024-7-4 16:42
非常感谢老师耐心的回答！
我还有一个问题，就是输出文件中有一个dsp.xvg文件，打开后看到这个文件就是 ...

没有线性段是因为一个case误差太大，所以才需要跑5-10组independent cases然后取平均得到线性段，注意最好用nvt跑，因为这个量对体积非常敏感。
作者
Author: jlyjlysjd    时间: 2024-7-5 13:58
二分音符 发表于 2024-7-4 22:44
没有线性段是因为一个case误差太大，所以才需要跑5-10组independent cases然后取平均得到线性段，注意最 ...

明白了 谢谢老师！
作者
Author: 二分音符    时间: 2024-7-5 17:31
本帖最后由 二分音符 于 2024-7-6 16:54 编辑
jlyjlysjd 发表于 2024-7-5 13:58
明白了 谢谢老师！

最近自己动手做了下，完整步骤：
1. 计算5-10个NVT，速度输出频率做到0.1ps （这个要自己测试，case可能要20个，频率可能要0.01ps）
2. gmx current -caf caf.xvg, 这个输出的就是EACF，注意归一化，把第二列每个值都除以t=0时的值，确保t=0时ACF=1。多个体系平均后得到比较平滑的曲线，这样你就能做出文献中的FIG.6a
3. 利用公式E1对ACF进行积分，积分符号里面的内容就是ACF，由于ACF会趋近于0，所以这个积分也会收敛。得到的就是电导率。（注意用caf.xvg中的原始数据积分并带入公式E1，归一化只是为了画图用的，不归一化也能画图，注意带上xvg中的单位）
作者
Author: jlyjlysjd    时间: 2024-7-5 19:18
二分音符 发表于 2024-7-5 17:31
最近自己动手做了下，完整步骤：
1. 计算5-10个NVT，速度输出频率做到0.1ps
2. gmx current -caf caf.x ...

谢谢老师！
作者
Author: 二分音符    时间: 2024-7-7 16:47
这篇工作和你想做的比较像，可以参考
https://www.sciencedirect.com/sc ... i/S0167732216332718
作者
Author: jlyjlysjd    时间: 2024-7-7 19:25
本帖最后由 jlyjlysjd 于 2024-7-7 19:33 编辑
二分音符 发表于 2024-7-7 16:47
这篇工作和你想做的比较像，可以参考
https://www.sciencedirect.com/science/article/pii/S016773221633 ...

是的老师，这篇文章中电导率随浓度的变化确实是和我现在做的很类似。
老师，就是您能帮我看看我输出的caf.xvg是否是正确的呢，我根据您的方法绘制出来以后，发现与图6(a)不相符。
下面是prd.mdp，输出日志相关的部分：
integrator           = md        
dt                     = 0.002     ; 2 fs
nsteps               = 5000000    ; 10.0 ns

nstenergy          = 2000
nstlog               = 2000
nstxout             = 2000
nstvout                = 100
nstfout                = 2000

之后连续模拟了5次，使用命令：gmx current -s prd.tpr -f prd.trr   -caf caf.xvg,  分别输出了5个caf.xvg()。

使用归一化后，并使用origin里面的合并和计算多条曲线的均值，画出了下面的一幅图：
(, 下载次数 Times of downloads: 40) 这个与文献的太不相符了，主要时间好像不太对，我的太长了。
将x轴方法，发现一点类似的规律。
(, 下载次数 Times of downloads: 52)
所以老师，请问问题是出在哪里呢，是不是mdp中时间尺度出现的问题呢？根据文献“[color=rgba(0, 0, 0, 0.85)]MD 获得的 ECACF 在几皮秒内快速衰减到零”。是不是需要控制整体的模拟时间在3ps左右呢？或者是不是我的caf.xvg文件输出的有问题呢？
学生愚钝，请老师赐教，谢谢老师！


caf.xvg:
(, 下载次数 Times of downloads: 45)
(, 下载次数 Times of downloads: 49)
作者
Author: 二分音符    时间: 2024-7-8 11:52
最重要的一点是nstvout改成5，输出间隔必须非常小。展示Fig.6a时只需要显示几个ps的数据，但是积分的时候得做到1-2ns。
我的mdp文件给你参考下
作者
Author: jlyjlysjd    时间: 2024-7-8 19:50
二分音符 发表于 2024-7-8 11:52
最重要的一点是nstvout改成5，输出间隔必须非常小。展示Fig.6a时只需要显示几个ps的数据，但是积分的时候得 ...

感谢老师的文件，经过修改mdp后，我总计运行了10次NVT后得到10组caf数据，均一化后得到如下数据：
(, 下载次数 Times of downloads: 44)
放置于origin后，使用曲线求均值后得到了如下的图： (, 下载次数 Times of downloads: 51)
虽然后面确实是趋近于0的，但是是在0上下波动，并没有单调趋近，所以这个曲线不知道正不正确。
然后使用积分工具对0-2ns进行积分面积计算后，得到如下数据：
(, 下载次数 Times of downloads: 54)
面积是：0.06043
通过同时,积分面积/3VkbT=0.06043/1.31E10=4.61E08
请问老师，这个电导率的计算哪里出现问题了呢？
是不是还是曲线的问题呢？


老师可以帮我看看我的NPT和NVT系综的mdp是否有问题呢？
非常感谢老师！
作者
Author: 二分音符    时间: 2024-7-9 10:46
注意单位，全部换成SI制。
你的ACF看起来还可以，再画一个积分值随时间变化的曲线（t时刻的函数值为ACF从0积分到t的面积）。
mdp看着没问题
作者
Author: jlyjlysjd    时间: 2024-7-9 16:34
本帖最后由 jlyjlysjd 于 2024-7-9 16:52 编辑
二分音符 发表于 2024-7-9 10:46
注意单位，全部换成SI制。
你的ACF看起来还可以，再画一个积分值随时间变化的曲线（t时刻的函数值为ACF从0 ...

好的老师，下面是ACF(t)对t积分的图：
(, 下载次数 Times of downloads: 48)

单位问题我计算如下：
[attach]94928[/attach]
所以需要乘以10^21。
通过gro文件获得体积：15.7514632602629 nm^3，Kb=1.380649 * 10^-23 ,T=298.15
(, 下载次数 Times of downloads: 51)
得到了比较合理的数据。
谢谢老师！

作者
Author: jlyjlysjd    时间: 2024-7-11 14:15
二分音符 发表于 2024-7-9 10:46
注意单位，全部换成SI制。
你的ACF看起来还可以，再画一个积分值随时间变化的曲线（t时刻的函数值为ACF从0 ...

老师，我在算新的一个浓度的时候发现积分结果为负值，这说明结果是错误的吗
作者
Author: 二分音符    时间: 2024-7-11 16:52
jlyjlysjd 发表于 2024-7-11 14:15
老师，我在算新的一个浓度的时候发现积分结果为负值，这说明结果是错误的吗

你可以检查一下之前的结果，不同case之间算出的电导率差多少，然后画出errorbar（根据误差分析的知识），然后和文献里的errorbar比较一下。
如果波动极大，那么说明case还是少了，要么加case，要么加大模拟体系。要是还降不下来，那这个思路就走不下去了，得再换方法。
作者
Author: jlyjlysjd    时间: 2024-7-11 18:28
本帖最后由 jlyjlysjd 于 2024-7-11 19:08 编辑
二分音符 发表于 2024-7-11 16:52
你可以检查一下之前的结果，不同case之间算出的电导率差多少，然后画出errorbar（根据误差分析的知识）， ...

老师，我做的是0.5M MgCl2溶液的电导率，刚查阅了一下文献库，发现好像并没有这样类似的文献。
有一篇中文的文献，拟合了电导率：
(, 下载次数 Times of downloads: 48)
我做了10次电导率的积分errorbar:
(, 下载次数 Times of downloads: 51)
这样是不是可以说明波动是比较大的。

其实从10组电导率的数据来看，发现波动还是很大的：
(, 下载次数 Times of downloads: 45)
但基本都是呈现负值情况，感觉再重复多次，是不是大概率还是负值呢？

作者
Author: 二分音符    时间: 2024-7-11 20:35
这种情况可能不适合用电流自相关算电导率了，因为数值误差过大，可能还是只能用NE方程计算。
作者
Author: 二分音符    时间: 2024-7-16 20:08
我又重新试了下，把步长改为0.2fs，每步输出速度可以收敛。会有小部分case的caf文件std.dev列显著大于其他case或者是负值，把这些舍掉。
这样做会导致一个case占据约2T硬盘空间，只能算完就删。
作者
Author: jlyjlysjd    时间: 2024-7-21 12:43
二分音符 发表于 2024-7-16 20:08
我又重新试了下，把步长改为0.2fs，每步输出速度可以收敛。会有小部分case的caf文件std.dev列显著大于其他c ...

老师。我做了10组数据，对每一组进行积分，得到了以下的面积数据：
1、-0.02303(std先正后负，714.059开始转为负，最后为-3.27208) 舍去
2、0.5655(std全为正，最后为80.2761)
3、-0.24452(std先正后负，但正的很少，大部分为负，从24.648开始转为负。最后为-34.7106) 舍去
4、-0.17958(std先正后负，386.723开始转为负。最后为-25.4915) 舍去
5、-0.18859(std先正后负，138.034开始转为负。最后为-26.7718)舍去
6、-0.22724(std先正后负，460.708开始转为负。最后为-32.2557)舍去
7、0.03704(std中间有一点为负，大部分为正，最后为5.25587)
8、0.2163(std全为正，逐渐增大，最后为30.6996)
9、-0.08554(std先正后负，34.024开始转为负，最后为-12.1424)舍去
10、0.08845(std全为正，最后为12.5571)
舍去数据后留下了:
2、0.5655(std全为正，最后为80.2761)
7、0.03704(std中间有一点为负，大部分为正，最后为5.25587)
8、0.2163(std全为正，逐渐增大，最后为30.6996)
10、0.08845(std全为正，最后为12.5571)
但是2、8偏离较大，是不是就只剩下：
7、0.03704(std中间有一点为负，大部分为正，最后为5.25587)
10、0.08845(std全为正，最后为12.5571)
这两组有效数据呢？
但是这样是不是就太少了数据，还要继续做模拟呢
作者
Author: jlyjlysjd    时间: 2024-7-23 14:06
本帖最后由 jlyjlysjd 于 2024-7-23 14:56 编辑
二分音符 发表于 2024-7-16 20:08
我又重新试了下，把步长改为0.2fs，每步输出速度可以收敛。会有小部分case的caf文件std.dev列显著大于其他c ...

老师，请问标准差大于多少就舍去呢
我做了20次：
(, 下载次数 Times of downloads: 48)
不太知道那些面积是有效面积（左边是面积，右边是最后2000ps的标准差）
作者
Author: 二分音符    时间: 2024-7-25 17:15
现在的数据波动还是太大了，可能还是没有物理意义。我不确定继续减小步长是否可行。
作者
Author: jlyjlysjd    时间: 2024-7-25 19:16
二分音符 发表于 2024-7-25 17:15
现在的数据波动还是太大了，可能还是没有物理意义。我不确定继续减小步长是否可行。

好的，谢谢老师!
作者
Author: 梵墨    时间: 2024-10-31 09:49
二分音符 发表于 2024-7-8 11:52
最重要的一点是nstvout改成5，输出间隔必须非常小。展示Fig.6a时只需要显示几个ps的数据，但是积分的时候得 ...

此mdp中的constraint用了all-bonds? 确定这样没问题么
作者
Author: 二分音符    时间: 2024-10-31 18:43
梵墨 发表于 2024-10-31 09:49
此mdp中的constraint用了all-bonds? 确定这样没问题么

我的体系测试过all-bonds和h-bonds算出来的性质都几乎没有差别，如果其他case有需要的话还是要验证。



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3