From 6227d83b747ff38ac71baa2cfbb0ffb7c83af3de Mon Sep 17 00:00:00 2001 From: Tom Joseph Date: Thu, 22 Jul 2021 14:00:06 -0400 Subject: [PATCH 1/9] Tolerate truncated and restarted NAMD fepout files --- src/alchemlyb/parsing/namd.py | 120 ++++++++++++++++++++++++++-------- 1 file changed, 93 insertions(+), 27 deletions(-) diff --git a/src/alchemlyb/parsing/namd.py b/src/alchemlyb/parsing/namd.py index e3bdb068..eab8452b 100644 --- a/src/alchemlyb/parsing/namd.py +++ b/src/alchemlyb/parsing/namd.py @@ -9,6 +9,63 @@ k_b = R_kJmol * kJ2kcal +def get_lambdas(fep_files): + """Retrieves all lambda values included in the FEP files provided. + + We have to do this in order to tolerate truncated and restarted fepout files. + The IDWS lambda is not present at the termination of the window, presumably + for backwards compatibility with ParseFEP and probably other things. + + For a given lambda1, there can be only one lambda_idws. + + Parameters + ---------- + fep_files: str or list of str + Path(s) to fepout files to extract dasta from. + + Returns + ------- + List of floats, or None if there is more than one lambda_idws for each lambda1. + """ + + lambda_fwd_map, lambda_bwd_map = {}, {} + + for fep_file in sorted(fep_files): + with anyopen(fep_file, 'r') as f: + for line in f: + l = line.strip().split() + if l[0] not in ['#NEW', '#Free']: + continue + + # We might not have a #NEW line so make the best guess + if l[0] == '#NEW': + lambda1, lambda2 = l[6], l[8] + lambda_idws = l[10] if 'LAMBDA_IDWS' in l else None + elif l[0] == '#Free': + lambda1, lambda2, lambda_idws = l[7], l[8], None + + # Make sure the lambda2 values are consistent + if lambda1 in lambda_fwd_map and lambda_fwd_map[lambda1] != lambda2: + print(f'fwd: lambda1 {lambda1} has lambda2 {lambda_fwd_map[lambda1]} instead of {lambda2}') + return None + + lambda_fwd_map[lambda1] = lambda2 + + # Make sure the lambda_idws values are consistent + if lambda_idws is not None: + if lambda1 in lambda_bwd_map and lambda_bwd_map[lambda1] != lambda_idws: + print(f'bwd: lambda1 {lambda1} has lambda_idws {lambda_bwd_map[lambda1]} instead of {lambda_idws}') + return None + lambda_bwd_map[lambda1] = lambda_idws + + all_lambdas = set() + all_lambdas.update(lambda_fwd_map.keys()) + all_lambdas.update(lambda_fwd_map.values()) + all_lambdas.update(lambda_bwd_map.keys()) + all_lambdas.update(lambda_bwd_map.values()) + return list(sorted(all_lambdas)) + + @_init_attrs def extract_u_nk(fep_files, T): """Return reduced potentials `u_nk` from NAMD fepout file(s). @@ -25,18 +82,6 @@ def extract_u_nk(fep_files, T): u_nk : DataFrame Potential energy for each alchemical state (k) for each frame (n). - Note - ---- - If the number of forward and backward samples in a given window are different, - the extra sample(s) will be discarded. This is typically zero or one sample. - - .. versionchanged:: 0.5.0 - The :mod:`scipy.constants` is used for parsers instead of - the constants used by the corresponding MD engine. - - Support for Interleaved Double-Wide Sampling files added. - - `fep_files` can now be a list of filenames. """ beta = 1/(k_b * T) @@ -57,51 +102,73 @@ def extract_u_nk(fep_files, T): fep_files = [fep_files] time = 0 + # Extract the lambda values only from the fepouts + all_lambdas = get_lambdas(fep_files) # open and get data from fep file. - for fep_file in fep_files: + # We sort the list of fep files in case some of them represent restarted windows. + # The assumption is that they make sense in lexicographic order. + for fep_file in sorted(fep_files): + lambda1_at_start, lambda2_at_start, lambda_idws_at_start = None, None, None with anyopen(fep_file, 'r') as f: + has_idws = False for line in f: l = line.strip().split() + # We don't know if IDWS was enabled just from the #Free line, and we might not have + # a #NEW line in this file, so we have to check for the existence of FepE_back lines + # We rely on short-circuit evaluation to avoid the string comparison most of the time + if has_idws is False and l[0] == 'FepE_back:': + has_idws = True # New window, get IDWS lambda if any + # We keep track of lambdas from the #NEW line and if they disagree with the #Free line + # within the same file, then complain. This can happen if truncated fepout files + # are presented in the wrong order. if l[0] == '#NEW': - if 'LAMBDA_IDWS' in l: - lambda_idws = l[10] - else: - lambda_idws = None + lambda1_at_start, lambda2_at_start = l[6], l[8] + lambda_idws_at_start = l[10] if 'LAMBDA_IDWS' in l else None # this line marks end of window; dump data into dataframe if '#Free' in l: - # extract lambda values for finished window # lambda1 = sampling lambda (row), lambda2 = comparison lambda (col) lambda1 = l[7] lambda2 = l[8] + lambda1_idx = all_lambdas.index(lambda1) + if has_idws is True and lambda1_idx > 0: + lambda_idws = all_lambdas[lambda1_idx - 1] + else: + lambda_idws = None + + # If the lambdas are not what we thought they would be, return None, ensuring the calculation + # fails. + if (lambda1, lambda2, lambda_idws) != (lambda1_at_start, lambda2_at_start, lambda_idws_at_start): + print("Error: Lambdas appear to have changed within the same fepout file", fep_file) + return None # convert last window's work and times values to np arrays win_de_arr = beta * np.asarray(win_de) win_ts_arr = np.asarray(win_ts) if lambda_idws is not None: - # Mimic classic DWS data - # Arbitrarily match up fwd and bwd comparison energies on the same times - # truncate extra samples from whichever array is longer + # Mimic classic DWS data + # Arbitrarily match up fwd and bwd comparison energies on the same times + # truncate extra samples from whichever array is longer win_de_back_arr = beta * np.asarray(win_de_back) n = min(len(win_de_back_arr), len(win_de_arr)) tempDF = pd.DataFrame({ - 'time': win_ts_arr[:n], - 'fep-lambda': np.full(n,lambda1), + 'time': win_ts_arr[:n], + 'fep-lambda': np.full(n,lambda1), lambda1: 0, - lambda2: win_de_arr[:n], - lambda_idws: win_de_back_arr[:n]}) + lambda2: win_de_arr[:n], + lambda_idws: win_de_back_arr[:n]}) else: # create dataframe of times and work values # this window's data goes in row LAMBDA1 and column LAMBDA2 tempDF = pd.DataFrame({ 'time': win_ts_arr, - 'fep-lambda': np.full(len(win_de_arr),lambda1), + 'fep-lambda': np.full(len(win_de_arr), lambda1), lambda1: 0, lambda2: win_de_arr}) @@ -137,7 +204,6 @@ def extract_u_nk(fep_files, T): tempDF = pd.DataFrame({ 'time': win_ts_arr, 'fep-lambda': lambda2}) - u_nk = pd.concat([u_nk, tempDF], sort=True) u_nk.set_index(['time','fep-lambda'], inplace=True) From c1e0d9c6cf3c437890ada84e5526235b452c4940 Mon Sep 17 00:00:00 2001 From: Tom Joseph Date: Thu, 22 Jul 2021 14:18:44 -0400 Subject: [PATCH 2/9] Better tolerate when fepout file has no #NEW line --- src/alchemlyb/parsing/namd.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/alchemlyb/parsing/namd.py b/src/alchemlyb/parsing/namd.py index eab8452b..cadc3f78 100644 --- a/src/alchemlyb/parsing/namd.py +++ b/src/alchemlyb/parsing/namd.py @@ -141,8 +141,10 @@ def extract_u_nk(fep_files, T): # If the lambdas are not what we thought they would be, return None, ensuring the calculation # fails. - if (lambda1, lambda2, lambda_idws) != (lambda1_at_start, lambda2_at_start, lambda_idws_at_start): + if lambda1_at_start is not None \ + and (lambda1, lambda2, lambda_idws) != (lambda1_at_start, lambda2_at_start, lambda_idws_at_start): print("Error: Lambdas appear to have changed within the same fepout file", fep_file) + print(f"{lambda1_at_start} {lambda2_at_start} {lambda_idws_at_start} => {lambda1} {lambda2} {lambda_idws}") return None # convert last window's work and times values to np arrays From d4620bea3dd3cf6661f55d155418ff0ee8a4190d Mon Sep 17 00:00:00 2001 From: Tom Joseph Date: Thu, 29 Jul 2021 15:00:16 -0400 Subject: [PATCH 3/9] Explicitly parse lambda value strings to float. And, cosmetic changes --- src/alchemlyb/parsing/namd.py | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/alchemlyb/parsing/namd.py b/src/alchemlyb/parsing/namd.py index cadc3f78..ed80fa8f 100644 --- a/src/alchemlyb/parsing/namd.py +++ b/src/alchemlyb/parsing/namd.py @@ -21,7 +21,7 @@ def get_lambdas(fep_files): Parameters ---------- fep_files: str or list of str - Path(s) to fepout files to extract dasta from. + Path(s) to fepout files to extract data from. Returns ------- @@ -39,14 +39,14 @@ def get_lambdas(fep_files): # We might not have a #NEW line so make the best guess if l[0] == '#NEW': - lambda1, lambda2 = l[6], l[8] - lambda_idws = l[10] if 'LAMBDA_IDWS' in l else None + lambda1, lambda2 = float(l[6]), float(l[8]) + lambda_idws = float(l[10]) if 'LAMBDA_IDWS' in l else None elif l[0] == '#Free': - lambda1, lambda2, lambda_idws = l[7], l[8], None + lambda1, lambda2, lambda_idws = float(l[7]), float(l[8]), None # Make sure the lambda2 values are consistent if lambda1 in lambda_fwd_map and lambda_fwd_map[lambda1] != lambda2: - print(f'fwd: lambda1 {lambda1} has lambda2 {lambda_fwd_map[lambda1]} instead of {lambda2}') + print(f'namd.py: get_lambdas: Error: fwd: lambda1 {lambda1} has lambda2 {lambda_fwd_map[lambda1]} but it should be {lambda2}') return None lambda_fwd_map[lambda1] = lambda2 @@ -54,7 +54,7 @@ def get_lambdas(fep_files): # Make sure the lambda_idws values are consistent if lambda_idws is not None: if lambda1 in lambda_bwd_map and lambda_bwd_map[lambda1] != lambda_idws: - print(f'bwd: lambda1 {lambda1} has lambda_idws {lambda_bwd_map[lambda1]} instead of {lambda_idws}') + print(f'namd.py: get_lambdas: Error: bwd: lambda1 {lambda1} has lambda_idws {lambda_bwd_map[lambda1]} but it should be {lambda_idws}') return None lambda_bwd_map[lambda1] = lambda_idws @@ -124,15 +124,16 @@ def extract_u_nk(fep_files, T): # within the same file, then complain. This can happen if truncated fepout files # are presented in the wrong order. if l[0] == '#NEW': - lambda1_at_start, lambda2_at_start = l[6], l[8] - lambda_idws_at_start = l[10] if 'LAMBDA_IDWS' in l else None + lambda1_at_start, lambda2_at_start = float(l[6]), float(l[8]) + lambda_idws_at_start = float(l[10]) if 'LAMBDA_IDWS' in l else None + has_idws = True if lambda_idws_at_start is not None else False # this line marks end of window; dump data into dataframe if '#Free' in l: # extract lambda values for finished window # lambda1 = sampling lambda (row), lambda2 = comparison lambda (col) - lambda1 = l[7] - lambda2 = l[8] + lambda1 = float(l[7]) + lambda2 = float(l[8]) lambda1_idx = all_lambdas.index(lambda1) if has_idws is True and lambda1_idx > 0: lambda_idws = all_lambdas[lambda1_idx - 1] @@ -143,8 +144,9 @@ def extract_u_nk(fep_files, T): # fails. if lambda1_at_start is not None \ and (lambda1, lambda2, lambda_idws) != (lambda1_at_start, lambda2_at_start, lambda_idws_at_start): - print("Error: Lambdas appear to have changed within the same fepout file", fep_file) - print(f"{lambda1_at_start} {lambda2_at_start} {lambda_idws_at_start} => {lambda1} {lambda2} {lambda_idws}") + print("namd.py: extract_u_nk: Error: Lambdas appear to have changed within the same namd fepout file", fep_file) + print(f"namd.py: extract_u_nk: Error: l1, l2, lidws: {lambda1_at_start}, {lambda2_at_start}, {lambda_idws_at_start} changed to {lambda1}, {lambda2}, {lambda_idws}") + print(f"namd.py: extract_u_nk: Error: fep_file = {fep_file}; has_idws = {has_idws}") return None # convert last window's work and times values to np arrays @@ -197,10 +199,10 @@ def extract_u_nk(fep_files, T): if '#STARTING' in l: parsing = True - if (len(win_de) != 0 or len(win_de_back) != 0): + if len(win_de) != 0 or len(win_de_back) != 0: print('Warning: trailing data without footer line (\"#Free energy...\"). Interrupted run?') - if (float(lambda2) == 1.0 or float(lambda2) == 0.0): + if lambda2 in (0.0, 1.0): # this excludes the IDWS case where a dataframe already exists for both endpoints # create last dataframe for fep-lambda at last LAMBDA2 tempDF = pd.DataFrame({ From 11b01779f01920d44d987b9767e0661d79000d96 Mon Sep 17 00:00:00 2001 From: Tom Joseph Date: Sat, 31 Jul 2021 10:52:50 -0400 Subject: [PATCH 4/9] More descriptive comments --- src/alchemlyb/parsing/namd.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/src/alchemlyb/parsing/namd.py b/src/alchemlyb/parsing/namd.py index ed80fa8f..d081b3d6 100644 --- a/src/alchemlyb/parsing/namd.py +++ b/src/alchemlyb/parsing/namd.py @@ -109,6 +109,8 @@ def extract_u_nk(fep_files, T): # The assumption is that they make sense in lexicographic order. for fep_file in sorted(fep_files): lambda1_at_start, lambda2_at_start, lambda_idws_at_start = None, None, None + # Note we have not set parsing=False because we could be continuing one window across + # more than one fepout file with anyopen(fep_file, 'r') as f: has_idws = False for line in f: @@ -129,7 +131,7 @@ def extract_u_nk(fep_files, T): has_idws = True if lambda_idws_at_start is not None else False # this line marks end of window; dump data into dataframe - if '#Free' in l: + if l[0] == '#Free': # extract lambda values for finished window # lambda1 = sampling lambda (row), lambda2 = comparison lambda (col) lambda1 = float(l[7]) @@ -141,22 +143,23 @@ def extract_u_nk(fep_files, T): lambda_idws = None # If the lambdas are not what we thought they would be, return None, ensuring the calculation - # fails. + # fails. This can happen if fepouts where one window spans multiple fepouts are processed out of order if lambda1_at_start is not None \ and (lambda1, lambda2, lambda_idws) != (lambda1_at_start, lambda2_at_start, lambda_idws_at_start): - print("namd.py: extract_u_nk: Error: Lambdas appear to have changed within the same namd fepout file", fep_file) + print("namd.py: extract_u_nk: Error: Lambdas changed unexpectedly while processing", fep_file) print(f"namd.py: extract_u_nk: Error: l1, l2, lidws: {lambda1_at_start}, {lambda2_at_start}, {lambda_idws_at_start} changed to {lambda1}, {lambda2}, {lambda_idws}") print(f"namd.py: extract_u_nk: Error: fep_file = {fep_file}; has_idws = {has_idws}") return None - # convert last window's work and times values to np arrays - win_de_arr = beta * np.asarray(win_de) - win_ts_arr = np.asarray(win_ts) + # As we are at the end of a window, convert last window's work and times values to np arrays + # (with energy unit kT since they were kcal/mol in the fepouts) + win_de_arr = beta * np.asarray(win_de) # dE values + win_ts_arr = np.asarray(win_ts) # timesteps if lambda_idws is not None: - # Mimic classic DWS data - # Arbitrarily match up fwd and bwd comparison energies on the same times - # truncate extra samples from whichever array is longer + # Mimic classic DWS data + # Arbitrarily match up fwd and bwd comparison energies on the same times + # truncate extra samples from whichever array is longer win_de_back_arr = beta * np.asarray(win_de_back) n = min(len(win_de_back_arr), len(win_de_arr)) @@ -195,7 +198,7 @@ def extract_u_nk(fep_files, T): win_de_back.append(float(l[6])) win_ts_back.append(float(l[1])) - # turn parsing on after line 'STARTING COLLECTION OF ENSEMBLE AVERAGE' + # Turn parsing on after line 'STARTING COLLECTION OF ENSEMBLE AVERAGE' if '#STARTING' in l: parsing = True From 1d8b0476be676591b73ed325d649365633e54708 Mon Sep 17 00:00:00 2001 From: Tom Joseph Date: Thu, 12 Aug 2021 14:56:45 -0400 Subject: [PATCH 5/9] Add jhenin and myself to AUTHORS --- AUTHORS | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/AUTHORS b/AUTHORS index c48f9fa2..60383a9b 100644 --- a/AUTHORS +++ b/AUTHORS @@ -35,4 +35,8 @@ Chronological list of authors - Hyungro Lee (@lee212) - Mohammad S. Barhaghi (@msoroush) 2020 - - Zhiyi Wu (@xiki-tempula) \ No newline at end of file + - Zhiyi Wu (@xiki-tempula) + +2021 + - Jérôme Hénin (@jhenin) + - Thomas T. Joseph (@ttjoseph) From 1cbc36feb3b121c221378bac7103154f48409ef6 Mon Sep 17 00:00:00 2001 From: Tom Joseph Date: Thu, 12 Aug 2021 15:07:36 -0400 Subject: [PATCH 6/9] Use logger instead of print() --- src/alchemlyb/parsing/namd.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/alchemlyb/parsing/namd.py b/src/alchemlyb/parsing/namd.py index d081b3d6..a801b3d5 100644 --- a/src/alchemlyb/parsing/namd.py +++ b/src/alchemlyb/parsing/namd.py @@ -3,13 +3,17 @@ """ import pandas as pd import numpy as np +import logging from .util import anyopen from . import _init_attrs from ..postprocessors.units import R_kJmol, kJ2kcal +logger = logging.getLogger("alchemlyb.parsers.NAMD") + k_b = R_kJmol * kJ2kcal -def get_lambdas(fep_files): + +def _get_lambdas(fep_files): """Retrieves all lambda values included in the FEP files provided. We have to do this in order to tolerate truncated and restarted fepout files. @@ -46,7 +50,7 @@ def get_lambdas(fep_files): # Make sure the lambda2 values are consistent if lambda1 in lambda_fwd_map and lambda_fwd_map[lambda1] != lambda2: - print(f'namd.py: get_lambdas: Error: fwd: lambda1 {lambda1} has lambda2 {lambda_fwd_map[lambda1]} but it should be {lambda2}') + logger.error(f'fwd: lambda1 {lambda1} has lambda2 {lambda_fwd_map[lambda1]} but it should be {lambda2}') return None lambda_fwd_map[lambda1] = lambda2 @@ -54,7 +58,7 @@ def get_lambdas(fep_files): # Make sure the lambda_idws values are consistent if lambda_idws is not None: if lambda1 in lambda_bwd_map and lambda_bwd_map[lambda1] != lambda_idws: - print(f'namd.py: get_lambdas: Error: bwd: lambda1 {lambda1} has lambda_idws {lambda_bwd_map[lambda1]} but it should be {lambda_idws}') + logger.error(f'namd.py: get_lambdas: Error: bwd: lambda1 {lambda1} has lambda_idws {lambda_bwd_map[lambda1]} but it should be {lambda_idws}') return None lambda_bwd_map[lambda1] = lambda_idws @@ -103,7 +107,7 @@ def extract_u_nk(fep_files, T): time = 0 # Extract the lambda values only from the fepouts - all_lambdas = get_lambdas(fep_files) + all_lambdas = _get_lambdas(fep_files) # open and get data from fep file. # We sort the list of fep files in case some of them represent restarted windows. # The assumption is that they make sense in lexicographic order. @@ -146,9 +150,9 @@ def extract_u_nk(fep_files, T): # fails. This can happen if fepouts where one window spans multiple fepouts are processed out of order if lambda1_at_start is not None \ and (lambda1, lambda2, lambda_idws) != (lambda1_at_start, lambda2_at_start, lambda_idws_at_start): - print("namd.py: extract_u_nk: Error: Lambdas changed unexpectedly while processing", fep_file) - print(f"namd.py: extract_u_nk: Error: l1, l2, lidws: {lambda1_at_start}, {lambda2_at_start}, {lambda_idws_at_start} changed to {lambda1}, {lambda2}, {lambda_idws}") - print(f"namd.py: extract_u_nk: Error: fep_file = {fep_file}; has_idws = {has_idws}") + logger.error("namd.py: extract_u_nk: Error: Lambdas changed unexpectedly while processing", fep_file) + logger.error(f"namd.py: extract_u_nk: Error: l1, l2, lidws: {lambda1_at_start}, {lambda2_at_start}, {lambda_idws_at_start} changed to {lambda1}, {lambda2}, {lambda_idws}") + logger.error(f"namd.py: extract_u_nk: Error: fep_file = {fep_file}; has_idws = {has_idws}") return None # As we are at the end of a window, convert last window's work and times values to np arrays @@ -203,7 +207,7 @@ def extract_u_nk(fep_files, T): parsing = True if len(win_de) != 0 or len(win_de_back) != 0: - print('Warning: trailing data without footer line (\"#Free energy...\"). Interrupted run?') + logger.warning('Trailing data without footer line (\"#Free energy...\"). Interrupted run?') if lambda2 in (0.0, 1.0): # this excludes the IDWS case where a dataframe already exists for both endpoints From da5d5b5453ab2f4798dd03c95faf183ec6813a7c Mon Sep 17 00:00:00 2001 From: Tom Joseph Date: Thu, 12 Aug 2021 15:17:23 -0400 Subject: [PATCH 7/9] Tidy up error/warning messages --- src/alchemlyb/parsing/namd.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/alchemlyb/parsing/namd.py b/src/alchemlyb/parsing/namd.py index a801b3d5..ccaa6e40 100644 --- a/src/alchemlyb/parsing/namd.py +++ b/src/alchemlyb/parsing/namd.py @@ -58,7 +58,7 @@ def _get_lambdas(fep_files): # Make sure the lambda_idws values are consistent if lambda_idws is not None: if lambda1 in lambda_bwd_map and lambda_bwd_map[lambda1] != lambda_idws: - logger.error(f'namd.py: get_lambdas: Error: bwd: lambda1 {lambda1} has lambda_idws {lambda_bwd_map[lambda1]} but it should be {lambda_idws}') + logger.error(f'bwd: lambda1 {lambda1} has lambda_idws {lambda_bwd_map[lambda1]} but it should be {lambda_idws}') return None lambda_bwd_map[lambda1] = lambda_idws @@ -150,9 +150,9 @@ def extract_u_nk(fep_files, T): # fails. This can happen if fepouts where one window spans multiple fepouts are processed out of order if lambda1_at_start is not None \ and (lambda1, lambda2, lambda_idws) != (lambda1_at_start, lambda2_at_start, lambda_idws_at_start): - logger.error("namd.py: extract_u_nk: Error: Lambdas changed unexpectedly while processing", fep_file) - logger.error(f"namd.py: extract_u_nk: Error: l1, l2, lidws: {lambda1_at_start}, {lambda2_at_start}, {lambda_idws_at_start} changed to {lambda1}, {lambda2}, {lambda_idws}") - logger.error(f"namd.py: extract_u_nk: Error: fep_file = {fep_file}; has_idws = {has_idws}") + logger.error("Lambdas changed unexpectedly while processing", fep_file) + logger.error(f"Error: l1, l2, lidws: {lambda1_at_start}, {lambda2_at_start}, {lambda_idws_at_start} changed to {lambda1}, {lambda2}, {lambda_idws}") + logger.error(f"Error: fep_file = {fep_file}; has_idws = {has_idws}") return None # As we are at the end of a window, convert last window's work and times values to np arrays From 1376eae1c2b7402e30a76769685df79a488ea963 Mon Sep 17 00:00:00 2001 From: Tom Joseph Date: Thu, 12 Aug 2021 15:24:24 -0400 Subject: [PATCH 8/9] Raise exceptions instead of returning a magic None --- src/alchemlyb/parsing/namd.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/alchemlyb/parsing/namd.py b/src/alchemlyb/parsing/namd.py index ccaa6e40..64f5d3b5 100644 --- a/src/alchemlyb/parsing/namd.py +++ b/src/alchemlyb/parsing/namd.py @@ -51,7 +51,7 @@ def _get_lambdas(fep_files): # Make sure the lambda2 values are consistent if lambda1 in lambda_fwd_map and lambda_fwd_map[lambda1] != lambda2: logger.error(f'fwd: lambda1 {lambda1} has lambda2 {lambda_fwd_map[lambda1]} but it should be {lambda2}') - return None + raise ValueError('Inconsistent lambda values') lambda_fwd_map[lambda1] = lambda2 @@ -59,7 +59,7 @@ def _get_lambdas(fep_files): if lambda_idws is not None: if lambda1 in lambda_bwd_map and lambda_bwd_map[lambda1] != lambda_idws: logger.error(f'bwd: lambda1 {lambda1} has lambda_idws {lambda_bwd_map[lambda1]} but it should be {lambda_idws}') - return None + raise ValueError('Inconsistent lambda values') lambda_bwd_map[lambda1] = lambda_idws all_lambdas = set() @@ -146,14 +146,13 @@ def extract_u_nk(fep_files, T): else: lambda_idws = None - # If the lambdas are not what we thought they would be, return None, ensuring the calculation + # If the lambdas are not what we thought they would be, raise an exception to ensure the calculation # fails. This can happen if fepouts where one window spans multiple fepouts are processed out of order if lambda1_at_start is not None \ and (lambda1, lambda2, lambda_idws) != (lambda1_at_start, lambda2_at_start, lambda_idws_at_start): logger.error("Lambdas changed unexpectedly while processing", fep_file) - logger.error(f"Error: l1, l2, lidws: {lambda1_at_start}, {lambda2_at_start}, {lambda_idws_at_start} changed to {lambda1}, {lambda2}, {lambda_idws}") - logger.error(f"Error: fep_file = {fep_file}; has_idws = {has_idws}") - return None + logger.error(f"l1, l2, lidws: {lambda1_at_start}, {lambda2_at_start}, {lambda_idws_at_start} changed to {lambda1}, {lambda2}, {lambda_idws}") + raise ValueError("Inconsistent lambda values") # As we are at the end of a window, convert last window's work and times values to np arrays # (with energy unit kT since they were kcal/mol in the fepouts) From cbf06989e789c865ce8424af388dbaa66175f91c Mon Sep 17 00:00:00 2001 From: Tom Joseph Date: Fri, 27 Aug 2021 11:28:39 -0400 Subject: [PATCH 9/9] Add test for restarted NAMD FEP run --- src/alchemlyb/tests/parsing/test_namd.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/alchemlyb/tests/parsing/test_namd.py b/src/alchemlyb/tests/parsing/test_namd.py index b41766c1..4fe784aa 100644 --- a/src/alchemlyb/tests/parsing/test_namd.py +++ b/src/alchemlyb/tests/parsing/test_namd.py @@ -6,6 +6,7 @@ from alchemlyb.parsing.namd import extract_u_nk from alchemtest.namd import load_tyr2ala from alchemtest.namd import load_idws +from alchemtest.namd import load_restarted @pytest.fixture(scope="module") @@ -34,3 +35,14 @@ def test_u_nk_idws(): assert u_nk.index.names == ['time', 'fep-lambda'] assert u_nk.shape == (29252, 11) + +def test_u_nk_restarted(): + """Test that u_nk has the correct form when extracted from an IDWS + FEP run that includes a termination and restart. + """ + + filenames = load_restarted()['data']['both'] + u_nk = extract_u_nk(filenames, T=300) + + assert u_nk.index.names == ['time', 'fep-lambda'] + assert u_nk.shape == (30061, 11)