update

example/hubbard.py

@@ -0,0 +1,102 @@
+from renormalizer.mps.backend import np, xp
+from renormalizer.model.op import Op
+from renormalizer.model.basis import BasisHalfSpin
+from renormalizer import Model, Mps, Mpo, optimize_mps
+from renormalizer.utils import EvolveConfig, EvolveMethod
+import logging
+
+logger = logging.getLogger("renormalizer")
+
+"""
+one dimensional Hubbard model with open boundary condition
+H = t \sum_{i=1}^{N-1} a_i^\dagger a_i+1 + a_i+1^\dagger a_i + U \sum_i n_{i\arrowup} n_{i\arrowdown}
+"""
+#------------------------------------------------------------------------
+# Jordan-Wigner transformation maps fermion problem into spin problem
+#
+# |0> => |alpha> and |1> => |beta >:
+#
+#    a_j   => Prod_{l=0}^{j-1}(sigma_z[l]) * sigma_+[j]
+#    a_j^+ => Prod_{l=0}^{j-1}(sigma_z[l]) * sigma_-[j]
+# j starts from 0 as in computer science convention to be consistent
+# with the following code.
+#------------------------------------------------------------------------
+
+nsites = 10
+t = -1
+U = 4
+# the ordering of the spin orbital is
+# 0up, 0down, 1up, 1down,...
+
+# the first number of the two-element list is the change of # of alpha
+# electrons, the second number is for beta electrons
+qn_dict_up = {"+": [-1, 0], "-": [1, 0], "Z": [0, 0]}
+qn_dict_do = {"+": [0, -1], "-": [0, 1], "Z": [0, 0]}
+
+ham_terms = []
+
+for i in range(2*(nsites-1)):
+    if i % 2 == 0:
+        qn1 = [qn_dict_up["Z"], qn_dict_up["+"], qn_dict_do["Z"], qn_dict_up["-"]]
+        qn2 = [qn_dict_up["Z"], qn_dict_up["-"], qn_dict_do["Z"], qn_dict_up["+"]]
+    else:
+        qn1 = [qn_dict_do["Z"], qn_dict_do["+"], qn_dict_up["Z"],
+                qn_dict_do["-"]]
+        qn2 = [qn_dict_do["Z"], qn_dict_do["-"], qn_dict_up["Z"],
+                qn_dict_do["+"]]
+
+    op1 = Op("Z + Z -", [i,i,i+1,i+2], factor=t, qn=qn1)
+    op2 = Op("Z - Z +", [i,i,i+1,i+2], factor=-t, qn=qn2)
+    ham_terms.extend([op1, op2])
+
+for i in range(0,2*nsites,2):
+    qn = [qn_dict_up["-"], qn_dict_up["+"], qn_dict_do["-"], qn_dict_do["+"]]
+    op = Op("- + - +", [i,i,i+1,i+1], factor=U, qn=qn)
+    ham_terms.append(op)
+
+basis = []
+for i in range(2*nsites):
+    if i % 2 == 0:
+        sigmaqn = np.array([[0, 0], [1, 0]])
+    else:
+        sigmaqn = np.array([[0, 0], [0, 1]])
+    basis.append(BasisHalfSpin(i, sigmaqn=sigmaqn))
+
+model = Model(basis, ham_terms)
+mpo = Mpo(model)
+logger.info(f"mpo_bond_dims:{mpo.bond_dims}")
+
+nelec = [5, 5]
+M = 100
+procedure = [[M, 0.4], [M, 0.2], [M, 0.1], [M, 0], [M, 0], [M,0], [M,0]]
+mps = Mps.random(model, nelec, M, percent=1.0)
+logger.info(f"initial mps: {mps}")
+
+# algorithm 1:  DMRG sweep
+mps.optimize_config.procedure = procedure
+mps.optimize_config.method = "2site"
+energies, mps = optimize_mps(mps.copy(), mpo)
+gs_e = min(energies)
+logger.info(f"lowest energy: {gs_e}")
+
+# algorithm 2: imaginary time propagation
+evolve_config = EvolveConfig(EvolveMethod.tdvp_ps,
+        adaptive=True,
+        guess_dt=1e-3/1j,
+        adaptive_rtol=5e-4,
+        ivp_solver="RK45"
+        )
+mps.evolve_config = evolve_config
+evolve_dt = 0.5/1j
+energy_old = 0
+istep = 0
+while True:
+    mps = mps.evolve(mpo, evolve_dt)
+    energy = mps.expectation(mpo)
+    logger.info(f"current mps: {mps}")
+    logger.info(f"istep={istep}, energy={energy}")
+    if np.abs(energy-energy_old) < 1e-5:
+        logger.info("converge!")
+        break
+    istep += 1
+    energy_old = energy

renormalizer/cv/tests/test_H_chain.py

@@ -23,7 +23,6 @@ def test_H_chain_LDOS():
 
     model = Model(basis, ham_terms)
     mpo = Mpo(model)
-
     nelec = [spatial_norbs//2, spatial_norbs//2]
     M = 50
     procedure = [[M, 0.4], [M, 0.2]] + [[M, 0],]*6
@@ -51,7 +50,10 @@ def photoelectron_operator(idx):
     #test_freq = np.linspace(0.25, 1.25, 100, endpoint=False).tolist()
     test_freq = np.linspace(0.25, 1.25, 20, endpoint=False).tolist()
     eta = 0.05
-    M = 10
+    # when M=10, the last site of random mps will sometimes of size (1,2,1)
+    # because of the qn conservation. I slightly enlarge it to make the test
+    # more robust.
+    M = 16
     procedure_cv = [0.4, 0.2] + [0]*6
     spectra = SpectraZtCV(model, None, M, eta, h_mpo=mpo, method="2site",
             procedure_cv=procedure_cv, b_mps=b_mps.scale(-eta), e0=mps_e)

renormalizer/lib/tests/test_krylov.py

@@ -3,6 +3,7 @@
 
 from renormalizer.mps.backend import xp
 from renormalizer.lib import expm_krylov
+from renormalizer.mps.matrix import asxp
 import pytest
 import numpy as np
 from scipy.linalg import expm, eigh
@@ -33,6 +34,6 @@ def test_expm(N, imag, block_size):
     # https://github.com/scipy/scipy/issues/18086
 
     w, x = eigh(a1)
-    res1 = x @ xp.diag(xp.exp(w)) @ x.conj().T @ v
+    res1 = x @ np.diag(np.exp(w)) @ x.conj().T @ v
     res2, _ = expm_krylov(lambda x: a2.dot(x), 1, xp.array(v), block_size)
     assert xp.allclose(res1, res2)

renormalizer/mps/lib.py

@@ -417,13 +417,18 @@ def block_select(basdic, qn, n):
 def compressed_sum(mps_list, batchsize=5, temp_m_trunc=None):
     assert len(mps_list) != 0
     mps_queue = deque(mps_list)
-    while len(mps_queue) != 1:
-        term_to_sum = []
-        for i in range(min(batchsize, len(mps_queue))):
-            term_to_sum.append(mps_queue.popleft())
-        s = _sum(term_to_sum, temp_m_trunc=temp_m_trunc)
-        mps_queue.append(s)
-    return mps_queue[0]
+    if len(mps_queue) > 1:
+        while len(mps_queue) != 1:
+            term_to_sum = []
+            for i in range(min(batchsize, len(mps_queue))):
+                term_to_sum.append(mps_queue.popleft())
+            s = _sum(term_to_sum, temp_m_trunc=temp_m_trunc)
+            mps_queue.append(s)
+        return mps_queue[0]
+    else:
+        new_mps = mps_list[0].canonicalise()
+        new_mps.compress(temp_m_trunc=temp_m_trunc)
+        return new_mps
 
 
 def _sum(mps_list, compress=True, temp_m_trunc=None):