mcodev31 / libmsym

molecular point group symmetry lib
MIT License
75 stars 32 forks source link

Error: "Point group has no primary axis for reorientation" #17

Closed ChangChunHe closed 3 years ago

ChangChunHe commented 4 years ago

C59N.txt I want to obtain the point group of the C59N structures I have uploaded, but an error came when I used ./msym_example.

error msg is:

Molecule has center of mass [-0.009490; 0.002012; 0.001649] and a radius of 3.562240
Error Error determining point group: Point group has no primary axis for reorientation

The version I used is the current mater branch, any reply will be appreciated.

mcodev31 commented 4 years ago

Hi Chang,

It's typically not possible to detect all symmetries with a single set of thresholds like in the example code, and some structures where the eigenvalues of the inertial tensor are very similar the thresholds may need changing. You can do that in code by supplying a msym_thresholds_t structure (can also be supplied to the example).

If it's not extremely important to use those exact numbers, you can change the N to a C, use msym_example and choose the Ih symmertry, then at the end of the example choose Y when asked if you want to symmetrize the coordinates. Then use the symmetrized coordinates and replace the C back to an N again.

I believe you would get Cs symmetry since there should be no proper rotations through an atom in C60 (don't quote me on that, it's been a while since I looked at it).

AlexB67 commented 3 years ago

Hello,

Thanks for this great contribution. I have a similar issue. I already allow for changing the thresholds in my code. 'till now, this has fixed that issue for cases I come across occasionally when that happens, that is, until today.

34 dichloro butene, to no avail will it succeed. unit angstrom C 0.210965397593 0.243244574037 0.407277464365 C 1.395665397593 -0.495455425963 -0.148022535635 C -0.957034602407 0.313544574037 -0.574522535635 Cl 0.735065397593 1.914144574037 0.803777464365 C 1.865765397593 -1.636855425963 0.368477464365 Cl -1.699034602407 -1.291255425963 -0.814522535635 H -0.120234602407 -0.192555425963 1.356377464365 H 1.889465397593 -0.083655425963 -1.025522535635 H -1.750334602407 0.964544574037 -0.193122535635 H -0.649434602407 0.687844574037 -1.556922535635 H 2.720865397593 -2.131955425963 -0.078122535635 H 1.406865397593 -2.097355425963 1.236577464365

I've tried altering the thresholds significantly by a factor of 10 (what is reasonable ?) libmsym always comes back with, Point group has no primary axis for reorientation"

psi4 symmetry detector works fine with this molecule without any tricks, just to double check. I am using the library via C++ btw, so no calls via python. structure looks okay in avogrado.

Any suggestions before I start delving into the inner guts of libmsym ?

Cheers.

ghutchis commented 3 years ago

@AlexB67 - what symmetry are you expecting?

I popped this into Avogadro2 and I don't see anything obvious. What does Psi4 assign it?

Screen Shot 2021-02-07 at 7 14 54 PM
AlexB67 commented 3 years ago

Hello, thank for the reply. C1, it's what i would expect .. correct ? C1 is what psi4 says it is. but in libmsym it ends up empty NULL (because it fails earlier on), I can in a hack sort of way work around it and assign C1 if it fails as a fall back, for now. For my purposes that's okay. Naturally I rather do it properly :)

ghutchis commented 3 years ago

Yes, it's C1. I guess your concern is the message? It's not so much of an error as an indication that the molecule has no symmetry.

AlexB67 commented 3 years ago

It is an error alright in that libmsym is failing somewhere. I had a quick look at the source earlier at lunch, but since I am not familiar with the code base, i would need to go through it with a debugger properly one day. By inspection not so easy :).

In any case, libmsym should return C1, not NULL for the point group, that is only because it is failing somewhere (I suspect) because no primary axis could be found.

Indeed, the msym error struct returns an error set somewhere in or around alignAxes (of top of my head, not at the sourcce code right now)

Cheers

mcodev31 commented 3 years ago

Hi Alex,

C1 resulting in an error is mostly out of laziness (it's intentional) since the missing primary axis results in problems with symmetrization of orbitals down the line. It's not so much that I want it to return an error as I wanted to make the code in some quantum calculcation software that uses this easier, and I took a shortcut, rather than having special handing of C1 in several places. Obviously it would be be better to implement those cases, but I haven't had any time to spend on this project in a long time, and no one has ever wanted to use the C1 group.

This patch will return C1, but keep in mind that some things may break when using that group

diff --git a/src/point_group.c b/src/point_group.c
index e0cfa13..ceca630 100644
--- a/src/point_group.c
+++ b/src/point_group.c
@@ -757,6 +757,7 @@ err:

 msym_error_t transformAxes(msym_point_group_type_t type, int n, msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t *sops, msym_thresholds_t *thresholds, double transform[3][3]){
     msym_error_t ret = MSYM_SUCCESS;
+   if(type == MSYM_POINT_GROUP_TYPE_Cn && n == 1) return ret;
     switch (type){
         case (MSYM_POINT_GROUP_TYPE_Ci)  :
         case (MSYM_POINT_GROUP_TYPE_K)   :
AlexB67 commented 3 years ago

Hello, thank you for the reply. I'll leave it as then, as in my own code I can work with it.

I only reported it as a suspected or unintended bug after looking at the code again, but you have now confirmed my thoughts. Funnily enough your thought of the proposed change is what I suspected I needed to do along those lines. Right now I do it after libmsym has returned with said error, so it boils down to the same thing for my needs.

Case closed as far as I am concerned. Many thanks,

Cheers