add charge density wave order parameter to example file

docs/source/graphene_example.md

@@ -26,7 +26,9 @@ import numpy as np
 import matplotlib.pyplot as plt
 import pymf.model as model
 import pymf.solvers as solvers
+import pymf.mf as mf
 import pymf.tb as tb
+import pymf.observables as observables
 import pymf.kwant_helper.kwant_examples as kwant_examples
 import pymf.kwant_helper.utils as kwant_utils
 ```
@@ -52,21 +54,124 @@ U=1
 V=0.1
 params = dict(U=U, V=V)
 h_int = kwant_utils.builder_to_tb(int_builder, params)
-model = model.Model(h_0, h_int, filling=2)
+_model = model.Model(h_0, h_int, filling=2)
 ```
 
 To start the meanfield calculation we also need a starting guess. We will use our random guess generator for this. It creates a random Hermitian hopping dictionary based on the hopping keys provided and the number of degrees of freedom specified. As we don't expect the mean-field solution to contain terms more than the hoppings from the interacting part, we can use the hopping keys from the interacting part. We will use the same number of degrees as freedom as both the non-interacting and interacting part, so that they match.
 
 ```{code-cell} ipython3
 guess = tb.utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
-mf_sol = solvers.solver(model, guess, nk=18)
+mf_sol = solvers.solver(_model, guess, nk=18, optimizer_kwargs={'M':0})
 full_sol = tb.tb.add_tb(h_0, mf_sol)
 ```
 
 After we have defined the guess, we feed it together with the model into the meanfield solver. The meanfield solver will return a hopping dictionary with the meanfield approximation. We can then add this solution to the non-interacting part to get the full solution. In order to get the solution, we specified the number of k-points to be used in the calculation. This refers to the k-grid used in the Brillouin zone for the density matrix.
 
-## Creating a phase diagram
+## Creating a phase diagram of the gap
 
-We can now create a phase diagram of the gap of the interacting solution. We will use the same hopping dictionary for the non-interacting part as before. We will vary the onsite Hubbard interactio $U$ strength from $0$ to $2$ and the nearest neighbor interaction strength $V$ from $0$ to $1.5$.
+We can now create a phase diagram of the gap of the interacting solution. In order to calculate the gap we first create a function which takes a hopping dictionary and a Fermi energy and returns the indirect gap. The gap is defined as the difference between the highest occupied and the lowest unoccupied energy level. We will use a dense k-grid to calculate the gap. In order to obtain the Hamiltonian on a dense k-grid, we use the `tb_to_khamvector` function from the `transforms` module.
 
 ```{code-cell} ipython3
+def compute_gap(h, fermi_energy=0, nk=100):
+    kham = tb.transforms.tb_to_khamvector(h, nk, ks=None)
+    vals = np.linalg.eigvalsh(kham)
+
+    emax = np.max(vals[vals <= fermi_energy])
+    emin = np.min(vals[vals > fermi_energy])
+    return np.abs(emin - emax)
+```
+
+Now that we can calculate the gap, we create a phase diagram of the gap as a function of the Hubbard interaction strength $U$ and the nearest neighbor interaction strength $V$. We vary the onsite Hubbard interactio $U$ strength from $0$ to $2$ and the nearest neighbor interaction strength $V$ from $0$ to $1.5$.
+
+```{code-cell} ipython3
+def compute_phase_diagram(Us, Vs, int_builder, h_0):
+    gap = []
+    mf_sols = []
+    for U in Us:
+      gap_U = []
+      mf_sols_U = []
+      for V in Vs:
+        params = dict(U=U, V=V)
+        h_int = kwant_utils.builder_to_tb(int_builder, params)
+        _model = model.Model(h_0, h_int, filling=2)
+        converged=False
+        while not converged:
+          guess = tb.utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
+          try:
+            mf_sol = solvers.solver(_model, guess, nk=18, optimizer_kwargs={'M':0})
+            converged=True
+          except:
+            converged=False
+        gap_U.append(compute_gap(tb.tb.add_tb(h_0, mf_sol), fermi_energy=0, nk=300))
+        mf_sols_U.append(mf_sol)
+      guess = None
+      gap.append(gap_U)
+      mf_sols.append(mf_sols_U)
+    return np.asarray(gap, dtype=float), np.asarray(mf_sols)
+```
+We chose to initialize a new guess for each $U$ value, but not for each $V$ value. Instead, for consecutive $V$ values we use the previous mean-field solution as a guess. We do this because the mean-field solution is expected to be smooth in the interaction strength and therefore by using an inspired guess we can speed up the calculation.
+
+We can now compute the phase diagram and plot it.
+
+```{code-cell} ipython3
+Us = np.linspace(0, 4, 5)
+Vs = np.linspace(0, 1.5, 5)
+gap, mf_sols = compute_phase_diagram(Us, Vs, int_builder, h_0)
+plt.imshow(gap.T, extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto')
+plt.colorbar()
+plt.xlabel('V')
+plt.ylabel('U')
+plt.title('Gap')
+plt.show()
+```
+
+This phase diagram has gap openings at the same places as shown in the [literature](https://arxiv.org/abs/1204.4531).
+
+## Order parameters
+
+We might also want to calculate order parameters of the mean field solution. From literature we know to expect a charge density wave (CDW) and a spin density wave (SDW) in the mean field solution in different regions of the phase diagram. Here we show how to create both the CDW and SDW order parameters and evaluate them for the mean field solution in the phase diagram.
+
+### Charge density wave
+
+We first start with the CDW order parameter. In CDW the spins on different sublattices have opposite signs. This means that our order parameter can be constructed as:
+
+```{code-cell} ipython3
+sz = np.array([[1, 0], [0, -1]])
+cdw_order_parameter = {}
+cdw_order_parameter[(0,0)] = np.kron(sz, np.eye(2))
+```
+
+We choose a point in the phase diagram where we expect there to be a CDW phase and calculate the expectation value with the CDW order parameter. In order to do this we first construct the density matrix from the mean field solution.
+
+```{code-cell} ipython3
+params = dict(U=0, V=2)
+h_int = kwant_utils.builder_to_tb(int_builder, params)
+_model = model.Model(h_0, h_int, filling=2)
+guess = tb.utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
+mf_sol = solvers.solver(_model, guess, nk=18, optimizer_kwargs={'M':0})
+full_sol = tb.tb.add_tb(h_0, mf_sol)
+
+rho, _ = mf.construct_density_matrix(full_sol, filling=2, nk=40)
+expectation_value = observables.expectation_value(rho, cdw_order_parameter)
+print(expectation_value)
+```
+
+We can also perform the same calculation over the complete phase diagram where we calculated the gap earlier:
+
+```{code-cell} ipython3
+expectation_value_list = []
+for mf_sol in mf_sols.flatten():
+    rho, _ = mf.construct_density_matrix(tb.tb.add_tb(h_0, mf_sol), filling=2, nk=40)
+    expectation_value = observables.expectation_value(rho, cdw_order_parameter)
+    expectation_value_list.append(expectation_value)
+```
+
+```{code-cell} ipython3
+expectation_value_list = np.asarray(expectation_value_list).reshape(mf_sols.shape)
+plt.imshow(np.abs(expectation_value_list.T.real), extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto')
+plt.colorbar()
+plt.xlabel('V')
+plt.ylabel('U')
+plt.title('Charge Density Wave Order Parameter')
+plt.show()
+```

examples/graphene_extended_hubbard.ipynb

@@ -2,32 +2,24 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 11,
    "id": "cb509096-42c6-4a45-8dc4-a8eed3116e67",
    "metadata": {
     "tags": []
    },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
-    "from codes.model import Model\n",
-    "from codes.solvers import solver\n",
-    "from codes import kwant_examples\n",
-    "from codes.kwant_helper import utils as kwant_utils\n",
-    "from codes.tb.tb import add_tb\n",
-    "from codes.tb import utils\n",
-    "from tqdm import tqdm\n",
-    "import xarray as xr\n"
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import pymf.model as model\n",
+    "import pymf.solvers as solvers\n",
+    "import pymf.mf as mf\n",
+    "import pymf.tb as tb\n",
+    "import pymf.observables as observables\n",
+    "import pymf.kwant_helper.kwant_examples as kwant_examples\n",
+    "import pymf.kwant_helper.utils as kwant_utils"
    ]
   },
   {
@@ -35,64 +27,231 @@
    "execution_count": 2,
    "id": "99f0e60c",
    "metadata": {},
+   "outputs": [],
+   "source": [
+    "graphene_builder, int_builder = kwant_examples.graphene_extended_hubbard()\n",
+    "h_0 = kwant_utils.builder_to_tb(graphene_builder)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 86,
+   "id": "9189a20f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "U=4\n",
+    "V=0\n",
+    "params = dict(U=U, V=V)\n",
+    "h_int = kwant_utils.builder_to_tb(int_builder, params)\n",
+    "_model = model.Model(h_0, h_int, filling=2)\n",
+    "guess = tb.utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))\n",
+    "mf_sol = solvers.solver(_model, guess, nk=18, optimizer_kwargs={'M':0})\n",
+    "full_sol = tb.tb.add_tb(h_0, mf_sol)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 87,
+   "id": "9d6ce961",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sx = np.array([[0, 1], [1, 0]])\n",
+    "sy = np.array([[0, -1j], [1j, 0]])\n",
+    "sz = np.array([[1, 0], [0, -1]])\n",
+    "\n",
+    "s_list = [sx, sy, sz]\n",
+    "rho, _ = mf.construct_density_matrix(full_sol, filling=2, nk=40)\n",
+    "\n",
+    "order_parameter_list = []\n",
+    "for s in s_list: \n",
+    "    order_parameter = {}\n",
+    "    order_parameter[(0,0)] =  np.kron(sz, s)\n",
+    "    order_parameter_list.append(order_parameter)\n",
+    "\n",
+    "expectation_value_list = []\n",
+    "for order_parameter in order_parameter_list:\n",
+    "    expectation_value = observables.expectation_value(rho, order_parameter)\n",
+    "    expectation_value_list.append(expectation_value)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 88,
+   "id": "d77b4bea",
+   "metadata": {},
    "outputs": [
     {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
-      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n"
-     ]
+     "data": {
+      "text/plain": [
+       "(1.804655008699922+0j)"
+      ]
+     },
+     "execution_count": 88,
+     "metadata": {},
+     "output_type": "execute_result"
     }
    ],
    "source": [
-    "# Create translationally-invariant `kwant.Builder`\n",
-    "graphene_builder, int_builder = kwant_examples.graphene_extended_hubbard()\n",
-    "h_0 = kwant_utils.builder_to_tb(graphene_builder)"
+    "np.sum(np.array(expectation_value_list)**2)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
-   "id": "f4d1bb07",
+   "execution_count": 64,
+   "id": "10dfae00",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#antiferromagnetic order parameter\n",
+    "anti_order_parameter = {}\n",
+    "anti_order_parameter[(0,0)] = np.kron(sz, sz)\n",
+    "expectation_value = observables.expectation_value(rho, anti_order_parameter)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7d3394a1",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 54,
+   "id": "260bd7c1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "cdw_order_parameter = {}\n",
+    "cdw_order_parameter[(0,0)] = np.kron(sz, np.eye(2))\n",
+    "\n",
+    "\n",
+    "sdw_order_parameter = {}\n",
+    "sdw_order_parameter[(0,0)] = np.kron(np.eye(2), sz)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 55,
+   "id": "345644dc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "expectation_value = observables.expectation_value(rho, cdw_order_parameter)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 56,
+   "id": "fe4b9564",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(-1.5833996732759665+0j)"
+      ]
+     },
+     "execution_count": 56,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "expectation_value"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "83e823ed",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{(0,\n",
+       "  0): array([[0., 0., 1., 0.],\n",
+       "        [0., 0., 0., 1.],\n",
+       "        [1., 0., 0., 0.],\n",
+       "        [0., 1., 0., 0.]]),\n",
+       " (0,\n",
+       "  -1): array([[0., 0., 1., 0.],\n",
+       "        [0., 0., 0., 1.],\n",
+       "        [0., 0., 0., 0.],\n",
+       "        [0., 0., 0., 0.]]),\n",
+       " (0,\n",
+       "  1): array([[0., 0., 0., 0.],\n",
+       "        [0., 0., 0., 0.],\n",
+       "        [1., 0., 0., 0.],\n",
+       "        [0., 1., 0., 0.]]),\n",
+       " (1,\n",
+       "  -1): array([[0., 0., 1., 0.],\n",
+       "        [0., 0., 0., 1.],\n",
+       "        [0., 0., 0., 0.],\n",
+       "        [0., 0., 0., 0.]]),\n",
+       " (-1,\n",
+       "  1): array([[0., 0., 0., 0.],\n",
+       "        [0., 0., 0., 0.],\n",
+       "        [1., 0., 0., 0.],\n",
+       "        [0., 1., 0., 0.]])}"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "h_0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e1894ac3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def compute_gap(h, fermi_energy=0, nk=100):\n",
+    "    kham = tb.transforms.tb_to_khamvector(h, nk, ks=None)\n",
+    "    vals = np.linalg.eigvalsh(kham)\n",
+    "\n",
+    "    emax = np.max(vals[vals <= fermi_energy])\n",
+    "    emin = np.min(vals[vals > fermi_energy])\n",
+    "    return np.abs(emin - emax)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d0092977",
    "metadata": {},
    "outputs": [],
    "source": [
-    "np.random.seed(5)\n",
     "def compute_phase_diagram(Us, Vs, int_builder, h_0): \n",
     "    gap = []\n",
-    "    for U in tqdm(Us): \n",
-    "        gap_U = []\n",
-    "        guess=None\n",
-    "        for V in Vs: \n",
-    "            params = dict(U=U, V=V)\n",
-    "            h_int = kwant_utils.builder_to_tb(int_builder, params)\n",
-    "            if guess==None:\n",
-    "                guess = utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))\n",
-    "            model = Model(h_0, h_int, filling=2)\n",
+    "    for U in Us: \n",
+    "      for V in Vs: \n",
+    "        params = dict(U=U, V=V)\n",
+    "        h_int = kwant_utils.builder_to_tb(int_builder, params)\n",
+    "        _model = model.Model(h_0, h_int, filling=2)\n",
     "\n",
-    "            mf_sol = solver(model, guess, nk=18)    \n",
-    "            gap_U.append(utils.compute_gap(add_tb(h_0, mf_sol), fermi_energy=0, nk=300))\n",
-    "            guess = None\n",
-    "        gap.append(gap_U)\n",
-    "    return np.asarray(gap, dtype=float)"
+    "        converged=False\n",
+    "        while not converged:\n",
+    "          guess = tb.utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))\n",
+    "          try:\n",
+    "            mf_sol = solvers.solver(_model, guess, nk=18, optimizer_kwargs={'M':0})\n",
+    "            converged=True\n",
+    "          except:\n",
+    "            converged=False\n",
+    "        gap.append(compute_gap(tb.tb.add_tb(h_0, mf_sol), fermi_energy=0, nk=300))\n",
+    "      guess = None\n",
+    "    return np.asarray(gap, dtype=float).reshape(len(Us), len(Vs))"
    ]
   },
   {
@@ -369,7 +528,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.6"
+   "version": "3.9.13"
   },
   "widgets": {
    "application/vnd.jupyter.widget-state+json": {