Continuing the discussion from Constant shift in MPS shift-and-invert algo:

I finally have time to come back to this problem. In short, I’m trying to variationally solve the following linear system:

(H - \lambda) |\psi\rangle = |\phi\rangle

Or that

```
o--o--o- -o--o--o--o--o--o <psi|
| | | | | | | | | | |
o--o--o--o--o--o--o--o--o--o--o H
| | | | | | | | | | |
o--o--o--o--o--o--o--o--o--o--o H
| | | | | | | | | | |
o--o--o- -o--o--o--o--o--o |psi>
=
o--o--o- -o--o--o--o--o--o <psi|
| | | | | | | | | | |
o--o--o--o--o--o--o--o--o--o--o H
| | | | | | | | | | |
o--o--o--o--o--o--o--o--o--o--o |phi>
```

(For now let’s assume \lambda=0, since I believe you can always add back the cross terms ( \lambda \langle \psi | H | \psi\rangle , \lambda^2\langle \psi|\psi\rangle, \lambda\langle \psi|\phi\rangle ) which doesn’t change the problem)

(Or there’s a way to add constant terms?)

The way I approach this, is that I defined a new position! function:

```
function position!(P::AbstractProjMPO, psi::MPS, phi::MPS, pos::Int)
L = makeL!(P, psi, phi, pos - 1)
R = makeR!(P, psi, phi, pos + nsite(P))
# start of the left end
if L == nothing
return phi[pos] * phi[pos + 1] * R
# start of the right end
elseif R == nothing
return L * phi[pos] * phi[pos + 1]
# in the middle
else
return L * phi[pos] * phi[pos + 1] * R
end
end
```

with appropriate modification to makeL! and makeR!

```
function _makeR!(P::AbstractProjMPO, psi::MPS, phi::MPS, k::Int)
"""everything were kept the same except"""
R = R * phi[rl - 1] * P.H[rl - 1] * dag(prime(psi[rl - 1]))
"""instead of R = R * psi[rl - 1] * P.H[rl - 1] * dag(prime(psi[rl - 1]))"""
end
```

Eventually I defined a hamiltonian ‘squared’:

```
H2 = contract(H', H; cutoff=1e-12)
H2 = replaceprime(H2, 2 => 1)
```

And I calculate the LHS and RHS and a guess

```
RHS = position!(PH, psi, phi, b)
position!(PH2, psi, b)
guess = psi[b] * psi[b + 1]
```

Put it through linsolve

```
vecs, conv = linsolve(
PH2,
RHS,
guess, 0, 1;
ishermitian=ishermitian,
tol=linsolve_tol,
krylovdim=linsolve_krylovdim,
maxiter=linsolve_maxiter,
)
```

This error initially occurred

`LoadError: DimensionMismatch("In scalar(T) or T[], ITensor T is not a scalar`

After I removed ‘guess’ from the linsolve arguments, another error

```
LoadError: The order of the ProjMPO-ITensor product P*v is not equal to the order of the ITensor v, this is probably due to an index mismatch.
Common reasons for this error:
(1) You are trying to multiply the ProjMPO with the 2-site wave-function at the wrong position.
(2) `orthogonalize!` was called, changing the MPS without updating the ProjMPO.```
```