diff --git a/user-doc/tutorials/a-trieste-4.txt b/user-doc/tutorials/a-trieste-4.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fe5e4e24361566d6e2a03009c766f9a527541e0b
--- /dev/null
+++ b/user-doc/tutorials/a-trieste-4.txt
@@ -0,0 +1,269 @@
+/**
+\page trieste-4 Trieste tutorial: Metadynamics simulations with PLUMED
+
+\section trieste-4-aims Aims
+
+The aim of this tutorial is to train the users to perform 
+metadynamics simulations with PLUMED, analyze the results, calculating free-energies as a function 
+of the collective variables used and estimating the associated error.
+
+\section trieste-4-objectives Objectives
+
+Once this tutorial is completed students will be able to:
+
+- Write the PLUMED input file to perform metadynamics simulations 
+- Calculate free energies from a metadynamics run
+- Assert the choice of the collective variable
+- Evaluate the convergence of metadynamics simulations
+- Compute the error associated to the reconstructed free energies 
+
+\section trieste-4-resources Resources
+
+The \tarball{trieste-4} for this project contains the following files:
+- diala.pdb : a PDB file for alanine dipeptide in vacuo
+- topol.tpr : GROMACS run file to run MD of alanine dipeptide
+- XXXX.py   : a python script to block analysis
+
+This tutorial has been tested on a pre-release version of version 2.4. However, it should not take advantage
+of 2.4-only features, thus should also work with version 2.3.
+
+\section trieste-4-intro Introduction
+
+We have seen that PLUMED can be used to compute collective variables. However, PLUMED
+is most often use to add forces on the collective variables. To this aim,
+we have implemented a variety of possible biases acting on collective variables.
+The complete documentation for
+all the biasing methods available in PLUMED can be found at the \ref Bias page.
+In the following we will see how to build an adaptive bias potential with metadynamics.
+Here you can find a brief recap of the metadynamics theory.
+
+\hidden{Summary of theory}
+
+In metadynamics, an external history-dependent bias potential is constructed in the space of 
+a few selected degrees of freedom \f$ \vec{s}({q}) \f$, generally called collective variables (CVs) \cite metad.
+This potential is built as a sum of Gaussians deposited along the trajectory in the CVs space:
+
+\f[
+V(\vec{s},t) = \sum_{ k \tau < t} W(k \tau)
+\exp\left(
+-\sum_{i=1}^{d} \frac{(s_i-s_i({q}(k \tau)))^2}{2\sigma_i^2}
+\right).
+\f]
+
+where \f$ \tau \f$ is the Gaussian deposition stride, 
+\f$ \sigma_i \f$ the width of the Gaussian for the ith CV, and \f$ W(k \tau) \f$ the
+height of the Gaussian. The effect of the metadynamics bias potential is to push the system away 
+from local minima into visiting new regions of the phase space. Furthermore, in the long
+time limit, the bias potential converges to minus the free energy as a function of the CVs:
+
+\f[
+V(\vec{s},t\rightarrow \infty) = -F(\vec{s}) + C.
+\f]
+
+In standard metadynamics, Gaussians of constant height are added for the entire course of a 
+simulation. As a result, the system is eventually pushed to explore high free-energy regions
+and the estimate of the free energy calculated from the bias potential oscillates around
+the real value. 
+In well-tempered metadynamics \cite Barducci:2008, the height of the Gaussian 
+is decreased with simulation time according to:
+
+\f[
+ W (k \tau ) = W_0 \exp \left( -\frac{V(\vec{s}({q}(k \tau)),k \tau)}{k_B\Delta T} \right ),
+\f]
+
+where \f$ W_0 \f$ is an initial Gaussian height, \f$ \Delta T \f$ an input parameter 
+with the dimension of a temperature, and \f$ k_B \f$ the Boltzmann constant. 
+With this rescaling of the Gaussian height, the bias potential smoothly converges in the long time limit,
+but it does not fully compensate the underlying free energy:
+
+\f[
+V(\vec{s},t\rightarrow \infty) = -\frac{\Delta T}{T+\Delta T}F(\vec{s}) + C.
+\f]
+
+where \f$ T \f$ is the temperature of the system.
+In the long time limit, the CVs thus sample an ensemble
+at a temperature \f$ T+\Delta T \f$ which is higher than the system temperature \f$ T \f$.
+The parameter \f$ \Delta T \f$ can be chosen to regulate the extent of free-energy exploration:
+ \f$ \Delta T = 0\f$ corresponds to standard MD, \f$ \Delta T \rightarrow \infty \f$ to standard
+metadynamics. In well-tempered metadynamics literature and in PLUMED, you will often encounter
+the term "biasfactor" which is the ratio between the temperature of the CVs (\f$ T+\Delta T \f$) 
+and the system temperature (\f$ T \f$):
+
+\f[
+\gamma = \frac{T+\Delta T}{T}.
+\f]
+
+The biasfactor should thus be carefully chosen in order for the relevant free-energy barriers to be crossed
+efficiently in the time scale of the simulation.
+ 
+Additional information can be found in the several review papers on metadynamics 
+\cite gerv-laio09review \cite WCMS:WCMS31 \cite WCMS:WCMS1103.
+
+\endhidden
+
+We will play with a toy system, alanine dipeptide simulated in vacuo using the AMBER99SB force field (see Fig. \ref trieste-4-ala-fig).
+This rather simple molecule is useful to benchmark data analysis and free-energy methods.
+This system is a nice example because it presents two metastable states separated by a high free-energy barrier.
+It is conventional use to characterize the two states in terms of Ramachandran dihedral angles, which are denoted with \f$ \Phi \f$ and \f$ \Psi \f$ in Fig. \ref trieste-4-transition-fig .
+
+\anchor trieste-4-ala-fig
+\image html belfast-2-ala.png "The molecule of the day: alanine dipeptide."
+
+\anchor trieste-4-transition-fig
+\image html belfast-2-transition.png "Two metastable states of alanine dipeptide are characterized by their Ramachandran dihedral angles."
+
+
+\section trieste-4-ex-1 Exercize 1: my first metadynamics calculation
+
+In this excercise we will setup and perform a well-tempered metadynamics run using the backbone dihedral \f$ \phi \f$
+as collective variable. During the calculation, we will also monitor the behavior of the other backbone dihedral \f$ \psi \f$.
+
+Here you can find a sample `plumed.dat` file that you can use as a template.
+Whenever you see an highlighted \highlight{FILL} string, this is a string that you should replace.
+
+\plumedfile
+# Compute the backbone dihedral angle phi, defined by atoms C-N-CA-C
+phi: TORSION ATOMS=__FILL__
+# Compute the backbone dihedral angle psi, defined by atoms N-CA-C-N
+psi: TORSION ATOMS=__FILL__
+
+# Activate well-tempered metadynamics in phi
+metad: __FILL__ ARG=phi...
+# Deposit a Gaussian every 500 time steps, with initial height equal
+# to 1.2 kJoule/mol, biasfactor equal to 10.0
+PACE=500 HEIGHT=1.2 BIASFACTOR=10.0
+# Gaussian width (sigma) should be chosen based on CV fluctuation in unbiased run
+SIGMA=__FILL__
+# Gaussians will be written to file and also stored on grid
+FILE=HILLS GRID_MIN=-pi GRID_MAX=pi
+...
+
+# Print the collective variables and the value of the bias potential on COLVAR file
+PRINT ARG=phi,psi,__FILL__ FILE=COLVAR STRIDE=10
+\endplumedfile
+
+Once your `plumed.dat` file is complete, you can run a XX-ns long metadynamics simulations with the following command
+\verbatim
+> gmx mdrun -s topol.tpr -plumed plumed.dat 
+\endverbatim
+
+During the metadynamics simulation, PLUMED will create two files, named COLVAR and HILLS.
+The COLVAR file contains all the information specified by the PRINT command, in this case
+the value of the CVs every 10 steps of simulation, along with the current value of the metadynamics bias potential. 
+We can use `gnuplot` to visualize the behavior of the CV during the simulation, as reported in the COLVAR file:
+
+\verbatim
+gnuplot> p "COLVAR" u 1:2
+\endverbatim
+
+\anchor trieste-4-phi-fig
+\image html munster-metad-phi.png "Time evolution of the CV phi during the first 2 ns of a metadynamics simulation of alanine dipeptide in vacuum."
+
+By inspecting Figure \ref trieste-4-phi-fig, we can see that the system is initialized in one of the two metastable
+states of alanine dipeptide. After a while (t=0.1 ns), the system is pushed
+by the metadynamics bias potential to visit the other local minimum. As the simulation continues,
+the bias potential fills the underlying free-energy landscape, and the system is able to diffuse in the
+entire phase space.
+
+The HILLS file contains a list of the Gaussians deposited along the simulation.
+If we give a look at the header of this file, we can find relevant information about its content:
+
+\verbatim
+#! FIELDS time phi psi sigma_phi sigma_psi height biasf
+#! SET multivariate false
+#! SET min_phi -pi
+#! SET max_phi pi
+#! SET min_psi -pi
+#! SET max_psi pi
+\endverbatim 
+
+The line starting with FIELDS tells us what is displayed in the various columns of the HILLS file:
+the time of the simulation, the value of phi and psi, the width of the Gaussian in phi and psi,
+the height of the Gaussian, and  the biasfactor.
+We can use the HILLS file to visualize the decrease of the Gaussian height during the simulation,
+according to the well-tempered recipe:
+
+\anchor trieste-4-phihills-fig
+\image html munster-metad-phihills.png "Time evolution of the Gaussian height."
+
+If we look carefully at the scale of the y-axis, we will notice that in the beginning the value
+of the Gaussian height is higher than the initial height specified in the input file, which should be 1.2 kJoule/mol.
+In fact, this column reports the height of the Gaussian rescaled by the pre-factor that
+in well-tempered metadynamics relates the bias potential to the free energy.
+
+
+
+One can estimate the free energy as a function of the metadynamics CVs directly from the metadynamics
+bias potential. In order to do so, the utility \ref sum_hills should be used to sum the Gaussians
+deposited during the simulation and stored in the HILLS file.  
+To calculate the free energy as a function of phi, it is sufficient to use the following command line:
+
+\verbatim
+plumed sum_hills --hills HILLS
+\endverbatim
+
+The command above generates a file called fes.dat in which the free-energy surface as function
+of phi is calculated on a regular grid. One can modify the default name for the free energy file,
+as well as the boundaries and bin size of the grid, by using the following options of \ref sum_hills :
+
+\verbatim
+--outfile - specify the outputfile for sumhills
+--min - the lower bounds for the grid
+--max - the upper bounds for the grid
+--bin - the number of bins for the grid
+--spacing - grid spacing, alternative to the number of bins
+\endverbatim 
+
+The result should look like this:
+
+\anchor trieste-4-metad-phifes-fig
+\image html munster-metad-phifes.png "Estimate of the free energy as a function of the dihedral phi from a 10ns-long well-tempered metadynamics simulation."
+
+To assess the convergence of a metadynamics simulation, one can calculate the estimate of the free energy as a function
+of simulation time. At convergence, the reconstructed profiles should be similar.
+The option --stride should be used to give an estimate of the free energy every N Gaussians deposited, and
+the option --mintozero can be used to align the profiles by setting the global minimum to zero.
+If we use the following command line:
+
+\verbatim
+plumed sum_hills --hills HILLS --stride 100 --mintozero
+\endverbatim
+
+one free energy is calculated every 100 Gaussians deposited, and the global minimum is set to zero in all profiles.
+The resulting plot should look like the following:
+
+\anchor trieste-4-metad-phifest-fig
+\image html munster-metad-phifest.png "Estimates of the free energy as a function of the dihedral phi calculated every 100 Gaussians deposited."
+
+These two qualitative observations:
+- the system is diffusing efficiently in the collective variable space 
+- the estimated free energy does not change significantly as a function of time 
+
+suggest that the simulation is converged. 
+
+\warning The fact that the Gaussian height is decreasing to zero should not be used as a measure of convergence
+of your metadynamics simulation!
+
+\note The two observations above are necessary, but qualitative conditions for convergence.
+A quantitative analysis of convergence can be obtained by proper error analisys of the
+recontructed free energy, as explained in the last exercise
+
+\section trieste-4-ex-2 Exercize 2: playing with collective variables
+
+\section trieste-4-ex-3 Exercize 3: quantifying the error in free-energy reconstructions
+ 
+
+\section trieste-4-conclusions Conclusions
+
+In summary, in this tutorial you should have learned how to use PLUMED to:
+- Manipulate atomic coordinates.
+- Compute collective variables.
+
+
+*/
+
+link: @subpage trieste-4
+
+description: This tutorial explains how to use PLUMED to run metadynamics simulations 
+
+additional-files: trieste-4