计算化学公社

标题: CotA 三铜中心截断模型约束优化长期卡在 Maximum Force [打印本页]

作者
Author: zenmefeishi 时间: yesterday 21:31
标题: CotA 三铜中心截断模型约束优化长期卡在 Maximum Force

各位老师好，最近在做 CotA multicopper oxidase 的 trinuclear copper center（TNC）QM 模型优化，目的是后续用于 MCPB.py 参数化并转到 GROMACS。现在遇到几何优化长期不收敛的问题，想请教一下这种多铜截断模型应该如何处理。

我的 TNC 模型来自 CotA 晶体结构，截取 TNC 铜中心及第一配位层残基，共127 个原子，包含 3 个 Cu 和多个 His 配体以及桥氧/含氧配体。

当前计算设置为：

Gaussian 16
UB3LYP/6-31G(d)
Charge = +4
Multiplicity = 4
NoSymm
Integral=UltraFine
SCF=(XQC,MaxCycle=512)

最初先对原始结构进行了：

Stable=Opt Guess=Mix

波函数稳定性检查正常：

The wavefunction is stable under the perturbations considered.

SCF 本身也没有明显问题，通常约 20–30 个 cycle 可得到 SCF Done。

问题主要出现在几何优化。

原始晶体结构附近三个 Cu–Cu 距离约为：

4.283 Å4.640 Å4.667 Å

一开始对截断 cluster 做完全自由优化时，其中一个 Cu 配位片段会整体漂离，Cu–Cu 距离最终达到约 7–8 Å，可能是由于截断模型缺失蛋白环境约束。

后来固定了 8 个位于截断 His 边界位置的 C 原子，只优化其余部分。目前采用：

#p UB3LYP/6-31G(d) Guess=Read
Opt=(Cartesian,MaxCycles=100,MaxStep=1)
SCF=(XQC,MaxCycle=512) Integral=UltraFine NoSymm

加边界约束后，Cu 第一配位层没有再解离，Cu–N/Cu–O 距离基本保持在约 1.8–2.1 Å，但三个 Cu–Cu 距离逐渐变为大约：3.98 Å5.58 Å5.32 Å

与原始结构相比，两条 Cu–Cu 距离分别增加了约 0.6–0.9 Å。

目前几何优化的另一个问题是 Maximum Force 长期降不下来。例如 restart 后前几步大致为：

Step Max Force RMS Force
1    0.003352    0.000748
2    0.003353    0.000735
3    0.003334    0.000667
4    0.003330    0.000499

Gaussian 默认收敛阈值为：

Maximum Force 0.000450RMS Force 0.000300

RMS Force 有下降，但 Maximum Force 一直卡在约 0.0033，主要残余梯度集中在个别配位 His 原子附近。能量总体仍在下降，SCF 每一步都正常收敛。

前一轮优化跑到 100 多步后，Berny approximate Hessian 曾出现较大的负 eigenvalue，随后出现：

Skip linear search -- no minimum in search direction.Steepest descent instead of Quadratic search.

然后我从较低能量、结构仍正常的一个中间几何重新开启了新的 constrained optimization。

目前主要有几个问题想请教：

对这种从蛋白质中截取出来的多铜 QM cluster，仅固定截断边界原子是否足够？还是应该增加第二配位层/扩大 QM cluster，或者采用其他约束方式？
原始 TNC 的 Cu–Cu 距离约为 4.28/4.64/4.67 Å，而 constrained QM 优化后稳定在约 3.98/5.58/5.32 Å。虽然 Cu–N/O 第一配位层仍完整，这种程度的 Cu–Cu 变化是否已经说明模型发生了不合理重排？
Maximum Force 长期约为 0.0033、而 SCF 和能量都正常时，继续优化是否还有意义？
对这种体系是否有必要尝试其他优化算法，或者重新选择约束原子？我的最终目的主要是获得适合 MCPB.py Hessian/force constant 参数化的几何。

大佬们是否有什么建议。谢谢！

作者
Author: sobereva 时间: yesterday 22:36
你的描述方式过于繁复抽象，写一堆文字，不如直接把你的gjf文件、输出文件（超过500KB就先压缩）上传，内行人一看就知道有没有硬伤、哪里需要修改、怎么计算最得当

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