Hi Dave and everyone,
Thanks to your advice, the code below worked to get corresponding bonds.
// mol1: ground truth loaded from SDF
// mol2: predicted molecule structure loaded from SDF
OBMol mol1, mol2;
// Get mapping between two molecules
OBQuery *query = CompileMoleculeQuery(&mol1);
OBIsomorphismMapper *mapper = OBIsomorphismMapper::GetInstance(query);
OBIsomorphismMapper::Mappings mappings;
mapper->MapAll(&mol2, mappings);
map<int, int> idxMapping; // mol1 index to mol2 index
for(size_t m = 0; m < mappings.size(); ++m) {
bool isCorrect = true;
for(size_t i = 0; i < mappings[m].size(); ++i) {
//obmol indices are 1-indexed while the mapper is zero indexed
unsigned int idx1 = mappings[m][i].first + 1;
unsigned int idx2 = mappings[m][i].second + 1;
idxMapping[idx1] = idx2;
if(mol1.GetAtom(idx1)->GetAtomicNum() !=
mol2.GetAtom(idx2)->GetAtomicNum()) {
idxMapping.clear();
isCorrect = false;
break;
}
}
if(isCorrect) break;
}
// Display the length of each bond
FOR_BONDS_OF_MOL(bond1, mol1) {
unsigned int begin1 = idxMapping[bond1->GetBeginAtomIdx()];
unsigned int end1 = idxMapping[bond1->GetEndAtomIdx()];
FOR_BONDS_OF_MOL(bond2, mol2) {
unsigned int begin2 = bond2->GetBeginAtomIdx();
unsigned int end2 = bond2->GetEndAtomIdx();
if((begin1 == begin2 && end1 == end2)) {
cout << bond1->GetLength() << "," << bond2->GetLength() << endl;
break;
}
}
}
I'm wondering why the following code using MapFirst to get idxMapping
did not work.
I assume that OBIsomorphism only considers connection, but is this
behavior intended?
// Get mapping between two molecules
OBQuery *query = CompileMoleculeQuery(&mol1);
OBIsomorphismMapper *mapper = OBIsomorphismMapper::GetInstance(query);
OBIsomorphismMapper::Mapping mapping;
mapper->MapFirst(&mol2, mapping);
map<int, int> idxMapping; // mol1 index to mol2 index
for(size_t i = 0; i < mapping.size(); ++i) {
//obmol indices are 1-indexed while the mapper is zero indexed
unsigned int idx1 = mapping[i].first + 1;
unsigned int idx2 = mapping[i].second + 1;
idxMapping[idx1] = idx2;
// this assertion fails in some molecules
assert(mol1.GetAtom(idx1)->GetAtomicNum() !=
mol2.GetAtom(idx2)->GetAtomicNum());
}
Thanks,
Naruki
2018年10月14日(日) 21:33 Koes, David <[email protected]>:
>
> Try using OBIsomorphismMapper to find tge right mapping of atoms which will
> implicitly give you the right matching for bonds. For example:
> https://github.com/dkoes/openbabel/blob/master/tools/obrms.cpp
>
> -Dave
>
> On Oct 14, 2018, at 4:01 AM, Naruki Yoshikawa <[email protected]>
> wrote:
>
> Hi everyone,
>
> I'm evaluating my fragment-based coordinate generation method.
> I want to measure the error of bond lengths and angles, but I have a
> problem in getting corresponding ones from two molecules.
>
> What I want to do is like this:
>
> // mol1: ground truth loaded from SDF
> // mol2: predicted molecule structure loaded from SDF
> OBMol mol1, mol2;
>
> // Display the length of each bond
> for (int idx = 0; idx < mol1.NumBonds(); ++idx) {
> OBBond *bond1 = mol1.GetBond(idx);
> OBBond *bond2 = mol2.GetBond(idx);
>
> // these assertion fails
> assert(bond1->GetBeginAtom()->GetAtomicNum()
> == bond2->GetBeginAtom()->GetAtomicNum());
> assert(bond1->GetEndAtom()->GetAtomicNum()
> == bond2->GetEndAtom()->GetAtomicNum());
>
> cout << bond1->GetLength() << "," << bond2->GetLength() << endl;
> }
>
> The problem is that the bonds are not ordered properly.
> When I execute the code above, the assertions fail.
> So n-th bond in mol1 was different from n-th bond in mol2.
> I also have a similar problem for angles.
>
> If I can get something like "canonical order of bonds", I can compare
> the length of each bond correctly.
> However, I don't know how to do it.
>
> Does anyone have a good idea to get corresponding bonds and angles?
>
> Thanks,
> Naruki
>
>
> _______________________________________________
> OpenBabel-Devel mailing list
> [email protected]
> https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.sourceforge.net%2Flists%2Flistinfo%2Fopenbabel-devel&data=02%7C01%7Cdkoes%40pitt.edu%7Ca4d3f190d0c74200797d08d631ab31a3%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C636751008682354300&sdata=6%2FgdvGOjbSNR5zoo20k4MPJVuZn3057ooWnoFpHo1rI%3D&reserved=0
_______________________________________________
OpenBabel-Devel mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/openbabel-devel
