KS-pies: Kohn-Sham Inversion Toolkit
KKohn-Sham Inversion Toolkit
Seungsoo Nam, Ryan J. McCarty, Hansol Park, and Eunji Sim ∗ Department of Chemistry, Yonsei University, 50 Yonsei-ro Seodaemun-gu, Seoul 03722, Korea Departments of Chemistry, University of California, Irvine, CA 92697, USA (Dated: August 21, 2020)
A Kohn-Sham (KS) inversion determines a KS potential and orbitals corresponding to a given electrondensity, a procedure which has applications in developing and evaluating functionals used in densityfunctional theory. Despite the utility of KS inversions, the application of these methods among theresearch community is disproportional small. We implement the KS inversion methods of Zhao-Morrison-Parr and Wu-Yang in a framework which simplifies analysis and conversion of the resulting potentialin real-space. Fully documented python scripts integrate with PySCF, a popular electronic structureprediction software, and alternative Fortran code is provided for computational hot spots.
1. PROGRAM SUMMARY
Program Title:
Kohn-Sham Python Inversion EvaluationSoftware
Developer’s respository link: https://doi.org/10.25351/V3.KS-PIES.2020
Licensing provisions:
Apache-2.0
Programming language:
Python3, Fortran90
Nature of problem:
Determining a KS potential from a givendensity requires a KS inversion, a method for mapping anelectron density to a KS potential which can reproduce thisdensity. A unique mapping between a KS potential and anelectron density does not exist within a localized basis sets.Practical calculation of KS potentials requires an approximatemethod accepting of previously determined electron densities.
Solution method:
We implement the theoretical methodologyof Zhao-Morrison-Parr [41] and Wu-Yang [39] in Python tocalculate KS potentials from provided densities. One-particledensity matrix obtained from PySCF calculation on the molecularsystem can be passed into our program to perform KS inversion.Our Wu-Yang implementation allows KS inversion of a user-defined quantum system. A utility module helps to prepareKS inversion inputs (density matrix pre-processing) or analyzinginversion output (evaluating KS potentials in real-space).
2. INTRODUCTION
Future progress in improving Kohn-Sham (KS) densityfunctional theory (DFT), one of the most popularcomputational techniques for materials and molecules,depends upon improved exchange-correlation (XC)functionals [6]. KS DFT assumes that non-interactingelectrons and a local multiplicative potential (i.e. KSpotential) approximate the electron density of the real ∗ [email protected] system’s interacting electrons. No systematic method forimproving XC functionals has been identified within theKS system. Instead, functional development draws upon abroad array of techniques, methods, and data sources forguidance on how to produce more accurate approximations.One resourceful method is the KS inversion, whichproduces a KS potential from a provided electron density.Insight from KS inversions have been used in functionaldevelopment since 1996 [33] and has seen a recent revivalin its use today [24, 25, 42]. It has also become a usefultool for studying DFT methods, such as time-dependentDFT [3, 21], density-corrected DFT [26], inter-molecularinteractions in partition-DFT [8, 23] and embedded-DFT[2, 13, 40], and even in symmetry-adapted perturbationtheory [4].A KS inversion can construct an exact potential froman exact electron density as defined by the one-to-onedensity-to-potential mapping stated in the Hohenberg-Kohntheorem [18]. However, in practice, the finite number oflocalized basis sets needed to expand KS orbitals destroysthe advantageous one-to-one mapping [17, 19, 30], andthe problem becomes ill-posed [19]. Although an exact KSinversion is no longer possible, several approximate methods[1, 7, 10, 14, 20, 29, 34, 38–41] have been proposed.Despite the developed theory and applications of KSinversions, there are no publicly available algorithms forroutinely preforming these calculations. In this work weimplement the most frequently cited KS inversion methods,Zhao-Morrison-Parr [41] (ZMP) and Wu-Yang [39] (WY),and provide our code under an open source licence. Wealso provide a utility module that helps simplify inversioncalculations by providing pre- and post- processing and real-space evaluation.Our python implementation is compatible with PySCF[31], and utilizes NumPy [37], and SciPy libraries[36]. A section of python code can alternatively call kspies fort , a compiled Fortran version which decreasescomputational cost.We include a theoretical summary, details on theimplementation, validation, and a discussion with examplesto highlight the simplicity of incorporating KS inversionsinto modern functional development workflows. Emphasison the theoretical background seeks to provide a simplified a r X i v : . [ phy s i c s . c o m p - ph ] A ug .2 Zhao-Morrison-Parr 3 BACKGROUND and practical entry point for theorists unfamiliar with themethodology.
3. BACKGROUND3.1. General Background on Kohn-Sham DFT andits Inverse
The conventional (forward) KS procedure solves a singleparticle equation {− ∇ + v S [ n ]( r ) } ψ i ( r ) = ε i ψ i ( r ) , (1) where ψ i and ε i are i -th KS orbitals and orbital energies, v S [ n ]( r ) is the KS potential, and n is the electron density.Typically v S [ n ]( r ) is written as v S [ n ]( r ) = v ext ( r ) + v H [ n ]( r ) + v XC [ n ]( r ) , (2) where v ext ( r ) is the external potential, v H [ n ]( r ) is theHartree potential v H [ n ]( r ) = (cid:90) n ( r (cid:48) ) (cid:107) r − r (cid:48) (cid:107) d r (cid:48) , (3) and v XC [ n ]( r ) is an XC potential, which is approximate inpractice. The electron density is determined by the sum ofthe occupied orbital densities as n ( r ) = N occ (cid:88) i | ψ i ( r ) | , (4) As a funcitonal of the density, v S [ n ]( r ) is used in Eq. 1 ina self consistent loop, until a converged electron density isdetermined.KS inversions operate in reverse, using a given density n (often refered as target density, n tar ( r ) ) to determine v S .Once determined, v S and Eq. 1 produce the KS orbitals andassociated eigenvalues. KS inversions has been developedfor use with input orbitals [12, 40] or wavefunctions [28],but we focus exclusively on density based methods [1, 7,10, 14, 20, 29, 34, 38, 39, 41] which includes ZMP [41] andWY [39]. The ZMP KS inversion method [41] minimizes anobjective self-repulsion functional C [ n λ ] = (cid:90) (cid:90) [ n λ ( r ) − n tar ( r )][ n λ ( r (cid:48) ) − n tar ( r (cid:48) )] (cid:107) r − r (cid:48) (cid:107) d r d r (cid:48) , (5) by solving a KS-like equation self-consistently under a givenLagrange multiplier λ {− ∇ + v S [ n tar , n λ ]( r ) } ψ λi ( r ) = ε λi ψ λi ( r ) . (6) A λ -dependent KS-potential v S [ n tar , n λ ]( r ) = v ext ( r ) + v H [ n tar ]( r )+ v g [ n tar ]( r ) + v λC [ n tar , n λ ]( r ) , (7) includes a guiding potential v g ( r ) and correction potential v λC [ n tar , n λ ]( r ) = λ (cid:90) n λ ( r (cid:48) ) − n tar ( r (cid:48) ) (cid:107) r − r (cid:48) (cid:107) d r (cid:48) . (8) In principle, as λ → ∞ , C [ n λ ] → and n λ → n tar . Only v λC depends on λ in Eq. 7 and accommodates all necessarypotential modification. If provided, additional potentialterms v ext ( r ) , v H ( r ) , and v g ( r ) accelerate the convergenceof v S with respect to λ .The guiding potential v g ( r ) mimics the XC potential,for which a variety of potentials can be used. Typically v g ( r ) is initially formulated to mimic the asymptotic decayof XC potential, − (1 /N ) v H ( r ) , where N is the number ofelectrons. We refer − (1 /N ) v H ( r ) as FAXC, the non-Hartreeportion of the Fermi-Amaldi potential [9]. In principle, anypotential can be used for v g ( r ) when the asymptotic decayof XC potential is not important [26]. DFTXC is used todenote standard DFT XC functionals used as the guidingpotential.For open-shell systems, ZMP is used as a spin-unrestricted formalism [32]. The correction potential isspin-dependent, and Eq. 8 is rewritten for α or β spin as v λC,σ [ n tar σ , n λσ ]( r ) = 2 λ (cid:90) n λσ ( r (cid:48) ) − n tar σ ( r (cid:48) ) (cid:107) r − r (cid:48) (cid:107) d r (cid:48) , (9) where σ denotes the spin index. The factor of 2 is requiredfor consistent results in closed shell systems for restrictedand unrestricted schemes. The guiding potential is spin-dependent for DFTXC potentials, but not for FAXC. The WY approach [39] maximizes an objective functional, W S [ { b t } ] = N/ (cid:88) i (cid:90) |∇ ψ b i ( r ) | d r + (cid:90) v bS [ n tar ]( r ) { n b ( r ) − n tar ( r ) } d r , (10) where v bS ( r ) is a KS-potential similar to Eq. 7 in ZMP,except with v C ( r ) represented as a linear combination ofpotential basis functions g t ( r ) , given as v b C ( r ) = (cid:88) t b t g t ( r ) , (11) making v bS ( r ) exclusively a functional of the target density.KS orbitals ψ b i ( r ) are determined by solving a KS-likeequation, {− ∇ + v bS [ n tar ]( r ) } ψ b i ( r ) = ε b i ψ b i ( r ) , (12)2 .1 General workflow of inversion procedures 4 IMPLEMENTATION that does not require an iterative self-consistent solution.The gradient and Hessian of W S with respect to { b } is givenin an analytical form ∂W S ∂b t = (cid:90) { n b ( r ) − n tar ( r ) } g t ( r ) d r , (13) ∂ W S ∂b t ∂b u = N occ (cid:88) i N vir (cid:88) a (cid:104) ψ b a | g t | ψ b i (cid:105)(cid:104) ψ b a | g u | ψ b i (cid:105) ε b i − ε b a , (14) where N occ ( N vir ) denotes the number of occupied (virtual)orbitals, simplifying maximization of Eq. 10.In a spin-unrestricted formalism, Eq. 10 becomes W S [ { b t } ] = 12 (cid:88) σ N σ (cid:88) i (cid:90) |∇ ψ b i,σ ( r ) | d r + (cid:90) { v ext ( r ) + v H [ n tar ]( r ) }{ n b ( r ) − n tar ( r ) } d r + (cid:88) σ (cid:90) { v g,σ [ n tar , n tar σ ]( r ) + v b C,σ ( r ) }× { n b σ ( r ) − n tar σ ( r ) } d r , (15) where σ denotes spin. Spin dependence of XC potentialsis reflected in the correction ( v C ) potential due to its spindependence on { b t } and the DFTXC guiding ( v g ) potential.Regularization is included to avoid oscillatory potentialsfrom unbalanced potential basis sets, especially near thenucleus [15]. The WY objective functional is regularizedwith W η S ( { b t } ) = W S ( { b t } ) + η (cid:107)∇ v b C ( r ) (cid:107) , (16) where η is a regularization strength hyperparamater, and thesmoothness of the correction potential is measured with (cid:107)∇ v b C ( r ) (cid:107) = (cid:90) v b C ( r ) ∇ v b C ( r ) d r . (17) Heaton-Burgess and Yang 2008 [16] describes how η can beselected.
4. IMPLEMENTATION4.1. General workflow of inversion procedures
Prior to KS inversion, PySCF is used to create a
Mole object and generate a density matrix from a targetdensity. The
Mole object contains standard detailsneeded for quantum chemical calculations, such as atomiccoordinates, number of electrons, and basis sets. KS-pies requires an atomic orbital basis set representation ofthe density matrix. Some PySCF calculations, such ascoupled cluster, output density matrices in molecular orbitalrepresentation, which can be converted to atomic orbital representation using the kspies.util.mo2ao function.The present implementation does not support periodicboundary conditions.KS-pies makes use of analytical functions within thePySCF integral library to convert Eq. 6 and 12 to solvablematrix equations. An exception is made for WY calculationswhen non-equivalent atomic and potential basis sets areused. A numerical integration of the three-center overlapintegral used in Eq. 10 is calculated with S ijt = (cid:90) φ i ( r ) φ j ( r ) g t ( r ) d r , (18) where φ i ( r ) is i -th orbital basis function. Although thisenables calculations of non-equivalent basis sets, numericalintegration of Eq. 18 adds a substantial computational cost.KS-pies constructs the Hartree and guiding potentialsfrom the target density. KS-pies then uses a method specificprocedure to optimize the correction potential. For ZMP, aself consistent calculation with Eq. 6 and a user provided λ is preformed. For WY, { b t } is adjusted to maximize Eq. 10.Several features reduce the computational cost ofdetermining the correction potential. For ZMP,recalculating the Hartree potential at each self consistentiteration can be accelerated using a density fittingprocedure. Density fitting results in minor difference in theinversion, while greatly reducing the computational cost.For WY, calculating Eq. 13 and Eq. 14 using a compiledFortran module kspies fort provides a substantial timesavings over the python implementation. A compiled binaryand the source code is provided.KS-pies manages data using python instance variables,storing the majority of results in memory. Data from oneinstance can be used as an initial guess for subsequentcalculations or analysis within KS-pies or PySCF. The KS potential in real-space provides valuable insightbeyond a target density calculation. For example, amajor motivation for performing a KS inversion is thevisualization of the exact KS potential which can showthe non-intuitive step structures present in some potentials[22, 35]. However, the use of atom-centered Gaussian basis-sets prevents the XC potentials produced by ZMP and WYfrom being directly converted into real-space representation.Real-space conversion of individual component of the XCpotentials provides a practical path forward. For WY,the correction potential can be converted directly usingEq. 11. Franchini et al. [11] method for converting Hartreepotentials can be applied to the guiding potentials and theZMP correlation potential.For a given system, the basis functions { φ i ( r ) } , real-spacefunction v ( r ) and real-space function matrix representation .2 Real-Space Potentials 5 VALIDATION AND PERFORMANCE V ij have the relation v ( r ) → V ij = (cid:90) φ i ( r ) v ( r ) φ j ( r ) d r , (19a) V ij → v ( r ) = (cid:88) ij φ i ( r ) V ij φ j ( r ) . (19b) Eq. 19 is only exact under the basis set limit, and inpractice, would result in large errors in the real-spacefunction v ( r ) . Real-space conversion of the invertedpotential requires a method other than Eq. 19b.Following Franchini et al. [11], the density can bedecomposed into one-center (i.e. atomic) contributions,and into different angular contributions as n ( r ) = N nuc (cid:88) i n i ( r ) ≈ N nuc (cid:88) i l max (cid:88) l (cid:88) m Z lm ( θ i , φ i ) s ilm ( r i ) , (20) where Z represents real spherical harmonics and s is a cubicspline interpolation of radial density at the i -th atoms radialgrid. The summation inside Eq. 20 can be used to calculatethe Hartree potential for each angular contribution of i -thatom with analytical form v H ,lm ( r i ) = 4 π l + 1 Z lm ( θ i , φ i ) × (cid:16) r l +1 i (cid:90) r i r (cid:48) l +2 s ilm ( r (cid:48) ) dr (cid:48) + r li (cid:90) ∞ r i s ilm ( r (cid:48) ) r (cid:48) l − dr (cid:48) (cid:17) . (21) We implement the Franchini et al. [11] approach forconverting the Hartree potential to real-space.Eq. 20 and 21 are used to calculate point-wise values ofFAXC guiding potentials and the ZMP correction potentials(Eq. 8). The 6.5 Utility section and associated figuresprovide an example of obtaining real-space representationusing this approach.The eval vxc function evaluates XC potential on userdefined grid points. Our implementation is based on thenumerical derivative. For Local Density Approximationfunctionals (LDA) we use, v LDAXC ( r ) = δE XC δn = d(cid:15) LDAXC dn ( r ) , (22) where (cid:15) LDAXC is the XC density of LDA, and can bedirectly obtained from pysef.dft.libxc.eval xc .For Generalized Gradient Approximation (GGA) functionalswe use, v GGAXC ( r ) = v n − {∇ n · ∇ v γ + v γ ∇ n } , (23) where v n = ∂ε GGAXC ∂n , v γ = ∂ε GGAXC ∂γ , and γ = ∇ n · ∇ n . Although v n and v γ are obtainable from pysef.dft.libxc.eval xc , ∇ v γ should be evaluated using a numerical derivative of v γ . For spin-polarizeddensities, Eq. 23 for α spin becomes v GGAXC ,α ( r ) = v n α − ∇ n α · ∇ v γ αα + v γ αα ∇ n α ) − ( ∇ v γ αβ · ∇ n β + v γ αβ ∇ n β ) , (24) and the formulation for β spin requires a trivial swappingof respective spins. Eq. 24 can be used to calculate point-wise value of DFT guiding potentials. The examples of eval vxc and related plots are provided in section 6.5Utility.
5. VALIDATION AND PERFORMANCE
Valid implementation is confirmed with accurate KSinversions of densities obtained from HF or correlatedwavefunction methods, in both restricted and unrestrictedschemes. Run time benchmarks are reported as wall time,and were performed an Intel(R) Xeon(R) Gold 6142 CPUusing 8 processors at 2.6 GHz. Calls to PySCF and kspies fort can take advantage of parallelization withOpenMP. In our KS inversion examples below, benzeneused the most memory, approximately 1.2 GB. The PySCFCCSD Calculation of molecular O used substantially morememory, however the inversion using KS-pies required lessthan that of benzene. Restricted and unrestricted inversionbenchmarks are included for convenience within the softwarerepository. T i m e ( s e c ) S C F It e r a t i on s λ d N ( m e ) I n t eg r a t ed den s i t y d i ff e r en c e Z M P s e l f -r epu l s i on f un c t i ona l C [ n λ ] λ = 8 λ = 32λ = 128 λ = 512 Angstroms0-0.05 0.05
C. Density error (target HF - ZMP)B. Computational performanceA. Analytical performance
10 2milielectronsDirectIterations TimeDensity FittingDirectDensity Fitting dN (me)
C[n λ ]
32 12
32 12 λ FIG. 1. ZMP results showing that density fitting or direct(without density fitting). For each λ value, level shift was set to 0 . × λ . (a) Analytical performance evaluatedby integrated density difference in milielectron (dN) andminimized self-repulsion functional value (Eq. 5). (b)Computational performance in terms of SCF iterations andtime needed for convergence. (c) Density difference betweentarget HF density and ZMP densities for increasing λ . Restricted inversion performance was validated usingZMP and WY on benzene ( R CC = 1 . ˚ A and R CH =4 USAGE . ˚ A ). The target density was generated with HF/cc-pVTZ. Resulting inversion potentials from ZMP and WYcan be used to accurately reproduce the target density.Furthermore, ZMP, qualitatively indicates that C (Eq. 8)is approaching as λ increases. kspies.zmp.RZMP (restricted ZMP) was used with aFAXC guiding potential and tested with and without densityfitting. Using the self-repulsion functional value C fromEq. 5 and the integrated density difference dN = (cid:82) | n λ ( r ) − n tar ( r ) | d r , we confirmed an expected decrease in C as λ increased. This is plotted in Figure 1a. As λ increases, dN also decreases. Our implementation shows almost noanalytical difference between with/without density fitting.Computationally, direct and density fitting methodswithin ZMP have comparable results, except that thedensity fitting method decreases computational cost.Figure 1b plots computational performance in terms ofSCF iterations and run time per λ ; highlighting negligibledifference in analytical and SCF cycles. Per SCF iteration,ZMP takes approximately 0.25 seconds with density fittingand 3 seconds without, requiring around 700 iterations toreach convergence on a small molecules.WY performance was evaluated using the same targetdensity as above and cc-pVTZ for the potential basis set.Optimization was complete after 8 iterations, taking lessthan 1 second on default settings. The maximum gradientelement was × − with dN = 170 . me. The dN agrees with the values from ZMP at λ =128, indicatingan accurately reproduced target density. Implementationverification is also incidentally discussed in section 6.4 wherea WY reproduces a target density from a user-definedharmonic potential Hamiltonian. To validate the unrestricted calculations of ZMP andWY, we use a coupled-cluster singles-and-doubles (CCSD)target density density of molecular oxygen ( R OO = 1 . ˚ A )obtained with UHF-UCCSD/cc-pVQZ. We used FAXCguiding potential for both cases.Benchmark values for ZMP with λ =2048 are C =1.10 × − and dN =5.75 me. A small and decreasing dN as λ is increased confirms the unrestricted ZMPimplementation. As a secondary validation benchmark,identical result are produced from spin-restricted andunrestricted ZMP calculation on closed-shell benzene.Benchmark values for WY with a cc-pVQZ potentialbasis set converges after 5 optimization steps, taking 0.07seconds, with maximum gradient element 3 × − , anddN=36.3 me, which agrees with ZMP dN at λ =128,verifying the WY implementation.
6. USAGE
To encourage the use of KS inversions in functionaldevelopment, use of kspies requires only a few inputs for simple use. We provide several examples that highlightthe relative simplicity of routine ZMP and WY, calculationswith and without real space conversions. from pyscf import gto, scfimport kspiesmol = gto.M(atom=’Ne’,basis=’aug-cc-pVQZ’)mf = scf.RHF(mol).run()P_tar = mf.make_rdm1()mz = kspies.zmp.RZMP(mol,P_tar)mz.zscf(8) converged SCF energy = -128.543755937285lambda= 8.00 niter: 12 gap= 0.6327865 dN=158.02 C= 4.63e-03
FIG. 2. Example inputs for calculating a ZMP KS inversion(top) and the terminal outputs (bottom). PySCF is used togenerate a Ne density which used in a ZMP KS inversion with λ = 8 ZMP calculations proceed by instantiating a kspies.zmp.RZMP object and then calling zscf(l) which generates a KS potential with the user specified l ,( λ ), value. The instance requires mol , a Mole objectwhich defines the basis set, and
P tar , an atomic orbitalrepresentation of the target density. Calculation results areprinted to the terminal, and a matrix representation formof the potential is stored as instance attributes. This canbe converted into real-space representation with the laterdiscussed kspies.util module.Figure 2 demonstrates generating a Ne HF density usingPySCF, generating a KS potential with kspies.zmp , andthe run outputs printed to the terminal. In the output, niter is the number of SCF iterations needed to reachconvergence, gap is the HOMO-LUMO gap in atomic unit, dN is the integrated density difference in millielectrons, and C is the minimized value from Eq. 5. In unrestricted ZMP,the returned C values accounts for both spins contributionsfrom both spins, C α + C β ) . .1 Basic use of ZMP 6 USAGE mz=kspies.zmp.RZMP(mol,P_tar)mz.diis_space=30mz.max_cycle=200mz.guide=’pbe’mz.conv_tol_dm=1e-10mz.conv_tol_diis=1e-7for l in [ 8, 32, 128, 512]:mz.level_shift=l*0.1mz.zscf(l)mz.level_shift=0.mz.zscf(512)P=mz.make_rdm1()print(’Total Energy:’,mf.energy_tot(P))print(mz.mo_energy[2:6])print(’Converged?’,mz.converged,mz.l) lambda= 8.00 niter: 14 gap= 0.6250675 dN=79.38 C= 4.69e-04lambda= 32.00 niter: 18 gap= 0.6531736 dN=31.10 C= 7.06e-05lambda= 128.00 niter: 29 gap= 0.6738522 dN=12.57 C= 9.41e-06lambda= 512.00 niter: 38 gap= 0.6896071 dN=6.10 C= 1.72e-06lambda= 512.00 niter: 12 gap= 0.6896071 dN=6.10 C= 1.72e-06Total Energy: -128.5434063648391 [ -0.63464901 -0.63464901 -0.634649010.05495812 ] Converged? True 512
FIG. 3. Use of ZMP with DIIS settings at multiple λ valuesusing the mol and P tar from Figure 2 (top). Terminaloutputs detaling ZMP potential results at the specified λ values and determined quantities form the final run (bottom). Additional options in ZMP include:(i) diis space (int): DIIS space size. Default is .(ii) max cycle (int): SCF maximum iteration. Defaultis .(iii) guide (character): Guiding potential. None sets v g ( r ) = 0 , ’faxc ’ sets v g ( r ) = − N v H ( r ) , or anyDFT XC functional defined in PySCF can be specified.(iv) level shift (float): Amount of level shift usedduring SCF cycles. Default is .(v) with df (bool): Use of density fitting. Default is False . (vi) conv tol dm (float): Density matrix convergencecriteria. Default is .(vii) conv tol diis (float): Convergence criteria ofDIIS error. Default is .Results from zscf(l) are stored as the following classattributes:(i) converged (boolean): If convergence criteria wasmet during SCF iterations(ii) dm (array): Density matrix(iii) mo coeff (array): Molecular orbital coefficient(iv) mo energy (array): Molecular orbital energy(v) mo occ (array): Orbital occupation numbersThe instance only stores one set of results; previous valueswill be overwritten and only the most recent calculationsresults are accessible.On the first zscf() SCF cycle,
P tar is used asthe initial density matrix guess. The guiding potential isconstructed from user specifications during the first call tothe function. Subsequent changes to the user specificationsfor the guiding potential after the first call are ignored.The direct inversion of the iterative subspace (DIIS)procedure [27] for convergence stabilization implementedin ks-pies is independent of the PySCF DIIS options.When specifying the DIIS space size, we recommend diis space to ≤ or greater to avoid singularityof matrix during DIIS extrapolation. DIIS convergencestabilization is not used when diis space ≤
1. ZMPoften fails to converge without DIIS, and users are cautionedagainst not using this feature, even for small λ .When specified, virtual orbital energies are increased by level shift , which can aid convergence. In general,setting level shift = 0.1 × λ is likely to be sufficientfor most systems. level shift is inactive when setto 0. When λ is larger than 10 ∼ level shift values of 1 ∼ level shift reduce the mixing betweenthe occupied and virtual orbits, which slows the orbitalrotation with the intent of assisting convergence. In somesystems, the initial ZMP iterations will significantly perturbsthe potential from a trajectory towards convergence.Increasing level shift to slow the orbital rotation asshow in Figure 3 can minimize or remove this effect.As level shift is an artificial value added to aid inconvergence, its contribution is ignored when calculatingproperties such as the HOMO-LUMO gap and orbitalenergies. Performing a KS inversion with kspies.wy requiresan input
Mole object and target density, as demonstratedin Figure 4. The potential basis ( pbas ) defaults to the .2 Basic use of WY 6 USAGE atomic orbital basis, and Sijt is integrated analytically.At the conclusion of Figure 4, the KS potential producedby maximizing W S is now accessible as an instance attribute. wy_instance=kspies.wy.RWY(mol,P_tar)wy_instance.run() FIG. 4. Minimum inputs for WY calculation using apredefined mol and target density from Figure 2.
Sijt=kspies.wy.numint_3c2b(mol,’aug-cc-pV5Z’)wy_i=kspies.wy.RWY(mol, P_tar, Sijt=Sijt,pbas=’aug-cc-pV5Z’)wy_i.method=’BFGS’wy_i.guide=’blyp’wy_i.tol=1e-7for eta in [ 1e-3, 1e-4, 1e-5, 1e-6 ]:wy_i.reg = etawy_i.run()gap=wy_i.mo_energy[5]-wy_i.mo_energy[4]v=wy_i.Dvb()print(f’eta= {eta:.1e} gap: {gap:.5f}v_grad: {v:.3f}’)Sijt calculation for n1:80 n2:127 basisfunctions: 0.0 mineta= 1.0e-03 gap: 0.65590 v_grad: 1.375eta= 1.0e-04 gap: 0.68991 v_grad: 4.781eta= 1.0e-05 gap: 0.70365 v_grad: 17.231eta= 1.0e-06 gap: 0.55057 v_grad: 93.973
FIG. 5. Example kspies.wy usage (upper panel) andterminal output (lower pannel) from PySCF defined mol , mf ,and P tar as created in Figure 2. The instance inherits thepotential basis aug-cc-pV5Z (127 basis functions) from the mol object, but uses the defined atomic orbital basis set aug-cc-pVQZ (80 basis functions for Be), requiring an three centeroverlap integral array (
Sijt , with shape (80, 80, 127)) toconvert between the two.
A second WY example in Figure 5 demonstratesadvanced options and non-equivalent basis function.The kspies.wy.numint 3c2b performs numericalintegration of Eq. 18, where the input is mol object andPySCF style basis function specification. The output
Sijt array is stored as an attribute of kspies.wy.RWY or kspies.wy.UWY to specify which three-center overlapintegral is used. The terminal output (lower panel, Figure 5)shows the numint 3c2b() output, the HOMO-LUMOgap, and the smoothness of correction potential (Eq. 17)for each eta .Figure 6 calls the instance from Figure 5 for informationon the run, the number of basis functions used in thepotential, and confirmation on convergence. The instance wy_i.info()print(len(wy_i.b))print(wy_i.converged)****Optimization Completed****after 350 iterationsfunc_value : -128.53550236max_grad : 0.00000007127True FIG. 6. Following Figure 5, information on the WY runconfirms a converged KS potential. overwrites itself after each calculation and the reportedvalues are for eta = 1e-6 .Options available in kspies.wy include:(i) method (string): Optimization algorithm used inSciPy. Default is ’trust-exact’.(ii) guide (string): Guiding potential. Default is ’faxc’.Usage is same as guide in ZMP.(iii) tol (float): Tolerance of the maximum gradientvalues used to determine optimization completion(Eq. 13). Default is 1e-6.(iv) reg (float): Potential regularization weight. ( η inEq. 16) Default is 0.The molecular orbital coefficients, energies, occupationnumbers, density matrices, and convergence are all storedas instance attributes. Optimized { b t } values are stored as’ b ’, and can be useful as an initial guess for subsequentcalculations using different regularization strength η , asdemonstrated in the for loop of Figure 5. Othernotable methods include: make rdm1() for generatinga density matrix in the same format and method asthe PySCF function of the same name, info() forprinting optimization results, and Dvb() for calculating thesmoothness of the correction potential (Eq. 17).
Regularized WY can compensate for issues resulting fromunbalanced potential basis sets [15]. Figure 7 providesexample code on a N molecule, and Figure 8 displaysexample plots used to select the regularization parameter. kspies.wy.RWY requires a Mole object, a target density,and a large even tempered Gaussian basis set for thepotential calculation, all of which can be produced usingPySCF. .3 Regularized WY 6 USAGE PBS=gto.expand_etbs([(0, 13, 2**-4 , 2),(1, 3 , 2**-2 , 2)])mw=ksPIES.wy.RWY(mol, P_tar, pbas=PBS)mw.tol=1e-7etas=[ 1/2**a for a in np.arange(5,27,1)]v=np.zeros(len(etas))W=np.zeros(len(etas))for i,eta in enumerate(etas):mw.reg=etamw.run()v[i]=mw.Dvb()W[i]=mw.Wsmw.reg=0.mw.run()Ws_fin=mw.Wsimport matplotlib.pyplot as pltfig,ax=plt.subplots(2)ax[0].scatter(np.log10(Ws_fin-W),np.log10(v))ax[1].scatter(np.log10(etas),v*etas/(Ws_fin-W))
FIG. 7. Regularized WY calculations require an additionalpotential basis set and eta parameter. Selection of eta using atrial and error process as plotted in Figure 8 is recommended. log || v c || ( ) log α| | v c || b ( ) W s − W s η a)b) FIG. 8. The L-curve (a) and its reciprocal derivative (b)used for selecting eta. The plots are outputs from Figure 7and displayed here next to the formulas used in the plot. Thehorizontal axis of (a) is log( W S − W η S ) and of (b) is log( η ). Theoptimal η value is the maximum of the reciprocal derivative,Eq. 25, approximately 10 − . as visible in (b). The optimal regularization parameter, η can be identifiedusing the L-curve method of [5], where the maximumvalue of the reciprocal derivative should be selected. As η decreases, W η S approaches W S and the curvature of theKS potential increases (i.e. increasing (cid:107)∇ v b C (cid:107) ). The needfor regularization arises because an unbalanced potentialbasis set may cause (cid:107)∇ v b C (cid:107) to increases sharply, but isnot guaranteed to correspond with a significant decreasesin W S − W η S (see Figure 8a). This implies that the increaseof potential curvature results from unphysical oscillation. According to [15], an optimal η gives a maximum reciprocalslope of the L-curve, which has the analytical form ( ∂ log (cid:0) (cid:107)∇ v b C (cid:107) (cid:1) ∂ log (cid:16) W S − W η S (cid:17) ) − = η (cid:107)∇ v b C (cid:107) W S − W η S , (25) as identified by Bulat et al. [5]. Using an N molecule withthe cc-pVDZ basis and a HF target density, the reciprocalderivative of the L-curve plotted for each η in Figure 8b),highlights the usefulness of this approach in selecting anoptimal value; η ≈ − . in the given example. KS-pies can calculate a WY KS potential for user-definedHamiltonians in the unrestricted and restricted formalism.A user must provide a mol object, density matrix of targetdensity
P tar , and
Sijt array to instantiate a user definedsystem, as well as the following instance attributes: Kineticenergy T , external potential V , overlap matrix S , kineticenergy matrix of the potential basis Tp , and a None override to prevent default use of a guiding potential. mw = kspies.wy.RWY(mol,P_tar,Sijt=Sijt)mw.T = Kinetic_matrixmw.Tp = Kinetic_matrix_potential_basismw.V = Potential_matrixmw.S = Overlap_matrixmw.guide = Nonemw.run()
FIG. 9. An example input for running WY with a user-definedHamiltonian, where the user has calculated all the necessaryvariables. See the online documentation for full examples.
Figure 9 shows calculation of a potential from a user-defined Hamiltonian of a harmonic potential ( v ext ( x ) =1 / x ) with four electrons (blue curve in Figure 10).Finite-difference HF with soft electron-electron repulsion w ( x , x ) = 1 / (cid:112) ( x − x ) + 0 . is used to generate thetarget density ( P tar ) within the domain x ∈ [ − , and a grid spacing of 0.02. KS-pies online documentationexpands upon the present example with information on howto obtain necessary inputs for user-defined Hamiltonians.The resulting WY KS potential is displayed in orange inFigure 10. The KS potential includes small oscillationsnear x = ± .5 Utility 6 USAGE Target density,
P_tar
WY calculated densityWY calculated KS potentialHarmonic potential, V Position (atomic units) P o t en t i a l ( a t o m i c un i t s ) FIG. 10. A user-defined harmonic potential Hamiltonian(blue) is used to generate a WY KS potential (orange) thatcalculates a density (red) in excellent agreement with the theinput target density (black). Densities are increased vertically10x for visibility. mol = gto.M(atom=’Ne’,basis=’aug-cc-pVTZ’)mf=scf.RHF(mol).run()dm_tar = mf.make_rdm1()coords = []for x in np.linspace(0, 5, 1001):coords.append((x, 0., 0.))coords = np.array(coords)import matplotlib.pyplotmz = zmp.RZMP(mol, dm_tar)mz.guide=’faxc’for l in [ 16, 128, 1024 ]:mz.level_shift = 0.1*lmz.zscf(l)dmxc =l*mz.dm-(l+1./mol.nelectron)*dm_tarvxc = util.eval_vh(mol, coords, dmxc )pyplot.plot(coords[:, 0], vxc, label =r’$\lambda$=’+str(l))pyplot.xlim(0, 5)pyplot.ylim(-9, 0)pyplot.legend()pyplot.show()
FIG. 11. An example of eval vh function for visualizing XCpotential obtained with ZMP. FIG. 12. The plot of exchange potential of Ne atom obtainedwith code in Figure 11. The vertical axis denotes the distancefrom Ne nucleus in Bohr, and and the horizontal axis denotesexchange potential obtained from the inversion of HF densityin atomic unit.
The usage of kspies.util.eval vh functionfor evaluating XC potentials (Figure 12) ispresented in Figure 11. Referring to the example, kspies.util.eval vh requires a
Mole object (mol),density matrix to calculate Hartree potential (dmxc),and specified xyz-coordinates (coords) which describe thepositions of grid points that are to be calculated. At agiven λ , the ZMP XC potential can be written as v XC ( r ) = − N v H [ n tar ]( r ) + λ ( v H [ n λ ]( r ) − v H [ n tar ]( r ))= v H [ λn λ − ( 1 N + λ ) n tar ]( r ) , (26) indicating that the XC potential obtained from ZMP at thespecific λ , is the Hartree potential of the density λn λ − ((1 /N ) + λ ) n tar . .5 Utility 7 CONCLUSION myWY=kspies.wy.RWY(mol,P_tar,pbas=’cc-pVQZ’)myWY.method=’BFGS’myWY.guide=’pbe’myWY.tol=1e-7from pyscf import dftao2=dft.numint.eval_ao(myWY.mol2,coords)vg=helper.eval_vxc(mol,P_tar,’pbe’,coords,delta=1e-8)for eta in [1e-2, 1e-3, 1e-4, 1e-5, 1e-6]:myWY.reg=etamyWY.run()vC=np.einsum(’t,rt->r’myWY.b,ao2)plt.plot(coords[:,0],vg+vC,label=str(eta)) FIG. 13. An example script for the usage of eval vxc function for visualizing XC potential obtained with WY.FIG. 14. Plot generate by the Figure 13 displaying theexchange potential produced by WY at several η . The verticalaxis denotes the distance from the Ne nucleus in Bohr, andthe horizontal axis denotes the exchange potential obtainedfrom the inversion of HF density in atomic units. The originis excluded since a numerical derivative is not possible there. mol and P tar are as defined in Figure 11.
An example of kspies.util.eval vxc is used inFigure 13 to create the visualization of the WY XC potential in Figure 14. The finite difference required for numericaldifferentiation of v γ in Eq. 23 and Eq. 24 is set with delta (atomic units) in eval vxc , and defaults to e − .Either ZMP or WY methods can beused with kspies.utils.eval vh and kspies.utils.eval vxc . Some useful possibilitiesbeyond our current examples include using ZMP withPBE XC guiding potential, to draw PBE XC potentialwith kspies.utils.eval vxc and ZMP correctionpotential with kspies.utils.eval vh . Visualizationwith kspies.utils is not limited to XC potentialobtained from KS inversion, but can be used independentlywith KS inversion. The kspies.utils.eval vxc function satisfies online feature requests unrelated to KSinversions.
7. CONCLUSION
KS-pies presents the first publicly available code forpreforming KS inversions of electron densities into KSpotentials. Since every KS inversion method is approximate,we implemented the two most cited inversion methods,ZMP and WY. It integrates with PySCF in an environmentfamiliar to the theoretical development community. Ourframework provides a starting point for the implementationof future KS inversion methods. This publication presentsthe theoretical context and examples which highlight thesimplicity of running KS inversions. Incorporating KSinversion methods to determine real-space potentials shouldbe beneficial for XC functional development and testing.With two implemented methods, users are able tochoose and compare results, leveraging advantages andeach method. ZMP requires many SCF iterations and iscomputationally intensive relative to WY, but the result ofinversion can be systematically improved by increasing λ .Alternatively, WY can preform inversions on user-definedHamiltonians, is computationally efficient. Potentialsproduced by both methods can be converted into real spacerepresentations using our software. Acknowledgements:
S.N., H.P., and E.S. are thankfulfor support from the National Research Foundation of Korea(NRF-2020R1A2C2007468 and NRF-2020R1A4A1017737).R.J.M is thankful for support from the University ofCalifornia President’s Postdoctoral Fellowship. [1] Carl-Olof Almbladh and Antonio Carlos Pedroza. Density-functional exchange-correlation potentials and orbitaleigenvalues for light atoms.
Phys. Rev. A. , 29(5):2322,1984.[2] Mojdeh Banafsheh and Tomasz Adam Wesolowski.Nonadditive kinetic potentials from inverted kohn–sham problem.
Int. J. Quantum Chem. , 118(1):e25410, 2018.[3] P Bleiziffer, A Heßelmann, CJ Umrigar, and AndreasG¨orling. Influence of the exchange-correlation potentialin methods based on time-dependent density-functionaltheory.
Phys. Rev. A. , 88(4):042513, 2013.[4] A Daniel Boese and Georg Jansen. Zmp-sapt: Dft-sapt CONCLUSION using ab initio densities.
J. Chem. Phys. , 150(15):154101,2019.[5] Felipe A Bulat, Tim Heaton-Burgess, Aron J Cohen,and Weitao Yang. Optimized effective potentials fromelectron densities in finite basis sets.
J. Chem. Phys. ,127(17):174101, 2007.[6] Kieron Burke. Perspective on density functional theory.
J.Chem. Phys. , 136(15):150901, 2012.[7] Timothy J Callow, Nektarios N Lathiotakis, and Nikitas IGidopoulos. Density-inversion method for the kohn–shampotential: Role of the screening density.
J. Chem. Phys. ,152(16):164114, 2020.[8] Peter Elliott, Kieron Burke, Morrel H Cohen, and AdamWasserman. Partition density-functional theory.
Phys. Rev.A. , 82(2):024501, 2010.[9] Enrico Fermi and Edoardo Amaldi. Le orbite infinito deglielementi.
Accad. Ital. Rome. , 6:117, 1934.[10] Kati Finzel, Paul W Ayers, and Patrick Bultinck. A simplealgorithm for the kohn–sham inversion problem applicableto general target densities.
Theor. Chem. Acc. , 137(3):30,2018.[11] Mirko Franchini, Pierre Herman Theodoor Philipsen, Erikvan Lenthe, and Lucas Visscher. Accurate coulombpotentials for periodic and molecular systems throughdensity fitting.
J. Chem. Theory Comput. , 10(5):1994–2004, 2014.[12] Alex P Gaiduk, Ilya G Ryabinkin, and Viktor N Staroverov.Removal of basis-set artifacts in kohn–sham potentialsrecovered from electron densities.
J. Chem. TheoryComput. , 9(9):3959–3964, 2013.[13] Jason D. Goodpaster, Nandini Ananth, Frederick R. Manby,and Thomas F. Miller. Exact nonadditive kinetic potentialsfor embedded density functional theory.
J. Chem. Phys. ,133(8):084103, 2010.[14] Andreas G¨orling. Kohn-sham potentials and wave functionsfrom electron densities.
Phys. Rev. A. , 46(7):3753, 1992.[15] Tim Heaton-Burgess, Felipe A Bulat, and Weitao Yang.Optimized effective potentials in finite basis sets.
Phys.Rev. Lett. , 98(25):256401, 2007.[16] Tim Heaton-Burgess and Weitao Yang. Optimized effectivepotentials from arbitrary basis sets.
J. Chem. Phys. ,129(19):194102, 2008.[17] So Hirata, Stanislav Ivanov, Ireneusz Grabowski, Rodney JBartlett, Kieron Burke, and James D Talman. Canoptimized effective potentials be determined uniquely?
J.Chem. Phys. , 115(4):1635–1649, 2001.[18] Pierre Hohenberg and Walter Kohn. Inhomogeneouselectron gas.
Phys. Rev. , 136(3B):B864, 1964.[19] Daniel S Jensen and Adam Wasserman. Numerical methodsfor the inverse problem of density functional theory.
Int. J.Quantum Chem. , 118(1):e25425, 2018.[20] Bikash Kanungo, Paul M Zimmerman, and Vikram Gavini.Exact exchange-correlation potentials from ground-stateelectron densities.
Nat. Commun. , 10(1):1–9, 2019.[21] Jaspreet Kaur, Egor Ospadov, and Viktor N Staroverov.What is the accuracy limit of adiabatic linear-responsetddft using exact exchange–correlation potentials andapproximate kernels?
J. Chem. Theory Comput. ,15(9):4956–4964, 2019.[22] Sviataslau V Kohut, Alexander M Polgar, and Viktor NStaroverov. Origin of the step structure of molecularexchange–correlation potentials.
Phys. Chem. Chem. Phys. ,18(31):20938–20944, 2016.[23] Jonathan Nafziger and Adam Wasserman. Density- based partitioning methods for ground-state molecularcalculations.
J. Phys. Chem. A , 118(36):7623–7639, 2014.[24] Ryo Nagai, Ryosuke Akashi, Shu Sasaki, and ShinjiTsuneyuki. Neural-network kohn-sham exchange-correlationpotential and its out-of-training transferability.
J. Chem.Phys. , 148(24):241737, 2018.[25] Tomoya Naito, Daisuke Ohashi, and Haozhao Liang.Improvement of functionals in density functional theoryby the inverse kohn–sham method and density functionalperturbation theory. 52(24):245003, 2019.[26] Seungsoo Nam, Suhwan Song, Eunji Sim, and KieronBurke. Measuring density-driven errors using kohn-shaminversion.
J. Chem. Theory Comput. , 16(8):5014–5023,2020.[27] P´eter Pulay. Convergence acceleration of iterativesequences. the case of scf iteration.
Chem. Phys. Lett. ,73(2):393–398, 1980.[28] Ilya G Ryabinkin, Sviataslau V Kohut, and Viktor NStaroverov. Reduction of electronic wave functions to kohn-sham effective potentials.
Phys. Rev. Lett. , 115(8):083001,2015.[29] Ilya G Ryabinkin and Viktor N Staroverov. Determinationof kohn–sham effective potentials from electron densitiesusing the differential virial theorem.
J. Chem. Phys. ,137(16):164113, 2012.[30] VN Staroverov, GE Scuseria, and ER Davidson. Optimizedeffective potentials yielding hartree-fock energies anddensities.
J. Chem. Phys. , 124(14):141103, 2006.[31] Qiming Sun, Timothy C. Berkelbach, Nick S. Blunt,George H. Booth, Sheng Guo, Zhendong Li, Junzi Liu,James D. McClain, Elvira R. Sayfutyarova, SandeepSharma, Sebastian Wouters, and Garnet Kin-Lic Chan.Pyscf: the python-based simulations of chemistryframework.
Wiley Interdiscip. Rev.: Comput. Mol. Sci. ,8(1):e1340, 2018.[32] David J Tozer, Nicholas C Handy, and William H Green.Exchange-correlation functionals from ab initio electrondensities.
Chem. Phys. Lett. , 273(3-4):183–194, 1997.[33] David J. Tozer, Victoria E. Ingamells, and Nicholas C.Handy. Exchangecorrelation potentials.
J. Chem. Phys. ,105(20):9200–9213, 1996.[34] R Van Leeuwen and EJ Baerends. Exchange-correlationpotential with correct asymptotic behavior.
Phys. Rev. A. ,49(4):2421, 1994.[35] Robert Van Leeuwen, Oleg Gritsenko, and Evert JanBaerends. Step structure in the atomic kohn-shampotential.
Z. Phys. D. , 33(4):229–238, 1995.[36] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, MattHaberland, Tyler Reddy, David Cournapeau, EvgeniBurovski, Pearu Peterson, Warren Weckesser, JonathanBright, St´efan J. van der Walt, Matthew Brett, JoshuaWilson, K. Jarrod Millman, Nikolay Mayorov, AndrewR. J. Nelson, Eric Jones, Robert Kern, Eric Larson,CJ Carey, ˙Ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman,Ian Henriksen, E. A. Quintero, Charles R Harris, Anne M.Archibald, Antˆonio H. Ribeiro, Fabian Pedregosa, Paulvan Mulbregt, and SciPy 1. 0 Contributors. SciPy1.0: Fundamental Algorithms for Scientific Computing inPython.
Nat. Methods , 17:261–272, 2020.[37] St´efan van der Walt, S Chris Colbert, and Gael Varoquaux.The numpy array: a structure for efficient numericalcomputation.
Comput. Sci. Eng. , 13(2):22–30, 2011.[38] Yue Wang and Robert G Parr. Construction of exact kohn- CONCLUSION sham orbitals from a given electron density.
Phys. Rev. A. ,47(3):R1591, 1993.[39] Qin Wu and Weitao Yang. A direct optimization methodfor calculating density functionals and exchange–correlationpotentials from electron densities.
J. Chem. Phys. ,118(6):2498–2509, 2003.[40] Xing Zhang and Emily A Carter. Kohn-sham potentials fromelectron densities using a matrix representation within finiteatomic orbital basis sets.
J. Chem. Phys. , 148(3):034105,2018.[41] Qingsheng Zhao, Robert C Morrison, and Robert GParr. From electron densities to kohn-sham kineticenergies, orbital energies, exchange-correlation potentials,and exchange-correlation energies.
Phys. Rev. A. ,50(3):2138, 1994.[42] Yi Zhou, Jiang Wu, Shuguang Chen, and GuanHua Chen.Toward the exact exchange–correlation potential: A three-dimensional convolutional neural network construct.
J.Phys. Chem. Lett. , 10(22):7264–7269, 2019., 10(22):7264–7269, 2019.