Basic concepts for plane wave pseudopotential calculation.

 

In a plane wave pseudopotential calculation, the single electron wavefunction y® is expanded using plane wave basis exp(iG*r):

 

here are the coefficients of the wavefunctions.  are the coefficients kept by the computer to represent the wavefunction, instead of the real space values. The selection of plane wave G is done using a sphere of Gc1, which is defined by Ecut1 as: , and the candidate G’s are the reciprocal vectors of the unit cell AL(3,3). As a result, the wavefunction y® is periodic to AL. Thus, a planewave  calculation is always for a periodic system.

 

The biggest advantage of planewave basis is its ability to perform the exact variational calculation based on a discrete numerical grid. Basically,

 

,

 

here, n1*n2*n3 are the real space grid points indexed by l, W is the unit cell volume. The above equation means that instead of doing a actual real space integrations, we can do a summation over a real space grid, the results are exactly the same. So, our numerical calculation based on the finite real space grid of n1*n2*n3 is an exact calculation, restricted only by the variational degree of freedom of the finite planewave basis for

the wavefunction y. This is a nice analytical property not shared by many other numerical methods.

 

Notice that,  even for , instead we have:

 

 

here, V(G) are the original Fourier components of V(r) and. Thus, the Fourier components beyond Gc2 do not have any contribution here.  In order for

the above integral=summation equation to be true, we need , and the

sphere of Gc2 being inside the Fourier box of n1*n2*n3 grid (i.e, the n1,n2,n3 should be large enough). This is represented graphically in the following. In summary, Ecut1 determines Ecut2 (usually by a factor of 4), and Ecut2 determines n1,n2,n3.

 

 

 

Ecut1,Ecut2, n1, n2, n3 are used as parameters in the PEtot calculation. Notice that, to save computation (in the expense of possible error),  sometime a value less than 4*Ecut1 is used for Ecut2 (e.g, 3*Ecut1, or even 2*Ecut1). It is Ecut2 which determines the n1,n2,n3. The current PEtot program requires that Gc2 being inside the Fourier box of n1,n2,n3.

 

Another advantage of planewave calculation is that the gradient operator is diagonal in G-space:

 

 

 As a result, we carry out different operations in different spaces (real-space, and G-space). This dual-space representation requires us to do frequent Fourier transform between them. This is done using the numerical Fast Fourier Transformation (FFT).

It is due to the FFT algorithm, the planewave calculation becomes attractive. Still, the FFT part is often the most time consuming part of a plane wave calculation.

 

 

The plane wave calculation requires that the wavefunction can be described byo the planewaves within Ecut1. In order to reduce Ecut1 to a reasonable value (thus reduce n1,n2,n3), the wavefunction should be smooth. While it is often smooth at the chemically important bonding area, near the nuclei, the valence wavefunctions often have a lot of wiggles. These wiggles must exist so that the valence wavefunctions are orthogonal to the deep level core wavefunctions. Besides, it is difficult to describe those chemically-not-important core wavefunctions using planewave basis. To overcome these difficulties, pseudopotentials are developed. Basically, using a pseudopotential, the core states will no longer exist, and the valence pseudo-wavefunctions becomes smooth near the nuclei. The pseudo-wavefunctions become the same as the original wavefunction for , and the eigen energies are not changed.  However, in order to have these good properties, it is necessary to have different pseudopotentials for s,p,d.. states, i.e, the pseudopotential is angular momentum dependent (non-local pseudopotential). This is illustrated in the following figure.

 

 

Mathematically, the angular dependent pseudopotential can be written as:

 

here  is an angular momentum projection operator. In computation, this can be evaluated using a full matrix: , analytical formula exists to write down the matrix elements. This is done in the old days for small system calculations. For large system calculations, the Hamiltonian matrix elements <G1|H|G2> are never explicitly written down. So, this because impractical. In such calculation (as in PEtot), the Kleinman-Bylander implementation of the nonlocal potential is used. Basically, we first take one l (s, or, p, or d) as the local potential, , then we can define a nonlocal part as . Note that this nonlocal part is error outside . Then the above nonlocal pseudopotential is approximated as:

 

 

here m is the azimuthal angular momentum,  is the solid angle spherical harmonic

function, and  are the atomic pseudowavefunctions for angular l. This Kleinman-Bylander formula can be carried out either in G-space (represent the projection function  in G-space), or real-space. In the real-space representation, the projection function can be cut off beyond a radius rcut (which is used as a parameter in PEtot).

As a result, for large system calculations, real-space implementation of the Kleinman-Bylander form is faster than G-space implementation.

 

 

For ultrasoft pseudopotential (D. Vanderbilt, Phys. Rev. B 41, 7892 (1990)), we have the following concepts. The total charge density is:

 

 

where the soft charge density is

 

 

In above, Q_mn^I(r) is a fixed form of additional charge density around each atom I. It is used to compensate the missed charge density by the ultrasoft wavefunction near each nuclei. The beta_n is the reference (projection) function for the nonlocal potential. The index m,n are the index for the reference functions.  Each (l,m) (angular momentum) might have two reference functions. Note the nonlocal pseudopotential operator in an ultrasoft pseudopotential Hamiltonian is:

 

in the total energy expression. While in the Single particle Kohn-Sham equation, it is:

 

 

where D_nm^I should be calculated selfconsistently as:

 

 

and the LDA total local potential V_L(r) is calculated using the LDA formula from the total charge density rho_L(r).  Finally, the wavefunctions satisfy the generalized orthonormality condition:

 

where .

 

Note, to solve the single particle Kohn-Sham equation, the real space wavefunction is expressed in a small grid: n1*n2*n3. The charge density rho(r), and potential V(r) are all expressed in this grid. However, since Q_nm^I(r) is sharp, we need to express rho_L(r) and calculate V_L(r) is a larger grid: n1L*n2L*n3L (usually n1L=2*n1, etc). After V_L(r) is calculated in the larger grid, V(r) is generated from V_L(r) using Fourier Space truncation technique (via 2 FFT). Note, Ecut2L corresponds to the n1L*n2L*n3L grid. For all these grids, the real space box is AL(3,3).