misc/dispersion2D-SVD.py

   1 #!/usr/bin/env python3
   2
   3
   4 import argparse, re, random, string
   5 import subprocess
   6 from scipy.constants import hbar, e as eV, pi, c
   7
   8 def make_action_sharedlist(opname, listname):
   9     class opAction(argparse.Action):
  10         def __call__(self, parser, args, values, option_string=None):
  11             if (not hasattr(args, listname)) or getattr(args, listname) is None:
  12                 setattr(args, listname, list())
  13             getattr(args,listname).append((opname, values))
  14     return opAction
  15
  16 parser = argparse.ArgumentParser()
  17 #TODO? použít type=argparse.FileType('r') ?
  18 parser.add_argument('--TMatrix', action='store', required=True, help='Path to TMatrix file')
  19 #parser.add_argument('--griddir', action='store', required=True, help='Path to the directory with precalculated translation operators')
  20 parser.add_argument('--output_prefix', action='store', required=True, help='Prefix to the pdf and/or npz output (will be appended frequency and hexside)')
  21 #sizepar = parser.add_mutually_exclusive_group(required=True)
  22 parser.add_argument('--hexside', action='store', type=float, required=True, help='Lattice hexagon size length')
  23 parser.add_argument('--output', action='store', help='Path to output PDF')
  24 parser.add_argument('--store_SVD', action='store_false', help='If specified without --SVD_output, it will save the data in a file named as the PDF output, but with .npz extension instead')
  25 parser.add_argument('--plot_TMatrix', action='store_true', help='Visualise TMatrix on the first page of the output')
  26 #parser.add_argument('--SVD_output', action='store', help='Path to output singular value decomposition result')
  27 parser.add_argument('--nSV', action='store', metavar='N', type=int, default=1, help='Store and draw N minimun singular values')
  28 parser.add_argument('--maxlayer', action='store', type=int, default=100, help='How far to sum the lattice points to obtain the dispersion')
  29 parser.add_argument('--scp_to', action='store', metavar='N', type=str, help='SCP the output files to a given destination')
  30 parser.add_argument('--background_permittivity', action='store', type=float, default=1., help='Background medium relative permittivity (default 1)')
  31 parser.add_argument('--eVfreq', action='store', required=True, type=float, help='Frequency in eV')
  32 parser.add_argument('--kdensity', action='store', type=int, default=33, help='Number of k-points per x-axis segment')
  33 parser.add_argument('--lMax', action='store', type=int, help='Override lMax from the TMatrix file')
  34 #TODO some more sophisticated x axis definitions
  35 parser.add_argument('--gaussian', action='store', type=float, metavar='σ', help='Use a gaussian envelope for weighting the interaction matrix contributions (depending on the distance), measured in unit cell lengths (?) FIxME).')
  36 popgrp=parser.add_argument_group(title='Operations')
  37 popgrp.add_argument('--tr', dest='ops', action=make_action_sharedlist('tr', 'ops'), default=list()) # the default value for dest can be set once
  38 popgrp.add_argument('--tr0', dest='ops', action=make_action_sharedlist('tr0', 'ops'))
  39 popgrp.add_argument('--tr1', dest='ops', action=make_action_sharedlist('tr1', 'ops'))
  40 popgrp.add_argument('--sym', dest='ops', action=make_action_sharedlist('sym', 'ops'))
  41 popgrp.add_argument('--sym0', dest='ops', action=make_action_sharedlist('sym0', 'ops'))
  42 popgrp.add_argument('--sym1', dest='ops', action=make_action_sharedlist('sym1', 'ops'))
  43 #popgrp.add_argument('--mult', dest='ops', nargs=3, metavar=('INCSPEC', 'SCATSPEC', 'MULTIPLIER'), action=make_action_sharedlist('mult', 'ops'))
  44 #popgrp.add_argument('--mult0', dest='ops', nargs=3, metavar=('INCSPEC', 'SCATSPEC', 'MULTIPLIER'), action=make_action_sharedlist('mult0', 'ops'))
  45 #popgrp.add_argument('--mult1', dest='ops', nargs=3, metavar=('INCSPEC', 'SCATSPEC', 'MULTIPLIER'), action=make_action_sharedlist('mult1', 'ops'))
  46 popgrp.add_argument('--multl', dest='ops', nargs=3, metavar=('INCL[,INCL,...]', 'SCATL[,SCATL,...]', 'MULTIPLIER'), action=make_action_sharedlist('multl', 'ops'))
  47 popgrp.add_argument('--multl0', dest='ops', nargs=3, metavar=('INCL[,INCL,...]', 'SCATL[,SCATL,...]', 'MULTIPLIER'), action=make_action_sharedlist('multl0', 'ops'))
  48 popgrp.add_argument('--multl1', dest='ops', nargs=3, metavar=('INCL[,INCL,...]', 'SCATL[,SCATL,...]', 'MULTIPLIER'), action=make_action_sharedlist('multl1', 'ops'))
  49 parser.add_argument('--frequency_multiplier', action='store', type=float, default=1., help='Multiplies the frequencies in the TMatrix file by a given factor.')
  50 # TODO enable more flexible per-sublattice specification
  51 pargs=parser.parse_args()
  52 print(pargs)
  53
  54 maxlayer=pargs.maxlayer
  55 hexside=pargs.hexside
  56 eVfreq = pargs.eVfreq
  57 freq = eVfreq*eV/hbar
  58
  59 TMatrix_file = pargs.TMatrix
  60 pdfout = pargs.output if pargs.output else (
  61         '%s_%dnm_%.4f.pdf' % (pargs.output_prefix,hexside/1e-9,eVfreq) if pargs.output_prefix else
  62         (''.join(random.choice(string.ascii_uppercase + string.digits) for _ in range(10)) + '.pdf'))
  63 print(pdfout)
  64 if(pargs.store_SVD):
  65     if re.search('.pdf$', pdfout):
  66         svdout = re.sub('.pdf$', r'.npz', pdfout)
  67     else:
  68         svdout = pdfout + '.npz'
  69 else:
  70      svdout = None
  71
  72
  73 epsilon_b = pargs.background_permittivity #2.3104
  74 gaussianSigma = pargs.gaussian if pargs.gaussian else None # hexside * 222 / 7
  75 interpfreqfactor = pargs.frequency_multiplier
  76 scp_dest = pargs.scp_to if pargs.scp_to else None
  77 kdensity = pargs.kdensity
  78 svn = pargs.nSV
  79
  80 # TODO multiplier operation definitions and parsing
  81 #factor13inc = 10
  82 #factor13scat=10
  83
  84 ops = list()
  85 opre = re.compile('(tr|sym|copy|multl|mult)(\d*)')
  86 for oparg in pargs.ops:
  87     opm = opre.match(oparg[0])
  88     if opm:
  89         ops.append(((opm.group(2),) if opm.group(2) else (0,1), opm.group(1), oparg[1]))
  90     else:
  91         raise # should not happen
  92 print(ops)
  93
  94 #ops = (
  95 #    # co, typ operace (symetrizace / transformace / kopie), specifikace (operace nebo zdroj),
  96 #    # co: 0, 1, (0,1), (0,), (1,), #NI: 'all'
  97 #    # typ operace: sym, tr, copy
  98 #    # specifikace:
  99 #    # sym, tr: 'σ_z', 'σ_y', 'C2'; sym: 'C3',
 100 #    # copy: 0, 1 (zdroj)
 101 #    ((0,1), 'sym', 'σ_z'),
 102 #    #((0,1), 'sym', 'σ_x'),
 103 #    #((0,1), 'sym', 'σ_y'),
 104 #    ((0,1), 'sym', 'C3'),
 105 #    ((1), 'tr', 'C2'),
 106 #
 107 #)
 108
 109 # -----------------finished basic CLI parsing (except for op arguments) ------------------
 110 import time
 111 begtime=time.time()
 112
 113 from matplotlib.path import Path
 114 import matplotlib.patches as patches
 115 import matplotlib.pyplot as plt
 116 import qpms
 117 import numpy as np
 118 import os, sys, warnings, math
 119 from matplotlib import pyplot as plt
 120 from matplotlib.backends.backend_pdf import PdfPages
 121 from scipy import interpolate
 122 nx = None
 123 s3 = math.sqrt(3)
 124
 125
 126
 127 pdf = PdfPages(pdfout)
 128
 129 # In[3]:
 130
 131 # specifikace T-matice zde
 132 cdn = c/ math.sqrt(epsilon_b)
 133 TMatrices_orig, freqs_orig, freqs_weirdunits_orig, lMaxTM = qpms.loadScuffTMatrices(TMatrix_file)
 134 lMax = lMaxTM
 135 if pargs.lMax:
 136     lMax = pargs.lMax if pargs.lMax else lMaxTM
 137 my, ny = qpms.get_mn_y(lMax)
 138 nelem = len(my)
 139 if pargs.lMax: #force commandline specified lMax
 140     TMatrices_orig = TMatrices_orig[...,0:nelem,:,0:nelem]
 141
 142 TMatrices = np.array(np.broadcast_to(TMatrices_orig[:,nx,:,:,:,:],(len(freqs_orig),2,2,nelem,2,nelem)) )
 143
 144 #TMatrices[:,:,:,:,:,ny==3] *= factor13inc
 145 #TMatrices[:,:,:,ny==3,:,:] *= factor13scat
 146 xfl = qpms.xflip_tyty(lMax)
 147 yfl = qpms.yflip_tyty(lMax)
 148 zfl = qpms.zflip_tyty(lMax)
 149 c2rot = qpms.apply_matrix_left(qpms.yflip_yy(3),qpms.xflip_yy(3),-1)
 150
 151 reCN = re.compile('(\d*)C(\d+)')
 152 #TODO C nekonečno
 153
 154 for op in ops:
 155     if op[0] == 'all':
 156         targets = (0,1)
 157     elif isinstance(op[0],int):
 158         targets = (op[0],)
 159     else:
 160         targets = op[0]
 161
 162     if op[1] == 'sym':
 163         mCN = reCN.match(op[2]) # Fuck van Rossum for not having assignments inside expressions
 164         if op[2] == 'σ_z':
 165             for t in targets:
 166                 TMatrices[:,t] = (TMatrices[:,t] + qpms.apply_ndmatrix_left(zfl,qpms.apply_ndmatrix_left(zfl, TMatrices[:,t], (-4,-3)),(-2,-1)))/2
 167         elif op[2] == 'σ_y':
 168             for t in targets:
 169                 TMatrices[:,t] = (TMatrices[:,t] + qpms.apply_ndmatrix_left(yfl,qpms.apply_ndmatrix_left(yfl, TMatrices[:,t], (-4,-3)),(-2,-1)))/2
 170         elif op[2] == 'σ_x':
 171             for t in targets:
 172                 TMatrices[:,t] = (TMatrices[:,t] + qpms.apply_ndmatrix_left(xfl,qpms.apply_ndmatrix_left(xfl, TMatrices[:,t], (-4,-3)),(-2,-1)))/2
 173         elif op[2] == 'C2': # special case of the latter
 174             for t in targets:
 175                 TMatrices[:,t] = (TMatrices[:,t] + qpms.apply_matrix_left(c2rot,qpms.apply_matrix_left(c2rot, TMatrices[:,t], -3),-1))/2
 176         elif mCN:
 177             rotN = int(mCN.group(2))
 178             TMatrix_contribs = np.empty((rotN,TMatrices.shape[0],2,nelem,2,nelem), dtype=np.complex_)
 179             for t in targets:
 180                 for i in range(rotN):
 181                     rotangle = 2*np.pi*i / rotN
 182                     rot =  qpms.WignerD_yy_fromvector(lMax,np.array([0,0,rotangle]))
 183                     rotinv = qpms.WignerD_yy_fromvector(lMax,np.array([0,0,-rotangle]))
 184                     TMatrix_contribs[i] = qpms.apply_matrix_left(rot,qpms.apply_matrix_left(rotinv, TMatrices[:,t], -3),-1)
 185                 TMatrices[:,t] = np.sum(TMatrix_contribs, axis=0) / rotN
 186         else:
 187             raise
 188     elif op[1] == 'tr':
 189         mCN = reCN.match(op[2]) # Fuck van Rossum for not having assignments inside expressions
 190         if op[2] == 'σ_z':
 191             for t in targets:
 192                 TMatrices[:,t] = qpms.apply_ndmatrix_left(zfl,qpms.apply_ndmatrix_left(zfl, TMatrices[:,t], (-4,-3)),(-2,-1))
 193         elif op[2] == 'σ_y':
 194             for t in targets:
 195                 TMatrices[:,t] = qpms.apply_ndmatrix_left(yfl,qpms.apply_ndmatrix_left(yfl, TMatrices[:,t], (-4,-3)),(-2,-1))
 196         elif op[2] == 'σ_x':
 197             for t in targets:
 198                 TMatrices[:,t] = qpms.apply_ndmatrix_left(xfl,qpms.apply_ndmatrix_left(xfl, TMatrices[:,t], (-4,-3)),(-2,-1))
 199         elif op[2] == 'C2':
 200             for t in targets:
 201                 TMatrices[:,t] = qpms.apply_matrix_left(c2rot,qpms.apply_matrix_left(c2rot, TMatrices[:,t], -3),-1)
 202         elif mCN:
 203             rotN = int(mCN.group(2))
 204             power = int(mCN.group(1)) if mCN.group(1) else 1
 205             TMatrix_contribs = np.empty((rotN,TMatrices.shape[0],2,nelem,2,nelem), dtype=np.complex_)
 206             for t in targets:
 207                 rotangle = 2*np.pi*power/rotN
 208                 rot = qpms.WignerD_yy_fromvector(lMax, np.array([0,0,rotangle]))
 209                 rotinv = qpms.WignerD_yy_fromvector(lMax, np.array([0,0,-rotangle]))
 210                 TMatrices[:,t] = qpms.apply_matrix_left(rot, qpms.apply_matrix_left(rotinv, TMatrices[:,t], -3),-1)
 211         else:
 212             raise
 213     elif op[1] == 'copy':
 214         raise # not implemented
 215     elif op[1] == 'mult':
 216         raise # not implemented
 217     elif op[1] == 'multl':
 218         incy = np.full((nelem,), False, dtype=bool)
 219         for incl in op[2][0].split(','):
 220             l = int(incl)
 221             incy += (l == ny)
 222         scaty = np.full((nelem,), False, dtype=bool)
 223         for scatl in op[2][1].split(','):
 224             l = int(scatl)
 225             scaty += (l == ny)
 226         for t in targets:
 227             TMatrices[np.ix_(np.arange(TMatrices.shape[0]),np.array([t]),np.array([0,1]),scaty,np.array([0,1]),incy)] *= float(op[2][2])
 228     else:
 229         raise #unknown operation; should not happen
 230
 231 TMatrices_interp = interpolate.interp1d(freqs_orig*interpfreqfactor, TMatrices, axis=0, kind='linear',fill_value="extrapolate")
 232
 233
 234
 235 # In[4]:
 236 if(pargs.plot_TMatrix):
 237     om = np.linspace(np.min(freqs_orig), np.max(freqs_orig),100)
 238     TMatrix0ip = np.reshape(TMatrices_interp(om)[:,0], (len(om), 2*nelem*2*nelem))
 239     f, axa = plt.subplots(2, 2, figsize=(15,15))
 240
 241     #print(TMatrices.shape)
 242     #plt.plot(om, TMatrices[:,0,0,0,0].imag,'r',om, TMatrices[:,0,0,0,0].real,'r--',om, TMatrices[:,0,2,0,2].imag,'b',om, TMatrices[:,0,2,0,2].real,'b--'))
 243
 244     ax = axa[0,0]
 245     ax2 = ax.twiny()
 246     ax2.set_xlim([ax.get_xlim()[0]/eV*hbar,ax.get_xlim()[1]/eV*hbar])
 247     ax.plot(
 248              om, TMatrix0ip[:,:].imag,'-',om, TMatrix0ip[:,:].real,'--',
 249     )
 250     ax = axa[0,1]
 251     ax2 = ax.twiny()
 252     ax2.set_xlim([ax.get_xlim()[0]/eV*hbar,ax.get_xlim()[1]/eV*hbar])
 253     ax.plot(
 254              om, abs(TMatrix0ip[:,:]),'-'
 255     )
 256     ax.set_yscale('log')
 257
 258     ax = axa[1,1]
 259     ax2 = ax.twiny()
 260     ax2.set_xlim([ax.get_xlim()[0]/eV*hbar,ax.get_xlim()[1]/eV*hbar])
 261     ax.plot(
 262              om, np.unwrap(np.angle(TMatrix0ip[:,:]),axis=0),'-'
 263     )
 264
 265     ax = axa[1,0]
 266     ax.text(0.5,0.5,str(pargs).replace(',',',\n'),horizontalalignment='center',verticalalignment='center',transform=ax.transAxes)
 267     pdf.savefig(f)
 268
 269
 270 # In[ ]:
 271
 272 '''
 273 #kdensity = 66 #defined from cl arguments
 274 bz_0 = np.array((0,0,0.,))
 275 bz_K1 = np.array((1.,0,0))*4*np.pi/3/hexside/s3
 276 bz_K2 = np.array((1./2.,s3/2,0))*4*np.pi/3/hexside/s3
 277 bz_M = np.array((3./4, s3/4,0))*4*np.pi/3/hexside/s3
 278 k0Mlist = bz_0 + (bz_M-bz_0) * np.linspace(0,1,kdensity)[:,nx]
 279 kMK1list = bz_M + (bz_K1-bz_M) * np.linspace(0,1,kdensity)[:,nx]
 280 kK10list = bz_K1 + (bz_0-bz_K1) * np.linspace(0,1,kdensity)[:,nx]
 281 k0K2list = bz_0 + (bz_K2-bz_0) * np.linspace(0,1,kdensity)[:,nx]
 282 kK2Mlist = bz_K2 + (bz_M-bz_K2) * np.linspace(0,1,kdensity)[:,nx]
 283 B1 = 2* bz_K1 - bz_K2
 284 B2 = 2* bz_K2 - bz_K1
 285 klist = np.concatenate((k0Mlist,kMK1list,kK10list,k0K2list,kK2Mlist), axis=0)
 286 kxmaplist = np.concatenate((np.array([0]),np.cumsum(np.linalg.norm(np.diff(klist, axis=0), axis=-1))))
 287 '''
 288 klist = qpms.generate_trianglepoints(kdensity, v3d=True, include_origin=True)*3*math.pi/(3*kdensity*hexside)
 289 TMatrices_om = TMatrices_interp(freq)
 290
 291 svdres = qpms.hexlattice_zsym_getSVD(lMax=lMax, TMatrices_om=TMatrices_om, epsilon_b=epsilon_b, hexside=hexside, maxlayer=maxlayer,
 292         omega=freq, klist=klist, gaussianSigma=gaussianSigma, onlyNmin=(0 if svdout else svn))
 293 if svdout:
 294    ((svUfullTElist, svSfullTElist, svVfullTElist), (svUfullTMlist, svSfullTMlist, svVfullTMlist)) = svdres
 295    (minsvElist, minsvTMlist) = (svSfullTElist[...,-svn:], svSfullTMlist[...,-svn:])
 296 else:
 297     minsvTElist, minsvTMlist = svdres
 298
 299
 300
 301 ''' The new pretty diffracted order drawing '''
 302 maxlayer_reciprocal=4
 303 cdn = c/ math.sqrt(epsilon_b)
 304 bz_0 = np.array((0,0,))
 305 bz_K1 = np.array((1.,0))*4*np.pi/3/hexside/s3
 306 bz_K2 = np.array((1./2.,s3/2))*4*np.pi/3/hexside/s3
 307 bz_M = np.array((3./4, s3/4))*4*np.pi/3/hexside/s3
 308
 309 # reciprocal lattice basis
 310 B1 = 2* bz_K1 - bz_K2
 311 B2 = 2* bz_K2 - bz_K1
 312
 313 if svdout:
 314     np.savez(svdout, omega = freq, klist = klist, bzpoints = np.array([bz_0, bz_K1, bz_K2, bz_M, B1, B2]),
 315                      uTE = svUfullTElist,
 316                      vTE = svVfullTElist,
 317                      sTE = svSfullTElist,
 318                      uTM = svUfullTMlist,
 319                      vTM = svVfullTMlist,
 320                      sTM = svSfullTMlist,
 321     )
 322
 323 k2density = 100
 324 k0Mlist = bz_0 + (bz_M-bz_0) * np.linspace(0,1,k2density)[:,nx]
 325 kMK1list = bz_M + (bz_K1-bz_M) * np.linspace(0,1,k2density)[:,nx]
 326 kK10list = bz_K1 + (bz_0-bz_K1) * np.linspace(0,1,k2density)[:,nx]
 327 k0K2list = bz_0 + (bz_K2-bz_0) * np.linspace(0,1,k2density)[:,nx]
 328 kK2Mlist = bz_K2 + (bz_M-bz_K2) * np.linspace(0,1,k2density)[:,nx]
 329 k2list = np.concatenate((k0Mlist,kMK1list,kK10list,k0K2list,kK2Mlist), axis=0)
 330 kxmaplist = np.concatenate((np.array([0]),np.cumsum(np.linalg.norm(np.diff(k2list, axis=0), axis=-1))))
 331
 332 centers2=qpms.generate_trianglepoints(maxlayer_reciprocal, v3d = False, include_origin= True)*4*np.pi/3/hexside
 333 rot90 = np.array([[0,-1],[1,0]])
 334 centers2=np.dot(centers2,rot90)
 335
 336 import matplotlib.pyplot as plt
 337 import matplotlib
 338 from matplotlib.path import Path
 339 import matplotlib.patches as patches
 340 cmap = matplotlib.cm.prism
 341 colormax = np.amax(np.linalg.norm(centers2,axis=0))
 342
 343
 344 # In[ ]:
 345 for minN in reversed(range(svn)):
 346     f, axes = plt.subplots(1,3, figsize=(20,4.8))
 347     ax = axes[0]
 348     sc = ax.scatter(klist[:,0], klist[:,1], c = np.clip(np.abs(minsvTElist[:,minN]),0,1), lw=0)
 349     for center in centers2:
 350         circle=plt.Circle((center[0],center[1]),omega/cdn, facecolor='none', edgecolor=cmap(np.linalg.norm(center)/colormax),lw=0.5)
 351         ax.add_artist(circle)
 352     verts = [(math.cos(math.pi*i/3)*4*np.pi/3/hexside/s3,math.sin(math.pi*i/3)*4*np.pi/3/hexside/s3) for i in range(6 +1)]
 353     codes = [Path.MOVETO,Path.LINETO,Path.LINETO,Path.LINETO,Path.LINETO,Path.LINETO,Path.CLOSEPOLY,]
 354     path = Path(verts, codes)
 355     patch = patches.PathPatch(path, facecolor='none', edgecolor='black',  lw=1)
 356     ax.add_patch(patch)
 357     ax.set_xticks([])
 358     ax.set_yticks([])
 359     ax.title.set_text('E in-plane ("TE")')
 360     f.colorbar(sc,ax=ax)
 361
 362
 363     ax = axes[1]
 364     sc = ax.scatter(klist[:,0], klist[:,1], c = np.clip(np.abs(minsvTMlist[:,minN]),0,1), lw=0)
 365     for center in centers2:
 366         circle=plt.Circle((center[0],center[1]),omega/cdn, facecolor='none', edgecolor=cmap(np.linalg.norm(center)/colormax),lw=0.5)
 367         ax.add_artist(circle)
 368     verts = [(math.cos(math.pi*i/3)*4*np.pi/3/hexside/s3,math.sin(math.pi*i/3)*4*np.pi/3/hexside/s3) for i in range(6 +1)]
 369     codes = [Path.MOVETO,Path.LINETO,Path.LINETO,Path.LINETO,Path.LINETO,Path.LINETO,Path.CLOSEPOLY,]
 370     path = Path(verts, codes)
 371     patch = patches.PathPatch(path, facecolor='none', edgecolor='black',  lw=1)
 372     ax.add_patch(patch)
 373     ax.set_xticks([])
 374     ax.set_yticks([])
 375     ax.title.set_text('E perpendicular ("TM")')
 376     f.colorbar(sc,ax=ax)
 377
 378     ax = axes[2]
 379     for center in centers2:
 380         ax.plot(kxmaplist, np.linalg.norm(k2list-center,axis=-1)*cdn, '-', color=cmap(np.linalg.norm(center)/colormax))
 381
 382     #ax.set_xlim([np.min(kxmlarr),np.max(kxmlarr)])
 383     #ax.set_ylim([np.min(omegalist),np.max(omegalist)])
 384     xticklist = [0, kxmaplist[len(k0Mlist)-1], kxmaplist[len(k0Mlist)+len(kMK1list)-1], kxmaplist[len(k0Mlist)+len(kMK1list)+len(kK10list)-1], kxmaplist[len(k0Mlist)+len(kMK1list)+len(kK10list)+len(k0K2list)-1], kxmaplist[len(k0Mlist)+len(kMK1list)+len(kK10list)+len(k0K2list)+len(kK2Mlist)-1]]
 385     ax.set_xticks(xticklist)
 386     for xt in xticklist:
 387         ax.axvline(xt, ls='dotted', lw=0.3,c='k')
 388     ax.set_xticklabels(['Γ', 'M', 'K', 'Γ', 'K\'','M'])
 389     ax.axhline(omega, c='black')
 390     ax.set_ylim([0,5e15])
 391     ax2 = ax.twinx()
 392     ax2.set_ylim([ax.get_ylim()[0]/eV*hbar,ax.get_ylim()[1]/eV*hbar])
 393
 394     pdf.savefig(f)
 395
 396 pdf.close()
 397
 398 if scp_dest:
 399     subprocess.run(['scp', pdfout, scp_dest])
 400     if svdout:
 401         subprocess.run(['scp', svdout, scp_dest])
 402
 403 print(time.strftime("%H.%M:%S",time.gmtime(time.time()-begtime)))