abacus 基本操作

文章目录

abacus 基本操作
- 下载和运行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, PseudopotentialLATTICE_CONSTANT
10.2            #Lattice constantLATTICE_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 3ATOMIC_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 format0        The Pseudo was generated with a Non-Relativistic Calculation0.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 NumberSi                   ElementNC                  Norm - Conserving pseudopotentialF                  Nonlinear Core CorrectionSLA  PZ   NOGX NOGC   PZ   Exchange-Correlation functional4.00000000000      Z valence0.00000000000      Total energy0.0000000  0.0000000 Suggested cutoff for wfc and rho1                  Max angular momentum component431                  Number of points in mesh2    2             Number of Wavefunctions, Number of ProjectorsWavefunctions         nl  l   occ3S  0  2.003P  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 __CUDAdouble *d_g2kin;
#endifpublic:// Calculate kinetic energyvoid 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-23std::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 __CUDAcudaMemcpy(this->d_g2kin, this->g2kin, GlobalC::kv.ngk[ik] * sizeof(double), cudaMemcpyHostToDevice);
#endifModuleBase::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,//0const int &iter,//1const int &ik,//0const double *precondition,//precondition[0] = 2.41424double &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-31void 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,//0const int nstart,//8const int n_band,//4const 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>的函数。