From 9bae718c02fddfe43f9f65ae43028de6bc803132 Mon Sep 17 00:00:00 2001 From: daniel Date: Mon, 7 Sep 2020 16:26:43 +0200 Subject: [PATCH] fix wham2d sampling regions with inf energy --- adaptiveumbrella/__init__.py | 2 +- adaptiveumbrella/runner.py | 37 +++++++++++++++++++++--------------- adaptiveumbrella/wham2d.py | 16 ++++++++-------- 3 files changed, 31 insertions(+), 24 deletions(-) diff --git a/adaptiveumbrella/__init__.py b/adaptiveumbrella/__init__.py index 71da874..484777b 100644 --- a/adaptiveumbrella/__init__.py +++ b/adaptiveumbrella/__init__.py @@ -3,4 +3,4 @@ from .runner import UmbrellaRunner from .wham2d import WHAM2DRunner __all__ = ['AdaptiveUmbrella'] -__version__ = "0.3.12" +__version__ = "0.3.13" diff --git a/adaptiveumbrella/runner.py b/adaptiveumbrella/runner.py index 47b0b09..707e35a 100755 --- a/adaptiveumbrella/runner.py +++ b/adaptiveumbrella/runner.py @@ -85,9 +85,7 @@ class UmbrellaRunner(): def _get_root_frames(self, pmf, frames, E_max): """ returns the index of all positions in the pmf where the energy is - smaller E_max""" - - # select positions of the pmf where E <= E_max and that have already been sampled (frames > 0) + smaller E_max and that have already been sampled""" selection = np.where((pmf <= E_max) & (frames > 0)) zipped = list(zip(*selection)) @@ -132,32 +130,41 @@ class UmbrellaRunner(): that have not an assigned energy yet, as well as their corresponding root frame in the format {new_frame1: root_frame1, new_frame2: root_frame2} """ - # find all neighboring frames and create a dict that associates them to the root frame with lowest energy + # find all neighboring frames and create a dict + # that associates them to the root frame with lowest energy + # new_frame -> corresponding root_frame new_frames = {} - for frame in root_frames: - neighbors = self._generate_neighbor_list(frame) + for root_frame in root_frames: - # remove neighbors if they are not valid (i.e not part of the pmf) + # all frames surrounding this root + neighbors = self._generate_neighbor_list(root_frame) + + # remove invalid (i.e not part of the pmf) neighbors = [n for n in neighbors if self.is_valid_frame(n)] - # for each neighbor, check if its already in the list and compare root frame energy - for n in neighbors: + # for every neighbor, we try + # to get the root frame energy + # see if its already with an associated root in new_frames + # if yes, we swap the root frame if the new one has lower energy + # if not, KeyError, so the current root is our best candidate so far + for neighbor in neighbors: try: - root_energy = pmf[frame] - old_root = new_frames[n] + root_energy = pmf[root_frame] + old_root = new_frames[neighbor] old_root_energy = pmf[old_root] if root_energy < old_root_energy: - new_frames[n] = frame + new_frames[neighbor] = root_frame except KeyError: - new_frames[n] = frame + new_frames[neighbor] = root_frame - - # remove already sampled frames (frames > 0) + # remove already sampled frames (frames with energy > 0) new_frames_list = list(new_frames.keys()) for idx in range(len(new_frames_list)): new_frame = new_frames_list[idx] if frames[new_frame] > 0: del(new_frames[new_frame]) + + # also remove all frames where the root has no energy return new_frames diff --git a/adaptiveumbrella/wham2d.py b/adaptiveumbrella/wham2d.py index 045a8ed..247379d 100644 --- a/adaptiveumbrella/wham2d.py +++ b/adaptiveumbrella/wham2d.py @@ -101,6 +101,8 @@ class WHAM2DRunner(UmbrellaRunner): def load_wham_pmf(self, wham_file): """ Load the new pmf into a pandas dataframe """ df = pd.read_csv(wham_file, delim_whitespace=True, names=['x', 'y', 'e', 'pro'], skiprows=1, index_col=None) + # if e is inf, that means this region is unsampled. we cannot really + # tell its energy and set it to NaN df = df.replace([np.inf, -np.inf], np.nan).dropna(subset=['e'], how='all') return df @@ -112,14 +114,12 @@ class WHAM2DRunner(UmbrellaRunner): lambdax, lambday = self._get_lambdas_for_index((x, y)) reduced = wham_pmf[ (abs(wham_pmf.x-lambdax) < self.cvs[0][2]) & (abs(wham_pmf.y-lambday) < self.cvs[1][2]) ] if len(reduced) == 0: - continue - - reduced['dist'] = reduced.apply(lambda row: np.linalg.norm((lambdax-row['x'], lambday-row['y'])), axis=1) - min = reduced[reduced.dist == reduced.dist.min()] - min_row = reduced[(reduced.x == min.x.iloc[0]) & (reduced.y == min.y.iloc[0])] - self.pmf[x, y] = min_row.e.iloc[0] - - + self.pmf[x,y] = np.inf + else: + reduced['dist'] = reduced.apply(lambda row: np.linalg.norm((lambdax-row['x'], lambday-row['y'])), axis=1) + min = reduced[reduced.dist == reduced.dist.min()] + min_row = reduced[(reduced.x == min.x.iloc[0]) & (reduced.y == min.y.iloc[0])] + self.pmf[x, y] = min_row.e.iloc[0] def calculate_new_pmf(self): metafile_path = self.create_metadata_file()