Real-space grid representation of momentum and kinetic energy operators for electronic structure calculations
We show that the central finite difference formula for the first and the second derivaxtive of a function can be derived, in the context of quantum mechanics, as matrix elements of the momentum and kinetic energy operators using, as a basis set, the discrete coordinate eigenkets defined on the uniform grid . Simple closed form expressions of the matrix elements are obtained starting from integrals involving the canonical commutation rule. A detailed analysis of the convergence toward the continuum limit with respect to both the grid spacing and the approximation order is presented. It is shown that the convergence from below of the eigenvalues in electronic structure calculations is an intrinsic feature of the finite difference method.
The last decade has witnessed a growing interest in electronic structure computational schemescol (2015); Beck (2000); Enkovaara et al. (2010); Michaud-Rioux et al. (2016); Hernández et al. (2007) based on the real-space finite difference representation of both the momentum and kinetic energy operators. The main advantage of the method when applied to either the Schrödinger or the DFT Kohn and Sham equationsKohn and Sham (1965), particularly within the pseudopotential density functional theoryChelikowsky et al. (1994), is the highly local structure of the relevant matrices. This feature allows for both a significant reduction of computer memory and for the use of efficient parallel algorithms Burdick et al. (2003) paving the way toward grid-based linear scaling methods. Ab initio electronic structureBeck (2000); Chelikowsky et al. (1994); Fattebert and Bernholc (2000); Castro et al. (2006); Andrade et al. (2015) and molecular dynamics codesSchmid (2004) have successfully been implemented and tested. Moreover, adaptive mesh refinementsModine et al. (1997) and multigrid strategiesBriggs et al. (1996) can be used for gaining local resolution only where it is needed. Finally, it is worth mentioning that the finite difference method is also a very good conceptualDatta (2005) and practicalFujimoto and Hirose (2003); Khomyakov and Brocks (2004) tool for the study of electric current flow in nanodevices.
Other electronic structure computational methods closely related to real space grids are finite elements schemesWhite et al. (1989), discrete variable representationsLiu et al. (2003), Lagrange meshesVarga et al. (2004), nonorthogonal generalized Wannier functionsSkylaris et al. (2005), and waveletsArias (1999); Genovese et al. (2008). All these methods have in common a real space grid whose points are associated to localized basis functionsBaye (2006). As such, it seems that there is not connection between this class of methods and finite differences. We shall show in the following that indeed there is a link when the continuum limit in the finite differences is approached.
For entering the core of the problem discussed here, let us briefly recall how the finite difference method works. The starting point is the choice of the way in which the kinetic energy operator is represented on a discretized space. Let us consider a one dimensional grid (the extension to two and three dimensional Cartesian grids is trivial) made of equally spaced points defined by with . The lowest-order approximation to the second derivative is obtained considering only three grid points and . With a Taylor series it is easy to see that the lowest order approximation of second derivative at is given by
Higher order expressions involving wave function values at can be obtained by means of specific algorithmsFornberg (1998). In general, a finite difference representation of the second derivative takes the form
where the integer number , which we define as the representation order, controls the accuracy. Whichever method is used for determining the coefficients , it is clear that the kinetic energy representation induced by Eq. (2) does not apparently have an explicit basis set. In other words, although the can formally be viewed as the entries of the kinetic energy matrix, it is not at all clear whether or not there exists an underlying basis set. This is, in our opinion, an interesting and fertile point that has not been discussed in the literature.
A further point is represented by the convergence properties toward the continuous space, that is, the limit of vanishing grid spacing. For instance, practical finite difference DFT electronic structure calculations show that the total energy converges from below with respect to both the grid spacing and the order Beck (2000); Maragakis et al. (2001); Skylaris et al. (2002). It is generally believed that this is an important limitation preventing from developing convergence schemes.
The core of this work is first of all to show that the finite difference method has discrete coordinate eigenkets as the underlying basis set. We obtain this result by explicitly constructing the matrix elements of both the linear momentum and kinetic energy operators. Interestingly, the unique starting point of the derivation is the use of integrals involving the canonical commutations between the operators and . In a way, we leave to the quantum mechanics rules the determination of the optimal matrix elements. We show that it is possible to obtain simple analytical expressions for the matrix elements so that their convergence towards the continuum limit and the connection with other real space methods can be discussed in detail. The theory is complemented with some explanatory numerical examples.
Ii The momentum operator
Let us associate to each grid point the ket which we assume to be an eigenket of the position operator . This means that
Any state defined on this discretized space can be written as
Although it is not strictly necessary for the present discussion, we can notice here that the numbers in Eq. (4) are the values of the wave function on the grid points .
A grid representation of the momentum operator is given by
where is a generic ket. The real numbers are the momentum matrix elements. For a given there are grid points over which the are different from zero. This means that the set of depends on . For indicating explicitly this dependence, we write, from now on, and call the representation order. Taking , from Eq. (5) we have
where . Since we want an Hermitian matrix representing , we must require
and this implies that . Additional conditions necessary for determining can be derived starting from integrals involving the canonical commutation
With continuum position eigenkets it is known that
where we have used Eq. (7) for limiting the sum to positive . Equation (11) is just enough for determining with . Indeed, one immediately gets which are exactly the central difference weights for the lowest order first derivative. This an interesting intermediate result. However, the issue now is to show how the most general case of can be handled, that is, how we can generate the necessary additional equations. Since multiple commutators of the type are all zero, we have
where we have defined
The momentum matrix elements are
The numerical values of Eq. (17) are identical to those obtained for the first derivative central finite difference weights using the algorithm discussed in Ref.Fornberg (1998). It is useful for the following to calculate the limit of Eq. (16) when . The result is extraordinarily simple
from which we immediately have
The actual meaning and implications of this result will be discussed in the next section.
Before closing this section, it may be useful to summarize the main points. The quantum mechanical ’rules’ expressed in Eq. (9) and (12) have guided us to the simple closed form expression (15) of the weights . Since our derivation has been built on the basis set , we have indirectly proved that this is the basis set of the finite difference scheme.
Iii The kinetic energy operator
We now derive an expression for the kinetic energy matrix elements following steps that are similar to what we have done for the linear momentum. The kinetic energy can be defined as
From Eq. (20) and the orthogonality of the set it is easy to see that
The coefficients in Eq. (21) can be related to taking the matrix elements of on the basis . The final result is
Eq. (23) gives a closed form expression for the kinetic energy valid to all orders . The advantage of having an analytical expression of the matrix elements is not only the simplicity with which they can be generated in an electronic structure computer code, but also the possibility of easily looking at their behavior when the representation order is very large. Indeed, from Eq. (18) it is easy to see that when Eq. (23) becomes
It should be noted that although Eq. (24) is strictly valid only for a grid made of an infinite number of points, it is nevertheless an interesting results because it establishes a connection between finite difference and pseudospectral methodsJordan and Mazziotti (2004); Baye (2006). It is worth recalling that in a pseudospectral method the wave function is written asBaye (2006)
where the basis has the cardinal property
Because of this property, Eq. (25) becomes
Although there is a similarity between Eqs. (27) and (4), it should be observed that while in the former the basis and the wave function are defined in a continuous space, in the latter they are defined only on the grid points. In a way we can say that Eq. (27) is the continuum space extension of Eq. (4).
Let us now calculate the matrix elements of the kinetic energy operator on the (also known as Whittaker cardinal) functions
It is quite interesting that the result of this calculationJordan and Mazziotti (2004); Baye (2006) is exactly Eq. (24). In a similar way, it is possible to show that the momentum matrix elements between functions are given by Eq. (19). We therefore see that the connection between the finite difference and the pseudospectral methods goes through the infinite order limit .
Iv The convergence towards the continuous space
For having more insights on the convergence properties with respect to it is necessary to give a sense to the limit on a space grid made of a finite number of points. As it is well known in solid state theory, an effective way of constructing an infinite lattice out of a finite one is that of taking periodic replicas. We can use exactly the same idea by taking periodic replicas of a grid of length . In this way the limit has a precise meaning.
On a replicated grid we need to impose periodic boundary conditions. We shall see in a moment that this choice is particularly useful because it allows the analytical calculation of both the momentum and the kinetic energy eigenvalues. Let us show this for the momentum matrix. The eigenvalue equation is
where we have indicated with the momentum eigenvectors on the grid. It is a simple matter to verify that Eq. (29) can be solved with
The wave vector is selected imposing the periodic boundary condition giving the well known result
The simple structure of Eq. (32) allows us to study in detail the convergence towards the continuum with respect to both the representation order and the grid spacing .
However, before doing that, let us also see what we obtain for the kinetic energy. From an eigenvalue equation similar to Eq. (29) we get
This is quite interesting: it tell us that with a fixed grid spacing , the linear momentum gets exactly the values one would expect for a free particle on a lattice of length with periodic boundary conditions. The type of convergence behind Eq. (34) can be better appreciated by computing the power series of Eq. (32) with respect to the grid spacing . A lengthy calculation in which we retain, for each order , the first non vanishing terms gives
It should be noted that since is always positive, Eq. (35) shows that the convergence to the exact momentum is from below.
where we have defined
The same conclusions we have drawn for the momentum hold in this case. In particular, it is worth stressing the convergence from below which plays an important role in electronic structure calculations.
To give a hint on how small the the relative error are, we show in Fig. 1 a plot of Eqs. (36) and (39) as a function of . It can be seen from this figure that quickly goes down so that one would expect that in a practical numerical calculation there should be no need of using values of larger than about . This is important because keeping a banded kinetic energy matrix is a prerequisite for a fast calculation.
V Numerical examples
The convergence properties discussed in the previous section are here verified on two very different numerical examples. The first one is the one dimensional Schrödinger equation with the potential for which an exact analytical solution is availableLandau (1976). Choosing and , it can be seen that the bound states are just three.
Using the kinetic energy matrix elements of Eq. (23), we have numerically calculated the three eigenvalues changing both and . The convergence towards the exact eigenvalues on changing is shown in Fig. 2 while that obtained on changing is shown in Fig. 3. Moreover, in the inset of Fig. 2 we also show the convergence of the eigenvalues sum. The data in Fig. 2 show that a monotonic convergence from below is recovered only for sufficiently large . On the contrary, as shown in Fig. 3, the convergence with respect to is fully monotonic. The difference is due to the way the potential contributes to the final result. In Fig. 2 we change the grid spacing so that the potential is sampled differently for each whereas in Fig. 3 we are only changing the representation order in the kinetic energy matrix. The convergence from below shown in Fig. 3 is consistent with the data of Fig. 1.
The second numerical example we have analyzed is that of a DFT calculation of the ground state of a C molecule. The calculations have been performed using the finite difference pseudopotential method developed by Chelikowsky et alChelikowsky et al. (1994). In Fig.4 we show the convergence of the total energies with respect to the representation order while in Fig.5 we show the convergence versus the grid spacing. In both the cases the convergence from below due to the kinetic energy is quite evident. Despite the enormous difference in the computational complexity with respect to the one dimensional case, the general trends are similar.
A last example is here reported to evaluate, in bulk silicon, the convergence properties of the total energy with respect to the grid spacing and the representation order. We used Octopus scientific code to perform DFT real space calculationsCastro et al. (2006).
A 2x2x2 Monkhorst Pack grid was adopted to perform integrations over the first Brillouin zone, together with a Perdew-Zunger LDA exchange-correlation functional, and a 8-atom cubic conventional cell, with a lattice constant fixed at a.u. The total energy was evaluated on changing the representation order and the grid spacing, each time redoing the selfconsistent procedure. The convergence of the total energy with respect to the grid spacing, for several values of , is reported in figure 6. The scaling is qualitatively similar to the example illustrated in the previous section, and the results discussed above are confirmed. Even in this case, the convergence is from below. The total energy is fully monotonic with respect to the representation order (at fixed grid spacing), and shows a striking similarity to figure 4 (for this reason the graph is not reported here). Instead, as reported in figure 6, the convergence is always from below, but with oscillating features, when the grid spacing decreases at a fixed value of . It is worth noticing that a higher representation order leads to a faster convergence with respect to the grid spacing. A tradeoff, in terms of computational time, has to be found between representation order and grid spacing. However, the use of a high value of is surely convenient. The default value chosen by the code for the representation order is .
With this work we have made an attempt of resolving two main misconceptions related to finite difference electronic structure calculations. The first is that it is not really true that a finite difference representation has not a basis set and the second is that the anti variational behavior is not a limitation of the method.
The proof that a basis set does exist has been a constructive one in the sense that we have derived explicit expressions for the momentum and kinetic energy matrix elements showing that they are identical to the central finite difference weights of the first and second derivatives.
The convergence from below has been fully characterized calculating the momentum and kinetic energy eigenvalues using periodic boundary conditions. We have derived for both the operators two formulas, Eqs. (36) and (39), that we think can be useful for developing convergence criteria. We hope to have made clear that the anti variational behavior is an intrinsic feature of finite difference method intimately connected to the fact for each couple of the integer number and one has a different representation of the operators.
An important point of this work is that the continuum momentum and kinetic energies can be obtained with two different limits. The first one is obtained taking the limit of large values of , the representation order, with a fixed grid spacing. This limit, taken on a finite grid, has a meaning only when periodic boundary conditions are imposed. The exact linear momentum and kinetic energy of a free particle on a lattice are recovered. The second limit consists in sending the grid spacing to zero while keeping unchanged. We have shown that also in this case the continuum momentum and kinetic energy are recovered.
A concluding remark may be that the above result may be useful for reviewing, from a different perspective, the quantum mechanics in a discretized space.
- Themed collection: Real space numerical grid methods in quantum chemistry, Phys. Chem. Chem. Phys. 2015, 17.
- Beck, T. L. Rev. Mod. Phys. 2000, 72, 1041–1080.
- Enkovaara et al. Journal of Physics: Condensed Matter 2010, 22, 253202.
- Michaud-Rioux, V.; Zhang, L.; Guo, H. Journal of Computational Physics 2016, 307, 593 – 613.
- Hernández, E. R.; Janecek, S.; Kaczmarski, M.; Krotscheck, E. Phys. Rev. B 2007, 75, 075108.
- Kohn, W.; Sham, L. J. Phys. Rev. 1965, 140, A1133–A1138.
- Chelikowsky, J. R.; Troullier, N.; Saad, Y. Phys. Rev. Lett. 1994, 72, 1240–1243.
- Burdick, W.; Saad, Y.; Kronik, L.; Vasiliev, I.; Jain, M.; Chelikowsky, J. Comput. Phys. Commun. 2003, 156, 22–42.
- Fattebert, J.-L.; Bernholc, J. Phys. Rev. B 2000, 62, 1713–1722.
- Castro, A.; Appel, H.; Oliveira, M.; Rozzi, C. A.; Andrade, X.; Lorenzen, F.; Marques, M. A. L.; Gross, E. K. U.; Rubio, A. Phys. Status Solidi B 2006, 243, 2465–2488.
- Andrade et al. Phys. Chem. Chem. Phys. 2015, 17, 31371–31396.
- Schmid, R. J. Comput. Chem 2004, 25, 799.
- Modine, N. A.; Zumbach, G.; Kaxiras, E. Phys. Rev. B 1997, 55, 10289–10301.
- Briggs, E. L.; Sullivan, D. J.; Bernholc, J. Phys. Rev. B 1996, 54, 14362–14375.
- Datta, S. Quantum transport atom to transistor; Cambridge University Press: Cambridge, UK, 2005.
- Fujimoto, Y.; Hirose, K. Phys. Rev. B 2003, 67, 195315.
- Khomyakov, P. A.; Brocks, G. Phys. Rev. B 2004, 70, 195402.
- White, S. R.; Wilkins, J. W.; Teter, M. P. Phys. Rev. B 1989, 39, 5819–5833.
- Liu, Y.; Yarne, D. A.; Tuckerman, M. E. Phys. Rev. B 2003, 68, 125110.
- Varga, K.; Zhang, Z.; Pantelides, S. T. Phys. Rev. Lett. 2004, 93, 176403.
- Skylaris, C.-K.; Haynes, P. D.; Mostofi, A. A.; Payne, M. C. J. Chem. Phys. 2005, 122, 084119.
- Arias, T. A. Rev. Mod. Phys. 1999, 71, 267–311.
- Genovese, L.; Neelov, A.; Goedecker, S.; Deutsch, T.; Ghasemi, S. A.; Willand, A.; Caliste, D.; Zilberberg, O.; Rayson, M.; Bergman, A.; Schneider, R. J. Chem. Phys. 2008, 129, 014109.
- Baye, D. Phys. Status Solidi B 2006, 243, 1095–1109.
- Fornberg, B. SIAM Rev. 1998, 40, 685–691.
- Maragakis, P.; Soler, J.; Kaxiras, E. Phys. Rev. B 2001, 64, 193101.
- Skylaris, C.-K.; Diéguez, O.; Haynes, P. D.; Payne, M. C. Phys. Rev. B 2002, 66, 073103.
- Jordan, D. K.; Mazziotti, D. A. J. Chem. Phys. 2004, 120, 574.
- Landau, L. D. Quantum Mechanics; Pergamon Press, Oxford, 1976.