xy-periodic lattice Beyn algorithm support in ScatteringSystem

Gives same results as newbeyn_unitcell 26d6e969, making it obsolete.


Former-commit-id: b1b1b1e2c11f60948efda237388bfdf9b6d689f5

misc/rectlat_simple_modes.py

@@ -0,0 +1,158 @@
+#!/usr/bin/env python3
+
+import math
+from qpms.argproc import ArgParser
+
+ap = ArgParser(['rectlattice2d', 'single_particle', 'single_lMax'])
+ap.add_argument("-k", nargs=2, type=float, required=True, help='k vector', metavar=('K_X', 'K_Y'))
+ap.add_argument("--kpi", action='store_true', help="Indicates that the k vector is given in natural units instead of SI, i.e. the arguments given by -k shall be automatically multiplied by pi / period (given by -p argument)")
+ap.add_argument("--rank-tol", type=float, required=False)
+ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)')
+ap.add_argument("-t", "--rank-tolerance", type=float, default=1e11)
+ap.add_argument("-c", "--min-candidates", type=int, default=1, help='always try at least this many eigenvalue candidates, even if their SVs in the rank tests are lower than rank_tolerance')
+ap.add_argument("-T", "--residual-tolerance", type=float, default=2.)
+ap.add_argument("-b", "--band-index", type=int, required=True, help="Argument's absolute value determines the empty lattice band order (counted from 1), -/+ determines whether the eigenvalues are searched below/above that empty lattice band.")
+ap.add_argument("-f", "--interval-factor", type=float, default=0.1)
+ap.add_argument("-N", type=int, default="150", help="Integration contour discretisation size")
+ap.add_argument("-i", "--imaginary-aspect-ratio",  type=float, default=1, help="Aspect ratio of the integration contour (Im/Re)")
+ap.add_argument("-O", "--plot-out", type=str, required=False, help="path to plot output (optional)")
+ap.add_argument("-P", "--plot", action='store_true', help="if -p not given, plot to a default path")
+
+a=ap.parse_args()
+
+px, py = a.period
+
+if a.kpi:
+    a.k[0] *= math.pi/px
+    a.k[1] *= math.pi/py
+
+particlestr = ("sph" if a.height is None else "cyl") + ("_r%gnm" % (a.radius*1e9))
+if a.height is not None: particlestr += "_h%gnm" % (a.height * 1e9)
+defaultprefix = "%s_p%gnmx%gnm_m%s_n%g_b%+d_k(%g_%g)um-1_L%d_cn%d" % (
+    particlestr, px*1e9, py*1e9, str(a.material), a.refractive_index, a.band_index, a.k[0]*1e-6, a.k[1]*1e-6, a.lMax, a.N)
+
+import logging
+logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO)
+
+import numpy as np
+import qpms
+from qpms.cybspec import BaseSpec
+from qpms.cytmatrices import CTMatrix
+from qpms.qpms_c import Particle, ScatteringSystem, empty_lattice_modes_xy
+from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude
+from qpms.constants import eV, hbar
+eh = eV/hbar
+
+def inside_ellipse(point_xy, centre_xy, halfaxes_xy):
+    x = point_xy[0] - centre_xy[0]
+    y = point_xy[1] - centre_xy[1]
+    ax = halfaxes_xy[0]
+    ay = halfaxes_xy[1]
+    return ((x/ax)**2 + (y/ay)**2) <= 1
+
+a1 = ap.direct_basis[0]
+a2 = ap.direct_basis[1]
+
+#Particle positions
+orig_x = [0]                                 
+orig_y = [0]            
+orig_xy = np.stack(np.meshgrid(orig_x,orig_y),axis=-1)
+
+if a.material in lorentz_drude:
+	emg = EpsMuGenerator(lorentz_drude[a.material])
+else: # constant refractive index
+        emg = EpsMuGenerator(EpsMu(a.material**2))
+
+beta = np.array(a.k)
+
+empty_freqs = empty_lattice_modes_xy(ap.background_epsmu, ap.reciprocal_basis2pi, np.array([0,0]), 1)
+empty_freqs = empty_lattice_modes_xy(ap.background_epsmu, ap.reciprocal_basis2pi, beta, (1+abs(a.band_index)) * empty_freqs[1])
+
+# make the frequencies in the list unique
+empty_freqs = list(empty_freqs)
+i = 0
+while i < len(empty_freqs) - 1:
+    if math.isclose(empty_freqs[i], empty_freqs[i+1], rel_tol=1e-13):
+        del empty_freqs[i+1]
+    else:
+        i += 1
+
+logging.info("Empty freqs: %s", str(empty_freqs))
+if a.band_index > 0:
+    top = empty_freqs[a.band_index]
+    bottom = empty_freqs[a.band_index - 1]
+    lebeta_om = bottom
+else: # a.band_index < 0
+    top = empty_freqs[abs(a.band_index) - 1]
+    bottom = empty_freqs[abs(a.band_index) - 2] if abs(a.band_index) > 1 else 0.
+    lebeta_om = top
+#print(top,bottom,lebeta_om)
+freqradius = .5 * (top - bottom) * a.interval_factor
+
+centfreq = bottom + freqradius if a.band_index > 0 else top - freqradius
+
+bspec = BaseSpec(lMax = a.lMax)
+pp = Particle(orig_xy[0][0], t = ap.tmgen, bspec=bspec)
+
+ss, ssw = ScatteringSystem.create([pp], ap.background_emg, centfreq, latticebasis = ap.direct_basis)
+
+if freqradius == 0:
+    raise ValueError("Integration contour radius is set to zero. Are you trying to look below the lowest empty lattice band at the gamma point?")
+
+freqradius *= (1-1e-13) # to not totally touch the singularities
+
+with qpms.pgsl_ignore_error(15):
+    res = ss.find_modes(centfreq, freqradius, freqradius * a.imaginary_aspect_ratio,
+            blochvector = a.k, contour_points = a.N, rank_tol = a.rank_tolerance,
+            res_tol = a.residual_tolerance, rank_min_sel = a.min_candidates)
+
+logging.info("Eigenfrequencies found: %s" % str(res['eigval']))
+
+res['inside_contour'] = inside_ellipse((res['eigval'].real, res['eigval'].imag),
+         (centfreq.real, centfreq.imag), (freqradius, freqradius * a.imaginary_aspect_ratio))
+
+res['refractive_index_internal'] = [emg(om).n for om in res['eigval']]
+
+#del res['omega']  If contour points are not needed...
+#del res['ImTW'] # not if dbg=false anyway
+outfile = defaultprefix + ".npz" if a.output is None else a.output
+np.savez(outfile, meta=vars(a), empty_freqs=np.array(empty_freqs), **res)
+logging.info("Saved to %s" % outfile)
+
+
+if a.plot or (a.plot_out is not None):
+    imcut = np.linspace(0, -freqradius)
+    recut1 = np.sqrt(lebeta_om**2+imcut**2) # incomplete Gamma-related cut
+    recut2 = np.sqrt((lebeta_om/2)**2-imcut**2) + lebeta_om/2 # odd-power-lilgamma-related cut
+
+    import matplotlib
+    matplotlib.use('pdf')
+    from matplotlib import pyplot as plt
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111,)
+    #ax.plot(res['omega'].real/eh, res['omega'].imag/eh*1e3, ':') #res['omega'] not implemented in ScatteringSystem
+    ax.add_artist(matplotlib.patches.Ellipse((centfreq.real/eh, centfreq.imag/eh*1e3),
+        2*freqradius/eh, 2*freqradius*a.imaginary_aspect_ratio/eh*1e3, fill=False,
+        ls=':'))
+    ax.scatter(x=res['eigval'].real/eh, y=res['eigval'].imag/eh*1e3 , c = res['inside_contour']
+                      )
+    ax.plot(recut1/eh, imcut/eh*1e3)
+    ax.plot(recut2/eh, imcut/eh*1e3)
+    for i,om in enumerate(res['eigval']):
+            ax.annotate(str(i), (om.real/eh, om.imag/eh*1e3))
+    xmin = np.amin(res['eigval'].real)/eh
+    xmax = np.amax(res['eigval'].real)/eh
+    xspan = xmax-xmin
+    ymin = np.amin(res['eigval'].imag)/eh*1e3
+    ymax = np.amax(res['eigval'].imag)/eh*1e3
+    yspan = ymax-ymin
+    ax.set_xlim([xmin-.1*xspan, xmax+.1*xspan])
+    ax.set_ylim([ymin-.1*yspan, ymax+.1*yspan])
+    ax.set_xlabel('$\hbar \Re \omega / \mathrm{eV}$')
+    ax.set_ylabel('$\hbar \Im \omega / \mathrm{meV}$')
+    plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
+    fig.savefig(plotfile)
+
+exit(0)
+

qpms/argproc.py

@@ -124,9 +124,9 @@ class ArgParser:
         if l != 2: raise ValueError('Two basis vectors must be specified (have %d)' % l)
         from .qpms_c import lll_reduce
         self.direct_basis = lll_reduce(self.args.basis_vector, delta=1.)
-        # import numpy as np
+        import numpy as np
-        # self.reciprocal_basis1 = np.linalg.inv(self.direct_basis)
+        self.reciprocal_basis1 = np.linalg.inv(self.direct_basis)
-        # self.reciprocal_basis2pi = 2 * np.pi * self.reciprocal_basis1
+        self.reciprocal_basis2pi = 2 * np.pi * self.reciprocal_basis1
     
     def _eval_rectlattice2d(self): # feature: rectlattice2d
         a = self.args
@@ -141,4 +141,6 @@ class ArgParser:
         import numpy as np
         a.basis_vectors = [(a.period[0], 0.), (0., a.period[1])]
         self.direct_basis = np.array(a.basis_vectors)
+        self.reciprocal_basis1 = np.linalg.inv(self.direct_basis)
+        self.reciprocal_basis2pi = 2 * np.pi * self.reciprocal_basis1
 

qpms/qpms_c.pyx

@@ -13,7 +13,7 @@ from .cybspec cimport BaseSpec
 from .cycommon cimport make_c_string
 from .cycommon import string_c2py, PointGroupClass
 from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator, TMatrixInterpolator
-from .cymaterials cimport EpsMuGenerator
+from .cymaterials cimport EpsMuGenerator, EpsMu
 from libc.stdlib cimport malloc, free, calloc
 import warnings
 
@@ -647,16 +647,27 @@ cdef class ScatteringSystem:
         return target_np
 
     def find_modes(self, cdouble omega_centre, double omega_rr, double omega_ri, iri = None,
+            blochvector = None,
             size_t contour_points = 20, double rank_tol = 1e-4, size_t rank_min_sel=1,
             double res_tol = 0):
         """
         Attempts to find the eigenvalues and eigenvectors using Beyn's algorithm.
         """
 
-        cdef beyn_result_t *res = qpms_scatsys_finite_find_eigenmodes(self.s, 
+        cdef beyn_result_t *res
-                self.iri_py2c(iri),
+        cdef double blochvector_c[3]
+        if self.lattice_dimension == 0:
+            assert(blochvector is None)
+            res = qpms_scatsys_finite_find_eigenmodes(self.s, self.iri_py2c(iri),
                 omega_centre, omega_rr, omega_ri, contour_points, 
                 rank_tol, rank_min_sel, res_tol)
+        else:
+            if iri is not None: raise NotImplementedError("Irreps decomposition not yet supported for periodic systems")
+            blochvector_c[0] = blochvector[0]
+            blochvector_c[1] = blochvector[1]
+            blochvector_c[2] = blochvector[2] if len(blochvector) > 2 else 0
+            res = qpms_scatsys_periodic_find_eigenmodes(self.s, blochvector_c,
+                    omega_centre, omega_rr, omega_ri, contour_points, rank_tol, rank_min_sel, res_tol)
         if res == NULL: raise RuntimeError
 
         cdef size_t neig = res[0].neig, i, j
@@ -694,10 +705,40 @@ cdef class ScatteringSystem:
                 'eigval_err':eigval_err,
                 'ranktest_SV':ranktest_SV,
                 'iri': iri,
+                'blochvector': blochvector
         }
 
         return retdict
 
+
+def empty_lattice_modes_xy(EpsMu epsmu, reciprocal_basis, wavevector, double maxomega):
+    '''Empty (2D, xy-plane) lattice mode (diffraction order) frequencies of a non-dispersive medium.
+
+    reciprocal_basis is of the "mutliplied by 2п" type.
+    '''
+    cdef double *omegas_c
+    cdef size_t n
+    cdef cart2_t k_, b1, b2
+    k_.x = wavevector[0]
+    k_.y = wavevector[1]
+    b1.x = reciprocal_basis[0][0]
+    b1.y = reciprocal_basis[0][1]
+    b2.x = reciprocal_basis[1][0]
+    b2.y = reciprocal_basis[1][1]
+    if(epsmu.n.imag != 0):
+        warnings.warn("Got complex refractive index", epsmu.n, "ignoring the imaginary part")
+    refindex = epsmu.n.real
+    n = qpms_emptylattice2_modes_maxfreq(&omegas_c, b1, b2, BASIS_RTOL,
+	    k_, GSL_CONST_MKSA_SPEED_OF_LIGHT / refindex, maxomega)
+    cdef np.ndarray[double, ndim=1] omegas = np.empty((n,), dtype=np.double)
+    cdef double[::1] omegas_v = omegas
+    cdef size_t i
+    for i in range(n):
+        omegas_v[i] = omegas_c[i]
+    free(omegas_c)
+    return omegas
+
+
 cdef class _ScatteringSystemAtOmegaK:
     '''
     Wrapper over the C qpms_scatsys_at_omega_k_t structure

qpms/qpms_cdefs.pxd

@@ -10,6 +10,9 @@ cdef extern from "gsl/gsl_errno.h":
     gsl_error_handler_t *gsl_set_error_handler(gsl_error_handler_t *new_handler)
     gsl_error_handler_t *gsl_set_error_handler_off();
 
+cdef extern from "gsl/gsl_const_mksa.h":
+    const double GSL_CONST_MKSA_SPEED_OF_LIGHT
+
 cdef extern from "qpms_types.h":
     cdef struct cart3_t:
         double x

qpms/scatsystem.c

@@ -2177,7 +2177,7 @@ beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
     complex double omega_centre, double omega_rr, double omega_ri, 
     size_t contour_npoints, 
     double rank_tol, size_t rank_sel_min, double res_tol) {
-  qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsys_finite_find_eigenmodes()");
+  qpms_ss_ensure_periodic_a(ss, "qpms_scatsys_finite_find_eigenmodes()");
   size_t n = ss->fecv_size; // matrix dimension
 
   beyn_contour_t *contour = beyn_contour_ellipse(omega_centre,