/
a-master-ISDD-2.txt
504 lines (388 loc) · 24 KB
/
a-master-ISDD-2.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
a-master-ISDD-2.txt/**
\page master-ISDD-2 Master ISDD tutorial 2024: Metadynamics simulations with PLUMED
\section master-ISDD-2-aims Aim
The aim of this tutorial is to train users to perform and analyze metadynamics simulations with PLUMED.
This tutorial has been prepared by Max Bonomi (adapting a lot of material from other tutorials) for
the <a href="http://isddteach.sdv.univ-paris-diderot.fr/fr/accueil.html">Master In Silico Drug Design</a>, held
at Université Paris Cité on April 8th, 2024.
\section master-ISDD-2-objectives Objectives
Once this tutorial is completed users will be able to:
- Write the PLUMED input file to perform metadynamics simulations.
- Calculate the free energy as a function of the metadynamics collective variables.
- Unbias metadynamics simulations.
- Estimate the error in the reconstructed free energies using block analysis.
- Assess the convergence of metadynamics simulations.
\section master-ISDD-2-resources Resources
The \tarball{master-ISDD-2} for this tutorial contains the following files:
- `diala.pdb`: a PDB file for alanine dipeptide in vacuo.
- `topol.tpr`: a GROMACS run (binary) file to perform MD simulations of alanine dipeptide.
- `do_block_fes.py`: a python script to perform error analysis of metadynamics simulations.
After dowloading the compressed archive to your local machine, you can unpack it using the following command:
\verbatim
tar xvzf master-ISDD-2.tar.gz
\endverbatim
Once unpacked, all the files can be found in the `master-ISDD-2` directory. To keep things clean,
it is recommended to run each exercise in a separate sub-directory that you can create inside `master-ISDD-2`.
\note This tutorial has been tested with PLUMED version 2.9.0 and GROMACS version 2020.7.
\section master-ISDD-2-intro Introduction
In the previous tutorial, we have seen that PLUMED can be used to compute collective variables (CVs) on a pre-calculated
trajectory. However, PLUMED is most often use to add forces on the CVs during a MD simulation, for example,
in order to accelerate sampling. To this aim, we have implemented a variety of possible biases acting on CVs.
The complete documentation for all the biasing methods available in PLUMED can be found at the \ref Bias page.
In the following we will learn how to use PLUMED to perform and analyze a metadynamics simulation.
Here you can find a brief recap of the metadynamics theory.
\hidden{Summary of theory}
In metadynamics, an external history-dependent bias potential is constructed in the space of
a few selected degrees of freedom \f$ \vec{s}({q})\f$, generally called collective variables (CVs) \cite metad.
This potential is built as a sum of Gaussian kernels deposited along the trajectory in the CVs space:
\f[
V(\vec{s},t) = \sum_{ k \tau < t} W(k \tau)
\exp\left(
-\sum_{i=1}^{d} \frac{(s_i-s_i({q}(k \tau)))^2}{2\sigma_i^2}
\right).
\f]
where \f$ \tau \f$ is the Gaussian deposition stride,
\f$ \sigma_i \f$ the width of the Gaussian for the \f$i\f$th CV, and \f$ W(k \tau) \f$ the
height of the Gaussian. The effect of the metadynamics bias potential is to push the system away
from local minima into visiting new regions of the phase space. Furthermore, in the long
time limit, the bias potential converges to minus the free energy as a function of the CVs:
\f[
V(\vec{s},t\rightarrow \infty) = -F(\vec{s}) + C.
\f]
In standard metadynamics, Gaussian kernels of constant height are added for the entire course of a
simulation. As a result, the system is eventually pushed to explore high free-energy regions
and the estimate of the free energy calculated from the bias potential oscillates around
the real value.
In well-tempered metadynamics \cite Barducci:2008, the height of the Gaussian
is decreased with simulation time according to:
\f[
W (k \tau ) = W_0 \exp \left( -\frac{V(\vec{s}({q}(k \tau)),k \tau)}{k_B\Delta T} \right ),
\f]
where \f$ W_0 \f$ is an initial Gaussian height, \f$ \Delta T \f$ an input parameter
with the dimension of a temperature, and \f$ k_B \f$ the Boltzmann constant.
With this rescaling of the Gaussian height, the bias potential smoothly converges in the long time limit,
but it does not fully compensate the underlying free energy:
\f[
V(\vec{s},t\rightarrow \infty) = -\frac{\Delta T}{T+\Delta T}F(\vec{s}) + C.
\f]
where \f$ T \f$ is the temperature of the system.
In the long time limit, the CVs thus sample an ensemble
at a temperature \f$ T+\Delta T \f$ which is higher than the system temperature \f$ T \f$.
The parameter \f$ \Delta T \f$ can be chosen to regulate the extent of free-energy exploration:
\f$ \Delta T = 0\f$ corresponds to standard MD, \f$ \Delta T \rightarrow\infty\f$ to standard
metadynamics. In well-tempered metadynamics literature and in PLUMED, you will often encounter
the term "bias factor" which is the ratio between the temperature of the CVs (\f$ T+\Delta T \f$)
and the system temperature (\f$ T \f$):
\f[
\gamma = \frac{T+\Delta T}{T}.
\f]
The bias factor should thus be carefully chosen in order for the relevant free-energy barriers to be crossed
efficiently in the time scale of the simulation.
Additional information can be found in the several review papers on metadynamics
\cite gerv-laio09review \cite WCMS:WCMS31 \cite WCMS:WCMS1103 \cite bussi2015free.
\endhidden
We will play with a toy system, alanine dipeptide simulated in vacuo using the AMBER99SB-ILDN
force field (see Fig. \ref master-ISDD-2-ala-fig).
This rather simple molecule is useful to understand data analysis and free-energy methods.
This system is a nice example because it presents two metastable states separated by a high free-energy barrier.
It is conventional use to characterize the two states in terms of Ramachandran dihedral angles, which are denoted with \f$ \phi \f$
(phi) and \f$ \psi \f$ (psi) in Fig. \ref master-ISDD-2-transition-fig.
\anchor master-ISDD-2-ala-fig
\image html belfast-2-ala.png "The molecule of the day: alanine dipeptide."
\anchor master-ISDD-2-transition-fig
\image html belfast-2-transition.png "Two metastable states of alanine dipeptide are characterized by their Ramachandran dihedral angles."
\section master-ISDD-2-ex Exercises
\subsection master-ISDD-2-ex-1 Exercise 1: My first metadynamics simulation
In this exercise we will setup and perform a well-tempered metadynamics run using the backbone dihedral \f$ \phi \f$
as collective variable. During the calculation, we will also monitor the behavior of the other backbone dihedral \f$ \psi \f$.
Here you can find a sample `plumed.dat` file that you can use as a template.
Whenever you see an highlighted \highlight{FILL} string, this is a string that you must replace.
\plumedfile
# Activate MOLINFO functionalities
MOLINFO STRUCTURE=diala.pdb
# Compute the backbone dihedral angle phi, defined by atoms C-N-CA-C
# you might want to use MOLINFO shortcuts
phi: TORSION ATOMS=__FILL__
# Compute the backbone dihedral angle psi, defined by atoms N-CA-C-N
# here also you might want to use MOLINFO shortcuts
psi: TORSION ATOMS=__FILL__
# Activate well-tempered metadynamics in phi
metad: __FILL__ ARG=__FILL__ ...
# Deposit a Gaussian every 500 time steps, with initial height equal to 1.2 kJ/mol
PACE=500 HEIGHT=1.2
# The bias factor should be wisely chosen
BIASFACTOR=__FILL__
# Gaussian width (sigma) should be chosen based on CV fluctuation in unbiased run
SIGMA=__FILL__
# Gaussians will be written to file and also stored on grid
FILE=HILLS GRID_MIN=-pi GRID_MAX=pi
...
# Print both collective variables on COLVAR file every 10 steps
PRINT ARG=__FILL__ FILE=COLVAR STRIDE=__FILL__
\endplumedfile
The syntax for the command \ref METAD is simple.
The directive is followed by a keyword `ARG` followed by the labels of the CVs
on which the metadynamics bias potential will act.
The keyword `PACE` determines the stride of Gaussian deposition in number of time steps,
while the keyword `HEIGHT` specifies the height of the Gaussian. For each CVs, one has
to specify the width of the Gaussian by using the keyword `SIGMA`. Gaussian will be written
to the file indicated by the keyword `FILE`.
In this example, the bias potential will be stored on a grid, whose boundaries are specified by the keywords `GRID_MIN` and `GRID_MAX`.
Notice that you can provide either the number of bins for every CV (`GRID_BIN`) or
the desired grid spacing (`GRID_SPACING`). In case you provide both PLUMED will use
the most conservative choice (highest number of bins) for each dimension.
In case you do not provide any information about bin size (neither `GRID_BIN` nor `GRID_SPACING`)
and if Gaussian width is fixed, PLUMED will use 1/5 of the Gaussian width as grid spacing.
This default choice should be reasonable for most applications.
Once your `plumed.dat` file is complete, you can run a 10-ns long metadynamics simulations with the following command:
\verbatim
gmx_mpi mdrun -s topol.tpr -nsteps 5000000 -plumed plumed.dat -ntomp 1
\endverbatim
During the metadynamics simulation, PLUMED will create two files, named `COLVAR` and `HILLS`.
The `COLVAR` file contains all the information specified by the \ref PRINT command, in this case
the value of the backbone dihedrals \f$ \phi \f$ and \f$ \psi \f$ every 10 steps of simulation.
We can use `gnuplot` to visualize the behavior of the metadynamics CV \f$ \phi \f$ during the simulation:
\verbatim
gnuplot> p "COLVAR" u 1:2
\endverbatim
\anchor master-ISDD-2-phi-fig
\image html munster-metad-phi.png "Time evolution of the metadynamics CV during the first 2 ns of a metadynamics simulation of alanine dipeptide in vacuum."
By inspecting Figure \ref master-ISDD-2-phi-fig, we can see that the system is initialized in one of the two metastable
states of alanine dipeptide. After a while (t=0.1 ns), the system is pushed
by the metadynamics bias potential to visit the other local minimum. As the simulation continues,
the bias potential fills the underlying free-energy landscape, and the system is able to diffuse in the
entire phase space.
The `HILLS` file contains a list of the Gaussian kernels deposited along the simulation.
If we give a look at the header of this file, we can find relevant information about its content:
\verbatim
#! FIELDS time phi sigma_phi height biasf
#! SET multivariate false
#! SET min_phi -pi
#! SET max_phi pi
\endverbatim
The line starting with `FIELDS` tells us what is displayed in the various columns of the `HILLS` file:
the simulation time, the instantaneous value of \f$ \phi \f$, the Gaussian width and height, and the bias factor.
We can use the `HILLS` file to visualize the decrease of the Gaussian height during the simulation,
according to the well-tempered recipe:
\anchor master-ISDD-2-phihills-fig
\image html munster-metad-phihills.png "Time evolution of the Gaussian height."
If we look carefully at the scale of the y-axis, we will notice that in the beginning the value
of the Gaussian height is higher than the initial height specified in the input file, which should be 1.2 kJ/mol.
In fact, this column reports the height of the Gaussian scaled by the pre-factor that
in well-tempered metadynamics relates the bias potential to the free energy.
\warning The fact that the Gaussian height is decreasing to zero should not be used as a measure of convergence
of your metadynamics simulation!
\subsection master-ISDD-2-ex-2 Exercise 2: Estimating the free energy as a function of the metadynamics CVs
One can estimate the free energy as a function of the metadynamics CVs directly from the metadynamics
bias potential. In order to do so, the utility \ref sum_hills can be used to sum the Gaussian kernels
deposited during the simulation and stored in the `HILLS` file.
To calculate the free energy as a function of \f$ \phi \f$, it is sufficient to use the following command line:
\verbatim
plumed sum_hills --hills HILLS
\endverbatim
The command above generates a file called `fes.dat` in which the free-energy surface as function
of \f$ \phi \f$ is calculated on a regular grid. One can modify the default name for the free-energy file,
as well as the boundaries and bin size of the grid, by using the following \ref sum_hills options:
\verbatim
--outfile - specify the outputfile for sumhills
--min - the lower bounds for the grid
--max - the upper bounds for the grid
--bin - the number of bins for the grid
--spacing - grid spacing, alternative to the number of bins
\endverbatim
The result should look like this:
\anchor master-ISDD-2-metad-phifes-fig
\image html munster-metad-phifes.png "Estimate of the free energy as a function of the dihedral phi from a 10ns-long well-tempered metadynamics simulation."
To give a preliminary assessment of the convergence of a metadynamics simulation, one can calculate the estimate of the free energy as a function
of simulation time. At convergence, the reconstructed profiles should be similar.
The \ref sum_hills option `--stride` should be used to give an estimate of the free energy every `N` Gaussian kernels deposited, and
the option `--mintozero` can be used to align the profiles by setting the global minimum to zero.
If we use the following command line:
\verbatim
plumed sum_hills --hills HILLS --stride 100 --mintozero
\endverbatim
one free energy is calculated every 100 Gaussian kernels deposited, and the global minimum is set to zero in all profiles.
The resulting plot should look like the following:
\anchor master-ISDD-2-metad-phifest-fig
\image html munster-metad-phifest.png "Estimates of the free energy as a function of the dihedral phi calculated every 100 Gaussian kernels deposited."
These two qualitative observations:
1. the system is diffusing rapidly in the entire CV space (Figure \ref master-ISDD-2-phi-fig)
2. the estimated free energy does not significantly change as a function of time (Figure \ref master-ISDD-2-metad-phifest-fig)
suggest that the simulation __might__ be converged.
\warning The two conditions listed above are necessary, but not sufficient to declare convergence.
For a quantitative analysis of the convergence of metadynamics simulations, please have a look below at \ref master-ISDD-2-ex-4.
\subsection master-ISDD-2-ex-3 Exercise 3: Reweighting (unbiasing) a metadynamics simulation
In the previous exercise we biased \f$\phi\f$ and computed the free energy as a function of
the same variable directly from the metadynamics bias potential using the \ref sum_hills utility.
However, in many cases you might decide which variable should be analyzed _after_
having performed a metadynamics simulation. For example,
you might want to calculate the free energy as a function of CVs other than those
biased during the metadynamics simulation, such as the dihedral \f$ \psi \f$.
At variance with standard MD simulations, you cannot simply calculate histograms of other variables directly from your metadynamics trajectory,
because the presence of the metadynamics bias potential has altered the statistical weight of each frame.
To remove the effect of this bias and thus be able to calculate properties of the system in the unbiased ensemble,
you must reweight (unbias) your simulation.
There are multiple ways to calculate the correct
statistical weight of each frame in your metadynamics trajectory and thus to reweight your simulation.
For example:
1. weights can be calculated by considering the time-dependence of the metadynamics bias
potential \cite Tiwary_jp504920s;
2. weights can be calculated using the metadynamics bias potential obtained at the end of the
simulation and assuming a constant bias during the entire course of the simulation \cite Branduardi:2012dl.
In this exercise we will use the second method, which resembles the umbrella-sampling reweighting approach.
In order to compute the weights we will use the \ref driver tool.
First of all, you need to prepare a `plumed_reweight.dat` file that is identical to the one you used
for running your metadynamics simulation except for a few modifications. First, you
need to add the keyword `RESTART=YES` to the \ref METAD command.
This will make this action behave as if PLUMED was restarting, i.e. PLUMED will
read from the `HILLS` file the Gaussians that have previously been accumulated.
Second, you need to set the Gaussian `HEIGHT` to zero and the `PACE` to a large number.
This will actually avoid adding new Gaussians (and even if they are added they will have
zero height). Finally, you need to modify the \ref PRINT statement so that you
write every frame and that, in addition to `phi` and `psi`,
you also write `metad.bias`. You might also want to change the name of the output file to `COLVAR_REWEIGHT`.
\plumedfile
# Activate MOLINFO functionalities
MOLINFO STRUCTURE=diala.pdb
__FILL__ # here goes the definitions of the phi and psi CVs
# Activate well-tempered metadynamics in phi
metad: __FILL__ ARG=__FILL__ ...
# Deposit a Gaussian every 10000000 time steps (never!), with initial height equal to 0.0 kJ/mol
PACE=10000000 HEIGHT=0.0 # <- this is the new stuff!
# The bias factor should be wisely chosen
BIASFACTOR=__FILL__
# Gaussian width (sigma) should be chosen based on CV fluctuation in unbiased run
SIGMA=__FILL__
# Gaussians will be written to file and also stored on grid
FILE=HILLS GRID_MIN=-pi GRID_MAX=pi
# Say that METAD should be restarting (= reading an existing HILLS file)
RESTART=YES # <- this is the new stuff!
...
# Print out the values of phi, psi and the metadynamics bias potential
# Make sure you print out the 3 variables in the specified order at every step
PRINT ARG=__FILL__ FILE=COLVAR_REWEIGHT STRIDE=__FILL__ # <- also change this one!
\endplumedfile
Then run the \ref driver tool using this command:
\verbatim
plumed driver --mf_xtc traj_comp.xtc --plumed plumed_reweight.dat --kt 2.494339
\endverbatim
Notice that you have to specify the value of \f$k_BT\f$ in energy units. While running your simulation
this information was communicated by the MD code.
As a result, PLUMED will produce a new `COLVAR_REWEIGHT` file with one additional column containing the
metadynamics bias potential \f$ V(s) \f$ calculated using all the Gaussians deposited along the entire trajectory.
The beginning of the file should look like this:
\verbatim
#! FIELDS time phi psi metad.bias
#! SET min_phi -pi
#! SET max_phi pi
#! SET min_psi -pi
#! SET max_psi pi
0.000000 -1.497988 0.273498 110.625670
1.000000 -1.449714 0.576594 110.873141
2.000000 -1.209587 0.831417 109.742353
3.000000 -1.475975 1.279726 110.752327
\endverbatim
You can easily obtain the weight \f$ w \f$ of each frame using the expression \f$w\propto\exp\left(\frac{V(s)}{k_BT}\right)\f$
(umbrella-sampling-like reweighting). At this point, you can read the `COLVAR_REWEIGHT` file using, for example, a python script
and compute a weighted histogram. Alternatively, if you want PLUMED to do the weighted histograms for you, you can add the following
lines at the end of the `plumed_reweight.dat` file:
\plumedfile
# Use the metadynamics bias as argument
as: REWEIGHT_BIAS ARG=__FILL__
# Calculate histograms of phi and psi dihedrals every 50 steps
# using the weights obtained from the metadynamics bias potentials (umbrella-sampling-like reweighting)
# Look at the manual to understand the parameters of the HISTOGRAM action!
hhphi: HISTOGRAM ARG=phi STRIDE=50 GRID_MIN=-pi GRID_MAX=pi GRID_BIN=50 BANDWIDTH=0.05 LOGWEIGHTS=as
hhpsi: HISTOGRAM ARG=psi STRIDE=50 GRID_MIN=-pi GRID_MAX=pi GRID_BIN=50 BANDWIDTH=0.05 LOGWEIGHTS=as
# Convert histograms h(s) to free energies F(s) = -kBT * log(h(s))
ffphi: CONVERT_TO_FES GRID=hhphi
ffpsi: CONVERT_TO_FES GRID=hhpsi
# Print out the free energies F(s) to file once the entire trajectory is processed
DUMPGRID GRID=ffphi FILE=ffphi.dat
DUMPGRID GRID=ffpsi FILE=ffpsi.dat
\endplumedfile
and plot the result using `gnuplot`:
\verbatim
gnuplot> p "ffphi.dat" u 1:2 w lp
gnuplot> p "ffpsi.dat" u 1:2 w lp
\endverbatim
You can now compare the free energies as a function of \f$ \phi \f$ calculated:
1. directly from the metadynamics bias potential using \ref sum_hills as done in \ref master-ISDD-2-ex-2;
2. using the reweighting procedure introduced in this exercise.
The results should be identical (see Fig. \ref master-ISDD-2-fescomp-fig).
\anchor master-ISDD-2-fescomp-fig
\image html master-ISDD-2-fescomp-fig.png "Comparison between the free energy as a function of the dihedral phi calculated from the metadynamics bias potential (bias) and by reweighting (rew)".
\subsection master-ISDD-2-ex-4 Exercise 4: Estimating the error in free energies using block-analysis
In the previous exercise, we calculated the _final_ bias \f$ V(s) \f$ on the entire metadynamics trajectory and we used
this quantity to calculate the correct statistical weight of each frame that we need to reweight the biased simulation.
In this exercise, we will see how this information can be used to calculate the error in
the reconstructed free energies and assess whether our simulation is converged or not.
Let's first calculate the un-biasing weights \f$w\propto\exp\left(\frac{V(s)}{k_BT}\right)\f$ from the `COLVAR_REWEIGHT`
file obtained at the end of \ref master-ISDD-2-ex-3:
\verbatim
# Find maximum value of bias to avoid numerical errors when calculating the un-biasing weights
bmax=`awk 'BEGIN{max=0.}{if($1!="#!" && $4>max)max=$4}END{print max}' COLVAR_REWEIGHT`
# Print phi values and un-biasing (un-normalized) weights
awk '{if($1!="#!") print $2,exp(($4-bmax)/kbt)}' kbt=2.494339 bmax=$bmax COLVAR_REWEIGHT > phi.weight
\endverbatim
If you inspect the `phi.weight` file, you will see that each line contains
the value of the dihedral \f$ \phi \f$ along with the corresponding (un-normalized) weight \f$ w \f$ for each frame of the
metadynamics trajectory:
\verbatim
0.907347 0.0400579
0.814296 0.0169656
1.118951 0.0651276
1.040781 0.0714174
1.218571 0.0344903
1.090823 0.0700568
1.130800 0.0622998
\endverbatim
At this point we can apply the block-analysis technique (for more info about the theory,
have a look at \ref trieste-2) to calculate the average free energy across the blocks
and the error as a function of block size. For your convenience, you can use the `do_block_fes.py` python
script to read the `phi.weight` file and produce the desired output.
We use a bash loop to test block sizes ranging from 1 to 1000:
\verbatim
# Arguments of do_block_fes.py
# - input file with CV value and weight for each frame of the trajectory: phi.weight
# - number of CVs: 1
# - CV range (min, max): (-3.141593, 3.141593)
# - # points in output free energy: 51
# - kBT (kJoule/mol): 2.494339
# - Block size: 1<=i<=1000 (every 10)
#
for i in `seq 1 10 1000`; do python3 do_block_fes.py phi.weight 1 -3.141593 3.141593 51 2.494339 $i; done
\endverbatim
For each value of block size `N`, you will obtain a separate `fes.N.dat` file, containing the value
of the \f$ \phi \f$ variable on a grid, the average free energy across the blocks with its associated error (in kJ/mol)
on each point of the grid:
\verbatim
-3.141593 23.184653 0.080659
-3.018393 17.264462 0.055181
-2.895194 13.360259 0.047751
-2.771994 10.772696 0.043548
-2.648794 9.403544 0.042022
\endverbatim
Finally, we can calculate the average error along each free-energy profile as a function of the block size:
\verbatim
for i in `seq 1 10 1000`; do a=`awk '{tot+=$3}END{print tot/NR}' fes.$i.dat`; echo $i $a; done > err.blocks
\endverbatim
and visualize it using `gnuplot`:
\verbatim
gnuplot> p "err.blocks" u 1:2 w lp
\endverbatim
As expected, the error increases with the block size until it reaches a plateau in correspondence of a dimension
of the block that exceeds the correlation between data points (Fig. \ref master-ISDD-2-block-phi).
\anchor master-ISDD-2-block-phi
\image html trieste-4-block-phi.png "Block analysis of a metadynamics simulation using phi as CV"
__What can we learn from this analysis about the convergence of the metadynamics simulation?__
\section master-ISDD-2-conclusions Conclusions
In summary, in this tutorial you should have learned how to use PLUMED to:
- Setup and run a metadynamics calculation.
- Compute free energies from the metadynamics bias potential using the \ref sum_hills utility.
- Reweight a metadynamics simulation.
- Calculate errors and assess convergence.
*/
link: @subpage master-ISDD-2
description: This tutorial explains how to use PLUMED to run metadynamics simulations
additional-files: master-ISDD-2