bool MoleculeAlleneStereo::_isAlleneCenter()

in core/indigo-core/molecule/src/molecule_allene_stereo.cpp [131:265]

• 100 lines of code
• 55 McCabe index (conditional complexity)

bool MoleculeAlleneStereo::_isAlleneCenter(BaseMolecule& mol, int idx, _Atom& atom, int* sensible_bonds_out)
{
    bool pure_h[4];

    if (!possibleCenter(mol, idx, atom.left, atom.right, atom.subst, pure_h))
        return false;

    int dirs[4] = {0, 0, 0, 0};
    Vec3f subst_vecs[4];
    int j, k;

    bool zero_bond_length = false;

    for (k = 0; k < 2; k++)
        if (atom.subst[k] >= 0)
        {
            dirs[k] = mol.getBondDirection2(atom.left, atom.subst[k]);
            subst_vecs[k].diff(mol.getAtomXyz(atom.subst[k]), mol.getAtomXyz(atom.left));
            if (!subst_vecs[k].normalize())
                zero_bond_length = true;
        }

    for (k = 2; k < 4; k++)
        if (atom.subst[k] >= 0)
        {
            dirs[k] = mol.getBondDirection2(atom.right, atom.subst[k]);
            subst_vecs[k].diff(mol.getAtomXyz(atom.subst[k]), mol.getAtomXyz(atom.right));
            if (!subst_vecs[k].normalize())
                zero_bond_length = true;
        }

    if (dirs[0] == 0 && dirs[1] == 0 && dirs[2] == 0 && dirs[3] == 0)
        return false; // no oriented bonds => no stereochemistry

    // check that they do not have the same orientation
    if (dirs[0] != 0 && dirs[0] == dirs[1] && dirs[0] != BOND_EITHER)
        return false;
    if (dirs[2] != 0 && dirs[2] == dirs[3] && dirs[2] != BOND_EITHER)
        return false;

    if (zero_bond_length)
        throw Error("zero bond length");

    Vec3f pos_center = mol.getAtomXyz(idx);
    Vec3f vec_left = mol.getAtomXyz(atom.left);
    Vec3f vec_right = mol.getAtomXyz(atom.right);

    vec_left.sub(pos_center);
    vec_right.sub(pos_center);

    if (!vec_left.normalize() || !vec_right.normalize())
        throw Error("zero bond length");

    // they should go in one line
    // 0.04 is equivalent to 16 degress because it is hard to draw a straight line accurately
    if (fabs(Vec3f::dot(vec_left, vec_right) + 1) > 0.04)
        return false;

    // check that if there are two left substituents, they do not lie on the same side
    if (atom.subst[1] != -1 && sameside(subst_vecs[0], subst_vecs[1], vec_left) != -1)
        return false;
    // the same check for the two right substituents
    if (atom.subst[3] != -1 && sameside(subst_vecs[2], subst_vecs[3], vec_right) != -1)
        return false;

    if (dirs[0] == BOND_EITHER || dirs[1] == BOND_EITHER || dirs[2] == BOND_EITHER || dirs[3] == BOND_EITHER)
        atom.parity = 3;
    else
    {
        if (dirs[0] == 0 && dirs[1] != 0)
            dirs[0] = 3 - dirs[1];
        if (dirs[2] == 0 && dirs[3] != 0)
            dirs[2] = 3 - dirs[3];

        int ss = sameside(subst_vecs[0], subst_vecs[2], vec_right);

        if (ss == 0)
            return false;

        if (dirs[0] == 0)
            dirs[0] = (ss == 1) ? 3 - dirs[2] : dirs[2];
        else if (dirs[2] == 0)
            dirs[2] = (ss == 1) ? 3 - dirs[0] : dirs[0];

        if ((ss == 1 && dirs[0] == dirs[2]) || (ss == -1 && dirs[0] != dirs[2]))
            return false; // square-planar configuration?

        if ((ss == 1 && dirs[0] == BOND_UP) || (ss == -1 && dirs[0] == BOND_DOWN))
            atom.parity = 1;
        else
            atom.parity = 2;
    }

    const Vertex& v_left = mol.getVertex(atom.left);
    const Vertex& v_right = mol.getVertex(atom.right);

    // mark bonds as sensible
    for (k = 0, j = v_left.neiBegin(); j != v_left.neiEnd(); j = v_left.neiNext(j))
    {
        int dir = mol.getBondDirection2(atom.left, v_left.neiVertex(j));
        if (dir != 0)
            sensible_bonds_out[v_left.neiEdge(j)] = 1;
    }
    for (k = 0, j = v_right.neiBegin(); j != v_right.neiEnd(); j = v_right.neiNext(j))
    {
        int dir = mol.getBondDirection2(atom.right, v_right.neiVertex(j));
        if (dir != 0)
            sensible_bonds_out[v_right.neiEdge(j)] = 1;
    }

    // "either" allene centers do not count
    if (atom.parity == 3)
        return false;

    Vec3f prod;

    prod.cross(vec_left, subst_vecs[0]);

    if (prod.z > 0)
        atom.parity = 3 - atom.parity;

    // move hydrogens from [0] and [2] to [1] and [3] respectively
    if (pure_h[0])
    {
        std::swap(atom.subst[0], atom.subst[1]);
        atom.parity = 3 - atom.parity;
    }
    if (pure_h[2])
    {
        std::swap(atom.subst[2], atom.subst[3]);
        atom.parity = 3 - atom.parity;
    }

    return true;
}