计算化学公社

标题: PM局域化优化到局部最优怎么办？ [打印本页]

作者
Author: Freeman 时间: 2024-2-16 11:37
标题: PM局域化优化到局部最优怎么办？
大家好。最近我需要手搓一个PM局域化的代码。目前用的是python的scipy.optimize.minimize函数来优化，参考DOI: 10.1002/jcc.23281的公式3、12～15（没有使用解析hessian，而使用数值的）。可是结果发现局域化优化到了局部最优而不是全局最优，具体情况就是如果体系存在对称性，相关的轨道都会混合在一起而非局域到一边。下面是两个例子，一个是4+2反应，一个是OTS，可以用multiwfn打开。
(, 下载次数 Times of downloads: 0)

(, 下载次数 Times of downloads: 1)
4+2反应体系是Cs对称的，每一个键局域轨道，都有上面一坨和下面一坨，但最优情况应该是上面一坨是一个轨道而下面一坨是另一个轨道。OTS情况好一些，毕竟整体不是对称的，但是局部对称性（一个甲基上的三个C-H）也会造成相似的情况，比如22和35轨道。

我一开始的猜测是，目前的情况是因为优化到了鞍点，于是多考虑了两点的数值hessian矩阵，优化算法也相应地改成了trust-ncg（带置信半径的牛顿共轭梯度法，具体算法我就不了解了），但是结果仍然如此。所以我觉得应该是优化到局部最优了。我又试过给初猜一个微扰来破坏对称性，结果虽然好了一些，但还是有很多轨道是分散成两坨的。例如还是这个4+2体系，43和45就明显混合在一起了。 (, 下载次数 Times of downloads: 0)

就我所知，至少社长是实现过PM局域化在multiwfn里的。毕竟PM局域化也挺悠久了，应该还有很多强人实现过。我想请教一下，像我说的这种情况该怎么办呢？

作者
Author: hebrewsnabla 时间: 2024-2-16 12:22
pyscf的实现 https://github.com/pyscf/pyscf/pull/1990

mokit的实现 https://gitlab.com/jxzou/mokit/- ... ref_type=heads#L400

作者
Author: sobereva 时间: 2024-2-16 16:27
Multiwfn用的PM原文的2*2旋转，基于你的ots_loc.mwfn按下文用Multiwfn做PM定域化后可以得到正常的C-H定域化轨道
Multiwfn的轨道定域化功能的使用以及与NBO、AdNDP分析的对比
http://sobereva.com/380（http://bbs.keinsci.com/thread-6053-1-1.html）

作者
Author: zjxitcc 时间: 2024-2-16 21:36
上面大家都说的很详细了，原理我就不讲了，讲两个解决办法：
（1）直接用MOKIT的loc功能。
（2）你自己实现PM局域化方法，要么用2*2旋转写，要么用二阶算法写（一阶算法容易掉到鞍点）。

作者
Author: Freeman 时间: 2024-2-18 18:08
本帖最后由 Freeman 于 2024-2-18 18:11 编辑

@hebrewsnabla @sobereva @zjxitcc

谢谢大家。我看了hebrewsnabla分享的pyscf的贴子，看到用梯度优化的方法确实往往优化到局部最优，而贴子里提供的方法好像就是在梯度优化前，做一些2*2旋转，作为初猜。所以无论如何，2*2旋转都是必要的，所以我现在打算实现2*2旋转的算法。

有关2*2旋转实现pm局域化的，我找到两篇文献，分别是ER（ (, 下载次数 Times of downloads: 1) ）和PM（ (, 下载次数 Times of downloads: 2) ）的原文。ER文章给出了2*2旋转的大致算法，而PM文章给出了一些2*2用于PM局域化的相关参数公式（Ast，Bst）。还有之前参考的梯度优化的文献SUSI（ (, 下载次数 Times of downloads: 1) ），给出了Mulliken和Lowdin电荷矩阵元的公式。

我试着写了2*2旋转的代码，但是运行出来的结果都是错的，所以想在这里复述以下这个算法的思路，请大家帮我看看有什么理解错误的地方。

1、循环遍历所有轨道对（i和j），找到使得Aij+(Aij+Bij)^2最大（ER eq26）的那对，其中Aij和Bij的定义是（PM eq29），Aij和Bij中的电荷矩阵元<i|P_A|j>可由（SUSI eq29-30）算得；
2、针对这对轨道ij，通过ER eq19，找到一个处于0-pi/2之间的alpha值，代入ER eq25混合ij两轨道成新的两个轨道，替换掉原有的ij轨道；
3、基于新的轨道，重新算Mulliken/Lowdin电荷矩阵；
4、重复1-3步直到轨道几乎不再变化。

感谢大家！

作者
Author: sobereva 时间: 2024-2-19 00:47
可以看Multiwfn源代码orbloc.f90，在里面搜!Carry out localization，从这开始到PM定域化结束也就150行左右，读起来不费劲，注释也清楚

作者
Author: Freeman 时间: 2024-2-19 13:50
本帖最后由 Freeman 于 2024-2-19 13:51 编辑

sobereva 发表于 2024-2-19 00:47
可以看Multiwfn源代码orbloc.f90，在里面搜!Carry out localization，从这开始到PM定域化结束也就150行左右 ...

谢谢。我看Multiwfn的代码，纠正了一个错误，结果就对了。前面5楼的步骤里，不需要找到ER eq26最大的轨道对，而是要对所有轨道对做第二步。至于为什么PM和ER局域化的步骤不一样，我也不知道，不管了。

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