Test coverage for tabulated angle potentials missing

[jngrad]

There is some refactoring going on in #2497, with new integration tests in 5c50bd6 covering the tabulated angle potential, among others.

The tabulated dihedral force appears to be zero on all particles.

Bonded angle calculations seem based on the external angle convention. It can be a source of confusion with tabulated angles if the user forgets to flip the energy and force arrays. In addition, the tabulated force drops rapidly when the angle is almost equal to 0 or pi within 5 decimal places. Here is a minimal example:

import espressomd
import numpy as np
system = espressomd.System(box_l=[10.0, 10.0, 10.0])
system.cell_system.skin = 0.4
system.time_step = .1
system.part.add(pos=[1,0,0])
system.part.add(pos=[0,0,0])
system.part.add(pos=[0,1,0])

# create a tabulated angle interaction, not flipping the arrays on purpose
tab_energy = [0.5 * (phi - np.pi / 3)**2 for phi in np.linspace(0, np.pi, 180)]
tab_force = 10 * np.abs(tab_energy)
angle_tabulated = espressomd.interactions.Tabulated(
    type='angle', energy=tab_energy, force=tab_force, min=0., max=np.pi)
system.bonded_inter.add(angle_tabulated)
system.part[1].add_bond((angle_tabulated, 0, 2))

epsilon = 0
#epsilon = 1e-4 # uncomment this line to erase the drop in force
a = []
for i, theta in enumerate(np.linspace(0 + epsilon, np.pi - epsilon, 180)):
    system.part[2].pos = [np.cos(theta), np.sin(theta), 0]
    system.integrator.run(recalc_forces=True, steps=0)
    E_sim = system.analysis.energy()["bonded"]
    F_sim = np.linalg.norm(system.part[0].f)
    j = -i - 1  # here use j = i  if the arrays were flipped
    E_ref = tab_energy[j]
    F_ref = tab_force[j]
    a.append((theta, E_sim, E_ref, F_sim, F_ref))

a = np.array(a)

from matplotlib import pyplot as plt
# plot the force
plt.plot(a[:,0], a[:,3], 'g')  # simulated force
plt.plot(a[:,0], a[:,4], 'b')  # expected force
plt.show()
# plot the energy
plt.plot(a[:,0], a[:,1], 'g')  # simulated energy
plt.plot(a[:,0], a[:,2], 'b')  # expected energy
plt.show()

Also, in /testsuite/python/interactions_bond_angle.py the magnitude of the bond angle force for the left and right particles decays rapidly in the range [0, phi0] regardless of the interaction (harmonic, tabulated, cosine, cossquare). Shouldn't the force be symmetric around phi0?

[RudolfWeeber]

Done