Introduction

In this chapter we consider the motion of nuclei in the classical limit. The laws of classical mechanics apply, and the nuclei move in a conservative potential, determined by the electrons in the Born–Oppenheimer approximation. The electrons are assumed to be in the ground state, and the energy is the ground-state solution of the time-independent Schrödinger equation, with all nuclear positions as parameters. This is similar to the assumptions of Car–Parrinello dynamics (see Section 6.3.1), but the derivation of the potential on the fly by quantum-mechanical methods is far too compute-intensive to be useful in general. In order to be able to treat large systems over reasonable time spans, a simple description of the potential energy surface is required to enable the simulation of motion on that surface. This is the first task: design a suitable force field from which the forces on each atom, as well as the total energy, can be efficiently computed, given the positions of each atom. Section 6.3 describes the principles behind force fields, and emphasizes the difficulties and insufficiencies of simple force field descriptions. But before considering force fields, we must define the system with its boundary conditions (Section 6.2). The way the interactions over the boundary of the simulated systems are treated is in fact part of the total potential energy description.

In this chapter a number of models for the thermodynamic properties of various phases will be described. The integral Gibbs energy will be used as the modeled thermodynamic property. The reason to model the Gibbs energy rather than any other thermodynamic function is that most experiments are done at constant temperature and pressure. From the Gibbs energy all other important quantities can be obtained according to Eqs. (2.12). Using the Gibbs energy means that the modeling is limited to a “mean-field” approximation. Thermodynamic calculations using “Monte Carlo” methods or “molecular dynamics” are outside the scope of this presentation, but these techniques can provide important information about the type of mean-field model to be selected.

A phase may sometimes have a particular physical or chemical feature that requires a special model in order for it to be described accurately. It is not uncommon that the mathematical expression for such a model may be identical to the expression derived to describe another physical feature. That simply means that the mathematical expression is more general than the physical model. Whenever such a generalized expression can be obtained, it will be called a formalism. A general formalism should be able to handle cases when various constituents added to a phase behave differently, for example some may dissolve interstitially or cause chemical ordering. Most of the models used in this book are special cases of the compound-energy formalism (CEF).

Introduction

Classical mechanics is not only an approximation of quantum mechanics, valid for heavy particles, but historically it also forms the basis on which quantum-mechanical notions are formed. We also need to be able to describe mechanics in generalized coordinates if we wish to treat constraints or introduce other ways to reduce the number of degrees of freedom. The basis for this is Lagrangian mechanics, from which the Hamiltonian description is derived. The latter is not only used in the Schrödinger equation, but forms also the framework in which (classical) statistical mechanics is developed. A background in Lagrangian and Hamiltonian mechanics is therefore required for many subjects treated in this book.

After the derivation of Lagrangian and Hamiltonian dynamics, we shall consider how constraints can be built in. The common type of constraint is a holonomic constraint that depends only on coordinates, such as a bond length constraint, or constraints between particles that make them behave as one rigid body. An example of a non-holonomic constraint is the total kinetic energy (to be kept at a constant value or at a prescribed time-dependent value).

We shall only give a concise review; for details the reader is referred to text books on classical mechanics, in particular to Landau and Lifschitz (1982) and to Goldstein et al. (2002).

There are several ways to introduce the principles of mechanics, leading to Newton's laws that express the equations of motion of a mechanical system.

Introduction

In the previous chapter we considered systems of interacting particles. They were treated as classical particles for which the interaction potential is known. We had to solve the classical equations of motion to simulate the behaviour of such a system at some nonzero temperature. Had we added frictional forces, the system would have evolved towards the ground state. In this chapter we discuss methods for simulating interacting atoms and molecules using quantum mechanical calculations. In fact, we consider the nuclei on a classical level but use quantum mechanics for the electronic degrees of freedom. Again, we can use this approach either to simulate a system of interacting particles at a finite temperature, or to find the ground state (minimum energy) configurations of solids and of molecules.

In Chapters 4 to 6 we studied methods for solving the electronic structure of molecular and solid state systems with a static configuration of nuclei (Born–Oppenheimer approximation). Knowledge of the electronic structure includes knowledge of the total energy. Therefore, by varying the positions of the nuclei, we can study the dependence of the total energy on these positions. The energy E(R1, R2, …, RN) as a function of the nuclear positions Ri is called the potential surface. As a simple example, consider the hydrogen molecule. We assume that the molecule is not rotating, so that the nuclear motion is a vibration along the molecular axis.

This is a book on computational methods used in theoretical physics research, with an emphasis on condensed matter applications.

Computational physics is concerned with performing computer calculations and simulations for solving physical problems. Although computer memory and processor performance have increased dramatically over the last two decades, most physical problems are too complicated to be solved without approximations to the physics, quite apart from the approximations inherent in any numerical method. Therefore, most calculations done in computational physics involve some degree of approximation. In this book, emphasis is on the derivation of algorithms and the implementation of these: it is a book which tells you how methods work, why they work, and what the approximations are. It does not contain extensive discussions on results obtained for large classes of different physical systems.

This book is not elementary: the reader should have a background in basic undergraduate physics and in computing. Some background in numerical analysis is also helpful. On the other hand, the topics discussed are not treated in a comprehensive way; rather, this book hopefully bridges the gap between more elementary texts by Koonin, Gould and Giordano, and specialised monographs and review papers on the applications described. The fact that a wide range of topics is included has the advantage that the many similarities in the methods used in seemingly very different fields could be highlighted. Many important topics and applications are however not considered in this book – the material presented obviously reflects my own expertise and interest.

Introduction

Here and in the following chapter we treat two different approaches to the many-electron problem: the Hartree–Fock theory and the density functional theory. Both theories are simplifications of the full problem of many electrons moving in a potential field. In fact, the physical systems we want to study, such as atoms, molecules and solids, consist not only of electrons but also of nuclei, and each of these particles moves in the field generated by the others.Afirst approximation is to consider the nuclei as being fixed, and to solve the Schrödinger equation for the electronic system in the field of the static nuclei. This approach, called the Born–Oppenheimer approximation, is justified by the nuclei being much heavier than the electrons so that they move at much slower speeds. It remains then to solve for the electronic structure.

The Hartree–Fock methodcan be viewed as a variational method in which the wave functions of the many-electron system have the form of an antisymmetrised product of one-electron wave functions (the antisymmetrisation is necessary because of the fermion character of the electrons). This restriction leads to an effective Schrödinger equation for the individual one-electron wave functions (called orbitals) with a potential determined by the orbitals occupied by the other electrons. This coupling between the orbitals via the potentials causes the resulting equations to become nonlinear in the orbitals, and the solution must be found iteratively in a self-consistency procedure.

In the previous chapter we encountered density functional theory (DFT) which is extensively used for calculating the electronic structure of periodic solids. Aside from DFT, carefully designed potentials often allow accurate electronic structures to be obtained by simply solving the Schrödinger equation without going through the self-consistency machinery of DFT. In both approaches it is necessary to solve the Schrödinger equation and the present chapter focuses on this problem, although some comments on implementing a DFT self-consistency loop will be made.

The large number of electrons contained in a macroscopic crystal prohibits a direct solution of the Schrödinger equation for such a system. Fortunately, the solid has periodic symmetry in the bulk, and this can be exploited to reduce the size of the problem significantly, using Bloch's theorem, which enables us to replace the problem of solving the Schrödinger equation for an infinite periodic solid by that of solving the Schrödinger equation in a unit cell with a series of different boundary conditions – the so-called Bloch boundary conditions. Having done this, there remains the problem that close to the nuclei the potential diverges, whereas it is weak when we are not too close to any of the nuclei (interstitial region).We can take advantage of the fact that the potential is approximately spherically symmetric close to the nuclei, but further away the periodicity of the crystal becomes noticeable. These two different symmetries render the solution of the Schrödinger equation in periodic solids difficult.

Introduction

It is not necessary to recall the dramatic increase in computer speed and the drop in cost of hardware over the last two decades. Today, anyone can buy a computer with which all of the programs in this book can be executed within a reasonable time – typically a few seconds to a few hours.

On the other hand, if there is one conclusion to be drawn from the enormous amount of research in computational physics, it should be that for most physical problems, a realistic treatment, one without severe approximations, is still not within reach. Quantum many-particle problems, for example, can only be treated if the correlations are treated in an approximate way (this does not hold for quantum Monte Carlo techniques, but there we suffer from minus-sign problems when treating fermions; see Chapter 12). It is easy to extend this list of examples.

Therefore the physical community always follows the developments in hardware and software with great interest. Developments in this area are so fast that if a particular type of machine were presented here as being today's state of the art, this statement would be outdated by the time the book is on the shelf. We therefore restrict ourselves here to a short account of some general principles of computer architecture and implications for software technology. The two main principles are pipelining and parallelism. Both concepts were developed a few decades ago, but pipelining became widespread in supercomputers from around 1980 and has found its way into most workstations, whereas parallelism has remained more restricted to the research community and to more expensive machines.

About numerical methods

In computational physics, many mathematical operations have to be carried out numerically. Techniques for doing this are studied in numerical analysis. For a large variety of numerical problems, commercial or public domain software packages are available, and on the whole these are recommended in preference to developing all the numerical software oneself, not only because this saves time but also because available routines are often coded by professionals and hence are hard to surpass in generality and efficiency. To avoid eventual problems with using existing numerical software, it is useful to have some background in numerical analysis. Moreover, knowledge of this field enables you to write codes for special problems for which no routine is available. This chapter reviews some numerical techniques which are often used in computational physics.

There are several ways of obtaining ‘canned’ routines. Commercially available libraries (such as NAG, IMSL) are of very high quality, but often the source code is not available, which might prevent your software from being portable. However, several libraries, such as the ones quoted, are available at many institutes and companies for various types of computers (‘platforms’), so that in practice this restriction is not so severe.

Via the internet, it is possible to obtain a wide variety of public domain routines; a particularly useful site is /www.netlib.org/. Most often these are provided in source code. Another cheap way of obtaining routines is by purchasing a book on numerical algorithms containing listings of source codes and, preferably, a CD (or an internet address) with these sources. A well-known book is Numerical Recipes by Press et al. [1].

Physics and computational physics

Solving a physical problem often amounts to solving an ordinary or partial differential equation. This is the case in classical mechanics, electrodynamics, quantum mechanics, fluid dynamics and so on. In statistical physics we must calculate sums or integrals over large numbers of degrees of freedom. Whatever type of problem we attack, it is very rare that analytical solutions are possible. In most cases we therefore resort to numerical calculations to obtain useful results. Computer performance has increased dramatically over the last few decades (see also Chapter 16) and we can solve complicated equations and evaluate large integrals in a reasonable amount of time.

Often we can apply numerical routines (found in software libraries for example) directly to the physical equations and obtain a solution. We shall see, however, that although computers have become very powerful, they are still unable to provide a solution to most problems without approximations to the physical equations. In this book, we shall focus on these approximations: that is, we shall concentrate on the development of computational methods (and also on their implementation into computer programs). In this introductory chapter we give a bird's-eye perspective of different fields of physics and the computational methods used to solve problems in these areas. We give examples of direct application of numerical methods but we also give brief and heuristic descriptions of the additional theoretical analysis and approximations necessary to obtain workable methods for more complicated problems which are described in more detail in the remainder of this book.

Introduction

Classical field theory enables us to calculate the behaviour of fields within the framework of classical mechanics. Examples of fields are elastic strings and sheets, and the electromagnetic field. Quantum field theory is an extension of ordinary quantum mechanics which not only describes extended media such as string and sheets, but which is also supposed to describe elementary particles. Furthermore, ordinary quantum many-particle systems in the grand canonical ensemble can be formulated as quantum field theories. Finally, classical statistical mechanics can be considered as a field theory, in particular when the classical statistical model is formulated on a lattice, such as the Ising model on a square lattice, discussed in Chapter 7.

In this chapter we shall describe various computational techniques that are used to extract numerical data from field theories. Renormalisation is a procedure without which field theories cannot be formulated consistently in continuous space-time. In computational physics, we formulate field theories usually on a lattice, thereby avoiding the problems inherent to a continuum formulation. Nevertheless, understanding the renormalisation concept is essential in lattice field theories in order to make the link to the real world. In particular, we want to make predictions about physical quantities (particle masses, interaction constants) which are independent of the lattice structure, and this is precisely where we need the renormalisation concept.

Quantum field theory is difficult. It does not belong to the standard repertoire of every physicist.

Introduction

In the previous chapter we saw that the experimental values of physical quantities of a many-particle system can be found as an ensemble average. Experimental systems are so large that it is impossible to determine this ensemble average by summing over all the accessible states in a computer. There exist essentially two methods for determining these physical quantities as statistical averages over a restricted set of states: the molecular dynamics and Monte Carlo methods. Imagine that we have a random sample of, say, 107 configurations of the system which are all compatible with the values of the system parameters. For such a large number we expect averages of physical quantities over the sample to be rather close to the ensemble average. It is unfortunately impossible to generate such a random sample; however, we can generate a sample consisting of a large number of configurations which are determined successively from each other and are hence correlated. This is done in the molecular dynamics and Monte Carlo methods. The latter will be described in Chapter 10.

Molecular dynamics is a widely used method for studying classical many-particle systems. It consists essentially of integrating the equations of motion of the system numerically. It can therefore be viewed as a simulation of the system as it develops over a period of time. The system moves in phase space along its physical trajectory as determined by the equations of motion, whereas in the Monte Carlo method it follows a (directed) random walk.

Introduction

In Chapters 8 and 10 we studied methods for simulating classical systems consisting of many interacting degrees of freedom. In these methods, a sequence of system configurations is generated, and from this sequence averages of physical quantities, given as functions of the degrees of freedom (positions and momenta, or spins), can be determined. These quantities are called mechanical quantities, and the expressions of their expectation values are called mechanical averages.

There exist, however, quantities that cannot be determined straightforwardly using these methods. These quantities include free energies and chemical potentials. The point is that these quantities are not given as a normalised average, which for mechanical quantities is replaced by an average over the configurations generated in the simulation. In the previous chapters we have seen that it is not straightforward to find free energies and chemical potentials using MC and MD methods (see Section 10.5).

In this chapter we discuss a method which enables us to find free energies for lattice spin models with very high accuracy; this is the transfer matrix method. This method calculates the free energy of a model defined on a strip of finite width and infinite length directly in terms of the largest eigenvalue of a large matrix, the transfer matrix. This matrix contains essentially the Boltzmann weights for adding an extra row of spins to the strip. In Chapter 12 we shall see that the transfer matrix is the analogue of the time evolution operator in quantum mechanics.