reorganize folders

a8cbb668 · Antonio Manesco · 38d18e7b · a8cbb668 · a8cbb668 · a8cbb668
Commit a8cbb668 authored 1 year ago by Antonio Manesco
--- a/codes/coarse_k_point_implementation.ipynb
+++ b/codes/coarse_k_point_implementation.ipynb
--- a/analysis/codes
+++ b/analysis/codes
+../codes/
\ No newline at end of file
--- a/codes/data/graphene_example.nc
+++ b/codes/data/graphene_example.nc
--- a/codes/graphene_test_extended_hubbard.ipynb
+++ b/codes/graphene_test_extended_hubbard.ipynb
--- a/codes/mean_field_demo.ipynb
+++ b/codes/mean_field_demo.ipynb
--- a/codes/scipy_optimizer.ipynb
+++ b/codes/scipy_optimizer.ipynb
@@ -308,11 +308,24 @@
  },
  {
   "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 32,
   "id": "e17fc96c-c463-4e1f-8250-c254d761b92a",
   "metadata": {},
   "outputs": [],
-   "source": []
+   "source": [
+    "import xarray as xr\n",
+    "gap_da = xr.DataArray(data=gap, coords=dict(Us=Us, Vs=Vs))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 33,
+   "id": "0cb395cd-84d1-49b4-89dd-da7a2d09c8d0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gap_da.to_netcdf('./data/graphene_example.nc')"
+   ]
  }
 ],
 "metadata": {

 %% Cell type:code id:cb509096-42c6-4a45-8dc4-a8eed3116e67 tags:

 ``` python
 import kwant
 import numpy as np
 import matplotlib.pyplot as plt
 import utils, hf
 from scipy.optimize import anderson
 from tqdm import tqdm

 s0 = np.identity(2)
 sz = np.diag([1, -1])

 norbs=2

 graphene = kwant.lattice.general(
    [[1, 0], [1 / 2, np.sqrt(3) / 2]], [[0, 0], [0, 1 / np.sqrt(3)]],
    norbs=norbs
 )
 a, b = graphene.sublattices

 # create bulk system
 bulk_graphene = kwant.Builder(kwant.TranslationalSymmetry(*graphene.prim_vecs))
 # add sublattice potential
 m0 = 0
 bulk_graphene[a.shape((lambda pos: True), (0, 0))] = m0 * sz
 bulk_graphene[b.shape((lambda pos: True), (0, 0))] = -m0 * sz
 # add hoppings between sublattices
 bulk_graphene[graphene.neighbors(1)] = s0

 # use kwant wraparound to sample bulk k-space
 wrapped_syst = kwant.wraparound.wraparound(bulk_graphene)
 wrapped_fsyst = kwant.wraparound.wraparound(bulk_graphene).finalized()
 ```

 %% Cell type:code id:9cc3b32d-404f-4bc5-a338-83571c9e4c4b tags:

 ``` python
 def func_onsite(site, U):
    return U * np.ones((2, 2))

 def func_hop(site1, site2, V):
    rij = np.linalg.norm(site1.pos - site2.pos)
    return V * np.ones((2, 2))

 def calculate_Hint(U, V, Uk, Vk):
    return U * Uk + V * Vk
 ```

 %% Cell type:code id:a341e0e5-330e-48d1-a20f-a0040688a9d7 tags:

 ``` python
 nk = 10
 # Generate coarse-grid k-points
 ks, dk = np.linspace(0, 2 * np.pi, nk, endpoint=False, retstep=True)
 # Generate Hamiltonian on a k-point grid
 hamiltonians_0 = utils.syst2hamiltonian(ks=ks, syst=wrapped_fsyst)
 ```

 %% Cell type:code id:d1ef154e-70bd-4f28-887f-72362d8533dd tags:

 ``` python
 Us = np.linspace(0, 4, 50)
 Vs = np.linspace(0, 1.5, 20)
 ```

 %% Cell type:code id:32b9e7c5-db12-44f9-930c-21e5494404b8 tags:

 ``` python
 _, deltas = utils.generate_scf_syst(
    max_neighbor=1, syst=wrapped_syst, lattice=graphene
 )
 deltas = np.asarray(deltas) #deltas are the hopping vecs
 deltas = np.unique(np.stack([*deltas, *-deltas]), axis=(0))

 def compute_gap(
    U,
    V,
    H_int,
    max_neighbor=1,
    lattice=graphene,
    filling=2,
    nk=12,
    tol=1e-5,
    norbs=norbs,
    nk_dense=30,
    mixing=0.5,
    order=1,
    guess=None
 ):
    # Generate coarse-grid k-points
    ks, dk = np.linspace(0, 2 * np.pi, nk, endpoint=False, retstep=True)
    # Generate Hamiltonian on a k-point grid
    hamiltonians_0 = utils.syst2hamiltonian(ks=ks, syst=wrapped_fsyst)
    # Generate guess on the same grid
    if guess is None:
        guess = guess = utils.generate_guess(nk, 2, deltas, ndof=4, scale=1)
    else:
        # guess *= 0.25
        guess += np.max(guess) * utils.generate_guess(nk, 2, deltas, ndof=4, scale=0.1)
        # guess = guess

    # Find groundstate Hamiltonian on the same grid
    hk = hf.find_groundstate_ham(
        H_int=H_int,
        nk=nk,
        filling=filling,
        hamiltonians_0=hamiltonians_0,
        tol=tol,
        guess=guess,
        mixing=mixing,
        order=order,
    )
    # Diagonalize groundstate Hamiltonian
    vals, vecs = np.linalg.eigh(hk)
    # Extract coarse-grid Fermi energy
    E_F = utils.get_fermi_energy(vals, 2)
    # Generate kwant system with k-grid Hamiltonian
    scf_syst = utils.hk2syst(deltas, hk, ks, dk, max_neighbor, norbs, lattice)
    # Generate dense-grid k-points
    ks_dense = np.linspace(0, 2 * np.pi, nk_dense, endpoint=False)
    # Compute groundstate Hamiltonian on a dense grid
    scf_ham = utils.syst2hamiltonian(
        ks=ks_dense, syst=scf_syst, params={}
    )
    # Diagonalize groundstate Hamiltonian
    vals, vecs = np.linalg.eigh(scf_ham)
    # Extract dense-grid Fermi energy
    E_F = utils.get_fermi_energy(vals, 2)

    gap = utils.calc_gap(vals, E_F)
    return gap, hk


 def compute_phase_diagram(Us, Vs, nk, tol, mixing, order):
    import qsymm
    import adaptive
    import utils, hf

    ks = np.linspace(0, 2 * np.pi, nk, endpoint=False)

    Uk = utils.potential2hamiltonian(
        syst=wrapped_syst,
        lattice=graphene,
        func_onsite=func_onsite,
        func_hop=func_hop,
        params=dict(U=1, V=0),
        ks=ks,
    )

    Vk = utils.potential2hamiltonian(
        syst=wrapped_syst,
        lattice=graphene,
        func_onsite=func_onsite,
        func_hop=func_hop,
        params=dict(U=0, V=1),
        ks=ks,
    )
    gap = []
    for U in tqdm(Us):
        guess = None
        gap_U = []
        for V in Vs:
            H_int = calculate_Hint(U, V, Uk, Vk)
            _gap, guess = compute_gap(
                U=U, V=V, H_int=H_int, nk=nk, tol=tol, mixing=mixing, order=order, guess=guess
            )
            gap_U.append(_gap)
        gap.append(gap_U)
    return np.asarray(gap, dtype=float)
 ```

 %% Cell type:code id:6a8c08a9-7e31-420b-b6b4-709abfb26793 tags:

 ``` python
 gap = compute_phase_diagram(Us, Vs, nk=15, tol=1e-5, mixing=0.01, order=10)
 ```

 %% Output

     22%|██▏       | 11/50 [02:33<07:34, 11.66s/it]/opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=4.86112e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=1.11881e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=6.6988e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=1.04075e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=1.03848e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=1.08015e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    /opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=1.07054e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
     32%|███▏      | 16/50 [03:37<06:44, 11.89s/it]/opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=6.43375e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
     34%|███▍      | 17/50 [03:58<08:05, 14.71s/it]/opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=7.90352e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
     54%|█████▍    | 27/50 [07:18<09:45, 25.48s/it]/opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=3.27211e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
     58%|█████▊    | 29/50 [09:18<15:24, 44.03s/it]/opt/conda/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:1074: LinAlgWarning: Ill-conditioned matrix (rcond=5.90338e-17): result may not be accurate.
      gamma = solve(self.a, df_f)
    100%|██████████| 50/50 [22:03<00:00, 26.47s/it]

 %% Cell type:code id:39edbf19 tags:

 ``` python
 plt.imshow(np.log10(gap).T, origin='lower', extent=(Us.min(), Us.max(), Vs.min(), Vs.max()), vmin=-1)
 plt.colorbar()
 ```

 %% Output

    <matplotlib.colorbar.Colorbar at 0x7f12eff193d0>



 %% Cell type:code id:27f9d0d8 tags:

 ``` python
 plt.imshow((gap).T, origin='lower', extent=(Us.min(), Us.max(), Vs.min(), Vs.max()), vmin=0)
 plt.colorbar()
 ```

 %% Output

    <matplotlib.colorbar.Colorbar at 0x7f12effcd650>



 %% Cell type:code id:e17fc96c-c463-4e1f-8250-c254d761b92a tags:

 ``` python
+import xarray as xr
+gap_da = xr.DataArray(data=gap, coords=dict(Us=Us, Vs=Vs))
+```
+
+%% Cell type:code id:0cb395cd-84d1-49b4-89dd-da7a2d09c8d0 tags:
+
+``` python
+gap_da.to_netcdf('./data/graphene_example.nc')
 ```