Instability of results as exponent approaches 3.

[kvhuguenin]

As discvoered by @E-Rum , the potentials obtained from the calculator seem to diverge as the exponent p-->3:

from torchpme.calculators import PMEPotential
import torch
for exp in [1.,2.,2.5,2.9,2.999,2.9999,2.999999,3.]:
    calculator = PMEPotential(
        exponent=exp,
        atomic_smearing=1.,
        mesh_spacing=0.05,
        interpolation_order=3,
        subtract_interior=True,
        full_neighbor_list=False,
    )
    potentials = calculator(
        positions=torch.rand(10, 3) * 10,
        charges=torch.rand(10,1) * 0.01,
        cell=torch.eye(3) * 10,
        neighbor_indices=torch.randint(0, 10, (15,2)),
        neighbor_distances=torch.rand(15),
    )
    print(min(abs(potentials)))

leading to Outputs:

tensor([0.0018])
tensor([0.0008])
tensor([0.0026])
tensor([0.0031])
tensor([0.3840])
tensor([2.7494])
tensor([210.4061])
tensor([nan])

Analytically, both the real and reciprocal space sums converge absolutely for any real number p>0, suggesting that this is due to a numerical instability. The obvious candidate in the current implementation is the incomplete gamma function gammaincc(a,x) = gammaincc((3-p)/2, x) whose first argument a needs to be positive for the torch (and also scipy) implementations. As a approaches zero, there could be numerical instabilities due to this.

[PicoCentauri]

Interesting. We can look into the torch code together next week if we can fix something there.

[ceriottm]

The ideal solution here would be to implement something different that also handles arbitrarily large (integer) exponents.