This first (well, zeroth) entry will be with no code, simply because it’s the first one and deserves some introduction. The real first post is up along with this one, so if you don’t feel like reading my rambling about how the establishment is doing it all wrong, jump to the second section. If you don’t want to read about profound philosophy of representation of function spaces in computers either (really?), then just skip this part all-together and go straight to next part.
What’s all this about?
The inspiration to this series comes from a course on topology in condensed matter, which I really like because it does not only talk about the physics, but actually does it and shows the results. It doesn’t just throw formulas at you, but uses them to show how the physical systems behave, at least on the theoretical level. But it’s not perfect. It focuses on, unsurprisingly, topology in condensed matter physics, which might not be something that interests everyone. And even if it is something that some of you would be interested in, the required amount of previous knowledge in condensed matter or quantum mechanics might be too high of a potential barrier to overcome. The second issue is that it obfuscates a large part of the computational algorithms and you have to do a lot of digging into the source code to find out how it works under the hood.
You could argue that there is a lot blogs, code repos and books out there to solve the inadequacies of an online topology course, and you would probably be right. If you know what you’re doing and know where to look or have the correct people around you to ask, you can build yourself from scratch to the point where you can comfortably do computational physics at or higher than the level of the course. But most people don’t fall into this category. The people who don’t know where and what to look for or who to ask (i.e. not higher undergrad or grad students working on computational topics) will go on the Internet and get flooded by hundreds of “I programmed a double-pendulum in Python” types of resources or resources that are too advanced for someone who doesn’t (plan to) do this for a living. It’s unfortunate that there’s this huge empty chasm between the absolute basics and science-grade code, which makes it hard to do anything cool.
This resource gap is where I want to find a niche and do something fun for me and useful for others. In this series, I want to show how to do computational physics that’s more than just drawing the trajectory of Earth around the Sun and slowly build up to some real, modern physics with no shortcuts.
So, what will you find in this series?
- Me bumbling around. I’m an experimental condensed matter physicist by trade but, back in college and early grad school, I was also working in the computational side of things. So you should expect things to mostly work, but they will not always be as efficient and the code might sometimes look ugly.
- A lot of computational quantum mechanics and statistical physics and maybe a tiny bit of particle physics. Things like simple astronomy and mechanics are covered by many other resources out there and I simply don’t know enough physics related to the more advanced stuff from these fields. Sorry, no computational general relativity or fluid dynamics (probably).
- Code that just works. Every single line of code will be shown and every function will be discussed. I don’t like hiding the computation in computational physics and I want you to see the process of going from textbook physics to a finished fancy figure that you could see in a published paper. Being able to see a worked example is a good way to overcome the learning curve and going from reading about exciting physics to doing it.
And what will you not find in my posts?
- A textbook on basic physics. I know I said that I don’t want this to become a collection of advanced physics topics, but there’s a limit to how much I can fit into this format. Through the text, I will do my best to reference freely available resources wherever possible, but don’t expect to understand everything I do here without knowledge of the basics. That doesn’t mean you won’t get anything out of this if you’re not a physics student or a full-grown physicist. You can always reach out to Physics Stack Exchange or Reddit AskPhysics where me and many others will help you if you have questions or just read the material that I link and you should be able to follow what’s going on.
- The same boring code and programming paradigms as everywhere else. If you were like me, your computational physics and numerics classes were taught in C, C++, Java or Python. Or Matlab (my condolences). We all had it hammered into our heads that object oriented programming is the One True Way and when you become a big boy scientist, C++ and Fortran will be the only acceptable languages. You will see that that’s not necessarily true anymore and I have to be the change I want to see in the world of scientific computing. If you still want to do this the conservative way, you easily can, because I will show all of the code and you will always have something to base your programs on.
So, who is this series for? Well, firstly me, I guess. I do this in my off time as a hobby. But if it were just that, I could just as easily keep all the code in private repo and not bother with all the text and derivations. It’s also for other people to read and learn. If you’re a physics student that wants to learn how to do computational physics, then you’re the target audience of this series. Similarly, if you’re a grown up physicist or material scientist that wants to see how the sausage is made, you should find this interesting. If you come from outside of physics and want to see and learn cool things, you absolutely can. Do expect, however, that some parts of my explanations might not be enough for you. Or, even worse, they might feel like they make sense to you but they actually don’t. Go to the sites I linked and ask if you don’t understand something. I’m not going to pretend that there are no stupid questions, but it’s always preferable to ask them rather just keeping a stupid question to yourself and not having an answer. We all were new to this at a certain point and no scientist will ever get mad at you or make fun of you not knowing something if you show genuine interest in a topic.
No snakes in this codebase
It shouldn’t surprise anyone that there will be a lot of code in a computational physics blog. And because there’s a lot of it, it’s worth discussing, at least for a little bit, what programming language I’ll be using and what are my preferred alternatives (and why Python isn’t any of them).
As I already mentioned, you will not find any of the workhorse languages like C++ or Fortran here. I sometimes have a weak moment and do small programs in Fortran but I find programming in C or C++ thoroughly unenjoyable and avoid it. At first, I was considering doing this series using OCaml because the balance of functional and imperative features is very pleasing and it has a NumPy/SciPy-like library specifically for scientific computing - the wonderful Owl library. I decided against following this route because the language has some rough edges and alien syntax that would distract from the main point. Although, I might do a bit here and there, where I’ll be showing how to do parallel algorithms on a GPU and that will be most definitely with Futhark, which is a ML-style language, just like OCaml. Another option I considered was the Rust programming language. It’s fast, has good multidimensional array library with associated linear algebra and the built-in tooling is infinitely better than C++. But it’s a compiled, strict and static language which makes the development of programs slow and the code way too verbose.
In the end, I settled for the Julia language. It’s a general purpose language with strong focus on scientific computing and numerics that’s trying to, at least seemingly, replace languages like R, Matlab and (heavily library-dependent) Python. And it’s being rather good at it: it’s very approachable, fast (with a few caveats) and getting it running is as easy as downloading a single app or setting up a free account on a Jupyter notebook server. It’s by no means as fun as coding in Scheme or Common Lisp, but it gets the job done. Be warned, though, that it’s not all happy-fun-times. Julia was really hyped a couple of years ago, early in it’s life cycle before the 1.0 release, which lead to a surprising explosion of libraries and community support. The dust has settled and the language has eventually matured and stabilized but at the expense of a bunch of breaking changes. Because of this, many of the resources you find on the web will just not work (at least 2/3 of the linear algebra tutorials and cookbooks tell you to use functions like eigs which no longer exist in the main linear algebra library) and some of the really interesting libraries just died before they could have been useful (e.g. seamless GPU support for arrays doesn’t work if you’re using OpenCL like me). C’est la vie.
After all this, we can finally get to the real stuff…
Computers and vector spaces or: How many bits does it take to store the concept of a sine function?
In an ideal world, this is where I would write a small book worth of concepts and derivations relevant to Hamiltonian mechanics and quantum mechanics. I will not do that. Not because I want to keep it secret from you but because it would be just too much to write (Griffiths has almost 470 pages just for introduction to quantum mechanics). Those of you who will read the following few paragraphs and feel completely lost, will have to either ask me or the Internet, or search for a book or pdf. There’s plenty of resources. And don’t lose hope, things will get a bit more down-to-earth once the programming part begins.
Before we get to the nitty-gritty of quantum mechanics, let’s have a bit of an informal discussion in context of classical mechanics. Let’s start with the concept of a state. A state is (as described by Wikipedia) a complete description of a system. It’s a roster of all the relevant quantities that describe an object. For example in orbital dynamics a state of a planet would be determined by giving the planet’s position and momentum. You can also have systems that have more than just one planet, in which case, the state is determined by enumerating all of the individual planet’s positions and momenta. In classical mechanics, all these properties are just well defined numbers.
The second pair of concepts is the phase space and the equations of motion. The phase space is a space of all the possible states of a system. If we take the example of a point mass in 1D (or a marble on a string, if you want), it’s phase space will be represented by a 2D plane with momentum as one coordinate and position as another. Equations of motion are equations that describe the time evolution of a system - given an initial state, they will tell you in what state the system will be in the future (physics is really just a posh form of clairvoyance). As the system evolves, it will draw out a trajectory in the phase space which, in case of a time-invariant system, will be contour of constant energy. These trajectories are solutions to the Hamilton’s equations:
\[ \frac{d\mathbf{p}}{dt} = -\frac{\partial\mathcal{H}}{\partial\mathbf{q}} \quad , \quad \frac{d\mathbf{q}}{dt} = \frac{\partial\mathcal{H}}{\partial\mathbf{p}}.\]
\(\mathcal{H}\) is the famous classical Hamiltonian and \(\mathbf{q}\)s and \(\mathbf{p}\)s are generalized coordinates and corresponding momenta. Generalized, because the coordinates need not be just Cartesian coordinates - they can be the angles of a pendulum, displacements in a vibrating sheet or many other things.
Now that we have all that, let’s ruin it by introducing quantum mechanics.
The concept of state still exists, somewhat. The state is now represented by a vector in vector space called a Hilbert space and this change has some interesting consequences. The relationship between a state and observables got subverted, where the state can no longer be defined by an enumeration of numbers corresponding to its canonical positions and momenta as in classical physics. Conceptually, momentum or position (or any other measurable quantity) are not just numbers anymore, but operators and the possible values of measurable quantities are the eigenvalues of these operators, i.e. the possible values \(a\) (numbers) of any physical quantity \(A\) of a system in state \(\psi\) are \(\hat{A} \psi = a\psi\).
Lastly, Newton’s equations of motion are not a thing anymore and the evolution of a system is governed by the Schrödinger equation:
\[i\hbar \frac{d}{dt} \left| \Psi(t) \right> = \hat{H} \left| \Psi(t)\right>\]
or the time-independent form:
\[ \hat{H} \left| \Psi \right> = E \left| \Psi \right>,\]
With \(\hat{H}\) the Hamiltonian operator and \(E\) it’s eigenvalues - the possible energies of the system.
Here we introduced the bra-ket notation commonly seen in physics. As I mentioned before, the state \(\Psi\) lives in a Hilbert space of wave functions and the angled braces notation is used for vectors in this space. Usually, a state \(\left| ... \right>\) (ket) is represented as a column vector and \(\left< ... \right|\) (bra) is it’s hermitian conjugate. An inner product of two vectors in this notations is written as \(\left< ... | ... \right>\) and the outer product as \(\left| ... \right> \left< ... \right|\).
The only problem is that Hilbert spaces and their vectors and operators are not something you can store in computer memory - it’s a fully abstract concept and my CPU doesn’t have registers for spherical harmonics or any other functions. How can we do simulations of such systems?
It’s easy, actually. When I tell you to make a program that calculates a position in real space, what will you do? It’s not like a CPU has registers for meters either (I’m betting that the CPU manufacturers are getting hard from the idea of region-locked hardware based on local preferred unit system). You’ll chose a basis which defines a good coordinate system and give me the position as a linear combination of the basis vectors - just numbers of how many of each there are between me and the position. What about operators? Well, consider what operators do: They transform vectors in the Hilbert space. Sounds suspiciously similar to something doesn’t it?
More formally, if we have some complete orthonormal basis \(\{\phi _i \}_i\), with \(\left< \phi _i | \phi _j \right> = \delta _{ij}\), an operator \(\hat{A}\) can be represented as:
\[ \hat{A} = \mathbf{1} \cdot \hat{A} \cdot \mathbf{1} = \left[ \sum_i \left| \phi_i \right> \left< \phi_i \right| \right] \cdot \hat{A} \cdot \left[ \sum_j \left| \phi_j \right> \left< \phi_j \right| \right] = \sum_{ij} \left<\phi_i\right| \hat{A} \left|\phi_j\right> \left| \phi_i \right> \left< \phi_j \right| = \sum_{ij} A_{ij} \left| \phi_i \right> \left< \phi_j \right|,\]
Where \(A_{ij}\) are elements of a complex matrix \(\mathbf{A}\), i.e. something that we can store in a computer memory and have the CPU do calculations with. You can do the same thing with a vector in the Hilbert space:
\[ \left| \psi \right> = \mathbf{1} \cdot \left| \psi \right> = \left[ \sum_i \left| \phi_i \right> \left< \phi_i \right| \right] \cdot \left| \psi \right> = \sum_i \left<\phi_i | \psi \right> \left| \phi_i \right> = \sum_i \psi_i \left|\phi_i\right>,\]
where \(\psi_i\) are components of a complex vector \(\mathbf{\psi}\). Just an array of complex numbers. Armed with this representation, we can finally write down computer-solvable problems, for example the time-independent Schrödinger equation:
\[ \sum_j H_{ij} \psi_j = E\psi_i . \]
It’s fascinating how dynamics of systems that might be fundamentally disconnected from the reality of our daily lives, described by mathematical objects which simply don’t exist in physical world can still be reduced to a simple algebraic problem solvable on a computer sitting on the desk in front of you. Unreasonable effectiveness of math indeed, Mr. Wigner.
That’s it for the this part. If you feel let down by reading the name of this series and then just seeing a bunch of text and equations with no real computation don’t worry. Other parts of the series won’t be like this. The next part is already out, as I wrote both of them back-to-back, but I would recommend to take a break. Stare at the last three equations. There is an important message being conveyed in the last couple of paragraphs and it would be a shame if it got lost before moving on to the number crunching.