Skip to content
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.

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

Add support for Foyer forcefields #517

Merged
merged 34 commits into from Jul 24, 2023
Merged

Add support for Foyer forcefields #517

merged 34 commits into from Jul 24, 2023

Conversation

aehogan
Copy link
Contributor

@aehogan aehogan commented Jul 17, 2023

Description

Add support for Foyer forcefields

Todos

  • Add ForceFieldSource and build system protocol
  • Add tests
  • Debugging

Questions

  • When running a test enthalpy of vaporization calculation I get NotImplementedError() in assign_parameters_gas

Status

  • Ready to go

@mattwthompson
Copy link
Member

Wow, this is a nice surprise!

Make sure you add foyer into the test environment.

Don't worry about the OpenEye parts of the matrix, but you can force the tests to run in those cases by adding if: always() to that step in the workflow.

@codecov
Copy link

codecov bot commented Jul 18, 2023

Codecov Report

Merging #517 (de3e374) into main (c6fb943) will decrease coverage by 2.48%.
The diff coverage is 93.82%.

Additional details and impacted files

@aehogan
Copy link
Contributor Author

aehogan commented Jul 18, 2023

@mattwthompson with the current code the tests run fine

assign_parameters = BuildFoyerSystem("assign_parameters")
assign_parameters.force_field_path = force_field_source_path
assign_parameters.coordinate_file_path = build_coordinates.coordinate_file_path
assign_parameters.substance = substance
assign_parameters.execute(directory)

but when I run a delta H vaporization simulation using my production script I get

a19a99f346ff46e3819f8fb530c67ba5|assign_parameters_gas failed to execute.

Traceback (most recent call last):
File "[...]/lib/python3.9/site-packages/openff/evaluator/workflow/protocols.py", line 1193, in _execute_protocol
protocol.execute(directory, available_resources)
File "[...]/lib/python3.9/site-packages/openff/evaluator/workflow/protocols.py", line 680, in execute
self._execute(directory, available_resources)
File "[...]/lib/python3.9/site-packages/openff/evaluator/protocols/forcefield.py", line 246, in _execute
raise NotImplementedError()
NotImplementedError

Would you know off hand why that is?

@mattwthompson
Copy link
Member

The traceback looks like it's not coming from a development/editable install - does from openff.evaluator import __file__; print(__file__) point to your clone of this repo?

@aehogan
Copy link
Contributor Author

aehogan commented Jul 18, 2023

Yes, I pip installed it to test it with the production script

openff/evaluator/forcefield/__init__.py Show resolved Hide resolved
openff/evaluator/forcefield/forcefield.py Show resolved Hide resolved
openff/evaluator/forcefield/forcefield.py Outdated Show resolved Hide resolved
openff/evaluator/protocols/forcefield.py Outdated Show resolved Hide resolved
openff/evaluator/protocols/forcefield.py Outdated Show resolved Hide resolved
openff/evaluator/protocols/forcefield.py Outdated Show resolved Hide resolved
openff/evaluator/tests/test_protocols/test_forcefield.py Outdated Show resolved Hide resolved
@mattwthompson
Copy link
Member

The error I'm seeing in tests:

openff.interchange.exceptions.UnsupportedCutoffMethodError: When using openmm.CustomNonbondedForce, vdW and electrostatics cutoff methods must agree on whether or not periodic boundary conditions should be used.OpenMM will throw an error. Found vdw method 1, and electrostatics method 0,

usually stems from the periodicity not being set. If Topology.box_vectors is set, which I assume it would be after the PACKMOL protocol finishes, then Interchange.box should also be set. Unless the error bubbles out from a system that is not periodic, in which case you're running up against a nasty confluence of corner cases in which the correct behavior is undefined (openforcefield/standards#7). We're hoping to resolve this within a month or so (we have a meeting in two weeks about it) but as a short-term stop-gap, you could just set the gas-phase calculations to have a large (say, 5 nm or greater) periodic box. If you choose this option please note it in a comment in the code so it can be reverted later.

Comment on lines 212 to 217
<HarmonicBondForce>
<Bond class1="OW_tip3p" class2="HW_tip3p" length="0.09572" k="1884.06"/>
</HarmonicBondForce>
<HarmonicAngleForce>
<Angle class1="HW_tip3p" class2="OW_tip3p" class3="HW_tip3p" angle="1.82421813" k="230.274"/>
</HarmonicAngleForce>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

TIP3P is supposed to be a fixed model, but I don't think Foyer support constraints in its force field schema? We always did that at simulation runtime, I think

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The test builds a rigid tip3p model HOWEVER this seems due to a short circuit in TemplateBuildSystem._execute that forces all water to be rigid tip3p with TemplateBuildSystem._build_tip3p_system. Is the intent to always override any water with tip3p in evaluator?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think that's the intent

Comment on lines 1101 to 1106
interchange = Interchange.from_foyer(topology=topology, force_field=force_field)

openmm_pdb_file = app.PDBFile(self.coordinate_file_path)
if openmm_pdb_file.topology.getPeriodicBoxVectors() is not None:
interchange.box = openmm_pdb_file.topology.getPeriodicBoxVectors()

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this mean Interchange.from_foyer isn't parsing the box vectors on the topology argument? If so, I think that's a bug

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

topology is not an argument to _parameterize_molecule, just the molecule itself which contains no box information, the better fix would be to pass the topology to _parameterize_molecule but I'm trying to make the smallest changes I can

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the topology being passed to interchange is from molecule.to_topology() not from TemplateBuildSystem._execute.topology (which isn't passed as an argument to _parameterize_molecule)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That makes sense, I was misreading the code then

@mattwthompson
Copy link
Member

@aehogan I think this is ready to go? Is there anything more you were looking to do, and have you been able to use this in your studies? If there aren't issues you've ran into, I think I'll give this another review and look to merge it into a release soon. How does that sound to you?

The changeset is smaller than I originally anticipated, but thinking through it again, it shouldn't be a surprise that it doesn't need to do much more than swap out the force field source.

@aehogan
Copy link
Contributor Author

aehogan commented Jul 20, 2023

@mattwthompson There are still one or two more issues with running a full delta H vaporization simulation I need to iron out.

Besides that I was going to change out the xml test with a methane so that part of the code is actually tested.

@mattwthompson
Copy link
Member

Sounds good, let me know when you think it's ready - no rush

@aehogan
Copy link
Contributor Author

aehogan commented Jul 24, 2023

@mattwthompson should be good now

Copy link
Member

@mattwthompson mattwthompson left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I once again trust this is working for simulations in your research; the code looks great but it's hard to say much comprehensively without battle-testing it since Evaluator relies so much on things that take a long time to test.

Could you also merge the main branch back in add a one-liner to the "Current development" section? Copying the name of the PR is plenty.

@mattwthompson
Copy link
Member

Awesome work!

@mattwthompson mattwthompson merged commit f6bdd5b into openforcefield:main Jul 24, 2023
8 of 15 checks passed
@aehogan
Copy link
Contributor Author

aehogan commented Jul 24, 2023

Thanks! Delta H vaporization, delta H mixing and densities have all run successfully now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants