Closed pbuslaev closed 1 year ago
Hey @pbuslaev, nice detective work again. The discrepancy here is much less exciting than you'd hope - the toolkit provides slightly different charges when it uses antechamber
vs. oechem
. More detail here:
This tag calculates partial charges using the default settings of the highest-priority cheminformatics toolkit that can perform AM1-BCC charge assignment. Currently, if the OpenEye toolkit is licensed and available, this will use QuacPac configured to generate charges using AM1-BCC ELF10 for each unique molecule in the topology. Otherwise RDKit will be used for initial conformer generation and the AmberTools antechamber/sqm software will be used for charge calculation.
The ELF10 variant provides slightly different charges for some molecules, but both are considered correct implementations of a SMIRNOFF force field. (This grates me a bit, personally, but electrostatics are notoriously difficult to get right and AM1-BCC vs. AM1-BCC ELF10 is just the tip of the iceberg of different partial charge methods that can be used to similar result.) Under the hood, these tools are called by ToolkitWrapper
s (OpenEyeToolkitWrapper
, AmberToolsToolkitWrapper
, RDKitToolkitWrapper
, ...).
The -15 number results from having access to OpenEye Toolkits (free for academics, and I usually use one in development) and the -6.8 number comes from using AmberTools/RDKit. The notebook you were looking at was probably rendered by a developer who had a license; we have another version floating around that reports the different number, which you got locally. To be a little more thorough, I checked myself what happens when OpenEyeToolkitWrapper
is available or not, and the results look familiar:
from openff.toolkit.utils.toolkit_registry import (
_toolkit_registry_manager,
RDKitToolkitWrapper,
OpenEyeToolkitWrapper,
AmberToolsToolkitWrapper,
ToolkitRegistry,
)
with _toolkit_registry_manager(
ToolkitRegistry(toolkit_precedence=[RDKitToolkitWrapper, AmberToolsToolkitWrapper])):
print(get_summary_data(forcefield.create_interchange(topology)))
with _toolkit_registry_manager(
ToolkitRegistry(toolkit_precedence=[OpenEyeToolkitWrapper])):
print(get_summary_data(forcefield.create_interchange(topology)))
Bond Angle Torsion Electrostatics vdW RBTorsion
OpenMM 0.233879 7.141818 25.580587 -6.840315 9.658657 NaN
Amber 0.233886 7.141670 25.580558 -6.840840 9.727382 NaN
GROMACS 0.233869 7.141836 25.580585 -6.858515 9.667051 0.0
LAMMPS 1.384295 7.141818 25.580587 -6.944298 9.630999 NaN
Bond Angle Torsion Electrostatics vdW RBTorsion
OpenMM 0.233879 7.141818 25.580587 -15.215644 9.658657 NaN
Amber 0.233886 7.141670 25.580558 -15.216371 9.727382 NaN
GROMACS 0.233869 7.141835 25.580582 -15.235316 9.667050 0.0
LAMMPS 1.384295 7.141818 25.580587 -15.301557 9.630999 NaN
Hi @mattwthompson
Thanks for the explanations. It makes perfect sense. Do you think it is worth adding a similar note (or providing additionally values for Amber AM1BCC charges in README for example) to the example notebook, so users are not confused by different energies?
Good idea @pbuslaev! https://github.com/openforcefield/openff-toolkit/pull/1671 has re-rendered the examples with only free (for everybody) and open source tools, and now includes the -6.8 number you expected it to.
I tried to reproduce the example notebook for using SMIRNOFF with GROMACS. I used
openff-toolkit
version0.14.0
and I used the structure file provided in the tutorial directory. The energies that I got areThe reported in the notebook values are:
The difference is only in the Electrostatics term. So I tried to check the assigned charges, by running
antechamber
on structures of cyclohexane and ethanol from1_cyclohexane_1_ethanol.pdb
separately. I compared those charges with charges assigned byopenff-toolkit
, by applyingI provided the conformer, to be sure that charges are calculated for the exact same conformations by
antechamber
andopenff-toolkit
. For both molecules (cyclohexane and ethanol) the charges computed with two different methods were identical, which tells me, that the difference in electrostatic energies is not due to the incorrect charge assignment.Next, I transformed interchange object created from
1_cyclohexane_1_ethanol.pdb
file into GROMACS format:and used default
.mdp
file, as implied by interchangeget_gromacs_energies
. Using created.mdp
,.gro
, and.top
files, I ranand checked the energies both manually and using routine from
openff.interchange
:_parse_gmx_energy
. The energies that I obtained manually, were the same that I obtained withget_gromacs_energies(interchange)
and thus, differed from those reported in example notebook. Could it be that the values reported are outdated?My conde environment is: