diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index bffd0a4..eae1737 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -700,7 +700,7 @@ cdef class BaseSpec: if 'norm' in kwargs.keys(): self.s.norm = kwargs['norm'] else: - self.s.norm = QPMS_NORMALISATION_POWER + self.s.norm = QPMS_NORMALISATION_POWER_CS # set the other metadata cdef qpms_l_t l self.s.lMax_L = -1 diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h index 33597ab..798b5d7 100644 --- a/qpms/qpms_error.h +++ b/qpms/qpms_error.h @@ -19,9 +19,9 @@ void qpms_pr_debug_at_flf(const char *filename, unsigned int linenum, //void qpms_error_at_line(const char *filename, unsigned int linenum, // const char *fmt, ...); -#define QPMS_CRASHING_MALLOC(pointer, size) {(pointer) = malloc(size); if(!pointer) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Allocation of %zd bytes for " #pointer " failed.", (size_t) (size));} +#define QPMS_CRASHING_MALLOC(pointer, size) {(pointer) = malloc(size); if(!pointer && (size)) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Allocation of %zd bytes for " #pointer " failed.", (size_t) (size));} -#define QPMS_CRASHING_REALLOC(pointer, size) {(pointer) = realloc(pointer, size); if(!pointer) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Rellocation of %zd bytes for " #pointer " failed.", (size_t) (size));} +#define QPMS_CRASHING_REALLOC(pointer, size) {(pointer) = realloc(pointer, size); if(!pointer && (size)) qpms_pr_debug_at_flf(__FILE__,__LINE__,__func__, "Rellocation of %zd bytes for " #pointer " failed.", (size_t) (size));} #define QPMS_WTF qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error.") diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 7a66196..4bf02ef 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -561,6 +561,10 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp current_orbit[0] = pi; ot_current.size = 1; ot_current.tmatrices[0] = ss->p[pi].tmatrix_id; +#ifdef DUMP_PARTICLE_POSITIONS + cart3_t pos = ss->p[pi].pos; + fprintf(stderr, "An orbit [%.4g, %.4g, %.4g] => ", pos.x, pos.y, pos.z); +#endif } for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi) { @@ -583,6 +587,10 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp qpms_particle_tid_t newparticle = {newpoint, new_tmi}; ss->p[ss->p_count] = newparticle; ++(ss->p_count); +#ifdef DUMP_PARTICLE_POSITIONS + if(new_orbit) + fprintf(stderr, "[%.4g, %.4g, %.4g] ", newpoint.x, newpoint.y, newpoint.z); +#endif } ss->p_sym_map[gmi + pi * sym->order] = pj; @@ -600,6 +608,9 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qp } } if (new_orbit) { // Now compare if the new orbit corresponds to some of the existing types. +#ifdef DUMP_PARTICLE_POSITIONS + fputc('\n', stderr); +#endif qpms_ss_oti_t oti; for(oti = 0; oti < ss->orbit_type_count; ++oti) if (orbit_types_equal(&ot_current, &(ss->orbit_types[oti]))) break; // HIT, orbit type already exists @@ -719,7 +730,7 @@ complex double *qpms_orbit_action_matrix(complex double *target, const qpms_gmi_t Row = ot->action[sym->order * Col + g]; for(size_t row = 0; row < bspec->n; ++row) for(size_t col = 0; col < bspec->n; ++col) - target[n*n*N*Row + n*Col + n*N*row + col] = tmp[row][col]; //CHECKCONJ + target[n*n*N*Row + n*Col + n*N*row + col] = conj(tmp[row][col]); //CHECKCONJ } #ifdef DUMP_ACTIONMATRIX fprintf(stderr,"%d: %s\n", @@ -855,6 +866,119 @@ complex double *qpms_orbit_irrep_basis(size_t *basis_size, return target; } +complex double *qpms_scatsys_irrep_transform_matrix(complex double *U, + const qpms_scatsys_t *ss, qpms_iri_t iri) { + const size_t packedlen = ss->saecv_sizes[iri]; + const size_t full_len = ss->fecv_size; + if (U == NULL) + QPMS_CRASHING_MALLOC(U,full_len * packedlen * sizeof(complex double)); + memset(U, 0, full_len * packedlen * sizeof(complex double)); + + size_t fullvec_offset = 0; + + for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { + const qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t; + const qpms_ss_orbit_type_t *const ot = ss->orbit_types + oti; + const qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn; + const qpms_ss_orbit_pi_t opi = ss->p_orbitinfo[pi].p; + // This is where the particle's orbit starts in the "packed" vector: + const size_t packed_orbit_offset = ss->saecv_ot_offsets[iri*ss->orbit_type_count + oti] + + osn * ot->irbase_sizes[iri]; + // Orbit coeff vector's full size: + const size_t orbit_fullsize = ot->size * ot->bspecn; + const size_t particle_fullsize = ot->bspecn; + const size_t orbit_packedsize = ot->irbase_sizes[iri]; + // This is the orbit-level matrix projecting the whole orbit onto the irrep. + const complex double *om = ot->irbases + ot->irbase_offsets[iri]; + + for (size_t prow = 0; prow < orbit_packedsize; ++prow) + for (size_t pcol = 0; pcol < ot->bspecn; ++pcol) + U[full_len * (packed_orbit_offset + prow) + (fullvec_offset + pcol)] + = om[orbit_fullsize * prow + (opi * ot->bspecn + pcol)]; + fullvec_offset += ot->bspecn; + } + + return U; +} + +complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed, + const complex double *orig_full, const qpms_scatsys_t *ss, + qpms_iri_t iri){ + const size_t packedlen = ss->saecv_sizes[iri]; + if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? + return target_packed; + const size_t full_len = ss->fecv_size; + if (target_packed == NULL) + target_packed = malloc(SQ(packedlen)*sizeof(complex double)); + if (target_packed == NULL) abort(); + memset(target_packed, 0, SQ(packedlen)*sizeof(complex double)); + + // Workspace for the intermediate matrix + complex double *tmp; + QPMS_CRASHING_MALLOC(tmp, full_len * packedlen * sizeof(complex double)); + + complex double *U = qpms_scatsys_irrep_transform_matrix(NULL, ss, iri); + + const complex double one = 1, zero = 0; + + // tmp = F U* + cblas_zgemm( + CblasRowMajor, CblasNoTrans, CblasConjTrans, + full_len /*M*/, packedlen /*N*/, full_len /*K*/, + &one /*alpha*/, orig_full/*A*/, full_len/*ldA*/, + U /*B*/, full_len/*ldB*/, + &zero /*beta*/, tmp /*C*/, packedlen /*LDC*/); + // target = U tmp + cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, + packedlen /*M*/, packedlen /*N*/, full_len /*K*/, + &one /*alpha*/, U/*A*/, full_len/*ldA*/, + tmp /*B*/, packedlen /*ldB*/, &zero /*beta*/, + target_packed /*C*/, packedlen /*ldC*/); + + free(tmp); + free(U); + return target_packed; +} + +/// Transforms a big "packed" matrix into the full basis (trivial implementation). +complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full, + const complex double *orig_packed, const qpms_scatsys_t *ss, + qpms_iri_t iri, bool add){ + const size_t packedlen = ss->saecv_sizes[iri]; + const size_t full_len = ss->fecv_size; + if (target_full == NULL) + target_full = malloc(SQ(full_len)*sizeof(complex double)); + if (target_full == NULL) abort(); + if(!add) memset(target_full, 0, SQ(full_len)*sizeof(complex double)); + + if(!packedlen) return target_full; // Empty irrep, do nothing. + + // Workspace for the intermediate matrix result + complex double *tmp; + QPMS_CRASHING_MALLOC(tmp, full_len * packedlen * sizeof(complex double)); + + complex double *U = qpms_scatsys_irrep_transform_matrix(NULL, ss, iri); + + const complex double one = 1, zero = 0; + + // tmp = P U + cblas_zgemm( + CblasRowMajor, CblasNoTrans, CblasNoTrans, + packedlen /*M*/, full_len /*N*/, packedlen /*K*/, + &one /*alpha*/, orig_packed/*A*/, packedlen/*ldA*/, + U /*B*/, full_len/*ldB*/, + &zero /*beta*/, tmp /*C*/, full_len /*LDC*/); + // target += U* tmp + cblas_zgemm(CblasRowMajor, CblasConjTrans, CblasNoTrans, + full_len /*M*/, full_len /*N*/, packedlen /*K*/, + &one /*alpha*/, U/*A*/, full_len/*ldA*/, + tmp /*B*/, full_len /*ldB*/, &one /*beta*/, + target_full /*C*/, full_len /*ldC*/); + free(tmp); + free(U); + return target_full; +} + complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed, const complex double *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri){ @@ -929,13 +1053,14 @@ complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed, } fullvec_offsetC += otC->bspecn; } - fullvec_offsetR += otR->bspecn; } + fullvec_offsetR += otR->bspecn; } free(tmp); return target_packed; } + /// Transforms a big "packed" matrix into the full basis. /** TODO doc */ complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full, diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index be503f3..439756c 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -426,6 +426,26 @@ qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const st /// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load. void qpms_scatsys_free(qpms_scatsys_t *s); +/// Creates a "full" transformation matrix U that takes a full vector and projects it onto an symmetry adapted basis. +/** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient. + * + * TODO doc about shape etc. + */ +complex double *qpms_scatsys_irrep_transform_matrix(complex double *target_U, + const qpms_scatsys_t *ss, qpms_iri_t iri); + +/// Projects a "big" matrix onto an irrep (slow reference implementation). +/** TODO doc */ +complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed, + const complex double *orig_full, const qpms_scatsys_t *ss, + qpms_iri_t iri); + +/// Transforms a big "packed" matrix into the full basis (slow reference implementation). +/** TODO doc */ +complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full, + const complex double *orig_packed, const qpms_scatsys_t *ss, + qpms_iri_t iri, bool add); + /// Projects a "big" matrix onto an irrep. /** TODO doc */ complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed, diff --git a/qpms/symmetries.py b/qpms/symmetries.py index 67fbb6f..d8cffc3 100644 --- a/qpms/symmetries.py +++ b/qpms/symmetries.py @@ -53,7 +53,7 @@ class SVWFPointGroupInfo: # only for point groups, coz in svwf_rep() I use I_tyt self.svwf_rep_gen_func = svwf_rep_gen_func self.irreps = dict() for irrepname, irrepgens in irrepgens_dict.items(): - is1d = isinstance(irrepgens[0], int) + is1d = isinstance(irrepgens[0], (int,float,complex)) irrepdim = 1 if is1d else irrepgens[0].shape[0] self.irreps[irrepname] = generate_grouprep(self.permgroup, 1 if is1d else np.eye(irrepdim), @@ -495,6 +495,24 @@ point_group_info = { # representation info of some useful point groups qpms.IRot3.identity(), ) ), + 'C2' : SVWFPointGroupInfo('C2', + # permutation group generators + (Permutation(0,1), # 180 deg rotation around z axis + ), + # dictionary with irrep generators + { + # Bradley, Cracknell p. 57; + 'A': (1,), + 'B': (-1,), + }, + # function that generates a tuple with svwf representation generators + lambda lMax : (qpms.zrotN_tyty(2, lMax),), + # quaternion rep generators + rep3d_gens = ( + qpms.IRot3.zrotN(2), + ) + + ), 'C2v' : SVWFPointGroupInfo('C2v', # permutation group generators (Permutation(0,1, size=4)(2,3), # x -> - x mirror operation (i.e. yz mirror plane) @@ -546,6 +564,24 @@ point_group_info = { # representation info of some useful point groups qpms.IRot3.zflip(), ) ), + 'C4' : SVWFPointGroupInfo('C4', + # permutation group generators + (Permutation(0,1,2,3, size=4),), #C4 rotation + # dictionary with irrep generators + { + # Bradley, Cracknell p. 58 + 'A': (1,), + 'B': (-1,), + '1E': (-1j,), + '2E': (1j,), + }, + # function that generates a tuple with svwf representation generators + lambda lMax : (qpms.zrotN_tyty(4, lMax), ), + # quaternion rep generators + rep3d_gens = ( + qpms.IRot3.zrotN(4), + ) + ), 'C4v' : SVWFPointGroupInfo('C4v', # permutation group generators (Permutation(0,1,2,3, size=4), #C4 rotation @@ -621,4 +657,41 @@ point_group_info = { # representation info of some useful point groups qpms.IRot3.zflip(), ) ), + 'x_and_z_flip': SVWFPointGroupInfo( + 'x_and_z_flip', + ( + Permutation(0,1, size=4), # x -> -x mirror op + Permutation(2,3, size=4), # z -> -z mirror op + ), + { + "P'": (1, 1), + "R'": (-1, 1), + "P''": (-1,-1), + "R''": (1, -1), + }, + lambda lMax : (qpms.xflip_tyty(lMax), qpms.zflip_tyty(lMax)), + rep3d_gens = ( + qpms.IRot3.xflip(), + qpms.IRot3.zflip(), + ) + + ), + 'y_and_z_flip': SVWFPointGroupInfo( + 'y_and_z_flip', + ( + Permutation(0,1, size=4), # y -> -y mirror op + Permutation(2,3, size=4), # z -> -z mirror op + ), + { + "P'": (1, 1), + "R'": (-1, 1), + "P''": (-1,-1), + "R''": (1, -1), + }, + lambda lMax : (qpms.yflip_tyty(lMax), qpms.zflip_tyty(lMax)), + rep3d_gens = ( + qpms.IRot3.yflip(), + qpms.IRot3.zflip(), + ) + ), } diff --git a/qpms/tmatrix_io.c b/qpms/tmatrix_io.c index 2108de9..ac77e26 100644 --- a/qpms/tmatrix_io.c +++ b/qpms/tmatrix_io.c @@ -33,9 +33,9 @@ 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) // CHECKME CORRECT? + if (bs->norm != QPMS_NORMALISATION_POWER_CS) // CHECKME CORRECT? qpms_pr_error_at_flf(__FILE__, __LINE__, __func__, - "Not implemented; only QPMS_NORMALISATION_POWER (CHECKME)" + "Not implemented; only QPMS_NORMALISATION_POWER_CS (CHECKME)" " norm supported right now."); int n_alloc = 128; // First chunk to allocate *n = 0; diff --git a/qpms/translations.c b/qpms/translations.c index db7af4d..4dbdbea 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -229,8 +229,14 @@ complex double qpms_trans_single_A(qpms_normalisation_t norm, #elif defined AN3 * ipow (nu - n) #endif +#ifdef AM1 + * ipow(-m+mu) +#endif //NNU #ifdef AM2 * min1pow(-m+mu) +#endif //NNU +#ifdef AM3 + * ipow(m-mu) #endif //NNU ; return presum * sum; @@ -409,8 +415,14 @@ complex double qpms_trans_single_B(qpms_normalisation_t norm, #elif defined AN3 * ipow (nu - n) #endif +#ifdef AM1 + * ipow(-m+mu) +#endif //NNU #ifdef AM2 * min1pow(-m+mu) +#endif //NNU +#ifdef AM3 + * ipow(m-mu) #endif //NNU ; @@ -559,8 +571,14 @@ static void qpms_trans_calculator_multipliers_A_general( #elif defined AN3 * ipow (nu - n) #endif +#ifdef AM1 + * ipow(-m+mu) +#endif //NNU #ifdef AM2 * min1pow(-m+mu) +#endif //NNU +#ifdef AM3 + * ipow(m-mu) #endif //NNU ; // FIXME I might not need complex here @@ -615,9 +633,15 @@ void qpms_trans_calculator_multipliers_B_general( #elif defined BN3 * ipow (nu - n) #endif +#ifdef BM1 + * ipow(-m+mu) +#endif #ifdef BM2 * min1pow(-m+mu) #endif +#ifdef BM3 + * ipow(m-mu) +#endif #ifdef BF1 * I #elif defined BF2 diff --git a/qpms/wigner.h b/qpms/wigner.h index d37c031..c98de1b 100644 --- a/qpms/wigner.h +++ b/qpms/wigner.h @@ -152,8 +152,8 @@ static inline cart3_t qpms_quat_rot_cart3(qpms_quat_t q, const cart3_t v) { //const qpms_quat_t qc = qpms_quat_normalise(qpms_quat_pow(q, -1)); // implementation of _pow wrong! const qpms_quat_t qc = qpms_quat_conj(q); const qpms_quat_t vv = qpms_quat_from_cart3(v); - return qpms_quat_to_cart3(qpms_quat_mult(qc, - qpms_quat_mult(vv, q))); + return qpms_quat_to_cart3(qpms_quat_mult(q, + qpms_quat_mult(vv, qc))); } /// Wigner D matrix element from a rotator quaternion for integer \a l. diff --git a/setup.py b/setup.py index 9993d9d..51696a0 100755 --- a/setup.py +++ b/setup.py @@ -31,11 +31,14 @@ qpms_c = Extension('qpms_c', 'qpms/error.c'], extra_compile_args=['-std=c99','-ggdb', '-O0', '-DQPMS_COMPILE_PYTHON_EXTENSIONS', # this is required + #'-DAN2', + #'-DBN2', #'-DQPMS_USE_OMP', '-DDISABLE_NDEBUG', # uncomment to enable assertions in the modules #'-fopenmp', ], - libraries=['gsl', 'lapacke', 'blas', 'gslcblas', #'omp' + libraries=['gsl', 'lapacke', 'cblas_LINUX', 'blas_LINUX', + '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 [] diff --git a/tests/fixproj-notes b/tests/fixproj-notes new file mode 100644 index 0000000..6b14a7e --- /dev/null +++ b/tests/fixproj-notes @@ -0,0 +1,4 @@ +c99 -g -DDUMP_ORBIT_ACTION -DDUMP_PROJECTORMATRIX -DDUMP_ACTIONMATRIX -I.. sss1.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c ../qpms/wigner.c -lm -lgsl -llapacke ~/repo/CBLAS/lib/cblas_LINUX.a ~/repo/BLAS-3.8.0/blas_LINUX.a + + +LD_LIBRARY_PATH=$HOME/repo/CBLAS/lib/ LDFLAGS=-L../CBLAS/lib/ python3 setup.py install --user diff --git a/tests/generate_groupcode.py b/tests/generate_groupcode.py index b69e675..846e765 100755 --- a/tests/generate_groupcode.py +++ b/tests/generate_groupcode.py @@ -2,7 +2,8 @@ from qpms.symmetries import point_group_info codestring = "#include \n" -for name, info in point_group_info.items(): +for name in sorted(point_group_info.keys()): + info = point_group_info[name] codestring += 'const qpms_finite_group_t QPMS_FINITE_GROUP_%s = ' %name codestring += info.generate_c_source() codestring += ";\n\n" diff --git a/tests/quaternion_O3.c b/tests/quaternion_O3.c index 7b880b3..da234eb 100644 --- a/tests/quaternion_O3.c +++ b/tests/quaternion_O3.c @@ -3,14 +3,18 @@ int main() { cart3_t v = {1, 2, 3}; - qpms_quat4d_t q4[4] = { + qpms_quat4d_t q4[8] = { {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, + {-1, 0, 0, 0}, + {0, -1, 0, 0}, + {0, 0, -1, 0}, + {0, 0, 0, -1}, }; printf("original: (%g, %g, %g)\n", v.x, v.y, v.z); - for(int i = 0; i < 4; ++i) { + for(int i = 0; i < 8; ++i) { for (int det = -1; det < 2; det += 2) { const qpms_quat4d_t qr = q4[i]; const qpms_quat_t qc = qpms_quat_2c_from_4d(qr); diff --git a/tests/scatsys.py b/tests/scatsys.py index 7206cd2..eedb12e 100644 --- a/tests/scatsys.py +++ b/tests/scatsys.py @@ -2,6 +2,8 @@ from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSyste from qpms.symmetries import point_group_info import numpy as np +np.random.seed(444) + sym = FinitePointGroup(point_group_info['D3h']) bspec2 = BaseSpec(lMax=2) bspec1 = BaseSpec(lMax=1) diff --git a/tests/ss_syms_packs.c b/tests/ss_syms_packs.c index 470fc5f..1af2c09 100644 --- a/tests/ss_syms_packs.c +++ b/tests/ss_syms_packs.c @@ -35,20 +35,27 @@ int main() for(qpms_l_t l = 1; l <= 1; ++l) for (qpms_m_t m = -l; m <= l; ++m) qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); - for(qpms_l_t l = 1; l <= 2; ++l) + for(qpms_l_t l = 1; l <= 1; ++l) for (qpms_m_t m = -l; m <= l; ++m) qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); #endif qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1); qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2); +#if 0 // Random diagonal T-matrices for(size_t i = 0; i < b1->n; ++i) t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1); for(size_t i = 0; i < b2->n; ++i) t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1); +#else + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = 1; + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = 1; +#endif - const cart3_t pp1 = {1, 2, 0}, pp2 = {1, 2, 3}, pp3 = {0.1, 2, 0}; + const cart3_t pp1 = {0, 0, 1}, pp2 = {0,0, 2}, pp3 = {0,0 , 0}; qpms_tmatrix_t * tmlist[] = {t1, t2}; qpms_particle_tid_t plist[] = {{pp1, 0}, {pp2, 0}, {pp3, 1}}; @@ -58,7 +65,7 @@ int main() protoss.p = plist; protoss.p_count=3; - qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D4h); + qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D3h); printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n", (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps, diff --git a/tests/sss1.c b/tests/sss1.c new file mode 100644 index 0000000..1e94b95 --- /dev/null +++ b/tests/sss1.c @@ -0,0 +1,105 @@ +// c99 -g -I.. ss_syms_packs.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c ../qpms/wigner.c -lm -lgsl -lblas -llapacke +typedef int qpms_gmi_t;// There is something wrong in the includes, apparently. +#include +#include +#include +#include +#include +#include +#include "staticgroups.h" + +const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h; +const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v; +const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g; +const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v; +const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h; +const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h; + +double uniform_random(double min, double max) { + double random_value = min + (max-min)*(double)rand()/RAND_MAX; + return random_value; +} + +int main() +{ + srand(666); +#if 0 + qpms_vswf_set_spec_t + *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS), + *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS); +#else + // Only electric waves + qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(), + *b2 = qpms_vswf_set_spec_init(); + b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS; + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -l; m <= l; ++m) + qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -l; m <= l; ++m) + qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); +#endif + qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1); + qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2); + +#if 0 + // Random diagonal T-matrices + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1); + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1); +#else + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = 1; + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = 1; +#endif + + const cart3_t pp1 = {1, 1, 0}; + qpms_tmatrix_t * tmlist[] = {t1, t2}; + qpms_particle_tid_t plist[] = {{pp1, 0}}; + + qpms_scatsys_t protoss; + protoss.tm = tmlist; + protoss.tm_count=2; + protoss.p = plist; + protoss.p_count=1; + + qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D2h); + + printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n", + (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps, + (int)ss->orbit_type_count); + + const double k = 1.7; + + complex double *S_full = qpms_scatsys_build_translation_matrix_full( + NULL, ss, k); + complex double *S_packed[ss->sym->nirreps]; + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) + S_packed[iri] = qpms_scatsys_irrep_pack_matrix(NULL, + S_full, ss, iri); + + complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix(NULL, + S_packed[0], ss, 0, false); + for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri) + qpms_scatsys_irrep_unpack_matrix(S_recfull, S_packed[iri], + ss, iri, true); + + double maxerr = 0; + for (size_t i = 0; i < ss->fecv_size; ++i) { + double err = cabs(S_full[i] - S_recfull[i]); + maxerr = (err > maxerr) ? err : maxerr; + } + + printf("maxerr: %lg\n", maxerr); + + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]); + free(S_full); + qpms_scatsys_free(ss); + qpms_tmatrix_free(t1); + qpms_tmatrix_free(t2); + qpms_vswf_set_spec_free(b1); + qpms_vswf_set_spec_free(b2); + return 0; +} diff --git a/tests/sss2.c b/tests/sss2.c new file mode 100644 index 0000000..54c959d --- /dev/null +++ b/tests/sss2.c @@ -0,0 +1,180 @@ +// c99 -g -DZLINE -DDAGRUP=D3h -DDUMP_PARTICLE_POSITIONS -DDUMP_ORBIT_ACTION -DDUMP_PROJECTORMATRIX -DDUMP_ACTIONMATRIX -I.. sss2.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c ../qpms/wigner.c -lm -lgsl -lblas -llapacke +typedef int qpms_gmi_t;// There is something wrong in the includes, apparently. +#include +#include +#include +#include +#include +#include +#include "staticgroups.h" + +const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h; +const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v; +const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g; +const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v; +const qpms_finite_group_t *C2 = &QPMS_FINITE_GROUP_C2; +const qpms_finite_group_t *C4 = &QPMS_FINITE_GROUP_C4; +const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h; +const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h; +const qpms_finite_group_t *x_and_z_flip = &QPMS_FINITE_GROUP_x_and_z_flip; +const qpms_finite_group_t *y_and_z_flip = &QPMS_FINITE_GROUP_y_and_z_flip; + +#ifndef DAGRUP +#define DAGRUP D4h +#endif + +double uniform_random(double min, double max) { + double random_value = min + (max-min)*(double)rand()/RAND_MAX; + return random_value; +} + +int main() +{ + srand(666); +#if 0 + qpms_vswf_set_spec_t + *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS), + *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS); +#else + // Only electric waves + qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(), + *b2 = qpms_vswf_set_spec_init(); + b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS; + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -0l; m <= l; m += 2) + qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -0l; m <= l; m += 2) + qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); +#endif + qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1); + qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2); + +#if 0 + // Random diagonal T-matrices + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1); + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1); +#else + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = 1; + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = 1; +#endif + +#ifdef YLINE + const cart3_t pp1 = {0, 1.1, 0}; + const cart3_t pp2 = {0, 1.4, 0}; +#elif defined XLINE + const cart3_t pp1 = {1.1, 0, 0}; + const cart3_t pp2 = {1.4, 0, 0}; +#elif defined ZLINE + const cart3_t pp1 = {0, 0, 1.1}; + const cart3_t pp2 = {0, 0, 1.4}; +#else + const cart3_t pp1 = {1.1, 1, 0}; + const cart3_t pp2 = {0, 1.4, 0}; +#endif + const cart3_t pp3 = {0, 0, 1}; + qpms_tmatrix_t * tmlist[] = {t1, t2}; + qpms_particle_tid_t plist[] = {{pp1,1}, {pp2, 0}, {pp3, 1}, + }; + + qpms_scatsys_t protoss; + protoss.tm = tmlist; + protoss.tm_count=sizeof(tmlist)/sizeof(qpms_tmatrix_t *); + protoss.p = plist; + protoss.p_count=sizeof(plist)/sizeof(qpms_particle_tid_t); + + qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, DAGRUP); + + printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n", + (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps, + (int)ss->orbit_type_count); + + const double k = 1.7; + + complex double *S_full = qpms_scatsys_build_translation_matrix_full( + NULL, ss, k); + { + const size_t full_len = ss->fecv_size; + size_t fullvec_offset_dest = 0; + for (qpms_ss_pi_t pdest = 0; pdest < ss->p_count; pdest++) { + size_t fullvec_offset_src = 0; + const size_t bspecn_dest = ss->tm[ss->p[pdest].tmatrix_id]->spec->n; + for (qpms_ss_pi_t psrc = 0; psrc < ss->p_count; psrc++) { + const size_t bspecn_src = ss->tm[ss->p[psrc].tmatrix_id]->spec->n; + fprintf(stderr, "Translation matrix element %d<-%d; (%g %g %g)<-(%g %g %g):\n", + (int)pdest, (int)psrc, ss->p[pdest].pos.x, ss->p[pdest].pos.y, ss->p[pdest].pos.z, + ss->p[psrc].pos.x, ss->p[psrc].pos.y, ss->p[psrc].pos.z); + + for(size_t row = 0; row < bspecn_dest; ++row) { + for(size_t col = 0; col < bspecn_src; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]), + cimag(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col])); + fputc('\n', stderr); + } + fullvec_offset_src += bspecn_src; + } + fullvec_offset_dest += bspecn_dest; + } + } + { + fputs("\n\n", stderr); + const size_t full_len = ss->fecv_size; + for (size_t row = 0 ; row < full_len; ++row) { + for (size_t col = 0 ; col < full_len; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * row + col]), cimag(S_full[full_len * row + col])); + fputc('\n', stderr); + } + } + + complex double *S_packed[ss->sym->nirreps]; + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) + S_packed[iri] = qpms_scatsys_irrep_pack_matrix(NULL, + S_full, ss, iri); + + complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix(NULL, + S_packed[0], ss, 0, false); + for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri) + qpms_scatsys_irrep_unpack_matrix(S_recfull, S_packed[iri], + ss, iri, true); + { + fputs("\n\n", stderr); + const size_t full_len = ss->fecv_size; + for (size_t row = 0 ; row < full_len; ++row) { + for (size_t col = 0 ; col < full_len; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_recfull[full_len * row + col]), cimag(S_recfull[full_len * row + col])); + fputc('\n', stderr); + } + } + + double maxerr = 0; + for (size_t i = 0; i < ss->fecv_size; ++i) { + double err = cabs(S_full[i] - S_recfull[i]); + maxerr = (err > maxerr) ? err : maxerr; + } + + + printf("maxerr: %lg\n", maxerr); + + fprintf(stderr, "pi\tpos\toti\tosn\tp\n"); + for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { + cart3_t pos = ss->p[pi].pos; + qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t; + qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn; + qpms_ss_orbit_pi_t p = ss->p_orbitinfo[pi].p; + fprintf(stderr, "%d\t(%.3g,%.3g,%.3g)\t%d\t%d\t%d\n", + (int)pi, pos.x, pos.y, pos.z, (int)oti, (int)osn, (int)p); + } + + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]); + free(S_full); + qpms_scatsys_free(ss); + qpms_tmatrix_free(t1); + qpms_tmatrix_free(t2); + qpms_vswf_set_spec_free(b1); + qpms_vswf_set_spec_free(b2); + return 0; +} diff --git a/tests/sss3.c b/tests/sss3.c new file mode 100644 index 0000000..9f07128 --- /dev/null +++ b/tests/sss3.c @@ -0,0 +1,236 @@ +// c99 -g -DNLINE -DDAGRUP=C4v -DDUMP_PARTICLE_POSITIONS -DDUMP_ORBIT_ACTION -DDUMP_PROJECTORMATRIX -DDUMP_ACTIONMATRIX -I.. sss3.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c ../qpms/wigner.c -lm -lgsl -llapacke ~/repo/CBLAS/lib/cblas_LINUX.a ~/repo/BLAS-3.8.0/blas_LINUX.a +typedef int qpms_gmi_t;// There is something wrong in the includes, apparently. +#include +#include +#include +#include +#include +#include +#include "staticgroups.h" + +const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h; +const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v; +const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g; +const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v; +const qpms_finite_group_t *C2 = &QPMS_FINITE_GROUP_C2; +const qpms_finite_group_t *C4 = &QPMS_FINITE_GROUP_C4; +const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h; +const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h; +const qpms_finite_group_t *x_and_z_flip = &QPMS_FINITE_GROUP_x_and_z_flip; +const qpms_finite_group_t *y_and_z_flip = &QPMS_FINITE_GROUP_y_and_z_flip; + +#ifndef DAGRUP +#define DAGRUP D4h +#endif + +double uniform_random(double min, double max) { + double random_value = min + (max-min)*(double)rand()/RAND_MAX; + return random_value; +} + +int main() +{ + srand(666); +#if 0 + qpms_vswf_set_spec_t + *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS), + *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS); +#else + // Only electric waves + qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(), + *b2 = qpms_vswf_set_spec_init(); + b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS; + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -0l; m <= l; m += 2) + qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -0l; m <= l; m += 2) + qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); +#endif + qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1); + qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2); + +#if 0 + // Random diagonal T-matrices + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1); + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1); +#else + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = 1; + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = 1; +#endif + +#ifdef YLINE + const cart3_t pp1 = {0, 1.1, 0}; + const cart3_t pp2 = {0, 1.4, 0}; +#elif defined XLINE + const cart3_t pp1 = {1.1, 0, 0}; + const cart3_t pp2 = {1.4, 0, 0}; +#elif defined ZLINE + const cart3_t pp1 = {0, 0, 1.1}; + const cart3_t pp2 = {0, 0, 1.4}; +#else + const cart3_t pp1 = {1.1, 1, 0}; + const cart3_t pp2 = {0, 1.4, 0}; +#endif + const cart3_t pp3 = {0, 0, 1}; + qpms_tmatrix_t * tmlist[] = {t1, t2}; + qpms_particle_tid_t plist[] = {{pp1,1}, {pp2, 0}, {pp3, 1}, + }; + + qpms_scatsys_t protoss; + protoss.tm = tmlist; + protoss.tm_count=sizeof(tmlist)/sizeof(qpms_tmatrix_t *); + protoss.p = plist; + protoss.p_count=sizeof(plist)/sizeof(qpms_particle_tid_t); + + qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, DAGRUP); + + printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n", + (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps, + (int)ss->orbit_type_count); + + fputs("Orbit projection matrices:\n", stderr); + for (qpms_ss_oti_t oti = 0; oti < ss->orbit_type_count; ++oti) { + fprintf(stderr, "Orbit type %d:\n", (int)oti); + const qpms_ss_orbit_type_t *ot = &(ss->orbit_types[oti]); + size_t row = 0; + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) { + assert(row*ot->size*ot->bspecn == ot->irbase_offsets[iri]); + fprintf(stderr, "---------- IR %d (%s) -------------\n", (int)iri, ss->sym->irreps[iri].name); + for (size_t irbi = 0; irbi < ot->irbase_sizes[iri]; ++irbi) { + for(qpms_ss_orbit_pi_t opi = 0; opi < ot->size; ++opi){ + fputs("| ", stderr); + for(size_t i = 0; i < ot->bspecn; ++i) { + const complex double elem = ot->irbases[row * ot->size * ot->bspecn + + opi * ot->bspecn + i]; + fprintf(stderr, "%+.3f%+.3fj ", creal(elem), cimag(elem)); + } + } + fputs("|\n", stderr); + ++row; + } + } + fputs("------------------------\n\n",stderr); + } + + const double k = 1.7; + + complex double *S_full = qpms_scatsys_build_translation_matrix_full( + NULL, ss, k); + { + const size_t full_len = ss->fecv_size; + size_t fullvec_offset_dest = 0; + for (qpms_ss_pi_t pdest = 0; pdest < ss->p_count; pdest++) { + size_t fullvec_offset_src = 0; + const size_t bspecn_dest = ss->tm[ss->p[pdest].tmatrix_id]->spec->n; + for (qpms_ss_pi_t psrc = 0; psrc < ss->p_count; psrc++) { + const size_t bspecn_src = ss->tm[ss->p[psrc].tmatrix_id]->spec->n; + fprintf(stderr, "Translation matrix element %d<-%d; (%g %g %g)<-(%g %g %g):\n", + (int)pdest, (int)psrc, ss->p[pdest].pos.x, ss->p[pdest].pos.y, ss->p[pdest].pos.z, + ss->p[psrc].pos.x, ss->p[psrc].pos.y, ss->p[psrc].pos.z); + + for(size_t row = 0; row < bspecn_dest; ++row) { + for(size_t col = 0; col < bspecn_src; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]), + cimag(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col])); + fputc('\n', stderr); + } + fullvec_offset_src += bspecn_src; + } + fullvec_offset_dest += bspecn_dest; + } + } + { + fputs("\n\n", stderr); + const size_t full_len = ss->fecv_size; + for (size_t row = 0 ; row < full_len; ++row) { + for (size_t col = 0 ; col < full_len; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * row + col]), cimag(S_full[full_len * row + col])); + fputc('\n', stderr); + } + } + + complex double *S_packed[ss->sym->nirreps]; + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) { + S_packed[iri] = qpms_scatsys_irrep_pack_matrix(NULL, + S_full, ss, iri); + fprintf(stderr, "--- Packed matrix for irrep %d (%s):\n", (int) iri, ss->sym->irreps[iri].name); + for (size_t row = 0; row < ss->saecv_sizes[iri]; ++row) { + for (size_t col = 0; col < ss->saecv_sizes[iri]; ++col) { + complex double elem = S_packed[iri][row * ss->saecv_sizes[iri] + col]; + fprintf(stderr, "%+.3f+%.3fj ", creal(elem), cimag(elem)); + } + fputc('\n', stderr); + } + } + { + complex double *S_partrecfull = qpms_scatsys_irrep_unpack_matrix(NULL, + S_packed[0], ss, 0, false); + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) { + qpms_scatsys_irrep_unpack_matrix(S_partrecfull, S_packed[iri], + ss, iri, false); + fprintf(stderr, "\nPartial reconstruction %d (%s):\n", (int)iri, ss->sym->irreps[iri].name); + const size_t full_len = ss->fecv_size; + for (size_t row = 0 ; row < full_len; ++row) { + for (size_t col = 0 ; col < full_len; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_partrecfull[full_len * row + col]), cimag(S_partrecfull[full_len * row + col])); + fputc('\n', stderr); + } + } + + double maxerr = 0; + for (size_t i = 0; i < ss->fecv_size; ++i) { + double err = cabs(S_full[i] - S_partrecfull[i]); + maxerr = (err > maxerr) ? err : maxerr; + } + free(S_partrecfull); + } + + + complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix(NULL, + S_packed[0], ss, 0, false); + for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri) + qpms_scatsys_irrep_unpack_matrix(S_recfull, S_packed[iri], + ss, iri, true); + { + fputs("\n\n", stderr); + const size_t full_len = ss->fecv_size; + for (size_t row = 0 ; row < full_len; ++row) { + for (size_t col = 0 ; col < full_len; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_recfull[full_len * row + col]), cimag(S_recfull[full_len * row + col])); + fputc('\n', stderr); + } + } + + double maxerr = 0; + for (size_t i = 0; i < ss->fecv_size; ++i) { + double err = cabs(S_full[i] - S_recfull[i]); + maxerr = (err > maxerr) ? err : maxerr; + } + + + printf("maxerr: %lg\n", maxerr); + + fprintf(stderr, "pi\tpos\toti\tosn\tp\n"); + for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { + cart3_t pos = ss->p[pi].pos; + qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t; + qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn; + qpms_ss_orbit_pi_t p = ss->p_orbitinfo[pi].p; + fprintf(stderr, "%d\t(%.3g,%.3g,%.3g)\t%d\t%d\t%d\n", + (int)pi, pos.x, pos.y, pos.z, (int)oti, (int)osn, (int)p); + } + + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]); + free(S_full); + qpms_scatsys_free(ss); + qpms_tmatrix_free(t1); + qpms_tmatrix_free(t2); + qpms_vswf_set_spec_free(b1); + qpms_vswf_set_spec_free(b2); + return 0; +} diff --git a/tests/sss4.c b/tests/sss4.c new file mode 100644 index 0000000..5a50ac8 --- /dev/null +++ b/tests/sss4.c @@ -0,0 +1,237 @@ +// c99 -g -DNLINE -DDAGRUP=C4v -DDUMP_PARTICLE_POSITIONS -DDUMP_ORBIT_ACTION -DDUMP_PROJECTORMATRIX -DDUMP_ACTIONMATRIX -I.. sss3.c staticgroups.c ../qpms/scatsystem.c ../qpms/vswf.c ../qpms/error.c ../qpms/translations.c ../qpms/symmetries.c ../qpms/legendre.c ../qpms/gaunt.c ../qpms/wigner.c -lm -lgsl -llapacke ~/repo/CBLAS/lib/cblas_LINUX.a ~/repo/BLAS-3.8.0/blas_LINUX.a +typedef int qpms_gmi_t;// There is something wrong in the includes, apparently. +#include +#include +#include +#include +#include +#include +#include "staticgroups.h" + +const qpms_finite_group_t *D3h = &QPMS_FINITE_GROUP_D3h; +const qpms_finite_group_t *C4v = &QPMS_FINITE_GROUP_C4v; +const qpms_finite_group_t *TRIVG = &QPMS_FINITE_GROUP_trivial_g; +const qpms_finite_group_t *C2v = &QPMS_FINITE_GROUP_C2v; +const qpms_finite_group_t *C2 = &QPMS_FINITE_GROUP_C2; +const qpms_finite_group_t *C4 = &QPMS_FINITE_GROUP_C4; +const qpms_finite_group_t *D2h = &QPMS_FINITE_GROUP_D2h; +const qpms_finite_group_t *D4h = &QPMS_FINITE_GROUP_D4h; +const qpms_finite_group_t *x_and_z_flip = &QPMS_FINITE_GROUP_x_and_z_flip; +const qpms_finite_group_t *y_and_z_flip = &QPMS_FINITE_GROUP_y_and_z_flip; + +#ifndef DAGRUP +#define DAGRUP D4h +#endif + +double uniform_random(double min, double max) { + double random_value = min + (max-min)*(double)rand()/RAND_MAX; + return random_value; +} + +int main() +{ + srand(666); +#if 1 + qpms_vswf_set_spec_t + *b1 = qpms_vswf_set_spec_from_lMax(1,QPMS_NORMALISATION_POWER_CS), + *b2 = qpms_vswf_set_spec_from_lMax(2,QPMS_NORMALISATION_POWER_CS); +#else + // Only electric waves + qpms_vswf_set_spec_t *b1 = qpms_vswf_set_spec_init(), + *b2 = qpms_vswf_set_spec_init(); + b1->norm = b2-> norm = QPMS_NORMALISATION_POWER_CS; + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -l; m <= l; m += 1) + qpms_vswf_set_spec_append(b1, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); + for(qpms_l_t l = 1; l <= 1; ++l) + for (qpms_m_t m = -l; m <= l; m += 1) + qpms_vswf_set_spec_append(b2, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, m, l)); +#endif + qpms_tmatrix_t *t1 = qpms_tmatrix_init(b1); + qpms_tmatrix_t *t2 = qpms_tmatrix_init(b2); + +#if 0 + // Random diagonal T-matrices + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = uniform_random(-1,1) + I*uniform_random(-1,1); + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = uniform_random(-1,1) + I*uniform_random(-1,1); +#else + for(size_t i = 0; i < b1->n; ++i) + t1->m[i + i*b1->n] = 1; + for(size_t i = 0; i < b2->n; ++i) + t2->m[i + i*b2->n] = 1; +#endif + +#ifdef YLINE + const cart3_t pp1 = {0, 1.1, 0}; + const cart3_t pp2 = {0, 1.4, 0}; +#elif defined XLINE + const cart3_t pp1 = {1.1, 0, 0}; + const cart3_t pp2 = {1.4, 0, 0}; +#elif defined ZLINE + const cart3_t pp1 = {0, 0, 1.1}; + const cart3_t pp2 = {0, 0, 1.4}; +#else + const cart3_t pp1 = {1.1, 1, 0}; + const cart3_t pp2 = {0, 1.4, 0}; + const cart3_t pp4 = {2, 1.4, 0}; +#endif + const cart3_t pp3 = {0, 0, 0}; + qpms_tmatrix_t * tmlist[] = {t1, t2}; + qpms_particle_tid_t plist[] = {{pp1,1}, {pp2, 0}, {pp3, 1}, {pp4,1} + }; + + qpms_scatsys_t protoss; + protoss.tm = tmlist; + protoss.tm_count=sizeof(tmlist)/sizeof(qpms_tmatrix_t *); + protoss.p = plist; + protoss.p_count=sizeof(plist)/sizeof(qpms_particle_tid_t); + + qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, DAGRUP); + + printf("p_count: %d, tm_count: %d, nirreps: %d, orbit_type_count: %d\n", + (int)ss->p_count, (int)ss->tm_count, (int)ss->sym->nirreps, + (int)ss->orbit_type_count); + + fputs("Orbit projection matrices:\n", stderr); + for (qpms_ss_oti_t oti = 0; oti < ss->orbit_type_count; ++oti) { + fprintf(stderr, "Orbit type %d:\n", (int)oti); + const qpms_ss_orbit_type_t *ot = &(ss->orbit_types[oti]); + size_t row = 0; + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) { + assert(row*ot->size*ot->bspecn == ot->irbase_offsets[iri]); + fprintf(stderr, "---------- IR %d (%s) -------------\n", (int)iri, ss->sym->irreps[iri].name); + for (size_t irbi = 0; irbi < ot->irbase_sizes[iri]; ++irbi) { + for(qpms_ss_orbit_pi_t opi = 0; opi < ot->size; ++opi){ + fputs("| ", stderr); + for(size_t i = 0; i < ot->bspecn; ++i) { + const complex double elem = ot->irbases[row * ot->size * ot->bspecn + + opi * ot->bspecn + i]; + fprintf(stderr, "%+.3f%+.3fj ", creal(elem), cimag(elem)); + } + } + fputs("|\n", stderr); + ++row; + } + } + fputs("------------------------\n\n",stderr); + } + + const double k = 1.7; + + complex double *S_full = qpms_scatsys_build_translation_matrix_full( + NULL, ss, k); + { + const size_t full_len = ss->fecv_size; + size_t fullvec_offset_dest = 0; + for (qpms_ss_pi_t pdest = 0; pdest < ss->p_count; pdest++) { + size_t fullvec_offset_src = 0; + const size_t bspecn_dest = ss->tm[ss->p[pdest].tmatrix_id]->spec->n; + for (qpms_ss_pi_t psrc = 0; psrc < ss->p_count; psrc++) { + const size_t bspecn_src = ss->tm[ss->p[psrc].tmatrix_id]->spec->n; + fprintf(stderr, "Translation matrix element %d<-%d; (%g %g %g)<-(%g %g %g):\n", + (int)pdest, (int)psrc, ss->p[pdest].pos.x, ss->p[pdest].pos.y, ss->p[pdest].pos.z, + ss->p[psrc].pos.x, ss->p[psrc].pos.y, ss->p[psrc].pos.z); + + for(size_t row = 0; row < bspecn_dest; ++row) { + for(size_t col = 0; col < bspecn_src; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col]), + cimag(S_full[full_len * (fullvec_offset_dest+row) + fullvec_offset_src+col])); + fputc('\n', stderr); + } + fullvec_offset_src += bspecn_src; + } + fullvec_offset_dest += bspecn_dest; + } + } + { + fputs("\n\n", stderr); + const size_t full_len = ss->fecv_size; + for (size_t row = 0 ; row < full_len; ++row) { + for (size_t col = 0 ; col < full_len; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_full[full_len * row + col]), cimag(S_full[full_len * row + col])); + fputc('\n', stderr); + } + } + + complex double *S_packed[ss->sym->nirreps]; + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) { + S_packed[iri] = qpms_scatsys_irrep_pack_matrix_stupid(NULL, + S_full, ss, iri); + fprintf(stderr, "--- Packed matrix for irrep %d (%s):\n", (int) iri, ss->sym->irreps[iri].name); + for (size_t row = 0; row < ss->saecv_sizes[iri]; ++row) { + for (size_t col = 0; col < ss->saecv_sizes[iri]; ++col) { + complex double elem = S_packed[iri][row * ss->saecv_sizes[iri] + col]; + fprintf(stderr, "%+.3f+%.3fj ", creal(elem), cimag(elem)); + } + fputc('\n', stderr); + } + } + { + complex double *S_partrecfull = qpms_scatsys_irrep_unpack_matrix_stupid(NULL, + S_packed[0], ss, 0, false); + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) { + qpms_scatsys_irrep_unpack_matrix_stupid(S_partrecfull, S_packed[iri], + ss, iri, false); + fprintf(stderr, "\nPartial reconstruction %d (%s):\n", (int)iri, ss->sym->irreps[iri].name); + const size_t full_len = ss->fecv_size; + for (size_t row = 0 ; row < full_len; ++row) { + for (size_t col = 0 ; col < full_len; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_partrecfull[full_len * row + col]), cimag(S_partrecfull[full_len * row + col])); + fputc('\n', stderr); + } + } + + double maxerr = 0; + for (size_t i = 0; i < ss->fecv_size; ++i) { + double err = cabs(S_full[i] - S_partrecfull[i]); + maxerr = (err > maxerr) ? err : maxerr; + } + free(S_partrecfull); + } + + + complex double *S_recfull = qpms_scatsys_irrep_unpack_matrix_stupid(NULL, + S_packed[0], ss, 0, false); + for (qpms_iri_t iri = 1; iri < ss->sym->nirreps; ++iri) + qpms_scatsys_irrep_unpack_matrix_stupid(S_recfull, S_packed[iri], + ss, iri, true); + { + fputs("\n\n", stderr); + const size_t full_len = ss->fecv_size; + for (size_t row = 0 ; row < full_len; ++row) { + for (size_t col = 0 ; col < full_len; ++col) + fprintf(stderr, "%+2.3f%+2.3fj ", creal(S_recfull[full_len * row + col]), cimag(S_recfull[full_len * row + col])); + fputc('\n', stderr); + } + } + + double maxerr = 0; + for (size_t i = 0; i < ss->fecv_size; ++i) { + double err = cabs(S_full[i] - S_recfull[i]); + maxerr = (err > maxerr) ? err : maxerr; + } + + + printf("maxerr: %lg\n", maxerr); + + fprintf(stderr, "pi\tpos\toti\tosn\tp\n"); + for(qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { + cart3_t pos = ss->p[pi].pos; + qpms_ss_oti_t oti = ss->p_orbitinfo[pi].t; + qpms_ss_osn_t osn = ss->p_orbitinfo[pi].osn; + qpms_ss_orbit_pi_t p = ss->p_orbitinfo[pi].p; + fprintf(stderr, "%d\t(%.3g,%.3g,%.3g)\t%d\t%d\t%d\n", + (int)pi, pos.x, pos.y, pos.z, (int)oti, (int)osn, (int)p); + } + + for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) free(S_packed[iri]); + free(S_full); + qpms_scatsys_free(ss); + qpms_tmatrix_free(t1); + qpms_tmatrix_free(t2); + qpms_vswf_set_spec_free(b1); + qpms_vswf_set_spec_free(b2); + return 0; +} diff --git a/tests/staticgroups.c b/tests/staticgroups.c index d9240e2..3040539 100644 --- a/tests/staticgroups.c +++ b/tests/staticgroups.c @@ -1,4 +1,42 @@ -#include "staticgroups.h" +#include +const qpms_finite_group_t QPMS_FINITE_GROUP_C2 = { + "C2", // name + 2, // order + 0, // idi + (qpms_gmi_t[]) { // mt + 0, 1, + 1, 0, + }, + (qpms_gmi_t[]) { // invi + 0, 1 + }, + (qpms_gmi_t[]) {1}, // gens + 1, // ngens + (qpms_permutation_t[]){ // permrep + "(1)", + "(0 1)", + }, + NULL, // elemlabels + 2, // permrep_nelem + (qpms_irot3_t[]) { // rep3d + {{1.0+0.0*I, 0.0+0.0*I}, 1}, + {{6.123233995736766e-17+1.0*I, 0.0+0.0*I}, 1}, + }, + 2, // nirreps + (struct qpms_finite_group_irrep_t[]) { // irreps + { + 1, // dim + "B", //name + (complex double []) {1, -1} // m + }, + { + 1, // dim + "A", //name + (complex double []) {1, 1} // m + }, + } // end of irreps +}; + const qpms_finite_group_t QPMS_FINITE_GROUP_C2v = { "C2v", // name 4, // order @@ -35,6 +73,11 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_C2v = { "A1", //name (complex double []) {1, 1, 1, 1} // m }, + { + 1, // dim + "A2", //name + (complex double []) {1, -1, 1, -1} // m + }, { 1, // dim "B2", //name @@ -45,10 +88,241 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_C2v = { "B1", //name (complex double []) {1, 1, -1, -1} // m }, + } // end of irreps +}; + +const qpms_finite_group_t QPMS_FINITE_GROUP_C4 = { + "C4", // name + 4, // order + 0, // idi + (qpms_gmi_t[]) { // mt + 0, 1, 2, 3, + 1, 2, 3, 0, + 2, 3, 0, 1, + 3, 0, 1, 2, + }, + (qpms_gmi_t[]) { // invi + 0, 3, 2, 1 + }, + (qpms_gmi_t[]) {1}, // gens + 1, // ngens + (qpms_permutation_t[]){ // permrep + "(3)", + "(0 1 2 3)", + "(0 2)(1 3)", + "(0 3 2 1)", + }, + NULL, // elemlabels + 4, // permrep_nelem + (qpms_irot3_t[]) { // rep3d + {{1.0+0.0*I, 0.0+0.0*I}, 1}, + {{0.7071067811865476+0.7071067811865475*I, 0.0+0.0*I}, 1}, + {{2.220446049250313e-16+1.0*I, 0.0+0.0*I}, 1}, + {{-0.7071067811865474+0.7071067811865477*I, 0.0+0.0*I}, 1}, + }, + 4, // nirreps + (struct qpms_finite_group_irrep_t[]) { // irreps + { + 1, // dim + "B", //name + (complex double []) {1, -1, 1, -1} // m + }, + { + 1, // dim + "2E", //name + (complex double []) {1, 1j, (-1+0j), (-0-1j)} // m + }, + { + 1, // dim + "A", //name + (complex double []) {1, 1, 1, 1} // m + }, + { + 1, // dim + "1E", //name + (complex double []) {1, -1j, (-1+0j), 1j} // m + }, + } // end of irreps +}; + +const qpms_finite_group_t QPMS_FINITE_GROUP_C4v = { + "C4v", // name + 8, // order + 0, // idi + (qpms_gmi_t[]) { // mt + 0, 1, 2, 3, 4, 5, 6, 7, + 1, 2, 3, 0, 7, 4, 5, 6, + 2, 3, 0, 1, 6, 7, 4, 5, + 3, 0, 1, 2, 5, 6, 7, 4, + 4, 5, 6, 7, 0, 1, 2, 3, + 5, 6, 7, 4, 3, 0, 1, 2, + 6, 7, 4, 5, 2, 3, 0, 1, + 7, 4, 5, 6, 1, 2, 3, 0, + }, + (qpms_gmi_t[]) { // invi + 0, 3, 2, 1, 4, 5, 6, 7 + }, + (qpms_gmi_t[]) {1, 7}, // gens + 2, // ngens + (qpms_permutation_t[]){ // permrep + "(3)", + "(0 1 2 3)", + "(0 2)(1 3)", + "(0 3 2 1)", + "(3)(0 2)", + "(0 3)(1 2)", + "(1 3)", + "(0 1)(2 3)", + }, + NULL, // elemlabels + 4, // permrep_nelem + (qpms_irot3_t[]) { // rep3d + {{1.0+0.0*I, 0.0+0.0*I}, 1}, + {{0.7071067811865476+0.7071067811865475*I, 0.0+0.0*I}, 1}, + {{2.220446049250313e-16+1.0*I, 0.0+0.0*I}, 1}, + {{-0.7071067811865474+0.7071067811865477*I, 0.0+0.0*I}, 1}, + {{0.0+0.0*I, 0.7071067811865477-0.7071067811865474*I}, -1}, + {{0.0+0.0*I, 1.0+2.220446049250313e-16*I}, -1}, + {{0.0+0.0*I, 0.7071067811865475+0.7071067811865476*I}, -1}, + {{0.0+0.0*I, 0.0+1.0*I}, -1}, + }, + 5, // nirreps + (struct qpms_finite_group_irrep_t[]) { // irreps + { + 1, // dim + "A1", //name + (complex double []) {1, 1, 1, 1, 1, 1, 1, 1} // m + }, { 1, // dim "A2", //name - (complex double []) {1, -1, 1, -1} // m + (complex double []) {1, 1, 1, 1, -1, -1, -1, -1} // m + }, + { + 2, // dim + "E", //name + (complex double []) { + // (3) + 1.0, 0.0, + 0.0, 1.0, + // (0 1 2 3) + 0.0, -1.0, + 1.0, 0.0, + // (0 2)(1 3) + -1.0, 0.0, + 0.0, -1.0, + // (0 3 2 1) + 0.0, 1.0, + -1.0, 0.0, + // (3)(0 2) + 0.0, 1.0, + 1.0, 0.0, + // (0 3)(1 2) + 1.0, 0.0, + 0.0, -1.0, + // (1 3) + 0.0, -1.0, + -1.0, 0.0, + // (0 1)(2 3) + -1.0, 0.0, + 0.0, 1.0, + } + }, + { + 1, // dim + "B2", //name + (complex double []) {1, -1, 1, -1, 1, -1, 1, -1} // m + }, + { + 1, // dim + "B1", //name + (complex double []) {1, -1, 1, -1, -1, 1, -1, 1} // m + }, + } // end of irreps +}; + +const qpms_finite_group_t QPMS_FINITE_GROUP_D2h = { + "D2h", // name + 8, // order + 0, // idi + (qpms_gmi_t[]) { // mt + 0, 1, 2, 3, 4, 5, 6, 7, + 1, 0, 3, 2, 5, 4, 7, 6, + 2, 3, 0, 1, 6, 7, 4, 5, + 3, 2, 1, 0, 7, 6, 5, 4, + 4, 5, 6, 7, 0, 1, 2, 3, + 5, 4, 7, 6, 1, 0, 3, 2, + 6, 7, 4, 5, 2, 3, 0, 1, + 7, 6, 5, 4, 3, 2, 1, 0, + }, + (qpms_gmi_t[]) { // invi + 0, 1, 2, 3, 4, 5, 6, 7 + }, + (qpms_gmi_t[]) {1, 3, 7}, // gens + 3, // ngens + (qpms_permutation_t[]){ // permrep + "(5)", + "(5)(0 1)(2 3)", + "(5)(0 2)(1 3)", + "(5)(0 3)(1 2)", + "(0 3)(1 2)(4 5)", + "(0 2)(1 3)(4 5)", + "(0 1)(2 3)(4 5)", + "(4 5)", + }, + NULL, // elemlabels + 6, // permrep_nelem + (qpms_irot3_t[]) { // rep3d + {{1.0+0.0*I, 0.0+0.0*I}, 1}, + {{0.0+0.0*I, 0.0+1.0*I}, -1}, + {{0.0+1.0*I, 0.0+0.0*I}, 1}, + {{0.0+0.0*I, 1.0+0.0*I}, -1}, + {{0.0+0.0*I, 0.0+1.0*I}, 1}, + {{-1.0+0.0*I, 0.0+0.0*I}, -1}, + {{0.0+0.0*I, -1.0+0.0*I}, 1}, + {{0.0+1.0*I, 0.0+0.0*I}, -1}, + }, + 8, // nirreps + (struct qpms_finite_group_irrep_t[]) { // irreps + { + 1, // dim + "A2\'", //name + (complex double []) {1, -1, 1, -1, -1, 1, -1, 1} // m + }, + { + 1, // dim + "B1\'", //name + (complex double []) {1, 1, -1, -1, -1, -1, 1, 1} // m + }, + { + 1, // dim + "A2\'\'", //name + (complex double []) {1, 1, 1, 1, -1, -1, -1, -1} // m + }, + { + 1, // dim + "B2\'", //name + (complex double []) {1, -1, -1, 1, 1, -1, -1, 1} // m + }, + { + 1, // dim + "A1\'", //name + (complex double []) {1, 1, 1, 1, 1, 1, 1, 1} // m + }, + { + 1, // dim + "A1\'\'", //name + (complex double []) {1, -1, 1, -1, 1, -1, 1, -1} // m + }, + { + 1, // dim + "B2\'\'", //name + (complex double []) {1, 1, -1, -1, 1, 1, -1, -1} // m + }, + { + 1, // dim + "B1\'\'", //name + (complex double []) {1, -1, -1, 1, -1, 1, 1, -1} // m }, } // end of irreps }; @@ -108,6 +382,11 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D3h = { }, 6, // nirreps (struct qpms_finite_group_irrep_t[]) { // irreps + { + 1, // dim + "A2\'", //name + (complex double []) {1, 1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1} // m + }, { 2, // dim "E\'", //name @@ -150,6 +429,21 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D3h = { 0.0, 0.9999999999999996, } }, + { + 1, // dim + "A2\'\'", //name + (complex double []) {1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1} // m + }, + { + 1, // dim + "A1\'", //name + (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} // m + }, + { + 1, // dim + "A1\'\'", //name + (complex double []) {1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1} // m + }, { 2, // dim "E\'\'", //name @@ -192,26 +486,6 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D3h = { 0.0, -0.9999999999999996, } }, - { - 1, // dim - "A2\'", //name - (complex double []) {1, 1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1} // m - }, - { - 1, // dim - "A1\'\'", //name - (complex double []) {1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1} // m - }, - { - 1, // dim - "A1\'", //name - (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} // m - }, - { - 1, // dim - "A2\'\'", //name - (complex double []) {1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1} // m - }, } // end of irreps }; @@ -284,8 +558,13 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D4h = { (struct qpms_finite_group_irrep_t[]) { // irreps { 1, // dim - "B1\'\'", //name - (complex double []) {1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1} // m + "A2\'", //name + (complex double []) {1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1} // m + }, + { + 1, // dim + "B1\'", //name + (complex double []) {1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1} // m }, { 2, // dim @@ -341,6 +620,36 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D4h = { 0.0, 1.0, } }, + { + 1, // dim + "A2\'\'", //name + (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1} // m + }, + { + 1, // dim + "B2\'", //name + (complex double []) {1, -1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1} // m + }, + { + 1, // dim + "A1\'", //name + (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} // m + }, + { + 1, // dim + "A1\'\'", //name + (complex double []) {1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1} // m + }, + { + 1, // dim + "B2\'\'", //name + (complex double []) {1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1} // m + }, + { + 1, // dim + "B1\'\'", //name + (complex double []) {1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1} // m + }, { 2, // dim "E\'\'", //name @@ -395,41 +704,6 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_D4h = { 0.0, -1.0, } }, - { - 1, // dim - "B1\'", //name - (complex double []) {1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1} // m - }, - { - 1, // dim - "A2\'\'", //name - (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1} // m - }, - { - 1, // dim - "B2\'", //name - (complex double []) {1, -1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1} // m - }, - { - 1, // dim - "A1\'\'", //name - (complex double []) {1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1} // m - }, - { - 1, // dim - "A2\'", //name - (complex double []) {1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1} // m - }, - { - 1, // dim - "B2\'\'", //name - (complex double []) {1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1} // m - }, - { - 1, // dim - "A1\'", //name - (complex double []) {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} // m - }, } // end of irreps }; @@ -463,184 +737,110 @@ const qpms_finite_group_t QPMS_FINITE_GROUP_trivial_g = { } // end of irreps }; -const qpms_finite_group_t QPMS_FINITE_GROUP_C4v = { - "C4v", // name - 8, // order +const qpms_finite_group_t QPMS_FINITE_GROUP_x_and_z_flip = { + "x_and_z_flip", // name + 4, // order 0, // idi (qpms_gmi_t[]) { // mt - 0, 1, 2, 3, 4, 5, 6, 7, - 1, 2, 3, 0, 7, 4, 5, 6, - 2, 3, 0, 1, 6, 7, 4, 5, - 3, 0, 1, 2, 5, 6, 7, 4, - 4, 5, 6, 7, 0, 1, 2, 3, - 5, 6, 7, 4, 3, 0, 1, 2, - 6, 7, 4, 5, 2, 3, 0, 1, - 7, 4, 5, 6, 1, 2, 3, 0, + 0, 1, 2, 3, + 1, 0, 3, 2, + 2, 3, 0, 1, + 3, 2, 1, 0, }, (qpms_gmi_t[]) { // invi - 0, 3, 2, 1, 4, 5, 6, 7 + 0, 1, 2, 3 }, - (qpms_gmi_t[]) {1, 7}, // gens + (qpms_gmi_t[]) {1, 3}, // gens 2, // ngens (qpms_permutation_t[]){ // permrep "(3)", - "(0 1 2 3)", - "(0 2)(1 3)", - "(0 3 2 1)", - "(3)(0 2)", - "(0 3)(1 2)", - "(1 3)", + "(3)(0 1)", "(0 1)(2 3)", + "(2 3)", }, NULL, // elemlabels 4, // permrep_nelem (qpms_irot3_t[]) { // rep3d {{1.0+0.0*I, 0.0+0.0*I}, 1}, - {{0.7071067811865476+0.7071067811865475*I, 0.0+0.0*I}, 1}, - {{2.220446049250313e-16+1.0*I, 0.0+0.0*I}, 1}, - {{-0.7071067811865474+0.7071067811865477*I, 0.0+0.0*I}, 1}, - {{0.0+0.0*I, 0.7071067811865477-0.7071067811865474*I}, -1}, - {{0.0+0.0*I, 1.0+2.220446049250313e-16*I}, -1}, - {{0.0+0.0*I, 0.7071067811865475+0.7071067811865476*I}, -1}, {{0.0+0.0*I, 0.0+1.0*I}, -1}, - }, - 5, // nirreps - (struct qpms_finite_group_irrep_t[]) { // irreps - { - 2, // dim - "E", //name - (complex double []) { - // (3) - 1.0, 0.0, - 0.0, 1.0, - // (0 1 2 3) - 0.0, -1.0, - 1.0, 0.0, - // (0 2)(1 3) - -1.0, 0.0, - 0.0, -1.0, - // (0 3 2 1) - 0.0, 1.0, - -1.0, 0.0, - // (3)(0 2) - 0.0, 1.0, - 1.0, 0.0, - // (0 3)(1 2) - 1.0, 0.0, - 0.0, -1.0, - // (1 3) - 0.0, -1.0, - -1.0, 0.0, - // (0 1)(2 3) - -1.0, 0.0, - 0.0, 1.0, - } - }, - { - 1, // dim - "A1", //name - (complex double []) {1, 1, 1, 1, 1, 1, 1, 1} // m - }, - { - 1, // dim - "B2", //name - (complex double []) {1, -1, 1, -1, 1, -1, 1, -1} // m - }, - { - 1, // dim - "B1", //name - (complex double []) {1, -1, 1, -1, -1, 1, -1, 1} // m - }, - { - 1, // dim - "A2", //name - (complex double []) {1, 1, 1, 1, -1, -1, -1, -1} // m - }, - } // end of irreps -}; - -const qpms_finite_group_t QPMS_FINITE_GROUP_D2h = { - "D2h", // name - 8, // order - 0, // idi - (qpms_gmi_t[]) { // mt - 0, 1, 2, 3, 4, 5, 6, 7, - 1, 0, 3, 2, 5, 4, 7, 6, - 2, 3, 0, 1, 6, 7, 4, 5, - 3, 2, 1, 0, 7, 6, 5, 4, - 4, 5, 6, 7, 0, 1, 2, 3, - 5, 4, 7, 6, 1, 0, 3, 2, - 6, 7, 4, 5, 2, 3, 0, 1, - 7, 6, 5, 4, 3, 2, 1, 0, - }, - (qpms_gmi_t[]) { // invi - 0, 1, 2, 3, 4, 5, 6, 7 - }, - (qpms_gmi_t[]) {1, 3, 7}, // gens - 3, // ngens - (qpms_permutation_t[]){ // permrep - "(5)", - "(5)(0 1)(2 3)", - "(5)(0 2)(1 3)", - "(5)(0 3)(1 2)", - "(0 3)(1 2)(4 5)", - "(0 2)(1 3)(4 5)", - "(0 1)(2 3)(4 5)", - "(4 5)", - }, - NULL, // elemlabels - 6, // permrep_nelem - (qpms_irot3_t[]) { // rep3d - {{1.0+0.0*I, 0.0+0.0*I}, 1}, - {{0.0+0.0*I, 0.0+1.0*I}, -1}, - {{0.0+1.0*I, 0.0+0.0*I}, 1}, - {{0.0+0.0*I, 1.0+0.0*I}, -1}, - {{0.0+0.0*I, 0.0+1.0*I}, 1}, - {{-1.0+0.0*I, 0.0+0.0*I}, -1}, {{0.0+0.0*I, -1.0+0.0*I}, 1}, {{0.0+1.0*I, 0.0+0.0*I}, -1}, }, - 8, // nirreps + 4, // nirreps (struct qpms_finite_group_irrep_t[]) { // irreps { 1, // dim - "B1\'\'", //name - (complex double []) {1, -1, -1, 1, -1, 1, 1, -1} // m + "P\'\'", //name + (complex double []) {1, -1, 1, -1} // m }, { 1, // dim - "A2\'\'", //name - (complex double []) {1, 1, 1, 1, -1, -1, -1, -1} // m + "P\'", //name + (complex double []) {1, 1, 1, 1} // m }, { 1, // dim - "B1\'", //name - (complex double []) {1, 1, -1, -1, -1, -1, 1, 1} // m + "R\'\'", //name + (complex double []) {1, 1, -1, -1} // m }, { 1, // dim - "A2\'", //name - (complex double []) {1, -1, 1, -1, -1, 1, -1, 1} // m + "R\'", //name + (complex double []) {1, -1, -1, 1} // m + }, + } // end of irreps +}; + +const qpms_finite_group_t QPMS_FINITE_GROUP_y_and_z_flip = { + "y_and_z_flip", // name + 4, // order + 0, // idi + (qpms_gmi_t[]) { // mt + 0, 1, 2, 3, + 1, 0, 3, 2, + 2, 3, 0, 1, + 3, 2, 1, 0, + }, + (qpms_gmi_t[]) { // invi + 0, 1, 2, 3 + }, + (qpms_gmi_t[]) {1, 3}, // gens + 2, // ngens + (qpms_permutation_t[]){ // permrep + "(3)", + "(3)(0 1)", + "(0 1)(2 3)", + "(2 3)", + }, + NULL, // elemlabels + 4, // permrep_nelem + (qpms_irot3_t[]) { // rep3d + {{1.0+0.0*I, 0.0+0.0*I}, 1}, + {{0.0+0.0*I, 1.0+0.0*I}, -1}, + {{0.0+0.0*I, 0.0+1.0*I}, 1}, + {{0.0+1.0*I, 0.0+0.0*I}, -1}, + }, + 4, // nirreps + (struct qpms_finite_group_irrep_t[]) { // irreps + { + 1, // dim + "P\'\'", //name + (complex double []) {1, -1, 1, -1} // m }, { 1, // dim - "A1\'\'", //name - (complex double []) {1, -1, 1, -1, 1, -1, 1, -1} // m + "P\'", //name + (complex double []) {1, 1, 1, 1} // m }, { 1, // dim - "B2\'", //name - (complex double []) {1, -1, -1, 1, 1, -1, -1, 1} // m + "R\'\'", //name + (complex double []) {1, 1, -1, -1} // m }, { 1, // dim - "B2\'\'", //name - (complex double []) {1, 1, -1, -1, 1, 1, -1, -1} // m - }, - { - 1, // dim - "A1\'", //name - (complex double []) {1, 1, 1, 1, 1, 1, 1, 1} // m + "R\'", //name + (complex double []) {1, -1, -1, 1} // m }, } // end of irreps }; diff --git a/tests/staticgroups.h b/tests/staticgroups.h index aac2f18..b6ac2bb 100644 --- a/tests/staticgroups.h +++ b/tests/staticgroups.h @@ -1,10 +1,14 @@ #ifndef STATICGROUPS_H #define STATICGROUPS_H #include +extern const qpms_finite_group_t QPMS_FINITE_GROUP_C2; extern const qpms_finite_group_t QPMS_FINITE_GROUP_C2v; extern const qpms_finite_group_t QPMS_FINITE_GROUP_trivial_g; +extern const qpms_finite_group_t QPMS_FINITE_GROUP_C4; extern const qpms_finite_group_t QPMS_FINITE_GROUP_C4v; extern const qpms_finite_group_t QPMS_FINITE_GROUP_D3h; extern const qpms_finite_group_t QPMS_FINITE_GROUP_D2h; extern const qpms_finite_group_t QPMS_FINITE_GROUP_D4h; +extern const qpms_finite_group_t QPMS_FINITE_GROUP_x_and_z_flip; +extern const qpms_finite_group_t QPMS_FINITE_GROUP_y_and_z_flip; #endif