Possible bug in half_cylinder_rotations

[Tristan Croll]

I call it a "possible bug" since it's only arisen for me since I copied the logic to a separate function (used for a position restraint implementation that avoids using fake atoms and pseudobonds, copied to a header file so I can call it directly from C++) and used it in double precision (code copied below). In ChimeraX's _geometry/cylinderrot.cpp the lines:

      if (d == 0)
	{ vx = vy = 0 ; vz = 1; }
      else
	{ vx /= d; vy /= d; vz /= d; }

seem to work just fine - but in double precision I found I had to make the modifications involving ALMOST_ZERO (currently defined as 1e-6) to avoid Numpy throwing up its hands over singular matrices when the bond length is zero (i.e. when the target positions are set to the atom coordinates).

template <typename T>
void bond_cylinder_transform_gl(T xyz0[3], T xyz1[3], T r, T *rot44)
{
    T bvec[3];
    for (size_t i=0; i<3; i++) {
        bvec[i] = xyz1[i]-xyz0[i];
    }
    T d = l2_norm_3d(bvec);
    if (d < ALMOST_ZERO) {
        bvec[0]=0; bvec[1]=0; bvec[2]=1;
        d = ALMOST_ZERO;
    } else {
        for (auto &b: bvec)
            b/=d;
    }
    T &c = bvec[2], c1;
    if (c <= -1) c1 = 0;
    else c1 = 1.0/(1.0+c);
    T wx = -bvec[1], wy = bvec[0];
    T cx = c1*wx, cy = c1*wy;
    T h = d*2;
    *rot44++ = r*(cx*wx + c);
    *rot44++ = r*cy*wx;
    *rot44++ = -r*wy;
    *rot44++ = 0;

    *rot44++ = r*cx*wy;
    *rot44++ = r*(cy*wy + c);
    *rot44++ = r*wx;
    *rot44++ = 0;

    *rot44++ = h*wy;
    *rot44++ = -h*wx;
    *rot44++ = h*c;
    *rot44++ = 0;

    *rot44++ = (xyz0[0]+xyz1[0])/2;
    *rot44++ = (xyz0[1]+xyz1[1])/2;
    *rot44++ = (xyz0[2]+xyz1[2])/2;
    *rot44++ = 1;
}

[goddard@…]

Avoiding divide by zero is tricky to do right.  I agree my test for equality to zero has problems although I have never seen it cause trouble in that ChimeraX code.  Unfortunately your fix has its own problems.  For instance if I decide to use units of meters for my molecules then all bond lengths are now smaller than your divide by zero cutoff of 1e-6 and all ChimeraX bonds would be drawn parallel to the z axis just because I switched units.  The vectors maintain their full precision scaled down by 1e-10.  So your test has introduced an assumption about the magnitudes of the vectors.  You may think that is far fetched but I once made an animation zooming from full human brain imaging down to molecular scale and units of meters would be very reasonable.  I think the only right solution is to use a serious vector math library normalization routine that has the code to handle loss of precision scenarios.  But I would just use your solution if it works for you.  But since a double has about 15 decimal digits precision I might lower your cutoff to 1e-10.

[Tom Goddard]

The resulting error is pretty easy to debug and I have no easy way to reduce the error potential so will ignore this until it becomes a problem.