Skip to content

Commit

Permalink
Default to None for charge_difference and check there are enough waters.
Browse files Browse the repository at this point in the history
  • Loading branch information
lohedges committed Aug 9, 2024
1 parent 5a56168 commit bfa8b99
Show file tree
Hide file tree
Showing 2 changed files with 23 additions and 4 deletions.
7 changes: 5 additions & 2 deletions src/somd2/config/_config.py
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ def __init__(
perturbable_constraint="h_bonds_not_heavy_perturbed",
include_constrained_energies=False,
dynamic_constraints=True,
charge_difference=0,
charge_difference=None,
com_reset_frequency=10,
minimise=True,
equilibration_time="0ps",
Expand Down Expand Up @@ -877,7 +877,10 @@ def charge_difference(self):
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")
try:
charge_difference = int(charge_difference)
except:
raise ValueError("'charge_difference' must be an integer")
self._charge_difference = charge_difference

@property
Expand Down
20 changes: 18 additions & 2 deletions src/somd2/runner/_runner.py
Original file line number Diff line number Diff line change
Expand Up @@ -183,14 +183,22 @@ def __init__(self, system, config):

# Make sure the charge difference matches the expected value
# from the config.
if self._config.charge_difference != charge_diff:
if (
self._config.charge_difference is not None
and self._config.charge_difference != charge_diff
):
_logger.warning(
f"The charge difference of {charge_diff} between the end states "
f"does not match the specified value of {self._config.charge_difference}"
)
else:
_logger.info(
f"There is a charge difference of {charge_diff} between the end states. "
f"Adding alchemical ions to keep the charge constant."
)

# The user value takes precedence.
if self._config.charge_difference != 0:
if self._config.charge_difference is not None:
charge_diff = self._config.charge_difference

# Create alchemical ions.
Expand Down Expand Up @@ -395,6 +403,14 @@ def _create_alchemical_ions(system, charge_diff):
# The number of waters to convert is the absolute charge difference.
num_waters = abs(charge_diff)

# Make sure there are enough waters to convert. The charge difference should
# never be this large, but it prevents a crash if it is.
if num_waters > len(system["water"].molecules()):
raise ValueError(
f"Insufficient waters to convert to ions. {num_waters} required, "
f"{len(system['water'].molecules())} available."
)

# Reference coordinates.
coords = system.molecules("property is_perturbable").coordinates()
coord_string = f"{coords[0].value()}, {coords[1].value()}, {coords[2].value()}"
Expand Down

0 comments on commit bfa8b99

Please sign in to comment.