Simulating molecular conductance using real-time density functional theory

We present real-time density functional calculations of finite-bias conductance in a polyacetylene molecular wire. Our approach is based on a novel, efficient method for numerically propagating the time-dependent Kohn-Sham equations in a Gaussian basis. Localized density constraints are used to create an appropriate chemical potential bias that, when released, causes charges to flow from one end of the molecule to the other, generating a current. Our numerical scheme is efficient enough that one is able to perform “brute force” conductance calculations by simply increasing the size of the electron reservoirs and propagating until a reasonable average current can be extracted. We demonstrate the feasibility of this approach on a simple polyacetylene wire. By varying the size of the finite leads and comparing to commonly used nonequilibrium Green’s function calculations, we show that reliable current-voltage curves can be obtained from a finite length of the molecular wire, even though the system never reaches a steady state. Our results indicate that it should be technically feasible to perform the same type of “brute force” simulations on molecular junctions, although it seems unlikely that a true steady state will ever be reached in these cases, due to the greater significance of current fluctuations at low transmittance.

Figure 1

(Color online) Predictor-corrector routine for the second-order Magnus integrator. The order row shows the time order (in $dt$) to which the matrices in the same column are correct to.

Figure 2

(Color online) Minimum wall time required to obtain a prescribed average absolute error in the final density matrix of methane ($\mathrm{B}3\mathrm{LYP}∕6–31{\mathrm{G}}^{*}$, $120\mathrm{a.u}$ of propagation) using various approximate propagators: second order Magnus (green dashed), fourth order Runge-Kutta (teal dot-dashed with squares), fourth order Magnus (red solid), and sixth order Magnus (blue dotted).

Figure 3

(Color online) Schematic of the source-wire-drain geometry used in the present simulations. The bias is applied to the left and right groups of atoms, which act as a source and drain for electrons, respectively. For different wire lengths (e.g., 50 carbons versus 100) the wire length is kept fixed and the size of the source and sink are varied.

Figure 4

(Color online) Initial density corresponding to a chemical potential bias in polyacetylene. Red indicates charge accumulation and green charge depletion relative to the unbiased ground state. At time $t=0$ the bias is removed and current flows from left to right.

Figure 5

(Color online) Transient current through the central four carbons in ${\mathrm{C}}_{50}{\mathrm{H}}_{52}$ at a series of different chemical potential biases. There is an increase in current as voltage is increased, along with large, persistent fluctuations in the current. The currents are converged with respect to time step and the apparent noise is a result of physical fluctuations in particle flow through the wire.

Figure 6

(Color online) Transient current through the central four carbons in ${\mathrm{C}}_{50}{\mathrm{H}}_{52}$ at a series of different chemical potential biases smoothed over a time window of width $\Delta t=0.36\mathrm{fs}$. The average currents are now move clearly visible. The slow decay of the current at later times results from the partial equilibration of the finite left and right leads.

Figure 7

(Color online) Maximum smoothed current through the central four carbons in ${\mathrm{C}}_{50}{\mathrm{H}}_{52}$ as a function of chemical potential bias (red pluses). For comparison, we also present the analogous result for the central carbons in ${\mathrm{C}}_{100}{\mathrm{H}}_{102}$ (green squares) demonstrating convergence of the calculation with respect to lead size. The blue line is a linear fit to the ${\mathrm{C}}_{50}{\mathrm{H}}_{52}$ data at low bias indicating that polyacetylene is an Ohmic resistor with a conductance of $\approx 0.8G_0$

Figure 8

(Color online) Maximum smoothed current through the central four carbons in ${\mathrm{C}}_{50}{\mathrm{H}}_{52}$ as a function of chemical potential bias using real-time TDDFT (red line) and an NEGF approach described in the text (blue line). The two calculations are nearly identical at low bias and differ somewhat at higher biases due to the lack of self-consistency in the NEGF results.

Figure 9

(Color online) Maximum smoothed current through the central four carbons in ${\mathrm{C}}_{50}{\mathrm{H}}_{52}$ as a function of chemical potential bias (red line) and voltage bias (green). The results are quite similar until $4\mathrm{V}$, at which point the bias is so large that the finite width of the valence band for polyacetylene causes the conductance to plateau in the voltage biased case.

Figure 10

(Color online) Transient current through the central four carbons in ${\mathrm{C}}_{100}{\mathrm{H}}_{102}$ at a series of different chemical potential biases. The current fluctuations previously observed with smaller reservoirs in ${\mathrm{C}}_{50}{\mathrm{H}}_{52}$ persist and are therefore not associated with a finite size effect.

Figure 11

(Color online) Statistical noise in the current through the central four carbons in ${\mathrm{C}}_{100}{\mathrm{H}}_{102}$. The data (squares) can be fit to a sub-Poissonian distribution (line).