Likely upcoming change to partial charge assignment using the OpenFF Toolkit with the RDKit+AmberTools backend

[j-wags]

Summary

We may change how charge assignment is done using the Open Force Field Toolkit with the RDKit+AmberTools backend. This would likely involve using a single “Electrostatically Least interacting Functional group” (ELF) conformer instead of a single random one, and using certain keywords when calling antechamber. These changes would correspond to a new version (0.4) of the ToolkitAM1BCC tag in the SMIRNOFF specification.

What this means for users

• Moving forward, OpenFF Toolkit users that rely on the RDKi+AmberTools backend are likely to see an improvement in parameter quality and reproducibility.
• Users of the OpenFF Toolkit with the OpenEye backend would see no change.
• New versions of the OpenFF Toolkit would only support the new behavior, and if a force field with the 0.3 ToolkitAM1BCC tag were loaded, a warning would be printed saying that the ToolkitAM1BCC tag is being auto-updated to version 0.4
• Users who want to reproduce the "old" (0.3) behavior of the ToolkitAM1BCC tag with the RDKit+AmberTools backend would need to install a legacy version of the OpenFF Toolkit with those backends.

What this means for developers

• There will likely be a SMIRNOFF specification update to the ToolkitAM1BCC tag to version 0.4, and new force fields will be published with this tag version. The new specification is very unlikely to change the desired behavior from the OpenEye backend. But it will add a suggestion that other backends use ELF conformer selection for AM1BCC charge assignment, and it may specify specific antechamber keywords or behaviors.

Background

All of OpenFF's main-line force fields include the ToolkitAM1BCC tag, which basically means "use AM1BCC to assign charges to small molecules in the topology". By default, this behavior explicitly ignores the geometr(ies) and partial charges, if any, associated with the input molecules. The reasoning behind this behavior is that SMIRNOFF parameter assignment only looks at the chemistry of the input molecule, and aims to assign identical parameters to chemically identical molecules, regardless of factors like input conformer or local environment. The current behavior in an RDKit+AmberTools build of the OpenFF Toolkit is:

• RDKit is called to generate one conformer of the molecule, using a constant set of keywords, and (currently) if the initial attempt fails a keyword argument is changed to try again before raising an error.
• Antechamber is called using a specific set of keywords to take the RDKit-generated conformer and perform an AM1 optimization and then return Mulliken charges from the optimized structure.

We've been shown preliminary findings from a user suggesting that this workflow can produce different charges for the same molecule under different circumstances, which in some cases leads to significant differences in simulation results. The current findings show that this is due to:

• RDKit generating different conformers
• antechamber optimizing the same geometry to different minimum energy structures

in different circumstances, which are thought to include:

• differences in availability of system resources at the time of calling these programs
• different system architectures
• different OSes/system libraries

The user hasn't finished and formally shared the study, but preliminarily suggested that two changes could improve the current behavior:

• Having RDKit generate a large number of conformers, and then using the "Electrostatically Least interacting Functional groups" (ELF) method to select a single conformer from the ensemble.
• Setting additional keywords to antechamber to reduce runtime variability, for example by explicitly specifying one matrix diagonalization routine to use instead of having the program automatically select one at runtime.

Impact on parameter reproducibility and quality

The OpenFF lead team discussed these potential changes in terms of their impact on two metrics:

Reproducibility

• with respect to the RDKit conformer generation step: The report that we received showed that the conformers generated during partial charge assignment differed more than we had expected. So the reproducibility of the existing behavior has room for improvement. In considering the new behavior, in the case that two different computers behave identically, generating a large number of conformers and selecting one using the ELF method will lead to the same conformer being used for charging. But in cases where two computers generate DIFFERENT ensembles, using the ELF method is likely to result in a similar conformer being selected from both ensembles, and similar partial charges emerging from both.
• with respect to the antechamber optimization and charge assignment step: Forcing the use of the same settings (such as diagonalization algorithm) by AmberTools is almost certain to improve the reproducibility of the resulting partial charges.

Quality

• with respect to the RDKit conformer generation step: Using the ELF method to select a conformer is expected to lead to systematic changes in partial charges, since the electrostatically least interacting conformer will experience less self-polarization than other possible conformers. However, our main-line force field training is performed using the OpenEye backend, which DOES use the ELF method during conformer generation in partial charge assignment. So by making the partial charge assignment performed by the RDKit+AmberTools backend more similar to that performed by the OpenEye backend, which is used in our force field training, this change is likely to improve the quality of the results that users on the RDKit+AmberTools backend get.
• with respect to the antechamber optimization and charge assignment step: Forcing the use of specific antechamber settings (such as diagonalization routine) may lead to different results in some cases, but these differences are unlikely to be systematically better or worse.

So in summary, in every area that we considered, the proposed changes are either expected to not have a systematic impact, or to systematically improve the quality and consistency of results. Therefore we are likely to approve+implement this user's specific recommendations once they present them to us in a final form.

SMIRNOFF specification update

This change would probably be accompanied by a change to the ToolkitAM1BCC description in the SMIRNOFF specification. This tag already specifies use of the ELF10 method when the OpenEye backend is used, but this would be extended to also specify the use of single-conformer ELF when other backends are used. Since substantial updates to the SMIRNOFF specification are recorded as section or whole-file version changes, the ToolkitAM1BCC section would be updated from version 0.3 to 0.4.

Process

This change will likely proceed in three stages:

• The user presents their final report to us (and we make it public)
• The change is proposed as a OFF EP and voted on by the SMIRNOFF committee
• The change is implemented in the software and in future force fields

We're happy to hear your thoughts, especially if we're overlooking factors that should cause us not to proceed with these changes. That said, we've also learned that OFF EPs can quickly balloon in scope, so please be understanding that we have limited resources and would prefer to consolidate gradual gains rather than fail to execute on grand plans.

[hjuinj]

Thanks for the investigation, it is interesting to know that the OE backend is using ELF. I always thought they do conformer ensemble averaged charges. I wonder how well the ELF approach handles with floppier chemical matters (e.g. going towards the beyond rule of 5 region), where I would have hoped an (diverse) conformer ensemble averaged partial charges assignment can better model such molecules.

[j-wags]

it is interesting to know that the OE backend is using ELF. I always thought they do conformer ensemble averaged charges.

They kinda do - OpenEye uses ELF10, which means "average the partial charges from up to 10 conformers, selected using the ELF method". There's some secret sauce in there, like the RMS cutoff between the conformers and how they choose whether to use 10 conformers or fewer.

The big reason we're not proposing to implement multi-conformer AM1BCC ELF in the RDKit/AmberTools backend in this proposal is just because it would take a really long time for antechamber/sqm to run for multiple conformers.