Better support virtual sites plugins with vdW interactions #71
Conversation
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #71 +/- ##
==========================================
+ Coverage 84.28% 85.34% +1.05%
==========================================
Files 4 6 +2
Lines 541 614 +73
==========================================
+ Hits 456 524 +68
- Misses 85 90 +5
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. |
|
@jthorton could you possibly share some scripts you've used to test what you were working on? I think I'd be guessing both what kinds of virtual sites you want to test (maybe as simple as water?) but more importantly I'd probably guess the wrong numerical parameters, or pull something from the release force field |
|
Thanks for coordinating all of this multi-branch fun Matt! I haven't managed to do any fits with these sites yet so I don't have any real values yet but I planned to use our evaluate water at distances test and move the oxygen parameters onto a v-site and check the energies are the same. Do we need to test the v-site exceptions again now they are handled by a plugin? |
|
Hmm, okay, that seems simple enough. Doesn't need to be optimized - just have a 4-site water model with the vdW interactions moved to the virtual site. What charges would you use for this models? It'd be useful, if only for testing purposes, to have something like a TIP4P variant with the charge moved onto the virtual site.
Yes - the number of corner cases this sort of thing brings out means any test we can think of should be added. It wasn't that long ago I thought virtual sites would be tractable because "oh well, at least it's just the geometry and charges" 🤣 |
You could use the water model from the dexp force field it has a vsite with charge. |
|
An upcoming 0.16.0 should include upstream changes that make this easier! |
|
This'll take a bit more work, right now the virtual sites aren't getting passed through to the |
|
Okay, I was doing some goofy stuff, here's a better (if artificial) water with the double exponential interaction shifted to the virtual site, overly golfed. The actual double exponential parameters aren't correctly set (somewhere in the OpenMM code) but I think the internal storage is okay. IIRC the 12-6 parameters are stored in a vdW handler, so this is consistent behavior. In [35]: y = ForceField(
...: "smirnoff_plugins/_tests/data/de-vs-1.0.1.offxml", load_plugins=True
...: ).create_interchange(Molecule.from_smiles("O").to_topology())
In [36]: y = ForceField(
...: "smirnoff_plugins/_tests/data/de-vs-1.0.1.offxml",
...: load_plugins=True,
...: ).create_interchange(Molecule.from_smiles("O").to_topology())
In [37]: y["DoubleExponential"].dict()
Out[37]:
{'type': 'DoubleExponential',
'is_plugin': True,
'expression': 'CombinedEpsilon*RepulsionFactor*RepulsionExp-CombinedEpsilon*AttractionFactor*AttractionExp;CombinedEpsilon=epsilon1*epsilon2;RepulsionExp=exp(-alpha*ExpDistance);AttractionExp=exp(-beta*ExpDistance);ExpDistance=r/CombinedR;CombinedR=r_min1+r_min2;',
'key_map': {TopologyKey with atom indices (0,): {'id': '[#1]-[#8X2H2+0:1]-[#1]',
'mult': None,
'associated_handler': 'DoubleExponential',
'bond_order': None,
'virtual_site_type': None,
'cosmetic_attributes': {}},
TopologyKey with atom indices (1,): {'id': '[#1:1]-[#8X2H2+0]-[#1]',
'mult': None,
'associated_handler': 'DoubleExponential',
'bond_order': None,
'virtual_site_type': None,
'cosmetic_attributes': {}},
TopologyKey with atom indices (2,): {'id': '[#1:1]-[#8X2H2+0]-[#1]',
'mult': None,
'associated_handler': 'DoubleExponential',
'bond_order': None,
'virtual_site_type': None,
'cosmetic_attributes': {}},
VirtualSiteKey with atom indices None: {'id': '[#1:2]-[#8X2H2+0:1]-[#1:3] EP once',
'mult': None,
'associated_handler': 'vdw',
'bond_order': None,
'virtual_site_type': None,
'cosmetic_attributes': {}}},
'potentials': {PotentialKey associated with handler 'DoubleExponential' with id '[#1]-[#8X2H2+0:1]-[#1]': {'parameters': {'r_min': 1 <Unit('angstrom')>,
'epsilon': 0.0 <Unit('kilocalorie / mole')>},
'map_key': None},
PotentialKey associated with handler 'DoubleExponential' with id '[#1:1]-[#8X2H2+0]-[#1]': {'parameters': {'r_min': 1 <Unit('angstrom')>,
'epsilon': 0.0 <Unit('kilocalorie / mole')>},
'map_key': None},
PotentialKey associated with handler 'vdw' with id '[#1:2]-[#8X2H2+0:1]-[#1:3] EP once': {'parameters': {'r_min': 3.538259263094 <Unit('angstrom')>,
'epsilon': 0.2130260271865 <Unit('kilocalorie / mole')>},
'map_key': None}},
'cutoff': 9.0 <Unit('angstrom')>,
'scale_12': 0.0,
'scale_13': 0.0,
'scale_14': 0.5,
'scale_15': 1.0,
'acts_as': 'vdW',
'periodic_method': 'cutoff',
'nonperiodic_method': 'no-cutoff',
'mixing_rule': '',
'switch_width': 1.0 <Unit('angstrom')>,
'alpha': 16.76632136164 <Unit('dimensionless')>,
'beta': 4.426755081718 <Unit('dimensionless')>}
In [38]: y["VirtualSites"].dict()
Out[38]:
{'type': 'VirtualSites',
'is_plugin': True,
'expression': '',
'key_map': {VirtualSiteKey with atom indices None: {'id': '[#1:2]-[#8X2H2+0:1]-[#1:3] EP once',
'mult': None,
'associated_handler': 'VirtualSites',
'bond_order': None,
'virtual_site_type': 'DivalentLonePair',
'cosmetic_attributes': {}}},
'potentials': {PotentialKey associated with handler 'VirtualSites' with id '[#1:2]-[#8X2H2+0:1]-[#1:3] EP once': {'parameters': {'distance': -0.010743 <Unit('nanometer')>,
'outOfPlaneAngle': 0.0 <Unit('degree')>,
'inPlaneAngle': None},
'map_key': None}},
'virtual_site_key_topology_index_map': {VirtualSiteKey with atom indices None: 3},
'exclusion_policy': 'parents',
'acts_as': 'VirtualSites'} |
|
I think this is another representation of the same problem; most stuff is in place, but I don't think the virtual site is actually getting its interaction(s) set |
|
I ought to fix this upstream first, since it's broken with or without plugins, and since fixing that should make this more straightforward: openforcefield/openff-interchange#982 |
|
But the example simulation starts okay looking okay (in my funky local setup): |
| if len(failures) > 0: | ||
| pytest.fail(f"failures at indices {failures} with energies {reported_energies}") |
There was a problem hiding this comment.
Getting this in CI:
E Failed: failures at indices [0, 1, 2] with energies [12694.954895019531, -35.317028284072876, -23.049033522605896]
I can't wrap my head around why these should be different with the interaction on the virtual site - do these reference values need to be re-computed or is something being set incorrectly here?
There was a problem hiding this comment.
Any idea what's going on here @jthorton?
mattwthompson
force-pushed
the
base-virtual-site-type
branch
from
50ac274
to
c90a696
Compare
Based off of #69 and meant to coordinate with openforcefield/openff-toolkit#1837