[gmx-developers] Use of invsqrt with tables

[David van der Spoel]

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)

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


-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se