I have been working with a Hubbard Diamond Chain model for some time now. I have done simulations both locally and in a cluster with simulations to find the ground State with iDMRG. I am now trying to calculate excited states but I am unable to do so. First, I put into context the simple calculation of the ground state:
Python: Select all
model_params = {
'L': 2, # Length of the chain
'bc_x':'periodic', # Bound conditions of the lattice
'bc_MPS':'infinite', # Bound conditions of the MPS
'mu':0,
't1': 1.,
'phi': np.pi/4, # Phase factor for the hopping term
't2': .8,
't3': 0.2, # Intercell coupling
'U': 0 # Hubbard parameter }
dmrg_params = {
'mixer': True,
'mixer_params': {
'amplitude': 1.e-5,
'decay': 1.2,
'disable_after': 30 },
'trunc_params': { 'svd_min': 1.e-10, },
'lanczos_params': {
'N_min': 5,
'N_max': 20 },
'chi_list': {
0: 9,
10: 49,
20: 128 },
'max_E_err': 1.e-10,
'max_S_err': 1.e-4,
'max_sweeps': 100,
'max_hours': 0.5
class DiamondChain(Lattice):
dim = 1 #: the dimension of the lattice
Lu = 4 #: the (expected) number of sites in the unit cell, ``len(unit_cell)``.
def __init__(self, L, site, **kwargs):
sites = [site] * 4 # Define the sites in the unit cell
basis = np.eye(1) # La base es la trivial
pos = np.array([[0.2],[0.4],[0.6],[0.8]]) # Equiespaciados en el intervalo [0,1], no pongo nada ni en 0 ni en 1 para no superponerlos con la siguiente celda
kwargs.setdefault('basis', basis) # Guardo en el diccionario kwargs la base
kwargs.setdefault('positions', pos) # Guardo en el diccionario kwargs las posiciones
NN = [(1, 0, np.array([0])), (0, 2, np.array([0])), (3, 1, np.array([0])),
(2, 3, np.array([0]))]
NNN = [(0,3,np.array([0])),(2,1,np.array([0]))]
intercell = [(3,0,np.array([1]))]
kwargs.setdefault('pairs', {}) # Creo pairs para guardar todos los vecinos.
kwargs['pairs'].setdefault('nearest_neighbors', NN)
kwargs['pairs'].setdefault('next_nearest_neighbors', NNN)
kwargs['pairs'].setdefault('intercell', intercell)
Lattice.__init__(self, [L], sites, **kwargs) # Llamo al constructor de la clase Lattice, con
def ordering(self, order):
if isinstance(order, str):
if order == "default":
priority = None
snake_winding = (False, False, False)
return get_order(self.shape, snake_winding, priority)
return super().ordering(order)
class HubbardDiamondChain_Theo(CouplingMPOModel):
default_lattice = DiamondChain
force_default_lattice = True
def init_sites(self, model_params):
cons_N = model_params.get('cons_N', 'N')
cons_Sz = model_params.get('cons_Sz', 'Sz')
site = SpinHalfFermionSite(cons_N=cons_N, cons_Sz=cons_Sz)
return site
def init_terms(self, model_params):
t1 = model_params.get('t1', 1.)
phi = model_params.get('phi',0.)
t1_complex = t1 * np.exp(1j * phi) # Complex hopping term
t1_complex_conj = t1 * np.exp(-1j * phi)
t2 = model_params.get('t2', 1.)
t3 = model_params.get('t3', 1.)
U = model_params.get('U', 0)
mu = model_params.get('mu', 0.)
for u in range(len(self.lat.unit_cell)):
self.add_onsite(mu, u, 'Ntot')
self.add_onsite(U, u, 'NuNd')
# The nearest neighbor terms are now complex, with a phase factor
for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
self.add_coupling(-t1_complex, u1, 'Cdu', u2, 'Cu', dx, plus_hc=True)
self.add_coupling(- t1_complex_conj , u1, 'Cdd', u2, 'Cd', dx, plus_hc=True)
# The other terms remain the same
for u1, u2, dx in self.lat.pairs['next_nearest_neighbors']:
self.add_coupling(-t2, u1, 'Cdu', u2, 'Cu', dx, plus_hc=True)
self.add_coupling(-t2, u1, 'Cdd', u2, 'Cd', dx, plus_hc=True)
for u1, u2, dx in self.lat.pairs['intercell']:
self.add_coupling(-t3, u1, 'Cdu', u2, 'Cu', dx, plus_hc=True)
self.add_coupling(-t3, u1, 'Cdd', u2, 'Cd', dx, plus_hc=True)
Python: Select all
pstate = ["up","up","down","down"] * (model_params['L'])
M = HubbardDiamondChain_Theo(model_params)
psi = MPS.from_product_state(M.lat.mps_sites(), pstate, bc=M.lat.bc_MPS)
eng = dmrg.TwoSiteDMRGEngine(psi, M, dmrg_params)
E0, psi = eng.run()
resume_data = eng.get_resume_data(sequential_simulations=True)
Python: Select all
init_env_data = resume_data['init_env_data'] # Get the initial environment data from the resume data
M_s = model.extract_segment(enlarge=10)
first, last = M_s.lat.segment_first_last
psi0_s = psi2.extract_segment(first, last)
E = dmrg.run(psi0_s, model, dmrg_params,orthogonal_to=psi2)
[issue]
Environments with segment b.c. need explicit environments!
[/issue]
The other thing I've tried is using OrthogonalExcitations :
Python: Select all
init_env_data = resume_data['init_env_data']
options = {'algorithm_params' : dmrg_params,
'model_params': model_params,
'model_class': M,
'segment_enlarge':10,
'init_env_data': init_env_data}
excite = OrthogonalExcitations(options, orthogonal_to=[psi.copy()])
excite.extract_segment_from_infinite(psi.copy(), M, resume_data=resume_data)
excite.run()
[issue]
raise KeyError("i = {0:d} out of bounds for finite MPO".format(i))
return i
KeyError: 'i = 159 out of bounds for finite MPO'
[/issue]