Finite-size scaling for quantum criticality using the finite-element method

Finite-size scaling for quantum criticality using the finite-element method

Edwin Antillon Departments of Physics, Purdue University, West Lafayette, IN 47907    Birgit Wehefritz-Kaufmann Departments of Physics, Purdue University, West Lafayette, IN 47907 Departments of Mathematics, Purdue University, West Lafayette, IN 47907    Sabre Departments of Physics, Purdue University, West Lafayette, IN 47907 Department of Chemistry, Purdue University, West Lafayette, IN 47907

Finite size scaling for the Schrödinger equation is a systematic approach to calculate the quantum critical parameters for a given Hamiltonian. This approach has been shown to give very accurate results for critical parameters by using a systematic expansion with global basis-type functions. Recently, the finite element method was shown to be a powerful numerical method for ab initio electronic structure calculations with a variable real-space resolution. In this work, we demonstrate how to obtain quantum critical parameters by combining the finite element method (FEM) with finite size scaling (FSS) using different ab initio approximations and exact formulations. The critical parameters could be atomic nuclear charges, internuclear distances, electron density, disorder, lattice structure, and external fields for stability of atomic, molecular systems and quantum phase transitions of extended systems. To illustrate the effectiveness of this approach we provide detailed calculations of applying FEM to approximate solutions for the two-electron atom with varying nuclear charge; these include Hartree-Fock, density functional theory under the local density approximation, and an “exact” formulation using FEM. We then use the FSS approach to determine its critical nuclear charge for stability; here, the size of the system is related to the number of elements used in the calculations. Results prove to be in good agreement with previous Slater-basis set calculations and demonstrate that it is possible to combine finite size scaling with the finite-element method by using ab initio calculations to obtain quantum critical parameters. The combined approach provides a promising first-principles approach to describe quantum phase transitions for materials and extended systems.

PACS numbers: 02.70.Dh 05.70.Jk 82.60.-s

I Introduction

The theory of finite size scaling (FSS) in statistical mechanics can provide us with numerical methods fisher (); widom (); barber (); privman (); cardy (); nightingale1 (); Peter1 (); Peter2 () capable of obtaining accurate results for infinite systems by studying the corresponding small systems. In the past, basis set approximation to many-body systems have been combined with FSS to obtain critical parameters, where calculations for these systems were done by expanding the wave function in Slater-type basis sets or Gaussian-type basis functions. However, a generalization of larger atomic and molecular systems proved itself difficult using these type of functions MSK (); PSK ().

More recently, we have combined FSS with the finite element method (FEM) in order to obtain critical parameters in a potential that admits a finite number of bound states and we obtained excellent results by using both finite-element and Slater-type basis functions AMWK (). For this potential we also established that the finite-difference (FD) method can be combined with FSS. The finite difference method, however, is not a variational method since it approximates the operators (i.e. the derivative) and therefore can result in errors of either sign TT2 (). As a result, we observed non monotonic behavior of the calculated critical parameters with respect to the number (N) of grid points used Daniel (). Nevertheless, for a large enough number of grid points N one can recover the correct asymptotic behavior in the large N limit. For this particular example we used FD to fourth order and in FEM we used -continuous basis.

In this paper, we present a method to combine FEM with FSS for electronic-structure calculations near regions of criticality. We note that any method where the wavefunction can be expanded in a complete basis (or set of grid points) should also work, provided that in the infinite basis limit () it becomes arbitrarily close to the the exact wavefunction. The main advantage of finite element methods, over other global basis methods, resides in the fact that matrices obtained by these numerical methods are often banded and tend to be sparse since variables are not coupled over arbitrarily large distances. This is ideal for describing extended systems since it facilitates parallelization pask (). Moreover, due to the polynomial nature of the basis set most integrals involved can be evaluated easily and accurately with a better convergence with respect to the basis set size.

We examine the two electron system with different levels of approximations to the ground state energy, by solving the system using FEM both exactly and using mean-field equations and demonstrate that this formalism is general and that it has potential for more complex systems by solving Hartree-Fock equations or Kohn-Sham equations in density-functional theory. The paper is organized as follows: In Sec. (II) we discuss FEM formalism, in Sec. (III) we discuss the FSS approach to quantum systems, in Sec. (III) and (IV) we combine FSS with FEM to compute calculations of critical parameters for two electron atoms, ,and, lastly in Sec. (V) we discuss the results and applications to extended systems.

Ii Finite Element Method in Electronic Structure Calculations

Computational methods in quantum chemistry typically make use of Slater or Gaussian basis sets to approximate solutions in atomic or molecular systems. The advantages of these methods are that some of the integrals involved in the calculations are analytical and very accurate results can be obtained for these integrals. However, there are problems in extrapolating to the complete basis limit. Plane wave methods, having a complete basis, present some difficulties if non-periodic boundary conditions are to be used such as in molecule calculations. A common trait of all these methods is that they employ global basis, which can result in some difficulties for the formulation of large problems. First, the matrix formulation becomes dense using global basis and therefore ill-conditioned for parallelism pask (). While for some complex objects, convergence issues may arise if there are large variations of potential terms in small regions which cannot be captured by the traditional Gaussian basis set that are centered on the nuclei.

The FEM presents an alternative solution to these problems. This method uses a variational method to solve differential equations by discretizing a continuous solution to a set of values in sub-domains, called elements, and employs a local basis expansion in terms of polynomials with variable real-space resolution, allowing for convergence to be controlled systematically. The locality of the bases in real space results in sparse and banded matrices, where the number of operations for matrix-vector multiplication can be reduced from for dense matrices to for sparse matrices, where N is the dimension of the matrix pask (). Some of the first applications of FEM in engineering date back to 1950, in problems such as elasticity and structural analysis in civil engineering and aeronautical engineering, respectively. In the 1970s, some applications to quantum mechanical problems appeared, but were limited to small systems due to the storage limitation FEM-Book1 (); FEM-Book2 () Recently, advances in high-performance computing make FEM a viable alternative to the traditional approaches of electronic-structure calculations FEM-Review1 (); FEM-Review2 (). The sparsity of the global matrices resulting from the FEM makes the parallel formulation more affordable. Alizadean et al, have recently applied a divide-and-conquer method to the FEM-Hartree-Fock approach for electronic calculations showing a facilitation on parallelization and reduced scaling for larger systems AHM ().

We briefly describe the FEM formulation for a given Hamiltonian to be used in later sections. Let us consider the following one-particle Hamiltonian:


where and are the energy eigenvalue and eigenvector respectively, and V is an arbitrary potential that will represent a sum of terms such as the electron-nuclear potential (), Hartree potential () and exchange-correlation potentials (). We can reformulate this problem into a weak form by taking an inner product of an arbitrary “test function” with Eq (1):


This relaxes the problem into a variational form where instead of finding an exact solution in the domain , we find a solution that satisfies Eq.(1) in an average sense everywhere. We then arrive at an equivalent equation by integrating by parts (in real space):


where is a boundary term and is the outward unit normal at r. In the absence of external fluxes we can employ Neumann boundary conditions []. Furthermore, in radial coordinates () we can write the surface term in the following form:


This surface term vanishes trivially at , while for , it can be removed if we choose a test function that vanishes at . Since Eq (3) is satisfied for any , we can choose and in this manner it reduces the problem to an eigenvalue problem. The wavefunction can be then discretized using the Garlekin approximation FEM-Book1 (); FEM-Book2 ():


where are piecewise-linear continuous functions (i.e. continuous but not necessarily smooth) and are the nodal degrees of freedom. Smoother bases can be constructed by using higher order polynomial ( continuous, ) and by extending the nodal space to include the derivative of the wavefunction at the node, viz, ; this ensues continuity of the wavefunction between the nodes. In the nodal basis, the generalized eigenvalue problem reads


where the above terms can be written in terms of local basis representation, viz

The finite element method takes advantage of the strict locality of its basis to yield matrices that are sparse and structured and due to the polynomial nature of the basis most integrals can be easily evaluated, as in global basis approaches. Moreover, since the method is variational, all errors are of the same sign, leading to monotonic convergence often from below pask2 ().

Iii Finite Size Scaling

We now discuss the role of finite size scaling in systems exhibiting critical behavior. Phase transitions are associated with singularities of the free energy that occur only in the thermodynamic limit yanglee1 (); yanglee2 (). Near criticality, fluctuations in the system become correlated over distances of the order of the correlation length , while for second order phase transitions, this correlation length diverges.

The theory of finite size scaling can be used to extract universal parameters describing the transition by studying systems of finite size widom (); barber (); privman (); cardy (); nightingale1 (); Peter1 (); Peter2 (). Any physical quantity () that would diverge in the thermodynamic limit, becomes now “rounded off” in the finite systems (of size L) by a scaling function cardy (); nightingale1 ()


where is a reduced critical parameter describing the phase transition and is the correlation length of the system in the thermodynamic limit.

While classical phase transitions are governed by the fluctuation of temperature that tend to minimize the free energy, at zero temperature, these fluctuations are driven by Heisenberg’s uncertainty principle sachdev (). A signature of a phase transition is non-analyticity of the ground state energy with respect to a parameter in the Hamiltonian. The non-analyticity can be the result of avoided crossing, or actual level crossing. In the first case, the configuration of the system that minimizes the energy maintains the energy gap () between the first excited energy state and the ground state as the parameter is varied. In the latter case, the energy levels cross so that the energy is now minimized by a configuration that previously corresponded to the first excited state of the system.

To pin down the transition point involves computing the energy accurately near the critical point. However, this can be quite a challenge since most real systems lack an analytical solution to the ground state energy and approximation methods, such as perturbation theory or mean-field methods (Hartree-Fock , DFT) which are viable options in stable systems, in general are not well defined or become ill-conditioned near critical points.

This is where the theory of finite size scaling can be recast in the context of the Hilbert space. Any given approximation to the ground state energy will result in a “rounding” effect from the “true” critical behavior, since any approximation to a wavefunction can be expanded in some basis and truncated to some order N; that is,


FSS has been used in quantum systems neirotti0 (); serra2 (); kais (); snk1 (); snk2 (); nsk (); qicun (); kais1 (); review (); adv (). In this approach, the finite size corresponds not to the spatial dimension but to the number of elements (N) in a complete basis set used to expand the exact eigenfunction of a given Hamiltonian dipole (); quadrupole (); qicun1 (); qicun2 (); Ferron-Last (). In the variational approach singularities in different mean values will occur only in the limit of infinite number basis function, viz.,


The FSS ansatz in the Hilbert space takes the following form (for expectation values of operators acting on the Hilbert space). There is a scaling function such that


where is a critical parameter in the Hamiltonian, is a quantum mechanical operator near criticality, and is a scaling function for each different operator but with a unique scaling exponent. At the critical point, the expectation value is related to as a power law, . In this manner the scaling function will remove the divergence as it is computed in the finite basis.

The asymptotic behavior of operators near the critical point, in this case the ground state energy, can be expressed using critical exponents similarly to classical phase transitions, viz


where is a universal critical exponent. In addition to the vanishing energy scale, a diverging length scale also occurs which is associated with the exponential decay of equal-time correlations in the ground state of the system sachdev (). For our case, we focus on the behavior of the wavefunction at threshold, where the correlation length is defined by:


In order to apply finite size scaling method to quantum system, we assume a Hamiltonian of the form


where is the coupling parameter showing a phase transition at . For potentials linear in ( ), Simon simon () showed that the critical exponent is equal to one if and only if the Hamiltonian has a normalizable eigenfunction with eigenvalue equal to zero. The existence or absence of a bound state at the critical point is related to the type of singularity in the energy. Using statistical mechanics terminology, we can associate a“first order phase transition” with the existence of a normalizable eigenfunction at the critical point. The absence of such a function could be related to a “continuous phase transition.” To obtain the critical parameters, we define the following function that will remove the N-dependence at a critical point:


For the energy operator , and using the customary Greek letter for the corresponding exponent we have


To arrive at an expression that isolates the critical exponent from numerical calculations, we make use of the Hellmann-Feynman theorem landau ()


At this point it becomes convenient to define another function that will be used in the following sections.


At the critical point this function is independent of and and takes the value of . Namely, for and any values of and we have


Iv Implementation

In the following, we illustrate the combined approach (FEM+FSS) for the two-electron atom, but emphasize that the method is general and can be applied to more complex atomic and molecular systems. We seek the minimum charge (critical charge) needed to bind both electrons to the nucleus with charge Z. It is well known that ignoring the electron-electron energy correlation will not predict the stability of chandr (), since the critical charge obtained by such assumption is . Stillinger (1974) Stillinger () took into account the e-e correlation by using a non-linear wavefunction in order to find the minimum energy as the nuclear charge Z was varied. He found a critical charge around . Baker et al (1990) Baker () performed a thorough 401-order perturbation calculation to resolve the controversy for the critical charge; they obtained a critical charge of . This is the value we use as reference. In the following sections we illustrate the mean-field equations and exact formulations for the ground state of the two-electron system, while in the next section we apply FSS to these approximations. Let us start with the Hamiltonian for the two-electron atoms; it is given by


where and are the coordinates of the two electrons referenced from the nuclei. We use Hartree units () throughout in our calculations.

Hartree Fock (HF) Approximation:

In the Hartree approximation the total wavefunction is approximated as a product of single particle wavefunctions is . We can arrive at a self-consistent form by multiplying the Hamiltonian by the complex conjugate of the aforementioned wavefunction and integrating one coordinate out, we arrive at


Using the FEM (see Sec. II) this expression can be reduced to a generalized eigenvalue problem . The Hamiltonian is written in terms of the single electron energy term (i.e. ) and a direct coulomb term defined below, viz


Since electrons are fermions, the product of the single particle states needs to be antisymmetric with respect to exchange of the coordinates on the electrons. This can be achieved by writing the wave function as a Slater determinant thijssen ():


Searching for the best Slater determinant self-consistently in the Hamiltonian [Eq.19] leads to the Hatree-Fock equations thijssen (); Schw (); GG ().


where the additional term on the right is known as the exchange Coulomb term, denoted by , and tends to repel states of the same (parallel) spin. We can write this term in the local basis form by looking at the direct Coulomb term and noting that the indices and similarly indices . The operator is obtained from by switching indices (), which effectively switches the coordinates and on two of the wavefunctions and makes this operator non-local.

The ground state energy (for electrons) is determined by solving for the ground state eigenfunction for the eigenvalue problem (Eq.23) and computing the expectation value below (to avoid double-counting terms),


Note that for our case, the exchange term does not arise since the ground state of the helium atom does not have parallel spins. In order to improve upon the Hartree-Fock approximation we can take into account the correlation energy from an approximation to the free electron gas due to Wigner. Here the “correlation” potential can be obtained from Refs. wigner1 (); wigner2 (). Note that the references cited use Rydbergs units, and reference wigner1 () additionally uses . In Hartrees the correlation potential is given by:


where is the local density parameter for the atom, that is


The last expression follows for helium since in the ground state the two electrons have opposite spin, but share the same spatial orbital. In this manner the correlation potential is included self-consistently and the density is obtained from the sum of the densities for each electron i.e., . Finally, the correlation energy is obtained as

and the value for the total energy is given by:


Density Functional Theory

In the Hartree-Fock approximation, the many-body wavefunction is written in the form of a Slater determinant while the correlation energy is treated separately. This can be done for small systems, but for larger systems this is not an easy task. The theory of Kohn and Sham developed in 1965 simplifies both the exchange and correlation effect of electrons by reducing their effect to a effective potential that is in principle exact. In their seminal paper KS (), they showed that to determine the ground state energy it is sufficient to determine the ground state electron density . By constructing a simple system of non-interacting electrons which yields the same density as the “real” fully interacting system, the ground state energy can be obtained.

The Kohn-Sham equations for two electrons are:


These equations are self-consistent in the density and differ from the Hartree-Fock equations only in the exchange-correlation potential which is a function of the total local density:


While the exact form of the exchange-correlation functional is not known, the simplest approximation is the local density approximation (LDA)


It follows from Eq.(30) that the potentials are given by:


and the ground state energy of the system is obtained from thijssen ():

LDA Total Energy
FEM FEM Perdue-Wang FEM FEM Exact
(N=100) (N=200) UG (); PW () (N=100) (N=200) UG ()
-2.813 661 -2.821 852 -2.834 455 -2.886 297 -2.904 925 -2.903 724
2.703 352 2.726 064 2.767 389 2.806 848 2.855 856 2.867 082
-6.538 681 -6.571 185 -6.624 884 -6.664 319 -6.737 783 -6.753 267
1.973 543 1.978 456 1.995 861 1.020 329 1.026 754 1.024 568
-0.850 814 -0.854 086 -0.861 740 0.000 000 0.000 000 0.000 000
-0.101 061 -0.101 092 -0.111 080 -0.049 135 -0.049 814 -0.042 107
-0.563 621 -0.565 805 -0.570 256 -0.957 555 -0.963 991 -0.903 724
Table 1: Comparison of various components to the Kohn-Sham and Hartree-Fock energies in Hartrees for helium. The reference values (Perdue-Wang, Exact) are correct to all digits quoted UG (); PW (), while the accuracy of the FEM estimates on is expected to be on the order of two digits (0.01) for N=100 and three digits (0.003) for N=200 for both LDA and Total Energy, refer to the discussion on the text for details.

The Kohn-Sham equations are in general simpler than HF equations, while including correlation effects beyond the HF approximation. Furthermore their computing time grows as a power of the number of the N-electrons, whereas the exact solution for the N-electron Schrodinger equations is exponential in N GG ().

Let us comment on the difference between the two approaches, even though both Hartree-Fock and Kohn-Sham equations have similar potential terms: [see Eqs.(23) and (29)]. In HF the energy is calculated using the wavefunction explicitly, whereas in DFT the energy is a functional of the density. For example, for the correlation energy we have in the HF approximation:

while for the DFT approximation we have:

Care should be exercised when computing the LDA exchange potential since this term is non-linear in the density and as a result the normalization does not trivially cancel out when solving the eigenvalue problem. For example, in this work we use the normalization so that a factor of is absorbed into the wavefunction . This in turn changes the density to , and therefore the exchange-potential in LDA becomes

In order to satisfy the boundary conditions as , we use a uniform mesh with a cutoff radius at . The radial equations where solved using -continuous polynomials, that is, and ; where the function interpolates the wavefunction from the left node to the right node and similarly for . We solved all the integrals numerically using 10-point Gaussian quadratures and in order to solve the generalized eigenvalue problem we have used a fortran subroutine EWEVGE EWEVGE ().

In Table 1 we show the energy components of the Kohn-Sham equations under LDA approximation and the Hartree-Fock equations using the Wigner correlation potential for the helium atom. We note that the reference values for the Total Energy (HF+Wigner correlation) are exact, therefore even in the infinite basis limit the expression will not converge to the value shown. We also point out that the Perdue-Wang values on the correlation energy are obtained by using the random-phase approximation (RPA) for the free electron gas PW (). Our convergence value is be limited by how well we can approximate this value. We have checked our level of convergence by increasing the number of elements to N=1000 in a linear mesh of size with -continuous basis (linear polynomials). We obtain a total energy of and correlation energy of , and when compared to the values by Perdue-Wang (see Table 1 ), the errors are and respectively. We see that the error in the total energy is dominated by the correlation term arising from using Wigner correlation rather than RPA description. However, taking this difference into account we can get an estimate for the error of all other terms, viz . In FEM, the error on integrals such as the energy is , where h is the element width and p is the order of the polynomial used. Thus, for our example using N=1000 in a uniform grid size with a.u. we have , and since we employ a linear basis, the error we expect is . In this manner, if we take into account the intrinsic error due to the different energy-correlation terms we can recover the level of convergence expected by using FEM.

Exact Formulation

Figure 1: example of an adaptive mesh GMSH () used in approximating the solution in three dimensions. In this case the number of elements is , while the size of the element progressively increases by a factor of 1.3 as is increased

We can compare the one-particle formalism above with an exact formulation of the system. In 1930 Breit showed that the Schrödinger equation for the two-electron atom becomes separable if one of the coordinate axes is aligned with one of the radius vector Breit (). For the s state, the exact wavefunction depending on six total variables can be reduced to only three variables: the positions of the two electrons and , and the angle between them , viz



We approximate the solution to the above Hamiltonian using FEM by discretizing each variable above independently as in Levin (). The boundary conditions implemented make the wavefunction vanish at a cutoff radius of a.u, and the angular range is between . Furthermore an adaptive grid was refined closer to the nuclei as shown in Fig.(1) to speed up calculations. In our implementation we use -continuous basis, that is Hermite polynomials and to interpolate the left and right nodal values of the wavefunction and additionally and to interpolate the derivatives of the wavefunction, where is the element width and . In this manner, the continuity of the wavefunction is enforced by using this basis. All the integrals involved in constructing the matrices were solved numerically using 10-point Gaussian quadratures, and in order to solve the generalized eigenvalue we used a sparse matrix solver in ARPACK++ ARPACK ().

For the mesh shown in Fig. 1, we obtain a ground state energy of using -continuous basis and using -continuous basis. The latter value is close to the energy obtained by Levin using the same approach Levin (), which is very close given that we perform all integrals numerically.

V Results and Discussion

In order to find the critical parameter we approximate the energy gap in the ground state numerically as

where is the binding energy obtained for two electrons and is binding energy for a single electron. The theory of finite size scaling is applied to a helium-like Hamiltonian of the form,


where and are the coordinates of the two electrons referenced from the nuclei, and Z is the charge of the nucleus, this will play the role of the critical parameter. For reference, we note that the transformation takes Eq. (38) into:


this form is often used in the literature to obtain the critical parameter (see Stillinger (); Baker ()). Below we outline a general procedure that can be used to combine FSS with numerical methods that expand a given wavefunction in the Hilbert space in any complete basis at a truncated level N.

General Procedure:

  1. Compute the expectation value of and for various values of N (number of elements). In this case the Hellmann-Feynman theorem (see Eq. 38) becomes useful since in this case is linear in and the derivative is quite simple .

  2. Construct the gamma function (see Eq. 17) for several values of vs . At criticality, the gamma function becomes independent of N and N’ and the curves will cross. In Fig.2 we show the gamma function we obtain by using different number of elements and how it begins to cross for various approximations to the energy.

  3. In order to go to the complete basis set limit , we have to let . The best choice is to keep minimal nightingale1 (); privman (). Then, we proceed as follows: We solve for (N-1, N , N+1). This provides us with two gamma curves: and , with minimum . The Crossing of these two curves gives us the pseudo-critical parameters and with a truncated basis set number N. Fig. 3 illustrates the N dependence of the crossing.

  4. In order to extrapolate the values and to the infinite N limit, we use the algorithm of Burlish and Stoer BS (). Fig.s 4 and 5 show their behavior by systematically increasing N for the mean-field approximations and exact formulation, respectively.

    (a)     (b)
    (c)     (d)
Figure 2: (Color online) Plot of (Eq.  17) obtained by using FSS method as a function of the charge Z for increasing number of elements (N) for: (a) Hartree-Fock only, (b) Total Energy (HF + correlation energy), (c) LDA approximation, and (d) exact formulation

While finding values and and extrapolating to the limit some numerical fluctuations can arise depending on the approximation used. We discuss some ways of dealing with these fluctuations. First, one can improve the basis used from -continuous polynomials to -continuous polynomials in order to obtain greater accuracy in the calculations, however, matrices double in size increasing the memory cost as well as the computing time. Second, one can construct a gamma function from a potential in the Hamiltonian that explicitly takes into account all variables. For example in the three dimensional (3D) case we found that it is better to use the transformed Hamiltonian (Eq.39) where the perturbation is given by rather than . This does not change the physics of the system since both expressions can be used to find the critical value: either or , which corresponds to one of the electrons going from a bound state into the continuum. The first perturbation takes into account all three variables when the expectation value is computed whereas the latter introduces numerical errors in the expectation value.

Figure 3: The crossing of and as N is increased using Total Energy approximation in HF. The crossing has been slightly enhanced in order to demonstrate the dependence of N.
HF   LDA Total Energy Exact Reference Baker (); PSK ()
critical charge
1.03114   0.92808 0.90946 0.91857 0.91103
critical exponent
0.99971   1.00011 1.00009 1.00082 1.00000
Table 2: Critical parameters extrapolated from Fig.(4) and Fig.(5)
          (a)           (b)           (c)
Figure 4: Extrapolated values for the critical exponent (top) and (bottom) for: (a) HF approximation, (b) LDA approximation, and (c) Total energy. The filled dots are the extrapolated values for the limit
Figure 5: Extrapolated values for the critical exponent and in the limit solving the system exactly in 3D. We construct the function using Hermite interpolation polynomials for two potentials: (denoted by ) and (denoted by ); the latter is obtained by making the transformation in the Hamiltonian. We checked that the fluctuations seen in the figure vanish by replacing the numerical approximation to the binding energy with the analytical form, , while the extrapolated value remains unchanged.

Last, another method is to increase the minimum difference used in finding a crossing between and until a systematic crossing is obtained for all values of N; in the LDA approximations we used , since the Gamma curves change very little as the number of elements increases (see Fig.2-c), whereas in the 3D formulation we used since in this case the number of elements increases as , and was kept fixed at three. In solving the self-consistent equations [ Eqs.(23) and (29)], we found that the region below became problematic when updating the self-consistent equations. This is because as the nuclear charge decreases, the wavefunction decay becomes slower with respect to r and the cutoff radius might lead to unphysical solutions (i.e. solutions which depend on ). To get around this problem, we found it is useful to perturb the wavefunction with the arithmetic average of the current and predicted wavefunctions, rather than just the predicted wavefunction. In this manner, we can avoid converging to a local minima that corresponds to an unphysical solution.

In phase transition theory it is well known that using mean field theory (for dimensionality ) does not give correct values near a critical point. This is indeed the case by using only the Hartree-Fock approximation (see Table 2 ), however we see that incorporating the exchange-correlation energy within the LDA approximation improves substantially the critical charge within 2% of the accepted value and within 1% using Total Energy approximation and exact formulation. It is interesting to note that all approaches give a consistent estimate for the critical exponent, regardless of the approximation used. This study demonstrates that Finite-element method is a viable tool that can be applied in combination with finite size scaling in order to obtain good estimates of critical parameters starting with ab initio calculations for electronic calculations such as DFT and HF and compare well with other methods that are deemed more precise.

We now check the consistency of our approach. In the theory of continuous phase transitions, a quantity diverging in the thermodynamic limit (), becomes a regular function for finite N, viz . We can check that this assumption is indeed satisfied for any value N. In the infinite basis limit we have the following form near criticality: (Eq.11). In this manner, can obtain the derivative of the energy near criticality given that we have a limiting expression for the energy, viz


In a finite basis the above expression becomes a power law in the number of basis with the following characteristic exponents.


Excluding the trivial case N=0, one can therefore use Taylor expansion around .


In this manner we obtain a universal scaling function . Plotting vs the argument of the function would confirm whether the scaling function depends on N or not. Fig. 6 shows that all quantities obtained from different N values do in fact collapse into one curve, and this confirms the scaling hypothesis for FSS with FEM (Eq.10). In order to make the data collapse, it should be emphasized that the values and used to produce Fig.6 correspond to those obtained in Table 2. If other values are used, the curves will not collapse and they will be distributed in a family of curves (N dependent) all deviating away from the unique curve shown in Fig.6. The correlation length exponent was varied until the best collapse curve was found. This occurred around for all implementations.

          (a)           (b)
          (c)           (d)
Figure 6: (Color online) Data collapse for various values of N using: (a) HF, (b) Total Energy, (c) LDA approximation, and (d) exact formulation. The value from this collapse study that give us the best correlation length exponent , it should be emphasized that the values and used in making these figures are the corresponding values for each approximation (i.e. those shown in Table  2)

The solution to a set of Hartree-Fock or Kohn-Sham equations for large number of atoms becomes the limiting step, and it is here that the FEM can extend the range of sizes that can be investigated by other current means. The efficiency of FEM in the context of large-scale calculations derives from the strict locality of the FEM basis in real space, thus minimizing the need for extensive interprocessor communications on massively parallel computational architectures pask (); pask2 (). The work of Pask et al pask (); pask2 () has formulated a real-space DFT approach for condensed matter applications without the need of using expensive Fourier transforms. Impressive results are obtained for GaAs calculations using this approach on DFT+FEM that systematically converges, as the number of elements is increased, to the values obtained using the conventional plane wave calculations, employing Fourier transforms and Ewald summations.

The scope of the dimension attainable with FEM has also been demonstrated by Sterne et al Sterne (). They have employed FEM to compute positron distributions and lifetimes for various defect structures. Because the distribution is non-electron-like (concentrating away from the nuclei rather than around them), the generality of the FEM basis makes it particularly well suited for such problems. In order to study defects they have modeled a 5488 atom cell; as far as we know, this is one of the largest such calculations reported by any method to date.

Another applications of FEM to large-scale ab initio calculations is the work of Tsuchida et al on molecular dynamics simulations of liquids TT (). They apply FEM study of liquid formamide (HCONH2) using a generalized gradient approximation to the exchange and correlation terms to improve the modeling of hydrogen bonds. The use of adaptive coordinate transformations helped to double the resolution of the basis at the C, N and O centers. The liquid was modeled by 100 molecules in a cubic supercell of side 35.5 a.u., for a total of 600 atoms. The scope of this simulation is of the order of the largest one performed by standard PW methods to date and demonstrates the capacity of FEM for such large-scale computations already.

A new promising approach developed recently by Mazziotti Mazziotti1 () is the large-scale semidefinite programming for many-electron quantum mechanics. In this approach, the energy of a many-electron quantum system can be approximated by a constrained optimization of the two-electron reduced density matrix (2-RDM) that is solvable in polynomial time by semidefinite programming (SDP). The developed SDP method for computing strongly correlated 2-RDMs is 10-20 times faster than previous methods Mazziotti2 (). This approach was illustrated for metal-to-insulator transition of H, with about variables.

These examples along with our current work, suggest that we can combine the Finite-element method with finite size scaling to obtain critical parameters in extended systems. One possibility is to calculate the phase transition at zero temperature in the 2-dimensional electron gas. The ground state energy of the system favors a crystallize state as the density of the electrons is lowered. Tanatar and Ceperey TC () performed a variational Monte Carlo simulations to the total energy using Hartree-Fock approach


They find the gas-to-crystal transition at . This system would be ideal to test FEM with finite size scaling for quantum criticality in extended systems.

Vi Conclusion

The theory of finite size scaling (FSS) has been demonstrated to be of great utility in analyzing numerical data of finite systems fisher (); widom (); barber (); privman (); cardy (); nightingale1 (); Peter1 (); Peter2 (). Moreover, it has been demonstrated that the theory also works for quantum systems, where the formalism applies to a complete basis expansion in the Hilbert space neirotti0 (); serra2 (); kais (); snk1 (); snk2 (); nsk (); qicun (); kais1 (); review (); adv (). In this study we implemented the finite-element method which uses a local basis expansion consisting of arbitrary polynomials in real space and with variable resolution. It was demonstrated that this method can be used to compute matrix elements (used in the finite size scaling formalism) accurately and with monotonic convergence with respect to the system size N. We have implemented the combined approach on different approximations to electronic calculations such as Hartree-Fock approximation (+ electron-electron correlation) and density-functional theory (DFT) under the local density approximation and found that even at the simplest level in DFT, the results we have obtained are in good agreement with other more accurate estimates such as large order perturbation theory Baker ().

The finite-element method has already shown great promise for extended systems since it uses local basis expansion to yield sparse matrices that lend the methodology easily to parallel architectures pask (); AHM (). Furtermore, the Khon-Sham equations in DFT provides a general method for extended system. Therefore we expect this combined approach can be applied to electron structure in atoms, molecules, and extended systems without loss of generality. Nevertheless, more studies need to be done to confirm the applicability of FSS with FEM for larger systems, where this combined approach can be used to shed some light into the nature of phase transitions.


We would like to acknowledge the financial support of ARO and B. W.-K. thankfully acknowledges support form the NSF under the grant PHY-0969689. We would also like to thank the referee for his helpful feedback and for pointing out the form for the convergence error in FEM.


  • (1) M.P. Nightingale, Physica 83A, 561 (1976).
  • (2) J.L. Cardy, Finite-Size Scaling, (Elsevier Science Publishers B.V. New York 1988).
  • (3) M.E. Fisher, in Critical Phenomena, Proceedings of the 51st enrico Fermi Summer School, Verena’s, Italy, edited by M.S. Green (Academic Press, New York, 1971); M.E. Fisher and M.N. Barber, Phys. Rev. Lett. 28, 1516 (1972)
  • (4) B. Widom, Critical Phenomena in Fundamental Problems in Statistical Mechanics Ed: E.G.D. Cohen, (Elsevier Publishing Company, NY 1975).
  • (5) M.N. Barber, Finite-size Scaling, in Phase Transitions and Critical Phenomena Vol. 8. C. Domb and J.L. Labors eds. (Academic Press, London 1983).
  • (6) V. Privman, Finite Size Scaling and Numerical Simulations of Statistical Systems, ( World Scientific,Singapore 1990).
  • (7) P. J. Reynolds, H. E. Stanley, and W. Klein, J. Phys. A 11, L199 (1978).
  • (8) P. J. Reynolds, H. E. Stanley, and W. Klein, Phys. Rev. B 21, 1980 (1223).
  • (9) W. Moy, S. Kais, and P. Serra, Molecular Physics 106, 203 (2008).
  • (10) J. P. Neirotti, P. Serra, and S. Kais, Phys. Rev. Lett. 79, 3142 (1997).
  • (11) E. Antillon, W. Moy, Q. Wei and S. Kais, J. Chem. Phys. 131, 104105 (2009).
  • (12) E. Tsuchida and M. Tsukada, Phys. Rev. B 54, 7602 (1996).
  • (13) D. Stiles (private communication).
  • (14) J.E. Pask and P.A. Sterne, Modelling Simul. Mater. Sci. Eng. 13, R71 (2005).
  • (15) L. R. Ram-Mohan, Finite Element and Boundary Element Applications in Quantum mechanics, (Oxford Uinv. Press, 2002).
  • (16) D.W. Pepper and J.C. Heinrich, the Finite element method, (Taylor and francis group, New York 2006).
  • (17) O. Taisuke; T. Masayuki, Comp. Phys. Commun., 182, 1245 (2011).
  • (18) J.E. Bylaska; M. Holst; J.H. Weare. J. Chem. Theor. Comp. 5, 937 ( 2009).
  • (19) R. Alizadegan, K. J. Hit, and T. J. Martinez J. Chem. Phys. 132, 034101 (2010).
  • (20) J.E Pask, B.M. Klein, C.Y. Fong, P.A Sterne, Phys Rev B 59 12352 (1999).
  • (21) C.N. Yang and T.D. Lee Phys. Rev. 87 404 (1952).
  • (22) T.D. Lee and C.N. Yang Phys. Rev. 87 410 (1952).
  • (23) S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, 1999).
  • (24) J.P. Neirotto, P. Serra, and S. Kais,Phys. Rev. Lett., 79, 3142 (1997).
  • (25) P. Serra, J.P. Neirotti, and S. Kais, Phys. Rev. Lett. 80, 5293 (1998).
  • (26) S. Kais, J.P. Neirotti, and P. Serra, Int. J. Mass Spectrometry 182/183, 23 (1999).
  • (27) P. Serra, J.P. Neirotti, and S. Kais, Phys. Rev. A 57, R1481 (1998).
  • (28) P. Serra, J.P. Neirotti, and S. Kais J. Phys. Chem. A, 102 47, 9518 (1998)
  • (29) J.P. Neirotti, P. Serra, and S. Kais, J. Chem. Phys. 108, 2765 (1998).
  • (30) Q. Shi and S. Kais, Mol. Phys. 98, 1485 (2000)
  • (31) S. Kais and Q Shi, Phys. Rev. A62,60502 (2000).
  • (32) S. Kais, and P. Serra, Int. Rev. Phys. Chem. 19, 97 (2000).
  • (33) S. Kais, and P. Serra, Adv. Chem. Phys., 125, 1 (2003).
  • (34) P. Serra, and S. Kais, Chem. Phys. Lett. 372, 205-209 (2003).
  • (35) A. Ferron, P. Serra and S. Kais, J. Chem. Phys. 120, 8412-8419 (2004).
  • (36) Q. Shi, and S. Kais, Mol. Phys. 98, 1485-1495 (2000)
  • (37) Q. Shi, and S. Kais, Int. J. Quantum Chem. 85, 307-314 (2001).
  • (38) A. Ferron, P. Serra, and S. Kais, J. Chem. Phys. 128, 044307 (2008).
  • (39) B. Simon J. Functional Analysis 25, 338 (1977).
  • (40) L.D. Landau and E.M. Lifshitz, Quantum Mechanics: Non-Relativistic Theory. Vol. 3 (3rd ed.) (Butterworth-Heinemann, London, 1981).
  • (41) S. Chandrasekhar, Astrophysical Journal, 100, 176 (1944).
  • (42) F. H. Stillinger and D. K. Stillinger, Phys. Rev. A 10, 1109 (1974).
  • (43) J. D. Baker, D. E. Freund, R. D. Hill, and J. D. Morgan, III. Phys. Rev. A 41, 1247 (1990)
  • (44) J.M. Thijssen, Computational Physics. (Cambride University Press 1999).
  • (45) W. Schweizer. Numerical Quantum Dynamics (Kluwer Academic Publishers, Dordrecht, 2001).
  • (46) G. Giuliani, G. Vignale, Quantum Theory of the Electron Liquid. (Cambride University Press 2005)
  • (47) H. Mitler, Phys. Rev. 99, 1835 (1955).
  • (48) C. W. Ufford,J.G. Thomas Phys. Rev. 133, A121 (1964).
  • (49) W. Kohn, L. Sham, Phys. Rev. 140, A1133 (1965).
  • (50) C. J. Umrigar and X. Gonze, Phys. Rev. A 50, 3827 (1994)
  • (51) J.P. Perdew, Y. Wang, Phys. Rev. B 45, 13244 (1992)
  • (52) D. Porezag,
  • (53) G. Breit, Phys. Rev. 35, 569 (1930)
  • (54) F. S. Levin and J. Shertzer, Phys. Rev. A 32, 3285 (1985).
  • (55) F. M. Gomes
  • (56) C. Geuzaine and J. Remacle.
  • (57) R. Bulirsch, J. Stoer, Num. Math, 6, 413 (1964)
  • (58) P.A. Sterne, J.E. Pask and B.M. Klein Appl. Surf. Sci. 149, 238 (1999).
  • (59) E. Tsuchida , J. Chem. Phys. 121, 4740 (2004).
  • (60) D. A. Mazziotti, Phys. Rev. Lett. 106, 083001 (2011).
  • (61) D. A. Mazziotti, Phys. Rev. Lett. 93, 213001 (2004).
  • (62) B. Tanatar and D. M. Ceperley, Phys. Rev. B 39, 5005 (1989).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description