I am working on a research project involving the 1D XXZ model, and I'm facing a persistent challenge when trying to run DMRG simulations with charge conservation. My primary goal is to calculate the <SxSx> correlation function in the critical phase (Δ <= 1) to extract a critical exponent.
For performance reasons, especially on larger chains (N > 60), it is essential for me to use the conserve='Sz' U(1) symmetry. However, this creates a fundamental conflict: guidance.
The Problem I'm Trying to Solve:
I need to calculate the ground state of the XXZ model efficiently and then compute correlation functions. The core conflict is:
To make the simulation computationally feasible for large systems, I must use conserve='Sz'.
To test my hypothesis, I must be able to measure the off-diagonal correlation function <SxSx>.
What I Have Tried So Far (The Debugging Journey):
My attempts to reconcile these two requirements have led me down a path of errors, which I believe point to a fundamental misunderstanding of the TeNPy API on my part.
- conserve=None: This approach works perfectly for measuring all observables, but it is computationally inviable. A single data point for N=40 takes hours on a powerful machine, making it unsuitable for a full campaign.
- conserve='Sz' (The Fast Way): This is incredibly fast, as expected. However, it immediately leads to a ValueError: ... doesn't have the operator 'Sx', because the default SpinHalfSite with U(1) symmetry does not define Sx.
My first attempt was to add Sx after model creation using site.add_op('Sx', 0.5 * (site.Sp + site.Sm)). This correctly failed with a ValueError: Arrays can't have different qtotal!. I now understand this is a physically illegal operation, as you cannot add operators with different quantum numbers.
My next attempts involved trying to access the site object from the model.lat attribute to modify it. This led to a series of AttributeErrors, such as 'Chain' object has no attribute 'sites' and, most recently, AttributeError: 'function' object has no attribute 'opnames' when trying to access the simulation becomes fast, but the default SpinHalfSite used by XXZChain no longer has the Sx operator, making it impossible to measure the very observable I'm interested in.
I have spent a significant amount of time trying to solve this, and I seem to be going in circles. I would be incredibly grateful for some guidance on the canonical way to solve this in TeNPy.
Summary of My Failed Attempts:
conserve=None: This approach is functionally correct. All operators are available, and I can measure everything. However, it is computationally inviable. Simulations for N=40 hang indefinitely on a powerful machine, as the Hilbert space is too large. This path is a dead end.
conserve='Sz' and Post-Creation add_op: This seemed like the most promising direction. My strategy was to create the XXZChain model (which defaults to conserve='Sz') and then manually add the Sx operator to the site afterward. This led to the following critical error:
ValueError: Arrays can't have different qtotal!
This error occurs on the line:
site.add_op('Sx', 0.5 * (site.Sp + site.Sm))
I now understand why this fails: Sp and Sm aremodel.lat.site.opnames. This shows I am fundamentally misunderstanding how to access and modify the site object within a pre-built model like XXZChain.
My Specific Question:
What is the canonical, correct pattern in TeNPy to construct a model for the XXZ chain that:
a) Uses conserve='Sz' for an efficient DMRG calculation.
b) Still allows for the measurement of off-diagonal correlation functions like <SxSx> on the resulting ground state MPS?
I feel I am missing one key piece of knowledge about the TeNPy API for handling this very common scenario.
Here is my latest (non-working) script:
This script uses the conserve='Sz' approach and fails with the AttributeError: 'function' object has no attribute 'opnamesnpcarrays with opposite charges (+1and-1), and they cannot be legally added together. My attempts to buildSx` from them were conceptually flawed.
Other API issues: My journey to this point has been plagued by various AttributeErrors (.sites vs .site) and TypeErrors ('ops' or 'add_ops' as keyword arguments), which I now believe are symptoms of me trying to force an API syntax that is incorrect for my TeNPy version.
My Specific Question:
What is the canonical, correct way in TeNPy (v1.0.6) to set up a XXZChain model that:
a) Uses conserve='Sz' for an efficient DMRG calculation.
b) Still allows for the subsequent'` because my attempt to modify the site is incorrect.
Python: Select all
import numpy as np
import pandas as pd
import time
import os
import json
from datetime import datetime
from tenpy.networks.mps import MPS
from tenpy.models.xxz_chain import XXZChain
from tenpy.algorithms import dmrg
import warnings
from scipy.optimize import curve_fit
from scipy.stats import linregress
# Imports necesarios para la solución correcta
from tenpy.networks.site import SpinHalfSite
from tenpy.models.lattice import Chain
# --- CONFIGURACIÓN FINAL (BASADA EN EL ÉXITO DE LA PRUEBA DE FUEGO) ---
CONFIG = {
"system_sizes": [20, 40, 60, 80],
"Delta_ranges": [
{"min": 0.5, "max": 1.5, "points": 21},
],
"dmrg_params": {
'chi_list': {0: 100, 5: 200, 15: 400},
'max_sweeps': 100,
'mixer': False,
# CORRECCIÓN DE SINTAXIS: Los parámetros de convergencia van al nivel principal
'max_E_err': 1.e-8,
'max_S_err': 1.e-7,
'verbose': 0
},
"correlation_params": {
"min_distance": 2,
"max_distance_ratio": 0.45,
},
"output_dir": "resultados_xxz_finalisimo",
"fitting_params": {
"min_r2_threshold": 0.90,
}
}
# --- FUNCIÓN DE CREACIÓN DEL MODELO (CON LA CORRECCIÓN DEFINITIVA) ---
def create_xxz_model(N, Delta):
"""
Crea el modelo XXZ usando la conservación de carga y asegurando
que los operadores Sx/Sy estén disponibles.
"""
model_params = {'L': N, 'Jz': Delta, 'conserve': 'Sz'}
model = XXZChain(model_params)
# LA CORRECCIÓN CLAVE: Acceder al sitio con .site (singular)
site = model.lat.site
if 'Sx' not in site.opnames:
site.add_op('Sx', 0.5 * (site.Sp + site.Sm))
site.add_op('Sy', -0.5j * (site.Sp - site.Sm))
return model
# --- FUNCIÓN DE AJUSTE "INTELIGENTE" ---
def fit_correlation_decay(distances, correlations, config, Delta):
op_to_fit = 'Sx' if Delta <= 1.0 else 'Sz'
valid_data = [(d, c) for d, c in zip(distances, correlations) if c is not None and c > 1e-12]
if len(valid_data) < 4: return {'fit_type': 'insufficient_data'}
r_data, corr_data = np.array([d for d,c in valid_data]), np.array([c for d,c in valid_data])
if op_to_fit == 'Sx':
try:
log_r, log_corr = np.log(r_data), np.log(corr_data)
slope, intercept, r_value, _, _ = linregress(log_r, log_corr)
if r_value**2 > config["fitting_params"]["min_r2_threshold"]:
return {'fit_type': 'power_law', 'alpha': -slope, 'R2': r_value**2}
except Exception: pass
else:
try:
def exponential_decay(r, A, xi): return np.abs(A) * np.exp(-r / xi)
popt, _ = curve_fit(exponential_decay, r_data, corr_data, p0=[1.0, 2.0])
r_squared = 1 - (np.sum((corr_data - exponential_decay(r_data, *popt))**2) / np.sum((corr_data - np.mean(corr_data))**2))
if r_squared > config["fitting_params"]["min_r2_threshold"]:
return {'fit_type': 'exponential', 'xi': popt[1], 'R2': r_squared}
except Exception: pass
return {'fit_type': 'no_good_fit'}
# --- FUNCIÓN DE SIMULACIÓN (CON LA LÓGICA CORRECTA) ---
def run_single_point(N, Delta, config):
start_time = time.time()
base_result = {'N': N, 'Delta': Delta, 'converged': False, 'fit_type': 'not_attempted', 'alpha': 'NO_FIT', 'xi': 'NO_FIT', 'R2': 'NO_FIT', 'operator_used': ''}
try:
model = create_xxz_model(N, Delta)
initial_state = ['up' if i % 2 == 0 else 'down' for i in range(N)]
psi = MPS.from_product_state(model.lat.mps_sites(), initial_state, bc='finite')
info = dmrg.run(psi, model, config["dmrg_params"])
converged = info.get('converged', False)
base_result.update({'converged': converged})
if converged:
op_to_fit = 'Sx' if Delta <= 1.0 else 'Sz'
ref_site = N // 2
max_dist = int(N * config["correlation_params"]["max_distance_ratio"])
sites2 = range(ref_site + config["correlation_params"]["min_distance"], ref_site + max_dist)
distances = np.array(sites2) - ref_site
corr_values = [abs(psi.correlation_function(op_to_fit, op_to_fit, sites1=[ref_site], sites2=[s2])[0]) for s2 in sites2]
fit_result = fit_correlation_decay(distances, corr_values, config, Delta)
base_result.update(fit_result)
base_result['operator_used'] = op_to_fit
except Exception as e:
print(f"\nERROR CRÍTICO (N={N}, Delta={Delta}): {e}\n")
import traceback
traceback.print_exc()
base_result['computation_time'] = time.time() - start_time
return base_result
# --- FUNCIÓN DE CAMPAÑA ---
def run_campaign(config):
print("=== CAMPAÑA DMRG XXZ (VERSIÓN FINAL CON CONSERVACIÓN) ===")
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
output_dir = f"{config['output_dir']}_{timestamp}"
os.makedirs(output_dir, exist_ok=True)
with open(os.path.join(output_dir, "config.json"), 'w') as f: json.dump(config, f, indent=2)
all_results = []
total_start_time = time.time()
for N in config["system_sizes"]:
print(f"\n=== Sistema N={N} ===")
all_deltas = sorted(list(set(np.concatenate([np.linspace(r["min"], r["max"], r["points"]) for r in config["Delta_ranges"]]))))
for i, Delta in enumerate(all_deltas):
print(f" [{i+1:2d}/{len(all_deltas)}] N={N}, Delta={Delta:.3f} -> Calculando...", end='', flush=True)
result = run_single_point(N, Delta, config)
all_results.append(result)
fit_type = result.get('fit_type', 'N/A'); op = result.get('operator_used', '')
status = '✓' if result.get('converged', False) else '✗'
comp_time = result.get('computation_time', 0.0)
alpha_val = result.get('alpha'); xi_val = result.get('xi'); r2_val = result.get('R2')
alpha_str = f"{alpha_val:.4f}" if isinstance(alpha_val, float) else str(alpha_val)
xi_str = f"{xi_val:.4f}" if isinstance(xi_val, float) else str(xi_val)
r2_str = f"{r2_val:.3f}" if isinstance(r2_val, float) else str(r2_val)
if fit_type == 'power_law': summary = f"fit=pow({op}), α={alpha_str}"
elif fit_type == 'exponential': summary = f"fit=exp({op}), ξ={xi_str}"
else: summary = f"fit={fit_type}"
print(f"\r [{i+1:2d}/{len(all_deltas)}] Delta={Delta:.3f} → {summary}, R²={r2_str} {status} ({comp_time:.1f}s)")
df = pd.DataFrame(all_results)
results_file = os.path.join(output_dir, "results_combined.csv")
df.to_csv(results_file, index=False)
print(f"\nResultados guardados en: {results_file}")
total_time = time.time() - total_start_time
print(f"\n=== CAMPAÑA COMPLETADA ===\nTiempo total: {total_time/60:.1f} minutos\nResultados en: {output_dir}")
# --- MAIN ---
if __name__ == "__main__":
warnings.filterwarnings("ignore", category=UserWarning)
run_campaign(CONFIG)
My Environment:
TeNPy Version: 1.0.6 (physics-tenpy)
Python Version: 3.9
OS: Ubuntu (running on AWS EC2)
Thank you for your time and expertise. Any guidance on the correct way to construct this model would be immensely appreciated.