Determination of the symmetry of normal modes seems uncertain

[physicien]

Hi. First of all, awesome code! I have long searched for a code capable of determining the symmetries of the normal modes.

I'm trying to use the code to assign the symmetries of the normal modes of vibration of the fullerene C60. The normal modes were calculated by DFT with ORCA.

Here is the hessian file containing the coordinates, symbols, normal modes and frequencies. FRQ_c60_f_PBE.hess.txt

Here is the parser used by my script to extract the data hess_parser.py.txt

And here is my python script.

#!/usr/bin/python3
from hess_parser import HessianData
from posym import SymmetryModes

# LOCAL VARIABLES
HESSPATH = './FRQ_c60_f_PBE.hess'
hessian = HessianData(HESSPATH)
gr='Ih'

coordinates = hessian.atoms
symbols = hessian.symbol
normal_modes = hessian.normal_modes[6:]
frequencies = hessian.freq[6:]

sym_modes_gs = SymmetryModes(group=gr, coordinates=coordinates, \
        modes=normal_modes, symbols=symbols)
for i in range(len(normal_modes)):
    print('Mode {:2}: {:8.3f} :'.format(i + 7, frequencies[i]), \
            sym_modes_gs.get_state_mode(i))

print('Total symmetry: ', sym_modes_gs)

The problem is that when assigning, there seems to be uncertainty for some modes. See for example modes 12, 13 and 14 in this extract of the result.

Mode  7:  257.019 : 0.2 Hg
Mode  8:  257.064 : 0.2 Hg
Mode  9:  257.076 : 0.2 Hg
Mode 10:  257.510 : 0.2 Hg
Mode 11:  257.532 : 0.2 Hg
Mode 12:  333.111 :  - 0.12 Ag - 0.06 T1g - 0.06 T2g + 0.12 Gg - 0.07 Au + 0.13 T1u + 0.13 T2u + 0.07 Gu
Mode 13:  333.165 :  - 0.12 Ag - 0.06 T1g - 0.06 T2g + 0.12 Gg - 0.07 Au + 0.13 T1u + 0.13 T2u + 0.07 Gu
Mode 14:  333.173 :  - 0.12 Ag - 0.06 T1g - 0.06 T2g + 0.12 Gg - 0.07 Au + 0.13 T1u + 0.13 T2u + 0.07 Gu
Mode 15:  343.869 : 0.25 Gu
Mode 16:  345.617 : 0.25 Gu
Mode 17:  345.702 : 0.25 Gu
Mode 18:  345.738 : 0.25 Gu
...

I tried to solve the problem by myself, but not finding the source of it I'm reporting it

Best, Emmanuel

[abelcarreras]

Hi Emmanuel,

I suspect it has something to do with the orientation of the molecule. This group has been never tested properly so I am happy you checked it and reported the issue. I will work on this issue next week.

[physicien]

Hi Abel,

I'm glad to help with testing of some point groups. I should report that I also tested the assignment of symmetries on molecules of groups D5d, D5h, D3d, D3h and D3. The molecules had between 60 and 150 atoms.

I had exactly the same problem with all these tests. In some cases, uncertain assignments strongly predominated.

If you need more tests, just contact me and I'll see what I can do.

Best, Emmanuel

[abelcarreras]

Hi Emmanuel,

I fixed I/Ih groups, and probable some others. There was some bugs in Rotation class and IR table. If you want to try I uploaded it to development branch. I tested D5h,D3, and D3h groups and seem to work now. Now I'm working on fixing the IR table generation for Dnd groups. If you do more tests and find other groups that fail would be great. I want to use this opportunity to fix them all.

Thanks for your help!

[abelcarreras]

I think I fixed most of the issues with the different groups. Let me know if you find more bugs.

[physicien]

Hi Abel,

Thanks for your quick work. I just installed version 0.23. Your fixes seem to work, but I believe I just found another underlying bug in my testing. I will therefore detail my tests for different groups and the errors observed:

I_h point group

I tested two data sets. The first one is the result of an analytical frequency calculation of C60 molecule in ORCA. The assignment is done perfectly. Here is the hessian file FRQ_c60_f_PBE.hess Here is the assignment c60-PBE_mode_sym.out

The second data set is the result of a Raman spectrum calculation for the same molecule once again in ORCA. Here is the hessian file RAM_c60-Ih_r.hess Here is the assignment c60-Ih_mode_sym.out

Observed errors are often in the form of a very small added term, such as Mode 34-36, or the mixing of the Gg and Gu from Mode 135 to Mode 142.

...
Mode 34:  519.952 : 0.31 T1u + 0.01 Hu
Mode 35:  519.998 : 0.31 T1u + 0.01 Hu
Mode 36:  520.018 : 0.31 T1u + 0.01 Hu
...
Mode 135: 1297.307 : 0.19 Gg + 0.06 Gu
Mode 136: 1297.735 : 0.17 Gg + 0.08 Gu
Mode 137: 1297.827 : 0.06 Gg + 0.19 Gu
Mode 138: 1297.953 : 0.08 Gg + 0.17 Gu
Mode 139: 1298.034 : 0.25 Gu
Mode 140: 1298.076 : 0.25 Gg
Mode 141: 1298.620 : 0.25 Gg
Mode 142: 1300.107 : 0.25 Gu
...

D_5d point group

I applied the same approach as for the previous point group, but here both the calculations gave assignments with a small term added. Here is the hessian file FRQ_c100_f_PBE.hess Here is the assignment c100-PBE_mode_sym.out

D_5h point group

Same as the previous point group.

D_3d point group

Same as the previous point group.

Hope this helps isolate the problem. Emmanuel

[abelcarreras]

About Ih point group I think the results may be correct. In the Raman spectrum calculation the mix symmetry values are not just small but the total symmetry returns an integer sum of IR. This may suggest that the molecule present some small distortion that mixes the symmetry of some particular modes.

D5d group seems very bad, I will check it.

Are D3d and D5h point groups also failing? This is strange, I tested these groups with Hexamethylbenzene and Pentaprismane and work perfectly.

[abelcarreras]

I fixed (hopefully) the issues with groups Dnd and Dnh. Now they should work.

[physicien]

Awesome! I'll redo the tests for D5d, D5h and the D3's, and I'll get back to you with the results.

[physicien]

Hi Abel,

I've been testing your fixes for four days and I still have strange results. I tested Ih, D5d, D5h, D3d and D3h point groups.

From now on, I'll consider that an assignment is good if

• the mix symmetry values are small
• the total symmetry returns an integer sum of IR

For the little story, I started by testing the point groups Ih and D5d. Regardless of the settings that I'll change later, Ih works. The case of D5d is more interesting. The assignment seems to work, except when the normal modes have been calculated with the PBE functional. In order to rule out any problem related to the numerical precision of my data, I systematically attempted the assignment for each group on normal modes calculated (FREQ) from the functionals PBE, r2SCAN and B3LYP (with the def2-TZVP basis set). I also tried to assign the symmetries of the normal modes taken from the Raman calculations (NUMFREQ) for some groups. Here is a summary of those results:

Gr  Method  Functional  Assignment
#############################################
C60-Ih
    FREQ
        PBE     GOOD
        r2SCAN      GOOD
        B3LYP       GOOD
    NUMFREQ
        PBE     GOOD

C70-D5h
    FREQ
        PBE     BAD
        r2SCAN      BAD
        B3LYP       BAD
    NUMFREQ
        PBE     BAD

C78-D3h
    FREQ
        PBE     BAD
        r2SCAN      BAD
        B3LYP       BAD

C96-D3d
    FREQ
        PBE     BAD
        r2SCAN      BAD
        B3LYP       BAD
    NUMFREQ
        PBE     BAD

C100-D5d
    FREQ
        PBE     BAD
        r2SCAN      GOOD
        B3LYP       GOOD
    NUMFREQ
        PBE     BAD
        r2SCAN      GOOD

About the bad assignments in D5d, as they only happen with PBE (so the least accurate calculations here) I wonder if it's not a numerical error problem. Correct me if I'm wrong, but the algorithm determines that there is a symmetry if it can permute a point with an operated point, and is this permutation determined by the distance between these two points (I suppose it must tend to 0?).

If so, what is the maximum allowable error on the distance to detect a permutation?

As for the other point groups, I could not say if the problem is the same, or if there is a problem specific to these groups. They all have non-integer coefficients in the total symmetry. I can share the assignments and hessian files with you if you need them.

Best, Emmanuel

---

EDIT

I was apparently using version 0.23. I redid all the assignments with 0.24 and

• the results are now correct for Ih and D5d
• the results are still bad for the other point groups

[abelcarreras]

Hi Emmanuel,

The way PoSym works is using the idea of continuous symmetry. This way there is no tolerance parameters. The summary is as follows: 1) The orientation of the molecule to the symmetry point group: This is done using only the atomic coordinates by applying each symmetry operation to the geometry and comparing the distances between the original geometry and operated geometry. This requires to find the best permutation (of atoms) that minimizes these distances. Each symmetry operation requires a different permutation. If the molecule has the symmetry of the group there should be a orientation for which each symmetry operation has a permutation that makes this distances zero. The average of these distances is given by the method measure_pos

https://github.com/abelcarreras/posym/blob/79468832c539987f7188c2450d38f4b3f8cdb36a/posym/__init__.py#L184

2) Once determined the orientation, for each symmetry operation the code determines the object overlap. This overlap is determined as the scalar product between the object original object and the object with the operation applied. In the case of Normal modes object, this is done by calculating the dot product of each mode eigenvector with itself with the operation (and permutation) applied. If perfect symmetry, for non degenerated modes, these overlaps will be 1 or -1 for all symmetry operations.

3) The overlaps of all symmetry operations can be seen as a vector in some kind of symmetry vector space expressed in the basis of symmetry operations. Using the symmetry tables of the corresponding group as a transformation matrix (change of basis) this vector is converted to the irreducible representations basis. If perfect symmetry, in this basis all the components of this vector will be zero, except for one. This non-zero one indicates the symmetry of the object (Normal mode in this case). If the object does not have the symmetry of the group you will find that several components of this vector are not zero, these are the mix symmetry cases. The sum of all the components of this vector (multiplied by its degeneracy) is always 1 by definition (basically because of the normalization of the original object). You can check this vector by accessing get_op_representation and get_ir_representation methods.

https://github.com/abelcarreras/posym/blob/79468832c539987f7188c2450d38f4b3f8cdb36a/posym/__init__.py#L56

https://github.com/abelcarreras/posym/blob/79468832c539987f7188c2450d38f4b3f8cdb36a/posym/__init__.py#L59

This scheme and machinery is the exactly the same used to determine the symmetry of all objects PoSym can analyze (orbitals/wavefunctions/density/etc..), the only difference is the way in which the overlaps are calculated.

[abelcarreras]

Regarding the groups that are not working properly yet, it would be very helpful if you can share the hessian files with me. I only tested these groups with a few small molecules and it seemed to work. I need non-working examples for debugging!

https://github.com/abelcarreras/posym/tree/master/test/molecules

Again, thanks a lot for your testing!

[physicien]

Thank you for the detailed explanation, it clarifies a lot.

Here is a selection of the point groups that I tested. I even found an example of D3 and it works. https://github.com/physicien/storage

I should also have examples of D2d and D2h point groups in the near future. I'll let you know of the result.

By the way, do you plan to add a "How to cite" section to your README.md (with a BibTeX entry)? It will eventually be useful for me to know how to reference PoSym.

[abelcarreras]

Hi Emmanuel,

I found that the issue is related to the point group orientation algorithm. If you define the orientation manually (using orientation_angles option) it seems to work. For example for C70:

sym_modes_gs = SymmetryModes(group=gr, coordinates=coordinates,
                             modes=normal_modes, symbols=symbols,
                             orientation_angles = [0.00000000e+00,  90, -1.13969995e-07])

These angles are the Euler angles that transform the original coordinates as this:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.transform.Rotation.from_euler.html

The tough part in determining the orientation is the fact that the function the I wrote to optimize the orientation has many local minima (roughly more atoms, more minima) and is not continuous everywhere (because of permutations). Therefore I do a preliminary systematic scan along the sphere to find a good starting point and then optimize using conjugate gradient assuming that the function will be continuous at least around the minima. This preliminary scan is slow (very slow for large molecules) but large molecules have more local minima and require finer scan grids to not get stuck in a local minima. These days I have tried some strategies to do this scan in a more clever way but none of them works in a reliable way for all examples I tested. I think the only solution will be to translate part of the optimization code to C for speed but it may take some time.

[abelcarreras]

I just updated PoSym to fix the orientation issue. Now all your examples work. Also the speed of the calculation is largely improved. About how to cite I will upload the code to zenodo.org and include the citation in the package.

[physicien]

Wow, I just tested v0.25 and it's definitely faster. Congratulations, all the problems I had are fixed. I'll close the post and open another one if I come across any other problematic groups.

Thanks again and great job.