Skip to content

Commit

Permalink
Merge pull request #695 from openforcefield/fix-694
Browse files Browse the repository at this point in the history 
FIX: Do not assume vdW handler exists
  • Loading branch information
@mattwthompson
mattwthompson authored May 1, 2023
2 parents c7b9c6d + 5f5f720 commit 83f0677
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 114 deletions.
2 changes: 2 additions & 0 deletions devtools/conda-envs/dev_env.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,3 +55,5 @@ dependencies:
- flake8
- snakeviz
- tuna
- pip:
- git+https://github.com/jthorton/de-forcefields
174 changes: 60 additions & 114 deletions openff/interchange/interop/openmm/_nonbonded.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,104 +92,49 @@ 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 "
"any electrostatics collection. This is a supported use case but should have been caught "
"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

Expand Down Expand Up @@ -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:
Expand Down

0 comments on commit 83f0677

Please sign in to comment.