CONUNDrum: A program for orbital-free density functional theory calculations

doi:10.1016/j.cpc.2020.107365

Computer Physics Communications

Volume 256, November 2020, 107365

https://doi.org/10.1016/j.cpc.2020.107365 Get rights and content

Abstract

We present a new code for energy minimization, structure relaxation and evaluation of bulk parameters in the framework of orbital-free density functional theory (OF-DFT). The implementation is based on solving the Euler–Lagrange equation on an equidistant real space grid on which density dependent variables and derivatives are computed. Some potential components are computed in Fourier space. The code is able to use semilocal and non-local kinetic energy functionals (KEF) as well as neural network based KEFs thus facilitating testing and development of emerging machine-learned KEFs. For semi-local and machine-learned KEFs the kinetic energy potentials are evaluated with real-space differentiation of the components, which are partial derivatives of the KE with respect to the electron density, its gradient and Laplacian.

Program summary

Program title: CONUNDrum.

CPC Library link to program files: http://dx.doi.org/10.17632/phnz2gg8mz.1

Licensing provision: GNU GPL v3

Programming language: C++

External routines: Fastest Fourier Transform in the West (FFTW) library (http://www.fftw.org/)

Nature of problem: Calculation of the electronic and structural properties of molecules and extended systems in the framework of the orbital-free density functional theory. Evaluation of the bulk parameters of solid compounds.

Solution method: High-order central finite-difference method and fast Fourier transform are used for calculation of different total energy components. Density optimization is performed with the steepest descent or the Polak–Ribière variant of the non-linear conjugate-gradient method with a line search procedure based on the Armijo condition. A numerical approach is used for structural optimization — the total energies with respect to small variations in lattice geometries are computed directly, with subsequent evaluation of the force components via a high-order central-finite difference method. The same numerical procedure is used for evaluation of bulk properties.

Restrictions: Local pseudopotentials.

Introduction

Ab initio modeling is a critically important part of computational material science. One of the most popular methods is the density functional theory (DFT) [1] that remains the only practical approach for ab initio modeling of many classes of materials. This is because in DFT, the energy is expressed as a functional of three-dimensional electron density rather than of a 3 $N_{e l e c}$ dimensional wavefunction (where $N_{e l e c}$ is the number of electrons). The Kohn–Sham (KS-DFT) orbital formulation [2] of DFT is implemented in a plethora of numerical codes and is widely used in applications. Its widespread use is justified by a well-defined ansatz of non-interacting particles, which facilitates the calculation of the kinetic energy leaving only exchange and correlation contributions to be approximated. However, the necessity to compute at least as many solutions of the Kohn–Sham equation as there are electrons leads to a near-cubic scaling of the CPU cost. While near-linear scaling formulations of KS-DFT exist [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], they introduce additional approximations and achieve advantages in scaling only for large calculations (where the sparsity of matrices due to finite basis function support becomes significant) for which the cost is still significant.

In Orbital-free (OF-DFT), one dispenses with the need to compute orbitals. Instead, the energy is computed as a function only of density-dependent quantities. A key problem of OF-DFT is expression of the kinetic energy – typically done based on the Kohn–Sham decomposition of the total energy – as a functional of the density (kinetic energy functional, KEF). Specifically, with semilocal KEFs (or KEDFs, kinetic energy density functionals) [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27], [28] one achieves near-linear scaling enabling routine large scale ab initio calculations. Even with non-local KEFs, systems with dozens of thousands of atoms are routinely computable on a desktop computer. Sufficiently accurate KEFs for applied simulations are still missing for major materials classes, and KEF development will continue to be an active area of research for the foreseeable future. Specifically, machine-learnt (ML) KEFs hold significant promise [29], [30]. Availability of OF-DFT software allowing for tests and development of new types of KEFs is therefore as important as the possibility to perform applied simulations with existing functionals.

Contrary to KS-DFT, the OF-DFT formulation is only implemented in several codes [31], [32], [33], [34]. This is not surprising, since only a limited set of problems can be successfully treated with an accuracy sufficient for applications with existing OF-DFT approximations, which include simplest metallic systems (Al, Li, Mg), their alloys and some semi-conductors [35], [36], [37], [38], [39], which limited software demand. As more and more powerful approximations are proposed, the need for OF-DFT software tools is increasing, specifically those capable of working with novel types of functionals (e.g. machine-learnt).

In terms of OF-DFT codes and capabilities easily accessible to the public, PROFESS [31] is a plane-wave based implementation relying on local pseudopotentials. The advantages of PROFESS include, among others, availability of multiple non-local KEFs and a molecular dynamics capability. Real-space approaches to DFT are promising and have recently been implemented in a number of electronic structure codes [40], [41], [42], [43], [44], [45], [46], [47]; they are also becoming a method of choice in some OF-DFT codes. As an example, the OF-DFT implementation in the publicly available GPAW code [34] permits using OF-DFT with the PAW formalism;however, the choice of KEFs is limited. Other groups proposed OF-DFT implementations which are not distributed [32], [33]. A finite-element formulation was suggested in Ref. [48], [49], where also an advantage of higher-order finite-element discretization over linear finite elements in terms of computational efficiency was shown. However, in the context of OF-DFT the corresponding implementations were made only for simplest KEFs.

In the present paper, we introduce a new OF-DFT code (CONUNDrum, C++ COde for NUmerical ab iNitio analysis withorbital-free Density functional theory), which is designed for careful treatment of semi-local KEDFs (and, in particular, to address problems of numerical stability in self-consistent minimization procedures we faced with existing codes), including the capability of using ML based KEDFs. In our implementation, we combine a real-space representation [32] of the density and finite-difference derivatives with a spectral representation for some operations (such as the Hartree energy). With this scheme, we can converge density-dependent quantities requiring higher-order derivatives and powers of the density (such as terms of the 4th order gradient expansion of the kinetic energy density). Our code explicitly permits using neural-network based KEFs. The code is designed to perform density and structure optimization and calculation of bulk properties of solid compounds and relies on periodic boundary conditions. We also kept in mind the use for testing of future KEDFs and therefore tried to make the implementation as simple as possible.

Section snippets

Theoretical background

In the context of DFT, the many-body electron problem is reformulated as a search for an electron density $n (r)$ in the three-dimensional space of electron coordinates r, that minimizes the total energy of the system. OF-DFT scheme assumes that this can be done without explicit reference to the KS one-electron wavefunctions (orbitals). The energy is a functional of density-dependent variables only. The total energy is considered as a sum of total electron energy $E_{t o t}^{e l e c} [n (r)]$ and the energy ofion–ion

General workflow of the density optimization procedure

We seek a solution to the unconstrained minimization problem applying the Lagrange multiplier method: $L [n (r)] = E_{t o t} [n (r)] - μ \{\int_{Ω} n (r) d r - N_{e l e c}\}, n (r) \geq 0,$ $\frac{δ L [n (r)]}{δ n (r)} = 0 .$ where the integration is over the volume of the system $Ω$ .

The chemical potential $μ$ ensures that the number of electrons $N_{e l e c}$ is conserved. Eq. (13) is solved self-consistently by varying a trial density until the norm of the potential (the Euclidean norm is used) is smaller than a predefined tolerance. It has been argued that the use of the electron

Numerical results

All calculations were performed on bulk systems of Al, Li and Mg. Initial cell geometries (conventional cells) were taken from experimental geometries (fcc Al: a $=$ 7.652 a.u.; bcc Li: a $=$ 6.633 a.u.; hcp Mg a $=$ 6.065 a.u., c $=$ 9.847 a.u.) [85]. For tested GGA-type KEDFs, the initial density state was precalculated with Wang–Teter (WT) KEF [57] (this is allowed by a special keyword in CONUNDrum), since the start from this state as initial yields a bit lower final total energies than the start

Conclusions

We developed a code for a numerical solution of theOF-DFT problems with enhanced possibilities for testing of new KE approximations. Special attention was paid to “problematic” semi-local KE functionals, which can be unstable during electron density minimization. Here, with the right choice of minimization parameters for these KE functionals the stability in density minimization procedure can be ensured up to 10 $^{- 3} eV$ per atom, in cell relaxation procedure up to 10 $^{- 3}$ a.u. for lattice parameters

Acknowledgment

This work was in part supported by the Ministry of Education of Singapore via an AcRF Tier 1 grant.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References (97)

HavuV. et al.
J. Comput. Phys.
(2009)
HineN.D.M. et al.
Comput. Phys. Comm.
(2009)
PerdewJ.P.
Phys. Lett. A
(1992)
HoG.S. et al.
Comput. Phys. Comm.
(2008)
MiW. et al.
Comput. Phys. Comm.
(2016)
ShaoX. et al.
Comput. Phys. Comm.
(2018)
GhoshS. et al.
Comput. Phys. Comm.
(2017)
GhoshS. et al.
Comput. Phys. Comm.
(2017)
Michaud-RiouxV. et al.
J. Comput. Phys.
(2016)
ValievM. et al.
Comput. Phys. Comm.
(2010)

MotamarriP. et al.

Comput. Phys. Comm.

(2020)

GaviniV. et al.

J. Mech. Phys. Solids

(2007)

MotamarriP. et al.

J. Comput. Phys.

(2012)

PulayP.

Chem. Phys. Lett.

(1980)

ToukmajiA.Y. et al.

Comput. Phys. Comm.

(1996)

ReshakA.H. et al.

J. Alloy Compd.

(2012)

ReshakA.H. et al.

Int. J. Electrochem. Sci.

(2013)

JamalM. et al.

Comput. Mater. Sci.

(2014)

WangJ. et al.

Int. J. Plast.

(2014)

SabraM.K. et al.

J. Alloy Compd.

(2017)

KarasievV.V. et al.

Adv. Quantum Chem.

(2015)

KarasievV.V. et al.

Phys. Rev. B

(2009)

EngelE. et al.

Phys. Rev. B

(1991)

HohenbergP. et al.

Phys. Rev.

(1964)

KohnW. et al.

Phys. Rev.

(1965)

GalliG. et al.

Phys. Rev. Lett.

(1992)

LiX.-P. et al.

Phys. Rev. B

(1993)

OrdejónP. et al.

Phys. Rev. B

(1995)

WangY. et al.

Phys. Rev. Lett.

(1995)

ArtachoE. et al.

Phys. Status Solidi b

(1999)

ScuseriaG.E.

J. Phys. Chem. A

(1999)

GoedeckerS.

Rev. Modern Phys.

(1999)

BowlerD.R. et al.

J. Phys. Condens. Matter.

(2002)

RudbergE. et al.

J. Chem. Theory Comput.

(2011)

MohrS. et al.

Phys. Chem. Chem. Phys.

(2015)

KirzhnitsD.A.

Sov. Phys. JRTP

(1957)

HodgesC.H.

Can. J. Phys.

(1973)

MurphyD.R.

Phys. Rev. A: At. Mol. Opt. Phys.

(1981)

Ou-YangH. et al.

Int. J. Quantum Chem.

(1991)

TranF. et al.

Int. J. Quantum Chem.

(2002)

KarasievV.V. et al.

J. Comput.-Aided Mater. Des.

(2006)

PerdewJ.P. et al.

Phys. Rev. B

(2007)

LaricchiaS. et al.

J. Chem. Theory Comput.

(2014)

KarasievV.V. et al.

Phys. Rev. B

(2013)

ConstantinL.A. et al.

J. Chem. Theory Comput.

(2017)

SmigaS. et al.

J. Chem. Phys.

(2017)

LuoK. et al.

Phys. Rev. B

(2018)

ConstantinL.A. et al.

J. Phys. Chem. Lett.

(2018)

Cited by (20)

2D MOF-enhanced SPR sensing platform: Facile and ultrasensitive detection of Sulfamethazine via supramolecular probe
2023, Journal of Hazardous Materials
Sulfamethazine (SMZ) is widely present in the environment and can cause severe allergic reactions and cancer in humans. Accurate and facile monitoring of SMZ is crucial for maintaining environmental safety, ecological balance, and human health. In this work, a real-time and label-free surface plasmon resonance (SPR) sensor was devised using a two-dimensional metal-organic framework with superior photoelectric performance as an SPR sensitizer. The supramolecular probe was incorporated at the sensing interface, allowing for the specific capture of SMZ from other analogous antibiotics through host-guest recognition. The intrinsic mechanism of the specific interaction of the supramolecular probe-SMZ was elucidated through the SPR selectivity test in combination with analysis by density functional theory, including p-π conjugation, size effect, electrostatic interaction, π-π stacking, and hydrophobic interaction. This method facilitates a facile and ultrasensitive detection of SMZ with a limit of detection of 75.54 pM. The accurate detection of SMZ in six environmental samples demonstrates the potential practical application of the sensor. Leveraging the specific recognition of supramolecular probes, this direct and simple approach offers a novel pathway for the development of novel SPR biosensors with outstanding sensitivity.
Highly accurate machine learning model for kinetic energy density functional
2021, Physics Letters, Section A: General, Atomic and Solid State Physics
Citation Excerpt :
Burke and coworkers original paper [18] was based on kernel regression and many followed the same approach [22–24]. Also, recently, neural network was applied [22,25,26]. For example, an interesting approach was recently proposed by Nagai et al. [22] to adopted generalized linear expansion to train neural network model.
Inspired by the remarkable ongoing progress of data-driven approaches, a very accurate predictive model is developed to estimate one-dimensional kinetic energy density functionals (KEDF) using Machine Learning (ML). Starting from possible analytical forms of kinetic energy density and by utilizing a variety of solvable models, a simple – yet highly – accurate linear regression model is statistically trained to estimate the kinetic energy as functionals of the density. The mean relative accuracy for even a small number of randomly generated potentials is found to be better than the standard KEDF (Thomas-Fermi (TF) and von Weizsäcker (vW)) by several orders of magnitudes. As more different potentials of model problems are mixed, the coefficients of the linear model significantly approach the known values of Thomas-Fermi and von Weizsäcker, suggesting the reliability of the statistical training approach. This work can provide an important step toward more accurate large-scale orbital free density functional theory (OF-DFT) calculations.
Nonlocal pseudopotential energy density functional for semiconductors
2024, Physical Review B
Machine learning of kinetic energy densities with target and feature smoothing: Better results with fewer training data
2023, Journal of Chemical Physics
Atoms, dimers, and nanoparticles from orbital-free density-potential-functional theory
2023, Physical Review A
Orbital-Free Density Functional Theory: An Attractive Electronic Structure Method for Large-Scale First-Principles Simulations
2023, Chemical Reviews

View all citing articles on Scopus

^☆: The review of this paper was arranged by Prof. D.P. Landau.

^☆☆: This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect (http://www.sciencedirect.com/science/journal/00104655)

View full text

CONUNDrum: A program for orbital-free density functional theory calculations☆,☆☆

Abstract

Program summary

Introduction

Section snippets

Theoretical background

General workflow of the density optimization procedure

Numerical results

Conclusions

Acknowledgment

Declaration of competing interest

J. Comput. Phys.

Comput. Phys. Comm.

Phys. Lett. A

Comput. Phys. Comm.

Comput. Phys. Comm.

Comput. Phys. Comm.

Comput. Phys. Comm.

Comput. Phys. Comm.

J. Comput. Phys.

Comput. Phys. Comm.

Comput. Phys. Comm.

J. Mech. Phys. Solids

J. Comput. Phys.

Chem. Phys. Lett.

Comput. Phys. Comm.

J. Alloy Compd.

Int. J. Electrochem. Sci.

Comput. Mater. Sci.

Int. J. Plast.

J. Alloy Compd.

Adv. Quantum Chem.

Phys. Rev. B

Phys. Rev. B

Phys. Rev.

Phys. Rev.

Phys. Rev. Lett.

Phys. Rev. B

Phys. Rev. B

Phys. Rev. Lett.

Phys. Status Solidi b

J. Phys. Chem. A

Rev. Modern Phys.

J. Phys. Condens. Matter.

J. Chem. Theory Comput.

Phys. Chem. Chem. Phys.

Sov. Phys. JRTP

Can. J. Phys.

Phys. Rev. A: At. Mol. Opt. Phys.

Int. J. Quantum Chem.

Int. J. Quantum Chem.

J. Comput.-Aided Mater. Des.

Phys. Rev. B

J. Chem. Theory Comput.

Phys. Rev. B

J. Chem. Theory Comput.

J. Chem. Phys.

Phys. Rev. B

J. Phys. Chem. Lett.