layout: post title: 力场拟合工具xff开发杂记 categories:
xff是一个用于拟合力场参数的在线工具, 可以生成GROMACS格式的拓扑文件. 目前它的功能尚简陋, 只能基于力常数矩阵(力常数为体系能量对坐标的二阶导数, 力常数矩阵常称为Hessian, 中文也有称海森矩阵, 或海森, 海塞的, 本文中称为海森, 以避免与力场参数中的力常数混淆)生成力场的简谐势参数. 后续可能会添加更多的功能, 也可能不会, 凭天意.
这篇文本是我开发xff工具过程中遇到或想到东西的杂乱记录, 算不得正式说明. 我暂时也没有时间将它爬梳成连贯有序, 逻辑自洽的正式文档, 那就采取这种所谓的论语体吧, 想到哪是哪, 想到哪说哪. 既然是杂记, 难免存在思维跳跃之处, 希望那只是因为我头脑活跃, 而不是因为痴呆健忘.
我不想把它梳理成比较正式的文档还有一个私心, 是我希望开发过程透明, 尽量展示而不是隐藏其中的一些细节, 以便有需要的人可以更清楚地了解实际开发过程是如何进行的, 并能从中学到些东西. 在正式发表的论文中, 无论公式推导还是实验步骤, 多是大混乱后清理出的道路, 看起来逻辑合理, 走起来干净通畅, 但实际过程往往远非如此. 数学王子高斯, 后人评价他像只狐狸, 会清理掉自己走过的痕迹, 单留下漂亮的结果让人赞叹而无法学习. 要学习, 只得去研究他未烧掉的草稿, 甚是可惜. 我当然没有狂妄到无知的地步, 只是在这里留下草稿, 凭需要者自取.
卢天(网名Sobereva)开发了sobtop, 可以基于海森矩阵拟合力场简谐势的力常数. 我试用了一下, 功能可以, 使用起来繁琐点但影响不大. 因为相应的源代码暂时还没有看到, 无法确知具体是如何实现的, 但猜测是调用了ambertools中的atomtype工具来获取分子中每个原子的GAFF原子类型, 自写或集成了LAPACK之类的线性代数库进行矩阵对角化, 此外还涉及一些UFF的功能, 自写的或是基于OpenBabel, RDKit之类工具.
邹文利(网名Zork, 翻译了高斯软件的手册, 维护着量子化学软件中文网, 算是我的前辈了, 我从读博到现在一直都在使用他翻译的高斯中文手册)以前做过内坐标力常数拟合类的研究, 比较理论化, 有论文发表, 代码也有发布. 我曾向他请教过相关的问题. 他的那些方法我不知道sobtop是否集成了.
我找来Seminario方法的相关论文看了看, 一篇是方法的原始文献, Calculation of intramolecular force fields from second-derivative tensors, Int. J. Quantum Chem. 60(7):1271-1277, 1996, 一篇是改进方法的论文, Harmonic Force Constants for Molecular Mechanics Force Fields via Hessian Matrix Projection, J. Chem. Theory Comput. 14(1):274-281, 2017. JCTC这篇论文我以前下载过, 只是当时没怎么在意, 因为我觉得这种方法局限性很大. 现在看来, 使用这种简单的方法生成力场参数还是有意义的, 至少可用于补充力场中一些缺失的参数, 或初步解决各种材料的力场参数化问题. 金属蛋白参数化常用的MCPB.py中听说也使用了这种方法, 但我没有去查证.
因为方法看起来不难实现, 所以我觉得自己也可以做个在线工具实现类似的功能. 既然是在线工具, 那就使用js吧; 既然使用了js, 那就不奢求速度了. 这种工具不大可能每个人运行成千上万次, 也不怎么需要自动运行, 所以速度不是很关键. 就把它当成个计算器, 需要用到的时候打开按两下, 最多等个几秒, 得到结果就够了. 也没有必要支持命令行, 敲一堆数字或代码才得到结果也省不了多少时间.
我做这个工具的另一个目的是想验证一下显示分子结构的ChemDoodle库能否指定成键. 如果可以, 那后续可以使用它做更多东西.
既然决定做这么个工具了, 那就先去研究下那两篇论文, 弄明白方法的原理.
JCTC论文的作者给了实现的matlab脚本和python脚本. 我测试了python脚本. 这个脚本支持高斯09的fchk文件, 如果使用高斯16的fchk文件, 需要删除其中的几项, 保证各项顺序与高斯09的一致.
弄明白方法的原理之后, 有了大致的实现流程: 读入海森矩阵, 抽取部分海森矩阵, 对角化部分海森矩阵获得其特征值和特征向量, 带入公式计算力常数.
先实现前面的部分: 读取高斯软件生成的fchk文件中的海森矩阵, 然后根据指定的原子编号抽取部分海森矩阵. 为此, 要弄清楚高斯fchk文件中海森矩阵的存储方式以及顺序. 因为整个海森矩阵是对称的, 所以一般只需要存储矩阵的一半, 上三角部分或下三角部分, 高斯手册中的相关说明为((F(J,I),J=1,I),I=1,3Natoms), 这样对于3x3矩阵, 各元素顺序为11 12 22 13 23 33, 相当于下三角部分. 要验证这个顺序是否正确, 或程序读取的数据是否正确, 可以对角化整个海森矩阵, 得到的特征值对应于对分子的振动频率, 结果应该与高斯给出的一致. 但这种验证方法比较麻烦, 而且后面也可以通过与已有计算结果对比进行验证, 所以就不使用这种麻烦的验证方式了.
抽取部分海森矩阵后需要对其进行对角化. 部分海森矩阵为3x3的实矩阵. 这种低维的矩阵, 特别是其维数又与我们所在空间的维数一致, 所以值得专门研究. 我以前开发分子尺寸大小的计算工具时用到过3x3矩阵的对角化. 论文Efficient numerical diagonalization of hermitian 3x3 matrices专门研究了3x3实对称矩阵的各种对角化方法, 发现其中最精确的是Jacobi旋转方法. 但这种方法只能用于实对称矩阵, 而部分海森矩阵一般不是对称的, 所以无法使用Jacobi旋转对角化方法, 而是需要使用非对称实矩阵的对角化方法.
那就看看对于3x3实矩阵, 有没有比较简单, 容易实现的对角化算法. 如果有的话, 就没有必要去折腾那些复杂的算法, LAPACK之类的东西了. 发现一篇论文Symbolic spectral decomposition of 3x3 matrices, 专门讨论了3x3实矩阵的解析对角化方法, 给出了具体详细的公式, 其中涉及的三次方程求根公式我以前也接触过, 比较熟悉, 所以我觉得基于这种方法对角化3x3实矩阵应该可行.
在对矩阵进行对角化之前, 我有个疑问, 是不是所有方阵都可以对角化? 答案是否定的, 只有半正定方阵才可以对角化. 物理上的理解, 或在我们这个问题背景下的理解是, 部分海森矩阵的半正定大致意味着所涉及两个原子的位置稍有变化时, 整个体系的能量变化要大于等于零. 因为所用的整个海森矩阵是在体系的平衡位置(稳定点, 极小点)计算的, 所以任何原子位置的变化都会导致体系能量升高. 这样看来, 从整个海森矩阵中抽取的部分海森矩阵也应该是半正定的, 从而可以对角化. 但我疑心严格来说, 却也未必总是这样, 只是我无法证明.
对角化涉及的另一个问题, 可对角化矩阵的特征值, 特征向量一定为实数么? 答案也是否定的. 矩阵的特征值可以是复数, 但复数特征值总是与其共轭值一起出现. 这样对于3x3实矩阵, 出现非实数特征值时, 其特征值只能是一个实数, 两个复数, 且两个复数相互共轭, 对应的特征向量情况类似. 这与三次方程根的存在情况一致.
数学上的细节总是很麻烦, 能否对角化, 是否存在复数特征值的问题就先不考虑了, 假定所处理的矩阵性质足够好, 直接对角化看看结果如何再说.
参考论文实现了其中的naive解析对角化方法. 测试发现对某些矩阵, 得到的特征值精度可以, 但特征向量总是与其他第三方结果不一致. 折腾了几天, 试了几种不同的特征向量计算方法, 有时结果能对上, 有时结果对不上. 似乎超出我的能力范围了. 这才认识到, 正如论文所指出的, naive方法误差容易传播, 对病态矩阵无法保证精度, 可能会导致所得结果, 特别是特征向量完全错误.
重新调研实矩阵对角化方法, 看到了Numerical reciples上的相关说明:
You have probably gathered by now that the solution of eigensystems is a fairly complicated business. It is. It is one of the few subjects covered in this book for which we do not recommend that you avoid canned routines. On the contrary, the purpose of this chapter is precisely to give you some appreciation of what is going on inside such canned routines, so that you can make intelligent choices about using them, and intelligent diagnoses when something goes wrong.
The algorithms for symmetric matrices given in the preceding sections are highly satisfactory in practice. By contrast, it is impossible to design equally satisfactory algorithms for the nonsymmetric case. There are two reasons for this. First, the eigenvalues of a nonsymmetric matrix can be very sensitive to small changes in the matrix elements. Second, the matrix itself can be defective, so that there is no complete set of eigenvectors. We emphasize that these difficulties are intrinsic properties of certain nonsymmetric matrices, and no numerical procedure can "cure" them. The best we can hope for are procedures that don't exacerbate such problems.
The presence of rounding error can only make the situation worse. With finiteprecision arithmetic, one cannot even design a foolproof algorithm to determine whether a given matrix is defective or not. Thus current algorithms generally try to find a complete set of eigenvectors and rely on the user to inspect the results. If any eigenvectors are almost parallel, the matrix is probably defective.
简而言之, 如果对实矩阵不做任何假定, 对角化是个很困难的问题, 而且这困难是本质性的, 和具体问题, 具体算法无关. 这也是数值方法书上讨论的多是实对称矩阵对角化的原因, 因为这种情况下矩阵非病态, 性质良好, 有可靠的对角化算法. 这是我一开始没有料到的, 还以为是我使用的算法有问题, 通过使用更精细的算法就可以避免结果不一致的问题.
在对非对称矩阵进行对角化时, 要注意检查结果的正确性. 下面一些工具可用于交叉验证结果:
使用一些测试矩阵, 对比了几种对角化程序的结果, 确定我的实现没有太大问题. 结果不一致是因为对有些矩阵其特征向量的计算非常困难. 在某些情况下, C++ Eigen库自带的对角化代码也无法得到正确的结果. 我没有深挖Eigen所用的对角化代码是否来自LAPACK, 也没有确认LAPCK能否给出一致的结果, 但猜测是可以的, 因为NumPy使用的对角化代码来自LAPACK, 在大多数情况下能够给出一致的结果. LAPACK对角化程序的使用, 可以参考实对称矩阵本征值问题的理论与实践.
我没有能力用js实现一个与LAPACK结果完全一致的矩阵对角化程序. 一种可行方法是用WebAssembly将C++代码翻译成js, 计算速度会下降, 但对简单应用不成问题. 网上对此有些实现可以参考.
最开始实现的naive解析对角化方法, 在矩阵病态的情况下给出的结果不可靠. 基于矩阵QR分解的对角化方法更可靠一些, Numerical reciples上的代码可以参考. 决定还是使用这种方法吧, 虽然实现比较麻烦, 结果也仍然无法保证完全与LAPACK一致.
对角化问题上遇到的挫折促使我思考, 使用Seminario方法计算力常数是不是一定要基于对角化, 有没有其他可替代或更好的方法? 在对非对称的部分海森矩阵进行对角化之前, 能不能简单粗暴地先将其对称化, 也就是取部分海森矩阵与其转置的平均值. 这样, 既避免了能否对角化, 特征值是否为实数这些问题, 也可以使用最精确的Jacobi旋转对角化方法得到可靠结果. 这引出了一个新的问题, 对称化矩阵的特征值和特征向量, 以及各个矩阵的特征值和特征向量之间有什么关系? 或更一般的, 矩阵和的特征值, 与矩阵特征值的和, 二者之间有何关系? 数学上对此并没有太多结论, 我能找到一点信息. 此外, 已知的是, 如果可进行对角化, 矩阵与其转置具有相同的特征值, 但特征向量不一定相同.
和对角化很类似的另一种矩阵分解是奇异值分解. 这种分解总是可行的, 而且很稳定. 那在Seminario方法中能否使用奇异值分解代替对角化呢? 对于实对称矩阵, 奇异值分解和对角化没有区别, 对一般的实矩阵, 奇异值与本征值之间关系如何? 对这个问题我也只找到一点说明What is the Relation between Eigenvalues & Singular Values?, 结论是一般而言矩阵的奇异值不大于其特征值.
直接对角化部分海森矩阵, 还存在另一个问题, 那就是计算的力常数会因原子编号顺序不同而不同, 通常差距较小. 为避免这种情况, JCTC论文作者的做法是对所有可能编号顺序, 如A-B, B-A分别计算, 然后再将所得力常数进行平均, 这样结果就与原子编号顺序无关了. 换用基于对称化或奇异值分解的方法可以自然地避免这个问题.
根据前面的说明, 无论是对称化后再对角化, 还是使用奇异值分解, 得到的特征值或奇异值不会大于直接对角化得到的特征值. 如果这样的话, 考虑到使用对角化方法所得的力常数常偏大, 采用这两种替代方法所得的力常数可能会小一点, 说不准更合理. 实际结果如何需要做测试才能确定, 但最多也只能得到些经验的结论.
解决了或说不再纠结矩阵对角化的问题, 剩下的就比较简单了, 照着论文中给出的公式实现即可. 没什么难点.
现在回顾评价下计算力常数的Seminario方法. 优点是概念简单, 实现容易(不考虑对角化中的各种问题), 但缺点/局限性也不少. 方法只适用于简单的简谐势力场, 二面角也只能是不常用简谐式, 不适用于柔性大的分子. 对成环体系方法也存在限制. 如前面说的, 所得力常数会不对称, 与原子编号顺序有关. 所给力常数的值常过大, 因为忽略了原子间的非键相互作用. 所以这种方法只适用于获得比较粗略的力场, 用于金属蛋白时效果应该比约束好些, 但具体好多少, 难说.
此外, 方法只能使用平衡位置的海森矩阵, 无法使用非平衡点的海森矩阵, 从而也就无法使用多个海森矩阵. 在非平衡位置, 部分海森矩阵很可能不是半正定的, 无法对角化, 或者特征值为复数, 这些都是棘手的问题. 当然, 也可算是方法的改进方向.
再者, 计算大体系的海森矩阵还是很费力的, 如果所用计算方法不支持解析的二阶导数, 那就更费力. 所以使用的计算水平不能太高, 体系也不能过大. 考虑到前面所说的那些限制, 我也不确定使用更高水平的方法是不是一定能得到更好的力场参数.
得到力场参数后, 如果使用它们对分子进行优化或模拟, 所得构型可能会与原本的平衡构型有差距. 其中一个原因是拟合时忽略了非键相互作用, 而力场计算时是要考虑的.
对某些对称性比较高的分子, 使用对称化的部分海森矩阵, 得到的结果可能更稳定, 所以我觉得还是值得为这种方法单独列一个选项, 虽然具体情况需要考虑更多. 另外, 使用奇异值分解我觉得也是一个可以尝试的方向, 理论上也能说通, 如果没有人试过, 值得研究.
对于二面角项, 原始论文中有其力常数的计算公式. 但一般力场中很少使用简谐势二面角, 更多是使用周期势, 所以这种刚性二面角的力常数意义不是很大. 即便是简谐势二面角, 我猜测其力常数也存在和键角类似的问题, 需要进行一定的修正, 具体如何修正, 值得考虑.
当三原子共线时, 键角为180度, 这种情况下键角力常数公式中的某些向量无法良定义, 所以不能直接使用. JCTC论文作者的matlab脚本中没有处理这种特殊情况, 其python脚本中的做法则是生成单位球面上的随机向量, 并计算各自对应的力常数, 最后再取平均值. 这种方法实现倒也不麻烦, 计算耗时也不是大问题, 只是我对其得到的结果甚为怀疑, 所以我索性没有实现. 真要使用这种方法, 或许也有解析结果, 因为空间中两个向量之间的夹角具有特殊的分布, 参考我以前的说明.
对于180度的键角, 采用线性势可能更好. 我以前查看GROMACS源代码时发现其中有一种线性键角势函数, 但文档, 网上对此没有说明. 我自己曾做过一些测试, 试图搞清楚是怎么回事. GROMACS 2022文档中添加了有关的说明, 来自2020年的论文 A potential for molecular simulation of compounds with linear moieties, J. Chem. Phys. 153(8)：084503, 2020. 使用这种势函数模拟CO₂之类的线性分子就很简单了. 当然, 至于采用线性键角势效果具体如何, 是不是更好, 值得做些测试进一步探讨. 但我没有时间, 也没有精力做这些测试并写成论文, 直接在工具中实现了这种方法. 这种线性键角势的力常数同样存在偏大的问题, 如何修正值得考虑.
由于Seminario方法的内在局限性, 改进某些地方或许有益, 但终归有限. 子夏曰: “虽小道, 必有可观者焉; 致远恐泥, 是以君子不为也.”对某些方法, 要有恰当的认知, 了解其边界所在, 简单改进一点有效果就够了, 没有必要拘泥于此.
比这种简单的方法更进一步的海森拟合方法, 可以参考Partial Hessian Fitting for Determining Force Constant Parameters in Molecular Mechanics, J. Comput. Chem. 37(26):2349-2359, 2016. 计算流程稍复杂一点, 但原理上更合理可靠. 这也是我曾考虑过的方法中的一种. 如果借助MD引擎计算力常数而不是自己计算, 那这种方法的实现还比较容易, 有能力的可以试试.
科学与传统学术的一个很大不同在于要始终保持前沿性, 而不是抱着已有的东西当成宝. 科学是求新知, 而不是博雅知. 科学要在state-of-the-art的基础上更进一步, 力求有所突破. 即便最终证明你更进的这一步是错的, 那也是进步必付的代价, 同样锻炼了你自己的思考能力, 证实了你不是单纯地在重复. 科学的孩子, 即便哭也要用自己的声音.
这也是为什么我一直反对将问题简单化的原因. 爱因斯坦说过, EVERYTHING SHOULD BE MADE AS SIMPLE AS POSSIBLE, BUT NOT SIMPLER. 简单化有助于系统地教授和学习, 但无助于培养科学求新能力. 现实很繁芜, 你要在其中找到属于自己的路, 走自己的路, 而不是人走得多了的路. 这样看来, 强求一致, 一刀切的思想对求新没有太多益处, 但也不要就此多演绎. 因为现实太过沉重, 不容你多做臆想.
见gmxtools.
Reset Bonding: 读入fchk文件后会自动根据原子间的距离确定成键. 成键距离标准基于原子的共价半径, 数据来自XTB代码, 引自Pyykko and Atsumi, Chem. Eur. J. 15, 2009, 188-197. 成键距离的容差为10%, 即如果两原子之间的距离小于两原子共价半径之和的1.1倍, 则认为它们之间成键. 对某些分子, 如果使用此标准确定的成键关系不正确, 可以自行在原子成键表中添加, 然后点击此按钮即可更新. 原子共价半径的实际大小可设置Model为covalent spheres查看.
Reset FFindex: 如果计算体系抽取自更大体系, 将得到的拓扑装配回原力场拓扑时要保证每个原子的编号与原来一致, 为此, 工具提供了指定每个原子在原体系中编号的功能, 输出拓扑中会使用这些编号. 这在处理金属蛋白时可能用的上. 出于同样的考虑, 也可以指定每个原子力场原子类型的功能. 具体的原子类型可以通过其他程序获知.
Scaling Factor: 用于对海森矩阵的值进行缩放, 力常数的所用的缩放值为其平方
Selected Atoms: 可以选择只在拓扑中只输出涉及某些原子的项, 方便装配回原力场. 选择基于从1开始的原子编号, 支持all和范围, 如all 1:5:2, 反选可使用-, 如-all -1 -2:5. 从左到右依次解析. 默认全选.
Diagonalization: QR直接对角化, Jacobi对称化后再对角化
c6h6.fchk来自JCTC作者所给示例, 频率缩放因子为0.957.
sf6.fchk来自sobtop所带示例, 频率缩放因子为1.