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

Code: Select all

CouplingMPOModel

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

Code: 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.

Code: 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.

Code: 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

Code: 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:

Code: 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)
            self.add_coupling(-Dy_p, u1, 'Sz', u2, 'Sx', dx)
            # 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).