I checked your program and it seems fine and the result should be 'almost' correct.

However, I would like to make some comments and corrections:

1) You are using wraparound and your system if uniform in the y-direction. No need to use a width W. A one unit cell length is enough.

This may seem at a first glance with no consequences, but it does change the definition of your ky such that two people using two different widths 

will have two different values of T(ky, E)  although the integrated value T(E) should be the same (if integrated on the correct interval in each case). 

2) Let us first be interested on the sum of T(ky, E) over some chosen continuous modes ky and discretized open modes in z-direction. Let us also put the same two leads at the same potential and put the L=1. This will help us in doing analytic calculation. 

The energy of each mode is:

 Em=6*t-2*t*cos(ky)-2*t*cos(m*pi/(H+1))  

There is no inter-mode scattering in your system and the transmission can be obtained using the Fisher-Lee formula as for a scatterer in 1D:

T(ky,E)=sum_m    1/(1.+(V0/Gamma)**2)     with Gamma=sqrt(4-(E-Em)**2).

This formula is tested in the below code.

If you change W>1, it will not work for the reason I mentioned before: some re normalization of the ky  (ky/W) and changes should be done probably. I tried to figure out how to do it but I failed. 

If you put V0 in one lead (as in your program) a quick calculation can give you the analytical result  in the same way as here (you will need Gamma1, Gamma2)

I hope this helps

Adel

############################

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from numpy import *
import kwant
import sys

#calculate the conductance for a sheet (Dim=2) or a Slab (Dim=3) that is translational invariant in y direction

a = 1.0              #lattice constant
V0 = 0.2            #potential step size

L = 1                #Length  x-direction  the length is a bit useless in this code
W = 1             #Width    y-direction
H = 4               #Hight     z-direction
t = 1                #hopping parameter

Dim=3          #dimension of the system
                    # for H=1, Dim=2
                    # for H>1, Dim=3
#params for integration and plotting
Nky=100             #points in ky space
NE=30                #points in energy
lat = kwant.lattice.cubic(a)


# Define the scattering region
# Infinite potential plane/Slab in y direction
sys = kwant.Builder(kwant.TranslationalSymmetry(lat.vec((0, W, 0))))


for i in range(L):
    for j in range(W):
        for k in range(H):
        # On-site Hamiltonian
            sys[lat(i, j, k)] = 2*Dim * t + V0                   #potential step located at the center between the leads

sys[lat.neighbors(1)] = -t

sys = kwant.wraparound.wraparound(sys)

#Left lead
leadL = kwant.Builder(kwant.TranslationalSymmetry((-a, 0, 0), lat.vec((0, W, 0))))  #trans inv in transport dir and ky
leadL[(lat(0,j,k) for j in range(W) for k in range(H))] = 2*Dim*t           #onsite terms in the left  lead
leadL[lat.neighbors(1)] = -t

#Right lead
leadR = kwant.Builder(kwant.TranslationalSymmetry(( a, 0, 0), lat.vec((0, W, 0))))
leadR[(lat(0,j,k) for j in range(W) for k in range(H))] = 2*Dim*t        #onsite terms in the right lead
leadR[lat.neighbors(1)] = -t

leadL = kwant.wraparound.wraparound(leadL, keep=0)  #keeping the translational invariance in x direction
leadR = kwant.wraparound.wraparound(leadR, keep=0)

sys.attach_lead(leadL)
sys.attach_lead(leadR)

kwant.plot(sys)
sys = sys.finalized()

# Calculation of the total transmission for a given Energy
ky_array = np.linspace(-pi,pi,Nky)
dky=ky_array[1]-ky_array[0]
energies = np.linspace(0.0, 5.0,NE)
Ttot=[]

for energy in energies:

    Tky_array=zeros(Nky)

    for i in arange(len(ky_array)):
        ky=ky_array[i]
        smatrix = kwant.smatrix(sys, energy, [ky])        
        Tky_array[i]=smatrix.transmission(1, 0)
    #Ttot.append((2*pi)/(W)*(sum(Tky_array)-0.5*Tky_array[0]-0.5*Tky_array[Nky-1])*dky) #cheap trapezodial rule to integrate over ky
    Ttot.append(sum(Tky_array)) #cheap trapezodial rule to integrate over ky


# Plot transmission
plt.plot(energies, Ttot)
plt.show()

def T(E,qy,W,H,V0,t=1):
    T0=0
    for m in range(1,H+1):
        Em=6*t-2*t*cos(qy/W)-2*t*cos(m*pi/(H+1))      
       
        if abs(E-Em)<=2: # condition for open  modes
             Gamma=sqrt(4-(E-Em)**2)
             T0+=1/(1.+(V0/Gamma)**2)  ##use Fisher Lee Formula T=Gamma G Gamma G^\dagger
    return T0



Trans=[sum(np.array([T(E,qy,W,H,V0,t=1) for qy in ky_array ] )) for E in energies]      
                   
                   
plt.plot(energies,Trans,'ro')
plt.plot(energies, Ttot)
plt.show()