diff --git a/src/adjmat/ClusterDiameter.cpp b/src/adjmat/ClusterDiameter.cpp
index ec26f2eeb90e42397f6fca68ed957532a1596581..8bebeab9289b82087748d6ba2ac6969087fe7d69 100644
--- a/src/adjmat/ClusterDiameter.cpp
+++ b/src/adjmat/ClusterDiameter.cpp
@@ -40,7 +40,7 @@ than 2.0 and if they are within 6.0 nm of each other.  Depth first search cluste
 between every pair of atoms that are within the largest of the clusters found is then calculated and the largest of these distances is output to a file named
 colvar.
 
-\verbatim
+\plumedfile
 # Calculate coordination numbers
 c1: COORDINATIONNUMBER SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
 # Select coordination numbers that are more than 2.0
@@ -52,7 +52,7 @@ dfs: DFSCLUSTERING MATRIX=mat
 clust1: CLUSTER_PROPERTIES CLUSTERS=dfs CLUSTER=1
 dia: CLUSTER_DIAMETER CLUSTERS=dfs CLUSTER=1
 PRINT ARG=dia FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/ClusterDistribution.cpp b/src/adjmat/ClusterDistribution.cpp
index 93f12fb04c16574535bd87bcc34e37bd91daf53c..afd7c956e523eb39e4120713f722c3d9fbf76e5e 100644
--- a/src/adjmat/ClusterDistribution.cpp
+++ b/src/adjmat/ClusterDistribution.cpp
@@ -44,7 +44,7 @@ algorithm on the corresponding graph. The number of componets in this graph that
 As discussed in \cite tribello-clustering this input was used to analyse the formation of a polycrystal of GeTe from amorphous
 GeTe.
 
-\verbatim
+\plumedfile
 q6: Q6 SPECIES=1-32768 SWITCH={GAUSSIAN D_0=5.29 R_0=0.01 D_MAX=5.3} LOWMEM
 lq6: LOCAL_Q6 SPECIES=q6 SWITCH={GAUSSIAN D_0=5.29 R_0=0.01 D_MAX=5.3} LOWMEM
 flq6: MFILTER_MORE DATA=lq6 SWITCH={GAUSSIAN D_0=0.19 R_0=0.01 D_MAX=0.2}
@@ -54,7 +54,7 @@ mat: CONTACT_MATRIX ATOMS=fcc SWITCH={GAUSSIAN D_0=3.59 R_0=0.01 D_MAX=3.6}
 dfs: DFSCLUSTERING MATRIX=mat
 nclust: CLUSTER_DISTRIBUTION CLUSTERS=dfs TRANSFORM={GAUSSIAN D_0=5.99 R_0=0.01 D_MAX=6.0} MORE_THAN={GAUSSIAN D_0=26.99 R_0=0.01 D_MAX=27}
 PRINT ARG=nclust.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/ClusterProperties.cpp b/src/adjmat/ClusterProperties.cpp
index 9a0d38e5be9d1460ff001f57992e5cd2a9507201..b1d891fe06f55fc5acdfc63ac8fd617d14c13e25 100644
--- a/src/adjmat/ClusterProperties.cpp
+++ b/src/adjmat/ClusterProperties.cpp
@@ -42,13 +42,13 @@ a graph.  This dfs action then finds the largest connected component in this gra
 numbers for the atoms in this largest connected component are then computed and this quantity is output to a colvar
 file.  The way this input can be used is described in detail in \cite tribello-clustering.
 
-\verbatim
+\plumedfile
 lq: COORDINATIONNUMBER SPECIES=1-100 SWITCH={CUBIC D_0=0.45  D_MAX=0.55} LOWMEM
 cm: CONTACT_MATRIX ATOMS=lq  SWITCH={CUBIC D_0=0.45  D_MAX=0.55}
 dfs: DFSCLUSTERING MATRIX=cm
 clust1: CLUSTER_PROPERTIES CLUSTERS=dfs CLUSTER=1 SUM
 PRINT ARG=clust1.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/ClusterSize.cpp b/src/adjmat/ClusterSize.cpp
index b0b22621d6b60a539054a095cf361a5ba0303988..4ae18bb491060052ff566ff3a686c16f9ef0dda4 100644
--- a/src/adjmat/ClusterSize.cpp
+++ b/src/adjmat/ClusterSize.cpp
@@ -39,7 +39,7 @@ The following input uses PLUMED to calculate a adjacency matrix that connects a
 than 2.0 and if they are within 6.0 nm of each other.  Depth first search clustering is used to find the connected components in this matrix and then
 the number of atoms in the largest cluster is found.  This quantity is then output to a file called colvar
 
-\verbatim
+\plumedfile
 # Calculate coordination numbers
 c1: COORDINATIONNUMBER SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
 # Select coordination numbers that are more than 2.0
@@ -51,7 +51,7 @@ dfs: DFSCLUSTERING MATRIX=mat
 clust1: CLUSTER_PROPERTIES CLUSTERS=dfs CLUSTER=1
 nat: CLUSTER_NATOMS CLUSTERS=dfs CLUSTER=1
 PRINT ARG=nat FILE=COLVAR
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/ClusterWithSurface.cpp b/src/adjmat/ClusterWithSurface.cpp
index ab472eee503791b75715cc02642baa29f6398874..1042accd051e4c0e25008f9df7fa580bc101cbfe 100644
--- a/src/adjmat/ClusterWithSurface.cpp
+++ b/src/adjmat/ClusterWithSurface.cpp
@@ -44,7 +44,7 @@ number of atoms with indices that are between 1 and 1996 and that are either in
 atoms within the the second largest cluster are then counted and this number of atoms is output to a file called size.  In addition the indices of the atoms
 that were counted are output to a file called dfs2.dat.
 
-\verbatim
+\plumedfile
 c1: COORDINATIONNUMBER SPECIES=1-1996 SWITCH={CUBIC D_0=0.34 D_MAX=0.38}
 cf: MFILTER_LESS DATA=c1 SWITCH={CUBIC D_0=13 D_MAX=13.5}
 mat: CONTACT_MATRIX ATOMS=cf SWITCH={CUBIC D_0=0.34 D_MAX=0.38}
@@ -53,7 +53,7 @@ clust2a: CLUSTER_WITHSURFACE CLUSTERS=dfs RCUT_SURF=0.3
 size2a: CLUSTER_NATOMS CLUSTERS=clust2a CLUSTER=2
 PRINT ARG=size2a FILE=size FMT=%8.4f
 OUTPUT_CLUSTER CLUSTERS=clust2a CLUSTER=2 FILE=dfs2.dat
-\endverbatim
+\endplumedfile
 
 
 */
diff --git a/src/adjmat/ContactAlignedMatrix.cpp b/src/adjmat/ContactAlignedMatrix.cpp
index e1bdfcdba0a18811235db2b8335482cfed924f32..d00542aea1abc7d94123f9f4ee2c4173bf70adef 100644
--- a/src/adjmat/ContactAlignedMatrix.cpp
+++ b/src/adjmat/ContactAlignedMatrix.cpp
@@ -57,12 +57,12 @@ is greater than 0.5.  The sum of the rows of this matrix are then computed.  The
 many of the molecules that are within the first coordination sphere of molecule \f$i\f$ have an orientation that is similar to that of
 molecule \f$i\f$.  We thus calculate the number of these "coordination numbers" that are greater than 1.0 and output this quantity to a file.
 
-\verbatim
+\plumedfile
 m1: MOLECULES MOL1=1,2 MOL2=3,4 MOL3=5,6 MOL4=7,8
 mat: ALIGNED_MATRIX ATOMS=m1 SWITCH={RATIONAL R_0=0.1} ORIENTATION_SWITCH={RATIONAL R_0=0.1 D_MAX=0.5}
 rr: ROWSUMS MATRIX=mat MORE_THAN={RATIONAL D_0=1.0 R_0=0.1}
 PRINT ARG=rr.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/ContactMatrix.cpp b/src/adjmat/ContactMatrix.cpp
index 1105a1aa230d4bf7f9437f8c74eb2cdabbc1d17f..e1a332a2af9a21f8ac2b7f703c27ea450ad9dcfc 100644
--- a/src/adjmat/ContactMatrix.cpp
+++ b/src/adjmat/ContactMatrix.cpp
@@ -48,11 +48,11 @@ The input shown below calculates a \f$6 \times 6\f$ matrix whose elements are eq
 of each other and which is zero otherwise.  The columns in this matrix are then summed so as to give the coordination number for each atom.
 The final quantity output in the colvar file is thus the average coordination number.
 
-\verbatim
+\plumedfile
 aa: CONTACT_MATRIX ATOMS=1-6 SWITCH={EXP D_0=0.2 R_0=0.1 D_MAX=0.66}
 COLUMNSUMS MATRIX=mat MEAN LABEL=csums
 PRINT ARG=csums.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/DFSClustering.cpp b/src/adjmat/DFSClustering.cpp
index af0dc327e8738a5590f72f12d68bdb404d66c738..23bf5c592caafc058baef4ba28a3bb17303c9ed2 100644
--- a/src/adjmat/DFSClustering.cpp
+++ b/src/adjmat/DFSClustering.cpp
@@ -53,13 +53,13 @@ a graph.  This dfs action then finds the largest connected component in this gra
 numbers for the atoms in this largest connected component are then computed and this quantity is output to a colvar
 file.  The way this input can be used is described in detail in \cite tribello-clustering.
 
-\verbatim
+\plumedfile
 lq: COORDINATIONNUMBER SPECIES=1-100 SWITCH={CUBIC D_0=0.45  D_MAX=0.55} LOWMEM
 cm: CONTACT_MATRIX ATOMS=lq  SWITCH={CUBIC D_0=0.45  D_MAX=0.55}
 dfs: DFSCLUSTERING MATRIX=cm
 clust1: CLUSTER_PROPERTIES CLUSTERS=dfs CLUSTER=1 SUM
 PRINT ARG=clust1.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/HbondMatrix.cpp b/src/adjmat/HbondMatrix.cpp
index f75e95754bb1697b4b1d474beb11e3f73ad72ba4..32e7989deee60972a4cad90ffe7008730109226f 100644
--- a/src/adjmat/HbondMatrix.cpp
+++ b/src/adjmat/HbondMatrix.cpp
@@ -66,13 +66,13 @@ on the number of hydrogen bonds each of the water molecules donates and accepts.
 five columns of data.  The first four of these columns are a label for the atom and the x, y and z position of the oxygen.  The last column is then
 the number of accepted/donated hydrogen bonds.
 
-\verbatim
+\plumedfile
 mat: HBOND_MATRIX ATOMS=1-192:3 HYDROGENS=2-192:3,3-192:3 SWITCH={RATIONAL R_0=3.20} HSWITCH={RATIONAL R_0=2.30} ASWITCH={RATIONAL R_0=0.167pi} SUM
 rsums: ROWSUMS MATRIX=mat MEAN
 csums: COLUMNSUMS MATRIX=mat MEAN
 DUMPMULTICOLVAR DATA=rsums FILE=donors.xyz
 DUMPMULTICOLVAR DATA=csums FILE=acceptors.x
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/MatrixColumnSums.cpp b/src/adjmat/MatrixColumnSums.cpp
index 7735023490ce6b93abc91d88cc1cfb47f75fd048..cae24fcec899d62d88cc879ead801480d41a35a5 100644
--- a/src/adjmat/MatrixColumnSums.cpp
+++ b/src/adjmat/MatrixColumnSums.cpp
@@ -43,21 +43,21 @@ tells you whether atoms \f$i\f$ and \f$j\f$ are within 1.0 nm of each other.  Th
 and the average value is computed.  As such the following input provides an alternative method for calculating the coordination numbers
 of atoms 1 to 10.
 
-\verbatim
+\plumedfile
 mat: CONTACT_MATRIX ATOMS=1-10 SWITCH={RATIONAL R_0=1.0}
 rsums: COLUMNSUMS MATRIX=mat MEAN
 PRINT ARG=rsums.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following input demonstrates another way that an average coordination number can be computed.  This input calculates the number of atoms
 with indices between 1 and 5 that are within the first coordination spheres of each of the atoms within indices between 6 and 15.  The average
 coordination number is then calculated from these fifteen coordination numbers and this quantity is output to a file.
 
-\verbatim
+\plumedfile
 mat2: CONTACT_MATRIX ATOMSA=1-5 ATOMSB=6-15 SWITCH={RATIONAL R_0=1.0}
 rsums: COLUMNSUMS MATRIX=mat2 MEAN
 PRINT ARG=rsums.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/MatrixRowSums.cpp b/src/adjmat/MatrixRowSums.cpp
index 038aa89a71fe36d25487433f621811ec705c3200..6ef15fb282c0ddabd29136e1a6ed55ba668e1cfd 100644
--- a/src/adjmat/MatrixRowSums.cpp
+++ b/src/adjmat/MatrixRowSums.cpp
@@ -43,21 +43,21 @@ tells you whether atoms \f$i\f$ and \f$j\f$ are within 1.0 nm of each other.  Th
 and the average value is computed.  As such the following input provides an alternative method for calculating the coordination numbers
 of atoms 1 to 10.
 
-\verbatim
+\plumedfile
 mat: CONTACT_MATRIX ATOMS=1-10 SWITCH={RATIONAL R_0=1.0}
 rsums: ROWSUMS MATRIX=mat MEAN
 PRINT ARG=rsums.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following input demonstrates another way that an average coordination number can be computed.  This input calculates the number of atoms
 with indices between 6 and 15 that are within the first coordination spheres of each of the atoms within indices between 1 and 5.  The average
 coordination number is then calculated from these five coordination numbers and this quantity is output to a file.
 
-\verbatim
+\plumedfile
 mat2: CONTACT_MATRIX ATOMSA=1-5 ATOMSB=6-15 SWITCH={RATIONAL R_0=1.0}
 rsums: ROWSUMS MATRIX=mat2 MEAN
 PRINT ARG=rsums.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/OutputCluster.cpp b/src/adjmat/OutputCluster.cpp
index c2b752f0d0821c70254e7da7e90f57588649f0cb..4dd430f3c71befba4c75b6a443e61b1698ab1ddc 100644
--- a/src/adjmat/OutputCluster.cpp
+++ b/src/adjmat/OutputCluster.cpp
@@ -47,13 +47,13 @@ that satisfy this criteria.  The DFS algorithm is then used to find the connecte
 in this matrix and the indices of the atoms in the largest connected component are then output
 to a file.
 
-\verbatim
+\plumedfile
 c1: COORDINATIONNUMBER SPECIES=1-1996 SWITCH={CUBIC D_0=0.34 D_MAX=0.38}
 cf: MFILTER_LESS DATA=c1 SWITCH={CUBIC D_0=13 D_MAX=13.5}
 mat: CONTACT_MATRIX ATOMS=cf SWITCH={CUBIC D_0=0.34 D_MAX=0.38}
 dfs: DFSCLUSTERING MATRIX=mat
 OUTPUT_CLUSTER CLUSTERS=dfs CLUSTER=1 FILE=dfs.dat
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/SMACMatrix.cpp b/src/adjmat/SMACMatrix.cpp
index ac433076bd9abef860031586015ac4810af0869b..a4e8288f49f5e12b0dfd07db3bb77eaeb1fe8882 100644
--- a/src/adjmat/SMACMatrix.cpp
+++ b/src/adjmat/SMACMatrix.cpp
@@ -49,7 +49,7 @@ molecules \f$i\f$ and \f$j\f$ are within 6 angstroms of each other and if the to
 of these molecules is close to 0 or \f$\pi\f$.  The various connected components of this matrix are determined using the
 \ref DFSCLUSTERING algorithm and then the size of the largest cluster of connectes molecules is output to a colvar file
 
-\verbatim
+\plumedfile
 UNITS LENGTH=A
 
 MOLECULES ...
@@ -70,7 +70,7 @@ SMAC_MATRIX ...
 dfs1: DFSCLUSTERING MATRIX=smacm
 cc2: CLUSTER_NATOMS CLUSTERS=dfs1 CLUSTER=1
 PRINT ARG=smac.*,cc1.*,cc2 FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/adjmat/Sprint.cpp b/src/adjmat/Sprint.cpp
index f7e0a0da8f5ff58c6fe3b2950c4e3156f6ac726e..3e90100e8e2d6c139b559a220e0edcb3d04e996e 100644
--- a/src/adjmat/Sprint.cpp
+++ b/src/adjmat/Sprint.cpp
@@ -46,17 +46,17 @@ This example input calculates the 7 SPRINT coordinates for a 7 atom cluster of L
 atoms and prints their values to a file.  In this input the SPRINT coordinates are calculated
 in the manner described in ?? so two atoms are adjacent if they are within a cutoff:
 
-\verbatim
+\plumedfile
 DENSITY SPECIES=1-7 LABEL=d1
 CONTACT_MATRIX ATOMS=d1 SWITCH={RATIONAL R_0=0.1} LABEL=mat
 SPRINT MATRIX=mat LABEL=ss
 PRINT ARG=ss.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 This example input calculates the 14 SPRINT coordinates foa a molecule composed of 7 hydrogen and
 7 carbon atoms.  Once again two atoms are adjacent if they are within a cutoff:
 
-\verbatim
+\plumedfile
 DENSITY SPECIES=1-7 LABEL=c
 DENSITY SPECIES=8-14 LABEL=h
 
@@ -71,7 +71,7 @@ CONTACT_MATRIX ...
 SPRINT MATRIX=mat LABEL=ss
 
 PRINT ARG=ss.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/analysis/Average.cpp b/src/analysis/Average.cpp
index c96b2556bbe11df0fd4ae0b4dcd7ac080918d8ca..9d0566da8e2a35c13cecdf769dfabab45e009741 100644
--- a/src/analysis/Average.cpp
+++ b/src/analysis/Average.cpp
@@ -51,11 +51,11 @@ and output this to a file called COLVAR.  In this example it is assumed that no
 on the system and that the weights, \f$w(t')\f$ in the formulae above can thus all be set equal
 to one.
 
-\verbatim
+\plumedfile
 d1: DISTANCE ATOMS=1,2
 d1a: AVERAGE ARG=d1
 PRINT ARG=d1a FILE=colvar STRIDE=100
-\endverbatim
+\endplumedfile
 
 The following example calculates the ensemble average for the torsional angle involving atoms 1, 2, 3 and 4.
 At variance with the previous example this quantity is periodic so the second formula in the above introduction
@@ -65,23 +65,23 @@ forgotten and the process of averaging is begun again.  The quantities output in
 block averages taken over the first 100 frames of the trajectory, the block average over the second 100 frames
 of trajectory and so on.
 
-\verbatim
+\plumedfile
 t1: TORSION ATOMS=1,2,3,4
 t1a: AVERAGE ARG=t1 CLEAR=100
 PRINT ARG=t1a FILE=colvar STRIDE=100
-\endverbatim
+\endplumedfile
 
 This third example incorporates a bias.  Notice that the effect the bias has on the ensemble average is removed by taking
 advantage of the \ref REWEIGHT_BIAS method.  The final ensemble averages output to the file are thus block ensemble averages for the
 unbiased canononical ensemble at a temperature of 300 K.
 
-\verbatim
+\plumedfile
 t1: TORSION ATOMS=1,2,3,4
 RESTRAINT ARG=t1 AT=pi KAPPA=100.
 ww: REWEIGHT_BIAS TEMP=300
 t1a: AVERAGE ARG=t1 LOGWEIGHTS=ww CLEAR=100
 PRINT ARG=t1a FILE=colvar STRIDE=100
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/analysis/ClassicalMultiDimensionalScaling.cpp b/src/analysis/ClassicalMultiDimensionalScaling.cpp
index 93bb25f1d485c7a82d1cbac211952c8678801e23..68ac6aeae304c4b06b38f5878d816ce4100fe3bf 100644
--- a/src/analysis/ClassicalMultiDimensionalScaling.cpp
+++ b/src/analysis/ClassicalMultiDimensionalScaling.cpp
@@ -55,7 +55,7 @@ can be found in the tutorial \ref belfast-3 and in the following <a href="https:
 The following command instructs plumed to construct a classical multidimensional scaling projection of a trajectory.
 The RMSD distance between atoms 1-256 have moved is used to measure the distances in the high-dimensional space.
 
-\verbatim
+\plumedfile
 CLASSICAL_MDS ...
   ATOMS=1-256
   METRIC=OPTIMAL-FAST
@@ -63,7 +63,7 @@ CLASSICAL_MDS ...
   NLOW_DIM=2
   OUTPUT_FILE=rmsd-embed
 ... CLASSICAL_MDS
-\endverbatim
+\endplumedfile
 
 The following section is for people who are interested in how this method works in detail. A solid understanding of this material is
 not necessary to use MDS.
diff --git a/src/analysis/Commit.cpp b/src/analysis/Commit.cpp
index 5cbacf4db1d10ef29c3d5ea8f2c45b612eaab869..d51d6b9e55f12af1600b6da880657438f589c341 100644
--- a/src/analysis/Commit.cpp
+++ b/src/analysis/Commit.cpp
@@ -37,7 +37,7 @@ The following input monitors two torsional angles during a simulation,
 defines two basins (A and B) as a function of the two torsions and
 stops the simulation when it falls in one of the two. In the log
 file will be shown the latest values for the CVs and the basin reached.
-\verbatim
+\plumedfile
 TORSION ATOMS=1,2,3,4 LABEL=r1
 TORSION ATOMS=2,3,4,5 LABEL=r2
 COMMITTOR ...
@@ -48,7 +48,7 @@ COMMITTOR ...
   BASIN_LL2=-0.15,-0.20
   BASIN_UL2=-0.25,-0.40
 ... COMMITTOR
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/analysis/Histogram.cpp b/src/analysis/Histogram.cpp
index b8a74c797d8197d4278653f5b77e76f562bc818d..16abb4a6dcb9ffa1791edd7a7721c8ca6dda1d3b 100644
--- a/src/analysis/Histogram.cpp
+++ b/src/analysis/Histogram.cpp
@@ -84,7 +84,7 @@ Additional material and examples can be also found in the tutorial \ref belfast-
 
 The following input monitors two torsional angles during a simulation
 and outputs a continuos histogram as a function of them at the end of the simulation.
-\verbatim
+\plumedfile
 TORSION ATOMS=1,2,3,4 LABEL=r1
 TORSION ATOMS=2,3,4,5 LABEL=r2
 HISTOGRAM ...
@@ -97,11 +97,11 @@ HISTOGRAM ...
 ... HISTOGRAM
 
 DUMPGRID GRID=hh FILE=histo
-\endverbatim
+\endplumedfile
 
 The following input monitors two torsional angles during a simulation
 and outputs a discrete histogram as a function of them at the end of the simulation.
-\verbatim
+\plumedfile
 TORSION ATOMS=1,2,3,4 LABEL=r1
 TORSION ATOMS=2,3,4,5 LABEL=r2
 HISTOGRAM ...
@@ -115,11 +115,11 @@ HISTOGRAM ...
 ... HISTOGRAM
 
 DUMPGRID GRID=hh FILE=histo
-\endverbatim
+\endplumedfile
 
 The following input monitors two torsional angles during a simulation
 and outputs the histogram accumulated thus far every 100000 steps.
-\verbatim
+\plumedfile
 TORSION ATOMS=1,2,3,4 LABEL=r1
 TORSION ATOMS=2,3,4,5 LABEL=r2
 HISTOGRAM ...
@@ -132,14 +132,14 @@ HISTOGRAM ...
 ... HISTOGRAM
 
 DUMPGRID GRID=hh FILE=histo STRIDE=100000
-\endverbatim
+\endplumedfile
 
 The following input monitors two torsional angles during a simulation
 and outputs a separate histogram for each 100000 steps worth of trajectory.
 Notice how the CLEAR keyword is used here and how it is not used in the
 previous example.
 
-\verbatim
+\plumedfile
 TORSION ATOMS=1,2,3,4 LABEL=r1
 TORSION ATOMS=2,3,4,5 LABEL=r2
 HISTOGRAM ...
@@ -153,7 +153,7 @@ HISTOGRAM ...
 ... HISTOGRAM
 
 DUMPGRID GRID=hh FILE=histo STRIDE=100000
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/analysis/PCA.cpp b/src/analysis/PCA.cpp
index 5b31a04faafd648f6280c5c1dc173e136dcfc82f..528c6290b73b8648e87d0e73a852db58e2a400eb 100644
--- a/src/analysis/PCA.cpp
+++ b/src/analysis/PCA.cpp
@@ -63,9 +63,9 @@ of the first 22 atoms.  The TYPE=OPTIMAL instruction ensures that translational
 The first two principal components will be output to a file called pca-comp.pdb.  Trajectory frames will be collected on every step and the PCA calculation
 will be performed at the end of the simulation.
 
-\verbatim
+\plumedfile
 PCA METRIC=OPTIMAL ATOMS=1-22 STRIDE=1 USE_ALL_DATA NLOW_DIM=2 OFILE=pca-comp.pdb
-\endverbatim
+\endplumedfile
 
 The following input instructs PLUMED to perform a principal component analysis in which the covariance matrix is calculated from chnages in the six distances
 seen in the previous lines.  Notice that here the TYPE=EUCLIDEAN keyword is used to indicate that no alighment has to be done when calculating the various
@@ -75,7 +75,7 @@ PCA analysis will be performed twice.  The REWEIGHT_BIAS keyword in this input t
 when calculating averages and covariances a reweighting should be performed based and each frames' weight in these calculations should be determined based on
 the current value of the instantaneous bias (see \ref REWEIGHT_BIAS).
 
-\verbatim
+\plumedfile
 d1: DISTANCE ATOMS=1,2
 d2: DISTANCE ATOMS=1,3
 d3: DISTANCE ATOMS=1,4
@@ -84,7 +84,7 @@ d5: DISTANCE ATOMS=2,4
 d6: DISTANCE ATOMS=3,4
 
 PCA ARG=d1,d2,d3,d4,d5,d6 METRIC=EUCLIDEAN STRIDE=5 RUN=1000 NLOW_DIM=2 REWEIGHT_BIAS OFILE=pca-comp.pdb
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/bias/ABMD.cpp b/src/bias/ABMD.cpp
index 5e8bef344e96a30c09e23d149f31709a06e9bc62..a76c5341c3b983cde34534ea31f8c57d18949113 100644
--- a/src/bias/ABMD.cpp
+++ b/src/bias/ABMD.cpp
@@ -70,13 +70,12 @@ an additional white noise acting on the minimum position of the bias.
 The following input sets up two biases, one on the distance between atoms 3 and 5
 and another on the distance between atoms 2 and 4. The two target values are defined
 using TO and the two strength using KAPPA. The total energy of the bias is printed.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 ABMD ARG=d1,d2 TO=1.0,1.5 KAPPA=5.0,5.0 LABEL=abmd
 PRINT ARG=abmd.bias,abmd.d1_min,abmd.d2_min
-\endverbatim
-(See also \ref DISTANCE and \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/bias/BiasValue.cpp b/src/bias/BiasValue.cpp
index 3e77621a55546b1510f126fbdd20e40f935efbfa..ff741a6c0c5c3e0e7b9b265fe11664d6442a7d10 100644
--- a/src/bias/BiasValue.cpp
+++ b/src/bias/BiasValue.cpp
@@ -42,18 +42,17 @@ to some collective variable then using the value of this function directly as a
 The following input tells plumed to use the value of the distance between atoms 3 and 5
 and the value of the distance between atoms 2 and 4 as biases.
 It then tells plumed to print the energy of the restraint
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=3,6 LABEL=d2
 BIASVALUE ARG=d1,d2 LABEL=b
 PRINT ARG=d1,d2,b.d1,b.d2
-\endverbatim
-(See also \ref DISTANCE and \ref PRINT).
+\endplumedfile
 
 Another thing one can do is asking one system to follow
 a circle in sin/cos according a  time dependence
 
-\verbatim
+\plumedfile
 t: TIME
 # this just print cos and sin of time
 cos: MATHEVAL ARG=t VAR=t FUNC=cos(t) PERIODIC=NO
@@ -72,9 +71,7 @@ vv1:  MATHEVAL ARG=mycos,mysin,cos,sin VAR=mc,ms,c,s  FUNC=100*((mc-c)^2+(ms-s)^
 cc: BIASVALUE ARG=vv1
 # some printout
 PRINT ARG=t,cos,sin,d.x,d.y,d.z,mycos,mysin,cc.bias.vv1 STRIDE=1 FILE=colvar FMT=%8.4f
-\endverbatim
-(see also \ref TIME, \ref MATHEVAL, \ref COM, \ref DISTANCE, and \ref PRINT).
-
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/bias/ExtendedLagrangian.cpp b/src/bias/ExtendedLagrangian.cpp
index 85de570dd5030e1c58e3a513d7747e7da5a2e2ed..ee3d3ada041452f28b612d2295b363c18ec80eca 100644
--- a/src/bias/ExtendedLagrangian.cpp
+++ b/src/bias/ExtendedLagrangian.cpp
@@ -88,28 +88,26 @@ and many shorter runs.
 
 The following input tells plumed to perform a metadynamics
 with an extended Lagrangian on two torsional angles.
-\verbatim
+\plumedfile
 phi: TORSION ATOMS=5,7,9,15
 psi: TORSION ATOMS=7,9,15,17
 ex: EXTENDED_LAGRANGIAN ARG=phi,psi KAPPA=20,20.0 TAU=0.1,0.1
 METAD ARG=ex.phi_fict,ex.psi_fict PACE=100 SIGMA=0.35,0.35 HEIGHT=0.1
 # monitor the two variables
 PRINT STRIDE=10 ARG=phi,psi,ex.phi_fict,ex.psi_fict FILE=COLVAR
-\endverbatim
-(See also \ref TORSION, \ref METAD, and \ref PRINT).
+\endplumedfile
 
 The following input tells plumed to perform a TAMD (or dAFED)
 calculation on two torsional angles, keeping the two variables
 at a fictitious temperature of 3000K with a Langevin thermostat
 with friction 10
-\verbatim
+\plumedfile
 phi: TORSION ATOMS=5,7,9,15
 psi: TORSION ATOMS=7,9,15,17
 ex: EXTENDED_LAGRANGIAN ARG=phi,psi KAPPA=20,20.0 TAU=0.1,0.1 FRICTION=10,10 TEMP=3000
 # monitor the two variables
 PRINT STRIDE=10 ARG=phi,psi,ex.phi_fict,ex.psi_fict FILE=COLVAR
-\endverbatim
-(See also \ref TORSION and \ref PRINT)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/bias/External.cpp b/src/bias/External.cpp
index 62bc54cf9252fc9102014483546b697d5b765816..7c714efae64d1c6568d29efa0e70b36fd7721353 100644
--- a/src/bias/External.cpp
+++ b/src/bias/External.cpp
@@ -39,11 +39,10 @@ Calculate a restraint that is defined on a grid that is read during start up
 \par Examples
 The following is an input for a calculation with an external potential that is
 defined in the file bias.dat and that acts on the distance between atoms 3 and 5.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 EXTERNAL ARG=d1 FILE=bias.dat LABEL=external
-\endverbatim
-(See also \ref DISTANCE \ref PRINT).
+\endplumedfile
 
 The header in the file bias.dat should read:
 \verbatim
@@ -61,11 +60,11 @@ with NOSPLINE you do not need to provide derivative information.
 You can also include grids that are a function of more than one collective
 variable.  For instance the following would be the input for an external
 potential acting on two torsional angles:
-\verbatim
+\plumedfile
 TORSION ATOMS=4,5,6,7 LABEL=t1
 TORSION ATOMS=6,7,8,9 LABEL=t2
 EXTERNAL ARG=t1,t2 FILE=bias.dat LABEL=ext
-\endverbatim
+\endplumedfile
 
 The header in the file bias.dat for this calculation would read:
 \verbatim
diff --git a/src/bias/LWalls.cpp b/src/bias/LWalls.cpp
index 33d7f7c04604a5a0217e918f4c200ff5b81e0d4b..26bb720f13ba11b7417230db587402eb8bc17e26 100644
--- a/src/bias/LWalls.cpp
+++ b/src/bias/LWalls.cpp
@@ -52,13 +52,13 @@ The following input tells plumed to add both a lower and an upper walls on the d
 between atoms 3 and 5 and the distance between atoms 2 and 4. The lower and upper limits
 are defined at different values. The strength of the walls is the same for the four cases.
 It also tells plumed to print the energy of the walls.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 UPPER_WALLS ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=uwall
 LOWER_WALLS ARG=d1,d2 AT=0.0,1.0 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=lwall
 PRINT ARG=uwall.bias,lwall.bias
-\endverbatim
+\endplumedfile
 (See also \ref DISTANCE and \ref PRINT).
 
 */
diff --git a/src/bias/MetaD.cpp b/src/bias/MetaD.cpp
index b586185705abc526b20c18ce567923982896535f..5023c38f2cfec80681d7ad36bb731d1b96607e6b 100644
--- a/src/bias/MetaD.cpp
+++ b/src/bias/MetaD.cpp
@@ -163,35 +163,35 @@ The following input is for a standard metadynamics calculation using as
 collective variables the distance between atoms 3 and 5
 and the distance between atoms 2 and 4. The value of the CVs and
 the metadynamics bias potential are written to the COLVAR file every 100 steps.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 METAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=restraint
 PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR
-\endverbatim
+\endplumedfile
 (See also \ref DISTANCE \ref PRINT).
 
 \par
 If you use adaptive Gaussians, with diffusion scheme where you use
 a Gaussian that should cover the space of 20 timesteps in collective variables.
 Note that in this case the histogram correction is needed when summing up hills.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 METAD ARG=d1,d2 SIGMA=20 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=DIFF
 PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR
-\endverbatim
+\endplumedfile
 
 \par
 If you use adaptive Gaussians, with geometrical scheme where you use
 a Gaussian that should cover the space of 0.05 nm in Cartesian space.
 Note that in this case the histogram correction is needed when summing up hills.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 METAD ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
 PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR
-\endverbatim
+\endplumedfile
 
 \par
 When using adaptive Gaussians you might want to limit how the hills width can change.
@@ -199,7 +199,7 @@ You can use SIGMA_MIN and SIGMA_MAX keywords.
 The sigmas should specified in terms of CV so you should use the CV units.
 Note that if you use a negative number, this means that the limit is not set.
 Note also that in this case the histogram correction is needed when summing up hills.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 METAD ...
@@ -207,7 +207,7 @@ METAD ...
   SIGMA_MIN=0.2,0.1 SIGMA_MAX=0.5,1.0
 ... METAD
 PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR
-\endverbatim
+\endplumedfile
 
 \par
 Multiple walkers can be also use as in  \cite multiplewalkers
@@ -215,7 +215,7 @@ These are enabled by setting the number of walker used, the id of the
 current walker which interprets the input file, the directory where the
 hills containing files resides, and the frequency to read the other walkers.
 Here is an example
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 METAD ...
    ARG=d1 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint
@@ -224,7 +224,7 @@ METAD ...
    WALKERS_DIR=../
    WALKERS_RSTRIDE=100
 ... METAD
-\endverbatim
+\endplumedfile
 where  WALKERS_N is the total number of walkers, WALKERS_ID is the
 id of the present walker (starting from 0 ) and the WALKERS_DIR is the directory
 where all the walkers are located. WALKERS_RSTRIDE is the number of step between
@@ -237,7 +237,7 @@ presented in \cite Tiwary_jp504920s as described above.
 This is enabled by using the keyword REWEIGHTING_NGRID where the grid used for
 the calculation is set. The number of grid points given in REWEIGHTING_NGRID
 should be equal or larger than the number of grid points given in GRID_BIN.
-\verbatim
+\plumedfile
 METAD ...
  LABEL=metad
  ARG=phi,psi SIGMA=0.20,0.20 HEIGHT=1.20 BIASFACTOR=5 TEMP=300.0 PACE=500
@@ -245,7 +245,7 @@ METAD ...
  REWEIGHTING_NGRID=150,150
  REWEIGHTING_NHILLS=20
 ... METAD
-\endverbatim
+\endplumedfile
 Here we have asked that the calculation is performed every 20 hills by using
 REWEIGHTING_NHILLS keyword. If this keyword is not given the calculation will
 by default be performed every 50 hills. The c(t) reweighting factor will be given
@@ -285,7 +285,7 @@ Alternatively, if you use a BIASFACTOR yout simulation will converge to a free
 energy that is a linear combination of the target free energy and of the intrinsic free energy
 determined by the original force field.
 
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 METAD ...
  LABEL=t1
@@ -295,7 +295,7 @@ METAD ...
 ... METAD
 
 PRINT ARG=d1,t1.bias STRIDE=100 FILE=COLVAR
-\endverbatim
+\endplumedfile
 
 The header in the file dist.dat for this calculation would read:
 
@@ -309,19 +309,19 @@ The header in the file dist.dat for this calculation would read:
 
 Notice that BIASFACTOR can also be chosen as equal to 1. In this case one will perform
 unbiased sampling. Instead of using HEIGHT, one should provide the TAU parameter.
-\verbatim
+\plumedfile
 d: DISTANCE ATOMS=3,5
 METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 BIASFACTOR=1.0
-\endverbatim
+\endplumedfile
 The HILLS file obtained will still work with `plumed sum_hills` so as to plot a free-energy.
 The case where this makes sense is probably that of RECT simulations.
 
 Regarding RECT simulations, you can also use the RECT keyword so as to avoid using multiple input files.
 For instance, a single input file will be
-\verbatim
+\plumedfile
 d: DISTANCE ATOMS=3,5
 METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 RECT=1.0,1.5,2.0,3.0
-\endverbatim
+\endplumedfile
 The number of elements in the RECT array should be equal to the number of replicas.
 
 
diff --git a/src/bias/Metainference.cpp b/src/bias/Metainference.cpp
index e4d3e8b94f0d63194082568859b63d1beafc5dd4..3953c436d1368922505817caa483ff033bc6ccb2 100644
--- a/src/bias/Metainference.cpp
+++ b/src/bias/Metainference.cpp
@@ -58,7 +58,7 @@ them and comparing them with a set of experimental values. RDCs are compared wit
 the experimental data but for a multiplication factor SCALE that is also sampled by
 MC on-the-fly
 
-\verbatim
+\plumedfile
 RDC ...
 LABEL=rdc
 SCALE=0.0001
@@ -83,12 +83,12 @@ LABEL=spe
 ... METAINFERENCE
 
 PRINT ARG=spe.bias FILE=BIAS STRIDE=1
-\endverbatim
+\endplumedfile
 
 in the following example instead of using one uncertainty parameter per data point we use
 a single uncertainty value in a long-tailed gaussian to take into account for outliers.
 
-\verbatim
+\plumedfile
 METAINFERENCE ...
 ARG=ardc.*
 NOISETYPE=OUTLIERS
@@ -99,9 +99,7 @@ SIGMA_MEAN=0.001
 TEMP=300
 LABEL=spe
 ... METAINFERENCE
-\endverbatim
-
-(See also \ref RDC and \ref ENSEMBLE).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/bias/MovingRestraint.cpp b/src/bias/MovingRestraint.cpp
index 13cafce97e63a71aa0a76718229986d4095963b4..d5317af1d9abb8a0e4e0771eede4b68e49e3681e 100644
--- a/src/bias/MovingRestraint.cpp
+++ b/src/bias/MovingRestraint.cpp
@@ -53,7 +53,7 @@ Additional material and examples can be also found in the tutorial \ref belfast-
 The following input is dragging the distance between atoms 2 and 4
 from 1 to 2 in the first 1000 steps, then back in the next 1000 steps.
 In the following 500 steps the restraint is progressively switched off.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=2,4 LABEL=d
 MOVINGRESTRAINT ...
   ARG=d
@@ -62,12 +62,12 @@ MOVINGRESTRAINT ...
   STEP2=2000 AT2=1.0
   STEP3=2500         KAPPA3=0.0
 ... MOVINGRESTRAINT
-\endverbatim
+\endplumedfile
 The following input is progressively building restraints
 distances between atoms 1 and 5 and between atoms 2 and 4
 in the first 1000 steps. Afterwards, the restraint is kept
 static.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=1,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 MOVINGRESTRAINT ...
@@ -75,10 +75,10 @@ MOVINGRESTRAINT ...
   STEP0=0    AT0=1.0,1.5 KAPPA0=0.0,0.0
   STEP1=1000 AT1=1.0,1.5 KAPPA1=1.0,1.0
 ... MOVINGRESTRAINT
-\endverbatim
+\endplumedfile
 The following input is progressively bringing atoms 1 and 2
 close to each other with an upper wall
-\verbatim
+\plumedfile
 DISTANCE ATOMS=1,2 LABEL=d1
 MOVINGRESTRAINT ...
   ARG=d1
@@ -86,7 +86,7 @@ MOVINGRESTRAINT ...
   STEP0=0    AT0=1.0 KAPPA0=10.0
   STEP1=1000 AT1=0.0
 ... MOVINGRESTRAINT
-\endverbatim
+\endplumedfile
 
 By default the Action is issuing some values which are
 the work on each degree of freedom, the center of the harmonic potential,
diff --git a/src/bias/PBMetaD.cpp b/src/bias/PBMetaD.cpp
index b50557f7d297f59f61da770a9ce711a77415f9ed..76c35154bd6b648132fb5396ca88fec4f694edea 100644
--- a/src/bias/PBMetaD.cpp
+++ b/src/bias/PBMetaD.cpp
@@ -150,18 +150,18 @@ The following input is for PBMetaD calculation using as
 collective variables the distance between atoms 3 and 5
 and the distance between atoms 2 and 4. The value of the CVs and
 the PBMetaD bias potential are written to the COLVAR file every 100 steps.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 PBMETAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=pb FILE=HILLS_d1,HILLS_d2
 PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
-\endverbatim
+\endplumedfile
 (See also \ref DISTANCE and \ref PRINT).
 
 \par
 If you use well-tempered metadynamics, you should specify a single biasfactor and initial
 Gaussian height.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 PBMETAD ...
@@ -170,11 +170,11 @@ PACE=500 BIASFACTOR=8 LABEL=pb
 FILE=HILLS_d1,HILLS_d2
 ... PBMETAD
 PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
-\endverbatim
+\endplumedfile
 
 \par
 The following input enables the MPI version of multiple-walkers.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 PBMETAD ...
@@ -184,7 +184,7 @@ FILE=HILLS_d1,HILLS_d2
 WALKERS_MPI
 ... PBMETAD
 PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
-\endverbatim
+\endplumedfile
 
 \par
 The disk version of multiple-walkers can be
@@ -192,7 +192,7 @@ enabled by setting the number of walker used, the id of the
 current walker which interprets the input file, the directory where the
 hills containing files resides, and the frequency to read the other walkers.
 Here is an example
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 PBMETAD ...
@@ -205,7 +205,7 @@ WALKERS_DIR=../
 WALKERS_RSTRIDE=100
 ... PBMETAD
 PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
-\endverbatim
+\endplumedfile
 where  WALKERS_N is the total number of walkers, WALKERS_ID is the
 id of the present walker (starting from 0 ) and the WALKERS_DIR is the directory
 where all the walkers are located. WALKERS_RSTRIDE is the number of step between
diff --git a/src/bias/Restraint.cpp b/src/bias/Restraint.cpp
index 8bb3df2027989545444b3157fceaf50333034ca0..85208d35ba1ad62377bcbf09bc662fb088261f1e 100644
--- a/src/bias/Restraint.cpp
+++ b/src/bias/Restraint.cpp
@@ -49,13 +49,12 @@ Additional material and examples can be also found in the tutorial \ref belfast-
 The following input tells plumed to restrain the distance between atoms 3 and 5
 and the distance between atoms 2 and 4, at different equilibrium
 values, and to print the energy of the restraint
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 RESTRAINT ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 LABEL=restraint
 PRINT ARG=restraint.bias
-\endverbatim
-(See also \ref DISTANCE and \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/bias/ReweightBias.cpp b/src/bias/ReweightBias.cpp
index d0b2d90a3e64af196838222f2656d920858fc278..60976e644d9d47bd2e931b9faddeaf68568eb049 100644
--- a/src/bias/ReweightBias.cpp
+++ b/src/bias/ReweightBias.cpp
@@ -43,7 +43,7 @@ restraint will have an effect on the region of phase space that will be sampled
 run using this variable.  Consequently, when the histogram as a function of the distance, \f$x\f$, is accumulated,
 we use reweighting into order to discount the effect of the bias from our final histogram.
 
-\verbatim
+\plumedfile
 x: DISTANCE ATOMS=1,2
 RESTRAINT ARG=x SLOPE=1.0 AT=0.0
 as: REWEIGHT_BIAS TEMP=300
@@ -59,7 +59,7 @@ HISTOGRAM ...
 ... HISTOGRAM
 
 DUMPGRID GRID=hB FILE=histoB STRIDE=1 FMT=%8.4f
-\endverbatim
+\endplumedfile
 
 
 */
diff --git a/src/bias/ReweightMetad.cpp b/src/bias/ReweightMetad.cpp
index 766ed40dfcd7b604dff26289fdf25bde58011875..42a4855c24d81e1b51232056c9682c836a05193b 100644
--- a/src/bias/ReweightMetad.cpp
+++ b/src/bias/ReweightMetad.cpp
@@ -39,7 +39,7 @@ we use reweighting into order to discount the effect of the bias from our final
 metadynamics instead.  Notice also that we have to specify how often we would like to calculate the c(t) reweighting
 factor and the grid over which we calculate c(t) in the input to the METAD command.
 
-\verbatim
+\plumedfile
 a: ANGLE ATOMS=1,2,3
 x: DISTANCE ATOMS=1,2
 METAD ARG=x PACE=100 SIGMA=0.1 HEIGHT=1.5 BIASFACTOR=5 GRID_MIN=0 GRID_MAX=10 GRID_BIN=100 REWEIGHTING_NGRID=100 REWEIGHTING_NHILLS=50
@@ -57,8 +57,7 @@ HISTOGRAM ...
 ... HISTOGRAM
 
 DUMPGRID GRID=hB FILE=histoB STRIDE=1 FMT=%8.4f
-
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/bias/UWalls.cpp b/src/bias/UWalls.cpp
index 0f096f11b1c8fa109252118ab5241c90a04ded4b..bfe4d8929ef2f3ae245594d70450a80517861f16 100644
--- a/src/bias/UWalls.cpp
+++ b/src/bias/UWalls.cpp
@@ -52,14 +52,13 @@ The following input tells plumed to add both a lower and an upper walls on the d
 between atoms 3 and 5 and the distance between atoms 2 and 4. The lower and upper limits
 are defined at different values. The strength of the walls is the same for the four cases.
 It also tells plumed to print the energy of the walls.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=3,5 LABEL=d1
 DISTANCE ATOMS=2,4 LABEL=d2
 UPPER_WALLS ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=uwall
 LOWER_WALLS ARG=d1,d2 AT=0.0,1.0 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=lwall
 PRINT ARG=uwall.bias,lwall.bias
-\endverbatim
-(See also \ref DISTANCE and \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/Angle.cpp b/src/colvar/Angle.cpp
index ed57d03ea56966effc65d70d6153d2a9d56f820e..285d6fa5050eca25969e490cbc96ab4199d34f2b 100644
--- a/src/colvar/Angle.cpp
+++ b/src/colvar/Angle.cpp
@@ -67,7 +67,7 @@ This command tells plumed to calculate the angle between the vector connecting a
 the vector connecting atom 2 to atom 3 and to print it on file COLVAR1. At the same time,
 the angle between vector connecting atom 1 to atom 2 and the vector connecting atom 3 to atom 4 is printed
 on file COLVAR2.
-\verbatim
+\plumedfile
 
 a: ANGLE ATOMS=1,2,3
 # equivalently one could state:
@@ -77,8 +77,7 @@ b: ANGLE ATOMS=1,2,3,4
 
 PRINT ARG=a FILE=COLVAR1
 PRINT ARG=b FILE=COLVAR2
-\endverbatim
-(see also \ref PRINT)
+\endplumedfile
 
 
 */
diff --git a/src/colvar/CS2Backbone.cpp b/src/colvar/CS2Backbone.cpp
index 040e0dfad8c6fdfff0da8c51a9254812f8dcceba..6c606e4e5caef102592328cb1fae021f21afff09 100644
--- a/src/colvar/CS2Backbone.cpp
+++ b/src/colvar/CS2Backbone.cpp
@@ -116,17 +116,17 @@ Additional material and examples can be also found in the tutorial \ref belfast-
 In this first example the chemical shifts are used to calculate a scoring function to be used
 in NMR driven Metadynamics \cite Granata:2013dk :
 
-\verbatim
+\plumedfile
 whole: GROUP ATOMS=2612-2514:-1,961-1:-1,2466-962:-1,2513-2467:-1
 WHOLEMOLECULES ENTITY0=whole
 cs: CS2BACKBONE ATOMS=1-2612 NRES=176 DATA=../data/ TEMPLATE=template.pdb CAMSHIFT NOPBC
 metad: METAD ARG=cs HEIGHT=0.5 SIGMA=0.1 PACE=200 BIASFACTOR=10
 PRINT ARG=cs,metad.bias FILE=COLVAR STRIDE=100
-\endverbatim
+\endplumedfile
 
 In this second example the chemical shifts are used as replica-averaged restrained as in \cite Camilloni:2012je \cite Camilloni:2013hs.
 
-\verbatim
+\plumedfile
 cs: CS2BACKBONE ATOMS=1-174 DATA=data/ NRES=13
 encs: ENSEMBLE ARG=(cs\.hn_.*),(cs\.nh_.*)
 stcs: STATS ARG=encs.* SQDEVSUM PARARG=(cs\.exphn_.*),(cs\.expnh_.*)
@@ -134,7 +134,7 @@ RESTRAINT ARG=stcs.sqdevsum AT=0 KAPPA=0 SLOPE=24
 
 PRINT ARG=(cs\.hn_.*),(cs\.nh_.*) FILE=RESTRAINT STRIDE=100
 
-\endverbatim
+\endplumedfile
 
 (See also \ref WHOLEMOLECULES, \ref STATS, \ref METAD, \ref RESTRAINT and \ref PRINT)
 
diff --git a/src/colvar/Cell.cpp b/src/colvar/Cell.cpp
index 5c59e9f49dbeae2521e1fb62204edd69617dc599..082fa237dbdba39aceb57dd087ef86f2ee71e370 100644
--- a/src/colvar/Cell.cpp
+++ b/src/colvar/Cell.cpp
@@ -36,14 +36,13 @@ Calculate the components of the simulation cell
 
 \par Examples
 The following input tells plumed to print the squared modulo of each of the three lattice vectors
-\verbatim
+\plumedfile
 cell: CELL
 aaa:    COMBINE ARG=cell.ax,cell.ay,cell.az POWERS=2,2,2 PERIODIC=NO
 bbb:    COMBINE ARG=cell.bx,cell.by,cell.bz POWERS=2,2,2 PERIODIC=NO
 ccc:    COMBINE ARG=cell.cx,cell.cy,cell.cz POWERS=2,2,2 PERIODIC=NO
 PRINT ARG=aaa,bbb,ccc
-\endverbatim
-(See also \ref COMBINE and \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/Constant.cpp b/src/colvar/Constant.cpp
index 49423413dfe307cd33c56cc063c8fd4a3f422e5d..247f9ac6d4ab5366692ac0b4eddada801d7c0691 100644
--- a/src/colvar/Constant.cpp
+++ b/src/colvar/Constant.cpp
@@ -44,21 +44,20 @@ The following input instructs plumed to compute the distance
 between atoms 1 and 2. If this distance is between 1.0 and 2.0, it is
 printed. If it is lower than 1.0 (larger than 2.0), 1.0 (2.0) is printed
 
-\verbatim
+\plumedfile
 cn: CONSTANT VALUES=1.0,2.0
 dis: DISTANCE ATOMS=1,2
 sss: SORT ARG=cn.v_0,dis,cn.v_1
 PRINT ARG=sss.2
-\endverbatim
-(See also \ref DISTANCE, \ref SORT, and \ref PRINT).
+\endplumedfile
 
 In case you want to pass a single value you can use VALUE:
-\verbatim
+\plumedfile
 cn: CONSTANT VALUE=1.0
 dis: DISTANCE ATOMS=1
 sss: SORT ARG=cn,dis
 PRINT ARG=sss.1
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/ContactMap.cpp b/src/colvar/ContactMap.cpp
index 9405f172f4fef5820e3c4bfb7389ba3cec06069a..9fe8f82c9289d32026da05cc7103d69da5e56d17 100644
--- a/src/colvar/ContactMap.cpp
+++ b/src/colvar/ContactMap.cpp
@@ -45,17 +45,17 @@ The following example calculates switching functions based on the distances betw
 1 and 2, 3 and 4 and 4 and 5. The values of these three switching functions are then output
 to a file named colvar.
 
-\verbatim
+\plumedfile
 CONTACTMAP ATOMS1=1,2 ATOMS2=3,4 ATOMS3=4,5 ATOMS4=5,6 SWITCH={RATIONAL R_0=1.5} LABEL=f1
 PRINT ARG=f1.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following example calculates the difference of the current contact map with respect
 to a reference provided. In this case REFERENCE is the fraction of contact that is formed
 (i.e. the distance between two atoms transformed with the SWITH), while R_0 is the contact
 distance. WEIGHT gives the relative weight of each contact to the final distance measure.
 
-\verbatim
+\plumedfile
 CONTACTMAP ...
 ATOMS1=1,2 REFERENCE1=0.1 WEIGHT1=0.5
 ATOMS2=3,4 REFERENCE2=0.5 WEIGHT2=1.0
@@ -67,7 +67,7 @@ CMDIST
 ... CONTACTMAP
 
 PRINT ARG=cmap FILE=colvar
-\endverbatim
+\endplumedfile
 
 The next example calculates calculates fraction of native contacts (Q)
 for Trp-cage mini-protein. R_0 is the distance at which the switch function is guaranteed to
@@ -80,7 +80,7 @@ WEIGHT is the 1/(number of contacts) giving equal weight to each contact.
 
 When using native contact Q switch function, please cite \cite best2013
 
-\verbatim
+\plumedfile
 # Full example available in regtest/basic/rt72/
 
 CONTACTMAP ...
@@ -97,9 +97,8 @@ SUM
 ... CONTACTMAP
 
 PRINT ARG=cmap FILE=colvar
-\endverbatim
-
-(See also \ref PRINT and \ref switchingfunction)
+\endplumedfile
+(See also \ref switchingfunction)
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/Coordination.cpp b/src/colvar/Coordination.cpp
index c1aa2da4f28f631d614b31615432f7c6c2d65871..fe0a8b08d332b6d71dc4e294fd3ebe80a9457568 100644
--- a/src/colvar/Coordination.cpp
+++ b/src/colvar/Coordination.cpp
@@ -62,19 +62,19 @@ so that they actually count as "zero".
 \par Examples
 
 The following example instructs plumed to calculate the total coordination number of the atoms in group 1-10 with the atoms in group 20-100.  For atoms 1-10 coordination numbers are calculated that count the number of atoms from the second group that are within 0.3 nm of the central atom.  A neighbour list is used to make this calculation faster, this neighbour list is updated every 100 steps.
-\verbatim
+\plumedfile
 COORDINATION GROUPA=1-10 GROUPB=20-100 R_0=0.3 NLIST NL_CUTOFF=0.5 NL_STRIDE=100
-\endverbatim
+\endplumedfile
 
 The following is a dummy example which should compute the value 0 because the self interaction
 of atom 1 is skipped. Notice that in plumed 2.0 "self interactions" were not skipped, and the
 same calculation should return 1.
-\verbatim
+\plumedfile
 c: COORDINATION GROUPA=1 GROUPB=1 R_0=0.3
 PRINT ARG=c STRIDE=10
-\endverbatim
+\endplumedfile
 
-\verbatim
+\plumedfile
 c1: COORDINATION GROUPA=1-10 GROUPB=1-10 R_0=0.3
 x: COORDINATION GROUPA=1-10 R_0=0.3
 c2: COMBINE ARG=x COEFFICIENTS=2
@@ -82,8 +82,7 @@ c2: COMBINE ARG=x COEFFICIENTS=2
 # since it runs on half of the pairs. Notice that to get the same result you
 # should double it
 PRINT ARG=c1,c2 STRIDE=10
-\endverbatim
-See also \ref PRINT and \ref COMBINE
+\endplumedfile
 
 
 
diff --git a/src/colvar/DHEnergy.cpp b/src/colvar/DHEnergy.cpp
index 23ed21d58dfdd50636fa2f4f1f6d35a461b42b1b..d13aeecdc5ae174ef36b45e7029061d8db15bded 100644
--- a/src/colvar/DHEnergy.cpp
+++ b/src/colvar/DHEnergy.cpp
@@ -56,12 +56,11 @@ Notice that if there are common atoms between GROUPA and GROUPB their interactio
 
 
 \par Examples
-\verbatim
+\plumedfile
 # this is printing the electrostatic interaction between two groups of atoms
 dh: DHENERGY GROUPA=1-10 GROUPB=11-20 EPSILON=80.0 I=0.1 TEMP=300.0
 PRINT ARG=dh
-\endverbatim
-(see also \ref PRINT)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/DRMSD.cpp b/src/colvar/DRMSD.cpp
index 5a3152ad452872f5115e30558eb709817a259e3e..49efb882ad8db179c4fe19e6545fd03dec9300c7 100644
--- a/src/colvar/DRMSD.cpp
+++ b/src/colvar/DRMSD.cpp
@@ -65,15 +65,15 @@ the positions of the atoms in the reference file and their instantaneous
 position. Only pairs of atoms whose distance in the reference structure is within
 0.1 and 0.8 nm are considered.
 
-\verbatim
+\plumedfile
 DRMSD REFERENCE=file.pdb LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8
-\endverbatim
+\endplumedfile
 
 The following tells plumed to calculate a DRMSD value for a pair of molecules.
 
-\verbatim
+\plumedfile
 DRMSD REFERENCE=file.pdb LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8 TYPE=INTER-DRMSD
-\endverbatim
+\endplumedfile
 
 In the input reference file (file.pdb) the atoms in each of the two molecules are separated by a TER
 command as shown below.
diff --git a/src/colvar/Dimer.cpp b/src/colvar/Dimer.cpp
index 827c843654c28946db8a0c5bae9618a821c39ccd..baead008e84b7b611fc71e5aed2f5cfa8f139417 100644
--- a/src/colvar/Dimer.cpp
+++ b/src/colvar/Dimer.cpp
@@ -85,43 +85,43 @@ the temperature of the system.
 
 This line tells Plumed to compute the Dimer interaction energy for every dimer in the system.
 
-\verbatim
+\plumedfile
 dim: DIMER TEMP=300 Q=0.5 ALLATOMS DSIGMA=0.002
-\endverbatim
+\endplumedfile
 
 If the simulation doesn't use virtual sites for the dimers centers of mass,
 Plumed has to know in order to determine correctly the total number of dimers from
 the total number of atoms:
-\verbatim
+\plumedfile
 dim: DIMER TEMP=300 Q=0.5 ALLATOMS DSIGMA=0.002 NOVSITES
-\endverbatim
+\endplumedfile
 
 The NOVSITES flag is not required if one provides the atom serial of each Dimer. This is
 defined as the atom serial of the first bead of the dimer and is thus a number between 1 and N.
 Along with the ATOMS list also the number N of lines describing the first beads has to be given.
 For example, the Dimer interaction energy of dimers 1,5,7 is:
-\verbatim
+\plumedfile
 dim: DIMER TEMP=300 Q=0.5 NATOMS=N ATOMS=1,5,7 DSIGMA=0.002
-\endverbatim
+\endplumedfile
 
 In a Replica Exchange simulation the keyword DSIGMA can be used in two ways:
 if a plumed.n.dat file is provided for each replica, then DSIGMA is passed as a single value,
 like in the previous examples, and each replica will read its own DSIGMA value. If
 a unique plumed.dat is given, DSIGMA has to be a list containing a value for each replica.
 For 4 replicas:
-\verbatim
+\plumedfile
 dim: DIMER TEMP=300 Q=0.5 NATOMS=N ATOMS=1,5,7 DSIGMA=0.002,0.002,0.004,0.01
-\endverbatim
+\endplumedfile
 
 
 \par Usage of the CV
 
 The dimer interaction is not coded in the driver program and has to be inserted
 in the hamiltonian of the system as a linear RESTRAINT (see \ref RESTRAINT):
-\verbatim
+\plumedfile
 dim: DIMER TEMP=300 Q=0.5 ALLATOMS DSIGMA=0.002
 RESTRAINT ARG=dim AT=0 KAPPA=0 SLOPE=1 LABEL=dimforces
-\endverbatim
+\endplumedfile
 
 In a replica exchange, Metadynamics (see \ref METAD) can be used on the Dimer CV to reduce
 the number of replicas. Just keep in mind that METAD SIGMA values should be tuned
diff --git a/src/colvar/Dipole.cpp b/src/colvar/Dipole.cpp
index 09a4a6580032fcc3a05bc04fa41d878e54c5bddd..3cffa3553e0bda7c637c26b7aa7e607918d1a5d1 100644
--- a/src/colvar/Dipole.cpp
+++ b/src/colvar/Dipole.cpp
@@ -37,11 +37,10 @@ Calculate the dipole moment for a group of atoms.
 \par Examples
 The following tells plumed to calculate the dipole of the group of atoms containing
 the atoms from 1-10 and print it every 5 steps
-\verbatim
+\plumedfile
 d: DIPOLE GROUP=1-10
 PRINT FILE=output STRIDE=5 ARG=5
-\endverbatim
-(see also \ref PRINT)
+\endplumedfile
 
 \attention
 If the total charge Q of the group in non zero, then a charge Q/N will be subtracted to every atom,
diff --git a/src/colvar/Distance.cpp b/src/colvar/Distance.cpp
index c3a981771540c28700c91715d06095fc097667f5..279961b7644759ed2012949b7c0b6037e88880fc 100644
--- a/src/colvar/Distance.cpp
+++ b/src/colvar/Distance.cpp
@@ -49,46 +49,44 @@ better to use SCALED_COMPONENTS.
 
 The following input tells plumed to print the distance between atoms 3 and 5,
 the distance between atoms 2 and 4 and the x component of the distance between atoms 2 and 4.
-\verbatim
+\plumedfile
 d1:  DISTANCE ATOMS=3,5
 d2:  DISTANCE ATOMS=2,4
 d2c: DISTANCE ATOMS=2,4 COMPONENTS
 PRINT ARG=d1,d2,d2c.x
-\endverbatim
-(See also \ref PRINT).
+\endplumedfile
 
 The following input computes the end-to-end distance for a polymer
 of 100 atoms and keeps it at a value around 5.
-\verbatim
+\plumedfile
 WHOLEMOLECULES ENTITY0=1-100
 e2e: DISTANCE ATOMS=1,100 NOPBC
 RESTRAINT ARG=e2e KAPPA=1 AT=5
-\endverbatim
-(See also \ref WHOLEMOLECULES and \ref RESTRAINT).
+\endplumedfile
 
 Notice that NOPBC is used
 to be sure that if the end-to-end distance is larger than half the simulation
 box the distance is compute properly. Also notice that, since many MD
 codes break molecules across cell boundary, it might be necessary to
 use the \ref WHOLEMOLECULES keyword (also notice that it should be
-_before_ distance). The list of atoms provided to WHOLEMOLECULES
+_before_ distance). The list of atoms provided to \ref WHOLEMOLECULES
 here contains all the atoms between 1 and 100. Strictly speaking, this
 is not necessary. If you know for sure that atoms with difference in
 the index say equal to 10 are _not_ going to be farther than half cell
 you can e.g. use
-\verbatim
+\plumedfile
 WHOLEMOLECULES ENTITY0=1,10,20,30,40,50,60,70,80,90,100
 e2e: DISTANCE ATOMS=1,100 NOPBC
 RESTRAINT ARG=e2e KAPPA=1 AT=5
-\endverbatim
-Just be sure that the ordered list provide to WHOLEMOLECULES has the following
+\endplumedfile
+Just be sure that the ordered list provide to \ref WHOLEMOLECULES has the following
 properties:
 - Consecutive atoms should be closer than half-cell throughout the entire simulation.
 - Atoms required later for the distance (e.g. 1 and 100) should be included in the list
 
 The following example shows how to take into account periodicity e.g.
 in z-component of a distance
-\verbatim
+\plumedfile
 # this is a center of mass of a large group
 c: COM ATOMS=1-100
 # this is the distance between atom 101 and the group
@@ -99,8 +97,7 @@ d: DISTANCE ATOMS=c,101 COMPONENTS
 dz: COMBINE ARG=d.z PERIODIC=-10,10
 # metadynamics on dd
 METAD ARG=dz SIGMA=0.1 HEIGHT=0.1 PACE=200
-\endverbatim
-(see also \ref COM, \ref COMBINE, and \ref METAD)
+\endplumedfile
 
 Using SCALED_COMPONENTS this problem should not arise because they are always periodic
 with domain (-0.5,+0.5).
diff --git a/src/colvar/ERMSD.cpp b/src/colvar/ERMSD.cpp
index 3602df0868737aa018969fa45be49a32715501b4..fb5ff6cada2e84ba9d984722777ea0800c038417 100644
--- a/src/colvar/ERMSD.cpp
+++ b/src/colvar/ERMSD.cpp
@@ -84,9 +84,9 @@ for a multi molecular system.
 Calculate the eRMSD from reference structure reference.pdb using the default cutoff (2.4). The list of residues involved in the calculation has to be specified. In this example, the eRMSD is calculated
 considering residues 1,2,3,4,5,6.
 
-\verbatim
+\plumedfile
 eRMSD1: ERMSD REFERENCE=reference.pdb ATOMS=@lcs-1,@lcs-2,@lcs-3,@lcs-4,@lcs-5,@lcs-6
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/Energy.cpp b/src/colvar/Energy.cpp
index 6f79eb883176bc996023dac26750b2286904c564..4b3f80736ef2a5d93163943604469821e9c4d9e5 100644
--- a/src/colvar/Energy.cpp
+++ b/src/colvar/Energy.cpp
@@ -48,11 +48,10 @@ using GROMACS with lambda replica exchange of with plumed-hrex branch.
 
 \par Examples
 The following input instructs plumed to print the energy of the system
-\verbatim
+\plumedfile
 ene: ENERGY
 PRINT ARG=ene
-\endverbatim
-(See also \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/Fake.cpp b/src/colvar/Fake.cpp
index 70bd7f278968c9192bba197ca87ee524a5eaecd0..20d70aed9ab5a93e06457ba15c75fc79f15e50ea 100644
--- a/src/colvar/Fake.cpp
+++ b/src/colvar/Fake.cpp
@@ -39,10 +39,9 @@ and just support input and period definition
 
 \par Examples
 
-\verbatim
+\plumedfile
 FAKE ATOMS=1 PERIODIC=-3.14,3.14   LABEL=d2
-\endverbatim
-(See also \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/FretEfficiency.cpp b/src/colvar/FretEfficiency.cpp
index 02600e10574cd08165d4f82ac203dabc910f778c..e6f1494529cf9ccf5df2d718e44c798d180f7e9f 100644
--- a/src/colvar/FretEfficiency.cpp
+++ b/src/colvar/FretEfficiency.cpp
@@ -51,22 +51,20 @@ boundary conditions. This behavior can be changed with the NOPBC flag.
 The following input tells plumed to print the FRET efficiencies
 calculated as a function of the distance between atoms 3 and 5 and
 the distance between atoms 2 and 4.
-\verbatim
+\plumedfile
 fe1:  FRET ATOMS=3,5 R0=5.5
 fe2:  FRET ATOMS=2,4 R0=5.5
 PRINT ARG=fe1,fe2
-\endverbatim
-(See also \ref PRINT).
+\endplumedfile
 
 The following input computes the FRET efficiency calculated on the
 terminal atoms of a polymer
 of 100 atoms and keeps it at a value around 0.5.
-\verbatim
+\plumedfile
 WHOLEMOLECULES ENTITY0=1-100
 fe: FRET ATOMS=1,100 R0=5.5 NOPBC
 RESTRAINT ARG=fe KAPPA=100 AT=0.5
-\endverbatim
-(See also \ref WHOLEMOLECULES and \ref RESTRAINT).
+\endplumedfile
 
 Notice that NOPBC is used
 to be sure that if the distance is larger than half the simulation
diff --git a/src/colvar/Gyration.cpp b/src/colvar/Gyration.cpp
index 791f6a3ed81b0f950c5f430218ee545887703201..3764ad404bef8faf64af02476b36038a9780d475 100644
--- a/src/colvar/Gyration.cpp
+++ b/src/colvar/Gyration.cpp
@@ -75,11 +75,10 @@ periodic image.
 
 The following input tells plumed to print the radius of gyration of the
 chain containing atoms 10 to 20.
-\verbatim
+\plumedfile
 GYRATION TYPE=RADIUS ATOMS=10-20 LABEL=rg
 PRINT ARG=rg STRIDE=1 FILE=colvar
-\endverbatim
-(See also \ref PRINT)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/Implicit.cpp b/src/colvar/Implicit.cpp
index ad53f243cbfd352570d46de5a569ebb354181ac3..b43b3ad604fdaf0f41136d639180fcc0d889cd45 100644
--- a/src/colvar/Implicit.cpp
+++ b/src/colvar/Implicit.cpp
@@ -58,7 +58,7 @@ where \f$\Delta G^\mathrm{free}_i\f$ is the solvation free energy of the isolate
 The output from this collective variable, the free energy of solvation, can be used with the \ref BIASVALUE keyword to provide implicit solvation to a system. All parameters are designed to be used with a modified CHARMM36 force field. It takes only non-hydrogen atoms as input, these can be conveniently specified using the \ref GROUP action with the NDX_GROUP parameter. To speed up the calculation, IMPLICIT internally uses a neighbourlist with a cutoff dependent on the type of atom (maximum of 1.95 nm). This cutoff can be extended further by using the NL_BUFFER keyword.
 
 \par Examples
-\verbatim
+\plumedfile
 MOLINFO MOLTYPE=protein STRUCTURE=peptide.pdb
 WHOLEMOLECULES ENTITY0=1-111
 
@@ -72,8 +72,7 @@ solv: IMPLICIT ATOMS=protein-h NL_STRIDE=10 NL_BUFFER=0.2
 bias: BIASVALUE ARG=solv
 
 PRINT ARG=solv FILE=SOLV
-\endverbatim
-(see also \ref PRINT, \ref GROUP, \ref MOLINFO, \ref WHOLEMOLECULES)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/Jcoupling.cpp b/src/colvar/Jcoupling.cpp
index d192b8bbd6e58c6a039a1796822861b606fba64a..8d8975f91b1e281616f5429546df8c0606a88e9b 100644
--- a/src/colvar/Jcoupling.cpp
+++ b/src/colvar/Jcoupling.cpp
@@ -62,7 +62,7 @@ In the following example we calculate the Ha-N J-coupling from a set of atoms in
 dihedral \f$\psi\f$ angles in the peptide backbone. We also add the experimental datapoints and compute
 the correlation and other measures and finally print the results.
 
-\verbatim
+\plumedfile
 
 MOLINFO MOLTYPE=protein STRUCTURE=peptide.pdb
 WHOLEMOLECULES ENTITY0=1-111
@@ -81,11 +81,7 @@ JCOUPLING ...
 jhanst: STATS ARG=(jhan\.j_.*) PARARG=(jhan\.exp_.*)
 
 PRINT ARG=jhanst.*,jhan.* FILE=COLVAR STRIDE=100
-
-ENDPLUMED
-
-\endverbatim
-(See also \ref PRINT, \ref STATS)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/MultiRMSD.cpp b/src/colvar/MultiRMSD.cpp
index 7656c232be81b77c7a7220a3babe19ed77aa217b..4cf4ff08bfcdf0424ad531bfa4ef7525524f59e4 100644
--- a/src/colvar/MultiRMSD.cpp
+++ b/src/colvar/MultiRMSD.cpp
@@ -99,17 +99,17 @@ The following tells plumed to calculate the RMSD distance between
 the positions of the atoms in the reference file and their instantaneous
 position.  The Kearseley algorithm for each of the domains.
 
-\verbatim
+\plumedfile
 MULTI-RMSD REFERENCE=file.pdb TYPE=MULTI-OPTIMAL
-\endverbatim
+\endplumedfile
 
 The following tells plumed to calculate the RMSD distance btween the positions of
 the atoms in the domains of reference the reference structure and their instantaneous
 positions.  Here distances are calculated using the \ref DRMSD measure.
 
-\verbatim
+\plumedfile
 MULTI-RMSD REFERENCE=file.pdb TYPE=MULTI-DRMSD
-\endverbatim
+\endplumedfile
 
 in this case it is possible to use the following DRMSD options in the pdb file using the REMARK syntax:
 \verbatim
diff --git a/src/colvar/NOE.cpp b/src/colvar/NOE.cpp
index 82d2b4372a6a885f014351df0ba950be6d50d8b3..b46bd65978621d4098392dd7b86db172f71370f8 100644
--- a/src/colvar/NOE.cpp
+++ b/src/colvar/NOE.cpp
@@ -52,7 +52,7 @@ In the following examples three noes are defined, the first is calculated based
 of atom 1-2 and 3-2; the second is defined by the distance 5-7 and the third by the distances
 4-15,4-16,8-15,8-16.
 
-\verbatim
+\plumedfile
 NOE ...
 GROUPA1=1,3 GROUPB1=2,2
 GROUPA2=5 GROUPB2=7
@@ -61,8 +61,7 @@ LABEL=noes
 ... NOE
 
 PRINT ARG=noes.* FILE=colvar
-\endverbatim
-(See also \ref PRINT)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/PCARMSD.cpp b/src/colvar/PCARMSD.cpp
index f5d1cf6f2cd4451d6313d68d545cb774d0dc41ca..02861e10963e49e08558e9ce991e7a34e695a866 100644
--- a/src/colvar/PCARMSD.cpp
+++ b/src/colvar/PCARMSD.cpp
@@ -57,9 +57,9 @@ Note that beta and occupancy values in the pdb are neglected and all the weights
 
 \par Examples
 
-\verbatim
+\plumedfile
 PCARMSD AVERAGE=file.pdb EIGENVECTORS=eigenvectors.pdb
-\endverbatim
+\endplumedfile
 
 The input is taken so to be compatible with the output you get from g_covar utility of gromacs (suitably adapted to have a pdb input format).
 
diff --git a/src/colvar/PRE.cpp b/src/colvar/PRE.cpp
index 906826b1f514fe8a04a073e0c31f998e4af06544..82c94eee23f3a4f10465b27cabc5b537ef209ba8 100644
--- a/src/colvar/PRE.cpp
+++ b/src/colvar/PRE.cpp
@@ -47,7 +47,7 @@ In the following example five PRE intensities are calculated using the distance
 oxigen of the spin label and the backbone hydrogens. Omega is the NMR frequency, RTWO the
 R2 for the hydrogens, INEPT of 8 ms for the experiment and a TAUC of 1.21 ns
 
-\verbatim
+\plumedfile
 PRE ...
 LABEL=HN_pre
 INEPT=8
@@ -63,7 +63,7 @@ GROUPA5=451 RTWO5=0.0086341843
 
 PRINT ARG=HN_pre.* FILE=PRE.dat STRIDE=1
 
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/PathMSD.cpp b/src/colvar/PathMSD.cpp
index cdfd96ee1d2a1f025a451e7e1df1889cff42afb5..17567f2d42b4001705bc34a895eb12de98db6fb1 100644
--- a/src/colvar/PathMSD.cpp
+++ b/src/colvar/PathMSD.cpp
@@ -42,10 +42,10 @@ in input ("sss" component) and the distance from them ("zzz" component).
 Here below is a case where you have defined three frames and you want to
 calculate the progress along the path and the distance from it in p1
 
-\verbatim
+\plumedfile
 p1: PATHMSD REFERENCE=file.pdb  LAMBDA=500.0 NEIGH_STRIDE=4 NEIGH_SIZE=8
 PRINT ARG=p1.sss,p1.zzz STRIDE=1 FILE=colvar FMT=%8.4f
-\endverbatim
+\endplumedfile
 
 note that NEIGH_STRIDE=4 NEIGH_SIZE=8 control the neighborlist parameter (optional but
 recommended for performance) and states that the neighbor list will be calculated every 4
diff --git a/src/colvar/Position.cpp b/src/colvar/Position.cpp
index a381e039e08212c6fac510c7e4b6ec5532cead40..b66ae5e6ca03a37de32414e7d5a442600ef36bb5 100644
--- a/src/colvar/Position.cpp
+++ b/src/colvar/Position.cpp
@@ -54,14 +54,12 @@ This can be done e.g. using \ref FIT_TO_TEMPLATE.
 
 \par Examples
 
-\verbatim
+\plumedfile
 # align to a template
 FIT_TO_TEMPLATE REFERENCE=ref.pdb
 p: POSITION ATOM=3
 PRINT ARG=p.x,p.y,p.z
-\endverbatim
-(see also \ref FIT_TO_TEMPLATE)
-
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/PropertyMap.cpp b/src/colvar/PropertyMap.cpp
index a8a21fef2057803070dd7dc3a46aee30ad62a7ac..eb1c27f0255a2b44d1482133ff2993e57cd3dc99 100644
--- a/src/colvar/PropertyMap.cpp
+++ b/src/colvar/PropertyMap.cpp
@@ -46,10 +46,10 @@ where the parameters \f$X_i\f$  and  \f$Y_i\f$ are provided in the input pdb (al
  \f$D_i(x)\f$  is the MSD after optimal alignment calculated on the pdb frames you input (see Kearsley).
 
 \par Examples
-\verbatim
+\plumedfile
 p3: PROPERTYMAP REFERENCE=../../trajectories/path_msd/allv.pdb PROPERTY=X,Y LAMBDA=69087 NEIGH_SIZE=8 NEIGH_STRIDE=4
 PRINT ARG=p3.X,p3.Y,p3.zzz STRIDE=1 FILE=colvar FMT=%8.4f
-\endverbatim
+\endplumedfile
 
 note that NEIGH_STRIDE=4 NEIGH_SIZE=8 control the neighborlist parameter (optional but
 recommended for performance) and states that the neighbor list will be calculated every 4
diff --git a/src/colvar/Puckering.cpp b/src/colvar/Puckering.cpp
index b90bdfbc4c732f6bb69932e7c4ddc7d68d17a0cf..66edb54c5d0e3fd1fa476cee61d43676782097c7 100644
--- a/src/colvar/Puckering.cpp
+++ b/src/colvar/Puckering.cpp
@@ -60,11 +60,11 @@ namespace colvar {
  \par Examples
 
  This input tells plumed to print the puckering phase angle of the 3rd nucleotide of a RNA molecule on file COLVAR.
- \verbatim
+ \plumedfile
  MOLINFO STRUCTURE=rna.pdb MOLTYPE=rna
  PUCKERING ATOMS=@sugar-3 LABEL=puck
  PRINT ARG=puck.phs FILE=COLVAR
- \endverbatim
+ \endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/RDC.cpp b/src/colvar/RDC.cpp
index 18db380141d62a7dc45c604a662c9451abb5e9bf..43dc4c1ee3f563989c03d597b7e74b007361cc7e 100644
--- a/src/colvar/RDC.cpp
+++ b/src/colvar/RDC.cpp
@@ -94,7 +94,7 @@ respect to a set of experimental data is calculated and restrained. In addition,
 and only for analysis purposes, the same RDCs are calculated using a Single Value
 Decomposition algorithm.
 
-\verbatim
+\plumedfile
 RDC ...
 GYROM=-72.5388
 SCALE=1.0
@@ -124,8 +124,7 @@ LABEL=svd
 
 PRINT ARG=nh.corr,rdce.bias FILE=colvar
 PRINT ARG=svd.* FILE=svd
-\endverbatim
-(See also \ref PRINT, \ref RESTRAINT)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/colvar/RMSD.cpp b/src/colvar/RMSD.cpp
index 8801dfeed200153fcd9368e4e72cbaff987774e5..ffdcefeebd071303979a9ac0cf5fa5e4821b1e8f 100644
--- a/src/colvar/RMSD.cpp
+++ b/src/colvar/RMSD.cpp
@@ -133,9 +133,9 @@ The following tells plumed to calculate the RMSD distance between
 the positions of the atoms in the reference file and their instantaneous
 position.  The Kearseley algorithm is used so this is done optimally.
 
-\verbatim
+\plumedfile
 RMSD REFERENCE=file.pdb TYPE=OPTIMAL
-\endverbatim
+\endplumedfile
 
 ...
 
diff --git a/src/colvar/Template.cpp b/src/colvar/Template.cpp
index 85cb95d27ec8837ed4b37a8cf64ec4508b6c0aa8..6fb3c58451c46dc18ee631315692c5920f2e1df4 100644
--- a/src/colvar/Template.cpp
+++ b/src/colvar/Template.cpp
@@ -40,11 +40,11 @@ This file provides a template for if you want to introduce a new CV.
 
 <!---You should put an example of how to use your CV here--->
 
-\verbatim
+\plumedfile
 # This should be a sample input.
 t: TEMPLATE ATOMS=1,2
 PRINT ARG=t STRIDE=100 FILE=COLVAR
-\endverbatim
+\endplumedfile
 <!---You should reference here the other actions used in this example--->
 (see also \ref PRINT)
 
diff --git a/src/colvar/Torsion.cpp b/src/colvar/Torsion.cpp
index 274d0565a424354ef9a8d88ec84edd6037e8859a..4b697b06976c3bdc7508a4961a2c6f4d58b8dc00 100644
--- a/src/colvar/Torsion.cpp
+++ b/src/colvar/Torsion.cpp
@@ -43,23 +43,23 @@ orthogonal to an axis.
 
 This input tells plumed to print the torsional angle between atoms 1, 2, 3 and 4
 on file COLVAR.
-\verbatim
+\plumedfile
 t: TORSION ATOMS=1,2,3,4
 # this is an alternative, equivalent, definition:
 # t: TORSION VECTOR1=2,1 AXIS=2,3 VECTOR2=3,4
 PRINT ARG=t FILE=COLVAR
-\endverbatim
+\endplumedfile
 
 If you are working with a protein you can specify the special named torsion angles \f$\phi\f$, \f$\psi\f$, \f$\omega\f$ and \f$\chi_1\f$
 by using TORSION in combination with the \ref MOLINFO command.  This can be done by using the following
 syntax.
 
-\verbatim
+\plumedfile
 MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
 t1: TORSION ATOMS=@phi-3
 t2: TORSION ATOMS=@psi-4
 PRINT ARG=t1,t2 FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 
 Here, \@phi-3 tells plumed that you would like to calculate the \f$\phi\f$ angle in the third residue of the protein.
 Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle of the 4th residue of the protein.
diff --git a/src/colvar/Volume.cpp b/src/colvar/Volume.cpp
index 993dc802af7e6da736f327a09a12a02774ed53cd..837f6ce75871ad84840e8ea1d5be7a20dee2bc77 100644
--- a/src/colvar/Volume.cpp
+++ b/src/colvar/Volume.cpp
@@ -36,11 +36,10 @@ Calculate the volume of the simulation box.
 
 \par Examples
 The following input tells plumed to print the volume of the system
-\verbatim
+\plumedfile
 vol: VOLUME
 PRINT ARG=vol
-\endverbatim
-(See also \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/crystallization/Fccubic.cpp b/src/crystallization/Fccubic.cpp
index fc1793ab39cea57d52603652660a35153aa5761e..47154f4c6f7d393c605f85a0a656b1da4717fe0c 100644
--- a/src/crystallization/Fccubic.cpp
+++ b/src/crystallization/Fccubic.cpp
@@ -63,10 +63,10 @@ so on.  Notice also that you can rotate the reference frame if you are using a n
 The following input calculates the FCCUBIC parameter for the 64 atoms in the system
 and then calculates and prints the average value for this quantity.
 
-\verbatim
+\plumedfile
 FCCUBIC SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5} MEAN LABEL=d
 PRINT ARG=d.* FILE=colv
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/crystallization/Gradient.cpp b/src/crystallization/Gradient.cpp
index 173d9e5adb0afbfe84d18ee5c22e13f4a089db57..e3af25856de604e00803a5a136a57895590263c7 100644
--- a/src/crystallization/Gradient.cpp
+++ b/src/crystallization/Gradient.cpp
@@ -39,22 +39,22 @@ The input below calculates the gradient of the density of atoms in the manner
 described in \cite fede-grad in order to detect whether or not atoms are distributed
 uniformly along the x-axis of the simulation cell.
 
-\verbatim
+\plumedfile
 d1: DENSITY SPECIES=1-50
 s1: GRADIENT ORIGIN=1 DATA=d1 DIR=x NBINS=4 SIGMA=1.0
 PRINT ARG=s1 FILE=colvar
-\endverbatim
+\endplumedfile
 
 The input below calculates the coordination numbers of the 50 atoms in the simulation cell.
 The gradient of this quantity is then evaluated in the manner described using the equation above
 to detect whether the average values of the coordination number are uniformly distributed along the
 x-axis of the simulation cell.
 
-\verbatim
+\plumedfile
 d2: COORDINATIONNUMBER SPECIES=1-50 SWITCH={RATIONAL R_0=2.0} MORE_THAN={EXP R_0=4.0}
 s2: GRADIENT ORIGIN=1 DATA=d2 DIR=x NBINS=4 SIGMA=1.0
 PRINT ARG=s2 FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/crystallization/InterMolecularTorsions.cpp b/src/crystallization/InterMolecularTorsions.cpp
index 3839e094d6793285e2f3909fe557b0d656562fa7..58294492f6fa982d8463e5cc476f5dde4a74080c 100644
--- a/src/crystallization/InterMolecularTorsions.cpp
+++ b/src/crystallization/InterMolecularTorsions.cpp
@@ -64,12 +64,12 @@ between the two molecules.  As such the torsional angles between molecules that
 histogram while the torsional angles between molecules that are far apart does not contribute to the histogram.  The histogram is
 averaged over the whole trajectory and output once all the trajectory frames have been read.
 
-\verbatim
+\plumedfile
 m1: MOLECULES MOL1=1,2 MOL2=3,4 MOL3=5,6 MOL4=7,8
 tt_p: INTERMOLECULARTORSIONS MOLS=m1 SWITCH={RATIONAL R_0=0.25 D_0=2.0 D_MAX=3.0}
 htt_p: HISTOGRAM DATA=tt_p GRID_MIN=-pi GRID_MAX=pi BANDWIDTH=0.1 GRID_BIN=200 STRIDE=1
 DUMPGRID GRID=htt_p FILE=myhist.out
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/crystallization/MoleculeOrientation.cpp b/src/crystallization/MoleculeOrientation.cpp
index 5bade7022c81bbdc6b6a391f4723cf704daee5e5..6001f6f095c25bd9d16b7acc0ab9f23c14252df1 100644
--- a/src/crystallization/MoleculeOrientation.cpp
+++ b/src/crystallization/MoleculeOrientation.cpp
@@ -39,10 +39,10 @@ The following input tells plumed to calculate the distances between two of the a
 This is done for the same set of atoms four different molecules and the average separation is then
 calculated.
 
-\verbatim
+\plumedfile
 MOLECULES MOL1=1,2 MOL2=3,4 MOL3=5,6 MOL4=7,8 MEAN LABEL=mm
 PRINT ARG=mm.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 
 */
diff --git a/src/crystallization/Q3.cpp b/src/crystallization/Q3.cpp
index 6c3aeedc2f6f8d2fbfe9606e4d37f50e144535af..038115a92158389807b110c24ea959e6785c31c7 100644
--- a/src/crystallization/Q3.cpp
+++ b/src/crystallization/Q3.cpp
@@ -60,28 +60,28 @@ the \f$q_{3}\f$ vectors on adjacent atoms.  More information on these variables
 The following command calculates the average Q3 parameter for the 64 atoms in a box of Lennard Jones and prints this
 quantity to a file called colvar:
 
-\verbatim
+\plumedfile
 Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN LABEL=q3
 PRINT ARG=q3.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following command calculates the histogram of Q3 parameters for the 64 atoms in a box of Lennard Jones and prints these
 quantities to a file called colvar:
 
-\verbatim
+\plumedfile
 Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=q3
 PRINT ARG=q3.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following command could be used to measure the Q3 paramters that describe the arrangement of chlorine ions around the
 sodium atoms in NaCl.  The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
 with the 64 Na\f$^+\f$ ions followed by the 64 Cl\f$-\f$ ions.  Once again the average Q3 paramter is calculated and output to a
 file called colvar
 
-\verbatim
+\plumedfile
 Q3 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN LABEL=q3
 PRINT ARG=q3.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
@@ -136,30 +136,30 @@ adjacent atoms is correlated.
 The following command calculates the average value of the LOCAL_Q3 parameter for the 64 Lennard Jones atoms in the system under study and prints this
 quantity to a file called colvar.
 
-\verbatim
+\plumedfile
 Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q3
 LOCAL_Q3 ARG=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq3
 PRINT ARG=lq3.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following input calculates the distribution of LOCAL_Q3 parameters at any given time and outputs this information to a file.
 
-\verbatim
+\plumedfile
 Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q3
 LOCAL_Q3 ARG=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq3
 PRINT ARG=lq3.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following calculates the LOCAL_Q3 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
 are done with those of all the other atoms in the system.  The final quantity is the average and is outputted to a file
 
-\verbatim
+\plumedfile
 Q3 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q3a
 Q3 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q3b
 
 LOCAL_Q3 ARG=q3a,q3b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w3
 PRINT ARG=w3.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/crystallization/Q4.cpp b/src/crystallization/Q4.cpp
index fa76b0d10e762b3b2efd30382c0ca09cfdc35670..a4f7d327b2545848609996a8d2fab50965a9616e 100644
--- a/src/crystallization/Q4.cpp
+++ b/src/crystallization/Q4.cpp
@@ -60,28 +60,28 @@ the \f$q_{4}\f$ vectors on adjacent atoms.  More information on these variables
 The following command calculates the average Q4 parameter for the 64 atoms in a box of Lennard Jones and prints this
 quantity to a file called colvar:
 
-\verbatim
+\plumedfile
 Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN LABEL=q4
 PRINT ARG=q4.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following command calculates the histogram of Q4 parameters for the 64 atoms in a box of Lennard Jones and prints these
 quantities to a file called colvar:
 
-\verbatim
+\plumedfile
 Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=q4
 PRINT ARG=q4.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following command could be used to measure the Q4 paramters that describe the arrangement of chlorine ions around the
 sodium atoms in NaCl.  The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
 with the 64 Na\f$^+\f$ ions followed by the 64 Cl\f$-\f$ ions.  Once again the average Q4 paramter is calculated and output to a
 file called colvar
 
-\verbatim
+\plumedfile
 Q4 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN LABEL=q4
 PRINT ARG=q4.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
@@ -136,30 +136,30 @@ adjacent atoms is correlated.
 The following command calculates the average value of the LOCAL_Q4 parameter for the 64 Lennard Jones atoms in the system under study and prints this
 quantity to a file called colvar.
 
-\verbatim
+\plumedfile
 Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q4
 LOCAL_Q4 ARG=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq4
 PRINT ARG=lq4.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following input calculates the distribution of LOCAL_Q4 parameters at any given time and outputs this information to a file.
 
-\verbatim
+\plumedfile
 Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q4
 LOCAL_Q4 ARG=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq4
 PRINT ARG=lq4.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following calculates the LOCAL_Q4 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
 are done with those of all the other atoms in the system.  The final quantity is the average and is outputted to a file
 
-\verbatim
+\plumedfile
 Q4 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q4a
 Q4 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q4b
 
 LOCAL_Q4 ARG=q4a,q4b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w4
 PRINT ARG=w4.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/crystallization/Q6.cpp b/src/crystallization/Q6.cpp
index 447b0e1cb42d6ec71c558e22be071f1773586d93..a848174ab88fbee83044c3488cbec611d3e8f190 100644
--- a/src/crystallization/Q6.cpp
+++ b/src/crystallization/Q6.cpp
@@ -60,28 +60,28 @@ the \f$q_{6}\f$ vectors on adjacent atoms.  More information on these variables
 The following command calculates the average Q6 parameter for the 64 atoms in a box of Lennard Jones and prints this
 quantity to a file called colvar:
 
-\verbatim
+\plumedfile
 Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN LABEL=q6
 PRINT ARG=q6.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following command calculates the histogram of Q6 parameters for the 64 atoms in a box of Lennard Jones and prints these
 quantities to a file called colvar:
 
-\verbatim
+\plumedfile
 Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=q6
 PRINT ARG=q6.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following command could be used to measure the Q6 paramters that describe the arrangement of chlorine ions around the
 sodium atoms in NaCl.  The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
 with the 64 Na\f$^+\f$ ions followed by the 64 Cl\f$-\f$ ions.  Once again the average Q6 paramter is calculated and output to a
 file called colvar
 
-\verbatim
+\plumedfile
 Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN LABEL=q6
 PRINT ARG=q6.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
@@ -136,30 +136,30 @@ adjacent atoms is correlated.
 The following command calculates the average value of the LOCAL_Q6 parameter for the 64 Lennard Jones atoms in the system under study and prints this
 quantity to a file called colvar.
 
-\verbatim
+\plumedfile
 Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q6
 LOCAL_Q6 ARG=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq6
 PRINT ARG=lq6.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following input calculates the distribution of LOCAL_Q6 parameters at any given time and outputs this information to a file.
 
-\verbatim
+\plumedfile
 Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q6
 LOCAL_Q6 ARG=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq6
 PRINT ARG=lq6.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following calculates the LOCAL_Q6 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
 are done with those of all the other atoms in the system.  The final quantity is the average and is outputted to a file
 
-\verbatim
+\plumedfile
 Q6 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q6a
 Q6 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q6b
 
 LOCAL_Q6 ARG=q4a,q4b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w4
 PRINT ARG=w6.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/crystallization/SMAC.cpp b/src/crystallization/SMAC.cpp
index ce83df73808be7d532e308f945b942590e0a5b68..3d7085a979e6979fe76d70b420e2f8fa8aad17cf 100644
--- a/src/crystallization/SMAC.cpp
+++ b/src/crystallization/SMAC.cpp
@@ -65,7 +65,7 @@ the indices of three atoms for each of the MOL keywords below we are telling PLU
 numbers to determine the orientation of the molecule that will ultimately be used when calculating the \f$\theta_{ij}\f$
 terms in the formula above.  The atom with the third index meanwhile is used when we calculate \f$r_{ij}\f$.
 
-\verbatim
+\plumedfile
 MOLECULES ...
 MOL1=9,10,9
 MOL2=89,90,89
@@ -86,13 +86,13 @@ SMAC ...
 ... SMAC
 
 PRINT ARG=s2.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 This second example works in a way that is very similar to the previous command.  Now, however,
 the orientation of the molecules is determined by finding the plane that contains the positions
 of three atoms.
 
-\verbatim
+\plumedfile
 PLANES ...
 MOL1=9,10,11
 MOL2=89,90,91
@@ -114,7 +114,7 @@ SMAC ...
 ... SMAC
 
 PRINT ARG=s2.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/crystallization/SimpleCubic.cpp b/src/crystallization/SimpleCubic.cpp
index 4d50cc7d59a408b9f138385d736d1d430b941512..2c15767c272aad5c6a60f595cccd39a86d1c60e8 100644
--- a/src/crystallization/SimpleCubic.cpp
+++ b/src/crystallization/SimpleCubic.cpp
@@ -55,15 +55,15 @@ so on.  Notice also that you can rotate the reference frame if you are using a n
 
 The following input tells plumed to calculate the simple cubic parameter for the atoms 1-100 with themselves.
 The mean value is then calculated.
-\verbatim
+\plumedfile
 SIMPLECUBIC SPECIES=1-100 R_0=1.0 MEAN
-\endverbatim
+\endplumedfile
 
 The following input tells plumed to look at the ways atoms 1-100 are within 3.0 are arranged about atoms
 from 101-110.  The number of simple cubic parameters that are greater than 0.8 is then output
-\verbatim
+\plumedfile
 SIMPLECUBIC SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=0.8 NN=6 MM=12 D_0=0}
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/crystallization/Tetrahedral.cpp b/src/crystallization/Tetrahedral.cpp
index 616cf766bf1bc8ebdbcbede38aa94a8d0d4cb4b0..87c6e33ae87d942c012716bdcfc272538f20dd58 100644
--- a/src/crystallization/Tetrahedral.cpp
+++ b/src/crystallization/Tetrahedral.cpp
@@ -54,20 +54,20 @@ when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zer
 The following command calculates the average value of the tetrahedrality parameter for a set of 64 atoms all of the same type
 and outputs this quantity to a file called colvar.
 
-\verbatim
+\plumedfile
 tt: TETRAHEDRAL SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
 PRINT ARG=tt.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following command calculates the number of tetrahedrality parameters that are greater than 0.8 in a set of 10 atoms.
 In this calculation it is assumed that there are two atom types A and B and that the first coordination sphere of the
 10 atoms of type A contains atoms of type B.  The formula above is thus calculated for ten different A atoms and within
 it the sum over \f$j\f$ runs over 40 atoms of type B that could be in the first coordination sphere.
 
-\verbatim
+\plumedfile
 tt: TETRAHEDRAL SPECIESA=1-10 SPECIESB=11-40 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=0.8}
 PRINT ARG=tt.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/eds/EDS.cpp b/src/eds/EDS.cpp
index e5e575abe167259d29c679cefae64ac164f8c5f7..e7e0e4331d1b9c30f2938ae1e1f73cd317f1233e 100644
--- a/src/eds/EDS.cpp
+++ b/src/eds/EDS.cpp
@@ -63,7 +63,7 @@ It is not possible to set the target value of the observable to zero with the de
 
 The following input for a harmonic oscillator of two beads will adaptively find a linear bias to change the mean and variance to the target values. The PRINT line shows how to access the value of the coupling constants.
 
-\verbatim
+\plumedfile
 dist: DISTANCE ATOMS=1,2
 # this is the squared of the distance
 dist2: COMBINE ARG=dist POWERS=2 PERIODIC=NO
@@ -71,19 +71,19 @@ dist2: COMBINE ARG=dist POWERS=2 PERIODIC=NO
 #bias mean and variance
 eds: EDS ARG=dist,dist2 CENTER=2.0,1.0 PERIOD=50000 TEMP=1.0
 PRINT ARG=dist,dist2,eds.dist_coupling,eds.dist2_coupling,eds.bias,eds.force2 FILE=colvars.dat STRIDE=100
-\endverbatim
+\endplumedfile
 
 Rather than trying to find the coupling constants adaptively, can ramp up to a constant value.
-\verbatim
+\plumedfile
 #ramp couplings from 0,0 to -1,1 over 50000 steps
 eds: EDS ARG=dist,dist2 CENTER=2.0,1.0 FIXED=-1,1 RAMP PERIOD=50000 TEMP=1.0
 
 #same as above, except starting at -0.5,0.5 rather than default of 0,0
 eds: EDS ARG=dist,dist2 CENTER=2.0,1.0 FIXED=-1,1 INIT=-0.5,0.5 RAMP PERIOD=50000 TEMP=1.0
-\endverbatim
+\endplumedfile
 
 A restart file can be added to dump information needed to restart/continue simulation using these parameters every STRIDE.
-\verbatim
+\plumedfile
 #add the option to write to a restart file
 eds: EDS ARG=dist,dist2 CENTER=2.0,1.0 PERIOD=50000 TEMP=1.0 OUT_RESTART=restart.dat
 
@@ -98,7 +98,7 @@ eds: EDS ARG=dist,dist2 CENTER=2.0,1.0 PERIOD=50000 TEMP=1.0 IN_RESTART=restart.
 
 #add the option to read in a previous restart file and continue the bias, but use the mean from the previous run as the starting point
 eds: EDS ARG=dist,dist2 CENTER=2.0,1.0 PERIOD=50000 TEMP=1.0 IN_RESTART=restart.dat MEAN
-\endverbatim
+\endplumedfile
 
 
 */
diff --git a/src/function/Combine.cpp b/src/function/Combine.cpp
index e76ed7b8bb2339ecdfb18fa5603dcd2eb1ec63e3..62f66e35d81b082325724ee919ba968e72e39d6c 100644
--- a/src/function/Combine.cpp
+++ b/src/function/Combine.cpp
@@ -52,23 +52,23 @@ is periodic.
 The following input tells plumed to print the distance between atoms 3 and 5
 its square (as computed from the x,y,z components) and the distance
 again as computed from the square root of the square.
-\verbatim
+\plumedfile
 DISTANCE LABEL=dist      ATOMS=3,5 COMPONENTS
 COMBINE  LABEL=distance2 ARG=dist.x,dist.y,dist.z POWERS=2,2,2 PERIODIC=NO
 COMBINE  LABEL=distance  ARG=distance2 POWERS=0.5 PERIODIC=NO
 PRINT ARG=distance,distance2
-\endverbatim
+\endplumedfile
 (See also \ref PRINT and \ref DISTANCE).
 
 The following input tells plumed to add a restraint on the
 cube of a dihedral angle. Notice that since the angle has a
 periodic domain
 -pi,pi its cube has a domain -pi**3,pi**3.
-\verbatim
+\plumedfile
 t: TORSION ATOMS=1,3,5,7
 c: COMBINE ARG=t POWERS=3 PERIODIC=-31.0062766802998,31.0062766802998
 RESTRAINT ARG=c KAPPA=10 AT=0
-\endverbatim
+\endplumedfile
 
 
 
diff --git a/src/function/Ensemble.cpp b/src/function/Ensemble.cpp
index 9976c25015d86e42ee3d5ea911d1c2d2c47029db..bc1e05fb28569a74cea4b7b471bae55b784fb950 100644
--- a/src/function/Ensemble.cpp
+++ b/src/function/Ensemble.cpp
@@ -38,12 +38,11 @@ Each collective variable is averaged separately and stored in a component labell
 \par Examples
 The following input tells plumed to calculate the distance between atoms 3 and 5
 and the average it over the available replicas.
-\verbatim
+\plumedfile
 dist: DISTANCE ATOMS=3,5
 ens: ENSEMBLE ARG=dist
 PRINT ARG=dist,ens.dist
-\endverbatim
-(See also \ref PRINT and \ref DISTANCE).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/function/FuncPathMSD.cpp b/src/function/FuncPathMSD.cpp
index f4d481d733f1ad840057c769254a3f098642cf85..25a87e7182ba39dcfe26eecf6808a97ad82e096d 100644
--- a/src/function/FuncPathMSD.cpp
+++ b/src/function/FuncPathMSD.cpp
@@ -49,17 +49,17 @@ It is a function of MSD that are obtained by the joint use of MSD variable and S
 Here below is a case where you have defined three frames and you want to
 calculate the progress alng the path and the distance from it in p1
 
-\verbatim
+\plumedfile
 t1: RMSD REFERENCE=frame_1.dat TYPE=OPTIMAL SQUARED
 t2: RMSD REFERENCE=frame_21.dat TYPE=OPTIMAL SQUARED
 t3: RMSD REFERENCE=frame_42.dat TYPE=OPTIMAL SQUARED
 p1: FUNCPATHMSD ARG=t1,t2,t3 LAMBDA=500.0
 PRINT ARG=t1,t2,t3,p1.s,p1.z STRIDE=1 FILE=colvar FMT=%8.4f
-\endverbatim
+\endplumedfile
 
 In this second example is shown how to define a PATH in the \ref CONTACTMAP space:
 
-\verbatim
+\plumedfile
 CONTACTMAP ...
 ATOMS1=1,2 REFERENCE1=0.1
 ATOMS2=3,4 REFERENCE2=0.5
@@ -92,7 +92,7 @@ CMDIST
 
 p1: FUNCPATHMSD ARG=c1,c2,c3 LAMBDA=500.0
 PRINT ARG=c1,c2,c3,p1.s,p1.z STRIDE=1 FILE=colvar FMT=%8.4f
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/function/LocalEnsemble.cpp b/src/function/LocalEnsemble.cpp
index 6e3bd96f689e6001e8e7e4b55bbda0c3a415c2f8..429a5a1436913f1a308dd700a67d8bb5903d742d 100644
--- a/src/function/LocalEnsemble.cpp
+++ b/src/function/LocalEnsemble.cpp
@@ -40,7 +40,7 @@ The following input tells plumed to calculate the chemical shifts for four
 different proteins in the same simulation box then average them, calcualated
 the sum of the squared deviation with respect to the experiemntal values and
 applies a linear restraint.
-\verbatim
+\plumedfile
 MOLINFO STRUCTURE=data/template.pdb
 
 chaina: GROUP ATOMS=1-1640
@@ -68,7 +68,7 @@ sthn: STATS ARG=(enshn\.csa\.hn_.*) PARARG=(csa\.exphn_.*) SQDEVSUM
 stnh: STATS ARG=(ensnh\.csa\.nh_.*) PARARG=(csa\.expnh_.*) SQDEVSUM
 
 res: RESTRAINT ARG=stca.*,stcb.*,stco.*,sthn.*,stnh.* AT=0.,0.,0.,0.,0. KAPPA=0.,0.,0.,0.,0 SLOPE=16.,16.,12.,24.,0.5
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/function/Matheval.cpp b/src/function/Matheval.cpp
index 59e6e9bf05007580e9ece7b1506a4c55103d3869..d21907e68eb6e92c6533d4755439390a8ee6af41 100644
--- a/src/function/Matheval.cpp
+++ b/src/function/Matheval.cpp
@@ -54,14 +54,14 @@ PLUMED has been linked to it
 
 The following input tells plumed to perform a metadynamics
 using as a CV the difference between two distances.
-\verbatim
+\plumedfile
 dAB: DISTANCE ARG=10,12
 dAC: DISTANCE ARG=10,15
 diff: MATHEVAL ARG=dAB,dAC FUNC=y-x PERIODIC=NO
 # notice: the previous line could be replaced with the following
 # diff: COMBINE ARG=dAB,dAC COEFFICIENTS=-1,1
 METAD ARG=diff WIDTH=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
-\endverbatim
+\endplumedfile
 (see also \ref DISTANCE, \ref COMBINE, and \ref METAD).
 Notice that forces applied to diff will be correctly propagated
 to atoms 10, 12, and 15.
@@ -75,7 +75,7 @@ The following input tells plumed to print the angle between vectors
 identified by atoms 1,2 and atoms 2,3
 its square (as computed from the x,y,z components) and the distance
 again as computed from the square root of the square.
-\verbatim
+\plumedfile
 DISTANCE LABEL=d1 ATOMS=1,2 COMPONENTS
 DISTANCE LABEL=d2 ATOMS=2,3 COMPONENTS
 MATHEVAL ...
@@ -86,7 +86,7 @@ MATHEVAL ...
   PERIODIC=NO
 ... MATHEVAL
 PRINT ARG=theta
-\endverbatim
+\endplumedfile
 (See also \ref PRINT and \ref DISTANCE).
 
 Notice that the matheval library implements a large number of functions (trigonometric, exp, log, etc).
@@ -96,13 +96,13 @@ a straightforward implementation of if clauses.
 
 For example, imagine that you want to implement a restraint that only acts when a
 distance is larger than 0.5. You can do it with
-\verbatim
+\plumedfile
 d: DISTANCE ATOMS=10,15
 m: MATHEVAL ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
 # check the function you are applying:
 PRINT ARG=d,n FILE=checkme
 RESTRAINT ARG=d AT=0.5 KAPPA=10.0
-\endverbatim
+\endplumedfile
 (see also \ref DISTANCE, \ref PRINT, and \ref RESTRAINT)
 
 The meaning of the function `0.5*step(0.5-x)+x*step(x-0.5)` is:
@@ -133,7 +133,7 @@ MATHEVAL can be used in combination with \ref DISTANCE to implement variants of
 DISTANCE keyword that were present in PLUMED 1.3 and that allowed to compute
 the distance of a point from a line defined by two other points, or the progression
 along that line.
-\verbatim
+\plumedfile
 # take center of atoms 1 to 10 as reference point 1
 p1: CENTER ATOMS=1-10
 # take center of atoms 11 to 20 as reference point 2
@@ -156,7 +156,7 @@ fromaxis: MATHEVAL ARG=d13,d23,d12,onaxis VAR=x,y,z,o FUNC=(0.5*(y^2+x^2)-o^2-0.
 
 PRINT ARG=onaxis,fromaxis
 
-\endverbatim
+\endplumedfile
 
 Notice that these equations have been used to combine \ref RMSD
 from different snapshots of a protein so as to define
diff --git a/src/function/Piecewise.cpp b/src/function/Piecewise.cpp
index e0f4bbd52ac4ec5b4eb51c3ab9de53f119a1b9cb..ebd886de09e4c77fb8f2685b8d5dc5d1d50d53b2 100644
--- a/src/function/Piecewise.cpp
+++ b/src/function/Piecewise.cpp
@@ -55,15 +55,14 @@ in a vector of values. Each value will be named as the name of the original
 argument with suffix _pfunc.
 
 \par Examples
-\verbatim
+\plumedfile
 dist1: DISTANCE ATOMS=1,10
 dist2: DISTANCE ATOMS=2,11
 
 pw: PIECEWISE POINT0=1,10 POINT1=1,PI POINT2=3,10 ARG=dist1
 ppww: PIECEWISE POINT0=1,10 POINT1=1,PI POINT2=3,10 ARG=dist1,dist2
 PRINT ARG=pw,ppww.dist1_pfunc,ppww.dist2_pfunc
-\endverbatim
-(See also \ref PRINT and \ref DISTANCE).
+\endplumedfile
 
 
 */
diff --git a/src/function/Sort.cpp b/src/function/Sort.cpp
index f1e9b2a71773a1ab8dc98195e654fc7d29dc9d2a..fc8e4ca233822da29806cd3b87e5809526235367 100644
--- a/src/function/Sort.cpp
+++ b/src/function/Sort.cpp
@@ -43,15 +43,14 @@ labelled <em>label</em>.1, the second lowest will be labelled <em>label</em>.2 a
 \par Examples
 The following input tells plumed to print the distance of the closest and of
 the farthest atoms to atom 1, chosen among atoms from 2 to 5
-\verbatim
+\plumedfile
 d12:  DISTANCE ATOMS=1,2
 d13:  DISTANCE ATOMS=1,3
 d14:  DISTANCE ATOMS=1,4
 d15:  DISTANCE ATOMS=1,5
 sort: SORT ARG=d12,d13,d14,d15
 PRINT ARG=sort.1,sort.4
-\endverbatim
-(See also \ref PRINT and \ref DISTANCE).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/function/Stats.cpp b/src/function/Stats.cpp
index 429442e6d833bc41ea63beb9283d3dd096be5b69..438d8f98b86cede7896a77e57e794ee926ad47a5 100644
--- a/src/function/Stats.cpp
+++ b/src/function/Stats.cpp
@@ -41,13 +41,13 @@ from other actions using PARARG (for example using experimental values from coll
 The following input tells plumed to print the distance between three couple of atoms
 and compare them with three reference distances.
 
-\verbatim
+\plumedfile
 d1: DISTANCE ATOMS=10,50
 d2: DISTANCE ATOMS=1,100
 d3: DISTANCE ATOMS=45,75
 st: STATS ARG=d1,d2,d3 PARAMETERS=1.5,4.0,2.0
 PRINT ARG=d1,d2,d3,st.*
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/function/Target.cpp b/src/function/Target.cpp
index 5ce91249563a68b9599abf0806d02b3aa7c95f65..1724cf3ff8263b4c40dfc68e4bd9e35ff7a4c3c0 100644
--- a/src/function/Target.cpp
+++ b/src/function/Target.cpp
@@ -75,12 +75,12 @@ specified in the input.
 The following input calculates the distance between a reference configuration and the instaneous position of the system in the trajectory.
 The position of the reference configuration is specified by providing the values of the distance between atoms 1 and 2 and atoms 3 and 4.
 
-\verbatim
+\plumedfile
 d1: DISTANCE ATOMS=1,2
 d2: DISTANCE ATOMS=3,4
 t1: TARGET REFERENCE=myref.pdb TYPE=EUCLIDEAN
 PRINT ARG=t1 FILE=colvar
-\endverbatim
+\endplumedfile
 
 The contents of the file containing the reference structure (myref.pdb) is shown below.  As you can see you must provide information on the
 labels of the CVs that are being used to define the position of the reference configuration in this file together with the values that these
diff --git a/src/generic/Debug.cpp b/src/generic/Debug.cpp
index 742b3b85420f4d77b99907274a3958945d3b1034..36654b0d90fdb901fcdfc15fa027b3f71112c05d 100644
--- a/src/generic/Debug.cpp
+++ b/src/generic/Debug.cpp
@@ -37,12 +37,12 @@ Can be used while debugging or optimizing plumed.
 
 \par Examples
 
-\verbatim
+\plumedfile
 # print detailed (action-by-action) timers at the end of simulation
 DEBUG DETAILED_TIMERS
 # dump every two steps which are the atoms required from the MD code
 DEBUG logRequestedAtoms STRIDE=2
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/generic/DumpAtoms.cpp b/src/generic/DumpAtoms.cpp
index d6ee377488219851af8eb407533666d13690548b..c144b5f0780082713a8a5946806822a0a4a2f06e 100644
--- a/src/generic/DumpAtoms.cpp
+++ b/src/generic/DumpAtoms.cpp
@@ -70,15 +70,14 @@ Notice that gro/xtc/trr files can only contain coordinates in nm.
 The following input instructs plumed to print out the positions of atoms
 1-10 together with the position of the center of mass of atoms 11-20 every
 10 steps to a file called file.xyz.
-\verbatim
+\plumedfile
 COM ATOMS=11-20 LABEL=c1
 DUMPATOMS STRIDE=10 FILE=file.xyz ATOMS=1-10,c1
-\endverbatim
-(see also \ref COM)
+\endplumedfile
 
 The following input is very similar but dumps a .gro (gromacs) file,
 which also contains atom and residue names.
-\verbatim
+\plumedfile
 # this is required to have proper atom names:
 MOLINFO STRUCTURE=reference.pdb
 # if omitted, atoms will have "X" name...
@@ -87,8 +86,7 @@ COM ATOMS=11-20 LABEL=c1
 DUMPATOMS STRIDE=10 FILE=file.gro ATOMS=1-10,c1
 # notice that last atom is a virtual one and will not have
 # a correct name in the resulting gro file
-\endverbatim
-(see also \ref COM and \ref MOLINFO)
+\endplumedfile
 
 
 */
diff --git a/src/generic/DumpDerivatives.cpp b/src/generic/DumpDerivatives.cpp
index c74343aca0983c529c9521c9d137b271c0c785ad..38c686d7b4aec532684f2ce6f0dcd40f4319520f 100644
--- a/src/generic/DumpDerivatives.cpp
+++ b/src/generic/DumpDerivatives.cpp
@@ -42,11 +42,11 @@ can be done by outputting the derivatives calculated analytically and numericall
 \par Examples
 The following input instructs plumed to write a file called deriv that contains both the
 analytical and numerical derivatives of the distance between atoms 1 and 2.
-\verbatim
+\plumedfile
 DISTANCE ATOM=1,2 LABEL=distance
 DISTANCE ATOM=1,2 LABEL=distanceN NUMERICAL_DERIVATIVES
 DUMPDERIVATIVES ARG=distance,distanceN STRIDE=1 FILE=deriv
-\endverbatim
+\endplumedfile
 
 (See also \ref DISTANCE)
 
diff --git a/src/generic/DumpForces.cpp b/src/generic/DumpForces.cpp
index f5e204d65b8297973a50099f9b74797885d788ac..a376de281a1cacc6b5a0b9dd7a0841aa10046309 100644
--- a/src/generic/DumpForces.cpp
+++ b/src/generic/DumpForces.cpp
@@ -43,12 +43,10 @@ by specifying more than one argument. You can control the buffering of output us
 \par Examples
 The following input instructs plumed to write a file called forces that contains
 the force acting on the distance between atoms 1 and 2.
-\verbatim
+\plumedfile
 DISTANCE ATOM=1,2 LABEL=distance
 DUMPFORCES ARG=distance STRIDE=1 FILE=forces
-\endverbatim
-
-(See also \ref DISTANCE)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/generic/DumpMassCharge.cpp b/src/generic/DumpMassCharge.cpp
index 4153bc83844773ab00f19ec811be5d704a9f4b3c..25224f3b3c5d643d9e721e5d215a5e874a8df80f 100644
--- a/src/generic/DumpMassCharge.cpp
+++ b/src/generic/DumpMassCharge.cpp
@@ -49,14 +49,13 @@ masses for all atoms are written.
 You can add the DUMPMASSCHARGE action at the end of the plumed.dat
 file that you use during an MD simulations:
 
-\verbatim
+\plumedfile
 c1: COM ATOMS=1-10
 c2: COM ATOMS=11-20
 PRINT ARG=c1,c2 FILE=colvar STRIDE=100
 
 DUMPMASSCHARGE FILE=mcfile
-\endverbatim
-(see also \ref COM and \ref PRINT)
+\endplumedfile
 
 In this way, you will be able to use the same masses while processing
 a trajectory from the \ref driver . To do so, you need to
@@ -67,11 +66,11 @@ plumed driver --mc mcfile --plumed plumed.dat --ixyz traj.xyz
 
 With the following input you can dump only the charges for a specific
 group.
-\verbatim
+\plumedfile
 solute_ions: GROUP ATOMS=1-121,200-2012
 DUMPATOMS FILE=traj.gro ATOMS=solute_ions STRIDE=100
 DUMPMASSCHARGE FILE=mcfile ATOMS=solute_ions
-\endverbatim
+\endplumedfile
 Notice however that if you want to process the charges
 with the driver (e.g. reading traj.gro) you have to fix atom
 numbers first, e.g. with the script
diff --git a/src/generic/DumpProjections.cpp b/src/generic/DumpProjections.cpp
index 3f0d404297be87078fe7e833f8b5efede814cde5..c3a0ca7509104403884c8ca07bcd296681d14507 100644
--- a/src/generic/DumpProjections.cpp
+++ b/src/generic/DumpProjections.cpp
@@ -39,12 +39,12 @@ Dump the derivatives with respect to the input parameters for one or more object
 Compute the distance between two groups and write on a file the
 derivatives of this distance with respect to all the atoms of the two groups
 
-\verbatim
+\plumedfile
 x1: CENTER ATOMS=1-10
 x2: CENTER ATOMS=11-20
 d: DISTANCE ATOMS=x1,x2
 DUMPPROJECTIONS ARG=d FILE=proj STRIDE=20
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/generic/EffectiveEnergyDrift.cpp b/src/generic/EffectiveEnergyDrift.cpp
index 1c3961390aa08a733668048356ffd9a75809a718..a9ca15e6c8e122d836072fb2924b7f5940b0d50c 100644
--- a/src/generic/EffectiveEnergyDrift.cpp
+++ b/src/generic/EffectiveEnergyDrift.cpp
@@ -55,18 +55,18 @@ Print the effective energy drift described in Ref \cite Ferrarotti2015
 This is to monitor the effective energy drift for a metadynamics
 simulation on the Debye-Huckel energy. Since this variable is very expensive,
 it could be conveniently computed every second step.
-\verbatim
+\plumedfile
 dh: DHENERGY GROUPA=1-10 GROUPB=11-20 EPSILON=80.0 I=0.1 TEMP=300.0
 METAD ARG=dh HEIGHT=0.5 SIGMA=0.1 PACE=500 STRIDE=2
 EFFECTIVE_ENERGY_DRIFT PRINT_STRIDE=100 FILE=eff
-\endverbatim
+\endplumedfile
 
 This is to monitor if a restraint is too stiff
-\verbatim
+\plumedfile
 d: DISTANCE ATOMS=10,20
 RESTRAINT ARG=d KAPPA=100000 AT=0.6
 EFFECTIVE_ENERGY_DRIFT PRINT_STRIDE=100 FILE=eff
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/generic/FitToTemplate.cpp b/src/generic/FitToTemplate.cpp
index 37e38333f64e5190b6f5b658a27ed920d93a4737..6fa379e1fd37dd7227a19f600d3f08673b0b6dda 100644
--- a/src/generic/FitToTemplate.cpp
+++ b/src/generic/FitToTemplate.cpp
@@ -78,13 +78,12 @@ this action is performed at every MD step.
 \par Examples
 
 Align the atomic position to a template then print them
-\verbatim
+\plumedfile
 # to see the effect, one could dump the atoms before alignment
 DUMPATOMS FILE=dump-before.xyz ATOMS=1-20
 FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=SIMPLE
 DUMPATOMS FILE=dump-after.xyz ATOMS=1-20
-\endverbatim
-(see also \ref DUMPATOMS)
+\endplumedfile
 
 
 
diff --git a/src/generic/Flush.cpp b/src/generic/Flush.cpp
index 78c00f87d5bbaf0169c36eaa3ca6e0f2ac6f26cd..71e9a06d3e144e991373ce501fa3d450205caec4 100644
--- a/src/generic/Flush.cpp
+++ b/src/generic/Flush.cpp
@@ -41,7 +41,7 @@ plumed input file, it will flush all the open files.
 
 \par Examples
 A command like this in the input will instruct plumed to flush all the output files every 100 steps
-\verbatim
+\plumedfile
 d1: DISTANCE ATOMS=1,10
 PRINT ARG=d1 STRIDE=5 FILE=colvar1
 
@@ -50,7 +50,7 @@ FLUSH STRIDE=100
 d2: DISTANCE ATOMS=2,11
 # also this print is flushed every 100 steps:
 PRINT ARG=d2 STRIDE=10 FILE=colvar2
-\endverbatim
+\endplumedfile
 (see also \ref DISTANCE and \ref PRINT).
 */
 //+ENDPLUMEDOC
diff --git a/src/generic/Group.cpp b/src/generic/Group.cpp
index 80393a9ce3924e73438f19113cb0212058a2842d..0f04cf75fdf9b93d6aee33cdb1370d601c2d0934 100644
--- a/src/generic/Group.cpp
+++ b/src/generic/Group.cpp
@@ -65,7 +65,7 @@ the \ref INCLUDE command so as to store long group definitions in a separate fil
 
 This command create a group of atoms containing atoms 1,4,7,11 and 14 (labeled 'o'), and another containing
 atoms 2,3,5,6,8,9,12,13 (labeled 'h'):
-\verbatim
+\plumedfile
 o: GROUP ATOMS=1,4,7,11,14
 h: GROUP ATOMS=2,3,5,6,8,9,12,13
 # compute the coordination among the two groups
@@ -75,45 +75,41 @@ c: COORDINATION GROUPA=o GROUPB=h R_0=0.3
 
 # print the coordination on file 'colvar'
 PRINT ARG=c FILE=colvar
-\endverbatim
-(see also \ref COORDINATION and \ref PRINT)
+\endplumedfile
 
 Groups can be conveniently stored in a separate file.
 E.g. one could create a file named 'groups.dat' which reads
-\verbatim
+\plumedfile
 o: GROUP ATOMS=1,4,7,11,14
 h: GROUP ATOMS=2,3,5,6,8,9,12,13
-\endverbatim
+\endplumedfile
 and then include it in the main 'plumed.dat' file
-\verbatim
+\plumedfile
 INCLUDE FILE=groups.dat
 # compute the coordination among the two groups
 c: COORDINATION GROUPA=o GROUPB=h R_0=0.3
 # print the coordination on file 'colvar'
 PRINT ARG=c FILE=colvar
-\endverbatim
-(see also \ref INCLUDE, \ref COORDINATION, and \ref PRINT).
+\endplumedfile
 The groups.dat file could be very long and include lists of thousand atoms without cluttering the main plumed.dat file.
 
 A GROMACS index file can also be imported
-\verbatim
+\plumedfile
 # import group named 'protein' from file index.ndx
 pro: GROUP NDX_FILE=index.ndx NDX_GROUP=protein
 # dump all the atoms of the protein on a trajectory file
 DUMPATOMS ATOMS=pro FILE=traj.gro
-\endverbatim
-(see also \ref DUMPATOMS)
+\endplumedfile
 
 A list can be edited with REMOVE
-\verbatim
+\plumedfile
 # take one atom every three
 ox: GROUP ATOMS=1-90:3
 # take the remaining atoms
 hy: GROUP ATOMS=1-90 REMOVE=ox
 DUMPATOMS ATOMS=ox FILE=ox.gro
 DUMPATOMS ATOMS=hy FILE=hy.gro
-\endverbatim
-(see also \ref DUMPATOMS)
+\endplumedfile
 
 
 */
diff --git a/src/generic/Include.cpp b/src/generic/Include.cpp
index 664c84ba67407f6f44055137f645a09b072940c4..30539f0bb05b5529a60dac5af9e219ad63e31a76 100644
--- a/src/generic/Include.cpp
+++ b/src/generic/Include.cpp
@@ -38,25 +38,25 @@ Useful to split very large plumed.dat files.
 \par Examples
 
 This input
-\verbatim
+\plumedfile
 c1: COM ATOMS=1-100
 c2: COM ATOMS=101-202
 d: DISTANCE ARG=c1,c2
 PRINT ARG=d
-\endverbatim
+\endplumedfile
 
 can be replaced with
-\verbatim
+\plumedfile
 INCLUDE FILE=pippo.dat
 d: DISTANCE ARG=c1,c2
 PRINT ARG=d
-\endverbatim
+\endplumedfile
 
 where the content of file pippo.dat is
-\verbatim
+\plumedfile
 c1: COM ATOMS=1-100
 c2: COM ATOMS=101-202
-\endverbatim
+\endplumedfile
 
 (see also \ref COM, \ref DISTANCE, and \ref PRINT).
 
diff --git a/src/generic/Print.cpp b/src/generic/Print.cpp
index 2a87d1f7bacef394bb5b5407db01311f8dfdfd2c..74245cf99a3c7ffc2437d2369d6486019e9e68a9 100644
--- a/src/generic/Print.cpp
+++ b/src/generic/Print.cpp
@@ -40,13 +40,12 @@ to different files.  You can control the buffering of output using the \subpage
 The following input instructs plumed to print the distance between atoms 3 and 5 on a file
 called COLVAR every 10 steps, and the distance and total energy on a file called COLVAR_ALL
 every 1000 steps.
-\verbatim
+\plumedfile
 DISTANCE ATOMS=2,5 LABEL=distance
 ENERGY             LABEL=energy
 PRINT ARG=distance          STRIDE=10   FILE=COLVAR
 PRINT ARG=distance,energy   STRIDE=1000 FILE=COLVAR_ALL
-\endverbatim
-(See also \ref DISTANCE and \ref ENERGY).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/generic/RandomExchanges.cpp b/src/generic/RandomExchanges.cpp
index be36873bccf097a317e99aaf1270792240a7c319..60d2071207366b50bd2b883df290f33b0172e7ba 100644
--- a/src/generic/RandomExchanges.cpp
+++ b/src/generic/RandomExchanges.cpp
@@ -44,25 +44,25 @@ metadynamics simulation using a different angle in each replica.
 Exchanges will be randomly tried between replicas 0-1, 0-2 and 1-2
 
 Here is plumed.0.dat
-\verbatim
+\plumedfile
 RANDOM_EXCHANGES
 t: TORSION ATOMS=1,2,3,4
 METAD ARG=t HEIGHT=0.1 PACE=100 SIGMA=0.3
-\endverbatim
+\endplumedfile
 
 Here is plumed.1.dat
-\verbatim
+\plumedfile
 RANDOM_EXCHANGES
 t: TORSION ATOMS=2,3,4,5
 METAD ARG=t HEIGHT=0.1 PACE=100 SIGMA=0.3
-\endverbatim
+\endplumedfile
 
 Here is plumed.2.dat
-\verbatim
+\plumedfile
 RANDOM_EXCHANGES
 t: TORSION ATOMS=3,4,5,6
 METAD ARG=t HEIGHT=0.1 PACE=100 SIGMA=0.3
-\endverbatim
+\endplumedfile
 
 \warning Multi replica simulations are presently only working with gromacs.
 
diff --git a/src/generic/Read.cpp b/src/generic/Read.cpp
index 88c279630021b47e6fd1e3e3ea4aa904e161ac69..63b8f8fc801157b83b460b1f0010bf777002c8d8 100644
--- a/src/generic/Read.cpp
+++ b/src/generic/Read.cpp
@@ -52,11 +52,11 @@ This input reads in data from a file called input_colvar.data that was generated
 in a calculation that involved PLUMED.  The first command reads in the data from the
 column headed phi1 while the second reads in the data from the column headed phi2.
 
-\verbatim
+\plumedfile
 rphi1:       READ FILE=input_colvar.data  VALUES=phi1
 rphi2:       READ FILE=input_colvar.data  VALUES=phi2
 PRINT ARG=rphi1,rphi2 STRIDE=500  FILE=output_colvar.data
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/generic/ResetCell.cpp b/src/generic/ResetCell.cpp
index fd14c1fe102f5a873c6ce5d4623df4e0b2b94100..5268d52de348b3271e0b076a072c32d480d09e69 100644
--- a/src/generic/ResetCell.cpp
+++ b/src/generic/ResetCell.cpp
@@ -71,14 +71,13 @@ this action is performed at every MD step.
 \par Examples
 
 Reset cell to be triangular after a rototranslational fit
-\verbatim
+\plumedfile
 DUMPATOMS FILE=dump-original.xyz ATOMS=1-20
 FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=OPTIMAL
 DUMPATOMS FILE=dump-fit.xyz ATOMS=1-20
 RESET_CELL TYPE=TRIANGULAR
 DUMPATOMS FILE=dump-reset.xyz ATOMS=1-20
-\endverbatim
-(see also \ref DUMPATOMS)
+\endplumedfile
 
 
 */
diff --git a/src/generic/Time.cpp b/src/generic/Time.cpp
index 3512cb53da5c7f1a1fac71517c3712a8b16aa115..9a821f0b65dcf0a1318d6bb19ed5f005527f340a 100644
--- a/src/generic/Time.cpp
+++ b/src/generic/Time.cpp
@@ -35,11 +35,10 @@ retrieve the time of the simulation to be used elsewere
 
 \par Examples
 
-\verbatim
+\plumedfile
 TIME            LABEL=t1
 PRINT ARG=t1
-\endverbatim
-(See also \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/generic/UpdateIf.cpp b/src/generic/UpdateIf.cpp
index 2d067fcf3db179ca6ce259b41a611dfd262908fc..696e66c20f4eb423b661ba3214fcf2b68ca8953e 100644
--- a/src/generic/UpdateIf.cpp
+++ b/src/generic/UpdateIf.cpp
@@ -62,7 +62,7 @@ can lead to unexpected results.
 \par Examples
 The following input instructs plumed dump all the snapshots where an atom is in touch with
 the solute.
-\verbatim
+\plumedfile
 solute: GROUP ATOMS=1-124
 coord: COORDINATION GROUPA=solute GROUPB=500 R_0=0.5
 
@@ -72,8 +72,7 @@ coord: COORDINATION GROUPA=solute GROUPB=500 R_0=0.5
 UPDATE_IF ARG=coord MORE_THAN=0.5
 DUMPATOMS ATOMS=solute,500 FILE=output.xyz
 UPDATE_IF ARG=coord END
-\endverbatim
-(See also \ref GROUP, \ref COORDINATION, and \ref DUMPATOMS)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/generic/WholeMolecules.cpp b/src/generic/WholeMolecules.cpp
index 12e9c82813bb10b53540f9598ff4a0304b9adb8d..e43c0aa07b571e2dfa9498146b5f2f4b24dfe7e3 100644
--- a/src/generic/WholeMolecules.cpp
+++ b/src/generic/WholeMolecules.cpp
@@ -74,30 +74,27 @@ to skip some atoms, provided consecute chosen atoms are close enough.
 This command instructs plumed to reconstruct the molecule containing atoms 1-20
 at every step of the calculation and dump them on a file.
 
-\verbatim
+\plumedfile
 # to see the effect, one could dump the atoms as they were before molecule reconstruction:
 # DUMPATOMS FILE=dump-broken.xyz ATOMS=1-20
 WHOLEMOLECULES ENTITY0=1-20
 DUMPATOMS FILE=dump.xyz ATOMS=1-20
-\endverbatim
-(see also \ref DUMPATOMS)
+\endplumedfile
 
 This command instructs plumed to reconstruct two molecules containing atoms 1-20 and 30-40
 
-\verbatim
+\plumedfile
 WHOLEMOLECULES ENTITY0=1-20 ENTITY1=30-40
 DUMPATOMS FILE=dump.xyz ATOMS=1-20,30-40
-\endverbatim
-(see also \ref DUMPATOMS)
+\endplumedfile
 
 This command instructs plumed to reconstruct the chain of backbone atoms in a
 protein
 
-\verbatim
+\plumedfile
 MOLINFO STRUCTURE=helix.pdb
 WHOLEMOLECULES RESIDUES=all MOLTYPE=protein
-\endverbatim
-(See also \ref MOLINFO)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/generic/WrapAround.cpp b/src/generic/WrapAround.cpp
index 946a060fdef55569d473b09559adf427d4da969b..6320ee4bbca6434e9dad48ffdf0a2ac2d9226dcb 100644
--- a/src/generic/WrapAround.cpp
+++ b/src/generic/WrapAround.cpp
@@ -69,20 +69,19 @@ consider the possibility of using the STRIDE keyword here (with great care).
 This command instructs plumed to move all the ions to their periodic image that is as close as possible to
 the rna group.
 
-\verbatim
+\plumedfile
 rna: GROUP ATOMS=1-100
 ions: GROUP ATOMS=101-110
 # first make the rna molecule whole
 WHOLEMOLECULES ENTITY0=rna
 WRAPAROUND ATOMS=ions AROUND=rna
 DUMPATOMS FILE=dump.xyz ATOMS=rna,ions
-\endverbatim
-(see also \ref WHOLEMOLECULES, \ref GROUP and \ref DUMPATOMS)
+\endplumedfile
 
 In case you want to do it during a simulation and you only care about wrapping the ions in
 the `dump.xyz` file, you can use the following:
 
-\verbatim
+\plumedfile
 # add some restraint that do not require molecules to be whole:
 a: TORSION ATOMS=1,2,10,11
 RESTRAINT ARG=a AT=0.0 KAPPA=5
@@ -98,15 +97,14 @@ ions: GROUP ATOMS=101-110
 WHOLEMOLECULES ENTITY0=rna STRIDE=100
 WRAPAROUND ATOMS=ions AROUND=rna STRIDE=100
 DUMPATOMS FILE=dump.xyz ATOMS=rna,ions STRIDE=100
-\endverbatim
-(see also \ref TORSION, \ref GROUP, \ref WHOLEMOLECULES and \ref DUMPATOMS)
+\endplumedfile
 
 Notice that if the biased variable requires a molecule to be whole, you might have to put
 just the \ref WHOLEMOLECULES command before computing that variable and leave the default STRIDE=1.
 
 This command instructs plumed to center all atoms around the center of mass of a solute molecule.
 
-\verbatim
+\plumedfile
 solute: GROUP ATOMS=1-100
 all: GROUP ATOMS=1-1000
 # center of the solute:
@@ -116,8 +114,7 @@ com: COM ATOMS=solute
 # notice that we wrap around a single atom. this should be fast
 WRAPAROUND ATOMS=all AROUND=com
 DUMPATOMS FILE=dump.xyz ATOMS=all
-\endverbatim
-(see also \ref GROUP \ref COM \ref DUMPATOMS)
+\endplumedfile
 
 Notice that whereas \ref WHOLEMOLECULES is designed to make molecules whole,
 \ref WRAPAROUND can easily break molecules. In the last example,
@@ -136,7 +133,7 @@ in the following examples all the water oxygens will be brought
 close to the solute, and all the hydrogens will be kept close
 to their related oxygen.
 
-\verbatim
+\plumedfile
 solute: GROUP ATOMS=1-100
 water: GROUP ATOMS=101-1000
 com: COM ATOMS=solute
@@ -145,8 +142,7 @@ WRAPAROUND ATOMS=solute AROUND=com
 # notice that we wrap around a single atom. this should be fast
 WRAPAROUND ATOMS=water AROUND=com GROUPBY=3
 DUMPATOMS FILE=dump.xyz ATOMS=solute,water
-\endverbatim
-(see also \ref GROUP \ref COM \ref DUMPATOMS)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/gridtools/ConvertToFES.cpp b/src/gridtools/ConvertToFES.cpp
index 4e3bbdf1705a63638b0ef08c1fe6658eda4b304a..661aa2bde96e2116adb664bc58c5156d97d0923c 100644
--- a/src/gridtools/ConvertToFES.cpp
+++ b/src/gridtools/ConvertToFES.cpp
@@ -46,12 +46,12 @@ and the HISTOGRAM action.  All the data within this trajectory is used in the co
 HISTOGRAM.  Finally, once all the data has been read in, the histogram is converted to a free energy
 using the formula above and the free energy is output to a file called fes.dat
 
-\verbatim
+\plumedfile
 x: DISTANCE ATOMS=1,2
 hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
 ff: CONVERT_TO_FES GRID=hA1 TEMP=300
 DUMPGRID GRID=ff FILE=fes.dat
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/gridtools/DumpCube.cpp b/src/gridtools/DumpCube.cpp
index 0fbc47f566575199de0688a5277fb391d96efcb2..f1df7406e75d598937d67bd313b0bdfcd3715ce4 100644
--- a/src/gridtools/DumpCube.cpp
+++ b/src/gridtools/DumpCube.cpp
@@ -46,14 +46,14 @@ between atoms 1 and 2, the distance between atom 1 and 3 and the angle between t
 all the kernels have been added the resulting histogram is output to a file called histoA1.cube.  This file has the
 Gaussian cube file format.  The histogram can thus be visualized using tools such as VMD.
 
-\verbatim
+\plumedfile
 x1: DISTANCE ATOMS=1,2
 x2: DISTANCE ATOMS=1,3
 x3: ANGLE ATOMS=1,2,3
 
 hA1: HISTOGRAM ARG=x1,x2,x3 GRID_MIN=0.0,0.0,0.0 GRID_MAX=3.0,3.0,3.0 GRID_BIN=10,10,10 BANDWIDTH=1.0,1.0,1.0
 DUMPCUBE GRID=hA1 FILE=histoA1.cube
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/gridtools/DumpGrid.cpp b/src/gridtools/DumpGrid.cpp
index 608742a8bfb48a1b9395b7906ac21363c42aa4e1..e3524d563d67aa96c091fb389984959385c6a9de 100644
--- a/src/gridtools/DumpGrid.cpp
+++ b/src/gridtools/DumpGrid.cpp
@@ -67,7 +67,7 @@ for y.  This block is then followed by a blank line again and this pattern conti
 
 The following input monitors two torsional angles during a simulation
 and outputs a continuos histogram as a function of them at the end of the simulation.
-\verbatim
+\plumedfile
 TORSION ATOMS=1,2,3,4 LABEL=r1
 TORSION ATOMS=2,3,4,5 LABEL=r2
 HISTOGRAM ...
@@ -80,11 +80,11 @@ HISTOGRAM ...
 ... HISTOGRAM
 
 DUMPGRID GRID=hh FILE=histo
-\endverbatim
+\endplumedfile
 
 The following input monitors two torsional angles during a simulation
 and outputs a discrete histogram as a function of them at the end of the simulation.
-\verbatim
+\plumedfile
 TORSION ATOMS=1,2,3,4 LABEL=r1
 TORSION ATOMS=2,3,4,5 LABEL=r2
 HISTOGRAM ...
@@ -98,11 +98,11 @@ HISTOGRAM ...
 ... HISTOGRAM
 
 DUMPGRID GRID=hh FILE=histo
-\endverbatim
+\endplumedfile
 
 The following input monitors two torsional angles during a simulation
 and outputs the histogram accumulated thus far every 100000 steps.
-\verbatim
+\plumedfile
 TORSION ATOMS=1,2,3,4 LABEL=r1
 TORSION ATOMS=2,3,4,5 LABEL=r2
 HISTOGRAM ...
@@ -115,14 +115,14 @@ HISTOGRAM ...
 ... HISTOGRAM
 
 DUMPGRID GRID=hh FILE=histo STRIDE=100000
-\endverbatim
+\endplumedfile
 
 The following input monitors two torsional angles during a simulation
 and outputs a separate histogram for each 100000 steps worth of trajectory.
 Notice how the CLEAR keyword is used here and how it is not used in the
 previous example.
 
-\verbatim
+\plumedfile
 TORSION ATOMS=1,2,3,4 LABEL=r1
 TORSION ATOMS=2,3,4,5 LABEL=r2
 HISTOGRAM ...
@@ -136,7 +136,7 @@ HISTOGRAM ...
 ... HISTOGRAM
 
 DUMPGRID GRID=hh FILE=histo STRIDE=100000
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/gridtools/FindContour.cpp b/src/gridtools/FindContour.cpp
index a58e4f650e4dddf7a05733ae729cb9d77ab700f5..1bae86d431cce76f9d7b98b24420bdcec41292b2 100644
--- a/src/gridtools/FindContour.cpp
+++ b/src/gridtools/FindContour.cpp
@@ -67,7 +67,7 @@ solidness at each point in the simulation cell.  The Willard-Chandler dividing s
 at which the value of this phase field is equal to 0.5.  This set of points is output to file called mycontour.dat.  A new contour
 is found on every single step for each frame that is read in.
 
-\verbatim
+\plumedfile
 UNITS NATURAL
 FCCUBIC ...
   SPECIES=1-96000 SWITCH={CUBIC D_0=1.2 D_MAX=1.5}
@@ -84,7 +84,7 @@ MULTICOLVARDENS ...
 ... MULTICOLVARDENS
 
 FIND_CONTOUR GRID=dens CONTOUR=0.5 FILE=mycontour.dat
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/gridtools/FindSphericalContour.cpp b/src/gridtools/FindSphericalContour.cpp
index a2ee2aa8627c0a1fb231f79111d936d5e5e25ffe..7d36e818c882e99fcba0085a6409093cf8d7b425 100644
--- a/src/gridtools/FindSphericalContour.cpp
+++ b/src/gridtools/FindSphericalContour.cpp
@@ -81,7 +81,7 @@ by exploiting the functionality within \ref CENTER_OF_MULTICOLVAR.  We can then
 cluster and the values of the coordination numbers of these atoms.  The final line in the input then finds the a set of points on the dividing surface that separates
 teh droplet from the surrounding gas.  The value of the phase field on this isocontour is equal to 0.75.
 
-\verbatim
+\plumedfile
 # Calculate coordination numbers
 c1: COORDINATIONNUMBER SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
 # Select coordination numbers that are more than 2.0
@@ -98,7 +98,7 @@ cent: CENTER_OF_MULTICOLVAR DATA=trans1
 dens: MULTICOLVARDENS DATA=trans1 ORIGIN=cent DIR=xyz NBINS=30,30,30 BANDWIDTH=2.0,2.0,2.0
 # Find the isocontour around the nucleus
 FIND_SPHERICAL_CONTOUR GRID=dens CONTOUR=0.85 INNER_RADIUS=10.0 OUTER_RADIUS=40.0 FILE=mysurface.xyz UNITS=A PRECISION=4 NPOINTS=100
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/gridtools/FourierTransform.cpp b/src/gridtools/FourierTransform.cpp
index 5a53b390dec72635b5c802995b12635942260a58..87a078336ebc5f7a4dfd2d7a8426a16d4be91386 100644
--- a/src/gridtools/FourierTransform.cpp
+++ b/src/gridtools/FourierTransform.cpp
@@ -57,9 +57,9 @@ The default values of these parameters are: \f$a=1\f$ and \f$b=1\f$.
 \par Examples
 The following example tells Plumed to compute the complex 2D 'backward' Discrete Fourier Transform by taking the data saved on a grid called 'density', and normalizing the output by \f$ \frac{1}{\sqrt{N_x\, N_y}}\f$, where \f$N_x\f$ and \f$N_y\f$ are the number of data on the grid (it can be the case that \f$N_x \neq N_y\f$):
 
-\verbatim
+\plumedfile
 FOURIER_TRANSFORM STRIDE=1 GRID=density FT_TYPE=complex FOURIER_PARAMETERS=0,-1 FILE=fourier.dat
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/gridtools/InterpolateGrid.cpp b/src/gridtools/InterpolateGrid.cpp
index e8f97c0d4abe25e86dff08372a69e63bdc5b62cd..7fbc1d65d334757ce8e712cc59345f6202eb878f 100644
--- a/src/gridtools/InterpolateGrid.cpp
+++ b/src/gridtools/InterpolateGrid.cpp
@@ -38,12 +38,12 @@ distance between atoms 1 and 2 using kernel density estimation.  During the calc
 are evaluated at 100 points on a uniform grid between 0.0 and 3.0.  Prior to outputting this function at the end of the
 simulation this function is interpolated onto a finer grid of 200 points between 0.0 and 3.0.
 
-\verbatim
+\plumedfile
 x: DISTANCE ATOMS=1,2
 hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
 ii: INTERPOLATE_GRID GRID=hA1 GRID_BIN=200
 DUMPGRID GRID=ii FILE=histo.dat
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/maketools/plumedcheck b/src/maketools/plumedcheck
index 120808decbabe905049f6c393a0cfbf14e68cfa3..5886ec086329d18d12c27278b206160a4c6c2869 100755
--- a/src/maketools/plumedcheck
+++ b/src/maketools/plumedcheck
@@ -354,6 +354,12 @@ BEGINFILE{
     if(in_plumed_doc && (in_plumed_doc in provide_examples) && match($0,"^ *\\\\verbatim *$")){
       provide_verbatim[in_plumed_doc]=1
     }
+# check if, besides the \par Examples text, there is a plumedfile command
+# now it is considered equivalent to a verbatim
+# later on we might make plumedfile compulsory instead of verbatim
+    if(in_plumed_doc && (in_plumed_doc in provide_examples) && match($0,"^ *\\\\plumedfile *$")){
+      provide_verbatim[in_plumed_doc]=1
+    }
 #print match($0,"^ *\\\\verbatim *$"),"X" $0 "X"
 
 # take note of actions or cltools that provide examples
diff --git a/src/manyrestraints/LWalls.cpp b/src/manyrestraints/LWalls.cpp
index 7c0ad569ecb1ba59a48b18616bae7292086ab711..c98405d04e40d3a232d052f43ac38f3df99662c3 100644
--- a/src/manyrestraints/LWalls.cpp
+++ b/src/manyrestraints/LWalls.cpp
@@ -45,10 +45,10 @@ The following set of commands can be used to stop any of the 800 atoms in group
 in the z direction from atom 34137.  This is done by adding a lower wall on the z-distance between all the atoms in
 group A and the position of 34137.
 
-\verbatim
+\plumedfile
 l: ZDISTANCES GROUPA=1-800 GROUPB=34137 NOPBC
 LWALLS DATA=l AT=2.46465 KAPPA=150.0 EXP=2 EPS=1 OFFSET=0 LABEL=lwall
-\endverbatim
+\endplumedfile
 
 
 */
diff --git a/src/manyrestraints/UWalls.cpp b/src/manyrestraints/UWalls.cpp
index f4c75f0d8976ac05aed11169d607fe16040c9b76..2c0f75367da490b1cba34e54896379e47a4eddd5 100644
--- a/src/manyrestraints/UWalls.cpp
+++ b/src/manyrestraints/UWalls.cpp
@@ -48,11 +48,11 @@ and the center of mass of the cluster.  These distances are then passed to the U
 a \ref UPPER_WALLS restraint on each of them and thereby prevents each of them from moving very far from the centre
 of mass of the cluster.
 
-\verbatim
+\plumedfile
 COM ATOMS=1-20 LABEL=c1
 DISTANCES GROUPA=c1 GROUPB=1-20 LABEL=d1
 UWALLS DATA=d1 AT=2.5 KAPPA=0.2 LABEL=sr
-\endverbatim
+\endplumedfile
 
 
 */
diff --git a/src/mapping/PCAVars.cpp b/src/mapping/PCAVars.cpp
index d029720ed843f89007cb6e344684bcd1ab06b1ec..20ae003016ae482155fc3d8ee446e8a7813fb829 100644
--- a/src/mapping/PCAVars.cpp
+++ b/src/mapping/PCAVars.cpp
@@ -76,10 +76,10 @@ frame are removed from these displacements.  The matrix \f$A\f$ and the referenc
 configuration \f$R^{ref}\f$ are specified in the pdb input file reference.pdb and the
 value of all projections (and the residual) are output to a file called colvar2.
 
-\verbatim
+\plumedfile
 PCAVARS REFERENCE=reference.pdb TYPE=OPTIMAL LABEL=pca2
 PRINT ARG=pca2.* FILE=colvar2
-\endverbatim
+\endplumedfile
 
 The reference configurations can be specified using a pdb file.  The first configuration that you provide is the reference configuration,
 which is refered to in the above as \f$X^{ref}\f$ subsequent configurations give the directions of row vectors that are contained in
diff --git a/src/mapping/Path.cpp b/src/mapping/Path.cpp
index 37db777d2e1677f2fa75198affee22c484dbd8e8..4808fad15dd52176845faa3afea24d097c0caad0 100644
--- a/src/mapping/Path.cpp
+++ b/src/mapping/Path.cpp
@@ -73,21 +73,21 @@ In the example below the path is defined using RMSD distance from frames.
 The reference frames in the path are defined in the pdb file.  In this frame
 each configuration in the path is separated by a line containing just the word END.
 
-\verbatim
+\plumedfile
 p1: PATH REFERENCE=file.pdb TYPE=OPTIMAL LAMBDA=500.0
 PRINT ARG=p1.sss,p1.zzz STRIDE=1 FILE=colvar FMT=%8.4f
-\endverbatim
+\endplumedfile
 
 In the example below the path is defined using the values of two torsional angles (t1 and t2).
 In addition, the \f$s\f$ and \f$z\f$ are calculated using the geometric expressions described
 above rather than the alegebraic expressions that are used by default.
 
-\verbatim
+\plumedfile
 t1: TORSION ATOMS=5,7,9,15
 t2: TORSION ATOMS=7,9,15,17
 pp: PATH TYPE=EUCLIDEAN REFERENCE=epath.pdb GPATH NOSPATH NOZPATH
 PRINT ARG=pp.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 Notice that the LAMBDA parameter is not required here as we are not calculating \f$s\f$ and \f$s\f$
 using the algebraic formulae defined earlier.  The positions of the frames in the path are defined
@@ -109,19 +109,19 @@ The following input instructs PLUMED to calculate the values of the path collect
 path are defined in the file all.pdb and all distances are measured using the OPTIMAL metric that is discussed in the manual
 page on \ref RMSD.
 
-\verbatim
+\plumedfile
 p2: PATH REFERENCE=all.pdb LAMBDA=69087
 PRINT ARG=p2.spath,p2.zpath STRIDE=1 FILE=colvar
-\endverbatim
+\endplumedfile
 
 If you wish to use collective variable values in the definition of your path you would use an input file with something like this:
 
-\verbatim
+\plumedfile
 d1: DISTANCE ATOMS=1,2
 d2: DISTANCE ATOMS=3,4a
 p2: PATH REFERENCE=mypath.pdb LAMBDA=2 TYPE=EUCLIDEAN
 PRINT ARG=p2.spath,p2.zpath STRIDE=1 FILE=colvar
-\endverbatim
+\endplumedfile
 
 The corresponding pdb file containing the  definitions of the frames in the path would then look like this:
 
diff --git a/src/mapping/PropertyMap.cpp b/src/mapping/PropertyMap.cpp
index 9da32b47bbad984703c06559560ec992115d27a5..91c96f4dd650734700561adadb19a868b2fd1ce1 100644
--- a/src/mapping/PropertyMap.cpp
+++ b/src/mapping/PropertyMap.cpp
@@ -47,10 +47,10 @@ that these properties take at a set of reference configurations and using the fo
 between the reference configurations and the instantaneous configurations are calculated using the OPTIMAL metric that is
 discussed at length in the manual pages on \ref RMSD.
 
-\verbatim
+\plumedfile
 p2: GPROPERTYMAP REFERENCE=allv.pdb PROPERTY=X,Y LAMBDA=69087
 PRINT ARG=p2.X,p2.Y,p2.zpath STRIDE=1 FILE=colvar
-\endverbatim
+\endplumedfile
 
 The additional input file for this calculation, which contains the reference frames and the values of X and Y at these reference
 points has the following format.
diff --git a/src/multicolvar/AlphaBeta.cpp b/src/multicolvar/AlphaBeta.cpp
index 01379c58374c04184de72ee675fef4c8b338b2d5..aea79fc2ea288a922e5897281e718f759902cfc1 100644
--- a/src/multicolvar/AlphaBeta.cpp
+++ b/src/multicolvar/AlphaBeta.cpp
@@ -49,7 +49,7 @@ The \f$\phi_i^{\textrm{Ref}}\f$ values are the user-specified reference values f
 
 The following provides an example of the input for an alpha beta similarity.
 
-\verbatim
+\plumedfile
 ALPHABETA ...
 ATOMS1=168,170,172,188 REFERENCE1=3.14
 ATOMS2=170,172,188,190 REFERENCE2=3.14
@@ -57,11 +57,11 @@ ATOMS3=188,190,192,230 REFERENCE3=3.14
 LABEL=ab
 ... ALPHABETA
 PRINT ARG=ab FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 
 Because all the reference values are the same we can calculate the same quantity using
 
-\verbatim
+\plumedfile
 ALPHABETA ...
 ATOMS1=168,170,172,188 REFERENCE=3.14
 ATOMS2=170,172,188,190
@@ -69,13 +69,13 @@ ATOMS3=188,190,192,230
 LABEL=ab
 ... ALPHABETA
 PRINT ARG=ab FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 
 Writing out the atoms involved in all the torsions in this way can be rather tedious. Thankfully if you are working with protein you
 can avoid this by using the \ref MOLINFO command.  PLUMED uses the pdb file that you provide to this command to learn
 about the topology of the protein molecule.  This means that you can specify torsion angles using the following syntax:
 
-\verbatim
+\plumedfile
 MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
 ALPHABETA ...
 ATOMS1=@phi-3 REFERENCE=3.14
@@ -84,7 +84,7 @@ ATOMS3=@phi-4
 LABEL=ab
 ... ALPHABETA
 PRINT ARG=ab FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 
 Here, \@phi-3 tells plumed that you would like to calculate the \f$\phi\f$ angle in the third residue of the protein.
 Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle of the 4th residue of the protein.
diff --git a/src/multicolvar/Angles.cpp b/src/multicolvar/Angles.cpp
index a42233cbd4de4fbdfef9942dc4aa674e6522256d..cdc062579991198e890a2aeed8675dd69cf33ebc 100644
--- a/src/multicolvar/Angles.cpp
+++ b/src/multicolvar/Angles.cpp
@@ -58,28 +58,28 @@ an atom / molecule \cite lj-recon.
 The following example instructs plumed to find the average of two angles and to
 print it to a file
 
-\verbatim
+\plumedfile
 ANGLES ATOMS1=1,2,3 ATOMS2=4,5,6 MEAN LABEL=a1
 PRINT ARG=a1.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following example tells plumed to calculate all angles involving
 at least one atom from GROUPA and two atoms from GROUPB in which the distances
 are less than 1.0. The number of angles between \f$\frac{\pi}{4}\f$ and
 \f$\frac{3\pi}{4}\f$ is then output
 
-\verbatim
+\plumedfile
 ANGLES GROUPA=1-10 GROUPB=11-100 BETWEEN={GAUSSIAN LOWER=0.25pi UPPER=0.75pi} SWITCH={GAUSSIAN R_0=1.0} LABEL=a1
 PRINT ARG=a1.between FILE=colvar
-\endverbatim
+\endplumedfile
 
 This final example instructs plumed to calculate all the angles in the first coordination
 spheres of the atoms. A discretized-normalized histogram of the distribution is then output
 
-\verbatim
+\plumedfile
 ANGLES GROUP=1-38 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=pi NBINS=20} SWITCH={GAUSSIAN R_0=1.0} LABEL=a1
 PRINT ARG=a1.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/Bridge.cpp b/src/multicolvar/Bridge.cpp
index 42198ffa1e26ddebc91b2a840dd0a40f92c22a22..53e76158b5f6339f3da4b666ff2bdfaba340ec91 100644
--- a/src/multicolvar/Bridge.cpp
+++ b/src/multicolvar/Bridge.cpp
@@ -51,10 +51,10 @@ The following example instructs plumed to calculate the number of water molecule
 that are bridging betweeen atoms 1-10 and atoms 11-20 and to print the value
 to a file
 
-\verbatim
+\plumedfile
 BRIDGE BRIDGING_ATOMS=100-200 GROUPA=1-10 GROUPB=11-20 LABEL=w1
 PRINT ARG=a1.mean FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/CenterOfMultiColvar.cpp b/src/multicolvar/CenterOfMultiColvar.cpp
index f8b5be66001a8defaef7b40d443d8d667ce2d7fe..0708bd4025863b288082bb56e17cee3a6e39bd83 100644
--- a/src/multicolvar/CenterOfMultiColvar.cpp
+++ b/src/multicolvar/CenterOfMultiColvar.cpp
@@ -55,10 +55,10 @@ As you want to calculate the position of the droplets you thus recognise that th
 numbers should have a high weight in the weighted average you are using to calculate the position of the droplet.
 You can thus calculate the position of the droplet using an input like the one shown below:
 
-\verbatim
+\plumedfile
 c1: COORDINATIONNUMBER SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5}
 cc: CENTER_OF_MULTICOLVAR DATA=c1
-\endverbatim
+\endplumedfile
 
 The first line here calclates the coordination numbers of all the atoms in the system.  The virtual atom then uses the values
 of the coordination numbers calculated by the action labelled c1 when it calculates the Berry Phase average described above.
@@ -66,11 +66,11 @@ of the coordination numbers calculated by the action labelled c1 when it calcula
 
 The above input is fine we can, however, refine this somewhat by making use of a multicolvar transform action as shown below:
 
-\verbatim
+\plumedfile
 c1: COORDINATIONNUMBER SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5}
 cf: MTRANSFORM_MORE DATA=c1 SWITCH={RATIONAL D_0=2.0 R_0=0.1} LOWMEM
 cc: CENTER_OF_MULTICOLVAR DATA=cf
-\endverbatim
+\endplumedfile
 
 This input once again calculates the coordination numbers of all the atoms in the system.  The middle line then transforms these
 coordinations numbers to numbers between 0 and 1.  Essentially any atom with a coordination number larger than 2.0 is given a weight
diff --git a/src/multicolvar/CoordinationNumbers.cpp b/src/multicolvar/CoordinationNumbers.cpp
index a3823b4887d687ed220cb85eee20efa29fed304b..38c26d348b9486b9c000c68e261a6395af081b3c 100644
--- a/src/multicolvar/CoordinationNumbers.cpp
+++ b/src/multicolvar/CoordinationNumbers.cpp
@@ -48,16 +48,16 @@ s = \frac{ 1 - \left(\frac{r-d_0}{r_0}\right)^n } { 1 - \left(\frac{r-d_0}{r_0}\
 
 The following input tells plumed to calculate the coordination numbers of atoms 1-100 with themselves.
 The minimum coordination number is then calculated.
-\verbatim
+\plumedfile
 COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 MIN={BETA=0.1}
-\endverbatim
+\endplumedfile
 
 The following input tells plumed to calculate how many atoms from 1-100 are within 3.0 of each of the atoms
 from 101-110.  In the first 101 is the central atom, in the second 102 is the central atom and so on.  The
 number of coordination numbers more than 6 is then computed.
-\verbatim
+\plumedfile
 COORDINATIONNUMBER SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=6.0 NN=6 MM=12 D_0=0}
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/Density.cpp b/src/multicolvar/Density.cpp
index 0351140a3f309a0f63d4ebdf387c73188761b676..598e16b30929e759fcc4e69d15ef93f7d13da9c6 100644
--- a/src/multicolvar/Density.cpp
+++ b/src/multicolvar/Density.cpp
@@ -40,11 +40,11 @@ the number of atoms in half the box.
 
 The following example calculates the number of atoms in one half of the simulation box.
 
-\verbatim
+\plumedfile
 DENSITY SPECIES=1-100 LABEL=d
 AROUND ARG=d XLOWER=0.0 XUPPER=0.5 LABEL=d1
 PRINT ARG=d1.* FILE=colvar1 FMT=%8.4f
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/DihedralCorrelation.cpp b/src/multicolvar/DihedralCorrelation.cpp
index e75bde668092da07212a91a7aa4641bd7a4f7406..a1f5ea87295680e3a7cd10a76addfefd4944a6da 100644
--- a/src/multicolvar/DihedralCorrelation.cpp
+++ b/src/multicolvar/DihedralCorrelation.cpp
@@ -48,14 +48,14 @@ where the \f$\phi_i\f$ and \f$\psi\f$ values and the instantaneous values for th
 
 The following provides an example input for the dihcor action
 
-\verbatim
+\plumedfile
 DIHCOR ...
   ATOMS1=1,2,3,4,5,6,7,8
   ATOMS2=5,6,7,8,9,10,11,12
   LABEL=dih
 ... DIHCOR
 PRINT ARG=dih FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 
 In the above input we are calculating the correation between the torsion angle involving atoms 1, 2, 3 and 4 and the torsion angle
 involving atoms 5, 6, 7 and 8.	This is then added to the correlation betwene the torsion angle involving atoms 5, 6, 7 and 8 and the
@@ -65,7 +65,7 @@ Writing out the atoms involved in all the torsions in this way can be rather ted
 can avoid this by using the \ref MOLINFO command.  PLUMED uses the pdb file that you provide to this command to learn
 about the topology of the protein molecule.  This means that you can specify torsion angles using the following syntax:
 
-\verbatim
+\plumedfile
 MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
 DIHCOR ...
 ATOMS1=@phi-3,@psi-3
@@ -73,7 +73,7 @@ ATOMS2=@psi-3,@phi-4
 ATOMS4=@phi-4,@psi-4
 ... DIHCOR
 PRINT ARG=dih FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 
 Here, \@phi-3 tells plumed that you would like to calculate the \f$\phi\f$ angle in the third residue of the protein.
 Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle of the 4th residue of the protein.
diff --git a/src/multicolvar/DistanceFromContour.cpp b/src/multicolvar/DistanceFromContour.cpp
index 214930392e0e2dffaab6c222a02c5043e50172b0..6f99e214a478ff53c3f3c5fea9ed3a059e09cec0 100644
--- a/src/multicolvar/DistanceFromContour.cpp
+++ b/src/multicolvar/DistanceFromContour.cpp
@@ -58,10 +58,10 @@ In this example atoms 2-100 are assumed to be concentraed along some part of the
 an interface between a liquid/solid and the vapour.  The quantity dc measures the distance between the
 surface at which the density of 2-100 atoms is equal to 0.2 and the position of the test particle atom 1.
 
-\verbatim
+\plumedfile
 dens: DENSITY SPECIES=2-100
 dc: DISTANCE_FROM_CONTOUR DATA=dens ATOM=1 BANDWIDTH=0.5,0.5,0.5 DIR=z CONTOUR=0.2
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/Distances.cpp b/src/multicolvar/Distances.cpp
index d5cb8fde0026618b2b4df7c5de24c3c9c38319d9..579668b8807085dc092ac1db6171b4abd2aa00c2 100644
--- a/src/multicolvar/Distances.cpp
+++ b/src/multicolvar/Distances.cpp
@@ -42,36 +42,36 @@ distances such as the minimum, the number less than a certain quantity and so on
 
 The following input tells plumed to calculate the distances between atoms 3 and 5 and between atoms 1 and 2 and to
 print the minimum for these two distances.
-\verbatim
+\plumedfile
 DISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 
 The following input tells plumed to calculate the distances between atoms 3 and 5 and between atoms 1 and 2
 and then to calculate the number of these distances that are less than 0.1 nm.  The number of distances
 less than 0.1nm is then printed to a file.
-\verbatim
+\plumedfile
 DISTANCES ATOMS1=3,5 ATOMS2=1,2 LABEL=d1 LESS_THAN={RATIONAL R_0=0.1}
 PRINT ARG=d1.lt0.1
-\endverbatim
+\endplumedfile
 (See also \ref PRINT \ref switchingfunction).
 
 The following input tells plumed to calculate all the distances between atoms 1, 2 and 3 (i.e. the distances between atoms
 1 and 2, atoms 1 and 3 and atoms 2 and 3).  The average of these distances is then calculated.
-\verbatim
+\plumedfile
 DISTANCES GROUP=1-3 MEAN LABEL=d1
 PRINT ARG=d1.mean
-\endverbatim
+\endplumedfile
 (See also \ref PRINT)
 
 The following input tells plumed to calculate all the distances between the atoms in GROUPA and the atoms in GROUPB.
 In other words the distances between atoms 1 and 2 and the distance between atoms 1 and 3.  The number of distances
 more than 0.1 is then printed to a file.
-\verbatim
+\plumedfile
 DISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1}
 PRINT ARG=d1.gt0.1
-\endverbatim
+\endplumedfile
 (See also \ref PRINT \ref switchingfunction)
 
 
@@ -79,10 +79,10 @@ PRINT ARG=d1.gt0.1
 
 To calculate and print the minimum distance between two groups of atoms you use the following commands
 
-\verbatim
+\plumedfile
 d1: DISTANCES GROUPA=1-10 GROUPB=11-20 MIN={BETA=500.}
 PRINT ARG=d1.min FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 (see \ref DISTANCES and \ref PRINT)
 
 In order to ensure differentiability the minimum is calculated using the following function:
@@ -102,7 +102,7 @@ allow you to calculate multiple functions of a distribution of simple collective
 can calculate the number of distances less than 1.0, the minimum distance, the number of distances more than
 2.0 and the number of distances between 1.0 and 2.0 by using the following command:
 
-\verbatim
+\plumedfile
 DISTANCES ...
  GROUPA=1-10 GROUPB=11-20
  LESS_THAN={RATIONAL R_0=1.0}
@@ -111,7 +111,7 @@ DISTANCES ...
  MIN={BETA=500.}
 ... DISTANCES
 PRINT ARG=d1.lessthan,d1.morethan,d1.between,d1.min FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 (see \ref DISTANCES and \ref PRINT)
 
 A calculation performed this way is fast because the expensive part of the calculation - the calculation of all the distances - is only
diff --git a/src/multicolvar/DumpMultiColvar.cpp b/src/multicolvar/DumpMultiColvar.cpp
index 27029a43690b7b658389a5b8bf465b2fe67fed25..fb08af5b4437eb54b8cb6e07d63e2380f13626ef 100644
--- a/src/multicolvar/DumpMultiColvar.cpp
+++ b/src/multicolvar/DumpMultiColvar.cpp
@@ -49,13 +49,13 @@ In this examples we calculate the distances between the  atoms of the first and
 group and we write them in the file MULTICOLVAR.xyz. For each couple it writes the
 coordinates of their geometric center and their distance.
 
-\verbatim
+\plumedfile
 pos:   GROUP ATOMS=220,221,235,236,247,248,438,439,450,451,534,535
 neg:   GROUP ATOMS=65,68,138,182,185,267,270,291,313,316,489,583,621,711
 DISTANCES GROUPA=pos GROUPB=neg LABEL=slt
 
 DUMPMULTICOLVAR DATA=slt FILE=MULTICOLVAR.xyz
-\endverbatim
+\endplumedfile
 
 (see also \ref DISTANCES)
 
diff --git a/src/multicolvar/FilterBetween.cpp b/src/multicolvar/FilterBetween.cpp
index d6a4187118d5c3a35407aa466f5b2900ebad281a..17ad0be9e5308560121d3d8ea4e78263c3f49cc2 100644
--- a/src/multicolvar/FilterBetween.cpp
+++ b/src/multicolvar/FilterBetween.cpp
@@ -59,22 +59,22 @@ In other words, you are calculating the mean for the transformed colvar.
 The following input gives an example of how a MTRANSFORM_BETWEEN action can be used to duplicate
 functionality that is elsehwere in PLUMED.
 
-\verbatim
+\plumedfile
 DISTANCES ...
  GROUPA=1-10 GROUPB=11-20
  LABEL=d1
 ... DISTANCES
 MTRANSFORM_BETWEEN DATA=d1 LOWER=1.0 UPPER=2.0 SMEAR=0.5
-\endverbatim
+\endplumedfile
 
 In this case you can achieve the same result by using:
 
-\verbatim
+\plumedfile
 DISTANCES ...
  GROUPA=1-10 GROUPB=11-20
  BETWEEN={GAUSSIAN LOWER=1.0 UPPER=2.0}
 ... DISTANCES
-\endverbatim
+\endplumedfile
 (see \ref DISTANCES)
 
 The advantage of MTRANSFORM_BETWEEN comes, however, if you want to use transformed colvars as input
@@ -112,10 +112,10 @@ One is thus calculating the mean for those colvars that are within the range of
 
 The example shown below calculates the mean for those distances that are between 0 and 3 nm in length
 
-\verbatim
+\plumedfile
 DISTANCES GROUPA=1 GROUPB=2-50 MEAN LABEL=d1
 MFILTER_BETWEEN DATA=d1 LOWER=0 UPPER=3.0 SMEAR=0.0001 MEAN LABEL=d4
-\endverbatim
+\endplumedfile
 
 More complicated things can be done by using the label of a filter as input to a new multicolvar as shown
 in the example below.  Here the coordination numbers of all atoms are computed.  The atoms with a coordination
@@ -123,11 +123,11 @@ number between 4 and 6 are then identified using the filter.  This reduced list
 to a second coordination number calculation.  This second coordination number thus measures the number of atoms
 4-6 coordinated atoms each of the 4-6 coordination atoms is bound to.
 
-\verbatim
+\plumedfile
 c1: COORDINATIONNUMBER SPECIES=1-150 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
 cf: MFILTER_BETWEEN DATA=c1 LOWER=4 UPPER=6 SMEAR=0.5 LOWMEM
 c2: COORDINATIONNUMBER SPECIES=cf SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0} MORE_THAN={RATIONAL D_0=2.0 R_0=0.1}
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/FilterLessThan.cpp b/src/multicolvar/FilterLessThan.cpp
index 6363bcbb430b7faf42f6d189d241bf9ff5b820a4..e65d5fcd1872a840a20050ff82521bd8b47199d0 100644
--- a/src/multicolvar/FilterLessThan.cpp
+++ b/src/multicolvar/FilterLessThan.cpp
@@ -51,22 +51,22 @@ In other words, you are calculating the mean for the transformed colvar.
 The following input gives an example of how a MTRANSFORM_LESS action can be used to duplicate
 functionality that is elsehwere in PLUMED.
 
-\verbatim
+\plumedfile
 DISTANCES ...
  GROUPA=1-10 GROUPB=11-20
  LABEL=d1
 ... DISTANCES
 MTRANSFORM_LESS DATA=d1 SWITCH={GAUSSIAN D_0=1.5 R_0=0.00001}
-\endverbatim
+\endplumedfile
 
 In this case you can achieve the same result by using:
 
-\verbatim
+\plumedfile
 DISTANCES ...
  GROUPA=1-10 GROUPB=11-20
  LESS_THAN={GAUSSIAN D_0=1.5 R_0=0.00001}
 ... DISTANCES
-\endverbatim
+\endplumedfile
 (see \ref DISTANCES)
 
 The advantage of MTRANSFORM_LESS comes, however, if you want to use transformed colvars as input
@@ -99,10 +99,10 @@ One is thus calculating the mean for those colvars that are less than the target
 
 The example shown below calculates the mean for those distances that less than 1.5 nm in length
 
-\verbatim
+\plumedfile
 DISTANCES GROUPA=1 GROUPB=2-50 MEAN LABEL=d1
 MFILTER_LESS DATA=d1 SWITCH={GAUSSIAN D_0=1.5 R_0=0.00001} MEAN LABEL=d4
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/FilterMoreThan.cpp b/src/multicolvar/FilterMoreThan.cpp
index 961c33d01627d1205d051c3e0430dd7032a2be33..50c1dc59b5228cfd0dccf01c038f1e07f7da0808 100644
--- a/src/multicolvar/FilterMoreThan.cpp
+++ b/src/multicolvar/FilterMoreThan.cpp
@@ -51,22 +51,22 @@ In other words, you are calculating the mean for the transformed colvar.
 The following input gives an example of how a MTRANSFORM_MORE action can be used to duplicate
 functionality that is elsehwere in PLUMED.
 
-\verbatim
+\plumedfile
 DISTANCES ...
  GROUPA=1-10 GROUPB=11-20
  LABEL=d1
 ... DISTANCES
 MTRANSFORM_MORE DATA=d1 SWITCH={GAUSSIAN D_0=1.5 R_0=0.00001}
-\endverbatim
+\endplumedfile
 
 In this case you can achieve the same result by using:
 
-\verbatim
+\plumedfile
 DISTANCES ...
  GROUPA=1-10 GROUPB=11-20
  MORE_THAN={GAUSSIAN D_0=1.5 R_0=0.00001}
 ... DISTANCES
-\endverbatim
+\endplumedfile
 (see \ref DISTANCES)
 
 The advantage of MTRANSFORM_MORE comes, however, if you want to use transformed colvars as input
@@ -104,10 +104,10 @@ One is thus calculating the mean for those colvars that are greater than the tar
 
 The example shown below calculates the mean for those distances that greater than 1.5 nm in length
 
-\verbatim
+\plumedfile
 DISTANCES GROUPA=1 GROUPB=2-50 MEAN LABEL=d1
 MFILTER_MORE DATA=d1 SWITCH={GAUSSIAN D_0=1.5 R_0=0.00001} MEAN LABEL=d4
-\endverbatim
+\endplumedfile
 
 More complicated things can be done by using the label of a filter as input to a new multicolvar as shown
 in the example below.  Here the coordination numbers of all atoms are computed.  The atoms with a coordination
@@ -115,11 +115,11 @@ number greater than 2 are then identified using the filter.  This reduced list o
 to a second coordination number calculation.  This second coordination number thus measures the number of
 two-coordinated atoms that each of the two-coordinated atoms is bound to.
 
-\verbatim
+\plumedfile
 1: COORDINATIONNUMBER SPECIES=1-150 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
 cf: MFILTER_MORE DATA=c1 SWITCH={RATIONAL D_0=2.0 R_0=0.1} LOWMEM
 c2: COORDINATIONNUMBER SPECIES=cf SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0} MORE_THAN={RATIONAL D_0=2.0 R_0=0.1}
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/InPlaneDistances.cpp b/src/multicolvar/InPlaneDistances.cpp
index d041fa1e0248477f15a1850e401b37115f116152..ef9f4a55f52d68932aec143af01146c508050e2f 100644
--- a/src/multicolvar/InPlaneDistances.cpp
+++ b/src/multicolvar/InPlaneDistances.cpp
@@ -53,10 +53,10 @@ Keywords such as MORE_THAN and LESS_THAN can then be used to calculate the numbe
 The following input can be used to calculate the number of atoms that have indices greater than 3 and less than 101 that
 are within a cylinder with a radius of 0.3 nm that has its long axis aligned with the vector connecting atoms 1 and 2.
 
-\verbatim
+\plumedfile
 d1: INPLANEDISTANCES VECTORSTART=1 VECTOREND=2 GROUP=3-100 LESS_THAN={RATIONAL D_0=0.2 R_0=0.1}
 PRINT ARG=d1.lessthan FILE=colvar
-\endverbatim
+\endplumedfile
 
 
 */
diff --git a/src/multicolvar/LocalAverage.cpp b/src/multicolvar/LocalAverage.cpp
index fdc3e09c4f71cefa57e26558666d4de08ce494d4..c8bae1d56ed6e70460c2fdb8762a505621646603 100644
--- a/src/multicolvar/LocalAverage.cpp
+++ b/src/multicolvar/LocalAverage.cpp
@@ -59,21 +59,21 @@ and so on.  You can also probe the value of these averaged variables in regions
 This example input calculates the coordination numbers for all the atoms in the system.  These coordination numbers are then averaged over
 spherical regions.  The number of averaged coordination numbers that are greater than 4 is then output to a file.
 
-\verbatim
+\plumedfile
 COORDINATIONNUMBER SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=d1
 LOCAL_AVERAGE ARG=d1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=4} LABEL=la
 PRINT ARG=la.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 This example input calculates the \f$q_4\f$ (see \ref Q4) vectors for each of the atoms in the system.  These vectors are then averaged
 component by component over a spherical region.  The average value for this quantity is then outputeed to a file.  This calculates the
 quantities that were used in the paper by Lechner and Dellago \cite dellago-q6
 
-\verbatim
+\plumedfile
 Q4 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=q4
 LOCAL_AVERAGE ARG=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=la
 PRINT ARG=la.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/MultiColvarDensity.cpp b/src/multicolvar/MultiColvarDensity.cpp
index 7ea3dc0096f3f8c2731cf885c94242efce8e3194..8e63b6e513a6043bc8e4edf8ee0a657a7142fbb5 100644
--- a/src/multicolvar/MultiColvarDensity.cpp
+++ b/src/multicolvar/MultiColvarDensity.cpp
@@ -58,11 +58,11 @@ The following example shows perhaps the simplest way in which this action can be
 input computes the density of atoms at each point on the grid and ouptuts this quantity to a file.  In
 other words this input instructs plumed to calculate \f$\rho(\mathbf{r}) = \sum_i K(\mathbf{r} - \mathbf{r}_i )\f$
 
-\verbatim
+\plumedfile
 dens: DENSITY SPECIES=1-100
 grid: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=xyz NBINS=100,100,100 BANDWIDTH=0.05,0.05,0.05 STRIDE=1
 DUMPGRID GRID=grid STRIDE=500 FILE=density
-\endverbatim
+\endplumedfile
 
 In the above example density is added to the grid on every step.  The PRINT_GRID instruction thus tells PLUMED to
 output the average density at each point on the grid every 500 steps of simulation.  Notice that the that grid output
@@ -72,11 +72,11 @@ of data separately you must use the CLEAR flag.
 This second example computes an order parameter (in this case \ref FCCUBIC) and constructs a phase field model
 for this order parameter using the equation above.
 
-\verbatim
+\plumedfile
 fcc: FCCUBIC SPECIES=1-5184 SWITCH={CUBIC D_0=1.2 D_MAX=1.5} ALPHA=27
 dens: MULTICOLVARDENS DATA=fcc ORIGIN=1 DIR=xyz NBINS=14,14,28 BANDWIDTH=1.0,1.0,1.0 STRIDE=1 CLEAR=1
 DUMPCUBE GRID=dens STRIDE=1 FILE=dens.cube
-\endverbatim
+\endplumedfile
 
 In this example the phase field model is computed and output to a file on every step of the simulation.  Furthermore,
 because the CLEAR=1 keyword is set on the MULTICOLVARDENS line each Gaussian cube file output is a phase field
diff --git a/src/multicolvar/NumberOfLinks.cpp b/src/multicolvar/NumberOfLinks.cpp
index b8257c5e4b822510b0c2b6df77f6162297efa430..c8850755bdd7c7e54a672b66485deb188960bd84 100644
--- a/src/multicolvar/NumberOfLinks.cpp
+++ b/src/multicolvar/NumberOfLinks.cpp
@@ -43,22 +43,22 @@ similar orientations.  The vectors on individual atoms could be Steinhardt param
 The following calculates how many bonds there are in a system containing 64 atoms and outputs
 this quantity to a file.
 
-\verbatim
+\plumedfile
 DENSITY SPECIES=1-64 LABEL=d1
 NLINKS ARG=d1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=dd
 PRINT ARG=dd FILE=colvar
-\endverbatim
+\endplumedfile
 
 The following calculates how many pairs of neighbouring atoms in a system containg 64 atoms have
 similar dispositions for the atoms in their coordination sphere.  This calculation uses the
 dot product of the Q6 vectors on adjacent atoms to measure whether or not two atoms have the same
 ``orientation"
 
-\verbatim
+\plumedfile
 Q6 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=q6
 NLINKS ARG=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=dd
 PRINT ARG=dd FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/Torsions.cpp b/src/multicolvar/Torsions.cpp
index 59d02cbfd72ae82b03f05604ab1b851fa281856d..4e3fcfbb2d897de06b98b95ab2f6c854502a2d75 100644
--- a/src/multicolvar/Torsions.cpp
+++ b/src/multicolvar/Torsions.cpp
@@ -40,7 +40,7 @@ Calculate whether or not a set of torsional angles are within a particular range
 
 The following provides an example of the input for the torsions command
 
-\verbatim
+\plumedfile
 TORSIONS ...
 ATOMS1=168,170,172,188
 ATOMS2=170,172,188,190
@@ -48,13 +48,13 @@ ATOMS3=188,190,192,230
 LABEL=ab
 ... TORSIONS
 PRINT ARG=ab.* FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 
 Writing out the atoms involved in all the torsions in this way can be rather tedious. Thankfully if you are working with protein you
 can avoid this by using the \ref MOLINFO command.  PLUMED uses the pdb file that you provide to this command to learn
 about the topology of the protein molecule.  This means that you can specify torsion angles using the following syntax:
 
-\verbatim
+\plumedfile
 MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
 TORSIONS ...
 ATOMS1=@phi-3
@@ -63,7 +63,7 @@ ATOMS3=@phi-4
 LABEL=ab
 ... TORSIONS
 PRINT ARG=ab FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 
 Here, \@phi-3 tells plumed that you would like to calculate the \f$\phi\f$ angle in the third residue of the protein.
 Similarly \@psi-4 tells plumed that you want to calculate the \f$\psi\f$ angle of the 4th residue of the protein.
diff --git a/src/multicolvar/VolumeAround.cpp b/src/multicolvar/VolumeAround.cpp
index bab79d8e0ab33d5893e78efe7fcf8202c19f8b09..d44f36d4b27857829b5a34877d8fa64a10c7e800 100644
--- a/src/multicolvar/VolumeAround.cpp
+++ b/src/multicolvar/VolumeAround.cpp
@@ -55,11 +55,11 @@ When AROUND is used with the \ref DENSITY action the number of atoms in the spec
 
 The following commands tell plumed to calculate the average coordination number for the atoms
 that have x (in fractional coordinates) within 2.0 nm of the com of mass c1. The final value will be labeled s.mean.
-\verbatim
+\plumedfile
 COM ATOMS=1-100 LABEL=c1
 COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 LABEL=c
 AROUND DATA=c ORIGIN=c1 XLOWER=-2.0 XUPPER=2.0 SIGMA=0.1 MEAN LABEL=s
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/VolumeCavity.cpp b/src/multicolvar/VolumeCavity.cpp
index af14a2be05491eae1e3d06f0ef3bff612286ef16..70d8254f323433ae025019d363ee46254bbabac5 100644
--- a/src/multicolvar/VolumeCavity.cpp
+++ b/src/multicolvar/VolumeCavity.cpp
@@ -86,19 +86,19 @@ described above and the resulting projections determine the \f$u'\f$, \f$v'\f$ a
 The following commands tell plumed to calculate the number of atoms in an ion chanel in a protein.
 The extent of the chanel is calculated from the positions of atoms 1, 4, 5 and 11. The final value will be labeled cav.
 
-\verbatim
+\plumedfile
 d1: DENSITY SPECIES=20-500
 CAVITY DATA=d1 ATOMS=1,4,5,11 SIGMA=0.1 LABEL=cav
-\endverbatim
+\endplumedfile
 
 The following command tells plumed to calculate the coordination numbers (with other water molecules) for the water
 molecules in the protein channel described above.  The average coordination number and the number of coordination
 numbers more than 4 is then calculated.  The values of these two quantities are given the labels cav.mean and cav.morethan
 
-\verbatim
+\plumedfile
 d1: COORDINATIONNUMBER SPECIES=20-500
 CAVITY DATA=d1 ATOMS=1,4,5,11 SIGMA=0.1 MEAN MORE_THAN={RATIONAL R_0=4} LABEL=cav
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/VolumeInCylinder.cpp b/src/multicolvar/VolumeInCylinder.cpp
index eb6bc611e4cb7af7b3b3a3f55ab9fb5cd6a8f407..3e8acfb226ac8fcd2392c4b28d576b4eb23c2e68 100644
--- a/src/multicolvar/VolumeInCylinder.cpp
+++ b/src/multicolvar/VolumeInCylinder.cpp
@@ -59,11 +59,11 @@ When INCYLINDER is used with the \ref DENSITY action the number of atoms in the
 The input below can be use to calculate the average coordination numbers for those atoms that are within a cylindrical tube
 of radius 1.5 nm that is centered on the position of atom 101 and that has its long axis parallel to the z-axis.
 
-\verbatim
+\plumedfile
 c1: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=0.1}
 d2: INCYLINDER ATOM=101 DATA=d1 DIRECTION=Z RADIUS={TANH R_0=1.5} SIGMA=0.1 LOWER=-0.1 UPPER=0.1 MEAN
 PRINT ARG=d2.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/VolumeInSphere.cpp b/src/multicolvar/VolumeInSphere.cpp
index 632af88df5e7d799fe0fd9735c8a3a939ade215b..937e067dd82d68fcaedcb9321c353b2bf0e7f3ec 100644
--- a/src/multicolvar/VolumeInSphere.cpp
+++ b/src/multicolvar/VolumeInSphere.cpp
@@ -58,11 +58,11 @@ When INCYLINDER is used with the \ref DENSITY action the number of atoms in the
 The input below can be use to calculate the average coordination numbers for those atoms that are within a sphere
 of radius 1.5 nm that is centered on the position of atom 101.
 
-\verbatim
+\plumedfile
 c1: COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=0.1}
 d2: INSPHERE ATOM=101 DATA=d1 RADIUS={TANH R_0=1.5} SIGMA=0.1 LOWER=-0.1 UPPER=0.1 MEAN
 PRINT ARG=d2.* FILE=colvar
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/VolumeTetrapore.cpp b/src/multicolvar/VolumeTetrapore.cpp
index 98fe915c33b7f55b5ee8b13ca983b2c12a97bcba..c7dcd66d04324e564adc15fd9a66a4faafb89787 100644
--- a/src/multicolvar/VolumeTetrapore.cpp
+++ b/src/multicolvar/VolumeTetrapore.cpp
@@ -94,19 +94,19 @@ This is in fact the only point of difference between these two actions.
 The following commands tell plumed to calculate the number of atom inside a tetrahedral cavity.  The extent of the tetrahedral
 cavity is calculated from the positions of atoms 1, 4, 5, and 11,  The final value will be labeled cav.
 
-\verbatim
+\plumedfile
 d1: DENSITY SPECIES=20-500
 TETRAHEDRALPORE DATA=d1 ATOMS=1,4,5,11 SIGMA=0.1 LABEL=cav
-\endverbatim
+\endplumedfile
 
 The following command tells plumed to calculate the coordination numbers (with other water molecules) for the water
 molecules in the tetrahedral cavity described above.  The average coordination number and the number of coordination
 numbers more than 4 is then calculated.  The values of these two quantities are given the labels cav.mean and cav.morethan
 
-\verbatim
+\plumedfile
 d1: COORDINATIONNUMBER SPECIES=20-500
 CAVITY DATA=d1 ATOMS=1,4,5,11 SIGMA=0.1 MEAN MORE_THAN={RATIONAL R_0=4} LABEL=cav
-\endverbatim
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/XAngle.cpp b/src/multicolvar/XAngle.cpp
index fbc2ab1d7104a9bb48620010dc2e633a5ce8d7e1..dd90c187f6f0fda16262a69a6df39e800316326c 100644
--- a/src/multicolvar/XAngle.cpp
+++ b/src/multicolvar/XAngle.cpp
@@ -41,10 +41,10 @@ Calculate the angles between the vector connecting two atoms and the x axis.
 
 The following input tells plumed to calculate the angles between the x-axis and the vector connecting atom 3 to atom 5 and between the x-axis
 and the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then
-\verbatim
+\plumedfile
 XANLGES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 */
 //+ENDPLUMEDOC
@@ -57,10 +57,10 @@ Calculate the angles between the vector connecting two atoms and the y axis.
 
 The following input tells plumed to calculate the angles between the y-axis and the vector connecting atom 3 to atom 5 and between the y-axis
 and the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then
-\verbatim
+\plumedfile
 YANLGES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 */
 //+ENDPLUMEDOC
@@ -73,10 +73,10 @@ Calculate the angles between the vector connecting two atoms and the z axis.
 
 The following input tells plumed to calculate the angles between the z-axis and the vector connecting atom 3 to atom 5 and between the z-axis
 and the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then
-\verbatim
+\plumedfile
 ZANLGES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 */
 //+ENDPLUMEDOC
diff --git a/src/multicolvar/XDistances.cpp b/src/multicolvar/XDistances.cpp
index 009b04de595c7f0a84f1566e07664769a34e1126..fe29a8b650baa8c14f47157ce4fa40f3e88191af 100644
--- a/src/multicolvar/XDistances.cpp
+++ b/src/multicolvar/XDistances.cpp
@@ -42,38 +42,38 @@ values such as the minimum, the number less than a certain quantity and so on.
 The following input tells plumed to calculate the x-component of the vector connecting atom 3 to atom 5 and
 the x-component of the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then
 printed
-\verbatim
+\plumedfile
 XDISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 
 
 The following input tells plumed to calculate the x-component of the vector connecting atom 3 to atom 5 and
 the x-component of the vector connecting atom 1 to atom 2.  The number of values that are
 less than 0.1nm is then printed to a file.
-\verbatim
+\plumedfile
 XDISTANCES ATOMS1=3,5 ATOMS2=1,2 LABEL=d1 LESS_THAN={RATIONAL R_0=0.1}
 PRINT ARG=d1.lt0.1
-\endverbatim
+\endplumedfile
 (See also \ref PRINT \ref switchingfunction).
 
 The following input tells plumed to calculate the x-components of all the distinct vectors that can be created
 between atoms 1, 2 and 3 (i.e. the vectors between atoms 1 and 2, atoms 1 and 3 and atoms 2 and 3).
 The average of these quantities is then calculated.
-\verbatim
+\plumedfile
 XDISTANCES GROUP=1-3 AVERAGE LABEL=d1
 PRINT ARG=d1.average
-\endverbatim
+\endplumedfile
 (See also \ref PRINT)
 
 The following input tells plumed to calculate all the vectors connecting the the atoms in GROUPA to the atoms in GROUPB.
 In other words the vector between atoms 1 and 2 and the vector between atoms 1 and 3.  The number of values
 more than 0.1 is then printed to a file.
-\verbatim
+\plumedfile
 XDISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1}
 PRINT ARG=d1.gt0.1
-\endverbatim
+\endplumedfile
 (See also \ref PRINT \ref switchingfunction)
 */
 //+ENDPLUMEDOC
@@ -89,38 +89,38 @@ values such as the minimum, the number less than a certain quantity and so on.
 The following input tells plumed to calculate the y-component of the vector connecting atom 3 to atom 5 and
 the y-component of the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then
 printed
-\verbatim
+\plumedfile
 YDISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 
 
 The following input tells plumed to calculate the y-component of the vector connecting atom 3 to atom 5 and
 the y-component of the vector connecting atom 1 to atom 2.  The number of values that are
 less than 0.1nm is then printed to a file.
-\verbatim
+\plumedfile
 YDISTANCES ATOMS1=3,5 ATOMS2=1,2 LABEL=d1 LESS_THAN={RATIONAL R_0=0.1}
 PRINT ARG=d1.lt0.1
-\endverbatim
+\endplumedfile
 (See also \ref PRINT \ref switchingfunction).
 
 The following input tells plumed to calculate the y-components of all the distinct vectors that can be created
 between atoms 1, 2 and 3 (i.e. the vectors between atoms 1 and 2, atoms 1 and 3 and atoms 2 and 3).
 The average of these quantities is then calculated.
-\verbatim
+\plumedfile
 YDISTANCES GROUP=1-3 AVERAGE LABEL=d1
 PRINT ARG=d1.average
-\endverbatim
+\endplumedfile
 (See also \ref PRINT)
 
 The following input tells plumed to calculate all the vectors connecting the the atoms in GROUPA to the atoms in GROUPB.
 In other words the vector between atoms 1 and 2 and the vector between atoms 1 and 3.  The number of values
 more than 0.1 is then printed to a file.
-\verbatim
+\plumedfile
 YDISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1}
 PRINT ARG=d1.gt0.1
-\endverbatim
+\endplumedfile
 (See also \ref PRINT \ref switchingfunction)
 
 */
@@ -137,38 +137,38 @@ values such as the minimum, the number less than a certain quantity and so on.
 The following input tells plumed to calculate the z-component of the vector connecting atom 3 to atom 5 and
 the z-component of the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then
 printed
-\verbatim
+\plumedfile
 ZDISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 
 
 The following input tells plumed to calculate the z-component of the vector connecting atom 3 to atom 5 and
 the z-component of the vector connecting atom 1 to atom 2.  The number of values that are
 less than 0.1nm is then printed to a file.
-\verbatim
+\plumedfile
 ZDISTANCES ATOMS1=3,5 ATOMS2=1,2 LABEL=d1 LESS_THAN={RATIONAL R_0=0.1}
 PRINT ARG=d1.lt0.1
-\endverbatim
+\endplumedfile
 (See also \ref PRINT \ref switchingfunction).
 
 The following input tells plumed to calculate the z-components of all the distinct vectors that can be created
 between atoms 1, 2 and 3 (i.e. the vectors between atoms 1 and 2, atoms 1 and 3 and atoms 2 and 3).
 The average of these quantities is then calculated.
-\verbatim
+\plumedfile
 ZDISTANCES GROUP=1-3 AVERAGE LABEL=d1
 PRINT ARG=d1.average
-\endverbatim
+\endplumedfile
 (See also \ref PRINT)
 
 The following input tells plumed to calculate all the vectors connecting the the atoms in GROUPA to the atoms in GROUPB.
 In other words the vector between atoms 1 and 2 and the vector between atoms 1 and 3.  The number of values
 more than 0.1 is then printed to a file.
-\verbatim
+\plumedfile
 ZDISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=0.1}
 PRINT ARG=d1.gt0.1
-\endverbatim
+\endplumedfile
 (See also \ref PRINT \ref switchingfunction)
 
 */
diff --git a/src/multicolvar/XYDistances.cpp b/src/multicolvar/XYDistances.cpp
index 767c64a86aafc268d23c4f0bce0c6f5bf80c1b4d..c0c8d4aa5272b49f7f1a05ab674fb1e3faa38f44 100644
--- a/src/multicolvar/XYDistances.cpp
+++ b/src/multicolvar/XYDistances.cpp
@@ -43,10 +43,10 @@ The following input tells plumed to calculate the projection of the length of th
 to atom 5 projected in the xy-plane and the projection of the length of the vector
 the vector connecting atom 1 to atom 2 in the xy-plane.  The minimum of these two quantities is then
 printed
-\verbatim
+\plumedfile
 XYDISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 
 */
@@ -64,10 +64,10 @@ The following input tells plumed to calculate the projection of the length of th
 to atom 5 projected in the xz-plane and the projection of the length of the vector
 the vector connecting atom 1 to atom 2 in the xz-plane.  The minimum of these two quantities is then
 printed
-\verbatim
+\plumedfile
 XZDISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 
 */
@@ -85,10 +85,10 @@ The following input tells plumed to calculate the projection of the length of th
 to atom 5 in the yz-plane and the projection of the length of the vector
 the vector connecting atom 1 to atom 2 in the yz-plane.  The minimum of these two quantities is then
 printed
-\verbatim
+\plumedfile
 YZDISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 
 */
diff --git a/src/multicolvar/XYTorsion.cpp b/src/multicolvar/XYTorsion.cpp
index fd512799a44f1e8121994fbb61b0ff8e3b708638..4e0804a508811722ae292890ea8c1b0ed1db84d4 100644
--- a/src/multicolvar/XYTorsion.cpp
+++ b/src/multicolvar/XYTorsion.cpp
@@ -41,10 +41,10 @@ Calculate the torsional angle around the x axis from the positive y direction.
 
 The following input tells plumed to calculate the angle around the x direction between the positive y-axis and the vector connecting atom 3 to atom 5 and
 the angle around the x direction between the positive y axis and the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then output
-\verbatim
+\plumedfile
 XYTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 */
 //+ENDPLUMEDOC
@@ -57,10 +57,10 @@ Calculate the torsional angle around the x axis from the positive z direction.
 
 The following input tells plumed to calculate the angle around the x direction between the positive z-axis and the vector connecting atom 3 to atom 5 and
 the angle around the x direction between the positive z direction and the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then output
-\verbatim
+\plumedfile
 XZTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 */
 //+ENDPLUMEDOC
@@ -73,10 +73,10 @@ Calculate the torsional angle around the y axis from the positive x direction.
 
 The following input tells plumed to calculate the angle around the y direction between the positive x-direction and the vector connecting atom 3 to atom 5 and
 the angle around the y direction between the positive x axis and the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then output
-\verbatim
+\plumedfile
 YXTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 */
 //+ENDPLUMEDOC
@@ -89,10 +89,10 @@ Calculate the torsional angle around the y axis from the positive z direction.
 
 The following input tells plumed to calculate the angle around the y direction between the positive z-direction and the vector connecting atom 3 to atom 5 and
 the angle around the y direction between the positive z direction and the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then output
-\verbatim
+\plumedfile
 YZTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 */
 //+ENDPLUMEDOC
@@ -105,10 +105,10 @@ Calculate the torsional angle around the z axis from the positive x direction.
 
 The following input tells plumed to calculate the angle around the z direction between the positive x-direction and the vector connecting atom 3 to atom 5 and
 the angle around the z direction between the positive x-direction and the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then output
-\verbatim
+\plumedfile
 ZXTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 */
 //+ENDPLUMEDOC
@@ -121,10 +121,10 @@ Calculate the torsional angle around the z axis from the positive y direction.
 
 The following input tells plumed to calculate the angle around the z direction between the positive y-axis and the vector connecting atom 3 to atom 5 and
 the angle around the z direction between the positive y axis and the vector connecting atom 1 to atom 2.  The minimum of these two quantities is then output
-\verbatim
+\plumedfile
 ZYTORSIONS ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
 PRINT ARG=d1.min
-\endverbatim
+\endplumedfile
 (See also \ref PRINT).
 */
 //+ENDPLUMEDOC
diff --git a/src/secondarystructure/AlphaRMSD.cpp b/src/secondarystructure/AlphaRMSD.cpp
index 6c9b49083ea11d2bc3f49863ab4984125333b8a6..ecf10eacb32de321aa4b63b509c050279713080b 100644
--- a/src/secondarystructure/AlphaRMSD.cpp
+++ b/src/secondarystructure/AlphaRMSD.cpp
@@ -68,11 +68,10 @@ anthing other than TYPE=DRMSD.  For more details as to how to do this see \ref W
 The following input calculates the number of six residue segments of
 protein that are in an alpha helical configuration.
 
-\verbatim
+\plumedfile
 MOLINFO STRUCTURE=helix.pdb
 ALPHARMSD RESIDUES=all TYPE=DRMSD LESS_THAN={RATIONAL R_0=0.08 NN=8 MM=12} LABEL=a
-\endverbatim
-(see also \ref MOLINFO)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/secondarystructure/AntibetaRMSD.cpp b/src/secondarystructure/AntibetaRMSD.cpp
index 6e60d54454e0d60860f6322d9e4e30e900938a7f..637e4c0b38522f58aec85f9860ac058964c0c6e9 100644
--- a/src/secondarystructure/AntibetaRMSD.cpp
+++ b/src/secondarystructure/AntibetaRMSD.cpp
@@ -70,11 +70,10 @@ anthing other than TYPE=DRMSD.  For more details as to how to do this see \ref W
 The following input calculates the number of six residue segments of
 protein that are in an antiparallel beta sheet configuration.
 
-\verbatim
+\plumedfile
 MOLINFO STRUCTURE=helix.pdb
 ANTIBETARMSD RESIDUES=all TYPE=DRMSD LESS_THAN={RATIONAL R_0=0.08 NN=8 MM=12} LABEL=a
-\endverbatim
-(see also \ref MOLINFO)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/secondarystructure/ParabetaRMSD.cpp b/src/secondarystructure/ParabetaRMSD.cpp
index fe87408643184d808adea19581f33671f7b78435..1a67ec8afcf6de744c7b99882055c7731615ddea 100644
--- a/src/secondarystructure/ParabetaRMSD.cpp
+++ b/src/secondarystructure/ParabetaRMSD.cpp
@@ -70,11 +70,10 @@ anthing other than TYPE=DRMSD.  For more details as to how to do this see \ref W
 The following input calculates the number of six residue segments of
 protein that are in an parallel beta sheet configuration.
 
-\verbatim
+\plumedfile
 MOLINFO STRUCTURE=helix.pdb
 PARABETARMSD RESIDUES=all TYPE=DRMSD LESS_THAN={RATIONAL R_0=0.08 NN=8 MM=12} LABEL=a
-\endverbatim
-(see also \ref MOLINFO)
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/setup/Load.cpp b/src/setup/Load.cpp
index 93f9d39b8831e7e18e926b3b61b0e01e66ce6b62..4002ed0f0ff6671d04b5f9d9aa7ef1068e13e027 100644
--- a/src/setup/Load.cpp
+++ b/src/setup/Load.cpp
@@ -43,9 +43,9 @@ If you have a shared object named extensions.so and want to
 use the functionalities implemented in it within PLUMED you can
 load it with the following syntax
 
-\verbatim
+\plumedfile
 LOAD FILE=extensions.so
-\endverbatim
+\endplumedfile
 
 As a more practical example, imagine that you want to make a
 small change to one collective variable that is already implemented
@@ -61,7 +61,7 @@ with different names. Then you can compile it into a shared object using
 This will generate a file `Distance2.so` (or `Distance2.dylib` on a mac)
 that can be loaded.
 Now you can use your new implementation with the following input
-\verbatim
+\plumedfile
 # load the new library
 LOAD FILE=Distance2.so
 # compute standard distance
@@ -70,11 +70,11 @@ d: DISTANCE ATOMS=1,10
 d2: DISTANCE2 ATOMS=1,10
 # print them on a file
 PRINT ARG=d,d2 FILE=compare-them
-\endverbatim
+\endplumedfile
 
 You can even skip the initial step and directly feed PLUMED
 with the `Distance2.cpp` file: it will be compiled on the fly.
-\verbatim
+\plumedfile
 # load the new definition
 # this is a cpp file so it will be compiled
 LOAD FILE=Distance2.cpp
@@ -84,7 +84,7 @@ d: DISTANCE ATOMS=1,10
 d2: DISTANCE2 ATOMS=1,10
 # print them on a file
 PRINT ARG=d,d2 FILE=compare-them
-\endverbatim
+\endplumedfile
 
 This will allow to make quick tests while developing your own
 variables. Of course, after your implementation is ready you might
diff --git a/src/setup/MolInfo.cpp b/src/setup/MolInfo.cpp
index 03a0012f97a8eaac8c5561aa7fa55f303f718282..24f2cb68922e8b89c84d8c28dca4784a3bfef40d 100644
--- a/src/setup/MolInfo.cpp
+++ b/src/setup/MolInfo.cpp
@@ -111,23 +111,21 @@ interpret terminal residue 1.
 In the following example the MOLINFO command is used to provide the information on which atoms
 are in the backbone of a protein to the ALPHARMSD CV.
 
-\verbatim
+\plumedfile
 MOLINFO STRUCTURE=reference.pdb
 ALPHARMSD RESIDUES=all TYPE=DRMSD LESS_THAN={RATIONAL R_0=0.08 NN=8 MM=12} LABEL=a
-\endverbatim
-(see also \ref ALPHARMSD)
+\endplumedfile
 
 The following example prints the distance corresponding to the hydrogen bonds
 in a GC Watson-Crick pair.
 
-\verbatim
+\plumedfile
 MOLINFO STRUCTURE=reference.pdb
 hb1: DISTANCE ATOMS=@N2-1,@O2-14
 hb2: DISTANCE ATOMS=@N1-1,@N3-14
 hb3: DISTANCE ATOMS=@O6-1,@N4-14
 PRINT ARG=hb1,hb2,hb3
-\endverbatim
-(see also \ref DISTANCE).
+\endplumedfile
 
 
 */
diff --git a/src/setup/Restart.cpp b/src/setup/Restart.cpp
index ffa4f46b1b4487f628f4dea4d9965778378369c1..4cfb03ef0364cd97e66388dd0ca503e05aaa03e9 100644
--- a/src/setup/Restart.cpp
+++ b/src/setup/Restart.cpp
@@ -51,31 +51,28 @@ and \ref PBMETAD and on some analysis action.
 \par Examples
 
 Using the following input:
-\verbatim
+\plumedfile
 d: DISTANCE ATOMS=1,2
 PRINT ARG=d FILE=out
-\endverbatim
-(See also \ref DISTANCE and \ref PRINT).
+\endplumedfile
 a new 'out' file will be created. If an old one is on the way, it will be automatically backed up.
 
 On the other hand, using the following input:
-\verbatim
+\plumedfile
 RESTART
 d: DISTANCE ATOMS=1,2
 PRINT ARG=d FILE=out
-\endverbatim
+\endplumedfile
 the file 'out' will be appended.
-(See also \ref DISTANCE and \ref PRINT).
 
 In the following case, file out1 will be backed up and file out2 will be concatenated
-\verbatim
+\plumedfile
 RESTART
 d1: DISTANCE ATOMS=1,2
 d2: DISTANCE ATOMS=1,2
 PRINT ARG=d1 FILE=out1 RESTART=NO
 PRINT ARG=d2 FILE=out2
-\endverbatim
-(See also \ref DISTANCE and \ref PRINT).
+\endplumedfile
 
 
 
diff --git a/src/setup/Units.cpp b/src/setup/Units.cpp
index 1b4225d24482ef09d48d08da5561b24f48b8e1cb..617e8c0a6f1a89bca9a736ef9c5833040438835f 100644
--- a/src/setup/Units.cpp
+++ b/src/setup/Units.cpp
@@ -44,16 +44,15 @@ the units. For example, trajectories written in .gro format (with \ref DUMPATOMS
 are going to be always in nm.
 
 \par Examples
-\verbatim
+\plumedfile
 # this is using nm - kj/mol - fs
 UNITS LENGTH=nm TIME=fs
-\endverbatim
+\endplumedfile
 If a number, x, is found, the new unit is equal to x (default units)
-\verbatim
+\plumedfile
 # this is using nm - kj/mol - fs
 UNITS LENGTH=nm TIME=0.001
-\endverbatim
-
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/vatom/COM.cpp b/src/vatom/COM.cpp
index 8eb60b64d6c6ecab0086c20f2a3d977e15a3fa41..d0cfd0b1feffc1f163e57309aca191192d738157 100644
--- a/src/vatom/COM.cpp
+++ b/src/vatom/COM.cpp
@@ -54,13 +54,12 @@ periodic image.
 
 The following input instructs plumed to print the distance between the
 center of mass for atoms 1,2,3,4,5,6,7 and that for atoms 15,20:
-\verbatim
+\plumedfile
 c1: COM ATOMS=1-7
 c2: COM ATOMS=15,20
 d1: DISTANCE ATOMS=c1,c2
 PRINT ARG=d1
-\endverbatim
-(See also \ref DISTANCE and \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/vatom/Center.cpp b/src/vatom/Center.cpp
index 7e4cdd440239f6dbe9ef99ced9da6d836017dacc..64eea7d7419381476acfa879d1640cd38fdc8d79 100644
--- a/src/vatom/Center.cpp
+++ b/src/vatom/Center.cpp
@@ -54,7 +54,7 @@ periodic image.
 
 \par Examples
 
-\verbatim
+\plumedfile
 # a point which is on the line connecting atoms 1 and 10, so that its distance
 # from 10 is twice its distance from 1:
 c1: CENTER ATOMS=1,1,10
@@ -67,8 +67,7 @@ c2: CENTER ATOMS=2,3,4,5 MASS
 d1: DISTANCE ATOMS=c1,c2
 
 PRINT ARG=d1
-\endverbatim
-(See also \ref DISTANCE, \ref COM and \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/src/vatom/FixedAtom.cpp b/src/vatom/FixedAtom.cpp
index ccf7eb52ca84063a68e269607643d8d4bc81c1b4..3baf40379a1c7c7f271c33e5563d16d7878ff8c6 100644
--- a/src/vatom/FixedAtom.cpp
+++ b/src/vatom/FixedAtom.cpp
@@ -51,24 +51,22 @@ then it is safe to add further fixed atoms without breaking translational invari
 
 The following input instructs plumed to compute the angle between
 distance of atoms 15 and 20 and the z axis and keeping it close to zero.
-\verbatim
+\plumedfile
 a: FIXEDATOM AT=0,0,0
 b: FIXEDATOM AT=0,0,1
 an: ANGLE ATOMS=a,b,15,20
 RESTRAINT ARG=an AT=0.0 KAPPA=100.0
-\endverbatim
-(See also \ref ANGLE and \ref RESTRAINT).
+\endplumedfile
 
 The following input instructs plumed to align a protein on a template
 and then compute the distance of one of its atom from the point
 (10,20,30).
-\verbatim
+\plumedfile
 FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=SIMPLE
 a: FIXEDATOM AT=10,20,30
 d: DISTANCE ARG=a,20
 PRINT ARG=d FILE=colvar
-\endverbatim
-(See also \ref FIT_TO_TEMPLATE and \ref DISTANCE).
+\endplumedfile
 
 
 */
diff --git a/src/vatom/Ghost.cpp b/src/vatom/Ghost.cpp
index b2c7837172b17f90901431676c4ddaf3c1c37f6e..7c57cee4ffad29f228ab766fd241a685065e0134 100644
--- a/src/vatom/Ghost.cpp
+++ b/src/vatom/Ghost.cpp
@@ -39,13 +39,12 @@ an atom list through the the label for the GHOST action that creates it.
 \par Examples
 The following input instructs plumed to print the distance between the
 ghost atom and the center of mass for atoms 15,20:
-\verbatim
+\plumedfile
 c1: GHOST ATOMS=1,5,10 COORDINATES=10.0,10.0,10.0
 c2: COM ATOMS=15,20
 d1: DISTANCE ATOMS=c1,c2
 PRINT ARG=d1
-\endverbatim
-(See also \ref DISTANCE and \ref PRINT).
+\endplumedfile
 
 */
 //+ENDPLUMEDOC
diff --git a/user-doc/Analysis.txt b/user-doc/Analysis.txt
index a5a0e55729a80ac63bba00f7ee2844e9b33694ed..ef3a005893acfb9f66e3f5bcfc294c8c73f96456 100644
--- a/user-doc/Analysis.txt
+++ b/user-doc/Analysis.txt
@@ -45,11 +45,11 @@ As an example the following set of commands instructs PLUMED to calculate the di
 atoms 1 and 2 for every 5th frame in the trajectory and to accumulate a histogram from this data
 which will be output every 100 steps (i.e. when 20 distances have been added to the histogram).
 
-\verbatim
+\plumedfile
 x: DISTANCE ATOMS=1,2
 h: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1 STRIDE=5
 DUMPGRID GRID=h FILE=histo STRIDE=100 
-\endverbatim
+\endplumedfile
 
 It is important to note when using commands such as the above the first frame in the trajectory is assumed 
 to be the initial configuration that was input to the MD code. It is thus ignored.  Furthermore, if you are 
@@ -57,11 +57,11 @@ running with driver and you would like to analyse the whole trajectory (without
 and then print the result you simply call \ref DUMPGRID (or any of the commands above) without a STRIDE 
 keyword as shown in the example below. 
 
-\verbatim
+\plumedfile
 x: DISTANCE ATOMS=1,2
 h: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1 STRIDE=5
 DUMPGRID GRID=h FILE=histo 
-\endverbatim
+\endplumedfile
 
 Please note that even with this calculation the first frame in the trajectory is ignored when computing the 
 histogram.
diff --git a/user-doc/Colvar.txt b/user-doc/Colvar.txt
index 7c6557796d1e5f042bd893987976443b627b175e..e3e20d2a1b87511159eede1ccbb0a094db31c703 100644
--- a/user-doc/Colvar.txt
+++ b/user-doc/Colvar.txt
@@ -92,9 +92,9 @@ Oftentimes the simplest way to specify the atoms involved is to use multiple ins
 i.e. ATOMS1, ATOMS2, ATOMS3,...  Separate instances of the quantity specified by NAME are then calculated for 
 each of the sets of atoms.  For example if the command issued contains the following:
 
-\verbatim
+\plumedfile
 DISTANCES ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
-\endverbatim
+\endplumedfile
 
 The distances between atoms 1 and 2, atoms 3 and 4, and atoms 5 and 6 are calculated. Obviously, generating 
 this sort of input is rather tedious so short cuts are also available many of the collective variables. 
@@ -163,11 +163,11 @@ prevent the cluster subliming.  Alternatively, you may wish to insist that a par
 coordination number greater than 2.  You can add these sorts of restraints by employing the following biases, which all act 
 on the set of collective variable values calculated by a multicolvar.  So for example the following set of commands:
 
-\verbatim
+\plumedfile
 COM ATOMS=1-20 LABEL=c1
 DISTANCES GROUPA=c1 GROUPB=1-20 LABEL=d1
 UWALLS DATA=d1 AT=2.5 KAPPA=0.2 LABEL=sr
-\endverbatim
+\endplumedfile
 
 creates the aforementioned set of restraints on the distances between the 20 atoms in a cluster and the center of mass of the cluster.
 
diff --git a/user-doc/Files.txt b/user-doc/Files.txt
index f360af5d28e4a86540994c23ef7ce69bd9e6d5ca..2ed8f21ee892fb39d4d1419bff7ddbdf3e3c0ee2 100644
--- a/user-doc/Files.txt
+++ b/user-doc/Files.txt
@@ -34,11 +34,10 @@ since your disk might fill up quickly with this setting.
 When running with multiple replicas (e.g., with GROMACS, -multi option) PLUMED adds the replica index as a suffix to
 all the files. The following command
 will thus print files named COLVAR.0, COLVAR.1, etc for the different replicas.
-\verbatim
+\plumedfile
 d: DISTANCE ATOMS=1,2
 PRINT ARG=d FILE=COLVAR
-\endverbatim
-(see also \ref DISTANCE and \ref PRINT).
+\endplumedfile
 
 When reading a file, PLUMED will try to add the suffix. If the file is not found, it will fall back to
 the name without suffix. The most important case is the reading of the plumed input file.
@@ -51,12 +50,11 @@ extension. Before PLUMED 2.2, the only recognized suffix was ".gz". Since 2.2, a
 less or equal to five letters is recognized.
 
 This means that using in a multireplica context an input such as
-\verbatim
+\plumedfile
 d: DISTANCE ATOMS=1,2
 PRINT ARG=d FILE=COLVAR.gz
 METAD ARG=d FILE=test.HILLS SIGMA=0.1 HEIGHT=0.1
-\endverbatim
-(see \ref DISTANCE, \ref PRINT, and \ref METAD)
+\endplumedfile
 PLUMED will write files named COLVAR.0.gz, COLVAR.1.gz, test.0.HILLS, test.1.HILLS, etc
 etc. This is useful since the preserved extension makes it easy
 to process the files later.
diff --git a/user-doc/Functions.txt b/user-doc/Functions.txt
index 8bf37c42db0e66f5e2727bb98b419edeb4f88217..53de40f365fef8acff528a4a3d0596dd259e7db4 100644
--- a/user-doc/Functions.txt
+++ b/user-doc/Functions.txt
@@ -23,7 +23,7 @@ is the smallest value and `B` is the largest value. In case the answer is that t
 values at the discontinuity are not always the same, then you cannot construct a variable that
 can be biased with PLUMED. Consider the following examples:
 
-\verbatim
+\plumedfile
 t: TORSION ATOMS=1,2,3,4
 # When atoms are moved, t could jump suddenly from -pi to +pi
 
@@ -39,7 +39,7 @@ d: DISTANCE ATOMS=1,10 COMPONENTS
 # make a new variable equal to d.z but with the correct periodicity
 dz: COMBINE ARG=d.z PERIODIC=-10,10
 # here we assumed the system is in a orthorhombic box with z side = 20
-\endverbatim
+\endplumedfile
 
 @FUNCTION@
 
diff --git a/user-doc/Group.txt b/user-doc/Group.txt
index 6365d9f7dd91b987cb4588a4fb6ab2a1b80610df..6edc0225b513c21370ca85c0f3712589ab624ce1 100644
--- a/user-doc/Group.txt
+++ b/user-doc/Group.txt
@@ -56,12 +56,12 @@ See http://en.wikipedia.org/wiki/Ramachandran_plot
 The following example shows how to use \ref MOLINFO with \ref TORSION to calculate the torsion angles phi and psi for the first and fourth residue
 of the protein:
 
-\verbatim
+\plumedfile
 MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
 t1: TORSION ATOMS=@phi-3
 t2: TORSION ATOMS=@psi-4
 PRINT ARG=t1,t2 FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 
 \subsection pbc Broken Molecules and PBC 
 
@@ -76,11 +76,11 @@ using the ALIGN_ATOMS keyword.  In PLUMED 2 the same effect can be achieved usin
 
 The following input computes the end-to-end distance for a polymer of 100 atoms and keeps it at a value around 5.
 
-\verbatim
+\plumedfile
 WHOLEMOLECULES ENTITY0=1-100
 e2e: DISTANCE ATOMS=1,100 NOPBC
 RESTRAINT ARG=e2e KAPPA=1 AT=5
-\endverbatim
+\endplumedfile
 
 Notice that NOPBC is used to be sure in \ref DISTANCE that if the end-to-end distance is larger than half the simulation box the distance 
 is compute properly. Also notice that, since many MD codes break molecules across cell boundary, it might be necessary to use the 
@@ -90,12 +90,12 @@ Notice that most expressions are invariant with respect to a change in the order
 but some of them depend on that order. E.g., with \ref WHOLEMOLECULES it could be useful to
 specify atom lists in a reversed order.
 
-\verbatim
+\plumedfile
 # to see the effect, one could dump the atoms as they were before molecule reconstruction:
 # DUMPATOMS FILE=dump-broken.xyz ATOMS=1-20
 WHOLEMOLECULES STRIDE=1 ENTITY0=1-20
 DUMPATOMS FILE=dump.xyz ATOMS=1-20
-\endverbatim
+\endplumedfile
 
 Notice that there are other ways to manipulate the coordinates stored within PLUMED:
 - Using the \subpage FIT_TO_TEMPLATE they can be aligned to a template structure.
@@ -116,10 +116,10 @@ To specify to a colvar that you want to use the position of a virtual atom to ca
 in your system you simply use the label for your virtual atom in place of the usual numerical index.  Virtual
 atoms and normal atoms can be mixed together in the input to colvars as shown below:
 
-\verbatim
+\plumedfile
 COM ATOMS=1,10 LABEL=com1
 DISTANCE ATOMS=11,com1
-\endverbatim 
+\endplumedfile
 
 If you don't want to calculate CVs from the virtual atom.  That is to say you just want to monitor the position of a virtual atom 
 (or any set of atoms) over the course of your trajectory you can do this using \ref DUMPATOMS.
diff --git a/user-doc/Misc.txt b/user-doc/Misc.txt
index 284a2a6e9bdfc65ce0b7e12f698bafd089cea19d..4fceb133df2803ed71979ab1277ea6c00c770dc7 100644
--- a/user-doc/Misc.txt
+++ b/user-doc/Misc.txt
@@ -19,12 +19,12 @@ particular simulation you might find it useful to put comments in your input fil
 comments can be added using a # sign.  On any given line everything after the # sign is ignored so 
 erm... yes add lines of comments or trailing comments to your hearts content as shown below (using Shakespeare is optional):
 
-\verbatim
+\plumedfile
 # This is the distance between two atoms:
 DISTANCE ATOM=1,2 LABEL=d1
 UPPER_WALLS ARG=d1 AT=3.0 KAPPA=3.0 LABEL=Snout # In this same interlude it doth befall.
 # That I, one Snout by name, present a wall.
-\endverbatim
+\endplumedfile
 (see \ref DISTANCE and \ref UPPER_WALLS)
 
 An alternative to including comments in this way is to use line starting ENDPLUMED.  Everything in the PLUMED input after this
@@ -37,14 +37,14 @@ We at PLUMED are aware of this fact and thus have provided a way of doing line c
 easier - aren't we kind?  Well no not really, we have to use this code too.  Anyway, you can do continuations by using the "..." syntax
 as this makes this: 
 
-\verbatim
+\plumedfile
 DISTANCES ATOMS1=1,300 ATOMS2=1,400 ATOMS3=1,500 LABEL=dist
-\endverbatim
+\endplumedfile
 (see \ref DISTANCES)
 
 equivalent to this:
 
-\verbatim
+\plumedfile
 DISTANCES ...
   LABEL=dist
 # we can also insert comments here
@@ -54,11 +54,11 @@ DISTANCES ...
 #empty lines are also allowed
 
 ... DISTANCES
-\endverbatim
+\endplumedfile
 
 Notice that the closing `...` is followed by the word `DISTANCES`. This is optional, but might be
 useful to find more easily which is the matching start of the statement. The following is equally correct
-\verbatim
+\plumedfile
 DISTANCES ...
   LABEL=dist
 # we can also insert comments here
@@ -68,12 +68,12 @@ DISTANCES ...
 #empty lines are also allowed
 
 ...
-\endverbatim
+\endplumedfile
 
 Notice that PLUMED makes a check that the word following the closing `...` is actually identical to
 the first word in the line with the first `...`. If not, it will throw an error.
 Also notice that you might put more than one word in the first line. E.g.
-\verbatim
+\plumedfile
 DISTANCES LABEL=dist ...
 # we can also insert comments here
   ATOMS1=1,300
@@ -81,9 +81,9 @@ DISTANCES LABEL=dist ...
   ATOMS2=1,400 ATOMS3=1,500
 #empty lines are also allowed
 ...
-\endverbatim
+\endplumedfile
 or, equivalently,
-\verbatim
+\plumedfile
 dist: DISTANCES ...
 # we can also insert comments here
   ATOMS1=1,300
@@ -91,7 +91,7 @@ dist: DISTANCES ...
   ATOMS2=1,400 ATOMS3=1,500
 #empty lines are also allowed
 ...  
-\endverbatim
+\endplumedfile
 
 
 \page VimSyntax Using VIM syntax file
@@ -171,11 +171,11 @@ Add to your `.vimrc` file the following commands
 :set modelines=5
 \endverbatim
 Then, at the beginning of your PLUMED input file, put the following comment:
-\verbatim
+\plumedfile
 # vim:ft=plumed
 d: DISTANCE ATOMS=1,2
 RESTRAINT ARG=d AT=0.0 KAPPA=1.0
-\endverbatim
+\endplumedfile
 Now, every time you open this file, you will see it highlighted.
 
 \par Syntax highlighting
@@ -220,27 +220,27 @@ Folded lines can be expanded with `zo` and folded with `zc`. Look at VIM documen
 learn more.
 In case you want to use this feature, we suggest you to put both label
 and action type on the first line of multi-line statements. E.g.
-\verbatim
+\plumedfile
 m: METAD ...
   ARG=d
   HEIGHT=1.0
   SIGMA=0.5
   PACE=100
 ...
-\endverbatim
+\endplumedfile
 will be folded to
 \verbatim
 +--  6 lines: m: METAD ...------------------------------------------------------
 \endverbatim
 and
-\verbatim
+\plumedfile
 METAD LABEL=m ...
   ARG=d
   HEIGHT=1.0
   SIGMA=0.5
   PACE=100
 ...
-\endverbatim
+\endplumedfile
 will be folded to
 \verbatim
 +--  6 lines: METAD LABEL=m ...-------------------------------------------------
@@ -293,9 +293,9 @@ Try to do the following. Enable plumed syntax:
 :set ft=plumed
 \endverbatim
 Then add the following line
-\verbatim
+\plumedfile
 DISTANCE
-\endverbatim
+\endplumedfile
 Now, in normal mode, go with the cursor on the `DISTANCE` line and type
 \verbatim
 :PHelp
@@ -365,38 +365,37 @@ vimrc file.
 
 If, for some reason, you want to spread your PLUMED input over a number of files you can use \subpage INCLUDE as shown below:
 
-\verbatim
+\plumedfile
 INCLUDE FILE=filename
-\endverbatim
+\endplumedfile
 
 So, for example, a single "plumed.dat" file:
 
-\verbatim
+\plumedfile
 DISTANCE ATOMS=0,1 LABEL=dist
 RESTRAINT ARG=dist
-\endverbatim
-(see \ref DISTANCE and \ref RESTRAINT)
+\endplumedfile
 
 could be split up into two files as shown below:
  
-\verbatim
+\plumedfile
 DISTANCE ATOMS=0,1 LABEL=dist
 INCLUDE FILE=toBeIncluded.dat
-\endverbatim
+\endplumedfile
 plus a "toBeIncluded.dat" file
-\verbatim
+\plumedfile
 RESTRAINT ARG=dist
-\endverbatim
+\endplumedfile
 
 However, when you do this it is important to recognise that \ref INCLUDE is a real directive that is only resolved
 after all the \ref comments have been stripped and the \ref ContinuationLines have been unrolled.  This means it
 is not possible to do things like:
 
-\verbatim
+\plumedfile
 # this is wrong:
 DISTANCE INCLUDE FILE=options.dat
 RESTRAINT ARG=dist
-\endverbatim
+\endplumedfile
 
 \page load Loading shared libraries
 
@@ -405,9 +404,9 @@ PLUMED libraries.  Alternatively, if you want to keep your code independent from
 so you can release it independely - we won't be offended), then you can create your own dynamic library.  To use this 
 in conjuction with PLUMED you can then load it at runtime by using the \subpage LOAD keyword as shown below:
 
-\verbatim
+\plumedfile
 LOAD FILE=library.so
-\endverbatim
+\endplumedfile
  
 N.B.  If your system uses a different suffix for dynamic libraries (e.g. macs use .dylib) then PLUMED will try to 
 automatically adjust the suffix accordingly.
diff --git a/user-doc/Performances.txt b/user-doc/Performances.txt
index 57bc2c66a554883c04f067b5a2173e32045378a6..00fc973acca246bc871164c6f777ee7e8008cc6e 100644
--- a/user-doc/Performances.txt
+++ b/user-doc/Performances.txt
@@ -105,12 +105,12 @@ a factor equal to STRIDE.
 This technique can allow your simulation to run faster if you need
 the apply a bias potential on some very expensive collective variable.
 Consider the following input:
-\verbatim
+\plumedfile
 c1: COM ATOMS=1-1000
 c2: COM ATOMS=1001-2000
 d:  DISTANCE ATOMS=c1,c2
 METAD ARG=d HEIGHT=1 SIGMA=0.1 BIASFACTOR=5 PACE=500
-\endverbatim
+\endplumedfile
 This performs a \ref METAD simulation biasing the distance between two
 centers of mass. Since computing these centers requires a lot of atoms
 to be imported from the MD engine, it could slow down significantly the
@@ -118,12 +118,12 @@ simulation. Notice that whereas the bias is changed every PACE=500 steps,
 it is applied every STRIDE step, where STRIDE=1 by default.
 The following input could lead to a significantly faster simulation at the price
 of a negligible systematic error
-\verbatim
+\plumedfile
 c1: COM ATOMS=1-1000
 c2: COM ATOMS=1001-2000
 d:  DISTANCE ATOMS=c1,c2
 METAD ARG=d HEIGHT=1 SIGMA=0.1 BIASFACTOR=5 PACE=500 STRIDE=2
-\endverbatim
+\endplumedfile
 Similarly, the STRIDE keyword can be used with other biases (e.g. \ref RESTRAINT).
 
 The technique is discussed in details here \cite Ferrarotti2015.
@@ -133,9 +133,9 @@ See also \subpage EFFECTIVE_ENERGY_DRIFT.
 
 Whenever you have a multicolvar action such as:
 
-\verbatim
+\plumedfile
 COORDINATIONNUMBER SPECIES=1-100 SWITCH={RATIONAL R_0=1. D_MAX=3.0} MORE_THAN={RATIONAL R_0=6.0 NN=6 MM=12 D_0=0}
-\endverbatim
+\endplumedfile
 
 You will get a collosal speedup by specifying the D_MAX keyword in all switching functions that act on distances.
 D_MAX tells PLUMED that the switching function is strictly zero if the distance is greater than this value.  As a result
diff --git a/user-doc/Regex.txt b/user-doc/Regex.txt
index b24524a0f8424daa1a03ab3853efe475f18549fa..93561f632200703b840f5a97f9f5530691341fdb 100644
--- a/user-doc/Regex.txt
+++ b/user-doc/Regex.txt
@@ -14,14 +14,14 @@ Regular expressions are enclosed in round braces and must not contain spaces (th
 names have no spaces indeed, so why use them?).
 
 As an example then command
-\verbatim
+\plumedfile
 d1: DISTANCE ATOMS=1,2 COMPONENTS
 PRINT ARG=(d1\.[xy])   STRIDE=100 FILE=colvar FMT=%8.4f
-\endverbatim
+\endplumedfile
 will cause both the d1.x and d1.y components of the DISTANCE action to be printed out in the order that they are created by plumed.
 The "." character must be escaped in order to interpret it as a literal ".". An unescaped dot is a wildcard which is matched by any character,
 So as an example
-\verbatim
+\plumedfile
 d1: DISTANCE ATOMS=1,2 COMPONENTS
 dxy: DISTANCE ATOMS=1,3
 
@@ -30,26 +30,26 @@ PRINT ARG=(d1.[xy])   STRIDE=100 FILE=colvar FMT=%8.4f
 
 # while this will match d1.x,d1.y only
 PRINT ARG=(d1\.[xy])   STRIDE=100 FILE=colvar FMT=%8.4f
-\endverbatim
+\endplumedfile
 
 You can include more than one regular expression by using comma separated regular expressions 
 
-\verbatim
+\plumedfile
 t1: TORSION ATOMS=5,7,9,15
 t2: TORSION ATOMS=7,9,15,17
 d1: DISTANCE ATOMS=7,17 COMPONENTS
 PRINT ARG=(d1\.[xy]),(t[0-9]) STRIDE=100 FILE=colvar FMT=%8.4f
-\endverbatim
+\endplumedfile
 
 (this selects t1,t2,d1.x and d2.x) Be aware that if you have overlapping selection they will be duplicated so it 
 a better alternative is to use the "or" operator "|". 
 
-\verbatim
+\plumedfile
 t1: TORSION ATOMS=5,7,9,15
 t2: TORSION ATOMS=7,9,15,17
 d1: DISTANCE ATOMS=7,17 COMPONENTS
 PRINT ARG=(d1\.[xy]|t[0-9]) STRIDE=100 FILE=colvar FMT=%8.4f
-\endverbatim
+\endplumedfile
 
 this selects the same set of arguments as the previous example.
 
diff --git a/user-doc/Syntax.txt b/user-doc/Syntax.txt
index 47c73c1f2f6066f5bc66a2e5e6a282e57baee995..277d339e9e98234513d0177008319d53701c42a8 100644
--- a/user-doc/Syntax.txt
+++ b/user-doc/Syntax.txt
@@ -25,18 +25,16 @@ is enclosed in curly braces (e.g. ATOMS={1 2 3 4}).  Please note that you can sp
 The most important of these keywords is the label keyword as it is only by using these labels that we can pass data 
 from one action to another.  As an example if you do:
 
-\verbatim
+\plumedfile
 DISTANCE ATOMS=1,2
-\endverbatim
-(see \ref DISTANCE)
+\endplumedfile
 
 Then PLUMED will do nothing other than read in your input file.  In contrast if you do:
 
-\verbatim
+\plumedfile
 DISTANCE ATOMS=1,2 LABEL=d1
 PRINT ARG=d1 FILE=colvar STRIDE=10
-\endverbatim
-(see \ref PRINT)
+\endplumedfile
 
 then PLUMED will print out the value of the distance between atoms 1 and 2 every 10 steps to the file colvar as you have told
 PLUMED to take the value calculated by the action d1 and to print it. You can use any character string to label your actions
@@ -46,10 +44,10 @@ code-generated groups of atoms and to give labels to any Actions for which the u
 Notice that if a word followed by a column is added at the beginning of the line (e.g. pippo:), PLUMED automatically
 removes it and adds an equivalent label (LABEL=pippo).
 Thus, a completely equivalent result can be obtained with the following shortcut:
-\verbatim
+\plumedfile
 d1: DISTANCE ATOMS=1,2
 PRINT ARG=d1 FILE=colvar STRIDE=10
-\endverbatim
+\endplumedfile
 
 Also notice that all the actions can be labeled, and that many actions besides normal collective variables can define
 one or more value, which can be then referred using the corresponding label.