abacus 基本操作

叶炜

2023-12-01

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);

};

$\left|\psi_{j}\right\rangle=H\left|\phi_{j}\right\rangle$ 等量计算的一些关键函数。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 是哈密顿量和（多个）特征向量的内积，即
$\left|\psi_{j}\right\rangle=H\left|\phi_{j}\right\rangle$
en 和 evc 是最小化能量泛函
$E\left[\phi_{i}\right]=\left\langle\phi_{i}|H| \phi_{j}\right\rangle$
得到的 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|\phi_i>$ 和 $S|\phi_i>$ 的函数。

abacus 基本操作