C implementation of the basic point group symmetries in uvswf basis.

Former-commit-id: c5559425bf62f741dce8f59a8dc603872ee701f8

qpms/indexing.h

@@ -51,6 +51,8 @@ static inline qpms_y_t qpms_lMax2nelem_sc(qpms_l_t lmax){
 	return lmax * ((qpms_y_t)lmax + 2) + 1;
 }
 
+// TODO maybe enable crashing / validity control by macro definitions...
+
 /// Conversion from VSWF type, order and degree to universal index.
 static inline qpms_uvswfi_t qpms_tmn2uvswfi(
 		qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n) {

qpms/symmetries.c

@@ -0,0 +1,177 @@
+#include "symmetries.h"
+#include "tiny_inlines.h"
+#include "indexing.h"
+
+// 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.
+  }
+}
+
+complex double *qpms_zflip_uvswi_dense(
+    complex double *target,
+    const qpms_vswf_set_spec_t *bspec) 
+{
+  check_norm_compat(bspec);
+  const size_t n = bspec->n;
+  if (target == NULL)
+    target = malloc(n * n * sizeof(complex double));
+  if (target == NULL) abort();
+
+  for (size_t row = 0; row < n; row++) {
+    qpms_vswf_type_t rt;
+    qpms_l_t rl;
+    qpms_m_t rm;
+    qpms_uvswfi2tmn(bspec->ilist[row], &rt, &rm, &rl);
+    for (size_t col = 0; col < n; col++) {
+      qpms_vswf_type_t ct;
+      qpms_l_t cl;
+      qpms_m_t cm;
+      if(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl)) abort();
+      if (rl == cl && rm == cm && rt == ct)
+        switch(rt) {
+          case QPMS_VSWF_ELECTRIC:
+          case QPMS_VSWF_LONGITUDINAL:
+            target[n*row + col] = min1pow(cm + cl);
+            break;
+          case QPMS_VSWF_MAGNETIC:
+            target[n*row + col] = -min1pow(cm + cl);
+            break;
+          default:
+            abort();
+        }
+      else target[n*row + col] = 0;
+    }
+  }
+  return target;
+}
+
+complex double *qpms_yflip_uvswi_dense(
+    complex double *target,
+    const qpms_vswf_set_spec_t *bspec) 
+{
+  check_norm_compat(bspec);
+  const size_t n = bspec->n;
+  if (target == NULL)
+    target = malloc(n * n * sizeof(complex double));
+  if (target == NULL) abort();
+
+  for (size_t row = 0; row < n; row++) {
+    qpms_vswf_type_t rt;
+    qpms_l_t rl;
+    qpms_m_t rm;
+    qpms_uvswfi2tmn(bspec->ilist[row], &rt, &rm, &rl);
+    for (size_t col = 0; col < n; col++) {
+      qpms_vswf_type_t ct;
+      qpms_l_t cl;
+      qpms_m_t cm;
+      if(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl)) abort();
+      if (rl == cl && rm == -cm && rt == ct)
+        switch(rt) {
+          case QPMS_VSWF_ELECTRIC:
+          case QPMS_VSWF_LONGITUDINAL:
+            target[n*row + col] = min1pow(rm);
+            break;
+          case QPMS_VSWF_MAGNETIC:
+            target[n*row + col] = -min1pow(rm);
+            break;
+          default:
+            abort();
+        }
+      else target[n*row + col] = 0;
+    }
+  }
+  return target;
+}
+
+complex double *qpms_xflip_uvswi_dense(
+    complex double *target,
+    const qpms_vswf_set_spec_t *bspec) 
+{
+  check_norm_compat(bspec);
+  const size_t n = bspec->n;
+  if (target == NULL)
+    target = malloc(n * n * sizeof(complex double));
+  if (target == NULL) abort();
+
+  for (size_t row = 0; row < n; row++) {
+    qpms_vswf_type_t rt;
+    qpms_l_t rl;
+    qpms_m_t rm;
+    qpms_uvswfi2tmn(bspec->ilist[row], &rt, &rm, &rl);
+    for (size_t col = 0; col < n; col++) {
+      qpms_vswf_type_t ct;
+      qpms_l_t cl;
+      qpms_m_t cm;
+      if(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl)) abort();
+      if (rl == cl && rm == -cm && rt == ct)
+        switch(rt) {
+          case QPMS_VSWF_ELECTRIC:
+          case QPMS_VSWF_LONGITUDINAL:
+            target[n*row + col] = 1;
+            break;
+          case QPMS_VSWF_MAGNETIC:
+            target[n*row + col] = -1;
+            break;
+          default:
+            abort();
+        }
+      else target[n*row + col] = 0;
+    }
+  }
+  return target;
+}
+
+// Dense matrix representation of a rotation around the z-axis
+complex double *qpms_zrot_uvswi_dense(
+                complex double *target, ///< If NULL, a new array is allocated.
+                const qpms_vswf_set_spec_t *bspec,
+                double phi ///< Rotation angle
+                )
+{
+  check_norm_compat(bspec);
+  const size_t n = bspec->n;
+  if (target == NULL)
+    target = malloc(n * n * sizeof(complex double));
+  if (target == NULL) abort();
+
+  for (size_t row = 0; row < n; row++) {
+    qpms_vswf_type_t rt;
+    qpms_l_t rl;
+    qpms_m_t rm;
+    qpms_uvswfi2tmn(bspec->ilist[row], &rt, &rm, &rl);
+    for (size_t col = 0; col < n; col++) {
+      qpms_vswf_type_t ct;
+      qpms_l_t cl;
+      qpms_m_t cm;
+      if(qpms_uvswfi2tmn(bspec->ilist[col], &ct, &cm, &cl)) abort();
+      if (rl == cl && rm == cm && rt == ct) // TODO COMPARE WITH PYTHON
+        target[n*row + col] = cexp(/* - ?*/I * rm * phi);
+      else target[n*row + col] = 0;
+    }
+  }
+  return target;
+}
+
+// Dense matrix representation of a "rational" rotation around the z-axis
+/* Just for convenience. Corresponds to the angle \f$ \phi = 2\piw/N \f$.
+ */
+complex double *qpms_zrot_rational_uvswi_dense(
+                complex double *target, ///< If NULL, a new array is allocated.
+                const qpms_vswf_set_spec_t *bspec,
+                int N,
+                int w
+) 
+{
+  double phi = 2 * M_PI * w / N;
+  return qpms_zrot_uvswi_dense(target, bspec, phi);
+}
+

qpms/symmetries.h

@@ -1,10 +1,10 @@
-#ifndef SYMMETRIES_H
+/*! \file symmetries.h
-#define SYMMETRIES_H
+ * \brief Functions providing point group operations operating on translation
-/* TODO.
+ * operators and T-matrices.
  *
  * Here will be functions providing point group operations
  * operating on translation operators and T-matrices
- * as in tmatrices.py.
+ * similar to tmatrices.py.
  * 
  * At least I want:
  * - Wigner D matrices
@@ -19,5 +19,37 @@
  *
  * 
  */
+#ifndef SYMMETRIES_H
+#define SYMMETRIES_H
+#include "vswf.h"
+#include <cblas.h>
+
+/// Dense matrix representation of the z coordinate sign flip operation (xy-plane mirroring).
+complex double *qpms_zflip_uvswi_dense(
+		complex double *target, ///< If NULL, a new array is allocated.
+		const qpms_vswf_set_spec_t *bspec);
+/// Dense matrix representation of the y coordinate sign flip operation (xz-plane mirroring).
+complex double *qpms_yflip_uvswi_dense(
+		complex double *target, ///< If NULL, a new array is allocated.
+		const qpms_vswf_set_spec_t *bspec);
+/// Dense matrix representation of the x coordinate sign flip operation (yz-plane mirroring).
+complex double *qpms_xflip_uvswi_dense(
+		complex double *target, ///< If NULL, a new array is allocated.
+		const qpms_vswf_set_spec_t *bspec);
+/// Dense matrix representation of a rotation around the z-axis
+complex double *qpms_zrot_uvswi_dense(
+		complex double *target, ///< If NULL, a new array is allocated.
+		const qpms_vswf_set_spec_t *bspec,
+		double phi ///< Rotation angle
+		);
+/// Dense matrix representation of a "rational" rotation around the z-axis
+/** Just for convenience. Corresponds to the angle \f$ \phi = 2\piw/N \f$.
+ */
+complex double *qpms_zrot_rational_uvswi_dense(
+		complex double *target, ///< If NULL, a new array is allocated.
+		const qpms_vswf_set_spec_t *bspec,
+		int N,
+		int w
+		);
 
 #endif // SYMMETRIES_H