in core/indigo-core/reaction/src/reaction_enumerator_state.cpp [1384:1640]
bool ReactionEnumeratorState::_attachFragments(Molecule& ready_product_out, Array<int>& ucfrag_mapping)
{
QS_DEF(Array<int>, frags2product_mapping);
_findFragments2ProductMapping(frags2product_mapping);
QS_DEF(QueryMolecule, product);
product.clear();
product.clone(_full_product, NULL, NULL);
QS_DEF(Molecule, uncleaned_fragments);
uncleaned_fragments.clear();
uncleaned_fragments.clone(_fragments, NULL, NULL);
QS_DEF(Molecule, mol_product);
mol_product.clear();
QS_DEF(Array<int>, mapping);
mapping.clear();
QS_DEF(Array<int>, all_forbidden_atoms);
all_forbidden_atoms.clear();
all_forbidden_atoms.clear_resize(product.vertexEnd() + _fragments.vertexCount());
all_forbidden_atoms.zerofill();
for (int i = product.vertexBegin(); i != product.vertexEnd(); i = product.vertexNext(i))
all_forbidden_atoms[i] = 1;
_buildMolProduct(_full_product, mol_product, uncleaned_fragments, all_forbidden_atoms, mapping);
_cleanFragments();
QS_DEF(Array<int>, frags_mapping);
frags_mapping.clear_resize(_fragments.vertexEnd());
frags_mapping.fffill();
mol_product.mergeWithMolecule(_fragments, &frags_mapping);
for (int i = _fragments.vertexBegin(); i < _fragments.vertexEnd(); i = _fragments.vertexNext(i))
if (i < _monomer_forbidden_atoms.size() && _monomer_forbidden_atoms[i])
all_forbidden_atoms[frags_mapping[i]] = _monomer_forbidden_atoms[i];
QS_DEF(Array<int>, product_mapping);
product_mapping.clear_resize(_full_product.vertexEnd());
for (int i = 0; i < _full_product.vertexEnd(); i++)
product_mapping[i] = i;
for (int i = _full_product.vertexBegin(); i != _full_product.vertexEnd(); i = _full_product.vertexNext(i))
{
if (_att_points[i].size() == 0)
continue;
const Vertex& pr_v = mol_product.getVertex(i);
QS_DEF(Array<int>, pr_neibours);
pr_neibours.clear();
for (int j = pr_v.neiBegin(); j != pr_v.neiEnd(); j = pr_v.neiNext(j))
pr_neibours.push(pr_v.neiVertex(j));
if (_is_rg_exist && (pr_neibours.size() == 2))
for (int j = 0; j < pr_neibours.size(); j++)
if (_product_aam_array[pr_neibours[j]] == 0)
throw Error("There are no AAM on RGroups attachment points");
if (_is_rg_exist)
{
if (pr_neibours.size() > 2)
throw Error("RGroup atom can't have more than two neighbors");
/* Setting the order of rgroup atom neighbors by AAM (less is first) */
if (pr_neibours.size() == 2)
{
if (((_product_aam_array.size() > pr_neibours[0]) && (_product_aam_array.size() > pr_neibours[1])) &&
_product_aam_array[pr_neibours[0]] > _product_aam_array[pr_neibours[1]])
{
int tmp = pr_neibours[0];
pr_neibours[0] = pr_neibours[1];
pr_neibours[1] = tmp;
}
}
}
bool is_valid = false;
if (is_transform && _att_points[i].size() != 0 && _checkValence(mol_product, frags_mapping[_att_points[i][0]]))
is_valid = true;
if (_is_rg_exist)
{
for (int j = 0; j < pr_neibours.size(); j++)
{
if (mol_product.findEdgeIndex(pr_neibours[j], frags_mapping[_att_points[i][j]]) != -1)
return false;
int atom_from = mapping[i];
int atom_to = frags_mapping[_att_points[i][j]];
if (mol_product.stereocenters.exists(atom_from) && mol_product.stereocenters.getPyramid(atom_from)[3] == -1)
return false;
if (mol_product.stereocenters.exists(atom_to) && mol_product.stereocenters.getPyramid(atom_to)[3] != -1)
return false;
mol_product.flipBond(pr_neibours[j], atom_from, atom_to);
// TODO:
// Check that corresponding R-group fragment in monomer has cis-trans bond
// and check that AAM mapping is specified for that.
// For example for reaction OC([*])=O>>OC([*])=O and monomer C\C=C\C(O)=O
// product shouldn't have should have cis-trans bonds because
// AAM is not specified on R-group atom neighbor
// Cis-trans bonds should be saved for such reaction: O[C:1]([*])=O>>O[C:1]([*])=O
}
mol_product.removeAtom(mapping[i]);
}
else
{
for (int j = 0; j < _att_points[i].size(); j++)
{
int mon_atom = frags_mapping[_att_points[i][j]];
int pr_atom = mapping[i];
const Vertex& mon_v = mol_product.getVertex(mon_atom);
// const Vertex& pr_v = mol_product.getVertex(pr_atom);
for (int k = mon_v.neiBegin(); k != mon_v.neiEnd(); k = mon_v.neiNext(k))
if (MoleculeCisTrans::isGeomStereoBond(mol_product, mon_v.neiEdge(k), NULL, false))
mol_product.cis_trans.setParity(mon_v.neiEdge(k), 0);
if (mol_product.stereocenters.exists(mon_atom))
mol_product.stereocenters.remove(mon_atom);
QS_DEF(Array<int>, neighbors);
neighbors.clear();
for (int k = mon_v.neiBegin(); k != mon_v.neiEnd(); k = mon_v.neiNext(k))
neighbors.push(mon_v.neiVertex(k));
for (int k = 0; k < neighbors.size(); k++)
if (mol_product.findEdgeIndex(neighbors[k], pr_atom) == -1)
mol_product.flipBond(neighbors[k], mon_atom, pr_atom);
frags_mapping[_att_points[i][j]] = pr_atom;
mol_product.removeAtom(mon_atom);
// if (mol_product.mergeAtoms(frags_mapping[_att_points[i][0]], mapping[i]) == -1)
// return false;
}
}
product_mapping[mapping[i]] = frags_mapping[_att_points[i][0]];
/*
if (is_transform && _att_points[i].size() != 0 && is_valid && !_checkValence(mol_product, mapping[i]))
{
_product_forbidden_atoms.copy(_monomer_forbidden_atoms);
_product_forbidden_atoms[_att_points[i][0]] = max_reuse_count;
return false;
}*/
/* Border stereocenters copying */
int nv_idx = 0;
for (int j = 0; j < _att_points[i].size(); j++)
{
if (uncleaned_fragments.stereocenters.exists(_att_points[i][j]) && !mol_product.stereocenters.exists(frags_mapping[_att_points[i][j]]))
{
int type, group, pyramid[4];
uncleaned_fragments.stereocenters.get(_att_points[i][j], type, group, pyramid);
int new_pyramid[4];
bool invalid_stereocenter = false;
for (int k = 0; k < 4; k++)
{
if (pyramid[k] == -1)
new_pyramid[k] = -1;
else if (!_is_needless_atom[pyramid[k]])
new_pyramid[k] = frags_mapping[pyramid[k]];
else if (frags2product_mapping[pyramid[k]] != -1)
{
new_pyramid[k] = frags2product_mapping[pyramid[k]];
}
else
{
invalid_stereocenter = true;
break;
}
}
if (!invalid_stereocenter)
mol_product.addStereocenters(frags_mapping[_att_points[i][j]], type, group, new_pyramid);
}
if (nv_idx == 2)
break;
}
}
/* Updating of cis-trans information on product & monomer's fragment border */
_completeCisTrans(mol_product, uncleaned_fragments, frags_mapping);
QS_DEF(Array<int>, out_mapping);
out_mapping.clear_resize(mol_product.vertexEnd());
ready_product_out.clone(mol_product, NULL, &out_mapping);
_product_forbidden_atoms.clear_resize(ready_product_out.vertexEnd());
_product_forbidden_atoms.zerofill();
QS_DEF(Array<int>, temp_orig_hydr);
temp_orig_hydr.clear();
if (is_transform)
{
for (int i = mol_product.vertexBegin(); i != mol_product.vertexEnd(); i = mol_product.vertexNext(i))
if (out_mapping[i] != -1 && all_forbidden_atoms[i])
_product_forbidden_atoms[out_mapping[i]] = all_forbidden_atoms[i];
for (int i = 0; i < _original_hydrogens.size(); i++)
{
int new_h_idx = frags_mapping[_original_hydrogens[i]];
if (new_h_idx == -1)
continue;
temp_orig_hydr.push(out_mapping[new_h_idx]);
}
_original_hydrogens.copy(temp_orig_hydr);
}
ucfrag_mapping.clear_resize(_fragments.vertexEnd());
ucfrag_mapping.fffill();
QS_DEF(Array<int>, old_mapping);
old_mapping.copy(_mapping);
_mapping.clear_resize(_fragments.vertexEnd());
_mapping.fffill();
for (int i = uncleaned_fragments.vertexBegin(); i != uncleaned_fragments.vertexEnd(); i++)
{
if (frags_mapping[i] != -1)
ucfrag_mapping[i] = frags_mapping[i];
else if (frags2product_mapping[i] != -1)
ucfrag_mapping[i] = frags2product_mapping[i];
else
continue;
ucfrag_mapping[i] = out_mapping[ucfrag_mapping[i]];
if (old_mapping.size() > 0)
{
int i_id = old_mapping.find(i);
if ((i_id != -1) && (i_id < _mapping.size()))
_mapping[i_id] = ucfrag_mapping[i];
}
else
_mapping[i] = ucfrag_mapping[i];
}
if (refine_proc)
return refine_proc(uncleaned_fragments, ready_product_out, ucfrag_mapping, userdata);
return true;
}