Hi Deepak,
I think I see the problem. The implementation of tebd includes optimizations for conserving U(1) symmetries, like particle number. If the on-site irreps of the symmetry are all one-dimensional (as they are for spins conserving magnetization of spinless particle models), then tebd uses a compressed mps format (which in the code means it uses the Gammasm instead of the Gammas). The kink state routine is set up for spins, which use this optimization, but this optimization does not apply for the spinful fermions. However, you can replace Gammasm by Gammas in AllocateGamLam, replace Gammasm by Gammas and remove the LabelLeft and LabelRight in the call to ProductStateMPD, and add
CALL ProductStateLabels(LabelLeft, LabelRight, carray,intDegFree=1)
after the line calling ProductStateMPD. This will use the correct, less compressed structure for the mps. The source code/syntax for all the above calls can be found in StateOps.f90, in the core directory. This should fix the issue.
-Michael