From a5b6f5ce348a8d3f17c5c4fa126afe55e6c461cb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 20 Feb 2019 00:58:49 +0000 Subject: [PATCH] C implementation of the basic point group symmetries in uvswf basis. Former-commit-id: c5559425bf62f741dce8f59a8dc603872ee701f8 --- qpms/indexing.h | 2 + qpms/symmetries.c | 177 ++++++++++++++++++++++++++++++++++++++++++++++ qpms/symmetries.h | 40 +++++++++-- 3 files changed, 215 insertions(+), 4 deletions(-) create mode 100644 qpms/symmetries.c diff --git a/qpms/indexing.h b/qpms/indexing.h index 96ae1de..267eddb 100644 --- a/qpms/indexing.h +++ b/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) { diff --git a/qpms/symmetries.c b/qpms/symmetries.c new file mode 100644 index 0000000..f2dd050 --- /dev/null +++ b/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); +} + diff --git a/qpms/symmetries.h b/qpms/symmetries.h index 1b7207b..8b64713 100644 --- a/qpms/symmetries.h +++ b/qpms/symmetries.h @@ -1,10 +1,10 @@ -#ifndef SYMMETRIES_H -#define SYMMETRIES_H -/* TODO. +/*! \file symmetries.h + * \brief Functions providing point group operations operating on translation + * 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 + +/// 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