Division by zero in distance calculation of sphere constraint

[RiccardoFrenner]

When placing a particle at the center of a sphere constraint, the distance between the particle and the constraint is calculated incorrectly, resulting in a force on the particle which contains NANs.

Here is a MWE:

import espressomd
import espressomd.shapes

PARTICLE_TYPE = 0
CONSTRAINT_TYPE = PARTICLE_TYPE

system = espressomd.System(box_l=3*[10], periodicity=3*[False])
system.time_step = 1
system.cell_system.skin = 0.1

sphere_shape = espressomd.shapes.Sphere(
    center=0.5*system.box_l,
    radius=0.25*system.box_l[0],
    direction=-1
)
sphere_constraint = espressomd.constraints.ShapeBasedConstraint(
    shape=sphere_shape, penetrable=True, particle_type=CONSTRAINT_TYPE)
system.constraints.add(sphere_constraint)

system.non_bonded_inter[PARTICLE_TYPE,
                        CONSTRAINT_TYPE].wca.set_params(epsilon=1., sigma=1.)

assert(system.constraints[0].penetrable == True)

system.part.add(pos=system.box_l/2, type=PARTICLE_TYPE)

print(system.constraints[0].min_dist())  # 2.5
system.integrator.run(0)
print(system.part[:].f)  # [[nan, nan, nan]]

system.integrator.run(1)  # ERROR: Constraint violated by particle 0 dist -nan

The scalar distance is actually calculated correctly, but the distance vector vec in Sphere.cpp contains NANs due to c_dist being zero.

[RudolfWeeber]

What should we do in this case? Return a vector with the correct length and a random orientation? (since the center is equally far from all points on the spherical surface)

[RiccardoFrenner]

I think it should be similar to the cylinder constraint, and therefore return a vector of the correct length and any direction.