in gsoc2022/seagrid-rich-client/molview/src/js/chem/chem/stereocenters.js [226:563]
chem.Stereocenters.prototype._buildOneCenter = function (atom_idx /*, int group, int type, const int *bond_orientations*/ )
{
var atom = this.molecule.atoms.get(atom_idx);
var nei_list = this.getNeighbors.call(this.context, atom_idx);
var degree = nei_list.length;
var implicit_degree = -1;
var stereocenter = {
group: 0, // = group;
type: 0, // = type;
pyramid: new Array(4)
};
var nei_idx = 0;
var edge_ids = new Array(4);
var last_atom_dir = 0;
var n_double_bonds = 0;
stereocenter.pyramid[0] = -1;
stereocenter.pyramid[1] = -1;
stereocenter.pyramid[2] = -1;
stereocenter.pyramid[3] = -1;
var n_pure_hydrogens = 0;
if(degree > 4)
throw new Error("Stereocenter with " + degree + " bonds are not supported");
nei_list.each(function (nei)
{
var nei_atom = this.molecule.atoms.get(nei.aid);
var bond = this.molecule.bonds.get(nei.bid);
edge_ids[nei_idx] = {
edge_idx: nei.bid,
nei_idx: nei.aid,
rank: nei.aid,
vec: util.Vec2.diff(nei_atom.pp, atom.pp).yComplement()
};
if(nei_atom.pureHydrogen())
{
n_pure_hydrogens++;
edge_ids[nei_idx].rank = 10000;
}
else if(nei_atom.label == 'H')
edge_ids[nei_idx].rank = 5000;
if(!edge_ids[nei_idx].vec.normalize())
throw new Error("Zero bond length");
if(bond.type == chem.Struct.BOND.TYPE.TRIPLE)
throw new Error("Non-single bonds not allowed near stereocenter");
else if(bond.type == chem.Struct.BOND.TYPE.AROMATIC)
throw new Error("Aromatic bonds not allowed near stereocenter");
else if(bond.type == chem.Struct.BOND.TYPE.DOUBLE)
n_double_bonds++;
nei_idx++;
}, this);
util.find(chem.Stereocenters.allowed_stereocenters, function (as)
{
if(as.elem == atom.label && as.charge == atom.charge &&
as.degree == degree && as.n_double_bonds == n_double_bonds)
{
implicit_degree = as.implicit_degree;
return true;
}
return false;
}, this);
if(implicit_degree == -1)
throw new Error("Unknown stereocenter configuration: " + atom.label + ", charge " + atom.charge + ", " + degree + " bonds (" + n_double_bonds + " double)");
if(degree == 4 && n_pure_hydrogens > 1)
throw new Error(n_pure_hydrogens + " hydrogens near stereocenter");
if(degree == 3 && implicit_degree == 4 && n_pure_hydrogens > 0)
throw new Error("Has hydrogen(s) besides implicit hydrogen near stereocenter");
/*
if (stereocenter.type == ATOM_ANY)
{
_stereocenters.insert(atom_idx, stereocenter);
return;
}
*/
if(degree == 4)
{
// sort by neighbor atom index (ascending)
if(edge_ids[0].rank > edge_ids[1].rank)
edge_ids.swap(0, 1);
if(edge_ids[1].rank > edge_ids[2].rank)
edge_ids.swap(1, 2);
if(edge_ids[2].rank > edge_ids[3].rank)
edge_ids.swap(2, 3);
if(edge_ids[1].rank > edge_ids[2].rank)
edge_ids.swap(1, 2);
if(edge_ids[0].rank > edge_ids[1].rank)
edge_ids.swap(0, 1);
if(edge_ids[1].rank > edge_ids[2].rank)
edge_ids.swap(1, 2);
var main1 = -1,
main2 = -1,
side1 = -1,
side2 = -1;
var main_dir = 0;
for(nei_idx = 0; nei_idx < 4; nei_idx++)
{
var stereo = this._getBondStereo(atom_idx, edge_ids[nei_idx].edge_idx);
if(stereo == chem.Struct.BOND.STEREO.UP || stereo == chem.Struct.BOND.STEREO.DOWN)
{
main1 = nei_idx;
main_dir = stereo;
break;
}
}
if(main1 == -1)
throw new Error("None of 4 bonds going from stereocenter is stereobond");
var xyz1, xyz2;
// find main2 as opposite to main1
if(main2 == -1)
{
xyz1 = chem.Stereocenters._xyzzy(edge_ids[main1].vec, edge_ids[(main1 + 1) % 4].vec, edge_ids[(main1 + 2) % 4].vec);
xyz2 = chem.Stereocenters._xyzzy(edge_ids[main1].vec, edge_ids[(main1 + 1) % 4].vec, edge_ids[(main1 + 3) % 4].vec);
if(xyz1 + xyz2 == 3 || xyz1 + xyz2 == 12)
{
main2 = (main1 + 1) % 4;
side1 = (main1 + 2) % 4;
side2 = (main1 + 3) % 4;
}
}
if(main2 == -1)
{
xyz1 = chem.Stereocenters._xyzzy(edge_ids[main1].vec, edge_ids[(main1 + 2) % 4].vec, edge_ids[(main1 + 1) % 4].vec);
xyz2 = chem.Stereocenters._xyzzy(edge_ids[main1].vec, edge_ids[(main1 + 2) % 4].vec, edge_ids[(main1 + 3) % 4].vec);
if(xyz1 + xyz2 == 3 || xyz1 + xyz2 == 12)
{
main2 = (main1 + 2) % 4;
side1 = (main1 + 1) % 4;
side2 = (main1 + 3) % 4;
}
}
if(main2 == -1)
{
xyz1 = chem.Stereocenters._xyzzy(edge_ids[main1].vec, edge_ids[(main1 + 3) % 4].vec, edge_ids[(main1 + 1) % 4].vec);
xyz2 = chem.Stereocenters._xyzzy(edge_ids[main1].vec, edge_ids[(main1 + 3) % 4].vec, edge_ids[(main1 + 2) % 4].vec);
if(xyz1 + xyz2 == 3 || xyz1 + xyz2 == 12)
{
main2 = (main1 + 3) % 4;
side1 = (main1 + 2) % 4;
side2 = (main1 + 1) % 4;
}
}
if(main2 == -1)
throw new Error("Internal error: can not find opposite bond");
if(main_dir == chem.Struct.BOND.STEREO.UP && this._getBondStereo(atom_idx, edge_ids[main2].edge_idx) == chem.Struct.BOND.STEREO.DOWN)
throw new Error("Stereo types of the opposite bonds mismatch");
if(main_dir == chem.Struct.BOND.STEREO.DOWN && this._getBondStereo(atom_idx, edge_ids[main2].edge_idx) == chem.Struct.BOND.STEREO.UP)
throw new Error("Stereo types of the opposite bonds mismatch");
if(main_dir == this._getBondStereo(atom_idx, edge_ids[side1].edge_idx))
throw new Error("Stereo types of non-opposite bonds match");
if(main_dir == this._getBondStereo(atom_idx, edge_ids[side2].edge_idx))
throw new Error("Stereo types of non-opposite bonds match");
if(main1 == 3 || main2 == 3)
last_atom_dir = main_dir;
else
last_atom_dir = (main_dir == chem.Struct.BOND.STEREO.UP ? chem.Struct.BOND.STEREO.DOWN : chem.Struct.BOND.STEREO.UP);
sign = chem.Stereocenters._sign(edge_ids[0].vec, edge_ids[1].vec, edge_ids[2].vec);
if((last_atom_dir == chem.Struct.BOND.STEREO.UP && sign > 0) ||
(last_atom_dir == chem.Struct.BOND.STEREO.DOWN && sign < 0))
{
stereocenter.pyramid[0] = edge_ids[0].nei_idx;
stereocenter.pyramid[1] = edge_ids[1].nei_idx;
stereocenter.pyramid[2] = edge_ids[2].nei_idx;
}
else
{
stereocenter.pyramid[0] = edge_ids[0].nei_idx;
stereocenter.pyramid[1] = edge_ids[2].nei_idx;
stereocenter.pyramid[2] = edge_ids[1].nei_idx;
}
stereocenter.pyramid[3] = edge_ids[3].nei_idx;
}
else if(degree == 3)
{
// sort by neighbor atom index (ascending)
if(edge_ids[0].rank > edge_ids[1].rank)
edge_ids.swap(0, 1);
if(edge_ids[1].rank > edge_ids[2].rank)
edge_ids.swap(1, 2);
if(edge_ids[0].rank > edge_ids[1].rank)
edge_ids.swap(0, 1);
var stereo0 = this._getBondStereo(atom_idx, edge_ids[0].edge_idx);
var stereo1 = this._getBondStereo(atom_idx, edge_ids[1].edge_idx);
var stereo2 = this._getBondStereo(atom_idx, edge_ids[2].edge_idx);
var n_up = 0,
n_down = 0;
n_up += ((stereo0 == chem.Struct.BOND.STEREO.UP) ? 1 : 0);
n_up += ((stereo1 == chem.Struct.BOND.STEREO.UP) ? 1 : 0);
n_up += ((stereo2 == chem.Struct.BOND.STEREO.UP) ? 1 : 0);
n_down += ((stereo0 == chem.Struct.BOND.STEREO.DOWN) ? 1 : 0);
n_down += ((stereo1 == chem.Struct.BOND.STEREO.DOWN) ? 1 : 0);
n_down += ((stereo2 == chem.Struct.BOND.STEREO.DOWN) ? 1 : 0);
if(implicit_degree == 4) // have implicit hydrogen
{
if(n_up == 3)
throw new Error("All 3 bonds up near stereoatom");
if(n_down == 3)
throw new Error("All 3 bonds down near stereoatom");
if(n_up == 0 && n_down == 0)
throw new Error("No up/down bonds near stereoatom -- indefinite case");
if(n_up == 1 && n_down == 1)
throw new Error("One bond up, one bond down -- indefinite case");
main_dir = 0;
if(n_up == 2)
last_atom_dir = chem.Struct.BOND.STEREO.DOWN;
else if(n_down == 2)
last_atom_dir = chem.Struct.BOND.STEREO.UP;
else
{
main1 = -1;
side1 = -1;
side2 = -1;
for(nei_idx = 0; nei_idx < 3; nei_idx++)
{
dir = this._getBondStereo(atom_idx, edge_ids[nei_idx].edge_idx);
if(dir == chem.Struct.BOND.STEREO.UP || dir == chem.Struct.BOND.STEREO.DOWN)
{
main1 = nei_idx;
main_dir = dir;
side1 = (nei_idx + 1) % 3;
side2 = (nei_idx + 2) % 3;
break;
}
}
if(main1 == -1)
throw new Error("Internal error: can not find up or down bond");
var xyz = chem.Stereocenters._xyzzy(edge_ids[side1].vec, edge_ids[side2].vec, edge_ids[main1].vec);
if(xyz == 3 || xyz == 4)
throw new Error("Degenerate case for 3 bonds near stereoatom");
if(xyz == 1)
last_atom_dir = main_dir;
else
last_atom_dir = (main_dir == chem.Struct.BOND.STEREO.UP ? chem.Struct.BOND.STEREO.DOWN : chem.Struct.BOND.STEREO.UP);
}
var sign = chem.Stereocenters._sign(edge_ids[0].vec, edge_ids[1].vec, edge_ids[2].vec);
if((last_atom_dir == chem.Struct.BOND.STEREO.UP && sign > 0) ||
(last_atom_dir == chem.Struct.BOND.STEREO.DOWN && sign < 0))
{
stereocenter.pyramid[0] = edge_ids[0].nei_idx;
stereocenter.pyramid[1] = edge_ids[1].nei_idx;
stereocenter.pyramid[2] = edge_ids[2].nei_idx;
}
else
{
stereocenter.pyramid[0] = edge_ids[0].nei_idx;
stereocenter.pyramid[1] = edge_ids[2].nei_idx;
stereocenter.pyramid[2] = edge_ids[1].nei_idx;
}
stereocenter.pyramid[3] = -1;
}
else // 3-connected P, N or S; no implicit hydrogens
{
var dir;
if(n_down > 0 && n_up > 0)
throw new Error("One bond up, one bond down -- indefinite case");
else if(n_down == 0 && n_up == 0)
throw new Error("No up-down bonds attached to stereocenter");
else if(n_up > 0)
dir = 1;
else
dir = -1;
if(chem.Stereocenters._xyzzy(edge_ids[0].vec, edge_ids[1].vec, edge_ids[2].vec) == 1 ||
chem.Stereocenters._xyzzy(edge_ids[0].vec, edge_ids[2].vec, edge_ids[1].vec) == 1 ||
chem.Stereocenters._xyzzy(edge_ids[2].vec, edge_ids[1].vec, edge_ids[0].vec) == 1)
// all bonds belong to the same half-plane
dir = -dir;
sign = chem.Stereocenters._sign(edge_ids[0].vec, edge_ids[1].vec, edge_ids[2].vec);
if(sign == dir)
{
stereocenter.pyramid[0] = edge_ids[0].nei_idx;
stereocenter.pyramid[1] = edge_ids[2].nei_idx;
stereocenter.pyramid[2] = edge_ids[1].nei_idx;
}
else
{
stereocenter.pyramid[0] = edge_ids[0].nei_idx;
stereocenter.pyramid[1] = edge_ids[1].nei_idx;
stereocenter.pyramid[2] = edge_ids[2].nei_idx;
}
stereocenter.pyramid[3] = -1;
}
}
this.atoms.set(atom_idx, stereocenter);
};