计算化学公社

标题: 请教一个DFT波函数二阶优化的问题 [打印本页]

作者
Author: Freeman 时间: 2025-1-10 14:28
标题: 请教一个DFT波函数二阶优化的问题
大家好。

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'呢？

感谢解答。

作者
Author: hebrewsnabla 时间: 2025-1-10 15:41
本帖最后由 hebrewsnabla 于 2025-1-10 16:06 编辑

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

作者
Author: Freeman 时间: 2025-1-10 18:00

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

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

作者
Author: hebrewsnabla 时间: 2025-1-10 19:26

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

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

作者
Author: Freeman 时间: 2025-1-10 19:33

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

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

作者
Author: hebrewsnabla 时间: 2025-1-10 20:01
本帖最后由 hebrewsnabla 于 2025-1-10 20:06 编辑

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的。

作者
Author: Freeman 时间: 2025-1-11 02:12
本帖最后由 Freeman 于 2025-1-11 02:17 编辑

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文献写的一样。
(, 下载次数 Times of downloads: 23)

第二步是要找（2.24）的working equation（我对“working equation”的理解是“不需要弯弯绕、照着抄就能写出代码的方程”）。我找到了这篇文献B（https://sci-hub.se/10.1016/0009-2614(96)00440-x）。B有两块地方对我很重要，其中一个是这里。
(, 下载次数 Times of downloads: 28)
在（12）里，又出现了上面（2.24）的“Exc的对两个自变量不同的电子密度的变分”，而（15）告诉我这个量（以及HF交换项）存放在了M张量里。

另一个重要的地方是这里。
(, 下载次数 Times of downloads: 26)
这里终于给出了M的具体形式。除了第一行是HF交换项，剩下的就是我想要的。

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

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

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

作者
Author: hebrewsnabla 时间: 2025-1-11 17:09
本帖最后由 hebrewsnabla 于 2025-1-11 18:47 编辑

1. 见下一楼。

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

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

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

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

作者
Author: Freeman 时间: 2025-1-11 19:56

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好像其实是能量对轨道旋转系数的导数。

作者
Author: hebrewsnabla 时间: 2025-1-11 20:03

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

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

作者
Author: Freeman 时间: 2025-1-11 22:56

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

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

作者
Author: wzkchem5 时间: 2025-1-12 04:18

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

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

作者
Author: Warm_Cloud 时间: 2025-1-13 15:34
先回答你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道理一样。

作者
Author: Warm_Cloud 时间: 2025-1-13 15:43
补充一下：
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楼的公式。

作者
Author: Warm_Cloud 时间: 2025-1-13 21:36
第一个问题参考Robert G Parr Weitao Yang的书的附录A，见附件。

作者
Author: wzkchem5 时间: 2025-1-13 21:45

Warm_Cloud 发表于 2025-1-13 08:43
补充一下：
TDDFT的工作方程：
|A B|(X)=omega|1 0|(X)

做UHF instability检测的时候，是需要用到A+B的三重态形式的。RHF instability要用A+B的单重态形式，complex instability用到A-B

作者
Author: hebrewsnabla 时间: 2025-1-13 22:21
本帖最后由 hebrewsnabla 于 2025-1-13 22:27 编辑

wzkchem5 发表于 2025-1-13 21:45
做UHF instability检测的时候，是需要用到A+B的三重态形式的。RHF instability要用A+B的单重态形式，comp ...

rhf to uhf instability，或者叫rhf external 是 triplet A+B. uhf internal 并不用 triplet，和uhf soscf是一样的。参见 doi:10.1063/1.434318.

作者
Author: wzkchem5 时间: 2025-1-13 22:27

hebrewsnabla 发表于 2025-1-13 15:21
rhf to uhf instability，或者叫rhf external 是 triplet A+B. uhf internal 并不用 triplet。参见 doi:1 ...

对，我就是想说RHF to UHF instability。
UHF internal的情况应该说是没有singlet、triplet的分别，也就是说用的A+B既不能说是triplet的，也不能说是singlet的，它的变分空间包括了RHF情况下的singlet和triplet的变分空间

作者
Author: Warm_Cloud 时间: 2025-1-14 09:10

wzkchem5 发表于 2025-1-13 21:45
做UHF instability检测的时候，是需要用到A+B的三重态形式的。RHF instability要用A+B的单重态形式，comp ...

对，三重态的忘了，用R做U是需要的。

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