求助：固液相转变的分子动力学模拟

qczgzly

    各位坛友好，我在做固液相转变的动力学模拟中遇到了困难，在此向大家求助。
    体系为一个COF，仅含C、H、N、O、S元素，该COF在常压1300K可熔化。该COF可吸附有机强酸如甲基磺酸（MSA，吸附比例为晶胞:MSA=1:3），或可通过分子间质子转移形成阴阳离子[COFH]+[MS]-。吸附甲基磺酸的COF常压下熔点仅为500K。希望通过模拟得到的结果为，COF在300K、1000K保持框架结构不变，在1700K收缩成无定型结构；[COFH]+[MS]-在300K保持框架结构不变，在1000K、1700K收缩成无定型结构。
    我目前已搭建好COF和[COFH]+[MS]-的周期性结构，力场通过sobtop调用GAFF得到，原子电荷为晶胞单元用H饱和后B3LYP/6-31G(d)计算esp电荷拟合所得，最后调整至晶胞单元整体电荷为整数。md在NPT系综下分别设置不同温度，其它参数设置如下：

1. integrator  = md
   
2. nsteps      = 4000000   ; 4 ns
   
3. dt          = 0.001
   
4. ; Output control
   
5. nstxout             = 40000
   
6. nstvout             = 400000
   
7. nstenergy           = 4000
   
8. nstlog              = 4000
   
9. nstxout-compressed  = 40000
   
10. compressed-x-grps   = System
    
11. ; Bond parameters
    
12. continuation            = yes
    
13. constraint_algorithm    = lincs
    
14. constraints             = h-bonds
    
15. lincs_iter              = 1
    
16. lincs_order             = 4
    
17. ; Neighborsearching
    
18. cutoff-scheme   = Verlet
    
19. ns_type         = grid
    
20. nstlist         = 40
    
21. rcoulomb        = 0.9
    
22. rvdw            = 0.9
    
23. ; Electrostatics
    
24. coulombtype     = PME
    
25. pme_order       = 4
    
26. fourierspacing  = 0.16
    
27. ; Temperature coupling is on
    
28. tcoupl      = V-rescale
    
29. tc-grps     = AZA       MS-    ; AZA is [COFH]+
    
30. tau_t       = 0.2       0.2
    
31. ref_t       = 1098      1098
    
32. pcoupl              = Parrinello-Rahman
    
33. pcoupltype          = isotropic
    
34. tau_p               = 6.0
    
35. ref_p               = 1.0
    
36. compressibility     = 4.5e-5
    
37. pbc     = xyz       ; 3-D PBC
    
38. periodic-molecules      = yes
    
39. DispCorr    = EnerPres
    
40. gen_vel     = no

复制代码

因为运行时警告体积变化超过1%，故将tau_p调大至6.0以避免报警。Berendsen压浴无法正常跑md，故将pcoupl设置为P-R。但在此设置条件下COF和[COFH]+[MS]-在不同温度点的表现均相同。
   我怀疑是体系结构的问题，于是又将周期性结构调整为超胞，边界用氢饱和，这样体系就变成了一个大分子，并在.mdp参数中将periodic-molecules调为no。但COF和[COFH]+[MS]-在不同温度点仍然跑出了相同的结果，或者都保持框架结构，或者都缩成一团。
   我想请教的问题是，我在体系和参数设置上是否存在问题？有一个教学教程是模拟升温过程中冰变成水，我感觉我做的体系与此相似，但实在找不到问题所在。还请大家多多提出宝贵意见或建议，谢谢！

sobereva

我不知道你的COF具体是什么结构。如果是纯共价相连的，靠GAFF力场描述熔化没戏，得用CP2K跑AIMD（DFT跑不动的话至少也用GFN1-xTB或SCC-DFTB）。

bobosiji

sobereva 发表于 2022-12-14 13:17
我不知道你的COF具体是什么结构。如果是纯共价相连的，靠GAFF力场描述熔化没戏，得用CP2K跑AIMD（DFT跑不动 ...

如果熔化不涉及共价键断裂（比如2D cof，层与层之间没有化学键），classical力场比如gaff也没戏么？

qczgzly

sobereva 发表于 2022-12-14 13:17
我不知道你的COF具体是什么结构。如果是纯共价相连的，靠GAFF力场描述熔化没戏，得用CP2K跑AIMD（DFT跑不动 ...

谢谢Sob！COF是纯共价键相连。那我试试AIMD。可否追问一下，为什么GAFF力场不能描述这种COF体系的熔化？
我在模拟中发现，如果将体系设置为周期性结构，则确实看不到熔化现象。如果将体系截断为超胞大分子，边界用氢饱和，那么高温下还是能看到熔化，即体系失去框架和层状结构。

sobereva

qczgzly 发表于 2022-12-17 20:35
谢谢Sob！COF是纯共价键相连。那我试试AIMD。可否追问一下，为什么GAFF力场不能描述这种COF体系的熔化？
...

本身GAFF的初衷主要描述的是一般有机体系在常温附近的特征
高温情况下很多变量可能涉及到偏离势能面极小点很远的地方，GAFF对这些区域可能描述得很不足。如果是牵扯到pi共轭，可能偏得更多

qczgzly

sobereva 发表于 2022-12-18 11:20
本身GAFF的初衷主要描述的是一般有机体系在常温附近的特征
高温情况下很多变量可能涉及到偏离势能面极小 ...

明白了。再次感谢Sob！