Where do the laws of quantum mechanics break down in simulations?

Question

As someone who holds a BA in physics I was somewhat scandalized when I began working with molecular simulations. It was a bit of a shock to discover that even the most detailed and computationally expensive simulations can't quantitatively reproduce the full behavior of water from first principles.

Previously, I had been under the impression that the basic laws of quantum mechanics were a solved problem (aside from gravity, which is usually assumed to be irrelevant at molecular scale). However, it seems that once you try to scale those laws up and apply them to anything larger or more complex than a hydrogen atom their predictive power begins to break down.

From a mathematics point of view, I understand that the wave functions quickly grow too complicated to solve and that approximations (such as Born-Oppenheimer) are required to make the wave functions more tractable. I also understand that those approximations introduce errors which propagate further and further as the time and spatial scales of the system under study increase.

What is the nature of the largest and most significant of these approximation errors? How can I gain an intuitive understanding of those errors? Most importantly, how can we move towards an ab-initio method that will allow us to accurately simulate whole molecules and populations of molecules? What are the biggest unsolved problems that are stopping people from developing these kinds of simulations?

Answer 1

As far as I'm aware, the most accurate methods for static calculations are Full Configuration Interaction with a fully relativistic four-component Dirac Hamiltonian and a "complete enough" basis set. I'm not an expert in this particular area, but from what I know of the method, solving it using a variational method (rather than a Monte-Carlo based method) scales shockingly badly, since I think the number of Slater determinants you have to include in your matrix scales something like $O(^{n_{orbs}}C_{n_e})$. (There's an article on the computational cost here.) The related Monte-Carlo methods and methods based off them using "walkers" and networks of determinants can give results more quickly, but as implied above, aren't variational. And are still hideously costly.

Approximations currently in practical use just for energies for more than two atoms include:

• Born Oppenheimer, as you say: this is almost never a problem unless your system involves hydrogen atoms tunneling, or unless you're very near a state crossing/avoided crossing. (See, for example, conical intersections.) Conceptually, there are non-adiabatic methods for the wavefunction/density, including CPMD, and there's also Path-Integral MD which can account for nuclear tunneling effects.
• Nonrelativistic calculations, and two-component approximations to the Dirac equation: you can get an exact two-component formulation of the Dirac equation, but more practically the Zeroth-Order Regular Approximation (see Lenthe et al, JChemPhys, 1993) or the Douglas-Kroll-Hess Hamiltonian (see Reiher, ComputMolSci, 2012) are commonly used, and often (probably usually) neglecting spin-orbit coupling.
• Basis sets and LCAO: basis sets aren't perfect, but you can always make them more complete.
• DFT functionals, which tend to attempt to provide a good enough attempt at the exchange and correlation without the computational cost of the more advanced methods below. (And which come in a few different levels of approximation. LDA is the entry-level one, GGA, metaGGA and including exact exchange go further than that, and including the RPA is still a pretty expensive and new-ish technique as far as I'm aware. There are also functionals which use differing techniques as a function of separation, and some which use vorticity which I think have application in magnetic or aromaticity studies.) (B3LYP, the functional some people love and some people love to hate, is a GGA including a percentage of exact exchange.)
• Configuration Interaction truncations: CIS, CISD, CISDT, CISD(T), CASSCF, RASSCF, etc. These are all approximations to CI which assume the most important excited determinants are the least excited ones.
• Multi-reference Configuration Interaction (truncations): Ditto, but with a few different starting reference states.
• Coupled-Cluster method: I don't pretend to properly understand how this works, but it obtains similar results to Configuration Interaction truncations with the benefit of size-consistency (i.e. $E(H_2) \times 2 = E((H_2)_2$ (at large separation)).

For dynamics, many of the approximations refer to things like the limited size of a tractable system, and practical timestep choice -- it's pretty standard stuff in the numerical time simulation field. There's also temperature maintenance (see Nose-Hoover or Langevin thermostats). This is mostly a set of statistical mechanics problems, though, as I understand it.

Anyway, if you're physics-minded, you can get a pretty good feel for what's neglected by looking at the formulations and papers about these methods: most commonly used methods will have at least one or two papers that aren't the original specification explaining their formulation and what it includes. Or you can just talk to people who use them. (People who study periodic systems with DFT are always muttering about what different functionals do and don't include and account for.) Very few of the methods have specific surprising omissions or failure modes. The most difficult problem appears to be proper treatment of electron correlation, and anything above the Hartree-Fock method, which doesn't account for it at all, is an attempt to include it.

As I understand it, getting to the accuracy of Full relativistic CI with complete basis sets is never going to be cheap without dramatically reinventing (or throwing away) the algorithms we currently use. (And for people saying that DFT is the solution to everything, I'm waiting for your pure density orbital-free formulations.)

There's also the issue that the more accurate you make your simulation by including more contributions and more complex formulations, the harder it is to actually do anything with. For example, spin orbit coupling is sometimes avoided solely because it makes everything more complicated to analyse (but sometimes also because it has negligable effect), and the canonical Hartree-Fock or Kohn-Sham orbitals can be pretty useful for understanding qualitative features of a system without layering on the additional output of more advanced methods.

(I hope some of this makes sense, it's probably a bit spotty. And I've probably missed someone's favourite approximation or niggle.)

Answer 2

The fundamental challenge of quantum mechanical calculations is that they do not scale very well—from what I recall, the current best-case scaling is approximately $O(N_e^{3.7})$, where $N_e$ is the number of electrons contained in the system. Thus, 13 water molecules will scale as having $N_e = 104$ electrons instead of just $N = 39$ atoms. (That's a factor of nearly 40.) For heavier atoms, the discrepancy becomes even greater.

The main issue will be that, in addition to increased computational horsepower, you will need to come up with better algorithms that can knock down the 3.7 exponent to something that is more manageable.

Answer 3

The problem is broadly equivalent to the difference between classical computers and quantum computers. Classical computers work on single values at once, as only one future/history is possible for one deterministic input. However, a quantum computer can operate on every possible input simultaneously, because it can be put in a superposition of all the possible states.

In the same way, a classical computer has to calculate every property individually, but the quantum system it is simulating has all the laws of the universe to calculate all the properties simultaneously.

The problem is exacerbated by the way we have to pass data almost serially through a CPU, or at most a few thousand CPUs. By contrast, the universe has a nearly unlimited set of simultaneous calculations going on at the same time.

Consider as an example 3 electrons in a box. A computer has to pick a timestep (first approximation), and keep recalculating the interactions of each electron with each other electron, via a limited number of CPUs. In reality, the electrons have an unknowable number of real and virtual exchange particles in transit, being absorbed and emitted, as a continuous process. Every particle and point in space has some interaction going on, which would need a computer to simulate.

Simulation is really the art of choosing your approximations and your algorithms to model the subject as well as possible with the resources you have available. If you want perfection, I'm afraid it's the mathematics of spherical chickens in vacuums; we can only perfectly simulate the very simple.

Answer 4

I don't know if the following helps, but for me it was very insightful to visualize the scaling behavior of quantum systems:

The main problem comes from the fact that the Hilbert space of quantum states grows exponentially with the number of particles. This can be seen very easily in discrete systems. Think of a couple of potential wells that are conected to each other, may just two: well 1 and well 2. Now add bosons (e.g., Rubidium 87, just as an example), at first only one. How many possible basis vectors are there?

• basis vector 1: boson in well 1
• basis vector 2: boson in well 2

They can be written like $\left|1,0 \right\rangle$ and $\left|0,1 \right\rangle$

Now suppose the boson can hop (or tunnel) from one well to the other. The Hamiltonian that describes the system can then be written is matrix notation as

$$ \hat{H}=\pmatrix{ \epsilon_{1} & t \\ t & \epsilon_{2}} $$

where $\epsilon_{1,2}$ are just the energies of the boson in well 1 and 2, respectively, and t the tunneling amplitude. The complete solution of this system, i.e., the solution containing all the information necessary to compute the system's state at any given point of time (given an initial condition), is given by the eigenstates and eigenvalues. The eigenstates are linear superpositions of the basis vectors (in this case $\left|1,0 \right\rangle$ and $\left|0,1 \right\rangle$).

This problem is so simple that it can be solved by hand.

Now suppose we have more potential wells and more bosons, e.g., in the case of four wells with two bosons there are 10 different possibilities to distribute the bosons among the wells. Then the Hamiltonian would have 10x10=100 elements and 10 eigenstates.

One can quickly see that the number of eigenstates is given by the binomial coefficient: $$ \text{number of eigenstates}=\pmatrix{\text{number of wells} + \text{number of bosons} - 1 \\ \text{number of bosons}} $$

So even for "just" ten bosons and ten different potential wells (a very small system), we'd have 92,378 eigenstates. The size of the Hamiltonian is then $92,378^2$ (approximately 8.5 billion elements). In a computer they'd occupy (depending on your system) about 70 gigabytes of RAM and is therefore probably impossible to solve on most computers.

Now let's assume we have a continuous system (i.e. no potential wells, but free space) and 13 water molecules (for simplicity I treat them each of them as a particle). Now in computer we can still model free space using many tiny potential wells (we discretize space... which is ok, as long as the relevant physics takes place on larger length scales then the discretization length). Let's say there are 100 different possible positions for each of the molecules in each of the x, y and z directions. So we end up with 100*100*100 = 1,000,000 little boxes. Then we'd have more than $2.7 \cdot 10^{53}$ basis vectors, the Hamiltonian would have almost $10^{107}$ elements, occupying so much space that we'd need all the particles from 10 million universes like ours just to encode that information.

Answer 5

One problem is that quantum mechanics suffers from the "curse of dimensionality": for most methods of solving PDEs, the number of basis functions needed to get a certain accuracy scales exponentially with the number of dimensions. Since an $n$-electron atom has $3n$ degrees of freedom, even relatively small systems require an enormous number of dimensions to simulate exactly.

Monte Carlo can be used to get around this problem, as the error scales like $\text{points}^{-\frac{1}{2}}$ regardless of the number of dimensions, but convergence is slow.

Density functional theory is another way to deal with this problem, but it's an approximation. It's a very good approximation in some cases, but in other cases it can be surprisingly bad.

Answer 6

I think a highly-accurate simulation of water was the topic of one of the very first and large simulations performed using the Jaguar supercomputer. You might want to look into this paper and their follow-up work (which, by the way, was a finalist for the Gordon-Bell prize in 2009):

"Liquid water: obtaining the right answer for the right reasons", Aprà, Rendell, Harrison, Tipparaju, deJong, Xantheas.

Answer 7

This problem is solved by Density Functinal Theory. The essence is replacing many body degrees of freedom by several fields one of them beeing the density of electrons. For a grand exposition see the nobel lecture of one of the founders of DFT: http://www.nobelprize.org/nobel_prizes/chemistry/laureates/1998/kohn-lecture.pdf