[gmx-developers] Use of invsqrt with tables

Berk Hess hess at kth.se
Wed Mar 23 14:42:55 CET 2016

```On 2016-03-23 14:19, David van der Spoel wrote:
> On 23/03/16 14:07, Erik Lindahl wrote:
>> Hi,
>>
>>
>>> On 23 Mar 2016, at 14:00, Berk Hess <hess at kth.se> wrote:
>>>>>
>>>>> But avoiding the inversion is not possible, since you need to
>>>>> multiply
>>>>> the scaling force by the normalized distance vector.
>>>> The force should be zero at r=0 for symmetry reasons and hence the
>>>> code should not crash.
>>> I know that. But that doesn't help in avoiding the division needed
>>> to get the normalized distance vector. You can avoid the division by
>>> exactly zero using a conditional. But you will still have issues for
>>> distances very close to zero. A solution I use for forces between
>>> excluded atoms in the non-bonded kernel, which can potentially
>>> overlap, is to add a very small value epsilon to r^2, such that
>>> invsqrt produces a large number that not close to max real.
>>> Multiplying this with r gives a value close to zero that's slightly
>>> off from 1, but the error will be less than real precision for
>>> potentials that are symmetric around 0.
>>
>> I have planned infrastructure for the new kernels that can handle
>> but my priority for the 2016 release is to get rid of group kernels,
>> so this might not make it…
>>
>
> The issue is the interaction between two Gaussian-smeared charges.
> Both the energy and the force go to zero smoothly. For distances close
> to zero we have to multiply the scalar force by the normalized
> distance vector. If we do the normalization first that should not give
> numerical problems, right?
You would still need to do the normalization conditional, which is not
good for performance. And you need to be very careful.
I prefer my solution: invr = sqrt(r^2 + 1e-38)
>
> A further alternative (tables only) would be to store F/r in the
> tables if that is well behaved (which it is for the Gaussian charges).
This would probably work, but you need to fill the table carefully.
However, you still need r for getting the table indices, so this doesn't
help much.

Berk
>
>> Cheers,
>>
>> Erik
>>
>
>

```