diff --git a/qpms/groups.h b/qpms/groups.h index 0ab132e..34f2c44 100644 --- a/qpms/groups.h +++ b/qpms/groups.h @@ -65,6 +65,21 @@ typedef struct qpms_finite_group_t { struct qpms_finite_group_irrep_t *irreps; ///< Irreducible representations of the group. } qpms_finite_group_t; +/// Group multiplication. +static inline qpms_gmi_t qpms_finite_group_mul(const qpms_finite_group_t *G, + qpms_gmi_t a, qpms_gmi_t b) { + assert(a < G->order); + assert(b < G->order); + return G->mt[G->order * a + b]; +} + +/// Group element inversion. +static inline qpms_gmi_t qpms_finite_group_inv(const qpms_finite_group_t *G, + qpms_gmi_t a) { + assert(a < G->order); + return G->invi[a]; +} + /// NOT IMPLEMENTED Get irrep index by name. qpms_iri_t qpms_finite_group_find_irrep_by_name(qpms_finite_group_t *G, char *name); diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index c917cd6..3614437 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -270,13 +270,34 @@ static bool orbit_types_equal(const qpms_ss_orbit_type_t *a, const qpms_ss_orbit return true; } -/// Add orbit type to a scattering system, updating the ss->otspace_end pointer accordingly +// Extend the action to all particles in orbit if only the action on the 0th +// particle has been filled. +static void extend_orbit_action(qpms_scatsys_t *ss, qpms_ss_orbit_type_t *ot) { + for(qpms_ss_orbit_pi_t src = 1; src < ot->size; ++src) { + // find any element g that sends 0 to src: + qpms_gmi_t g; + for (g = 0; g < ss->sym->order; ++g) + if (ot->action[g] == src) break; + assert(g < ss->sym->order); + // invg sends src to 0 + qpms_gmi_t invg = qpms_finite_group_inv(ss->sym, g); + for (qpms_gmi_t f = 0; f < ss->sym->order; ++f) + // if f sends 0 to dest, then f * invg sends src to dest + ot->action[src * ss->sym->order + + qpms_finite_group_mul(ss->sym,f,invg)] = ot->action[f]; + } +} + +//Add orbit type to a scattering system, updating the ss->otspace_end pointer accordingly static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_current) { qpms_ss_orbit_type_t * const ot_new = & (ss->orbit_types[ss->orbit_type_count]); ot_new->size = ot_current->size; - const size_t actionsiz = sizeof(ot_current->action[0]) * ss->sym->order; + const size_t actionsiz = sizeof(ot_current->action[0]) * ot_current->size + * ss->sym->order; ot_new->action = (void *) (ss->otspace_end); memcpy(ot_new->action, ot_current->action, actionsiz); + // N.B. we copied mostly garbage ^^^, most of it is initialized just now: + extend_orbit_action(ss, ot_new); ss->otspace_end += actionsiz; const size_t tmsiz = sizeof(ot_current->tmatrices[0]) * ot_current->size; ot_new->tmatrices = (void *) (ss->otspace_end); @@ -384,12 +405,12 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp ss->orbit_types = calloc(ss->p_count, sizeof(qpms_ss_orbit_type_t)); ss->otspace_end = ss->otspace = malloc( // reallocate later - (sizeof(qpms_ss_orbit_pi_t) * sym->order + (sizeof(qpms_ss_orbit_pi_t) * sym->order * sym->order + sizeof(qpms_ss_tmi_t) * sym->order) * ss->p_count); // Workspace for the orbit type determination qpms_ss_orbit_type_t ot_current; - qpms_ss_orbit_pi_t ot_current_action[sym->order]; + qpms_ss_orbit_pi_t ot_current_action[sym->order * sym->order]; qpms_ss_tmi_t ot_current_tmatrices[sym->order]; qpms_ss_pi_t current_orbit[sym->order]; @@ -499,3 +520,44 @@ void qpms_scatsys_free(qpms_scatsys_t *ss) { free(ss); } + + +// (copypasta from symmetries.c) +// TODO at some point, maybe support also other norms. +// (perhaps just use qpms_normalisation_t_factor() at the right places) +static inline void check_norm_compat(const qpms_vswf_set_spec_t *s) +{ + switch (qpms_normalisation_t_normonly(s->norm)) { + case QPMS_NORMALISATION_POWER: + break; + case QPMS_NORMALISATION_SPHARM: + break; + default: + abort(); // Only SPHARM and POWER norms are supported right now. + } +} + +#if 0 +complex double *qpms_orbit_matrix_action(complex double *target, + const qpms_ss_orbit_type_t *ot, const qpms_vswf_set_spec_t *bspec, + const qpms_finite_group_t *sym, const qpms_gmi_t g) { + assert(sym); assert(g < sym->order); + assert(sym->rep3d); + assert(ot); assert(ot->size > 0); + check_norm_compat(bspec); + const size_t n = bspec->n; + const qpms_gmi_t N = ot->size; + if (target == NULL) + target = malloc(n*n*N*N*sizeof(complex double)); + if (target == NULL) abort(); + memset(target, 0, n*n*N*N*sizeof(complex double)); + complex double tmp[n][n]; // this is the 'per-particle' action + qpms_irot3_uvswfi_dense(tmp[0], bspec, sym->rep3d[g]); + for(qpms_gmi_t Col = 0; Col < ot->size; ++Col) { + // Row is the 'destination' of the symmetry operation, Col is the 'source' + qpms_gmi_t Row = ot->action + +#endif + + + diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 208113e..3547d2c 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -224,13 +224,14 @@ typedef qpms_ss_tmi_t qpms_ss_oti_t; ///< Auxilliary type used for labeling orbi * */ typedef struct qpms_ss_orbit_type_t { - qpms_gmi_t size; ///< Size of the orbit (a divisor of the group order). - /// Action of the group elements onto the "first" element in this orbit. - /** Its size is sym->order and its values lie between 0 and \a this.size − 1. - * - * The corresponding stabilizer {\a g} of the i-th particle on the orbit - * is given by action[i] = g. + qpms_ss_orbit_pi_t size; ///< Size of the orbit (a divisor of the group order). + /// Action of the group elements onto the elements in this orbit. + /** Its size is sym->order * this.size + * and its values lie between 0 and \a this.size − 1. * + * Action of the group element g onto the pi-th particle + * is given by action[g + pi*sym->order]. + * */ qpms_ss_orbit_pi_t *action; /// T-matrix IDs of the particles on this orbit (type). @@ -243,6 +244,25 @@ typedef struct qpms_ss_orbit_type_t { qpms_ss_tmi_t *tmatrices; } qpms_ss_orbit_type_t; +/// Construct a "full matrix action" of a point group element for an orbit type. +/** TODO detailed doc */ +complex double *qpms_orbit_matrix_action( + /// Target array. If NULL, a new one is allocated. + /** The size of the array is (orbit->size * bspec->n)**2 + * (it makes sense to assume all the T-matrices share their spec). + */ + complex double *target, + /// The orbit (type). + const qpms_ss_orbit_type_t *orbit, + /// Base spec of the t-matrices (we don't know it from orbit, as it has + /// only T-matrix indices. + const qpms_vswf_set_spec_t *bspec; + /// The symmetry group used to generate the orbit (must have rep3d filled). + const qpms_finite_group_t *sym, + /// The index of the operation in sym to represent. + const qpms_gmi_t g); + + /// Auxillary type used in qpms_scatsys_t that identifies the particle's orbit and its id inside that orbit. typedef struct qpms_ss_particle_orbitinfo { qpms_ss_oti_t t; ///< Orbit type.