Generating symmetrized harmonics for a given point group

phockett commented 3 years ago

Hi,

Possibly a silly question (but hopefully not) - is it possible to use libmsym to generate symmetrized harmonics for a single-point type expansion^? I think this should be very similar to multi-point SALC generation, and I tried to do it based on the example code (in python) by setting a single dummy core system with multiple Ylms, but couldn't get it to work. Specifically, I tried something like the following (and a few variations on basis set and point group):

# Set single element
elements = [msym.Element(name = "C", coordinates = [0.0, 0.0,0.0])]

# Test basis for arb l,m,n...
def set_basis(element, l, m, n, name):
    basis_function = msym.RealSphericalHarmonic(element = element, l=l, m=m, n=n, name = name)
    element.basis_functions = [basis_function]
    return basis_function

basis_functions = []
for l in range(0,5):
    for m in range(-l,l+1):
        basis_functions.extend([set_basis(e, l=l, m=m, n=l, name=f"{l}s{m}") for e in elements])

point_group = "C2v"

ctx = msym.Context(elements = elements, basis_functions = basis_functions)
ctx.point_group = point_group
ctx.symmetrize_elements()
ctx.salcs

This fails at SALC generation with Error: 'Invalid basis functions': 'Number of generated SALC wavefunctions (0) does not match orbital basis (25)', which I suspect is due to the single core (although I may have done something else silly, since I'm new to the library).

Is this possible, or am I barking up the wrong tree?

Thanks in advance (and thanks for the great library too!).

^ Specifically I'm working on problems in photoionization, where the continuum wavefunctions are described by a set of symmetrized harmonics, for the given point group of the ionizing system.

mcodev31 commented 3 years ago

Hi,

The problem is the principal quantum number n=l.

Note that these are real spherical harmonics so the are linear combinations of the magnetic quantum number (e.g. 2pz = 0, 2px = 1, 2py = -1). The coefficients won't be that interesting, since each orbital will simply belong to a symmetry species or be part of a degenerate pair (or higher dim).

This page would probably be useful if you just want to look things up rather than generate: Character Tables for Point Groups used in Chemistry

This should print something reasonably useful for C3v (I rarely use python so you'll have to excuse the code)

n = 5
basis_functions = []
for l in range(0,n):
    for m in range(-l,l+1):
        basis_functions.extend([set_basis(e, l=l, m=m, n=n, name=f"{l}s{m}") for e in elements])

ctx = msym.Context(elements = elements, basis_functions = basis_functions, point_group = "C3v")
for ss in ctx.subrepresentation_spaces:
    print(ctx.character_table.symmetry_species[ss.symmetry_species].name)
    for salcix, salc in enumerate(ss.salcs):
        print(f"salc {salcix}")
        for pfix, pf in enumerate(salc.partner_functions):
            print(f"\tpartern function {pfix}")
            for ix, v in enumerate(pf):
                if (abs(v) > 100*sys.float_info.epsilon):
                    print("\t\t{0}*{1} ".format(v,salc.basis_functions[ix].name));

phockett commented 2 years ago

Hi @mcodev31,

Many thanks for the suggestion (and apologies for the slow response), this does seem to be what I am after, or at least pretty close^! In general it's both the symmetry classification I'm after and the coefficients (although I agree they're rather trivial for the C3v example case), so a basic look-up for Symmetry of Spherical Harmonics is not what I need (although, obviously, is often pretty useful). I also need to go to relatively high L potentially (usually not more than 10 or so, but possibly as high as 30 - 40 in some cases), and didn't want to spend time reinventing the wheel here - so it's great that libmsym has this covered. :+1: :smile:

^ I still need to work out some sign and normalisation conventions, but otherwise all good I think.

As a more representative example for anyone interested, here's the output from your code snippet for point_group = "Td" (with n=5), which can be compared with literature values.

A1
salc 0
    partern function 0
        1.0*0s0 
salc 1
    partern function 0
        1.0*3s-2 
salc 2
    partern function 0
        0.7637626158259733*4s0 
        0.6454972243679029*4s4 
A2
E
salc 0
    partern function 0
        0.5000000000000003*2s0 
        0.8660254037844385*2s2 
    partern function 1
        -0.8660254037844385*2s0 
        0.5000000000000001*2s2 
salc 1
    partern function 0
        0.32274861218395134*4s0 
        -0.8660254037844388*4s2 
        -0.3818813079129865*4s4 
    partern function 1
        -0.5590169943749477*4s0 
        -0.49999999999999983*4s2 
        0.6614378277661476*4s4 
T1
salc 0
    partern function 0
        0.790569415042095*3s1 
        0.6123724356957945*3s3 
    partern function 1
        -1.0*3s2 
    partern function 2
        -0.6123724356957947*3s-3 
        0.7905694150420948*3s-1 
salc 1
    partern function 0
        0.3535533905932735*4s-3 
        0.9354143466934854*4s-1 
    partern function 1
        -1.0*4s-4 
    partern function 2
        0.9354143466934853*4s1 
        -0.35355339059327395*4s3 
T2
salc 0
    partern function 0
        1.0*1s1 
    partern function 1
        -1.0*1s0 
    partern function 2
        1.0*1s-1 
salc 1
    partern function 0
        1.0*2s-1 
    partern function 1
        -1.0*2s-2 
    partern function 2
        1.0*2s1 
salc 2
    partern function 0
        0.6123724356957945*3s1 
        -0.790569415042095*3s3 
    partern function 1
        1.0*3s0 
    partern function 2
        0.7905694150420947*3s-3 
        0.6123724356957946*3s-1 
salc 3
    partern function 0
        0.9354143466934854*4s-3 
        -0.3535533905932735*4s-1 
    partern function 1
        -1.0*4s-2 
    partern function 2
        -0.3535533905932739*4s1 
        -0.9354143466934854*4s3

And the literature results (quick test for A1, l=4, m=0,+/-4 case only):

# From Boiko et. al. (2006), table VI, should have coeffs sqrt(5/24) and sqrt(7/12)
# Looks OK aside from extra sqrt(2) norm.
print(np.sqrt(7/12)) 
print(np.sqrt(5/24)*np.sqrt(2))

0.7637626158259734
0.6454972243679029

# Chandra (1987) Appendix 3 gives...
# Again OK aside from extra sqrt(2)
print(np.sqrt(7)/(2*np.sqrt(3)))
print(np.sqrt(5)/(2*np.sqrt(6))*np.sqrt(2))

0.7637626158259734
0.645497224367903

Refs:

Altmann, S.L. and Bradley, C.J. (1963) ‘On the Symmetries of Spherical Harmonics’, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 255(1054), pp. 199–215. doi:10.1098/rsta.1963.0002.
Boiko, D.L., Féron, P. and Besnard, P. (2006) ‘Simple method to obtain symmetry harmonics of point groups’, The Journal of Chemical Physics, 125(9), p. 094104. doi:10.1063/1.2338040.
Chandra, N. (1987) ‘Photoelectron spectroscopic studies of polyatomic molecules. I. Theory’, Journal of Physics B: Atomic and Molecular Physics, 20(14), pp. 3405–3415. doi:10.1088/0022-3700/20/14/013.

phockett commented 2 years ago

In case it's of any interest @mcodev31 (or to anyone else), I've now wrapped your solutions for symmetrized harmonic generation along with some basic handling routines with SHtools, Pandas and Xarray. Still a little crude, but hopefully will be useful to some, see here for details.

mcodev31 commented 2 years ago

That looks great! Thanks for the info.

mcodev31 / libmsym

Generating symmetrized harmonics for a given point group #21