There’s a couple of links on measuring correlations in ITensor.
http://itensor.org/docs.cgi?vers=cppv3&page=formulas/correlator_mps
http://itensor.org/docs.cgi?vers=cppv3&page=formulas/spinless_correlator_mps
Those recipes work when i < j. What should I do to measure Cdagup * Cup or Cdagup*Cdagdn*Cup*Cdn for Spinful fermions and i=j (on the same site)?
Good news: we have recently added a function correlationMatrix to the C++ version of ITensor that measures two-point correlations like \langle c^\dagger_i c_j \rangle for you, including handling details of fermions. (There is a similar function correlation_matrix for the Julia version too.) You can read about how to use it here:
http://itensor.org/docs.cgi?vers=cppv3&page=classes/mps_mpo_algs
Regarding implementing it yourself, if you want to work out how the i > j case relates to te i < j case you can use the following kind of logic. For i < j:
So you can just compute the correlator for i < j (thought of as a matrix) and fill in the other i > j values by taking the conjugate transpose.
For the diagonal i = j values you can just measure the n_i = c^\dagger_i c_i operator on each site.
For four-operator terms, it’s more complicated of course, and you’ll need to have a more complex code and logic about how you advance the tensors from the left and right, but the basic pattern is not terribly different from the two-operator case. Of course there is also the matter of the Jordan-Wigner string or “F” operators which becomes a little more complicated in the four-fermion case.
Thanks a lot. I’m using C++ version of iTensor. Now I’ve come with a function like that:
complex<double> Correlation(MPS &psi, string opname1, int i, string opname2, int j, const SiteSet &sites)
{
// if (j <= i)
// cout << "Error in Correlation: i=" << i << " j=" << j << endl, exit(0);
ITensor op_i, op_j;
for (int o = 1; o <= 2; o++)
{
ITensor &op_s = (o == 1) ? op_i : op_j;
string name_s = (o == 1) ? opname1 : opname2;
int s = (o == 1) ? i : j;
if (name_s == "P") // Pair annihilation operator
{
op_s = op(sites, "Cup", s);
op_s *= prime(op(sites, "Cdn", s));
op_s.mapPrime(2, 1);
}
else if (name_s == "Pdag") // Pair creation operator
{
op_s = op(sites, "Cdagdn", s);
op_s *= prime(op(sites, "Cdagup", s));
op_s.mapPrime(2, 1);
}
else //Other single-site operator defined in ITEnsor for the FermionSite local Hilbert space
{
op_s = op(sites, name_s, s);
}
}
// Are the two operators "fermionic" (i.e. anticommmute) ?
const bool IsFermionic1 = (opname1 == "Cdagup" || opname1 == "Cdagdn" || opname1 == "Cup" || opname1 == "Cdn");
const bool IsFermionic2 = (opname2 == "Cdagup" || opname2 == "Cdagdn" || opname2 == "Cup" || opname2 == "Cdn");
if (IsFermionic1 != IsFermionic2)
cout << "Error: mixing a fermionic and a non-fermion operator." << endl, exit(1);
psi.position(i); //'gauge' the MPS to site i
//Create the <psi| ('bra') from |psi> ('ket')
auto psidag = dag(psi);
//Prime the link indices to make them distinct from
//the original ket links
psidag.prime("Link");
//below we will assume j > i
//index linking i-1 to i:
auto li_1 = leftLinkIndex(psi, i);
auto C = prime(psi(i), li_1) * op_i;
C *= prime(psidag(i), "Site");
// if (i!=j){
if (IsFermionic1)
{
for (int k = i + 1; k < j; ++k)
{
C *= psi(k);
C *= op(sites, "F", k); // Jordan-Wigner fermion 'string' operator F =(−1)^{n}
C *= prime(psidag(k), "Site");
}
}
else
{
for (int k = i + 1; k < j; ++k)
{
C *= psi(k);
C *= psidag(k);
}
}
// }
//index linking j to j+1:
auto lj = rightLinkIndex(psi, j);
C *= prime(psi(j), lj) * op_j;
C *= prime(psidag(j), "Site");
auto result = eltC(C);
return result;
}
Did I get that right that for c^{\dag}_{i,\uparrow}c_{i,\uparrow} I should just calculate
psi.position(i);
auto psidag = dag(psi);
psidag.prime("Link");
auto li_1 = leftLinkIndex(psi, i);
auto C = prime(psi(i), li_1) * op_i;
C *= prime(psidag(i), "Site");
auto lj = rightLinkIndex(psi, j);
C *= prime(psi(j), lj) * op_j;
C *= prime(psidag(j), "Site");
auto result = eltC(C);
right? For i = j it seems like there’s no need in Jordan-Wigner “strings”. The same thing should apply for “pair” creation/annihilation operators c^{\dag}_{i,\uparrow}c^{\dag}_{I,\downarrow} (pair creation) and c_{i,\uparrow}c_{I,\downarrow}(pair annihilation) as pairs are bosons.
If so, I’m getting an error:
Wrong number of IndexVals passed to elt/eltC (expected 0, got 2)
Great that you came up with this function. Not to steal your thunder, but I had forgotten that actually there is a new function called correlationMatrix in the C++ version of ITensor that will compute two-point correlation functions of MPS for you, including handling fermion string.
The documentation on it is here: http://itensor.org/docs.cgi?vers=cppv3&page=classes/mps_mpo_algs
and you can use it like this
auto C = correlationMatrix(psi,sites,"Cdag","C");
As to why you are getting the “Wrong number of IndexVals” error, it means that the final ITensor which you thought was a scalar actually isn’t, and has one or more indices left over on it. So it’s an error with your priming or indexing logic somewhere. I would suggest printing out the various ITensors, starting with C to see if they have the indices on them that you are expecting.
Thanks a lot! It solves half of my problem. Is that possible to input operators other than pre-defined to correlationMatrix? Say, if I have bosonic operators p_i = c_{i,\uparrow}*c_{i,\downarrow} and p^{\dag}_j = c^{\dag}_{j, \uparrow}*c^{\dag}_{j, \downarrow}
Unfortunately no, you can only input predefined operators into correlationMatrix. But it’s not too hard to define your own operators by modifying the site set (I presume you’re using the Spinless site set?), or by making a custom site set based on one of the included ones.
Another option is to take the code from within correlationMatrix and use it as a starting point for your own code. Especially if you sites are all of the same type, then it would be easy to just replace the calls to the op function within correlationMatrix with a function that makes custom ITensors that encode the operators you want.
Finally, maybe you will get tired of us saying this, but all of this is much easier to do in the Julia version of ITensor where you can make custom operators in a few lines of code without having to define a site set, or modifying any ITensor library code or anything like that.
I’m using Electron siteset. Thanks a lot, I got it. May I have an example of the easiest way (you told, that it should be in Julia), that I could implement custom operators , say, p^{\dag}_i = c^{\dag}_i(spin-up)* c^{\dag}_i(spin-down) and p_i = c_i(spin-up)* c_i(spin-down) to plug it the. into correlationMatrix function. I don’t quite get how to contract indices or deal with Jordan/wigner “strings” especially in case of ladder lattice, say, 2 x N or 3 x N
Good question: the way to implement custom operators is to define overloads of the ITensors.op method. For example, if you want to define the operator c^\dagger_\uparrow c^\dagger_\downarrow and give it the name "CduCdd" you write the following code:
function ITensors.op(::OpName"CduCdd", ::SiteType"Electron")
return [
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0
]
end
where I believe I have the matrix elements correct (the only non-zero one being c^\dagger_\uparrow c^\dagger_\downarrow |0\rangle = |\!\uparrow\downarrow\rangle).
Because this operator is bosonic (made up of an even number of c or c^\dagger operators) you don’t have to do anything involving Jordan-Wigner string for it. You do have to think about whether the operator definition has the correct signs for the non-zero matrix elements when acting on a single Electron site, but that’s the only part involving fermions you have to think about.
If you do want to also define an operator which consists of an odd number of creation or annihilation operators, or which acts on multiple sites, please let me know and we can discuss further about how to put in the Jordan-Wigner string correctly for those kind of operators.
Finally, I would always encourage you to test custom operators that you make this way, by preparing specific states as MPS then computing matrix elements of your operators (sandwiching them between two MPS) or inputting them into the expect or correlation_matrix functions for some cases where you know what the answer should be.
Thanks a lot, Miles. I’m on the way to rewriting my code (now C++) in Julia so that I could try that.
Did I get that right that then there hermitian conjugate (which is c_{\uparrow}c
_{\downarrow} should just be a transposed matrix:
function ITensors.op(::OpName"CuCd", ::SiteType"Electron")
return [
0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
]
end
Yes I just checked it and I believe your definition is correct as being the conjugate of your other operator. However I don’t think the name is correct nor the identification of it with c_\uparrow c_\downarrow. Here’s what I worked out:
where note the following two things:
If you do want to also define the operator c_\uparrow c_\downarrow then it is the negative of c_\downarrow c_\uparrow by the rules of fermion anticommutation.
Lastly, I should have mentioned this earlier but ITensor recognizes the following kind of operator input strings: “A * B” where A and B are two valid operator names such as “Cdagup” or “Cup” etc. So you should also be able to input things like “Cdagup * Cdagdn” into OpSum and into correlation_matrix.
Thank you so much. I’m getting some strange results as I try to use
auto C = correlationMatrix(psi,sites,"Cup*Cdn", "Cdagdn*Cdagup");
Cup*Cdn, for instance, on site s should look like:
op_s = op(sites, “Cup”, s);
op_s *= prime(op(sites, “Cdn”, s));
op_s.mapPrime(2, 1);
I’m not quite sure that it is so for correlation_matrix.
It works just fine for standard operators like:
auto C = correlationMatrix(psi,sites,"Cdagup","Cup");
I’m using correlations implementation in ALPS package for verification. Unfortunately results for auto C = correlationMatrix(psi,sites,"Cup*Cdn", "Cdagdn*Cdagup"); do not match and I’m trying to find out why. How could we check that correlationMatrix works correctly for “CupCdn", "CdagdnCdagup” or other pairly multiplied operators?
I don’t know if it will solve the issue you’re seeing, but I think you have the order reversed about multiplying the operators (your op_s code). When writing “Cup * Cdn” what is meant is that “Cdn” acts on the state first, then “Cup”, so like “Cup * Cdn” acting on a state \Psi means c_\uparrow c_\downarrow |\Psi\rangle.
If you draw the tensor diagrams for you piece of code, you’ll see that “Cup” has an unprimed index coming out, meaning it can connect to the site index of the state, whereas “Cup” is “stacked” on top (or to the left, depending on how you draw it) so it acts second.
Another way you could check whether you are getting the right results, instead of comparing to ALPS, the simplest thing would just be to prepare an input state psi for which you know what the answer should be by some analytical calculation. Like it could be a wavefunction you got from solving a non-interacting fermion problem. Do you know of a good example like that?
ITensor Website
·
Twitter @ITensorLib
© 2022 ITensor Collaboration