diff --git a/devtools/conda-envs/dev_env.yaml b/devtools/conda-envs/dev_env.yaml index c8faba4b6..a48edd4b1 100644 --- a/devtools/conda-envs/dev_env.yaml +++ b/devtools/conda-envs/dev_env.yaml @@ -55,3 +55,5 @@ dependencies: - flake8 - snakeviz - tuna + - pip: + - git+https://github.com/jthorton/de-forcefields diff --git a/openff/interchange/interop/openmm/_nonbonded.py b/openff/interchange/interop/openmm/_nonbonded.py index f96a18a39..7b15a1077 100644 --- a/openff/interchange/interop/openmm/_nonbonded.py +++ b/openff/interchange/interop/openmm/_nonbonded.py @@ -92,73 +92,28 @@ def _process_nonbonded_forces( ) # TODO: Process ElectrostaticsHandler.exception_potential - # TODO: Improve logic to be less reliant on these names - if "vdW" in interchange.collections or "Electrostatics" in interchange.collections: - _data = _prepare_input_data(interchange) + has_vdw = False - if combine_nonbonded_forces: - _func = _create_single_nonbonded_force - else: - _func = _create_multiple_nonbonded_forces - - _func( - _data, - interchange, - system, - molecule_virtual_site_map, - openff_openmm_particle_map, - ) - - elif "Buckingham-6" in interchange.collections: - if has_virtual_sites: - raise UnsupportedExportError( - "Virtual sites with Buckingham-6 potential not supported. If this use case is important to you, " - "please raise an issue describing the functionality you wish to see.", - ) - - buck = interchange["Buckingham-6"] - - non_bonded_force = openmm.CustomNonbondedForce( - "A * exp(-B * r) - C * r ^ -6; A = sqrt(A1 * A2); B = 2 / (1 / B1 + 1 / B2); C = sqrt(C1 * C2)", - ) - non_bonded_force.addPerParticleParameter("A") - non_bonded_force.addPerParticleParameter("B") - non_bonded_force.addPerParticleParameter("C") - system.addForce(non_bonded_force) - - for molecule in interchange.topology.molecules: - for _ in molecule.atoms: - non_bonded_force.addParticle([0.0, 0.0, 0.0]) - - if interchange.box is None: - non_bonded_force.setNonbondedMethod(openmm.NonbondedForce.NoCutoff) - else: - non_bonded_force.setNonbondedMethod(openmm.NonbondedForce.CutoffPeriodic) - non_bonded_force.setCutoffDistance(buck.cutoff * unit.angstrom) - - for top_key, pot_key in buck.key_map.items(): - atom_idx = top_key.atom_indices[0] - - # TODO: Add electrostatics - params = buck.potentials[pot_key].parameters - a = to_openmm_quantity(params["A"]) - b = to_openmm_quantity(params["B"]) - c = to_openmm_quantity(params["C"]) - non_bonded_force.setParticleParameters(atom_idx, [a, b, c]) - - openff_openmm_particle_map - - else: - # Here we assume there are no vdW interactions in any collections + for name, collection in interchange.collections.items(): + if name == "vdW": + has_vdw = True + break + if collection.is_plugin: + # TODO: Here is where to detect an electrostatics plugin, if one ever exists + if collection.acts_as == "vdW": + has_vdw = True + break + if not has_vdw: if has_virtual_sites: raise UnsupportedExportError( - "Virtual sites with no vdW handler not currently supported. If this use case is important to you, " - "please raise an issue describing the functionality you wish to see.", + "Virtual sites with no vdW handler not currently supported. If this use case is " + "important to you, please raise an issue describing the functionality you wish to " + "see.", ) try: - electrostatics = interchange["Electrostatics"] + interchange["Electrostatics"] except LookupError: raise InternalInconsistencyError( "In a confused state, could not find any vdW interactions but also failed to find " @@ -166,30 +121,20 @@ def _process_nonbonded_forces( "earlier in this function. Please file an issue with a minimal reproducing example.", ) - electrostatics_method = ( - electrostatics.periodic_potential if electrostatics else None - ) + _data = _prepare_input_data(interchange) - non_bonded_force = openmm.NonbondedForce() - system.addForce(non_bonded_force) - - for molecule in interchange.topology.molecules: - for _ in molecule.atoms: - non_bonded_force.addParticle(0.0, 1.0, 0.0) + if combine_nonbonded_forces: + _func = _create_single_nonbonded_force + else: + _func = _create_multiple_nonbonded_forces - if electrostatics_method in ["Coulomb", None]: - non_bonded_force.setNonbondedMethod(openmm.NonbondedForce.NoCutoff) - # TODO: Would setting the dispersion correction here have any impact? - elif electrostatics_method == _PME: - non_bonded_force.setNonbondedMethod(openmm.NonbondedForce.LJPME) - non_bonded_force.setEwaldErrorTolerance(1.0e-4) - else: - raise UnsupportedCutoffMethodError( - f"Found no vdW interactions but an electrostatics method of {electrostatics_method}. " - "This is either unsupported or ambiguous. If you believe this exception has been raised " - "in error, please file an issue with a minimally reproducing example and your motivation " - "for this use case.", - ) + _func( + _data, + interchange, + system, + molecule_virtual_site_map, + openff_openmm_particle_map, + ) return openff_openmm_particle_map @@ -700,46 +645,47 @@ def _create_multiple_nonbonded_forces( openmm_pairs.append(openmm_indices) - for i in range(electrostatics_force.getNumExceptions()): - (p1, p2, _, _, _) = electrostatics_force.getExceptionParameters(i) - - if (p1, p2) in openmm_pairs or (p2, p1) in openmm_pairs: - if vdw_force is not None: - if data["vdw_collection"].is_plugin: - # Since we fed in in r_min1, epsilon1, ..., r_min2, epsilon2, ... - # each as individual parameters, we need to prepare a list of - # length 2 * len(potential_parameters) to this constructor - parameters1 = vdw_force.getParticleParameters(p1) - parameters2 = vdw_force.getParticleParameters(p2) - vdw_14_force.addBond(p1, p2, [*parameters1, *parameters2]) + if electrostatics_force is not None: + for i in range(electrostatics_force.getNumExceptions()): + (p1, p2, _, _, _) = electrostatics_force.getExceptionParameters(i) + + if (p1, p2) in openmm_pairs or (p2, p1) in openmm_pairs: + if vdw_force is not None: + if data["vdw_collection"].is_plugin: + # Since we fed in in r_min1, epsilon1, ..., r_min2, epsilon2, ... + # each as individual parameters, we need to prepare a list of + # length 2 * len(potential_parameters) to this constructor + parameters1 = vdw_force.getParticleParameters(p1) + parameters2 = vdw_force.getParticleParameters(p2) + vdw_14_force.addBond(p1, p2, [*parameters1, *parameters2]) - else: - # Look up the vdW parameters for each particle - sig1, eps1 = vdw_force.getParticleParameters(p1) - sig2, eps2 = vdw_force.getParticleParameters(p2) + else: + # Look up the vdW parameters for each particle + sig1, eps1 = vdw_force.getParticleParameters(p1) + sig2, eps2 = vdw_force.getParticleParameters(p2) - # manually compute ... - sig_14 = (sig1 + sig2) * 0.5 - eps_14 = (eps1 * eps2) ** 0.5 * vdw_14 + # manually compute ... + sig_14 = (sig1 + sig2) * 0.5 + eps_14 = (eps1 * eps2) ** 0.5 * vdw_14 - # ... and set the 1-4 interactions - vdw_14_force.addBond(p1, p2, [sig_14, eps_14]) + # ... and set the 1-4 interactions + vdw_14_force.addBond(p1, p2, [sig_14, eps_14]) - # Look up the partial charges for each particle - q1 = electrostatics_force.getParticleParameters(p1)[0] - q2 = electrostatics_force.getParticleParameters(p2)[0] + # Look up the partial charges for each particle + q1 = electrostatics_force.getParticleParameters(p1)[0] + q2 = electrostatics_force.getParticleParameters(p2)[0] - # manually compute ... - qq = q1 * q2 * coul_14 + # manually compute ... + qq = q1 * q2 * coul_14 - # ... and set the 1-4 interactions - coul_14_force.addBond(p1, p2, [qq]) + # ... and set the 1-4 interactions + coul_14_force.addBond(p1, p2, [qq]) - if vdw_force is not None: - vdw_force.addExclusion(p1, p2) + if vdw_force is not None: + vdw_force.addExclusion(p1, p2) - if electrostatics_force is not None: - electrostatics_force.setExceptionParameters(i, p1, p2, 0.0, 0.0, 0.0) + if electrostatics_force is not None: + electrostatics_force.setExceptionParameters(i, p1, p2, 0.0, 0.0, 0.0) for force in [vdw_force, electrostatics_force, vdw_14_force, coul_14_force]: if force is not None: