From d63dd2cae2f29cef2ef52d761c0e656e9fec1ad9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Mar 2019 07:01:10 +0000 Subject: [PATCH 01/16] Fix projections - work in progress Former-commit-id: f2f6d3d762fc7c2d728af0cf31c1229ed8db1b1c --- qpms/qpms_c.pyx | 2 +- qpms/tmatrix_io.c | 4 +- setup.py | 5 +- tests/fixproj-notes | 4 ++ tests/quaternion_O3.c | 8 +++- tests/scatsys.py | 2 + tests/ss_syms_packs.c | 13 ++++-- tests/sss1.c | 105 ++++++++++++++++++++++++++++++++++++++++++ tests/sss2.c | 105 ++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 239 insertions(+), 9 deletions(-) create mode 100644 tests/fixproj-notes create mode 100644 tests/sss1.c create mode 100644 tests/sss2.c 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/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/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/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..0e9b458 --- /dev/null +++ b/tests/sss2.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 = {0, 0, 1}; + 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; +} From b97f3e6e2eeb2296f8dbe366ee75f1603cb3b857 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Mar 2019 15:07:04 +0200 Subject: [PATCH 02/16] Dump the whole translation matrix... The translation matrix has non-zero diagonal elements (which is obviously wrong). Former-commit-id: 91fbec3412af62799eec2784ddd5d4e43c30d7a4 --- tests/sss2.c | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/tests/sss2.c b/tests/sss2.c index 0e9b458..ac4b258 100644 --- a/tests/sss2.c +++ b/tests/sss2.c @@ -23,7 +23,7 @@ double uniform_random(double min, double max) { int main() { srand(666); -#if 0 +#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); @@ -55,7 +55,7 @@ int main() t2->m[i + i*b2->n] = 1; #endif - const cart3_t pp1 = {0, 0, 1}; + const cart3_t pp1 = {1, 2, 0}; qpms_tmatrix_t * tmlist[] = {t1, t2}; qpms_particle_tid_t plist[] = {{pp1, 0}}; @@ -75,6 +75,39 @@ int main() 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, From 8efd582daf522b45dc42f43cb582a4fc9963ee1c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Mar 2019 15:43:28 +0200 Subject: [PATCH 03/16] Fix offset incrementation on the diagonal. Solves the translation matrix projection discrepancy. Former-commit-id: 7bb7dbba47aa5ac3df32177098c096a86d2879fd --- qpms/scatsystem.c | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 07f70a4..7a66196 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1120,13 +1120,14 @@ complex double *qpms_scatsys_build_translation_matrix_full( const cart3_t posR = ss->p[piR].pos; size_t fullvec_offsetC = 0; for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { - if(piC == piR) continue; // The diagonal will be dealt with later. const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec; - const cart3_t posC = ss->p[piC].pos; - QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, - target + fullvec_offsetR*full_len + fullvec_offsetC, - bspecR, full_len, bspecC, 1, - k, posR, posC)); + if(piC != piR) { // The diagonal will be dealt with later. + const cart3_t posC = ss->p[piC].pos; + QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, + target + fullvec_offsetR*full_len + fullvec_offsetC, + bspecR, full_len, bspecC, 1, + k, posR, posC)); + } fullvec_offsetC += bspecC->n; } assert(fullvec_offsetC = full_len); @@ -1162,19 +1163,20 @@ complex double *qpms_scatsys_build_modeproblem_matrix_full( // dest particle T-matrix const complex double *tmmR = ss->tm[ss->p[piR].tmatrix_id]->m; for(qpms_ss_pi_t piC = 0; piC < ss->p_count; ++piC) { - if(piC == piR) continue; // The diagonal will be dealt with later. const qpms_vswf_set_spec_t *bspecC = ss->tm[ss->p[piC].tmatrix_id]->spec; - const cart3_t posC = ss->p[piC].pos; - QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, - tmp, // tmp is S(piR<-piC) - bspecR, bspecC->n, bspecC, 1, - k, posR, posC)); - cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, + if(piC != piR) { // The diagonal will be dealt with later. + const cart3_t posC = ss->p[piC].pos; + QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(ss->c, + tmp, // tmp is S(piR<-piC) + bspecR, bspecC->n, bspecC, 1, + k, posR, posC)); + cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, &one/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/, tmp/*b*/, bspecC->n/*ldb*/, &zero/*beta*/, target + fullvec_offsetR*full_len + fullvec_offsetC /*c*/, full_len /*ldc*/); + } fullvec_offsetC += bspecC->n; } fullvec_offsetR += bspecR->n; From 9a592b4682329788023eb5a86b84d21ca55f8c1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Mar 2019 17:17:49 +0200 Subject: [PATCH 04/16] Test with more particles; Apparetly, the result are now correct for particles lying strictly on the z-axis. particles strictly Former-commit-id: 56c968d3415a01519437479c6f4524a0cd10dc2f --- tests/sss2.c | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/tests/sss2.c b/tests/sss2.c index ac4b258..37666e1 100644 --- a/tests/sss2.c +++ b/tests/sss2.c @@ -55,17 +55,19 @@ int main() t2->m[i + i*b2->n] = 1; #endif - const cart3_t pp1 = {1, 2, 0}; + const cart3_t pp1 = {0, 0, 1}; + const cart3_t pp2 = {0, 0, 1.4}; + const cart3_t pp3 = {0, 0, 0}; qpms_tmatrix_t * tmlist[] = {t1, t2}; - qpms_particle_tid_t plist[] = {{pp1, 0}}; + qpms_particle_tid_t plist[] = {{pp1, 0}, {pp2, 1}, {pp3, 1}}; qpms_scatsys_t protoss; protoss.tm = tmlist; - protoss.tm_count=2; + protoss.tm_count=sizeof(tmlist)/sizeof(qpms_tmatrix_t *); protoss.p = plist; - protoss.p_count=1; + protoss.p_count=sizeof(plist)/sizeof(qpms_particle_tid_t); - qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, D2h); + qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, C4v); 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, @@ -92,10 +94,10 @@ int main() 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; + fullvec_offset_src += bspecn_src; } + fullvec_offset_dest += bspecn_dest; } } { From a2e61ad67adb9d6e39e0ac0b1b0f3b350d9ace93 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Mar 2019 18:32:38 +0200 Subject: [PATCH 05/16] Particle position dump Former-commit-id: 75cfdfff52531da4d9e52a162774fe59d2595a1b --- qpms/scatsystem.c | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 7a66196..341b838 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 From 6f66ddc8450e737f2525cbf85198007270ab0904 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Mar 2019 21:22:31 +0200 Subject: [PATCH 06/16] Additional symmetries Former-commit-id: dce9d43c0bfbce7166e6365f40591fbf5f9b4873 --- qpms/symmetries.py | 75 ++++- tests/generate_groupcode.py | 3 +- tests/staticgroups.c | 596 ++++++++++++++++++++++++------------ tests/staticgroups.h | 4 + 4 files changed, 478 insertions(+), 200 deletions(-) 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/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/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 From 208e05d8112f300a9a6651a8c678cd64d4e445c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Mar 2019 21:47:01 +0200 Subject: [PATCH 07/16] Macro-driven sss2.c test Former-commit-id: 493a2c0c5a1149f57752241077d57985cfe9ca69 --- tests/sss2.c | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/tests/sss2.c b/tests/sss2.c index 37666e1..2791d80 100644 --- a/tests/sss2.c +++ b/tests/sss2.c @@ -1,4 +1,4 @@ -// 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 +// 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 @@ -12,8 +12,16 @@ 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; @@ -23,7 +31,7 @@ double uniform_random(double min, double max) { int main() { srand(666); -#if 1 +#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); @@ -55,8 +63,16 @@ int main() t2->m[i + i*b2->n] = 1; #endif - const cart3_t pp1 = {0, 0, 1}; +#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}; +#else + const cart3_t pp1 = {0, 0, 1.1}; const cart3_t pp2 = {0, 0, 1.4}; +#endif const cart3_t pp3 = {0, 0, 0}; qpms_tmatrix_t * tmlist[] = {t1, t2}; qpms_particle_tid_t plist[] = {{pp1, 0}, {pp2, 1}, {pp3, 1}}; @@ -67,7 +83,7 @@ int main() protoss.p = plist; protoss.p_count=sizeof(plist)/sizeof(qpms_particle_tid_t); - qpms_scatsys_t *ss = qpms_scatsys_apply_symmetry(&protoss, C4v); + 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, From 02ec8804d0834dd3cbee827dabcb74d81a3feb43 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Mar 2019 22:22:02 +0200 Subject: [PATCH 08/16] Fix quaternion conjugation order on vector rotation. Now it does not crash on C4v and the projectors seem to be indeed projectors, but the translation operator decomposition is still wrong. Former-commit-id: 6273c5a76b4e424f77b12514df209b995999b24b --- qpms/wigner.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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. From 8cba5396842c9b83ce119647307ac2b7f9c8a15e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Mar 2019 23:03:26 +0200 Subject: [PATCH 09/16] =?UTF-8?q?Tenhle=20p=C5=99=C3=ADpad=20rozb=C3=ADj?= =?UTF-8?q?=C3=AD=20rozborku=20sborku=20na=20diagon=C3=A1le.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit (Přitom prohození první a třetí částice v zadání dává správný výsledek.) Former-commit-id: 0da6fd465c44e803d15164a98dfd6896ff97d363 --- tests/sss3.c | 180 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 180 insertions(+) create mode 100644 tests/sss3.c diff --git a/tests/sss3.c b/tests/sss3.c new file mode 100644 index 0000000..e2120a9 --- /dev/null +++ b/tests/sss3.c @@ -0,0 +1,180 @@ +// 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); + + 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; +} From 6f219b2c847e24c2e35c4b5db09feaa821bcff7a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 13 Mar 2019 08:57:23 +0000 Subject: [PATCH 10/16] Dump orbit projection matrices in sss3.c Former-commit-id: 57814624995124945704e8a1bfc52273845633af --- tests/sss3.c | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/tests/sss3.c b/tests/sss3.c index e2120a9..23cedcd 100644 --- a/tests/sss3.c +++ b/tests/sss3.c @@ -93,6 +93,30 @@ int main() (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( From aee4d2db2a7d5a81427badfb801f7645237e900e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 13 Mar 2019 11:24:09 +0200 Subject: [PATCH 11/16] sss3 dump packed matrices Former-commit-id: 5b2f8a4ebe9829d63546ae91e24c25ec5966fe72 --- tests/sss3.c | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/tests/sss3.c b/tests/sss3.c index 23cedcd..deac10a 100644 --- a/tests/sss3.c +++ b/tests/sss3.c @@ -155,9 +155,18 @@ int main() } complex double *S_packed[ss->sym->nirreps]; - for (qpms_iri_t iri = 0; iri < ss->sym->nirreps; ++iri) + 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_recfull = qpms_scatsys_irrep_unpack_matrix(NULL, S_packed[0], ss, 0, false); From cc61a64313060824e64d41a3992f4b6c59191553 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 13 Mar 2019 11:36:52 +0200 Subject: [PATCH 12/16] Dump parts of the full matrix. Former-commit-id: 53b84dc031fab778b089226be8b26d89e6f6a4a0 --- tests/sss3.c | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/tests/sss3.c b/tests/sss3.c index deac10a..9f07128 100644 --- a/tests/sss3.c +++ b/tests/sss3.c @@ -167,6 +167,29 @@ int main() 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); From d11db980f98bab10c470e02b0fff1d781d332a76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 13 Mar 2019 12:28:40 +0200 Subject: [PATCH 13/16] A slow but working reference implementation of matrix (un)packing. Former-commit-id: 5f85c4983a8754ebae03980242b896a79954a73e --- qpms/qpms_error.h | 4 +- qpms/scatsystem.c | 114 ++++++++++++++++++++++++++++++++++++++++++++++ qpms/scatsystem.h | 20 ++++++++ 3 files changed, 136 insertions(+), 2 deletions(-) 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 341b838..600dada 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -866,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){ @@ -947,6 +1060,7 @@ complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed, 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, From 07aed6feb725f43858c456fe7e35a915c1bf3bb4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 13 Mar 2019 13:24:44 +0200 Subject: [PATCH 14/16] Fix offset increment in matrix packing. Now it works for a z-dipole system. Former-commit-id: 3e4aa9d4736ecc18e561d0ba5c2c19b0a027e37c --- qpms/scatsystem.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 600dada..87a5994 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1053,8 +1053,8 @@ 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; From ddfd4954dca9090311adfde54614087bb1cba9ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 13 Mar 2019 15:59:05 +0200 Subject: [PATCH 15/16] Fix complex conjugation; now all available point groups seem to work. Former-commit-id: 764155cddfd30162b01c6552cea8718d55a26fd6 --- qpms/scatsystem.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 87a5994..4bf02ef 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -730,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", From 12da9f3b5a7701f271817f5bae79ed5552d685a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 13 Mar 2019 16:18:47 +0200 Subject: [PATCH 16/16] Various tests and macro conditions. Former-commit-id: ed964c6284d98ed315a442be48b2dab3d7685222 --- qpms/translations.c | 24 +++++ tests/sss2.c | 34 ++++++- tests/sss4.c | 237 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 290 insertions(+), 5 deletions(-) create mode 100644 tests/sss4.c 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/tests/sss2.c b/tests/sss2.c index 2791d80..54c959d 100644 --- a/tests/sss2.c +++ b/tests/sss2.c @@ -41,10 +41,10 @@ int main() *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) + 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 = -l; m <= l; ++m) + 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); @@ -69,13 +69,17 @@ int main() #elif defined XLINE const cart3_t pp1 = {1.1, 0, 0}; const cart3_t pp2 = {1.4, 0, 0}; -#else +#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, 0}; + const cart3_t pp3 = {0, 0, 1}; qpms_tmatrix_t * tmlist[] = {t1, t2}; - qpms_particle_tid_t plist[] = {{pp1, 0}, {pp2, 1}, {pp3, 1}}; + qpms_particle_tid_t plist[] = {{pp1,1}, {pp2, 0}, {pp3, 1}, + }; qpms_scatsys_t protoss; protoss.tm = tmlist; @@ -136,6 +140,15 @@ int main() 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) { @@ -143,8 +156,19 @@ int main() 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); 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; +}