1.) The toric code model acutally didn't allow for that, you would have had to override the
method. I've just updated it in 4cc12c0ab9f38cc970da7e5338a1ee3210964997
to allow torus-boundary conditions with
bc_MPS="finite", bc_x="periodic", bc_y="periodic"
You can see for Ly=4, that you need bond dimension 64 for the torus compared to bond dimension 8 for the infinite cylinder....
Also, I've now implemented the tenpy.models.lattice.IrregularLattice
, in case you want to remove boundary sites not included into couplings for the finite case.
2.) A gapless system formally always requires infinite MPS bond dimension for the ground state, even in 1D!
Start off with small circumferences. Threading a flux through the cylinder can cause a gap opening, where DMRG converges better, but of course that also changes the physics! Same holds for tuning to a different phase by changeing other parameters (e.g., Jz > Jx + Jy in your case, or adding a magnetic field) where there is a gap.
For the Kitaev Honeycomb model, I'd guess that the best order is
(but I didn't try it...).
You might want to do a scaling analysis with the MPS bond dimension, i.e. run DMRG multiple times at increasing bond and look how results change. This kind of analysis still allows to infere quite a lot of information, even if you never get the fully converged ground state, c.f. arXiv:0812.2903
Changing the initial product state will not change what the required MPS bond dimension is (unless you change the charge sector...), the most it can do is converge in less sweeps if the guess is good.
As I said before, there is no general "magic trick" to make it compuationally easy; if there were, we wouldn't need to develop quantum computers....
For this particular model, of course, you can do all calculuations analytically, making the "compuations" trivial