1 import Numeric
   2 from ASE.Transport import TransmissionTools
   3 from ASE.Transport.ScatteringHamiltonianMatrix import ScatteringHamiltonianMatrix
   4 from ASE.Transport.LeadHamiltonianMatrix import LeftLeadHamiltonianMatrix,RightLeadHamiltonianMatrix
   5 
   6 # Construct a lead Hamiltonian matrix describing 2 principal layers of 
   7 # a 1d TB chain with up to second-nearest neighbor hopping                
   8 
   9 H_lead=Numeric.zeros([4,4],Numeric.Complex)
  10 
  11 # On-site energies are zero
  12 for i in range(4):
  13     H_lead[i,i]=0.0
  14 # Nearest neighbor hopping is -1.0
  15 for i in range(3):
  16     H_lead[i,i+1]=-1.0
  17     H_lead[i+1,i]=-1.0
  18 # Next-nearest neighbor hopping is 0.2
  19 for i in range(2):
  20     H_lead[i,i+2]=0.2
  21     H_lead[i+2,i]=0.2
  22 
  23 # Construct the Hamiltonian of the scattering region including 1 principal 
  24 # layer of the lead. The scattering region consists of 2 sites at zero        
  25 # energy and with external coupling of 0.2 and internal coupling og 0.8.      
  26 
  27 H_scat=Numeric.zeros([6,6],Numeric.Complex)
  28 
  29 # Principal layers
  30 H_scat[:2,:2]=H_lead[:2,:2]
  31 H_scat[-2:,-2:]=H_lead[:2,:2]
  32 # Scattering region
  33 H_scat[2,2]=0.0
  34 H_scat[3,3]=0.0
  35 H_scat[2,3]=0.8
  36 H_scat[3,2]=0.8
  37 # External coupling
  38 H_scat[1,2]=0.2
  39 H_scat[2,1]=0.2
  40 H_scat[3,4]=0.2
  41 H_scat[4,3]=0.2
  42 
  43 
  44 # Initializing classes to contain scattering and lead Hamiltonians
  45 # Principal layer is two sites and is the same for left and right leads 
  46 prin=2
  47 hmat=ScatteringHamiltonianMatrix(leftprincipallayer=prin,rightprincipallayer=prin,hamiltonianmatrix=H_scat)
  48 left=LeftLeadHamiltonianMatrix(principallayer=prin,hamiltonianmatrix=H_lead)
  49 right=RightLeadHamiltonianMatrix(principallayer=prin,hamiltonianmatrix=H_lead)
  50 
  51 # Calculate transmission function at energies
  52 en=Numeric.arange(-3,3,0.02)
  53 
  54 #####################################################################################
  55 # The following illustrates how to manipulate the Hamiltonian for analysis          #
  56 # purposes. In particular it is shown how to (i) sub-diagonalize a subspace to      #
  57 # obtain renormalized energy levels in a certain region. (ii) How to cut the        #
  58 # coupling to selected orbitals e.g. some of the renormalized energy levels.        # 
  59 #####################################################################################
  60 
  61 # Sub-diagonalize the subspace of the two sites to obtain two (renormalized) energy levels
  62 from ASE.Utilities.Wannier import HamiltonianTools
  63 H=hmat.GetHamiltonianMatrix()
  64 # H has dimension 6x6 corresponding to the scattering region (two levels).
  65 # The two levels thus have indices 2 and 3
  66 states=[2,3]
  67 # Sub-diagonalize the subspace of the two sites to obtain two (renormalized) energy levels
  68 H2,U,eig=HamiltonianTools.SubDiagonalize(H,states)
  69 print "The renormalized molecular energy levels",eig
  70 # Cut all couplings to one of the two levels. This level is now effectively removed 
  71 # from the calculation, and only the other level will show up as a peak in the transmission.
  72 H3=HamiltonianTools.CutCoupling(H2,[3])
  73 hmat.SetHamiltonianMatrix(H3)
  74 
  75 # Initialize Transmission object
  76 trans=TransmissionTools.TransmissionTool(hmat,left,right,en,prin)
  77 trans.InitGreenFunctions()
  78 
  79 # We want to calculate: 
  80 # the transmission function and the spectral functions of all sites in scattering region
  81 trans.UpdateTransmission()
  82 trans.UpdateSpectralFunctions()
  83 
  84 # Run calculation
  85 t=trans.UpdateToNetCDFFile('TBchain_trans.nc')