TDVP time evolution with interaction / onsite state probability

How do I use this algorithm? What does that parameter do?
Post Reply
Jirawat S
Posts: 4
Joined: 01 Mar 2022, 14:12

TDVP time evolution with interaction / onsite state probability

Post by Jirawat S »

Hi,

I am trying to simulate a time evolution of driving 2 spin-1 system with external field by using TDVP method. Starting from ground state, when the coupling term did not present, the system has Rabi oscillation in both spin. However, once I added the coupling term in form of (Sz.Sz), I noticed that probability of the first spin to be in ground state (\(P|0>\)) is just flat to around zero. This doesn't happened to the second spin.

I believed there was a bug in TDVPEngine, but I also got the same result with ExpMPOEvolution method. What should I do with this?

1st spin
c9916ae6-9d34-4521-ab0d-b85bde4c3b2c.png
c9916ae6-9d34-4521-ab0d-b85bde4c3b2c.png (21.86 KiB) Viewed 192 times
2nd spin
e700f2fd-6632-4d02-a26c-ebda8bdc5601.png
e700f2fd-6632-4d02-a26c-ebda8bdc5601.png (35.88 KiB) Viewed 192 times
This is the code that I use.

Code: Select all

import numpy as np
from tqdm import tqdm
import matplotlib.pyplot as plt
from tenpy.networks.site import SpinSite
from tenpy.models.lattice import Chain
from tenpy.networks.mps import MPS
from TenPy_Model import TenPy_Dipole
from tenpy.algorithms import tdvp
from tenpy.models.model import CouplingModel, NearestNeighborModel, MPOModel

def get_probability(psi, site_idx:int):
    coeff = psi.get_B(site_idx).to_ndarray().reshape(-1)
    return (np.conjugate(coeff) * coeff).real

class TenPy_Dipole(CouplingModel, NearestNeighborModel, MPOModel):
    def __init__(self, lattice):
        self.lattice = lattice

        CouplingModel.__init__(self, self.lattice)
        self.add_onsite(2870., 0, 'Sz Sz')  #ZFS
        self.add_onsite(185., 0, 'Sz')     #Zeeman
        self.add_onsite(-2685., 0, 'Sz Sz')  #MW
        self.add_onsite(2*np.pi*15.0, 0, 'Sx')       #MW
        self.add_coupling_term(1.0, 0, 1, 'Sz', 'Sz', plus_hc=False)
        MPOModel.__init__(self, self.lattice, self.calc_H_MPO())

#Parameter for time evolution.
T = 50e-3 #(us)
dt = 1e-3 #(us)
steps = round(T/dt)
cutoff = 1.e-10
tdvp_params = {
        'start_time': 0,
        'N_steps':1,
        'dt': dt,
        'trunc_params': {
            'chi_max': 10,
            'svd_min': 1.e-10,
            'trunc_cut': None
        }}

# predefined charges and Sp,Sm,Sz operators
S = 1.0        #Total spin number.
L = 2           #Number of defects.
site = SpinSite(S=S, conserve='None')
lat = Chain(L, site, bc_MPS='finite')
state = ["0.0","0.0"]
psi = MPS.from_product_state(lat.mps_sites(), state, lat.bc_MPS)

which_site = 0

times = []
Pg, Pm, Pp = [], [], [] #Population 
for i in tqdm(range(steps)):
    #Getting probability in |0> state of spin 1.
    prob_array = get_probability(psi,which_site)
    Pg.append(prob_array[1])
    Pm.append(prob_array[0])
    Pp.append(prob_array[2])

    Model = TenPy_Dipole(lat)
    engine = tdvp.TDVPEngine(psi, Model,options=tdvp_params) 
    engine.run()
    times.append(dt*i)
    psi = engine.psi

fig = plt.figure()
plt.plot(times, Pg, label='P(|0>)')
plt.plot(times, Pm, label='P(|-1>)')
plt.plot(times, Pp, label='P(|+1>)')
plt.legend()
plt.show()
Jirawat S
Posts: 4
Joined: 01 Mar 2022, 14:12

Re: TDVP time evolution with interaction.

Post by Jirawat S »

Okay, after searching for an old post again, I've found an answer. So I think it could be beneficial to the others who has the same problem as me.

The reason of the problem is that I had been using incorrect way to calculate the state probability. One can find the correct method in this post
User avatar
Johannes
Site Admin
Posts: 316
Joined: 21 Jul 2018, 12:52
Location: UC Berkeley

Re: TDVP time evolution with interaction / onsite state probability

Post by Johannes »

Indeed, the get_probability is not working; it even gives you wrong dimensions of the array when you go beyond 2 sites.

Using the onsite reduced density matrix is as pointed out in the post you linked to is indeed one way.

Another and maybe even clearer solution would be to define projection operators, e.g. \(P_{+1} = \ket{+1}\bra{+1} \), and measure the expectation values of those.

Expanding your example:

Code: Select all

import numpy as np
from tqdm import tqdm
import matplotlib.pyplot as plt
from tenpy.networks.site import SpinSite
from tenpy.models.lattice import Chain
from tenpy.networks.mps import MPS
#  from TenPy_Model import TenPy_Dipole
from tenpy.algorithms import tdvp
from tenpy.models.model import CouplingMPOModel, NearestNeighborModel


class TenPy_Dipole(CouplingMPOModel, NearestNeighborModel):

    def init_sites(self, model_params):
        S = 1.0        #Total spin number.
        site = SpinSite(S=S, conserve='None')
        for label, opname in zip(['-1.0', '0.0', '1.0'], ['Pm', 'P0', 'Pp']):
            diag = np.zeros(3) #  diagonal
            diag[site.state_index(label)] = 1.
            op = np.diag(diag)  # 3x3 numpy array
            site.add_op(opname, op)
        return site

    def init_terms(self, model_params):
        self.add_onsite(2870., 0, 'Sz Sz')  #ZFS
        self.add_onsite(185., 0, 'Sz')     #Zeeman
        self.add_onsite(-2685., 0, 'Sz Sz')  #MW
        self.add_onsite(2*np.pi*15.0, 0, 'Sx')       #MW
        self.add_coupling_term(1.0, 0, 1, 'Sz', 'Sz', plus_hc=False)


#Parameter for time evolution.
T = 50e-3 #(us)
dt = 1e-3 #(us)
steps = round(T/dt)
cutoff = 1.e-10
tdvp_params = {
        'start_time': 0,
        'N_steps':1,
        'dt': dt,
        'trunc_params': {
            'chi_max': 10,
            'svd_min': 1.e-10,
            'trunc_cut': None
        }}

model_params = {
    'L': 2,
    'lattice': 'Chain',
    'bc_MPS': 'finite',
}

Model = TenPy_Dipole(model_params)
state = ["0.0","0.0"]
psi = MPS.from_product_state(Model.lat.mps_sites(), state, Model.lat.bc_MPS)
engine = tdvp.TDVPEngine(psi, Model,options=tdvp_params)

times = []
P0, Pm, Pp = [], [], [] #Population
for i in tqdm(range(steps)):
    #Getting probability in |0> state of spin 1.

    P0.append(psi.expectation_value('P0'))  # two values, one for each site!
    Pp.append(psi.expectation_value('Pp'))
    Pm.append(psi.expectation_value('Pm'))

    engine.run()  # evolve by dt * N_steps
    times.append(engine.evolved_time) # same as dt*i)
    #  psi = engine.psi


P0 = np.array(P0)
Pm = np.array(Pm)
Pp = np.array(Pp)

fig = plt.figure()
for i in range(psi.L):
    linestyle = '--' if i == 1 else '-'
    plt.plot(times, P0[:, i], linestyle=linestyle, label=f'P(|0>), site {i:d}')
    plt.plot(times, Pm[:, i], linestyle=linestyle, label=f'P(|-1>), site {i:d}')
    plt.plot(times, Pp[:, i], linestyle=linestyle, label=f'P(|+1>), site {i:d}')
plt.legend()
plt.show()

Post Reply