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)