From 7ed8d94327deab15a5552d8a46dba4a348965ae9 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 8 Mar 2024 14:36:04 +0000 Subject: [PATCH 1/4] Flag perturbable atom as light by max mass and check constraints. --- src/somd2/runner/_dynamics.py | 10 +++---- src/somd2/runner/_runner.py | 52 ++++++++++++++++++++++++++++++++--- 2 files changed, 52 insertions(+), 10 deletions(-) diff --git a/src/somd2/runner/_dynamics.py b/src/somd2/runner/_dynamics.py index 44a3ac2..985a9b6 100644 --- a/src/somd2/runner/_dynamics.py +++ b/src/somd2/runner/_dynamics.py @@ -158,11 +158,6 @@ def _setup_dynamics(self, equilibration=False): else: pressure = None - try: - map = self._config._extra_args - except: - map = None - self._dyn = self._system.dynamics( temperature=self._config.temperature, pressure=pressure, @@ -186,7 +181,7 @@ def _setup_dynamics(self, equilibration=False): swap_end_states=self._config.swap_end_states, com_reset_frequency=self._config.com_reset_frequency, vacuum=not self._has_space, - map=map, + map=self._config._extra_args, ) def _minimisation(self, lambda_min=None, perturbable_constraint="none"): @@ -200,6 +195,7 @@ def _minimisation(self, lambda_min=None, perturbable_constraint="none"): Lambda value at which to run minimisation, if None run at pre-set lambda_val. """ + if lambda_min is None: _logger.info(f"Minimising at {_lam_sym} = {self._lambda_val}") try: @@ -210,6 +206,7 @@ def _minimisation(self, lambda_min=None, perturbable_constraint="none"): lambda_value=self._lambda_val, platform=self._config.platform, vacuum=not self._has_space, + constraint=self._config.constraint, perturbable_constraint=perturbable_constraint, include_constrained_energies=self._config.include_constrained_energies, dynamic_constraints=self._config.dynamic_constraints, @@ -230,6 +227,7 @@ def _minimisation(self, lambda_min=None, perturbable_constraint="none"): lambda_value=lambda_min, platform=self._config.platform, vacuum=not self._has_space, + constraint=self._config.constraint, perturbable_constraint=perturbable_constraint, include_constrained_energies=self._config.include_constrained_energies, dynamic_constraints=self._config.dynamic_constraints, diff --git a/src/somd2/runner/_runner.py b/src/somd2/runner/_runner.py index 0ab4f34..6b13381 100644 --- a/src/somd2/runner/_runner.py +++ b/src/somd2/runner/_runner.py @@ -94,6 +94,7 @@ def __init__(self, system, config): if not isinstance(config, _Config): raise TypeError("'config' must be of type 'somd2.config.Config'") self._config = config + self._config._extra_args = {} # We're running in SOMD1 compatibility mode. if self._config.somd1_compatibility: @@ -138,9 +139,21 @@ def __init__(self, system, config): "Unable to convert water topology to AMBER format for SOMD1 compatibility." ) + # Only check for light atoms by the maxium end state mass if running + # in SOMD1 compatibility mode. + if self._config.somd1_compatibility: + self._config._extra_args["check_for_h_by_max_mass"] = True + self._config._extra_args["check_for_h_by_mass"] = False + self._config._extra_args["check_for_h_by_mass"] = False + self._config._extra_args["check_for_h_by_element"] = False + self._config._extra_args["check_for_h_by_ambertype"] = False + # Check for a periodic space. self._check_space() + # Check the end state contraints. + self._check_end_state_constraints() + # Set the lambda values. self._lambda_values = [ round(i / (self._config.num_lambda - 1), 5) @@ -235,6 +248,39 @@ def _check_space(self): ) self._config.cutoff_type = "rf" + def _check_end_state_constraints(self): + """ + Internal function to check whether the constrants are the same at the two + end states. + """ + + # Assume a single perturbable molecule for now. + pert_mol = self._system.molecules("property is_perturbable")[0] + + # Create a dynamics object. + d = pert_mol.dynamics( + constraint=self._config.constraint, + perturbable_constraint=self._config.perturbable_constraint, + platform="cpu", + map=self._config._extra_args, + ) + + # Get the constraints at lambda = 0. + constraints0 = d.get_constraints() + + # Update to lambda = 1. + d.set_lambda(1) + + # Get the constraints at lambda = 1. + constraints1 = d.get_constraints() + + # Check for equivalence. + for c0, c1 in zip(constraints0, constraints1): + if c0 != c1: + _logger.info( + f"Constraints are at not the same at {_lam_sym} = 0 and {_lam_sym} = 1." + ) + def _check_directory(self): """ Find the name of the file from which simulations will be started. @@ -657,17 +703,15 @@ def run(self): threads_per_worker = ( self._config.max_threads // self._config.num_lambda ) - self._config._extra_args = {"threads": threads_per_worker} + self._config._extra_args["threads"] = threads_per_worker # (Multi-)GPU platform. elif self._is_gpu: self.max_workers = len(self._gpu_pool) - self._config._extra_args = {} # All other platforms. else: self._max_workers = 1 - self._config._extra_args = {} import concurrent.futures as _futures @@ -696,7 +740,7 @@ def run(self): # Serial configuration. elif self._config.num_lambda is not None: if self._config.platform == "cpu": - self._config.extra_args = {"threads": self._config.max_threads} + self._config._extra_args = {"threads": self._config.max_threads} self._lambda_values = [ round(i / (self._config.num_lambda - 1), 5) for i in range(0, self._config.num_lambda) From b57619e9efd4547ecc23b5f249a325948b31cc2e Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 Mar 2024 09:08:17 +0000 Subject: [PATCH 2/4] Check constraints for all perturbable molecules. --- src/somd2/runner/_runner.py | 46 +++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 22 deletions(-) diff --git a/src/somd2/runner/_runner.py b/src/somd2/runner/_runner.py index 6b13381..bac3f1b 100644 --- a/src/somd2/runner/_runner.py +++ b/src/somd2/runner/_runner.py @@ -254,32 +254,34 @@ def _check_end_state_constraints(self): end states. """ - # Assume a single perturbable molecule for now. - pert_mol = self._system.molecules("property is_perturbable")[0] - - # Create a dynamics object. - d = pert_mol.dynamics( - constraint=self._config.constraint, - perturbable_constraint=self._config.perturbable_constraint, - platform="cpu", - map=self._config._extra_args, - ) + # Find all perturbable molecules in the system.. + pert_mols = self._system.molecules("property is_perturbable") + + # Check constraints at lambda = 0 and lambda = 1 for each perturbable molecule. + for mol in pert_mols: + # Create a dynamics object. + d = mol.dynamics( + constraint=self._config.constraint, + perturbable_constraint=self._config.perturbable_constraint, + platform="cpu", + map=self._config._extra_args, + ) - # Get the constraints at lambda = 0. - constraints0 = d.get_constraints() + # Get the constraints at lambda = 0. + constraints0 = d.get_constraints() - # Update to lambda = 1. - d.set_lambda(1) + # Update to lambda = 1. + d.set_lambda(1) - # Get the constraints at lambda = 1. - constraints1 = d.get_constraints() + # Get the constraints at lambda = 1. + constraints1 = d.get_constraints() - # Check for equivalence. - for c0, c1 in zip(constraints0, constraints1): - if c0 != c1: - _logger.info( - f"Constraints are at not the same at {_lam_sym} = 0 and {_lam_sym} = 1." - ) + # Check for equivalence. + for c0, c1 in zip(constraints0, constraints1): + if c0 != c1: + _logger.info( + f"Constraints are at not the same at {_lam_sym} = 0 and {_lam_sym} = 1." + ) def _check_directory(self): """ From 00886a816f13b5df941ea21317a27e76648b1514 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 Mar 2024 09:14:09 +0000 Subject: [PATCH 3/4] Handle "None" as a command-line option for pressure. --- src/somd2/config/_config.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py index 0718629..fa04d06 100644 --- a/src/somd2/config/_config.py +++ b/src/somd2/config/_config.py @@ -432,6 +432,10 @@ def pressure(self, pressure): try: p = _sr.u(pressure) except: + # Handle special case of pressure = "none" + if pressure.lower().replace(" ", "") == "none": + self._pressure = None + return raise ValueError( f"Unable to parse 'pressure' as a Sire GeneralUnit: {pressure}" ) From d3bc06fd38c82b079d72e7ebdccf82beae437d19 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 Mar 2024 09:32:50 +0000 Subject: [PATCH 4/4] Expose barostat_frequency kwarg. --- src/somd2/config/_config.py | 19 +++++++++++++++++++ src/somd2/runner/_dynamics.py | 1 + 2 files changed, 20 insertions(+) diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py index fa04d06..e24e05c 100644 --- a/src/somd2/config/_config.py +++ b/src/somd2/config/_config.py @@ -74,6 +74,7 @@ def __init__( timestep="4fs", temperature="300K", pressure="1 atm", + barostat_frequency=25, integrator="langevin_middle", cutoff_type="pme", cutoff="7.5A", @@ -125,6 +126,9 @@ def __init__( pressure: str Simulation pressure. (Simulations will run in the NVT ensemble unless a pressure is specified.) + barostat_frequency: int + The number of integration steps between barostat updates. + integrator: str Integrator to use for simulation. @@ -249,6 +253,7 @@ def __init__( self.runtime = runtime self.temperature = temperature self.pressure = pressure + self.barostat_frequency = barostat_frequency self.integrator = integrator self.cutoff_type = cutoff_type self.cutoff = cutoff @@ -447,6 +452,20 @@ def pressure(self, pressure): else: self._pressure = pressure + @property + def barostat_frequency(self): + return self._barostat_frequency + + @barostat_frequency.setter + def barostat_frequency(self, barostat_frequency): + if not isinstance(barostat_frequency, int): + raise TypeError("'barostat_frequency' must be of type 'int'") + + if barostat_frequency <= 0: + raise ValueError("'barostat_frequency' must be a positive integer") + + self._barostat_frequency = barostat_frequency + @property def integrator(self): return self._integrator diff --git a/src/somd2/runner/_dynamics.py b/src/somd2/runner/_dynamics.py index 985a9b6..6c97993 100644 --- a/src/somd2/runner/_dynamics.py +++ b/src/somd2/runner/_dynamics.py @@ -161,6 +161,7 @@ def _setup_dynamics(self, equilibration=False): self._dyn = self._system.dynamics( temperature=self._config.temperature, pressure=pressure, + barostat_frequency=self._config.barostat_frequency, timestep=( self._config.equilibration_timestep if equilibration