Edit: I have made some progress with this but I'm on holiday at the moment. I will add something in the next couple of weeks.

Edit 2: I came up with the following code to implement the transfer matrix, which can then be used with arnoldi. However it runs quite slowly for larger unit cells and bond lengths. Please note that I've only tested this code briefly and it should be used with caution

```
#ifndef __TRANSFER_H_
#define __TRANSFER_H_
#include "itensor/all.h"
using namespace itensor;
using namespace std;
struct TMatrix {
ITensor left, right;
ITensor T;
int N;
Index link1, link2;
Index left1, left2, right1, right2;
TMatrix(MPS psi1, MPS psi2) {
N = length(psi1);
// Use primes to distinguish some indices
psi1.position(1);
psi2.position(1);
link1 = commonIndex(psi1.A(1), psi1.A(N), {"Link"});
link2 = commonIndex(psi2.A(1), psi2.A(N), {"Link"});
left1 = link1;
left1.prime(1);
left2 = link2;
left2.prime(2);
right1 = link1;
right1.prime(3);
right2 = link2;
right2.prime(4);
T = dag(psi1.A(1)) * delta(link1, left1);
T *= psi2.A(1) * delta(link2, left2);
auto s1 = findIndex(psi1.A(1), "Site");
auto s2 = findIndex(psi2.A(1), "Site");
if (s1 != s2)
T *= delta(findIndex(psi1.A(1), "Site"), findIndex(psi2.A(1), "Site"));
for (int i = 2; i <= N; i++) {
ITensor z = dag(psi1.A(i));
if (s1 != s2)
z *= delta(findIndex(psi1.A(i), "Site"), findIndex(psi2.A(i), "Site"));
if (i == N)
z *= delta(link1, right1);
T *= z;
z = psi2.A(i);
if (i == N)
z *= delta(link2, right2);
T *= z;
}
}
ITensor make_guess() { return randomITensor(right1, right2); }
size_t size() { return dim(right1) * dim(right2); }
void product(ITensor &other, ITensor &result) {
result = other * T;
result *= delta(right1, left1);
result *= delta(right2, left2);
}
};
#endif // __TRANSFER_H_
```

Dear Miles/Matt/ITensor,

I've been using your iDMRG example code (idmrg.h etc) and I've adapted it to some models I'm interested in. My goal is to use the fidelity and fidelity per site (and their susceptibilities) to examine the phase diagram. For a finite system, the fidelity is simply the overlap of the ground states at two different values of some parameter of the Hamiltonian. I'm able to calculate the "overlap" of two iMPS's with some code adapted from your examples - this required some modifications to make the site indices match if they were created from different SiteSet instances. I believe this is appropriate for the fidelity but I'm not sure, and it offers no way to calculate the fidelity per site. I've added this code below.

In McCulloch's, "Infinite size density matrix renormalization group, revisited" (section D. Fidelity), he suggests a method using the eigenvalues of a transfer matrix Eq. (39). I'm having difficulty understanding this section, in particular what 'E' denotes in Eq. (39). I've gone through the references but I haven't been able to decipher it. However, McCulloch seems to imply the procedure is fairly simple.

So, I was wondering if you have any suggestions on how to calculate the fidelity (per site) from iMPS's.

Cheers,

Alex Henry

```
Complex ioverlap(MPS &psi1, MPS &psi2, bool same) {
if (same) {
// if psi1 and psi2 already have the same site indices
auto O = dag(psi1(1)) * psi2(1);
for (auto n : range1(2, length(psi2)))
O = O * dag(psi1(n)) * psi2(n);
auto ovrlap = real(eltC(O));
return ovrlap;
}
auto O = dag(psi1(1)) * psi2(1);
auto l1 = findIndex(dag(psi1(1)), "Site");
auto l2 = findIndex(psi2(1), "Site");
O = O * delta(l1, l2);
for (auto n : range1(2, length(psi2))) {
auto a = dag(psi1(n));
auto b = psi2(n);
l1 = findIndex(a, "Site");
l2 = findIndex(b, "Site");
auto link = delta(l1, l2);
a = a * link;
O = O * a * b;
}
auto overlap = real(eltC(O));
return overlap;
}
```