stablization_green_function.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.insert(0, \"./admf\")\n",
    "from admf import utils\n",
    "from collections import namedtuple\n",
    "import numpy as np\n",
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "basis = namedtuple(\"basis\", [\"x\", \"y\", \"sub\", \"spin\"])\n",
    "\n",
    "\n",
    "def generate_lattice(nx, ny):\n",
    "    for i in range(nx):\n",
    "        for j in range(ny):\n",
    "            yield basis(i, j, 0, 0)\n",
    "            yield basis(i, j, 0, 1)\n",
    "            yield basis(i, j, 1, 0)\n",
    "            yield basis(i, j, 1, 1)\n",
    "\n",
    "\n",
    "def nn(t):\n",
    "    if t.sub == 0:  # A sublattice\n",
    "        yield basis((t.x - 1), (t.y - 1), 1, t.spin)\n",
    "        yield basis(t.x, (t.y - 1), 1, t.spin)\n",
    "        yield basis(t.x, t.y, 1, t.spin)\n",
    "    elif t.sub == 1:  # B sublattice\n",
    "        yield basis(t.x, (t.y + 1), 0, t.spin)\n",
    "        yield basis((t.x + 1), (t.y + 1), 0, t.spin)\n",
    "        yield basis(t.x, t.y, 0, t.spin)\n",
    "\n",
    "\n",
    "pnn = utils.pbc(nn)\n",
    "\n",
    "\n",
    "def real_position(t):\n",
    "    return (\n",
    "        np.sqrt(3) * t.x - np.sqrt(3) / 2.0 * t.y,\n",
    "        -3 / 2 * t.y if t.sub == 0 else -3 / 2 * t.y - 1,\n",
    "    )\n",
    "\n",
    "\n",
    "def rashba(t, nx, ny, lmbd=1):\n",
    "    if t.spin == 0:\n",
    "        sigma = np.array([1, -1j, 0])\n",
    "    else:\n",
    "        sigma = np.array([1, 1j, 0])\n",
    "    for site in nn(t):\n",
    "        psite = utils.site_mod(site, {\"x\": nx, \"y\": ny})\n",
    "        dsite = utils.spin_flip(psite)\n",
    "        dx, dy = real_position(site)\n",
    "        ux, uy = real_position(t)\n",
    "        d = np.array([dx - ux, dy - uy, 0.0])\n",
    "        cr = np.cross(sigma, d)[2]\n",
    "        cr = cr / np.linalg.norm(cr)\n",
    "        yield dsite, 1j * lmbd * cr\n",
    "\n",
    "\n",
    "nx = 10\n",
    "ny = 10\n",
    "loc, _ = utils.loc_index(generate_lattice(nx, ny))\n",
    "uloc, _ = utils.loc_index(generate_lattice(nx, ny), lambda t: t.spin == 0)\n",
    "nloc = uloc\n",
    "\n",
    "hsize = len(loc)\n",
    "K, RS = utils.generate_np_zeros(2, [hsize, hsize])\n",
    "for site in loc:\n",
    "    for hopsite in pnn(site, {\"x\": nx, \"y\": ny}):\n",
    "        K[loc[site], loc[hopsite]] = 1\n",
    "    for rashbasite, lam in rashba(site, nx, ny):\n",
    "        RS[loc[site], loc[rashbasite]] = lam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta= 10\n",
    "dt = 0.1\n",
    "l = int(beta/dt)\n",
    "u = 5.\n",
    "lbd = np.arccosh(np.exp(u*dt/2))\n",
    "expK = linalg.expm(dt*(K+RS))\n",
    "V = [np.zeros([len(loc)]) for _ in range(int(l))]\n",
    "B = []\n",
    "for p in range(l):\n",
    "    s = 2*np.random.choice(2, size=[nx*ny*2])-1\n",
    "    for i, j in loc.items():\n",
    "        V[p][j]  = s[int(j/2)]*(-1)**(i.spin)*u\n",
    "    B.append([email protected](np.exp(lbd*V[p])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def decomp(b):\n",
    "    u = np.eye(b.shape[0])\n",
    "    s = np.eye(b.shape[0])\n",
    "    v = b\n",
    "    return u, s, v\n",
    "\n",
    "def decompsvd(b):\n",
    "    u, s, v = np.linalg.svd(b)\n",
    "    return u, np.diag(s), v\n",
    "\n",
    "def decompsdd(b):\n",
    "    u, s, v = linalg.svd(b, lapack_driver=\"gesdd\")\n",
    "    return u, np.diag(s), v\n",
    "\n",
    "def decompqr(b, epsilon=1e-20):\n",
    "    q, r = np.linalg.qr(b)\n",
    "    d = np.diagonal(r)\n",
    "    di = 1/(d+epsilon)\n",
    "    di = np.diag(di)\n",
    "    d = np.diag(d)\n",
    "    return q, d, di@r\n",
    "\n",
    "def decompqrp(a, epsilon=1e-20):\n",
    "    q, r, p = linalg.qr(a, pivoting=True)\n",
    "    pm = np.zeros_like(a)\n",
    "    for i in range(a.shape[0]):\n",
    "        pm[i,p[i]] = 1\n",
    "    d = np.diagonal(r)\n",
    "    di = 1/(d+epsilon)\n",
    "    di = np.diag(di)\n",
    "    d = np.diag(d)\n",
    "    return q, d, di@r@pm\n",
    "\n",
    "def ipbbb(blist, decomp_func=None):\n",
    "    if decomp_func is None:\n",
    "        decomp_func = decomp\n",
    "    assert len(blist)>1\n",
    "    u, s, v = decomp(blist[-1])\n",
    "    for i in range(1, len(blist)):\n",
    "        c = blist[len(blist)-1-i]@u@s\n",
    "        u, s, v1 = decomp_func(c)\n",
    "        v = v1@v\n",
    "    # now we have u,s,v\n",
    "    s = np.diagonal(s)\n",
    "    sbi = np.array([1/d if np.abs(d)>1 else 1 for  d in s])\n",
    "    sbi = np.diag(sbi)\n",
    "    ss = np.array([d if np.abs(d)<1 else 1 for  d in s])\n",
    "    ss = np.diag(ss)\n",
    "    return np.linalg.inv([email protected]().T+ss@v)@([email protected]().T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "qprr = ipbbb(B, decomp_func=decompqr)\n",
    "diff = (ipbbb(B, decomp_func=decompqrp)-qprr)/np.abs(qprr) #ipbbb(B, decomp_func=decompqr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[ 0.00327751-0.07621173j,  0.23255039+0.86073044j,\n",
       "         -1.09127022+0.38705612j, ..., -0.99984207+0.00999757j,\n",
       "          0.68929645+0.48991378j,  0.31080438-0.94581675j],\n",
       "        [ 0.61068991+1.16481611j, -0.93011994-0.23367805j,\n",
       "          0.1795815 -0.35238937j, ...,  1.63284139-0.64540915j,\n",
       "         -0.85272525+0.40671605j,  0.80246317+0.16163303j],\n",
       "        [ 0.32620772+0.94536451j, -0.69073084+0.53413475j,\n",
       "         -0.43050278+0.73567789j, ..., -0.23464304-0.9723987j ,\n",
       "          0.54961247+0.51750164j,  0.47679298+0.87899484j],\n",
       "        ...,\n",
       "        [-2.75035903+0.72614842j,  0.98565085-0.26911691j,\n",
       "          0.75486867+0.59059479j, ..., -0.38206473+0.98490994j,\n",
       "          0.44419098-0.86624267j, -2.60086089-0.88171352j],\n",
       "        [ 0.07048176+0.01052458j, -0.90702722-0.03783568j,\n",
       "         -1.25981596+0.23289341j, ...,  1.2629737 -0.78354849j,\n",
       "         -0.65694429+0.54969214j,  0.89856469+0.2001996j ],\n",
       "        [-0.97452756-0.33171337j,  0.63220686+0.63930236j,\n",
       "          0.30182817+0.9993272j , ..., -0.95364993+0.5139401j ,\n",
       "          0.85257681+0.08520754j, -0.0133771 -0.0107126j ]]),\n",
       " 1.792867175421553)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff, np.mean(np.abs(diff)) # 20*20 beta=10 u=5 somehow limit "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-3.31800416e-06-2.43204422e-06j,  1.07604692e-08+6.91113510e-08j,\n",
       "         -1.85638989e-07+9.54964823e-08j, ...,\n",
       "          2.04886583e-07-1.85569368e-07j, -2.20339915e-07+3.66557481e-07j,\n",
       "          4.12637342e-08-7.98851427e-08j],\n",
       "        [-3.52849323e-07-1.61343408e-07j,  3.74116111e-09-8.19471210e-09j,\n",
       "         -1.01368446e-04-4.60879058e-05j, ...,\n",
       "         -3.98001360e-08+3.86474491e-07j,  8.58753627e-08+7.86400984e-08j,\n",
       "         -6.53330394e-06-1.14839258e-06j],\n",
       "        [-8.09229067e-09-4.25138724e-07j,  7.61287637e-05+1.05034578e-04j,\n",
       "         -1.30726799e-07+5.73249756e-07j, ...,\n",
       "          1.43239758e-08-1.87375640e-07j,  2.20046370e-06-1.54310598e-06j,\n",
       "         -1.09899831e-05+1.50880080e-05j],\n",
       "        ...,\n",
       "        [ 5.37209299e-07+2.22343549e-06j,  8.69493899e-08-1.06047557e-07j,\n",
       "          4.21118837e-07-9.31489547e-07j, ...,\n",
       "         -2.24120219e-07-6.79408078e-07j,  1.02939737e-06+3.89361215e-07j,\n",
       "         -1.64997515e-08+9.47170600e-09j],\n",
       "        [ 1.23971820e-06+2.90745017e-07j, -8.13779117e-08-2.59132089e-07j,\n",
       "         -1.67473522e-06-2.78830782e-06j, ...,\n",
       "         -5.13326999e-08-1.06007799e-07j,  7.79926911e-08-2.30592409e-08j,\n",
       "          2.33846465e-08+7.20346250e-09j],\n",
       "        [ 5.65984950e-07-1.08178808e-06j, -4.43792191e-07+4.04380344e-08j,\n",
       "         -1.06144914e-06-2.01351299e-07j, ...,\n",
       "         -4.38212715e-08+3.48440005e-08j,  1.20185501e-07-3.66768239e-07j,\n",
       "          3.56915195e-09+4.75803666e-09j]]), 8.4692823384249e-05)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff, np.mean(np.abs(diff)) # 10*10 beta=10 u=5 somehow limit "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}