fix solvent atom index shifting

This commit is contained in:
Josh Horton
2026-01-21 11:06:25 +00:00
parent cbac853ee7
commit 36b623d918
5 changed files with 305 additions and 5 deletions

View File

@@ -755,6 +755,7 @@ def _scale_angles_and_torsions(htf: HybridTopologyFactory, scale_factor: float =
softened_hybrid_system = openmm.System()
# add all the particles
logger.info("Copying particles and constraints to new hybrid system.")
print("Copying particles and constraints to new hybrid system.")
for i in range(htf.hybrid_system.getNumParticles()):
softened_hybrid_system.addParticle(htf.hybrid_system.getParticleMass(i))
# add all constraints
@@ -799,7 +800,9 @@ def _scale_angles_and_torsions(htf: HybridTopologyFactory, scale_factor: float =
softened_harmonic_angle_force = openmm.HarmonicAngleForce()
softened_hybrid_system.addForce(softened_harmonic_angle_force)
logger.info("Processing dummy-core junction angles for softening.")
print("Processing dummy-core junction angles for softening.")
logger.info("Adding softened angles to core_angle_force.")
print("Adding softened angles to core_angle_force.")
default_hybrid_angle_force = hybrid_forces["standard_angle_force"]
softened_custom_angle_force = new_hybrid_forces["CustomAngleForce"]
for i in range(default_hybrid_angle_force.getNumAngles()):
@@ -812,12 +815,14 @@ def _scale_angles_and_torsions(htf: HybridTopologyFactory, scale_factor: float =
# add the term to the interpolated custom angle force
new_k = k * scale_factor
logger.info(f"Softening angle {angle} at lambda=0: original k = {k}, new k = {new_k}")
print(f"Softening angle {angle} at lambda=0: original k = {k}, new k = {new_k}")
softened_custom_angle_force.addAngle(p1, p2, p3, [theta_eq, new_k, theta_eq, k])
elif 1 <= len(dummy_old_atoms.intersection(angle)) < 3:
# if we match an old unique atom the angle must be softened at lambda = 1
# add the term to the interpolated custom angle force
new_k = k * scale_factor
logger.info(f"Softening angle {angle} at lambda=1: original k = {k}, new k = {new_k}")
print(f"Softening angle {angle} at lambda=1: original k = {k}, new k = {new_k}")
softened_custom_angle_force.addAngle(p1, p2, p3, [theta_eq, k, theta_eq, new_k])
else:
# the term does not involve any dummy atoms, so we can just copy it
@@ -825,7 +830,9 @@ def _scale_angles_and_torsions(htf: HybridTopologyFactory, scale_factor: float =
# process torsions
logger.info("Processing dummy-core junction torsions for softening.")
print("Processing dummy-core junction torsions for softening.")
logger.info("Adding softened torsions to core_torsion_force.")
print("Adding softened torsions to core_torsion_force.")
default_hybrid_torsion_force = hybrid_forces["unique_atom_torsion_force"]
softened_custom_torsion_force = new_hybrid_forces["CustomTorsionForce"]
for i in range(default_hybrid_torsion_force.getNumTorsions()):
@@ -838,6 +845,7 @@ def _scale_angles_and_torsions(htf: HybridTopologyFactory, scale_factor: float =
# add the term to the interpolated custom torsion force
new_k = k * scale_factor
logger.info(f"Softening torsion {torsion} at lambda=0: original k = {k}, new k = {new_k}")
print(f"Softening torsion {torsion} at lambda=0: original k = {k}, new k = {new_k}")
softened_custom_torsion_force.addTorsion(p1, p2, p3, p4,
[periodicity, phase,
new_k, periodicity,
@@ -847,6 +855,7 @@ def _scale_angles_and_torsions(htf: HybridTopologyFactory, scale_factor: float =
# add the term to the interpolated custom torsion force
new_k = k * scale_factor
logger.info(f"Softening torsion {torsion} at lambda=1: original k = {k}, new k = {new_k}")
print(f"Softening torsion {torsion} at lambda=1: original k = {k}, new k = {new_k}")
softened_custom_torsion_force.addTorsion(p1, p2, p3, p4,
[periodicity, phase,
k, periodicity,
@@ -885,6 +894,10 @@ def _load_ghostly_corrections(ghostly_output_string: str) -> dict:
def _shift_ghostly_correction_indices(htf: HybridTopologyFactory, corrections: dict) -> dict:
"""Shift the atom indices in the ghostly corrections to account for the solvent environment in the HTF."""
core_atoms = htf._atom_classes["core_atoms"]
unique_old_atoms = htf._atom_classes["unique_old_atoms"]
# get the maximum index of the core and unique old atoms all stateB atoms are after this
max_state_a_atom = max(core_atoms.union(unique_old_atoms))
# we need the shift the indices of the atoms in the corrections to account for the solvent system
# this involves shifting the unique new atoms to be after the last water atom
@@ -892,7 +905,11 @@ def _shift_ghostly_correction_indices(htf: HybridTopologyFactory, corrections: d
if len(htf._atom_classes["environment_atoms"]) > 0:
print("shifting ghostly correction atom indices to account for solvent environment.")
last_water_atom = max(htf._atom_classes["environment_atoms"])
print("last water atom", last_water_atom)
last_core_atom = max(htf._atom_classes["core_atoms"])
print("last core atom", last_core_atom)
print("old unique atoms", htf._atom_classes["unique_old_atoms"])
print("new unique atoms", htf._atom_classes["unique_new_atoms"])
shifted_corrections = {}
for lambda_key in corrections.keys():
for correction_type in corrections[lambda_key].keys():
@@ -900,10 +917,10 @@ def _shift_ghostly_correction_indices(htf: HybridTopologyFactory, corrections: d
shifted_list = []
for angle in corrections[lambda_key][correction_type]:
shifted_angle = tuple(
atom_idx + (last_water_atom - last_core_atom)
atom_idx + (last_water_atom - max_state_a_atom)
# only shift if not a core atom and greater than last core atom
# else this is a unique old atom and should not be shifted
if atom_idx not in core_atoms and atom_idx > last_core_atom else atom_idx
if atom_idx not in core_atoms and atom_idx > max_state_a_atom else atom_idx
for atom_idx in angle
)
shifted_list.append(shifted_angle)
@@ -912,8 +929,8 @@ def _shift_ghostly_correction_indices(htf: HybridTopologyFactory, corrections: d
new_dict = {}
for angle, params in corrections[lambda_key][correction_type].items():
shifted_angle = tuple(
atom_idx + (last_water_atom - last_core_atom)
if atom_idx not in core_atoms and atom_idx > last_core_atom else atom_idx
atom_idx + (last_water_atom - max_state_a_atom)
if atom_idx not in core_atoms and atom_idx > max_state_a_atom else atom_idx
for atom_idx in angle
)
new_dict[shifted_angle] = params
@@ -974,6 +991,7 @@ def _apply_ghostly_corrections(htf: HybridTopologyFactory, corrections: dict) ->
p1, p2, p3, theta_eq, k = old_hybrid_angle_force.getAngleParameters(i)
# check if we have one ghost atom for this angle
angle = (p1, p2, p3)
print(angle)
if 1<= len(dummy_old_atoms.intersection(angle)) < 3 or 1<= len(dummy_new_atoms.intersection(angle)) < 3:
angle_reversed = (p3, p2, p1)
# set up containers for the end state values

View File

@@ -0,0 +1,92 @@
ejm_46
RDKit 3D
36 38 0 0 0 0 0 0 0 0999 V2000
-3.2903 22.5712 -29.5240 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.9192 21.7173 -28.5965 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.2025 21.1503 -27.5518 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.2582 25.0245 -31.6039 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.1053 23.0704 -30.6507 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.9452 24.4099 -32.6722 C 0 0 0 0 0 0 0 0 0 0 0 0
-6.7996 25.1874 -33.4514 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.9269 22.8784 -29.3601 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.8593 21.4510 -27.4045 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.2252 22.3093 -28.3025 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.6128 21.4071 -28.7312 Cl 0 0 0 0 0 0 0 0 0 0 0 0
-4.4756 22.2608 -31.5073 O 0 0 0 0 0 0 0 0 0 0 0 0
-4.4499 24.3695 -30.6718 N 0 0 0 0 0 0 0 0 0 0 0 0
-6.9355 26.4911 -33.2090 N 0 0 0 0 0 0 0 0 0 0 0 0
-6.3008 27.1207 -32.2194 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.4625 26.3854 -31.3778 C 0 0 0 0 0 0 0 0 0 0 0 0
-6.4935 28.5024 -32.0877 N 0 0 0 0 0 0 0 0 0 0 0 0
-5.7410 29.3927 -31.4053 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.7130 29.1087 -30.7881 O 0 0 0 0 0 0 0 0 0 0 0 0
-6.2261 30.8526 -31.4064 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.9465 31.6762 -30.1512 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.1392 31.9250 -31.3944 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.0986 23.9892 -30.4092 Cl 0 0 0 0 0 0 0 0 0 0 0 0
-3.6932 20.4911 -26.8557 H 0 0 0 0 0 0 0 0 0 0 0 0
-5.8817 23.3629 -32.9126 H 0 0 0 0 0 0 0 0 0 0 0 0
-7.3630 24.7668 -34.2730 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.2990 21.0136 -26.5902 H 0 0 0 0 0 0 0 0 0 0 0 0
-0.1903 22.5719 -28.1610 H 0 0 0 0 0 0 0 0 0 0 0 0
-4.0744 24.9351 -29.9290 H 0 0 0 0 0 0 0 0 0 0 0 0
-4.9676 26.8422 -30.5384 H 0 0 0 0 0 0 0 0 0 0 0 0
-7.3482 28.8436 -32.5018 H 0 0 0 0 0 0 0 0 0 0 0 0
-7.1314 31.0429 -31.9814 H 0 0 0 0 0 0 0 0 0 0 0 0
-5.4329 31.2354 -29.2953 H 0 0 0 0 0 0 0 0 0 0 0 0
-6.6880 32.4155 -29.8451 H 0 0 0 0 0 0 0 0 0 0 0 0
-5.3071 32.8596 -31.9267 H 0 0 0 0 0 0 0 0 0 0 0 0
-4.0875 31.6318 -31.4083 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 2 0
1 5 1 0
1 8 1 0
2 3 1 0
2 11 1 0
3 9 2 0
3 24 1 0
4 6 2 0
4 13 1 0
4 16 1 0
5 12 2 0
5 13 1 0
6 7 1 0
6 25 1 0
7 14 2 0
7 26 1 0
8 10 2 0
8 23 1 0
9 10 1 0
9 27 1 0
10 28 1 0
13 29 1 0
14 15 1 0
15 16 2 0
15 17 1 0
16 30 1 0
17 18 1 0
17 31 1 0
18 19 2 0
18 20 1 0
20 21 1 0
20 22 1 0
20 32 1 0
21 22 1 0
21 33 1 0
21 34 1 0
22 35 1 0
22 36 1 0
M END
> <atom.dprop.PartialCharge> (1)
-0.14360000000000001 0.052400000000000009 -0.129 0.12659999999999999 0.69169999999999998 -0.30330000000000001 0.44419999999999998 0.052400000000000009 -0.095999999999999988 -0.129
-0.059899999999999995 -0.55010000000000003 -0.46010000000000001 -0.73599999999999999 0.55020000000000002 -0.32029999999999997 -0.5474 0.70009999999999994 -0.59309999999999996
-0.22969999999999999 -0.10389999999999999 -0.10389999999999999 -0.059899999999999995 0.157 0.183 0.025100000000000011 0.14799999999999999 0.157 0.33050000000000002 0.17799999999999999
0.33350000000000002 0.09870000000000001 0.084200000000000011 0.084200000000000011 0.084200000000000011 0.084200000000000011
> <ofe-name> (1)
ejm_46
> <r_exp_dg> (1)
-11.31
$$$$

View File

@@ -0,0 +1,92 @@
jmc_27
RDKit 3D
36 38 0 0 0 0 0 0 0 0999 V2000
-3.2903 22.5712 -29.5240 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.9192 21.7173 -28.5965 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.2025 21.1503 -27.5518 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.2582 25.0245 -31.6039 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.1053 23.0704 -30.6507 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.9452 24.4099 -32.6722 C 0 0 0 0 0 0 0 0 0 0 0 0
-6.7996 25.1874 -33.4514 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.9269 22.8784 -29.3601 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.8593 21.4510 -27.4045 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.2252 22.3093 -28.3025 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.6128 21.4071 -28.7312 Cl 0 0 0 0 0 0 0 0 0 0 0 0
-4.4756 22.2608 -31.5073 O 0 0 0 0 0 0 0 0 0 0 0 0
-4.4499 24.3695 -30.6718 N 0 0 0 0 0 0 0 0 0 0 0 0
-6.9355 26.4911 -33.2090 N 0 0 0 0 0 0 0 0 0 0 0 0
-6.3008 27.1207 -32.2194 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.4625 26.3854 -31.3778 C 0 0 0 0 0 0 0 0 0 0 0 0
-6.5061 28.4999 -32.0789 N 0 0 0 0 0 0 0 0 0 0 0 0
-5.7409 29.4041 -31.4257 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.6551 29.1451 -30.9087 O 0 0 0 0 0 0 0 0 0 0 0 0
-6.3480 30.8113 -31.2740 C 0 0 1 0 0 0 0 0 0 0 0 0
-5.9277 31.5796 -30.0153 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.4977 32.0811 -31.3560 C 0 0 1 0 0 0 0 0 0 0 0 0
-1.0986 23.9892 -30.4092 Cl 0 0 0 0 0 0 0 0 0 0 0 0
-3.6931 20.4910 -26.8557 H 0 0 0 0 0 0 0 0 0 0 0 0
-5.8816 23.3628 -32.9125 H 0 0 0 0 0 0 0 0 0 0 0 0
-7.3629 24.7669 -34.2730 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.2990 21.0136 -26.5902 H 0 0 0 0 0 0 0 0 0 0 0 0
-0.1903 22.5719 -28.1611 H 0 0 0 0 0 0 0 0 0 0 0 0
-4.0740 24.9351 -29.9293 H 0 0 0 0 0 0 0 0 0 0 0 0
-4.9684 26.8537 -30.5430 H 0 0 0 0 0 0 0 0 0 0 0 0
-7.4039 28.8141 -32.4185 H 0 0 0 0 0 0 0 0 0 0 0 0
-7.3665 30.8855 -31.6383 H 0 0 0 0 0 0 0 0 0 0 0 0
-5.1588 31.1817 -29.3505 H 0 0 0 0 0 0 0 0 0 0 0 0
-6.6706 32.1636 -29.4737 H 0 0 0 0 0 0 0 0 0 0 0 0
-5.9673 33.0048 -31.6947 H 0 0 0 0 0 0 0 0 0 0 0 0
-3.7666 32.1420 -31.7726 Cl 0 0 0 0 0 0 0 0 0 0 0 0
1 2 2 0
1 5 1 0
1 8 1 0
2 3 1 0
2 11 1 0
3 9 2 0
3 24 1 0
4 6 2 0
4 13 1 0
4 16 1 0
5 12 2 0
5 13 1 0
6 7 1 0
6 25 1 0
7 14 2 0
7 26 1 0
8 10 2 0
8 23 1 0
9 10 1 0
9 27 1 0
10 28 1 0
13 29 1 0
14 15 1 0
15 16 2 0
15 17 1 0
16 30 1 0
17 18 1 0
17 31 1 0
18 19 2 0
18 20 1 0
20 21 1 0
20 22 1 0
20 32 1 6
21 22 1 0
21 33 1 0
21 34 1 0
22 35 1 6
22 36 1 0
M END
> <atom.dprop.PartialCharge> (1)
-0.14357222222222224 0.052427777777777768 -0.12897222222222224 0.12662777777777776 0.69172777777777772 -0.30227222222222222 0.44422777777777778 0.052427777777777768 -0.095972222222222237
-0.12897222222222224 -0.06037222222222223 -0.5500722222222223 -0.46107222222222222 -0.73097222222222225 0.54822777777777776 -0.31927222222222218 -0.53937222222222225 0.70012777777777768
-0.58207222222222221 -0.21167222222222223 -0.092372222222222231 -0.02427222222222223 -0.06037222222222223 0.15702777777777777 0.18302777777777776 0.026127777777777771 0.14802777777777776
0.15702777777777777 0.33052777777777781 0.17702777777777776 0.33852777777777782 0.11472777777777776 0.095727777777777759 0.095727777777777759 0.10972777777777777 -0.11737222222222224
> <ofe-name> (1)
jmc_27
> <r_exp_dg> (1)
-11.279999999999999
$$$$

View File

@@ -440,4 +440,76 @@ def chloroethane_to_ethane_mapping_ghostly(chloroethane_to_ethane_mapping, chlor
componentA=chloro_ghostly,
componentB=chloroethane_to_ethane_mapping.componentB,
componentA_to_componentB=chloroethane_to_ethane_mapping.componentA_to_componentB
)
)
@pytest.fixture(scope="module")
def tyk2_ghostly_modifications():
return {
"lambda_0": {
"removed_angles": [],
"removed_dihedrals": [
"33,20,21,36",
"19,20,21,36",
"17,19,21,36",
"31,19,21,36",
"20,19,21,36",
"32,20,21,36"
],
"stiffened_angles": [],
"softened_angles": {
"19,21,36": {
"k": 5.0,
"theta0": 2.222322474687164
},
"34,21,36": {
"k": 5.0,
"theta0": 2.000809911756968
},
"20,21,36": {
"k": 5.0,
"theta0": 1.6181044310934367
}
}
},
"lambda_1": {
"removed_angles": [],
"removed_dihedrals": [
"33,20,21,35",
"19,20,21,35",
"17,19,21,35",
"31,19,21,35",
"20,19,21,35",
"32,20,21,35"
],
"stiffened_angles": [],
"softened_angles": {
"19,21,35": {
"k": 5.0,
"theta0": 2.2915349600475916
},
"34,21,35": {
"k": 5.0,
"theta0": 1.6334317078992815
},
"20,21,35": {
"k": 5.0,
"theta0": 2.395734337620435
}
}
}
}
@pytest.fixture(scope="module")
def tyk2_ghostly_mapping(tyk2_ghostly_modifications):
"""Return a mapping from tyk2 ligand ejm_31 to ejm_32 with ghosted corrections stored on componentA."""
with resources.as_file(resources.files("openfe.tests.data.htf")) as f:
jmc_27 = openfe.SmallMoleculeComponent.from_sdf_file(f / "jmc_27.sdf")
ejm_46 = openfe.SmallMoleculeComponent.from_sdf_file(f / "ejm_46.sdf")
from kartograf import KartografAtomMapper
# we need to use the exact mappings we use in the hybrid topology so only map Hs to Hs
mapper = KartografAtomMapper(map_hydrogens_on_hydrogens_only=True)
tyk2_ghostly = jmc_27.copy_with_replacements(**{"molprops": {"ghostly_corrections": json.dumps(tyk2_ghostly_modifications)}})
return next(mapper.suggest_mappings(tyk2_ghostly, ejm_46))

View File

@@ -797,6 +797,32 @@ def test_apply_ghostly_corrections_torsions_chloro(chloroethane_to_ethane_mappin
assert k == ethane_torsion.k[term_index].m_as(unit.kilojoule_per_mole) * omm_unit.kilojoule_per_mole
def test_tyk2_ghostly_corrections(tyk2_ghostly_mapping, solv_settings, tmpdir, vac_settings):
"""Test that the Tyk2 ligand transformation with ghostly corrections runs without error and applies corrections."""
state_a_dict = {"ligand": tyk2_ghostly_mapping.componentA}
state_b_dict = {"ligand": tyk2_ghostly_mapping.componentB}
# settings = vac_settings
settings = solv_settings
state_a_dict["solvent"] = openfe.SolventComponent()
state_b_dict["solvent"] = openfe.SolventComponent()
protocol = openmm_rfe.RelativeHybridTopologyProtocol(settings=settings)
# create DAG from protocol and take first (and only) work unit from within
dag = protocol.create(
stateA=openfe.ChemicalSystem(state_a_dict),
stateB=openfe.ChemicalSystem(state_b_dict),
mapping=tyk2_ghostly_mapping,
)
dag_unit = list(dag.protocol_units)[0]
with tmpdir.as_cwd():
sampler = dag_unit.run(dry=True)["debug"]["sampler"]
htf = sampler._hybrid_factory
def test_dry_run_gaff_vacuum(
benzene_vacuum_system, toluene_vacuum_system, benzene_to_toluene_mapping, vac_settings, tmpdir
):