Open source FCI code for quantum dots and effective interactions
aa r X i v : . [ c ond - m a t . m e s - h a ll ] O c t Open source FCI code for quantum dots and effective interactions
Simen Kvaal ∗ Centre of Mathematics for Applications, University of Oslo, N-0316 Oslo, Norway (Dated: October 15. 2008)We describe
OpenFCI , an open source implementation of the full configuration-interactionmethod (FCI) for two-dimensional quantum dots with optional use of effective renormalized in-teractions. The code is written in C++ and is available under the Gnu General Public License. Thecode and core libraries are well documented and structured in a way such that customizations andgeneralizations to other systems and numerical methods are easy tasks. As examples we provide amatrix element tabulation program and an implementation of a simple model from nuclear physics,in addition to the quantum dot application itself.
PACS numbers: 02.60.-x, 95.75.Pq, 73.21.La, 24.10.CnKeywords: full configuration interaction; open source; C++; quantum dots; effective interactions
I. INTRODUCTION
Quantum dots, nanometre-scale semiconductor devicesconfining a varying number of electrons, have been stud-ied intensely in the last two decades. Quantum dotsare fabricated using essentially macroscopic tools, forexample etching techniques, but the resulting confine-ment allows for quantum mechanical behaviour of theelectrons. Many of the parameters are directly control-lable, thereby justifying the term “artificial atoms” or“designer atoms”. These considerations explain the im-mense research activity on these systems. For a generalintroduction, see Ref. [1] and references therein.A very common model is that of a parabolic quantumdot, in which N electrons are confined in an isotropicharmonic oscillator potential in d spatial dimensions,where d is determined by the semi-conductor environ-ment. Electronic structure calculations on the parabolicdot and similar systems are often carried out using thefull configuration-interaction method (FCI), also calledexact diagonalization [1]. The Hamiltonian is thenprojected onto a finite-dimensional subspace of the N -electron Hilbert space and diagonalized. Care is taken inorder to exploit dynamical and discrete symmetries of theexact problem, such as conservation of angular momen-tum and total electron spin, in order to block-diagonalizethe Hamiltonian matrix and reduce the computationalcomplexity.In this article, we describe OpenFCI , a recently de-veloped open source C++ code implementing the FCImethod for quantum dots [2]. The code has a genericframework in the shape of library functions, thereby al-lowing easy customization and extension to other sys-tems and methods, e.g., three-dimensional quantum dotsor the nuclear no-core shell model.
OpenFCI implements a renormalization of the two-body interactions, a technique widely used in nuclear no-core shell model calculations. This allows for accelerated ∗ Electronic address: [email protected] convergence with respect to Slater determinant basis size[3, 4]. To the author’s knowledge, no other available codeprovides such effective interactions for quantum dot sys-tems. The code can be easily modified to create effectiveinteractions for almost many-body problem using a har-monic oscillator basis.The code is developed in a Linux environment usingGnu C++, and is readily portable to other environmentsand compilers. The Fortran 77 libraries
Lapack and
Arpack are required, but as these are available on a widerange of platforms, portability should not be affected.
OpenFCI is released under the Gnu General Public Li-cense (Gnu GPL) [5] and is documented using Doxygen[6]. As an open source project, the code can freely beused and modified.The article is organized as follows: In Section II, theFCI method is introduced in the context of the parabolicquantum dot, where we also discuss the reduction of theHamiltonian matrix by means of commuting operatorsand configurational state functions. In Section III we dis-cuss the effective two-body interaction. As the techniqueis likely to be unfamiliar to most readers outside the nu-clear physics community, this is done in some detail. InSection IV we discuss the organization and use of
Open-FCI . We also give some results from example runs, andin particular an analytically solvable non-trivial modeldue to Johnson and Payne is considered [7], where theonly modification of the parabolic quantum dot is theinteraction. Finally, we conclude our article in SectionV.Two appendices have been provided, Appendix A de-tailing the heavily-used centre-of-mass transformationand Appendix B discussing the exact numerical solutionof the two-electron quantum dot needed for the effectiveinteraction scheme.
II. FCI METHODA. Hamiltonian in occupation number formalism
We consider N electrons trapped in an isotropic har-monic oscillator potential in d spatial dimensions. Theelectrons interact via the Coulomb potential given by U ( r ij ) = λ/r ij , where r ij = k ~r i − ~r j k is the inter-particledistance and λ is a constant. The quantum dot Hamil-tonian then reads H := N X i =1 H ( i ) + N X i 440 eV · nm, and ω = ~ /m ∗ a , a beingthe trap size and length unit, and m ∗ being the effectiveelectron mass. Typical values for GaAs quantum dotsare ǫ = 12 . m ∗ = 0 . 067 electron masses, and a = 20nm, yielding λ = 2 . ~ ω , in thiscase ~ ω = 2 . 84 meV.Choosing a complete set { φ α ( x ) } α ∈ A of single-particleorbitals (where x = ( ~r, s ) denotes both spatial and spindegrees of freedom, and α = ( a, σ ) denotes both genericspatial quantum numbers a and spin projection quan-tum numbers σ = ± H can be written in occupationnumber form as H = X a,b X σ h ab a † a,σ a b,σ + 12 X abcd X στ u abcd a † a,σ a † b,τ a d,τ a c,σ , (3)where a † α ( a α ) creates (destroys) a particle in the orbital φ α ( x ). These operators obey the usual anti-commutationrelations { a α , a † β } = δ α,β , { a α , a β } = 0 . (4)For a review of second quantization and occupation num-ber formalism, see for example Ref. [8]. The single-particle orbitals are chosen on the form φ ( a,σ ) ( x ) := ϕ a ( ~r ) χ σ ( s ) , where { ϕ a ( ~r ) } are spinless orbitals and χ σ ( s ) = δ σ,s arespinor basis functions corresponding to the eigenstates ofthe spin-projection operator S z with eigenvalues σ/ { ϕ a ( ~r ) } a ∈ A are denumerable, we may choose an ordering on the set A , such that A can in fact be identified witha range of integers, A ≅ { , , , · · · , L/ } . In most abinitio systems L is infinite, since the Hilbert space isinfinite-dimensional. Similarly, α = ( a, +1) is identifiedwith even integers, and α = ( a, − 1) with odd integers,creating an ordering of the single-particle orbitals φ α ( x ),and α is identified with an integer 0 ≤ I ( α ) ≤ L .The single-particle matrix elements h ab and the two-particle elements u abcd are defined by h ab := h ϕ a | H | ϕ b i = Z ϕ a ( ~r ) H ϕ b ( ~r ) d d r, and u abcd := h ϕ a ϕ b | U ( ~r ) | ϕ c ϕ d i = λ Z ϕ a ( ~r ) ϕ b ( ~r ) 1 r ϕ c ( ~r ) ϕ d ( ~r ) d d r d d r , (5)respectively.The spatial orbitals ϕ a ( ~r ) are usually chosen as eigen-functions of H , so that h ab = δ a,b ǫ a .The basis functions for N -particle Hilbert space areSlater determinants | Φ α ,α , ··· ,α N i defined by | Φ α , ··· α N i := a † α a † α · · · a † α N |−i , where |−i is the zero-particle vacuum. In terms of single-particle orbitals, the spatial representation isΦ α , ··· ,α N ( x , · · · , x N ) = 1 √ N ! X p ∈ S N ( − ) | p | N Y i =1 φ α p ( i ) ( x i ) , where S N is the group of permutations of N symbols.The Slater determinants are anti-symmetric with respectto permutations of both x i and α i , so that the orbitalnumbers α i must all be distinct to give a nonzero func-tion. Each orbital is then occupied by at most one parti-cle. Moreover, for a given set { α i } Ni =1 of orbitals, one cancreate N ! distinct Slater determinants that are linearlydependent. In order to remove this ambiguity, we chooseonly orbital numbers such that I ( α i ) < I ( α j ) whenever i < j .It follows, that there is a natural one-to-one corre-spondence between Slater determinants with N parti-cles and integers b whose binary representations have N bits set. (If | A | = L < ∞ , the integers are limited to0 ≤ b < L .) Each bit position k corresponds to anorbital φ α ( x ) through k = I ( α ), and the bit is set ifthe orbital is occupied. Creating and destroying parti-cles in | Φ α , ··· ,α N i simply amounts to setting or clearingbits (possibly obtaining the zero-vector if a particle isdestroyed or created twice in the same orbital), keepingtrack of the possible sign change arising from bringingthe set { α i } on ordered form using Eqn. (4). Note thatthe vacuum |−i corresponds to b = 0, which is not thezero vector, but the single state with zero particles. − − − − R m n = 0 n = 0 n = 0 n = 0 n = 0 n = 0 n = 0 n = 0 n = 1 n = 1 n = 1 n = 1 n = 1 n = 2 ... n = 0 FIG. 1: Structure of single-particle orbitals of the two-dimensional harmonic oscillator. Angular momentum andshell number/energy on axes, and nodal quantum number n at each orbital. Orbital n = 0, m = − R = 3 isoccupied by two electrons for illustration. B. Model spaces The FCI calculations are done in a finite-dimensionalsubspace P of the N -particle Hilbert space, called themodel space. The model space has a basis B of Slater de-terminants, and P has the orthogonal projector P givenby P := X | Φ b i∈B | Φ b ih Φ b | . (6)The configuration-interaction method in general nowamounts to diagonalizing (in the sense of finding a fewof the lowest eigenvalues of) the, in general, large andsparse matrix P HP . The only approximation we havemade is the truncation of the N -particle Hilbert space.The model space P is seen to be a function of the singleparticle orbitals ϕ a ( ~r ), whom we choose to be the eigen-functions of H , i.e., harmonic oscillator eigenfunctions.These may be given on several equivalent forms, but it isconvenient to utilize rotational symmetry of H to createeigenfunctions of the projection of the angular momen-tum L z . In d = 2 dimensions we obtain the Fock-Darwinorbitals defined in polar coordinates by ϕ n,m ( r, θ ) = 1 √ π e imθ r | m | ˜ L | m | n ( r ) e − r / . (7)Here ˜ L kn ( x ) = ( − n [ n ! / ( n + | m | )!] / L kn ( x ) is the nor-malized generalized Laguerre polynomial. The factor( − n is for convenience, see Appendix A 1. The har-monic oscillator energy is 2 n + | m | + 1 and the eigenvalueof L z = − i∂/∂ θ is m . All eigenfunctions with the sameenergy 2 n + | m | + 1 =: R + 1 span a single-particle shell .The single-particle orbitals are illustrated in Fig. 1.For a Slater determinant | Φ α , ··· ,α N i , we have N X i =1 H ( i ) | Φ α , ··· ,α N i = E α , ··· ,α N | Φ α , ··· ,α N i with E α , ··· ,α N := N X i =1 ( R i + 1) , where R i = 2 n i + | m i | , and N X i =1 L z ( i ) | Φ α , ··· ,α N i = M | Φ α , ··· ,α N i , where M = P Ni =1 m i .To complete our definition of P , we let B = B R = ( | Φ α , ··· ,α N i : N X i =1 R i ≤ R ) , (8)where R is called the energy cut, for obvious reasons. As R → ∞ , the whole Hilbert space is spanned, and theeigenpairs of P HP converge to those of H . C. Configurational state functions and blockdiagonality In order to reduce the complexity of the computa-tions, we need to exploit symmetries of H . First of all,[ H, L z ] = 0, and it is obvious that also [ H, S z ] = 0, wherethe spin projection operator S z is given by S z := 12 X a,σ σa † a,σ a a,σ . The Slater determinants are eigenvectors of both L z and S z with eigenvalues M and s z = P Ni =1 σ i / 2, respectively.We obtain a natural splitting of the model space P intosubspaces with constant angular momentum M and spinprojection s z , viz, P = M M,s z P M,s z , P = X M X s z P M,s z . The diagonalization of H can thus be done within eachspace P M,s z separately, amounting to diagonalizing indi-vidual blocks P M,s z HP M,s z .The Hamiltonian (3) also commutes with total electronspin S , [ S , S z ] = 0, given by S := S z + 12 ( S + S − + S − S + ) , with S ± := X a a † a ± a a ∓ , so that a common basis for S z and S would lead to evensmaller matrix blocks.The eigenvalues of S are on the form s ( s + 1), where0 ≤ s ≤ N is an odd (even) integer for odd (even) N . For a joint eigenfunction of S z and S , called a con-figurational state function (CSF), | s z | ≤ s . The Slaterdeterminants are, however, not eigenfunctions of S , butsuch can be constructed by taking linear combinationsof a small number Slater determinants. For details onthis algorithm, see Ref. [9]. Suffice it to say here, that S only couples Slater determinants with identical sets ofdoubly occupied orbitals (meaning that φ ( a, +) and φ ( a, − ) are both occupied, as in Figure 1) and singly occupiedorbitals (meaning that only one of φ ( a, +) and φ ( a, − ) areoccupied). It is easy to see that S does not couple Slaterdeterminants in P M,s z to another P M ′ ,s ′ z . Thus, we ob-tain the splitting P M,s z = M s P M,s z ,s . We stress that all the mentioned operators commute witheach other, viz, [ H, Ω i ] = [Ω i , Ω j ] = 0 , with Ω i ∈ { L z , S z , S } . If a modified problem breaks,say, rotational symmetry, such that [ H, L z ] = 0, we maystill split the model space into to the eigenspaces of S z and S . D. Matrix elements of Coulomb interaction The remaining ingredient in the FCI method is theCoulomb matrix elements u abcd defined in Eqn. (5). Thesecan be calculated by first expanding L kn ( x ) in powers of x using L kn ( x ) ≡ n X m =0 ( − m ( n + k )!( n − m )!( k + m )! m ! x m , and evaluating the resulting integral term-by-term by an-alytical methods [10]. The resulting expression is a seven-fold nested sum, which can be quite time-consuming, es-pecially if a large number of Fock-Darwin orbitals occursin the basis B . Moreover, the terms are fractions of fac-torials with alternating signs, which is a potential sourceof loss of numerical precision.We therefore opt for a more indirect approach, giv-ing a procedure applicable to a wide range of potentials U ( r ) in addition to the Coulomb potential. Moreover,it can be generalized to arbitrary spatial dimensions d .The approach is based on directly transforming the prod-uct functions ϕ a ( ~r ) ϕ b ( ~r ) to the centre-of-mass system,where the interaction U ( r ) only acts on the relativecoordinate, and then transforming back to the lab sys-tem. This reduces the computational cost to a doublynested sum, as well as the pre-computation of the centre-of-mass transformation and the relative coordinate inter-action matrix. Both can be done exactly using Gaussianquadrature. The transformations to and from the cen-tre of mass frame are unitary transformations, which arestable and will not magnify round-off errors. In Appendix A we provide the details of the centre-of-mass transformation. One then obtains the follow-ing prescription for the interaction matrix elements u abcd :Let a = ( µ , ν ), b = ( µ , ν ), c = ( µ , ν ), and d =( µ , ν ) be the circular quantum number equivalents ofthe usual polar coordinate quantum numbers n i and m i .Due to conservation of angular momentum, we assume m + m = m + m ; otherwise, the matrix element u abcd = 0. Define M = µ + µ , M ′ = µ + µ , N = ν + ν ,and N ′ = ν + ν . Since u abcd is linear in λ , we set λ = 1without loss of generality. Now, u abcd = M X p = p T ( M ) p,µ T ( M ′ ) p ′ ,µ N X q = q T ( N ) q,ν T ( N ′ ) q ′ ,ν C | p − q | n,n + s , (9)where n = min( p, q ), s = M ′ − M , p ′ = p + M ′ − M , q ′ = q + N ′ − N . Moreover, p = max( M ′ − M, 0) and q = max( N ′ − N, T ( N ) are centre-of-mass transformation coeffi-cients defined in Appendix A, while the relative coor-dinate interaction matrix elements C | m | n,n ′ , n, n ′ ≥ 0, aredefined by C | m | n,n ′ := h ϕ n,m ( r, θ ) | U ( √ r ) | ϕ n ′ ,m ( r, θ ) i = 2 Z ∞ r | m | ˜ L | m | n ( r ) ˜ L | m | n ′ ( r ) U ( √ r ) e − r r d r. (10)Depending on U ( r ), the integral is best computed us-ing generalized half-range Hermite quadrature (see Ap-pendix B and Ref.[11]) or Gauss-Hermite quadrature.Weights and abscissa for quadratures are convenientlycomputed using the Golub-Welsch algorithm [12], whichonly depends on the ability to compute the coefficientsof the three-term recursion relation for the polynomialclass in question, as well as diagonalizing a symmetrictri-diagonal matrix.Let p ( r ) be a polynomial, and let α < β benon-negative constants. Then U ( r ) = r α p ( r ) e − βr (11)admit exact evaluations using generalized half-rangeGauss-Hermite quadrature. The Coulomb potential,Gaussian potentials, and the parabolic interaction − λr / α = 1 and p ( r ) = q ( r ) (i.e., an evenpolynomial), the integral is more convenient to eval-uate using standard Gauss-Hermite quadrature. TheCoulomb interaction falls into this class.Of course, one may let p ( r ) be a non-polynomial func-tion as well and still obtain very good results, as longas p ( r ) is well approximated with a polynomial, e.g., issmooth. III. EFFECTIVE INTERACTIONSA. Motivation The FCI calculations converge relatively slowly asfunction of the model space parameter R [3], as the er-ror ∆ E in the eigenvalue behaves like o ( R − k ) in general,where k = O (1). This behaviour comes from the singularnature of the Coulomb interaction.In Ref. [3], numerical results using an effective interac-tion were presented. This method is widely used in no-core shell model calculations in nuclear physics, wherethe nucleon-nucleon interaction is basically unknown buthighly singular [4]. This so-called sub-cluster effectiveinteraction scheme replaces the Coulomb interaction (oranother interaction) U ( r ij ) = λ/r ij with a renormalizedinteraction ˜ U ( i, j ) obtained by a unitary transformationof the two-body Hamiltonian that decouples the modelspace P and its complement [13]. Therefore, the two-body problem becomes exact in a finite number of har-monic oscillator shells. Loosely speaking, the effectiveinteraction incorporates information about the interac-tion’s action outside the model space. In general, ˜ U ( i, j )is non-linear in λ and not a local potential.Using the renormalized ˜ U ( i, j ), the many-body systemdoes not become exact, of course, but ˜ U ( i, j ) will per-form better than the bare interaction in this setting aswell. To the author’s knowledge, there exists no rigorousmathematical treatment with respect to this, but it hasnevertheless enjoyed great success in the nuclear physicscommunity [4, 14, 15], and our numerical experimentsunambiguously demonstrate that the convergence of theFCI method is indeed improved drastically [3], especiallyfor N ≤ U ( i, j ) is very small compared to the remaining calcula-tions. B. Unitary transformation of two-bodyHamiltonian We now describe the unitary transformation of thetwo-body Hamiltonian (i.e., Eqn. (1) or (3) with N = 2)that de-couples P and its complement. This approachdates back as far as 1929, when Van Vleck introducedsuch a generic unitary transformation to de-couple themodel space to first order in the interaction [16, 17].Let P be given by Eqn. (6), and let D = dim( P ). Theidea is to find a unitary transformation H = Z † HZ of H such that (1 − P ) H P = 0 , i.e., H is block diagonal. This implies that H eff definedby H eff := P H P has eigenvalues identical to D of those of the full operator H . Since D is finite, H eff is called an effective Hamilto-nian.Selecting Z is equivalent to selecting a set of effectiveeigenpairs { ( E k , | Ψ eff k i ) } Dk =1 , where E k is an eigenvalue of H and {| Ψ eff k i} Dk =1 ⊂ P are the effective eigenvectors; anorthonormal basis for P . It is clear that Z is not unique,since there are many ways to pick D eigenvalues of H ,and for each such selection any unitary D × D matrixwould yield an eigenvector set.However, some choices are more natural than others,since the eigenvectors and eigenvalues are usually con-tinuous functions of λ . We then select the D eigenvalues E k ( λ ) that develop adiabatically from λ = 0. For thecorresponding effective eigenvectors | Ψ eff k ( λ ) i , we choosethe orthonormal set that minimizes the distance to theexact eigenvectors {| Ψ k i} Dk =1 , i.e., {| Ψ eff k i} Dk =1 := argmin {| Ψ ′ k i} Dk =1 D X k =1 k| Ψ k i − | Ψ ′ k ik , (12)where the minimization is taken over orthonormal setsonly. The effective eigenvectors also turn out to be con-tinuous functions of λ , so H eff will also be continuous.Let U is the D × D matrix whose columns contain P | Ψ k i in the chosen basis, and let V be the correspond-ing matrix containing | Ψ eff k i . Clearly, V is unitary, while U only approximately so. Equation (12) can then bewritten V := argmin U ′ trace[( U − U ′ )( U − U ′ ) † ] , (13)where the minimum is taken over all unitary matrices. If U has singular value decomposition given by U = X Σ Y † , (14)the solution V is given by V := XY † . (15)If E = diag( E , · · · , E D ) is the diagonal matrix whoseelements are the chosen eigenvalues, we have H eff = V EV † . See Ref. [13] for a thorough discussion of the above pre-scription for H eff .Having computed the two-body H eff , we define the ef-fective interaction ˜ U (1 , 2) by˜ U (1 , 2) := H eff − P X i =1 H ( i ) P, which gives meaning solely in the model space. In secondquantization,˜ U (1 , 2) := 12 X abcd X στ ˜ u abcd a † aσ a † bτ a dτ a cσ , and the N -body H eff becomes (cf. Eqn. (1)) H eff = N X i =1 H ( i ) + N X i 2) wascomputed in a two-body energy cut space with parame-ter R , H eff is well-defined on the many-body model spacewith the same cut R . C. A comment concerning the choice of modelspace The two-body problem is classically integrable, i.e.,there exists 2 d − i , such thattheir quantum mechanical observables commute with H and each other, viz,[ H, Ω i ] = [Ω i , Ω j ] = 0 , for all i, j. Indeed, the centre-of-mass harmonic oscillator H C de-fined in Eqn. (18) below and the corresponding centre-of-mass angular momentum provides two constants, whiletotal angular momentum L z provides a third.Using the model space P defined by an energy cut, wehave [ P, Ω i ] = 0as well, which is equivalent [13] to[ H eff , Ω i ] = 0 , (17)so that H eff is integrable as well. In particular, ˜ U (1 , 2) isblock-diagonal with respect to Ω i .If we consider the commonly encountered model space P ′ defined by the Slater determinant basis B ′ given by B ′ := {| Φ α , ··· ,α N i : max( R i ) ≤ R } instead of Eqn. (8), we will have[ P ′ , H C ] = 0 , as is easily verified. Indeed, P ′ is not an invariant sub-space of the centre-of-mass transformation T defined inAppendix A. Thus, [ H ′ eff , H C ] = 0, so that the centre-of-mass energy no longer is a constant of motion! Thesymmetry-breaking of the effective Hamiltonian in this case is problematic, since in the limit λ → 0, the ex-act eigenfunctions that develop adiabatically are not alleither in the model space or in the complement. The adi-abatic continuation of the eigenpairs starting out in P isthus not well-defined.We comment, that the model space P ′ is often usedin both no-core shell model calculations and quantumdot calculations, but the effective interaction becomes,in fact, ill-behaved in this case. D. Solution of the two-body problem What remains for the effective interaction, is the com-putation of the exact eigenpairs { ( E k , | Ψ k i ) } Dk =1 . Wemust also solve the problem of following eigenpairs adia-batically from λ = 0.For the two-body Coulomb problem, analytical solu-tions are available only for very special values for λ [18].These are useless for our purpose, so we must use numer-ical methods.A direct application of the FCI method using Fock-Darwin orbitals with a large R ′ > R will converge slowly,and there is no device in the method for following eigen-values adiabatically. As the eigenvalues may cross, select-ing, e.g., the lowest eigenvalues will not work in general.For the two-body problem, the Pauli principle leads toa symmetric spatial wave function for the singlet s = 0spin state, and an anti-symmetric wave function for thetriplet s = 1 spin states. For the spatial part, we exploitthe integrability of the system as follows. Define centre-of-mass coordinates by ~R := 1 √ ~r + ~r )and ~r := 1 √ ~r − ~r ) . Using these coordinates, the two-body Hamiltonian be-comes H = H ( ~R ) + h H ( ~r ) + U ( √ r ; λ ) i =: H C + H rel (18)where r = √ r := √ k ~r k . We have introduced theparameter λ explicitly in the potential in this equation. H is clearly separable, and the centre-of-mass coordinateHamiltonian H C is a trivial harmonic oscillator, while therelative coordinate Hamiltonian can be written as H rel := − ∇ + 12 r + U ( √ r ; λ ) , where in polar coordinates ~r = ( r cos θ, r sin θ ) we have ∇ = 1 r ∂∂r r ∂∂r + ∂ ∂θ . Applying separation of variables again, the eigenfunc-tions of H rel can be written ψ n,m ( ~r ) := e imθ √ π u n,m ( r )where n is the nodal quantum number. u n,m ( r ) satisfies K | m | u n,m ( r ) = µ n,m u n,m ( r ) (19)where K | m | := − r ∂∂r r ∂∂r + m r + 12 r + U ( √ r ; λ ) . (20)Equation (19) is an eigenvalue problem in the Hilbertspace L ([0 , ∞ ) , r d r ), where the measure r d r is inducedby the polar coordinate transformation. Although it isnatural to try and solve the radial problem using Fock-Darwin orbitals, this will converge slowly. The solutionto this problem is to use a radial basis of generalized half-range Hermite functions [11]. In Appendix B this is laidout in some detail.Equation (19) is a one-dimensional equation, so therewill be no degeneracy in the eigenvalues µ m,n for fixed m . In particular, the eigenvalues as function of the inter-action strength λ will not cross, and will be continuousfunctions of λ . We thus have µ m,n < µ m,n +1 for all n ,where n is the nodal quantum number.At λ = 0 we regain the harmonic oscillator eigen-values 2 n + | m | + 1. Correspondingly, the eigenfunc-tions ψ m,n ( r, θ ) approaches the Fock-Darwin orbitals ϕ m,n ( r, θ ), i.e., the harmonic oscillator eigenfunctions.For the radial part,lim λ =0 u m,n ( r ) = g | m | n ( r ) := √ r | m | ˜ L | m | n ( r ) e − r / . Reintroducing spin, the full eigenfunctions Ψ =Ψ n ,m ,n ,m are on the formΨ( x , x ) = ϕ n ,m ( ~R ) e im θ √ π u n ,m ( r ) χ s,s z , where s = 0 for odd m , and s = 1 for even m , and | s z | ≤ s is an integer.Let R i = 2 n i + | m i | be the shell numbers for the centre-of-mass coordinate and relative coordinate, respectively.The eigenvalue E = E n ,m ,n ,m is E = R + 1 + µ n ,m , with limit E −→ λ → R + R + 2 , which is the harmonic oscillator eigenvalue.As the centre-of-mass coordinate transformation con-serves harmonic oscillator energy, at λ = 0, the eigen-functions that are in the model space are exactly those obeying R + R ≤ R . Turning on the interaction adi-abatically, the eigenpairs we must choose for the effec-tive Hamiltonian at a given λ are exactly those with R + R ≤ R .The model-space projection P Ψ needed in Eqns. (12)and (13) is now given by P Ψ( x , x ) = ϕ n ,m ( ~R ) e im θ √ π h ˜ P | m | R − R u n ,m ( r ) i χ s,s z , where˜ P | m | R := ¯ n X n =0 | g | m | n ih g | m | n | , ¯ n = (cid:22) R − | m | (cid:23) , (21)where ⌊ x ⌋ is the integer part of x . This operator thusprojects onto the ¯ n + 1 first radial basis functions withgiven | m | .Due to Eqn. (17), the unitary operator Z can be de-composed into its action on blocks defined by tuples of n , m and m [13]. The minimization (12) can thenbe applied on block-per-block basis as well. Each sub-problem is equivalent to the calculation of an effectiveHamiltonian K eff of the radial problem for a given m and ¯ n .To this end, let ¯ n and m = m be given. Let U be the(¯ n + 1) × (¯ n + 1) matrix whose elements are given by U n,k = h g | m | n | u m,n i , ≤ n, k ≤ ¯ n, i.e., the model space projections of the exact eigenvec-tors with the lowest eigenvalues. Let U = X Σ Y † be thesingular value decomposition, and let V = XY † . Then, K eff = V diag( E , · · · , E ¯ n ) V † and˜ C ¯ n, | m | := K eff − diag( | m | +1 , | m | +1 , · · · , n + | m | +1)is the ( n , m , m )-block of the effective interaction. Ifwe return to Eqn. (9), the effective interaction matrixelements ˜ u abcd are now given by replacing the matrix ele-ments C | p − q | n,n + s by the matrix elements ˜ C ¯ n, | p − q | n,n + s , where¯ n = (cid:22) R − R − | m | (cid:23) , R = N + M − ( p + q ) , where N , M , p and q are defined immediately afterEqn. (9). IV. CODE ORGANIZATION AND USEA. Overview The main program is called qdot , and processes atextual configuration file with problem parameters be-fore proceeding with the diagonalization of the Hamilto-nian. Eventually, it writes the resulting data to a Mat-lab/Gnu Octave compatible script for further process-ing.As a C++ library as well as stand-alone application, OpenFCI is organized in several namespaces, whichlogically separate independent units. There are threemain namespaces: manybody , gauss , and quantumdot .Put simply, manybody provides generic tools for many-body calculations, such as occupation number formal-ism, Slater determinants and CSFs, while gauss providestools for orthogonal polynomials and Gaussian quadra-ture. These namespaces are independent of each other,and are in no way dependent on the particular quantumdot model. On the other hand, quantumdot synthesizeselements from the two former into a quantum dot FCIlibrary. In qdot , the main work is thus processing ofthe configuration file.Two other namespaces are also defined, being simple sparse and simple dense , which are, respec-tively, simple implementations of sparse and dense ma-trices suitable for our needs. We will not go into detailsin the present article.It should be clear that extending and customizing qdot is a relatively easy task. The application qdot is provided as a tool with a minimum of functionality,and the interested will almost certainly desire to furtherdevelop this small application.In order to help with getting started on such tasks,some stand-alone demonstration applications are pro-vided, all based on the core classes and functions. Theseinclude an interaction matrix element tabulator tabu-late , and a simple program pairing for studying thewell-known pairing Hamiltonian [19], which we will notdiscuss further here. Finally, there is a small interac-tive console-based Slater determinant demonstration pro-gram slater demo as well. These applications will alsoserve as indicators of the flexibility of OpenFCI . OpenFCI does not yet support parallel computationon clusters of computers, using for example the Mes-sage Passing Interface [20]. Future versions will almostcertainly be parallelized, but the present version in factcompetes with parallel implementations of the standardFCI method with respect to convergence due to the ef-fective interaction implemented, see Sec. IV C. The sim-ple structure of OpenFCI also allows users with less re-sources to compile and run the code. B. Core functionality The manybody namespace currently contains four mainclasses: Slater , CsfMachine , NChooseKBitset , and MatrixMachine . These will probably form the backboneof any manybody computation with OpenFCI .The class Slater provides Slater determinants, cre-ation and annihilation operators, and so on. It is basedon the standard template library’s (STL) bitset class,which provides generic bit set manipulations. The class NChooseKBitset provides means for generating sets of k objects out of n possible represented as bit-patterns, i.e.,bit patterns corresponding to Slater determinants in the basis B or B ′ . This results in a STL vector CsfMachine is a tool for converting a basis ofSlater determinants into a basis of configurational statefunctions. These are represented as vector MatrixMachine , which is a template class, andgenerates a sparse matrix P AP of an operator A , where P projects onto the basis. It also handles bases of pureSlater determinants as they are trivially dealt within the CSF framework. The template parameter to MatrixMachine is a class that should provide the matrixelements h αβ , u αβγδ , etc, of the generic operator given by A = X αβ h αβ a † α a β + 12 X αβγδ u αβγδ a † α a † β a δ a γ + 13! X αβ ··· v αβγδǫζ a † α a † β a † γ a ζ a ǫ a δ . Notice, that the indices are generic orbitals, and not as-sumed to be on the form ( a, σ ) as in Eqn. (3).Currently, only one-, two-, and three-body operatorsare implemented. The reason is, that the matrix ele-ments are not computed by directly applying the sumof creation- and annihilation operators to Slater deter-minants, since this approach, however natural, is veryinefficient. Instead, we apply Wick’s theorem directly [8]on the matrix elements known to be not identically zero.In the gauss namespace, several functions are de-fined which computes sequences of orthogonal polyno-mials via recurrence relations and weights and abscissafor Gaussian quadratures based on these. The lat-ter is done using the Golub-Welsh algorithm, whichonly depends on being able to compute the coefficientsof the recurrence relation [12]. The most importantfunctions are perhaps computeLaguerrePolys() and computeGenHalfGaussHermite() , which computes a se-quence of generalized Laguerre polynomials evaluated ata given set of points and quadrature rules for generalizedhalf-range Hermite functions, respectively.Finally, the quantumdot namespace defines classes andfunctions that combined define the quantum dot prob-lem. The class RadialPotential encapsulates poten-tials on the form (11). It also computes effective inter-action blocks ˜ C ¯ n, | m | . The class QdotHilbertSpace pro-vides means for generating the bases B and B ′ , utilizingconservation of angular momentum, using a fast, custommade algorithm independent of NChooseKBitset . Theclass QdotFci sews everything together and is basicallya complete solver for the FCI method with effective in-teractions. FIG. 2: Simple configuration file for qdot TABLE I: Some ground state eigenvalues produced by qdot for N = 3 , λ = 2. Both the bare and theeffective interaction are used N = 3, M = 0, s = N = 4, M = s = 0 R E E , eff E E , eff . . . . . . . . . . . . . . . . . . . . . . . . C. Sample runs A basic configuration file for qdot is shown in Fig. 2.Varying the parameters lambda , R , S and the number ofparticles A [23], and running qdot each time, we pro-duce a table of ground state energies, shown in TableI. By changing the parameter use veff we turn on andoff the effective interaction. The corresponding effectiveinteraction ground states are also shown in the table.Notice, that with the effective interaction we obtain thesame precision as the bare interaction, but with muchsmaller model spaces. This indicates that OpenFCI canproduce results that in fact compete with parallel im-plementations of the standard FCI method, even in itspresent serial form.In Table II we compare the ground state energies re-ported in Ref. [9] with the alternative model space P ′ tothe corresponding values produced by qdot , also using P ′ . This Table also appears in Ref. [3], and serves as a check of the validity of the calculations.By uncommenting the lines following the definition of lambda , we override the default Coulomb interaction,and produce a configuration file for the analytically solv-able model given by Johnson and Payne [7], where theCoulomb interaction is replaced by the parabolic inter-action U ( r ) = − λr . If λ is sufficiently small, all the eigenvalues of this modelare on the form E j,k = 1 + j + ( k + N − √ − N λ, j, k ≥ . (22)Since the potential is smooth, the eigenfunctions are allsmooth, implying exponential convergence with respectto R [3]. We therefore expect very accurate eigenvalueseven with moderate R . In Table III we show the firsteigenvalues along with the error computed for N = 4electrons with λ = 1 / 8. The computations are done inthe M = s = s z = 0 model space with R = 10 and R = 15. Some duplicates exist, and they are includedfor illustration purposes. It is evident, that the eigen-values become very accurate with increasing R ; a clearindication of the correctness of the implementation. V. CONCLUSION AND OUTLOOK We have presented OpenFCI , an open source full con-figuration interaction implementation for quantum dotsand similar systems. OpenFCI also implements a renor-malized effective interaction widely used in nuclear no-core shell model calculations, and we demonstrated thatsuch interactions are indeed useful in the quantum dotcalculations as well. OpenFCI is easy to extend and adapt. Possible appli-cations are computations on systems with more generalsymmetry-breaking geometries and in d = 3 spatial di-mensions. Also, a generalization of the CSF part of thecode to handle isobaric spin would allow us to handlenuclear systems.There is one more symmetry of the Hamiltonian H thatcan be exploited, namely that of conservation of centre-of-mass motion, which would further reduce the blocksizes of the matrices. We exploited this symmetry for theeffective interaction, but it is a fact that it is a symmetryfor the full Hamiltonian as well. Using the energy cutmodel space P we may take care of this symmetry in away similar to the CSF treatment [21].As mentioned, we have not parallelized the code atthe time of writing, but it is not difficult to do so. Afuture version will almost certainly provide parallelizedexecutables, for example using the Message Passing In-terface [20].0 TABLE II: Comparison of current code and Ref. [9], taken from Ref. [3] R = 5 R = 6 R = 7 N λ M s Current Ref. 9 Current Ref. 9 Current Ref. 92 1 0 0 3 . . . . . . . . . . . . . . . . . . . . . . 046 11 . . . . . 055 11 . . . . . 691 23 . . . . . 870 23 . . . . . 15 21 . . 13 21 . . 134 0 5 29 . . 44 29 . . 31 29 . . R = 10 R = 15 E ∆ E E ∆ E . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − . · − APPENDIX A: CENTRE OF MASSTRANSFORMATION1. Cartesian coordinates In this appendix, we derive the centre-of-mass (COM)transformation utilized in Eqn. (9) for the interactionmatrix elements u abcd .The one-dimensional harmonic oscillator (HO) Hamil-tonian ( p x + x ) / φ n ( x ) = (2 n n ! √ π ) − / H n ( x ) e − x / = ( n !) − / A nx φ ( x ) , (A1)where A x := ( x − ip x ) / √ x -coordinate, and where φ ( x ) = π − / exp( − x / n + 1 / H in x and x is found to have eigenfunctions on theform Φ n ,n ( x , x ) := φ n ( x ) φ n ( x ) and eigenvalues n + n + 1. Define the raising operators A x i := ( x i − ip x i ) / √ 2, so thatΦ n ,n ( ~x ) := ( n ! n !) − / A n x A n x Φ , ( ~x ) , (A2)where the non-degenerate ground state is given byΦ , ( x , x ) = 1 √ π e − ( x + x ) / . (A3)Note that this can equally well describe two (distinguish-able and spinless) particles in one dimension.To this end, we introduce normalized COM frame co-ordinates by (cid:20) ξ ξ (cid:21) = 1 √ (cid:20) − (cid:21) (cid:20) x x (cid:21) =: F (cid:20) x x (cid:21) (A4)The matrix F is symmetric and orthogonal, i.e., F T F = F = 1, transforming a set of Cartesian coordinates intoanother. The operator H is invariant under this trans-formation, so the eigenfunctions have the same form withrespect to these coordinates, viz,Φ ′ n ,n ( ξ , ξ ) := φ n ( ξ ) φ n ( ξ )= ( n ! n !) − / A n ξ A n ξ Φ ′ , , (A5)where A ξ i = ( ξ i − ip ′ i ) / √ p ′ i are the corre-sponding momentum components.Define the operator T by T ψ ( x , x ) := ψ ( ξ , ξ ) = ψ (cid:18) x + x √ , x − x √ (cid:19) , (A6)so that T Φ n ,n := Φ ′ n ,n . (A7)1Since T maps eigenfunctions in the two frames onto eachother, T must be a unitary operator, and the invarianceof H under the coordinate transformation is the same as[ H , T ] = 0, i.e., that energy is conserved. This in turnmeans that T is block diagonal with respect to each shell R = n + n , viz, T Φ R − n ,n = R X n =0 h Φ R − n,n | Φ ′ R − n ,n i Φ R − n,n =: R X n =0 T ( R ) n ,n Φ R − n,n , (A8)where T ( R ) is the ( R + 1) × ( R + 1) transformation ma-trix within shell R . It is real, symmetric, and orthogo-nal. Numerically, the matrix elements are convenientlycomputed using two-dimensional Gauss-Hermite quadra-ture of sufficiently high order, producing exact matrixelements.In a two-dimensional setting, the two-particle har-monic oscillator becomes a 4-dimensional oscillator. Let ~r i = ( x i , y i ), i = 1 , 2, be the particles’ coordinates, andlet a = ( m , n ) and b = ( m , n ) to compress the nota-tion a little. An eigenfunction is now on the formΨ a,b ( ~r , ~r ) := Φ a ( ~r )Φ b ( ~r )= CA m x A n y A m x A n y Ψ , ( ~r , ~r ) (A9)where C = ( m ! n ! m ! n !) − / .The COM coordinate transformation now acts in the x and y directions separately, viz, F acts on x i and y i to yield the COM coordinates ξ i and η i : [ ξ , ξ ] T = F [ x , x ] T and [ η , η ] T = F [ y , y ] T . The induced oper-ator T again conserves energy. Let M = m + m and N = n + n . It is readily verifiable that the COM frametransformation becomes T Ψ M − m ,N − n ,m ,n := Ψ ′ M − m ,N − n ,m ,n = M X p =0 T ( M ) m ,p N X q =0 T ( N ) n ,p Ψ M − p,N − q,p,q . (A10)Note that the shell number is R = N + M , which isconserved by T . 2. Centre of mass transformation for Fock-Darwinorbitals Consider a Fock-Darwin orbital ϕ n,m ( ~r ) in shell R =2 n + | m | with energy R + 1. It is straightforward butsomewhat tedious to show that these can be written interms of co-called circular raising operators B + and B − [22] defined by (cid:20) B + B − (cid:21) := 1 √ (cid:20) i − i (cid:21) (cid:20) A x A y (cid:21) . (A11) Letting µ = n + max(0 , m ) and ν = n + max(0 , − m )(which gives µ, ν ≥ 0) one obtains ϕ n,m ( ~r ) = ( µ ! ν !) − / B µ + B ν − Φ , ( ~r ) , (A12)which should be compared with Eqn. (A2). Moreover, R = µ + ν and m = µ − ν , giving energy and angularmomentum, respectively. We comment that this is thereason for the non-standard factor ( − n in the normal-ization of the Fock-Darwin orbitals in Eqn. (7).Let a two-particle HO state be given by˜Ψ µ ,ν ,µ ,ν := ϕ n ,m ( ~r ) ϕ n ,m ( ~r ) , = CB µ B ν − B µ B ν − Ψ , ( ~r , ~r ) , (A13)where µ i = n i + max(0 , m i ) and ν i = n i + max(0 , − m i ),and where C = ( µ ! ν ! µ ! ν !) − / . We will now provethat, in fact, when applying the centre-of-mass transfor-mation to Eqn. (A13), we obtain an expression on thesame form as Eqn. (A10) viz, T ˜Ψ M − µ ,N − ν ,µ ,ν := ˜Ψ ′ M − µ ,N − ν ,µ ,ν = M X p =0 T ( M ) µ ,p N X q =0 T ( N ) ν ,q ˜Ψ M − p,N − q,p,q , (A14)where M := µ + µ and N := ν + ν .To this end, we return to the raising operators A ξ i and A η i , and express them in terms of A x i and A y i . By usingEqn. (A4), we obtain (cid:20) A ξ A ξ (cid:21) = F (cid:20) A x A x (cid:21) , (A15)and similarly for A η i in terms of A y i . In terms of theraising operators, the COM transformation becomesΨ ′ m ,n ,m ,n = C ( A x + A x ) m ( A y + A y ) n × ( A x − A x ) m ( A y − A y ) n Ψ , , (A16)where C = (2 n + n + m + m n ! n ! m ! m !) − / , (A17)and we have used Eqn. (A9), but in the analogous COMcase. Expanding the powers using the binomial formula(and the fact that the raising operators commute), weobtain a linear combination of the individual eigenfunc-tions, which must be identical to Eqn. (A10).Let the COM circular ladder operators be defined by (cid:20) B ′ j + B ′ j − (cid:21) := 1 √ (cid:20) i − i (cid:21) (cid:20) A ξ j A η j (cid:21) . (A18)Using Eqn. (A15), we obtain that the circular raisingoperators transform in the same way as the Cartesianoperators when going to the COM frame, i.e., (cid:20) B ′ B ′ (cid:21) = F (cid:20) B B (cid:21) , (A19)2and similarly for B ′ i − and terms of B i − . UsingEqn. (A13) in the COM case, we obtain˜Ψ ′ µ ,ν ,µ ,ν = C ( B + B ) µ ( B − + B − ) ν × ( B − B ) µ ( B − − B − ) ν Ψ , , (A20)with C = (2 µ + ν + µ + ν µ ! ν ! µ ! ν !) − / . (A21)Eqn. (A20) is on the same form as Eqn. (A16). Again,by expanding the powers using the binomial formula (andthat the raising operators commute), we obtain a linearcombination with coefficients identical to those of the ex-pansion of Eqn. (A16). It then follows that Eqn. (A14)holds. APPENDIX B: NUMERICAL TREATMENT OFRADIAL PROBLEM We now briefly discuss the numerical method used forsolving the radial problem (19), i.e., the eigenvalue prob-lem for the operator K | m | defined in Eqn. (20). This is aneigenproblem in the Hilbert space L ([0 , ∞ ) , r d r ), wherethe measure r d r is induced by the polar coordinate trans-formation. The inner product on this space is thus givenby h f | g i = Z ∞ f ( r ) g ( r ) r d r. (B1)Let the Fock-Darwin orbitals be given by ϕ n,m ( r, θ ) = e imθ √ π g | m | n ( r ) , (B2)with radial part g | m | n ( r ) := √ L | m | n ( r ) r | m | exp( − r / 2) (B3)Thus, h g | m | n | g | m | n ′ i = δ n,n ′ , and these functions form an orthonormal sequence in L ([0 , ∞ ) , r d r ) for fixed | m | .In the electronic case, U ( √ r ) = λ/ √ r has a singular-ity at r = 0 which gives rise to a cusp in the eigenfunction u m,n ( r ) at r = 0, or in one of its derivatives. Away from r = 0, the eigenfunction is smooth. These considera-tions are also true for more general potentials smooth for r > K with respect to thetruncated basis { g | m | n ( r ) } ¯ nn =0 will give eigenpairs converg-ing slowly with respect to increasing ¯ n due to the non-smoothness of u n,m ( r ) at r = 0. This is easily seen forthe m = 0 case and λ = 1, which has the exact groundstate u , ( r ) = (cid:18) r + 1 √ (cid:19) e − r / , a polynomial of odd degree multiplied by a Gaussian.The cusp at r = 0 is evident. However, g n ( r ) = √ L n ( r ) e − r / , which are all even polynomials. It is clear, that ¯ n mustbe large to resolve the cusp of u , ( r ).The eigenproblem is best solved using a basis of gener-alized half-range Hermite functions f j ( r ) [11], which willresolve the cusp nicely. These functions are defined by f j ( r ) := P j ( r ) exp( − r / , where P j ( r ) are the orthonormal polynomials definedby Gram-Schmidt orthogonalization of the monomials r k with respect to the weight function r exp( − r ). Thus, h f j | f j ′ i = Z ∞ f j ( r ) f j ′ ( r ) r d r = δ j,j ′ . The fundamental difference between f j ( r ) and g | m | n ( r )is that the latter contains only even (odd) powers of r for even (odd) | m | . Both sets constitute orthonormalbases, but f j ( r ) will in general have better approximationproperties.Moreover, since deg( r | m | L | m | n ( r )) = 2 n + | m | , g | m | n ( r ) = n + | m | X j =0 h f j | g | m | n i f j ( r ) (B4)gives the Fock-Darwin orbitals as a finite linear combi-nation of the generalized half-range Hermite functions,while the converse is not possible.Computing the matrix of K with respect to { f j ( r ) } ¯ jj =0 and diagonalizing will give eigenpairs converging expo-nentially fast with respect to increasing ¯ j . The resultingeigenfunctions’ expansion in g | m | n are readily computedusing Eqn. (B4), whose coefficients h f j | g | m | n i can be com-puted numerically exactly using Gaussian quadrature in-duced by P J ( r ), for J sufficiently large.The basis size ¯ j to use in the diagonalization dependson how many eigenfunctions ¯ n we desire. We adjust ¯ j semi-empirically, noting that 2¯ n + | m | is sufficient to re-solve g | m | ¯ n , and assuming that the exact eigenfunctionsare dominated by the latter. We then add a fixed num-ber j to get ¯ j = 2¯ n + | m | + j , and numerical experimentsconfirm that this produces eigenvalues that indeed haveconverged within desired precision. ACKNOWLEDGMENTS