diff --git a/openfe/protocols/openmm_rfe/_rfe_utils/relative.py b/openfe/protocols/openmm_rfe/_rfe_utils/relative.py index db4002ace..2e151e945 100644 --- a/openfe/protocols/openmm_rfe/_rfe_utils/relative.py +++ b/openfe/protocols/openmm_rfe/_rfe_utils/relative.py @@ -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 diff --git a/openfe/protocols/openmm_rfe/equil_rfe_methods.py b/openfe/protocols/openmm_rfe/equil_rfe_methods.py index d2f95f2c3..daea4d3f9 100644 --- a/openfe/protocols/openmm_rfe/equil_rfe_methods.py +++ b/openfe/protocols/openmm_rfe/equil_rfe_methods.py @@ -777,22 +777,41 @@ def run(self, *, dry=False, verbose=True, # 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( + 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( @@ -800,7 +819,7 @@ def run(self, *, dry=False, verbose=True, 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, ) @@ -815,7 +834,7 @@ def run(self, *, dry=False, verbose=True, # 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'])], @@ -834,7 +853,7 @@ def run(self, *, dry=False, verbose=True, # 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'])], diff --git a/openfe/tests/protocols/test_openmm_equil_rfe_protocols.py b/openfe/tests/protocols/test_openmm_equil_rfe_protocols.py index 513962355..b9efd8640 100644 --- a/openfe/tests/protocols/test_openmm_equil_rfe_protocols.py +++ b/openfe/tests/protocols/test_openmm_equil_rfe_protocols.py @@ -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 @@ -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] @@ -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)