Computational Physics (Lecture 18)

1. Computational Physics(Lecture 18) PHY4370

2. Molecular dynamics simulations • Most physical systems are collections of interacting objects. • a drop of water contains more than 1022water molecules. • a galaxy is a collection of millions and millions of stars. • no analytical solution that can be found for an interacting system with more than two objects. • We can solve the problem of a two-body system, such as the Earth–Sun system, analytically, but not a three-body system, such as the Moon–Earth–Sun system.

3. The situation is similar in quantum mechanics, • one can obtain the energy levels of the hydrogen atom • (one electron and one proton) analytically, but not those the helium atom (two electrons and a nucleus) • Numerical techniques are needed to study a system of a large number of interacting objects, or the so-called many-body system.

4. a distinction between three-body systems such as the Moon–Earth–Sun system and a more complicated system, such as a drop of water. • Statistical mechanics has to be applied to the latter.

5. General behavior of a classical system • MD solves the dynamics of a classical many-body system described by the Hamiltonian

6. The methods for solving Newton’s equation discussed early Chapters can be used to solve the above equation set. • However, those methods are not as practical as MD, in terms of the speed and accuracy of the computation and given the statistical nature of large systems.

7. There are several ways to simulate a many-body system. • Most simulations are done either through a stochastic process, such as the Monte Carlo simulation • or through a deterministic process, such as a molecular dynamics simulation. • Some numerical simulations are performed in a hybridized form of both, • for example, Langevindynamics, and Brownian dynamics.

8. Another issue is the distribution function of the system. • In statistical mechanics, each special environment is dealt with by way of a special ensemble.

9. For an isolated system, we use the microcanonical ensemble, • which assumes a constant total energy, number of particles, and volume. • A system in good contact with a thermal bath is dealt with using the canonical ensemble • which assumes a constant temperature, number of particles, and volume (or pressure). • Forany given ensemble, the system is described by a probability function W(R, P), • which is in general a function of phase space, consisting of all coordinates and momenta of the particles R = (r1, r2, . . . , rN) and P = (p1, p2, . . . , pN), and other quantities, such as temperature, total particle number of the system, and so forth.

10. For the canonical ensemble, we have: • where T is the temperature of the system, kB is the Boltzmann constant, and N is a normalization constant.

11. We can separate the position dependence and momentum dependence in W(R, P) if they are not coupled in H. • Any average of the momentum-dependent quantity becomes quite simple because of the quadratic behavior of the momentum in H. • So we concentrate on the position dependence here. • The statistical average of a physical quantity A(R, P) is then given by:

12. if the system is ergodic: that is, every possible state is accessed with an equal probability. • Because molecular dynamics simulations are deterministic in nature, almost all physical quantities are obtained through time averages. • Sometimesthe average over all the particles is also needed to characterize the system. • Forexample, the average kinetic energy of the system can be obtained from any ensemble average, and the result is given by the partition theorem

14. Basic methods for many-body systems

15. where ρ is the average density from the points r = 0 and r.

16. If the angular distribution is not the information needed, we can take the angular average to obtain the radial distribution function:

17. The behavior of the pair-distribution function can provide a lot of information regarding the translational nature of the particles in the system. • For example, a solid structure would have a pair-distribution function with sharp peaks at the distances of nearest neighbors, next nearest neighbors, and so forth. • If the system is a liquid, the pair-distribution still has some broad peaks at the average distances of nearest neighbors, next nearest neighbors, and so forth, but the feature fades away after several peaks.

19. Similarly, we can obtain the radial distribution function numerically. • We need to measure the radius r from the position of a specific particle r_i, and then the • radial distribution function g(r ) is the probability of another particle’s showing up at a distance r .

20. Numerically, we have

21. The dynamics of the system can be measured from the displacement of the particles in the system. • We can evaluate the time dependence of the mean-square displacement of all the particles,

22. The very first issue in numerical simulations for a bulk system is howtoextend a finite simulation box to model the nearly infinite system. • A common practice is to use a periodic boundary condition, that is, to approximate an infinite system by piling up identical simulation boxes periodically. • A periodic boundary condition removes the conservation of the angular momentum of the simulated system (particles in one simulation box), but still preserves the translational symmetry of the center of mass. • So the temperature is related to the average kinetic energy byE_K = 3/2(N − 1) k_BT, • where the factor (N − 1) is due to the removal of the rotation around the center of mass.

23. how to include the interactions among the particlesin different simulation boxes • If the interaction is a short-range interaction, • one can truncate it at a cut-off length r_c. • The interaction V(rc) has to be small enough that the truncation does not affect the simulation results significantly. • Atypical simulation box usually has much larger dimensions than r_c. • For a threedimensionalcubic box with sides of length L, the total interaction potential can be evaluated with many fewer summations than N!/2, the number of possible pairs in the system. • For example, if we have L/2 > rc, and if |xi j |, |yi j |, and|zij | are all smaller than L/2, we can use Vi j = V(r_ij ); • otherwise, we use the corresponding point in the neighboring box. • In order to avoid a finite jump at the truncation, one can always shift the interaction to V(r ) − V(r_c) to make sure that it is zero at the truncation.

24. which can be evaluated quite easily, because at every time step the force fi j =−∇V(ri j ) is calculated for each particle pair.

25. The Verlet algorithm • Hamilton’s equations are equivalent to Newton’s equation • To simplify the notation, we can rewrite:

26. The Verlet algorithm can be started if the first two positions R0 and R1 of the particles are given. • However, in practice, only the initial position R0 and initial velocity V0 are given. • Therefore, we need to figure out R1 before we can start the recursion. • A common practice is to treat the force during the first time interval [0, τ] as a constant, and then to apply the kinematic equation to obtain G0 is the acceleration vector evaluated at the initial configuration R0.

27. Of course, the position R1 can be improved by carrying out the Taylor expansion to higher-order terms if the accuracy in the first two points is critical. • We can also replace G0 with the average (G0 + G1)/2, with G1 evaluated at R1.This procedure can be iterated several times before starting the algorithm for the velocity V1 and the next position R2.

28. The Verlet algorithm has advantages and disadvantages. • It preserves the time reversibility that is one of the important properties of Newton’s equation. • Therounding error may eventually destroy this time symmetry. • The error in the velocity is two orders of magnitude higher than the error in the position. • Inmany applications, we may only need information about the positions of the • particles, and the Verlet algorithm yields very high accuracy for the position. • If the velocity is not needed, we can totally ignore the evaluation of the velocity, since the evaluation of the position does not depend on the velocity at each time step. • The biggest disadvantage of the Verlet algorithm is that the velocity is evaluated one time step behind the position! • However, this lag can be removed if the velocity is evaluated directly from the force. • A two-point formula would yield: V_k+1 = V_k + τ G_k + O(τ^ 2).

29. Note that the evaluation of the position still has the same accuracy because the velocity is now updated according to Eq. (8.41), which provides the cancelation of the third-order term in the new position.