openforcefield / openff-toolkit

The Open Forcefield Toolkit provides implementations of the SMIRNOFF format, parameterization engine, and other tools. Documentation available at http://open-forcefield-toolkit.readthedocs.io
http://openforcefield.org
MIT License
313 stars 92 forks source link

VdW Sigma Incorrectly Calculated #807

Closed SimonBoothroyd closed 3 years ago

SimonBoothroyd commented 3 years ago

Describe the bug

From version 0.8.1 of the toolkit the sigma parameter is incorrectly calculated as rmin_half / 2 ** (1/6) when it should be computed as 2 * rmin_half / 2 ** (1/6), i.e. it will be off by a factor of two. The reverse is also true when computing rmin_half from sigma.

This bug was introduced in #750.

As, due to the same PR, sigma will now never be None, its value will always be used when creating an OpenMM system. Hence, any force field which defines parameters in terms of rmin_half (all of the OpenFF production force fields) will be incorrectly applied to systems, and the science likely invalid.

I discovered this bug as all of my simulations NaN'd after upgrading to the latest toolkit version so hopefully this should be obvious when it affects users, although users just doing an energy minimisation or a single point evaluation will likely not notice this bug even though it will be present.

To Reproduce

    molecule = Molecule.from_smiles("C")
    force_field = ForceField("openff-1.2.0.offxml")

    system = force_field.create_openmm_system(molecule.to_topology())

    with open(f"system-{openforcefield.__version__}.xml", "w") as file:
        file.write(XmlSerializer.serialize(system))

Output

0.8.0:

<Particles>
    <Particle eps=".4577296" q="-.10868000239133835" sig=".3399669508423535"/>
    <Particle eps=".06568879999999999" q=".027170000597834587" sig=".2649532787749369"/>
    <Particle eps=".06568879999999999" q=".027170000597834587" sig=".2649532787749369"/>
    <Particle eps=".06568879999999999" q=".027170000597834587" sig=".2649532787749369"/>
    <Particle eps=".06568879999999999" q=".027170000597834587" sig=".2649532787749369"/>
</Particles>

0.8.1:

<Particles>
    <Particle eps=".4577296" q="-.10868000239133835" sig=".16998347542117676"/>
    <Particle eps=".06568879999999999" q=".027170000597834587" sig=".13247663938746845"/>
    <Particle eps=".06568879999999999" q=".027170000597834587" sig=".13247663938746845"/>
    <Particle eps=".06568879999999999" q=".027170000597834587" sig=".13247663938746845"/>
    <Particle eps=".06568879999999999" q=".027170000597834587" sig=".13247663938746845"/>
</Particles>

Computing environment (please complete the following information):

Additional context Add any other context about the problem here.

mattwthompson commented 3 years ago

Ideally we should look into adding regression tests so that human eyes alone (even if multiple pairs of them) aren't the only defense against bugs like this. Maybe the current benchmarking effort can be made into a small suite of tests that run periodically?

j-wags commented 3 years ago

Thanks for catching this, @SimonBoothroyd. I agree that this is an error, and we should make a release immediately.