Direct products in python

[phockett]

Hi @mcodev31,

Hopefully a quick (and possibly stupid...) question: is there a simply way to compute direct products with the python bindings?

I feel like this should be simple, but I can't quite get a handle on it, so if you have an easy way to do it with libmsym, or pointers for how else to go about it (since it should be a well-solved problem) that'd be great - I just can't find a nice implementation in python and don't really want to work from scratch unless absolutely necessary! And I hopefully haven't missed anything obvious... my only real experience of this in practice is quite old school, with pen and paper and using direct product tables.

I've played around a bit, and I think I also can get something close to what I need in the general case by doing a suitable multiplication of my basis functions and then pulling out the irreps (basically as per the wavefunction symmetrization case in the libmsym python demo script), but this isn't very elegant.

For background: I'm interested in computing multiple direct products to determine allowed cases in photoionization problems (see, e,g. https://epsproc.readthedocs.io/en/latest/ePS_ePSproc_tutorial/ePS_tutorial_080520.html#Symmetry-selection-rules-in-photoionization).

Thanks, as ever, for any help you can provide!

[mcodev31]

Hi @phockett,

There is no direct product accessible from python (there is internally in c, but I don't think it generates what you are after). Correct me if I'm wrong but you are looking for something like B₁xB₂=A₂ in C_4v (and looking for the cases where that has a component in A₁)? In that case you should be able to calculate that directly from the character table, then just take the "dot product" with A₁ (just 1s). Note however that the character table is not the character of each operation, but has the general structure of a character table (e.g. 2C₄ rather than C₄ and C₄³) so the "dot" product is not just a normal dot product. Then check if the result is orthogonal to A₁

[mcodev31]

Also note that since libmsym currently only works with RSH certain point groups will have reducible representations where the characters otherwise would have been complex.

[mcodev31]

Hopefully this is useful:

import libmsym as msym, numpy as np

def prod(character_table,s1,s2):
    r = np.zeros(len(character_table.symmetry_operations));
    for i in range(len(character_table.symmetry_operations)):
        r[i] = character_table.table[s1][i]*character_table.table[s2][i]
    return r;

def dot(character_table,s,p):
    r = 0
    l = 0
    for i in range(len(character_table.symmetry_operations)):
        l += character_table.class_count[i];
        r += p[i]*character_table.table[s][i]*character_table.class_count[i]
    return int(round(r))//l;

with msym.Context(elements=[msym.Element(name = "H", coordinates = [.0,.0,.0])], basis_functions=[], point_group="C4v") as ctx:
    for (i,s1) in enumerate(ctx.character_table.symmetry_species):
        for (j,s2) in enumerate(ctx.character_table.symmetry_species):
            p = prod(ctx.character_table,i,j);
            print(s1.name+'x'+s2.name+'=',end='')
            s = ''
            for (k,s3) in enumerate(ctx.character_table.symmetry_species):
                d = dot(ctx.character_table,k,p)
                if(d != 0):
                    if s != '':
                        s += '+'
                    s += str(d) + s3.name
            print(s)

[phockett]

Yes, great, that's exactly what I was trying to get too - but I was struggling to get from the standard group theory formulas to numerical implementation. Thanks very much for the help, you've almost made it look easy ;)

I plan to wrap this into a larger routine for the photoionization case I mentioned above, so will update here when I have something useful to share back.

[mcodev31]

Great, always fun to see how this is used =)

[phockett]

Still a bit of a work-in-progress (as is ever the case!), and messy (as is also ever the case!), but I now have a few wrappers/extended/modified routines in place.

In case you are interested:

• Demo doc for the routines. This shows use for the photoionization case I mentioned above, specifically determining dipole symmetries and allowed photoionization symmetries, with associated matrix elements, based on the libmsym routines.
• Direct product routines, basically as per your comment above plus extensions to multiple-product cases and some vectorization: https://github.com/phockett/PEMtk/blob/master/pemtk/sym/_directProduct.py (these are wrapped for my class, but should be generally useful).
• Routines to handle dipole symmetries and associated selection rules: https://github.com/phockett/PEMtk/blob/master/pemtk/sym/_dipoleTerms.py (probably not very generally useful, these are mainly for working with photoionization matrix elements).

And thanks again for the help - I wouldn't have got most of this in place without libmsym, plus your additional comments!

[mcodev31]

Looks very nice