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

Non-bonded interactions set_params() is broken #4569

Closed
jngrad opened this issue Sep 12, 2022 · 0 comments · Fixed by #4558
Closed

Non-bonded interactions set_params() is broken #4569

jngrad opened this issue Sep 12, 2022 · 0 comments · Fixed by #4558
Assignees
Labels

Comments

@jngrad
Copy link
Member

jngrad commented Sep 12, 2022

The non-bonded interactions set_params() method has undocumented behavior that can lead to unphysical potentials.

Once a non-bonded interaction is activated with valid parameters using set_params(), the second call to set_params() silently changes behavior: required parameters become optional, and missing arguments are populated using the values currently in the core, completely bypassing the default_params() method.

This unexpected behavior can compromise simulation scripts in two ways:

  1. The user intends to reset an interaction with only the required parameters, expecting the default parameters to populate the missing non-required parameters: the old parameters are silently recycled
  2. The user discovers this undocumented behavior and intends to use it to update one parameter of a LennardJonesInteraction or GenericLennardJonesInteraction potential that was previously set up with shift="auto": the shift parameter is not automatically recalculated and the potential has a jump in energy.

Here is a MWE to reproduce the bug in the second scenario:

import espressomd
import numpy as np
import matplotlib.pyplot as plt

system = espressomd.System(box_l=3 * [10.])
system.time_step = .1

def sample():
    start_pos = [0., 0., 0.]
    step = np.array([0.01 / 600, 0., 0.])
    system.part.clear()
    p0, p1 = system.part.add(pos=[start_pos, start_pos + 600 * step], type=[0, 0])
    xval = []
    yval = []
    for _ in np.arange(180):
        p1.pos = p1.pos + step
        xval.append(np.linalg.norm(p1.pos - p0.pos))
        yval.append(system.analysis.energy()["non_bonded"])
    return np.array(xval), np.array(yval)

system.non_bonded_inter[0, 0].lennard_jones.set_params(epsilon=1.92, sigma=0.05, cutoff=0.11, shift="auto")
# trigger bug: forget to use shift="auto", so the shift isn't recalculated for the new cutoff
system.non_bonded_inter[0, 0].lennard_jones.set_params(cutoff=0.012)
xval1, yval1 = sample()
# fix bug: manually add shift="auto"
system.non_bonded_inter[0, 0].lennard_jones.set_params(cutoff=0.012, shift="auto")
xval2, yval2 = sample()

plt.plot(xval2, yval2, label="set_params(cutoff=0.012, shift=\"auto\")")
plt.plot(xval1, yval1, label="set_params(cutoff=0.012)")
plt.xlabel("Inter-particle distance")
plt.ylabel("Energy")
plt.legend()
plt.show()

Plot of the MWE showing the energy jump for a LJ potential

@jngrad jngrad added the Bug label Sep 12, 2022
@jngrad jngrad added this to the Espresso 4.3.0 milestone Sep 12, 2022
@jngrad jngrad self-assigned this Sep 12, 2022
@kodiakhq kodiakhq bot closed this as completed in #4558 Sep 20, 2022
kodiakhq bot added a commit that referenced this issue Sep 20, 2022
Fixes #4555 and fixes #4569

Description of changes:
- rewrite non-bonded interactions as script interface objects
   - convert Cython file `interaction.pyx` to Python file `interaction.py`
   - 1600 lines of Cython code were removed
- bugfix: the `set_params()` method now raises an error if required arguments are not provided
   - scripts may silently break when a potential is set and later updated in a way that is affected by the bug described in #4569
- new feature: possibility to deactivate/reset non-bonded interactions
- reduce code duplication in the script interface of bonded interactions
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

1 participant