Note: theoretical framework of COHP may have strong basis set dependence

[kirk0830]

Details

COHP is actually a Mulliken population analysis-like way decomposing DOS:
$E^{\mathrm{band}}=\sum_{n\mathbf{k}}{w\left( \mathbf{k}\right)f_{n\mathbf{k}}\epsilon_{n\mathbf{k}}}=\sum_{\mathbf{k}}{w\left( \mathbf{k}\right)}\int_{-\infty}^{\epsilon_F}{\mathrm{d}\epsilon\sum_n{f_{n\mathbf{k}}\epsilon_{n\mathbf{k}}\delta \left( \epsilon _{n\mathbf{k}}-\epsilon \right)}}$

rewritten the eigenenergy as:

$\epsilon_{n\mathbf{k}}=\sum_{Ii,Jj}{c_{Ii,n}^{*}\left(\mathbf{k} \right)H_{Ii,Ji}\left(\mathbf{k}\right) c_{Jj,n}\left(\mathbf{k}\right)}$

Band energy is then in the following form:

$E^{\mathrm{band}}=\sum_{n\mathbf{k}}{w\left( \mathbf{k} \right) f_{n\mathbf{k}}\sum_{Ii,Jj}{c_{Ii,n}^{*}\left( \mathbf{k} \right) H_{Ii,Ji}\left( \mathbf{k} \right) c_{Jj,n}\left( \mathbf{k} \right)}}$

Rearrange, then COHP term will emerge:

$\sum_{IJ}{\left( \sum_{n\mathbf{k}}{w\left( \mathbf{k} \right) f_{n\mathbf{k}}\sum_{ij}{c_{Ii,n}^{*}\left( \mathbf{k} \right) H_{Ii,Ji}\left( \mathbf{k} \right) c_{Jj,n}\left( \mathbf{k} \right)}} \right)}$

Define COHP between atom I and atom J is:

$COHP_{IJ}\left( \epsilon \right) =\sum_{ij}{COHP_{ij}\left( \epsilon \right)}$

indices i and j are those of orbitals belong to I and J, respectively,

$COHP_{ij}\left( \epsilon \right) =\sum_{n\mathbf{k}}{w\left( \mathbf{k} \right) f_{n\mathbf{k}}\delta \left( \epsilon_{n\mathbf{k}}-\epsilon \right) c_{Ii,n}^{*}\left( \mathbf{k} \right) H_{Ii,Ji}\left( \mathbf{k} \right) c_{Jj,n}\left( \mathbf{k} \right)}$

.
I have quickly implement one version of COHP directly based on ABACUS lcao:

To calculate $\Re \left[ c_{Ii,n}^{*}\left( \mathbf{k} \right) H_{Ii,Ji}\left( \mathbf{k} \right) c_{Jj,n}\left( \mathbf{k} \right) \right] $

def cal_COHPmatskIJ_e_ij(Hk, Sk, Ck, atomI_orbs, atomJ_orbs, mode="COHP"):
    """Calculate the energy-resolved COHP or COOP matrix for all bands, every selected orbital pair of the atom I and J, and a k-point.
    
    Args: 
    Hk: Hamiltonian matrix for a k-point, H(k),
    Sk: Overlap matrix for a k-point, S(k),
    Ck: Wavefunction set of one k-point, arranged as C(k)[i, n], where i is the local index and n is the band index.
    atomI_orbs: indices of the orbitals of the atom I,
    atomJ_orbs: indices of the orbitals of the atom J.
    mode: "COHP" or "COOP".

    Returns:
    np.ndarray: energy-resolved COHP or COOP matrix for a k-point.
    """
    nrows, ncols = Hk.shape
    nlocal, nbands = Ck.shape
    assert nrows == ncols
    assert nrows == nlocal
    # allocate memory for value
    value = [np.zeros((len(atomI_orbs), len(atomJ_orbs)), dtype=np.float64) for i in range(nbands)]

    for ib in range(nbands):
        for i_inval, i_inorb in enumerate(atomI_orbs):
            for j_inval, j_inorb in enumerate(atomJ_orbs):
                power = Sk[i_inorb, j_inorb] if mode == "COOP" else Hk[i_inorb, j_inorb]
                value[ib][i_inval, j_inval] = (Ck[i_inorb, ib].conj()*power*Ck[j_inorb, ib]).real
    return value

Then calculate $\sum_{ij}{\Re \left[ c_{Ii,n}^{*}\left( \mathbf{k} \right) H_{Ii,Ji}\left( \mathbf{k} \right) c_{Jj,n}\left( \mathbf{k} \right) \right]}$:

def cal_COHPvalskIJ_e(Hk, Sk, Ek, Ck, atomI_orbs, atomJ_orbs, mode="COHP"):
    """Compute the sum_{ij} over c*{Ii,n}(k)H{IiJj}(k)c{Jj,n}(k), 
    where I, J are atom indexes and i, j are orbitals indexes.
    n is the band index and this quantity corresponds to an eigenenergy.

    Args:
        Hk (np.ndarray): Hamiltonian matrix for a k-point.
        Ck (np.ndarray): Wavefunction set of one k-point.
        atomI_orbs (list): indices of the orbitals of the atom I.
        atomJ_orbs (list): indices of the orbitals of the atom J.

    Returns:
        tuple: eigenenergies and the sum_{ij} on c*{Ii,n}(k)H{IiJj}(k)c{Jj,n}(k) for each band.
    """
    
    nlocal, nband = Ck.shape
    matskIJ_ij_e = cal_COHPmatskIJ_e_ij(Hk, Sk, Ck, atomI_orbs, atomJ_orbs, mode=mode)
    # dimension assertation
    assert len(Ek) == nband
    assert len(matskIJ_ij_e) == nband
    for ib in range(nband):
        assert matskIJ_ij_e[ib].shape == (len(atomI_orbs), len(atomJ_orbs))
    # compute the sum_{ij} over c*{Ii,n}(k)H{IiJj}(k)c{Jj,n}(k), i.e., for atom-pair IJ
    valskIJ_e = [np.sum(matskIJ_ij_e[i]) for i in range(nband)]
    return Ek, valskIJ_e

Then $w\left( \mathbf{k} \right) \sum_{ij}{\Re \left[ c_{Ii,n}^{*}\left( \mathbf{k} \right) H_{Ii,Ji}\left( \mathbf{k} \right) c_{Jj,n}\left( \mathbf{k} \right) \right]}$, which means summation over kpoints should consider their weights due to symmetry.

def cal_COHPvalsIJ_e(Hks, Sks, Eks, Cks, wk = None, 
                     atomI_orbs = None, atomJ_orbs = None,
                     mode: str = "COHP"):
    nks = len(Hks)
    nbands = len(Eks[0])

    wk = [1/nks for i in range(nks)] if wk is None else wk

    Evals = []
    COHPvalsIJ_e = []
    for ik in range(nks):
        Evalsk, COHPvalskIJ_e = cal_COHPvalskIJ_e(Hk=Hks[ik], Sk=Sks[ik], Ek=Eks[ik], Ck=Cks[ik], 
                                                  atomI_orbs=atomI_orbs, atomJ_orbs=atomJ_orbs, mode=mode)
        # dimension assertation
        assert len(Evalsk) == nbands
        assert len(COHPvalskIJ_e) == nbands
        # add k-point weight
        COHPvalskIJ_e = [wk[ik]*COHPvalskIJ_e[i] for i in range(nbands)]
        # add to the global list
        Evals += Evalsk.tolist()
        COHPvalsIJ_e += COHPvalskIJ_e
    # dimension assertation
    assert len(Evals) == nks*nbands
    assert len(COHPvalsIJ_e) == nks*nbands
    # sort COHPvalsIJ_e according to the Evals
    Evals, COHPvalsIJ_e = zip(*sorted(zip(Evals, COHPvalsIJ_e)))
    # energy unit conversion
    Evals = np.array([rao.unit_conversion(e, "Ry", "eV") for e in Evals])
    return dos_integral(Evals, COHPvalsIJ_e)

Then after zero padding, the COHP between atomI and atomJ is obtained.

However, this result is totally different with what LOBSTER calculated. LOBSTER uses STO or GTO->STO, and the COHP looks similar with Mulliken population analysis, therefore it is not basis set independent. I doubt about the whether COHP is directly useful for ABACUS LCAO with numerical orbitals.

Task list for Issue attackers (only for developers)

• Reproduce the performance issue on a similar system or environment.
• Identify the specific section of the code causing the performance issue.
• Investigate the issue and determine the root cause.
• Research best practices and potential solutions for the identified performance issue.
• Implement the chosen solution to address the performance issue.
• Test the implemented solution to ensure it improves performance without introducing new issues.
• Optimize the solution if necessary, considering trade-offs between performance and other factors (e.g., code complexity, readability, maintainability).
• Review and incorporate any relevant feedback from users or developers.
• Merge the improved solution into the main codebase and notify the issue reporter.

[WHUweiqingzhou]

Could you explain a little about your next plan?

[kirk0830]

Could you explain a little about your next plan?

@WHUweiqingzhou I am afraid that simply declaring the difference between COHP result of ABACUS numerical atomic orbitals and LMTO or sth. is because of basis-dependence of Mulliken would not be convincing to all users. Therefore a LCAO_IN_PW run produces wavefunction in pw representation, which is also compatible with LOBSTER, will be a promising choice.

I will try to connect a hdf5 lib to print out pw wavefunction.