misc/lat2d_modes.py

   1 #!/usr/bin/env python3
   2
   3 import math
   4 from qpms.argproc import ArgParser, sfloat, annotate_pdf_metadata
   5
   6 ap = ArgParser(['const_real_background', 'lattice2d', 'multi_particle']) # TODO general analytical background
   7
   8 ap.add_argument("-k", nargs=2, type=sfloat, required=True, help='k vector', metavar=('K_X', 'K_Y'))
   9 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)")
  10 ap.add_argument("--rank-tol", type=float, required=False)
  11 ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)')
  12 ap.add_argument("-t", "--rank-tolerance", type=float, default=1e11)
  13 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')
  14 ap.add_argument("-T", "--residual-tolerance", type=float, default=2.)
  15 ap.add_argument("-N", type=int, default="150", help="Integration contour discretisation size")
  16
  17
  18 #TODO alternative specification of the contour by center and half-axes
  19 dospec = ap.add_argument_group("Eigenvalue search area by diffracted order specification", "Specification of eigenvalue search area by diffracted order number (requires constant real refractive index for background): the integration contour 'touches' the empty lattice band specified by -b, and its axis lying on the real axis reaches '-f'-way to the next diffractive order")
  20 dospec.add_argument("-d", "--band-index", type=int, 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.", required=True)
  21 dospec.add_argument("-f", "--interval-factor", type=float, default=0.1, help="Relative length of the integration ellipse axis w.r.t. the interval between two empty lattice bands; this should be be less than 1.") #TODO check
  22 dospec.add_argument("-i", "--imaginary-aspect-ratio", type=float, default=1., help="Aspect ratio of the integration ellipse (Im/Re); this should not exceed 1/interval_factor.")
  23
  24 ap.add_argument("-P", "--plot", action='store_true', help="if -p not given, plot to a default path")
  25 ap.add_argument("-O", "--plot-out", type=str, required=False, help="path to plot output (optional)")
  26
  27 a = ap.parse_args()
  28
  29 a1 = ap.direct_basis[0]
  30 a2 = ap.direct_basis[1]
  31
  32 particlestr = "modes" # TODO particle string specifier or some hash, do this in argproc.py
  33 defaultprefix = "%s_basis%gnm_%gnm__%gnm_%gnm_n%g_b%+d_k(%g_%g)um-1_cn%d" % (
  34             particlestr, a1[0]*1e9, a1[1]*1e9, a2[0]*1e9, a2[1]*1e9, a.refractive_index, a.band_index, a.k[0]*1e-6, a.k[1]*1e-6, a.N)
  35
  36
  37 import logging
  38 logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO)
  39
  40 import numpy as np
  41 import qpms
  42 from qpms.cybspec import BaseSpec
  43 from qpms.cytmatrices import CTMatrix
  44 from qpms.qpms_c import Particle, ScatteringSystem, empty_lattice_modes_xy
  45 from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude
  46 from qpms.constants import eV, hbar
  47 eh = eV/hbar
  48
  49 def inside_ellipse(point_xy, centre_xy, halfaxes_xy):
  50     x = point_xy[0] - centre_xy[0]
  51     y = point_xy[1] - centre_xy[1]
  52     ax = halfaxes_xy[0]
  53     ay = halfaxes_xy[1]
  54     return ((x/ax)**2 + (y/ay)**2) <= 1
  55
  56 beta = np.array(a.k)
  57
  58 if True: # TODO alternative specification of the contour by center and half-axes
  59     empty_freqs = empty_lattice_modes_xy(ap.background_epsmu, ap.reciprocal_basis2pi, np.array([0,0]), 1)
  60     empty_freqs = empty_lattice_modes_xy(ap.background_epsmu, ap.reciprocal_basis2pi, beta, (1+abs(a.band_index)) * empty_freqs[1])
  61
  62     # make the frequencies in the list unique
  63     empty_freqs = list(empty_freqs)
  64     i = 0
  65     while i < len(empty_freqs) - 1:
  66         if math.isclose(empty_freqs[i], empty_freqs[i+1], rel_tol=1e-13):
  67             del empty_freqs[i+1]
  68         else:
  69             i += 1
  70
  71     logging.info("Empty freqs: %s", str(empty_freqs))
  72     logging.info("Empty freqs (eV): %s", str([ff / eh for ff in empty_freqs]))
  73     if a.band_index > 0:
  74         top = empty_freqs[a.band_index]
  75         bottom = empty_freqs[a.band_index - 1]
  76         lebeta_om = bottom
  77     else: # a.band_index < 0
  78         top = empty_freqs[abs(a.band_index) - 1]
  79         bottom = empty_freqs[abs(a.band_index) - 2] if abs(a.band_index) > 1 else 0.
  80         lebeta_om = top
  81     #print(top,bottom,lebeta_om)
  82     freqradius = .5 * (top - bottom) * a.interval_factor
  83
  84     centfreq = bottom + freqradius if a.band_index > 0 else top - freqradius
  85     if freqradius == 0:
  86             raise ValueError("Integration contour radius is set to zero. Are you trying to look below the lowest empty lattice band at the gamma point?")
  87
  88     freqradius *= (1-1e-13) # to not totally touch the singularities
  89
  90 logging.info("Direct lattice basis: %s" % str(ap.direct_basis))
  91 logging.info("Reciprocal lattice basis: %s" % str(ap.reciprocal_basis2pi))
  92
  93 ss, ssw = ScatteringSystem.create(ap.get_particles(), ap.background_emg, centfreq, latticebasis=ap.direct_basis)
  94
  95 logging.info("Finding eigenvalues around %s (= %s eV)" % (str(centfreq), str(centfreq/eh)))
  96 logging.info("Real half-axis %s (= %s eV)" % (str(freqradius), str(freqradius/eh)))
  97 logging.info("Imaginary half-axis %s (= %s eV)" % (str(freqradius*a.imaginary_aspect_ratio), str(freqradius*a.imaginary_aspect_ratio/eh)))
  98
  99 with qpms.pgsl_ignore_error(15):
 100     res = ss.find_modes(centfreq, freqradius, freqradius * a.imaginary_aspect_ratio,
 101             blochvector = a.k, contour_points = a.N, rank_tol = a.rank_tolerance,
 102             res_tol = a.residual_tolerance, rank_min_sel = a.min_candidates)
 103
 104 logging.info("Eigenfrequencies found: %s" % str(res['eigval']))
 105 logging.info("Eigenfrequencies found (eV): %s" % str(res['eigval'] / eh))
 106
 107 res['inside_contour'] = inside_ellipse((res['eigval'].real, res['eigval'].imag),
 108          (centfreq.real, centfreq.imag), (freqradius, freqradius * a.imaginary_aspect_ratio))
 109
 110 #res['refractive_index_internal'] = [emg(om).n for om in res['eigval']]
 111
 112 #del res['omega']  If contour points are not needed...
 113 #del res['ImTW'] # not if dbg=false anyway
 114 outfile = defaultprefix + ".npz" if a.output is None else a.output
 115 np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, empty_freqs=np.array(empty_freqs),
 116         ss_positions=ss.positions, ss_fullvec_poffsets=ss.fullvec_poffsets,
 117         ss_fullvec_psizes=ss.fullvec_psizes,
 118         ss_bspecs_flat = np.concatenate(ss.bspecs),
 119         ss_lattice_basis=ss.lattice_basis, ss_reciprocal_basis = ss.reciprocal_basis,
 120         **res)
 121 logging.info("Saved to %s" % outfile)
 122
 123
 124 if a.plot or (a.plot_out is not None):
 125     if len(res['eigval']) == 0:
 126         logging.info("No eigenvalues found; nothing to plot")
 127         exit(1)
 128     imcut = np.linspace(0, -freqradius)
 129     recut1 = np.sqrt(lebeta_om**2+imcut**2) # incomplete Gamma-related cut
 130     recut2 = np.sqrt((lebeta_om/2)**2-imcut**2) + lebeta_om/2 # odd-power-lilgamma-related cut
 131
 132     import matplotlib
 133     matplotlib.use('pdf')
 134     from matplotlib import pyplot as plt
 135     from matplotlib.backends.backend_pdf import PdfPages
 136
 137     fig = plt.figure()
 138     ax = fig.add_subplot(111,)
 139     #ax.plot(res['omega'].real/eh, res['omega'].imag/eh*1e3, ':') #res['omega'] not implemented in ScatteringSystem
 140     ax.add_artist(matplotlib.patches.Ellipse((centfreq.real/eh, centfreq.imag/eh*1e3),
 141         2*freqradius/eh, 2*freqradius*a.imaginary_aspect_ratio/eh*1e3, fill=False,
 142         ls=':'))
 143     ax.scatter(x=res['eigval'].real/eh, y=res['eigval'].imag/eh*1e3 , c = res['inside_contour']
 144                       )
 145     ax.plot(recut1/eh, imcut/eh*1e3)
 146     ax.plot(recut2/eh, imcut/eh*1e3)
 147     for i,om in enumerate(res['eigval']):
 148             ax.annotate(str(i), (om.real/eh, om.imag/eh*1e3))
 149     xmin = np.amin(res['eigval'].real)/eh
 150     xmax = np.amax(res['eigval'].real)/eh
 151     xspan = xmax-xmin
 152     ymin = np.amin(res['eigval'].imag)/eh*1e3
 153     ymax = np.amax(res['eigval'].imag)/eh*1e3
 154     yspan = ymax-ymin
 155     ax.set_xlim([xmin-.1*xspan, xmax+.1*xspan])
 156     ax.set_ylim([ymin-.1*yspan, ymax+.1*yspan])
 157     ax.set_xlabel('$\hbar \Re \omega / \mathrm{eV}$')
 158     ax.set_ylabel('$\hbar \Im \omega / \mathrm{meV}$')
 159
 160     plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
 161     with PdfPages(plotfile) as pdf:
 162         pdf.savefig(fig)
 163         annotate_pdf_metadata(pdf, scriptname='lat2d_modes.py')
 164
 165 exit(0)
 166