diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e31e431 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +venv +.idea +plot + diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 index 0000000..6ca207e --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,122 @@ +Creative Commons Legal Code + +CC0 1.0 Universal + + CREATIVE COMMONS CORPORATION IS NOT A LAW FIRM AND DOES NOT PROVIDE + LEGAL SERVICES. DISTRIBUTION OF THIS DOCUMENT DOES NOT CREATE AN + ATTORNEY-CLIENT RELATIONSHIP. CREATIVE COMMONS PROVIDES THIS + INFORMATION ON AN "AS-IS" BASIS. CREATIVE COMMONS MAKES NO WARRANTIES + REGARDING THE USE OF THIS DOCUMENT OR THE INFORMATION OR WORKS + PROVIDED HEREUNDER, AND DISCLAIMS LIABILITY FOR DAMAGES RESULTING FROM + THE USE OF THIS DOCUMENT OR THE INFORMATION OR WORKS PROVIDED + HEREUNDER. + +Statement of Purpose + +The laws of most jurisdictions throughout the world automatically confer +exclusive Copyright and Related Rights (defined below) upon the creator +and subsequent owner(s) (each and all, an "owner") of an original work of +authorship and/or a database (each, a "Work"). + +Certain owners wish to permanently relinquish those rights to a Work for +the purpose of contributing to a commons of creative, cultural and +scientific works ("Commons") that the public can reliably and without fear +of later claims of infringement build upon, modify, incorporate in other +works, reuse and redistribute as freely as possible in any form whatsoever +and for any purposes, including without limitation commercial purposes. +These owners may contribute to the Commons to promote the ideal of a free +culture and the further production of creative, cultural and scientific +works, or to gain reputation or greater distribution for their Work in +part through the use and efforts of others. + +For these and/or other purposes and motivations, and without any +expectation of additional consideration or compensation, the person +associating CC0 with a Work (the "Affirmer"), to the extent that he or she +is an owner of Copyright and Related Rights in the Work, voluntarily +elects to apply CC0 to the Work and publicly distribute the Work under its +terms, with knowledge of his or her Copyright and Related Rights in the +Work and the meaning and intended legal effect of CC0 on those rights. + +1. Copyright and Related Rights. A Work made available under CC0 may be +protected by copyright and related or neighboring rights ("Copyright and +Related Rights"). Copyright and Related Rights include, but are not +limited to, the following: + + i. the right to reproduce, adapt, distribute, perform, display, + communicate, and translate a Work; + ii. moral rights retained by the original author(s) and/or performer(s); +iii. publicity and privacy rights pertaining to a person's image or + likeness depicted in a Work; + iv. rights protecting against unfair competition in regards to a Work, + subject to the limitations in paragraph 4(a), below; + v. rights protecting the extraction, dissemination, use and reuse of data + in a Work; + vi. database rights (such as those arising under Directive 96/9/EC of the + European Parliament and of the Council of 11 March 1996 on the legal + protection of databases, and under any national implementation + thereof, including any amended or successor version of such + directive); and +vii. other similar, equivalent or corresponding rights throughout the + world based on applicable law or treaty, and any national + implementations thereof. + +2. Waiver. To the greatest extent permitted by, but not in contravention +of, applicable law, Affirmer hereby overtly, fully, permanently, +irrevocably and unconditionally waives, abandons, and surrenders all of +Affirmer's Copyright and Related Rights and associated claims and causes +of action, whether now known or unknown (including existing as well as +future claims and causes of action), in the Work (i) in all territories +worldwide, (ii) for the maximum duration provided by applicable law or +treaty (including future time extensions), (iii) in any current or future +medium and for any number of copies, and (iv) for any purpose whatsoever, +including without limitation commercial, advertising or promotional +purposes (the "Waiver"). Affirmer makes the Waiver for the benefit of each +member of the public at large and to the detriment of Affirmer's heirs and +successors, fully intending that such Waiver shall not be subject to +revocation, rescission, cancellation, termination, or any other legal or +equitable action to disrupt the quiet enjoyment of the Work by the public +as contemplated by Affirmer's express Statement of Purpose. + +3. Public License Fallback. Should any part of the Waiver for any reason +be judged legally invalid or ineffective under applicable law, then the +Waiver shall be preserved to the maximum extent permitted taking into +account Affirmer's express Statement of Purpose. In addition, to the +extent the Waiver is so judged Affirmer hereby grants to each affected +person a royalty-free, non transferable, non sublicensable, non exclusive, +irrevocable and unconditional license to exercise Affirmer's Copyright and +Related Rights in the Work (i) in all territories worldwide, (ii) for the +maximum duration provided by applicable law or treaty (including future +time extensions), (iii) in any current or future medium and for any number +of copies, and (iv) for any purpose whatsoever, including without +limitation commercial, advertising or promotional purposes (the +"License"). The License shall be deemed effective as of the date CC0 was +applied by Affirmer to the Work. Should any part of the License for any +reason be judged legally invalid or ineffective under applicable law, such +partial invalidity or ineffectiveness shall not invalidate the remainder +of the License, and in such case Affirmer hereby affirms that he or she +will not (i) exercise any of his or her remaining Copyright and Related +Rights in the Work or (ii) assert any associated claims and causes of +action with respect to the Work, in either case contrary to Affirmer's +express Statement of Purpose. + +4. Limitations and Disclaimers. + + a. No trademark or patent rights held by Affirmer are waived, abandoned, + surrendered, licensed or otherwise affected by this document. + b. Affirmer offers the Work as-is and makes no representations or + warranties of any kind concerning the Work, express, implied, + statutory or otherwise, including without limitation warranties of + title, merchantability, fitness for a particular purpose, non + infringement, or the absence of latent or other defects, accuracy, or + the present or absence of errors, whether or not discoverable, all to + the greatest extent permissible under applicable law. + c. Affirmer disclaims responsibility for clearing rights of other persons + that may apply to the Work or any use thereof, including without + limitation any person's Copyright and Related Rights in the Work. + Further, Affirmer disclaims responsibility for obtaining any necessary + consents, permissions or other rights required for any use of the + Work. + d. Affirmer understands and acknowledges that Creative Commons is not a + party to this document and has no duty or obligation with respect to + this CC0 or use of the Work. + diff --git a/README.md b/README.md new file mode 100644 index 0000000..a504f2e --- /dev/null +++ b/README.md @@ -0,0 +1,15 @@ +# Python module for adaptive umbrella sampling + +This module can be used to perform adaptive umbrella sampling of a multi-dimensional potential of mean force. The +algorithm involves:: + +1) calculate the free energy landscape +2) Among existing windows, select windows with E < E_max +3) For each selected window, generate 3^N-1 neighbor windows +4) Sample new windows, then go to 1) or stop if no new windows can be found + +For more details about the algorithm: See + +Self-Learning Adaptive Umbrella Sampling Method for the Determination of Free Energy Landscapes in Multiple Dimensions (Wojtas-Niziurski†, Meng, Roux, Bernèche, 2013) +[(https://doi.org/10.1021/ct300978b)](https://doi.org/10.1021/ct300978b) + diff --git a/adaptiveumbrella/__init__.py b/adaptiveumbrella/__init__.py new file mode 100644 index 0000000..a4a144c --- /dev/null +++ b/adaptiveumbrella/__init__.py @@ -0,0 +1,5 @@ +from __future__ import absolute_import +from .runner import UmbrellaRunner + +__all__ = ['AdaptiveUmbrella'] +__version__ = "0.1.0" \ No newline at end of file diff --git a/adaptive_umbrella.py b/adaptiveumbrella/runner.py similarity index 81% rename from adaptive_umbrella.py rename to adaptiveumbrella/runner.py index 188ae06..ee84a25 100755 --- a/adaptive_umbrella.py +++ b/adaptiveumbrella/runner.py @@ -216,99 +216,6 @@ class UmbrellaRunner(): """ This can be implemented to hook into the simulation cycle before the simulation runs """ pass -class WHAM2DRunner(UmbrellaRunner): - """ Umbrella runner implementation that uses wham-2d to perform - the pmf calculation. - - Attributes: - WHAM_EXEC: path to wham executeable - - """ - - def __init__(self): - UmbrellaRunner.__init__(self) - self.WHAM_EXEC = 'wham-2d' - - def calculate_new_pmf(self): - import os - from shutil import copyfile - - simulation_dir = "simulations" - print("Collecting sampling data from simulations folder") - - # collect COLVARs - wham_dir = "WHAM/" - if not os.path.exists(wham_dir): - os.makedirs(wham_dir) - - for folder in os.listdir(simulation_dir): - src = os.path.join(simulation_dir, folder, "COLVAR") - dst = os.path.join(wham_dir, folder + ".xvg") - copyfile(src, dst) - - # create metadata file - metadata_file = os.path.join(wham_dir, "{}_metadata.dat".format(self.num_iterations)) - fc_x = 100 - fc_y = 100 - with open(metadata_file, 'w') as out: - for f in os.listdir(simulation_dir): - prefix, x, y = f.split("_") - out.write("WHAM/{}.xvg {} {} {} {}\n".format(f, x, y, fc_x, fc_y)) - - # run WHAM2d - print("Running WHAM-2d") - wham_output = os.path.join(wham_dir, "{}_freeenergy.dat".format(self.num_iterations)) - periodicity_x = "pi" - 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] - - 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, - px=periodicity_x, - min_x=min_x, - max_x=max_x, - frames_x=frames_x, - py=periodicity_y, - min_y=min_y, - max_y=max_y, - frames_y=frames_y, - tol=tolerance, - metafile=metadata_file, - outfile=wham_output - ) - print(cmd) - os.system(cmd) - - # read wham to new pmf - return self.read_pmf(wham_output) - - def read_pmf(self, pmf_path): - import pandas as pd - print("Update PMF from WHAM") - df = pd.read_csv(pmf_path, delim_whitespace=True, names=['x', 'y', 'e', 'pro'], skiprows=1, - index_col=None) - df = df.replace([np.inf, -np.inf], np.nan).dropna(subset=['e'], how='all') - new_pmf = deepcopy(self.pmf) - for x in range(new_pmf.shape[0]): - for y in range(new_pmf.shape[1]): - lambdax, lambday = self._get_lambdas_for_index((x, y)) - x_selection = (df.x - lambdax).abs() < 0.01 - y_selection = (df.y - lambday).abs() < 0.01 - selected_energies = df[(x_selection) & (y_selection)].e - if len(selected_energies) == 0: - new_pmf[x, y] = -1 - else: - new_pmf[x, y] = selected_energies.iloc[0] - - return new_pmf - - - class UmbrellaRunnerTest(unittest.TestCase): def test_init_pmf_3d(self): diff --git a/example.py b/example.py deleted file mode 100644 index 945f071..0000000 --- a/example.py +++ /dev/null @@ -1,59 +0,0 @@ -from adaptive_umbrella import WHAM2DRunner - -from copy import deepcopy -import os -import numpy as np -import matplotlib.pyplot as plt - - -class MyUmbrellaRunner(WHAM2DRunner): - def after_run_hook(self): - filename = "pmf_{}.pdf".format(self.num_iterations) - print("Writing new pmf to {}".format(filename)) - pmf_to_plot = deepcopy(self.pmf.T) - pmf_to_plot[pmf_to_plot < 0] = None - plt.figure() - plt.imshow(pmf_to_plot, origin="bottom", cmap='jet') - 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 ] - - plt.xticks(tick_positions, tick_labels) - plt.yticks(tick_positions, tick_labels) - - - cb = plt.colorbar(pad=0.1) - cb.set_label("kJ/mol") - plt.savefig(filename) - os.system("cp {} {}".format(filename, "pmf_current.pdf")) - - - def simulate_frames(self, lambdas, frames): - print("{} new simulations:".format(len(lambdas))) - counter = 0 - - threads = [] - for f in lambdas: - counter += 1 - if os.path.exists("sim/sim_{}_{}/COLVAR".format(*f)): - print("{}) Skipping lambdas={}/{}: COLVAR exists".format(counter, *f)) - continue - - print("{}) Simulate lambda1={}, lambda2={}".format(counter, *f)) - command = "bash sim.sh {} {} 2>&1 > run.log".format(*f) - # print("Running {}".format(command)) - os.system(command) - -runner = MyUmbrellaRunner() -runner.WHAM_EXEC = "/opt/wham/wham-2d/wham-2d" -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 = 100 - -runner.run() diff --git a/examples/example.py b/examples/example.py new file mode 100644 index 0000000..0ad6329 --- /dev/null +++ b/examples/example.py @@ -0,0 +1,152 @@ +import os +from copy import deepcopy + +import matplotlib.pyplot as plt +import numpy as np + +from adaptiveumbrella.runner import UmbrellaRunner + + +class WHAM2DRunner(UmbrellaRunner): + """ Umbrella runner implementation that uses wham-2d to perform + the pmf calculation. + + Attributes: + WHAM_EXEC: path to wham executeable + + """ + + def __init__(self): + UmbrellaRunner.__init__(self) + self.WHAM_EXEC = 'wham-2d' + + def calculate_new_pmf(self): + import os + from shutil import copyfile + + simulation_dir = "simulations" + print("Collecting sampling data from simulations folder") + + # collect COLVARs + wham_dir = "WHAM/" + if not os.path.exists(wham_dir): + os.makedirs(wham_dir) + + for folder in os.listdir(simulation_dir): + src = os.path.join(simulation_dir, folder, "COLVAR") + dst = os.path.join(wham_dir, folder + ".xvg") + copyfile(src, dst) + + # create metadata file + metadata_file = os.path.join(wham_dir, "{}_metadata.dat".format(self.num_iterations)) + fc_x = 100 + fc_y = 100 + with open(metadata_file, 'w') as out: + for f in os.listdir(simulation_dir): + prefix, x, y = f.split("_") + out.write("WHAM/{}.xvg {} {} {} {}\n".format(f, x, y, fc_x, fc_y)) + + # run WHAM2d + print("Running WHAM-2d") + wham_output = os.path.join(wham_dir, "{}_freeenergy.dat".format(self.num_iterations)) + periodicity_x = "pi" + 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] + + 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, + px=periodicity_x, + min_x=min_x, + max_x=max_x, + frames_x=frames_x, + py=periodicity_y, + min_y=min_y, + max_y=max_y, + frames_y=frames_y, + tol=tolerance, + metafile=metadata_file, + outfile=wham_output + ) + print(cmd) + os.system(cmd) + + # read wham to new pmf + return self.read_pmf(wham_output) + + def read_pmf(self, pmf_path): + import pandas as pd + print("Update PMF from WHAM") + df = pd.read_csv(pmf_path, delim_whitespace=True, names=['x', 'y', 'e', 'pro'], skiprows=1, + index_col=None) + df = df.replace([np.inf, -np.inf], np.nan).dropna(subset=['e'], how='all') + new_pmf = deepcopy(self.pmf) + for x in range(new_pmf.shape[0]): + for y in range(new_pmf.shape[1]): + lambdax, lambday = self._get_lambdas_for_index((x, y)) + x_selection = (df.x - lambdax).abs() < 0.01 + y_selection = (df.y - lambday).abs() < 0.01 + selected_energies = df[(x_selection) & (y_selection)].e + if len(selected_energies) == 0: + new_pmf[x, y] = -1 + else: + new_pmf[x, y] = selected_energies.iloc[0] + + return new_pmf + + +class MyUmbrellaRunner(WHAM2DRunner): + def after_run_hook(self): + filename = "pmf_{}.pdf".format(self.num_iterations) + print("Writing new pmf to {}".format(filename)) + pmf_to_plot = deepcopy(self.pmf.T) + pmf_to_plot[pmf_to_plot < 0] = None + plt.figure() + plt.imshow(pmf_to_plot, origin="bottom", cmap='jet') + 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 ] + + plt.xticks(tick_positions, tick_labels) + plt.yticks(tick_positions, tick_labels) + + + cb = plt.colorbar(pad=0.1) + cb.set_label("kJ/mol") + plt.savefig(filename) + os.system("cp {} {}".format(filename, "pmf_current.pdf")) + + + def simulate_frames(self, lambdas, frames): + print("{} new simulations:".format(len(lambdas))) + counter = 0 + + threads = [] + for f in lambdas: + counter += 1 + if os.path.exists("sim/sim_{}_{}/COLVAR".format(*f)): + print("{}) Skipping lambdas={}/{}: COLVAR exists".format(counter, *f)) + continue + + print("{}) Simulate lambda1={}, lambda2={}".format(counter, *f)) + command = "bash sim.sh {} {} 2>&1 > run.log".format(*f) + # print("Running {}".format(command)) + os.system(command) + +runner = MyUmbrellaRunner() +runner.WHAM_EXEC = "/opt/wham/wham-2d/wham-2d" +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 = 100 + +runner.run() diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..cd9df90 --- /dev/null +++ b/setup.py @@ -0,0 +1,32 @@ +#!/usr/bin/env python3 +try: + from setuptools import setup +except ImportError: + from distutils.core import setup + +LONG_DESCRIPTION = """`adaptiveumbrella` is a python implementation of the +[self-learning adaptive umbrella sampling algorithm (Wojtas-Niziurski†, Meng, Roux, Bernèche, 2013)] +(https://pubs.acs.org/doi/abs/10.1021/ct300978b) technique. It allows the calculation of a multidimensional potential +of mean force while automatically exploring the phase space. +""" + +# Parse the version from the fiona module. +with open('adaptiveumbrella/__init__.py') as f: + for line in f: + if line.find("__version__") >= 0: + version = line.split("=")[1].strip() + version = version.strip('"') + version = version.strip("'") + break + +setup( + name='adaptiveumbrella', + version=version, + description='Adaptive umbrella sampling in Python', + license='CC0', + author='Daniel Bauer', + author_email='bauer@cbs.tu-darmstadt.de', + url='https://github.com/danijoo/adaptiveumbrella', + long_description=LONG_DESCRIPTION, + packages=['adaptiveumbrella'], +install_requires=['numpy']) \ No newline at end of file