miaolab20 / pyreweighting

PyReweighting

A toolkit of python scripts "PyReweighting" is provided to facilitate the reweighting of accelerated molecular dynamics (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.

Usage

• Update dir_codes in the .sh files to the folder where you save the .py scripts

• Example command lines for using the .sh scripts:

./reweight-1d.sh $Emax $cutoff $binx $data $T

./reweight-2d.sh $Emax $cutoff $binx $biny $data $T

./reweight-3d.sh $Emax $cutoff $binx $biny $binz $data $T

Make sure set temperate for $T, which may have been ignored in previous calculations.

Tutorial

• 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.001987300)}' > 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

./run.sh

NOTES

(1) Maclaurin series "pmf-2D-Phi_Psi-reweight-MC-order10-disc6.png" is equivalent to cumulant expansion on the 1st order "pmf-2D-Phi_Psi-reweight-CE1.xvg"

(2) Check out cumulant expansion to the 2nd order "pmf-2D-Phi_Psi-reweight-CE2.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) For aMD, 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.

(5) GaMD has been developed to reduce the energetic noise in simulations of large systems as noted in (4). Accurate energetic reweighting can be obtained from GaMD simulations of large biomolecular systems (e.g., with millions of atoms) using cumulant expansion to the 2nd order, provided that standard deviation of the boost potential is <= ~ 9 kcal/mol and a total of roughly 0.5 - 1 million frames are saved in the GaMD simulations.

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): 2677-2689.

Miao Y, Bhattarai, A and Wang, J (2020) Ligand Gaussian accelerated molecular dynamics (LiGaMD): Characterization of ligand binding thermodynamics and kinetics, Journal of Chemical Theory and Computation, 16(9): 5526–5547.

pyreweighting's Issues

Extreme PMF values for reweighting DBMD

Hello,
Im trying to use the PyReweighting-2D.py to get the reweighted psi-phi PMF from a DeepBoosted MD simulation. However, the energy values I'm getting are crazy high like 100 Kcal/mol or more.

This is the command Im using to get the input weights file:
awk '!/^#/{if(NR%1==0) print ($8+$7)/(0.001987*300)" " $2 " " ($8+$7)}' ../gamd-0.log > weights.dat

And this is the one I'm using to run the script:
./pyreweighting/reweight-2d.sh 20 10 6 6 coi1_phi-2_psi-2.txt 300

The files gamd-0.log and coi1_phi-2_psi-2.txt are attached

Do you have any idea what's going on?

Thanks

gamd-0.log
coi1_phi-1_psi-1.txt

Some trouble in NAMD input file when I use pyreweighting

PyReweighting test example is outdated, disagrees with latest version of python2 and python3 scripts

On the Miao Lab website that discusses reweighting, there is a set of scripts under the header "Test Example", this provides a tar.xz archive of files from 2013-2014. This test example runs through reweighting using an example Phi_Psi dat file and then shows the expected output XVG and PNG heatmap files from the reweighted Phi/Psi angles.

This issue is specifically about the 2nd order cumulant expansion, not sure if this applies to the other options

However, the Py-Reweighting-2D.py file used to make these XVG and PNG script is outdated, and thus provides completely different results if we try the latest version of the python2 or python3 PyReweighting-2D.py scripts on the same Phi_Psi dat file.

Here is the output given in the test example, if you use the old python script given in the tar archive:

python PyReweighting-2D.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

DATA LOADED: Phi_Psi
dV all: avg = 4.16023073679 std = 2.13785379582
pmf_min-c1 = -8.44361460879
pmf_min-c2 = -15.5787296211
pmf_min-c3 = -304.297343506

Here is the output, with the latest version of the python2 reweighting script. (here, we run into a bug if we try to explicitly set the dimY variable, so I just left that out of the command and let it use the default value of [-180, 180] since that was exactly what I was trying to explicitly set).

python PyReweighting-2D.py -cutoff 10 -input Phi_Psi -Xdim -180 180 -discX 6 -discY 6 -Emax 20 -job amdweight_CE -weight weights.dat | tee -a reweight_variable.log

DATA LOADED: Phi_Psi
bigdata: False
Reweighting with big matrix
hist_max: 684
dV all: avg = 4.16023073679 std = 2.13785379582
pmf_min-c1 = -8.44361460879
pmf_min-c2 = -20.7355467231
pmf_min-c3 = -48.1616340098

Here is the output, with the latest version of the python3 reweighting script. (Ran into a bug "AttributeError: module 'scipy.misc' has no attribute 'factorial'", so changed that scipy.misc.factorial to scipy.math.factorial)

python PyReweighting-2D.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

DATA LOADED: Phi_Psi
dV all: avg = 4.160230736794477 std = 2.1378537958222306
pmf_min-c1 = -8.443614608787296
pmf_min-c2 = -20.73554672305713
pmf_min-c3 = -48.16163400979573

I have also attached the XVG files for the 2nd order cumulant expansion reweighting. The latest python2 and python3 XVG files are pretty much identical except the python3 is rounded to 3 decimal places, and the test example script XVG output is completely different from the python2/python3 XVG output.

pmf-c2-Phi_Psi-PYTHON-3-VERSION.txt
pmf-c2-Phi_Psi-PYTHON-2-VERSION.txt
pmf-c2-Phi_Psi-TEST_EXAMPLE.txt

TL;DR: Looks like the latest versions of the python2 and python3 reweighting scripts are identical, but they both are different from the test example script output, on the same exact data. I suspect that the reweighting algorithm was changed sometime between when the test example was uploaded to the website and when the latest versions of the python2 and python3 scripts were uploaded.

I think it would be good to update that test example folder with the latest Python3 code (and update the latest python3 code to use scipy.math.factorial instead of scipy.misc.factorial). My mentor and I were trying to follow that test example and fully re-capitulate the reweighted heatmap from the Phi/Psi data with the new scripts, and I was very puzzled as to why my results were so different. Updating the test example may help new GaMD users understand the reweighting better with the more modern version of the code.

P.S. I appreciate from the bottom of my heart the hard work that has gone into the PPI-GaMD, Pep-GaMD, Lig-GaMD, and GaMD innovations. It is an incredibly beautiful idea to get around the earlier inconveniences of aMD reweighting. I've learned so much about MD and enhanced sampling from the Miao Lab papers, I'm truly thankful for all your hard work.

GaMD with 4 fs and HMR

Dear,

I want do ask if is possible perform GaMD with HMR (4fs) in NAMD2 o NAMD3.

best.