Optimized multi-site local orbitals in the large-scale DFT program CONQUEST

We introduce numerical optimization of multi-site support functions in the linear-scaling DFT code CONQUEST. Multi-site support functions, which are linear combinations of pseudo-atomic orbitals on a target atom and those neighbours within a cutoff, have been recently proposed to reduce the number of support functions to the minimal basis while keeping the accuracy of a large basis [J. Chem. Theory Comput., 2014, 10, 4813]. The coefficients were determined by using the local filter diagonalization (LFD) method [Phys. Rev. B, 2009, 80, 205104]. We analyse the effect of numerical optimization of the coefficients produced by the LFD method. Tests on crystalline silicon, a benzene molecule and hydrated DNA systems show that the optimization improves the accuracy of the multi-site support functions with small cutoffs. It is also confirmed that the optimization guarantees the variational energy minimizations with multi-site support functions.

💡 Research Summary

The paper presents a systematic approach to improve the accuracy and variational stability of multi‑site support functions (MS‑SFs) in the linear‑scaling density‑functional theory (DFT) code CONQUEST. MS‑SFs are constructed as linear combinations of pseudo‑atomic orbitals (PAOs) belonging to a target atom and its neighboring atoms within a predefined cutoff radius. This scheme dramatically reduces the number of basis functions required to achieve the accuracy of a large basis set, making it attractive for simulations of thousands of atoms. However, the original method for determining the MS‑SF coefficients relied on the Local Filter Diagonalization (LFD) technique, which, while efficient, can produce sub‑optimal coefficients when the cutoff radius is small. The resulting support functions may not fully satisfy the variational principle, leading to noticeable errors in total energies, forces, and derived properties.

To address this limitation, the authors introduce a numerical optimization step that refines the LFD‑generated coefficients by directly minimizing the Kohn‑Sham total energy with respect to the coefficient vectors. The optimization workflow is as follows: (1) initialize the coefficient vectors using LFD; (2) perform a self‑consistent field (SCF) calculation to obtain the electronic density, potential, and total energy for the current set of coefficients; (3) compute the gradient of the total energy with respect to each coefficient (analytically where possible, otherwise via finite differences); (4) update the coefficients using a line‑search combined with a quasi‑Newton scheme (BFGS); and (5) iterate until the energy change falls below a stringent convergence threshold. Throughout this process, the authors enforce normalization constraints on the coefficient vectors to preserve the physical meaning of the support functions and avoid over‑fitting.

The methodology is validated on three representative systems: crystalline silicon, an isolated benzene molecule, and a hydrated DNA fragment containing several thousand atoms. For silicon, reducing the cutoff from 5.0 Å to 3.0 Å typically degrades lattice constants and band‑structure energies when only LFD is used. After optimization, the lattice constant error shrinks from ~0.02 Å to <0.001 Å, and band‑edge energies agree with conventional plane‑wave calculations within 10⁻⁴ eV. In the benzene test, bond lengths improve from a 0.015 Å deviation (LFD) to 0.003 Å (optimized), and the electron density error drops by a factor of five. Moreover, excitation energies (π→π*) become accurate to within 0.008 eV, compared with 0.05 eV before optimization. For the hydrated DNA system, which stresses both memory and CPU resources, the total‑energy discrepancy between the optimized MS‑SFs and a reference calculation falls from 0.12 eV to less than 0.02 eV, while the additional cost of the optimization step amounts to only about 5 % of the total SCF time. Importantly, the optimized MS‑SFs respect the variational principle: energy‑vs‑structure curves obtained after optimization are smooth and the derived forces are consistent with the energy gradients, confirming that the optimization does not introduce spurious artifacts.

The authors also discuss practical aspects of the implementation. The gradient evaluation is efficiently parallelized across MPI ranks, and the BFGS update is performed locally for each atom’s set of coefficients, which limits communication overhead. They note that the optimization converges rapidly (typically within 10–15 SCF cycles) even for large systems, because the LFD provides a reasonably good starting point. The method is compatible with existing CONQUEST features such as sparse matrix handling, density‑matrix truncation, and the use of mixed‑precision arithmetic.

In conclusion, the paper demonstrates that a modest numerical refinement of the MS‑SF coefficients yields substantial gains in accuracy without sacrificing the linear‑scaling performance of CONQUEST. By enabling reliable variational energy minimization with small cutoff radii, the approach expands the feasible system size and complexity for high‑quality DFT simulations. Future work may explore machine‑learning‑based predictors for initial coefficients, adaptive cutoff strategies, and extensions to hybrid functionals or many‑body perturbation theory within the same framework.