Ab Initio Molecular Dynamics

Ab Initio Molecular Dynamics

(Parte 6 de 7)

2.8.4 Wavelets

Similar to using generalized plane waves is the idea to exploit the powerful multiscale–properties of wavelets. Since this approach requires an extensive introductory discussion (see e.g. Ref. 242 for a gentle introduction) and since it seems still quite far from being used in large–scale electronic structure calculations the interested reader is referred to original papers 134,674,699,652,241,25 and the general wavelet literature cited therein. Wavelet–based methods allow intrinsically to exploit multiple length scales without introducing Pulay forces and can be efficiently handled by fast wavelet transforms. In addition, they are also a powerful route to linear scaling or “order–N” methods 453,243 as first demonstrated in Ref. 241 with the calculation of the Hartree potential for an all–electron uranium dimer.

2.8.5 Mixed and Augmented Basis Sets

Localized Gaussian basis functions on the one hand and plane waves on the other hand are certainly two extreme cases. There has been a tremendous effort to combine such localized and originless basis functions in order to exploit their mutual strengths. This resulted in a rich collection of mixed and augmented basis sets with very specific implementation requirements. This topic will not be covered here and the interested reader is referred to Refs. 75,654,498,370,371 and references given therein for some recent implementations used in conjunction with ab initio molecular dynamics.

2.8.6 Wannier Functions

An alternative to the plane wave basis set in the framework of periodic calculations in solid–state theory are Wannier functions, see for instance Sect. 10 in Ref. 27. These functions are formally obtained from a unitary transformation of the Bloch orbitals Eq. (114) and have the advantage that they can be exponentially localized under certain circumstances. The so–called maximally localized generalized Wannier functions 413 are the periodic analogues of Boys’ localized orbitals defined for isolated systems. Recently the usefulness of Wannier functions for numerical purposes was advocated by several groups, see Refs. 339,184,413,10 and references given therein.

2.8.7 Real Space Grids

A quite different approach is to leave conventional basis set approaches altogether and to resort to real–space methods where continuous space is replaced by a discrete space r → rp. This entails that the derivative operator or the entire energy expression has to be discretized in some way. The high–order central–finite difference approach leads to the expression n =−N Cn ψi(rp

+ nxh,rp , rp

n =−N Cn ψi(rp

, rp + nyh,rp

n =−N Cn ψi(rp for the Laplacian which is correct up to the order h2N+2. Here, h is the uniform grid spacing and {Cn} are known expansion coefficients that depend on the selected order 130. Within this scheme, not only the grid spacing h but also the order are disposable parameters that can be optimized for a particular calculation. Note that the discretization points in continuous space can also be considered to constitute a sort of “finite basis set” – despite different statements in the literature – and that the “infinite basis set limit” is reached as h → 0 for N fixed. A variation on the theme are Mehrstellen schemes where the discretization of the entire differential equation and not only of the derivative operator is optimized 89.

The first real–space approach devised for ab initio molecular dynamics was based on the lowest–order finite–difference approximation in conjunction with a equally– spaced cubic mesh in real space 109. A variety of other implementations of more sophisticated real–space methods followed and include e.g. non–uniform meshes, multigrid acceleration, different discretization techniques, and finite–element methods 686,61,39,130,131,632,633,431,634. Among the chief advantages of the real–space methods is that linear scaling approaches 453,243 can be implemented in a natural way and that the multiple–length scale problem can be coped with by adapting the grid. However, the extension to such non–uniform meshes induces the (in)famous Pulay forces (see Sect. 2.5) if the mesh moves as the nuclei move.

3 Basic Techniques: Implementation within the CPMD Code

3.1 Introduction and Basic Definitions

This section discusses the implementation of the plane wave–pseudopotential molecular dynamics method within the CPMD computer code 142. It concentrates on the basics leaving advanced methods to later chapters. In addition all formulas are for the non-spin polarized case. This allows to show the essential features of a plane wave code as well as the reasons for its high performance in detail. The implementation of other versions of the presented algorithms and of the more advanced techniques in Sect. 4 is in most cases very similar.

There are many reviews on the pseudopotential plane wave method alone or in connection with the Car–Parrinello algorithm. Older articles 312,157,487,591 as well as the book by Singh 578 concentrate on the electronic structure part. Other reviews 513,472,223,224 present the plane wave method in connection with the molecular dynamics technique.

3.1.1 Unit Cell and Plane Wave Basis

The unit cell of a periodically repeated system is defined by the Bravais lattice vectors a1, a2, and a3. The Bravais vectors can be combined into a three by three

Further, scaled coordinates s are introduced that are related to r via h r = hs . (105)

Distances in scaled coordinates are related to distances in real coordinates by the metric tensor G = hth

where [·]NINT denotes the nearest integer value. The coordinates rpbc will be always within the box centered around the origin of the coordinate system. Recip- rocal lattice vectors bi are defined as bi · aj = 2pi δij (108) and can also be arranged to a three by three matrix

Plane waves build a complete and orthonormal basis with the above periodicity (see also the section on plane waves in Sect. 2.8)

with the reciprocal space vectors where g = [i,j,k] is a triple of integer values. A periodic function can be expanded in this basis

where ψ(r) and ψ(G) are related by a three-dimensional Fourier transform. The direct lattice vectors L connect equivalent points in different cells.

3.1.2 Plane Wave Expansions

The Kohn–Sham potential (see Eq. (82)) of a periodic system exhibits the same periodicity as the direct lattice

where k is a vector in the first Brillouin zone. The functions ui(r,k) have the periodicity of the direct lattice

The index i runs over all states and the states have an occupation fi(k) associated with them. The periodic functions ui(r,k) are now expanded in the plane wave basis

and the Kohn–Sham orbitals are

where ci(G,k) are complex numbers. With this expansion the density can also be expanded into a plane wave basis

where the sum over G vectors in Eq. (119) expands over double the range given by the wavefunction expansion. This is one of the main advantages of the plane wave basis. Whereas for atomic orbital basis sets the number of functions needed to describe the density grows quadratically with the size of the system, there is only a linear dependence for plane waves.

3.1.3 K–Points and Cutoffs

In actual calculations the infinite sums over G vectors and cells has to be truncated. Furthermore, we have to approximate the integral over the Brillouin zone by a finite sum over special k–points ∫

where wk are the weights of the integration points. Schemes on how to choose the integration points efficiently are available in the literature 30,123,435 where also an overview 179 on the use of k–points in the calculation of the electronic structure of solids can be found.

The truncation of the plane wave basis rests on the fact that the Kohn–Sham potential V KS(G) converges rapidly with increasing modulus of G. For this reason, at each k–point, only G vectors with a kinetic energy lower than a given maximum cutoff

are included in the basis. With this choice of the basis the precision of the calculation within the approximations of density functional theory is controlled by one parameter Ecut only.

The number of plane waves for a given cutoff depends on the unit cell and the k–point. An estimate for the size of the basis at the center of the Brillouin zone is

where Ecut is in Hartree units. The basis set needed to describe the density calculated from the Kohn-Sham orbitals has a corresponding cutoff that is four times the cutoff of the orbitals. The number of plane waves needed at a given density cutoff is therefore eight times the number of plane waves needed for the orbitals.

It is a common approximation in density functional theory calculations 536,169 to use approximate electronic densities. Instead of using the full description, the density is expanded in an auxiliary basis. An incomplete plane wave basis can be considered as an auxiliary basis with special properties 371. Because of the filter property of plane waves the new density is an optimal approximation to the true density. No additional difficulties in calculations of the energy or forces appear. The only point to control is, if the accuracy of the calculation is still sufficient.

Finally, sums over all unit cells in real space have to be truncated. The only term in the final energy expression with such a sum is the real space part of the Ewald sum (see Sect. 3.2). This term is not a major contribution to the workload in a density functional calculation, that is the cutoff can be set rather generously.

3.1.4 Real Space Grid

A function given as a finite linear combination of plane waves can also be defined as a set of functional values on a equally spaced grid in real space. The sampling theorem (see e.g. Ref. 492) gives the maximal grid spacing that still allows to hold the same information as the expansion coefficients of the plane waves. The real space sampling points R are defined where N is a diagonal matrix with the entries 1/Ns and q is a vector of integers ranging from 0 to Ns − 1 (s = x, y, z). To fulfill the sampling theorem Ns has to be bigger than 2max(gs) + 1. To be able to use fast Fourier techniques, Ns must be decomposable into small prime numbers (typically 2, 3, and 5). In applications the smallest number Ns that fulfills the above requirements is chosen. A periodic function can be calculated at the real space grid points

Nx igxqx]exp

Ny igyqy]exp

extended over all indices in the cube −gmaxsgmaxs . The functions f(R) and

The function f(G) is zero outside the cutoff region and the sum over g can be 46 f(G) are related by three–dimensional Fourier transforms

The Fourier transforms are defined by l=0 fG jklexp[ i




(Parte 6 de 7)