氧化铝表面羟基O-H键在NPxyLzT下异常快速旋转报错

chieko

在做一个氧化铝固液界面的模拟，现在先只做羟基化氧化铝的模型。力场用CLAYFF，使用Sobtop生成拓扑文件，再跑em，nvt和npt平衡。最后产生相也是npt。em和nvt目前看都正常，但是npt总是报错，提示表面羟基的O-H键bonds that rotated more than 30 degrees。具体相关信息如下：
1、模型：在MS中切割氧化铝0001表面，暴露O，加H，上下加真空层，用castep优化下结构。
2、用Sobtop基于CLAYFF力场建立拓扑。原子类型、各种力场参数我都反复检查了，没有问题。
3、em：不约束hbonds，emtol  = 50.0，完成后发现固体扁了约2 A；（不知道正常与否）
4、NVT：约束hbonds，步长1 fs，v-rescale热浴，经过100 ps从0 K退火到298K，保持100 ps，全程200 ps。轨迹正常，固体形态维持，羟基H规律振动。
5、NPxyLzT：xy方向berendsen控压1 bar，z方向维持盒子尺寸（这是文献中的做法，我也试了各向同性控压，结果一样），步长1 fs，约束hbonds（这里和nvt我也试了不约束，因为步长不大，但是也一样报错），ref_p=5.0。接着nvt的继续算，报错。提示表面羟基的O-H键bonds that rotated more than 30 degrees。查看轨迹，模拟几乎刚开始（0.1 ps左右），部分表面的O-H键就开始快速旋转，并且这些羟基上O原子的位置也不能很好保持。

我反复确认过力场参数，应该没有问题（我也算过不加羟基的氧化铝表面，npt没有任何问题）。nvt也没有问题，npt下各向异性控压的mdp设置我觉得也没有问题也给盒子填满水计算，还是不行。

这么排查了一周多，始终无法解决，希望Sob老师和大家给点意见，谢谢！
所有结构、拓扑、mdp和运行产生的文件都打包在压缩包中了，terminal的报错在Typescript中。再次感谢各位！

nuanwen

你好！建立拓扑文件的时候，其它LJ_param，bonded_param等需要修改吗？能看一下私信吗

biogon

你用sobtop怎么搞的clayff的拓扑文件

chieko

Lacrimosa 发表于 2023-3-13 21:09
从你轨迹来看在崩溃之前表面的铝/氧原子都挤到一起了，所以羟基才会旋转。问题可能出在结构或者控压参数 ...

非常感谢您的建议！我明天来试试排查一下，感谢！

Lacrimosa

chieko 发表于 2023-3-11 22:06
文献中有类似的模型，是这么设置的，并没有说具体原因；我这里就当是先重复别人的模拟试一试了。
在我看 ...

从你轨迹来看在崩溃之前表面的铝/氧原子都挤到一起了，所以羟基才会旋转。问题可能出在结构或者控压参数的设置上。
你可以先检查一下你的初始结构，氧化铝的三个角度应该是90/90/120，而你的结构中参数都偏了一点，xy方向的结构不同了，使用xy方向耦合的控压方式是否合理？试试把晶胞正交化，做成矩形盒子再尝试。
另外需要确认一下compressibility设置的数值是否合理

chieko

Lacrimosa 发表于 2023-3-11 10:28
为什么要 xy方向控压 z方向维持尺寸? 不应该是反过来的吗？

文献中有类似的模型，是这么设置的，并没有说具体原因；我这里就当是先重复别人的模拟试一试了。
在我看来，对于我要做的东西，只NPT也是没有问题的，暂时并没有一定要固定z方向尺寸的理由。
当然如问题中所述，我也试了在NPT下的计算，会出现一样的问题。

Lacrimosa

为什么要 xy方向控压 z方向维持尺寸? 不应该是反过来的吗？