diff --git a/qpms/groups.h b/qpms/groups.h index b1198b7..6bfa851 100644 --- a/qpms/groups.h +++ b/qpms/groups.h @@ -19,7 +19,7 @@ typedef const char * qpms_permutation_t; /// To be used only in qpms_finite_group_t struct qpms_finite_group_irrep_t { int dim; ///< Irrep dimension. - char name[]; ///< Irrep label. + char *name; ///< Irrep label. /// Irrep matrix data. /** The r-th row, c-th column of the representation of the i'th element is retrieved as * m[i * dim * dim + r * dim + c] @@ -49,12 +49,12 @@ typedef struct qpms_finite_group_t { qpms_gmi_t *invi; ///< Group elem inverse indices. qpms_gmi_t *gens; ///< A canonical set of group generators. int ngens; ///< Number of the generators in gens; - qpms_permutation_t permrep[]; ///< Permutation representations of the elements. + qpms_permutation_t *permrep; ///< Permutation representations of the elements. char **elemlabels; ///< Optional human readable labels for the group elements. int permrep_nelem; ///< Number of the elements over which the permutation representation acts. - struct qpms_irot3_t rep3d[]; ///< The quaternion representation of a 3D point group (if applicable). + struct qpms_irot3_t *rep3d; ///< The quaternion representation of a 3D point group (if applicable). int nirreps; ///< How many irreps does the group have - struct qpms_finite_group_irrep_t irreps[]; ///< Irreducible representations of the group. + struct qpms_finite_group_irrep_t *irreps; ///< Irreducible representations of the group. } qpms_finite_group_t; diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index f197262..9d20fea 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -185,5 +185,18 @@ cdef extern from "translations.h": char *phi_data, np.npy_intp *phi_shape, np.npy_intp *phi_strides, char *r_ge_d_data, np.npy_intp *phi_shape, np.npy_intp *phi_strides) nogil +cdef extern from "scatsystem.h": + struct qpms_tmatrix_t: + qpms_vswf_set_spec_t *spec + cdouble *m + int owns_m # FIXME in fact bool + struct qpms_particle_t: + cart3_t pos + const qpms_tmatrix_t *tmatrix + struct qpms_tmatrix_interpolator_t: + const qpms_vswf_set_spec_t *bspec + struct qpms_scatsys_t: + pass # TODO + diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 91873c9..872cbe6 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -2,8 +2,19 @@ #include #include "scatsystem.h" #include "indexing.h" +#include "vswf.h" +#include "groups.h" +#include "symmetries.h" #include +#include +#include +#include "vectors.h" +#include "wigner.h" + +#define QPMS_SCATSYS_LEN_RTOL 1e-13 +#define QPMS_SCATSYS_TMATRIX_ATOL 1e-14 +#define QPMS_SCATSYS_TMATRIX_RTOL 1e-12 qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) { qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t)); if (!t) abort(); @@ -30,6 +41,28 @@ void qpms_tmatrix_free(qpms_tmatrix_t *t){ free(t); } +qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace( + qpms_tmatrix_t *T, + const complex double *M + ) +{ + //qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec); + const size_t n = T->spec->n; + complex double tmp[n][n]; + // tmp = M T + const complex double one = 1, zero = 0; + cblas_zgemm(CblasRowMajor, + CblasNoTrans, + CblasNoTrans, + n, n, n, &one, M, n, T->m, n, &zero, tmp, n); + // t->m = tmp M* = M T M* + cblas_zgemm(CblasRowMajor, + CblasNoTrans, + CblasConjTrans, + n, n, n, &one, tmp, n, M, n, &zero, T->m, n); + return T; +} + qpms_tmatrix_t *qpms_tmatrix_apply_symop( const qpms_tmatrix_t *T, const complex double *M @@ -129,7 +162,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(qpms_tmatrix_t *T, int N) { return T; } -bool qpms_tmatrix_isnear(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B, +bool qpms_tmatrix_isclose(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B, const double rtol, const double atol) { if (!qpms_vswf_set_spec_isidentical(A->spec, B->spec)) @@ -218,3 +251,146 @@ qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t return t; } + + +// ------------ Stupid implementation of qpms_scatsys_apply_symmetry() ------------- + +#define MIN(x,y) (((x)<(y))?(x):(y)) +#define MAX(x,y) (((x)>(y))?(x):(y)) + +qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym) { + // TODO check data sanity + + // First, determine the rough radius of the array + double lenscale = 0; + { + double minx = +INFINITY, miny = +INFINITY, minz = +INFINITY; + double maxx = -INFINITY, maxy = -INFINITY, maxz = -INFINITY; + for (ssize_t i = 0; i < orig->p_count; ++i) { + minx = MIN(minx,orig->p[i].pos.x); + miny = MIN(miny,orig->p[i].pos.y); + minz = MIN(minz,orig->p[i].pos.z); + maxx = MAX(maxx,orig->p[i].pos.x); + maxy = MAX(maxy,orig->p[i].pos.y); + maxz = MAX(maxz,orig->p[i].pos.z); + } + lenscale = ((maxx-minx)+(maxy-miny)+(maxz-minz)) / 3; + } + + // Second, check that there are no duplicit positions in the input system. + for (ssize_t i = 0; i < orig->p_count; ++i) + for (ssize_t j = 0; j < i; ++j) + assert(!cart3_isclose(orig->p[i].pos, orig->p[j].pos, 0, QPMS_SCATSYS_LEN_RTOL * lenscale)); + + // Allocate T-matrix and particle arrays + qpms_scatsys_t *ss = malloc(sizeof(qpms_scatsys_t)); + ss->lenscale = lenscale; + + ss->tm_capacity = sym->order * orig->tm_count; + ss->tm = malloc(ss->tm_capacity * sizeof(qpms_tmatrix_t *)); + + ss->p_capacity = sym->order * orig->p_count; + ss->p = malloc(ss->p_capacity * sizeof(qpms_particle_tid_t)); + + // Copy T-matrices; checking for duplicitie + + ptrdiff_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities + ss->tm_count = 0; + for (ptrdiff_t i = 0; i < orig->tm_count; ++i) { + bool found_dupl = false; + ptrdiff_t j; + for (j = 0; j < ss->tm_count; ++j) + if (qpms_tmatrix_isclose(orig->tm[i], ss->tm[j], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL)) { + break; + } + if (j == ss->tm_count) { // duplicity not found, copy the t-matrix + ss->tm[i] = qpms_tmatrix_copy(orig->tm[i]); + ++(ss->tm_count); + } + tm_dupl_remap[i] = j; + } + + // Copy particles, remapping the t-matrix indices + for (ptrdiff_t i = 0; i < orig->p_count; ++(i)) { + ss->p[i] = orig->p[i]; + ss->p[i].tmatrix_id = tm_dupl_remap[ss->p[i].tmatrix_id]; + } + ss->p_count = orig->p_count; + + // allocate t-matrix symmetry map + ss->tm_sym_map = malloc(sizeof(ptrdiff_t) * sym->order * sym->order * ss->tm_count); + + // Extend the T-matrices list by the symmetry operations + for (ptrdiff_t tmi = 0; tmi < ss->tm_count; ++tmi) + for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi){ + const size_t d = ss->tm[tmi]->spec->n; + complex double M[d][d]; // transformation matrix + qpms_irot3_uvswfi_dense(M[0], ss->tm[tmi]->spec, sym->rep3d[gmi]); + qpms_tmatrix_t *transformed = qpms_tmatrix_apply_symop(ss->tm[tmi], M[0]); + ptrdiff_t tmj; + for (tmj = 0; tmj < ss->tm_count; ++tmj) + if (qpms_tmatrix_isclose(transformed, ss->tm[tmj], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL)) + break; + if (tmj < ss->tm_count) { // HIT, transformed T-matrix already exists + qpms_tmatrix_free(transformed); + } else { // MISS, save the matrix and increment the count + ss->tm[ss->tm_count] = transformed; + ++(ss->tm_count); + } + ss->tm_sym_map[gmi + tmi * sym->order] = tmj; // In any case, tmj now indexes the correct transformed matrix + } + // Possibly free some space using the new ss->tm_count instead of (old) ss->tm_count*sym->order + ss->tm_sym_map = realloc(ss->tm_sym_map, sizeof(ptrdiff_t) * sym->order * ss->tm_count); + // tm could be realloc'd as well, but those are just pointers, not likely many. + + // allocate particle symmetry map + ss->p_sym_map = malloc(sizeof(ptrdiff_t) * sym->order * sym->order * ss->p_count); + // Extend the particle list by the symmetry operations, check that particles mapped by symmetry ops on themselves + // have the correct symmetry + // TODO this could be sped up to O(npart * log(npart)); let's see later if it is needed. + for (ptrdiff_t pi = 0; pi < ss->p_count; ++pi) + for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi) { + cart3_t newpoint = qpms_irot3_apply_cart3(sym->rep3d[gmi], ss->p[pi].pos); + ptrdiff_t new_tmi = ss->tm_sym_map[gmi + ss->p[pi].tmatrix_id * sym->order]; // transformed t-matrix index + ptrdiff_t pj; + for (pj = 0; pj < ss->p_count; ++pj) + if (cart3_isclose(newpoint, ss->p[pj].pos, 0, ss->lenscale * QPMS_SCATSYS_LEN_RTOL)) { + if (new_tmi != ss->p[pj].tmatrix_id) + abort(); // The particle is mapped to another (or itself) without having the required symmetry. TODO give some error message. + break; + } + if (pj < ss->p_count) { // HIT, the particle is transformed to an "existing" one. + ; + } else { // MISS, the symmetry transforms the particle to a new location, so add it. + qpms_particle_tid_t newparticle = {newpoint, new_tmi}; + ss->p[ss->p_count] = newparticle; + ++(ss->p_count); + } + ss->p_sym_map[gmi + pi * sym->order] = pi; + } + // Possibly free some space using the new ss->p_count instead of (old) ss->p_count*sym->order + ss->p_sym_map = realloc(ss->p_sym_map, sizeof(ptrdiff_t) * sym->order * ss->p_count); + ss->p = realloc(ss->p, sizeof(qpms_particle_tid_t) * ss->p_count); ss->p_capacity = ss->p_count; + + ss->sym = sym; + return ss; +} + + + + + + + + + + + + + + + + + + + diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 34c32c3..651176b 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -52,7 +52,7 @@ typedef struct qpms_particle_t { * This function actually checks for identical vswf specs. * TODO define constants with "default" atol, rtol for this function. */ -bool qpms_tmatrix_isnear(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2, +bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2, const double rtol, const double atol); /// Creates a T-matrix from another matrix and a symmetry operation. @@ -66,6 +66,17 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop( const complex double *M ///< the symmetry op matrix in row-major format ); +/// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data. +/** The symmetry operation is expected to be a unitary (square) + * matrix \a M with the same + * VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then + * correspond to CHECKME \f[ T' = MTM^\dagger \f] + */ +qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace( + qpms_tmatrix_t *T, ///< the original T-matrix + const complex double *M ///< the symmetry op matrix in row-major format + ); + /// Symmetrizes a T-matrix with an involution symmetry operation. /** The symmetry operation is expected to be a unitary (square) * matrix \a M with the same @@ -154,7 +165,7 @@ extern const gsl_interp_type * gsl_interp_akima_periodic; extern const gsl_interp_type * gsl_interp_steffen; */ -// struct gsl_interp_accel; // use if lookup shows to be too slow +// struct gsl_interp_accel; // use if lookup proves to be too slow typedef struct qpms_tmatrix_interpolator_t { const qpms_vswf_set_spec_t *bspec; //bool owns_bspec; @@ -177,8 +188,42 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num //, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix. /*, ...? */); + +/// A particle, defined by its T-matrix INDEX and position, to be used in qpms_scatsys_t. +typedef struct qpms_particle_tid_t { + // Does it make sense to ever use other than cartesian coords for this? + cart3_t pos; ///< Particle position in cartesian coordinates. + ptrdiff_t tmatrix_id; ///< T-matrix index +} qpms_particle_tid_t; + +struct qpms_finite_group_t; + typedef struct qpms_scatsys_t { - ; + // TODO does bspec belong here? + qpms_tmatrix_t **tm; ///< T-matrices in the system + ptrdiff_t tm_count; ///< Number of all different T-matrices + ptrdiff_t tm_capacity; ///< Capacity of tm[]. + qpms_particle_tid_t *p; ///< Particles. + ptrdiff_t p_count; ///< Size of particles array. + ptrdiff_t p_capacity; ///< Capacity of p[]. + + //TODO the index types do not need to be so big. + struct qpms_finite_group_t *sym; ///< Symmetry group of the array + ptrdiff_t *p_sym_map; ///< Which particles map onto which by the symmetry ops. + ///< p_sym_map[idi + pi * sym->order] gives the index of pi-th particle under the idi'th sym op. + ptrdiff_t *tm_sym_map; ///< Which t-matrices map onto which by the symmetry ops. Lookup by tm_sum_map[idi + tmi * sym->order]. + + // TODO shifted origin of the symmetry group etc. + // TODO some indices for fast operations here. + // private + double lenscale; // radius of the array, used as a relative tolerance measure } qpms_scatsys_t; +/// Creates a new scatsys by applying a symmetry group, copying particles if needed. +/** In fact, it copies everything except the vswf set specs, so keep them alive until scatsys is destroyed. + */ +qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym); +/// Destroys the result of qpms_scatsys_apply_symmetry. +void qpms_scatsys_free(qpms_scatsys_t *s); + #endif //QPMS_SCATSYSTEM_H diff --git a/qpms/vectors.h b/qpms/vectors.h index ece28fb..64bf664 100644 --- a/qpms/vectors.h +++ b/qpms/vectors.h @@ -122,6 +122,14 @@ static inline cart3_t cart3_scale(const double c, const cart3_t v) { return res; } +static inline double cart3_dist(const cart3_t a, const cart3_t b) { + return cart3norm(cart3_substract(a,b)); +} + +static inline bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) { + return cart3_dist(a,b) <= atol + rtol * (cart3norm(b) + cart3norm(a)) * .5; +} + static inline ccart3_t ccart3_scale(const complex double c, const ccart3_t v) { ccart3_t res = {c * v.x, c * v.y, c * v.z}; return res; @@ -161,6 +169,7 @@ static inline double cart3_dot(const cart3_t a, const cart3_t b) { return a.x * b.x + a.y * b.y + a.z * b.z; } + // equivalent to sph_loccart2cart in qpms_p.py static inline ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at) { const double st = sin(at.theta);