Chapter 1
Introduction
Context: We start by introducing the concept of the potential energy and the interatomic force. Those are the central ingredients to the molecular dynamics simulation method.
Additional resources:
- Chapter 1 of Allen & Tildesley, Computer Simulation of Liquids
- Chapters 1 and 2 of Interatomic potentials: Achievements and challenges
1.1 Structure of matter at the atomic scale
All matter is build out of quark and leptons, or perhaps even smaller particles, but for the sake of modeling the real material world the atom is the fundamental unit. Atoms can be described by nuclei and electrons or through “coarse-grained” models that ignore the fact that there are electrons. Both types of models are useful for describing materials, and the latter ones will be extensively used in this class.
Atoms in solids can arrange in different configurations that are called crystals when there is long-ranged order or glasses when there is not. (All solid matter typically has short-ranged order that is determined by the chemical bonds between atoms.) Atoms in solids are immobile and self-diffusion is limited. Conversely, liquids and gases are disordered (like glasses) but have mobile constituent atoms. Macroscopic object typically contain a lot of atoms – on the order of Avogadro’s constant \(N_A\approx 6\times 10^{23}\). The atomic-scale simulation techniques discussed in this class can at the time of this writing (the 2020s) casually treat on order of \(\sim 10^7\) atoms. This can be pushed to \(10^{13}\) Eckhardt et al. (2013); Tchipev et al. (2019) if you use the biggest computers available to us (in combination with computationally-efficient models). Of course, this boundary moves towards larger systems as computer technology evolves.
Nowadays, we can even observe matter at atomic scales and “see” individual atoms. The collaborative research center 103 has produced an extremely instructive video on the structure of specific type of alloys, dubbed “superalloy”, that is used in e.g. turbine blades. Enjoy the ride from the blade to the atom. This class is about modeling matter at the smallest scales that you see in this video.
https://www.youtube.com/embed/wYHch5QIWTQ
1.2 Interatomic forces and the potential energy
Atoms interact via forces. As Feynman put it in his famous lectures on physics, the fundamental truth about man’s understanding of the physical world is “that all things are made of atoms – little particles that move around in perpetual motion, attracting each other when they are a little distance apart, but repelling upon being squeezed into one another” Feynman et al. (1964). Indeed this statement summarizes the essence of the molecular dynamics simulation method.
As the simplest example why atoms attract each other, let us consider the example of simple salt, e.g. Na-Cl that we all have sitting in our kitchen. Na-Cl in its solid form is an ionic crystal. Na atoms have approximately a charge of \(q_{\text {Na}}=+1|e|\), where \(e\) is the electron charge, and Cl atoms have a charge of approximately \(q_{\text {Na}}=-1|e|\). The Coulomb interaction between these atoms is a fundamental force of nature. Basic physical principles tell us, that the interaction energy between a Na and a Cl atom is given by \begin {equation} V_{\text {Coulomb}}(r;q_{\text {Na}},q_{\text {Cl}}) = \frac {1}{4\pi \varepsilon _0} \frac {q_{\text {Na}} q_{\text {Cl}}}{r}. \end {equation} We also know that this energy is pair-wise additive, allowing us to write down the Coulomb interaction energy for Na-Cl consisting of \(N\) atoms, \begin {equation} E_{\text {Coulomb}}(\{\vec {r}_i\}) = \sum _{i=1}^N \sum _{j=i+1}^N V_{\text {Coulomb}}(r_{ij};q_i, q_j) = \frac {1}{4\pi \varepsilon _0} \sum _{i=1}^N \sum _{j=i+1}^N \frac {q_i q_j}{r_{ij}} \label {eqn:coulomb} \end {equation} where \(q_i\) is the charge on atom \(i\) and \(r_{ij}=|\vec {r}_i - \vec {r}_j|\) the distance between atom \(i\) and \(j\). Note that we have introduced — in passing — a central quantity of the molecular dynamics method, the atomic positions \(\vec {r}_i\) and Eq. \eqref{eqn:coulomb} indicates that the interaction energy depends on the positions of all atoms.
The Coulomb interaction has a singularity at \(r\to 0\). The attractive force between opposite charges becomes infinitely large. The salt crystal does not collapse because atoms are, as Feynman puts it, “repelling upon being squeezed into one another”. While the attraction between our Na and Cl atoms are described by a fundamental force of nature, it is more difficult to understand the origin of this repulsion. Hand-wavingly, it goes back to the fact that electrons are Fermions and electrons from the electron shells of Na and Cl therefore cannot exist at the same point in space (and the same spin state). This is the Pauli exclusion principle, and the resulting repulsive force is called Pauli repulsion.
Different models for the Pauli repulsion exist. While the Coulomb interaction is a fundamental force of nature, these models are approximations to the true quantum physics that is the origin of the repulsive form. Two common forms are exponential repulsion, \begin {equation} E_{\text {rep,exp}}(\{\vec {r}_i\}) = \sum _{i=1}^N \sum _{j=i+1}^N A e^{-r/\rho }, \end {equation} or an algebraic repulsion of the form \begin {equation} E_{\text {rep,12}}(\{\vec {r}_i\}) = \sum _{i=1}^N \sum _{j=i+1}^N A r^{-12}. \end {equation} Note that \(A\) and \(\rho \) are parameters, that need to be somehow determined. This can be done with the help of either experimental data or first-principles calculations, that treat the electrons explicitly. These parameters depend on the atom types under consideration and, in contrast to the parameter that show up in the Coulomb interaction (the permittivity \(\varepsilon _0\)), they are not universal.
For our Na-Cl model, we combine Coulomb interaction with an exponential repulsion, to give the total energy \begin {equation} E_{\text {pot}}(\{\vec {r}_i\}) = \sum _{i=1}^N \sum _{j=i+1}^N \left (A_{ij} e^{-r_{ij}/\rho _{ij}} + \frac {1}{4\pi \varepsilon _0} \frac {q_i q_j}{r_{ij}}\right ). \label {eqn:NaCl} \end {equation} This energy is called the potential energy and is the central property of an atomic-scale model. With Eq. \eqref{eqn:NaCl}, we have also encountered our first atomic-scale model for a real material, called the Born-Mayer potential. Potentials that can be decomposed as Eq. \eqref{eqn:NaCl} into pair-wise terms are called pair potentials. They are often written as \begin {equation} E_{\text {pot}}(\{\vec {r}_i\}) = \sum _{i=1}^N \sum _{j=i+1}^N V(r_{ij}), \end {equation} with \begin {equation} V(r_{ij}) = A_{ij} e^{-r_{ij}/\rho _{ij}} + \frac {1}{4\pi \varepsilon _0} \frac {q_i q_j}{r_{ij}} \end {equation} for the above potential. The quantity \(V(r_{ij})\) is called the pair interaction energy.
The most famous pair-potential is likely the Lennard-Jones potential. (Unlike Born and Mayer, Lennard-Jones is a single person.) Its pair interaction energy is given by \begin {equation} V(r_{ij}) = 4\varepsilon \left [\left (\frac {\sigma }{r}\right )^{12} - \left (\frac {\sigma }{r}\right )^6\right ]. \end {equation} The repulsive term \(\propto r^{-12}\) is one of the models for Pauli repulsion discussed above. The attractive term \(\propto r^{-6}\) arises from London dispersion interactions. Dispersion forces exist between all atoms, even uncharged molecules or noble gases. They are widely employed for the nonbonded portion of what are called valence force-fields. Simple Lennard-Jones systems are often used to study generic phenomena found in real materials, e.g. the glass transition or plasticity of amorphous materials. However, there are limitations to pair potentials, and more sophisticated potential energy models have been developed over the past decades to address shortcomings of simple potentials. In particular, the assumption of pair additive intrinsic to the equations of this chapter is often not valid. We will discuss how to lift these limitations in Chapter ??.
Note: A repulsive term of the form \(r^{-12}\) is advantageous from a simulation point of view, since it is faster to compute than an exponential. This has helped popularize the Lennard-Jones potential in the early days of molecular dynamics simulations.
Writing the potential energy of a system of particles allows the derivation of the forces acting on these particles: \begin {equation} \vec {f}_k = -\frac {\partial }{\partial \vec {r}_k} E_{\text {pot}}(\{\vec {r}_i\}). \label {eq:forces} \end {equation} These forces are the essential ingredient to molecular dynamics, as they determine the motion of the atoms, in accordance to Newton’s second law.
The potential energy itself describes what is called the potential energy landscape. The potential energy landscape depends on \(3N\) degrees of freedom (as compared to the landscape we experience while walking, which depends on \(2\) degrees of freedom); it is therefore an object that is complex to visualize. Simplifying some of its aspects is the core of molecular statics. For example, it is often important to identify the ground-state of a system; this is the most stable configuration of a material and has the lowest possible potential energy. There is usually some crystal that is lower in energy than the energy of a glass with the same stoichiometry. Yet, in many real-world engineering applications, the materials are not in their crystalline ground-state: the most common material we encounter with this property may be window glass. In molecular statics, we therefore seek to enumerate those local minima of the potential energy landscape and the energy barriers between them.
Since the dynamics of a molecular system is determined by the forces acting on the individual atoms, we only need to specify the potential energy up to a constant, which disappears in the derivative Eq. \eqref{eq:forces}. We can therefore measure the potential energy with respect to any reference configuration. This reference configuration is often the atomized state of the material, where all the constituent atoms sit individually in vacuum and are not interacting with each other. If this reference situation is assigned the energy \(0\), then the potential energy is generally negative, because if it was positive the system would spontaneously atomize. (Remember, any physical system evolves to a state of lower energy.)
Bibliography
W. Eckhardt, A. Heinecke, R. Bader, M. Brehm, N. Hammer, H. Huber, H.-G. Kleinhenz, J. Vrabec, H. Hasse, M. Horsch, M. Bernreuther, C. W. Glass, C. Niethammer, A. Bode, and H.-J. Bungartz. 591 TFLOPS multi-trillion particles simulation on SuperMUC. In Supercomputing, pages 1–12. Springer Berlin Heidelberg, 2013. URL https://doi.org/10.1007/978-3-642-38750-0_1.
R. P. Feynman, R. B. Leighton, and M. Sands. The Feynman Lectures on Physics. Addison-Wesley, New York, 1964.
N. Tchipev, S. Seckler, M. Heinen, J. Vrabec, F. Gratl, M. Horsch, M. Bernreuther, C. W. Glass, C. Niethammer, N. Hammer, B. Krischok, M. Resch, D. Kranzlmüller, H. Hasse, H.-J. Bungartz, and P. Neumann. TweTriS: Twenty trillion-atom simulation. Int. J. High Perform. Comput. Appl., 33(5):838–854, Sept. 2019. URL https://doi.org/10.1177/1094342018819741.