Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

VdW Sigma Incorrectly Calculated #807

Closed
SimonBoothroyd opened this issue Jan 8, 2021 · 2 comments · Fixed by #808
Closed

VdW Sigma Incorrectly Calculated #807

SimonBoothroyd opened this issue Jan 8, 2021 · 2 comments · Fixed by #808

Comments

@SimonBoothroyd
Copy link
Contributor

SimonBoothroyd commented Jan 8, 2021

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):

  • Operating system: OSX
  • Output of running conda list

Additional context
Add any other context about the problem here.

@mattwthompson
Copy link
Member

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
Copy link
Member

j-wags commented Jan 8, 2021

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants