Tip
All input files can be downloaded: Files.
Metadynamics for Calculating Free Energy Surface
This tutorial will lead you step by step to do molecular dynamics (MD) simulations with enhanced sampling using Qbics. We recommend you read Standard Molecular Dynamics Simulations first.
Enhanced samping is a method to sample rare events in MD simulations, mostly the overcome of high barrier, including difficult conformational change and chemical reactions. This is not a black-box method and you need some experience to select collective variables (CVs) and sampling methods. In Qbics, enhanced sampling is implemented with Plumed. All CVs and sampling methods in Plumed can be used. They can be combined with (TSO-)DFT, xTB, NDDO, CHARMM, and QM/MM in Qbics. In this tutorial, we will focus on metadynamics, although other methods like OPES can also be used.
Example: Conformation Change of Dialaline
Here, we want to calculate the free energy surface (FES) of 2 dihedral angles (phi and psi) for dialaline with CHARMM force field, see below:
We will use metadynamics, which will add a bias potential to the place it has visited:
Here, \(\mathbf{s}\) is the collective variables, \(t\) is the time, \(\tau\) is the time step, \(W(k \tau)\) is the height of the bias potential, \(\sigma_i\) is the width of the bias potential, \(k\) is the frequency of adding bias potential.
The input file for Qbics is almost same as standard MD simulations:
1charmm
2 parameters par_all22_prot.prm # Parameter files. There can be several ones.
3 topology diala.psf # Topology file.
4 scaling14 1.0 # Scaling factor of 1-4 interactions.
5 rcutoff 12.0 # The cutoff radius for nonbonded interactions.
6 rswitch 10.0
7end
8
9mol
10 diala.pdb
11end
12
13md
14 type NVT # Type: NVE, NVT
15 dt 0.001 # The time step in ps.
16 temp 300. # The temperature in Kelvin.
17 temp_nhc_length 5 # The length of temperature Nose-Hoover chain.
18 temp_nhc_tau 0.2 # The time constant of temperature Nose-Hoover chain in ps.
19 num_steps 10000000 # The number of steps: 10 ns
20 output_freq 100000 # The output frequency: 0.1 ns
21 plumed_fn plumed.dat # Plumed settings file.
22end
23
24task
25 md charmm
26end
The only difference is the addition of plumed_fn in md...end, which is the file that contains the settings for Plumed. With plumed_fn added, enhanced sampling will be automatically activated; if you delete this line, an ordinary MD simulation will be performed.
We can see that a 100 ns NVT MD will be carried out with enhanced sampling set in plumed.dat:
1# Compute the backbone dihedral angle phi, defined by atoms C-N-CA-C
2phi: TORSION ATOMS=5,7,9,15
3psi: TORSION ATOMS=7,9,15,17
4
5metad: METAD ...
6 ARG=phi,psi
7 SIGMA=0.20,0.20
8 HEIGHT=0.0836799999612436
9 #BIASFACTOR=5
10 #TEMP=300.0
11 PACE=100
12 FILE=HILLS
13 GRID_MIN=-pi,-pi
14 GRID_MAX=pi,pi
15 #GRID_BIN=150,150
16...
17
18# Print both collective variables on COLVAR file every 10 steps
19PRINT ARG=phi,psi FILE=COLVAR STRIDE=10
More details about Plumed can refer to Plumed documentation. Here, we give some explainations:
Line 1: A comment.
Line 2: Defines a CV named
phi. TheTORSIONkeyword specifies that this CV is a dihedral angle (the angle between two planes formed by four atoms).ATOMS=5,7,9,15lists the indices (5, 7, 9, 15) of the four atoms in the simulation system that form the phi dihedral angle. See the figure below.Line 3: Similar to Line 2.
Line 5:
metad: METAD ...Starts the definition of a metadynamics simulation block namedmetad(theMETADkeyword initializes metadynamics). The...is a line continuation marker, indicating the configuration for this metadynamics block continues on the next lines.Line 6: Specifies the CVs to be used in metadynamics. Here,
ARG=phi,psimeans the metadynamics simulation will bias the system based on the two dihedral anglesphiandpsi. This is \(\mathbf{s} = (\phi, \psi)\) in the equation above.Line 7: Sets the width (Gaussian width) of the bias potentials added in metadynamics. The two values
0.20,0.20correspond to the \(\sigma_i\) forphiandpsirespectively (units are typically radians for angles), controlling how “spread out” each added Gaussian bias is.Line 8: Defines the height of each Gaussian bias potential added during metadynamics. The value
0.0836799999612436(units are typically in energy units like kilojoules per mole, kJ/mol) determines how strong each individual bias contribution is. This is \(W(k \tau)\) in the equation above.Line 9: Set the
BIASFACTOR(a parameter sometimes used to scale the bias potential) to 5. It is currently commented out.Line 10: Specify the simulation temperature (
TEMP) as300.0Kelvin (K) for metadynamics calculations if uncommented. It is currently commented out.Line 11: Sets the frequency at which Gaussian bias potentials are added, which is \(k\) in the equation above.
PACE=100means a new Gaussian bias is added every 100 MD steps. In our case, 100*0.001 = 0.1 ps.Line 12: Specifies the output file where information about the added Gaussian biases (called “hills” in metadynamics) is stored. The file named
HILLSwill record details like the position, height, and width of each added hill for later analysis.Line 13: Defines the minimum bounds of the grid used for analyzing or visualizing the metadynamics FES. The values
-pi,-pi(in radians) set the lower limit for thephiandpsiCVs respectively (equivalent to -180 degrees).Line 14: Similar to Line 13, but defines the maximum bounds of the FES grid.
Line 15: Set the number of grid bins (resolution) for the FES.
GRID_BIN=150,150would divide thephiandpsiranges into 150 bins each (resulting in a 150x150 grid). It is currently commented out.Line 16: Ends the line continuation for the
metadmetadynamics block, closing the configuration for this section.Line 19: Defines an output section using the
PRINTkeyword.ARG=phi,psispecifies the CVs to print (phiandpsi),FILE=COLVARsets the output file name toCOLVAR, andSTRIDE=10means the CV values are written toCOLVARevery 10 MD steps.
Attention
In Plumed, the default units for length, degree, and energy is nm, radian, and kJ/mol, not Angstrom, degree, and kcal/mol used in Qbics.
With this input files, we can run a 100-ns metadynamics simulation for dialaline.
After the simulation, a file named HILLS will be generated. It saves all \(V(\mathbf{s},t)\) during simulation, with which we can generate free energy data (FES) using the following command:
plumed sum_hills --hills HILLS
The executable of Plumed can be found in plumed/bin in the distribution of Qbics. This will generate the file fes.dat. With the following Python script, we can plot the FES \(F(\phi, \psi)\):
1import numpy as np
2import matplotlib.pyplot as plt
3from matplotlib import cm
4
5def plot_contour_from_file(filename):
6 # Read file.
7 data = np.loadtxt(filename)
8 x = data[:, 0]
9 y = data[:, 1]
10 z = data[:, 2]
11 x_unique = np.unique(x)
12 y_unique = np.unique(y)
13 # Create grids.
14 X, Y = np.meshgrid(x_unique, y_unique)
15 Z = np.zeros_like(X)
16 for i, xi in enumerate(x_unique):
17 for j, yi in enumerate(y_unique):
18 mask = (x == xi) & (y == yi)
19 if np.any(mask):
20 Z[j, i] = z[mask][0]
21 X *= 180/3.1415926 # To degree
22 Y *= 180/3.1415926 # To degree
23 Z *= 0.23900574 # To kcal/mol
24 # Plot.
25 plt.figure(figsize=(10, 8))
26 contourf = plt.contourf(X, Y, Z, levels=25, cmap='jet')
27 contour = plt.contour(X, Y, Z, levels=25, colors='black', linewidths=1)
28 # Style.
29 cbar = plt.colorbar(contourf)
30 cbar.set_label('$F(\phi, \psi)$ (kcal/mol)', fontsize=12)
31 plt.xlabel('$\phi$ (degree)', fontsize=12)
32 plt.ylabel('$\psi$ (degree)', fontsize=12)
33 plt.title('Free Energy Surface $F(\phi, \psi)$', fontsize=14)
34 plt.tight_layout()
35 # Show.
36 plt.show()
37
38if __name__ == "__main__":
39 plot_contour_from_file('fes.dat')
You can find several free energy minima and maxima on the FES.
Example: Proton Transfer of Malonaldehyde
Here we consider the following proton transfer reaction of malonaldehyde:
In you just do a standard MD, you may never see a protron transfering to another oxygen atom. To sample this rare event, we will use metadynamics. The intial structure is shown below. A nature set of CVs to describe the proton transfer is the 2 OH distances, but in this case, we want to plot a one-dimensional FES, so we will use the difference between the 2 distances \(d\equivd_1-d_2\):
The input file is:
1mol
2 H 2.60849612 -0.72353966 0.29921536
3 O 1.96497994 -0.12074426 -0.17095304
4 C 2.66761718 0.78172105 -0.81146959
5 C 4.02335932 0.83968127 -0.80715851
6 C 4.77421680 -0.12520070 -0.05714631
7 O 4.27344987 -1.03185159 0.59780685
8 H 2.06272652 1.49828618 -1.35825644
9 H 4.54176858 1.60740187 -1.35230520
10 H 5.86993567 -0.01026417 -0.09028312
11end
12
13md
14 type NVT # Type: NVE, NVT
15 dt 0.001 # The time step in ps.
16 temp 300. # The temperature in Kelvin.
17 temp_nhc_length 5 # The length of temperature Nose-Hoover chain.
18 temp_nhc_tau 0.2 # The time constant of temperature Nose-Hoover chain in ps.
19 num_steps 1000000 # The number of steps: 1 ns
20 output_freq 1000 # The output frequency: 1 ps
21 plumed_fn plumed.dat
22end
23
24task
25 md xtb
26end
So, this will be a 10-ns NVT metadynamics using xTB, since only quantum chemistry can describe chemical reactions. xTB is fast, but less accurate than DFT.
The Plumed setting is shown below:
1# Compute the OH distances
2d1: DISTANCE ATOMS=1,2
3d2: DISTANCE ATOMS=1,6
4d21: COMBINE ARG=d1,d2 COEFFICIENTS=1,-1 PERIODIC=NO
5
6metad: METAD ...
7 ARG=d21
8 SIGMA=0.05
9 HEIGHT=0.1
10 PACE=100
11 FILE=HILLS
12 GRID_MIN=-0.5
13 GRID_MAX=+0.5
14...
15
16# Print both collective variables on COLVAR file every 10 steps
17PRINT ARG=d1,d2,d21 FILE=COLVAR STRIDE=10
Most settings have already been explained in the last exmaple. Just a few points:
Line 2: Define the CV as distance between atom 1 and 2 named
d1using the keywordDISTANCE.Line 3: Similar to Line 2, but define the CV as distance between atom 1 and 6 named
d2.Line 4:
COMBINEis used to combine the two CVs, whose operation isd1`times1+d2times-1.
Others should be self-explanatory. Now you can run the calculation.
Before calculating the FES, we can visualize the trajectory using malonaldehyde-md-traj.xyz in Qbics-MolStar:
We can see several proton transfer events between the 2 oxygen atoms. Without bias potentials from metadynamics, this is very difficult to happen.
Now, calculate the FES using the following command:
plumed sum_hills --hills HILLS
This will generate the file fes.dat. With the following Python script, we can plot the FES \(F(d)\):
1
Below is the FES plot:
Although there are many local minima, we focus on the transfer region, which is [-1,1] Angstrom:
We show some representative structures in the picuture.
Fiinally, there are 3 points:
The FES should be symmetric, but it is not. This is because the calculation is not completely converged.
The xTB FES is not accurate. One should use DFT to get a more accurate FES in production research.
One can add some restraints in
plumed.datto let \(d\) stay in [-1,1] during the metadynamics, but this needs an unbiasing/reweighting process to modify the FES, which is not discussed in this tutorial and can reter to Plumed documentation.