From 74fb432d31905a43b2d2c1d9b88a9f4bbfd0a8b6 Mon Sep 17 00:00:00 2001 From: Thomas Christensen Date: Sun, 6 Jun 2021 14:23:16 -0400 Subject: [PATCH] Symmetry-transformed overlap integrals (#134) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * WIP: Initial commit for calculating symmetry transformed overlaps between Bloch states of H * Symmetry operations act on the full wave, not just the Bloch part; fix that * Minor comment improvements and an error message for detW ≠ 1 * generalize so that nontrivial mu can be used, and also allow passing in the d-field instead of the b-field * swap a few cnumbers for scalar_complex to hopefully avoid copying things unnecessarily * Fix bug in vector transformation and account for the fact that vector fields are in a Cartesian basis * Add band-index function interfaces to transformed_overlap without invoking get_bfield/get_dfield and add a clearer method description: - compute_symmetry (band function) - compute_symmetries (all bands) Also minor revisions to comments and micro-refactoring. * Fix erroneous operation of {W|w}^-1 on coordinate p - The translation component of {W|w}^-1 contains a factor of W^-1 - Previously, nonsymmorphic operations were consequently wrongly implemented * Remove MPI guards/warnings (#112 has been merged) * throw when used with mpbi - also add a comment to remind ourselves why the current implementation doesn't work with mpbi + a likely explanation for the origin of the issue * put guard-rails back on MPI: doesn't seem to work as is. - add a note describing the problem * improve consistency of comments * add Scheme documentation * comment nits * add tests and improve `get_bloch_field_point_` function parameter signature * Update mpb/fields.c Revert use of `static` for C89 compatibility. Co-authored-by: Steven G. Johnson Co-authored-by: Steven G. Johnson --- doc/docs/Scheme_User_Interface.md | 25 ++++ examples/check.ctl | 41 ++++++ mpb/Makefile.am | 2 +- mpb/fields.c | 31 ++--- mpb/mpb.h | 1 + mpb/mpb.scm.in | 7 + mpb/transform.c | 210 ++++++++++++++++++++++++++++++ src/matrices/scalar.h | 8 ++ 8 files changed, 307 insertions(+), 18 deletions(-) create mode 100644 mpb/transform.c diff --git a/doc/docs/Scheme_User_Interface.md b/doc/docs/Scheme_User_Interface.md index 10d44909..1b19fa5c 100644 --- a/doc/docs/Scheme_User_Interface.md +++ b/doc/docs/Scheme_User_Interface.md @@ -571,6 +571,31 @@ Returns a list of the group-velocity components (units of *c*) in the given *dir            As above, but returns the group velocity or component thereof only for band `which-band`. +### Band symmetry + +MPB can compute the "characters" of a symmetry operation $g$ applied to a field at a particular point, i.e. MPB can compute the overlap integrals: + +```math +(\mathbf{F}|g\mathbf{F}') = \int \mathbf{F}(\mathbf{r})^\dagger [g\mathbf{F}'(\mathbf{r})] \, \mathrm{d}\mathbf{r} = \int \mathbf{F}(\mathbf{r})^\dagger (g\mathbf{F}')(g^{-1}\mathbf{r}) \, \mathrm{d}\mathbf{r} +``` + +where **F** denotes either the **E**- or **H**-field and **F**′ denotes the associated **D**- or **B**-field (including Bloch phases) and integration is over the unit cell. +The space group operation $g$ may include both rotational ($W$) and translational ($w$) parts, and are specified by the associated matrix-column pair $g = \{W,w\}$ acting on points as $\{W,w\}\mathbf{r} = W\mathbf{r} + w$. +$W$ and $w$ should be supplied as `matrix3x3` and `vector3` variables in the basis of the lattice (listings of space group operations are available e.g. on the [Bilbao Crystallographic Server](https://www.cryst.ehu.es/)). +(Note that $g$ transforms both the vectorial components of **F**′ as well as its coordinates: if **F**′ refers to a **B**-field, an additional factor of $\mathop{\mathrm{det}}W$ is incorporated to account for its pseudovectorial nature.) + +Given the character tables associated with the little groups of a considered photonic crystal, this functionality can e.g. be used to infer which irreducible representation (irrep) a given band (or band grouping) transforms as across the Brillouin zone. + +**`(compute-symmetry which_band W w)`** +Compute the $(\mathbf{H}|g\mathbf{B})$ character for `which-band` under a symmetry operation with rotation `W` and translation `w`, returning a complex number. + +**`(compute-symmetries W w)`** +Equivalent to `compute-symmetry`, but returns characters for *all* computed bands, returning a list of complex numbers. + +**`(transformed_overlap W w)`** +Sets **F**′ to `curfield` (must be either a **D**- or **B**-field) and computes the associated character under the symmetry operation specified by `W` and `w`, returning a complex number. +Usually, it will be more convenient to call `compute-symmetry` (which wraps `transformed_overlap` with an initialization of `curfield` via `get-bfield`). + Field Manipulation ------------------ diff --git a/examples/check.ctl b/examples/check.ctl index b60948aa..4a6c789a 100644 --- a/examples/check.ctl +++ b/examples/check.ctl @@ -303,6 +303,47 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +(print + "**************************************************************************\n" + " Test case: symmetry transformed overlaps & inversion/mirror eigenvalues.\n" + "**************************************************************************\n" +) + +(set! resolution 16) +(set! num-bands 6) +(set! default-material air) +; define a simple geometry with inversion and z-mirror +(set! geometry-lattice + (make lattice (size 1 1 1) (basis1 1 0 0) (basis2 0 1 0) (basis3 0 0 1))) +(set! geometry + (list (make sphere (center 0 0 0) (radius 0.25) (material (make dielectric (epsilon 13) ))))) + +; compare inversion eigenvalues at k = [1,1,1]/2 against precomputed values +(set! k-points (list (vector3 0.5 0.5 0.5))) ; little group includes inversion +(define W (matrix3x3 (vector3 -1 0 0) (vector3 0 -1 0) (vector3 0 0 -1))) ; inversion as operation {W|w} +(define w (vector3 0 0 0)) + +(run) +(define symeigs-inv (compute-symmetries W w)) +(check-almost-equal (map real-part symeigs-inv) '(+1 +1 +1 -1 -1 -1)) +(check-almost-equal (map imag-part symeigs-inv) '( 0 0 0 0 0 0)) + +; compare with run-zeven and run-zodd at k = 0 [must be 0 due to https://github.com/NanoComp/mpb/issues/114] +(set! k-points (list (vector3 0 0 0))) +(define mz (matrix3x3 (vector3 1 0 0) (vector3 0 1 0) (vector3 0 0 -1))) + +(run-zeven) ; even z-parity check +(define symeigs-zeven (compute-symmetries mz w)) +(check-almost-equal (map real-part symeigs-zeven) '(+1 +1 +1 +1 +1 +1)) +(check-almost-equal (map imag-part symeigs-zeven) '( 0 0 0 0 0 0)) + +(run-zodd) ; odd z-parity check +(define symeigs-zodd (compute-symmetries mz w)) +(check-almost-equal (map real-part symeigs-zodd) '(-1 -1 -1 -1 -1 -1)) +(check-almost-equal (map imag-part symeigs-zodd) '( 0 0 0 0 0 0)) + +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + (display-eigensolver-stats) (print "Relative error ranged from " min-err " to " max-err ", with a mean of " (/ sum-err num-err) "\n") diff --git a/mpb/Makefile.am b/mpb/Makefile.am index 27cf00ea..b8d85c54 100644 --- a/mpb/Makefile.am +++ b/mpb/Makefile.am @@ -20,7 +20,7 @@ EXTRA_DIST = mpb.scm.in epsilon.c mu.c nodist_pkgdata_DATA = $(SPECIFICATION_FILE) -MY_SOURCES = medium.c epsilon_file.c field-smob.c fields.c \ +MY_SOURCES = transform.c medium.c epsilon_file.c field-smob.c fields.c \ material_grid.c material_grid_opt.c matrix-smob.c mpb.c field-smob.h matrix-smob.h mpb.h my-smob.h MY_LIBS = $(top_builddir)/src/matrixio/libmatrixio.a $(top_builddir)/src/libmpb@MPB_SUFFIX@.la $(NLOPT_LIB) -lctl $(GUILE_LIBS) diff --git a/mpb/fields.c b/mpb/fields.c index 72ff978d..c154f11c 100644 --- a/mpb/fields.c +++ b/mpb/fields.c @@ -706,20 +706,23 @@ number get_energy_point(vector3 p) return f_interp_val(p, mdata, (real *) curfield, 1, 0); } -cvector3 get_bloch_field_point(vector3 p) +void get_bloch_field_point_(scalar_complex field[3], vector3 p) { - scalar_complex field[3]; - cvector3 F; - CHECK(curfield && strchr("dhbecv", curfield_type), "field must be must be loaded before get-*field*-point"); field[0] = f_interp_cval(p, mdata, &curfield[0].re, 6); field[1] = f_interp_cval(p, mdata, &curfield[1].re, 6); field[2] = f_interp_cval(p, mdata, &curfield[2].re, 6); - F.x = cscalar2cnumber(field[0]); - F.y = cscalar2cnumber(field[1]); - F.z = cscalar2cnumber(field[2]); - return F; +} + +cvector3 get_bloch_field_point(vector3 p) +{ + scalar_complex field[3]; + cvector3 F; + + get_bloch_field_point_(field, p); + + return cscalar32cvector3(field); } cvector3 get_field_point(vector3 p) @@ -727,11 +730,7 @@ cvector3 get_field_point(vector3 p) scalar_complex field[3], phase; cvector3 F; - CHECK(curfield && strchr("dhbecv", curfield_type), - "field must be must be loaded before get-*field*-point"); - field[0] = f_interp_cval(p, mdata, &curfield[0].re, 6); - field[1] = f_interp_cval(p, mdata, &curfield[1].re, 6); - field[2] = f_interp_cval(p, mdata, &curfield[2].re, 6); + get_bloch_field_point_(field, p); if (curfield_type != 'v') { double phase_phi = TWOPI * @@ -744,10 +743,8 @@ cvector3 get_field_point(vector3 p) CASSIGN_MULT(field[2], field[2], phase); } - F.x = cscalar2cnumber(field[0]); - F.y = cscalar2cnumber(field[1]); - F.z = cscalar2cnumber(field[2]); - return F; + F = cscalar32cvector3(field); + return F; } cnumber get_bloch_cscalar_point(vector3 p) diff --git a/mpb/mpb.h b/mpb/mpb.h index ce47e88e..2e2187c8 100644 --- a/mpb/mpb.h +++ b/mpb/mpb.h @@ -74,6 +74,7 @@ extern int kpoint_index; extern void compute_field_squared(void); void get_efield(integer which_band); real mean_medium_from_matrix(const symmetric_matrix *eps_inv); +void get_bloch_field_point_(scalar_complex *field, vector3 p); /**************************************************************************/ diff --git a/mpb/mpb.scm.in b/mpb/mpb.scm.in index b4a951bd..364c45e0 100644 --- a/mpb/mpb.scm.in +++ b/mpb/mpb.scm.in @@ -211,6 +211,13 @@ (define-external-function compute-energy-in-object-list false false 'number (make-list-type 'geometric-object)) +(define-external-function transformed-overlap false false + 'cnumber 'matrix3x3 'vector3) +(define-external-function compute-symmetries false false + (make-list-type 'cnumber) 'matrix3x3 'vector3) +(define-external-function compute-symmetry false false + 'cnumber 'integer 'matrix3x3 'vector3) + (define-external-function output-field-to-file false false no-return-value 'integer 'string) diff --git a/mpb/transform.c b/mpb/transform.c new file mode 100644 index 00000000..a84f605d --- /dev/null +++ b/mpb/transform.c @@ -0,0 +1,210 @@ +/* Copyright (C) 1999-2014 Massachusetts Institute of Technology. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + + +#include +#include +#include + +#include "config.h" +#include +#include +#include + +#include "mpb.h" +#include +#include +#include + + + +/* If `curfield` is the real-space D-field (of band-index i), computes the overlap + * ∫ Eᵢ†(r){W|w}Dᵢ(r) dr = ∫ Eᵢ†(r)(WDᵢ)({W|w}⁻¹r) dr, + * for a symmetry operation {W|w} with rotation W and translation w; each specified + * in the lattice basis. The vector fields Eᵢ and Dᵢ include Bloch phases. + * If `curfield` is the real-space B-field, the overlap + * ∫ Hᵢ†(r){W|w}Bᵢ(r) dr = det(W) ∫ Hᵢ†(r)(WBᵢ)({W|w}⁻¹r) dr, + * is computed instead. Note that a factor det(W) is then included since B & H are + * pseudovectors. As a result, the computed symmetry expectation values are + * independent of whether the D- or B-field is used. + * No other choices for `curfield` are allowed: to set `curfield` to the real-space + * B- or D-field use the `get-bfield` and `get-dfield` Scheme functions or the + * `get_bfield` and `get_dfield` in C functions. + * Usually, it will be more convenient to use the wrapper `compute_symmetry(i, W, w)` + * which defaults to the B-field (since μ = 1 usually) and takes a band-index `i`. + */ +cnumber transformed_overlap(matrix3x3 W, vector3 w) +{ + int n1, n2, n3; + real s1, s2, s3, c1, c2, c3; + cnumber integral = {0,0}, integral_sum; + number detW; + vector3 kvector = cur_kvector; + matrix3x3 invW, Wc; + + if (!curfield || !strchr("db", curfield_type)) { + mpi_one_fprintf(stderr, "curfield must refer to either a B- or D-field\n"); + return integral; + } + + #ifndef SCALAR_COMPLEX + CHECK(0, "transformed_overlap(..) is not yet implemented for mpbi"); + /* NOTE: Not completely sure why the current implementation does not work for mpbi + * (i.e for assumed-inversion): one clue is that running this with mpbi and the + * error-check off gives symmetry eigenvalues whose norm are ≈(ni÷2+1)/ni (where + * ni=n1=n2=n3) instead of near-unity (as they should be). This suggests we are not + * counting the number of grid points correctly somehow: I think the root issue is + * that we use the LOOP_XYZ macro, which has special handling for mbpi (i.e. only + * "visits" the "inversion-reduced" part of the unitcell): but here, we actually + * really want to (or at least assume to) visit _all_ the points in the unitcell. + */ + #endif + #ifdef HAVE_MPI + CHECK(0, "transformed_overlap(..) is not yet implemented for MPI"); + /* NOTE: It seems there's some racey stuff going on with the MPI implementation, + * unfortunately, so it doesn't end giving consistent (or even meaningful) results. + * The issue could be that both `LOOP_XYZ` is distributed over workers _and_ that + * `get_bloch_field_point_` also is (via the interpolation). I'm imagining that such + * a naive "nested parallelism" doesn't jive here. + * Long story short: disable it for now. + */ + #endif + + /* prepare before looping ... */ + n1 = mdata->nx; n2 = mdata->ny; n3 = mdata->nz; + + s1 = geometry_lattice.size.x / n1; /* pixel spacings */ + s2 = geometry_lattice.size.y / n2; + s3 = geometry_lattice.size.z / n3; + c1 = n1 <= 1 ? 0 : geometry_lattice.size.x * 0.5; /* offsets (negative) */ + c2 = n2 <= 1 ? 0 : geometry_lattice.size.y * 0.5; + c3 = n3 <= 1 ? 0 : geometry_lattice.size.z * 0.5; + + detW = matrix3x3_determinant(W); + if (fabs(fabs(detW) - 1.0) > 1e-12) { + mpi_one_fprintf(stderr, "valid symmetry operations {W|w} must have |det(W)| = 1.0\n"); + return integral; + } + invW = matrix3x3_inverse(W); + /* W is specified in the lattice basis, but field *vectors* evaluated in real-space + * are in a Cartesian basis: when we transform the vector components, we must account + * for this difference. We transform W to a Cartesian basis Wc = RWR⁻¹ (with + * R = geometry_lattice.basis, a matrix w/ columns of Cartesian basis vectors) and + * then use Wc to transform vector fields. + */ + Wc = matrix3x3_mult(matrix3x3_mult(geometry_lattice.basis, W), + matrix3x3_inverse(geometry_lattice.basis)); + + /* hoist rescalings outside the loop (maybe licm takes care of it, but do it anyway) */ + kvector.x *= TWOPI/geometry_lattice.size.x; + kvector.y *= TWOPI/geometry_lattice.size.y; + kvector.z *= TWOPI/geometry_lattice.size.z; + + /* loop over coordinates (introduces int vars `i1`, `i2`, `i3`, `xyz_index`) */ + LOOP_XYZ(mdata) { /* implies two opening braces `{{` */ + vector3 p, pt; + scalar_complex F[3], Ft[3], Ftemp[3], integrand, phase; + double deltaphi; + + /* current lattice coordinate */ + p.x = i1 * s1 - c1; + p.y = i2 * s2 - c2; + p.z = i3 * s3 - c3; + + /* transformed coordinate pt = {W|w}⁻¹p = W⁻¹(p-w) since {W|w}⁻¹={W⁻¹|-W⁻¹w} */ + pt = matrix3x3_vector3_mult(invW, vector3_minus(p, w)); + + /* Bloch field value at transformed coordinate pt: interpolation is needed to + * ensure generality in the case of fractional translations. */ + get_bloch_field_point_(Ftemp, pt); /* assign `Ftemp` to field at `pt` (without eⁱᵏʳ factor) */ + + /* define `Ft` as the vector components of `Ftemp` transformed by `Wc`; we just + * write out the matrix-product manually here, for both real & imag parts */ + Ft[0].re = Wc.c0.x*Ftemp[0].re + Wc.c1.x*Ftemp[1].re + Wc.c2.x*Ftemp[2].re; + Ft[0].im = Wc.c0.x*Ftemp[0].im + Wc.c1.x*Ftemp[1].im + Wc.c2.x*Ftemp[2].im; + Ft[1].re = Wc.c0.y*Ftemp[0].re + Wc.c1.y*Ftemp[1].re + Wc.c2.y*Ftemp[2].re; + Ft[1].im = Wc.c0.y*Ftemp[0].im + Wc.c1.y*Ftemp[1].im + Wc.c2.y*Ftemp[2].im; + Ft[2].re = Wc.c0.z*Ftemp[0].re + Wc.c1.z*Ftemp[1].re + Wc.c2.z*Ftemp[2].re; + Ft[2].im = Wc.c0.z*Ftemp[0].im + Wc.c1.z*Ftemp[1].im + Wc.c2.z*Ftemp[2].im; + + /* get the Bloch field value at current point `p` (without eⁱᵏʳ factor). + * We multiply the input field `F` (either B or D-field) with μ⁻¹ or ε⁻¹ to get + * H- or E-fields, as the relevant overlap is ⟨F|Ft⟩ = ⟨H|Bt⟩ or ⟨E|Dt⟩, with + * t-postscript denoting a field transformed by {W|w}. Here, we essentially + * adapt some boiler-plate code from compute_field_energy_internal in fields.c */ + if (curfield_type == 'd') { + assign_symmatrix_vector(F, mdata->eps_inv[xyz_index], curfield+3*xyz_index); + } + else if (curfield_type == 'b' && mdata->mu_inv != NULL) { + assign_symmatrix_vector(F, mdata->mu_inv[xyz_index], curfield+3*xyz_index); + } + else { + F[0] = curfield[3*xyz_index]; + F[1] = curfield[3*xyz_index+1]; + F[2] = curfield[3*xyz_index+2]; + } + + /* inner product of F and Ft={W|w}F in Bloch form */ + CASSIGN_CONJ_MULT(integrand, F[0], Ft[0]); /* add adjoint(F)*Ft to integrand */ + CACCUMULATE_SUM_CONJ_MULT(integrand, F[1], Ft[1]); + CACCUMULATE_SUM_CONJ_MULT(integrand, F[2], Ft[2]); + + /* so far, Bloch phases have been excluded; we include them here. It saves two + * trigonometric operations if we do them just jointly for `p` and `pt`. */ + deltaphi = (kvector.x*(pt.x-p.x) + kvector.y*(pt.y-p.y) + kvector.z*(pt.z-p.z)); + CASSIGN_SCALAR(phase, cos(deltaphi), sin(deltaphi)); + + /* finally, add integrand-contribution to integral, with Bloch phases included */ + integral.re += CSCALAR_MULT_RE(integrand, phase); + integral.im += CSCALAR_MULT_IM(integrand, phase); + }}} + + integral.re *= Vol / H.N; + integral.im *= Vol / H.N; + + mpi_allreduce(&integral, &integral_sum, 2, number, + MPI_DOUBLE, MPI_SUM, mpb_comm); + + if (curfield_type == 'b') { /* H & B are pseudovectors => transform includes det(W) */ + integral_sum.re *= detW; + integral_sum.im *= detW; + } + + return integral_sum; +} + +cnumber compute_symmetry(int which_band, matrix3x3 W, vector3 w) +{ + cnumber symval; + get_bfield(which_band); + symval = transformed_overlap(W, w); + + return symval; +} + +cnumber_list compute_symmetries(matrix3x3 W, vector3 w) +{ + cnumber_list symvals; + int ib; + symvals.num_items = num_bands; + CHK_MALLOC(symvals.items, cnumber, num_bands); + + for (ib = 0; ib < num_bands; ++ib){ + symvals.items[ib] = compute_symmetry(ib+1, W, w); + } + return symvals; +} \ No newline at end of file diff --git a/src/matrices/scalar.h b/src/matrices/scalar.h index dbe21d06..9749453a 100644 --- a/src/matrices/scalar.h +++ b/src/matrices/scalar.h @@ -63,6 +63,14 @@ typedef struct { bbbb_re * cccc_im + bbbb_im * cccc_re); \ } +/* a = conj(b) * c */ +#define CASSIGN_CONJ_MULT(a, b, c) { \ + real bbbb_re = (b).re, bbbb_im = (b).im; \ + real cccc_re = (c).re, cccc_im = (c).im; \ + CASSIGN_SCALAR(a, bbbb_re * cccc_re + bbbb_im * cccc_im, \ + bbbb_re * cccc_im - bbbb_im * cccc_re); \ +} + /* a = b / c = b * conj(c) / |c|^2 */ #define CASSIGN_DIV(a, b, c) { \ scalar_complex aaaa_tmp; real aaaa_tmp_norm; \