change the optimizer kwargs to make convergence behave better

06d1a482 · Johanna Zijderveld · 6d791558 · 06d1a482 · 06d1a482
Commit 06d1a482 authored May 7, 2024 by Johanna Zijderveld
--- a/docs/source/hubbard_1d.md
+++ b/docs/source/hubbard_1d.md
@@ -91,7 +91,7 @@ def compute_sol(U, h_0, nk, filling=2):
    }
    guess = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
    full_model = pymf.Model(h_0, h_int, filling)
-    mf_sol = pymf.solver(full_model, guess, nk=nk)
+    mf_sol = pymf.solver(full_model, guess, nk=nk, optimizer_kwargs={"M":0})
    return pymf.add_tb(h_0, mf_sol)

--- a/examples/1d_hubbard_totalenergy.ipynb
+++ b/examples/1d_hubbard_totalenergy.ipynb
@@ -124,7 +124,7 @@
    "    }\n",
    "    guess = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))\n",
    "    full_model = pymf.Model(h_0, h_int, filling)\n",
-    "    mf_sol = pymf.solver(full_model, guess, nk=nk)\n",
+    "    mf_sol = pymf.solver(full_model, guess, nk=nk, optimizer_kwargs={\"M\":0})\n",
    "    return pymf.add_tb(h_0, mf_sol)\n",
    "\n",
    "\n",

 %% Cell type:code id:cb509096-42c6-4a45-8dc4-a8eed3116e67 tags:
 ``` python
 import numpy as np
 import matplotlib.pyplot as plt
 import pymf
 from tqdm import tqdm
 ```
 %% Cell type:markdown id:396d935c-146e-438c-878b-04ed70aa6d63 tags:
 To simulate infinite systems, we provide the corresponding tight-binding model.
 We exemplify this construction by computing the ground state of an infinite spinful chain with onsite interactions.
 Because the ground state is an antiferromagnet, so we must build a two-atom cell. We name the two sublattices, $A$ and $B$. The Hamiltonian in is:
 $$
 H_0 = \sum_i c_{i, B}^{\dagger}c_{i, A} + c_{i, A}^{\dagger}c_{i+1, B} + h.c.
 $$
 We write down the spinful by simply taking $H_0(k) \otimes \mathbb{1}$.
 %% Cell type:code id:5529408c-fb7f-4732-9a17-97b0718dab29 tags:
 ``` python
 hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))
 h_0 = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}
 ```
 %% Cell type:markdown id:050f5435-6699-44bb-b31c-8ef3fa2264d4 tags:
 To build the tight-binding model, we need to generate a Hamiltonian on a k-point and the corresponding hopping vectors to generate a guess. We then verify the spectrum and see that the bands indeed consistent of two bands due to the Brillouin zone folding.
 %% Cell type:code id:b39a2976-7c35-4670-83ef-12157bd3fc0e tags:
 ``` python
 # Set number of k-points
 nk = 100
 ks = np.linspace(0, 2*np.pi, nk, endpoint=False)
 hamiltonians_0 = pymf.tb_to_khamvector(h_0, nk, ks=ks)
 # Perform diagonalization
 vals, vecs = np.linalg.eigh(hamiltonians_0)
 # Plot data
 plt.plot(ks, vals, c="k")
 plt.xticks([0, np.pi, 2 * np.pi], ["$0$", "$\pi$", "$2\pi$"])
 plt.xlim(0, 2 * np.pi)
 plt.ylabel("$E - E_F$")
 plt.xlabel("$k / a$")
 plt.show()
 ```
 %% Output
 %% Cell type:markdown id:6ec53b08-053b-4aad-87a6-525dd7f61687 tags:
 Here, in the workflow to find the ground state, we use a helper function to build the initial guess. because we don't need a dense k-point grid in the self-consistent loop, we compute the spectrum later on a denser k-point grid.
 %% Cell type:markdown id:dc59e440-1289-4735-9ae8-b04d0d13c94a tags:
 Finally, we compute the eigen0alues for a set of Ualues of $U$. For this case, since the interaction is onsite only, the interaction matrix is simply
 $$
 H_{int} =
 \left(\begin{array}{cccc}
    U & U & 0 & 0\\
    U & U & 0 & 0\\
    0 & 0 & U & U\\
    0 & 0 & U & U
 \end{array}\right)~.
 $$
 %% Cell type:code id:32b9e7c5-db12-44f9-930c-21e5494404b8 tags:
 ``` python
 def compute_sol(U, h_0, nk, filling=2):
    h_int = {
        (0,): U * np.kron(np.eye(2), np.ones((2, 2))),
    }
    guess = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
    full_model = pymf.Model(h_0, h_int, filling)
-    mf_sol = pymf.solver(full_model, guess, nk=nk)
+    mf_sol = pymf.solver(full_model, guess, nk=nk, optimizer_kwargs={"M":0})
    return pymf.add_tb(h_0, mf_sol)
 def compute_gap_and_vals(full_sol, nk_dense, fermi_energy=0):
    h_kgrid = pymf.tb_to_khamvector(full_sol, nk_dense)
    vals = np.linalg.eigvalsh(h_kgrid)
    emax = np.max(vals[vals <= fermi_energy])
    emin = np.min(vals[vals > fermi_energy])
    return np.abs(emin - emax), vals
 def compute_phase_diagram(
    Us,
    nk,
    nk_dense,
 ):
    gaps = []
    vals = []
    for U in tqdm(Us):
        full_sol = compute_sol(U, h_0, nk)
        gap, _vals = compute_gap_and_vals(full_sol, nk_dense)
        gaps.append(gap)
        vals.append(_vals)
    return np.asarray(gaps, dtype=float), np.asarray(vals)
 ```
 %% Cell type:code id:6a8c08a9-7e31-420b-b6b4-709abfb26793 tags:
 ``` python
 # Interaction strengths
 Us = np.linspace(0.5, 10, 20, endpoint=True)
 nk, nk_dense = 40, 100
 gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)
 ```
 %% Output
     75%|███████▌  | 15/20 [00:02<00:00,  5.07it/s]/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=9.07651e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=8.54968e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=8.17567e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=8.88621e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=9.35564e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=9.73285e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
     85%|████████▌ | 17/20 [00:06<00:01,  2.47it/s]
    ---------------------------------------------------------------------------
    NoConvergence                             Traceback (most recent call last)
    /var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/1169016709.py in <module>
          2 Us = np.linspace(0.5, 10, 20, endpoint=True)
          3 nk, nk_dense = 40, 100
    ----> 4 gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)
    /var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/3682270308.py in compute_phase_diagram(Us, nk, nk_dense)
         26     vals = []
         27     for U in tqdm(Us):
    ---> 28         full_sol = compute_sol(U, h_0, nk)
         29         gap, _vals = compute_gap_and_vals(full_sol, nk_dense)
         30         gaps.append(gap)
    /var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/3682270308.py in compute_sol(U, h_0, nk, filling)
          5     guess = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
          6     full_model = pymf.Model(h_0, h_int, filling)
    ----> 7     mf_sol = pymf.solver(full_model, guess, nk=nk)
          8     return pymf.add_tb(h_0, mf_sol)
          9
    ~/kwant-scf/pymf/solvers.py in solver(Model, mf_guess, nk, optimizer, optimizer_kwargs)
         60     f = partial(cost, Model=Model, nk=nk)
         61     result = rparams_to_tb(
    ---> 62         optimizer(f, mf_params, **optimizer_kwargs), list(Model.h_int), shape
         63     )
         64     fermi = calculate_fermi_energy(add_tb(Model.h_0, result), Model.filling, nk=nk)
    ~/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py in anderson(F, xin, iter, alpha, w0, M, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, **kw)
    ~/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py in nonlin_solve(F, x0, jacobian, iter, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, full_output, raise_exception)
        239     else:
        240         if raise_exception:
    --> 241             raise NoConvergence(_array_like(x, x0))
        242         else:
        243             status = 2
    NoConvergence: [ 1.86715628e+00  4.47287535e-02 -1.23356829e-04  1.20430729e-03
  4.47287535e-02  2.63543104e+00  2.63234161e-03  3.71702031e-03
 -1.23356829e-04  2.63234161e-03  1.86002862e+00  4.51677226e-02
  1.20430729e-03  3.71702031e-03  4.51677226e-02  2.64245378e+00
  0.00000000e+00  1.61192759e-01 -1.51104882e-03  1.47113079e-03
 -1.61192759e-01  0.00000000e+00  5.03808388e-03  1.99828384e-03
  1.51104882e-03 -5.03808388e-03  0.00000000e+00  1.64032886e-01
 -1.47113079e-03 -1.99828384e-03 -1.64032886e-01  0.00000000e+00]
 %% Cell type:code id:e17fc96c-c463-4e1f-8250-c254d761b92a tags:
 ``` python
 import xarray as xr
 ds = xr.Dataset(
    data_vars=dict(vals=(["Us", "ks", "n"], vals), gap=(["Us"], gap)),
    coords=dict(
        Us=Us,
        ks=np.linspace(0, 2 * np.pi, nk_dense),
        n=np.arange(vals.shape[-1])
    ),
 )
 ```
 %% Cell type:markdown id:5a87dcc1-208b-4602-abad-a870037ec95f tags:
 We observe that as the interaction strength increases, a gap opens due to the antiferromagnetic ordering.
 %% Cell type:code id:50f02d06 tags:
 ``` python
 plt.scatter(Us, gap)
 ```
 %% Output
    <matplotlib.collections.PathCollection at 0x7fe0a154a130>
 %% Cell type:code id:868cf368-45a0-465e-b042-6182ff8b6998 tags:
 ``` python
 ds.vals.plot.scatter(x="ks", hue="Us", ec=None, s=5)
 plt.axhline(0, ls="--", c="k")
 plt.xticks([0, np.pi, 2 * np.pi], ["$0$", "$\pi$", "$2\pi$"])
 plt.xlim(0, 2 * np.pi)
 plt.ylabel("$E - E_F$")
 plt.xlabel("$k / a$")
 plt.show()
 ```
 %% Output
    ---------------------------------------------------------------------------
    AttributeError                            Traceback (most recent call last)
    /var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/3492086229.py in <module>
    ----> 1 ds.vals.plot.scatter(x="ks", hue="Us", ec=None, s=5)
          2 plt.axhline(0, ls="--", c="k")
          3 plt.xticks([0, np.pi, 2 * np.pi], ["$0$", "$\pi$", "$2\pi$"])
          4 plt.xlim(0, 2 * np.pi)
          5 plt.ylabel("$E - E_F$")
    AttributeError: '_PlotMethods' object has no attribute 'scatter'
 %% Cell type:markdown id:0761ed33-c1bb-4f12-be65-cb68629f58b9 tags:
 The Hartree-Fock dispersion should follow (see [these notes](https://www.cond-mat.de/events/correl11/manuscript/Lechermann.pdf))
 $$
 \epsilon_{HF}^{\sigma}(\mathbf{k}) = \epsilon(\mathbf{k}) + U \left(\frac{n}{2} + \sigma m\right)
 $$
 where $m=(\langle n_{i\uparrow} \rangle - \langle n_{i\downarrow} \rangle) / 2$ is the magnetization per atom and $n = \sum_i \langle n_i \rangle$ is the total number of atoms per cell. Thus, for the antiferromagnetic groundstate, $m=1/2$ and $n=2$. The gap thus should be $\Delta=U$. And we can confirm it indeed follows the expected trend.
 %% Cell type:code id:ac2eb725-f3bd-4d5b-a509-85d0d0071958 tags:
 ``` python
 ds.gap.plot()
 plt.plot(ds.Us, ds.Us)
 plt.show()
 ```
 %% Output
 %% Cell type:markdown id:06e0d356-558e-40e3-8287-d7d2e0bee8cd tags:
 We can also fit
 %% Cell type:code id:5499ea62-cf1b-4a13-8191-ebb73ea38704 tags:
 ``` python
 ds.gap.polyfit(dim="Us", deg=1).polyfit_coefficients[0].data
 ```
 %% Output
    array(0.9997852)
 %% Cell type:code id:0cb395cd-84d1-49b4-89dd-da7a2d09c8d0 tags:
 ``` python
 ds.to_netcdf("./data/1d_hubbard_example.nc")
 ```
 %% Cell type:code id:ce428241 tags:
 ``` python
 ```{code-cell} ipython3
 def compute_phase_diagram(
    Us,
    nk,
    nk_dense,
    filling=2,
 ):
    gap = []
    vals = []
    for U in tqdm(Us):
        # onsite interactions
        guess = utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
        full_model = Model(h_0, h_int, filling)
        mf_sol = solver(full_model, guess, nk=nk)
        hkfunc = transforms.tb_to_kfunc(add_tb(h_0, mf_sol))
        ks_dense = np.linspace(0, 2 * np.pi, nk_dense, endpoint=False)
        hkarray = np.array([hkfunc(kx) for kx in ks_dense])
        _vals = np.linalg.eigvalsh(hkarray)
        _gap = (utils.compute_gap(add_tb(h_0, mf_sol), fermi_energy=0, nk=nk_dense))
        gap.append(_gap)
        vals.append(_vals)
    return np.asarray(gap, dtype=float), np.asarray(vals)
 import xarray as xr
 ds = xr.Dataset(
    data_vars=dict(vals=(["Us", "ks", "n"], vals), gap=(["Us"], gap)),
    coords=dict(
        Us=Us,
        ks=np.linspace(0, 2 * np.pi, nk_dense),
        n=np.arange(vals.shape[-1])
    ),
 )
 # Interaction strengths
 Us = np.linspace(0.5, 10, 20, endpoint=True)
 nk, nk_dense = 40, 100
 gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)
 ds.vals.plot.scatter(x="ks", hue="Us", ec=None, s=5)
 plt.axhline(0, ls="--", c="k")
 plt.xticks([0, np.pi, 2 * np.pi], ["$0$", "$\pi$", "$2\pi$"])
 plt.xlim(0, 2 * np.pi)
 plt.ylabel("$E - E_F$")
 plt.xlabel("$k / a$")
 plt.show()
 ```
 The Hartree-Fock dispersion should follow (see [these notes](https://www.cond-mat.de/events/correl11/manuscript/Lechermann.pdf))
 $$
 \epsilon_{HF}^{\sigma}(\mathbf{k}) = \epsilon(\mathbf{k}) + U \left(\frac{n}{2} + \sigma m\right)
 $$
 where $m=(\langle n_{i\uparrow} \rangle - \langle n_{i\downarrow} \rangle) / 2$ is the magnetization per atom and $n = \sum_i \langle n_i \rangle$ is the total number of atoms per cell. Thus, for the antiferromagnetic groundstate, $m=1/2$ and $n=2$. The gap thus should be $\Delta=U$. And we can confirm it indeed follows the expected trend.
 ```