diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 4537269..4b21ebd 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -11,6 +11,7 @@ #include "vectors.h" #include "wigner.h" #include +#include "qpms_error.h" #define QPMS_SCATSYS_LEN_RTOL 1e-13 @@ -86,6 +87,71 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop( return t; } +qpms_errno_t qpms_symmetrise_tmdata_irot3arr( + complex double *tmdata, const size_t tmcount, + const qpms_vswf_set_spec_t *bspec, + const size_t n_symops, const qpms_irot3_t *symops) { + const size_t n = bspec->n; + qpms_tmatrix_t *tmcopy = qpms_tmatrix_init(bspec); + complex double *symop_matrices = malloc(n*n*sizeof(complex double) * n_symops); + if(!symop_matrices) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "malloc() failed."); + for (size_t i = 0; i < n_symops; ++i) + qpms_irot3_uvswfi_dense(symop_matrices + i*n*n, bspec, symops[i]); + complex double tmp[n][n]; + const complex double one = 1, zero = 0; + for (size_t tmi = 0; tmi < tmcount; ++tmi) { + // Move the data in tmcopy; we will then write the sum directly into tmdata. + memcpy(tmcopy->m, tmdata+n*n*tmi, n*n*sizeof(complex double)); + memset(tmdata+n*n*tmi, 0, n*n*sizeof(complex double)); + for (size_t i = 0; i < n_symops; ++i) { + const complex double *const M = symop_matrices + i*n*n; + // tmp = M T + cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, + n, n, n, &one, M, n, tmcopy->m, n, &zero, tmp, n); + // tmdata[...] += tmp M* = M T M* + cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans, + n, n, n, &one, tmp, n, M, n, &one, tmdata + tmi*n*n, n); + } + for (size_t ii = 0; ii < n*n; ++ii) + tmdata[n*n*tmi+ii] /= n_symops; + } + free(symop_matrices); + qpms_tmatrix_free(tmcopy); + return QPMS_SUCCESS; +} + +qpms_errno_t qpms_symmetrise_tmdata_finite_group( + complex double *tmdata, const size_t tmcount, + const qpms_vswf_set_spec_t *bspec, + const qpms_finite_group_t *pointgroup) { + if (!(pointgroup->rep3d)) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "This function requires pointgroup->rep3d to be set correctly!"); + return qpms_symmetrise_tmdata_irot3arr(tmdata, tmcount, bspec, + pointgroup->order, pointgroup->rep3d); +} + +qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace( + qpms_tmatrix_t *T, + size_t n_symops, + const qpms_irot3_t *symops + ) { + if(qpms_symmetrise_tmdata_irot3arr(T->m, 1, + T->spec, n_symops, symops) != QPMS_SUCCESS) + return NULL; + else return T; +} + +qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace( + qpms_tmatrix_t *T, + const qpms_finite_group_t *pointgroup + ) { + if(qpms_symmetrise_tmdata_finite_group(T->m, 1, + T->spec, pointgroup) != QPMS_SUCCESS) + return NULL; + else return T; +} + qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace( qpms_tmatrix_t *T, const complex double *M diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 4d6e540..4805eb4 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -1,5 +1,12 @@ /*! \file scatsystem.h * \brief Modern interface for finite lattice calculations, including symmetries. + * + * N.B. Only "reasonably normalised" waves are supported now in most of the + * functions defined here, i.e. those that can be rotated by the usual + * Wigner matrices, i.e. the "power" or "spharm" -normalised ones. + * + * TODO FIXME check whether Condon-Shortley phase can have some nasty influence + * here; I fear that yes. */ #ifndef QPMS_SCATSYSTEM_H #define QPMS_SCATSYSTEM_H @@ -27,6 +34,9 @@ typedef struct qpms_tmatrix_t { bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free() } qpms_tmatrix_t; +struct qpms_finite_group_t; +typedef struct qpms_finite_group_t qpms_finite_group_t; + /// Returns a pointer to the beginning of the T-matrix row \a rowno. static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){ return t->m + (t->spec->n * rowno); @@ -162,7 +172,6 @@ qpms_errno_t qpms_load_scuff_tmatrix( qpms_tmatrix_t **tmatrices_array, complex double **tmdata ///< The T-matrices raw contents ); - /// Loads a scuff-tmatrix generated file. /** A simple wrapper over qpms_read_scuff_tmatrix() that needs a * path instead of open FILE. @@ -185,20 +194,53 @@ qpms_errno_t qpms_read_scuff_tmatrix( complex double ** tmdata ); -/// Loads scuff-tmatrix generated files. -/** - * freqs, freqs_su, tmatrices_array and tmdata arrays are allocated by this function - * and have to be freed by the caller after use. - * The contents of tmatrices_array is NOT supposed to be freed element per element. +/// In-place application of point group elements on raw T-matrix data. +/** \a tmdata can be e.g. obtained by qpms_load_scuff_tmatrix(). + * The \a symops array should always contain all elements of a finite + * point (sub)group, including the identity operation. + * + * TODO more doc. */ -qpms_errno_t qpms_load_scuff_tmatrix( - const char *path, ///< Path to the TMatrix file - const qpms_vswf_set_spec_t *bspec, ///< VSWF set spec - size_t *n, ///< Number of successfully loaded t-matrices - double **freqs, ///< Frequencies in SI units - double **freqs_su, ///< Frequencies in SCUFF units (optional) - qpms_tmatrix_t **tmatrices_array, ///< The resulting T-matrices. - complex double **tmdata ///< The t-matrices raw contents +qpms_errno_t qpms_symmetrise_tmdata_irot3arr( + complex double *tmdata, const size_t tmcount, + const qpms_vswf_set_spec_t *bspec, + size_t n_symops, + const qpms_irot3_t *symops + ); + +/// In-place application of a point group on raw T-matrix data. +/** This does the same as qpms_symmetrise_tmdata_irot3arr(), + * but takes a valid finite point group as an argument. + * + * TODO more doc. + */ +qpms_errno_t qpms_symmetrise_tmdata_finite_group( + complex double *tmdata, const size_t tmcount, + const qpms_vswf_set_spec_t *bspec, + const qpms_finite_group_t *pointgroup + ); + +/// In-place application of point group elements on a T-matrix. +/** The \a symops array should always contain all elements of a finite + * point (sub)group, including the identity operation. + * + * TODO more doc. + */ +qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace( + qpms_tmatrix_t *T, + size_t n_symops, + const qpms_irot3_t *symops + ); + +/// In-place application of point group elements on a T-matrix. +/** This does the same as qpms_tmatrix_symmetrise_irot3arr(), + * but takes a valid finite point group as an argument. + * + * TODO more doc. + */ +qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace( + qpms_tmatrix_t *T, + const qpms_finite_group_t *pointgroup ); /* Fuck this, include the whole @@ -244,7 +286,6 @@ typedef struct qpms_particle_tid_t { qpms_ss_tmi_t tmatrix_id; ///< T-matrix index } qpms_particle_tid_t; -struct qpms_finite_group_t; typedef qpms_gmi_t qpms_ss_orbit_pi_t; ///< Auxilliary type used in qpms_ss_orbit_type_t for labeling particles inside orbits. typedef qpms_ss_tmi_t qpms_ss_oti_t; ///< Auxilliary type used for labeling orbit types. diff --git a/qpms/tmatrix_io.c b/qpms/tmatrix_io.c index 2ce341d..e9972ae 100644 --- a/qpms/tmatrix_io.c +++ b/qpms/tmatrix_io.c @@ -29,6 +29,10 @@ qpms_errno_t qpms_read_scuff_tmatrix( if (!(freqs && n && tmdata)) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, "freqs, n, and tmdata are mandatory arguments and must not be NULL."); + if (bs->norm != QPMS_NORMALISATION_POWER_CS) // CHECKME CORRECT? + qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)" + " norm supported right now."); int n_alloc = 128; // First chunk to allocate *n = 0; *freqs = malloc(n_alloc * sizeof(double));