SMARTS for hexa-coordinated compounds

[summeraz]

Are we able to define SMARTS for atoms attached to a hexa-coordinated atom with the SMARTS functionality we currently have? The compound of interest in this case is PF6. If all angles within the molecule are to be considered, then three fluorine atom types will be necessary. All Fx-P-Fx angles should be 180 degrees and all Fx-P-Fy angles should be 90 degrees. However, because the atom types are dependent on orientation rather than bonding, I'm not sure we're able to assign atom types with our current SMARTS functionality.

One possible solution to this involves adding the NOT operator "!" into our SMARTS logic. I'm then envisioning an XML that would look something like this:

<Type name="P" class="P" element="P"  mass="30.9738"  def="P" desc="phosphorous"/>
<Type name="F1" class="F1" element="F"  mass="18.9984"  def="FP([!%F1])([!%F1])([!%F1])([!%F1])F" desc="fluorine type 1"/>
<Type name="F2" class="F2" element="F"  mass="18.9984"  def="FP([!%F1;!%F2])([!%F1;!%F2])(F1)(F1)F" desc="fluorine type 2" overrides="F1"/>
<Type name="F3" class="F3" element="F"  mass="18.9984"  def="FP(F1)(F1)(F2)(F2)F" desc="fluorine type 3" overrides="F1,F2"/>

This should atom type two fluorines each as F1, F2, and F3. The geometry wouldn't be guaranteed, but a simple energy minimization would fix that.

[davidlmobley]

I'm just spectating in general, since I'm not currently involved in this project (but am an interested onlooker). However, I can't resist chiming in on this one. I'd encourage you guys to get in as much of SMARTS as you can - which means you should be able to use "not", when appropriate.

I'm not totally sure I understand the problem here, though -- but just to raise the question, is this a problem that would be more easily solved with SMIRKS, where you have the ability to tag an atom to refer to it? e.g. [#8:1]~[#6:2]~[#9:3] would be a oxygen-carbon-fluorine, where the carbon is tagged as "atom 1" and the fluorine as "atom 3".

[summeraz]

My explanation of the problem was not very clear, but basically compounds featuring a hexa-coordinated atom are an interesting edge case where two different angles are present between atom triplets that would normally have the same atom types (in the case of PF6, there are F-P-F angles of both 180 degrees and 90 degrees). If all fluorines are given the same atom type (as would normally be the case due to all fluorines featuring the same connectivity/local chemistry), Foyer would be unable to choose which F-P-F angle to assign. The workaround is to define three fluorine atom types (F1, F2, F3) such that F-P-F angles of like fluorine types (e.g. F1-P-F1) are 180 degrees, and those of unlike fluorine types (e.g. F1-P-F2) are 90 degrees. However, because there are no differences in the local chemistry for each fluorine, we run into a problem with how to define the SMARTS string for each of the F1, F2, and F3 atom types.

I think adding the "not" operator into Foyer and defining the SMARTS as I showed above will get around this problem. The plan to add this operator (and more of the SMARTS functionality) into Foyer has already been in the works.

Having the ability to assign a tag to individual atoms as you've suggested, following the SMIRKS format would also solve this problem I think. The issue here is where the tag would live within the data structure. This could possibly be added as an attribute to an mBuild Compound, but internally would need to be able to be passed on to both ParmEd Structures and OpenMM Topologies.

[ctk3b]

I'd encourage you guys to get in as much of SMARTS as you can - which means you should be able to use "not", when appropriate.

That's absolutely high priority. It's hard to tell people that we're SMARTS based but you can only use a tiny subset :)

is this a problem that would be more easily solved with SMIRKS, where you have the ability to tag an atom to refer to it?

I don't think so. As @summeraz explained, the weirdness here comes from the fact that the force field angle types are written in a geometry dependent way - a pure topology doesn't contain enough information to directly resolve whether a given F-C-F angle is supposed to be the 90 or 180 degree term.

[ctk3b]

@summeraz I think I understand what you were saying the other day now. I suspect your idea could work. I am curious how trivially an energy minimization would resolve this situation though. We should be getting more SMARTS functionality built in early this week so we can give it a whirl.

[davidlmobley]

Ah, OK, I understand the issue now. You just need to end up with three fluorines of each type (and yes, the minimization should sort it out, er, hopefully. One would need to check.) So it seems like you'd want to just come in at the first pass and make them all (type 1) and then make type 2 be any fluorine which is connected to a P which already has at least three of (type 1). There should be a variety of well-formed strings which will let you do this. But I agree, you'll likely need the "not" operator.

This is a bit odd, though, as to handle this you need to make your SMARTS strings reference atom types, not chemistry, right?? To me this seems like the language is telling you that you're committing an "unnatural act".

I wonder if there is a way to get the geometry you want without doing this. (Though, I suppose your perspective is that you're trying to put the existing force field into this format, not find a more natural way of getting the same geometry even if it means modifying the force field...)

[ctk3b]

Though, I suppose your perspective is that you're trying to put the existing force field into this format, not find a more natural way of getting the same geometry even if it means modifying the force field...

This is the crux of the problem. There are a few weird edge cases in force fields because they were not developed with considerations like this in mind. Something that OpenFF will hopefully solve but doesn't immediately allow us to use existing force fields 😕