Skip to content

Apply setAmberWater unconditionally to ensure fully rigid water constraints #163

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lohedges merged 1 commit into from
Jun 26, 2026
Merged

Apply setAmberWater unconditionally to ensure fully rigid water constraints #163

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 39 additions & 48 deletions src/somd2/runner/_base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -266,54 +266,6 @@ def __init__(self, system, config):
)
self._system = _sr.morph.link_to_reference(self._system)

# Next, swap the water topology so that it is in AMBER format.

try:
waters = self._system["water"]
except:
waters = []

if len(waters) > 0:
from sire.legacy.IO import isAmberWater as _isAmberWater
from sire.legacy.IO import setAmberWater as _setAmberWater

if not _isAmberWater(waters[0]):
num_atoms = waters[0].num_atoms()

if num_atoms == 3:
# Here we assume TIP3P for any 3-point water model.
model = "tip3p"
elif num_atoms == 4:
# Check for OPC water.
try:
if (
waters[0]
.search("element Xx")
.atoms()[0]
.charge()
.value()
< -1.1
):
model = "opc"
else:
model = "tip4p"
except:
model = "tip4p"
elif num_atoms == 5:
model = "tip5p"
try:
self._system = _System(
_setAmberWater(self._system._system, model)
)
_logger.info(
"Converting water topology to AMBER format for SOMD1 compatibility."
)
except Exception as e:
_logger.error(
"Unable to convert water topology to AMBER format for SOMD1 compatibility."
)
raise e

# Ghost atoms are considered light when adding bond constraints.
self._config._extra_args["ghosts_are_light"] = True

Expand All @@ -334,6 +286,45 @@ def __init__(self, system, config):
_logger.error(msg)
raise RuntimeError(msg)

# Convert water topology to AMBER format if not already done. AMBER
# format adds an explicit H-H bond, giving fully rigid water (O-H and
# H-H constraints) rather than just O-H constraints under h_bonds.
try:
waters = self._system["water"]
except:
waters = []

if len(waters) > 0:
from sire.legacy.IO import isAmberWater as _isAmberWater
from sire.legacy.IO import setAmberWater as _setAmberWater

if not _isAmberWater(waters[0]):
num_atoms = waters[0].num_atoms()

if num_atoms == 3:
model = "tip3p"
elif num_atoms == 4:
try:
if (
waters[0].search("element Xx").atoms()[0].charge().value()
< -1.1
):
model = "opc"
else:
model = "tip4p"
except:
model = "tip4p"
elif num_atoms == 5:
model = "tip5p"
try:
self._system = _System(_setAmberWater(self._system._system, model))
_logger.info(
f"Converting water topology to AMBER {model.upper()} format."
)
except Exception as e:
_logger.error("Unable to convert water topology to AMBER format.")
raise e

# Check the end state constraints.
self._check_end_state_constraints()

Expand Down
Loading