id	summary	reporter	owner	description	type	status	priority	milestone	component	version	resolution	keywords	cc	blockedby	blocking	notify_on_close	platform	project
1033	Possible bug in half_cylinder_rotations	Tristan Croll	Tom Goddard	"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;
}
}}}"	defect	closed	minor		Core		wontfix						all	ChimeraX