| blob | 08adb4836d443a406427f16aef229700b1834b68 |
1 {
2 "cells": [
3 {
4 "cell_type": "code",
5 "execution_count": 1,
6 "metadata": {},
7 "outputs": [],
8 "source": [
9 "%matplotlib inline\n",
10 "import re\n",
11 "import numpy as np\n",
12 "import matplotlib\n",
13 "from matplotlib import pyplot as plt\n",
14 "from scipy.constants import hbar, e as eV, c\n",
15 "eh = eV/hbar\n",
16 "import glob\n",
17 "def ri(z): return (z.real, z.imag)\n",
18 "#m = re.compile(r\"([^_]+)_r([0-9.]+)nm_\")\n",
19 "#removek = re.compile(r\"(k\\([^)]+\\)um-1_)\")\n",
20 "remover = re.compile(r\"r[0-9.]+nm_\")\n",
21 "\n",
22 "\n",
23 "markerdict = {\n",
24 " 4: \"3\",\n",
25 " -4: \"4\",\n",
26 " 3: \"^\",\n",
27 " -3: \"v\",\n",
28 " -2: 'x',\n",
29 " 2: '+',\n",
30 " 1: 's',\n",
31 " -1: 'd',\n",
32 "}\n",
33 "\n",
34 "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
35 "colors = prop_cycle.by_key()['color']\n",
36 "colordict = {i: colors[(i+1)] for i in range(-4,8)}\n",
37 "\n",
38 "def markerfun(b):\n",
39 " if b in markerdict.keys():\n",
40 " return markerdict[b]\n",
41 " else: return 'X'\n",
42 "\n",
43 "def colorfun(b):\n",
44 " if (b+1) in colordict.keys():\n",
45 " return colordict[b+1]\n",
46 " else: return colordict[0]"
47 ]
48 },
49 {
50 "cell_type": "code",
51 "execution_count": 2,
52 "metadata": {},
53 "outputs": [],
54 "source": [
55 "allfiles=glob.glob('*sph*k(0_0)*.npz')\n",
56 "allgraphs=dict()\n",
57 "for f in allfiles:\n",
58 " base = remover.sub('', f)\n",
59 " if base in allgraphs.keys():\n",
60 " allgraphs[base] += 1\n",
61 " else:\n",
62 " allgraphs[base] = 1\n",
63 "for k in sorted(allgraphs.keys()):\n",
64 " print(k, allgraphs[k])"
65 ]
66 },
67 {
68 "cell_type": "code",
69 "execution_count": 3,
70 "metadata": {},
71 "outputs": [
72 {
73 "ename": "FileNotFoundError",
74 "evalue": "[Errno 2] No such file or directory: 'projectors_D4h_lMax1.npz'",
75 "output_type": "error",
76 "traceback": [
77 "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
78 "\u001b[0;31mFileNotFoundError\u001b[0m Traceback (most recent call last)",
79 "\u001b[0;32m<ipython-input-3-0c266089be08>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mlMaxes\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mlMax\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mlMax\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mlMax\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlMaxes\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mproj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'projectors_D4h_lMax%d.npz'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mlMax\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 6\u001b[0m \u001b[0mirlabels\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msorted\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mproj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkeys\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0mproj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mproj\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mf\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mirlabels\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
80 "\u001b[0;32m~/.local/lib/python3.7/site-packages/numpy-1.17.3-py3.7-linux-x86_64.egg/numpy/lib/npyio.py\u001b[0m in \u001b[0;36mload\u001b[0;34m(file, mmap_mode, allow_pickle, fix_imports, encoding)\u001b[0m\n\u001b[1;32m 426\u001b[0m \u001b[0mown_fid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 427\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 428\u001b[0;31m \u001b[0mfid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mos_fspath\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfile\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"rb\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 429\u001b[0m \u001b[0mown_fid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 430\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
81 "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'projectors_D4h_lMax1.npz'"
82 ]
83 }
84 ],
85 "source": [
86 "projectors = dict()\n",
87 "projectors_list = dict()\n",
88 "lMaxes = [lMax for lMax in range(1,6)]\n",
89 "for lMax in lMaxes:\n",
90 " proj = np.load('projectors_D4h_lMax%d.npz' % lMax)\n",
91 " irlabels = sorted(proj.keys())\n",
92 " proj = {f: proj[f] for f in irlabels}\n",
93 " proj_list = [proj[irlabels[i]] for i in range(len(proj))]\n",
94 " projectors[lMax] = proj\n",
95 " projectors_list[lMax] = proj_list\n",
96 "globpattern = '*sph_r*_p580nmx580nm_mAu_n1.52_b?2_k(0_0)um-1_L?_cn???.npz'\n",
97 "filenames=glob.glob(globpattern)\n",
98 "plotfilename = 'collected_' + globpattern.replace('*', 'XXX').replace('?', 'X').replace('npz','pdf')\n",
99 "print(filenames[:4], plotfilename)\n"
100 ]
101 },
102 {
103 "cell_type": "code",
104 "execution_count": 41,
105 "metadata": {},
106 "outputs": [],
107 "source": [
108 "#projectors\n",
109 "#glob.glob('cyl_r100nm*L3*3100.npz')\n",
110 "#glob.glob('sph_r100*m5*.npz')\n",
111 "#dat['meta'][()],list(dat.keys())"
112 ]
113 },
114 {
115 "cell_type": "code",
116 "execution_count": null,
117 "metadata": {},
118 "outputs": [],
119 "source": []
120 },
121 {
122 "cell_type": "code",
123 "execution_count": 42,
124 "metadata": {},
125 "outputs": [
126 {
127 "name": "stdout",
128 "output_type": "stream",
129 "text": [
130 "inpure result detected [1. 0.99999999 1. 0.97991334 0.99999996 0.9999989\n",
131 " 0.99999983 0.99999966 0.99999322 0.99999721 0.99999653] [3.28735741e-04 2.66532534e-05 2.47011478e-05 1.45012420e-01\n",
132 " 2.44785416e-04 7.05405359e-04 1.60203586e-03 1.71245137e-03\n",
133 " 1.03244480e-02 9.18732728e-03 1.18651583e-02]\n",
134 "inpure result detected [1. 1. 0.99999998 0.99999999 0.99999996 0.96608887\n",
135 " 0.99999852 0.99999397 0.99998951 0.99999912 0.99982435] [2.66223026e-04 2.12357147e-05 3.54211968e-05 1.06651057e-04\n",
136 " 2.79595790e-04 2.41939163e-01 2.17645058e-03 3.41541473e-03\n",
137 " 1.14507609e-02 1.49639498e-02 2.33483138e-02]\n",
138 "inpure result detected [1. 1. 0.92521572 1. 0.99999627 0.99990293\n",
139 " 0.99946049] [1.59712906e-05 3.60193407e-05 2.48341492e-01 1.21848930e-03\n",
140 " 3.81805601e-03 2.42649228e-02 2.99534246e-02]\n",
141 "inpure result detected [1. 1. 0.99999998 0.99999961 0.93267685 0.99999964\n",
142 " 0.99999822 0.99921774 0.99995547 0.99997301] [5.22490396e-04 3.01556792e-05 4.88795563e-05 6.29703960e-04\n",
143 " 2.34414238e-01 3.72766210e-03 4.72444059e-03 7.62106094e-02\n",
144 " 6.32796684e-02 5.63231562e-02]\n"
145 ]
146 }
147 ],
148 "source": [
149 "plotdata = {}\n",
150 "for file in filenames:\n",
151 " dat = np.load(file, allow_pickle=True)\n",
152 " kx = dat['meta'][()]['k'][0]\n",
153 " radius = dat['meta'][()]['radius']\n",
154 " b = dat['meta'][()]['band_index']\n",
155 " eigvals = dat['eigval']\n",
156 " lMax = dat['meta'][()]['lMax']\n",
157 " residuals = dat['residuals']\n",
158 " ef =dat['empty_freqs']\n",
159 " eigvecs = dat['eigvec']\n",
160 " irweights = []\n",
161 " #for proj in projectors_list[lMax]:\n",
162 " # try:\n",
163 " # irweights.append(np.linalg.norm(np.tensordot(proj, eigvecs, axes=(-1, -1)), axis=0,ord=2) if len(proj) != 0 else np.zeros((len(eigvecs),)))\n",
164 " # except ValueError as err:\n",
165 " # print(proj, len(proj))\n",
166 " # raise err\n",
167 " irweights = np.array(irweights)\n",
168 " #print(irweights)\n",
169 " irweights = np.array([np.linalg.norm(np.tensordot(proj, eigvecs, axes=(-1, -1)), axis=0,ord=2) if len(proj) != 0 else np.zeros((len(eigvecs),)) for proj in projectors_list[lMax]]).T\n",
170 " irclass = np.argmax(irweights, axis=-1)\n",
171 " purities = np.amax(irweights, axis=-1)\n",
172 " if (np.any(purities < 0.98)):\n",
173 " print(\"inpure result detected\", purities, residuals)\n",
174 " #print(purities)\n",
175 " \n",
176 " #for i in range(len(residuals)): \n",
177 " # if residuals[i] < 0.01:\n",
178 " # vec = eigvecs[i]\n",
179 " # for irlabel, proj in projectors.items():\n",
180 " # print(irlabel, np.linalg.norm(np.dot(proj, vec))) #maybe some conj() here?\n",
181 " # print('--->', irlabels[irclass[i]])\n",
182 "\n",
183 " \n",
184 " plotdata[(lMax,radius)] = (eigvals, residuals, b, ef, irclass,)\n"
185 ]
186 },
187 {
188 "cell_type": "code",
189 "execution_count": null,
190 "metadata": {},
191 "outputs": [],
192 "source": [
193 "fig = plt.figure(figsize=(15,6))\n",
194 "axesR = {}\n",
195 "axesI = {}\n",
196 "for i, lMax in enumerate(lMaxes):\n",
197 " axesR[lMax] = fig.add_subplot(2,len(lMaxes),i+1)\n",
198 " axesR[lMax].set_xlim([50,300])\n",
199 " axesR[lMax].set_ylim([1.25, ef[1]/eh])\n",
200 " axesI[lMax] = fig.add_subplot(2,len(lMaxes),len(lMaxes)+i+1)\n",
201 " axesI[lMax].set_xlim([50,300])\n",
202 " axesI[lMax].set_ylim([-60, 30])\n",
203 " axesR[lMax].set_title('$l_\\max = %d $' % lMax) \n",
204 " axesR[lMax].tick_params(labelbottom=False) \n",
205 " if i == len(lMaxes)//2:\n",
206 " axesI[lMax].set_xlabel(\"Particle base radius / nm\")\n",
207 " if i == 0:\n",
208 " axesR[lMax].set_ylabel('$\\hbar \\Re \\omega / \\mathrm{eV}$')\n",
209 " axesI[lMax].set_ylabel('$\\hbar \\Im \\omega / \\mathrm{meV}$')\n",
210 " else:\n",
211 " axesR[lMax].tick_params(labelleft=False) \n",
212 " axesI[lMax].tick_params(labelleft=False) \n",
213 "\n",
214 "res_thr = 0.005\n",
215 "\n",
216 "ir_labeled=set()\n",
217 "if True:\n",
218 " for (lMax, radius), (eigvals, residuals, b, ef, irclass) in plotdata.items():\n",
219 " for i, (e, res, iri) in enumerate(zip(eigvals, residuals, irclass)):\n",
220 " #if i == 0:\n",
221 " if res < res_thr:# and e.real < 2.14e15:\n",
222 " if iri in ir_labeled: \n",
223 " label=None\n",
224 " else:\n",
225 " ir_labeled.add(iri)\n",
226 " label=irlabels[iri]\n",
227 " axesR[lMax].plot(radius*1e9, e.real/eh, \n",
228 " marker='.',\n",
229 " #marker=markerfun(b),\n",
230 " ms=4, #c=colorfun(b)\n",
231 " c=matplotlib.cm.hsv(iri/9),\n",
232 " #c = colorfun(iri),\n",
233 " label=label,\n",
234 " )\n",
235 " axesI[lMax].plot(radius*1e9, e.imag/eh*1000, \n",
236 " #marker='x', \n",
237 " #c=colorfun(b), \n",
238 " c=matplotlib.cm.hsv(iri/9),#colorfun(iri),\n",
239 " marker='.', #markerfun(b),\n",
240 " ms=4,\n",
241 " #label=label\n",
242 " )\n",
243 "fig.legend(title=\"Irrep\", loc=\"center right\")\n",
244 "#fig.suptitle('$l_\\mathrm{max}=%d$, residual threshold = %g' % (lMax, res_thr) )\n",
245 "fig.savefig(plotfilename)\n",
246 "fig.savefig(plotfilename.replace('pdf', 'png'))\n",
247 "print(plotfilename)"
248 ]
249 },
250 {
251 "cell_type": "code",
252 "execution_count": 55,
253 "metadata": {},
254 "outputs": [
255 {
256 "data": {
257 "text/plain": [
258 "array([0. , 1.40635433, 1.98888536, 2.81270865, 3.14470387,\n",
259 " 3.97777072, 4.21906298, 4.44728287])"
260 ]
261 },
262 "execution_count": 55,
263 "metadata": {},
264 "output_type": "execute_result"
265 }
266 ],
267 "source": [
268 "ef / eh"
269 ]
270 },
271 {
272 "cell_type": "code",
273 "execution_count": null,
274 "metadata": {},
275 "outputs": [],
276 "source": []
277 }
278 ],
279 "metadata": {
280 "kernelspec": {
281 "display_name": "Python 3",
282 "language": "python",
283 "name": "python3"
284 },
285 "language_info": {
286 "codemirror_mode": {
287 "name": "ipython",
288 "version": 3
289 },
290 "file_extension": ".py",
291 "mimetype": "text/x-python",
292 "name": "python",
293 "nbconvert_exporter": "python",
294 "pygments_lexer": "ipython3",
295 "version": "3.7.6"
296 }
297 },
298 "nbformat": 4,
299 "nbformat_minor": 4
300 }