CP2K模拟LiPF6电解液体系AIMD温度不合理

aipacino

各位老师好：
我是用packmol生成的2mol/kgLiPF6盐＋CCOC:CCOCC=1:1电解液盒子，包含553个原子，用multiwfn生成的cp2k.inp输入文件，先进行了结构优化（没有优化cell），用优化后的结构做AIMD模拟了2000步，理论方法是PBE，热浴选择CSVR, TIMECON 100。温度设定298.15K，模拟结果首先下降到200K以下，最后稳定在230K附近，请问是什么原因不能保持在我的设定温度？

另外想请教一下，我的目的是想探究Li离子在不同溶剂组成的电解液中溶剂化结构，Li离子迁移速率，电导率等性质，请问AIMD的结果是否/如何体现这些性质？我初始结构是不是应该把Li原子放在电解液盒子外边？

新手入门AIMD，如有描述不当请老师们见谅！

Uus/pMeC6H4-/キ

没开色散校正的PBE泛函并不足以正确描述溶液里的相互作用，参考http://sobereva.com/413，所以动力学行为当然会有问题。另外NVE系综根本没法控制温度，要在298.15 K下平衡得用NVT系综（热浴就用CSVR）。除此之外，仅模拟2 ps或许并不足以使溶液达到平衡态，参考http://sobereva.com/627，应当继续延长。另外这种体系可以开启OT加快模拟的速度，参考http://sobereva.com/665，而且PRINT_LEVEL设为LOW就够了。

aipacino

Uus/pMeC6H4-/キ 发表于 2025-2-24 22:31
没开色散校正的PBE泛函并不足以正确描述溶液里的相互作用，参考http://sobereva.com/413，所以动力学行为当 ...

谢谢老师。
如果我想尽量大限度地提高模拟电解质溶液的准确性，应该用什么理论方法呢？
Print_level是我后续工作的需要MEDIUM及以上。
还有个问题就是我初始构型的Li离子应该放在哪里呀？

Uus/pMeC6H4-/キ

aipacino 发表于 2025-2-25 13:03
谢谢老师。
如果我想尽量大限度地提高模拟电解质溶液的准确性，应该用什么理论方法呢？
Print_level是 ...

（看来也是准备用CP2K的数据作为训练集构建机器学习势了吧）
这种体系的分子动力学就用PBE泛函加上D3(BJ)色散校正的理论方法即可，注意Multiwfn生成输入文件时用选项3设置色散校正。一开始锂离子和阴离子、溶剂分子等一样随机放在盒子里就行，反正跑一段足够长时间的模拟（平衡相）让体系达到平衡之后就没区别了，之后的轨迹（产生相）再用来研究溶剂化结构、计算离子迁移率和电导率等。

aipacino

Uus/pMeC6H4-/キ 发表于 2025-2-25 15:33
（看来也是准备用CP2K的数据作为训练集构建机器学习势了吧）
这种体系的分子动力学就用PBE泛函加上D3(BJ ...

哈哈哈是的老师
我已经加上D3校正和OT了
但是我看到别人的文献里初始构型是把锂原子和溶剂、阴离子严格分开，如下图（图片来自doi:10.1021/jacs.4c09027）。

我的理解是 AIMD模拟这个时候开始的Li原子从电极迁移到溶剂里直到形成稳定的溶剂化结构，统计这个过程中的Li离子轨迹计算迁移率，还有溶剂化能这些。
如果已经平衡了，再进行AIMD这些结构不是就不会发生大的变化了吗？还是说产生相我的输入文件会做相应的改动？
谢谢老师耐心解答！

Uus/pMeC6H4-/キ

aipacino 发表于 2025-2-25 16:44
哈哈哈是的老师
我已经加上D3校正和OT了
但是我看到别人的文献里初始构型是把锂原子和溶剂、阴离子严格 ...

啥，前面完全没提电极的事啊。你要研究的到底是哪个体系的锂离子：(1)溶解了锂盐的电解液，属于单一的均匀液相；(2)新鲜锂盐电解液和锂金属电极接触，除了(1)以外还形成固液两相界面；(3)实际的金属锂电池的负极，除(2)外还受到外电位限制？这几种体系模拟起来难度差别可大了。

(1)里单纯的锂盐电解液可以按不发生化学反应处理，还需要假定温度不很高、外电场均匀且不很强，等等。(2)里锂盐电解液和锂金属电极会发生化学反应生成固体电解质界面层（SEI），涉及一大堆溶剂和锂盐发生化学反应的产物。(3)更加复杂了，参考http://bbs.keinsci.com/thread-50714-1-1.html的讨论，虽然说的是计算反应能量变化但同样涉及电化学。

电解液体相里“溶剂化结构、离子迁移率、电导率”等性质一般对应上面说的(1)，这个做普通的分子动力学是可以处理的。至于(2)和(3)的情形，涉及的变化太多，想用分子动力学合理描述并不简单。

aipacino

Uus/pMeC6H4-/キ 发表于 2025-2-25 17:24
啥，前面完全没提电极的事啊。你要研究的到底是哪个体系的锂离子：(1)溶解了锂盐的电解液，属于单一的均 ...

老师不好意思 前几天不被允许发帖
我想了一下，目前应该是符合您说的（1），所以不用单独做电极的模型，只做电解液的模型
我目前在298k跑了3ps的AIMD，然后拿这个数据训练MLP，但是喂给LAMMPS模拟温度异常奇高，所以不知道是我训练数据的问题还是训练设置的问题还是LAMMPS输入文件的问题
训练数据方面，您说“一开始锂离子和阴离子、溶剂分子等一样随机放在盒子里就行，反正跑一段足够长时间的模拟（平衡相）让体系达到平衡之后就没区别了”，但是我看到有帖子说这样有可能陷入局部最优的构型，让我先升温退火后再进行接下来的工作，您对这个说法怎么看？
我现在的想法是先用lammps做升温退火再长时间预平衡，再做AIMD

Uus/pMeC6H4-/キ

aipacino 发表于 2025-3-3 21:42
老师不好意思 前几天不被允许发帖
我想了一下，目前应该是符合您说的（1），所以不用单独做电极的模型， ...

常温下3 ps的AIMD轨迹数据作为训练集太少了，不足以从构型中得到相互作用势的完整描述，后续模拟很容易出现原子的异常团聚导致模拟崩溃（参考这段介绍）。我也在LAMMPS模拟的轨迹中观察到过这种现象，“温度奇高”算是体系出问题的症状之一而不是原因本身。解决这个问题需要扩充数据集，包括但不限于升温加快构型变化、延长模拟时间等（上面那个链接后面讲的DP-GEN的原理也很重要，但我还没用过这个，无法就此提供帮助）。另外训练NNP的过程中也注意关注能量和受力的RMSE变化，训练结束时NNP的预测和DFT的结果应当差异足够小。

“陷入局部最优构型”说的是几何优化找势能面上极小点的静态问题，但分子动力学本来就是遍历势能面的动态采样过程，一般热浴温度下溶液体系里的能量涨落交换频繁，离局部极小点远着呢（你看到的怕不是那个通过分子动力学得到无定形晶体结构的例子，与当前目的无关）。有关模拟流程，仔细阅读2楼链接的627号博文以及我4楼的最后一句话。

aipacino

Uus/pMeC6H4-/キ 发表于 2025-3-3 22:40
常温下3 ps的AIMD轨迹数据作为训练集太少了，不足以从构型中得到相互作用势的完整描述，后续模拟很容易出 ...

谢谢老师耐心解惑
因为我刚上手AIMD，但是对机器学习还是比较了解，相比于扩充数据集，我更担心的是我现有的3ps数据集的质量问题，我应该如何评估运行的AIMD结果的好坏呢？
比如那个陷入局部最优的可能，我用VMD看了我的构型变化，感觉很有道理，下面附上某教学视频截图和我的构象快照。
一共五个阴离子，4个都集中在了盒子的一侧，是否不太合理（我感觉和这个视频截图类似），3ps的动力学期间位置变化也不大，估计再延长也不会有什么变化。

Uus/pMeC6H4-/キ

aipacino 发表于 2025-3-4 12:30
谢谢老师耐心解惑
因为我刚上手AIMD，但是对机器学习还是比较了解，相比于扩充数据集，我更担心的是我现 ...

单看这最后一帧没感觉有什么异常（不用管某些奇特的成键显示，这种几何关系仅供参考，见http://sobereva.com/534），瞬时的聚集不足以说明模拟全程的趋势。本来分子大一点就运动得慢，模拟时间又不够长的话那自然没法充分遍历势能面、观察到很大位置变化，所以我说了要延长模拟时间。

要说数据质量的话，&MGRID下的CUTOFF倒是有点低了。截断能不够高的话原子能量及受力的计算精度不够，有可能导致虚假的动力学行为，而且受力可是训练集里至关重要的数据。严格的做法是测试一下能量和受力随截断能收敛的情况（参考https://manual.cp2k.org/trunk/methods/dft/cutoff.html，注意根据超算等不同运行环境随机应变）；若不想做测试，鉴于体系里氧和氟对截断能的要求高一些，也得提升到至少400 Ry（参考https://www.cp2k.org/faq:cutoff）。

aipacino

Uus/pMeC6H4-/キ 发表于 2025-3-4 13:49
单看这最后一帧没感觉有什么异常（不用管某些奇特的成键显示，这种几何关系仅供参考，见http://sobereva. ...

谢谢老师
我已经在延长时间继续运行了，看起来只有多做测试才能知道性能
等我先用延长时间的数据训练之后再把截断能提高重新做