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, implement the method, run the simulations described below, and present what you have learned in a 10-minute talk at your final defense. Build your talk around the following guiding principles:

  1. Which advanced method have you chosen, how does the algorithm work, and what is the physical motivation behind it?

  2. What is the problem with standard simulations that the method attempts to solve?

  3. What are your results? Show the figures from your results notebook, with error bars, and compare them with what you expected.

  4. What are the limitations of the method?

  5. 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. Build on your own code from Projects I and II, in the same repository; do not start from scratch. Work on a new branch named final, created from your project2 branch: git checkout -b final project2 (the first time you push it: git push -u origin final), or with the Codespaces interface as described in the README. Submit your results as a short notebook with the required figures; captions are enough, no written summary is needed.

How you are assessed. The final project is 50% of your grade, entirely through an in-person defense after the break, about 20 minutes long and without AI tools:

  1. A presentation: up to 10 slides in 10 minutes, answering the five questions above.

  2. An oral exam on your own code (about 10 minutes), in the same format as the oral exams during the semester: explain a part of your implementation and why you implemented it that way; read a different version of that part and explain what it does, how it differs from yours, and what would happen if you used it; and answer a question on the theory behind your method or on one of your figures.

Before the defense, your runPIMD or runMetaD method is checked automatically (see your option below). Bonus points are given when the bonus checks pass.

Option A - Path Integral Molecular Dynamics / Monte Carlo#

The gist: Throughout this class, we have only discussed classical dynamics and classical statistical mechanics. 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:

  1. M.E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulations [F1], Chapter 12.

  2. 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\), with \(\hbar \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, keeping \(\beta \hbar \omega_0 / P\) fixed (i.e., increasing the temperature by a factor of two allows you to halve the number of beads).

Figures to show:

  1. The average energy against \(\beta \hbar \omega_0 = 1, \dots, 6\), with error bars from independent runs, together with the exact curve.

  2. The energy against the number of beads at the lowest temperature (convergence with \(P\)).

What to deliver: add a method with exactly this signature to your Simulation class:

def runPIMD(self, betahw, P, nsteps, seed, method="MD"):
    """Return (E_mean, E_error): the average total energy and its statistical error, in meV."""

It simulates one Ar atom in the one-dimensional harmonic trap above (\(\hbar \omega_0 = 50\) meV) at \(\beta \hbar \omega_0 =\) betahw, with P beads, nsteps steps and random seed seed. method is "MD" or, for the bonus, "MC". How you implement it (thermostat, energy estimator, how you estimate the error, etc.) is up to you, and you will explain your choices in the talk. Use this method to produce the figures in your results notebook.

Before the defense, your method will be checked automatically. We will call it as

sim = Simulation(dt=1e-15, L=1.0, ftype="Harm")
E, dE = sim.runPIMD(betahw=..., P=..., nsteps=100000, seed=1)

for several temperatures from the list above, with numbers of beads that we choose large enough to be converged, and compare E with the exact result you derive from the partition function. We expect agreement within about three times your error estimate dE. So your error bar must be honest: neither too large nor too small.

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: implement method="MC" in runPIMD as well; it should give the same energies as method="MD" within the errors.

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 (see the discussion of energy conservation in the Molecular dynamics chapter). Since every time step requires evaluating the forces on all atoms, the computational cost limits typical MD simulations 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:

  1. G. Bussi, A. Laio, Using metadynamics to explore complex free-energy landscapes [F2].

  2. A. Valsson, P. Tiwary, and M. Parrinello, Enhancing Important Fluctuations: Rare Events and Metadynamics from a Conceptual Viewpoint [F3].

Simulation details: Perform Metadynamics simulations at \(298\) K for an Ar atom in the double well potential energy surface:

\[ V(x) = A x^4 - Bx^2, \]

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. Accumulate the bias on a grid, for computational efficiency.

Figures to show:

  1. \(x(t)\) of the unbiased run: no transition.

  2. \(x(t)\) of a Metadynamics run, showing the transitions, together with the bias at the end of the run.

  3. The free energy from the bias, averaged over independent runs with error bars, against the exact \(V(x)/k_B T\) after shifting by a constant.

What to deliver: add a method with exactly this signature to your Simulation class:

def runMetaD(self, nsteps, height, width, pace, seed, biasfactor=None):
    """Return (x, F): a grid of positions in Angstrom, covering at least -1.5 to 1.5 Angstrom,
    and the free energy on it, in units of kB*T, with min(F) = 0."""

It simulates one Ar atom in the double well above at \(298\) K, starting in the left well (\(x = -1\) Angstrom), for nsteps steps with random seed seed. A Gaussian of height height (in units of \(k_B T\)) and width width (in Angstrom) is deposited every pace steps; height=0 is an unbiased run. biasfactor is for the bonus (well-tempered Metadynamics). Implement the double-well potential as a new force method, evalDoubleWell(A, B), called through evalForce like the other potentials. How you implement the rest (thermostat, how you reconstruct the free energy, etc.) is up to you, and you will explain your choices in the talk.

Before the defense, your method will be checked automatically. We will call it as

sim = Simulation(dt=1e-15, L=1.0, ftype="DoubleWell")
x, F = sim.runMetaD(nsteps=..., height=..., width=..., pace=..., seed=s)   # several seeds s

with parameters of our choice, and compare the average of F over the seeds with the exact free energy you expect from the potential, over the wells and the barrier, after shifting it by a constant. The deviation should be small: about the size of the differences between your independent runs.

Bonus: Implement the well-tempered version of Metadynamics through the biasfactor argument of runMetaD, 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.

References#

F1

Mark Tuckerman. Statistical Mechanics: Theory and Molecular Simulation. Oxford Graduate Texts. Oxford University Press, Oxford, New York, February 2010. ISBN 978-0-19-852526-4.

F2

Giovanni Bussi and Alessandro Laio. Using metadynamics to explore complex free-energy landscapes. Nature Reviews Physics, 2(4):200–212, April 2020. URL: https://doi.org/10.1038/s42254-020-0153-0, doi:10.1038/s42254-020-0153-0.

F3

Omar Valsson, Pratyush Tiwary, and Michele Parrinello. Enhancing Important Fluctuations: Rare Events and Metadynamics from a Conceptual Viewpoint. Annual Review of Physical Chemistry, 67(1):159–184, 2016. _eprint: https://doi.org/10.1146/annurev-physchem-040215-112229. URL: https://doi.org/10.1146/annurev-physchem-040215-112229, doi:10.1146/annurev-physchem-040215-112229.