Skip to content

Instantly share code, notes, and snippets.

@oroszl

oroszl/sisil.ipynb

Created December 13, 2017 11:55
Show Gist options
  • Save oroszl/55298ba276aaa07cdcbaf2295aa4edbb to your computer and use it in GitHub Desktop.
Save oroszl/55298ba276aaa07cdcbaf2295aa4edbb to your computer and use it in GitHub Desktop.
Download ZIP
bismutahn
Display the source blob
Display the rendered blob
Raw
sisil.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"%pylab inline\n",
"import scipy.linalg as sl\n",
"import sisl\n",
"import tqdm # this is just here for a nice progress bar"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"dat=sisl.get_sile('bi_mono_forwilson_SOC.nc')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"dh=dat.read_hamiltonian()\n",
"ds=dat.read_supercell()\n",
"dg=dat.read_geometry()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
" 0%| | 0/300 [00:00<?, ?it/s]/usr/local/lib/python3.4/dist-packages/scipy/sparse/compressed.py:774: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.\n",
" SparseEfficiencyWarning)\n",
"100%|██████████| 300/300 [00:57<00:00, 5.34it/s]\n"
]
}
],
"source": [
"kran=linspace(0,0.5,300)\n",
"qran=linspace(0,1,7,endpoint=False)\n",
"theta=[]\n",
"for k in tqdm.tqdm(kran):\n",
" W=matrix(eye(2*dh.no))\n",
" for q in qran:\n",
" SK=matrix(sl.fractional_matrix_power(dh.Sk(k=(k,q,0)).toarray(),-0.5))\n",
" HK=SK*matrix(dh.Hk(k=(k,q,0)).toarray())*SK\n",
" w,v=eigh(HK)\n",
" P=zeros_like(HK);\n",
" for i in arange(len(w)):\n",
" if w[i]<0.001:\n",
" P=P+(v[:,i]*v[:,i].H)\n",
" W=W*P\n",
" ev=eigvals(W)\n",
" theta.append((angle(ev[abs(ev)>1e-10])))\n",
"theta=vstack(theta)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAETCAYAAADKy1riAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+Q1Hed5/Hnu5sB3JAsDiGDpeHYi0YLgh4Lq84yXIbV\nNa57BBC39PZ2qfP20mFI9jwkTky2yr2ru5jNGClzlwDT8bxd6qygJcuP0eRQa5kA0tkSSFYcvPPM\nFqJbYkYGTCYkMD/e90dPT3p6ume+3dPd3293vx5VU0n3fOf7/Xzome/r+/58vj/M3REREQkiFnYD\nRESkdig0REQkMIWGiIgEptAQEZHAFBoiIhKYQkNERAJTaIiISGAKDRERCUyhISIigc0KuwHlduON\nN/qSJUtK+tlXX32V6667rrwNijj1uTGoz41hJn0+derUr9x94XTL1V1oLFmyhJMnT5b0s729vbS3\nt5e3QRGnPjcG9bkxzKTPZvbTIMtpeEpERAJTaIiISGAKDRERCUyhISIigSk0REQkMIWGiIgEptDI\n0tfXx8MPP0wqlQq7KSIikVR312mUKpVKsX37doaGhojFYjzxxBMkEomwmyUiEimRDg0zmwscBeYA\ns4GD7v7ZSmyrt7eXa9eu4e6Mjo7S0dEBoOAQEckS9eGpq8Dvuft7gHcDa81sTSU21N7eTiz2xj/H\n6Ogo9957r4aqRESyRDo0PG1w7GUTEAcuVWJbra2tfOpTn5oQHMPDw+zZs6cSmxMRqUmRDg0AM4ub\n2QvAS0Cvu/+wUttat24du3btIh6PA+DuPPnkkySTyUptUkSkppi7h92GQMxsPnAY+Ky7H8n5XgJI\nALS0tKzcu3dvSdsYHBxk3rx57Nixg56enuz1s23bNtatW1dy+6Mq0+dGoj43BvW5OGvXrj3l7qum\nXdDda+YL+BzwmamWWblypZfqyJEj7u5+4sQJnzVrlgPjX7FYzLu7u0ted1Rl+txI1OfGoD4XBzjp\nAfbDkR6eMrOFYxUGZvYm4PeBFyq93dbWVp544olJE+MdHR0aqhKRhhbpU26BtwB/Y2Yx0vMv/8vd\nv1ONDWdOte3o6GB0dBRAp+KKSMOLdGi4+w+AFWFtX8EhIjJRpIenoiCRSLBr1y4NVYmIoNAIRMEh\nIpKm0AhIwSEiotAoioJDRBqdQqNICg4RaWQKjRIoOESkUSk0SqTgEJFGFOnrNKJuqus4XnzxRebP\nn097ezutra1hNlNEpGwUGjNUKDi6urowM+LxuJ4CKCJ1Q8NTZZBvqArSN4McHh7WkJWI1A2FRplk\ngqOpqQkzm/C90dFRtmzZwsaNG/UkQBGpaQqNMkokEjz77LM89NBDdHZ2Tqg83J0DBw6wZs0aVR0i\nUrM0p1Fmra2t4xPft9xyy4S5DoCRkRG2bNnCM888Q2dnpybJRaSmKDQqKDP5vXXrVkZGRsbfz1Qd\nPT097Ny5U5PkDSSVSrFnzx4uXLhQ1e329/ezcOHCCe8tWrSIFStW8Pzzz5elPYsWLWLz5s06EKpz\nCo0KSyQSLF++nK6uLg4ePJh5AiGgqqOWlbLzHxgY4Pjx4xMqz3qTTCZpa2ujubm57OtWKEWDQqMK\nWltb2b9/P8lksmDVcejQIe68806FR8gyYQAUPApvhJ1/qUZHRzl69GjF1l9MKPX397N8+fJpqymF\nUXEUGlU0VdUxOjo6Hh5tbW0sXbpUv8hllgmEM2fOTBqqAYVBLSg2lL73ve8FWq5SFVI9BpJl77jq\nwapVq/zkyZMl/Wxvby/t7e3lbVAB+aqOXLFYrOLVRzX7XA2Fho2iFAixWGxGO6iBgQHOnz/P66+/\nHvhnrl69ypw5cya819zczIoVKzh9+jSXLl0quh1z585l8eLFNDc3R+rfN2pm+nnnyv38sz/HK1eu\n8OCDD5Y0T2pmp9x91bTLKTTeUO0daCqVoquri0OHDk35x5b5patE9VGLoRGlYFiwYAFNTU2Bl8/e\n0ZYiajvn7B1iKWEW1LVr1xgYGCj7eutVd3d30cGh0ChBWDvQoOEBE/9Iy1H6RjU0wg6GJUuWsHr1\n6rxH4UNDQ1y8eLGi2xeZiQ996EMcPny4qJ8JGhqa04iAzER5Zkd59uzZgjvG3DHd3LHYWhtDzRcO\nUTiaPnfuHOfOnQtt+yIzsWnTpoqtW6ERIdkXBgatPvJNDBaa1As7UHIDIgrhUKuKHScv93Ua+uzC\nlxkaLdecRlAansoSxaGaINVHMXJ3Nvl2JkFMtcPJN7Zdj2PS5Z7gDKqU8K/E73ZYFyoGFfSU21oz\n1ec/k8+5LoanzOxmYA/QAjiQdPfHwm1VdeVWH5k/0lKP9Cp9Hn2tmG6HP91ReNhVWxRk/25GURQP\nAutBpEMDGAa2u/tpM7seOGVm33H3s2E3LAy5f6Qa7kkr9oh/ZGSEBx54INI7PJGoinRouPsvgF+M\n/f8rZvYj4K1AQ4ZGrnxHemGfdTRTxQZAqUM1CgyR0tTMnIaZLQGOAre5+8s530sACYCWlpaVe/fu\nLWkbg4ODzJs3b2YNjbC+vj4OHz48YW5hZGSEeDwOwMsvv0x/fz9Xr15leHiYV155ZUbbu/7665k1\n643jktmzZ3PTTTdxww035F2+ubmZO+64g2XLls1ou9Op9885H/W5Mcykz2vXrq2f6zTMbB7wLPCQ\nu//tVMvW20R4JaVSKR5++GHi8fiMK5FyXz9SSY32OYP63CgafiIcwMyagH3AV6cLDJlaOSbSs69O\nz0wUA5EOCREpn0iHhqWfm/o/gB+5+46w21OLZnrKbi1VECJSeZEODWA18KfAGTN7Yey9B9396RDb\nFHkzCQqFhIhMJdKh4e7HAQu7HbWg1KAwM9asWaOQEJFAIh0aMr1ibnYIkyuJ2267jXvuuacKLRWR\neqDQqEHFVhVT3Vq9t7e3wq0VkXqi0KghpdxCXU8AFJFyUmjUgCg8rElEBBQakVZMWFT6sbAiIqDQ\niKQgYaGqQkTCoNCIkKBhoapCRMKi0IiIZDLJ1q1bGRkZyft9hYWIRIFCI2SZ6uLgwYPku3mkwkJE\nokShEZLphqIUFiISRQqNEEw1FGVmrF+/XmEhIpGk0KiyZDJJR0dH3uoiHo+zc+dOEolECC0TEZle\nLOwGNIpUKsXGjRvZsmXLpMCIxWJs2LCBY8eOKTBEJNJUaVRBoeEoDUWJSK1RaFRYoeGoWCzGrl27\nVFmISE3R8FQFFQqMeDyuwBCRmqRKo8xSqRS9vb1cvnyZRx99dEJgaDhKRGqdQqOMkskk9957L8PD\nw5Mu1NNwlIjUA4VGmUx1Kq0CQ0TqheY0yqBQYJgZs2bNUmCISN1QpTFD+QIjFotx3333MX/+fNrb\n2zV/ISJ1Q6FRokI3GtRQlIjUs8iHhpl9BfhXwEvuflvY7YHCF+spMESk3tXCnMZfAx8OuxEZmeEo\nBYaINKLIh4a7HwUGwm4H6GI9EZHIh0ZU5AsMM9ONBkWkoVi+p8VFjZktAb5ZaE7DzBJAAqClpWXl\n3r17S9rO4OAg8+bNm/BeX18fTz31FCdOnJgw4W1mbNu2jXXr1pW0rajI1+d6pz43BvW5OGvXrj3l\n7qumXdDdI/8FLAF+GGTZlStXeqmOHDky4XV3d7fH43EHJnzFYjHv7u4ueTtRktvnRqA+Nwb1uTjA\nSQ+wj4382VNh0d1pRUQmi/ychpk9BaSAd5rZz83szyq9TU14i4jkF/lKw93/dbW21dfXx2OPPTbp\ngj3dnVZEJC3yoVEt999/P1/4whd0d1oRkSkoNEgPR3V1dU16X4EhIjJR5Oc0qmHfvn2T3jMzBYaI\nSA6FBrBp06ZJ77k7zzzzDKlUKoQWiYhEk0IDSCQSdHZ2Tnr/wIEDrFmzhmQyGUKrRESiR6Ex5pFH\nHuHTn/40sdjEf5KRkRG2bNnCxo0bVXWISMNTaGRZt24du3btIh6PT3jf3VV1iIig0JgkkUhw7Ngx\nNmzYgJlN+J6qDhFpdAqNPFpbW9m/fz+7d+9W1SEikkWhMQVVHSIiEyk0phGk6mhra1N4iEhDUGgE\nNFXVMTo6qiErEWkICo0iTFV1gIasRKT+KTRKkF115F7XkRmyuv3229m4cSMbN26ko6NDISIidUGh\nUaJM1XH8+PG8Q1ZDQ0McOHCAAwcOsHv3bs17iEhdUGjM0HRDVhmZeY+2tjZuv/12VR8iUpMUGmWS\nGbLasmULGzZsoKmpKe9yo6OjHD16dLz6UICISC3R8zTKqLW1dfzJfqlUij179nD27FmOHz8+6dGx\n8EaAHD16lGQySVtbG83NzSxatIjNmzfrKYEiJUqlUuzYsYOvfe1rrFixgueff54LFy6E3axABgYG\nOH/+PK+//nqg5Zubm8f7+Oqrr/Lggw9W9JEOCo0KyQ2Qrq4uDh06lDc84I0AyVCISLlkdqCPPfbY\nhPcXLVoU6R1qsTvPjGvXrjEwMFChVkXPhQsXOHv27Pjru+++G6BiwaHQqILMvEeQ6iNjqhDJUJiE\nI/M5RnFHm2tgYGDa3zWpP/v27VNo1INih6+y5YZIRr4wyZYbLLk7vP7+fhYuXFh0X6J+lDqVIH3O\n17+BgQF+/OMf11x/pfHke7BcuZi7V2zlYVi1apWfPHmy6J9LJpN8/vOf5+rVqzNuQ2aM8fTp01y6\ndGna5YeGhrhy5QojIyO4O0NDQzNuQ66mpibcneHh4bKvW0TKb8GCBQVPqMlWrjkNMzvl7qumW27a\nSsPMbgLWAIuAa8DPgLPufr7oVkVUMpkcHwcsh9wxxiioRBCJVEPQnWe2uXPncv3117N69eqaq4hn\nMuzc29tLe3t7+RuVpWBomNks4L8D/548p+aa2c+AbwH/092LP7QPyMw+DDwGxIEvu/tflXsb+/bt\nK/cqRSJvyZIlrF69OnBFXAlmRqHRjrlz57J48eKSTwapxg60EU1Vafxn4G7gPLAfuAjMBT4IvBd4\nC9ABbDGzp4F7yl19mFkceAL4feDnwPfN7JC7l/UwftOmTXz7298u5ypFIu/cuXOcO3cu7GZMKbt9\n083f5erv72f58uU1V2kUq9onxBSc0zCz88Cvgfe7+6tZ7/8l8DngN4E7gE8CfwBcAv7Q3f++bI0z\nawX+k7vfMfb6AQB3f7jQz9TinEYlZc+XAMTjcd70pjcxe/ZsAK5evcqcOXOKXm+U+lisIH3O17/s\nI99i6SwmqaRYLMbNN98MEOqcxkLgq9mBkc3dB4F9wD4zez/wFPAtM1vu7r8ousX5vZX0HErGz4H3\nlWndEyQSCW699daGK2cbsYQPq89hnqqbOWNM4VWfRkdH+elPfwqEe53GOWBJkJW4+3Nm1g68QLoK\n6Zhpw4phZgkgAdDS0kJvb29J6xkcHCz5Z2uV+lxdH//4x+nr6+Pw4cNVvQBtZGSE/v5+AJYtW0Z/\nfz/uzksvvVRwTkFq15NPPsmtt95akXVPFRpPAQ+a2fvd/bnpVuTuPzWzrwN/WLbWwT8BN2e9ftvY\ne7nbTgJJSA9PlXIUWc/DU5nhmHxnoJRjeOry5cs1teMpdnjqtddeK3lYKpeO9KUa7rrrropV01PN\nabwJOEl6p90JPOnuI5k5DXefdEtXM9sBbHH33yhL49JncP0Y+ADpsPg+8Mfu3lfoZ0qZ0yj3Kbci\nUjmxWCzQhHgjTYTfcMMN9PT0cOXKlfDmNNz9NTP7APBN0mcw/YWZPUX6aD/fBm8D/gR4sejWFm7D\nsJndCxwmfcrtV6YKjFLplFtpVPF4nJtvvrno+ztV0lQnHBRzplCjzdc98sgj4V6nAeDuF8bOYPoM\ncN/YlwOY2Q9IH/1fJX3h30rSO/bt5Wyguz8NPF3OdebSKbfSqEZGRiJ52u25c+cKXtR34MCBSe/l\nC5paqTRq7R5ygW8jYmZzgY8CHwHagMU5i/wI+K/u/lRZW1ikWjzlthq3ERGRaCvlyneYuL+pxvBU\nyfeeMrPfIH1K7GzggrtfLGlFZVZqaEB1y9lib1hYyHTju9PdWFA3LMwvu3/PPfccL7zwQpVaJzJz\n3d3dRQdH2e49VYi7XwH+X6k/34hmEhT5wqEcZW2jjftCaX2upduhZwwMDHDs2LGaOrNNykO3Rq9x\nQR7ClCs7JGptzLMeZd/WvpZkwu7MmTOTqquoVovFPnxpaGiIixcjMdARGZW8NbpCo4KKCQuFhFRC\nJuzqvaLMVwnWykR4qU8ozCjXnEZQCo0yS6VS9Pb2cvnyZb74xS+O3/Mpn0xQLF26VCEhMgP5KsF6\nD8p8Qj/lVoLLVBU9PT2Mjo4WHEdWUIhILVNozFDQIahYLMadd95JZ2engkJEapZCo0RBwsLMiMVi\nrFu3TmEhInVBoVGCZDLJ1q1bC85XxONxtm/fzvz582lvb1dYiEjdUGgUIVNdHDx4MO+chYagRKTe\nKTQCmqq6UFiISKNQaExjqurCzFi/fr3CQkQahkJjClNVF/F4nJ07d1b0IhoRkaiJhd2AqEomk3R0\ndEwKDDNjw4YNHDt2TIEhIg1HlUaOqYajVF2ISKNTaGTp6enhsccey1tdaO5CREShMe7+++9nx44d\nk96PxWLs2rVL1YWICJrTANLzF11dXZPej8fjCgwRkSwKDdIPLMllZmzfvl2BISKSRaFB/geWuDuP\nPvooyWQyhBaJiESTQgNIJBJ0dnZOen90dJSOjg4Fh4jIGIXGmEceeYTHH3+cDRs2YGbj74+OjrJl\nyxY2btxIKpUKsYUiIuGLbGiY2R+ZWZ+ZjZrZqmpsc9myZezfv5/du3cTi73xT+PuHDhwgDVr1qjq\nEJGGFtnQAH4IfBQ4Wu0NJxIJdu3aNSE4AEZGRjRcJSINLbLXabj7j4AJQ0XVlDlrKvfeU5l5juxl\nREQahRV6lnVUmFkvcJ+7n5ximQSQAGhpaVm5d+/ekrY1ODjIvHnzJrzX19fHU089xYkTJybcVsTM\n2LZtG+vWrStpW1GRr8/1Tn1uDOpzcdauXXvK3aefCnD30L6A75Iehsr9Wp+1TC+wKug6V65c6aU6\ncuRIwe91d3d7LBZzYPzLzHzDhg1+4sSJkrcZtqn6XK/U58agPhcHOOkB9rGhzmm4+wfd/bY8XwfD\nbFc++eY5XBPkItJgojwRHjmaIBeRRhfZ0DCzjWb2c6AV+JaZHQ67TfBGcMTj8Qnv60JAEWkEkQ0N\nd9/v7m9z9znu3uLud4TdpoxEIsGxY8fyXgio4BCRehbZU26jrrW1lf37948/4W90dBR4IzhefPFF\n5s+fT3t7u57BISJ1Q6ExQ5lrNXKDo6urCzMjHo/zxBNP6JoOEakLkR2eqiWFJsjdneHhYQ1ZiUjd\nUGiUSSY4mpqaJl3FrrkOEakXCo0ySiQSPPvsszz00EN0dnZOqDx0t1wRqQea0yiz1tbW8YnvW265\nZcJcR+ZiwJ6eHnbu3Kl5DhGpOao0KkgXA4pIvVFoVNhUFwNquEpEao1CowoKXQyoe1eJSK1RaFRJ\n5mLA3KcCQnq4SlWHiNQChUaVFRquylQdbW1tCg8RiSyFRggKDVdBeq5D4SEiUaXQCEn2cFVu1QEK\nDxGJJoVGyLKrjty5DlB4iEi0KDQiIFN1HD9+XOEhIpGm0IgQhYeIRJ1CI4KKDY/bb7+djo4OBYiI\nVJxCI8KChsfRo0fZvXu3AkREKk6hUQOChAcoQESk8hQaNSRoeIACREQqQ7dGr0GZ8EilUuzZs4ez\nZ89y/Pjx8Vuw58oEyNGjR0kmk7S1tdHc3MyiRYu47bbbaG9vr24HRKRmKTRqWPazO4oNkAwz4+tf\n//p4iGzevHl8nSIiuSIbGmb2BWAdcA14Efiku18Ot1XRVUqAQPqeV9khkluJKEREJFtkQwP4DvCA\nuw+b2SPAA8D9IbepJpQaIDC5ElGIiEi2yIaGu3876+VzwMfCakstyxcgFy5cYGBgYEYhsnTpUlas\nWMHFixdpb29XkIg0iMiGRo5/B3wt7EbUuuwAgXSIPPzww8Tj8aJDJBMkZkYsFmP16tU0NzcDqCIR\nqWPm7uFt3Oy7wKI83/oLdz84tsxfAKuAj3qBxppZAkgAtLS0rNy7d29J7RkcHGTevHkl/Wwt6Ovr\n4/DhwwwMDIy/NzIyMn6X3Zdffpn+/n6uXr3K8PAwr7zyyoy2d/311zNrVvq4ZPbs2dx0003ccMMN\neZdtbm7mjjvuYNmyZTPaZhD1/jnnoz43hpn0ee3atafcfdV0y4UaGtMxs38L3A18wN2vBPmZVatW\n+cmTJ0vaXm9vb82ffpo9BJUtaCURplgsNj5/EkSpFU09fM7FUp8bw0z6bGaBQiOyw1Nm9mGgE7g9\naGA0mtyAqIVgmEru/EkQ2RP1QejaFJGZiWxoAI8Dc4DvjD3d7jl33xJuk8JX7NlQ9a6UoMm+NqWQ\nRYsWsWLFCp5//vlJVZvmbKSRRTY03P3tYbchCko542kquUNA/f39LFy4sOj1TLVTHRgY4Pz587z+\n+usAXLt2bcI8Sthyr00pRbEVzkzVQlAVGhoNS39/P8uXLy/4e1provI7EOk5jVLUy5xGKpWiq6uL\nQ4cOFR0SheYG8v3SVavP+XYotT6cVm3Fzvlky3dwMFXwF0ufZfUsWLCApqam8dfNzc2sWLGC06dP\nc+XKFR588EESiUTR6635OY1GVMrQU+6OJCpHI7lyT/fNUJgEV8pQnNSfixcvTnh94cIFzp49O/76\n7rvvBigpOIJQaERAsVVFJiiWLl0ayYAoRjFhAtUJlAULFtDS0jJ+9Hbp0qXx7w0NDU36oxWJmn37\n9ik06lHQsMiuJqJaSZRboTCBylcnFy9e5OLFixOO3grJHSqYyty5c1m8eHHRw0u1VHmZGS0tLRVZ\nd9TmxqJs06ZNFVu3QiMkyWSSrVu3MjIyUnCZWCzGnXfeSWdnZ92HRDGKrU6gcjveYquO8+fPFz0v\n0dzcTFtb24STC4p19epV5syZM2m9+aqpmfjlL39Z8xPOlTCT+aiM3BNMMso1pxGUQqPKMtXFwYMH\nyXcSQj0NPVXbVNUJvBEqZ86cmTQpXK2j+SjNS+SOhTe6cuzY86nm6EA1TmxRaFTRVNWFqorKy4RK\noT+s7EplulOKa2W4qB4F3bkXc8ptowz7loNCowqmqi7MjPXr1yssImC6SiVbsdckNFLQVOqIHYrb\nuUfpFPp6otCosKmqi3g8zs6dOys6/iiVUUzAZIR58Vulr9PIXqeO2OubQqOCkskkHR0dk44uVV00\nplKCplx01C3lotAos1QqRW9vL5cvX+bRRx+dFBiqLkSklik0yiiZTHLvvfcyPDysuQsRqUsKjTIp\nNBQF6YnBXbt2qboQkZoXC7sB9WCquYtZs2YpMESkbqjSmKF8gRGLxbjvvvuYP38+7e3tGo4Skbqh\n0JiBQoGhykJE6pWGp0qkwBCRRqTQKIECQ0QalUKjSAoMEWlkCo0iKDBEpNEpNAJSYIiIKDQCUWCI\niKQpNKahwBAReUNkQ8PM/ouZ/cDM/sHM/s7MFle7DQoMEZGJIhsawBfc/d3u/h7gAPCX1dy4AkNE\nZLLIhoa7v5z18jrgYrW2nUqluOeeexQYIiI5LPcW3lFiZg8Bm4HXgPe5+6UCyyWABEBLS8vKvXv3\nlrS9wcFB5s2bx44dO+jp6cleP9u2bWPdunUlrTfKMn1uJOpzY1Cfi7N27dpT7r5q2gXdPbQv4LvA\nD/N8rc9Z7gHgr4Osc+XKlV6qI0eOeHd3t8fjcQcc8Hg87t3d3SWvM+qOHDkSdhOqTn1uDOpzcYCT\nHmAfG+oNC939gwEX/SrwTCXbAtDT08OXvvSl8WEpM+Ouu+7SkJSIyJjIzmmY2TuyXq4HXqjk9u6/\n/3527NgxYR5j1qxZbN68uZKbFRGpKVG+Nfpfmdk7gRHgH4GOSm0omUzS1dU14b1YLMbjjz+uZ2GI\niGSJbGi4+6ZqbWvfvn0TXpuZzpQSEckjssNT1bRp08R8+sxnPqPAEBHJI7KVRjVlAuLJJ5/UxLeI\nyBQUGmMSiQS33nor7e3tYTdFRCSyNDwlIiKBKTRERCQwhYaIiASm0BARkcAUGiIiEphCQ0REAov0\nrdFLYWb9wE9L/PEbgV+VsTm1QH1uDOpzY5hJn/+Zuy+cbqG6C42ZMLOTHuR+8nVEfW4M6nNjqEaf\nNTwlIiKBKTRERCQwhcZEybAbEAL1uTGoz42h4n3WnIaIiASmSkNERAJryNAwsw+b2f81s5+Y2Wfz\nfN/M7L+Nff8HZvbbYbSznAL0+V1mljKzq2Z2XxhtLLcAff43Y5/vGTM7YWbvCaOd5RKgv+vH+vuC\nmZ02sw+E0c5ymq7PWcv9jpkNm9nHqtm+SgjwObeb2a/HPucXzOxzZW2AuzfUFxAHXgT+OTAb+Adg\nac4yHwGeAQx4P/D3Ybe7Cn2+Cfgd4CHgvrDbXKU+/y7w5rH//4Na/pwD9ncebwxJvxt4Mex2V7rP\nWcv9HfA08LGw212Fz7kd+Gal2tCIlcZ7gZ+4+z+6+zVgL7A+Z5n1wB5Pew6Yb2ZvqXZDy2jaPrv7\nS+7+fWAojAZWQJA+n3D3S2MvnwPeVuU2llOQ/g762F4FuA64WOU2lluQv2WAPwf2AS9Vs3EVErTP\nFdOIofFW4GdZr38+9l6xy9SSeutPEMX2+c9IV5e1KlB/zWyjmf0f4H8D/6FKbauUaftsZm8FNgK7\nqtiuSgr6e/27Y0ORz5jZsnI2QE/uk4ZnZmtJh0Zb2G2pNHffD+w3s38J7DGzd7n7aNjtqqAvAfe7\n+6iZhd2WajkNLHb3QTP7CHAAeEe5Vt6IlcY/ATdnvX7b2HvFLlNL6q0/QQTqs5m9G/gysN7da3m4\npqjP2N2Pkj5oXFDhdlVSkD6vAvaa2TngY8BOM9tQneZVxLR9dveX3X1w7P+fBprM7MZyNaARQ+P7\nwDvM7LfMbDbwCeBQzjKHgM1jZ1G9H/i1u/+i2g0toyB9rjfT9tnMFgN/C/ypu/84hDaWU5D+vt3G\nDrfHzgg0d++vflPLZto+u/tvufsSd18CfAPY6u4Hqt/UsgnyOS/K+pzfS3o/X7YDooYbnnL3YTO7\nFzhM+kwkLFhMAAACKklEQVSEr7h7n5ltGfv+btJnWXwE+AlwBfhkWO0thyB9NrNFwEngBmDUzP4j\n6bMyXg6t4TMQ8HP+HOkj7Z1jf2PDXqM3uAvY302kD4aGgFdJ73BqVsA+15WAff4Y0GFmw8BrwCey\nToCYMV0RLiIigTXi8JSIiJRIoSEiIoEpNEREJDCFhoiIBKbQEBGRwBQaIhVkZtvMzM3sj8Nui0g5\nKDREKmvl2H9PhtoKkTLRdRoiFWRmZ0nfUG5+OS+wEgmLKg2RCjGz64B3As8rMKReKDREKudfkP4b\nO5X9ppm92cwOjs11fMnMmsJpnkjxGu7eUyJVlHlM8HhomNn7gK8Bbwb+yN2/EUbDREqlSkOkcjKT\n4KcAzOzTwDHgZWCVAkNqkSoNkcr5beAV4FdmdhC4E/gboMPdXwu1ZSIl0tlTIhVgZnNJB8ZLwDBw\nE/Dn7v7lUBsmMkManhKpjPeQruTnAIuBbygwpB4oNEQqIzMJvh34JvAnYw+2EqlpCg2Rysi+EvwT\npCfDv2hmHw2vSSIzpzkNkQows9PAu4Dr3X1k7HG6z5Ge2/g9d38u1AaKlEiVhkiZmdls4DbgjLuP\nALj7BdLPnb8K9JjZ20NsokjJFBoi5bccaAKez37T3c8CHwV+E3jGzG4MoW0iM6LhKRERCUyVhoiI\nBKbQEBGRwBQaIiISmEJDREQCU2iIiEhgCg0REQlMoSEiIoEpNEREJDCFhoiIBKbQEBGRwP4/f73k\nhpQAODwAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb88b171470>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot(kran,theta,'k.');\n",
"grid()\n",
"ylabel(r'$\\theta$',fontsize=20)\n",
"xlabel(r'$k$',fontsize=20);"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def Skd(dh, k=(0, 0, 0)):\n",
" \n",
" dtype = np.complex128 \n",
" k = np.asarray(k, np.float64)\n",
" k.shape = (-1,)\n",
" # Setup the quantity for this k-point\n",
" Vf = dh.tocsr(dh.S_idx).todense()\n",
" # number of orbitals\n",
" no = dh.no\n",
" V = zeros((no, no), dtype=dtype)\n",
" # Get the reciprocal lattice vectors dotted with k\n",
" kr = dot(dh.rcell, k)\n",
" for si, isc in dh.sc:\n",
" phase = exp(-1j * dot(kr, dot(dh.cell, isc)))\n",
" V += Vf[:, si * no:(si + 1) * no] * phase\n",
"\n",
" del Vf\n",
" return V\n",
"\n",
"def Hkd(dh, k=(0, 0, 0)):\n",
" \n",
" k = np.asarray(k, np.float64)\n",
" k.shape = (-1,)\n",
"\n",
" # number of orbitals\n",
" no = dh.no\n",
" dtype = np.complex128\n",
" # sparse matrix dimension (2 * self.no)\n",
" v = [dh.tocsr(i).todense() for i in range(len(dh.spin))]\n",
" V11 = zeros((no, no), dtype=dtype)\n",
" V22 = zeros_like(V11)\n",
" V21 = zeros_like(V11)\n",
" V12 = zeros_like(V11)\n",
"\n",
" # Get the reciprocal lattice vectors dotted with k\n",
" kr = dot(dh.rcell, k)\n",
" for si, isc in dh.sc:\n",
" phase = exp(-1j * dot(kr, dot(dh.cell, isc)))\n",
" sl = slice(si*no, (si+1) * no, None)\n",
"\n",
" # diagonal elements\n",
" V11 += (v[0][:, sl] + 1j * v[4][:, sl]) * phase\n",
" V22 += (v[1][:, sl] + 1j * v[5][:, sl]) * phase\n",
"\n",
" # lower off-diagonal elements\n",
" V21 += (v[2][:, sl] - 1j * v[3][:, sl]) * phase\n",
"\n",
" # upper off-diagonal elements\n",
" V12 += (v[6][:, sl] + 1j * v[7][:, sl]) * phase\n",
"\n",
" del v\n",
"\n",
" V = zeros((len(dh), len(dh)), dtype=dtype)\n",
" V[::2, ::2] = V11\n",
" V[1::2, 1::2] = V22\n",
" V[1::2, ::2] = V21\n",
" V[::2, 1::2] = V12\n",
"\n",
" del V11, V22, V21, V12\n",
"\n",
" return V\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 300/300 [00:27<00:00, 10.93it/s]\n"
]
}
],
"source": [
"kran=linspace(0,0.5,300)\n",
"qran=linspace(0,1,7,endpoint=False)\n",
"theta=[]\n",
"for k in tqdm.tqdm(kran):\n",
" W=matrix(eye(2*dh.no))\n",
" for q in qran:\n",
" SK=matrix(sl.fractional_matrix_power(Skd(dh,k=(k,q,0)),-0.5))\n",
" SK=kron(SK,eye(2))\n",
" HK=SK*matrix(Hkd(dh,k=(k,q,0)))*SK\n",
" w,v=eigh(HK)\n",
" P=zeros_like(HK);\n",
" for i in arange(len(w)):\n",
" if w[i]<0.001:\n",
" P=P+(v[:,i]*v[:,i].H)\n",
" W=W*P\n",
" ev=eigvals(W)\n",
" theta.append((angle(ev[abs(ev)>1e-10])))\n",
"theta=vstack(theta)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAETCAYAAADKy1riAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+Q1Hed5/Hnu5sB3JAsDiGDpeHYi0YLgh4Lq84yXIbV\nNa57BBC39PZ2qfP20mFI9jwkTky2yr2ru5jNGClzlwDT8bxd6qygJcuP0eRQa5kA0tkSSFYcvPPM\nFqJbYkYGTCYkMD/e90dPT3p6ume+3dPd3293vx5VU0n3fOf7/Xzome/r+/58vj/M3REREQkiFnYD\nRESkdig0REQkMIWGiIgEptAQEZHAFBoiIhKYQkNERAJTaIiISGAKDRERCUyhISIigc0KuwHlduON\nN/qSJUtK+tlXX32V6667rrwNijj1uTGoz41hJn0+derUr9x94XTL1V1oLFmyhJMnT5b0s729vbS3\nt5e3QRGnPjcG9bkxzKTPZvbTIMtpeEpERAJTaIiISGAKDRERCUyhISIigSk0REQkMIWGiIgEptDI\n0tfXx8MPP0wqlQq7KSIikVR312mUKpVKsX37doaGhojFYjzxxBMkEomwmyUiEimRDg0zmwscBeYA\ns4GD7v7ZSmyrt7eXa9eu4e6Mjo7S0dEBoOAQEckS9eGpq8Dvuft7gHcDa81sTSU21N7eTiz2xj/H\n6Ogo9957r4aqRESyRDo0PG1w7GUTEAcuVWJbra2tfOpTn5oQHMPDw+zZs6cSmxMRqUmRDg0AM4ub\n2QvAS0Cvu/+wUttat24du3btIh6PA+DuPPnkkySTyUptUkSkppi7h92GQMxsPnAY+Ky7H8n5XgJI\nALS0tKzcu3dvSdsYHBxk3rx57Nixg56enuz1s23bNtatW1dy+6Mq0+dGoj43BvW5OGvXrj3l7qum\nXdDda+YL+BzwmamWWblypZfqyJEj7u5+4sQJnzVrlgPjX7FYzLu7u0ted1Rl+txI1OfGoD4XBzjp\nAfbDkR6eMrOFYxUGZvYm4PeBFyq93dbWVp544olJE+MdHR0aqhKRhhbpU26BtwB/Y2Yx0vMv/8vd\nv1ONDWdOte3o6GB0dBRAp+KKSMOLdGi4+w+AFWFtX8EhIjJRpIenoiCRSLBr1y4NVYmIoNAIRMEh\nIpKm0AhIwSEiotAoioJDRBqdQqNICg4RaWQKjRIoOESkUSk0SqTgEJFGFOnrNKJuqus4XnzxRebP\nn097ezutra1hNlNEpGwUGjNUKDi6urowM+LxuJ4CKCJ1Q8NTZZBvqArSN4McHh7WkJWI1A2FRplk\ngqOpqQkzm/C90dFRtmzZwsaNG/UkQBGpaQqNMkokEjz77LM89NBDdHZ2Tqg83J0DBw6wZs0aVR0i\nUrM0p1Fmra2t4xPft9xyy4S5DoCRkRG2bNnCM888Q2dnpybJRaSmKDQqKDP5vXXrVkZGRsbfz1Qd\nPT097Ny5U5PkDSSVSrFnzx4uXLhQ1e329/ezcOHCCe8tWrSIFStW8Pzzz5elPYsWLWLz5s06EKpz\nCo0KSyQSLF++nK6uLg4ePJh5AiGgqqOWlbLzHxgY4Pjx4xMqz3qTTCZpa2ujubm57OtWKEWDQqMK\nWltb2b9/P8lksmDVcejQIe68806FR8gyYQAUPApvhJ1/qUZHRzl69GjF1l9MKPX397N8+fJpqymF\nUXEUGlU0VdUxOjo6Hh5tbW0sXbpUv8hllgmEM2fOTBqqAYVBLSg2lL73ve8FWq5SFVI9BpJl77jq\nwapVq/zkyZMl/Wxvby/t7e3lbVAB+aqOXLFYrOLVRzX7XA2Fho2iFAixWGxGO6iBgQHOnz/P66+/\nHvhnrl69ypw5cya819zczIoVKzh9+jSXLl0quh1z585l8eLFNDc3R+rfN2pm+nnnyv38sz/HK1eu\n8OCDD5Y0T2pmp9x91bTLKTTeUO0daCqVoquri0OHDk35x5b5patE9VGLoRGlYFiwYAFNTU2Bl8/e\n0ZYiajvn7B1iKWEW1LVr1xgYGCj7eutVd3d30cGh0ChBWDvQoOEBE/9Iy1H6RjU0wg6GJUuWsHr1\n6rxH4UNDQ1y8eLGi2xeZiQ996EMcPny4qJ8JGhqa04iAzER5Zkd59uzZgjvG3DHd3LHYWhtDzRcO\nUTiaPnfuHOfOnQtt+yIzsWnTpoqtW6ERIdkXBgatPvJNDBaa1As7UHIDIgrhUKuKHScv93Ua+uzC\nlxkaLdecRlAansoSxaGaINVHMXJ3Nvl2JkFMtcPJN7Zdj2PS5Z7gDKqU8K/E73ZYFyoGFfSU21oz\n1ec/k8+5LoanzOxmYA/QAjiQdPfHwm1VdeVWH5k/0lKP9Cp9Hn2tmG6HP91ReNhVWxRk/25GURQP\nAutBpEMDGAa2u/tpM7seOGVm33H3s2E3LAy5f6Qa7kkr9oh/ZGSEBx54INI7PJGoinRouPsvgF+M\n/f8rZvYj4K1AQ4ZGrnxHemGfdTRTxQZAqUM1CgyR0tTMnIaZLQGOAre5+8s530sACYCWlpaVe/fu\nLWkbg4ODzJs3b2YNjbC+vj4OHz48YW5hZGSEeDwOwMsvv0x/fz9Xr15leHiYV155ZUbbu/7665k1\n643jktmzZ3PTTTdxww035F2+ubmZO+64g2XLls1ou9Op9885H/W5Mcykz2vXrq2f6zTMbB7wLPCQ\nu//tVMvW20R4JaVSKR5++GHi8fiMK5FyXz9SSY32OYP63CgafiIcwMyagH3AV6cLDJlaOSbSs69O\nz0wUA5EOCREpn0iHhqWfm/o/gB+5+46w21OLZnrKbi1VECJSeZEODWA18KfAGTN7Yey9B9396RDb\nFHkzCQqFhIhMJdKh4e7HAQu7HbWg1KAwM9asWaOQEJFAIh0aMr1ibnYIkyuJ2267jXvuuacKLRWR\neqDQqEHFVhVT3Vq9t7e3wq0VkXqi0KghpdxCXU8AFJFyUmjUgCg8rElEBBQakVZMWFT6sbAiIqDQ\niKQgYaGqQkTCoNCIkKBhoapCRMKi0IiIZDLJ1q1bGRkZyft9hYWIRIFCI2SZ6uLgwYPku3mkwkJE\nokShEZLphqIUFiISRQqNEEw1FGVmrF+/XmEhIpGk0KiyZDJJR0dH3uoiHo+zc+dOEolECC0TEZle\nLOwGNIpUKsXGjRvZsmXLpMCIxWJs2LCBY8eOKTBEJNJUaVRBoeEoDUWJSK1RaFRYoeGoWCzGrl27\nVFmISE3R8FQFFQqMeDyuwBCRmqRKo8xSqRS9vb1cvnyZRx99dEJgaDhKRGqdQqOMkskk9957L8PD\nw5Mu1NNwlIjUA4VGmUx1Kq0CQ0TqheY0yqBQYJgZs2bNUmCISN1QpTFD+QIjFotx3333MX/+fNrb\n2zV/ISJ1Q6FRokI3GtRQlIjUs8iHhpl9BfhXwEvuflvY7YHCF+spMESk3tXCnMZfAx8OuxEZmeEo\nBYaINKLIh4a7HwUGwm4H6GI9EZHIh0ZU5AsMM9ONBkWkoVi+p8VFjZktAb5ZaE7DzBJAAqClpWXl\n3r17S9rO4OAg8+bNm/BeX18fTz31FCdOnJgw4W1mbNu2jXXr1pW0rajI1+d6pz43BvW5OGvXrj3l\n7qumXdDdI/8FLAF+GGTZlStXeqmOHDky4XV3d7fH43EHJnzFYjHv7u4ueTtRktvnRqA+Nwb1uTjA\nSQ+wj4382VNh0d1pRUQmi/ychpk9BaSAd5rZz83szyq9TU14i4jkF/lKw93/dbW21dfXx2OPPTbp\ngj3dnVZEJC3yoVEt999/P1/4whd0d1oRkSkoNEgPR3V1dU16X4EhIjJR5Oc0qmHfvn2T3jMzBYaI\nSA6FBrBp06ZJ77k7zzzzDKlUKoQWiYhEk0IDSCQSdHZ2Tnr/wIEDrFmzhmQyGUKrRESiR6Ex5pFH\nHuHTn/40sdjEf5KRkRG2bNnCxo0bVXWISMNTaGRZt24du3btIh6PT3jf3VV1iIig0JgkkUhw7Ngx\nNmzYgJlN+J6qDhFpdAqNPFpbW9m/fz+7d+9W1SEikkWhMQVVHSIiEyk0phGk6mhra1N4iEhDUGgE\nNFXVMTo6qiErEWkICo0iTFV1gIasRKT+KTRKkF115F7XkRmyuv3229m4cSMbN26ko6NDISIidUGh\nUaJM1XH8+PG8Q1ZDQ0McOHCAAwcOsHv3bs17iEhdUGjM0HRDVhmZeY+2tjZuv/12VR8iUpMUGmWS\nGbLasmULGzZsoKmpKe9yo6OjHD16dLz6UICISC3R8zTKqLW1dfzJfqlUij179nD27FmOHz8+6dGx\n8EaAHD16lGQySVtbG83NzSxatIjNmzfrKYEiJUqlUuzYsYOvfe1rrFixgueff54LFy6E3axABgYG\nOH/+PK+//nqg5Zubm8f7+Oqrr/Lggw9W9JEOCo0KyQ2Qrq4uDh06lDc84I0AyVCISLlkdqCPPfbY\nhPcXLVoU6R1qsTvPjGvXrjEwMFChVkXPhQsXOHv27Pjru+++G6BiwaHQqILMvEeQ6iNjqhDJUJiE\nI/M5RnFHm2tgYGDa3zWpP/v27VNo1INih6+y5YZIRr4wyZYbLLk7vP7+fhYuXFh0X6J+lDqVIH3O\n17+BgQF+/OMf11x/pfHke7BcuZi7V2zlYVi1apWfPHmy6J9LJpN8/vOf5+rVqzNuQ2aM8fTp01y6\ndGna5YeGhrhy5QojIyO4O0NDQzNuQ66mpibcneHh4bKvW0TKb8GCBQVPqMlWrjkNMzvl7qumW27a\nSsPMbgLWAIuAa8DPgLPufr7oVkVUMpkcHwcsh9wxxiioRBCJVEPQnWe2uXPncv3117N69eqaq4hn\nMuzc29tLe3t7+RuVpWBomNks4L8D/548p+aa2c+AbwH/092LP7QPyMw+DDwGxIEvu/tflXsb+/bt\nK/cqRSJvyZIlrF69OnBFXAlmRqHRjrlz57J48eKSTwapxg60EU1Vafxn4G7gPLAfuAjMBT4IvBd4\nC9ABbDGzp4F7yl19mFkceAL4feDnwPfN7JC7l/UwftOmTXz7298u5ypFIu/cuXOcO3cu7GZMKbt9\n083f5erv72f58uU1V2kUq9onxBSc0zCz88Cvgfe7+6tZ7/8l8DngN4E7gE8CfwBcAv7Q3f++bI0z\nawX+k7vfMfb6AQB3f7jQz9TinEYlZc+XAMTjcd70pjcxe/ZsAK5evcqcOXOKXm+U+lisIH3O17/s\nI99i6SwmqaRYLMbNN98MEOqcxkLgq9mBkc3dB4F9wD4zez/wFPAtM1vu7r8ousX5vZX0HErGz4H3\nlWndEyQSCW699daGK2cbsYQPq89hnqqbOWNM4VWfRkdH+elPfwqEe53GOWBJkJW4+3Nm1g68QLoK\n6Zhpw4phZgkgAdDS0kJvb29J6xkcHCz5Z2uV+lxdH//4x+nr6+Pw4cNVvQBtZGSE/v5+AJYtW0Z/\nfz/uzksvvVRwTkFq15NPPsmtt95akXVPFRpPAQ+a2fvd/bnpVuTuPzWzrwN/WLbWwT8BN2e9ftvY\ne7nbTgJJSA9PlXIUWc/DU5nhmHxnoJRjeOry5cs1teMpdnjqtddeK3lYKpeO9KUa7rrrropV01PN\nabwJOEl6p90JPOnuI5k5DXefdEtXM9sBbHH33yhL49JncP0Y+ADpsPg+8Mfu3lfoZ0qZ0yj3Kbci\nUjmxWCzQhHgjTYTfcMMN9PT0cOXKlfDmNNz9NTP7APBN0mcw/YWZPUX6aD/fBm8D/gR4sejWFm7D\nsJndCxwmfcrtV6YKjFLplFtpVPF4nJtvvrno+ztV0lQnHBRzplCjzdc98sgj4V6nAeDuF8bOYPoM\ncN/YlwOY2Q9IH/1fJX3h30rSO/bt5Wyguz8NPF3OdebSKbfSqEZGRiJ52u25c+cKXtR34MCBSe/l\nC5paqTRq7R5ygW8jYmZzgY8CHwHagMU5i/wI+K/u/lRZW1ikWjzlthq3ERGRaCvlyneYuL+pxvBU\nyfeeMrPfIH1K7GzggrtfLGlFZVZqaEB1y9lib1hYyHTju9PdWFA3LMwvu3/PPfccL7zwQpVaJzJz\n3d3dRQdH2e49VYi7XwH+X6k/34hmEhT5wqEcZW2jjftCaX2upduhZwwMDHDs2LGaOrNNykO3Rq9x\nQR7ClCs7JGptzLMeZd/WvpZkwu7MmTOTqquoVovFPnxpaGiIixcjMdARGZW8NbpCo4KKCQuFhFRC\nJuzqvaLMVwnWykR4qU8ozCjXnEZQCo0yS6VS9Pb2cvnyZb74xS+O3/Mpn0xQLF26VCEhMgP5KsF6\nD8p8Qj/lVoLLVBU9PT2Mjo4WHEdWUIhILVNozFDQIahYLMadd95JZ2engkJEapZCo0RBwsLMiMVi\nrFu3TmEhInVBoVGCZDLJ1q1bC85XxONxtm/fzvz582lvb1dYiEjdUGgUIVNdHDx4MO+chYagRKTe\nKTQCmqq6UFiISKNQaExjqurCzFi/fr3CQkQahkJjClNVF/F4nJ07d1b0IhoRkaiJhd2AqEomk3R0\ndEwKDDNjw4YNHDt2TIEhIg1HlUaOqYajVF2ISKNTaGTp6enhsccey1tdaO5CREShMe7+++9nx44d\nk96PxWLs2rVL1YWICJrTANLzF11dXZPej8fjCgwRkSwKDdIPLMllZmzfvl2BISKSRaFB/geWuDuP\nPvooyWQyhBaJiESTQgNIJBJ0dnZOen90dJSOjg4Fh4jIGIXGmEceeYTHH3+cDRs2YGbj74+OjrJl\nyxY2btxIKpUKsYUiIuGLbGiY2R+ZWZ+ZjZrZqmpsc9myZezfv5/du3cTi73xT+PuHDhwgDVr1qjq\nEJGGFtnQAH4IfBQ4Wu0NJxIJdu3aNSE4AEZGRjRcJSINLbLXabj7j4AJQ0XVlDlrKvfeU5l5juxl\nREQahRV6lnVUmFkvcJ+7n5ximQSQAGhpaVm5d+/ekrY1ODjIvHnzJrzX19fHU089xYkTJybcVsTM\n2LZtG+vWrStpW1GRr8/1Tn1uDOpzcdauXXvK3aefCnD30L6A75Iehsr9Wp+1TC+wKug6V65c6aU6\ncuRIwe91d3d7LBZzYPzLzHzDhg1+4sSJkrcZtqn6XK/U58agPhcHOOkB9rGhzmm4+wfd/bY8XwfD\nbFc++eY5XBPkItJgojwRHjmaIBeRRhfZ0DCzjWb2c6AV+JaZHQ67TfBGcMTj8Qnv60JAEWkEkQ0N\nd9/v7m9z9znu3uLud4TdpoxEIsGxY8fyXgio4BCRehbZU26jrrW1lf37948/4W90dBR4IzhefPFF\n5s+fT3t7u57BISJ1Q6ExQ5lrNXKDo6urCzMjHo/zxBNP6JoOEakLkR2eqiWFJsjdneHhYQ1ZiUjd\nUGiUSSY4mpqaJl3FrrkOEakXCo0ySiQSPPvsszz00EN0dnZOqDx0t1wRqQea0yiz1tbW8YnvW265\nZcJcR+ZiwJ6eHnbu3Kl5DhGpOao0KkgXA4pIvVFoVNhUFwNquEpEao1CowoKXQyoe1eJSK1RaFRJ\n5mLA3KcCQnq4SlWHiNQChUaVFRquylQdbW1tCg8RiSyFRggKDVdBeq5D4SEiUaXQCEn2cFVu1QEK\nDxGJJoVGyLKrjty5DlB4iEi0KDQiIFN1HD9+XOEhIpGm0IgQhYeIRJ1CI4KKDY/bb7+djo4OBYiI\nVJxCI8KChsfRo0fZvXu3AkREKk6hUQOChAcoQESk8hQaNSRoeIACREQqQ7dGr0GZ8EilUuzZs4ez\nZ89y/Pjx8Vuw58oEyNGjR0kmk7S1tdHc3MyiRYu47bbbaG9vr24HRKRmKTRqWPazO4oNkAwz4+tf\n//p4iGzevHl8nSIiuSIbGmb2BWAdcA14Efiku18Ot1XRVUqAQPqeV9khkluJKEREJFtkQwP4DvCA\nuw+b2SPAA8D9IbepJpQaIDC5ElGIiEi2yIaGu3876+VzwMfCakstyxcgFy5cYGBgYEYhsnTpUlas\nWMHFixdpb29XkIg0iMiGRo5/B3wt7EbUuuwAgXSIPPzww8Tj8aJDJBMkZkYsFmP16tU0NzcDqCIR\nqWPm7uFt3Oy7wKI83/oLdz84tsxfAKuAj3qBxppZAkgAtLS0rNy7d29J7RkcHGTevHkl/Wwt6Ovr\n4/DhwwwMDIy/NzIyMn6X3Zdffpn+/n6uXr3K8PAwr7zyyoy2d/311zNrVvq4ZPbs2dx0003ccMMN\neZdtbm7mjjvuYNmyZTPaZhD1/jnnoz43hpn0ee3atafcfdV0y4UaGtMxs38L3A18wN2vBPmZVatW\n+cmTJ0vaXm9vb82ffpo9BJUtaCURplgsNj5/EkSpFU09fM7FUp8bw0z6bGaBQiOyw1Nm9mGgE7g9\naGA0mtyAqIVgmEru/EkQ2RP1QejaFJGZiWxoAI8Dc4DvjD3d7jl33xJuk8JX7NlQ9a6UoMm+NqWQ\nRYsWsWLFCp5//vlJVZvmbKSRRTY03P3tYbchCko542kquUNA/f39LFy4sOj1TLVTHRgY4Pz587z+\n+usAXLt2bcI8Sthyr00pRbEVzkzVQlAVGhoNS39/P8uXLy/4e1provI7EOk5jVLUy5xGKpWiq6uL\nQ4cOFR0SheYG8v3SVavP+XYotT6cVm3Fzvlky3dwMFXwF0ufZfUsWLCApqam8dfNzc2sWLGC06dP\nc+XKFR588EESiUTR6635OY1GVMrQU+6OJCpHI7lyT/fNUJgEV8pQnNSfixcvTnh94cIFzp49O/76\n7rvvBigpOIJQaERAsVVFJiiWLl0ayYAoRjFhAtUJlAULFtDS0jJ+9Hbp0qXx7w0NDU36oxWJmn37\n9ik06lHQsMiuJqJaSZRboTCBylcnFy9e5OLFixOO3grJHSqYyty5c1m8eHHRw0u1VHmZGS0tLRVZ\nd9TmxqJs06ZNFVu3QiMkyWSSrVu3MjIyUnCZWCzGnXfeSWdnZ92HRDGKrU6gcjveYquO8+fPFz0v\n0dzcTFtb24STC4p19epV5syZM2m9+aqpmfjlL39Z8xPOlTCT+aiM3BNMMso1pxGUQqPKMtXFwYMH\nyXcSQj0NPVXbVNUJvBEqZ86cmTQpXK2j+SjNS+SOhTe6cuzY86nm6EA1TmxRaFTRVNWFqorKy4RK\noT+s7EplulOKa2W4qB4F3bkXc8ptowz7loNCowqmqi7MjPXr1yssImC6SiVbsdckNFLQVOqIHYrb\nuUfpFPp6otCosKmqi3g8zs6dOys6/iiVUUzAZIR58Vulr9PIXqeO2OubQqOCkskkHR0dk44uVV00\nplKCplx01C3lotAos1QqRW9vL5cvX+bRRx+dFBiqLkSklik0yiiZTHLvvfcyPDysuQsRqUsKjTIp\nNBQF6YnBXbt2qboQkZoXC7sB9WCquYtZs2YpMESkbqjSmKF8gRGLxbjvvvuYP38+7e3tGo4Skbqh\n0JiBQoGhykJE6pWGp0qkwBCRRqTQKIECQ0QalUKjSAoMEWlkCo0iKDBEpNEpNAJSYIiIKDQCUWCI\niKQpNKahwBAReUNkQ8PM/ouZ/cDM/sHM/s7MFle7DQoMEZGJIhsawBfc/d3u/h7gAPCX1dy4AkNE\nZLLIhoa7v5z18jrgYrW2nUqluOeeexQYIiI5LPcW3lFiZg8Bm4HXgPe5+6UCyyWABEBLS8vKvXv3\nlrS9wcFB5s2bx44dO+jp6cleP9u2bWPdunUlrTfKMn1uJOpzY1Cfi7N27dpT7r5q2gXdPbQv4LvA\nD/N8rc9Z7gHgr4Osc+XKlV6qI0eOeHd3t8fjcQcc8Hg87t3d3SWvM+qOHDkSdhOqTn1uDOpzcYCT\nHmAfG+oNC939gwEX/SrwTCXbAtDT08OXvvSl8WEpM+Ouu+7SkJSIyJjIzmmY2TuyXq4HXqjk9u6/\n/3527NgxYR5j1qxZbN68uZKbFRGpKVG+Nfpfmdk7gRHgH4GOSm0omUzS1dU14b1YLMbjjz+uZ2GI\niGSJbGi4+6ZqbWvfvn0TXpuZzpQSEckjssNT1bRp08R8+sxnPqPAEBHJI7KVRjVlAuLJJ5/UxLeI\nyBQUGmMSiQS33nor7e3tYTdFRCSyNDwlIiKBKTRERCQwhYaIiASm0BARkcAUGiIiEphCQ0REAov0\nrdFLYWb9wE9L/PEbgV+VsTm1QH1uDOpzY5hJn/+Zuy+cbqG6C42ZMLOTHuR+8nVEfW4M6nNjqEaf\nNTwlIiKBKTRERCQwhcZEybAbEAL1uTGoz42h4n3WnIaIiASmSkNERAJryNAwsw+b2f81s5+Y2Wfz\nfN/M7L+Nff8HZvbbYbSznAL0+V1mljKzq2Z2XxhtLLcAff43Y5/vGTM7YWbvCaOd5RKgv+vH+vuC\nmZ02sw+E0c5ymq7PWcv9jpkNm9nHqtm+SgjwObeb2a/HPucXzOxzZW2AuzfUFxAHXgT+OTAb+Adg\nac4yHwGeAQx4P/D3Ybe7Cn2+Cfgd4CHgvrDbXKU+/y7w5rH//4Na/pwD9ncebwxJvxt4Mex2V7rP\nWcv9HfA08LGw212Fz7kd+Gal2tCIlcZ7gZ+4+z+6+zVgL7A+Z5n1wB5Pew6Yb2ZvqXZDy2jaPrv7\nS+7+fWAojAZWQJA+n3D3S2MvnwPeVuU2llOQ/g762F4FuA64WOU2lluQv2WAPwf2AS9Vs3EVErTP\nFdOIofFW4GdZr38+9l6xy9SSeutPEMX2+c9IV5e1KlB/zWyjmf0f4H8D/6FKbauUaftsZm8FNgK7\nqtiuSgr6e/27Y0ORz5jZsnI2QE/uk4ZnZmtJh0Zb2G2pNHffD+w3s38J7DGzd7n7aNjtqqAvAfe7\n+6iZhd2WajkNLHb3QTP7CHAAeEe5Vt6IlcY/ATdnvX7b2HvFLlNL6q0/QQTqs5m9G/gysN7da3m4\npqjP2N2Pkj5oXFDhdlVSkD6vAvaa2TngY8BOM9tQneZVxLR9dveX3X1w7P+fBprM7MZyNaARQ+P7\nwDvM7LfMbDbwCeBQzjKHgM1jZ1G9H/i1u/+i2g0toyB9rjfT9tnMFgN/C/ypu/84hDaWU5D+vt3G\nDrfHzgg0d++vflPLZto+u/tvufsSd18CfAPY6u4Hqt/UsgnyOS/K+pzfS3o/X7YDooYbnnL3YTO7\nFzhM+kwkLFhMAAACKklEQVSEr7h7n5ltGfv+btJnWXwE+AlwBfhkWO0thyB9NrNFwEngBmDUzP4j\n6bMyXg6t4TMQ8HP+HOkj7Z1jf2PDXqM3uAvY302kD4aGgFdJ73BqVsA+15WAff4Y0GFmw8BrwCey\nToCYMV0RLiIigTXi8JSIiJRIoSEiIoEpNEREJDCFhoiIBKbQEBGRwBQaIhVkZtvMzM3sj8Nui0g5\nKDREKmvl2H9PhtoKkTLRdRoiFWRmZ0nfUG5+OS+wEgmLKg2RCjGz64B3As8rMKReKDREKudfkP4b\nO5X9ppm92cwOjs11fMnMmsJpnkjxGu7eUyJVlHlM8HhomNn7gK8Bbwb+yN2/EUbDREqlSkOkcjKT\n4KcAzOzTwDHgZWCVAkNqkSoNkcr5beAV4FdmdhC4E/gboMPdXwu1ZSIl0tlTIhVgZnNJB8ZLwDBw\nE/Dn7v7lUBsmMkManhKpjPeQruTnAIuBbygwpB4oNEQqIzMJvh34JvAnYw+2EqlpCg2Rysi+EvwT\npCfDv2hmHw2vSSIzpzkNkQows9PAu4Dr3X1k7HG6z5Ge2/g9d38u1AaKlEiVhkiZmdls4DbgjLuP\nALj7BdLPnb8K9JjZ20NsokjJFBoi5bccaAKez37T3c8CHwV+E3jGzG4MoW0iM6LhKRERCUyVhoiI\nBKbQEBGRwBQaIiISmEJDREQCU2iIiEhgCg0REQlMoSEiIoEpNEREJDCFhoiIBKbQEBGRwP4/f73k\nhpQAODwAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb88a804c88>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot(kran,theta,'k.');\n",
"grid()\n",
"ylabel(r'$\\theta$',fontsize=20)\n",
"xlabel(r'$k$',fontsize=20);"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4.52 ms ± 9.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
]
}
],
"source": [
"%%timeit\n",
"\n",
"Hkd(dh)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/usr/local/lib/python3.4/dist-packages/scipy/sparse/compressed.py:774: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.\n",
" SparseEfficiencyWarning)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"5.61 ms ± 41.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
]
}
],
"source": [
"%%timeit\n",
"\n",
"dh.Hk()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"gnubands=loadtxt('bi_mono_forwilson_SOC.gnubands')"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/usr/local/lib/python3.4/dist-packages/scipy/sparse/compressed.py:774: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.\n",
" SparseEfficiencyWarning)\n"
]
},
{
"data": {
"text/plain": [
"(-5, 5)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGAhJREFUeJzt3WuMZGldx/Hfv6vAyy4E4gyXsMw0Go0iyG0kCkRBIrvM\n4K4mvMBLjMSkdruDwUSzO0Ci+0acXQ0XY3pg3aASjKuJRNadkZFLuBhcYUYWlkskC07DQgizoG5c\nXpDu+ftiumZrqutUnVN1znM7308yyXT1qaqnz3meXz31nOc8x9xdAIByrMUuAACgXQQ7ABSGYAeA\nwhDsAFAYgh0ACkOwA0BhCHYAKAzBDgCFIdgBoDDDGG964MABX19fj/HWAJCtc+fOPeTuBxdtFyXY\n19fXdfbs2RhvDQDZMrPtOtsxFAMAhSHYAaAwBDsAFIZgB4DCEOwAUBiCHQAKQ7ADQGEIdgAoDMEO\nAIUh2AGgMAQ7ABSGYAeAwhDsAFAYgh0ACkOwA0BhCHYAKExrwW5mAzP7tJnd09ZrAgCaa7PH/npJ\nX2zx9QAAS2gl2M3sGknHJN3ZxusBAJbXVo/9bZJulnSxpdcDACxp5WA3s1dJ+pa7n1uw3cjMzprZ\n2QsXLqz6tgCACm302F8s6XozOy/pLkm/YGbvmd7I3e9w9yPufuTgwYMtvC0AYJaVg93d3+Du17j7\nuqTXSPqwu//GyiUDACyFeewAUJhhmy/m7h+R9JE2XxMA0Aw9dgAoDMEOAIUh2AGgMAQ7ABSGYAeA\nwhDsAFAYgh0ACkOwA0BhCHYAKAzBjs5tbm5qOBzKzDQcDrW5uRm7SEDRWl1SAJi2ubmpkydPXv55\nd3f38s9bW1uxigUULZtg39zc1B133KHd3V0NBgONRqN9wbB+/NRSr33+xLE2iogp68dPafsd75z5\nu5PveKdOP37/fudYtKdJe0hlv9OG22HuHvxNjxw54mfPnq29/XSvb+yq5x7VgWuv/Frf9ADPqkhU\nkuVM78vzJ47JzCq3n657VY2a47HYKvtuUZh2uf/baH99asNmds7djyzcLodgHw6H2t3d3ff4YDDQ\nzs5Om0W7opKUWjnaNCvMJ6167PrUaJsItV+W7UE30Xa5S27DRQV7k15fW6YrR52hoD6p23iqvm1t\nbGwstf/6HPSLPkRxpRL3V1HBHrLHPm39+Ck9dGZLj9x3et/vlg2nnC3TG+ryQ7Hk4Zs+f4i1rZRe\nfFHB3navr6mYHyypyKlhtBn2ob6plfwBlZKc6vEsRQW7FK6BzRJjKGharL8/94YwVmesePrva6ND\nkePMlD7ItV4XF+wxVfXYZWs6fPPdl3/sqoKE/sZS4tjkIrMCePv26yW/uH/jqeO+SB/2X65yq+sE\ne4vqBGuXswfaCpgmUq/gIaTwTQ3hhJgBNLZs+yLYW9b3oaA+4twKUlM32FkrpqatrS3t7OzI3bWz\nsxN0NsxgMGj0ONoxGo0aPQ6kgmDPAAETx9bWljY2Ni5/gA4Gg15OcUV+GIrJBBdIAWCMHQAKwxg7\nAPQUwV4TN4tYHvsOCCub9dhj4mYRy2PfAeH1Yox91QsPqi4Q6tt85mX247yLq/zijKt5gRnavngo\n1wvwijt5umhWSJc3C5h3gdDhW+7JtpIs0sbqgov23Syl7k/MF+qGHzmvmllUsNe9g1JXB2feFYjX\n/P77On//kNpeO6Pp1ZusctgfqQRsTguCFRXssS/tbrpWTOqVY5auyt/VCok57uNSrTJMktJxzKEN\nFxXsKayVUvcCoRwqx7RxmbtcnbLNi6sI+rhK3/8pt+Gigj12j30ZKVeOsRzKWEduS6/mpvQgr5Ji\n+ygq2GPfQWkVXfeGl5FihW1LX0OoCyXXkyZSasPBgt3Mni7p3ZKeLMkl3eHub5/3nNzuoLSqlBpI\nSpU0hFWDPnS9i33ymG8/s6XShkMG+1MlPdXd/8PMHifpnKRfdvcvVD2nr2vFxAzVVCpmbE2Cvq1v\nik1PLk6Xp+7zlzmufMNpJnY7ijYUY2bvk/Tn7v6Bqm36GuxS+IoRuyKmbl5otnnnqq73fSkzU3IR\nq5MWJdjNbF3SxyQ9y90frtquz8E+FiJw+zbs0rYUZmMhXTE6TcGD3cyulvRRSX/k7u+d8fuRpJEk\nHTp06AXb29utvG/u2q4cjJG2J8fZWAgvZMAHDXYze4ykeySdcfe3LNqeHvt+iwK5yZIKhHk7cp6N\nhfBCdKpCnjw1SX8t6Tvu/rt1nkOwzzddQR46s6VH7ju9b7tQSyr0Wc6zsRDPvHMeK65bFSzYXyLp\n45LulzQ+0/RGd9+fRHsI9mYYEgDKsOp5r7rBvvJ67O7+r5KqzzJhZbNCfd7jANIU6ls1d1DKwGAw\naPQ4gH4j2DMwGo0aPQ6g37g1XgbGJ+s4iQegjiwWAQMA1D95ylAMABSGYAeAwhDsAFAYgh0ACkOw\nA0BhCHYAKAzz2IFC1bn5BgvHlYlgBwqwzL1S14+fYrnnQnGBEpCpNtf/JuDzEO2ep3UQ7MByur6Z\nAwGfNoIdKESM2x1yv9w0BVuPHUC7Zo2Xhw7Y8fsR8Hmixw4kIOWbkBPu6aDHDiQu5TCfRO89P/TY\ngYByCfMqhHtcnDwFEpB7kM/CzJl4CHagBZubm43uXJXCic9Q6L2HR7ADK9rc3NTJkyf3PX7Vc4/q\nwLWblc/rU9DRew+LYG9Z054b0lZnHZXt26+X/OK+xweDgXZ2drooVrbovYdBsLeoque2sbFBuCdu\nXoAvCiEzq/xdjHaTOnrv3Ssu2GP2mG1tMLPnJlvT4Zvv3vdwjpW6Tg92UpO/MdSxW2YhrHmGw6F2\nd3f3PU6Pfb5UA77LOh5KUcEeqsdcdeC3b3tV5XOm918uJ89WKeeiBjL5Ol0cu1V64U3wTW01MZdC\nqLJKHU+hHRcV7FU9pzo95jY+pVfpuaU23a3r3tT031s1Ti1bk198dJ+m2pvi3Eo7mh5faf+FUU2f\n16ZUvoUUFexNxjq7+KRtq+cWs3LEeu95x+7wLfdc8XPsDz2kI7UO0VjsgC8q2FMY62yz5xaycsSu\niCkcO6BtETtK5QR7qWOdXU4Rix3oY6UeO0AK386KCnap3LHOLipGanOKSz12wFioNldcsJeujYBP\npZcO9FGI9lc32Nc6eXc0dv7EsbmzADY3NzUcDmVmGg6H2tx89JL2yZsST74OgHAWteGQ6LEnaPqT\nv86aJYQ5kJYuhmeCDsWY2XWS3i5pIOlOdz8xb3uCvZ5xxag7FxxAetqcgh0s2M1sIOlLkn5R0oOS\nPiXpV939C1XPIdibYc0SAFLYMfYXSnrA3b/i7t+TdJekG1p4XewZDAaNHgfQb20E+9MkfW3i5wf3\nHkNLRqNRo8cB9FuwWTFmNjKzs2Z29sKFC6HetghbW1va2Ni43EMfDAZc4AOgUhtj7D8r6VZ3v3bv\n5zdIkrv/cdVzGGMHgOZCjrF/StKPmtkzzOyxkl4jaf+SiwCAIFYOdnffkfQ6SWckfVHS37v751d9\nXVxp3gVKADCJC5QiaXJl2kNntvTIfaf3Pb7opspjXLwEtC/GOvGsFZOgZdeYrro13zI3+mjyvgCu\ntGwbbmt9eYI9EW0c0DYvUGKhMKCZlG76QbBH1maArq2tzQxwM9PFizOWGqghpcoKpCjFThDBHkFX\nYdn1kgIpVmAghi47PG20s7rBPlzq1XGFLoMxxPKfVTf/JuTRFyHacMj2RLDXNH0XoO9/9rVXzEjp\n8vZ2g8Fg5n1DZWtaP36q1feeXE86tTsxAW3ruiMTqw0xFFNDjPt2TlaIee9/+vHHLm/XZTm6fA8g\npFDnl7Jfj72p3IJ9OBzO7DHXmW64jFkVYt59Q0P1Cgh55Cxk/e2qTRLsLVk/fkrbt72q8vdt779l\nK0TIr3wEPHISur522RYJ9hVMV4RQPfZVK0To8TwCHimLUT+7boME+xKqKkKIMfa2KkSMkzUEPFIR\n8/qMEG2PYG+gTjDNG+Nu6/3bqhAxZ7MQ8oghdr0L1eYI9hpiV4bJMnQxdzZmsE7u26MPn+rsQxH9\nVnIbnoVgnyOFyjBZji6nW8XuNV/9vGMzV6bkDlBYRV/a8DSCfUpqa6OEnKIY828NPVUUZUsl0KU4\nQ54sKbAnpYowFrJCnD9xLGq4z7xids7jwCyptePUr8ouNthTqwhjMSpEzHBftByClNbxQVpSbMc5\n1Nvigj3FijAWs0LECvfRaDR7quhNN2prr0w5NBSEQxteXTFj7ClXBimdChGjHHWmiqZ2DgTh0YYX\nK+7kaVU4pF4ZpDQqxKTUyjMth2OK9uRwvFNpM0UFe9WVn+ObOcfe2fOkUiGmpVquSTk0eCwvl+Ob\nUlspKthznTKXUoWYJfXyjeUSAHhU3RvEpH48U2sjRQV717eGa1tOQZRaxZ2Hcfg0VYV47scnxbZR\nVLDn1GNPsTIsknOZpbzKnbtSQ3xSynWrqAuUqqbMjUajCKWplmNASo9Ohcyp/NyntXt9CPFpObWB\nebLosUvdrq64qpKCJeeKzVDN8voY4tNyqPtFDcWkLIfK0FQpf1MbQZ9yh2IZ805q5n68l5VTx4xg\nb9msBt71jaRjKyXgpeV6pDFuYt4WArye3Oo4wd6iefPo/+/T9aZ15Sq3it/Eoil527dfL/nFfY+n\ndNKeAF9OTr30SQR7i3KaldOVkgO+yrxptodvuaf267R1u8OuXr9vcq7LBHuLcptH35VceznLauMD\nve6FOov0YX93rYT6W9R0x5jWj5+SbK3yK3mfjBvDeGpkro2jrjam2Za+j3JQQqA3tRa7AKmanNe9\ncdONM7dJbR59KNMBX6qtrS1tbGxc/gAfDAZZnDjFoyaHXfoS6hJDMTPNGoMrbdpbW/rYG0L6Sq2X\njLEvIeXKkPoHS8r7Dv1Rej0MEuxm9ieSfknS9yR9WdJr3f1/Fj0vtWBPvTLkNJ869X2JMvWl3oUK\n9ldI+rC775jZbZLk7rcsel4qwZ56ZRiXr2o+tWxNh2++O+mynz9xLPlvG8hXLm142rJlDTIrxt3/\nZeLHeyW9epXXCyXFyjCvAthtM0Jduhz2s54b++8av//VzzumR+47ffnx3d3dy98+CHcsK7c2HFpr\nY+xm9k+S/s7d37No21g99tQqQ921TJrOp05pMSwu7kLbUrrAKHRba63HbmYflPSUGb96k7u/b2+b\nN0nakfQ3c15nJGkkSYcOHVr0tq1KKdCXqQhN51NPL2kbsyHMCvV5jwNVUm3Hscsyy8o9djP7LUk3\nSnq5u3+3znNC9dhT2vmrlmXVcepY+6Kqx57y+QGkJZV2nEI5Qp08vU7SWyT9vLtfqPu8EMGeyte1\nFCrDpNDlmTejZ7w6ZqiyIC8ptZ1U8iRUsD8g6fskfXvvoXvd/aZFz+sy2FM5AFJaZZkWstEs+raR\nUgNGfCnVh5TKIvXwAqWUQjS1yjAP+w2pSOn4p1SWScUFe1WvL7UDkFJQ1pXqPpTSKA+6ldrxTrkN\nFxXs8250ceDazSQOQGqVcxmpVegS9imqpXh8U2sD04oK9tTnQqdeGZpI8W9JMQAwX91VP1M5nrnU\nsaKCPeUbXaQYhG1I8e/KpfH1UYpXP9eVYl2vUlSwp9pjz6lCLCPVvy+lK2v7LvcP21TreJWigj2F\n1Q0nT97K1nTVc64r/kbWUvoVP/dgyVUJ+z31uj1LUcEuxV2PPIUPlphyaAAlBE3qStnHOf8dxQV7\nTKkOBYWUS2NgmKYbOXy415H730Gwtyjlk7eh5dQwCPnV5XS8FynhbyHYW2Rrg5k3uuhTj31Sjg1k\n2ZDv601CcvmGVleOdXYWgr0F48rw0JmtK24WMdaXMfZZcm4odW+I0MdzK6UFupR3XZ1GsK9oujL0\ntec2T0kNZlbYL7olYQih9m2JgS6VVUclgn0lpVWGLpW8r2KfW5l39WYb+7vUMJfKrZcE+5JKrRBd\nK3G/pTobatl7a+Z8dWgTJdbFMYJ9CSVXiBBK2385jbHntjZLV0qrg9MI9gZK/koaWmkNi3Mr+Sit\n7s1CsNfUh8oQGvsUIfWpY0aw10AAdadPjQ3x9K0N1w32tRCFSVHfKkRo508cu7xvr37eMQ2HQ5mZ\nhsOhNjc3I5cOJaANVxvGLkBo9CTDOvrwKZ2cuLhrd3f38glJxqqxLEJ9vl4NxVAZwkt1yiDy1Pc2\nXHcophc9dnrp8cwK9XmPA1X6HupNFB/sVIa4BoPB7BC3Na0fP8VxwUJ0zJorNtipDGkYjUazL/K5\n6Uadlgh3VKINL6+4YKcypGV8grTqIp/146f4VhXQsssRhEadWE1RJ0+pDPniA7ldTRYQW7QcQcjj\nQT2Yr7gLlOZd2k1lKAfHcnld3DGq6xUmZ70Px71aUcFetRjTVc89qgPXXrrYhcpQFhp6PbFu/9dW\nL5/j3ExRwc5c6P7ivqX7pb5P6q40OZZa+VNWVLDHvuEB0tDX3l1f1lHHYkVdoFQ1F3owGEQoDWKZ\nDLPUe62rIMixqiyCvWou9Gg0ilAapGA65LsKwy7XYw91YhL9k8VQjMQND9BMG+O8y9xBifFldKmo\nMXagS1VhvH379ZJf3P8LW9Phm++ufD3CGl0JGuxm9nuS/lTSQXd/aNH2BDtywEl7pCbYjTbM7OmS\nXiHpq6u+FpCSqpPznLRH6tq4g9JbJd0siS4MilJ1cp6T9kjdSrNizOwGSV9398/M+9oK5GjRAmZA\nqhaOsZvZByU9Zcav3iTpjZJe4e7/a2bnJR2pGmM3s5GkkSQdOnToBdvb26uUGwB6p/OTp2b2bEkf\nkvTdvYeukfQNSS9092/Oey4nTwGguc6vPHX3+yU9aeINz2tOjx0AEEYbJ08BAAlpbUkBd19v67UA\nAMujxw4AhSHYAaAwBDsAFIZgB4DCEOwAUBiCHQAKQ7ADQGEIdgAoDMEOAIUh2AGgMAQ7ABSGYAeA\nwhDsAFAYgh0ACkOwA0BhCHYAKAzBDgCFWfpm1iu9qdkFSdtLPv2ApFzuq0pZ25dLOSXK2pVcytpF\nOQ+7+8FFG0UJ9lWY2dk6d+lOAWVtXy7llChrV3Ipa8xyMhQDAIUh2AGgMDkG+x2xC9AAZW1fLuWU\nKGtXcilrtHJmN8YOAJgvxx47AGCOZIPdzK4zs/80swfM7PiM35uZ/dne7z9rZs+PUc69siwq66/v\nlfF+M/uEmT0nxXJObPfTZrZjZq8OWb6pMiwsq5m91MzuM7PPm9lHQ5dxohyLjv8BM3u/mX1mr6yv\njVTOd5nZt8zscxW/T6lNLSprEm1qryxzyzqxXbh25e7J/ZM0kPRlST8s6bGSPiPpmVPbHJX0z5JM\n0s9I+veEy/oiSU/c+/8rY5S1TjkntvuwpNOSXp3wPn2CpC9IOrT385MSLuutkm7b+/9BSd+R9NgI\nZf05Sc+X9LmK3yfRpmqWNXqbqlvWiXoSrF2l2mN/oaQH3P0r7v49SXdJumFqmxskvdsvuVfSE8zs\nqaELqhpldfdPuPt/7/14r6RrApdRqrdPJel3JP2DpG+FLNyUOmX9NUnvdfevSpK7xypvnbJ+U9Lj\nzMwkXa1Lwb4TtpiSu39s772rpNKmFpY1kTY1Lsui/SoFblepBvvTJH1t4ucH9x5ruk0ITcvx27rU\nKwptYTnN7GmSfkXSyYDlmqXOPv0xSU80s4+Y2Tkz+81gpbtSnbL+haRnSvqGpPslvd7dL4YpXiOp\ntKmmYrWpWmK0q2GoN4JkZi/TpUr4kthlqfA2Sbe4+8VLncukDSW9QNLLJf2ApH8zs3vd/UtxizXT\nGyR9VtLLJP2IpA+Y2cfd/eG4xcpfBm1KitCuUg32r0t6+sTP1+w91nSbEGqVw8x+StKdkl7p7t8O\nVLZJdcp5RNJde5XvgKSjZrbj7v8YpoiX1Snrg5K+7e6PSHrEzD4m6TmSQgd7nbK+WNKb/dJg6wNm\n9l+SflzSJ8MUsbZU2lQtCbSpusK3q1gnHBacjBhK+oqkZ+jRE1I/ObXNMV15oueTCZf1kKQHJL0o\n5X06tf1fKd7J0zr79CckfWhv2x+U9DlJz0q0rG+VdOve/5+sS2F5INK+XVf1Cckk2lTNskZvU3XL\nOrVdkHaVZI/d3XfM7HWSzujS2eR3ufvnzeymvd+/Q5fOLh/VpYP7XUlRppDVLOsfSPohSVt7n9o7\nHnhxoJrlTEKdsrr7F83s/bo0xHFR0p3uPne6WayySnqzpL80s8/q0nmtW9w9+OqEZva3kl4q6YCZ\nPSjpDyU9ZqKcSbQpqVZZo7epBmUNX6a9TxEAQCFSnRUDAFgSwQ4AhSHYAaAwBDsAFIZgB4DCEOwA\nUBiCHQAKQ7ADQGH+H1mBvbHglYncAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb88a732d68>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot(gnubands[:,0],gnubands[:,1]-dat.variables['Ef'][0]*sisl.unit_convert('Ry','eV'),',')\n",
"plot([gnubands[0,0]],[dh.eigh(k=(1/3,-1/3,0))],'ko');\n",
"plot([gnubands[500,0]],[dh.eigh(k=(1/6,-1/6,0))],'ko');\n",
"plot([gnubands[1000,0]],[dh.eigh(k=(0,0,0))],'ko');\n",
"plot([gnubands[1500,0]],[dh.eigh(k=(1/4,0,0))],'ko');\n",
"plot([gnubands[2000,0]],[dh.eigh(k=(0.5,0,0))],'ko');\n",
"ylim(-5,5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```latex\n",
"BandLinesScale ReciprocalLatticeVectors\n",
"%block BandLines\n",
"1\t0.333333333 -0.3333333 0.0000000 #K\n",
"1000\t0.000000000 0.00000000 0.0000000 #G\n",
"1000 \t0.500000000 0.00000000 0.0000000 #M\n",
"1000\t0.333333333 -0.3333333 0.0000000 #K\n",
"1000\t0.000000000 0.00000000 0.0000000 #G\n",
"%endblock BandLines\n",
"\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from scipy.sparse import diags"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def Skd_MK2(dh, k=(0, 0, 0)):\n",
" dtype = np.complex128\n",
" k = np.asarray(k, np.float64)\n",
" k.shape = (-1,)\n",
" # sparse matrix dimension (dh.no)\n",
" V = zeros((len(dh), len(dh)), dtype=dtype)\n",
" # Get the reciprocal lattice vectors dotted with k\n",
" kr = dot(dh.rcell, k)\n",
" # Calculate all phases\n",
" phases = exp(-1j * dot(kr, dot(dh.cell, dh.sc.sc_off.T)))\n",
" # Now create offsets\n",
" offsets = - np.arange(len(phases)) * dh.no\n",
" dia = diags(phases, offsets, shape=(dh.shape[1], dh.shape[0]), dtype=dtype).todense()\n",
" #dia = diags(phases, offsets, shape=(dh.shape[1], dh.shape[0]), dtype=dtype).todense()\n",
" V = (dh.tocsr(dh.S_idx)).todense().dot(dia)\n",
" return kron(V,eye(2))\n",
" \n",
"def Hkd_MK2(dh, k=(0, 0, 0)):\n",
" dtype = np.complex128\n",
" k = np.asarray(k, np.float64)\n",
" k.shape = (-1,)\n",
" # sparse matrix dimension (2 * dh.no)\n",
" V = zeros((len(dh), len(dh)), dtype=dtype)\n",
" # Get the reciprocal lattice vectors dotted with k\n",
" kr = np.dot(dh.rcell, k)\n",
" # Calculate all phases\n",
" phases = exp(-1j * dot(kr, dot(dh.cell, dh.sc.sc_off.T)))\n",
" # Now create offsets\n",
" offsets = - arange(len(phases)) * dh.no\n",
" dia = diags(phases, offsets, shape=(dh.shape[1], dh.shape[0]), dtype=dtype).todense()\n",
" V[::2, ::2] = (dh.tocsr(0) + 1j * dh.tocsr(4)).todense().dot(dia)\n",
" V[1::2, 1::2] = (dh.tocsr(1) + 1j * dh.tocsr(5)).todense().dot(dia)\n",
" V[1::2, ::2] = (dh.tocsr(2) - 1j * dh.tocsr(3)).todense().dot(dia)\n",
" V[::2, 1::2] = (dh.tocsr(6) + 1j * dh.tocsr(7)).todense().dot(dia)\n",
" return V"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/usr/local/lib/python3.4/dist-packages/scipy/sparse/compressed.py:774: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.\n",
" SparseEfficiencyWarning)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"2 ms ± 18.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
]
}
],
"source": [
"%%timeit\n",
"dh.Sk()"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"958 µs ± 24.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
]
}
],
"source": [
"%%timeit\n",
"Skd_MK2(dh)"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"797 µs ± 2.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
]
}
],
"source": [
"%%timeit\n",
"kron(Skd(dh),eye(2))"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 300/300 [00:47<00:00, 6.32it/s]\n"
]
}
],
"source": [
"kran=linspace(0,0.5,300)\n",
"qran=linspace(0,1,7,endpoint=False)\n",
"theta=[]\n",
"for k in tqdm.tqdm(kran):\n",
" W=matrix(eye(2*dh.no))\n",
" for q in qran:\n",
" SK=matrix(sl.fractional_matrix_power(Skd_MK2(dh,k=(k,q,0)),-0.5))\n",
" #SK=kron(SK,eye(2))\n",
" HK=SK*matrix(Hkd_MK2(dh,k=(k,q,0)))*SK\n",
" w,v=eigh(HK)\n",
" P=zeros_like(HK);\n",
" for i in arange(len(w)):\n",
" if w[i]<0.001:\n",
" P=P+(v[:,i]*v[:,i].H)\n",
" W=W*P\n",
" ev=eigvals(W)\n",
" theta.append((angle(ev[abs(ev)>1e-10])))\n",
"theta=vstack(theta)"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAETCAYAAADKy1riAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+Q1Hed5/Hnu5sB3JAsDiGDpeHYi0YLgh4Lq84yXIbV\nNa57BBC39PZ2qfP20mFI9jwkTky2yr2ru5jNGClzlwDT8bxd6qygJcuP0eRQa5kA0tkSSFYcvPPM\nFqJbYkYGTCYkMD/e90dPT3p6ume+3dPd3293vx5VU0n3fOf7/Xzome/r+/58vj/M3REREQkiFnYD\nRESkdig0REQkMIWGiIgEptAQEZHAFBoiIhKYQkNERAJTaIiISGAKDRERCUyhISIigc0KuwHlduON\nN/qSJUtK+tlXX32V6667rrwNijj1uTGoz41hJn0+derUr9x94XTL1V1oLFmyhJMnT5b0s729vbS3\nt5e3QRGnPjcG9bkxzKTPZvbTIMtpeEpERAJTaIiISGAKDRERCUyhISIigSk0REQkMIWGiIgEptDI\n0tfXx8MPP0wqlQq7KSIikVR312mUKpVKsX37doaGhojFYjzxxBMkEomwmyUiEimRDg0zmwscBeYA\ns4GD7v7ZSmyrt7eXa9eu4e6Mjo7S0dEBoOAQEckS9eGpq8Dvuft7gHcDa81sTSU21N7eTiz2xj/H\n6Ogo9957r4aqRESyRDo0PG1w7GUTEAcuVWJbra2tfOpTn5oQHMPDw+zZs6cSmxMRqUmRDg0AM4ub\n2QvAS0Cvu/+wUttat24du3btIh6PA+DuPPnkkySTyUptUkSkppi7h92GQMxsPnAY+Ky7H8n5XgJI\nALS0tKzcu3dvSdsYHBxk3rx57Nixg56enuz1s23bNtatW1dy+6Mq0+dGoj43BvW5OGvXrj3l7qum\nXdDda+YL+BzwmamWWblypZfqyJEj7u5+4sQJnzVrlgPjX7FYzLu7u0ted1Rl+txI1OfGoD4XBzjp\nAfbDkR6eMrOFYxUGZvYm4PeBFyq93dbWVp544olJE+MdHR0aqhKRhhbpU26BtwB/Y2Yx0vMv/8vd\nv1ONDWdOte3o6GB0dBRAp+KKSMOLdGi4+w+AFWFtX8EhIjJRpIenoiCRSLBr1y4NVYmIoNAIRMEh\nIpKm0AhIwSEiotAoioJDRBqdQqNICg4RaWQKjRIoOESkUSk0SqTgEJFGFOnrNKJuqus4XnzxRebP\nn097ezutra1hNlNEpGwUGjNUKDi6urowM+LxuJ4CKCJ1Q8NTZZBvqArSN4McHh7WkJWI1A2FRplk\ngqOpqQkzm/C90dFRtmzZwsaNG/UkQBGpaQqNMkokEjz77LM89NBDdHZ2Tqg83J0DBw6wZs0aVR0i\nUrM0p1Fmra2t4xPft9xyy4S5DoCRkRG2bNnCM888Q2dnpybJRaSmKDQqKDP5vXXrVkZGRsbfz1Qd\nPT097Ny5U5PkDSSVSrFnzx4uXLhQ1e329/ezcOHCCe8tWrSIFStW8Pzzz5elPYsWLWLz5s06EKpz\nCo0KSyQSLF++nK6uLg4ePJh5AiGgqqOWlbLzHxgY4Pjx4xMqz3qTTCZpa2ujubm57OtWKEWDQqMK\nWltb2b9/P8lksmDVcejQIe68806FR8gyYQAUPApvhJ1/qUZHRzl69GjF1l9MKPX397N8+fJpqymF\nUXEUGlU0VdUxOjo6Hh5tbW0sXbpUv8hllgmEM2fOTBqqAYVBLSg2lL73ve8FWq5SFVI9BpJl77jq\nwapVq/zkyZMl/Wxvby/t7e3lbVAB+aqOXLFYrOLVRzX7XA2Fho2iFAixWGxGO6iBgQHOnz/P66+/\nHvhnrl69ypw5cya819zczIoVKzh9+jSXLl0quh1z585l8eLFNDc3R+rfN2pm+nnnyv38sz/HK1eu\n8OCDD5Y0T2pmp9x91bTLKTTeUO0daCqVoquri0OHDk35x5b5patE9VGLoRGlYFiwYAFNTU2Bl8/e\n0ZYiajvn7B1iKWEW1LVr1xgYGCj7eutVd3d30cGh0ChBWDvQoOEBE/9Iy1H6RjU0wg6GJUuWsHr1\n6rxH4UNDQ1y8eLGi2xeZiQ996EMcPny4qJ8JGhqa04iAzER5Zkd59uzZgjvG3DHd3LHYWhtDzRcO\nUTiaPnfuHOfOnQtt+yIzsWnTpoqtW6ERIdkXBgatPvJNDBaa1As7UHIDIgrhUKuKHScv93Ua+uzC\nlxkaLdecRlAansoSxaGaINVHMXJ3Nvl2JkFMtcPJN7Zdj2PS5Z7gDKqU8K/E73ZYFyoGFfSU21oz\n1ec/k8+5LoanzOxmYA/QAjiQdPfHwm1VdeVWH5k/0lKP9Cp9Hn2tmG6HP91ReNhVWxRk/25GURQP\nAutBpEMDGAa2u/tpM7seOGVm33H3s2E3LAy5f6Qa7kkr9oh/ZGSEBx54INI7PJGoinRouPsvgF+M\n/f8rZvYj4K1AQ4ZGrnxHemGfdTRTxQZAqUM1CgyR0tTMnIaZLQGOAre5+8s530sACYCWlpaVe/fu\nLWkbg4ODzJs3b2YNjbC+vj4OHz48YW5hZGSEeDwOwMsvv0x/fz9Xr15leHiYV155ZUbbu/7665k1\n643jktmzZ3PTTTdxww035F2+ubmZO+64g2XLls1ou9Op9885H/W5Mcykz2vXrq2f6zTMbB7wLPCQ\nu//tVMvW20R4JaVSKR5++GHi8fiMK5FyXz9SSY32OYP63CgafiIcwMyagH3AV6cLDJlaOSbSs69O\nz0wUA5EOCREpn0iHhqWfm/o/gB+5+46w21OLZnrKbi1VECJSeZEODWA18KfAGTN7Yey9B9396RDb\nFHkzCQqFhIhMJdKh4e7HAQu7HbWg1KAwM9asWaOQEJFAIh0aMr1ibnYIkyuJ2267jXvuuacKLRWR\neqDQqEHFVhVT3Vq9t7e3wq0VkXqi0KghpdxCXU8AFJFyUmjUgCg8rElEBBQakVZMWFT6sbAiIqDQ\niKQgYaGqQkTCoNCIkKBhoapCRMKi0IiIZDLJ1q1bGRkZyft9hYWIRIFCI2SZ6uLgwYPku3mkwkJE\nokShEZLphqIUFiISRQqNEEw1FGVmrF+/XmEhIpGk0KiyZDJJR0dH3uoiHo+zc+dOEolECC0TEZle\nLOwGNIpUKsXGjRvZsmXLpMCIxWJs2LCBY8eOKTBEJNJUaVRBoeEoDUWJSK1RaFRYoeGoWCzGrl27\nVFmISE3R8FQFFQqMeDyuwBCRmqRKo8xSqRS9vb1cvnyZRx99dEJgaDhKRGqdQqOMkskk9957L8PD\nw5Mu1NNwlIjUA4VGmUx1Kq0CQ0TqheY0yqBQYJgZs2bNUmCISN1QpTFD+QIjFotx3333MX/+fNrb\n2zV/ISJ1Q6FRokI3GtRQlIjUs8iHhpl9BfhXwEvuflvY7YHCF+spMESk3tXCnMZfAx8OuxEZmeEo\nBYaINKLIh4a7HwUGwm4H6GI9EZHIh0ZU5AsMM9ONBkWkoVi+p8VFjZktAb5ZaE7DzBJAAqClpWXl\n3r17S9rO4OAg8+bNm/BeX18fTz31FCdOnJgw4W1mbNu2jXXr1pW0rajI1+d6pz43BvW5OGvXrj3l\n7qumXdDdI/8FLAF+GGTZlStXeqmOHDky4XV3d7fH43EHJnzFYjHv7u4ueTtRktvnRqA+Nwb1uTjA\nSQ+wj4382VNh0d1pRUQmi/ychpk9BaSAd5rZz83szyq9TU14i4jkF/lKw93/dbW21dfXx2OPPTbp\ngj3dnVZEJC3yoVEt999/P1/4whd0d1oRkSkoNEgPR3V1dU16X4EhIjJR5Oc0qmHfvn2T3jMzBYaI\nSA6FBrBp06ZJ77k7zzzzDKlUKoQWiYhEk0IDSCQSdHZ2Tnr/wIEDrFmzhmQyGUKrRESiR6Ex5pFH\nHuHTn/40sdjEf5KRkRG2bNnCxo0bVXWISMNTaGRZt24du3btIh6PT3jf3VV1iIig0JgkkUhw7Ngx\nNmzYgJlN+J6qDhFpdAqNPFpbW9m/fz+7d+9W1SEikkWhMQVVHSIiEyk0phGk6mhra1N4iEhDUGgE\nNFXVMTo6qiErEWkICo0iTFV1gIasRKT+KTRKkF115F7XkRmyuv3229m4cSMbN26ko6NDISIidUGh\nUaJM1XH8+PG8Q1ZDQ0McOHCAAwcOsHv3bs17iEhdUGjM0HRDVhmZeY+2tjZuv/12VR8iUpMUGmWS\nGbLasmULGzZsoKmpKe9yo6OjHD16dLz6UICISC3R8zTKqLW1dfzJfqlUij179nD27FmOHz8+6dGx\n8EaAHD16lGQySVtbG83NzSxatIjNmzfrKYEiJUqlUuzYsYOvfe1rrFixgueff54LFy6E3axABgYG\nOH/+PK+//nqg5Zubm8f7+Oqrr/Lggw9W9JEOCo0KyQ2Qrq4uDh06lDc84I0AyVCISLlkdqCPPfbY\nhPcXLVoU6R1qsTvPjGvXrjEwMFChVkXPhQsXOHv27Pjru+++G6BiwaHQqILMvEeQ6iNjqhDJUJiE\nI/M5RnFHm2tgYGDa3zWpP/v27VNo1INih6+y5YZIRr4wyZYbLLk7vP7+fhYuXFh0X6J+lDqVIH3O\n17+BgQF+/OMf11x/pfHke7BcuZi7V2zlYVi1apWfPHmy6J9LJpN8/vOf5+rVqzNuQ2aM8fTp01y6\ndGna5YeGhrhy5QojIyO4O0NDQzNuQ66mpibcneHh4bKvW0TKb8GCBQVPqMlWrjkNMzvl7qumW27a\nSsPMbgLWAIuAa8DPgLPufr7oVkVUMpkcHwcsh9wxxiioRBCJVEPQnWe2uXPncv3117N69eqaq4hn\nMuzc29tLe3t7+RuVpWBomNks4L8D/548p+aa2c+AbwH/092LP7QPyMw+DDwGxIEvu/tflXsb+/bt\nK/cqRSJvyZIlrF69OnBFXAlmRqHRjrlz57J48eKSTwapxg60EU1Vafxn4G7gPLAfuAjMBT4IvBd4\nC9ABbDGzp4F7yl19mFkceAL4feDnwPfN7JC7l/UwftOmTXz7298u5ypFIu/cuXOcO3cu7GZMKbt9\n083f5erv72f58uU1V2kUq9onxBSc0zCz88Cvgfe7+6tZ7/8l8DngN4E7gE8CfwBcAv7Q3f++bI0z\nawX+k7vfMfb6AQB3f7jQz9TinEYlZc+XAMTjcd70pjcxe/ZsAK5evcqcOXOKXm+U+lisIH3O17/s\nI99i6SwmqaRYLMbNN98MEOqcxkLgq9mBkc3dB4F9wD4zez/wFPAtM1vu7r8ousX5vZX0HErGz4H3\nlWndEyQSCW699daGK2cbsYQPq89hnqqbOWNM4VWfRkdH+elPfwqEe53GOWBJkJW4+3Nm1g68QLoK\n6Zhpw4phZgkgAdDS0kJvb29J6xkcHCz5Z2uV+lxdH//4x+nr6+Pw4cNVvQBtZGSE/v5+AJYtW0Z/\nfz/uzksvvVRwTkFq15NPPsmtt95akXVPFRpPAQ+a2fvd/bnpVuTuPzWzrwN/WLbWwT8BN2e9ftvY\ne7nbTgJJSA9PlXIUWc/DU5nhmHxnoJRjeOry5cs1teMpdnjqtddeK3lYKpeO9KUa7rrrropV01PN\nabwJOEl6p90JPOnuI5k5DXefdEtXM9sBbHH33yhL49JncP0Y+ADpsPg+8Mfu3lfoZ0qZ0yj3Kbci\nUjmxWCzQhHgjTYTfcMMN9PT0cOXKlfDmNNz9NTP7APBN0mcw/YWZPUX6aD/fBm8D/gR4sejWFm7D\nsJndCxwmfcrtV6YKjFLplFtpVPF4nJtvvrno+ztV0lQnHBRzplCjzdc98sgj4V6nAeDuF8bOYPoM\ncN/YlwOY2Q9IH/1fJX3h30rSO/bt5Wyguz8NPF3OdebSKbfSqEZGRiJ52u25c+cKXtR34MCBSe/l\nC5paqTRq7R5ygW8jYmZzgY8CHwHagMU5i/wI+K/u/lRZW1ikWjzlthq3ERGRaCvlyneYuL+pxvBU\nyfeeMrPfIH1K7GzggrtfLGlFZVZqaEB1y9lib1hYyHTju9PdWFA3LMwvu3/PPfccL7zwQpVaJzJz\n3d3dRQdH2e49VYi7XwH+X6k/34hmEhT5wqEcZW2jjftCaX2upduhZwwMDHDs2LGaOrNNykO3Rq9x\nQR7ClCs7JGptzLMeZd/WvpZkwu7MmTOTqquoVovFPnxpaGiIixcjMdARGZW8NbpCo4KKCQuFhFRC\nJuzqvaLMVwnWykR4qU8ozCjXnEZQCo0yS6VS9Pb2cvnyZb74xS+O3/Mpn0xQLF26VCEhMgP5KsF6\nD8p8Qj/lVoLLVBU9PT2Mjo4WHEdWUIhILVNozFDQIahYLMadd95JZ2engkJEapZCo0RBwsLMiMVi\nrFu3TmEhInVBoVGCZDLJ1q1bC85XxONxtm/fzvz582lvb1dYiEjdUGgUIVNdHDx4MO+chYagRKTe\nKTQCmqq6UFiISKNQaExjqurCzFi/fr3CQkQahkJjClNVF/F4nJ07d1b0IhoRkaiJhd2AqEomk3R0\ndEwKDDNjw4YNHDt2TIEhIg1HlUaOqYajVF2ISKNTaGTp6enhsccey1tdaO5CREShMe7+++9nx44d\nk96PxWLs2rVL1YWICJrTANLzF11dXZPej8fjCgwRkSwKDdIPLMllZmzfvl2BISKSRaFB/geWuDuP\nPvooyWQyhBaJiESTQgNIJBJ0dnZOen90dJSOjg4Fh4jIGIXGmEceeYTHH3+cDRs2YGbj74+OjrJl\nyxY2btxIKpUKsYUiIuGLbGiY2R+ZWZ+ZjZrZqmpsc9myZezfv5/du3cTi73xT+PuHDhwgDVr1qjq\nEJGGFtnQAH4IfBQ4Wu0NJxIJdu3aNSE4AEZGRjRcJSINLbLXabj7j4AJQ0XVlDlrKvfeU5l5juxl\nREQahRV6lnVUmFkvcJ+7n5ximQSQAGhpaVm5d+/ekrY1ODjIvHnzJrzX19fHU089xYkTJybcVsTM\n2LZtG+vWrStpW1GRr8/1Tn1uDOpzcdauXXvK3aefCnD30L6A75Iehsr9Wp+1TC+wKug6V65c6aU6\ncuRIwe91d3d7LBZzYPzLzHzDhg1+4sSJkrcZtqn6XK/U58agPhcHOOkB9rGhzmm4+wfd/bY8XwfD\nbFc++eY5XBPkItJgojwRHjmaIBeRRhfZ0DCzjWb2c6AV+JaZHQ67TfBGcMTj8Qnv60JAEWkEkQ0N\nd9/v7m9z9znu3uLud4TdpoxEIsGxY8fyXgio4BCRehbZU26jrrW1lf37948/4W90dBR4IzhefPFF\n5s+fT3t7u57BISJ1Q6ExQ5lrNXKDo6urCzMjHo/zxBNP6JoOEakLkR2eqiWFJsjdneHhYQ1ZiUjd\nUGiUSSY4mpqaJl3FrrkOEakXCo0ySiQSPPvsszz00EN0dnZOqDx0t1wRqQea0yiz1tbW8YnvW265\nZcJcR+ZiwJ6eHnbu3Kl5DhGpOao0KkgXA4pIvVFoVNhUFwNquEpEao1CowoKXQyoe1eJSK1RaFRJ\n5mLA3KcCQnq4SlWHiNQChUaVFRquylQdbW1tCg8RiSyFRggKDVdBeq5D4SEiUaXQCEn2cFVu1QEK\nDxGJJoVGyLKrjty5DlB4iEi0KDQiIFN1HD9+XOEhIpGm0IgQhYeIRJ1CI4KKDY/bb7+djo4OBYiI\nVJxCI8KChsfRo0fZvXu3AkREKk6hUQOChAcoQESk8hQaNSRoeIACREQqQ7dGr0GZ8EilUuzZs4ez\nZ89y/Pjx8Vuw58oEyNGjR0kmk7S1tdHc3MyiRYu47bbbaG9vr24HRKRmKTRqWPazO4oNkAwz4+tf\n//p4iGzevHl8nSIiuSIbGmb2BWAdcA14Efiku18Ot1XRVUqAQPqeV9khkluJKEREJFtkQwP4DvCA\nuw+b2SPAA8D9IbepJpQaIDC5ElGIiEi2yIaGu3876+VzwMfCakstyxcgFy5cYGBgYEYhsnTpUlas\nWMHFixdpb29XkIg0iMiGRo5/B3wt7EbUuuwAgXSIPPzww8Tj8aJDJBMkZkYsFmP16tU0NzcDqCIR\nqWPm7uFt3Oy7wKI83/oLdz84tsxfAKuAj3qBxppZAkgAtLS0rNy7d29J7RkcHGTevHkl/Wwt6Ovr\n4/DhwwwMDIy/NzIyMn6X3Zdffpn+/n6uXr3K8PAwr7zyyoy2d/311zNrVvq4ZPbs2dx0003ccMMN\neZdtbm7mjjvuYNmyZTPaZhD1/jnnoz43hpn0ee3atafcfdV0y4UaGtMxs38L3A18wN2vBPmZVatW\n+cmTJ0vaXm9vb82ffpo9BJUtaCURplgsNj5/EkSpFU09fM7FUp8bw0z6bGaBQiOyw1Nm9mGgE7g9\naGA0mtyAqIVgmEru/EkQ2RP1QejaFJGZiWxoAI8Dc4DvjD3d7jl33xJuk8JX7NlQ9a6UoMm+NqWQ\nRYsWsWLFCp5//vlJVZvmbKSRRTY03P3tYbchCko542kquUNA/f39LFy4sOj1TLVTHRgY4Pz587z+\n+usAXLt2bcI8Sthyr00pRbEVzkzVQlAVGhoNS39/P8uXLy/4e1provI7EOk5jVLUy5xGKpWiq6uL\nQ4cOFR0SheYG8v3SVavP+XYotT6cVm3Fzvlky3dwMFXwF0ufZfUsWLCApqam8dfNzc2sWLGC06dP\nc+XKFR588EESiUTR6635OY1GVMrQU+6OJCpHI7lyT/fNUJgEV8pQnNSfixcvTnh94cIFzp49O/76\n7rvvBigpOIJQaERAsVVFJiiWLl0ayYAoRjFhAtUJlAULFtDS0jJ+9Hbp0qXx7w0NDU36oxWJmn37\n9ik06lHQsMiuJqJaSZRboTCBylcnFy9e5OLFixOO3grJHSqYyty5c1m8eHHRw0u1VHmZGS0tLRVZ\nd9TmxqJs06ZNFVu3QiMkyWSSrVu3MjIyUnCZWCzGnXfeSWdnZ92HRDGKrU6gcjveYquO8+fPFz0v\n0dzcTFtb24STC4p19epV5syZM2m9+aqpmfjlL39Z8xPOlTCT+aiM3BNMMso1pxGUQqPKMtXFwYMH\nyXcSQj0NPVXbVNUJvBEqZ86cmTQpXK2j+SjNS+SOhTe6cuzY86nm6EA1TmxRaFTRVNWFqorKy4RK\noT+s7EplulOKa2W4qB4F3bkXc8ptowz7loNCowqmqi7MjPXr1yssImC6SiVbsdckNFLQVOqIHYrb\nuUfpFPp6otCosKmqi3g8zs6dOys6/iiVUUzAZIR58Vulr9PIXqeO2OubQqOCkskkHR0dk44uVV00\nplKCplx01C3lotAos1QqRW9vL5cvX+bRRx+dFBiqLkSklik0yiiZTHLvvfcyPDysuQsRqUsKjTIp\nNBQF6YnBXbt2qboQkZoXC7sB9WCquYtZs2YpMESkbqjSmKF8gRGLxbjvvvuYP38+7e3tGo4Skbqh\n0JiBQoGhykJE6pWGp0qkwBCRRqTQKIECQ0QalUKjSAoMEWlkCo0iKDBEpNEpNAJSYIiIKDQCUWCI\niKQpNKahwBAReUNkQ8PM/ouZ/cDM/sHM/s7MFle7DQoMEZGJIhsawBfc/d3u/h7gAPCX1dy4AkNE\nZLLIhoa7v5z18jrgYrW2nUqluOeeexQYIiI5LPcW3lFiZg8Bm4HXgPe5+6UCyyWABEBLS8vKvXv3\nlrS9wcFB5s2bx44dO+jp6cleP9u2bWPdunUlrTfKMn1uJOpzY1Cfi7N27dpT7r5q2gXdPbQv4LvA\nD/N8rc9Z7gHgr4Osc+XKlV6qI0eOeHd3t8fjcQcc8Hg87t3d3SWvM+qOHDkSdhOqTn1uDOpzcYCT\nHmAfG+oNC939gwEX/SrwTCXbAtDT08OXvvSl8WEpM+Ouu+7SkJSIyJjIzmmY2TuyXq4HXqjk9u6/\n/3527NgxYR5j1qxZbN68uZKbFRGpKVG+Nfpfmdk7gRHgH4GOSm0omUzS1dU14b1YLMbjjz+uZ2GI\niGSJbGi4+6ZqbWvfvn0TXpuZzpQSEckjssNT1bRp08R8+sxnPqPAEBHJI7KVRjVlAuLJJ5/UxLeI\nyBQUGmMSiQS33nor7e3tYTdFRCSyNDwlIiKBKTRERCQwhYaIiASm0BARkcAUGiIiEphCQ0REAov0\nrdFLYWb9wE9L/PEbgV+VsTm1QH1uDOpzY5hJn/+Zuy+cbqG6C42ZMLOTHuR+8nVEfW4M6nNjqEaf\nNTwlIiKBKTRERCQwhcZEybAbEAL1uTGoz42h4n3WnIaIiASmSkNERAJryNAwsw+b2f81s5+Y2Wfz\nfN/M7L+Nff8HZvbbYbSznAL0+V1mljKzq2Z2XxhtLLcAff43Y5/vGTM7YWbvCaOd5RKgv+vH+vuC\nmZ02sw+E0c5ymq7PWcv9jpkNm9nHqtm+SgjwObeb2a/HPucXzOxzZW2AuzfUFxAHXgT+OTAb+Adg\nac4yHwGeAQx4P/D3Ybe7Cn2+Cfgd4CHgvrDbXKU+/y7w5rH//4Na/pwD9ncebwxJvxt4Mex2V7rP\nWcv9HfA08LGw212Fz7kd+Gal2tCIlcZ7gZ+4+z+6+zVgL7A+Z5n1wB5Pew6Yb2ZvqXZDy2jaPrv7\nS+7+fWAojAZWQJA+n3D3S2MvnwPeVuU2llOQ/g762F4FuA64WOU2lluQv2WAPwf2AS9Vs3EVErTP\nFdOIofFW4GdZr38+9l6xy9SSeutPEMX2+c9IV5e1KlB/zWyjmf0f4H8D/6FKbauUaftsZm8FNgK7\nqtiuSgr6e/27Y0ORz5jZsnI2QE/uk4ZnZmtJh0Zb2G2pNHffD+w3s38J7DGzd7n7aNjtqqAvAfe7\n+6iZhd2WajkNLHb3QTP7CHAAeEe5Vt6IlcY/ATdnvX7b2HvFLlNL6q0/QQTqs5m9G/gysN7da3m4\npqjP2N2Pkj5oXFDhdlVSkD6vAvaa2TngY8BOM9tQneZVxLR9dveX3X1w7P+fBprM7MZyNaARQ+P7\nwDvM7LfMbDbwCeBQzjKHgM1jZ1G9H/i1u/+i2g0toyB9rjfT9tnMFgN/C/ypu/84hDaWU5D+vt3G\nDrfHzgg0d++vflPLZto+u/tvufsSd18CfAPY6u4Hqt/UsgnyOS/K+pzfS3o/X7YDooYbnnL3YTO7\nFzhM+kwkLFhMAAACKklEQVSEr7h7n5ltGfv+btJnWXwE+AlwBfhkWO0thyB9NrNFwEngBmDUzP4j\n6bMyXg6t4TMQ8HP+HOkj7Z1jf2PDXqM3uAvY302kD4aGgFdJ73BqVsA+15WAff4Y0GFmw8BrwCey\nToCYMV0RLiIigTXi8JSIiJRIoSEiIoEpNEREJDCFhoiIBKbQEBGRwBQaIhVkZtvMzM3sj8Nui0g5\nKDREKmvl2H9PhtoKkTLRdRoiFWRmZ0nfUG5+OS+wEgmLKg2RCjGz64B3As8rMKReKDREKudfkP4b\nO5X9ppm92cwOjs11fMnMmsJpnkjxGu7eUyJVlHlM8HhomNn7gK8Bbwb+yN2/EUbDREqlSkOkcjKT\n4KcAzOzTwDHgZWCVAkNqkSoNkcr5beAV4FdmdhC4E/gboMPdXwu1ZSIl0tlTIhVgZnNJB8ZLwDBw\nE/Dn7v7lUBsmMkManhKpjPeQruTnAIuBbygwpB4oNEQqIzMJvh34JvAnYw+2EqlpCg2Rysi+EvwT\npCfDv2hmHw2vSSIzpzkNkQows9PAu4Dr3X1k7HG6z5Ge2/g9d38u1AaKlEiVhkiZmdls4DbgjLuP\nALj7BdLPnb8K9JjZ20NsokjJFBoi5bccaAKez37T3c8CHwV+E3jGzG4MoW0iM6LhKRERCUyVhoiI\nBKbQEBGRwBQaIiISmEJDREQCU2iIiEhgCg0REQlMoSEiIoEpNEREJDCFhoiIBKbQEBGRwP4/f73k\nhpQAODwAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb88a7bd080>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot(kran,theta,'k.');\n",
"grid()\n",
"ylabel(r'$\\theta$',fontsize=20)\n",
"xlabel(r'$k$',fontsize=20);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"kr = dot(dh.rcell, (0,0,0))\n",
"phases = exp(-1j * dot(kr, dot(dh.cell, dh.sc.sc_off.T))) \n",
"offsets = - np.arange(len(phases)) * dh.no\n",
"dia = diags(phases, offsets, shape=(dh.shape[1], dh.shape[0]), dtype=complex128)\n",
"V = (dh.tocsr(dh.S_idx)).todense().dot(dia)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 0.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 1.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.],\n",
" [ 2.]])"
]
},
"execution_count": 103,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"transpose(kron([arange(3)],ones((dh.no))))*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from scipy.sparse.linalg import "
]
},
{
"cell_type": "code",
"execution_count": 90,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-18.15016975+0.j, -14.84379374+0.j, 1.48413093+0.j,\n",
" -11.34543432+0.j, 1.48425660+0.j, 1.81764548+0.j,\n",
" -4.39790634+0.j, 1.81770364+0.j, 11.43780014+0.j,\n",
" 11.71093358+0.j, 5.45676251+0.j, 11.71151583+0.j,\n",
" 11.43778311+0.j, -31.73407295+0.j, -28.55014816+0.j,\n",
" 4.35919994+0.j, 6.46936701+0.j, 4.35885880+0.j,\n",
" 5.91290510+0.j, 5.79933095+0.j, 5.91252297+0.j,\n",
" 10.22921143+0.j, 14.20176378+0.j, 13.02755238+0.j,\n",
" 14.20169975+0.j, 10.22830190+0.j])"
]
},
"execution_count": 90,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dh.tocsr(0)*kron(phases,ones(dh.no))"
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x7fb888053e80>"
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAOfCAYAAADfPeYuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGCtJREFUeJzt3V+s7Wl91/HPt8MwOFBxRprJKaCAGTTYIOhJU23TkBAt\nVuPgzQSSmtGQjBfYUGOi0Jv2hoSa2rQX2mQEdBqRZkJRJk3TkY5tqjGhPfzJ8GeEmVRGhs4fCipU\nkukAjxdnUU9Pz5nzmb32Pmuts1+vm733Wmuv/cyC886z9v7m+c1aKwBX8h27XgBwGMQCqIgFUBEL\noCIWQEUsgMrexWJm3jAzn52Zh2fm7btezzOZmc/PzCdn5hMzc27X67nQzLx3Zp6cmU9dcNvNM/Ph\nmXlo8/GmXa7x2y6z1p+cmS9uXttPzMwP73KN3zYzL52ZX5+Zz8zMp2fmbZvb9+61fYa1Hum1nX2a\ns5iZ65J8LslfT/Jokt9O8ua11md2urDLmJnPJzm71vq9Xa/lYjPzg0l+P8kvrLW+Z3PbP0/ylbXW\nuzYhvmmt9c92uc7Nui611p9M8vtrrZ/e5douNjNnkpxZa31sZr4zyUeTvDHJ38+evbbPsNbbc4TX\ndt92Ft+b5OG11u+stf4gyS8muW3HazpIa63fTPKVi26+Lcndm8/vzvn/4+zcZda6l9Zaj621Prb5\n/GtJHkzy4uzha/sMaz2SfYvFi5N84YKvH80W/3FXwUryazPz0Zm5c9eLKdyy1nps8/njSW7Z5WIK\nPzozD2zepux8W3+xmXlZktcm+Uj2/LW9aK3JEV7bfYvFofmBtdZrkvzNJG/dbKcPwjr//nN/3oP+\ncT+f5BVJXpPksST/YrfL+aNm5gVJfinJj621vnrhffv22l5irUd6bfctFl9M8tILvn7J5ra9tNb6\n4ubjk0n+Q86/jdpnT2zex377/eyTO17PZa21nlhrfXOt9a0k/zp79NrOzPU5/4/vfWutD25u3svX\n9lJrPepru2+x+O0kt87My2fmuUnelOTeHa/pkmbm+ZtfGmVmnp/kbyT51DN/187dm+SOzed3JPnQ\nDtfyjL79D2/j72ZPXtuZmSTvSfLgWutnLrhr717by631qK/tXv01JEk2f8b52STXJXnvWuudO17S\nJc3MK3J+N5Ekz0ny7/dprTPz/iSvS/KiJE8k+Ykk/zHJPUn+TJJHkty+1tr5LxYvs9bX5fw2eSX5\nfJJ/eMHvBHZmZn4gyX9J8skk39rc/OM5/7uAvXptn2Gtb84RXtu9iwWwn/btbQiwp8QCqIgFUBEL\noCIWQGVvY3Eg49NJDmeth7LOxFpPyjZr3dtYJDmY/wFyOGs9lHUm1npSDjsWM/Oru14DnFbtv78T\nG8qamTck+bmcn8R891rrXZd77Au/8zvWra947h+57Utf/ma+609f94dff+6BG09kncfh6TyV63PD\nrpdxRYeyzsRaT8ql1vq1/K+vrrVeeKXvfc5JLGhziM2/zAWH2MzMvZc7xObWVzw3v3XfSy911x/6\noe9+zbGvE0h+bX3goeZxJ/U2xCE2cI05qVgc2iE2wBXs7BecM3PnzJybmXNf+vI3d7UMoHRSsbji\nITZrrbvWWmfXWmcv/EUmsJ9OKhYHc4gN0DmRv4astb4xM/8oyX35/4fYfPokfhZwdZxILJJkrfUr\nSX6leeznHrjxin8ave93P/GM9/vTKpysvZjgBPafWAAVsQAqYgFUxAKoiAVQEQugcmJzFsdt2zmM\n5jmAy7OzACpiAVTEAqiIBVARC6AiFkBFLICKWACVgxnKupJm4MoBOnB0dhZARSyAilgAFbEAKmIB\nVMQCqIgFULlm5iwaDtCBo7OzACpiAVTEAqiIBVARC6AiFkBFLICKWACVUzWUdSXHcYBO+zxwaOws\ngIpYABWxACpiAVTEAqiIBVARC6AiFkDFUNaz5MpnnFZ2FkBFLICKWAAVsQAqYgFUxAKoiAVQMWdx\nAlz5jGuRnQVQEQugIhZARSyAilgAFbEAKmIBVMQCqBjK2gFXPuMQ2VkAFbEAKmIBVMQCqIgFUBEL\noCIWQMWcxZ5yMSP2jZ0FUBELoCIWQEUsgIpYABWxACpiAVTEAqgYyjpgrnzG1WRnAVTEAqiIBVAR\nC6AiFkBFLICKWAAVcxbXMAfocJzsLICKWAAVsQAqYgFUxAKoiAVQEQugIhZAxVDWKecAHVp2FkBF\nLICKWAAVsQAqYgFUxAKoiAVQEQugYiiLZ3Qcp221z8N+s7MAKmIBVMQCqIgFUBELoCIWQEUsgIo5\nC7bmymeng50FUBELoCIWQEUsgIpYABWxACpiAVTEAqgYyuKqcOWzw2dnAVTEAqiIBVARC6AiFkBF\nLICKWAAVcxbsBQfo7D87C6AiFkBFLICKWAAVsQAqYgFUxAKoiAVQMZTFwXCAzm7ZWQAVsQAqYgFU\nxAKoiAVQEQugIhZARSyAiqEsrhnHcdpW+zynkZ0FUBELoCIWQEUsgIpYABWxACpiAVTMWXCquPLZ\n0dlZABWxACpiAVTEAqiIBVARC6AiFkBFLICKoSy4iCufXdqRdxYz89KZ+fWZ+czMfHpm3ra5/eaZ\n+fDMPLT5eNPxLRfYlW3ehnwjyT9Za70qyfcleevMvCrJ25Pcv9a6Ncn9m6+BA3fkWKy1HltrfWzz\n+deSPJjkxUluS3L35mF3J3njtosEdu9YfmcxMy9L8tokH0lyy1rrsc1djye55TLfc2eSO5Pkebnx\nOJYBnKCt/xoyMy9I8ktJfmyt9dUL71trrSTrUt+31rprrXV2rXX2+tyw7TKAE7ZVLGbm+pwPxfvW\nWh/c3PzEzJzZ3H8myZPbLRHYB9v8NWSSvCfJg2utn7ngrnuT3LH5/I4kHzr68oB9sc3vLL4/yd9L\n8smZ+fYfnn88ybuS3DMzb0nySJLbt1si7JfTejGjI8dirfVfk8xl7n79UZ8X2E/GvYGKWAAVsQAq\nYgFUxAKoiAVQEQug4vAbOAHX4pXP7CyAilgAFbEAKmIBVMQCqIgFUBELoGLOAnbk0C5mZGcBVMQC\nqIgFUBELoCIWQEUsgIpYABWxACqGsmBP7dsBOnYWQEUsgIpYABWxACpiAVTEAqiIBVARC6BiKAsO\n2HGctnXdme5n2VkAFbEAKmIBVMQCqIgFUBELoCIWQGUv5ixe+eqv5777rt4hHnBadP9uHq6ey84C\nqIgFUBELoCIWQEUsgIpYABWxACpiAVT2Yijrcw/cuPUhHoa24GTZWQAVsQAqYgFUxAKoiAVQEQug\nIhZAZS/mLBrHcTEVsxhwdHYWQEUsgIpYABWxACpiAVTEAqiIBVARC6ByMENZV9IMXDlAB47OzgKo\niAVQEQugIhZARSyAilgAFbEAKmIBVK6ZoayG07bg6OwsgIpYABWxACpiAVTEAqiIBVARC6ByquYs\nruQ4DtBpnwcOjZ0FUBELoCIWQEUsgIpYABWxACpiAVTEAqgYynqWXPmM08rOAqiIBVARC6AiFkBF\nLICKWAAVsQAq5ixOgIsZcS2yswAqYgFUxAKoiAVQEQugIhZARSyAilgAFUNZO+DKZxwiOwugIhZA\nRSyAilgAFbEAKmIBVMQCqIgFUDGUtadc+Yx9Y2cBVMQCqIgFUBELoCIWQEUsgIpYABVzFgfMlc+4\nmuwsgIpYABWxACpiAVTEAqiIBVARC6AiFkDFUNY1zAE6HCc7C6AiFkBFLICKWAAVsQAqYgFUxAKo\nmLM45RygQ8vOAqiIBVARC6AiFkBFLICKWAAVsQAqYgFUDGXxjI7jAJ32edhvdhZARSyAilgAFbEA\nKmIBVMQCqIgFUDFnwdZczOh0sLMAKlvHYmaum5mPz8wvb76+eWY+PDMPbT7etP0ygV07jp3F25I8\neMHXb09y/1rr1iT3b74GDtxWsZiZlyT5W0nefcHNtyW5e/P53UneuM3PAPbDtjuLn03yT5N864Lb\nbllrPbb5/PEkt1zqG2fmzpk5NzPnns5TWy4DOGlHjsXM/O0kT661Pnq5x6y1VpJ1mfvuWmudXWud\nvT43HHUZwFWyzZ9Ovz/J35mZH07yvCR/cmb+XZInZubMWuuxmTmT5MnjWCiwW0feWay13rHWesla\n62VJ3pTkP6+1fiTJvUnu2DzsjiQf2nqVwM6dxFDWu5LcMzNvSfJIkttP4GdwYFz57PAdSyzWWr+R\n5Dc2n385yeuP43mB/WGCE6iIBVARC6AiFkBFLICKWAAVsQAqTspiLzhta//ZWQAVsQAqYgFUxAKo\niAVQEQugIhZAxZwFB8MBOrtlZwFUxAKoiAVQEQugIhZARSyAilgAFbEAKoayuGYcxwE67fOcRnYW\nQEUsgIpYABWxACpiAVTEAqiIBVAxZ8Gp4mJGR2dnAVTEAqiIBVARC6AiFkBFLICKWAAVsQAqhrLg\nIq58dml2FkBFLICKWAAVsQAqYgFUxAKoiAVQEQugYigLnqXTeuUzOwugIhZARSyAilgAFbEAKmIB\nVMQCqJizgBNwLV75zM4CqIgFUBELoCIWQEUsgIpYABWxACpiAVQMZcGOHNqVz+wsgIpYABWxACpi\nAVTEAqiIBVARC6BizgL21L4doGNnAVTEAqiIBVARC6AiFkBFLICKWAAVsQAqhrLggB3HATrXnel+\nlp0FUBELoCIWQEUsgIpYABWxACpiAVT2Ys7ila/+eu677+od4gGnRffv5uHquewsgIpYABWxACpi\nAVTEAqiIBVARC6AiFkBlL4ayPvfAjVsf4mFoC06WnQVQEQugIhZARSyAilgAFbEAKmIBVMQCqOzF\nUFbjOK68ZHALjs7OAqiIBVARC6AiFkBFLICKWAAVsQAqBzNncSXNDIUDdODo7CyAilgAFbEAKmIB\nVMQCqIgFUBELoCIWQOWaGcpqOEAHjs7OAqiIBVARC6AiFkBFLICKWAAVsQAqp2rO4kqO4wCd9nng\n0NhZABWxACpiAVTEAqiIBVARC6AiFkBFLICKoaxnyZXPOK3sLICKWAAVsQAqYgFUxAKoiAVQEQug\nIhZAxVDWCXDlM65FdhZARSyAilgAFbEAKmIBVMQCqIgFUDFnsQOufMYhsrMAKlvFYmb+1Mx8YGb+\n+8w8ODN/dWZunpkPz8xDm483Hddigd3Zdmfxc0l+da31F5L8pSQPJnl7kvvXWrcmuX/zNXDgjhyL\nmXlhkh9M8p4kWWv9wVrrfye5Lcndm4fdneSN2y4S2L1tdhYvT/KlJP9mZj4+M++emecnuWWt9djm\nMY8nueVS3zwzd87MuZk593Se2mIZwNWwTSyek+QvJ/n5tdZrk/zfXPSWY621kqxLffNa66611tm1\n1tnrc8MWywCuhm1i8WiSR9daH9l8/YGcj8cTM3MmSTYfn9xuicA+OHIs1lqPJ/nCzPz5zU2vT/KZ\nJPcmuWNz2x1JPrTVCoG9sO1Q1o8med/MPDfJ7yT5BzkfoHtm5i1JHkly+5Y/41Ry5TP2zVaxWGt9\nIsnZS9z1+m2eF9g/JjiBilgAFbEAKmIBVMQCqIgFUHH4zQFzMSOuJjsLoCIWQEUsgIpYABWxACpi\nAVTEAqiIBVAxlHUNc4AOx8nOAqiIBVARC6AiFkBFLICKWAAVsQAqYgFUDGWdck7bomVnAVTEAqiI\nBVARC6AiFkBFLICKWAAVcxY8o+M4QKd9HvabnQVQEQugIhZARSyAilgAFbEAKmIBVMQCqBjKYmuu\nfHY62FkAFbEAKmIBVMQCqIgFUBELoCIWQMWcBVeFixkdPjsLoCIWQEUsgIpYABWxACpiAVTEAqiI\nBVAxlMVecIDO/rOzACpiAVTEAqiIBVARC6AiFkBFLICKOQsOhgN0dsvOAqiIBVARC6AiFkBFLICK\nWAAVsQAqYgFUDGVxzTiOA3Ta5zmN7CyAilgAFbEAKmIBVMQCqIgFUBELoCIWQMVQFqeKK58dnZ0F\nUBELoCIWQEUsgIpYABWxACpiAVTMWcBFXPns0uwsgIpYABWxACpiAVTEAqiIBVARC6AiFkDFUBY8\nS6f1ymd2FkBFLICKWAAVsQAqYgFUxAKoiAVQMWcBJ+BavJiRnQVQEQugIhZARSyAilgAFbEAKmIB\nVMQCqBjKgh05tCuf2VkAFbEAKmIBVMQCqIgFUBELoCIWQEUsgIqhLNhT+3balp0FUBELoCIWQEUs\ngIpYABWxACpiAVTMWcABO44DdK470/0sOwugIhZARSyAilgAFbEAKmIBVMQCqIgFUNmLoaxXvvrr\nue++q3eIB5wW3b+bh6vnsrMAKmIBVMQCqIgFUBELoCIWQEUsgMpezFl87oEbtz7EwxwGnCw7C6Ai\nFkBFLICKWAAVsQAqYgFUxAKoiAVQ2YuhrMZxXHnJ4BYc3VY7i5l5x8x8ZmY+NTPvn5nnzczNM/Ph\nmXlo8/Gm41ossDtHjsXMvCzJnUn+ylrre5Jcl+RNSd6e5P611q1J7t98DRy4bXYWX03ydJI/MTPP\nSXJjkt9NcluSuzePuTvJG7daIbAXjhyLtdZXkvx0kv+Z5LEk/2et9Z+S3LLWemzzsMeT3LL1KoGd\n2+ZtyJ9L8o+TvDzJdyd5/sz8yIWPWWutJOsy33/nzJybmXNP56mjLgO4SrZ5G3I2yX9ba31prfV0\nkg8m+WtJnpiZM0my+fjkpb55rXXXWuvsWuvs9blhi2UAV8M2sfhsku+bmRtnZpK8PsmDSe5Ncsfm\nMXck+dB2SwT2wZHnLNZan5iZX0hyLsm3knw8yV1JXpDknpl5S5JHktx+HAu9kmaGwgE6cHRbDWWt\ntX4qyU9ddPNTOb/LAK4hxr2BilgAFbEAKmIBVMQCqIgFUBELoHIwh98cBwfowNHZWQAVsQAqYgFU\nxAKoiAVQEQugIhZARSyAyqkayrqS4zhtq30eODR2FkBFLICKWAAVsQAqYgFUxAKoiAVQMWfxLLny\nGaeVnQVQEQugIhZARSyAilgAFbEAKmIBVMQCqBjKOgGufMa1yM4CqIgFUBELoCIWQEUsgIpYABWx\nACrmLHbAxYw4RHYWQEUsgIpYABWxACpiAVTEAqiIBVARC6BiKGtPufIZ+8bOAqiIBVARC6AiFkBF\nLICKWAAVsQAqYgFUDGUdMFc+42qyswAqYgFUxAKoiAVQEQugIhZARSyAijmLa5gDdDhOdhZARSyA\nilgAFbEAKmIBVMQCqIgFUBELoGIo65RzgA4tOwugIhZARSyAilgAFbEAKmIBVMQCqJiz4BkdxwE6\n7fOw3+wsgIpYABWxACpiAVTEAqiIBVARC6AiFkDFUBZbc+Wz08HOAqiIBVARC6AiFkBFLICKWAAV\nsQAqYgFUDGVxVbjy2eGzswAqYgFUxAKoiAVQEQugIhZARSyAijkL9oIDdPafnQVQEQugIhZARSyA\nilgAFbEAKmIBVMQCqBjK4mA4QGe37CyAilgAFbEAKmIBVMQCqIgFUBELoGLOgmvGcRyg0z7PaWRn\nAVTEAqiIBVARC6AiFkBFLICKWAAVsQAqhrI4VVz57OjsLICKWAAVsQAqYgFUxAKoiAVQEQugYs4C\nLuJiRpdmZwFUxAKoiAVQEQugIhZARSyAilgAFbEAKoay4Fk6rVc+s7MAKmIBVMQCqIgFUBELoCIW\nQEUsgIpYABVDWXACrsUrn11xZzEz752ZJ2fmUxfcdvPMfHhmHtp8vOmC+94xMw/PzGdn5odOauHA\n1dW8Dfm3Sd5w0W1vT3L/WuvWJPdvvs7MvCrJm5L8xc33/KuZue7YVgvszBVjsdb6zSRfuejm25Lc\nvfn87iRvvOD2X1xrPbXW+h9JHk7yvce0VmCHjvoLzlvWWo9tPn88yS2bz1+c5AsXPO7RzW1/zMzc\nOTPnZubc03nqiMsArpat/xqy1lpJ1hG+76611tm11tnrc8O2ywBO2FFj8cTMnEmSzccnN7d/MclL\nL3jcSza3AQfuqLG4N8kdm8/vSPKhC25/08zcMDMvT3Jrkt/abonAPrjinMXMvD/J65K8aGYeTfIT\nSd6V5J6ZeUuSR5LcniRrrU/PzD1JPpPkG0neutb65gmtHQ7aoV357IqxWGu9+TJ3vf4yj39nkndu\nsyhg/xj3BipiAVTEAqiIBVARC6AiFkBFLICKw29gT+3bATp2FkBFLICKWAAVsQAqYgFUxAKoiAVQ\nMWcBB+w4DtC57kz3s+wsgIpYABWxACpiAVTEAqiIBVARC6AiFkBlzl/XeMeLmPlSzl/Z7EIvSvJ7\nO1jOURzKWg9lnYm1npRLrfXPrrW+60rfuBexuJSZObfWOrvrdTQOZa2Hss7EWk/KNmv1NgSoiAVQ\n2edY3LXrBTwLh7LWQ1lnYq0n5chr3dvfWQD7ZZ93FsAeEQugIhZARSyAilgAlf8HV9T1ZC7FyYoA\nAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fb88a296438>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"matshow(abs(dia.todense())[:2*52,:])"
]
},
{
"cell_type": "code",
"execution_count": 114,
"metadata": {},
"outputs": [],
"source": [
"A=rand(7,3)\n",
"B=arange(3).reshape(1,3)"
]
},
{
"cell_type": "code",
"execution_count": 115,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0. , 0.12041409, 1.87048978],\n",
" [ 0. , 0.48879153, 1.47919414],\n",
" [ 0. , 0.81748285, 1.06503386],\n",
" [ 0. , 0.87313179, 0.80500327],\n",
" [ 0. , 0.9514684 , 1.81707067],\n",
" [ 0. , 0.92036652, 1.35496943],\n",
" [ 0. , 0.74585022, 1.20294771]])"
]
},
"execution_count": 115,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"B*A"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.4.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment