From 40f7b161ca04c9edf85dfb4c2c000dcbb6436caa Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 17 Jul 2024 11:06:53 +0100 Subject: [PATCH 01/11] Save both end state topologies for trajectory reconstruction. --- src/somd2/runner/_dynamics.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/somd2/runner/_dynamics.py b/src/somd2/runner/_dynamics.py index e7eeabb..53aa671 100644 --- a/src/somd2/runner/_dynamics.py +++ b/src/somd2/runner/_dynamics.py @@ -151,7 +151,8 @@ def increment_filename(base_filename, suffix): raise ValueError("lambda_value not in lambda_array") lam = f"{lambda_value:.5f}" filenames = {} - filenames["topology"] = "system.prm7" + filenames["topology0"] = "system0.prm7" + filenames["topology1"] = "system1.prm7" filenames["checkpoint"] = f"checkpoint_{lam}.s3" filenames["energy_traj"] = f"energy_traj_{lam}.parquet" filenames["trajectory"] = f"traj_{lam}.dcd" @@ -307,8 +308,12 @@ def _run(self, lambda_minimisation=None, is_restart=False): # Save the system topology to a PRM7 file that can be used to load the # trajectory. - topology = str(self._config.output_directory / self._filenames["topology"]) - sr.save(self._system, topology) + topology0 = str(self._config.output_directory / self._filenames["topology0"]) + topology1 = str(self._config.output_directory / self._filenames["topology1"]) + mols0 = sr.morph.link_to_reference(self._system) + mols1 = sr.morph.link_to_perturbed(self._system) + sr.save(mols0, topology0) + sr.save(mols1, topology1) def generate_lam_vals(lambda_base, increment): """Generate lambda values for a given lambda_base and increment""" @@ -560,7 +565,7 @@ def generate_lam_vals(lambda_base, increment): traj_chunks = [f"{traj_filename}.bak"] + traj_chunks # Load the topology and chunked trajectory files. - system = sr.load([topology] + traj_chunks) + system = sr.load([topology0] + traj_chunks) # Save the final trajectory to a single file. traj_filename = str( From 6dd4958073a01ab3e08f89b776fe2899f561d9b3 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 17 Jul 2024 11:14:42 +0100 Subject: [PATCH 02/11] Mamba is now the default with MiniForge. --- .github/workflows/main.yaml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/.github/workflows/main.yaml b/.github/workflows/main.yaml index de7dcef..32c05f4 100644 --- a/.github/workflows/main.yaml +++ b/.github/workflows/main.yaml @@ -55,12 +55,10 @@ jobs: activate-environment: somd2 environment-file: environment.yaml miniforge-version: latest - miniforge-variant: Mambaforge - use-mamba: true run-post: ${{ matrix.platform.name != 'windows' }} # - name: Install pytest - run: mamba install pytest + run: conda install pytest # - name: Install the package run: pip install . From 8a005c53655bf2746d64436f0b7fa34783ae5cdf Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 17 Jul 2024 11:14:58 +0100 Subject: [PATCH 03/11] Add support for alchemical ions. --- src/somd2/config/_config.py | 19 ++++++++ src/somd2/runner/_runner.py | 90 +++++++++++++++++++++++++++++++++++++ 2 files changed, 109 insertions(+) diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py index 51062f1..fbb733e 100644 --- a/src/somd2/config/_config.py +++ b/src/somd2/config/_config.py @@ -104,6 +104,7 @@ def __init__( perturbable_constraint="h_bonds_not_heavy_perturbed", include_constrained_energies=False, dynamic_constraints=True, + charge_difference=0, com_reset_frequency=10, minimise=True, equilibration_time="0ps", @@ -211,6 +212,12 @@ def __init__( to the value of r0 at that lambda value. If this is False, then the constraint is set based on the current length. + charge_difference: int + The charge difference between the two end states. (Perturbed minus + reference.) If specified, then a number of alchemical ions will be + added to the system to neutralise the charge difference, i.e. make + perturbed state charge the same as the reference state. + com_reset_frequency: int Frequency at which to reset the centre of mass of the system. @@ -315,6 +322,7 @@ def __init__( self.perturbable_constraint = perturbable_constraint self.include_constrained_energies = include_constrained_energies self.dynamic_constraints = dynamic_constraints + self.charge_difference = charge_difference self.com_reset_frequency = com_reset_frequency self.minimise = minimise self.equilibration_time = equilibration_time @@ -861,6 +869,17 @@ def dynamic_constraints(self, dynamic_constraints): raise ValueError("'dynamic_constraints' must be of type 'bool'") self._dynamic_constraints = dynamic_constraints + @property + def charge_difference(self): + return self._charge_difference + + @charge_difference.setter + def charge_difference(self, charge_difference): + if charge_difference is not None: + if not isinstance(charge_difference, int): + raise ValueError("'charge_difference' must be an integer") + self._charge_difference = charge_difference + @property def com_reset_frequency(self): return self._com_reset_frequency diff --git a/src/somd2/runner/_runner.py b/src/somd2/runner/_runner.py index 3ca0f05..8e33516 100644 --- a/src/somd2/runner/_runner.py +++ b/src/somd2/runner/_runner.py @@ -173,6 +173,28 @@ def __init__(self, system, config): # Check the end state contraints. self._check_end_state_constraints() + # Get the charge difference between the two end states. + charge_diff = self._get_charge_difference() + + # Make sure the difference is integer valued to 5 decimal places. + if not round(charge_diff, 4).is_integer(): + _logger.warning("Charge difference between end states is not an integer.") + charge_diff = round(charge_diff, 4) + + # Make sure the charge difference matches the expected value + # from the config. + if self._config.charge_difference != charge_diff: + _logger.warning( + f"The charge difference of {charge_diff} between the end states " + f"does not match the expected value of {self._config.charge_difference}. " + "Please specify the 'charge_difference' if you wish to maintain " + "charge neutrality." + ) + + # Create alchemical ions. + if self._config.charge_difference != 0: + self._create_alchemical_ions(self._system, self._config.charge_difference) + # Set the lambda values. if self._config.lambda_values: self._lambda_values = self._config.lambda_values @@ -318,6 +340,74 @@ def _check_end_state_constraints(self): ) break + def _get_charge_difference(self): + """ + Internal function to check the charge difference between the two end states. + """ + + reference = _morph.link_to_reference(self._system).charge().value() + perturbed = _morph.link_to_perturbed(self._system).charge().value() + + return perturbed - reference + + def _create_alchemical_ions(self, system, charge_diff): + """ + Internal function to create alchemical ions to maintain charge neutrality. + """ + + from sire.legacy.IO import createChlorineIon as _createChlorineIon + from sire.legacy.IO import createSodiumIon as _createSodiumIon + + # The number of waters to convert is the absolute charge difference. + num_waters = abs(charge_diff) + + # Reference coordinates. + coords = system.molecules("property is_perturbable").coordinates() + coord_string = f"{coords[0].value()}, {coords[1].value()}, {coords[2].value()}" + + # Find the furthest N waters from the perturbable molecule. + waters = self._system[ + f"furthest {num_waters} waters from {coord_string}" + ].molecules() + + # Determine the water model. + if waters[0].num_atoms() == 3: + model = "tip3p" + elif waters[0].num_atoms() == 4: + model = "tip4p" + elif waters[0].num_atoms() == 5: + # Note that AMBER has no ion model for tip5p. + model = "tip4p" + + # Store the molecule numbers for the system. + numbers = self._system.numbers() + + # Create the ions. + for water in waters: + # Create an ion. This is to adjust the charge of the perturbed state + # to match that of the reference. + if charge_diff > 0: + ion = _createChlorineIon(water["element O"].coordinates(), model) + ion_str = "Cl-" + else: + ion = _createSodiumIon(water["element O"].coordinates(), model) + ion_str = "Na+" + + # Create a perturbable molecule: water --> ion. + merged = _morph.merge(water, ion, map={"as_new_molecule": False}) + + # Update the system. + self._system.update(merged) + + # Get the index of the perturbed water. + index = numbers.index(water.number()) + + # Log that the water was perturbed. + _logger.info( + f"Water at molecule index {index} will be perturbed to " + f"{ion_str} to preserve charge neutrality." + ) + def _check_directory(self): """ Find the name of the file from which simulations will be started. From 0a7c83837655e4efd3b3ed27ae54aabcbbc60eea Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 17 Jul 2024 11:31:21 +0100 Subject: [PATCH 04/11] Convert the charge difference to an integer. --- src/somd2/runner/_runner.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/somd2/runner/_runner.py b/src/somd2/runner/_runner.py index 8e33516..78e58de 100644 --- a/src/somd2/runner/_runner.py +++ b/src/somd2/runner/_runner.py @@ -179,7 +179,7 @@ def __init__(self, system, config): # Make sure the difference is integer valued to 5 decimal places. if not round(charge_diff, 4).is_integer(): _logger.warning("Charge difference between end states is not an integer.") - charge_diff = round(charge_diff, 4) + charge_diff = int(round(charge_diff, 4)) # Make sure the charge difference matches the expected value # from the config. From 6b7095b3878497363cdfb5c48fafdd978c53575e Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 17 Jul 2024 12:16:49 +0100 Subject: [PATCH 05/11] Rename the reference topology file for trajectory loading. --- tests/runner/test_restart.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/runner/test_restart.py b/tests/runner/test_restart.py index a0e509f..5f9639d 100644 --- a/tests/runner/test_restart.py +++ b/tests/runner/test_restart.py @@ -45,7 +45,7 @@ def test_restart(mols, request): # Load the trajectory. traj_1 = sr.load( - [str(Path(tmpdir) / "system.prm7"), str(Path(tmpdir) / "traj_0.00000.dcd")] + [str(Path(tmpdir) / "system0.prm7"), str(Path(tmpdir) / "traj_0.00000.dcd")] ) # Check that both config and lambda have been written @@ -89,7 +89,7 @@ def test_restart(mols, request): # Reload the trajectory. traj_2 = sr.load( - [str(Path(tmpdir) / "system.prm7"), str(Path(tmpdir) / "traj_0.00000.dcd")] + [str(Path(tmpdir) / "system0.prm7"), str(Path(tmpdir) / "traj_0.00000.dcd")] ) # Check that the trajectory is twice as long as the first. From fa2c8f35d46882b1e4d7a526d1b9b989d91d1cb1 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 17 Jul 2024 12:27:02 +0100 Subject: [PATCH 06/11] Update Python version. --- .github/workflows/main.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/main.yaml b/.github/workflows/main.yaml index 32c05f4..9547046 100644 --- a/.github/workflows/main.yaml +++ b/.github/workflows/main.yaml @@ -28,7 +28,7 @@ jobs: max-parallel: 9 fail-fast: false matrix: - python-version: ["3.11"] + python-version: ["3.12"] platform: - { name: "windows", os: "windows-latest", shell: "pwsh" } - { name: "linux", os: "ubuntu-latest", shell: "bash -l {0}" } From 7012ab2d387bbe534c1e49c29112e612939eeefd Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 17 Jul 2024 15:11:42 +0100 Subject: [PATCH 07/11] Convert to staticmethods so they can be tested. --- src/somd2/runner/_runner.py | 52 ++++++++++++++++++++++++++++++------- 1 file changed, 42 insertions(+), 10 deletions(-) diff --git a/src/somd2/runner/_runner.py b/src/somd2/runner/_runner.py index 78e58de..8850c53 100644 --- a/src/somd2/runner/_runner.py +++ b/src/somd2/runner/_runner.py @@ -174,7 +174,7 @@ def __init__(self, system, config): self._check_end_state_constraints() # Get the charge difference between the two end states. - charge_diff = self._get_charge_difference() + charge_diff = self._get_charge_difference(self._system) # Make sure the difference is integer valued to 5 decimal places. if not round(charge_diff, 4).is_integer(): @@ -340,24 +340,56 @@ def _check_end_state_constraints(self): ) break - def _get_charge_difference(self): + @staticmethod + def _get_charge_difference(system): """ Internal function to check the charge difference between the two end states. + + Parameters + ---------- + + system: :class: `System ` + The system to be perturbed. + + Returns + ------- + + charge_diff: int + The charge difference between the perturbed and reference states. """ - reference = _morph.link_to_reference(self._system).charge().value() - perturbed = _morph.link_to_perturbed(self._system).charge().value() + reference = _morph.link_to_reference(system).charge().value() + perturbed = _morph.link_to_perturbed(system).charge().value() return perturbed - reference - def _create_alchemical_ions(self, system, charge_diff): + @staticmethod + def _create_alchemical_ions(system, charge_diff): """ Internal function to create alchemical ions to maintain charge neutrality. + + Parameters + ---------- + + system: :class: `System ` + The system to be perturbed. + + charge_diff: int + The charge difference between perturbed and reference states. + + Returns + ------- + + system: :class: `System ` + The perturbed system with alchemical ions added. """ from sire.legacy.IO import createChlorineIon as _createChlorineIon from sire.legacy.IO import createSodiumIon as _createSodiumIon + # Clone the system. + system = system.clone() + # The number of waters to convert is the absolute charge difference. num_waters = abs(charge_diff) @@ -366,9 +398,7 @@ def _create_alchemical_ions(self, system, charge_diff): coord_string = f"{coords[0].value()}, {coords[1].value()}, {coords[2].value()}" # Find the furthest N waters from the perturbable molecule. - waters = self._system[ - f"furthest {num_waters} waters from {coord_string}" - ].molecules() + waters = system[f"furthest {num_waters} waters from {coord_string}"].molecules() # Determine the water model. if waters[0].num_atoms() == 3: @@ -380,7 +410,7 @@ def _create_alchemical_ions(self, system, charge_diff): model = "tip4p" # Store the molecule numbers for the system. - numbers = self._system.numbers() + numbers = system.numbers() # Create the ions. for water in waters: @@ -397,7 +427,7 @@ def _create_alchemical_ions(self, system, charge_diff): merged = _morph.merge(water, ion, map={"as_new_molecule": False}) # Update the system. - self._system.update(merged) + system.update(merged) # Get the index of the perturbed water. index = numbers.index(water.number()) @@ -408,6 +438,8 @@ def _create_alchemical_ions(self, system, charge_diff): f"{ion_str} to preserve charge neutrality." ) + return system + def _check_directory(self): """ Find the name of the file from which simulations will be started. From f12716a29f06b11ea732462fd5e2792dd88ae195 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 17 Jul 2024 15:23:14 +0100 Subject: [PATCH 08/11] Add test for creation of alchemical ions. --- tests/runner/test_alchemical_ions.py | 31 ++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) create mode 100644 tests/runner/test_alchemical_ions.py diff --git a/tests/runner/test_alchemical_ions.py b/tests/runner/test_alchemical_ions.py new file mode 100644 index 0000000..c808406 --- /dev/null +++ b/tests/runner/test_alchemical_ions.py @@ -0,0 +1,31 @@ +import math + +from somd2.runner import Runner + + +def test_alchemical_ions(ethane_methanol): + """Ensure that alchemical ions are added correctly.""" + + # Clone the system. + mols = ethane_methanol.clone() + + # Helper function to return the charge difference between the end states. + def charge_difference(mols): + from sire import morph + + reference = morph.link_to_reference(mols) + perturbed = morph.link_to_perturbed(mols) + + return (perturbed.charge() - reference.charge()).value() + + # Add 10 Cl- ions. + new_mols = Runner._create_alchemical_ions(mols, 10) + + # Make sure the charge difference is correct. + assert math.isclose(charge_difference(new_mols), -10.0, rel_tol=1e-6) + + # Add 10 Na+ ions. + new_mols = Runner._create_alchemical_ions(mols, -10) + + # Make sure the charge difference is correct. + assert math.isclose(charge_difference(new_mols), 10.0, rel_tol=1e-6) From cc2df101451c1037c9517a7619b2fb4982928d66 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 17 Jul 2024 16:09:49 +0100 Subject: [PATCH 09/11] Use staticmethod for calculating the charge difference. [ci skip] --- tests/runner/test_alchemical_ions.py | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/tests/runner/test_alchemical_ions.py b/tests/runner/test_alchemical_ions.py index c808406..47e1b07 100644 --- a/tests/runner/test_alchemical_ions.py +++ b/tests/runner/test_alchemical_ions.py @@ -9,23 +9,14 @@ def test_alchemical_ions(ethane_methanol): # Clone the system. mols = ethane_methanol.clone() - # Helper function to return the charge difference between the end states. - def charge_difference(mols): - from sire import morph - - reference = morph.link_to_reference(mols) - perturbed = morph.link_to_perturbed(mols) - - return (perturbed.charge() - reference.charge()).value() - # Add 10 Cl- ions. new_mols = Runner._create_alchemical_ions(mols, 10) # Make sure the charge difference is correct. - assert math.isclose(charge_difference(new_mols), -10.0, rel_tol=1e-6) + assert math.isclose(Runner._get_charge_difference(new_mols), -10.0, rel_tol=1e-6) # Add 10 Na+ ions. new_mols = Runner._create_alchemical_ions(mols, -10) # Make sure the charge difference is correct. - assert math.isclose(charge_difference(new_mols), 10.0, rel_tol=1e-6) + assert math.isclose(Runner._get_charge_difference(new_mols), 10.0, rel_tol=1e-6) From 6feceb9c42916c04f7e4e5363822ba66d1eb611a Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 19 Jul 2024 09:17:13 +0100 Subject: [PATCH 10/11] Refer to keeping charge constant rather than neutralising. [ci skip] --- src/somd2/config/_config.py | 2 +- src/somd2/runner/_runner.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py index fbb733e..5cc1e96 100644 --- a/src/somd2/config/_config.py +++ b/src/somd2/config/_config.py @@ -215,7 +215,7 @@ def __init__( charge_difference: int The charge difference between the two end states. (Perturbed minus reference.) If specified, then a number of alchemical ions will be - added to the system to neutralise the charge difference, i.e. make + added to the system to keep the charge constant, i.e. make the perturbed state charge the same as the reference state. com_reset_frequency: int diff --git a/src/somd2/runner/_runner.py b/src/somd2/runner/_runner.py index 8850c53..e86bcce 100644 --- a/src/somd2/runner/_runner.py +++ b/src/somd2/runner/_runner.py @@ -187,8 +187,8 @@ def __init__(self, system, config): _logger.warning( f"The charge difference of {charge_diff} between the end states " f"does not match the expected value of {self._config.charge_difference}. " - "Please specify the 'charge_difference' if you wish to maintain " - "charge neutrality." + "Please specify the 'charge_difference' if you wish to keep the charge " + "constant." ) # Create alchemical ions. @@ -366,7 +366,7 @@ def _get_charge_difference(system): @staticmethod def _create_alchemical_ions(system, charge_diff): """ - Internal function to create alchemical ions to maintain charge neutrality. + Internal function to create alchemical ions to maintain a constant charge. Parameters ---------- @@ -435,7 +435,7 @@ def _create_alchemical_ions(system, charge_diff): # Log that the water was perturbed. _logger.info( f"Water at molecule index {index} will be perturbed to " - f"{ion_str} to preserve charge neutrality." + f"{ion_str} to keep charge constant." ) return system From ab86a7358d5b4af8bf273072c3bb7e928bbf6878 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 19 Jul 2024 12:09:44 +0100 Subject: [PATCH 11/11] Update system with return value. --- src/somd2/runner/_runner.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/somd2/runner/_runner.py b/src/somd2/runner/_runner.py index e86bcce..a4b616e 100644 --- a/src/somd2/runner/_runner.py +++ b/src/somd2/runner/_runner.py @@ -193,7 +193,9 @@ def __init__(self, system, config): # Create alchemical ions. if self._config.charge_difference != 0: - self._create_alchemical_ions(self._system, self._config.charge_difference) + self._system = self._create_alchemical_ions( + self._system, self._config.charge_difference + ) # Set the lambda values. if self._config.lambda_values: