 |
 |
 |
| |
|
|
Once this tutorial is completed students will learn to:
The tarball for this project contains the following folders:
We consider the association/dissociation of NaCl in aqueous solution. The dissociation barrier is expected to be around 2.5 \(k_\mathrm{B} T\). One interesting aspect of the ion dissociation problem is that collective solvent motions play an important role in the transition. This problem has been considered in the original metadynamics paper [64] and also in reference [14] . We will use the potential developed in ref. [84] for NaCl and TIP3P water with parameters corrected to be used with long-range Coulomb solvers [91]. The system contains 1 Na, 1 Cl, and 106 water molecules (total 320 atoms).
We first perform a standard MD simulation and control the distance Na-Cl. All the files needed for this example are contained in the folder Example1 . The distance Na-Cl can be calculated in Plumed 2 using:
d1: DISTANCE ATOMS=319,320
The coordination number of Na with respect to O in water will also be calculated for later use. This variable will represent the collective motion of the solvent.
COORDINATION ... GROUPA=319 GROUPB=1-318:3 SWITCH={RATIONAL R_0=0.315 D_MAX=0.5 NN=12 MM=24} NLIST NL_CUTOFF=0.55 NL_STRIDE=10 LABEL=coord ... COORDINATION
To run LAMMPS you can use the run.sh script:
#!/bin/bash
############################################################################
# Definition of variables
############################################################################
EXE=lmp_mpi
totalCores=2
############################################################################
mpirun -np ${totalCores} ${EXE} < start.lmp > out.lmpThis command runs LAMMPS using 2 MPI threads. The use of partitions will be discussed when using multiple walkers
Once the simulation is launched, the so called COLVAR file is written. In this case it contains the following:
#! FIELDS time d1 coord 0.200000 0.568067 5.506808 0.400000 0.500148 4.994588 0.600000 0.449778 4.931140 0.800000 0.528272 5.105816 1.000000 0.474371 5.089863 1.200000 0.430620 5.091551 1.400000 0.470374 4.993886 1.600000 0.458768 4.940097 1.800000 0.471886 4.952868 2.000000 0.489058 4.897593 . . .
If you plot the time (column 1) vs the distance (column 2), for instance in gnuplot:
pl "./COLVAR" u 1:2 w lp,
you will see that the ion pair is stuck in the dissociated state during the 1 ns simulation. It is unable to cross the \(\sim5k_\mathrm{B}T\) barrier located at a distance of approximately 0.4 nm. You can also observe this behavior in the trajectory using VMD:
vmd out.dcd -psf nacl.psf
The trajectory has been saved in unwrapped format in order to avoid bonds stretching from one side to the box to the other due to periodic boundary conditions. In VMD we can wrap the atoms without breaking the bonds and show the box using the commands:
pbc wrap -compound res -all pbc box
You can play with different visualization styles and options that VMD has. Therefore, if we want the system to go back and forth between the associated and dissociated state, we will need enhanced sampling.
We now construct a bias potential \( V(\mathbf{s}) \) on the distance Na-Cl using well-tempered metadynamics. The files for this example are contained in the directory Example2. As argument for the construction of the potential we will use the distance Na-Cl (label d1). We choose a gaussian height of 1 kJ/mol which is slightly less than 0.5 \(k_\mathrm{B} T \). The gaussian width is 0.02 nm, in the same order of the features in the FES. A rule of thumb for choosing the gaussian width is to use the standard deviation of the unbiased fluctuations of the CV. The bias factor is set to 5 since the largest barrier in the FES is expected to be roughly 5 \(k_\mathrm{B}T\). Once the metadynamics simulation is converged, the bias will be (up to an arbitrary constant):
\[ V(\mathbf{s})= - \left ( 1- \frac{1}{\gamma} \right ) F(\mathbf{s}) \]
and therefore the system will evolve under an effective free energy:
\[ \tilde{F}(\mathbf{s})=F(\mathbf{s})+V(\mathbf{s})= \frac{F(\mathbf{s})}{\gamma}, \]
that is to say, the largest barrier will be of around 1 \(k_\mathrm{B} T\). The input is:
d1: DISTANCE ATOMS=319,320 METAD ... LABEL=metad ARG=d1 SIGMA=0.02 HEIGHT=1. BIASFACTOR=5 TEMP=300.0 PACE=500 GRID_MIN=0.2 GRID_MAX=1.0 GRID_BIN=300 CALC_RCT ... METAD
Here the CALC_RCT keyword turns on the calculation of the time dependent constant \( c(t) \) that we will use below when reweighting the simulations.
We will also limit the exploration of the CV space by introducing an upper wall bias \( V_{wall}(\mathbf{s}) \):
\[ V_{wall}(s) = \kappa (s-s_0)^2 \mathrm{\: if \:} s>s_0 \mathrm{ \: and \: 0 \: otherwise}. \]
The wall will focus the sampling in the most interesting region of the free energy surface. The effect of this bias potential will have to be corrected later in order to calculate ensemble averages. The syntax in Plumed 2 is:
d1: DISTANCE ATOMS=319,320 UPPER_WALLS ... ARG=d1 AT=0.6 KAPPA=2000.0 EXP=2 EPS=1 OFFSET=0. LABEL=uwall ... UPPER_WALLS
You can run the simulation with the run.sh script as done in the previous example.
It is possible to try different bias factors to check the effect that it has on the trajectory and the effective FES.
In principle the cost of a metadynamics simulation should increase as the simulation progresses due to the need of adding an increasing number of Gaussian kernels to calculate the bias. However, since a grid is used to build up the bias this effect is not observed. You can check what happens if you do not use the GRID_* keywords. Remember that the bins in the grid should be small enough to describe the changes in the bias created by the Gaussian kernels. Normally a bin of size \( \sigma / 5 \) (with \( \sigma \) the gaussian width) is small enough.
One way to identify if a well tempered metadynamics simulation has converged is observing the estimated free energy surface at different times. The FES is estimated by using the relation (again up to an arbitrary constant):
\[ F(\mathbf{s}) = - \left ( \frac{\gamma}{\gamma-1} \right ) V(\mathbf{s},t) \]
and the bias potential is calculated as a sum of Gaussian kernels. This can be done with Plumed 2 using for instance:
plumed sum_hills --hills ../HILLS --min 0.1 --max 0.8 --bin 300 --stride 100
Most of the flags are self explanatory. The flag –stride 100 will result in the FES been written every 100 Gaussian kernels, i.e. 100 ps. Inside the folder FES_calculation you will find a script run.sh that executes the sum_hills command and a gnuplot script plot.gpi that can be used typing:.
gnuplot plot.gpi
After roughly 3 ns the free energy surface does not change significantly with time except for an immaterial constant \( c(t) \) that grows in time. This is in line with well-tempered metadynamics asymptotic behavior:
\[ V(\mathbf{s},t)= - \left ( 1- \frac{1}{\gamma} \right ) F(\mathbf{s}) + c(t). \]
The behavior of \( c(t) \) will be studied with greater detail later.
It should be stressed that we are actually not calculating the free energy \( F(\mathbf{s}) \) but rather \( F(\mathbf{s})+V_{wall}(\mathbf{s}) \). For this reason for distances higher than 0.6 nm the free energy increases sharply.
An alternative way to observe the evolution of the bias is plotting the final free energy plus the instantaneous bias. The scripts for this example can be found in the folder Bias_calculation. In this case we observe only the first 200 ps of the simulation. The (negative) bias can be calculated using:
plumed sum_hills --hills ../HILLS --min 0.1 --max 0.8 --bin 300 --stride 10 --negbias
This plot illustrates clearly how the bias is constructed to progressively "fill" the FES.
It is also possible to track convergence by controlling the evolution of some quantity connected to the free energy surface. In this case we will calculate the dissociation barrier, e.g. the height of the barrier that separates the associated state from the dissociated one. The scripts for this example are found in the folder Barrier_calculation. Using the python script calculate_barrier.py we can compute the barrier, for instance:
The script can be executed using:
python barrier_calculation.py > barrier.txt
and the results can be plotted using the gnuplot script plot.gpi. After roughly 4 ns the barrier stabilizes around 3 and 3.5 \(k_\mathrm{B} T\).
It is important to stress that it is only possible to calculate the free energy difference between two points if the system has gone back and forth between these points several times. This applies both for the calculation of a barrier and the difference in free energy between two basins. It is also important to understand that none of the free energy methods described in this series of tutorials will be able to calculate free energies of regions that have not been sampled, i.e. visited.
Reweighting the simulation on the same CV that was used for biasing can also be used as a test of convergence. We will show that in the next section.
We first reweight the simulation on the distance Na-Cl, the same CV used for biasing. This reweighting is useful to check convergence and to have an estimate of the free energy that does not rely on using kernels. For instance if some features of the FES could not be captured by the kernels, the reweighting procedure will show them. This will become clearer in the VES tutorial. The scripts to perform this calculation are found in the folder ReweightDistance. In metadynamics quasi-stationary limit the weight assigned to a given configuration is [103] :
\[ w(\mathbf{R}) \propto e^{\beta ( V(\mathbf{s},t) - c(t) )}. \]
By plotting time (column 1) versus \( c(t) \) (column 6) using the COLVAR file, the importance of taking \( c(t) \) into account becomes clear. \( c(t) \) keeps growing even after long times, reflecting the approximately rigid shift of the bias with time. Normally the first part of the trajectory is not used for reweighting since during this period the simulation has not reached the quasi-stationary limit. In this case we can discard the first 2 or 3 ns of simulation. To disregard the first 3 ns of simulation we can use sed to delete the first 15000 lines from the COLVAR file:
sed '2,15000d' ../COLVAR > COLVAR
We then use Plumed to calculate two histograms, one taking into account the wall bias and the other one neglecting it. The weights for the reweighting involving only the metadynamics bias have already been discussed while the weights considering both biases are:
\[ w(\mathbf{R}) \propto e^{\beta ( V(\mathbf{s},t) - c(t) + V_{wall}(\mathbf{s}) )}. \]
The header of the COLVAR file should be:
#! FIELDS time d1 coord metad.bias metad.rbias metad.rct metad.work uwall.bias
The input script for Plumed is:
# Read COLVAR file distance: READ FILE=COLVAR IGNORE_TIME VALUES=d1 metad: READ FILE=COLVAR IGNORE_TIME VALUES=metad.rbias uwall: READ FILE=COLVAR IGNORE_TIME VALUES=uwall.bias # Define weights weights1: REWEIGHT_METAD TEMP=300 weights2: REWEIGHT_BIAS TEMP=300 ARG=metad.rbias,uwall.bias # Calculate histograms HISTOGRAM ... ARG=distance GRID_MIN=0.2 GRID_MAX=0.8 GRID_BIN=100 BANDWIDTH=0.002 LOGWEIGHTS=weights1 LABEL=hh1 ... HISTOGRAM HISTOGRAM ... ARG=distance GRID_MIN=0.2 GRID_MAX=0.8 GRID_BIN=100 BANDWIDTH=0.002 LOGWEIGHTS=weights2 LABEL=hh2 ... HISTOGRAM # Print histograms to file DUMPGRID GRID=hh1 FILE=histo FMT=%24.16e DUMPGRID GRID=hh2 FILE=histo_wall FMT=%24.16e
This example can be run with:
plumed --no-mpi driver --plumed plumed.dat --noatoms > plumed.out
and will generate the files histo and histo_wall. The histograms represent the probability \(p(\mathbf{s})\) of observing a given value of the CV \( \mathbf{s} \) . From the histograms the FES can be calculated using:
\[ \beta F(\mathbf{s})= - \log p(\mathbf{s}) \]
and therefore we plot the FES in gnuplot using for instance:
pl "./histo" u 1:(-log($2)) w lp
The next plot compares the estimations of the FES from sum_hills, reweighting with metadynamics bias, and reweighting using both the metadynamics bias and the upper wall bias You will find a gnuplot script plot.gpi to make this plot inside the ReweightDistance folder.
We can obtain important information of the system by reweighting on 2 CVs: The distance Na-Cl and the coordination of Na with O. This reweighting is similar to the one already done and the files that you will need are located in the ReweightBoth folder. Additionally the COLVAR file with the omitted first steps is required. The plot of the FES as a function of these 2 CVs provides important information of the association/dissociation mechanism. In the dissociated state, Na can have a coordination of 5 or 6, though it is more likely to find a coordination number of 6. However, in order to associate Na must have a coordination with O of 5. In the associated state Na can have a coordination of 3, 4 or 5. The transition state is characterized by a coordination number of ~5.
As an exercise, you can write the input files for a simulation in which a bias potential is constructed on the coordination Na-0, i.e. the solvent degree of freedom. You can use the same gaussian height as before and \( \sigma = 0.1 \).
You will find that the exploration of the CV space is not efficient. The reason is that there is a slow degree of freedom that it is not being biased: the distance Na-Cl. Furthermore you can see in the 2 CV reweighting that the coordination Na-O shows significant overlap between the associated and dissociated states.
Bear in mind that this is a rather trivial example since the existing barriers are relatively low. Real problems in materials science usually involve large barriers and are not as forgiving as this example; a bad CV may lead to huge hysteresis and problems in convergence.
We will now construct a bias potential on both CVs. We have already calculated the FES as a function of both CVs through reweighting. In this example the FES will be calculated using the metadynamics bias potential. You can use the input files from Example2.tar and changed the plumed.dat file. To construct a 2 dimensional bias with metadynamics use the following input:
__FILL__ METAD ... LABEL=metad ARG=d1,coord SIGMA=0.02,0.1 HEIGHT=1. BIASFACTOR=5 TEMP=300.0 PACE=500 GRID_MIN=0.15,2. GRID_MAX=0.9,9. GRID_BIN=400,400 CALC_RCT ... METAD
Once that the simulation is completed you can run plumed sum_hills to calculate the FES:
plumed sum_hills --hills HILLS --mintozero
and plot the results using the following lines in gnuplot:
set pm3d map set zr [0:15] spl "fes.dat" u 1:2:3
Some valuable tools for metadynamics simulations will be discussed in the VES tutorial. These include:
1.8.14