diff --git a/benchmark/evaluation/__init__.py b/benchmark/evaluation/__init__.py index edb7fee..cf1f97c 100644 --- a/benchmark/evaluation/__init__.py +++ b/benchmark/evaluation/__init__.py @@ -1,7 +1,7 @@ -from analysis import estimate_nonequilibrium_free_energy, compute_free_energy -from entropy import estimate_entropy, estimate_marginal_entropies +from .analysis import estimate_nonequilibrium_free_energy, compute_free_energy +from .entropy import estimate_entropy, estimate_marginal_entropies __all__ = ["estimate_nonequilibrium_free_energy", "compute_free_energy", "estimate_entropy", "estimate_marginal_entropies" - ] \ No newline at end of file + ] diff --git a/benchmark/evaluation/analysis.py b/benchmark/evaluation/analysis.py index a75d433..209a622 100644 --- a/benchmark/evaluation/analysis.py +++ b/benchmark/evaluation/analysis.py @@ -1,6 +1,6 @@ import numpy as np from pickle import load -from entropy import estimate_marginal_entropies, estimate_entropy +from .entropy import estimate_marginal_entropies, estimate_entropy import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt @@ -58,7 +58,7 @@ def compute_free_energy(xv, potential, kinetic_energy, beta): if __name__ == "__main__": print("reading and re-analyzing results") name = "alanine_unconstrained_null_results.pkl" - with open(name, "r") as f: results = load(f) + with open(name, "rb") as f: results = load(f) for scheme in results.keys(): print(scheme) W_shads_F, W_shads_R, DeltaF_neq, sq_uncertainty, W_midpoints = results[scheme] diff --git a/benchmark/experiments/baoab_vs_vvvr.py b/benchmark/experiments/baoab_vs_vvvr.py index bfee845..a63bcf8 100644 --- a/benchmark/experiments/baoab_vs_vvvr.py +++ b/benchmark/experiments/baoab_vs_vvvr.py @@ -6,24 +6,26 @@ import pickle import numpy as np from simtk import unit -from benchmark.plotting import plot_scheme_comparison +#from baoab_vs_aboba_analysis import plot_results from benchmark.evaluation.analysis import estimate_nonequilibrium_free_energy +from benchmark.plotting import plot, savefig if __name__ == "__main__": n_protocol_samples, protocol_length = 1000, 50 - system_name = "alanine_unconstrained" - equilibrium_simulator = benchmark.testsystems.alanine_unconstrained - target_filename = os.path.join(DATA_PATH, "scheme_comparison_{}.pkl".format(system_name)) + system_name = "alanine_constrained" + equilibrium_simulator = benchmark.testsystems.alanine_constrained + target_filename = os.path.join(DATA_PATH, "baoab_vs_aboba_{}.pkl".format(system_name)) - schemes = {"BAO": "V R O", "BABO": "V R V O", - "ABO": "R V O", "ABAO": "R V R O"} - timesteps = np.linspace(0.1, 1.5, 5) + schemes = {"BAOAB": "V R O O R V", "VVVR": "O V R R V O"} + timesteps = np.linspace(0.1, 3.0, 10) noneq_simulators = {} for name, scheme in schemes.items(): for timestep in timesteps: noneq_simulators[(name, timestep)] = NonequilibriumSimulator(equilibrium_simulator, LangevinSplittingIntegrator( - splitting=scheme, timestep=timestep * unit.femtosecond)) + splitting=scheme, + timestep=timestep * unit.femtosecond, + collision_rate=1.0/unit.picoseconds)) # need to catch "Exception: Particle coordinate is nan" results = {} @@ -40,7 +42,7 @@ DeltaF_neq, squared_uncertainty = estimate_nonequilibrium_free_energy(W_shads_F, W_shads_R) print("\t{:.3f} +/- {:.3f}".format(DeltaF_neq, np.sqrt(squared_uncertainty))) - with open(target_filename, "w") as f: + with open(target_filename, "wb") as f: pickle.dump(results, f) - plot_scheme_comparison(target_filename, system_name) \ No newline at end of file + # TODO: Plot stuff diff --git a/benchmark/experiments/potential_energy_distributions.py b/benchmark/experiments/potential_energy_distributions.py index 698f990..2598b12 100644 --- a/benchmark/experiments/potential_energy_distributions.py +++ b/benchmark/experiments/potential_energy_distributions.py @@ -3,7 +3,7 @@ import matplotlib -from code.testsystems import system_params +from benchmark.testsystems import system_params matplotlib.use("Agg") import matplotlib.pyplot as plt @@ -13,7 +13,7 @@ from simtk import unit import numpy as np from tqdm import tqdm -from code.utils import stderr +from benchmark.utilities import stderr histstyle = {"alpha":0.3, "histtype":"stepfilled", "bins":50 @@ -98,4 +98,4 @@ def get_samples(simulation): from pickle import dump - with open("constrained_energies.pkl", "w") as f: dump(results, f) \ No newline at end of file + with open("constrained_energies.pkl", "wb") as f: dump(results, f) diff --git a/benchmark/experiments/quartic_kl_validation.py b/benchmark/experiments/quartic_kl_validation.py index 76afc93..4338fec 100644 --- a/benchmark/experiments/quartic_kl_validation.py +++ b/benchmark/experiments/quartic_kl_validation.py @@ -58,7 +58,7 @@ DeltaF_neq, squared_uncertainty = estimate_nonequilibrium_free_energy(W_shads_F, W_shads_R) print("\t{:.5f} +/- {:.5f}".format(DeltaF_neq, np.sqrt(squared_uncertainty))) - with open(target_filename, "w") as f: + with open(target_filename, "wb") as f: pickle.dump(results, f) - plot_scheme_comparison(target_filename, system_name) \ No newline at end of file + plot_scheme_comparison(target_filename, system_name) diff --git a/benchmark/experiments/quartic_protocol_length_dependence.py b/benchmark/experiments/quartic_protocol_length_dependence.py index e1c41be..c071462 100644 --- a/benchmark/experiments/quartic_protocol_length_dependence.py +++ b/benchmark/experiments/quartic_protocol_length_dependence.py @@ -48,7 +48,7 @@ errors.append(np.sqrt(squared_uncertainty)) print("\t{:.5f} +/- {:.5f}".format(DeltaF_neq, np.sqrt(squared_uncertainty))) - with open(data_filename, "w") as f: + with open(data_filename, "wb") as f: pickle.dump(results, f) plt.figure() diff --git a/benchmark/experiments/replicate_paper_plot.py b/benchmark/experiments/replicate_paper_plot.py index 63401bb..a8a128c 100644 --- a/benchmark/experiments/replicate_paper_plot.py +++ b/benchmark/experiments/replicate_paper_plot.py @@ -5,12 +5,12 @@ from benchmark.experiments.benchmark import get_equilibrium_samples, null_midpoint_operator from benchmark.integrators.integrators import LangevinSplittingIntegrator -from code.testsystems import system_params +from benchmark.testsystems import system_params matplotlib.use("Agg") import matplotlib.pyplot as plt -from code.timestep_search import sweep_over_timesteps +from benchmark.experiments.timestep_search import sweep_over_timesteps if __name__ == "__main__": unconstrained_timesteps_to_try = [0.5, 0.75, 1, 1.25, 1.5] diff --git a/benchmark/experiments/self_averaging.py b/benchmark/experiments/self_averaging.py index 794930f..8573c36 100644 --- a/benchmark/experiments/self_averaging.py +++ b/benchmark/experiments/self_averaging.py @@ -51,5 +51,5 @@ def oscillator_factory(coupling_strength=0.0): results[coupling_strength] = W_shads - with open(target_filename, "w") as f: + with open(target_filename, "wb") as f: pickle.dump(results, f) diff --git a/benchmark/experiments/benchmark.py b/benchmark/experiments/tools.py similarity index 95% rename from benchmark/experiments/benchmark.py rename to benchmark/experiments/tools.py index 2335810..e058564 100644 --- a/benchmark/experiments/benchmark.py +++ b/benchmark/experiments/tools.py @@ -2,8 +2,8 @@ from simtk import unit from simtk.openmm import app -from code.integrators import LangevinSplittingIntegrator -from code.testsystems import system_params +from benchmark.integrators import LangevinSplittingIntegrator +from benchmark.testsystems import system_params W_unit = unit.kilojoule_per_mole from benchmark.evaluation.analysis import estimate_nonequilibrium_free_energy @@ -11,7 +11,9 @@ from openmmtools.integrators import GHMCIntegrator from tqdm import tqdm -from code.utils import plot, get_total_energy, get_summary_string, measure_shadow_work +from benchmark.plotting import plot +from benchmark.utilities import get_total_energy, get_summary_string +from benchmark.tests.measure_shadow_work import measure_shadow_work def randomization_midpoint_operator(simulation, temperature): @@ -93,7 +95,7 @@ def collect_and_save_results(schemes, simulation_factory, equilibrium_samples, M=M, midpoint_operator=midpoint_operator, temperature=temperature) print("\n".join((scheme, get_summary_string(results[scheme], linebreaks=True)))) - with open("{}_results.pkl".format(name), "w") as f: dump(results, f) + with open("{}_results.pkl".format(name), "wb") as f: dump(results, f) return results diff --git a/benchmark/experiments/waterbox_test.py b/benchmark/experiments/waterbox_test.py index 561f597..5371177 100644 --- a/benchmark/experiments/waterbox_test.py +++ b/benchmark/experiments/waterbox_test.py @@ -12,11 +12,11 @@ from simtk.openmm import app from tqdm import tqdm -from benchmark.experiments.benchmark import estimate_nonequilibrium_free_energy -from benchmark.experiments.benchmark import get_equilibrium_samples, \ +from benchmark.experiments.tools import estimate_nonequilibrium_free_energy +from benchmark.experiments.tools import get_equilibrium_samples, \ randomization_midpoint_operator, apply_protocol from benchmark.integrators.integrators import LangevinSplittingIntegrator -from code.testsystems import system_params +from benchmark.testsystems import system_params matplotlib.use("Agg") import matplotlib.pyplot as plt @@ -133,4 +133,4 @@ def plot_deltaF_neq(schemes, Ms, midpoint_operator): plot_deltaF_neq(schemes, Ms, lambda simulation: randomization_midpoint_operator(simulation, temperature)) plt.savefig("../figures/{}_conf_DeltaFs_over_Ms.jpg".format(name)) - plt.close() \ No newline at end of file + plt.close() diff --git a/benchmark/integrators/__init__.py b/benchmark/integrators/__init__.py index 55ca8ab..b15f872 100644 --- a/benchmark/integrators/__init__.py +++ b/benchmark/integrators/__init__.py @@ -1,9 +1,9 @@ -from ghmc import CustomizableGHMC -from langevin import LangevinSplittingIntegrator -from numba_integrators import baoab_factory, vvvr_factory, metropolis_hastings_factory, aboba_factory -from mts_utilities import condense_splitting, generate_sequential_BAOAB_string, generate_all_BAOAB_permutation_strings +from .ghmc import CustomizableGHMC +from .langevin import LangevinSplittingIntegrator +from .numba_integrators import baoab_factory, vvvr_factory, metropolis_hastings_factory, aboba_factory +from .mts_utilities import condense_splitting, generate_sequential_BAOAB_string, generate_all_BAOAB_permutation_strings, generate_solvent_solute_splitting_string __all__ = ["CustomizableGHMC", "LangevinSplittingIntegrator", "baoab_factory", "vvvr_factory", "aboba_factory", "metropolis_hastings_factory", "condense_splitting", "generate_sequential_BAOAB_string", "generate_all_BAOAB_permutation_strings" - ] \ No newline at end of file + ] diff --git a/benchmark/integrators/langevin.py b/benchmark/integrators/langevin.py index f11bf14..24fd6f2 100644 --- a/benchmark/integrators/langevin.py +++ b/benchmark/integrators/langevin.py @@ -122,6 +122,9 @@ def __init__(self, # Define substep functions def R_step(): + self.addConstrainPositions() # TODO: Constrain initial step only? + self.addConstrainVelocities() # TODO: Constrain initial step only? + if measure_shadow_work: self.addComputeGlobal("old_pe", "energy") self.addComputeSum("old_ke", kinetic_energy) @@ -147,6 +150,9 @@ def V_step(fg): Force group to use in this substep. "" means all forces, "0" means force-group 0, etc. """ + self.addConstrainPositions() # TODO: Constrain initial step only? + self.addConstrainVelocities() # TODO: Constrain initial step only? + if measure_shadow_work: self.addComputeSum("old_ke", kinetic_energy) @@ -163,6 +169,9 @@ def V_step(fg): self.addComputeGlobal("shadow_work", "shadow_work + (new_ke - old_ke)") def O_step(): + self.addConstrainPositions() # TODO: Constrain initial step only? + self.addConstrainVelocities() # TODO: Constrain initial step only? + if measure_heat: self.addComputeSum("old_ke", kinetic_energy) diff --git a/benchmark/plotting/__init__.py b/benchmark/plotting/__init__.py index 5278512..6b0e938 100644 --- a/benchmark/plotting/__init__.py +++ b/benchmark/plotting/__init__.py @@ -1,3 +1,3 @@ -from plotting_utilities import generate_figure_filename, plot_scheme_comparison +from .plotting_utilities import savefig, plot, generate_figure_filename, plot_scheme_comparison -__all__ = ["generate_figure_filename", "plot_scheme_comparison"] \ No newline at end of file +__all__ = ["generate_figure_filename", "plot_scheme_comparison"] diff --git a/benchmark/plotting/plotting_utilities.py b/benchmark/plotting/plotting_utilities.py index f6f49d7..a227433 100644 --- a/benchmark/plotting/plotting_utilities.py +++ b/benchmark/plotting/plotting_utilities.py @@ -6,6 +6,7 @@ import pickle import os from benchmark import FIGURE_PATH +from benchmark.utilities import unpack_trajs figure_directory = FIGURE_PATH figure_format = ".jpg" @@ -129,7 +130,7 @@ def plot_scheme_comparison(target_filename, name): and one for the i.e. results[marginal][(name, timestep)] = W_shads_F, W_shads_R """ - with open(target_filename, "r") as f: + with open(target_filename, "rb") as f: results = pickle.load(f) def get_schemes(result_dict): @@ -172,4 +173,4 @@ def plot_curves(results): plt.savefig(generate_figure_filename("scheme_comparison_{}.jpg".format(name)), dpi=300) plt.close() - plot_curves(results) \ No newline at end of file + plot_curves(results) diff --git a/benchmark/serialization/serialize_integrators.py b/benchmark/serialization/serialize_integrators.py index 7c7d544..4be433f 100644 --- a/benchmark/serialization/serialize_integrators.py +++ b/benchmark/serialization/serialize_integrators.py @@ -25,7 +25,7 @@ integrator = LangevinSplittingIntegrator(scheme) # Export integrator to XML - with open(os.path.join(DATA_PATH, "serialized_integrators/{}.xml".format(scheme), "w")) as f: + with open(os.path.join(DATA_PATH, "serialized_integrators/{}.xml".format(scheme), "wb")) as f: f.writelines(XmlSerializer.serialize(integrator)) # Also pretty-print (using `step_type_dict`) @@ -36,5 +36,5 @@ print(readable_lines) # Save pretty-printed results to .txt - with open(os.path.join(DATA_PATH, "serialized_integrators/{}.txt".format(scheme), "w")) as f: + with open(os.path.join(DATA_PATH, "serialized_integrators/{}.txt".format(scheme), "wb")) as f: f.writelines(readable_lines) diff --git a/benchmark/serialization/serialize_testsystems.py b/benchmark/serialization/serialize_testsystems.py index 9472ebe..043e749 100644 --- a/benchmark/serialization/serialize_testsystems.py +++ b/benchmark/serialization/serialize_testsystems.py @@ -8,5 +8,5 @@ top, sys, pos = system_params[name]["loader"](constrained) c = "constrained" if not constrained: c = "unconstrained" - with open("serialized_testsystems/{}_{}.xml".format(name, c), "w") as f: + with open("serialized_testsystems/{}_{}.xml".format(name, c), "wb") as f: f.writelines(XmlSerializer.serialize(sys)) diff --git a/tests/__init__.py b/benchmark/tests/__init__.py similarity index 100% rename from tests/__init__.py rename to benchmark/tests/__init__.py diff --git a/tests/measure_shadow_work.py b/benchmark/tests/measure_shadow_work.py similarity index 99% rename from tests/measure_shadow_work.py rename to benchmark/tests/measure_shadow_work.py index 9cdaad9..c8cb3f4 100644 --- a/tests/measure_shadow_work.py +++ b/benchmark/tests/measure_shadow_work.py @@ -89,4 +89,4 @@ def measure_shadow_work(simulation, n_steps): elif ("W_shad" in global_variable_names): return measure_shadow_work_via_W_shad(simulation, n_steps) else: - raise (RuntimeError("Simulation doesn't support shadow work computation")) \ No newline at end of file + raise (RuntimeError("Simulation doesn't support shadow work computation")) diff --git a/tests/substep_test.py b/benchmark/tests/substep_test.py similarity index 69% rename from tests/substep_test.py rename to benchmark/tests/substep_test.py index 1bcf2f6..80f2b33 100644 --- a/tests/substep_test.py +++ b/benchmark/tests/substep_test.py @@ -1,5 +1,5 @@ -from tests import compare_substep_energy_changes, simulation_factory -from code.utils import generate_solvent_solute_splitting_string +from .tests import compare_substep_energy_changes, simulation_factory +from benchmark.integrators import generate_solvent_solute_splitting_string for constrained in [True, False]: if constrained: print("\n\nWith constraints\n") @@ -12,4 +12,4 @@ ]: simulation = simulation_factory(scheme) print("Scheme: {}".format(scheme)) - compare_substep_energy_changes(simulation, n_steps=10) \ No newline at end of file + compare_substep_energy_changes(simulation, n_steps=10) diff --git a/tests/tests.py b/benchmark/tests/tests.py similarity index 94% rename from tests/tests.py rename to benchmark/tests/tests.py index c31a5e6..b6b54ae 100644 --- a/tests/tests.py +++ b/benchmark/tests/tests.py @@ -1,11 +1,13 @@ import numpy as np -from benchmark.integrators.integrators import LangevinSplittingIntegrator +from benchmark.integrators import LangevinSplittingIntegrator from openmmtools.testsystems import AlanineDipeptideVacuum from simtk import unit from simtk.openmm import app W_unit = unit.kilojoule_per_mole -from code.utils import configure_platform, get_total_energy, strip_unit, generate_solvent_solute_splitting_string +from benchmark.testsystems.configuration import configure_platform +from benchmark.utilities import get_total_energy, strip_unit +from benchmark.integrators import generate_solvent_solute_splitting_string from benchmark.testsystems.testsystems import load_alanine @@ -124,4 +126,4 @@ def check_W_shad_consistency(energy_changes): simulation.context.setPositions(testsystem.positions) simulation.context.setVelocitiesToTemperature(temperature) simulation.step(100) - check_W_shad_consistency(record_energy_changes(simulation, W_shad_name="shadow_work")) \ No newline at end of file + check_W_shad_consistency(record_energy_changes(simulation, W_shad_name="shadow_work")) diff --git a/benchmark/testsystems/__init__.py b/benchmark/testsystems/__init__.py index 5711157..ab22773 100644 --- a/benchmark/testsystems/__init__.py +++ b/benchmark/testsystems/__init__.py @@ -1,9 +1,11 @@ -from alanine_dipeptide import alanine_constrained, alanine_unconstrained, solvated_alanine_unconstrained -from waterbox import flexible_waterbox, waterbox_constrained -from coupled_power_oscillators import coupled_power_oscillators -from low_dimensional_systems import quartic, NumbaNonequilibriumSimulator -from bookkeepers import EquilibriumSimulator, NonequilibriumSimulator +from .alanine_dipeptide import alanine_constrained, alanine_unconstrained +#from .alanine_dipeptide import solvated_alanine_constrained, solvated_alanine_unconstrained # TODO: Uncomment when this is working +from .waterbox import flexible_waterbox, waterbox_constrained +from .coupled_power_oscillators import coupled_power_oscillators +from .low_dimensional_systems import quartic, NumbaNonequilibriumSimulator +from .bookkeepers import EquilibriumSimulator, NonequilibriumSimulator +from .testsystems import system_params __all__ = ["alanine_constrained", "alanine_unconstrained", "solvated_alanine_unconstrained", "flexible_waterbox", "waterbox_constrained", "coupled_power_oscillators", - "EquilibriumSimulator", "NonequilibriumSimulator", "quartic", "NumbaNonequilibriumSimulator"] \ No newline at end of file + "EquilibriumSimulator", "NonequilibriumSimulator", "quartic", "NumbaNonequilibriumSimulator"] diff --git a/benchmark/testsystems/alanine_dipeptide.py b/benchmark/testsystems/alanine_dipeptide.py index 701c2cd..ac44aad 100644 --- a/benchmark/testsystems/alanine_dipeptide.py +++ b/benchmark/testsystems/alanine_dipeptide.py @@ -4,7 +4,7 @@ from openmmtools.testsystems import CustomExternalForcesTestSystem, AlanineDipeptideVacuum, WaterBox, AlanineDipeptideExplicit, SrcImplicit from simtk.openmm import app from simtk import unit -from configuration import configure_platform +from benchmark.testsystems.configuration import configure_platform from benchmark.utilities import keep_only_some_forces @@ -34,27 +34,35 @@ def load_solvated_alanine(constrained=True): temperature = 298 * unit.kelvin -from bookkeepers import EquilibriumSimulator +from benchmark.testsystems.bookkeepers import EquilibriumSimulator top, sys, pos = load_alanine(constrained=True) alanine_constrained = EquilibriumSimulator(platform=configure_platform("Reference"), topology=top, system=sys, positions=pos, temperature=temperature, - ghmc_timestep=2.5 * unit.femtosecond, - burn_in_length=1000, n_samples=10000, - thinning_interval=100, name="alanine_constrained") + ghmc_timestep=1.0 * unit.femtosecond, + burn_in_length=50000, n_samples=1000, + thinning_interval=10000, name="alanine_constrained") top, sys, pos = load_alanine(constrained=False) alanine_unconstrained = EquilibriumSimulator(platform=configure_platform("Reference"), topology=top, system=sys, positions=pos, temperature=temperature, - ghmc_timestep=2.0 * unit.femtosecond, - burn_in_length=1000, n_samples=10000, - thinning_interval=100, name="alanine_unconstrained") + ghmc_timestep=1.0 * unit.femtosecond, + burn_in_length=50000, n_samples=1000, + thinning_interval=10000, name="alanine_unconstrained") -top, sys, pos = load_solvated_alanine(constrained=False) -solvated_alanine_unconstrained = EquilibriumSimulator(platform=configure_platform("CPU"), - topology=top, system=sys, positions=pos, - temperature=temperature, - ghmc_timestep=0.25 * unit.femtosecond, - burn_in_length=1000, n_samples=1000, - thinning_interval=5, name="solvated_alanine_unconstrained") \ No newline at end of file +#top, sys, pos = load_solvated_alanine(constrained=False) +#solvated_alanine_unconstrained = EquilibriumSimulator(platform=configure_platform("CPU"), +# topology=top, system=sys, positions=pos, +# temperature=temperature, +# ghmc_timestep=0.25 * unit.femtosecond, +# burn_in_length=200000, n_samples=1000, +# thinning_interval=200000, name="solvated_alanine_unconstrained") +# +#top, sys, pos = load_solvated_alanine(constrained=True) +#solvated_alanine_constrained = EquilibriumSimulator(platform=configure_platform("CPU"), +# topology=top, system=sys, positions=pos, +# temperature=temperature, +# ghmc_timestep=0.25 * unit.femtosecond, +# burn_in_length=200000, n_samples=1000, +# thinning_interval=200000, name="solvated_alanine_constrained") diff --git a/benchmark/testsystems/bookkeepers.py b/benchmark/testsystems/bookkeepers.py index 55ef056..0f3b1b5 100644 --- a/benchmark/testsystems/bookkeepers.py +++ b/benchmark/testsystems/bookkeepers.py @@ -88,6 +88,7 @@ def reset_ghmc_statistics(self): def collect_equilibrium_samples(self): """Collect equilibrium samples, return as (n_samples, n_atoms, 3) numpy array""" + print("Collecting equilibrium samples for '%s'..." % self.name) # Minimize energy by gradient descent print("Minimizing...") minimizer = GradientDescentMinimizationIntegrator() @@ -184,6 +185,11 @@ def accumulate_shadow_work(self, x_0, v_0, n_steps): set_positions(self.simulation, x_0) set_velocities(self.simulation, v_0) + # Apply position and velocity constraints. + tol = self.simulation.context.getIntegrator().getConstraintTolerance() + self.simulation.context.applyConstraints(tol) + self.simulation.context.applyVelocityConstraints(tol) + E_0 = get_energy() Q_0 = get_heat() @@ -221,4 +227,4 @@ def collect_protocol_samples(self, n_protocol_samples, protocol_length, marginal if __name__ == "__main__": - pass \ No newline at end of file + pass diff --git a/benchmark/testsystems/configuration.py b/benchmark/testsystems/configuration.py index 39c2c61..8a0c87b 100644 --- a/benchmark/testsystems/configuration.py +++ b/benchmark/testsystems/configuration.py @@ -1,19 +1,55 @@ import simtk.openmm as mm +# Check if we support the OpenCL platform in mixed precision +_opencl_supports_mixed_precision = True +from openmmtools.testsystems import HarmonicOscillator +try: + integrator = mm.VerletIntegrator(1.0) + testsystem = HarmonicOscillator() + platform = mm.Platform.getPlatformByName('OpenCL') + properties = {'OpenCLPrecision' : 'mixed'} + context = mm.Context(testsystem.system, integrator, platform, properties) + del context, properties, testsytem, integrator +except Exception as e: + print('OpenCL unavailable or does not support mixed precision') + _opencl_supports_mixed_precision = False -def configure_platform(platform_name='Reference'): - """Set precision, etc...""" - if platform_name.upper() == 'Reference'.upper(): - platform = mm.Platform.getPlatformByName('Reference') - elif platform_name.upper() == "CPU": - platform = mm.Platform.getPlatformByName("CPU") - elif platform_name.upper() == 'OpenCL'.upper(): - platform = mm.Platform.getPlatformByName('OpenCL') - platform.setPropertyDefaultValue('OpenCLPrecision', 'mixed') - elif platform_name.upper() == 'CUDA'.upper(): - platform = mm.Platform.getPlatformByName('CUDA') - platform.setPropertyDefaultValue('CUDAPrecision', 'double') - else: - raise(ValueError("Invalid platform name")) - return platform \ No newline at end of file +def configure_platform(platform_name='Reference', fallback_platform_name='CPU'): + """ + Retrieve the requested platform with platform-appropriate precision settings. + + platform_name : str, optional, default='Reference' + The requested platform name + fallback_platform_name : str, optional, default='CPU' + If the requested platform cannot be provided, the fallback platform will be provided. + + Returns + ------- + platform : simtk.openmm.Platform + The requested platform with precision configured appropriately, + or the fallback platform if this is not available. + + """ + try: + if platform_name.upper() == 'Reference'.upper(): + platform = mm.Platform.getPlatformByName('Reference') + elif platform_name.upper() == "CPU": + platform = mm.Platform.getPlatformByName("CPU") + elif platform_name.upper() == 'OpenCL'.upper(): + if not _opencl_supports_mixed_precision: + # Return CPU platform if OpenCL doesn't support mixed precision + return mm.Platform.getPlatformByName(fallback_platform_name) + platform = mm.Platform.getPlatformByName('OpenCL') + platform.setPropertyDefaultValue('OpenCLPrecision', 'mixed') + elif platform_name.upper() == 'CUDA'.upper(): + platform = mm.Platform.getPlatformByName('CUDA') + platform.setPropertyDefaultValue('CUDAPrecision', 'double') + else: + raise(ValueError("Invalid platform name")) + except: + # Return CPU platform if we can't provide requested platofrm + print("Warning: Returning '%s' platform instead of requested platform '%s'" % (fallback_platform_name, platform_name)) + return mm.Platform.getPlatformByName(fallback_platform_name) + + return platform diff --git a/benchmark/testsystems/coupled_power_oscillators.py b/benchmark/testsystems/coupled_power_oscillators.py index 31c50c6..5d06367 100644 --- a/benchmark/testsystems/coupled_power_oscillators.py +++ b/benchmark/testsystems/coupled_power_oscillators.py @@ -6,8 +6,8 @@ from openmmtools.testsystems import TestSystem -from bookkeepers import EquilibriumSimulator -from configuration import configure_platform +from benchmark.testsystems.bookkeepers import EquilibriumSimulator +from benchmark.testsystems.configuration import configure_platform class CoupledPowerOscillators(TestSystem): @@ -112,7 +112,7 @@ def __init__(self, nx=3, ny=3, nz=3, # Add these HarmonicBondForces to the system force = openmm.HarmonicBondForce() for bond in bonds: - force.addBond(bond[0], bond[1], grid_spacing, coupling_strength) + force.addBond(int(bond[0]), int(bond[1]), grid_spacing, coupling_strength) system.addForce(force) # Create topology. @@ -136,4 +136,4 @@ def __init__(self, nx=3, ny=3, nz=3, temperature=temperature, ghmc_timestep=1.0 * unit.femtosecond, burn_in_length=1000, n_samples=1000, - thinning_interval=10, name="coupled_power_oscillators") \ No newline at end of file + thinning_interval=10, name="coupled_power_oscillators") diff --git a/benchmark/testsystems/low_dimensional_systems.py b/benchmark/testsystems/low_dimensional_systems.py index cc88afc..ddcd64c 100644 --- a/benchmark/testsystems/low_dimensional_systems.py +++ b/benchmark/testsystems/low_dimensional_systems.py @@ -3,7 +3,7 @@ import os from tqdm import tqdm -from bookkeepers import BookkeepingSimulator, NonequilibriumSimulator +from benchmark.testsystems.bookkeepers import BookkeepingSimulator, NonequilibriumSimulator from benchmark.integrators import metropolis_hastings_factory from benchmark import DATA_PATH @@ -163,4 +163,4 @@ def collect_protocol_samples(self, n_protocol_samples, protocol_length, marginal W_shads_R.append(self.accumulate_shadow_work(x_1, v_1, protocol_length)) - return np.array(W_shads_F), np.array(W_shads_R) \ No newline at end of file + return np.array(W_shads_F), np.array(W_shads_R) diff --git a/benchmark/testsystems/testsystems.py b/benchmark/testsystems/testsystems.py index 1091029..8f69c4a 100644 --- a/benchmark/testsystems/testsystems.py +++ b/benchmark/testsystems/testsystems.py @@ -2,7 +2,7 @@ from openmmtools.testsystems import SrcImplicit, DHFRExplicit from simtk.openmm import app from simtk import unit -from configuration import configure_platform +from benchmark.testsystems.configuration import configure_platform from benchmark.utilities import keep_only_some_forces @@ -80,14 +80,19 @@ def get_src_explicit_test_system(temperature): } harmonic_oscillator_params = simple_params.copy() +from .low_dimensional_systems import load_harmonic_oscillator harmonic_oscillator_params["loader"] = load_harmonic_oscillator quartic_params = simple_params.copy() +from .low_dimensional_systems import load_quartic_potential quartic_params["loader"] = load_quartic_potential mts_params = simple_params.copy() +from .low_dimensional_systems import load_mts_test mts_params["loader"] = load_mts_test +from .waterbox import load_waterbox +from .alanine_dipeptide import load_alanine system_params = { "harmonic_oscillator": harmonic_oscillator_params, "quartic_potential": quartic_params, @@ -115,4 +120,4 @@ def get_src_explicit_test_system(temperature): "collision_rate": 91 / unit.picoseconds, } } -# TODO: Add Waterbox, AlanineExplicit EquilibriumSimulators \ No newline at end of file +# TODO: Add Waterbox, AlanineExplicit EquilibriumSimulators diff --git a/benchmark/testsystems/waterbox.py b/benchmark/testsystems/waterbox.py index 4987c4d..e66f2f5 100644 --- a/benchmark/testsystems/waterbox.py +++ b/benchmark/testsystems/waterbox.py @@ -2,7 +2,7 @@ from openmmtools.testsystems import CustomExternalForcesTestSystem, AlanineDipeptideVacuum, WaterBox, AlanineDipeptideExplicit, SrcImplicit from simtk.openmm import app from simtk import unit -from configuration import configure_platform +from benchmark.testsystems.configuration import configure_platform from benchmark.utilities import keep_only_some_forces from benchmark import DATA_PATH @@ -14,7 +14,7 @@ def load_waterbox(constrained=True): temperature = 298 * unit.kelvin -from bookkeepers import EquilibriumSimulator +from benchmark.testsystems.bookkeepers import EquilibriumSimulator top, sys, pos = load_waterbox(constrained=True) waterbox_constrained = EquilibriumSimulator(platform=configure_platform("OpenCL"), topology=top, system=sys, positions=pos, diff --git a/benchmark/utilities/__init__.py b/benchmark/utilities/__init__.py index e1ec41c..cd98316 100644 --- a/benchmark/utilities/__init__.py +++ b/benchmark/utilities/__init__.py @@ -1,10 +1,10 @@ -from utils import strip_unit, stderr,\ - summarize, print_array +from .utils import strip_unit, stderr,\ + summarize, print_array, get_summary_string, unpack_trajs -from openmm_utilities import keep_only_some_forces,\ +from .openmm_utilities import keep_only_some_forces,\ get_total_energy, get_positions, get_velocities,\ set_positions, set_velocities __all__ = ["strip_unit", "get_total_energy", "stderr", "summarize", "get_positions", "get_velocities", "keep_only_some_forces", - "set_positions", "set_velocities", "print_array"] \ No newline at end of file + "set_positions", "set_velocities", "print_array"] diff --git a/data/waterbox.npy b/data/waterbox.npy deleted file mode 100644 index d5b20dd..0000000 Binary files a/data/waterbox.npy and /dev/null differ diff --git a/references/Efficient molecular dynamics using geodesic integration and solute-solvent splitting.pdf b/references/Efficient molecular dynamics using geodesic integration and solute-solvent splitting.pdf new file mode 100644 index 0000000..de0a60d Binary files /dev/null and b/references/Efficient molecular dynamics using geodesic integration and solute-solvent splitting.pdf differ