Extended `NEGF` class

[sofiasanz]

With this commit now one can solve open systems with periodicity along certain direction(s) that aren't the transport direction.

[sofiasanz]

This relates to issue #94.

[sofiasanz]

You might consider some way of having this in files. You'll probably end up having problems with memory when you do multi-electrode wide stuff. So somekind of wrapper class that can either "Buffer" the self-energies, "Calculate on the fly", or "Store in files". Just an idea.

First step would be to move these things into a container, once that is there, changing what a container does should be much simpler.

Thank you for the tip!! Will look into that :)

[sofiasanz]

I tested this extension on graphene with the following script and outputs:

import sisl
import numpy as np
from hubbard import HubbardHamiltonian, density, NEGF, plot

# Set U for the whole calculation
U = 3.0
kT = 0.025

# Build graphene
graphene = sisl.geom.graphene()
H_elec = sisl.Hamiltonian(graphene, spin='polarized')
H_elec[0, 1] = 2.7
H_elec[0, 1, (-1, 0)] = 2.7
H_elec[0, 1, (0, -1)] = 2.7
H_elec[1, 0] = 2.7
H_elec[1, 0, (1, 0)] = 2.7
H_elec[1, 0, (0, 1)] = 2.7

# Hubbard Hamiltonian of elecs
MFH_elec = HubbardHamiltonian(H_elec, U=U, nkpt=[102, 102, 1], kT=kT)
MFH_elec.random_density()

# Converge Electrode Hamiltonians
dn = MFH_elec.converge(density.calc_n, mixer=sisl.mixing.PulayMixer(weight=.7, history=7), tol=1e-10)
MFH_elec.H.write('MFH_elec.nc')

# Central region is a repetition of the electrodes along the periodic direction (1st axis) 
# and not periodic along the transport direction (2nd axis)
HC = H_elec.tile(3, axis=1)
HC.set_nsc([3, 1, 1])

# Map electrodes in the device region
elec_indx = [range(len(H_elec)), range(-len(H_elec), 0)]

# MFH object
MFH_HC = HubbardHamiltonian(HC.H, n=np.tile(MFH_elec.n, 3), nkpt=[100,1,1], U=U, kT=kT)

# First create NEGF object
negf = NEGF(MFH_HC, [(MFH_elec, '-B'), (MFH_elec, '+B')], elec_indx)
# Converge using Green's function method to obtain the densities
dn = MFH_HC.converge(negf.calc_n_open, steps=1, mixer=sisl.mixing.PulayMixer(weight=1., history=7), 
                             tol=1e-6, print_info=True)
print('Nup, Ndn: ', MFH_HC.n.sum(axis=1))

MFH_HC.H.write('MFH_HC.nc')

Output:

   HubbardHamiltonian: converge towards tol=1.00e-06
   1 iterations completed: 5.7087738092320706e-08 -38.95883679505442
   found solution in 1 iterations
Nup, Ndn:  [3.00000034 3.00000034]

Bandstructure of the resulting MFH object:

p = plot.Bandstructure(MFH_elec, bz=[([0., 0., 0.], r'$\Gamma$'), ([2/3, 1/3, 0.], r'K'), 
                                          ([0.5, 0.5, 0.], r'M'), ([0, 0., 0.], r'$\Gamma$')], figsize=(6,4))

Transmission probability output from TBtrans calculated with 100 k-points along the periodic direction (1st axis):

With input file:

SystemLabel device

TBT.HS MFH_HC.nc

%block TBT.Contours
  window
%endblock

%block TBT.Contour.window
 part line
   from -10.0 eV to 10.0 eV
    delta 0.01 eV
     method mid-rule
%endblock

%block TBT.k
100 0 0 0.
0 1 0 0.
0 0 1 0.
%endblock

%block TS.Elecs
  HGNRA
  HGNRB
%endblock TS.Elecs

%block TS.Elec.HGNRA
  TSHS MFH_elec.nc
  semi-inf-dir -a2
  elec-pos begin 1
%endblock TS.Elec.HGNRA

%block TS.Elec.HGNRB
  TSHS MFH_elec.nc
  semi-inf-dir +a2
  elec-pos end -1
%endblock TS.Elec.HGNRB

[sofiasanz]

I think I'll just merge this branch for the moment and try to program it better afterwards :)