diff --git a/adaptiveumbrella/runner.py b/adaptiveumbrella/runner.py index 4c2321c..bd4e892 100755 --- a/adaptiveumbrella/runner.py +++ b/adaptiveumbrella/runner.py @@ -150,11 +150,11 @@ class UmbrellaRunner(): new_frames = self._get_new_frames(self.pmf, self.sample_list, root_frames) self.num_iterations = 0 - self.E = self.E_min - + # outer main loop: increase E and calculate PMF until E > E_max while True: - + self.E = self.E_min + # stop if max iterations is reached self.num_iterations += 1 if(self.max_iterations > 0 and self.num_iterations > self.max_iterations): @@ -165,14 +165,14 @@ class UmbrellaRunner(): if len(np.where(self.pmf < 0)) == 0: print("Every window of the PMF appears to be sampled.") return - + print("~~~~~~~~~~~~~~~ Iteration {}/{} ~~~~~~~~~~~~~~~~".format(self.num_iterations, self.max_iterations)) lambdas = dict([(self._get_lambdas_for_index(x), self._get_lambdas_for_index(y)) for x,y in new_frames.items()]) self.pre_run_hook() print("Running simulations") self.simulate_frames(lambdas, new_frames) - + print("Calculating new PMF") self.pmf = self.calculate_new_pmf() @@ -180,11 +180,11 @@ class UmbrellaRunner(): for new_frame in new_frames.keys(): self.sample_list[new_frame] = self.num_iterations self.after_run_hook() - + while self.E <= self.E_max: root_frames = self._get_root_frames(self.pmf, self.sample_list, self.E) new_frames = self._get_new_frames(self.pmf, self.sample_list, root_frames) - + if len(new_frames) == 0: self.E += self.E_incr print("Max energy increased to {} (max={})".format(self.E, self.E_max)) diff --git a/examples/example.py b/examples/example.py index 7fbee68..70a786b 100644 --- a/examples/example.py +++ b/examples/example.py @@ -55,10 +55,14 @@ class WHAM2DRunner(UmbrellaRunner): periodicity_y = "pi" tolerance = 0.1 frames_x, frames_y = 1002, 1002 - min_x = self.cvs[0][0] - max_x = self.cvs[0][1] - min_y = self.cvs[1][0] - max_y = self.cvs[1][1] + # min_x = self.cvs[0][0] + # max_x = self.cvs[0][1] + # min_y = self.cvs[1][0] + # max_y = self.cvs[1][1] + min_x = -3 + max_x = 3 + min_y = -3 + max_y = 3 cmd = "{exec} Px={px} {min_x} {max_x} {frames_x} Py={py} {min_y} {max_y} {frames_y} {tol} 298 0 {metafile} {outfile} 0".format( exec=self.WHAM_EXEC, @@ -102,31 +106,42 @@ class WHAM2DRunner(UmbrellaRunner): class MyUmbrellaRunner(WHAM2DRunner): + + cum_frames = [0] + def after_run_hook(self): filename = "tmp/pmf_{}.pdf".format(self.num_iterations) print("Writing new pmf to {}".format(filename)) - pmf_to_plot = deepcopy(self.pmf.T) + pmf_to_plot = deepcopy(self.pmf) pmf_to_plot[pmf_to_plot < 0] = None - frames_to_plot = deepcopy(self.sample_list.T) - frames_to_plot[frames_to_plot == 0] = None + pmf_to_plot[self.sample_list == 0] = None + pmf_to_plot = pmf_to_plot.T - fig, (ax0, ax1) = plt.subplots(ncols=2, sharey=True) - im = ax0.imshow(pmf_to_plot, origin="bottom", cmap='jet') - cb = fig.colorbar(im, ax=ax0) + self.cum_frames.append(len(self.sample_list[self.sample_list > 0])) + + fig, (ax0, ax1) = plt.subplots(ncols=2) + im = ax0.imshow(pmf_to_plot, origin='lower', cmap='jet') + cb = fig.colorbar(im, ax=ax0, orientation='horizontal', pad=0.15) cb.set_label("kJ/mol") - im2 = ax1.imshow(frames_to_plot, origin="bottom") - cb2 = fig.colorbar(im2, ax=ax1, ticks=np.arange(0, self.max_iterations)) - cb.set_label("cycle") - # ticks = [(x,x) for x in [-3, -2, -1, 0, 1, 2, 3]] - # tick_positions = [ self._get_index_for_lambdas(x)[0] for x in ticks ] - # tick_labels = [ str(x[0]) for x in ticks ] - # ax1.set_xticks(tick_positions, tick_labels) - # ax1.set_yticks(tick_positions, tick_labels) - # ax2.set_xticks(tick_positions, tick_labels) - # ax2.set_yticks(tick_positions, tick_labels) + ax1.plot(self.cum_frames, linewidth=0.5, marker="o", color='black') + ax1.set_xlabel("Cycles") + ax1.set_ylabel("Number of umbrella Windows") - plt.savefig(filename) + ticks = [(x,x) for x in [-3, -2, -1, 0, 1, 2, 3]] + tick_positions = [ self._get_index_for_lambdas(x)[0] for x in ticks ] + tick_labels = [ str(x[0]) for x in ticks ] + ax0.set_yticks(tick_positions) + ax0.set_yticklabels(tick_labels) + ax0.set_xticks(tick_positions) + ax0.set_xticklabels(tick_labels) + ax0.set_ylabel("$\phi$") + ax0.set_xlabel("$\psi$") + + fig.subplots_adjust(wspace=.5) + + + plt.savefig(filename, bbox_inches='tight', dpi=200) os.system("cp {} {}".format(filename, "tmp/pmf_current.pdf")) @@ -151,10 +166,10 @@ runner.cvs = np.array([ (-3, 3, 0.2), (-3, 3, 0.2), ]) -runner.cvs_init = (1.4, -1.4) -runner.E_min = 10 -runner.E_max = 100 -runner.E_incr = 10 -runner.max_iterations = 1 +runner.cvs_init = (-1.4, 1.4) +runner.E_min = 2 +runner.E_max = 50 +runner.E_incr = 5 +runner.max_iterations = 30 runner.run()