Large-scale ab initio simulations based on systematically improvable atomic basis

Pengfei Li; Xiaohui Liu; Mohan Chen; Peize Lin; Xinguo Ren; Lin Lin; Chao Yang; Lixin He

Source-linked AI summary

Large-scale ab initio simulations based on systematically improvable atomic basis

Pengfei Li, Xiaohui Liu, Mohan Chen, Peize Lin, Xinguo Ren, Lin Lin, Chao Yang, Lixin He

arXiv:1503.00097v1Feb 28, 2015cond-mat.mtrl-sci

TL;DR

The paper presents ABACUS, a density-functional-theory package using numerical atomic basis sets and evaluates its hierarchical CGH basis sets across finite and extended systems. The results support systematic convergence toward plane-wave accuracy, while DZP provides a practical balance between accuracy and computational load.

Problem
The paper examines whether hierarchical numerical atomic basis sets can provide reliable electronic-structure simulations for both finite and extended systems.
Method
ABACUS combines density functional theory with numerical atomic basis sets, including hierarchically generated CGH orbitals and systematically improvable basis sets.
Results
The hierarchical CGH basis sets converge systematically toward plane-wave accuracy, and the DZP basis set balances accuracy with computational load.
Takeaways & Limitations
DZP can be used in production simulations, while CGH orbitals are reliable for finite and extended systems.

Abstract

from arXiv · show

We present a first-principles computer code package (ABACUS) that is based on density functional theory and numerical atomic basis sets. Theoretical foundations and numerical techniques used in the code are described, with focus on the accuracy and transferability of the hierarchical atomic basis sets as generated using a scheme proposed by Chen, Guo and He [J. Phys.:Condens. Matter \textbf{22}, 445501 (2010)]. Benchmark results are presented for a variety of systems include molecules, solids, surfaces, and defects. All results show that the ABACUS package with its associated atomic basis sets is an efficient and reliable tool for simulating both small and large-scale materials.

4 Department of Mathematics, University of California, Berkeley and Computational Research Division,

ABACUS is a DFT package built around numerical atomic orbitals and systematically improvable CGH basis sets, with an alternative plane-wave basis for benchmarking. Benchmarks across molecules, solids, surfaces, and defects support its accuracy, transferability, and production use.

Numerical atomic orbitals provide compact, localized bases that can support linear-scaling or other algorithms with better-than-O(N^3) behavior.
ABACUS is a first-principles DFT package developed around numerical atomic orbitals and the CGH procedure for generating systematically improvable basis sets.
At the LDA and GGA levels, ABACUS supports typical electronic-structure calculations, structure relaxations, and molecular dynamics.
The CGH orbitals are evaluated for molecules, solids, surfaces, and defects across a broader range of elements, including alkali, transition-metal, group VI, and group VII elements.
The benchmarks characterize ABACUS with CGH orbitals as reliable for finite and extended systems, while DZP offers a compromise between accuracy and computational cost suitable for production calculations.
The package offers both numerical atomic and plane-wave basis sets, enabling direct accuracy and consistency checks in benchmark calculations.

B. Systematically improvable atomic basis sets

ABACUS uses the CGH scheme to generate hierarchical numerical atomic orbitals that can be systematically improved and optimized for accuracy and transferability. The basis construction combines flexible spherical-Bessel expansions, spillage minimization, and kinetic-energy regularization against selected reference systems.

ABACUS adopts the CGH scheme to generate systematically improvable, optimized atomic basis sets.
Basis construction: Radial functions are expanded in spherical Bessel functions with a finite cutoff, whose number is controlled by a kinetic-energy cutoff.
Optimization: The coefficients are optimized by minimizing spillage between the atomic basis and wave functions from selected plane-wave reference systems.
Optimization: Simulated annealing determines the coefficients, while kinetic-energy minimization suppresses unphysical oscillations that can reduce basis-set transferability.
Transferability: Users can choose target systems, angular momenta, and radial multiplicities, producing a hierarchy with systematic convergence behavior toward the plane-wave reference.
Performance: The generated atomic orbitals show excellent accuracy and transferability across various systems.

C. Hamiltonian and overlap matrices construction

ABACUS constructs Hamiltonian and overlap matrices using efficient two-center and grid-integral techniques. Locality makes the matrices sparse, while grid-integral costs scale linearly with system size and can be parallelized.

Locality: Only matrix elements between overlapping atomic orbitals are evaluated, producing sparse matrices and O(N) scaling in the number of integrals.
Matrix construction: Hamiltonian components are computed using two-center integrals for kinetic, nonlocal pseudopotential, and overlap terms, alongside grid integration for local potentials.
Two-center integrals: Two-center integrals separate into radial and angular parts, with radial values tabulated and interpolated efficiently over orbital distances.
Grid integrals: Local potentials are evaluated on a uniform real-space grid, combining local pseudopotential, Hartree, and exchange-correlation contributions.
Computational cost: Grid integrals are among the most time-consuming parts of atomic-orbital algorithms, but their computational effort scales linearly and is readily parallelized.

D. Kohn-Sham equation solvers

ABACUS uses PEXSI as an alternative to diagonalization for solving Kohn-Sham problems in large systems. PEXSI exploits sparse Hamiltonian and overlap matrices, with cost scaling from O(N) in quasi-1D systems to O(N2) in three-dimensional bulk systems.

PEXSI method: PEXSI expands the density matrix through a pole expansion and computes only selected Green’s-function elements needed for the electron density.
PEXSI method: The pole expansion requires a number of terms proportional to log(β∆E), making Fermi-operator expansion efficient.
Scaling: PEXSI costs O(N) for quasi-1D systems, O(N1.5) for quasi-two-dimensional systems, and O(N2) for three-dimensional bulk systems.
Scaling: The favorable scaling depends on sparse Hamiltonian and overlap matrices rather than an assumption about density-matrix localization.
Capabilities: PEXSI obtains electron density, free energy, forces, and density-of-states quantities without computing eigenvalues or eigenvectors.

E. Total energy and force calculations

ABACUS evaluates total energies by combining Kohn–Sham and ion–ion contributions, while its force formalism includes Feynman–Hellmann, Pulay, and non-orthogonal terms. The package supports structural relaxation and systematically improvable localized atomic bases benchmarked against plane-wave results.

Forces include Feynman–Hellmann contributions, Pulay forces from changing atomic orbitals, and non-orthogonal-basis terms.
Etot is decomposed into the Kohn–Sham electronic energy and the Coulomb interaction energy between ions.
The force implementation combines reciprocal-space treatments, two-center integrals, and real-space grid integrals for different terms.
BFGS and conjugate-gradient algorithms are implemented for structural relaxation by searching local minima on the potential-energy surface.
CGH orbitals use hierarchical SZ, DZ, DZP, TZDP, and QZTP bases whose convergence is tested against plane-wave references.

A. Eggbox effect

The eggbox effect is an artificial energy ripple caused by finite real-space-grid integration errors, complicating force calculations. In ABACUS, reciprocal-space confinement of basis functions keeps the effect small enough for practical calculations and structural relaxation without correction.

The eggbox effect is artificial ground-state energy rippling caused by evaluating Hamiltonian integrals on a finite uniform real-space grid.
Suppressing high-energy Fourier components of localized orbitals allows a less dense real-space grid while preserving locality.
ABACUS automatically confines basis functions in reciprocal space below a chosen energy cutoff during construction.
For Si2 and O2, total-energy oscillations remain within 1 meV and force oscillations within 1 meV/Å.
Mn2 has the largest force oscillation, but it remains within 10 meV/Å, and structural relaxations require no further eggbox correction.

1. Bond lengths

Hierarchical atomic basis sets reproduce molecular bond lengths increasingly closely as their size grows, with polarization functions improving accuracy. At QZTP, mean absolute errors are 0.004 Å for LDA and 0.003 Å for PBE relative to plane-wave results.

Equilibrium bond lengths systematically approach plane-wave results as the atomic-orbital basis increases from SZ to QZTP.
0.018 Å is the DZP mean absolute error for both LDA and PBE bond-length calculations.
0.004 Å for LDA and 0.003 Å for PBE are the QZTP bond-length mean absolute errors.
Alkali-metal bonding is described satisfactorily for Na2 and LiH when large 10–12 Bohr cutoff radii are used.
Converged LDA bond lengths are systematically smaller than experiment, whereas converged PBE bond lengths are systematically larger.

4. Weak interaction energy

For weakly interacting S22 dimers, increasing the atomic-orbital basis systematically approaches plane-wave interaction energies and generally moves results closer to CCSD(T) references. The tests support ABACUS for van der Waals interactions and extend transferability checks from molecules to solids.

The S22 benchmark contains 22 dimers spanning hydrogen bonding, dispersion, and mixed interaction types.
Interaction energy is defined as the dimer energy minus the energies of its two fully relaxed monomers.
Increasing the number of atomic orbitals systematically approaches plane-wave results and averages closer to CCSD(T) references.
The ABACUS atomic-basis implementation with PBE-D2 is reported as suitable for describing van der Waals forces.
For crystalline solids, lattice-constant errors reach 0.02 Å with DZP and 0.01 Å with TZ(DP), confirming accurate structural descriptions.

2. Cohesive energies

The hierarchical LCAO basis sets systematically approach plane-wave accuracy for solid-state properties, with cohesive energies converging especially well. DZP and larger bases provide accurate bulk-modulus results, while DZP also reproduces Si(100) surface energetics and geometries against plane-wave references.

Cohesive energies: Cohesive-energy MAE is 0.04 eV at DZP for solids, compared with 0.09 eV for molecular atomization energies.At TZDP and QZTP levels, the solid cohesive-energy MAE is 0.01 eV, twice smaller than the molecular value.
Bulk moduli: Bulk-modulus MAE decreases from 21.2 GPa with SZ to 5.8 GPa with DZ and 2.1 GPa with DZP.Bases larger than DZP reduce the MAE to around 1.0 GPa.
Bulk properties: The basis-set quality is reported as equally good for transition-metal and main-group elements.This supports the transferability of the tested hierarchical basis sets across the solid systems considered.
Si(100) surface reconstruction: ABACUS and QE produce the same energy ordering for the three Si(100) reconstructions relative to the ideal surface.ABACUS with PW basis gives almost identical results to QE, while DZP yields energy lowerings of 1.481, 0.078, and 0.070 eV per Si dimer.
Si(100) surface reconstruction: For Si(100), ABACUS DZP reproduces plane-wave surface properties accurately, with bond lengths agreeing within 0.001 Å and bond angles within 0.4 degrees at most.The paper concludes that DZP demonstrates reliable surface calculations for all three tested reconstructions.

E. N defect in bulk GaAs

Large LCAO supercells enable first-principles study of dilute GaAs1−xNx, revealing continuous band-gap reduction with increasing nitrogen concentration and band-gap closure above x = 0.125. The approach reaches concentrations as low as x = 0.002, beyond earlier computational studies.

Band-gap behavior: GaAs1−xNx exhibits an unusually large bowing parameter of about 16 eV, producing pronounced band-gap narrowing even for x < 0.015.The paper contrasts this with most isovalent semiconductor alloys, whose bowing parameters are only a fraction of an eV.
Band-gap behavior: The calculated GaAs1−xNx band gap continuously decreases as x increases, reproducing the reported experimental trend.The LDA band gap of bulk GaAs is 1.13 eV in these calculations.
Band-gap behavior: Band-gap closure occurs when the nitrogen concentration exceeds x = 0.125.The exact curve may vary with supercell shape, but the general trend is expected to remain unchanged.

2. Formation energy of GaAs1−xNx

The formation-energy analysis evaluates nitrogen defects in GaAs1−xNx under As-rich and As-poor conditions using ABACUS LCAO bases and comparisons with plane-wave calculations. It reproduces the composition-dependent trends, while finite-temperature entropy and neglected defect interactions limit the phase-stability interpretation.

2. Formation energy of GaAs1−xNx: The formation-energy expression treats nitrogen impurities as point defects and neglects interactions between neighboring nitrogen defects.
2. Formation energy of GaAs1−xNx: Under As-rich conditions, nitrogen defect formation energies are positive for all compositions, indicating very low probability of nitrogen doping into GaAs.
2. Formation energy of GaAs1−xNx: Under As-poor conditions, formation energies become negative at small nitrogen concentrations, consistent with alloy formation only near composition endpoints.
2. Formation energy of GaAs1−xNx: ABACUS DZP formation-energy curves reproduce the same concentration dependence as ABACUS plane-wave calculations for both chemical conditions.
2. Formation energy of GaAs1−xNx: DZP formation energies are underestimated by 0.2 to 0.3 eV relative to plane-wave results because the DZP N2 energy differs appreciably from the plane-wave reference.
2. Formation energy of GaAs1−xNx: Increasing nitrogen from DZP to TZDP improves formation-energy differences to within 0.1 eV of plane-wave results for all cases.
IV. SUMMARY: The paper concludes that ABACUS LCAO bases can provide reliable defect calculations while supporting large-scale simulations with plane waves and localized orbitals.
2. Formation energy of GaAs1−xNx: A thorough phase-stability treatment would require finite-temperature and entropy effects, which are outside the present work.

Appendix A: Appendix

The appendix describes efficient matrix, potential, and force evaluation techniques for ABACUS using two-center, grid, reciprocal-space, and interpolation methods.

Matrix elements: Overlap and kinetic-energy matrix elements are evaluated through tabulated two-center quantities and interpolation for arbitrary interatomic distances.
Matrix elements: Atomic-orbital overlap construction uses angular-momentum decomposition, radial Fourier transforms, and tabulated distance-dependent integrals.
Local potentials: The local potential combines local ionic, Hartree, and exchange-correlation contributions on a grid.
Local potentials: Plane waves and FFTs evaluate long-ranged local and Hartree potentials efficiently by transforming between reciprocal and real space.
Forces: Atomic-orbital forces are decomposed into Feynman-Hellmann, Pulay, nonorthogonality, and Ewald contributions.
Forces: Nonlocal pseudopotential force terms use the same two-center integral technique as atomic-orbital matrix elements.
Forces: Pulay-force derivatives combine numerical radial-orbital calculations with analytic real spherical harmonics and grid integrals.

3. The force arises from the fact that atomic orbitals are not orthogonal,

The section describes force-related terms arising from nonorthogonal atomic orbitals, including overlap-matrix derivatives and orbital separations. It also notes analytical or numerical treatments for angular and radial components.

The force expression includes the derivative of the overlap matrix with respect to atomic coordinates.The surrounding method text identifies this term as an element of the energy density matrix.
The displacement D is defined from atomic coordinates and lattice translation as D = R + τν − τµ.Here τµ and τν denote the coordinates of orbitals µ and ν, respectively.
The distance between two orbitals is used in the further expansion of the force-related term.
Multiplying Y_lm(D̂) by D^l makes the angular function analytical at the origin.
The relevant term can be evaluated numerically by interpolation or analytically with real spherical harmonic functions.