At the
end of 2009-06, I tried to understand Davidson method for the first
time. But I almost got nothing,except Paolo
Giannozzi's paper
at
http://www.democritos.it/pipermail/pw_forum/2007-November/007584.html.
At than time I wrote a blog for it at 2009-06-24. Now the following
content covers that one, therefore I delete the old
blog:)
This is my
second try. I spend almost two
days to get some progress in understanding of Davidson
diagonalization. I am not interested in Davidson method itself. But
without knowing it, you can not understand the real meaning of
parameter
"nband".
Part
one:
Suppose the digenpair (λ, u) is exact solution of Au = λu, and
(λk~, uk~) is the Ritz approximate eigenpair of Au = λu.
Let
replace λ by the
currect approximate λk~ ,
M is
called suitable preconditioners.
From part one, you can see the general idea of Davidson.
Apparently, it can be viewed as a
perturbation method.
Part
two:
Be aware that the vector v1 in K
subspace is n*1, not m*1. And
y is m*1. Usually, n
>>
m In QE, n is
determined by parameter "ecutwfc". That indicates, although
n is very large, m can be constrained in a small
number, which represents parameter "nband".
The basic idea of Davison method is easy to understand. Suppose we
have a subspace K of dimension k, for which we have an
orthonormal basis {v1, v2, ..., vk}. We definite
Vk as the n by k matrix with column
vi, then the matrix of A reduced to a small matrix of
Hk, where Hk = conj(Vk)*A*Vk. The
dominant eigenpair (λ~,y) of Hk is compute, and the Ritz
approximate eigenvector u=Vk*y.
Now, you can
undertand what it indicates in the <QUANTUM
ESPRESSO: a modular and open-source software project for quantum
simulations of materials>
and
in Paolo Giannozzi's paper at
http://www.democritos.it/pipermail/pw_forum/2007-November/007584.html:
以下为个人体会，如果错误，概不负责。更期待高人指正；）
在计算中指定了截止动能，对于一个固定的k点，也就意味着指定了有多少个G（=2pi/a的整数倍）参与傅立叶变换（假设为M个）。而G的个数就是相应哈密顿量的维度，M*M维。通常M是非常大的，而我们只需要知道较低的几个被占据的能级，即可求出总能等一些列刚兴趣的物理量。因此直接对M*M维的哈密顿量进行对角化是不合适的。所以需要采用迭代的方法来求解几个低能的本征值和本征能量即可。
其实迭代方法在矩阵对角化中是非常常用的数值计算方法，他可以通过迭代的方式只求出某一个或者几个本征值。Davidson是其中一个非常出名的，因为他在量子力学中取得了很大的成功，这个成功来源于量子力学的哈密顿矩阵通常都是稀疏的，就是非对角元很多都是0，或者接近0。
当nband被指定后，本征方程的就被给定。。其中N应该是等于nband的。此时N远小于有ecutwfc决定的G的个数，即H的维度M*M。而波函数psi则是由M个平面波的叠加。此时的情况正好跟Part
two中的第一张图对应。在这里，可以对此方程的左乘波函数psi(j)，将得到一个N*N维的哈密顿量，然后将它对角化，求得初始猜测波函数对应的本征值。
然后按照求得对初始猜测波函数的修正（方程左边的波函数跟最左边的那个是正交的）。这样的修正波函数作为独立的基矢，加入到初始猜测的N个波函数上，构成N+1维的子空间。如果对N个初始猜测波函数均采用这种办法，则得到2*N维的子空间。然后，在2*N维的子空间中求得H矩阵，并对其对角化，比较前N个本征值是否有变化，判定是否收敛。
补记：
一：
Need to diagonalize a matrix of size
N_PW*N_PW
• NPW >> Nb = number of bands
required = Ne/2 or a little more (for metals).
• OK to obtain lowest few eigenvalues.
• Exact diagonalization is expensive!
• Use iterative diagonalizers that recast
diagonalization as a minimization problem.
(2009-10-07)
二：
Today I bought a
book<量子化学中的计算方法>。At 140 page, there
give some explanation about Davidson method.
将“然后按照......判定是否收敛"修正如下：
然后对感兴趣的第k个波函数按照(将其中的i全部改为k)求得对初始猜测波函数的修正（方程左边的波函数跟最右边的那个是正交的）。这样的修正波函数作为独立的基矢，加入到初始猜测的N个波函数上，构成N+1维的子空间。并对其对角化，并判定第k个本征值是否收敛。(这段的说明并不一定比以前的那段更正确)
The last sentence of
<量子化学中的计算方法> about Davidson is
“Pulay等指出，为了节省计算机内存，如果子空间大小超过某一设定的值，比如50时，可以只保留最新的两个特征矢量x1和x2，将子空间的大小将为2，重新开始迭代。数值计算表明，迭代的收敛性质基本保持不变”(2009-10-27)
I've
never tried computing quite that many bands (only gone up to about
1200 bands) or with that many atoms per cell. However, I've never
not been able to get the bands to converge using both 'cg'
diagonalization and also with increasing the cutoff. You mentioned
that you have tried increasing the cutoff, but perhaps you haven't
gone high enough?
Here
is a quote from Paolo that might be relevant in your
situation:
the algorithm used in PWscf is based on the assumption that the number of Kohn-Sham states (aka bands) is much smaller than the number of basis functions (i.e. plane waves). If this
assumption doesn't hold (and it doesn't if you calculate 500 bands in fcc Si unless you use a very large cutoff) you may run into trouble. Moreover iterative diagonalization becomes
slower and more memory-consuming than conventional diagonalization. If you really need to do that, you need to modify the code(2009-12-02)
This
explanation are copied from VASP
instruction:
One should chose NBANDS so that a considerable number of empty
bands is included in the calculation. As a minimum
we require one empty band. VASP will give a
warning, if this is not the case. NBANDS is also
important from a technical point of view: In iterative
matrix-diagonalization schemes eigenvectors
close to the top of the calculated number of
vectors converge much slower than the lowest eigenvectors. This
might result in a significant performance loss if
not enough empty bands are included in the calculation. Therefore
we recommend to set NBANDS to NELECT/2 + NIONS/2,
this is also the default setting of the makeparam utility and of
VASP.4.X. This setting is safe in most cases. In
some cases, it is also possible to decrease the number of
additional bands to NIONS/4 for large
systems without performance loss, but on the other
hand transition metals do require a much larger number of empty
bands (up to 2*NIONS).
To check this parameter perform several calculations for a
fixed potential (ICHARG=12) with an increasing number
of bands (e.g. starting from NELECT/2 + NIONS/2).
An accuracy of 10−6 should be obtained in 10-15 iterations. Mind
that the RMM-DIIS scheme (IALGO=48) is more
sensible to the number of bands than the default CG algorithm
(IALGO=8).(2009-03-30)
加载中，请稍候......
加载中…
(2009-08-08 15:32:24)