The Hypernetted Chain Equations for Periodic Systems.

Starting from the general inhomogeneous Fermi hypernetted chain equations, the equations for periodic systems are derived by simple Fourier transform. It is shown how the symmetry reduces the size of the involved quantities. First results for a one-dimensional (1D) model system are presented. The results allow a reliable estimation of the numerical demand even for realistic 3D systems, such as solids. It is shown that treatment of this systems is feasible with moderate computational resources.

Introduction

Jastrow correlated trial wave functions, i.e., a wave function which additionally to a Slater determinant includes two-particle correlations, are widely used in quantum Monte Carlo (MC) techniques [1, 2]. Less well known are the analytic methods to calculate expectation values with this wave function [3-6]. The Fermi hypernetted chain method (FHNC) is an effective scheme to calculate a large class of cluster diagrams. The optimal correlation function is obtained by minimization of the energy expectation value. A drawback of the method, compared to MC techniques, is that certain diagrams, i.e., the elementary diagrams, are not covered by the scheme. Therefore, the obtained quantities and matrix elements are only approximations. However, the advantage of the method is that it is generally numerically less demanding than MC methods. The analytic representation of the wave functions further allows to deal with excited states [7, 8].

An additional advantage of the used formulation of the FHNC method is that the functional form of the method allows a parameter-free optimization of the two-particle correlation function.

The FHNC method has been first developed for homogeneous systems, but has been generalized to inhomogeneous systems by Krotscheck et al. [9-12]. In order to keep the numerical demand low, one has to utilize the symmetries of the system. This has been done for slab [13] and for spherical geometries.

In this paper a possible discretization for periodic systems is presented which reduces the numerical demand considerably and leads to a simple form of the Euler equation. A different discretization and coordinates are proposed by Krotscheck [14]. Although it allows a further reduction in numerical demand by reducing the resolution of the center of mass coordinate, it sacrifices the simplicity of the relevant equations and reduces the resolution and consequently the accuracy of the result.

This work is related to an implementation of Fantoni and Schmidt [15], which formulates and solves the FHNC equations in a periodic box, but with a uniform density. The presented method extends this approach by explicitly including a non-constant density.

The Inhomogeneous HNC Equations

Here the generalization to the general inhomogeneous case of the FHNC equations in the formulation of Kallio and Piilo [16] is given. The reason for this is that results are shown for a 1D model system and the simplified FHNC of Krotscheck doesn’t work in 1D.

The inhomogeneous theory yields coupled equations for the one-particle and two-particle correlations. The one-particle equation is a generalized Hartree–Fock equation (gHF) as given in [12], and the two-particle equation is actually a set of equations given below. These two equations are coupled and need to be solved self-consistently. However, it seems to be a reasonable approximation for certain systems to omit the coupling [13]. This is the justification for using a model for the single-particle states in the next section and not solving the gHF.

The inhomogeneous Euler–Lagrange equations can be written aswith the definition of the convolutionand where S is the, to be determined, static structure function. is the particle hole irreducible interaction explicitly given below. The one-body operator is defined aswhere is the one-body density obtained from the gHF.

The induced interaction iswhere the from the previous iteration is used. (In order to start the iteration process, a reasonable guess for is needed.)

With these ingredients and by using the pair distribution function , the new particle-hole irreducible interaction is computedwhere is short for the gradient with respect to and v is the interaction potential of the particles in the system. Note that here no convolutions are involved, i.e., the functions are simply multiplied. The Fermi potential is defined aswithThe free static structure function is given bywhere l is the non-interacting density matrix as result of the gHF and is the spin degeneracy and is the free pair distribution function analog of .

In order to obtain the boson version of the equations has to be set to zero.

Solution of the Euler equation: Equation (1) is solved by considering the eigenvalue problemThe are orthonormalized according toThe static structure factor is obtained byand its inverse

The Equations for the Periodic System and Results for a 1D Model

Definition of the Fourier Transform

In order to describe a periodic system a unit cell, spanned by the lattice vectors , where d is the dimensionality of the system, is defined. One-body quantities like the density are periodic: . With the translation vector . A similar relation also holds for two-body quantities, e.g., the pair distribution function: . This symmetry is used to constrain the first coordinate to the basic unit cell. This is denoted by a bar: . Further, we may decompose the second coordinate: . In order to describe pair properties properly we have to choose a cutoff in the translations. Therefore, we restrict the integers to the interval ,1 which we call now crystal on which we impose periodic boundary conditions. (Although this is non-physical, it becomes a good approximation if N is large enough.) This crystal consists of unit cells. The unit cell is discretized with points in each direction, so the resolution is . A periodic function can be expressed as a sum of plane waveswhere are the reciprocal lattice vectors, with . For the discretized function in with integer, the reciprocal lattice sum (13) is limited by .

By utilizing the periodicity of the system the Fourier transform is defined asand its inverseWith this definition the convolution as defined above becomes in momentum spaceWith the above definitions the Euler equation is written in reciprocal space asThe induced interaction is also calculated most effectively in reciprocal space, by Fourier transform or Eq. (4). However, the particle hole irreducible interaction is calculated in r-space, like in the homogeneous case.

The most time-consuming part is the solution of the Euler equation. The Fourier transform reduces the eigenvalue problem Eq. (9) considerably since appears as a parameter:so the computation timescales are linear with the number of q-vectors.

Results for a 1D Model System

The described method is applied to a 1D model system similar to that described by Asgari [17]. The interaction potential is derived from the electron gas in a homogeneous trap where b characterizes the thickness of the wire. Since the single-particle equation is not solved, a simple model for the single-particle basis is usedThese single-particle states result in a sinusoidal modulated density with the amplitude determined by . is the reciprocal lattice vector and a the length of the unit cell.

Results of this model are plotted in Fig. 1 with the parameters and the Wigner–Seitz radius and Å the Bohr radius. In order to converge the result a real space resolution of was sufficient for the homogeneous and the inhomogeneous implementation. This results in reciprocal lattice points, and the number of q-points needed is which determine the size of the large cell. It has been verified that the implementation reproduces the FHNC/0 result of [17] in the homogeneous limit. In Fig. 1 the deviation of the inhomogeneous pair distribution function from the homogeneous on is clearly seen. Further, the minimum in at follows approximately a local density approximation.

Fig. 1

Pair distribution function of an inhomogeneous system is shown at different points g(0, r), g(1.5, r) and g(3, r) (solid lines, distinguished by the minimum). Also shown are the density (dashed), the non-interacting pair distribution function (dashed dotted) and the difference of the inhomogeneous and homogeneous result (short dashed)

The number of G vectors for a realistic 3D system is now estimated. Sodium is used as an example. The lattice constant of the primitive unit cell is , where for sodium. From the treatment of the homogeneous electron gas and the results for the 1D model system it is seen that points are sufficient for that length. It is even possible to reduce this number further. This results in a eigenvalue problem for each q-vector, which is numerically tractable. Only q-vectors in the irreducible Brillouin zone need to be calculated, which further reduces the numerical load.

Conclusions

The FHNC equations for periodic systems have been derived by simple Fourier transform of the general inhomogeneous equations. This considerably reduces the numerical demand so that realistic systems become numerically feasible. Further reduction in the numerical demand is possible, but only at the cost of simplicity of the theory.

First results for a 1D system have been shown. There is increasing interest in the theoretical description of 1D systems, e.g., [18-20]. The presented method could extend existing methods by applying it to more realistic systems. Also more fundamental questions could be addressed, like the extension of the local density approximation to pair quantities could be investigated. That is how to combine results of the homogeneous system for different densities to obtain an estimation of the inhomogeneous result. Work in that direction is in progress.