Extend the role of qpms_ss_orbit_type_t::action.

Former-commit-id: 8be18f9d502689673aaf4e44806c93fc462e2014

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);
 

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
+
+
+

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).
+	qpms_ss_orbit_pi_t size; ///< Size of the orbit (a divisor of the group order).
-	/// Action of the group elements onto the "first" element in this orbit.
+	/// Action of the group elements onto the elements in this orbit.
-	/** Its size is sym->order and its values lie between 0 and \a this.size − 1.
+	/** Its size is sym->order * this.size
-	 *  
+	 *  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.
 	 *
+	 *  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.