Skip to content

Commit

Permalink
pympb ports of of MPB's compute_symmetries (#1930)
Browse files Browse the repository at this point in the history
* pympb ports of of `transform_overlap`, `compute_symmetry`, & `compute_symmetries`

- corresponding MPB PRs: NanoComp/mpb#134 and NanoComp/mpb#100

* remove `!SCALAR_COMPLEX` and `HAVE_MPI` related code

* remove overlooked stale MPI-specific code
  • Loading branch information
thchr authored Feb 11, 2022
1 parent f359a04 commit ab5fe8c
Show file tree
Hide file tree
Showing 4 changed files with 254 additions and 13 deletions.
201 changes: 188 additions & 13 deletions libpympb/pympb.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,15 @@ void matrix3x3_to_arr(mpb_real arr[3][3], matrix3x3 m) {

cnumber cscalar2cnumber(scalar_complex cs) { return make_cnumber(CSCALAR_RE(cs), CSCALAR_IM(cs)); }

cvector3 cscalar32cvector3(const scalar_complex *cs)
{
cvector3 v;
v.x = cscalar2cnumber(cs[0]);
v.y = cscalar2cnumber(cs[1]);
v.z = cscalar2cnumber(cs[2]);
return v;
}

// Return a string describing the current parity, used for frequency and filename
// prefixes
const char *parity_string(maxwell_data *d) {
Expand Down Expand Up @@ -1701,14 +1710,16 @@ std::vector<mpb_real> mode_solver::get_output_k() {

mpb_real mode_solver::get_val(int ix, int iy, int iz, int nx, int ny, int nz, int last_dim_size,
mpb_real *data, int stride, int conjugate) {
// #ifndef SCALAR_COMPLEX
// CHECK(0, "get-*-point not yet implemented for mpbi!");
// #endif
// #ifdef HAVE_MPI
// CHECK(0, "get-*-point not yet implemented for MPI!");
// #else
// CHECK(0, "get-*-point not yet implemented for MPI!");
// #endif
(void)nx;
(void)last_dim_size;
(void)conjugate;
return data[(((ix * ny) + iy) * nz + iz) * stride];
// #endif
}

mpb_real mode_solver::interp_val(vector3 p, int nx, int ny, int nz, int last_dim_size,
Expand Down Expand Up @@ -1819,18 +1830,20 @@ mpb_real mode_solver::get_energy_point(vector3 p) {
return f_interp_val(p, mdata, (mpb_real *)curfield, 1, 0);
}

cvector3 mode_solver::get_bloch_field_point(vector3 p) {
scalar_complex field[3];
cvector3 F;

void mode_solver::get_bloch_field_point_(scalar_complex field[3], vector3 p) {
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]);
}

cvector3 mode_solver::get_bloch_field_point(vector3 p) {
scalar_complex field[3];
cvector3 F;

get_bloch_field_point_(field, p);
F = cscalar32cvector3(field);
return F;
}

Expand Down Expand Up @@ -1858,9 +1871,7 @@ cvector3 mode_solver::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]);
F = cscalar32cvector3(field);

return F;
}
Expand Down Expand Up @@ -2680,4 +2691,168 @@ bool with_hermitian_epsilon() {
return false;
#endif
}


/* --- port of mpb's mpb/transform.c --- */

/* 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` Python functions (_without_
* the Bloch included) or the `get_bfield` and `get_dfield` 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 mode_solver::transformed_overlap(matrix3x3 W, vector3 w)
{
int n1, n2, n3;
mpb_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)) {
meep::master_fprintf(stderr, "The D or B field must be loaded first.\n");
return integral;
}

/* Python interface of MPB does not run under mpbi or MPI */
// #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) {
meep::master_fprintf(stderr, "valid symmetry operations {W|w} must have |det(W)| = 1\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 == 0 ? 1e-20 : geometry_lattice.size.x);
kvector.y *= TWOPI/(geometry_lattice.size.y == 0 ? 1e-20 : geometry_lattice.size.y);
kvector.z *= TWOPI/(geometry_lattice.size.z == 0 ? 1e-20 : 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]);

/* include Bloch phases */
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));

/* add integrand-contribution to integral */
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 mode_solver::compute_symmetry(int which_band, matrix3x3 W, vector3 w) {
cnumber symval;
get_bfield(which_band); // _without_ Bloch phase
symval = transformed_overlap(W, w);

return symval;
}

} // namespace meep_mpb
4 changes: 4 additions & 0 deletions libpympb/pympb.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -175,6 +175,7 @@ struct mode_solver {
cmatrix3x3 get_epsilon_inverse_tensor_point(vector3 p);
mpb_real get_energy_point(vector3 p);
cvector3 get_field_point(vector3 p);
void get_bloch_field_point_(scalar_complex field[3], vector3 p);
cvector3 get_bloch_field_point(vector3 p);

void multiply_bloch_phase(std::complex<double> *cdata = NULL);
Expand All @@ -195,6 +196,9 @@ struct mode_solver {

vector3 get_dominant_planewave(int band);

cnumber transformed_overlap(matrix3x3 W, vector3 w);
cnumber compute_symmetry(int which_band, matrix3x3 W, vector3 w);

private:
int kpoint_index;
scalar_complex *curfield;
Expand Down
12 changes: 12 additions & 0 deletions python/solver.py
Original file line number Diff line number Diff line change
Expand Up @@ -1126,6 +1126,18 @@ def try_all_and_repeat(k):
def get_dominant_planewave(self, band):
return self.mode_solver.get_dominant_planewave(band)

def transformed_overlap(self, W, w):
return self.mode_solver.transformed_overlap(W, w)

def compute_symmetry(self, band, W, w):
return self.mode_solver.compute_symmetry(band, W, w)

def compute_symmetries(self, W, w):
symvals = []
for band in range(1, self.num_bands+1):
symvals.append(self.mode_solver.compute_symmetry(band, W, w))
return symvals


# Predefined output functions (functions of the band index), for passing to `run`

Expand Down
50 changes: 50 additions & 0 deletions python/tests/test_mpb.py
Original file line number Diff line number Diff line change
Expand Up @@ -1419,6 +1419,56 @@ def test_fractional_lattice(self):

self.check_band_range_data(expected_brd, ms.band_range_data)

def test_symmetry_overlap(self):
ms = mpb.ModeSolver(
num_bands=6,
k_points = [mp.Vector3(0.5, 0.5, 0.5),],
geometry = [mp.Sphere(0.25, material=mp.Medium(epsilon=13))],
geometry_lattice = mp.Lattice(size=mp.Vector3(1, 1, 1)),
resolution = 32,
filename_prefix = self.filename_prefix,
deterministic = True,
tolerance = 1e-12
)
ms.run()

# inversion symmetry
W = mp.Matrix([-1,0,0], [0,-1,0], [0,0,-1])
w = mp.Vector3(0,0,0)
symeigs = ms.compute_symmetries(W, w)
symeigs_expected = [1,1,1,-1,-1,-1]
for i in range(0,5):
self.assertAlmostEqual(symeigs[i].real, symeigs_expected[i], places=3)
self.assertAlmostEqual(symeigs[i].imag, 0, places=3)

# check that transformed_overlap agrees when called manually, with a
# "preloaded" bfield
for i in range(0,5):
ms.get_bfield(i+1, bloch_phase=False)
symeig_bfield = ms.transformed_overlap(W, w)
self.assertAlmostEqual(symeigs[i].real, symeig_bfield.real, places=3)
self.assertAlmostEqual(symeigs[i].imag, symeig_bfield.imag, places=3)

ms.get_dfield(i+1, bloch_phase=False)
symeig_dfield = ms.transformed_overlap(W, w)
self.assertAlmostEqual(symeigs[i].real, symeig_dfield.real, places=3)
self.assertAlmostEqual(symeigs[i].imag, symeig_dfield.imag, places=3)

# mirror symmetry: compare with run_zeven & run_zodd
ms.k_points = [mp.Vector3(0,0,0)] # must be at Γ cf. https://github.com/NanoComp/mpb/issues/114
Wz = mp.Matrix([1,0,0], [0,1,0], [0,0,-1])

ms.run_zeven()
symeigs_zeven = ms.compute_symmetries(Wz, w)
for i in range(0,5):
self.assertAlmostEqual(symeigs_zeven[i].real, 1, places=3)
self.assertAlmostEqual(symeigs_zeven[i].imag, 0, places=3)

ms.run_zodd()
symeigs_zodd = ms.compute_symmetries(Wz, w)
for i in range(0,5):
self.assertAlmostEqual(symeigs_zodd[i].real, -1, places=3)
self.assertAlmostEqual(symeigs_zodd[i].imag, 0, places=3)

if __name__ == '__main__':
unittest.main()

0 comments on commit ab5fe8c

Please sign in to comment.