How to model anisotropic exchange J, single-ion anisotropy and DMI for DMRG?

[DiracString]

Hi,

Currently, I am able to run DMRG calculations for the

Python: Select all

CouplingMPOModel

with isotropic exchange J parameters J1 and J2 for NN, and NNN using:

Python: Select all

def init_terms(self, model_params):
		J1 = model_params.get('J1', 1.0)
		J2 = model_params.get('J2', 0.0)
		for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
			self.add_coupling(J1, u1, coupling_axis, u2, coupling_axis, dx)
		for u1, u2, dx in self.lat.pairs['next_nearest_neighbors']:
			self.add_coupling(J2, u1, coupling_axis, u2, coupling_axis, dx)

However, I now want to include anisotropic terms like Jxx, Jyy, Jzz, Jxz, Jzx, Jyz, Jzx, ...(including antisymmetric); Dx, Dy, Dz; Axy, Axz, Azz-Axx, etc for anisotropic exchange J, DMI D and single-ion anisotropy A. I want to do this for NN, NNN and NNNN couplings; using the spin-1/2 Heisenberg model for a 2D triangular lattice.

Would someone mind clarifying how I could go about setting this up in TeNPy?
Apologies if this is a newbie question, I am still getting familiar with this.

[Johannes]

Not sure what your question is - you already have a loop over the coupling_axis, that you didn't include in your sample code.
Just read out different model_params and explicitly use those instead of the plain J1/J2 for the various terms, e.g.

Python: Select all

def init_terms(self, model_params):
    Jxx = model_params.get("Jxx", 0.)
    for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
        self.add_coupling(Jxx, u1, 'Sx', u2, 'Sx', dx)

Note that you need to be a little bit careful about symmetries, e.g. if you want to use Sz conservation. In that case, you need to 1) write it in terms of Sp and Sm instead of Sx, Sy, and 2) enable that conservation in the init_sites()
See the source code of tenpy.models.spins.SpinModel for a more detailed example; in fact that model already supports the Jx, Jy, Jz terms.
So you just need to add more terms for D, A and further range couplings.

[DiracString]

@Johannes, thank you! That source code is very helpful. I think I understand how to deal with anisotropic J and single-ion anisotropy. However, I am lost regarding incorporating DMI. I have the DMI vector D = (Dx,Dy,Dz), and this vector is different depending on the spin site chosen. I have a triangular lattice with P3m1's symmetry, and D is 120degree rotated between each consecutive nearest-neighbor pair. Do you have any advice on how I could go about incorporating this? The corresponding Hamiltonian term is like: D . (Si x Sj), where we have the cross product.
Thank you for your patience!

EDIT: Would it simply be something like writing out the cross product explicitly? [Dx (SySz-SzSy), -Dy(SxSz-SzSx), Dz(SxSy-SySx)]? In this case, I am not sure how to incorporate the direction (I am guessing it has to do with the 'dx' term below... but I am not too sure how to exactly do this for my triangular lattice.

Python: Select all

    Dx = model_params.get("Dx", 0.)
    Dy = model_params.get("Dy", 0.)
    for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
        self.add_coupling(Dx, u1, 'Sy', u2, 'Sz', dx)
        self.add_coupling(Dx, u1, 'Sz', u2, 'Sy', dx)
        
        self.add_coupling(-Dy, u1, 'Sx', u2, 'Sz', dx)
        self.add_coupling(-Dy, u1, 'Sz', u2, 'Sx', dx)

[Johannes]

Yes, for the DM interaction, you need to write out the cross product explicitly and add the Dx/Dy/Dz terms separately. Note that you got the signs in the cross product wrong, c.f. the code below.

Again, I'm not sure if I understood correctly - it seems you want the D to depend on which pairs of NN you couple, and the D should be e.g. parallel or perpendicular to that direction.
To acchive that, you should reflect this in your code: instead of just looping over the nearest-neighbor pairs without adjusting the D, specialize your model to the triangular model, and hard-code that the D depends on which hopping you have.
To specialize to the Triangular model, you can simply set the default_lattice and force_default_lattice. Then you know that they NN pairs are as defined in the source code of the tenpy.models.lattice.Triangular, i.e. as

Python: Select all

        NN = [(0, 0, np.array([1, 0])), (0, 0, np.array([-1, 1])), (0, 0, np.array([0, -1]))]

Note that those couplings have a direction going counter-clockwise, with angles of 30° to the x axis, then 150°, and finally 270°.
Your code could hence look like this:

Python: Select all

class MyModel(CouplingMPOModel):
    default_lattice = 'Triangular'
    force_default_lattice = True
    ... # def init_sites()
    def init_terms(self, model_params):
        Dx = model_params.get("Dx", 0.)
        Dy = model_params.get("Dy", 0.)
        for (u1, u2, dx), theta in zip(self.lat.pairs['nearest_neighbors'], 
                                       np.array([0., 1., 2.])*2.*np.pi/3. + np.pi/6.):
            rot_U = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
            Dx_p, Dy_p = rot_U @ np.array([Dx, Dy])
            self.add_coupling(Dx_p, u1, 'Sy', u2, 'Sz', dx)
            self.add_coupling(-Dx_p, u1, 'Sz', u2, 'Sy', dx)
            self.add_coupling(Dy_p, u1, 'Sx', u2, 'Sz', dx)   # EDIT: wrong sign here, see correct version below
            self.add_coupling(-Dy_p, u1, 'Sz', u2, 'Sx', dx)  # EDIT: wrong sign here, see correct version below
            # similar for Dz, if also needed

If you wanted to be very fancy/general, you could even get the angle of the coupling direction numerically from the (u1, u2, dx) and the lattice unit_cell and basis vectors, with a code very similar to distance, changing only the last line of that function to calculate the angle instead of the distance of the given vector(s).

[DiracString]

Thank you for your response @Johannes ! Would you mind clarifying the following?

1) Why do you add np.pi/6 to the array of angles?

2) If I specify anisotropic exchange \(J.S_i.S_j\) (Jxx,Jxy,Jxz,Jyx,Jyy,Jyz,Jzx,Jzy,Jzz) via add_coupling(), and my system has DMI too, will I need to specify DMI D separately (using \(D. S_i \times S_j\)) ? Or will the DMI be automatically treated from the anisotropic J values in TenPy, since they are related by: Dz = (Jxy - Jyx)/2, etc (since components of the DMI vector are related to the antisymmetric part of J)?

3) Regarding the incorrect cross product sign, it was my understanding that if I have D=(Dx,Dy,Dz), S1=(s1x,s1y,s1z) and S2=(s2x,s2y,s2z), then D.(S1 x S2) would be: Dx (s1y s2z - s1z s2y) + Dy (s1z s2x - s1x s2z) + Dz (s1x s2y - s1y s2x). Wouldn't this reverse the signs of the Dy_p in your code? Apologies if I got something wrong here.

Thank you very much for your time and patience!

[Johannes]

1) If you look at the imange of the tenpy.models.lattice.Triangular lattice, you'll see that it's 30 degrees rotatet compared to your image - that's the pi/6 offsett here.

2) If you explicitly have Jxy, Jyx etc, then these will contribute to the very same "terms" Sx_i Sy_j etc that the DMI interactions contribute to, yes. The model will simply add the prefactors for all of the terms you added with `add_coupling`, if multiple of them had the same operators. If they add up to 0, they will not explicitly appear in the MPO.

3) You are of course right, there's a sign error in the cross product I wrote, it should read

Python: Select all

class MyModel(CouplingMPOModel):
    default_lattice = 'Triangular'
    force_default_lattice = True
    ... # def init_sites()
    def init_terms(self, model_params):
        Dx = model_params.get("Dx", 0.)
        Dy = model_params.get("Dy", 0.)
        for (u1, u2, dx), theta in zip(self.lat.pairs['nearest_neighbors'], 
                                       np.array([0., 1., 2.])*2.*np.pi/3. + np.pi/6.):
            rot_U = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
            Dx_p, Dy_p = rot_U @ np.array([Dx, Dy])
            self.add_coupling(Dx_p, u1, 'Sy', u2, 'Sz', dx)
            self.add_coupling(-Dx_p, u1, 'Sz', u2, 'Sy', dx)
            self.add_coupling(Dy_p, u1, 'Sz', u2, 'Sx', dx)  # corrected order!!!
            self.add_coupling(-Dy_p, u1, 'Sx', u2, 'Sz', dx) # corrected order!!!
            # similar for Dz, if also needed