请教一个DFT波函数二阶优化的问题

Freeman

大家好。

TRAH-SCF的文章（https://doi.org/10.1063/5.0040798）的公式A10给出了交换相关势矩阵Vxc沿着轨道旋转系数kappa变化的方向导数（大概是这样吧？）。我对它有两点疑问。

一、这个公式要对格点求二次积分（dr和dr'），岂不是会非常非常非常慢吗？但是依我平时的经验，每一圈TRAH-SCF的micro-step的耗时和一圈普通SCF的耗时，差不多。

二、exc对rho(r)和rho(r')的二阶变分是怎么求出来的。我知道各种交换相关库可以算exc对rho的高阶导数，但分母的rho的自变量应该都是r，而不能一个是r另一个是r'。另外，公式里分子上的exc的自变量是r还是r'呢？

感谢解答。

hebrewsnabla

这是一种比较原始的写法，Exc是总的xc能量泛函，它对rho(r),rho(r')的二阶导数最终会转换成泛函核对rho(r)的二阶导数。

Freeman

hebrewsnabla 发表于 2025-1-10 15:41
这是一种比较原始的写法，Exc是总的xc能量泛函，它对rho(r),rho(r')的二阶导数最终会转换成泛函核对rho(r) ...

恐怕不是这样。同一篇文章的公式A8就是普通的LDA交换相关势矩阵的计算公式，里面也出现了exc。既然exc出现在那个位置，就说明exc的意思是格点上的单粒子交换相关能量，而不是总的交换相关能。另外，如果要表达总的交换相关能，一般都写作大写的Exc而不是小写的exc。

hebrewsnabla

Freeman 发表于 2025-1-10 18:00
恐怕不是这样。同一篇文章的公式A8就是普通的LDA交换相关势矩阵的计算公式，里面也出现了exc。既然exc出 ...

啊，这个影响不大，exc对rho(r),rho(r')的二阶导数最终也是转换成泛函核对rho(r)的二阶导数。

Freeman

hebrewsnabla 发表于 2025-1-10 19:26
啊，这个影响不大，exc对rho(r),rho(r')的二阶导数最终也是转换成泛函核对rho(r)的二阶导数。

原来如此。请问你有啥相关资料吗？我最近也在实现一个DFT的二阶优化，但苦于找不到资料，不知道交换相关势的导数该咋写。TRAH文章的附录是我见过最详细的了，但还不够让我整明白。自己推又太坐牢了。

hebrewsnabla

Freeman 发表于 2025-1-10 19:33
原来如此。请问你有啥相关资料吗？我最近也在实现一个DFT的二阶优化，但苦于找不到资料，不知道交换相关 ...

我觉得他这个写法很不好，完全不是working equation。我不太使用交换相关势的导数这种说法，感觉orbital hessian更常用一些。其中的核心部分fxc matrix可以参考 https://github.com/pyscf/pyscf/b ... dft/numint.py#L1418

fxc这个概念本身是通用的，所以早期tddft文献推过之后可能后面的soscf、stability文献就不推了，或者并不专门讲dft部分。而tddft文献里面是一定要处理fxc的。

Freeman

hebrewsnabla 发表于 2025-1-10 20:01
我觉得他这个写法很不好，完全不是working equation。我不太使用交换相关势的导数这种说法，感觉orbital  ...

谢谢提醒。我之前不知道还可以从stability和tddft方面的文献里找资料。不过即使是这些文献，也没有明确地告诉别人哪个公式是orbital hessian。我只能通过一些线索，猜测哪些是我需要的量。

我首先找到了这篇文献A（https://sci-hub.se/10.1063/1.5044202）。A有个两个公式。（2.23）是Fock矩阵关于密度矩阵的导数，间接就可以认为是orbital hessian了。很明显，（2.23）前两项是库伦和HF交换项，第三项是XC项。（2.24）是XC项的具体形式，正好和TRAH文献写的一样。


第二步是要找（2.24）的working equation（我对“working equation”的理解是“不需要弯弯绕、照着抄就能写出代码的方程”）。我找到了这篇文献B（https://sci-hub.se/10.1016/0009-2614(96)00440-x）。B有两块地方对我很重要，其中一个是这里。

在（12）里，又出现了上面（2.24）的“Exc的对两个自变量不同的电子密度的变分”，而（15）告诉我这个量（以及HF交换项）存放在了M张量里。

另一个重要的地方是这里。

这里终于给出了M的具体形式。除了第一行是HF交换项，剩下的就是我想要的。

不过我到这里还是有一些疑问：

一、我该用（24）（25）中的哪一个做闭壳层二阶SCF呢？
B是TDDFT最早期的文献之一。B指出，如果你的参考波函数是闭壳层的，那么你就可以只算单重态的激发态（此时用24）或只算三重态的激发态（此时用25）。这是不是意味着，如果我的波函数是闭壳层，我就该用（24）来算orbital hessian呢？只有当我想检查波函数是否有RHF->UHF不稳定性时，或者当我要做RO三重态的二阶优化，才需要（25）？

二、如果我的波函数是非限制自旋（unrestricted）的，该怎么办？文献B并没有给出UHF对应的M张量。

hebrewsnabla

1. 见下一楼。

2. uks就不区分什么 singlet fxc, triplet fxc了，那么soscf的fxc也是一样的。

关于几种 orbital hessian的理解，或许 10.1063/1.434318 有一定帮助，虽然这是hf的。

最后，需要提醒一下，一般来说我们不直接算 orbital hessian, 而是算 H_ia,jb和一个向量的积。

hebrewsnabla

https://doi.org/10.1063/1.471637 appendix 里面的A应该就是soscf和rks internal stability的orbital hessian。这个名字有点迷惑性，实际上相当于tddft的1A+1B，可以与你的文献B对照（L，M相当于A，B）。

Freeman

hebrewsnabla 发表于 2025-1-11 18:34
https://doi.org/10.1063/1.471637 appendix 里面的A应该就是soscf和rks internal stability的orbital hess ...

谢谢。在我正在研究的scheme里，就是需要Fock矩阵对密度矩阵的导数，而不是能量对轨道旋转系数的导数。（当然了，知道了其中一个，也能推导出另一个；只是如果有现成的前者，对我而言比较省事儿。）因此，我上面发的文献A的第一张图的\partial F / \partial P的形式其实正是我想要的。我猜，文献B大概也是\partial F / \partial P这个意思吧。你发的文献没直接说Appendix的A和B究竟是谁关于谁的导数，不过该文献里涉及了轨道旋转系数，比如公式13、14、21-25。这就让我有些担心，A和B好像其实是能量对轨道旋转系数的导数。

hebrewsnabla

Freeman 发表于 2025-1-11 19:56
谢谢。在我正在研究的scheme里，就是需要Fock矩阵对密度矩阵的导数，而不是能量对轨道旋转系数的导数。（ ...

文献B和https://doi.org/10.1063/1.471637 是同一个人写的，都是对轨道旋转系数的导数。

Freeman

hebrewsnabla 发表于 2025-1-11 20:03
文献B和https://doi.org/10.1063/1.471637 是同一个人写的，都是对轨道旋转系数的导数。

真是个悲伤的故事。不过我有点儿预感，对轨道旋转系数的导数可能跟对密度矩阵的导数是一样的。譬如，能量对轨道旋转系数的一阶导数是MO基的Fock矩阵，而能量对正交化AO基的密度矩阵的一阶导数则是正交化AO基的Fock矩阵，两者都是Fock矩阵，只是差个基变换。二阶导数会不会也有类似的规律？我打算先直接试试。谢谢啦。

wzkchem5

Freeman 发表于 2025-1-11 15:56
真是个悲伤的故事。不过我有点儿预感，对轨道旋转系数的导数可能跟对密度矩阵的导数是一样的。譬如，能量 ...

对，到二阶导为止都有明确的对应关系

Warm_Cloud

先回答你1楼的第一个问题，你发现了很多文献中的Exc的二阶变分项包含对dr和dr'的积分，但是实际计算的时候只有对dr的积分（你在7楼贴的图片），也就是积分从6维成了3维，原因参考https://www.annualreviews.org/co ... em.55.091602.094449中的(39)式。

另外从编程角度告诉你求所需矩阵元的方法：一阶：dE/dD_{uv}（Fock矩阵元），二阶：d2E/dD_{uv}^2（TDDFT和orbital hessian里面的fxc部分），三阶：d3E/dD_{uv}^3（TDDFT一阶导数中的gxc部分），D_{uv}是密度矩阵的矩阵元。举个例子：你要求GGA部分的fxc，那就是dExc[rho,|∇rho|^2]/dD_{uv}=∂Exc[rho,|∇rho|^2]/∂rho*∂rho/∂dD_{uv}+∂Exc[rho,|∇rho|^2]/∂|∇rho|^2*∂|∇rho|^2/∂dD_{uv}
其中∂Exc[rho,|∇rho|^2]/∂rho和∂Exc[rho,|∇rho|^2]/∂|∇rho|^2就是泛函库给你的结果，∂rho/∂dD_{uv}=u(r)v(r)就是两个基函数的基了，∂|∇rho|^2/∂dD_{uv}你可根据我前面说的自己推一下。meta-GGA道理一样。

Warm_Cloud

补充一下：
TDDFT的工作方程：
|A B|(X)=omega|1  0|(X)
|B A|(Y)              |0 -1|(Y)
TDA的：
AX=omegaX
波函数稳定性检测（也就是对角化orbital hessian）：
(A+B)X=omegaX

其中三个情况的A,B一样，区别是前两个还有三重态的情况，第三个只用到单重态的形式。
单重态：∂/∂a和∂/∂b，三重态：∂/∂a和∂/∂(-b)，利用这个规则你可以得到7楼的公式。