对社长博文《密度泛函计算中的格点积分方法》的一点疑问

Freeman

大家好。最近我在研读社长的《密度泛函计算中的格点积分方法》，有如下疑问。
单中心球极坐标积分∫W_i(r)*F(r)dr，其中W_i(r)是Becke划分下原子i对于格点r的空间权重，F是被积函数。由中心i生发出来的格点，其径向和角度权重可以由高斯积分多项式和lebedev方法给出；但是因为积分范围是所有格点，而且Becke划分是模糊的（fuzzy），积分范围内不可避免还有别的中心（例如j）生发出来的格点。后者肯定不会正好落在本中心的某个径向或角度节点上，那么它们的径向和角度权重怎么知道呢？

hebrewsnabla

先生成每个原子的格点及其权重，然后合并，合并后的权重是原来的权重和becke划分的权重的乘积。

Freeman

hebrewsnabla 发表于 2023-8-16 12:31
先生成每个原子的格点及其权重，然后合并，合并后的权重是原来的权重和becke划分的权重的乘积。

谢谢回复。您能细说一下“合并”吗？我的理解是，中心i给出的径向和角度的权重只是以中心i做积分时所使用的。如果中心i生成的某个格点落到相邻的中心j范围去了，那么以中心j做积分时，这个格点的径向和角度的权重是多少呢？应该不能再用中心i给出的权重了吧？

hebrewsnabla

Freeman 发表于 2023-8-16 12:52
谢谢回复。您能细说一下“合并”吗？我的理解是，中心i给出的径向和角度的权重只是以中心i做积分时所使 ...

直接合并。因为becke划分是fuzzy的，所以并不存在“中心j的范围”。径向和角度权重本来是多少还是多少。只有不fuzzy的划分才存在某个原子的范围。

可以看看sob博文里的代码，或者看这个 https://github.com/pyscf/pyscf/blob/master/pyscf/dft/gen_grid.py

Freeman

好的，那就忽略“范围”这几个字。
这么说吧，以中心i做积分I_i=∫W_i(r)*F(r)dr=∑[(r,Ω) within all grids]W_i(r,Ω)*F(r,Ω)*w_i_r(r)*w_i_a(Ω)
其中w_i_r(r)和w_i_a(Ω)是中心i给出的该格点的径向和角度权重。
您说“径向和角度权重本来是多少还是多少”，意思就是任何一个格点(r,Ω)，对于所有中心i，w_i_r和w_i_a都相等，于是
I_i=∑[(r,Ω) within all grids]W_i(r,Ω)*F(r,Ω)*w_r(r)*w_a(Ω)
除了第一项W_i(r,Ω)，被积函数的其他部分都和i没关系了。
这会带来一个问题：已知∑[i within all centres]W_i(r,Ω)=1，那么
∑I_i=∑[(r,Ω) within all grids]{∑[i within all centres]W_i(r,Ω)}*F(r,Ω)*w_r(r)*w_a(Ω)=∑[r within all grids]F(r,Ω)*w_r(r)*w_a(Ω)
这个结果和W_i(r,Ω)无关了，说明无论我们怎么划分原子，结果都是相等的。这不可能对。
所以我想，大概对于所有中心i，w_i_r(r,Ω)和w_i_a(r,Ω)是不相等的，这样才能让不同的原子区域导致不同的结果。对于中心i给出的格点，w_i_r(r,Ω)和w_i_a(r,Ω)很好计算，可以由高斯积分多项式和lebedev给出；可是其他中心j给出的格点呢？

Freeman

Freeman 发表于 2023-8-16 16:32
好的，那就忽略“范围”这几个字。
这么说吧，以中心i做积分I_i=∫W_i(r)*F(r)dr=∑[(r,Ω) within all gr ...

我想了很久，想到了一个路子，可以让结果取决于划分方法，希望是正常dft处理格点积分的方式。
中心i的积分I_i=∑[(r,Ω) within all grids]W_i(r,Ω)*F(r,Ω)*w_r(r)*w_a(Ω)的求和范围应该不是所有格点，而是从本中心生发出来的格点。对于中心i的某个格点(r,Ω)，既要算W_i(r,Ω)，又要算W_j(r,Ω)[j!=i]，但是W_j(r,Ω)单纯只是为了求和，将W_i(r,Ω)归一化用的，并不显式地出现在被积函数中。

hebrewsnabla

意思就是任何一个格点(r,Ω)，对于所有中心i，w_i_r和w_i_a都相等，

首先求一个全空间积分并不是 sum_k sum_i (假设k是合并后格点的编号），而是sum_i sum_{k 属于 i}, 或者换句话说只要 sum_k就完了，所以一个格点k如果在生成原子格点时是属于原子i的话，他只有w_i(k)，而其他的w_j(k)是不存在的。

卡开发发

Freeman 发表于 2023-8-16 16:51
我想了很久，想到了一个路子，可以让结果取决于划分方法，希望是正常dft处理格点积分的方式。
中心i的积 ...

1、你这么理解好了：
初始给定
（1）积分I=∫F(r)dr
（2）引入划分F_i(r)=F(r)*w_i(r)
（3）划分函数∑_i w_i(r)=1
在（1）积分插入1用（3）展开，得
I=∫∑_i w_i(r)F(r)dr=∑_i ∫w_i(r)F(r)dr，若定义I_i=∫w_i(r)F(r)dr，则I=∑_i I_i。
如果我们这时候只看I_i=∫w_i(r)F(r)这个关系，可以扔掉脚标，然后再把w_i(r)F(r)看成一个整体f_i(r)，即
I_i=∫f_i(r)dr，这时候我们单独选一个以f_i(r)中心为参考的网格对f_i(r)进行离散，并加权求和求I_i，此时和f_j(r)完全没啥关系，最后单独求和I=∑_i I_i。注：如果F(r)本身是先由g_i(r)叠加出来的，经过划分f_i(r)会有g_j(r)的贡献，但并不影响上述求积的形式，权重是随着分割f_i(r)走的。
2、

说明无论我们怎么划分原子，结果都是相等的。这不可能对。

（虽然我不确定是不是完全充要的条件）保证∑_i w_i(r)=1，其实任何划分都是可以的，之所以积分I的微小差异是因为不同划分方法可能对径向网格或角度lm收敛快慢不同所导致的，这是数值问题。对于单独的I_i，因为划分后f_i(r)不一样，各个区域贡献大小肯定也会不同，另外就是不同区域覆盖的原子不同，一般为了节约计算量，对不同原子选取的径向网格和角度网格数目有可能也会不同。Becke划分或Hirshfeld划分或者其他的方式都是可以的（例如DMol3当中的划分为了局域会选择w_A=(ρ_A/r_A)^2这种形式，参考此处11~14页），他们收敛的极限应该都相同。也就是说，原理上也允许你把空间划分成一大堆方块，每个方块上比方说也可以均匀离散相关区域的被积函数，然后各自进行求积，总的结果按道理不会不同，但是一般为什么不这么做，你仔细思考应该就能知道这个划分的意义到底是什么，反过来如果不做划分会怎样？

sobereva

全空间积分的时候依次循环各个原子（中心），将每个中心的积分结果加和
对每个中心积分时，循环由它向四周发散出来的点。对这些点，只需要计算：
(1)被积函数在这些点上的数值
(2)当前这个中心在这些点上的权重函数值（Becke方法）
(3)单中心积分权重（径向&角度）在这些点上的数值
对每个点的以上三项相乘，循环加和各个点即可。

积分当前中心格点时不用管其它中心。计算(2)的时候虽然也会顺带产生其它中心在这些点上的权重函数值，但不会利用到。也不用考虑什么其它中心的格点落入到当前中心的格点分布范围内的问题，因为那些中心由其Becke权重函数模糊式截断出来的范围函数，靠它们自己的单中心积分就已经很好地积分了。积分当前中心的时候也不会积分到该属于其它原子的部分（因为权重函数就只利用了当前中心的）。