Skip to content

Commit

Permalink
WIP: trying to make NAMD parser return energies in kT
Browse files Browse the repository at this point in the history
(see #75)
  • Loading branch information
orbeckst committed Jul 17, 2019
1 parent 1696d94 commit e2a078b
Show file tree
Hide file tree
Showing 3 changed files with 23 additions and 9 deletions.
21 changes: 16 additions & 5 deletions src/alchemlyb/parsing/namd.py
Original file line number Diff line number Diff line change
Expand Up @@ -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 = []
Expand All @@ -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
Expand All @@ -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])
Expand Down
2 changes: 1 addition & 1 deletion src/alchemlyb/tests/parsing/test_namd.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
9 changes: 6 additions & 3 deletions src/alchemlyb/tests/test_fep_estimators.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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),
Expand All @@ -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)
Expand Down

0 comments on commit e2a078b

Please sign in to comment.