How does Metadynamics calculate the potential mean force (PMF) or the free energy profile of a reaction coordinate?

Metadynamics works by throwing sands to a valley, in the dark. An illustration of "throwing sands in the dark" :) First of all, what is the free energy surface (FES), in comparison with the less obscure potential energy surface? A system with N atoms has 3N coordinates ([math]\textbf{q}[/math]) in the configurational space. At a constant temperature, the canonical configurational partition function is: [math]Z=\int d\textbf{q} \rho(\textbf{q}) [/math], with the density of each state [math]\rho(\textbf{q})=e^{-\beta V(q)}[/math] Imagine how difficult it is to analyze or visualize the density map in a 3N dimensional space! To alleviate this difficulty, the reaction coordinates (or sometimes called the "collective variables") were introduced. Given a reaction coordinate [math]S_1(\textbf{q})[/math] which is a function of [math]\textbf{q}[/math], its density distribution is: [math]\rho(S_1)=\int d\textbf{q} \rho(\textbf{q})\delta(S_1-S_1(\textbf{q}))[/math]. This can be useful because, most dimensions of [math]\textbf{q}[/math] do not contribute to the transitions that we are interested in, but form an enormous number of shallow local minima. By tactically choosing the function  [math]S_1(\textbf{q})[/math], the important degrees of freedoms can be singled out. (How to pick good reaction coordinates is another very difficult topic.) Selecting reaction coordinates is equivalent to obtaining a coarse-grained description of the system. FES for a cluster of 38 atoms (with two reaction coordinates) Having defined a reaction coordinate [math]S_1(\textbf{q})[/math] , the concept of a free energy can be introduced: [math]F(S_1)=-\frac{1}{\beta}ln(\rho(S_1))[/math], and the extension to a set of multiple dimensional reaction coordinates is straight forward. Right now we'll focus on one-dimension. Now, the problem becomes how to sample the FES of the reaction coordinate. We have already simply the problem tremendously at the moment. The original problem is to sample the fuzzy 3N-dimentional potential energy surface, now we only need to work on one dimension. So, suppose we now start a molecular dynamics (MD) simulation to sample the free energy profile of [math]S_1[/math]. What problem can we encounter here? Here is the problem that ALWAYS bothers the atomistic simulation people: the time span for MD simulations is so short (~ns), and the system does not have enough time to go through the transition states. The system can be trapped in a "valley", and the FES is not sufficiently sampled. This is when Michele Parrinello and Alessandro Laio (and some others) came up with a brilliant idea: to add a biased potential ([math]B(S_1)[/math]) to "push" the system out of the valley. But the FES is unknown a priori, how do we know where are the valleys? The system is totally in the dark! Here comes the other ingenious idea: we know the system is in a valley, if the system tends to stay there. Combine these two key elements, we have reached the essence of  Metadynamics. For each a few MD steps, a repulsive Gaussian bias  potential (B) of height w and width σ is added to the potential energy  surface: [math]V(\textbf{q}) \leftarrow V(\textbf{q}) + B(S_1(V(\textbf{q})),t)[/math]. This  Gaussian bias is centered on positions in collective variables space  that the system is visiting. The repulsive bias acts as “sand”, which  fills the valleys on the FES. The "sands" accumulates, and encourages  the system to go through the barriers and explore the new spaces. Finally, by adding up all the bias potential [math]B(S_1(V(\textbf{q})),t)[/math] at each step, one can the NEGATIVE value of the FES. This is equivalent to measure how deep is the valley by counting how many sands were thrown in.