Initialising MPS or MPO manually

Question

Dear ITensor team,

this is a kind of follow-up of my previous question on manually constructing MPOs in ITensor v3 c++. It appears to me that it is not well documented the general way to construct generic MPS, and it would be really useful to have it written somewhere. :)

Say one wants to generate an MPS

psi = A^sigma1{b1} A{b1,b2}^sigma2 A{b2,b3}^sigma_3 ....

with b1,b2 indices of dimension = dim and sigma local site indices and A{bi,bi+1}^sigmai matrices we know. The main problem to me is to assign correct indices that are correctly contracted and properly represent the MPS (or MPO, it's the same). The way I do it is the following:

psi = MPS(sites)
auto comInd = Index(dim,"Link"); // define a common index  

for(int n = 1; n <= N; n ++){ // run over the chain of size N

  auto sigma = sites(n);

    if(n==1){

        auto j = Index(dim,"Link");
            auto A = ITensor(sigma,i);
    // here set the values of the matrix elements of A        
           psi.set(n,A);
           comInd = j; // set the common index to the right index of A 
}

    if(n>1 && n< N){
        auto i = comInd;
        auto j = Index(dim,"Link");
       auto A = ITensor(sigma,i,j);
    // here set the values of the matrix elements of A        
       psi.set(n,A);
       comInd = j; // set the common index to the right index of A 
}


    if(n==N){
    auto i=comInd;
auto A = ITensor(sigma,i);
    // here set the values of the matrix elements of A        
    psi.set(n,A);
 }


 } // end of for over sites

I would like to know if this is the correct way to do it or if there are other ways or simpler methods!

Thanks a lot!
All the best,
Jacopo.

Answer 1

Hi Jacopo,
This is a good question, but the short answer is that basically, yes, this is the way to make a more general or arbitrary MPS with custom entries in ITensor.

You can of course come up with some small variations on this approach using existing ITensor features, such as making just a two-index ITensor which is like the "matrix" in the name "matrix product state" then outer-producting on the physical index afterward, or variations like that which might be suitable for diffeent cases, but your code is a good approach too and the most general.

One thing I would recommend though is making all of the bond or link indices ahead of time and storing them in a std::vector. Then it's simpler to just grab the ones you need in each part of the code rather than having a logic about saving the previous index to a variable outside the loop, which could introduce a subtle bug if done wrongly. (Not a bug to do with ITensor, but just to do with that style of programming.)

Please let us know if you either have a suggestion for an improvement to the ITensor interface or features which could make this task easier, or if there's something you find you can't do this way or is excessively hard to do this way.

(Noting that we are aware it would be nice if one could "slice" into an ITensor i.e. access whole subarrays or subtensors at a time. We are planning this feature to go in this year, though in the Julia language version of ITensor.)

Miles