Tracking atom index (ID) for equivalence sets

[ghutchis]

I implemented display of equivalent sets of atoms in Avogadro v2:

In the process, I realized that once you get the sets, there's no way to relate the atoms in the group back to the original molecule, except by element and coordinates.

Right now, I'm doing that, but it would be great to have an atom id (index) that travels through the context. Is there an easy way to do that?

  const msym_equivalence_set_t* smes = &m_es[group];
  const msym_element_t* a = smes->elements[atomInGroup];
  if (a == NULL)
    return;

  unsigned int length = m_molecule->atomCount();
  for (Index i = 0; i < length; ++i) {
    qDebug() << "checking atom" << i << " for " << a->n;
    m_molecule->setAtomSelected(i, false);
    if (m_molecule->atomicNumbers()[i] != a->n)
      continue;

    Vector3 ipos = m_molecule->atomPositions3d()[i];
    if (fabs(a->v[0] - (ipos[0] - m_cm[0])) < 0.05 &&
        fabs(a->v[0] - (ipos[0] - m_cm[0])) < 0.05 &&
        fabs(a->v[0] - (ipos[0] - m_cm[0])) < 0.05) {
      m_molecule->setAtomSelected(i, true);
      qDebug() << " got one!";
    }
  }

In other words, what I'd really like to do is:

  const msym_equivalence_set_t* smes = &m_es[group];
  const msym_element_t* a = smes->elements[atomInGroup];
  const int index = a->id;

[mcodev31]

Hi

Each element has an id which is intended for cases like these.

    typedef struct _msym_element {
        void *id;                               // custom identifier
        double m;                               // Mass
        double v[3];                            // Position
        int n;                                  // Nuclear charge
        char name[4];                           // Name
    } msym_element_t;

Or is this not what you meant?

[ghutchis]

I did see the id field, but not an example on how it's used. Am I supposed to supply those when handing off the molecule? I worry because in equivalence_set.c, there's this:

libmsym/src/equivalence_set.c

Line 87 in c994703

ge[gel].id = NULL;

At the outset, it seems like the id fields are cleared.

[mcodev31]

Yes, you supply it in each element when doingmsymSetElements.
That code only runs from msymGenerateElements in which case we want to remove the id for newly generated elements.
I.e. if you want to generate a new molecule from an atom and a group the new elements will not have an id (but the original will).

[mcodev31]

Not been doing much chemistry lately, but just downloaded Avogadro 2 from Flathub. It looks great! I can see that atom selection is working. I like the subgroup selection and rendering as well. Perhaps we should look at rendering some orbitals at some point. Would probably be a good educational feature to be able to quickly view (example) SALCs for some n.

[ghutchis]

That could be interesting (e.g., going from s or p orbitals on a particular atom to the SALC).

[mcodev31]

Shouldn't be too difficult assuming one could add an orbital to an atom (or equivalence set).
You could also show the symmetry species of all orbitals in e.g. a fchk, and/or vibrational modes.

[mcodev31]

I made a (very quick and ugly) hack to the Avogadro symmetry code to see if I could extract the symmetry species of the orbitals in an fchk file for CH2 (C2v) and basis C=1s2s2p, H=1s. It spits out this:

0.992799 A1
0.876721 A1
0.815339 B1
1.00216 A1
1 B2
1.8554 A1
1.701 B1

Which matches 4A1 + 0A2 + 2B1 + 1B2 expected by libmsym using that basis (code above does not print species with coefficient < epsilon (double))
So one could imagine adding a tab listing all orbitals and their symmetry species and perhaps another with the character table.
SALC generation should just reasonably similar, just in reverse.
Perhaps we should continue this conversation by mail.

[ghutchis]

This would be great, particularly since many calculations run in low symmetry.

I'd be happy to talk over email.

[mcodev31]

I've sent you a couple of patches (hopefully to the right email).