[gmx-developers] Use of invsqrt with tables

[Berk Hess]

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.

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