Skip to content

[DNM] Test a self self transformation with different charges #1116

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

Draft
IAlibay wants to merge 2 commits into
base: main
Choose a base branch
from
Draft

[DNM] Test a self self transformation with different charges #1116

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ac04e21
Add a test for a self self transformation with different charges
IAlibay Feb 5, 2025
0efb505
A really bad fix
IAlibay Feb 5, 2025
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
2 changes: 2 additions & 0 deletions openfe/protocols/openmm_rfe/_rfe_utils/relative.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1602,8 +1602,10 @@ def _check_indices(idx1, idx2):
# Get the parameters in the new and old systems:
old_index = hybrid_to_old_map[particle_index]
[charge_old, sigma_old, epsilon_old] = old_system_nonbonded_force.getParticleParameters(old_index)
print("core old", particle_index, [charge_old, sigma_old, epsilon_old])
new_index = hybrid_to_new_map[particle_index]
[charge_new, sigma_new, epsilon_new] = new_system_nonbonded_force.getParticleParameters(new_index)
print("core new", new_index, [charge_new, sigma_new, epsilon_new])

# Add the particle to the custom forces, interpolating between
# the two parameters; add steric params and zero electrostatics
Expand Down
39 changes: 29 additions & 10 deletions openfe/protocols/openmm_rfe/equil_rfe_methods.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -777,30 +777,49 @@
# Block out oechem backend in system_generator calls to avoid
# any issues with smiles roundtripping between rdkit and oechem
with without_oechem_backend():
system_generator = system_creation.get_system_generator(
system_generatorA = system_creation.get_system_generator(
forcefield_settings=forcefield_settings,
integrator_settings=integrator_settings,
thermo_settings=thermo_settings,
cache=ffcache,
cache=None,
has_solvent=solvent_comp is not None,
)

system_generatorB = system_creation.get_system_generator(
forcefield_settings=forcefield_settings,
integrator_settings=integrator_settings,
thermo_settings=thermo_settings,
cache=None,
has_solvent=solvent_comp is not None,
)

# c. force the creation of parameters
# This is necessary because we need to have the FF templates
# registered ahead of solvating the system.
for smc, mol in chain(off_small_mols['stateA'],
off_small_mols['stateB'],
off_small_mols['both']):
system_generator.create_system(mol.to_topology().to_openmm(),
molecules=[mol])
# stateA
system_generatorA.create_system(
off_small_mols['stateA'][0][1].to_topology().to_openmm(),
molecules=[off_small_mols['stateA'][0][1]]
)
system_generatorB.create_system(
off_small_mols['stateB'][0][1].to_topology().to_openmm(),
molecules=[off_small_mols['stateB'][0][1]]
)

for smc, mol in off_small_mols['both']:
for sys_gen in [system_generatorA, system_generatorB]:
sys_gen.create_system(

Check warning on line 811 in openfe/protocols/openmm_rfe/equil_rfe_methods.py

View check run for this annotation

Codecov / codecov/patch

openfe/protocols/openmm_rfe/equil_rfe_methods.py#L810-L811 

Added lines #L810 - L811 were not covered by tests
mol.to_topology().to_openmm(),
molecules=[mol]
)

# c. get OpenMM Modeller + a dictionary of resids for each component
stateA_modeller, comp_resids = system_creation.get_omm_modeller(
protein_comp=protein_comp,
solvent_comp=solvent_comp,
small_mols=dict(chain(off_small_mols['stateA'],
off_small_mols['both'])),
omm_forcefield=system_generator.forcefield,
omm_forcefield=system_generatorA.forcefield,
solvent_settings=solvation_settings,
)

Expand All @@ -815,7 +834,7 @@
# Block out oechem backend in system_generator calls to avoid
# any issues with smiles roundtripping between rdkit and oechem
with without_oechem_backend():
stateA_system = system_generator.create_system(
stateA_system = system_generatorA.create_system(
stateA_modeller.topology,
molecules=[m for _, m in chain(off_small_mols['stateA'],
off_small_mols['both'])],
Expand All @@ -834,7 +853,7 @@
# Block out oechem backend in system_generator calls to avoid
# any issues with smiles roundtripping between rdkit and oechem
with without_oechem_backend():
stateB_system = system_generator.create_system(
stateB_system = system_generatorB.create_system(
stateB_topology,
molecules=[m for _, m in chain(off_small_mols['stateB'],
off_small_mols['both'])],
Expand Down
46 changes: 26 additions & 20 deletions openfe/tests/protocols/test_openmm_equil_rfe_protocols.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -683,8 +683,11 @@ def test_dry_run_charge_backends(
np.testing.assert_allclose(c, ref, rtol=1e-4)


@pytest.mark.flaky(reruns=3) # bad minimisation can happen
def test_dry_run_user_charges(benzene_modifications, tmpdir):
#@pytest.mark.flaky(reruns=3) # bad minimisation can happen
@pytest.mark.parametrize('nameA, nameB', [
['benzene', 'toluene'], ['benzene', 'benzene'],
])
def test_dry_run_user_charges(benzene_modifications, nameA, nameB, tmpdir):
"""
Create a hybrid system with a set of fictitious user supplied charges
and ensure that they are properly passed through to the constructed
Expand Down Expand Up @@ -719,27 +722,29 @@ def check_propchgs(smc, charge_array):
np.testing.assert_allclose(prop_chgs, charge_array.m)

# Create new smc with overriden charges
benzene_offmol = benzene_modifications['benzene'].to_openff()
toluene_offmol = benzene_modifications['toluene'].to_openff()
benzene_rand_chg = assign_fictitious_charges(benzene_offmol)
toluene_rand_chg = assign_fictitious_charges(toluene_offmol)
benzene_offmol.partial_charges = benzene_rand_chg
toluene_offmol.partial_charges = toluene_rand_chg
benzene_smc = openfe.SmallMoleculeComponent.from_openff(benzene_offmol)
toluene_smc = openfe.SmallMoleculeComponent.from_openff(toluene_offmol)
molA_offmol = benzene_modifications[nameA].to_openff()
molB_offmol = benzene_modifications[nameB].to_openff()
molA_rand_chg = assign_fictitious_charges(molA_offmol)
molB_rand_chg = assign_fictitious_charges(molB_offmol)
molA_offmol.partial_charges = molA_rand_chg
molB_offmol.partial_charges = molB_rand_chg
molA_smc = openfe.SmallMoleculeComponent.from_openff(molA_offmol)
molB_smc = openfe.SmallMoleculeComponent.from_openff(molB_offmol)

# Check that the new smcs have the new overriden charges
check_propchgs(benzene_smc, benzene_rand_chg)
check_propchgs(toluene_smc, toluene_rand_chg)
check_propchgs(molA_smc, molA_rand_chg)
check_propchgs(molB_smc, molB_rand_chg)

# Create new mapping
mapper = openfe.setup.LomapAtomMapper(element_change=False)
mapping = next(mapper.suggest_mappings(benzene_smc, toluene_smc))
mapping = next(mapper.suggest_mappings(molA_smc, molB_smc))
print(molA_smc.to_openff().partial_charges)
print(molB_smc.to_openff().partial_charges)

# create DAG from protocol and take first (and only) work unit from within
dag = protocol.create(
stateA=openfe.ChemicalSystem({'l': benzene_smc, }),
stateB=openfe.ChemicalSystem({'l': toluene_smc, }),
stateA=openfe.ChemicalSystem({'l': molA_smc, }),
stateB=openfe.ChemicalSystem({'l': molB_smc, }),
mapping=mapping,
)
dag_unit = list(dag.protocol_units)[0]
Expand Down Expand Up @@ -792,21 +797,22 @@ def check_propchgs(smc, charge_array):
# offset (c_offsets) is equal to -(molA particle charge)
if i in htf._atom_classes['unique_old_atoms']:
idx = htf._hybrid_to_old_map[i]
np.testing.assert_allclose(c, benzene_rand_chg[idx])
np.testing.assert_allclose(c_offsets[i], -benzene_rand_chg[idx])
np.testing.assert_allclose(c, molA_rand_chg[idx])
np.testing.assert_allclose(c_offsets[i], -molA_rand_chg[idx])
# particle charge (c) is equal to 0
# offset (c_offsets) is equal to molB particle charge
elif i in htf._atom_classes['unique_new_atoms']:
idx = htf._hybrid_to_new_map[i]
np.testing.assert_allclose(c, 0 * unit.elementary_charge)
np.testing.assert_allclose(c_offsets[i], toluene_rand_chg[idx])
np.testing.assert_allclose(c_offsets[i], molB_rand_chg[idx])
# particle charge (c) is equal to molA particle charge
# offset (c_offsets) is equal to difference between molB and molA
elif i in htf._atom_classes['core_atoms']:
old_i = htf._hybrid_to_old_map[i]
new_i = htf._hybrid_to_new_map[i]
c_exp = toluene_rand_chg[new_i] - benzene_rand_chg[old_i]
np.testing.assert_allclose(c, benzene_rand_chg[old_i])
c_exp = molB_rand_chg[new_i] - molA_rand_chg[old_i]
print(c, c_exp)
np.testing.assert_allclose(c, molA_rand_chg[old_i])
np.testing.assert_allclose(c_offsets[i], c_exp)


Expand Down