TeNPy Forum

Python: Select all

# --- Imports ---
import numpy as np
import matplotlib.pyplot as plt
from tenpy.models.model import CouplingMPOModel
from tenpy.networks.site import SpinHalfSite
from tenpy.models.lattice import Chain
from tenpy.networks.mps import MPS
from tenpy.algorithms import dmrg

# --- Custom J1–J2 Heisenberg Model ---
class J1J2Model(CouplingMPOModel):
    def init_sites(self, model_params):
        conserve = model_params.get('conserve', None)
        return SpinHalfSite(conserve=conserve)

    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, 'Sz', u2, 'Sz', dx)
            self.add_coupling(0.5 * J1, u1, 'Sp', u2, 'Sm', dx, plus_hc=True)

        for u1, u2, dx in self.lat.pairs['next_nearest_neighbors']:
            self.add_coupling(J2, u1, 'Sz', u2, 'Sz', dx)
            self.add_coupling(0.5 * J2, u1, 'Sp', u2, 'Sm', dx, plus_hc=True)

# --- Parameters ---
L = 20
J1 = 1.0
J2 = 0.0

# --- Lattice with PBC ---
lat = Chain(L, site=SpinHalfSite(conserve=None), bc='periodic')

# --- Model ---
model_params = {'lattice': lat, 'J1': J1, 'J2': J2}
model = J1J2Model(model_params)

# --- Initial MPS ---
init_state = ['up', 'down'] * (L // 2)
psi = MPS.from_product_state(model.lat.mps_sites(), init_state, bc='finite')

# --- Run DMRG ---
dmrg_params = {
    'mixer': True,
    'max_E_err': 1e-10,
    'trunc_params': {'chi_max': 800, 'svd_min': 1e-10},
    'max_sweeps': 20,
}
info = dmrg.run(psi, model, dmrg_params)

# --- Print and rescale energy ---
E_tp = info['E']
E_tp_scaled = 4.0*E_tp   # match NetKet convention
print(f"Ground state energy (TeNPy):        {E_tp:.12f}")
print(f"Rescaled energy to match NetKet ED: {E_tp_scaled:.12f}")


# --- Site-local spin operators ---
Sx = [model.lat.mps_sites()[i].Sx for i in range(L)]
Sy = [model.lat.mps_sites()[i].Sy for i in range(L)]
Sz = [model.lat.mps_sites()[i].Sz for i in range(L)]




# --- Translation-averaged spin-spin correlator ---
def avg_corr(psi, op_list, L, r):
    vals = []
    for R in range(L):
        i = R
        j = (R + r) % L  # periodic wrap-around
        # Get site index directly from the lattice
        val = psi.expectation_value_term([[i, op_list[i].site.index], [j, op_list[j].site.index]])  
        vals.append(val)
    return np.mean(vals)

def isotropic_corr(psi, Sx, Sy, Sz, L, lat):  # Added lat as an argument
    Cxx = [avg_corr(psi, Sx, L, r) for r in range(L)]
    Cyy = [avg_corr(psi, Sy, L, r) for r in range(L)]
    Czz = [avg_corr(psi, Sz, L, r) for r in range(L)]
    return (np.array(Cxx) + np.array(Cyy) + np.array(Czz)) / 3.0

# --- Compute and plot ---
C_iso = isotropic_corr(psi, Sx, Sy, Sz, L, lat)  # Pass lat to isotropic_corr
r = np.arange(len(C_iso))


plt.figure(figsize=(8, 5))
plt.plot(r, C_iso, 'o-', label='DMRG (PBC): Isotropic $C(r)$')
plt.xlabel('r')
plt.ylabel('C(r)')
plt.title(f'DMRG Isotropic Spin-Spin Correlation $C(r)$ (L = {L}, J2 = {J2})')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

print("\nIsotropic spin-spin correlation C(r):")
for ri, ci in enumerate(C_iso):
    print(f"r = {ri:2d}, C(r) = {ci:.12f}")

Code: Select all

'Array' object has no attribute 'site'

Code: Select all

val = psi.expectation_value_term([[i, op_list[i].site.index], [j, op_list[j].site.index]])

Python: Select all

Sx = [model.lat.mps_sites().Sx for i in range(L)]

Python: Select all

model.lat.mps_sites()