LAMMPS writer can't export TIP3P water with constraints #1141
Comments
|
The constraints issue is genuine and a little surprising, but probably just not something we've happened to do yet. I'll add it to our triage process and use this issue as a stub for updates. I suspect the Omega warning you're seeing is from your call to |
|
@mattwthompson see comment in #1144, I've closed that issue following your clarification. That still does leave the LAMMPS constraints issue pretty seriously open however, can you provide some insight into why the constraint styles for TIP3P would be inconsistent? Stranger still, this Exception doesn't appear to crop up for other small molecules using a non-tip3p forcefield, but is instead replaced by an "unsupported electrostatics" error: ...: import logging
...: logging.basicConfig(level=logging.INFO)
...:
...: from openff.toolkit import Molecule, ForceField
...:
...: offmol = Molecule.from_smiles('CCO')
...: offmol.generate_conformers(n_conformers=1)
...:
...: ff = ForceField('openff-2.0.0.offxml')
...: inc = ff.create_interchange(offmol.to_topology())
...:
...: inc.to_lammps(prefix='ethanol')
INFO:rdkit:Enabling RDKit 2023.09.6 jupyter extensions
Warning: Cannot perform Hydrogen sampling with GPU-Omega: GPU-Omega disabled.
INFO:openff.toolkit.typing.engines.smirnoff.parameters:Attempting to up-convert vdW section from 0.3 to 0.4
INFO:openff.toolkit.typing.engines.smirnoff.parameters:Successfully up-converted vdW section from 0.3 to 0.4. `method="cutoff"` is now split into `periodic_method="cutoff"` and `nonperiodic_method="no-cutoff"`.
INFO:openff.toolkit.typing.engines.smirnoff.parameters:Attempting to up-convert Electrostatics section from 0.3 to 0.4
INFO:openff.toolkit.typing.engines.smirnoff.parameters:Successfully up-converted Electrostatics section from 0.3 to 0.4. `method="PME"` is now split into `periodic_potential="Ewald3D-ConductingBoundary"`, `nonperiodic_potential="Coulomb"`, and `exception_potential="Coulomb"`
------------------------------------------------------------------------------------------------------------------------------------------------------
UnsupportedExportError Traceback (most recent call last)
Cell In[1], line 12
9 ff = ForceField('openff-2.0.0.offxml')
10 inc = ff.create_interchange(offmol.to_topology())
---> 12 inc.to_lammps(prefix='ethanol')
File ~/miniconda3/envs/nrel-polymers/lib/python3.11/site-packages/openff/interchange/components/interchange.py:474, in Interchange.to_lammps(self, prefix)
472 datafile_path = prefix + ".lmp"
473 self.to_lammps_datafile(datafile_path)
--> 474 self.to_lammps_input(
475 prefix + "_pointenergy.in",
476 datafile_path,
477 )
File ~/miniconda3/envs/nrel-polymers/lib/python3.11/site-packages/openff/interchange/components/interchange.py:511, in Interchange.to_lammps_input(self, file_path, data_file)
508 data_file = Path(file_path).with_suffix(".lmp")
510 mdconfig = MDConfig.from_interchange(self)
--> 511 mdconfig.write_lammps_input(self, str(file_path), data_file=str(data_file))
File ~/miniconda3/envs/nrel-polymers/lib/python3.11/site-packages/openff/interchange/components/mdconfig.py:346, in MDConfig.write_lammps_input(self, interchange, input_file, data_file)
344 lmp.write(f"pair_style lj/cut/coul/cut {vdw_cutoff} {coul_cutoff}\n")
345 else:
--> 346 raise UnsupportedExportError(
347 f"Unsupported electrostatics method {self.coul_method}",
348 )
350 if self.mixing_rule == "lorentz-berthelot":
351 lmp.write("pair_modify mix arithmetic tail yes\n\n")
UnsupportedExportError: Unsupported electrostatics method Coulomb |
|
The topology = molecule.to_topology()
topology.box_vectors = Quantity([5, 5, 5], "nanometer")I've split that off to #1145 Back to the constraints - there's a lot that can be done with that code. It doesn't seem to be smart enough to process geometries from constrain parameters without looking up into the bond/angle parameters. I used to avoid special casing for water but that's probably unavoidable now. Some of the annotations are also wrong. TIP3P constraints are handled differently than Sage's general "let's constrain every H-to-heavy-atom bond" constraints. Sage ships with TIP3P, so you should see the same behavior with either force field. I applied this diff to make the error message more useful: diff --git a/openff/interchange/components/mdconfig.py b/openff/interchange/components/mdconfig.py
index ce4c3966..bc98e730 100644
--- a/openff/interchange/components/mdconfig.py
+++ b/openff/interchange/components/mdconfig.py
@@ -235,7 +235,8 @@ class MDConfig(_BaseModel):
if len(constraint_styles.difference({"Bonds", "Angles"})) > 0:
raise NotImplementedError(
- "Found unsupported constraints case in LAMMPS input writer.",
+ "Found unsupported constraints case in LAMMPS input writer. "
+ f"Constraint style(s) found: {constraint_styles}",
)
constrained_bond_smirks = {from openff.interchange._tests import ForceField, Topology, get_test_file_path
water_dimer = Topology.from_pdb(get_test_file_path("water-dimer.pdb"))
def try_lammps_export(force_field: str):
try:
ForceField(force_field).create_interchange(water_dimer).to_lammps('f')
except NotImplementedError as e:
print(f"failed with {force_field}: {e}")
try_lammps_export("tip3p.offxml")
# failed with tip3p.offxml: Found unsupported constraints case in LAMMPS input writer. Constraint style(s) found: {'Constraints'}
try_lammps_export("openff-2.0.0.offxml")
# failed with openff-2.0.0.offxml: Found unsupported constraints case in LAMMPS input writer. Constraint style(s) found: {'Constraints'}
# no TIP3P parameters - water treated as a small molecule - runs without error
try_lammps_export("openff-1.0.0.offxml")I think if we can get a water dimer to export sensibly that should handle most other cases. I'm a little surprised we haven't seen this yet (or I'm just forgetting it ...) |
|
@timbernat could you remind me / point me towards the relevant docs on what a LAMMPS used would expect here? The simplest case of a water molecule or two should be fine, but we will need to worry about constrained/rigid water in the same box as ligands or biopolymers with bonds between heavy atoms and hydrogens constrained. That's what's in Sage out of the box. The five "Shake X" sections here seems like overkill to my intuition? https://docs.lammps.org/Howto_tip3p.html |
mattwthompson
added
feedback needed
|
The difficulty is in assigning constraint information without bond and angle parameters, since TIP3P doesn't have harmonic bond or angle parameters. I'm not sure what LAMMPS users expect here. |
|
Looping in @DariaLazar who is an expert LAMMPS user, could you provide some insight here? |
Description
Interchange instances created from a Topology containing any number of water molecules (including just 1) and from the SMIRNOFF TIP3P force field cannot be written to LAMMPS input/data files using
Interchange.to_lammps(), "citing unsupported constraints in the LAMMPS writer".Incidentally (unless I've misunderstood what's going on), Interchange creation does NOT use TIP3P library charges supplied by the FF, but rather assigns explicit AM1-BCC charges. This is evidenced by the telltaleResolved in #1144Warning: Cannot perform Hydrogen sampling with GPU-Omega: GPU-Omega disabled.log message generated by the OpenEye charge model when callingForceField.create_interchange(...).Reproduction
Output
Software versions
conda list:The text was updated successfully, but these errors were encountered: