Example update and fixes for new SALC generation

More extensive example code Fix new SALC generation problems with primary axis New API got get eq set by element Change tex threshold Make setup.py require python 3
mcodev31 · Nov 8, 2015 · 4adbebf · 4adbebf
1 parent c150835
commit 4adbebf
Show file tree

Hide file tree

Showing 8 changed files with 399 additions and 169 deletions.
diff --git a/bindings/python/setup.py b/bindings/python/setup.py
@@ -9,10 +9,15 @@
 #
 
 from distutils.core import setup
+import sys
 
+if (not sys.version_info[0] is 3):
+    sys.exit('libmsym module requires python 3')
+
 setup(name='libmsym',
       version='0.2',
       description = 'libmsym python binding',
+      license='MIT',
       author='Marcus Johansson',
       author_email='mcodev31@gmail.com',
       url='https://github.com/mcodev31/libmsym',

diff --git a/examples/example.c b/examples/example.c
diff --git a/examples/tex_character_table.c b/examples/tex_character_table.c
@@ -121,7 +121,7 @@ void characterTableToTex(FILE *fp, msym_point_group_type_t type, int n, const ch
 #define M_PI 3.14159265358979323846264338327950288419716939937510582
 #endif
 
-#define CHARACTER_EQUAL 1.0e-9
+#define CHARACTER_EQUAL 1.0e-7
 
 void characterToTex(FILE *fp,int n, const msym_character_table_t *ct, int i, int j, int mode){
     double (*table)[ct->d] = (double (*)[ct->d]) ct->table;

diff --git a/src/basis_function.c b/src/basis_function.c
@@ -51,7 +51,7 @@ msym_error_t basisFunctionFromQuantumNumbers(int n, int l, int m, msym_basis_fun
         }
         default : {
             char t = l <= 20 ? ('f' - 3 + l + (l >= 7) + (l >= 12) + (l >= 14)) : '?';
-            char *d = (m < 0 ? "-" : "+");
+            char *d = (m == 0 ? "" : m < 0 ? "-" : "+");
             snprintf(bf->name, sizeof(bf->name), "%d%c%d%s",n,t,abs(m),d);
         }
     }
@@ -80,7 +80,7 @@ msym_error_t basisFunctionFromName(char *name, msym_basis_function_t *bf){
             break;
         }
         default :
-            if(cl < 'd' || cl == 'e' || cl == 'j' ||cl > 'z') goto err;
+            if(cl < 'd' || cl == 'e' || cl == 'j' || cl > 'z') goto err;
 
             l = cl - 'd' + 2 - (cl > 'e') - (cl > 'j') - (cl > 'p') - (cl > 's');
 

diff --git a/src/context.c b/src/context.c
@@ -265,6 +265,18 @@ msym_error_t msymGetSubgroups(msym_context ctx, int *sgl, const msym_subgroup_t
         int sgmax = numberOfSubgroups(ctx->pg);
         if(MSYM_SUCCESS != (ret = findPermutationSubgroups(ctx->pg->order, ctx->pg->perm, sgmax, ctx->pg->sops, &gsgl, &gsg))) goto err;
 
+        if(isLinearSubgroup(ctx->pg)){
+            gsg = realloc(gsg, (gsgl+1)*sizeof(*gsg));
+            memset(&gsg[gsgl], 0, sizeof(*gsg));
+            gsg[gsgl].n = ctx->pg->n;
+            gsg[gsgl].order = ctx->pg->order;
+            gsg[gsgl].sops = calloc(gsg[gsgl].order, sizeof(*gsg[gsgl].sops));
+            for(int i = 0;i < ctx->pg->order;i++){
+                gsg[gsgl].sops[i] = &ctx->pg->sops[i];
+            }
+            gsgl++;
+        }
+
         ctx->sg = gsg;
         ctx->sgl = gsgl;
 
@@ -663,6 +675,30 @@ msym_error_t ctxGetEquivalenceSets(msym_context ctx, int *esl, msym_equivalence_
     return ret;
 }
 
+msym_error_t msymGetEquivalenceSetByElement(msym_context ctx, msym_element_t *element, const msym_equivalence_set_t **es){
+    msym_error_t ret = MSYM_SUCCESS;
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; goto err;}
+    if(ctx->es == NULL) {ret = MSYM_INVALID_EQUIVALENCE_SET;goto err;}
+
+    if(element >= ctx->ext.set_elements_ptr && element < ctx->ext.set_elements_ptr + ctx->elementsl){
+        element = element - ctx->ext.set_elements_ptr + ctx->ext.elements;
+    } else if(!(element >= ctx->ext.elements && element < ctx->ext.elements + ctx->elementsl)){
+        ret = MSYM_INVALID_ELEMENTS;
+        msymSetErrorDetails("Element not within [%p,%p) or [%p,%p) but is at %p",ctx->ext.set_elements_ptr,ctx->ext.set_elements_ptr + ctx->elementsl,ctx->ext.elements,ctx->ext.elements + ctx->elementsl, element);
+        goto err;
+    }
+
+    msym_equivalence_set_t **eesmap = NULL;
+
+    if(MSYM_SUCCESS != (ret = ctxGetExternalElementEquivalenceSetMap(ctx, &eesmap))) goto err;
+
+    *es = eesmap[element - ctx->ext.elements];
+
+err:
+    return ret;
+
+}
+
 msym_error_t ctxGetExternalEquivalenceSets(msym_context ctx, int *esl, msym_equivalence_set_t **es){
     msym_error_t ret = MSYM_SUCCESS;
 

diff --git a/src/msym.h b/src/msym.h
@@ -135,7 +135,7 @@ extern "C" {
     } msym_partner_function_t;
 
     typedef struct _msym_salc {
-        int d;              // dimension of space (same as msym_character_table_t.s[msym_subspace_t.s].d)
+        int d;              // dimension of space (same as msym_character_table_t.s[msym_subrepresentation_space_t.s].d)
         int fl;             // number of basis functions
         void *pf;           // partner functions double[d][fl]
         msym_basis_function_t **f;
@@ -180,6 +180,7 @@ extern "C" {
     msym_error_t MSYM_EXPORT msymSelectSubgroup(msym_context ctx, const msym_subgroup_t *subgroup);
     msym_error_t MSYM_EXPORT msymGetSymmetryOperations(msym_context ctx, int *sopsl, const msym_symmetry_operation_t **sops);
     msym_error_t MSYM_EXPORT msymGetEquivalenceSets(msym_context ctx, int *l, const msym_equivalence_set_t **es);
+    msym_error_t MSYM_EXPORT msymGetEquivalenceSetByElement(msym_context ctx, msym_element_t *element, const msym_equivalence_set_t **es);
     msym_error_t MSYM_EXPORT msymGetSubrepresentationSpaces(msym_context ctx, int *l, const msym_subrepresentation_space_t **srs);
     msym_error_t MSYM_EXPORT msymGetCharacterTable(msym_context ctx, const msym_character_table_t **ct);
 

diff --git a/src/point_group.c b/src/point_group.c
@@ -897,12 +897,13 @@ msym_error_t reduceLinearPointGroup(msym_point_group_t *pg, int n, msym_threshol
 
 }
 
+/* Really should split the orientation and class from the symmetry operation
+ * and move it to the group so we don't have to regenerate */
 msym_error_t pointGroupFromSubgroup(const msym_subgroup_t *sg, msym_thresholds_t *thresholds, msym_point_group_t **opg){
     msym_error_t ret = MSYM_SUCCESS;
     *opg = calloc(1,sizeof(msym_point_group_t));
     msym_point_group_t *pg = *opg;
     pg->type = sg->type;
-    pg->primary = sg->primary;
     pg->n = sg->n;
     pg->sops = malloc(sizeof(msym_symmetry_operation_t[sg->order]));
     memcpy(pg->name,sg->name,sizeof(pg->name));
@@ -941,6 +942,15 @@ msym_error_t pointGroupFromSubgroup(const msym_subgroup_t *sg, msym_thresholds_t
     minv(pg->transform, T);
 
     for(int i = 0; i < pg->order;i++){
+        // vector and orientation may not be equal
+        if(NULL != sg->primary &&
+           NULL == pg->primary &&
+           pg->sops[i].type == sg->primary->type &&
+           pg->sops[i].order == sg->primary->order &&
+           pg->sops[i].power == sg->primary->power){
+
+            pg->primary = &pg->sops[i];
+        }
         mvmul(pg->sops[i].v,T,pg->sops[i].v);
     }
 

diff --git a/src/subspace.c b/src/subspace.c
@@ -666,11 +666,13 @@ msym_error_t generateSplittingOperation(msym_point_group_t *pg, msym_permutation
         case MSYM_POINT_GROUP_TYPE_Dn  :
         case MSYM_POINT_GROUP_TYPE_Dnh :
         case MSYM_POINT_GROUP_TYPE_Dnd : {
-            sop = pg->primary;
-            if(NULL == sop){
-                ret = MSYM_SUBSPACE_ERROR;
-                msymSetErrorDetails("Cannot determine splitting operation for point group %s", pg->name);
-                goto err;
+            if(pg->n > 2){
+                sop = pg->primary;
+                if(NULL == sop){
+                    ret = MSYM_SUBSPACE_ERROR;
+                    msymSetErrorDetails("Cannot determine splitting operation for point group %s", pg->name);
+                    goto err;
+                }
             }
 
             break;
@@ -895,7 +897,6 @@ msym_error_t symmetrySpeciesComponents(msym_point_group_t *pg, int srsl, msym_su
 
 
 msym_error_t generateSubrepresentationSpaces(msym_point_group_t *pg, int sgl, const msym_subgroup_t sg[sgl], int esl, msym_equivalence_set_t *es, msym_permutation_t **perm, int basisl, msym_basis_function_t basis[basisl], msym_element_t *elements, msym_equivalence_set_t **esmap, msym_thresholds_t *thresholds, int *osrsl, msym_subrepresentation_space_t **osrs, msym_basis_function_t ***osrsbf, int **ospan){
-    printf("NEW\n");
     msym_error_t ret = MSYM_SUCCESS;
     msym_character_table_t *ct = pg->ct;
     int lmax = -1, nmax = 0, eslmax = 0;
@@ -1341,7 +1342,7 @@ msym_error_t generateSubrepresentationSpaces(msym_point_group_t *pg, int sgl, co
 
 
 
-msym_error_t generateSubrepresentationSpacesOld(msym_point_group_t *pg, int sgl, const msym_subgroup_t sg[sgl], int esl, msym_equivalence_set_t *es, msym_permutation_t **perm, int basisl, msym_basis_function_t basis[basisl], msym_element_t *elements, msym_equivalence_set_t **esmap, msym_thresholds_t *thresholds, int *osrsl, msym_subrepresentation_space_t **osrs, msym_basis_function_t ***osrsbf, int **ospan){
+msym_error_t generateSubrepresentationSpacesLowMem(msym_point_group_t *pg, int sgl, const msym_subgroup_t sg[sgl], int esl, msym_equivalence_set_t *es, msym_permutation_t **perm, int basisl, msym_basis_function_t basis[basisl], msym_element_t *elements, msym_equivalence_set_t **esmap, msym_thresholds_t *thresholds, int *osrsl, msym_subrepresentation_space_t **osrs, msym_basis_function_t ***osrsbf, int **ospan){
     msym_error_t ret = MSYM_SUCCESS;
     msym_character_table_t *ct = pg->ct;
     int lmax = -1, nmax = 0;
@@ -1747,7 +1748,6 @@ msym_error_t generateSubrepresentationSpacesOld(msym_point_group_t *pg, int sgl,
 
                             memset(found,0,sizeof(int[ct->s[sk].d][mdim]));
                             for(int s = 0;s < pg->order ;s++){
-                                printf("loop\n");
                                 permutationMatrix(&perm[i][s], mperm);
                                 kron(d, mperm, ld, lst[s], dd, dscal);
                                 for(int dim = 0;dim < mdim;dim++){