Hi Alex,

The trick is to add the Peierls phase after you discretized your Hamiltonian, I do it in the following way:

def apply_peierls_to_template(template, xyz_offset=(0, 0, 0)):
    template = deepcopy(template)  # Needed because kwant.Builder is mutable
    x0, y0, z0 = xyz_offset
    lat = template.lattice
    a = np.max(lat.prim_vecs)  # lattice contant

    def phase(site1, site2, B_x, B_y, B_z, orbital, e, hbar):
        x, y, z = site1.tag
        direction = site2.tag - site1.tag
        A = [B_y * (z - z0) - B_z * (y - y0), 0, B_x * (y - y0)]
        A = np.dot(A, direction) * a**2 * 1e-18 * e / hbar
        phase = np.exp(-1j * A)

        if lat.norbs == 2:  # No PH degrees of freedom
            return phase
        elif lat.norbs == 4:
            return np.array([phase, phase.conj(), phase, phase.conj()],
                            dtype='complex128')

    for (site1, site2), hop in template.hopping_value_pairs():
        template[site1, site2] = combine(hop, phase, operator.mul, 2)

    return template

You can apply this function to your template builder. In short, it does this combine(hop, phase, operator.mul, 2), which is creating a new function from the phase function and the hopping function and combining the arguments.

You would probably need to adjust the form of the vector potential, the units (this is using nm) and change the form of the phase function.

I have attached the combine.py script you will need to combine the functions.

Please let me know whether it works or not.

Best, Bas