diff --git a/TODOs.rst b/TODOs.rst index c4d6a27..d85be8a 100644 --- a/TODOs.rst +++ b/TODOs.rst @@ -1,5 +1,13 @@ -TODO label description -====================== +URGENT +====== + +- Check exact normalisation convention of scuff-tmatrix output. +- Check whether the Condon-Shortley phase affects the form of Wigner matrices. +- The xflip, yflip and possible i-factor problem. + + +TODO label description (obsolete!) +================================== LMAXVAR - besides to int, the function should support an iterable for the lMax argument or similar, and in such case the output should have the according additional dimension. In the end, I want diff --git a/qpms/error.c b/qpms/error.c index 6027190..305cb11 100644 --- a/qpms/error.c +++ b/qpms/error.c @@ -43,6 +43,20 @@ void qpms_pr_error_at_flf(const char *filename, unsigned int linenum, fflush(stderr); abort(); } + +void qpms_pr_debug_at_flf(const char *filename, unsigned int linenum, + const char *func, + const char *fmt, ...) { + fprintf(stderr, "%s:%u, %s: ", filename, linenum, func); + va_list ap; + va_start(ap, fmt); + vfprintf(stderr, fmt, ap); + va_end(ap); + fputc('\n', stderr); + fflush(stderr); +} + + //void qpms_error_at_line(const char *filename, unsigned int linenum, // const char *fmt, ...); diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h index 45f88aa..70aa7f9 100644 --- a/qpms/qpms_error.h +++ b/qpms/qpms_error.h @@ -13,6 +13,9 @@ void qpms_pr_error_at_flf(const char *filename, unsigned int linenum, const char *func, const char *fmt, ...); +void qpms_pr_debug_at_flf(const char *filename, unsigned int linenum, + const char *func, + const char *fmt, ...); //void qpms_error_at_line(const char *filename, unsigned int linenum, // const char *fmt, ...); diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 4b21ebd..900d5c8 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1,5 +1,6 @@ #include #include +#include #include "scatsystem.h" #include "indexing.h" #include "vswf.h" @@ -13,7 +14,7 @@ #include #include "qpms_error.h" - +#define SQ(x) ((x)*(x)) #define QPMS_SCATSYS_LEN_RTOL 1e-13 #define QPMS_SCATSYS_TMATRIX_ATOL 1e-14 #define QPMS_SCATSYS_TMATRIX_RTOL 1e-12 @@ -358,6 +359,10 @@ static void extend_orbit_action(qpms_scatsys_t *ss, qpms_ss_orbit_type_t *ot) { static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_current) { qpms_ss_orbit_type_t * const ot_new = & (ss->orbit_types[ss->orbit_type_count]); ot_new->size = ot_current->size; + + const qpms_vswf_set_spec_t *bspec = ss->tm[ot_current->tmatrices[0]]->spec; + const size_t bspecn = bspec->n; + const size_t actionsiz = sizeof(ot_current->action[0]) * ot_current->size * ss->sym->order; ot_new->action = (void *) (ss->otspace_end); @@ -365,10 +370,36 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu // N.B. we copied mostly garbage ^^^, most of it is initialized just now: extend_orbit_action(ss, ot_new); ss->otspace_end += actionsiz; + const size_t tmsiz = sizeof(ot_current->tmatrices[0]) * ot_current->size; ot_new->tmatrices = (void *) (ss->otspace_end); memcpy(ot_new->tmatrices, ot_current->tmatrices, tmsiz); ss->otspace_end += tmsiz; + + const size_t irbase_sizes_siz = sizeof(ot_new->irbase_sizes[0]) * ss->sym->nirreps; + ot_new->irbase_sizes = (void *) (ss->otspace_end); + ss->otspace_end += irbase_sizes_siz; + ot_new->irbase_cumsizes = (void *) (ss->otspace_end); + ss->otspace_end += irbase_sizes_siz; + + const size_t irbases_siz = sizeof(ot_new->irbases[0]) * SQ(ot_new->size * bspecn); + ot_new->irbases = (void *) (ss->otspace_end); + ss->otspace_end += irbases_siz; + + size_t lastbs, bs_cumsum = 0; + for(qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) { + qpms_orbit_irrep_basis(&lastbs, + ot_new->irbases + bs_cumsum*ot_new->size*bspecn, + ot_new, bspec, ss->sym, iri); + ot_new->irbase_sizes[iri] = lastbs; + bs_cumsum += lastbs; + ot_new->irbase_cumsizes[iri] = bs_cumsum; + } + if(bs_cumsum != ot_new->size * bspecn) + qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "The cumulative size of the symmetry-adapted bases is wrong; " + "expected %d = %d * %d, got %d.", + ot_new->size * bspecn, ot_new->size, bspecn, bs_cumsum); } @@ -414,6 +445,8 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp } // Copy T-matrices; checking for duplicities + + size_t max_bspecn = 0; // We'll need it later.for memory alloc estimates. qpms_ss_tmi_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities ss->tm_count = 0; @@ -428,6 +461,7 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp ++(ss->tm_count); } tm_dupl_remap[i] = j; + max_bspecn = MAX(ss->tm[i]->spec->n, max_bspecn); } // Copy particles, remapping the t-matrix indices @@ -472,7 +506,10 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp ss->otspace_end = ss->otspace = malloc( // reallocate later (sizeof(qpms_ss_orbit_pi_t) * sym->order * sym->order - + sizeof(qpms_ss_tmi_t) * sym->order) * ss->p_count); + + sizeof(qpms_ss_tmi_t) * sym->order + + 2*sizeof(size_t) * sym->nirreps + + sizeof(complex double) * SQ(sym->order * max_bspecn)) * ss->p_count + ); // Workspace for the orbit type determination qpms_ss_orbit_type_t ot_current; @@ -502,7 +539,11 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp 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. + qpms_pr_error("The %d. particle with coords (%lg, %lg, %lg) " + "is mapped to %d. another (or itself) with cords (%lg, %lg, %lg) " + "without having the required symmetry", (int)pi, + ss->p[pi].pos.x, ss->p[pi].pos.y, ss->p[pi].pos.z, + (int)pj, ss->p[pj].pos.x, ss->p[pj].pos.y, ss->p[pj].pos.z); break; } if (pj < ss->p_count) { // HIT, the particle is transformed to an "existing" one. @@ -557,6 +598,8 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp for (size_t oi = 0; oi < ss->orbit_type_count; ++oi) { ss->orbit_types[oi].action = (void *)(((char *) (ss->orbit_types[oi].action)) + shift); ss->orbit_types[oi].tmatrices = (void *)(((char *) (ss->orbit_types[oi].tmatrices)) + shift); + ss->orbit_types[oi].irbase_sizes = (void *)(((char *) (ss->orbit_types[oi].irbase_sizes)) + shift); + ss->orbit_types[oi].irbases = (void *)(((char *) (ss->orbit_types[oi].irbases)) + shift); } } @@ -660,4 +703,65 @@ complex double *qpms_orbit_irrep_projector_matrix(complex double *target, return target; } +#define SVD_ATOL 1e-8 + +complex double *qpms_orbit_irrep_basis(size_t *basis_size, + complex double *target, + const qpms_ss_orbit_type_t *ot, const qpms_vswf_set_spec_t *bspec, + const qpms_finite_group_t *sym, const qpms_iri_t iri) { + assert(sym); + assert(sym->rep3d); + assert(ot); assert(ot->size > 0); + assert(iri < sym->nirreps); assert(sym->irreps); + check_norm_compat(bspec); + const size_t n = bspec->n; + const qpms_gmi_t N = ot->size; + if (target == NULL) + target = malloc(n*n*N*N*sizeof(complex double)); + if (target == NULL) abort(); + memset(target, 0, n*n*N*N*sizeof(complex double)); + + // Get the projector (also workspace for right sg. vect.) + complex double *projector = qpms_orbit_irrep_projector_matrix(NULL, + ot, bspec, sym, iri); + if(!projector) abort(); + // Workspace for the right singular vectors. + complex double *V_H = malloc(n*n*N*N*sizeof(complex double)); + if(!V_H) abort(); + // THIS SHOULD NOT BE NECESSARY + complex double *U = malloc(n*n*N*N*sizeof(complex double)); + if(!U) abort(); + double *s = malloc(n*N*sizeof(double)); if(!s) abort(); + + int info = LAPACKE_zgesdd(LAPACK_ROW_MAJOR, + 'A', // jobz; overwrite projector with left sg.vec. and write right into V_H + n*N /* m */, n*N /* n */, projector /* a */, n*N /* lda */, + s /* s */, U /* u */, n*N /* ldu, irrelev. */, V_H /* vt */, + n*N /* ldvt */); + if (info) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "Something went wrong with the SVD."); + + + size_t bs; + for(bs = 0; bs < n*N; ++bs) { + qpms_pr_debug_at_flf(__FILE__, __LINE__, __func__, + "%d. irrep, %zd. SV: %.16lf", (int) iri, bs, s[bs]); + if(s[bs] > 1 + SVD_ATOL) qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "%zd. SV too large: %.16lf", bs, s[bs]); + if(s[bs] > SVD_ATOL && fabs(1-s[bs]) > SVD_ATOL) + qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, + "%zd. SV in the 'wrong' interval: %.16lf", bs, s[bs]); + if(s[bs] < SVD_ATOL) break; + } + + memcpy(target, V_H, bs*N*n*n*sizeof(complex double)); + if(basis_size) *basis_size = bs; + + free(U); + free(V_H); + free(projector); + return target; +} + + diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 3ce6b37..188e071 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -332,13 +332,34 @@ typedef struct qpms_ss_orbit_type_t { * The size of this array is \a size. */ qpms_ss_tmi_t *tmatrices; + /// Sizes of the per-orbit irrep bases. + /** + * The order of the irreps corresponds to the order in \a ss->sym->irreps. + * The size of this array is (obviously) \a ss->sym->nirreps. + * + * TODO different type? + * TODO doc. + */ + size_t *irbase_sizes; + /// Cumulative sums of irbase_sizes. + size_t *irbase_cumsizes; + /// Per-orbit irreducible representation orthonormal bases. + /** This also defines the unitary operator that transforms the orbital excitation coefficients + * in the symmetry-adapted basis. + * + * The size is (\a this->size * \a this->tmatrices[0].spec->n)**2. + * + * TODO doc. + */ + complex double *irbases; + } qpms_ss_orbit_type_t; /// Auxillary type used in qpms_scatsys_t that identifies the particle's orbit and its id inside that orbit. typedef struct qpms_ss_particle_orbitinfo { qpms_ss_oti_t t; ///< Orbit type. #define QPMS_SS_P_ORBITINFO_UNDEF (-1) ///< This labels that the particle has not yet been assigned to an orbit. - qpms_ss_orbit_pi_t p; ///< "Order" of the particle inside that orbit type. + qpms_ss_orbit_pi_t p; ///< Order (sija, ei rankki) of the particle inside that orbit type. } qpms_ss_particle_orbitinfo_t; @@ -435,6 +456,26 @@ complex double *qpms_orbit_irrep_projector_matrix( /// The index of the irreducible representation of sym. const qpms_iri_t iri); +/// TODO DOC!!!!! +complex double *qpms_orbit_irrep_basis( + /// Here theh size of theh basis shall be saved, + size_t *basis_size, + /// Target array. If NULL, a new one is allocated. + /** The size of the array is basis_size * (orbit->size * bspec->n) + * (it makes sense to assume all the T-matrices share their spec). + */ + complex double *target, + /// The orbit (type). + const qpms_ss_orbit_type_t *orbit, + /// Base spec of the t-matrices (we don't know it from orbit, as it has + /// only T-matrix indices. + const qpms_vswf_set_spec_t *bspec, + /// The symmetry group used to generate the orbit (must have rep3d filled). + const struct qpms_finite_group_t *sym, + /// The index of the irreducible representation of sym. + const qpms_iri_t iri); + + #if 0 // Abstract types that describe T-matrix/particle/scatsystem symmetries // To be implemented later. See also the thoughts in the beginning of groups.h. diff --git a/setup.py b/setup.py index 068187d..9993d9d 100755 --- a/setup.py +++ b/setup.py @@ -35,7 +35,7 @@ qpms_c = Extension('qpms_c', '-DDISABLE_NDEBUG', # uncomment to enable assertions in the modules #'-fopenmp', ], - libraries=['gsl', 'blas', 'gslcblas', #'omp' + libraries=['gsl', 'lapacke', 'blas', 'gslcblas', #'omp' # TODO resolve the problem with openblas (missing gotoblas symbol) and preferable use other blas library ], runtime_library_dirs=os.environ['LD_LIBRARY_PATH'].split(':') if 'LD_LIBRARY_PATH' in os.environ else []