Magnetization from dmrg

[HamidArianZad]

What is your issue?

Dear quimb team,

I am interested to compute the magnetization of 1D and 2D models by implementing dmrg method. How can I do that? I searched throughout the quimb website but I could not find any direct resource to obtain the local magnetization and total magnetization of a desired 1D or 2D model by using dmrg method.

thanks for your support.

[jcmgray]

Hi @HamidArianZad. MPSs do have the method magnetization: https://quimb.readthedocs.io/en/latest/_autosummary/quimb.tensor.tensor_1d.html?highlight=magnetization#quimb.tensor.tensor_1d.MatrixProductState.magnetization.

Maybe this tutorial with usage demonstrated in the docs would be helpful too: https://quimb.readthedocs.io/en/latest/examples/ex_TEBD_evo.html?highlight=magnetization#Non-translationally-invariant-Hamiltonians.

For 2D, you would have to specify the local operators yourself and use compute_local_expectation: https://quimb.readthedocs.io/en/latest/_autosummary/quimb.tensor.tensor_2d.html?highlight=compute_local_expectation#quimb.tensor.tensor_2d.TensorNetwork2DVector.compute_local_expectation

Might these cover your use cases?

[HamidArianZad]

I tried following code to reproduce the Hamiltonian of a 2D square lattice of L = 4*5. Then I applied DMRG2 method to gain the ground-state energy and magnetization of the model.

---

def ham_heis_2D(m, n, j=1.0, bz=0.0, cyclic=False, sparse=True):

      dims = [[2] * m] * n  # shape (n, m)
    # generate tuple of all site coordinates
    sites = tuple(itertools.product(range(n), range(m)))

    # generate neighbouring pairs of coordinates
    def gen_pairs():
        for i, j, in sites:
            above, right = (i + 1) % n, (j + 1) % m
            # ignore wraparound coordinates if not cyclic
            if cyclic or above != 0:
                yield ((i, j), (above, j))
            if cyclic or right != 0:
                yield ((i, j), (i, right))

    # generate all pairs of coordinates and directions
    pairs_ss = tuple(itertools.product(gen_pairs(), 'xyz'))

    # make XX, YY and ZZ interaction from pair_s
    #     e.g. arg ([(3, 4), (3, 5)], 'z')
    def interactions(pair_s):
        pair, s = pair_s
        Sxyz = spin_operator(s, sparse=True)
        return ikron([j * Sxyz, Sxyz], dims, inds=pair)

    # function to make Z field at ``site``
    def fields(site):
        Sz = spin_operator('z', sparse=True)
        return ikron(bz * Sz, dims, inds=[site])

    # combine all terms
    all_terms = itertools.chain(map(interactions, pairs_ss),
                                map(fields, sites) if bz != 0.0 else ())
    H = sum(all_terms)

    # can improve speed of e.g. eigensolving if known to be real
    if isreal(H):
        H = H.real

    if not sparse:
        H = qarray(H.A)

    return H

n, m = 4, 5
dims = [[2] * m] * n

L = m * n

H = ham_heis_2D(m, n, cyclic=True)

dmrg = DMRG2(H, L)

dmrg.solve(max_sweeps=2, verbosity=1, cutoffs=1e-6)

---

But I get following error:

---

Traceback (most recent call last): File "/home/hamid/PycharmProjects/Hamid Python Project/quimb Project/Mag_DMRG_2D.py", line 64, in dmrg = DMRG2(H, L) File "/home/hamid/anaconda3/lib/python3.9/site-packages/quimb/tensor/tensor_dmrg.py", line 1080, in init super().init(ham, bond_dims=bond_dims, cutoffs=cutoffs, File "/home/hamid/anaconda3/lib/python3.9/site-packages/quimb/tensor/tensor_dmrg.py", line 541, in init self.L = ham.L File "/home/hamid/anaconda3/lib/python3.9/site-packages/scipy/sparse/_base.py", line 771, in getattr raise AttributeError(attr + " not found") AttributeError: L not found

---

Could you please guide me to resolve this problem?

[jcmgray]

Hi @HamidArianZad, the problem is that the DMRG algorithm requires its Hamiltonian in MPO form, see e.g. the docstring here. quimb does not yet have the functionality to generate MPO forms of 2D/arbitrary hamiltonians, but it will at some point soon (or you could try constructing the MPO yourself).

For now you would have to use a native 2D / PEPS algorithm, which should scale better for non-'thin' systems too.

[HamidArianZad]

Hello,

Thank you very much!

I tried to run below tutorial:

import quimb as qu
import quimb.tensor as qtn
from matplotlib import pyplot as plt

Lx = 4
Ly = 4
D = 2
max_bond = 64
chi = 64
ctm = 100
seed = 999

peps = qtn.PEPS.rand(Lx=Lx, Ly=Ly, bond_dim=D, seed=seed)

norm = peps.H & peps

norm.contract(all, optimize='auto-hq')

norm.contract_boundary(max_bond=max_bond)

peps.add_tag('KET')
pepsH = peps.conj().retag({'KET': 'BRA'})
norm = pepsH & peps

norm.contract_boundary_(max_bond=max_bond, layer_tags=['KET', 'BRA'], around=((2, 2), (2, 3)))

H2 = qu.ham_heis(2)

coo_a = (1, 1)
coo_b = (1, 2)

peps.compute_local_expectation(
    {(coo_a, coo_b): H2},
    max_bond=max_bond,
    normalized=True,
)

terms = {
    (coo_a, coo_b): H2
    for coo_a, coo_b in peps.gen_bond_coos()
}

peps.compute_local_expectation(
    terms,
    max_bond=max_bond,
    normalized=True,
)

ham = qtn.LocalHam2D(Lx, Ly, H2=H2)

H2 = {None: qu.ham_heis(2)}

H1 = {}

for i in range(Lx):
    for j in range(Ly):

        # add next nearest neighbor interactions
        if (i + 1 < Lx) and (j - 1 >= 0):
            H2[(i, j), (i + 1, j - 1)] = 0.5 * qu.ham_heis(2)
        if (i + 1 < Lx) and (j + 1 < Ly):
            H2[(i, j), (i + 1, j + 1)] = 0.5 * qu.ham_heis(2)

        # add a random field
        H1[i, j] = qu.randn() * qu.spin_operator('Z')

ham_nn_r = qtn.LocalHam2D(Lx, Ly, H2=H2, H1=H1)

ham = qtn.LocalHam2D(Lx, Ly, H2=qu.ham_heis(2))

energy_exact = qu.groundenergy(qu.ham_heis_2D(Lx, Ly, sparse=True)) / (Lx * Ly)

print("energy_exact = ", energy_exact, "\n")

psi0 = qtn.PEPS.rand(Lx, Ly, bond_dim=D, seed=seed, dtype=float)

su = qtn.SimpleUpdate(
    psi0,
    ham,
    chi=chi,  # boundary contraction bond dim for computing energy
    compute_energy_every=20,
    compute_energy_per_site=True,
    keep_best=True,
)

for tau in [0.1, 0.01]:
    su.evolve(ctm, tau=tau)

print("\n su.best: ", su.best)
print("su.energies: ", su.energies)

psi0 = su.best['state'].copy()

# psi0 = qtn.PEPS.rand(Lx, Ly, bond_dim=4)

fu = qtn.FullUpdate(
    psi0=psi0,
    ham=ham,
    # chi again is the boundary contraction max_bond
    # now used for the envs as well as any energy calc
    chi=chi,
    # we thus can cheaply compute the energy at every step
    compute_energy_every=1,
    compute_energy_per_site=True,
    keep_best=True,
)

fu.evolve(50, tau=0.02)

fu.fit_opts['als_enforce_pos'] = True

def loss(psi, terms):
    # the following functions simply scale the various tensors
    #     for the sake of numerical stability
    psi.balance_bonds_()
    psi.equalize_norms_(1.0)

    # then we just compute the energy of all the terms
    return psi.compute_local_expectation(
        terms,
        max_bond=max_bond,
        cutoff=0.0,
        normalized=True
    ) / (Lx * Ly)

tnopt = qtn.TNOptimizer(
    # initial TN to optimize
    psi0,
    # the function to minimize
    loss_fn=loss,
    # constant TNs, tensors, arrays
    loss_constants={'terms': ham.terms},
    # the library that computes the gradient
    autodiff_backend='jax',
    # the scipy optimizer that makes use of the gradient
    optimizer='L-BFGS-B',
)

psi_opt = tnopt.optimize(seed)

print("psi_opt = ", psi_opt)

But I get below warning:

energy_exact = -0.5743254415745596

n=100, tau=0.1, energy~-0.544667: 100%|##########| 100/100 [00:01<00:00, 64.62it/s] n=200, tau=0.01, energy~-0.544182: 100%|##########| 100/100 [00:01<00:00, 70.41it/s]

su.best: {'energy': -0.5446668257748429, 'state': <PEPS(tensors=16, indices=40, Lx=4, Ly=4, max_bond=2)>, 'it': 100} su.energies: [-0.01394248455141442, -0.42462802975791175, -0.5438133939333119, -0.5445453694088561, -0.5445222409923282, -0.5446668257748429, -0.5444968907069688, -0.5443706273752631, -0.5442786430312206, -0.5442215405191778, -0.5441819975379001] n=50, tau=0.02, energy~-0.544600: 100%|##########| 50/50 [00:08<00:00, 5.76it/s]

WARNING:jax._src.lib.xla_bridge:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)

0%| | 0/999 [00:00<?, ?it/s]

---

I tried several ways to remove this warning and I spent much time to optimize the results of the SU and FU methods with jax under gpu, but I failed.

I installed below dependencies for quimb and gpu usages:

1. Anaconda 3 + virtual environment + python=3.9 + quimb 1.4.0
2. jax=0.3.25 + jaxlib=0.3.25
3. Cudatoolkit = 11.2.2 + cuda-nvcc = 11.8.89 (+ cudatoolkit dev 11.3.1) + cupy=9.1.0

My laptop has below gpu config: Nvidia GeForce RTX 3060 6GB GDDR6 & Integrated graphics card: Intel Iris Xe Graphics HD

+-----------------------------------------------------------------------------+
| NVIDIA-SMI 515.65.01    Driver Version: 515.65.01    CUDA Version: 11.7     |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|                               |                      |               MIG M. |
|===============================+======================+======================|
|   0  NVIDIA GeForce ...  Off  | 00000000:01:00.0 Off |                  N/A |
| N/A   42C    P0    N/A /  N/A |      5MiB /  6144MiB |      0%      Default |
|                               |                      |                  N/A |
+-------------------------------+----------------------+----------------------+

+-----------------------------------------------------------------------------+
| Processes:                                                                  |
|  GPU   GI   CI        PID   Type   Process name                  GPU Memory |
|        ID   ID                                                   Usage      |
|=============================================================================|
|    0   N/A  N/A      1848      G   /usr/lib/xorg/Xorg                  4MiB |
+-----------------------------------------------------------------------------+

Could you please also address this problem?

[jcmgray]

I'm afraid that seems to be jax installation issue unrelated to quimb, so you might have to find support elsewhere.

My only suggestion is that you appear to have at least 4 different cuda versions installed (11.2, 11.8, 11.3, 11.7), which are probably conflicting. I'd setup a fresh conda environment and install jax in that via the recommended pip install --upgrade "jax[cuda]" -f https://storage.googleapis.com/jax-releases/jax_cuda_releases.html.