use numerical stability parameter for calculatiung d_ij
Currently the numerical stability parameter epsilon is only used in PaiNNMixing for calculating the norm and not for calculating the norm dir_ij.
For training force field models this is not necessary, because force field datasets do not contain any structures where two atoms come very close. However, in generative models (e.g. diffusion) it can happen that two atoms are located at the same position.
I suggest to add the same numerical stability parameter epsilon as for calculating the norm in PaiNNMixing:
105 mu_Vn = torch.sqrt(torch.sum(mu_V**2, dim=-2, keepdim=True) + self.epsilon)
I also added the argument return_vector_representation in PaiNN to make it consistent with SO3Net.
two intermediate commits should probably be squashed