WIP: trying to make NAMD parser return energies in kT

(see #75)

src/alchemlyb/parsing/namd.py

@@ -5,20 +5,31 @@
 import numpy as np
 from .util import anyopen
 
-def extract_u_nk(fep_file):
+# TODO: perhaps move constants elsewhere?
+# these are the units we need for dealing with NAMD, which uses
+# kcal/mol for energies  http://www.ks.uiuc.edu/Research/namd/2.13/ug/node12.html#SECTION00062200000000000000
+# (kB in kcal/molK)
+k_b = 1.9872041e-3
+
+
+def extract_u_nk(fep_file, T):
     """Return reduced potentials `u_nk` from NAMD fepout file.
 
     Parameters
     ----------
     fep_file : str
         Path to fepout file to extract data from.
+    T : float
+        Temperature in Kelvin at which the simulation was sampled.
 
     Returns
     -------
     u_nk : DataFrame
         Potential energy for each alchemical state (k) for each frame (n).
 
     """
+    beta = 1/(k_b * T)
+
     # lists to get timesteps and work values of each window
     win_ts = []
     win_de = []
@@ -40,7 +51,7 @@ def extract_u_nk(fep_file):
         if '#Free' in l:
 
             # convert last window's work and timestep values to np arrays
-            win_de_arr = np.asarray(win_de)
+            win_de_arr = beta * np.asarray(win_de)
             win_ts_arr = np.asarray(win_ts)
 
             # extract lambda values for finished window
@@ -51,10 +62,10 @@ def extract_u_nk(fep_file):
             # create dataframe of timestep and work values
             # this window's data goes in row LAMBDA1 and column LAMBDA2
             tempDF = pd.DataFrame({
-                'timestep':win_ts_arr,
+                'timestep': win_ts_arr,
                 'fep-lambda': np.full(len(win_de_arr),lambda1),
-                lambda1:0,
-                lambda2:win_de_arr})
+                lambda1: 0,
+                lambda2: win_de_arr})
 
             # join the new window's df to existing df
             u_nk = pd.concat([u_nk, tempDF])

src/alchemlyb/tests/parsing/test_namd.py

@@ -18,7 +18,7 @@ def test_u_nk(dataset, direction, shape):
     """Test that u_nk has the correct form when extracted from files.
     """
     for filename in dataset['data'][direction]:
-        u_nk = extract_u_nk(filename)
+        u_nk = extract_u_nk(filename, T=300)
 
         assert u_nk.index.names == ['timestep', 'fep-lambda']
         assert u_nk.shape == shape

src/alchemlyb/tests/test_fep_estimators.py

@@ -86,8 +86,8 @@ def amber_bace_example_complex_vdw():
 
 def namd_tyr2ala():
     dataset = alchemtest.namd.load_tyr2ala()
-    u_nk1 = namd.extract_u_nk(dataset['data']['forward'][0])
-    u_nk2 = namd.extract_u_nk(dataset['data']['backward'][0])
+    u_nk1 = namd.extract_u_nk(dataset['data']['forward'][0], T=300)
+    u_nk2 = namd.extract_u_nk(dataset['data']['backward'][0], T=300)
 
     # combine dataframes of fwd and rev directions
     u_nk1.replace(0, np.nan, inplace=True)
@@ -139,6 +139,9 @@ class TestBAR(FEPestimatorMixin):
     """
     cls = BAR
 
+    T = 300
+    kT_NAMD = namd.k_b * T
+
     @pytest.mark.parametrize('X_delta_f', [
         (gmx_benzene_coul_u_nk(), 3.044, 0.01640),
         (gmx_benzene_vdw_u_nk(), -3.033, 0.03438),
@@ -149,7 +152,7 @@ class TestBAR(FEPestimatorMixin):
         (gmx_water_particle_with_potential_energy(), -11.724, 0.064964),
         (gmx_water_particle_without_energy(), -11.660, 0.064914),
         (amber_bace_example_complex_vdw(), 2.37846, 0.050899),
-        (namd_tyr2ala(), 6.031269829, 0.069813058),
+        (namd_tyr2ala(), 6.031269829/kT_NAMD, 0.069813058/kT_NAMD),
     ])
     def test_bar(self, X_delta_f):
         self.compare_delta_f(X_delta_f)