DFT 平面波方法
变分法
考虑薛定谔方程，
H
^
ψ
=
E
ψ
\hat{H} \psi=E \psi
H^ψ=Eψ
这里的波函数我们一般不能精确得到，所以我们一般找一个数学上可以处理的函数去逼近，
ϕ
≈
ψ
\phi \approx \psi
ϕ≈ψ
特征值理论告诉我们，极小化能量泛函
E
~
=
⟨
ϕ
∣
H
^
∣
ϕ
⟩
⟨
ϕ
∣
ϕ
⟩
\tilde{E}=\frac{\langle\phi|\hat{H}| \phi\rangle}{\langle\phi \mid \phi\rangle}
E~=⟨ϕ∣ϕ⟩⟨ϕ∣H^∣ϕ⟩
便可以得到最小的特征值
E
0
E_0
E0（基态能量），此时的
ϕ
\phi
ϕ 也就是对应的特征函数。
如果我们把逼近波函数
ϕ
\phi
ϕ 找一组基底展开，
ϕ
(
x
⃗
)
=
∑
i
=
1
N
c
i
χ
i
(
x
⃗
)
\phi(\vec{x})=\sum_{i=1}^{N} c_{i} \chi_{i}(\vec{x})
ϕ(x)=i=1∑Nciχi(x)
在毕竟波函数单位化约束
⟨
ϕ
∣
ϕ
⟩
=
1
\langle\phi \mid \phi\rangle=1
⟨ϕ∣ϕ⟩=1 的条件下，我们使用拉格朗日乘子法，可以得到广义代数特征值问题，
H
⋅
C
=
λ
⋅
S
⋅
C
\mathbf{H} \cdot \mathbf{C}=\lambda \cdot \mathbf{S} \cdot \mathbf{C}
H⋅C=λ⋅S⋅C
其中，
H
n
m
=
⟨
χ
n
∣
H
^
∣
χ
m
⟩
S
n
m
=
⟨
χ
n
∣
χ
m
⟩
\begin{aligned} H_{n m} &=\left\langle\chi_{n}|\hat{H}| \chi_{m}\right\rangle \\ S_{n m} &=\left\langle\chi_{n} \mid \chi_{m}\right\rangle \end{aligned}
HnmSnm=⟨χn∣H^∣χm⟩=⟨χn∣χm⟩
类比有限元方法，其中的
S
n
m
S_{nm}
Snm 是质量矩阵，
H
n
m
H_{nm}
Hnm 可以看成是瑞利泛函意义下的“刚度”矩阵。
DFT 平面波方法
求解 KS 方程
[
T
^
s
+
V
e
f
f
]
ϕ
i
=
ϵ
i
ϕ
i
\left[\hat{T}_{s}+V_{e f f}\right] \phi_{i}=\epsilon_{i} \phi_{i}
[T^s+Veff]ϕi=ϵiϕi
对波函数进行平面波展开，并对有效势做傅里叶展开，和有限元方法类似，使用平面波基底做测试函数进行变分，忽略一堆推导，最后可以得到一个代数特征值问题，
∑
m
H
m
′
m
(
k
⃗
)
c
i
,
m
(
k
⃗
)
=
ϵ
i
(
k
⃗
)
c
i
,
m
′
(
k
⃗
)
\sum_{m} H_{m^{\prime} m}(\vec{k}) c_{i, m}(\vec{k})=\epsilon_{i}(\vec{k}) c_{i, m^{\prime}}(\vec{k})
m∑Hm′m(k)ci,m(k)=ϵi(k)ci,m′(k)
其中，
H
m
′
m
(
k
⃗
)
=
1
2
∣
k
⃗
+
G
⃗
m
∣
2
δ
m
′
m
+
V
e
f
f
(
G
⃗
m
−
G
⃗
m
′
)
H_{m^{\prime} m}(\vec{k})=\frac{1}{2}\left|\vec{k}+\vec{G}_{m}\right|^{2} \delta_{m^{\prime} m}+V_{e f f}\left(\vec{G}_{m}-\vec{G}_{m^{\prime}}\right)
Hm′m(k)=21∣∣∣k+Gm∣∣∣2δm′m+Veff(Gm−Gm′)
第一部分是动能项，比较好理解。第二部分可以进一步细化。略去推导，我们直接给出表达式。
- Hartree 势
V H ( G ⃗ ) = 4 π n ( G ⃗ ) G 2 V_{H}(\vec{G})=4 \pi \frac{n(\vec{G})}{G^{2}} VH(G)=4πG2n(G)
其中的 n n n 表示电子密度。 - 交换关联势
V x c ( G ⃗ ) = ∑ G ⃗ ′ n x c ( G ⃗ − G ⃗ ′ ) d ϵ x c d n ( G ⃗ ′ ) + ϵ x c ( G ⃗ ) V_{x c}(\vec{G})=\sum_{\vec{G}^{\prime}} n_{x c}\left(\vec{G}-\vec{G}^{\prime}\right) \frac{d \epsilon_{x c}}{d n}\left(\vec{G}^{\prime}\right)+\epsilon_{x c}(\vec{G}) Vxc(G)=G′∑nxc(G−G′)dndϵxc(G′)+ϵxc(G) - 外势（球谐）
V κ ( G ⃗ ) = V κ ( ∣ G ⃗ ∣ ) ∫ 0 ∞ d r r 2 j 0 ( ∣ G ⃗ ∣ r ) V κ ( r ) V^{\kappa}(\vec{G})=V^{\kappa}(|\vec{G}|) \int_{0}^{\infty} d r r^{2} j_{0}(|\vec{G}| r) V^{\kappa}(r) Vκ(G)=Vκ(∣G∣)∫0∞drr2j0(∣G∣r)Vκ(r)
对于平面波基底，电子密度可以表达为
n
i
k
⃗
(
G
⃗
)
=
1
Ω
∑
m
c
i
,
m
∗
(
k
⃗
)
c
i
,
m
′
′
(
k
⃗
)
n_{i \vec{k}}(\vec{G})=\frac{1}{\Omega} \sum_{m} c_{i, m}^{*}(\vec{k}) c_{i, m^{\prime \prime}}(\vec{k})
nik(G)=Ω1m∑ci,m∗(k)ci,m′′(k)
其中，
m
′
′
m^{\prime \prime}
m′′ 是向量
G
⃗
m
′
′
=
G
⃗
m
+
G
⃗
\vec{G}_{m^{\prime \prime}}=\vec{G}_{m}+\vec{G}
Gm′′=Gm+G 的指标。
截断能的选择拌住了倒格矢量，
1
2
∣
G
⃗
∣
2
<
E
c
u
t
\frac{1}{2}|\vec{G}|^{2}<E_{c u t}
21∣G∣2<Ecut
DFT 平面波迭代解法
一般来说，按上述的方式，以 m ′ m' m′ 为行标， m m m 为列标，去组装 H H H 矩阵，最后在 SCF 迭代过程中求解代数特征值问题即可。但是，这样做有个问题，即如果我们要达到相对较高的精度，就必须选取足够个数的基底，这使得矩阵规模很大，存不下，也算不动。事实上，我们需要的只是前几个特征值和特征向量，这时候的对角化求特征值特征向量，并不需要把整个矩阵组装出来。所用的方法就是子空间迭代，即先给一个特征向量初猜，慢慢地产生一系列向量，最后张成向量空间。比如说，我们可以用 CG 方法。
KS 哈密顿量的对角化是是任何电子结构系统 scf 求解的主要步骤。CG 方法的过程如下：
- 对每个能带，给定一个测试特征对： { ∣ ϕ i ( n ) ⟩ , ε i } \left\{\left|\phi_{i}^{(n)}\right\rangle, \varepsilon_{i}\right\} {∣∣∣ϕi(n)⟩,εi}
- 在约束条件
⟨ ϕ i ∣ S ∣ ϕ j ⟩ = δ i j , ∀ j ≤ i \left\langle\phi_{i}|S| \phi_{j}\right\rangle=\delta_{i j}, \quad \forall j \leq i ⟨ϕi∣S∣ϕj⟩=δij,∀j≤i
下，使用预条件 CG 方法极小化单粒子能量
E ( ∣ ϕ i ⟩ ) = ⟨ ϕ i ∣ H K S ∣ ϕ j ⟩ E\left(\left|\phi_{i}\right\rangle\right)=\left\langle\phi_{i}\left|H_{K S}\right| \phi_{j}\right\rangle E(∣ϕi⟩)=⟨ϕi∣HKS∣ϕj⟩ - 重复上述过程，直到完成。
abacus 中采用的是 band-by-band 的 CG 方法。
预条件 CG 也可以采用 PPCG：
abacus 源码阅读
我们把上述的迭代解法对应到 abacus 源码当中。
void Electrons::c_bands(const int &istep)
c_bands 就是计算哈密顿量并对角化的一个函数。正常情况下，构建哈密顿量就是建立一个 H 矩阵，然后把里面的矩阵元算出来，最后把 H 放进对角化函数，哈密顿量很大，不能全部储存。而且我们并不需要对角化得到所有本征值，只要最低的几个就可以了。所以，程序中使用 CG 迭代求解。CG 算法只需要能量泛函，所以单独储存的哈密顿量是不需要的，只要储存一个向量 ∣ ψ j ⟩ = H ∣ ϕ j ⟩ \left|\psi_{j}\right\rangle=H\left|\phi_{j}\right\rangle ∣ψj⟩=H∣ϕj⟩。对于动能项 g2kin ，因为平面波基底下下动能算符是对角矩阵，所以只要存一个 npw 维度的向量。我们把 c_bands 打开，对于其中的每个 k 点操作，进行阅读。
for (int ig = 0; ig < GlobalC::wf.npw; ig++)
{
h_diag[ig]
= 1 + GlobalC::wf.g2kin[ig] + sqrt(1 + (GlobalC::wf.g2kin[ig] - 1) * (GlobalC::wf.g2kin[ig] - 1));
if (GlobalV::NPOL == 2)
h_diag[ig + GlobalC::wf.npwx] = h_diag[ig];
}
这一步是构建 CG 算法所需要的预条件矩阵 h_diag。
GlobalC::hm.diagH_pw(istep, this->iter, ik, h_diag, avg_iter_k);
这一步是 CG 算法的核心，我们可以把它打开。其中，它核心的循环体是
if ( iter > 1 || istep > 1 || ntry > 0)//iter will grow up, when iter==2, carry out this branch
{
//cout<<"ls to here"<<endl;exit(0);//to here when iter==2, iter==1 at beginning
//cout<<"ls out:"<<"ik="<<ik0<<" ik="<<ik<<endl;exit(0);
this->diagH_subspace(
ik,
GlobalV::NBANDS,
GlobalV::NBANDS,
GlobalC::wf.evc[ik0],//0
GlobalC::wf.evc[ik0],//0
GlobalC::wf.ekb[ik]);//0
avg_iter += 1.0;
}
当迭代次数不是第一次的时候，需要先执行一次 diagH_subspace 的对角化，再执行 CG 的对角化。
if(GlobalV::NPOL==1)
{
//cout<<"ls to here"<<endl;exit(0);//to here
cg.diag(GlobalC::wf.evc[ik0], GlobalC::wf.ekb[ik], GlobalC::kv.ngk[ik], GlobalC::wf.npwx,
GlobalV::NBANDS, precondition, GlobalV::PW_DIAG_THR,
GlobalV::PW_DIAG_NMAX, reorder, notconv, avg);
}
这一步是 CG 对角化的核心，采用 ccgdiagg 算法，我们把它打开。
void Diago_CG::diag
(
ModuleBase::ComplexMatrix &phi,
double *e,
const int &dim,
const int &dmx,
const int &n_band,
const double *precondition,
const double &eps,
const int &maxter,
const bool &reorder,
int & notconv,
double &avg_iter
)
这里的 phi 是个矩阵，每行表示 了 ∣ ψ j ⟩ = H ∣ ϕ j ⟩ \left|\psi_{j}\right\rangle=H\left|\phi_{j}\right\rangle ∣ψj⟩=H∣ϕj⟩，e 是这些特征向量对应的特征值。采用的是预条件的共轭梯度算法，Band-by-band 算法最小化内存的使用。h_1psi 和 s_1psi分别计算 H|phi> 和 S|phi>。我们将其打开，看对每个特征对的获取操作。
void Hamilt_PW::s_1psi
(
const int dim,
const std::complex<double> *psi,
std::complex<double> *spsi
)
{
for (int i=0; i<dim; i++)
{
spsi[i] = psi[i];
}
return;
}
考虑 S 是个单位矩阵，这相当于就做了一个赋值操作。
void Diago_CG::schmit_orth
(
const int& dim,
const int& dmx,
const int& m, //end
const ModuleBase::ComplexMatrix &psi,
std::complex<double> *sphi,
std::complex<double> *psi_m
)
正交化操作，使得 psi_m 正交于 psi(start) ~ psi(m-1)，即
ψ
(
m
)
=
ψ
(
m
)
−
∑
i
<
m
<
ψ
(
i
)
∣
S
∣
ψ
(
m
)
>
ψ
(
i
)
\psi(m) = \psi(m) - \sum_{i < m} < \psi(i) | S | \psi(m) >\psi(i)
ψ(m)=ψ(m)−i<m∑<ψ(i)∣S∣ψ(m)>ψ(i)
void Hamilt_PW::h_1psi( const int npw_in, const std::complex < double> *psi,
std::complex<double> *hpsi, std::complex < double> *spsi)
{
this->h_psi(psi, hpsi);
for (int i=0;i<npw_in;i++)
{
spsi[i] = psi[i];
}
return;
}
计算了 hpsi 和 spsi。hpsi 的计算使用，
void Hamilt_PW::h_psi(const std::complex<double> *psi_in, std::complex<double> *hpsi, const int m)
其中包含了动能的计算：
//------------------------------------
//(1) the kinetical energy.
//------------------------------------
std::complex<double> *tmhpsi;
const std::complex<double> *tmpsi_in;
if(GlobalV::T_IN_H)
{
tmhpsi = hpsi;
tmpsi_in = psi_in;
for(int ib = 0 ; ib < m; ++ib)
{
for(ig = 0;ig < GlobalC::wf.npw; ++ig)
{
tmhpsi[ig] = GlobalC::wf.g2kin[ig] * tmpsi_in[ig];
}
if(GlobalV::NSPIN==4){
for(ig=GlobalC::wf.npw; ig < GlobalC::wf.npwx; ++ig)
{
tmhpsi[ig] = 0;
}
tmhpsi +=GlobalC::wf.npwx;
tmpsi_in += GlobalC::wf.npwx;
for (ig = 0;ig < GlobalC::wf.npw ;++ig)
{
tmhpsi[ig] = GlobalC::wf.g2kin[ig] * tmpsi_in[ig];
}
for(ig=GlobalC::wf.npw; ig < GlobalC::wf.npwx; ++ig)
{
tmhpsi[ig] =0;
}
}
tmhpsi += GlobalC::wf.npwx;
tmpsi_in += GlobalC::wf.npwx;
}
}
包含动能矩阵乘以 phi。phi 被拉成一条向量，每一段都乘以动能矩阵，以相同的方式放到 tmhphi 中。
局部势能的计算：
//------------------------------------
//(2) the local potential.
//-----------------------------------
ModuleBase::timer::tick("Hamilt_PW","vloc");
if(GlobalV::VL_IN_H)
{
tmhpsi = hpsi;
tmpsi_in = psi_in;
for(int ib = 0 ; ib < m; ++ib)
{
if(GlobalV::NSPIN!=4){
ModuleBase::GlobalFunc::ZEROS( GlobalC::UFFT.porter, GlobalC::pw.nrxx);
GlobalC::UFFT.RoundTrip( tmpsi_in, GlobalC::pot.vr_eff1, GR_index, GlobalC::UFFT.porter );
for (j = 0;j < GlobalC::wf.npw;j++)
{
tmhpsi[j] += GlobalC::UFFT.porter[ GR_index[j] ];
}
}
else
{
std::complex<double>* porter1 = new std::complex<double>[GlobalC::pw.nrxx];
ModuleBase::GlobalFunc::ZEROS( GlobalC::UFFT.porter, GlobalC::pw.nrxx);
ModuleBase::GlobalFunc::ZEROS( porter1, GlobalC::pw.nrxx);
for (int ig=0; ig< GlobalC::wf.npw; ig++)
{
GlobalC::UFFT.porter[ GR_index[ig] ] = tmpsi_in[ig];
porter1[ GR_index[ig] ] = tmpsi_in[ig + GlobalC::wf.npwx];
}
// (2) fft to real space and doing things.
GlobalC::pw.FFT_wfc.FFT3D( GlobalC::UFFT.porter, 1);
GlobalC::pw.FFT_wfc.FFT3D( porter1, 1);
std::complex<double> sup,sdown;
for (int ir=0; ir< GlobalC::pw.nrxx; ir++)
{
sup = GlobalC::UFFT.porter[ir] * (GlobalC::pot.vr_eff(0,ir) + GlobalC::pot.vr_eff(3,ir)) +
porter1[ir] * (GlobalC::pot.vr_eff(1,ir) - std::complex<double>(0.0,1.0) * GlobalC::pot.vr_eff(2,ir));
sdown = porter1[ir] * (GlobalC::pot.vr_eff(0,ir) - GlobalC::pot.vr_eff(3,ir)) +
GlobalC::UFFT.porter[ir] * (GlobalC::pot.vr_eff(1,ir) + std::complex<double>(0.0,1.0) * GlobalC::pot.vr_eff(2,ir));
GlobalC::UFFT.porter[ir] = sup;
porter1[ir] = sdown;
}
// (3) fft back to G space.
GlobalC::pw.FFT_wfc.FFT3D( GlobalC::UFFT.porter, -1);
GlobalC::pw.FFT_wfc.FFT3D( porter1, -1);
for (j = 0;j < GlobalC::wf.npw;j++)
{
tmhpsi[j] += GlobalC::UFFT.porter[ GR_index[j] ];
}
for (j = 0;j < GlobalC::wf.npw;j++ )
{
tmhpsi[j+GlobalC::wf.npwx] += porter1[ GR_index[j] ];
}
delete[] porter1;
}
tmhpsi += dmax;
tmpsi_in += dmax;
}
}
ModuleBase::timer::tick("Hamilt_PW","vloc");
把 phi 通过 FFT 放到实空间，和 V 内积之后，再通过 fft 到把乘积结果放倒空间。
非局部赝势的计算：
//------------------------------------
// (3) the nonlocal pseudopotential.
//------------------------------------
ModuleBase::timer::tick("Hamilt_PW","vnl");
if(GlobalV::VNL_IN_H)
{
if ( GlobalC::ppcell.nkb > 0)
{
//<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
//qianrui optimize 2021-3-31
int nkb=GlobalC::ppcell.nkb;
ModuleBase::ComplexMatrix becp(GlobalV::NPOL * m, nkb, false);
char transa = 'C';
char transb = 'N';
if(m==1 && GlobalV::NPOL==1)
{
int inc = 1;
zgemv_(&transa, &GlobalC::wf.npw, &nkb, &ModuleBase::ONE, GlobalC::ppcell.vkb.c, &GlobalC::wf.npwx, psi_in, &inc, &ModuleBase::ZERO, becp.c, &inc);
}
else
{
int npm = GlobalV::NPOL * m;
zgemm_(&transa,&transb,&nkb,&npm,&GlobalC::wf.npw,&ModuleBase::ONE,GlobalC::ppcell.vkb.c,&GlobalC::wf.npwx,psi_in,&GlobalC::wf.npwx,&ModuleBase::ZERO,becp.c,&nkb);
}
Parallel_Reduce::reduce_complex_double_pool( becp.c, nkb * GlobalV::NPOL * m);
this->add_nonlocal_pp(hpsi, becp.c, m);
}
}
ModuleBase::timer::tick("Hamilt_PW","vnl");
metaGGA 部分的计算：
//------------------------------------
// (4) the metaGGA part
//------------------------------------
ModuleBase::timer::tick("Hamilt_PW","meta");
if(XC_Functional::get_func_type() == 3)
{
tmhpsi = hpsi;
tmpsi_in = psi_in;
for(int ib = 0; ib < m; ++ib)
{
for(int j=0; j<3; j++)
{
ModuleBase::GlobalFunc::ZEROS( GlobalC::UFFT.porter, GlobalC::pw.nrxx);
for (int ig = 0;ig < GlobalC::kv.ngk[GlobalV::CURRENT_K] ; ig++)
{
double fact = GlobalC::pw.get_GPlusK_cartesian_projection(GlobalV::CURRENT_K,GlobalC::wf.igk(GlobalV::CURRENT_K,ig),j) * GlobalC::ucell.tpiba;
GlobalC::UFFT.porter[ GR_index[ig] ] = tmpsi_in[ig] * complex<double>(0.0,fact);
}
GlobalC::pw.FFT_wfc.FFT3D(GlobalC::UFFT.porter, 1);
for (int ir = 0; ir < GlobalC::pw.nrxx; ir++)
{
GlobalC::UFFT.porter[ir] = GlobalC::UFFT.porter[ir] * GlobalC::pot.vofk(GlobalV::CURRENT_SPIN,ir);
}
GlobalC::pw.FFT_wfc.FFT3D(GlobalC::UFFT.porter, -1);
for (int ig = 0;ig < GlobalC::kv.ngk[GlobalV::CURRENT_K] ; ig++)
{
double fact = GlobalC::pw.get_GPlusK_cartesian_projection(GlobalV::CURRENT_K,GlobalC::wf.igk(GlobalV::CURRENT_K,ig),j) * GlobalC::ucell.tpiba;
tmhpsi[ig] = tmhpsi[ig] - complex<double>(0.0,fact) * GlobalC::UFFT.porter[ GR_index[ig] ];
}
}//x,y,z directions
}
}
ModuleBase::timer::tick("Hamilt_PW","meta");
我们接着往下看。
e[m] = ModuleBase::GlobalFunc::ddot_real(dim, phi_m, hphi );
这是计算瑞利泛函极小，得到特征值的一个初猜。
for (iter = 0;iter < maxter;iter++)
{
this->calculate_gradient( precondition, dim, hphi, sphi, g, pphi );
this->orthogonal_gradient( dim,dmx, g, scg, lagrange, phi, m);
this->calculate_gamma_cg( iter, dim, precondition, g, scg,
g0, cg, gg_last, cg_norm, theta, phi_m);// scg used as sg
converged = this->update_psi( dim, cg_norm, theta, pphi, cg, scg, phi_m ,
e[m], eps, hphi, sphi); // pphi is used as hcg
if ( converged ) break;
}//end iter
这是 CG 迭代的一个核心。我们可以把它打开。
void Diago_CG::calculate_gradient(
const double* precondition, const int dim,
const std::complex<double> *hpsi, const std::complex<double> *spsi,
std::complex<double> *g, std::complex<double> *ppsi)
这个子函数是计算梯度的子函数。首先，它计算了预条件矩阵的左乘，即 PH|psi> 和 PS|psi>。其次更新由瑞利泛函得到的特征值
λ
=
<
ψ
∣
S
P
H
∣
ψ
>
<
ψ
∣
S
P
S
∣
ψ
>
\lambda =\frac{<\psi|SPH|\psi >}{<\psi|SPS|\psi >}
λ=<ψ∣SPS∣ψ><ψ∣SPH∣ψ>
最后更新梯度：
g
=
P
H
∣
ψ
>
−
<
ψ
∣
S
P
H
∣
ψ
>
<
ψ
∣
S
P
S
∣
ψ
>
∗
P
S
∣
ψ
>
g = PH|\psi>-\frac{<\psi|SPH|\psi>}{<\psi|SPS|\psi>}*PS|\psi>
g=PH∣ψ>−<ψ∣SPS∣ψ><ψ∣SPH∣ψ>∗PS∣ψ>
void Diago_CG::orthogonal_gradient( const int &dim, const int &dmx,
std::complex<double> *g, std::complex<double> *sg, std::complex<double> *lagrange,
const ModuleBase::ComplexMatrix &eigenfunction, const int m)
这个子函数做的是梯度的正交化，首先把上一步正交化得到的梯度 g 左乘 S，得到 S|g>。其次把 S|g> 左乘特征函数矩阵，得到拉格朗日乘子。最后正交化 |g> 和 |Sg> 对于 0 到 m-1 所有的态，即形如 ∣ g m ⟩ = ∣ g ⟩ − ∑ n = 0 m − 1 ∣ Ψ n ⟩ ⟨ Ψ n ∣ S ∣ g ⟩ \left|g_{m}\right\rangle=|g\rangle-\sum_{n=0}^{m-1}\left|\Psi_{n}\right\rangle\left\langle\Psi_{n}|S| g\right\rangle ∣gm⟩=∣g⟩−∑n=0m−1∣Ψn⟩⟨Ψn∣S∣g⟩。
void Diago_CG::calculate_gamma_cg(
const int iter,
const int dim,
const double *precondition,
const std::complex<double> *g,
const std::complex<double> *sg,
std::complex<double> *psg,
std::complex<double> *cg,
double &gg_last,
const double &cg_norm,
const double &theta,
const std::complex<double> *psi_m)
这个函数体开始正式地跑 CG 算法。第一步，若非第 0 步迭代，那么更新梯度在 S 意义下的内积，gg_inter = <g|psg>。第二步，更新 psg，即 |psg> = P | Sg >。第三步，更新当前的梯度内积，gg_now = < g|P|sg > = < g|psg >。第四步，如果是第 0 次迭代，那么更新 cg 方向为梯度方向，|cg> = |g>。否则，我们执行标准的共轭梯度步，
c
g
n
+
1
=
γ
c
g
n
+
g
n
\mathrm{cg}^{n+1}=\gamma \mathrm{cg}^{n}+\mathrm{g}^{n}
cgn+1=γcgn+gn
γ
=
⟨
g
n
+
1
∣
P
∣
g
n
+
1
⟩
−
⟨
g
n
+
1
∣
P
∣
g
n
⟩
⟨
g
n
∣
P
∣
g
n
⟩
\gamma=\frac{\left\langle g^{n+1}|\mathrm{P}| g^{n+1}\right\rangle-\left\langle g^{n+1}|\mathrm{P}| g^{n}\right\rangle}{\left\langle g^{n}|\mathrm{P}| g^{n}\right\rangle}
γ=⟨gn∣P∣gn⟩⟨gn+1∣P∣gn+1⟩−⟨gn+1∣P∣gn⟩
第五步，对共轭方向做一个矫正。
cg = cg - gamma * cg_norm * sin(theta)* psi_m;
bool Diago_CG::update_psi(
const int dim,
double &cg_norm,
double &theta,
std::complex<double> *hcg,
const std::complex<double> *cg,
std::complex<double> *scg,
std::complex<double> *psi_m ,
double &eigenvalue,
const double &threshold,
std::complex<double> *hpsi,
std::complex<double> *sphi)
这一步是利用得到的共轭方向，确定步长，进行迭代更新。终止条件为 cg 的模足够小。
A
0
=
2
Re
⟨
Ψ
m
n
∣
H
^
∣
cg
n
⟩
B
0
=
⟨
cg
n
∣
H
^
∣
cg
n
⟩
\begin{gathered} A_{0}=2 \operatorname{Re}\left\langle\Psi_{m}^{n}|\hat{\mathrm{H}}| \operatorname{cg}^{n}\right\rangle \\ B_{0}=\left\langle\operatorname{cg}^{n}|\hat{H}| \operatorname{cg}^{n}\right\rangle \end{gathered}
A0=2Re⟨Ψmn∣H^∣cgn⟩B0=⟨cgn∣H^∣cgn⟩
θ
=
1
2
tan
−
1
a
0
e
0
−
b
0
\theta=\frac{1}{2} \tan ^{-1} \frac{a_{0}}{e_{0}-b_{0}}
θ=21tan−1e0−b0a0
∣
Ψ
m
n
+
1
⟩
=
cos
(
θ
)
∣
Ψ
m
n
⟩
+
sin
(
θ
)
∣
c
g
n
⟩
H
∣
Ψ
m
n
+
1
⟩
=
cos
(
θ
)
H
∣
Ψ
m
n
⟩
+
sin
(
θ
)
H
∣
c
g
n
⟩
S
∣
Ψ
m
n
+
1
⟩
=
cos
(
θ
)
S
∣
Ψ
m
n
⟩
+
sin
(
θ
)
S
∣
c
g
n
⟩
\begin{aligned} \left|\Psi_{m}^{n+1}\right\rangle &=\cos (\theta)\left|\Psi_{m}^{n}\right\rangle+\sin (\theta)\left|\mathrm{cg}^{n}\right\rangle \\ \mathrm{H}\left|\Psi_{m}^{n+1}\right\rangle &=\cos (\theta) \mathrm{H}\left|\Psi_{m}^{n}\right\rangle+\sin (\theta) \mathrm{H}\left|\mathrm{cg}^{n}\right\rangle \\ \mathrm{S}\left|\Psi_{m}^{n+1}\right\rangle &=\cos (\theta) \mathrm{S}\left|\Psi_{m}^{n}\right\rangle+\sin (\theta) \mathrm{S}\left|\mathrm{cg}^{n}\right\rangle \end{aligned}
∣∣Ψmn+1⟩H∣∣Ψmn+1⟩S∣∣Ψmn+1⟩=cos(θ)∣Ψmn⟩+sin(θ)∣cgn⟩=cos(θ)H∣Ψmn⟩+sin(θ)H∣cgn⟩=cos(θ)S∣Ψmn⟩+sin(θ)S∣cgn⟩
for (int i = 0;i < dim;i++)
{
phi(m, i) = phi_m[i];
}
更新求出来的特征向量，放回到 phi 矩阵中。
if (m > 0 && reorder)
{
//cout<<" ls out: to here"<<endl;exit(0);
ModuleBase::GlobalFunc::NOTE("reorder bands!");
if (e[m]-e[m-1]<-2.0*eps)
{
// if the last calculated eigenvalue is not the largest...
int i=0;
for (i=m-2; i>= 0; i--)
{
if (e[m]-e[i]>2.0*eps) break;
}
i++;
moved++;
// last calculated eigenvalue should be in the i-th position: reorder
double e0 = e[m];
//dcopy(phi, m, pphi);
for (int ig=0;ig<dim;ig++)
{
pphi[ig]=phi(m,ig);
}
for (int j = m;j >= i + 1;j--)
{
e[j]=e[j-1];
for (int ig=0;ig<dim;ig++)
{
phi(j,ig) = phi(j-1,ig);
}
}
e[i] = e0;
//dcopy(pphi, phi, i);
for (int ig=0;ig<dim;ig++)
{
phi(i,ig) = pphi[ig];
}
// this procedure should be good if only a few inversions occur,
// extremely inefficient if eigenvectors are often in bad order
// (but this should not happen)
} // endif
} //end reorder
这一部分处理特征值和特征向量的重排序。我们先不去管他。
扫一扫
2020
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
