Hamiltonian MPO of different kinds of spins

[underway93]

Hello.

I'm trying to build an MPO for a spin chain of different spins.
For example, system is 1D chain of 4 spin-1 sites and 2 spin-1/2 sites.

Python: Select all

from tenpy.networks.site import SpinHalfSite, SpinSite

N_s1 = 3
N_half = 2
sites_s1 = [SpinSite(1, conserve=None) for _ in range(N_s1)]
sites_half = [SpinHalfSite(conserve=None) for _ in range(N_half)]
sites = sites_s1 + sites_half

And the Hamiltonian is as follows:
\(H_{s1} = \sum_{i \neq j} \vec{S_{i}} \cdot J_{ij} \cdot \vec{S_{j}} + \sum_{i} (D S_{zi}^2 - h S_{zi})\)
\(H_{h} = \sum_{k \neq l} \vec{I_{k}} \cdot J_{kl} \cdot \vec{I_{l}} - h \sum_{k} I_{zk} \)
\(H_{int} = \sum_{i, k} \vec{S_{i}} \cdot A_{ik} \cdot \vec{I_k} \)
\(H_{total} = H_{s1} + H_{h} + H_{int} \)

where \(\vec{S}\) and \(\vec{I}\) are spin operators for the spin-1 and spin-1/2 respectively, \(J\) is the dipole-dipole interaction coefficient, \(A\) is the hyperfine interaction coefficient, and \(h\) is the external magnetic field.

Does anyone know how to build the Hamiltonian and corresponding MPO with TeNPy?

[Johannes]

I would recommend to use the tenpy.models.model.CouplingMPOModel, and override the methods as detailed in the model user guide.
Since you have all-to-all interactions, you don't really have a chain, so you can use a TrivialLattice to write this down - but that implies that you have to manually sum over all the interactions on different sites. If you have two coupled "chains" of equal length and local interactions only, you might consider switching to a Ladder instead, and have implicit sums over nearest neighbors with add_coupling instead of the more explicit add_coupling_term.
(You're using spins here, but for those coming across this thinking about fermions: the add_coupling_term does *not* handle Jordan Wigner strings for you, you need to manually take care of that via the op_string argument!)


Something along the following lines should work:

Python: Select all

class MyModel(CouplingMPOModel):
    def init_sites(self, model_params):
        # TODO: could optionally enable Sz conservation here...
        return (SpinHalfSite(None), SpinSite(S=1., conserve=None))

    def init_lattice(self, model_params):
        N_one = model_params.get("N_one", 3)
        N_half = model_params.get("N_half", 2)
        S_half, S_one = self.init_sites(model_params)
        return tenpy.models.lattice.TrivialLattice([S_one] * N_one + [S_half] * N_half)

    def init_terms(self, model_params):
        # read out parameters
        N_one = model_params["N_one"]
        N_half = model_params["N_half"]
        # i = [0, ... N_one - 1] are spin one sites,
        # i = [N_one, ... N_one + N_half] are spin half sites
        J_one_ij = model_params.get("J_one_ij", np.zeros((N_one, N_one)))
        J_half_ij = model_params.get("J_half_ij", np.zeros((N_half, N_half)))
        D = model_params.get("D", 0.)
        h = model_params.get("D", 0.)
        # H_one terms
        for i in range(N_one):
            for j in range(N_one):
                for op in ['Sx', 'Sy', 'Sz']:
                    self.add_coupling_term(J_ij[i,j], i, j, op, op)
            self.add_onsite_term(D, i, 'Sz Sz')  # op='Sz Sz' is Sz^2
            self.add_onsite_term(-h, i, 'Sz')
        # H_half
        J_half_ij = model_params.get("J_ij", np.zeros((N_half, N_half)))
        for i in range(N_half):
            for j in range(N_half):
                for op in ['Sx', 'Sy', 'Sz']:
                    self.add_coupling_term(J_half_ij[i,j], N_one + i, N_one + j, op, op)
            self.add_onsite_term(-h, N_one + i, 'Sz')
        # H_int
        A_ij = model_params.get("A_ij", np.zeros((N_one, N_half)))
        for i in range(N_one):
            for j in range(N_half):
                for op in ['Sx', 'Sy', 'Sz']:
                    self.add_coupling_term(A_ij[i,j], i, N_one + j, op, op)

[underway93]

Thank you for your guidance.

I successfully have built an MPO. But I want an MPO with open boundary.

Python: Select all

import tenpy
from tenpy.models.model import CouplingMPOModel
from tenpy.networks.site import SpinHalfSite, SpinSite
import numpy as np
from utils import compute_J


class MyModel(CouplingMPOModel):
    def init_sites(self, model_params):
        # TODO: could optionally enable Sz conservation here...
        return (SpinHalfSite(None), SpinSite(S=1.0, conserve=None))

    def init_lattice(self, model_params):
        N_one = model_params.get("N_one", 3)
        N_half = model_params.get("N_half", 2)
        bc_MPS = model_params.get("bc_MPS", "finite")
        S_half, S_one = self.init_sites(model_params)
        lattice = tenpy.models.lattice.TrivialLattice(
            [S_one] * N_one + [S_half] * N_half
        )
        lattice.bc_MPS = "finite"
        self.bc_MPS = "open"
        self.bc_x = "open"
        return lattice

    def init_terms(self, model_params):
        # read out parameters
        S_name_arr = ["Sx", "Sy", "Sz"]
        N_one = model_params["N_one"]
        N_half = model_params["N_half"]
        N = N_one + N_half
        D = model_params.get("D", [0.0] * 3)  # D: x, y, z (1d array of length 3)
        B0 = model_params.get("B0", 0.0)
        gamma_one = model_params["gamma_one"]
        gamma_half = model_params["gamma_half"]
        pos_one = model_params["pos_one"]
        pos_half = model_params["pos_half"]
        pos_one = [np.array(x) for x in pos_one]
        pos_half = [np.array(x) for x in pos_half]
        # H_one terms
        for i in range(N_one - 1):
            for j in range(i + 1, N_one):
                pos1 = pos_one[i]
                pos2 = pos_one[j]
                J_one_ij = compute_J(pos1, pos2, gamma_one, gamma_one)
                for k in range(3):
                    for l in range(3):
                        op1 = S_name_arr[k]
                        op2 = S_name_arr[l]
                        self.add_coupling_term(J_one_ij[k, l], i, j, op1, op2)
            for s in range(3):
                self.add_onsite_term(
                    2.0 * np.pi * D[s], i, " ".join([S_name_arr[s]] * 2)
                )
                # self.add_onsite_term(D, i, "Sz Sz")  # op='Sz Sz' is Sz^2
            self.add_onsite_term(-B0 * gamma_one, i, "Sz")
        for s in range(3):
            self.add_onsite_term(
                2.0 * np.pi * D[s], N_one - 1, " ".join([S_name_arr[s]] * 2)
            )
            # self.add_onsite_term(D, i, "Sz Sz")  # op='Sz Sz' is Sz^2
        self.add_onsite_term(-B0 * gamma_one, N_one - 1, "Sz")

        # H_half
        # J_half_ij = model_params.get("J_ij", np.zeros((N_half, N_half)))
        for i in range(N_half - 1):
            for j in range(i + 1, N_half):
                pos1 = pos_half[i]
                pos2 = pos_half[j]
                J_half_ij = compute_J(pos1, pos2, gamma_half, gamma_half)
                for k in range(3):
                    for l in range(3):
                        op1 = S_name_arr[k]
                        op2 = S_name_arr[l]
                        self.add_coupling_term(
                            J_half_ij[k, l], N_one + i, N_one + j, op1, op2
                        )
            self.add_onsite_term(-B0 * gamma_half, N_one + i, "Sz")
        self.add_onsite_term(-B0 * gamma_half, N - 1, "Sz")

        # H_int
        for i in range(N_one):
            for j in range(N_half):
                pos1 = pos_one[i]
                pos2 = pos_half[j]
                A_ij = compute_J(pos1, pos2, gamma_one, gamma_half)
                for k in range(3):
                    for l in range(3):
                        op1 = S_name_arr[k]
                        op2 = S_name_arr[l]
                        self.add_coupling_term(A_ij[k, l], i, N_one + j, op1, op2)

Running the code above, I get the following output.

Code: Select all

>>> model
<MPO.MyModel object at 0x13f89a880>
>>> model.H_MPO.get_W(0)
<npc.Array shape=(2, 5, 3, 3) labels=['wL', 'wR', 'p', 'p*']>
>>> model.H_MPO.get_W(-1)
<npc.Array shape=(23, 2, 2, 2) labels=['wL', 'wR', 'p', 'p*']>

I want the bond dimension of 'wL' for MPO.get_W(0) is 1 and 'wR' for MPO.get_W(-1) is also 1.

[Simulation Description]
I want to perform real-time evolution of the all-to-all Hamiltonian composed of \(N\) spin-one and \(M\) spin-half system.
We perform the MPS time-evolution code for the 3 spin-one and 2 spin-half system with two-site TDVP engine.
Maximum bond dimension is set to 30.
The result is quite different from the exact time evolution result (Exact Diagonalization).

I guess that this is due to the connection of the first site to the last site in the MPO.

[Johannes]

The fact that you get bond dimension 2 on the boundaries is for technical reasons that we add the "IdL" (ready) and "IdR" (final) states of the MPO finite state machine on each bonds. If you don't want that, overwrite calc_H_MPO in the model by copy-pasting the current implementation, but adding the argument insert_all_id=False in his line

Python: Select all

H_MPO_graph = mpo.MPOGraph.from_terms((ot, ct, edt), self.lat.mps_sites(), self.lat.bc_MPS, insert_all_id=False)

However, having bond dimension 2 at the boundaries does not affect the validty of the MPO or how TDVP is done...

To debug your model, try to read the output of

Python: Select all

print(model.all_coupling_terms().to_TermList())
print(model.all_onsite_terms().to_TermList())

and see whether you got the hamiltonian terms as you expect (with the correct prefactors).

Only run time evolution once you're confident that H is right. For a small system with 3+2 sites, you can also use tenpy.algorithms.exact_diag.ExactDiag to see whether your model implementation agrees with whatever other ED you were running.
Once that works, try to find the sources of error for the TDVP - did you choose small enough time steps and large enough bond dimensions? Bond dimension 30 is small (typically simulations use several 100s), but for the 3+2 site sytem exact anyways, since the Hilbert space is so small in that case.