From 213853e4078d00396e95679c5397ac017cff9716 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Sat, 29 Feb 2020 01:35:12 +0200
Subject: [PATCH] WIP cython support for periodic lattices

Former-commit-id: f679f63941ec8dd597a5ccd140682de889e21807
---
 qpms/qpms_c.pyx     | 23 ++++++++++++++++++++++-
 qpms/qpms_cdefs.pxd | 19 +++++++++++++++++++
 qpms/scatsystem.h   |  2 +-
 3 files changed, 42 insertions(+), 2 deletions(-)

diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx
index f99095f..1143a15 100644
--- a/qpms/qpms_c.pyx
+++ b/qpms/qpms_c.pyx
@@ -244,6 +244,16 @@ cdef class FinitePointGroup:
             self.G = <qpms_finite_group_t *>0
             self.owns_data = False
 
+    property order: # LPTODO might instead be __len__() if iterable at some point
+        def __get__(self):
+            return self.G[0].order
+
+    @staticmethod 
+    def TRIVIAL(): # quite ugly
+        from .symmetries import point_group_info
+        return FinitePointGroup(point_group_info['trivial_g'])
+
+
 cdef class FinitePointGroupElement:
     '''TODO'''
     cdef readonly FinitePointGroup G
@@ -360,7 +370,12 @@ cdef class ScatteringSystem:
         #TODO is there a way to disable the constructor outside this module?
 
     @staticmethod # We don't have any "standard" constructor for this right now
-    def create(particles, medium, FinitePointGroup sym, cdouble omega): # TODO tolerances
+    def create(particles, medium, cdouble omega, FinitePointGroup sym = FinitePointGroup.TRIVIAL(),
+            latticebasis = None): # TODO tolerances
+
+        if latticebasis is not None and sym.order != 1:
+            raise NotImplementedError("Periodic systems don't currently support nontrivial point group symmetries")
+
         # These we are going to construct
         cdef ScatteringSystem self
         cdef _ScatteringSystemAtOmega pyssw
@@ -419,6 +434,12 @@ cdef class ScatteringSystem:
                 tm_derived_key = (tmg_key, None) # TODO unique representation of p.p.op instead of None
                 orig.p[pi].pos = p.cval().pos
                 orig.p[pi].tmatrix_id = tmindices[tm_derived_key]
+            if latticebasis is not None: # periodic system
+                assert(len(latticebasis) <= 3 and len(latticebasis) > 0)
+                orig.lattice_dimension = len(latticebasis)
+                for d in range(len(latticebasis)):
+                    orig.per.lattice_basis[d] = {'x' : latticebasis[d][0], 'y' : latticebasis[d][1], 'z' : latticebasis[d][2]}
+            else: orig.lattice_dimension = 0
             ssw = qpms_scatsys_apply_symmetry(&orig, sym.rawpointer(), omega, &QPMS_TOLERANCE_DEFAULT)
             ss = ssw[0].ss
         finally:
diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd
index a492820..1f39ebe 100644
--- a/qpms/qpms_cdefs.pxd
+++ b/qpms/qpms_cdefs.pxd
@@ -529,7 +529,11 @@ cdef extern from "scatsystem.h":
     struct qpms_ss_derived_tmatrix_t:
         qpms_ss_tmgi_t tmgi
         qpms_tmatrix_operation_t op
+    struct qpms_scatsys_periodic_info_t:
+        cart3_t lattice_basis[3]
+        #etc.
     struct qpms_scatsys_t:
+        int lattice_dimension
         qpms_epsmu_generator_t medium
         qpms_tmatrix_function_t *tmg
         qpms_ss_tmgi_t tmg_count
@@ -541,6 +545,9 @@ cdef extern from "scatsystem.h":
         size_t fecv_size
         size_t *saecv_sizes
         const qpms_finite_group_t *sym
+        qpms_scatsys_periodic_info_t per
+
+        # per[] and other stuff not currently needed in cython
     void qpms_scatsys_free(qpms_scatsys_t *s)
     qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path) #NI
     qpms_scatsys_t *qpms_scatsys_load(char *path) #NI
@@ -583,6 +590,7 @@ cdef extern from "scatsystem.h":
             const qpms_scatsys_at_omega_t *ssw)
     struct qpms_ss_LU:
         const qpms_scatsys_at_omega_t *ssw
+        const qpms_scatsys_at_omega_k_t *sswk
         bint full
         qpms_iri_t iri
         cdouble *a
@@ -602,6 +610,17 @@ cdef extern from "scatsystem.h":
     beyn_result_t *qpms_scatsys_finite_find_eigenmodes(const qpms_scatsys_t *ss, qpms_iri_t iri,
             cdouble omega_centre, double omega_rr, double omega_ri, size_t contour_npoints, 
             double rank_tol, size_t rank_sel_min, double res_tol)
+    # periodic-related funs
+    struct qpms_scatsys_at_omega_k_t:
+        const qpms_scatsys_at_omega_t *ssw
+        double k[3]
+    cdouble *qpms_scatsyswk_build_modeproblem_motrix_full(cdouble *target, const qpms_scatsys_at_omega_k_t *sswk)
+    cdouble *qpms_scatsys_periodic_build_translation_matrix_full(cdouble *target, const qpms_scatsys_t *ss, cdouble wavenumber, const cart3_t *wavevector)
+    cdouble *qpms_scatsyswk_build_translation_matrix_full(cdouble *target, const qpms_scatsys_at_omega_k_t *sswk)
+    qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(cdouble *target, int *target_piv, const qpms_scatsys_at_omega_t *sswk)
+    beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(const qpms_scatsys_t *ss, const double *k,
+            cdouble omega_centre, double omega_rr, double omega_ri, size_t contour_npoints, 
+            double rank_tol, size_t rank_sel_min, double res_tol)
 
 cdef extern from "ewald.h":
     struct qpms_csf_result:
diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h
index d13bac5..a83b1ad 100644
--- a/qpms/scatsystem.h
+++ b/qpms/scatsystem.h
@@ -158,7 +158,7 @@ struct qpms_epsmu_generator_t;
 /// Common "class" for system of scatterers, both periodic and non-periodic.
 /**
  * Infinite periodic structures (those with \a lattice_dimension > 0) 
- * have the \a per element allocated and filled.
+ * have the \a per filled.
  * These are ignored for finite systems (lattice_dimension == 0).
  */
 typedef struct qpms_scatsys_t {