PM局域化优化到局部最优怎么办？

Freeman · 发表于 Post on 2024-2-16 11:37:56

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

42_loc.mwfn (2.33 MB, 下载次数 Times of downloads: 0)

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

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

42_loc_new.mwfn (2.33 MB, 下载次数 Times of downloads: 0)

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

hebrewsnabla · 发表于 Post on 2024-2-16 12:22:47

pyscf的实现 https://github.com/pyscf/pyscf/pull/1990

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

sobereva · 发表于 Post on 2024-2-16 16:27:08

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）

zjxitcc · 发表于 Post on 2024-2-16 21:36:01

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

Freeman · 发表于 Post on 2024-2-18 18:08:43

本帖最后由 Freeman 于 2024-2-18 18:11 编辑

@hebrewsnabla @sobereva @zjxitcc

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

有关2*2旋转实现pm局域化的，我找到两篇文献，分别是ER（

RevModPhys.35.457.pdf (809.55 KB, 下载次数 Times of downloads: 1) ）和PM（

1.456588.pdf (1.08 MB, 下载次数 Times of downloads: 2) ）的原文。ER文章给出了2*2旋转的大致算法，而PM文章给出了一些2*2用于PM局域化的相关参数公式（Ast，Bst）。还有之前参考的梯度优化的文献SUSI（

ct400793q.pdf (1.08 MB, 下载次数 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步直到轨道几乎不再变化。

感谢大家！

sobereva · 发表于 Post on 2024-2-19 00:47:08

可以看Multiwfn源代码orbloc.f90，在里面搜!Carry out localization，从这开始到PM定域化结束也就150行左右，读起来不费劲，注释也清楚

Freeman · 发表于 Post on 2024-2-19 13:50:04

本帖最后由 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局域化的步骤不一样，我也不知道，不管了。

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[算法与编程] PM局域化优化到局部最优怎么办？

评分 Rate

评分 Rate

评分 Rate