计算化学公社

标题: 已解决_Gromacs 中分子动力学模拟聚合物中分子链几乎没有运动 [打印本页]

作者
Author: sarphuart 时间: 2025-5-1 20:35
标题: 已解决_Gromacs 中分子动力学模拟聚合物中分子链几乎没有运动
Gromacs 分子动力学模拟聚合物中分子链无变化。研究不同温度对单个分子聚合物链构象的影响，系综为 NPT，一个温度一次分子动力学。聚合物为 30 个重复单元的 N，N-异丙基丙烯酰胺（该分子链大于 305 K 分子链会卷缩，小于305 K 会舒展），聚合物使用 OPLS-AA/L 力场，电荷使用 XTB-RESP 。水力场使用 TIP4P/ICE，水分子数大概 23000，使用 solvate 填充水。目前已做温度为 298 K，305K，310 K 的 30 ns 的分子动力学模拟，三次结果都是分子链略微抖动，尝试延长模拟时间到 100 ns，分子链也没有大的变化。
聚合物动力学参考文献： https://doi.org/10.1039/C8CP00537K
Tip4p/ice itp 文件来源：https://zhuanlan.zhihu.com/p/75785937
具体操作步骤为
gmx editconf -f UNK.gro -o box.gro -c -d 1.8 -bt cubic # 将聚合物放入盒子
gmx solvate -cp box.gro -cs TIP4P.gro -p UNK.top -o solvated.gro # 加水
gmx grompp -f minim.mdp -c solvated.gro -p UNK.top -o em.tpr -maxwarn 1 # 能量最小化
gmx mdrun -deffnm em -v
gmx grompp -f nvt.mdp -c em.gro -p UNK.top -o nvt.tpr # nvt 预平衡
gmx mdrun -deffnm nvt -v
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -p UNK.top -o npt.tpr -maxwarn 1 # npt 预平衡
gmx mdrun -deffnm npt -v
gmx grompp -f md.mdp -c npt.gro -p UNK.top -o md.tpr # md 计算
gmx mdrun -deffnm md -v
本人第一次模拟，请各位老师指教，感谢。
以下是所有模拟过程中所需文件：

作者
Author: sobereva 时间: 2025-5-1 22:19
若非要用OPLS-AA/L力场，应当结合1.2*CM5原子电荷，看
计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本
http://sobereva.com/585（http://bbs.keinsci.com/thread-21462-1-1.html）
ORCA结合Multiwfn计算RESP、RESP2和1.2*CM5原子电荷的懒人脚本
http://sobereva.com/637（http://bbs.keinsci.com/thread-28178-1-1.html）

RESP、RESP2(0.5)电荷适合结合GAFF力场，拓扑文件可用http://sobereva.com/soft/Sobtop产生。而基于xtb程序的波函数算的RESP电荷质量太烂，强烈不建议使用，怎么算质量像样的RESP电荷看：
RESP拟合静电势电荷的原理以及在Multiwfn中的计算
http://sobereva.com/441（http://bbs.keinsci.com/thread-10880-1-1.html）
计算RESP原子电荷的超级懒人脚本（一行命令就算出结果）
http://sobereva.com/476（http://bbs.keinsci.com/thread-12858-1-1.html）
RESP2原子电荷的思想以及在Multiwfn中的计算
http://sobereva.com/531（http://bbs.keinsci.com/thread-16190-1-1.html）
ORCA结合Multiwfn计算RESP、RESP2和1.2*CM5原子电荷的懒人脚本
http://sobereva.com/637（http://bbs.keinsci.com/thread-28178-1-1.html）

你的mdp里全都设的是ref_t = 280。
而且倘若你已知比如实验上305K才是相变温度，模拟用的温度至少应当比这高10度，否则可能体系可能处于高度有序的晶态而缺乏流动性。而且用的力场、原子电荷越不理想，模拟时的相变温度可能比实验偏差得越大。

也没给体系完整结构文件，没法判断结构是否存在硬伤。

用TIP4P/ICE完全莫名其妙。又不是模拟冰，用专门给模拟冰点附近的这个水模型干嘛，而且四点水模型耗时还高。用SPC/E或OPC3才是优先应当考虑的选择。不要在没有水模型常识的情况下乱用水模型。PS：北京科音分子动力学与GROMACS培训班（http://www.keinsci.com/KGMX）里专有一节“溶剂模型” 极为全面介绍了所有水模型，并且详细讲解了如何正确计算。

作者
Author: sarphuart 时间: 2025-5-7 16:57

sobereva 发表于 2025-5-1 22:19
若非要用OPLS-AA/L力场，应当结合1.2*CM5原子电荷，看
计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒 ...

非常感谢 sob 老师的回复，按照您的建议，我这几天对模拟体系进行了调整，目前动力学行为已经初步呈现出预期的构象变化趋势。不过在进一步分析过程中，仍有以下几个概念性问题希望向您请教：
1. 关于回旋半径的波动问题：在使用 OPLS-AA/L 力场结合 SPCE 水模型进行模拟时，我发现所得回旋半径（Rg）数据的时间序列波动较大。请问这是否可能是由于 mdp 文件中的参数（如温控方式、积分步长或约束设置）不合理所致？或者这是否反映出 OPLS-AA/L 力场本身在描述该聚合物体系构象柔性方面存在局限？
2. 关于“完整结构文件”的理解：您提到需检查“体系完整结构文件”是否存在问题。我理解这里所指的应是 .gro 和 .top 文件在 GROMACS 模拟准备过程中可能存在的结构残缺或拓扑拼接错误。但由于我的模拟流程从初始构建、能量最小化到动力学阶段均可顺利完成，尚未出现报错，因此这个“完整结构文件”的具体是指哪些不需要依赖gmx自带的力场库多出来的类似于 ATP.itp 文件吗
3. 关于盐溶液体系下的力场适用性：后续我计划模拟聚合物在含氯化钙溶液中的构象变化行为。请问在此类含多价离子的体系中，OPLS-AA/L 结合 1.2 × CM5 原子电荷的方案是否依然适用？是否还需要针对 Ca²⁺ 离子进行额外参数优化或验证，特别是在考虑与聚合物羧基等极性基团的相互作用时
最后，附上当前回旋半径分析结果的波动趋势图

作者
Author: sobereva 时间: 2025-5-8 00:12

sarphuart 发表于 2025-5-7 16:57
非常感谢 sob 老师的回复，按照您的建议，我这几天对模拟体系进行了调整，目前动力学行为已经初步呈现出 ...

1 可能只是偶然性，恰好跑到了某种构象对应的势阱里。可以结合轨迹动画弄清楚情况

2 我就是指gro文件。别人没有结构文件、看不到实际模拟体系具体是什么样，就只能猜着回答，缺乏针对性

3 适用，但只用于描述聚合物。Ca2+用OPLS-AA目录下ions.itp里自带的参数，或者其它离子力场

作者
Author: sarphuart 时间: 2025-5-8 10:02
学生明白了，谢谢 sob 老师

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