Page 1 of 1

Kinetic energy

Posted: 09 Mar 2020, 12:57
by 011235
Hi! I think I have an easy question but not sure how to do it. I am running iDMRG algorithm to find the ground state of the system I am interested in. Apart from the total energy density, is there a built-in function in Tenpy to obtain the contribution to the total energy of each of the terms in the Hamiltonian separately? E.g., if H= T + V, what is the correct way to use the defined model to access eg the kinetic energy T?

My example code:

M = FermionModel(model_params)
psi = MPS.from_product_state(M.lat.mps_sites(), product_state, bc=M.lat.bc_MPS)
result = dmrg.run(psi, M, dmrg_params)
Energy = result['E']

Thanks!

Re: Kinetic energy

Posted: 11 Mar 2020, 15:11
by Leon
Hi, and welcome to the forum! In Tenpy, you can compute expectation values of any operator you'd like using psi.expectation_value(). You can read the documentation here: tenpy.networks.mps.MPS

Re: Kinetic energy

Posted: 11 Mar 2020, 16:36
by Johannes
I think the question aims at how to get the terms contained in the kinetic energy.
Actually, there's some functionality in the tenpy.models.model.CouplingModel, which is not explicitly used in most of the pre-defined models, though, and a bit hidden in the documentation.
When one calls add_coupling, it takes a keyword argument category.
The attribute coupling_terms of the model is simply a dictionary mapping the category string to instances of tenpy.networks.terms.CouplingTerms.
You basically want to take a subset of all the CouplingTerms and build an MPO from it to measure the expectation value of the MPO.

Full example:

Code: Select all

import tenpy
from tenpy.models.xxz_chain import XXZChain
from tenpy.networks.terms import CouplingTerms, OnsiteTerms
from tenpy.networks.mps import MPS
from tenpy.networks import mpo

L = 6
M = XXZChain({'L': L, 'bc_MPS': 'finite'})
print("coupling term categories:", M.coupling_terms.keys())
ct = M.coupling_terms['Sp_i Sm_j']
print("contained terms: ")
print(ct.to_TermList())
ot = OnsiteTerms(L)

T_MPO_graph = mpo.MPOGraph.from_terms(ot, ct, M.lat.mps_sites(), M.lat.bc_MPS)
T_MPO = T_MPO_graph.build_MPO()
T_MPO.max_range = ct.max_range()
# define psi as nearest-neighbor singlets
psi = MPS.from_singlets(M.lat.unit_cell[0], L, [(0, 1), (2, 3), (4, 5)])
print("<psi|T_MPO|psi> = ", T_MPO.expectation_value(psi))
assert abs(T_MPO.expectation_value(psi) - (-1.5)) < 1.e-13