[gmx-developers] Use of invsqrt with tables

[Berk Hess]

On 03/23/2016 04:25 PM, David van der Spoel wrote:
> On 23/03/16 14:42, Berk Hess wrote:
>> 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
>>>> much of this with masked (SIMD) operations instead of conditionals,
>>>> 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.
> But then one can use simply r = sqrt(r2)
You're right.
But then we need separate tables for the energy and force, since we 
can't use the /r trick for the combined cubic or quadratic spline table.
I still prefer my solution of:
invr = sqrt(r^2 + GMX_REAL_MIN)
This does not modify the results for float r down to 1e-15, which should 
be fine for all purposes. For smaller r it results in a invr of 1e-19, 
which is much larger than GMX_FLOAT_MIN, so we can safely do arithmetic 
with it.

Cheers,

Berk