Reweighting of Accelerated Molecular Dynamics (aMD) Simulations
Background
Accelerated molecular dynamics (aMD) is an enhanced sampling technique that works by flattening the molecular potential energy surface, often by adding a nonnegative boost potential when the system potential is lower than a reference energy. The boost potential, ΔV decreases the energy barriers and thus accelerates transitions between different lowenergy states. With this, aMD is able to sample distinct biomolecular conformations and rare barriercrossing events that are not accessible to conventional molecular dynamics (cMD) simulations. AMD has been successfully applied to a number of biological systems and hundredsofnanosecond aMD simulations have been shown to capture millisecondtimescale events in both globular and membrane proteins.
In addition, it is appealing to reweight aMD trajectories to recover canonical ensemble and the original free energy profile of functional biomolecules. Here, we provide a toolkit of python scripts "PyReweighting" to facilitate the reweighting of aMD simulations. PyReweighting implements a list of commonly used reweighting methods, including (1) exponential average that reweights trajectory frames by the Boltzmann factor of the boost potential and then calculates the ensemble average for each bin, (2) Maclaurin series expansion that approximates the exponential Boltzmann factor, and (3) cumulant expansion that expresses the reweighting factor as summation of boost potential cumulants.
Downloads
PyReweighting1D.py
PyReweighting2D.py
Test Example
testPyReweighting.tgz
Uncompress the file: tar xvzf testPyReweighting.tgz
Tutorial
The following tutorial can also be downloaded as a batch script: run.sh
####################################################
# Prepare input file "weights.dat" in the following format:
# Column 1: dV in units of kbT; column 2: timestep; column 3: dV in units of kcal/mol
# For AMBER14:
# awk 'NR%1==0' amd.log  awk '{print ($8+$7)/(0.001987*300)" " $2 " " ($8+$7)}' > weights.dat
# For AMBER12:
# awk 'NR%1==0' amd.log  awk '{print ($8+$7)" " $3 " " ($8+$7)*(0.001987*300)}' > weights.dat
# For NAMD simulation:
# grep "ACCELERATED MD" namd.log  awk 'NR%1==0'  awk '{print $6/(0.001987*300)" " $4 " " $6 " "$8}' > weights.dat
####################################################
# 1D data
# Prepare input data file "Psi.dat" in one column, e.g., a dihedral angle Psi
# ptraj can be used for AMBER simulations
# Reweighting using cumulant expansion
python PyReweighting1D.py input Psi.dat cutoff 10 Xdim 180 180 disc 6 Emax 20 job amdweight_CE weight weights.dat  tee a reweight_variable.log
mv v pmfc1Psi.dat.xvg pmfPsireweightCE1.xvg
mv v pmfc2Psi.dat.xvg pmfPsireweightCE2.xvg
mv v pmfc3Psi.dat.xvg pmfPsireweightCE3.xvg
# Reweighting using Maclaurin series expansion
python PyReweighting1D.py input Psi.dat disc 6 Emax 20 job amdweight_MC order 10 weight weights.dat  tee a reweight_variable.log
mv v pmfPsi.dat.xvg pmfPsireweightMCorder10.xvg
# Reweighting using exponential average
python PyReweighting1D.py input Psi.dat disc 6 Emax 20 job amdweight weight weights.dat  tee a reweight_variable.log
mv v pmfPsi.dat.xvg pmfPsireweight.xvg
# NOTE: Check out cumulant expansion to the 2nd order "pmfPsireweightCE2.xvg"; normally it gives the most accurate result!
####################################################
# 2D data
# Prepare input data file "Phi_Psi" in two columns
# ptraj can be used for AMBER simulations
# Reweighting using cumulant expansion
python PyReweighting2D.py cutoff 10 input Phi_Psi Xdim 180 180 discX 6 Ydim 180 180 discY 6 Emax 20 job amdweight_CE weight weights.dat  tee a reweight_variable.log
mv v pmfc1Phi_Psi.xvg pmf2DPhi_PsireweightCE1.xvg
mv v pmfc2Phi_Psi.xvg pmf2DPhi_PsireweightCE2.xvg
mv v pmfc3Phi_Psi.xvg pmf2DPhi_PsireweightCE3.xvg
mv v 2D_Free_energy_surface.png pmf2DPhi_PsireweightCE2.png
# Reweighting using Maclaurin series expansion
python PyReweighting2D.py input Phi_Psi Emax 100 discX 6 discY 6 job amdweight_MC order 10 weight weights.dat  tee a reweight_variable.log
mv v pmfPhi_Psi.xvg pmf2DPhi_PsireweightMCorder10disc6.xvg
mv v 2D_Free_energy_surface.png pmf2DPhi_PsireweightMCorder10disc6.png
# Reweighting using exponential average
python PyReweighting2D.py input Phi_Psi Emax 20 discX 6 discY 6 job amdweight weight weights.dat  tee a reweight_variable.log
mv v pmfPhi_Psi.xvg pmf2DPhi_Psireweight.xvg
mv v 2D_Free_energy_surface.png pmf2DPhi_Psireweight.png
NOTES:
1) Maclaurin series "pmf2DPhi_PsireweightMCorder10disc6.png" is equivalent to cumulant expansion on the 1st order "pmf2DPhi_PsireweightCE1.xvg"
2) Check out cumulant expansion to the 2nd order "pmf2DPhi_PsireweightCE2.png"; normally it gives the most accurate result!
3) The above python scripts work for any kind of reaction coordinates, e.g., atom distance, RMSD, or Principal Component Analysis (PCA) modes. You just need to change the default parameters "Xdim 180 180 discX 6 Ydim 180 180 discY 6" to the dimension and bin size of the corresponding reaction coordinates for reweighting.
4) The current reweighting scheme using cumulant expansion to the 2nd order is limited to aMD simulations of small systems, e.g., proteins with 10  40 residues. For larger proteins with more than 100 residues, the energetic noise would be too high for accurate reweighting. Further research has been focused on reducing such energetic noise in simulation of large proteins.
Citation
Miao Y, Sinko W, Pierce L, Bucher D, Walker RC, McCammon JA (2014) Improved reweighting of accelerated molecular dynamics simulations for free energy calculation. J Chemical Theory and Computation, 10(7): 26772689.
Related resources
Scaled Molecular Dynamics (scaled MD) with populationbased reweighting:
http://scaledmd.ucsd.edu/
◇ Last updated:
Dec 30, 2014
◇
