abacus 基本操作
下载和运行abacus
设置代理
export HTTP_PROXY=http://sys-proxy-rd-relay.byted.org:8118
export HTTPS_PROXY=http://sys-proxy-rd-relay.byted.org:8118
克隆
git clone https://github.com/deepmodeling/abacus-develop.git
运行方式，在参数文件夹下执行
abacus所在路径/abacus
输入输出
INPUT 文件
INPUT_PARAMETERS
#Parameters (General)
pseudo_dir ./
ntype 1
nbands 4
#Parameters (Accuracy)
ecutwfc 50
scf_nmax 20
symmetry 1
ntype：在原包内有多少种类型的原子。
nbands：要被计算的能带数量。
ecutwfc：平面波的截断能。
scf_nmax：迭代过程的最大迭代次数。默认 40。
symmetry：取值为 0 或者 1。如果是 1，将用对称分析来决定布拉菲点阵的类型。
STRU
#This is the atom file containing all the information
#about the lattice structure.
ATOMIC_SPECIES
Si 1.000 Si.pz-vbc.UPF #Element, Mass, Pseudopotential
LATTICE_CONSTANT
10.2 #Lattice constant
LATTICE_VECTORS
0.5 0.5 0.0 #Lattice vector 1
0.5 0.0 0.5 #Lattice vector 2
0.0 0.5 0.5 #Lattice vector 3
ATOMIC_POSITIONS
Cartesian #Cartesian(Unit is LATTICE_CONSTANT)
Si #Name of element
0.0 #Magnetic for this element.
2 #Number of atoms
0.00 0.00 0.00 0 0 0 #x,y,z, move_x, move_y, move_z
0.25 0.25 0.25 1 1 1
第 5 行，指定标签、质量、赝势文件名。
第 8 行，晶格缩放因子。
第 11~13 行，晶格矢量。
第 16 行，方向，笛卡尔或者直接坐标。
KPT 文件
K_POINTS
0
Gamma
4 4 4 0 0 0
第 2 行表示 k 点总数，0 表示自动。
第 3 行选择 Monkhorst-Pack 方法， Gamma 或者 MP。
第 4 行前三个数，沿着倒格矢的分割。后三个数网格的移动。
UPF 文件
<PP_INFO>
Generated using unknown code
Author: Von Barth-Car (<1984)
Info: automatically converted from PWSCF format
0 The Pseudo was generated with a Non-Relativistic Calculation
0.00000000000E+00 Local Potential cutoff radius
nl pn l occ Rcut Rcut US E pseu
3S 0 0 2.00 0.00000000000 0.00000000000 0.00000000000
3P 0 1 2.00 0.00000000000 0.00000000000 0.00000000000
</PP_INFO>
<PP_HEADER>
0 Version Number
Si Element
NC Norm - Conserving pseudopotential
F Nonlinear Core Correction
SLA PZ NOGX NOGC PZ Exchange-Correlation functional
4.00000000000 Z valence
0.00000000000 Total energy
0.0000000 0.0000000 Suggested cutoff for wfc and rho
1 Max angular momentum component
431 Number of points in mesh
2 2 Number of Wavefunctions, Number of Projectors
Wavefunctions nl l occ
3S 0 2.00
3P 1 2.00
</PP_HEADER>
……
……
……
定义赝势相关的参数。
源码解读
class WF_igk
{
public:
WF_igk();
~WF_igk();
//---------------------------------------------------
// npwx: max npw
// npw
// igk: [nks, npw_max]
// g2kin: [npwx],kinetic energy for current k point
//---------------------------------------------------
int npwx;
int npw;
ModuleBase::IntArray igk;
double *g2kin;
#ifdef __CUDA
double *d_g2kin;
#endif
public:
// Calculate kinetic energy
void ekin(const int ik);
double* get_qvec_cartesian(const int &ik);
ModuleBase::Vector3<double> get_1qvec_cartesian(const int ik,const int ig)const;
std::complex<double>* get_sk(const int ik, const int it, const int ia)const;
// pengfei 2016-11-23
std::complex<double>* get_skq(int ik, int it, int ia, ModuleBase::Vector3<double> q);
};
这个模块包含了动能相关，g2kin 是动能向量。其中 ekin 是动能计算函数。
//--------------------------------------------------------
// Compute kinetic energy for each k-point
//--------------------------------------------------------
void WF_igk::ekin(const int ik)
{
ModuleBase::timer::tick("WF_igk", "ekin");
ModuleBase::GlobalFunc::ZEROS(g2kin, this->npwx);
for (int ig = 0; ig < GlobalC::kv.ngk[ik]; ig++)
{
//--------------------------------------------------------
// EXPLAIN : Put the correct units on the kinetic energy
//--------------------------------------------------------
this->g2kin[ig] = GlobalC::pw.get_GPlusK_cartesian(ik, this->igk(ik, ig)).norm2() * GlobalC::ucell.tpiba2;
}
#ifdef __CUDA
cudaMemcpy(this->d_g2kin, this->g2kin, GlobalC::kv.ngk[ik] * sizeof(double), cudaMemcpyHostToDevice);
#endif
ModuleBase::timer::tick("WF_igk", "ekin");
return;
}
动能向量的具体计算函数。
void Electrons::c_bands(const int &istep)
这是电子能带计算的一个核心函数。包含了动能项的计算，计算哈密顿量和特征向量的乘积，通过最小化能量泛函得到特征值和特征向量。h_diag 为预条件矩阵。
else if (precondition_type==2)
{
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];
}
}
这是 c_bands 中构造某种预条件子部分，其中利用了构造出来的 g2_kin 动能向量。
void Hamilt::diagH_pw(
const int &istep,//0
const int &iter,//1
const int &ik,//0
const double *precondition,//precondition[0] = 2.41424
double &avg_iter)//0
对角分解是能带计算的核心步骤，iter 表示迭代次数，ik 表示计算第几个 k 点。precondition 表示预条件数组。
if(GlobalV::KS_SOLVER=="cg")
{
if ( iter > 1 || istep > 1 || ntry > 0)
{
this->diagH_subspace(
ik,
GlobalV::NBANDS,
GlobalV::NBANDS,
GlobalC::wf.evc[ik0],
GlobalC::wf.evc[ik0],
GlobalC::wf.ekb[ik]);
avg_iter += 1.0;
}
Diago_CG cg(&GlobalC::hm.hpw);
bool reorder = true;
if(GlobalV::NPOL==1)
{
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 ks 方程解法器执行的代码。可以看得出来，当 iter 小于等于 1 的时候，开始的那段代码不能执行，等于 2 的时候就执行。这是因为，在平面波方法中，对角化哈密顿量的时候，首先用 diagH_subspace 去对角化，然后再用 cg 方法。
wf.ekb[ik] 存的是第 ik 个 k 点对应的 n_band 个特征值。wf.evc[ik0] 原来存的是哈密顿量和特征向量的内积，后来存的是第 ik0 个 k 点计算得到的 n_band 个特征向量矩阵。
class Hamilt_PW
{
public:
void diagH_subspace(const int ik,
const int nstart,
const int nbnd,
const ModuleBase::ComplexMatrix &psi,
ModuleBase::ComplexMatrix &evc,
double *en);
void h_1psi(
const int npw,
const std::complex<double> *psi1d,
std::complex<double> *hpsi,
std::complex<double> *spsi);
void h_psi(
const std::complex<double> *psi,
std::complex<double> *hpsi,
const int m = 1); // qianrui add a default parameter 2021-3-31
void s_1psi(
const int npw,
const std::complex < double> *psi,
std::complex < double> *spsi);
};
∣ ψ j ⟩ = H ∣ ϕ j ⟩ \left|\psi_{j}\right\rangle=H\left|\phi_{j}\right\rangle ∣ψj⟩=H∣ϕj⟩ 等量计算的一些关键函数。h_1psi 调用的核心本质是 h_psi，动能项什么的其实也是在 h_psi 中进行了计算。利用 h_1psi 计算得到 hpsi 和 spsi。s_1psi 计算 sphi = S|psi(m)>。
//----------------------------------------------------------------------
// Hamiltonian diagonalization in the subspace spanned
// by nstart states psi (atomic or random wavefunctions).
// Produces on output n_band eigenvectors (n_band <= nstart) in evc.
//----------------------------------------------------------------------
void Hamilt_PW::diagH_subspace(
const int ik,//0
const int nstart,//8
const int n_band,//4
const ModuleBase::ComplexMatrix &psi,
ModuleBase::ComplexMatrix &evc,
double *en)
这个模块的作用是通过对角化得到特征分解。ik 为 k 点指标。n_band 为能带个数，也就是求得的前几个特征值个数。psi 是哈密顿量和（多个）特征向量的内积，即
∣
ψ
j
⟩
=
H
∣
ϕ
j
⟩
\left|\psi_{j}\right\rangle=H\left|\phi_{j}\right\rangle
∣ψj⟩=H∣ϕj⟩
en 和 evc 是最小化能量泛函
E
[
ϕ
i
]
=
⟨
ϕ
i
∣
H
∣
ϕ
j
⟩
E\left[\phi_{i}\right]=\left\langle\phi_{i}|H| \phi_{j}\right\rangle
E[ϕi]=⟨ϕi∣H∣ϕj⟩
得到的 n_band 个特征值和特征向量。
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
)
//-------------------------------------------------------------------
// "poor man" iterative diagonalization of a std::complex hermitian matrix
// through preconditioned conjugate gradient algorithm
// Band-by-band algorithm with minimal use of memory
// Calls h_1phi and s_1phi to calculate H|phi> and S|phi>
// Works for generalized eigenvalue problem (US pseudopotentials) as well
//-------------------------------------------------------------------
使用 CG 算法（ccgdiagg）进行对角化。其中 h_1phi 和 s_1phi 是计算 H ∣ ϕ i > H|\phi_i> H∣ϕi> 和 S ∣ ϕ i > S|\phi_i> S∣ϕi>的函数。
扫一扫
502
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
