请教一下，如何正确实现盒子的整体拉伸和蠕变模拟？

光荣道彭于晏

体系是PMMA聚合物盒子，盒子包含20条不包含侧链的分子链，每条分子链包含40个重复单元残基，共12040个原子。力场是GAFF，原子电荷是RESP2(0.5) 。GMX版本是2022.4（之前一直使用的2020版本在英伟达CUDA更新为12.0版本之后无法完成编译了，我暂时也解决不了），系统是RockyLinux 8 。
根据分子动力学培训班和初级量化班的课程，还有sob老师的建议，使用GMX的deform功能进行盒子的整体拉伸。以298.15K下的盒子拉伸和蠕变为例：
盒子初始状态是一个在298.15K下已经完全弛豫平衡的盒子，然后设置拉伸模拟的mdp文件：

1. md-elongation-298K-longterm-1.mdp
   
2. 
   
3. define =
   
4. integrator = md
   
5. dt         = 0.002  ; ps
   
6. nsteps     = 250000000 ; 500000ps 500ns
   
7. comm-grps  = system
   
8. energygrps =
   
9. ;
   
10. nstxout = 0
    
11. nstvout = 0
    
12. nstfout = 0
    
13. nstlog  = 0
    
14. nstenergy = 500
    
15. nstxout-compressed = 5000
    
16. compressed-x-grps  = system
    
17. ;
    
18. pbc = xyz
    
19. cutoff-scheme = Verlet
    
20. coulombtype   = PME
    
21. rcoulomb      = 1.0
    
22. vdwtype       = cut-off
    
23. rvdw          = 1.0
    
24. DispCorr      = EnerPres
    
25. ;
    
26. Tcoupl  = V-rescale
    
27. tau_t   = 0.2
    
28. tc_grps = system
    
29. ref_t   = 298.15
    
30. ;
    
31. Pcoupl     = Berendsen
    
32. pcoupltype = anisotropic
    
33. tau_p = 0.5
    
34. ref_p = 1.01325 1.01325 1.01325 1.01325 1.01325 1.01325
    
35. compressibility = 0 0 0 4.5e-5 4.5e-5 4.5e-5
    
36. deform       = 0 0 2.5e-6 0 0 0  ; 以2.5e-6nm/ps(2.5e-3nm/ns)速度压缩
    
37. ;
    
38. freezegrps  =
    
39. freezedim   =
    
40. constraints = hbonds

复制代码

deform功能通过控压部分的参数实现，目前这套参数是摸索出来的，跑着不会报错的，盒子可以按照设定在某个轴方向上拉伸到指定长度。但是，会有几个问题

a)  grompp时会出现两条warning

1. WARNING 1 [file md-elongation-298K-1.mdp, line 39]:
   
2.   All off-diagonal reference pressures are non-zero. Are you sure you want
   
3.   to apply a threefold shear stress?
   
4. 所有非对角线参考压力都是非零的。 您确定要施加三倍剪切应力吗？(谷歌直译）
   
5. 
   
6. 
   
7. WARNING 2 [file md-elongation-298K-1.mdp]:
   
8.   The Berendsen barostat does not generate any strictly correct ensemble,
   
9.   and should not be used for new production simulations (in our opinion).
   
10.   For isotropic scaling we would recommend the C-rescale barostat that also
    
11.   ensures fast relaxation without oscillations, and for anisotropic scaling
    
12.   you likely want to use the Parrinello-Rahman barostat.
    
13. 
    
14. Berendsen 恒压器不会生成任何严格正确的合奏，不应用于新的生产模拟（我们认为）。 对于各向同性缩放，我们建议使用 C-rescale 恒压器，它也可以确保快速松弛而不会振荡，对于各向异性缩放，您可能需要使用 Parrinello-Rahman 恒压器。（谷歌直译）

复制代码

b) 盒子在拉伸过程中，仅在z轴拉伸方向发生了变化，盒子在x轴和y轴方向的尺寸没有发生变化，按照正常的拉伸现象，这两个方向上的尺寸应该随拉伸逐渐缩小

c) 在所有的拉伸和蠕变动力学模拟中，盒子都会出不同程度的偏离矩形形态的情况，盒子的各个平面都变成了平行四边形。这种情况还有一定的随机性，完全相同的参数和初始帧，做两次模拟，盒子偏离矩形形态的程度不一样。（蠕变的mdp设置是在拉伸的基础上，仅把控压部分的参数改为deform= 0 0 0 0 0 0 ，即拉伸速度为0，如下）

1. Pcoupl     = Berendsen
   
2. pcoupltype = anisotropic
   
3. tau_p = 0.5
   
4. ref_p = 1.01325 1.01325 1.01325 1.01325 1.01325 1.01325
   
5. compressibility = 0 0 0 4.5e-5 4.5e-5 4.5e-5
   
6. deform       = 0 0 0 0 0 0  ; 以0nm/ps(0nm/ns)速度拉伸

复制代码

d) 拉伸之后，盒子内部出现了明显的空隙，分子链不是均匀地分布在盒子中——拉伸之后盒子左侧分子链非常稀疏，而右侧有一大团分子链，如图

我的问题是：
1，现有的参数应该如何进行修改，才能正确完成盒子在常压下的拉伸和蠕变过程？尤其是目前只有保证compressibility = 0 0 0 4.5e-5 4.5e-5 4.5e-5，才能保证模拟正常能跑完（虽然能正常跑完不报错，也不代表参数正确）。
2，在模拟过程中，是否有必要维持盒子保持严格的矩形形态，如何维持？
3，grompp时出现的warning是否正常，用-maxwarn直接忽略是否合理？
4，在经历长达500ns的模拟后，盒子出现的分子链分布不均匀的情况，是否是因为拉伸速度太快（弛豫不够充分）？是不是与拉伸参数没有正确设置导致的x轴和y轴方向没有随着z轴拉伸而收缩有关？

光荣道彭于晏

拉伸或蠕变之后，盒子偏离矩形形态的轨迹截图

光荣道彭于晏

@sobereva  麻烦老师了！

beibei319

为什么要设置剪切应力呀？一般设置x，y，z正应力就好了吧。以及提示用Parrinello-Rahman压浴就改成Parrinello-Rahman呗（Berendsen压浴多用于最初体系建立时，系统平衡后多用Parrinello-Rahman来模拟）
贴的图片是我做蠕变的时用的npt文件，希望能给你帮助

光荣道彭于晏

beibei319 发表于 2023-5-23 15:04
为什么要设置剪切应力呀？一般设置x，y，z正应力就好了吧。以及提示用Parrinello-Rahman压浴就改成Parrinel ...

感谢您提供的mdp参数。
1. 剪应力其实不是我要设置的，sob老师建议我使用deform功能对盒子进行整体的压缩或者拉伸，而在设置deform的mdp时，我参考了GROMACS论坛上的一个帖子，发帖人和我一样并不清楚相关控压参数是如何对盒子产生作用的，我把他的参数拿过来修改了一下，在忽略grompp的提醒之后，是可以拉伸的。当然，每次grompp都会有与剪应力相关的一条提醒，我一直也不理解到底是哪部分参数导致的，只能先保证模拟不报错。
2. sob老师课上也确实是讲了Berendsen压浴用于体系构建，平衡之后使用PR压浴，我在做拉伸或者蠕变时使用的初始模型就是用PR压浴在某个温度下平衡过的。但是，考虑到deform功能会导致体系内部压力出现剧烈变化，而PR压浴容易产生震荡，所以才会考虑用Berendsen。
3. 不知道您有没有做过拉伸或者压缩盒子的模拟，有没有相关的mdp参数，是不是方便供我参考一下，非常感谢！

beibei319

图片里这两行就是涉及拉伸率以及压力的参数。六个数据前三个是正应力，后三个是剪切力。对比你的mdp文件，压缩率一栏仅剪切方向设置了参数【compressibility = 0 0 0 4.5e-5 4.5e-5 4.5e-5】可以考虑把这行代码改为compressibility = 4.5e-5。
我给的mdp文件是做蠕变时的mdp文件，没有使用deform，考虑到蠕变解析是单轴恒定应力，因此如图片所示，x,y轴设置为大气压，z轴设置100bar的拉力。
做拉伸解析时我用了deform，你的代码按warning提示改正后应该可以正常运行，即compressibility = 4.5e-5，改用PR压浴。（PR压浴虽会产生震荡，但如果模拟能正常进行，结果也合理的情况下可以无视震荡）

开心の土豆雷

你好，我能问一下是怎么判断这个体系达到了弛豫平衡了吗？

foging

beibei319 发表于 2023-5-23 15:04
**** 作者被禁止或删除 内容自动屏蔽 ****

你好，为什么ref_p = 1 1 -100,这个-100是怎么来的？

beibei319

foging 发表于 2024-3-11 02:12
你好，为什么ref_p = 1 1 -100,这个-100是怎么来的？

我个人设置的解析条件…z方向施加100atm的拉力。

webridging

之所以在z轴发生变化，是因为z轴有deform且compressibility为0，compressibility为0那变形后控压就不会去压他吧？
同理，x，y没变化，是因为他们无deform且compressibility为0，compressibility为0控压就不会去压他吧？
至于形态发生剪切就好解释吧，是因为他们无deform且compressibility不为0，compressibility不为0控压就会去切他吧？
我是这样理解的，不知道对不对，希望大家指正。总体而言，手册意思是先给一个deform，变形之后，压浴根据compressibility缩放相应指标使得体系向相应压强平衡态发展，尽管deform是非平衡态。
手册Deformation can be used together with semiisotropic or anisotropic pressure coupling when the appropriate compressibilities are set to zero. The diagonal elements can be used to strain a solid. The off-diagonal elements can be used to shear a solid or a liquid.
For systems with interfaces, semi-isotropic scaling can be useful. In this case, the 𝑥/𝑦-directions are scaled isotropically and the 𝑧 direction is scaled independently. The compressibility in the 𝑥/𝑦 or 𝑧-direction can be set to zero, to scale only in the other direction(s).
anisotropic Same as before, but 6 values are needed for xx, yy, zz, xy/yx, xz/zx and yz/zy components, respectively. When the off-diagonal compressibilities are set to zero, a rectangular box will stay rectangular.
关于pullcode的When the distance between two groups is changed continuously, work is applied to the system, which means that the system is no longer in equilibrium. Although in the limit of very slow pulling the system is again in equilibrium, for many systems this limit is not reachable within reasonable computational time. However, one can use the Jarzynski relation 135 (page 516) to obtain the equilibrium free-energy difference ∆𝐺 between two distances from many non-equilibrium simulations: