Implement Library Charges

[davidlmobley]

Major design decision:

Currently, the charge method is handled on a "by system" basis; that is, as an argument to createSystem. So, all the components of a System are charged via the same mechanism, such as AM1-BCC, and we have no mechanism to flag specific components of a system for special treatment. Often, this is what we want, but water poses challenges (see #24 ) -- typically, for existing water models, we want to assign library charges and only charge our solute (or non-water solvent) atoms via normal charging schemes like AM1-BCC.

How shall we achieve this? @jchodera , @bannanc , @andrizzi, thoughts? Some ideas:

• Create a new "special section" in the SMIRNOFF format (kind of analogous to the BondChargeCorrections section) which allows SMIRKS patterns to specify use of library charges for specific groups
• Allow chargeMethod to be a list rather than a string, where a different charge method can be used for different components of the system

Other possibilities?

[davidlmobley]

See #24 (comment) for possible solutions to this. Probably the best is modifying the nonbonded force section to allow charge overrides, I think.

[jchodera]

I think the addition of allowing an additional charge attribute to override the atomic charge is the simplest approach for handling special water models to enable comparisons with existing models.

The real question here is what are the cases outside of special water models in which we would want to override charge handling? We already have a very flexible approach to molecule-based library charges (by providing OEMol molecules with the desired partial charges) as well as an approach to have createSystem compute and apply AM1-BCC-like charges to all molecules. The gaping hole in this capability is charging of polymeric systems or very large small molecules where AM1-BCC-like methods would be impractical, but are there really use cases beyond special water models where the ability to do something more than what we have now would be useful?

[davidlmobley]

I can't think of that many cases offhand aside from the issues you mention.

I have done work where I had sugars in a protein where I wanted to use GLYCAM parameters for the sugars (with their own charges) and "normal" parameters for the protein and ligand and co-solvents. But in that case, to avoid having to create my own GLYCAM parameterization engine, I presumably would have just created a ParmEd system for the GLYCAM components and merged them with the remainder of the system. So I don't necessarily see this as a use case.

But, there are other co-solvents one might sometimes want to use which are somewhat analogous to waters, I believe. For example, presumably someone has custom TRIS, or ... We probably want a mechanism for people to be able to say, "OK, do your normal thing on all of the rest of the stuff, but treat these components in this special way..." (Or, is that opening too big a can of worms?)

Another case is probably ions, but the issue there is less about charges and more about what other parameters one wants to use.

[jchodera]

But, there are other co-solvents one might sometimes want to use which are somewhat analogous to waters, I believe. For example, presumably someone has custom TRIS, or ... We probably want a mechanism for people to be able to say, "OK, do your normal thing on all of the rest of the stuff, but treat these components in this special way..." (Or, is that opening too big a can of worms?)

Specifying specific forcefields for complex buffer molecules like TRIS are harder to specify in ffxml format because the SMARTS strings will have to be ultra-specific to match individual atoms, bonds, angles, and torsions. It can still work, but it will become a huge pain to generate these manually at some point. I'm not sure how to fix that in a SMIRNOFF-friendly way, however. We could extend the syntax significantly to allow us to specify what amounts to the specific force terms applied to a specific molecule, but I'm not sure if the use case is so compelling that it would ever be worth the effort.

Another case is probably ions, but the issue there is less about charges and more about what other parameters one wants to use.

Supporting ion models is going to be important in both the short term (borrowing ions from elsewhere) and long term (allowing flexibility in parameterization and use of ion models). We should make sure that SMIRNOFF is able to easily let us bring in ion models from AMBER right now, and hopefully also support multisite ion models.

[jchodera]

Following #86, here's how library charges might work:

<LibraryCharges charge_unit="elementary_charge">
   <!-- match an alanine residue -->
   <Residue name="ALA" smirks="[$([NX3H:1](C)(C))][CX4H:2]([CH3X4:3])[CX3:4](=[OX1:5])([N])" charge1=0.2 charge2=-0.2 charge3=0.4 charge4=-0.1 charge5=-0.3>
   ...
</LibraryCharges>

This is just an example (and I neglected to use explicit hydrogens, though we would want to use them), but it illustrates how we might be able to specify a number of SMARTS/SMIRKS matches that could be applied to provide library charges for biopolymers (or other small molecules) prior to subjecting the rest of the uncharged molecules (or molecule fragments) in the system to an AM1-BCC-like scheme.

[davidlmobley]

That looks good to me, @jchodera , and is what I'd basically been thinking in the interim.

[jchodera]

Great. It should be very easy to add support for this.