python script to calculate the free energy #100
Comments
|
Hi @qasimpars, and welcome to alchemlyb! You can write a script like the following to get out what you're looking for: #!/usr/bin/env python
import pandas as pd
from alchemlyb.parsing.gmx import extract_dHdl, extract_u_nk
from alchemlyb.preprocessing import statistical_inefficiency
from alchemlyb.estimators import BAR, MBAR, TI
if __name__ == '__main__':
import glob
import os
import argparse
parser = argparse.ArgumentParser(description="Perform a simple alchemical analysis")
parser.add_argument('datadir', type=str,
help="directory of XVG files")
parser.add_argument('-T', type=float,
help="Temperature (K) simulations performed at")
args = parser.parse_args()
# get a list of files from the given directory
files = glob.glob(os.path.join(args.datadir, '*.xvg*'))
# extract and subsample dHdl using statistical inefficiency
print("Gathering dHdl from {} files".format(len(files)))
dHdl = [extract_dHdl(i, T=args.T) for i in files]
print("Subsampling dHdl on statistical inefficiency")
dHdl = pd.concat([statistical_inefficiency(i, i.iloc[:, 0]) for i in dHdl])
# extract and subsample u_nk using statistical inefficiency
print("Gathering u_nk from {} files".format(len(files)))
u_nk = [extract_u_nk(i, T=args.T) for i in files]
print("Subsampling u_nk on statistical inefficiency")
u_nk = pd.concat([statistical_inefficiency(i, i.iloc[:, 0]) for i in u_nk])
# fit dHdl estimators
print("Fitting TI on dHdl")
ti = TI().fit(dHdl)
# fit u_nk estimators
print("Fitting BAR on u_nk")
bar = BAR().fit(u_nk)
print("Fitting MBAR on u_nk")
mbar = MBAR().fit(u_nk)
print("====== Results ======")
# print free energy estimates between lowest lambda and highest lambda
# these will all be in units of kT!
# NOTE: for BAR we would need to perform bootstrap sampling to get
# uncertainties of the endpoints
print("TI: {} +/- {} kT".format(ti.delta_f_.iloc[0, -1], ti.d_delta_f_.iloc[0, -1]))
print("BAR: {} +/- {} kT".format(bar.delta_f_.iloc[0, -1], "unknown"))
print("MBAR: {} +/- {} kT".format(mbar.delta_f_.iloc[0, -1], mbar.d_delta_f_.iloc[0, -1]))You can then call this script on your directory with, e.g.: Please note that currently we don't give reliable estimates of the uncertainty for BAR, so I have left those out of the output. We are working to address this by adding an implementation of bootstrap sampling (#80). You'll want to make sure your files have the necessary data to support both |
|
Hi @dotsdl, Thanks for your reply and script. They have been installed properly. Then I ran the script: But the script gives the below error:
It gives the following error about alchemlyb when I delete
Do you know how I can fix those errors? My OS is ubuntu 18.04. I reinstalled the alchemlyb and pandas but it didn't help. Edit: Now the script still gives the below error:
|
|
Hi @qasimpars, It may be that your and paste the output of each command here? |
|
If I recall on Ubuntu 18.04 |
|
Hi @dotsdl, I used the
|
|
Ah, I believe you are running into #95. We need to get #98 finished and merged soon to address this. If you are feeling adventurous, you can clone the repository and switch to the |
|
Hi, I did what you said. It seems that the script works but it gives the below error. Do you know how I can fix the error about MBAR?
Edit: |
|
It would be really good if the following flags are implemented in the script also: |
|
Hmmm, I think it may be worth trying the contents of the quick script above out in a Jupyter notebook. This way, you can run individual components and have a look at the intermediates, such as what As for adding more flags, our eventual goal is to create a replacement for |
|
Hello @dotsdl , Edit:
|
|
I hope to close this issue with #114 |
|
I am closing this issue (in the light of discussion #159 )
|
Hi,
Is there an example python script to calculate the free energy with MBAR, BAR and TI from the dhdl.xvg files of GROMACS via alchemlyb?
Please note that all the dhdl.xvg files obtained from GROMACS are in a directory named dHdl_files.
Hope that you will help me.
Best,
The text was updated successfully, but these errors were encountered: