Rounding error in acos function

[Carlos Bueno]

Hi,
I was working with the CustomHbond force using an acos function to calculate the angle between two vectors. In order to compute the angle between two vectors I used the function:

t3 = acos(cost3lim);
cost3 = sin(t1)*sin(t2)*cos(phi)-cos(t1)*cos(t2);

This function shouldn't be more than 1 or less than -1 but I noticed that I was getting the following error frequently during my simulation:

Particle coordinate is nan

Which may have happened because of a rounding error. I tried to add a limit to this function using -1 and 1 but it didn't fix the problem.

t3 = acos(cost3lim);
cost3lim = min(max(cost3,-1),1);
cost3 = sin(t1)*sin(t2)*cos(phi)-cos(t1)*cos(t2);

Nevertheless restricting the limit further to 0.9999 fixes it:

t3 = acos(cost3lim);
cost3lim = min(max(cost3,-0.9999),0.9999);
cost3 = sin(t1)*sin(t2)*cos(phi)-cos(t1)*cos(t2);

This problem seems to be depend on the force used (CustomNonBonded vs CustomHbond) and the platform (CPU vs OpenCL). I include a small example that reproduces the problem using 3 types of forces:

acos2.zip

(3.6 KiB) Downloaded 1 time

"CustomNonbondedForce2": Only contains the acos function
nonbondedForce=simtk.openmm.CustomNonbondedForce(f"acos({lim})")


"CustomNonbondedForce": Contains the acos function and the min/max
nonbondedForce=simtk.openmm.CustomNonbondedForce(f"""acos(cost3lim);
cost3lim = min(max(cost3,-{lim}),{lim});
cost3=4*(r-1)^2-2;""")
"CustomHbondForce": Contains the acos function and the min/max and uses CustomHbond instead of CustomNonbondedForce
nonbondedForce=simtk.openmm.CustomHbondForce(f"""acos(cost3lim);
cost3lim = min(max(cost3,-{lim}),{lim});
cost3=4*(r-1)^2-2;
r=distance(a1,d1)""")

I include the output with annotations here:

Running: Simulation in platform CPU using force CustomNonbondedForce2 with limit 0.9999
(Completed)

Running: Simulation in platform CPU using force CustomNonbondedForce with limit 0.9999
(Completed)

Running: Simulation in platform CPU using force CustomHbondForce with limit 0.9999
(Completed)

Running: Simulation in platform OpenCL using force CustomNonbondedForce2 with limit 0.9999
(Completed)

Running: Simulation in platform OpenCL using force CustomNonbondedForce with limit 0.9999
(Completed)

Running: Simulation in platform OpenCL using force CustomHbondForce with limit 0.9999
(Completed)

Running: Simulation in platform Reference using force CustomNonbondedForce2 with limit 0.9999
(Completed)

Running: Simulation in platform Reference using force CustomNonbondedForce with limit 0.9999
(Completed)

Running: Simulation in platform Reference using force CustomHbondForce with limit 0.9999
(Completed)

Running: Simulation in platform CPU using force CustomNonbondedForce2 with limit 0.99999999
(Completed)

Running: Simulation in platform CPU using force CustomNonbondedForce with limit 0.99999999
(Completed)

Running: Simulation in platform CPU using force CustomHbondForce with limit 0.99999999
(Completed)

Running: Simulation in platform OpenCL using force CustomNonbondedForce2 with limit 0.99999999
(Completed)

Running: Simulation in platform OpenCL using force CustomNonbondedForce with limit 0.99999999
(Completed)

Running: Simulation in platform OpenCL using force CustomHbondForce with limit 0.99999999
=====Failed=====
Particle position is NaN
[*]CustomHBond force fails with a different limit

Running: Simulation in platform Reference using force CustomNonbondedForce2 with limit 0.99999999
(Completed)

Running: Simulation in platform Reference using force CustomNonbondedForce with limit 0.99999999
(Completed)

Running: Simulation in platform Reference using force CustomHbondForce with limit 0.99999999
(Completed)

Running: Simulation in platform CPU using force CustomNonbondedForce2 with limit 1
(Completed)

Running: Simulation in platform CPU using force CustomNonbondedForce with limit 1
=====Failed=====
Particle coordinate is nan
[*]The problem seems to occur after including the max/min functions

Running: Simulation in platform CPU using force CustomHbondForce with limit 1
=====Failed=====
Particle coordinate is nan

Running: Simulation in platform OpenCL using force CustomNonbondedForce2 with limit 1
(Completed)

Running: Simulation in platform OpenCL using force CustomNonbondedForce with limit 1
(Completed)
[*]The problem on CustomNonbondedForce doesn't happen using OpenCL

Running: Simulation in platform OpenCL using force CustomHbondForce with limit 1
=====Failed=====
Particle coordinate is nan

Running: Simulation in platform Reference using force CustomNonbondedForce2 with limit 1
(Completed)

Running: Simulation in platform Reference using force CustomNonbondedForce with limit 1
[*] Simulation hangs when using Reference

[Peter Eastman]

Welcome to the exciting world of floating point numerics, where numbers are never exactly what they appear to be! To understand what's going on, look at a graph of acos(x). Notice that it becomes vertical at 1 and -1. More precisely, its derivative is -1/sqrt(1-x^2). That becomes infinite at x=1 and x=-1. So if you don't want bad things to happen while computing the forces, don't let it reach those points! Even if it doesn't quite reach those points but just gets very close, you'll get a catastrophic loss of accuracy.

Notice, though, that this doesn't mean your potential function is singular. The angle itself is well defined and finite everywhere except when the length of one of the two edges is zero. It's just that you're calculating it in a way that introduces singularities. Could you just use the built in angle() function? Here's the actual code it uses to compute the angle, which carefully works around the singularities.

Code: Select all

real computeAngle(real4 vec1, real4 vec2) {
    real dotProduct = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
    real cosine = dotProduct*RSQRT(vec1.w*vec2.w);
    real angle;
    if (cosine > 0.99f || cosine < -0.99f) {
        // We're close to the singularity in acos(), so take the cross product and use asin() instead.

        real4 crossProduct = cross(vec1, vec2);
        real scale = vec1.w*vec2.w;
        angle = asin(SQRT(dot(crossProduct, crossProduct)/scale));
        if (cosine < 0.0f)
            angle = M_PI-angle;
    }
    else
       angle = acos(cosine);
    return angle;
}

[Carlos Bueno]

Thank you.
Is there a way to incorporate that computeAngle function on a CustomHbond force expresion using the python API?

[Peter Eastman]

See the documentation at http://docs.openmm.org/latest/api-pytho ... HbondForce. Just use the function angle(...) to compute the angle between any three specified particles.

[Carlos Bueno]

Hi Peter,
Sorry for the late response. I tried to use the compute Angle function, but I don't seem to be able to reach it from the python API or the XML setup. I want to compute the angle between 2 vectors defined by 4 particles, so the angle(...) function from customHBondForce doesn't meet my needs. The only way to compute the angle between vectors seems to be to define the angle as:

t3 = acos(cost3lim);
cost3lim = min(max(cost3,-0.9999),0.9999);
cost3 = sin(t1)*sin(t2)*cos(phi)-cos(t1)*cos(t2);
t1 = angle(d2,d1,a1);
t2 = angle(d1,a1,a2);
phi = dihedral(d2,d1,a1,a2);

,where the two vectors are d2->d1 and a2->a1. But that introduces the singularity in the derivative as you mentioned (unless I use the cutoff of 0.9999). Should I try to expose the computeAngle function to the python API? I would like to have this function available:

t3 = vector_angle(d2,d1,a1,a2);

[Peter Eastman]

Could you describe exactly what the potential function you want is? Don't worry about the details of how to compute it. Just describe the potential in the simplest way possible. I think I'm getting lost in the trees, and that's keeping me from seeing the forest.

[Carlos Bueno]

Yes. I'm implementing the Cross-Stacking term of the 3SPN.2C potential. The Cross-Stacking potential is a function of 5 particles: B0, S1, B1, S2 and B2, where B0, S1, and B1 are part of one chain and S2 and B2 are part of another chain (B stands for base and S for Sugar). The potential is a function of the angle between the vector S1->B1 and the vector S2->B2(t_3), the angle between the particles S2, B2 and B0(t_CS), and the distance between B1 and B2(r_ij)

Ucstk=f(t_3)*f(t_CS)*f(r_ij)

It is described on this paper:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3808442/