+1 vote
asked by (190 points)
edited by

Hello!

I'd like to ask a question about constructing the MPO operator being the product of local operators. In my problem I consider two types of bosons "a" and "b" with total number @@N = N_a + N_b@@ conserved, but @@N_a - N_b@@ is not conserved. Total number of sites is @@M = N@@.

I model my system by creating an artificial lattice with @@M = 2N@@ sites, and each two lattice sites form a unit cell where the first site corresponds to boson @@a@@ and second to boson @@b@@.

On such a lattice I can define spin operators @@S^x_i, S^y_i, S^z_i@@ as follows:

auto tmp_Sx = AutoMPO(sites);
auto tmp_Sy = AutoMPO(sites);
auto tmp_Sz = AutoMPO(sites);

for (int j = 1; j < M_sites; j+=2){
  tmp_Sx +=  0.5, "Adag",j,"A", j+1; 
  tmp_Sx +=  0.5, "Adag",j+1,"A", j; 

  tmp_Sy +=  0.5/Cplx_i, "Adag",j,"A", j+1;
  tmp_Sy += -0.5/Cplx_i, "Adag",j+1,"A", j;

  tmp_Sz +=  0.5, "N", j;
  tmp_Sz += -0.5, "N", j+1;
}
auto Sx = toMPO(tmp_Sx);
auto Sy = toMPO(tmp_Sy);
auto Sz = toMPO(tmp_Sz);

I'm interested in calculating expectation value of the @@D = \prod_j d_j@@
operator, where @@d_j = S^y_{j} + i S^z_{j}@@. I've tried the following construction (let's assume N = 2 particles):

vector<MPO> d_vec(N);  
int j = 0;
for(int i = 1; i<M; i+=2){ 
  auto tmp_d_i = AutoMPO(sites);
  tmp_d_i +=  0.5/Cplx_i, "Adag",i,  "A", i+1;
  tmp_d_i += -0.5/Cplx_i, "Adag",i+1,"A", i;
  tmp_d_i +=  0.5*Cplx_i, "N", i;
  tmp_d_i += -0.5*Cplx_i, "N", i+1;
  d_vec[j] = toMPO(tmp_d_i);    
  j = j+1;
}

auto D = nmultMPO(d_vec[1],prime(d_vec[0]), "Method=",nmultMPO_method,"MaxDim",nmultMPO_MaxDim,"Cutoff",nmultMPO_cutoff});

However, it doesn't work (I have exact-diagonalization results to compare).

Can someone point me to what I'm missing? Thanks in advance!

Best,
Marcin

1 Answer

0 votes
answered by (70.1k points)

Hi Marcin,
Thanks for the question. It's helpful to see the structure of your operator, which as I understand is a product of single-site operators. So then actually you should not have to use an MPO at all, or any machinery to sum operators or multiply MPOs together.

(1) Instead you can just apply each operator onto each tensor of your MPS to compute the product D|psi>. If you call the resulting (unnormalized) state |phi> = D|psi> you can then compute the overlap <phi|psi> to obtain the expectation value <psi|D|psi>.

(2) Alternatively you could think of D as an MPO, but which has a bond dimension of 1, so it is just a product of operators. In that approach, AutoMPO is not really the best tool to use to make the MPO. Instead it would be better to make the MPO "by hand" by setting each MPO tensor one by one.

I would recommend approach #1 above. To perform the action of D onto |psi> you could write code similar to in this "code formula" page to apply each operator d_j onto your MPS:
http://itensor.org/docs.cgi?vers=cppv3&page=formulas/mps_onesite_op
By looping over all j from 1 to N you will end up with the state |phi> = D|psi>.
Then you can just use the function inner to compute inner(phi,psi) to get your expectation value.

Hope that helps - we can discuss more if you have questions.

Miles

commented by (190 points)
edited by
Hi Miles,

thanks for your reply! In fact, the local operator @@d\_j@@ is a two-site operator acting on neighboring sites, so I applied an approach from http://itensor.org/docs.cgi?vers=cppv3&page=formulas/gate, i.e.

    //|phi> = D|psi> = \prod_{j} d_j |psi>, d_j acts on {j,j+1} sites
    auto phi = psi;
    for (int j = 1; j < M_sites;  j+=2) {        
    //application of two-site gate to MPS: http://itensor.org/docs.cgi?vers=cppv3&page=formulas/gate
    phi.position(j);  
    auto d_j = 0.5/Cplx_i*( sites.op("Adag",j)*sites.op("A",j+1) - sites.op("Adag",j+1)*sites.op("A",j));
    d_j = d_j + 0.5*Cplx_i*( sites.op("N",j)*sites.op("Id",j+1) - sites.op("N",j+1)*sites.op("Id",j) );
    auto phi_j_j_plus_1 = phi(j)*phi(j+1);              
    phi_j_j_plus_1 *= d_j;
    phi_j_j_plus_1.noPrime();
    auto [U,S,V] = svd(phi_j_j_plus_1,inds(phi(j)),{"Cutoff=",1E-8});
    phi.set(j,U);
    phi.set(j+1,S*V);      
    }   

Where "Id" operator is defined by me in boson.h file as

    if(opname == "Id" || opname == "id"){
      for(auto n : range1(1+maxOcc)){
        Op.set(s=n,sP=n,1);
       }
    }

Is there something missing from my side?

Thanks for help
Marcin
commented by (70.1k points)
Hi Marcin, that's helpful info. You could see why I thought d_j might act only on site j. Apart from the code itself, I would say to definitely be sure that you have taken into account that the operators on overlapping sites may not commute, so to make sure that you are applying them in an ordering that gives the result you want on paper. I'm sure you've thought of that but just mentioning for completeness.

I think you mentioned a bug in an earlier version of your comment? Hope you have found it now since I don't see a mention of that anymore.

But overall, yes this looks like a great approach for doing what you want.

Best,
Miles
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.

Categories

...