vumps_and_plane_wave.py

on github (download).

"""Excitations of the transverse field Ising model.

This example uses VUMPS to find the ground state of the transverse field Ising model and uses the
plane wave excitation ansatz to compute the first excited states.
"""
# Copyright (C) TeNPy Developers, Apache license

import numpy as np

from tenpy.algorithms import plane_wave_excitation, vumps
from tenpy.models.tf_ising import TFIChain
from tenpy.networks import mps


def tfi_vumps(g=1.5):
    model_params = dict(L=2, J=1.0, g=g, bc_MPS='infinite', conserve='parity')
    M = TFIChain(model_params)
    psi = mps.MPS.from_product_state(M.lat.mps_sites(), [0, 0], bc='infinite')
    vumps_pars = {
        'combine': False,
        'N_sweeps_check': 1,
        'mixer': False,
        'trunc_params': {
            'chi_max': 32,
            'svd_min': 1.0e-14,
        },
        'min_sweeps': 30,
        'max_sweeps': 50,
        'max_split_err': 1e-8,  # different criteria than DMRG
        'max_E_err': 1.0e-12,
        'max_S_err': 1.0e-8,
    }

    eng = vumps.TwoSiteVUMPSEngine(psi, M, vumps_pars)
    E, psi = eng.run()
    print(f'ground state energy: {E:.5f}')

    # to directly continue with plane wave excitations, return the uniform MPS
    uniform_psi = eng.psi

    return E, uniform_psi, M


def tfi_excitations(psi_gs, M):
    pw_params = {
        'lanczos_params': {
            'N_max': 50,
        },
    }
    eng_pw = plane_wave_excitation.PlaneWaveExcitationEngine(psi_gs, M, pw_params)

    momenta = np.arange(0, np.pi, np.pi / 8)  # compute for some momenta
    qtotal_change = [1]  # look for excitations in other parity sector
    num_ev = 1  # we only compute the lowest dispersion mode

    dispersions = []
    for p in momenta:
        Es, _, _ = eng_pw.run(p, qtotal_change, num_ev=num_ev)
        dispersions.append(Es)
        print(f'excitation energy for momentum {p / np.pi:.2f} Pi: {Es[0]:.5f}')
    return momenta, np.array(dispersions)


def tfi_dispersion(k, g):
    # exact dispersion for two site unit cell
    return np.min(
        [2 * np.sqrt(g**2 - 2 * g * np.cos(k) + 1), 2 * np.sqrt(g**2 - 2 * g * np.cos(k + np.pi) + 1)], axis=0
    )


if __name__ == '__main__':
    E, psi, M = tfi_vumps()
    momenta, dispersions = tfi_excitations(psi, M)

    # plot and compare to exact results
    import matplotlib.pyplot as plt

    plt.plot(momenta, dispersions, 'x', label='plane wave ansatz')
    plt.plot(np.arange(0, np.pi, 0.1), tfi_dispersion(np.arange(0, np.pi, 0.1), 1.5), ':', label='exact', c='black')
    plt.legend()
    plt.xlabel('momentum')
    plt.ylabel('excitation energy')
    plt.show()