MSD计算离子电导率的问题

二分音符

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中的单位）

jlyjlysjd

二分音符 发表于 2024-7-5 17:31
最近自己动手做了下，完整步骤：
1. 计算5-10个NVT，速度输出频率做到0.1ps
2. gmx current -caf caf.x ...

谢谢老师！

二分音符

这篇工作和你想做的比较像，可以参考
https://www.sciencedirect.com/sc ... i/S0167732216332718

jlyjlysjd

二分音符 发表于 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里面的合并和计算多条曲线的均值，画出了下面的一幅图：
这个与文献的太不相符了，主要时间好像不太对，我的太长了。
将x轴方法，发现一点类似的规律。

所以老师，请问问题是出在哪里呢，是不是mdp中时间尺度出现的问题呢？根据文献“[color=rgba(0, 0, 0, 0.85)]MD 获得的 ECACF 在几皮秒内快速衰减到零”。是不是需要控制整体的模拟时间在3ps左右呢？或者是不是我的caf.xvg文件输出的有问题呢？
学生愚钝，请老师赐教，谢谢老师！


caf.xvg:

二分音符

最重要的一点是nstvout改成5，输出间隔必须非常小。展示Fig.6a时只需要显示几个ps的数据，但是积分的时候得做到1-2ns。
我的mdp文件给你参考下

jlyjlysjd

二分音符 发表于 2024-7-8 11:52
最重要的一点是nstvout改成5，输出间隔必须非常小。展示Fig.6a时只需要显示几个ps的数据，但是积分的时候得 ...

感谢老师的文件，经过修改mdp后，我总计运行了10次NVT后得到10组caf数据，均一化后得到如下数据：

放置于origin后，使用曲线求均值后得到了如下的图：
虽然后面确实是趋近于0的，但是是在0上下波动，并没有单调趋近，所以这个曲线不知道正不正确。
然后使用积分工具对0-2ns进行积分面积计算后，得到如下数据：

面积是：0.06043
通过同时,积分面积/3VkbT=0.06043/1.31E10=4.61E08
请问老师，这个电导率的计算哪里出现问题了呢？
是不是还是曲线的问题呢？


老师可以帮我看看我的NPT和NVT系综的mdp是否有问题呢？
非常感谢老师！

二分音符

注意单位，全部换成SI制。
你的ACF看起来还可以，再画一个积分值随时间变化的曲线（t时刻的函数值为ACF从0积分到t的面积）。
mdp看着没问题

jlyjlysjd

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

好的老师，下面是ACF(t)对t积分的图：


单位问题我计算如下：

所以需要乘以10^21。
通过gro文件获得体积：15.7514632602629 nm^3，Kb=1.380649 * 10^-23 ,T=298.15

得到了比较合理的数据。
谢谢老师！

jlyjlysjd

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

老师，我在算新的一个浓度的时候发现积分结果为负值，这说明结果是错误的吗

二分音符

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

你可以检查一下之前的结果，不同case之间算出的电导率差多少，然后画出errorbar（根据误差分析的知识），然后和文献里的errorbar比较一下。
如果波动极大，那么说明case还是少了，要么加case，要么加大模拟体系。要是还降不下来，那这个思路就走不下去了，得再换方法。

jlyjlysjd

二分音符 发表于 2024-7-11 16:52
你可以检查一下之前的结果，不同case之间算出的电导率差多少，然后画出errorbar（根据误差分析的知识）， ...

老师，我做的是0.5M MgCl2溶液的电导率，刚查阅了一下文献库，发现好像并没有这样类似的文献。
有一篇中文的文献，拟合了电导率：

我做了10次电导率的积分errorbar:

这样是不是可以说明波动是比较大的。

其实从10组电导率的数据来看，发现波动还是很大的：

但基本都是呈现负值情况，感觉再重复多次，是不是大概率还是负值呢？

二分音符

这种情况可能不适合用电流自相关算电导率了，因为数值误差过大，可能还是只能用NE方程计算。

二分音符

我又重新试了下，把步长改为0.2fs，每步输出速度可以收敛。会有小部分case的caf文件std.dev列显著大于其他case或者是负值，把这些舍掉。
这样做会导致一个case占据约2T硬盘空间，只能算完就删。

jlyjlysjd

二分音符 发表于 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)
这两组有效数据呢？
但是这样是不是就太少了数据，还要继续做模拟呢

jlyjlysjd

二分音符 发表于 2024-7-16 20:08
我又重新试了下，把步长改为0.2fs，每步输出速度可以收敛。会有小部分case的caf文件std.dev列显著大于其他c ...

老师，请问标准差大于多少就舍去呢
我做了20次：

不太知道那些面积是有效面积（左边是面积，右边是最后2000ps的标准差）