A primer on Smeagol Víctor García Suárez
A primer on Smeagol http://www.smeagol.tcd.ie Víctor García Suárez Outline 1) Introduction 2) Theory 3) How to run the code 4) Simple examples 5) Some calculations 1) Introduction Introduction to Smeagol S pin and M olecular E lectronics in an A tomically G enerated O rbital L andscape Why Smeagol Need for smaller electronic devices. Atomic limit - Faster - Cheaper - More compact - Similar features as today electronic elements (rectification, NDR, etc). General features of Smeagol - Density functional theory (DFT) First-principles code based on localized functions: Siesta, Fireball, etc + - Non-equilibrium Green’s functions Calculation of the density matrix, transmission and current under finite biases Smeagol Characteristics and Capabilities Smeagol and spintronics: exploit the spin degree of freedom - Spin polarized - Non collinear - Spin orbit Calculations of both extended and isolated systems: G and k-point calculations Calculations of more than 100 atoms Parallelized 2) Theory The philosophy behind a Smeagol calculation V Low resistance Large resistance Low resistance Direction of electronic transport (z-axis) Left bank (reservoir) Left electrode Junction Right bank (reservoir) Extended molecule Right electrode Leads and extended molecule Left electrode Molecule Right electrode HL [nL] – e mL mL m – eV/2 No m No Fermi distribution HR [nR] – e mR mR = m + eV/2 Bulk left electrode Left reservoir Extended molecule Non-equilibrium Bulk right electrode Right reservoir H LM 0 H L eV / 2 H H ML H M [nM ] H MR 0 H RM H R eV / 2 H M [nM ] ? DF ( E ) f ( E ) Procedure to compute nM when the distribution function is not Fermi nM (r ) 2 ψ ( r i ) f (εi ) i Non-equilibrium Green’s function formalism Calculation of the leads properties - Unit cell that is repeated along the transport direction (z) - Use of k-points. Necessary to converge the density of states H1 H0 H1 H0 H1 H0 H1 H1 H0 H0 H1 H0 H1 H0 H0 - At the end of the calculation: * Hamiltonian and overlap matrices (H0, S0 and H1, S1) → Surface GF * Density matrix * Fermi energy Self-energies and G matrices The self-energies are calculated with the couplings and the surface GF G0R is the retarded Green’s function of the leads, calculated using a semianalytic formula 0R ˆ ˆ Hˆ ] ˆR [( E i )Sˆ Hˆ ] G ( E )[( E i ) S L/R M L/R L/R L/R M The G matrix contains information on the coupling between the extended molecule and the leads It is important in the calculation of the transmission ˆG ( E) i[ˆ R ( E) ˆ R ( E)] L/R R/L R/L Calculation of the surface GF Since the GF depends on energy it is necessary to calculate k(E) from the block vectors instead of E(k) → Solve the inverse secular equation [ K 1e ik K 0 K1e ]k 0 ; K i H i ESi ik This involves obtaininig the inverse of the K1 matrix, which can not always be inverted → Identify the singularities (GSVD) and get rid of them (decimation) Once the bulk GF is constructed the surface GF is obtained by applying the appropriate boundary conditions The surface GF should vanish at z1 → Add a wave function to the bulk GF z-1 z0 z1 Calculation of the extended molecule - Includes molecule or central scattering region + part of the leads - The first and last part of the EM must coincide with the unit cells of each of the leads (buffer layers). This implies that: * The same general parameters as in the bulk calculation have to be used in the unit cell: temperature, mesh, perpendicular k-points, spin-polarization, etc * The same particular parameters for the leads have to be used in the buffer layers and rest of the leads in the EM: basis set, atomic coordinates, etc. Buffer layers Buffer layers - The Hamiltonian, overlap and DM of the buffer layers is substituted by those of the leads - The buffer layers ensure that: * The electronic structure at the beginning and end of the EM is that of the leads * (forced) Convergence between the electronic strucrure of the leads near the scattering region and deep into the leads * Absence of spurious effects in case of surfaces TEST: infinite system Energy, charge and transmission T(E) E Equal leads - In case of equal leads it is possible to make the calculation of the EM periodic to avoid the presence of surfaces. This does not mean that the system is periodic This does not mean that the system is periodic No k-points are necessary Energy mismatch between the bulk and EM calculations In general the energy origin in a calculation of an infinite system is arbitrary It is necessary to determine such mismatch and correct it. Otherwise the Fermi energy can be incorrectly defined and the system can win or loose charge The correction is made at certain positions of the bulk slices It is necessary to calculate the bulk Hartree potential X X Calculation of the DM Siesta: Solve eigenvalue problem. Order-N or diagonalization * ˆ ˆ Hc ESc ˆ cˆn cˆn f n n Smeagol: * Semi-infinite leads * Non-equilibrium charge distribution 1 ˆ ˆ G ( E ) ˆ dE G ( E ) 2i Density matrix in equilibrium In equilibrium it is only necessary to know the retarded Green’s function Built with the Hamiltonian, overlap and self-energies R R R 1 ˆ ˆ ˆ ˆ ˆ G ( E ) [( E i ) S H R ( E ) L ( E )] The density matrix is obtained by integrating along an energy axis 1 R ˆ ˆ dE[Im G ( E )] f ( E m ) Complex energy contour The lesser Green’s is calculated on an energy contour in equlibrium Three parts: imaginary circle, imaginary line and Fermi poles R ˆ ˆ ˆ dE G ( E ) dE G ( E ) 2ik T z Im[ G ( zn )] C n Fermi function poles zn i(2n 1)kT Image from Atomistix www.quantumwise.com Brandbyge et al. Phys. Rev. B 65, 165401 (2002) Density matrix out of equilibrium Out of equilibrium it the GF is not analytic inside the contour The calculation of the DM is divided in two parts ˆ ˆ L/R ˆ V ˆ L/R 1 R ˆ dE[Im G ( E )] f ( E m L/R ) C 1 ˆ V 2 Rˆ R ˆ ˆ dE [G GR/LG ]( E )[ f L/R f R/L ]( E ) Out of equilibrium voltage profile The Hartree potential is defined up to a constant and a linear ramp (solution of the Poisson equation) A linear ramp related to the bias voltage is added to help the convergence out of equilibrium +eV/2 mL mR V -eV/2 Inside the Siesta version of Smeagol Initial guess Calculate effective potential Smeagol Calculate the DM using NEGF ˆ 1 dE Gˆ ( E ) 2i Smeagol Smeagol No Compute electron density Yes Self-consistent? Output quantities Transmission and current Calculation of the transmission and current The transmission is calculated at the end of the self-consistent cycle It is possible to simplify it due to the small size of the G matrices Rˆ ˆR ˆ ˆ T ( E ) Tr [GLG GR G ]( E ) Rˆ ˆR ˆ ˆ Tr [G G G G ] L LR R RL The current is calculated by integrating the transmission e I dE T ( E )[ f ( E m L ) f ( E m R )]( E ) h Single level coupled to wide band leads Wide band leads have a constant density of states at the Fermi level e1 After coupling the level the onsite energy (e0) is renormalized by the real part of the self-energy (e1) The imaginary part of the self energy give the inverse of the lifetime (width of the Breit-Wigner resonance) ˆ RL/R i G 2 2 e1 e 0 2 1 G G R (E) T (E) E e 1 iG ( E e1 ) 2 G 2 3) How to run the code Calculation of the leads First run before calculating the transport properties of the extended molecule - Include two variables in the input file: BulkTransport T . It specifies wehter or nor the bulk parameters are written BulkLead LR . Left (L), right (R) or both (LR) leads * At the end of the calculation three or four files are generated: - bulklft.DAT and bulkrgt.DAT: contain the label of the system and basic information such as the Fermi energy, temperature, etc. - SystemLabel.HST: contains the Hamiltonian and overlap matrices - SystemLabel.DM: contains the density matrix Example of bulk file SystemName Au SystemLabel Au NumberOfAtoms 2 NumberOfSpecies 1 %block ChemicalSpeciesLabel 1 79 Au %endblock ChemicalSpeciesLabel %block PAO.Basis Au 1 n=6 0 1 5.0 %endblock PAO.Basis %block Ps.lmax Au 1 %endblock Ps.lmax LatticeConstant 1.00 Ang %block LatticeVectors 10.00 0.00 0.00 0.00 10.00 0.00 0.00 0.00 4.08 %endblock LatticeVectors AtomicCoordinatesFormat Ang %block AtomicCoordinatesAndAtomicSpecies 0.00 0.00 0.00 1 Au 1 0.00 0.00 2.04 1 Au 2 %endblock AtomicCoordinatesAndAtomicSpecies %block kgrid_Monkhorst_Pack 1 0 0 0.0 0 1 0 0.0 0 0 100 0.0 %endblock kgrid_Monkhorst_Pack xc.functional GGA xc.authors PBE MeshCutoff 200. Ry MaxSCFIterations 10000 DM.MixingWeight 0.1 DM.NumberPulay 8 DM.Tolerance 1.d-4 SolutionMethod diagon ElectronicTemperature 150 K SaveElectrostaticPotential T BandLinesScale pi/a %block BandLines 1 0.00 0.00 0.00 200 0.00 0.00 0.50 %endblock Bandlines BulkTransport T BulkLead LR DM.UseSaveDM T Calculation of the extended molecule. General Run the extended molecule file in the same directory that contains the bulk files - Variable to define the transport calculation: EMTransport T . Performs the transport calculation - Variables related to the energy contour: NEnergReal 500 . Number of points along the real axis (out of equilibrium) NEnergImCircle 90 . Number of points in the imaginary circle NEnergImLine 30 . Number of points in the imaginary line Npoles 10 . Number of poles of the Fermi function Delta 1.0d-4 . Small imaginary part of the Green’s Function EnergyLowestBound -10.0 Ry . Energy of the lowest bound of the EC Nslices 1 . Number of slices that are substituted by the bulk H and S Calculation of the extended molecule. Out of equilibrium - Variables related to the bias voltage (out of equilibrium): VInitial 0.0 . Initial value of the bias voltage VFinal 2.0 . Final value of the bias voltage NIVPoints AtomLeftVcte AtomRightVcte 10 . Number of points where the bias is going to be applied 9 . Left position where the bias ramp starts 36 . Right position where the bias ramp ends At each bias point the electronic structure is converged. Optionally, it is also possible to relax the atomic coordinates When the electronic structure converges the transmission and current at that bias point are calculated Calculation of the extended molecule. Transmission and VH - Variables related to the calculation of the transmission NTransmPoints 800 . Number of energy points where the transmission is going to be calculated InitTransmRange -8.0 eV . Initial value of the transmission range FinalTransmRange 2.0 eV . Final value of the transmission range - Variables related to the energy mismatch between bulk and EM HartreeLeadsLeft HartreeLeadsRight 2.040 Ang . Left position where the correction is applied 10.200 Ang . Right position where the correction is applied HartreeLeadsBottom -1.711 eV . Value of the Hartree potential at a certain point in the leads Example of extended molecule file SystemName Au.em SystemLabel Au.em NumberOfAtoms 6 NumberOfSpecies 1 %block ChemicalSpeciesLabel 1 79 Au %endblock ChemicalSpeciesLabel %block PAO.Basis Au 1 n=6 0 1 5.0 %endblock PAO.Basis %block Ps.lmax Au 1 %endblock Ps.lmax LatticeConstant 1.00 Ang %block LatticeVectors 10.00 0.00 0.00 0.00 10.00 0.00 0.00 0.00 12.24 %endblock LatticeVectors AtomicCoordinatesFormat Ang %block AtomicCoordinatesAndAtomicSpecies 0.00 0.00 0.00 1 Au 1 0.00 0.00 2.04 1 Au 2 0.00 0.00 4.08 1 Au 3 0.00 0.00 6.12 1 Au 4 0.00 0.00 8.16 1 Au 5 0.00 0.00 10.20 1 Au 6 %endblock AtomicCoordinatesAndAtomicSpecies ... EMTransport T NEnergReal 500 NEnergImCircle 60 NEnergImLine 30 NPoles 10 Delta 2.d-4 EnergLowestBound -8.d0 Ry NSlices 1 VInitial 0.d0 eV VFinal 2.d0 eV NIVPoints 10 AtomLeftVCte 2 AtomRightVCte 7 TrCoefficients T NTransmPoints 800 InitTransmRange -14.0 eV FinalTransmRange 8.0 eV PeriodicTransp T UseLeadsGF F HartreeLeadsLeft 2.040 Ang HartreeLeadsRight 10.200 Ang HartreeLeadsBottom -1.711 eV #%block SaveBiasSteps #012 #%endblock SaveBiasSteps DM.UseSaveDM T 4) Simple examples Infinite atomic chain. Equilibrium Band structure and transmission Two atoms in the unit cell → two crossing bands in the Brillouin zone One single channel at every energy Perfect squared transmission Infinite atomic chain. Out of equilibrium Out of equiilibrium transmission Starts to disappear at the edges, where the bands of both leads mismatch Current Ohmic regime Diatomic molecule. Equilibrium I Effect of separating the atoms from the leads Changing the coupling configuration Diatomic molecule. Equilibrium II Effect of changing the intramolecular distance Decrease the distance Increase the distance Diatomic molecule. Out of equilibrium Bonding Antibonding Different behaviour of the bonding and antibonding orbitals Bias-dependent transmission Negative differential resistance (NDR) 4) Some calculations Magnetoresistance effects in nickel chains Properties of this system: - Highest magnetic moments in the middle of the chain - The spins invert close to the leads Spin-polarized Non-collinear - Abrupt change (collinear) of the magnetization Symmetric, parallel Symmetric, antiparallel Aymmetric, antiparallel Metallocenes inside carbon nanotubes Chains of metallocenes inside CNTs Reduction of the magnetoresistance due to charge transfer from the metallocene CoCp CoCp + CNT CNT Magnetoresistance effect Even-odd effect in monoatomic chains Different conductance depending on the number of atoms in the chain The oscillations depend on the type of contact configuration Gold Sodium Conductance of H2 molecules couples to Pt or Pd leads Two possible configurations: parallel or perpendicular to the transport direction In case of Pd the H2 molecule can go inside the bulk and contacts Continuous line: Pt Dashed line: Pd Oscillations in Pt chains Chains with a number of atoms between 1 and 5 Structural oscillations due to changes in the levels at the Fermi level as the chain is stretched from a zigzag to a linear configuration 2 atoms 3 atoms I-V calculations of polyynes between gold leads Two types of molecules with different coupling atoms S Two types of NDR due to a different evolution of the resonances with bias N (a) (b) Porphyrins between gold leads Very large calculations with more than 500 atoms Evolution of the conductance with the number of porphyrin units and the angle between them Fin
© Copyright 2024