DMRG on 2D triangular lattice for J1-J2 spin-1/2 model

[DiracString]

Hi,

I started using TeNPy just a few days ago, and I appreciate all the work being done on it. Please forgive me for the following long, beginner questions, but I wasn't sure where else to go to for help.

1) I want to use DMRG to find the magnetic ground state on a triangular lattice for a spin-1/2 model with nearest-neighbor and next-nearest neighbor interactions (as seen in this Wikipedia page). I wrote the following code, but am not sure whether it is correct (when visualized, the lattice looks reasonable, but I am not sure about what other sanity checks I can perform). Do you have any advice on what I could do? Also, notes on general TeNPy-specific good practices are welcome (since I am still getting used to TeNPy).

Python: Select all

import numpy as np 
import matplotlib.pyplot as plt
from tenpy.models.spins import SpinModel
from tenpy.networks.mps import MPS
from tenpy.algorithms import dmrg
from tenpy.models.lattice import Triangular
from tenpy.networks.site import SpinHalfSite
from tenpy.models.model import CouplingMPOModel
    
bc_xx = 'periodic'; bc_yy = 'periodic'; bc_MPSS = 'infinite'; # boundary conditions
coupling_axis = 'Sz' # adds exchange interaction along Sz axis

lx = 3; ly=3; tot_sites = lx*ly; # lattice
initial_state = ['up'] * (tot_sites) # initial state
j1= 1; j2 = 0.5; # exchange interactions

class MyModel(CouplingMPOModel):
	def init_sites(self, model_params):
		conserve = model_params.get('conserve', None) 
		site = SpinHalfSite(conserve=conserve)
		return site
	def init_lattice(self, model_params):
		Lx = model_params.get('Lx', 3)
		Ly = model_params.get('Ly', 3)
		bc_x = model_params.get('bc_x', bc_xx)
		bc_y = model_params.get('bc_y', bc_yy)
		bc_MPS = model_params.get('bc_MPS', bc_MPSS)
		lattice = Triangular(Lx, Ly, site=self.init_sites(model_params), bc=[bc_x, bc_y], bc_MPS=bc_MPS)
		return lattice
	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)

dmrg_params = {
    'mixer': True,  
    'trunc_params': {
        'chi_max': 10000,
        'svd_min': 1.e-10,
    },
    'max_E_err': 1.e-10,
    'verbose': 1,
}

model_params = {
    'conserve': None,
    'Lx': lx,
    'Ly': ly,
    'J1': j1,
    'J2': j2,
    'bc_MPS': bc_MPSS, 
    'bc_x': bc_xx, 
    'bc_y': bc_yy 
    }
    
M = MyModel(model_params)
lattice = M.init_lattice(model_params)
psi = MPS.from_product_state(lattice.mps_sites(), initial_state, bc=lattice.bc_MPS)

eng = dmrg.TwoSiteDMRGEngine(psi, M, dmrg_params)
info = eng.run() # the main work; modifies psi in place

E = info[0]; psi = info[1];
X = psi.expectation_value("Sx"); Y = psi.expectation_value("Sy"); Z = psi.expectation_value("Sz");

2) I want to test convergence of the ground state with respect to the number of states. However, almost all simulations I run (even up to 30 x 30 size, for both finite, semi-finite and infinite systems), use at most a maximum of chi = 6. Am I right in saying that this is the number of states used (bond dimension)? If so, is there a way for me to fix the number of states used, for me to do this type of convergence test? That is, I want to calculate the ground state energy for each number of states. Although I set chi_max to 10000 above, it doesn't go as anywhere close. I also tried chi_list, but to no avail. Any advice? Or just generally sanity checks I could use to make sure I am starting with a good ground state?

3) Finally, I want to play around and see whether I can use these simulations to check the type of magnetic ground state I have. For example, stripe, 120 AFM, FM, etc. Such as those in the image below (source):

Can anyone provide me a direction on how I could go about using TeNPy to classify simulation systems into those states? I naively thought that I could just calculate the expectation values of the spin operators at each site as I do with X, Y and Z in my code above. And then plot those as vectors (X,Y,Z) to visualize the spin. But, for the range of values I tested, I got only non-zero Z (and only values of 0.5 or -0.5, when I expected at least some cases to be 120 AFM, etc. Am I right in saying that it's more complicated than just getting the spin operator expectation values? If so, I would appreciate any directions I could use to start exploring. This paper uses TeNPy to achieve a similar end result. However, they use a so-called spin structure factor and momentum-space considerations to do so. I don't quite understand it yet (I am thinking I should perform a Fourier transform). But I would like to know what exactly I must do to intuitively classify these ground states.

I apologize for the long post, but I sincerely appreciate your time. Thanks again for all the work that goes into development!

[Johannes]

You need to be aware of symmetries in your system. The SU(2)-symmetric heisenberg interactions in particular preserve the total Sz spin, and the ferromagnetic all-up state (a) in your picture is an exact eigenstate. Hence, when you start from this state, DMRG just doens't update it - that's why your measurments of Z are all 0.5.

To distinguish (a) from (b) and (c), you need to find the actual ground state - depending on the sign and strength of the parameters, this might be the FM eigenstate you start with - but if you have AFM coupling, it's likely rather in the Sz=0 sector. To find the ground state in that sector, start with a state in that sector, e.g. by initializing a Neel state of type (b), and run DMRG from there. You'll see that the bond dimension goes way beyond the 6 you had so far (for the FM, the bond dimension is even just 1 if you don't enable the "mixer"), and that simulation need waaayyy longer to run. 30x30 is totally out of reach to converge, the required bond dimension grows exponentially with Ly, and in this frustated system quite quickly, it's already quite hard to converge the Ly=6 they present in the paper!

For 3), again, you should realize that a single <X> or <Y> is zero if you conserve Sz. Instead, you need to look at correlation functions like <X_i X_j> to distinguish (b) and(c) type states. Indeed, if you Fourier transform those, you get the static structure factor, which will then show peaks at the specific momenta of the states, i.e. at the M point for (b) and at the K/K' points for (c), respectively.