Thank you so much for the detailed answers, and I am sorry about the barrage of questions! :)
Your answers help a lot. I think I understand the spin version of things. I had some follow up questions about the Kitaev chain.
1. As you mentioned, the correct way to do fermion DMRG is by imposing fermion parity conservation. Thus, I could obtain the fermionic ground states in both the even and odd parity sectors. And then if the energies of both these ground states turn out to be the same (up to some DMRG resolution errors which I would need to estimate), I could conclude that the Hamiltonian is in a topologically non-trivial phase.
2. But I am also interested in another way of diagnosing the topological characteristics of the Hamiltonian. This second method involves adding appropriate tunneling terms which break the ground state degeneracy. One example of such a term would be the bond term connecting the first and last site which imposes periodic/ anti-periodic boundary conditions. I am interested in this method because it allows me to identify the operators which represent edge modes. However, the caveat is that I cannot impose fermion parity conservation anymore. Because if I do, then the imposed fermion parity conservation (and not the tunneling term) dictates which ground state is picked. So, I have a question about how to make this second method work.
Here is what I have noticed so far: When I turn off fermion parity conservation and add appropriate tunneling terms, DMRG picks up a definite parity ground state only if the bond dimension of the initial state psi0 is large enough (and the size of the required bond dimension grows with system size). For example, for a Kitaev chain with 100 sites, I observe that psi0 should have a bond dimension of at least 50. If the bond dimension is any smaller, then DMRG picks up a uniform superposition of the ground states with different parities (even if I increase the number of sweeps/ maximum bond dimension in intermediate sweeps). I don't understand why DMRG chooses ground states in this manner or how this required bond dimension scales with system size. Basically, how is DMRG dealing with these tunneling terms, and how should I set up the DMRG to ensure that the algorithm registers these terms?
Also, thanks for letting me know about Juan's suggestion. I will give it a shot and get back to you in case I face difficulties.