JDFTx: software for joint density-functional theory

Density-functional theory (DFT) has revolutionized computational prediction of atomic-scale properties from first principles in physics, chemistry and materials science. Continuing development of new methods is necessary for accurate predictions of new classes of materials and properties, and for connecting to nano- and mesoscale properties using coarse-grained theories. JDFTx is a fully-featured open-source electronic DFT software designed specifically to facilitate rapid development of new theories, models and algorithms. Using an algebraic formulation as an abstraction layer, compact C++11 code automatically performs well on diverse hardware including GPUs. This code hosts the development of joint density-functional theory (JDFT) that combines electronic DFT with classical DFT and continuum models of liquids for first-principles calculations of solvated and electrochemical systems. In addition, the modular nature of the code makes it easy to extend and interface with, facilitating the development of multi-scale toolkits that connect to ab initio calculations, e.g. photo-excited carrier dynamics combining electron and phonon calculations with electromagnetic simulations.


Motivation and significance
Density functional theory (DFT) enables computational prediction of material properties and chemical reactions starting from a quantum mechanical description of the electrons. DFT codes are now widely used to understand and design new materials from first principles through the prediction of electronic properties, structures and dynamics of molecules, solids and surfaces. Such studies commonly employ proprietary software such as GAUSSIAN [1] and VASP [2] as well as open-source software such as Quantum Espresso [3], ABINIT [4] and Qbox [5], to name just a few. 1 However, DFT offers limited accuracy for certain classes of materials and properties [6], and is extremely computationally expensive for amorphous materials, liquids and nanostructures [7]. The study of new systems, as soon as they become technologically and scientifically relevant, requires continual development of new methods to improve Listing 1: JDFTx code for a simplified solvation model with local dielectric screening. The extended DFT++ algebraic formulation [11,12] expresses the physical model in almost the same language as the defining equations (1) and (2), while parallel implementations of the operators make this code, written only once, run with multiple CPU threads (using pthreads) or on NVIDIA GPUs (using CUDA).
class LinearPCM : public LinearSolvable<ScalarFieldTilde> { ScalarField epsilon; //inhomogeneous dielectric public: virtual ScalarFieldTilde hessian(const ScalarFieldTilde& phiTilde) const { return (-1./(4*M_PI)) * divergence(J(epsilon * I(gradient(phiTilde)))); } double getAdiel(const ScalarFieldTilde& rhoTilde) const { solve(rhoTilde); //solve hessian(state) = rhoTilde return 0.5 * dot(state -coulomb(rhoTilde), O(rhoTilde)); } }; tion 2 presents the overall design of JDFTx using the algebraic formulation of DFT [11,12] to separate implementation into physics, algorithm and hardware layers, enabling rapid development of high-performance code that is easy to use. It also outlines commonly used features of the code, some of which are illustrated in more detail with examples in section 3. Finally, section 4 highlights new methods that have already been developed using JDFTx, including a hierarchy of JDFT models for the electronic structure of solvated systems and a toolkit for photo-excited carrier dynamics with ab initio electron and phonon properties, as well as key applications of these methods.

Software description
The core functionality of any electronic DFT software includes the calculation of ground-state electron densities, energies and forces within the Kohn-Sham DFT formalism [13], given a list of atoms and their positions. This facilitates prediction of structure and dynamics of materials, evaluation of reaction pathways and chemical kinetics, as well as determination of phase equilibria and stability.
Due to the nature of quantum-mechanical simulations of matter, DFT calculations become increasingly expensive with the number of atoms and electrons involved; computational complexity ranges from O(N 3 ) (with a smaller prefactor) to O(N ) (with a much larger prefactor), depending on the implementation. A brute force approach to nano and mesoscale systems with several thousands to millions of atoms is therefore not practical; it is instead logical to develop multi-scale theories for the properties of interest while still incorporating DFT electronic structure where appropriate.
JDFTx is an open-source DFT software designed specifically with the goals of coupling electronic DFT with coarsegrained theories to bridge atomic and system length scales, and of facilitating the rapid development of new classes of such combined theories. It implements electronic DFT in the plane-wave basis, which is best suited for periodic systems such as solids and solid surfaces, but is also applicable to molecular systems. A key functionality of JDFTx beyond standard electronic DFT codes is the modeling of liquids using classical DFT [12], and JDFT calculations of electronic structure in liquid environments by combining electronic DFT with classical DFT or simpler solvation models. Section 2.1 presents a birds-eye view of the code architecture along with a code example to illustrate the ease of developing new features, after which section 2.2 outlines the key features of the code.

Software Architecture
JDFTx achieves its goal of code simplicity and extensibility by using the 'DFT++' algebraic formulation of electronic DFT [11] and its generalization to classical DFT and JDFT [12]. This algebraic formulation cleanly separates the code into physics, algorithm and computational layers. Theories and algorithms are expressed concisely at a high-level of abstraction in the top layers of the code, whereas performance optimizations and support for specialized hardware such as GPUs are handled in the lower layers.
We illustrate this clean separation with an example of a simplified solvation model defined by −4πρ( r) = ∇ · ( ( r)∇φ( r)), and (1) Here, the liquid is treated as an inhomogeneous dielectric ( r) which interacts with the charge density of the electronic system, ρ( r). The net electrostatic potential φ( r) satisfies the modified Poisson equation (1), and the electrostatic solvation energy A diel is the difference between the dielectric-screened and unscreened electrostatic self energies of ρ( r) (2), whereK is the unscreened Coulomb operator. This is the essence of most solvation models used with DFT [14][15][16][17]; we have only skipped the determination of ( r) from atom positions or electron densities, and additional non-electrostatic correction terms in A diel for brevity. Regardless, solving the Poisson equation above is the most complex and time-consuming operation in these solvation models.
Listing 1 shows the implementation of this model in JDFTx. Class LinearPCM derives from a templated abstract base class LinearSolvable, which implements the Conjugate Gradients (CG) algorithm [18] on arbitrary vector spaces, instantiated in this case for scalar fields in reciprocal space, ScalarFieldTilde. The equation to be solved is defined by the virtual function hessian, whose one-line implementation can be recognized as (1) at a moment's glance. Note that the gradient (∇) and divergence (∇·) operators apply in reciprocal space (where they are diagonal), while the operators I and J Fourier transform from reciprocal to real space and vice versa [ [11]] (using Fast Fourier Transforms). The function getAdiel first calls function solve from base class LinearSolvable to solve (1) for φ( r) stored in a member variable state of the base class, and the second line evaluates A diel using (2). The integral is evaluated as a dot product with the overlap operator O, and coulomb implements the unscreened Coulomb operatorK, both of which are diagonal in the plane-wave basis [ [11]].
The algebraic formulation enables Listing 1 to resemble (1) and (2) as closely as possible, streamlining the translation of physics into code. In addition, the linear algebra representation allows the implementation of all involved operators (eg. I,J, *, gradient etc.) to be optimized under the hood for different hardware. In particular, all operators in JDFTx are already implemented using both pthreads for multi-core CPUs (Central Processing Units) and CUDA for NVIDIA GPUs. This division of labor allows all physics code in JDFTx to be implemented only once, yet run natively on all supported hardware configurations. Consequently, JDFTx could compute entirely on GPUs since its inception in 2012 (a first for plane-wave DFT).
Additionally, the use of C++11 features such as rvalue references and smart pointers in the implementation of data structures (eg. ScalarField) and their operators helps minimize memory overhead. For example, in hessian() in Listing 1, I(), epsilon * and J() automatically operate in-place because their inputs are temporary return values of other functions. Fig. 1 shows the organization of the key classes and files of the JDFTx codebase, containing 275 source files with approximately 60000 lines of code; this fairly compact implementation for a DFT code is possible because of the algebraic framework discussed above. Most of the functionality is compiled into a dynamic link library libjdftx, which is accessed through light-weight executables jdftx for most calculations, phonon for phonon dispersion and electron-phonon interaction calculations, and wannier for generation of maximally-localized Wannier functions [19] and ab initio tight binding models. This organization makes the procedure straightforward for other DFT codes to leverage JDFTx features, especially JDFT and related solvation models, by linking directly to libjdftx.
The core code (Fig. 1) implements the basic data structures, operators and algorithms used by the rest of JDFTx. ManagedMemory handles data transfers between CPUs and GPUs (automatically, when necessary), and is the base class for ScalarField s and electronic wavefunctions in ColumnBundle. Minimize and Pulay provide algorithms for variational minimization, including linear and nonlinear CG [18], L-BFGS [20]) and self-consistency by Pulay mixing [21]. GridInfo describes the plane-wave grid and Fourier transforms, while Coulomb implements the Coulomb kernel in various dimensions [22].
The fluid code is tied together by the abstract base class FluidSolver, which is implemented by FluidMixure containing several FluidComponents to provide the classical DFTs [12,23] used for JDFT; see Ref. [12] for more details on this framework to implement atomically-detailed classical DFT for molecular fluids in terms of IdealGas representations and Fex (excess) functionals. Alternately, a hierarchy of solvation models with linear (LinearPCM [24]), nonlinear (NonlinearPCM [16]) and even nonlocal response (SaLSA [25]) derive from base class PCM.
The electronic code contains the standard Kohn-Sham electronic DFT implementation (with class names derived from  the original implementation of DFT++ [11]). Here, ElecInfo and ElecVars contain electronic occupations, wavefunctions, density, potential etc. and the functions relating them. SpeciesInfo handles electron-ion interactions (pseudopotentials), ExCorr implements exchange-correlation functionals, and Dump handles outputs and post-processing. All this functionality (including fluids) is tied together into the container class Everything for convenience.
Finally, the commands code ( Fig. 1) provides an objectoriented interface to define commands and parse input files. All executables use this module to provide a consistent input file syntax that supports environment variable substitution and modular input files; phonon and wannier support all jdftx commands in addition to those specific to phonon dispersion and Wannier function calculations respectively. Complete description of the classes, files and connections between them is available in the API documentation generated automatically by Doxygen [26] on the JDFTx website, http://jdftx.org. Table 1 lists a selection of features available in JDFTx. It supports the full range of functionality found in all major electronic DFT software. It has built-in support for several semilocal [27][28][29], meta-GGA [30] and EXXhybrid [31][32][33] exchange-correlation functions, with additional functionals available through LibXC [34]. DFT+U [35] and DFT-D2 [36] pair potential dispersion corrections allow for handling localized electrons and van der Waals interactions respectively. JDFTx supports norm-conserving and ultrasoft pseudopotentials in the FHI, UPF and USPP formats, which allows easy interoperability with QE [3] and ABINIT [4]. It automatically installs two well-tested, open-source pseudopotential libraries, GBRV (ultrasoft) [37] and SG15 (norm-conserving) [38], enabling out-of-thebox calculations for most elements. JDFTx supports interactions with custom external potentials and fields, and allows accurate calculations of systems of any dimensionality: molecules (0D), wires (1D), slabs / 2D materials and bulk (3D), using truncated Coulomb interactions [22].

Software Functionalities
Importantly, JDFTx implements two distinct classes of algorithms for electronic DFT, variational minimization of the (analytically-continued) total energy [39,40] and the self-consistent field (SCF) iteration method [41]. Variational minimization is stable and guaranteed to converge, while SCF (the default method available in all DFT codes) is less stable in general, but faster when it converges well. JDFTx also uniquely implements grand-canonical DFT [42], where electron number adjusts automatically at a fixed electron chemical potential, which correctly describes the behavior of electrochemical systems. SCF convergence can be problematic in this mode, and for many solvated systems in general; variational minimization is therefore essential to have as an alternative.
JDFTx specializes in electronic structure calculations incorporating a continuum description of the environment, with a the range of techniques available for including fluids and solvation effects. Full JDFT calculations include a detailed classical DFT description of the liquid that captures atomic-scale structure in a number of solvents [12,23,43]. JDFTx also includes a hierarchy of solvation models for computationally inexpensive treatment of solvation effects, ranging from the nonlocal solvation model SaLSA [25] which captures atomic-scale liquid structure at a linear-response level, through models that capture nonlinear dielectric and ionic response [16], to the simplest linear-response models [24]. These linear response models include GLSSA13 [16] (later ported to VASP as VASPsol [17]), SCCS [15] (the solvation model available in QE) and CANDLE [44]. Table 2 compares the accuracy of these solvation models for a standard set of neutral solutes, cations and anions in water, and shows that CANDLE achieves the best accuracy for aqueous solvation of charged species by explicitly capturing charge asymmetry in solvation.   [19,46] and transform Hamiltonians and matrix elements into an ab initio tight-binding model. JDFTx is also interfaced with several commonly-used visualization software packages, and with the Atomistic Simulation Environment [47] for features including Nudged Elastic Band (NEB) [48] barrier calculations and alternate molecular dynamics methods. It provides solvation functionality to other electronic structure software through interfaces, such as for Quantum Monte-Carlo (QMC) simulations in CASINO [49].

Illustrative Example
JDFTx is used exactly like most plane-wave DFT software, with an input file that describes the atomic geometry, pseudopotentials, functionals, and other similar options alongside the operations to be performed. For example, Listing 2 shows an input file for a formate ion on a Pt(111) surface, modeled as a 3 layer slab in a 2×2 supercell, with the slab normal along the third lattice vector. The first section selects the GBRV pseudopotentials [37] for all elements (using the wildcard '$ID') and corresponding kinetic energy cutoffs (all energies in Hartrees (E h )) for wavefunctions and charge densities. The second section specifies the lattice geometry, slab-mode Coulomb truncation to isolate periodic images along z, ionic positions in lattice coordinates, and requests 10 iterations of ionic geometry optimization. The final 0 on the Pt atom positions constrains their positions completely, the 1 allows the rest The third section of Listing 2 selects Brillouin zone sampling, smearing, and a unique feature of JDFTx: grand canonical DFT [42] at a fixed electron chemical potential µ = −0.16 E h . The fourth section selects the CANDLE solvation model [44] for water with 1M (mol/L) Na+ and F-ions. The final section requests output every ionic step of ionic positions, electron density and bound charge density in the fluid, with files named using the pattern 'test.*'. The organization above is only for illustration; commands may appear in any order, may be split over multiple input files with 'include' statements, and may invoke environment variable substitution for ease of scripting. See http://jdftx.org/Commands.html for more details.
After saving Listing 2 to 'test.in' (as one possible naming convention), mpirun -n 4 jdftx -c8 -i test.in -o test.out   runs the code using 4 MPI processes with 8 threads each, logging output to 'test.out'. After completion, createXSF test.out test.xsf n creates file 'test.xsf' containing the structure and electron density n( r), which visualized in VESTA [50] yields Figure 2(a). The same procedure with 'nbound' instead of 'n' yields Figure 2(b) that shows the bound charge density in the fluid. As expected, the formate electron density (green) is greatest on the O atoms, surrounded by a positive (blue) fluid bound charge. A small negative (red) bound charge appears next to the H, but at this potential the Pt surface is negatively charged and mostly surrounded by positive charge in the fluid. (The raw output from JDFTx uses the opposite electron-is-positive sign convention rather than the usual convention employed in this discussion.) Figure 3 illustrates the performance of JDFTx using the above example for both (a) conventional vacuum DFT cal-culations (performed automatically before solvated calculations) and (b) fixed-potential solvated calculations, on varying numbers of both CPUs and GPUs. JDFTx implements MPI parallelization over symmetry-reduced kpoints and pthreads / CUDA parallelization over everything else. It exhibits almost linear MPI scaling when the number of processes is a factor of this k-point count (20 in above example). Each K80 GPU delivers approximately 3× performance compared to each 8-core Xeon approximately (25× compared to each core) for this problem size.

Impact
JDFTx simultaneously targets two partially-overlapping research communities: developers of new DFT methods, models, and algorithms, and DFT practitioners who take advantage of these new methods. It emphasizes ease of development as well as use, realized in practice with highlymodular code that expresses physics in almost the same language as the theoretical derivations, using the algebraic formulation of DFT and an efficient C++11 abstraction, as discussed above. Using this framework, we have rapidly implemented in a few years a complete feature set comparable to that of proprietary codes that have been in development for decades, making JDFTx now usable as a general-purpose DFT software.
JDFTx has facilitated the rapid development of new methods for a number of applications, most notably for DFT calculations of solvated and electrochemical systems using JDFT [10] and solvation models. JDFTx served as the medium for development of liquid free energy functionals (classical DFT functionals) [12,23,43] to realize JDFT, as well as a hierarchy of increasingly accurate solvation models including simple linear solvation models [16,24,51], nonlinear models [16,52,53] and models incorporating nonlocality: SaLSA [25] and CANDLE [44]. In fact, one of the simplest models developed using JDFTx, GLSSA13 [16], was later ported to the proprietary DFT code VASP as VASPsol [17], the predominant solvation option for that code. Solvation from JDFTx can be used for quantum Monte Carlo simulations through an interface with the CASINO code [49,54]. JDFTx recently enabled the development of grand-canonical DFT algorithms [42], making it possible to realistically simulate electrochemical processes by allowing number of electrons in the calculation to adjust automatically at fixed potential. More generally, JDFTx has also been used in electronic structure method development for exact exchange [22], dielectric matrices [55], X-ray measurements [56],and electron-phonon interactions [57].

Conclusions
We present JDFTx, a general-purpose open-source planewave DFT software, with a particularly rich feature set for solvated and electrochemical DFT calculations, and with an emphasis on ease of development and use. In its first five years, JDFTx has enabled rapid development of joint density-functional theory (JDFT) and a hierarchy of efficient and accurate solvation models, which are gradually being implemented or interfaced with proprietary and other open-source DFT codes. With these methods, JDFTx has already been widely applied to electrochemical problems in catalysis, energy conversion and energy storage.
The design of JDFTx allows short, easy-to-read code to perform well on various hardware architectures, making it ideally suited to rapidly prototype new methods, and then test and optimize them. All functionality is exposed as a library, libjdftx, that other software can link to for direct interfacing. Future interfaces of JDFTx with other opensource software for electronic structure calculations, direct computation of experimental observables, electromagnetic simulations, phase-field methods, and more will drive widespread applications of multi-scale techniques that start at the electronic scale.