Bucket Sampling

Introduction

In this tutorial we will learn to use bucket sampling (BS) for efficient sampling of high dimensional free energy landscape. This tutorial follows the method used in the paper: Javed, R., Kapakayala, A. B., & Nair, N. N., Buckets instead of umbrellas for enhanced sampling and free energy calculations. Journal of Chemical Theory and Computation, 20(19), 2024. It is recommended to read this paper for the details of parameters used for the simulation.

We study the free energy surface of alanine dipeptide in vacuo as a function of two backbone dihedral angles (Φ, Ψ) using BS. Here, we are assuming that users are familiar with GROMACS and can run MD simulation using GROMACS-PLUMED interface. The first thing we need to do is build an initial structure using ff14SB force-field and relax this at 300K so that we can use this as the input for our simulation.  

This tutorial consist of two parts:

• Part 1: BS Simulation

• Part 2: Reweighting

Here, we will run the simulation for six windows. The full range for our interest is from -3.14 rad to 3.14 rad along Φ. Thus, we divide the full region into 6 buckets  each spanning a region of 1 rad. Thus, we have the η_L and η_U values:

1. BS Simulation

a) Determination of the position of the WALL potentials for every bucket

Before starting the BS simulations, it is essential to estimate the collective variable (CV) diffusion length, λ, along which the bucket bias is to be applied. This step ensures that overly rigid wall potentials do not adversely affect the simulation or introduce artifacts.

Steps for computing λ:

1. Preliminary Metadynamics Simulation:
   Conduct short, well-tempered metadynamics (WT-MetaD) simulation starting with the equilibrium structure using the same set of parameters that we plan to use for the bucket sampling. WT-MetaD will be applied on both Φ, Ψ CVs. This simulation will help capture the natural diffusion behavior of the system across the chosen CV.

2. Diffusion Analysis:
   Use the script compute_lambda.py to calculate the diffusion coefficient from each independent simulation. Before running this Python code, you need to edit the following:

CV_COL = 1
cv_is_torsion= True # Make it False if not a torsion

 ETA_L=np.array([-3.14, -1.05, 1.05]) #\eta_L

 ETA_U=np.array([-1.05,  1.05, 3.14]) #\eta_u

CV_COL is the index of the CV along which the bucket potential is applied.
cv_is_torsion=True if the bucket-CV is a torsion.
ETA_L & ETA_U: List all the values of η_L and η_U are to be entered.

The code will compute the ξ_L and ξ_U values and print them. ξ_L and ξ_U are the positions where the wall potentials are to be positioned:

For the alanine dipeptide in water simulation, we get λ=0.4 rad; thus, we eliminate 0.2 rad from the right and left regions to consider any possibility of boundary effects. See the table below.

The main output files of these simulations which are required for analysis are:

1. COLVAR :  Collective Variables files for different windows

2. HILLS : List of added bias is stored in this file for different windows

2. Reweighting

The second part of our calculation is to use probability generated by BS simulation to get the free energy surface. But in our simulation we have applied MTD bias in all the 6 windows which will give biased probability distribution. We need to get unbiased probability distribution by reweighting and then patching all the windows using mean force based reweighting technique (Pal et al., 2021).

For reweighting, use the Py-bucket-reweight.py code. The manual for the reweighting code is here (PDF). 

NOTE: For the reconstruction part, only the η region for each window needs to be taken.

Use gnuplot to plot free_energy_2D.dat to get the free energy as below.