Final project
Contents
Final project¶
Congratulations! You have successfully finished your first two projects, implementing your very own MD and Metropolis MC codes. You have also learned some of the theory underlying both methods, including analytical mechanics, statistical mechanics, and random processes. This is a highly non-trivial achievement.
Now, to the final project that will determine \(50\%\) of your grade. This time, I will not teach you the theory myself. Instead, I will briefly describe two advanced simulation methods that you can choose from, and refer you to the relevant reading material. Your job is to read the material and prepare a short summary, no longer than 5 pages, explaining what you have learned. Focus on the following guiding principles:
Which advanced method have you chosen to implement?
Explain what is the problem with standard simulations the method attempts to solve.
How does the algorithm work? Explain the physical motivation behind it.
Discuss the limitations of the method.
Can you think of ways to improve upon it?
Then, you need to implement your advanced method and run a simple simulation, as described below. You need to create a repository with the source code of your simulation. The summary should be submitted in .ipynb format. I recommend using some of the code we developed in class, and the repository, but you don’t have to.
Option A - Path Integral Molecular Dynamics / Monte Carlo¶
The gist: Throughout this class, we have only discussed classical dynamics and classical statistical mehcanics. What about quantum mechanics? Can any of the methods we learned be utilized for quantum systems? The answer came from the work of the famous physicist Richard P. Feynman.
He developed a new formulation of quantum mechanics using a tool called path integrals. His formulation is entirely equivalent to the formulations of quantum mechanics you already know, like the Schrodinger equation, but it is especially illuminating in connecting quantum mechanics with classical mechanics. It essentially explains how the concept of the classical action appears also in quantum systems. But, more concretely to our purposes, it also shows how one can use classical MD or MC simulations to evaluate thermal (NVT) expectation values of quantum systems. Mind blown!
This is possible because there is a mapping of each quantum system to a larger system of classical particles which has the same partition function. This is achieved by replacing every quantum particle with a classical ring polymer composed of \(P\) beads which are copies of the original particle that are connected with harmonic springs. Use the following two sources to read all about it:
Reading material:
M.E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulations [1], Chapter 12.
M. Ceriotti, V, Kapil, Y. Litman, T. Markland and M. Rossi, Booklet of the CECAM Flagship School on Path Integral Quantum Mechanics. Sections 1-2, available on the moodle.
Simulation details: Evaluate the average total energy of a quantum particle (Ar atom) in a one-dimensional Harmonic trap of frequency \(\omega_0 = 50 meV\) as a function of inverse temperature \(\beta=(k_BT)^{-1}\). Use values of \(\beta \hbar \omega_0 = 1,2,3,4,5,6\). To estimate the statistical uncertainty of your results run several independent simulations at every temperature. Show that the energy converges to the correct ground state energy as you decrease the temperature (increase \(\beta\)). Compare the results at all temperatures with the analytical results that can be derived from the exact expression for the partition function. Don’t forget to show the convergence of your simulation with respect to the number of beads. You can do this only for the lowest temperature and then just scale the number of beads with temperature (i.e., decreasing the temperature by a factor of two leads to the same decrease in the number of beads).
Bonus: You can use either MD or MC simulations but PIMD is described in more detail in Tuckerman. You will get bonus points if you do both!
Option B - Metadynamics¶
The gist: MD simulations are an extremely powerful tool that can be thought of as a virtual microscope capturing the motion of every atom in a solid, liquid, gas, etc. The temporal resolution is determined by the time step such that the total energy is conserved. A typical time step is of the order of \(1 fs\), which is limited by the highest frequency normal modes of the system. The algorithms for time propagation accumulate a numerical error, which depends on the time step. Therefore, MD simulations are usually stable up to about \(10^9\) steps. This results in a total simulation time of the order of \(\sim 1 \mu s\).
Many fascinating chemical, physical, and biological processes occur on a much longer time scale. Important examples include protein folding, crystal nucleation, and growth, or chemical reactions with a high energy barrier. Standard MD simulations cannot describe such processes, which we call rare events.
Metadynamics is an algorithm that allows for enhanced sampling of such rare events by adding a bias potential to the system, which is built as the simulation progresses, pushing it from one state to another. It also provides a way of evaluating the free energy surface of the process without the added bias. This is achieved in practice by depositing a Gaussian every 100-500 steps during the simulation. You can read all about it in the following two review papers:
Reading material:
G. Bussi, A. Laio, Using metadynamics to explore complex free-energy landscapes [3].
A. Valsson, P. Tiwary, and M. Parrinello, Enhancing Important Fluctuations: Rare Events and Metadynamics from a Conceptual Viewpoint [4].
Simulation details: Perform Metadynamics simulations at \(298 K\) for a Ar atom in the double well potential energy surface:
where \(A=4.11 \cdot 10^{20} \, J m^{-4}, B = 8.22 \,J m^{-2}\).
First show that this potential leads to the transition from one well to another being a rare event. Then, use the position of the particle as a collective variable to perform a Metadynamics simulation biasing the system to go from one well to the next. Compare the Metadynamics simulations with the unbiased simulations using the same potential. Finally, obtain the unbiased free-energy surface of the system, using one of the approaches described in the review, and compare it with your expectation based on the potential above. Evaluate the error on your free energy surface estimate by running several independent trajectories and evaluating their standard deviation. Note that it is highly recommended to use a grid for accumulating the bias for computational efficiency, but it is not mandatory.
Bonus: Implement the well-tempered version of Metadynamics, and compare its performance against standard metadynamics. If you do implement well-tempered Metadynamics, note that you will have to use a grid to calculate \(c(t)\), the reweighting factor, to obtain the unbiased free energy.