Spin transport in the phenalenyl radical throughout molecular dynamics trajectories

Abstract

The phenalenyl family of radicals has previously been identified as either possessing high spin filter efficiency (SFE) or excellent conductance, with these behaviors being tunable by the application of bias voltage. However, the persistence of these traits in the presence of thermally-induced structural fluctuations has yet to be determined. In this study, we employ a combination of ab initio molecular dynamics (AIMD) and the non-equilibrium Green's function technique combined with density functional theory (NEGF-DFT) to show how spin and charge transport behaviors evolve in a phenalenyl (PLY) molecular junction under thermal motion. To this end, we performed 5000 fs AIMD trajectories on Au-PLY-Au molecular junctions at temperatures of 300 K, 500 K, and 700 K, calculating charge and spin transport properties every 100 fs over these simulations. Remarkably, the molecule remained bonded to the Au electrodes at temperatures up to 700 K. Similarly, the SFE of the junction was stable against fluctuations caused by thermal-induced motion, with the SFE remaining above 80% for the majority of the three trajectories. Although there were a few instances of SFE dropping below 80%, strong SFE was quickly recovered within 200 fs. Additionally, the conductance of the molecular junction broadly increased with the temperature of the AIMD trajectory, with the 700 K trajectory producing the highest set of conductance values. Through the combination of AIMD and NEGF-DFT techniques, we show the robustness of spin filtering behavior in the PLY against thermally-induced structural changes in molecular junctions and that higher temperature can result in higher conducting configurations.

---

1. Introduction

Typically, spintronic devices are constructed out of sheets of thin magnetic and non-magnetic materials, through which architectures such as a spin-transfer torque structure can be created.^1–3 However, these alternating magnetic and non-magnetic structures require careful processing and painstaking engineering. Organic molecules, in contrast, are easy to manufacture, possess a degree of purity unmatched in the inorganic world, have finely tunable electronic structures which are altered via well-known chemical reactions, and can be prepared in a wide variety of types and combinations.⁴ Most importantly, since organic molecules are made from elements of the upper rows of the periodic table (i.e., C, H, S, O), the spin–orbit and hyperfine interactions – the primary methods through which spins interact with the environment – are weak, leading to long lifetimes for electron spins.^4,5 Evidently, organic molecules provide an excellent testbed through which the interactions between spin and current can be further explored.

Organic spin valves achieve significant spin polarization through the use of ferromagnetic electrodes and typically operate in low-temperature conditions.^4,6 The spin polarization is determined by the overall combination of the organic molecule and the electrode, where the spin current originates from the spin-polarized density of states in the electrodes. Naturally, the question of imparting a spin-polarized density of states in the junction through means other than ferromagnetic electrodes arises. Radicals, a class of organic molecules containing one or more unpaired electrons, inherently possess a spin-split density of states which may translate to in-junction behavior and impart a spin-split density of states to the overall molecular junction.^4,6–10 As such, radicals are advantageous in that they do not require ferromagnetic electrodes, external magnetic fields, or transition metal elements to achieve significant spin polarization. For example, the phenalenyl class of radicals has previously been studied for use as a spin filter; it was revealed that phenalenyl (PLY-2SH) exhibits a high degree of spin filtering that can be turned on or off by application of bias voltage.⁶ However, that study was performed on a single, geometrically optimized (minimum-energy) structure. While focusing on a minimum-energy junction provides a representative junction geometry, and is hence a good approximation, this technique discounts the role thermal-induced motion could play on the molecular junction.

At room temperature, thermal fluctuations change the conformation of the molecular junction and induce variations in the geometry (compared with the optimized structure) that could change charge and spin transport characteristics as the junction evolves in time. Therefore, an ensemble of structures representing the various configurations visited by the junction must be generated to properly evaluate transport behavior in thermal environments.^11,12 This study investigates the room-temperature junction behavior of the PLY radical, and to that end we performed a series of brute-force calculations to determine the robustness of spin filtering in the PLY radical with –SH linker groups against thermal-induced fluctuations. We combined ab initio molecular dynamics (AIMD) trajectories with the non-equilibrium Green's function technique as implemented in density functional theory (NEGF-DFT) to calculate transport behavior over an ensemble of structures created by thermal-induced motion in the junction. Notably the molecule remained attached to the electrode tips at temperatures of 300 K, 500 K, and 700 K, and the spin filtering efficiency remained high throughout all trajectories. Additionally, the distribution of conductance values over the AIMD trajectories trended to greater conductance values at higher temperatures, with the highest conductance configurations visited in the 700 K trajectory.

2. Computational methods

The nature of this study necessitated the use of different methods to compute the electronic and physical structure of the PLY molecule, generate the AIMD trajectories of the molecule in-junction, and perform transport calculations. Most notably, the B3LYP exchange–correlation functional was utilized in the ORCA code to evaluate the isolated PLY molecule, while the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional was used for AIMD and transport calculations performed in VASP and NanoDCAL, respectively. This approach has been extensively used in other studies,^6,9,13,14 although subsequent calculations utilizing PBE were performed on isolated PLY to mitigate the discrepancy between the different methods.

2.1. Density functional theory

Calculations to determine the physical and electronic structure of PLY-2SH were performed using the ORCA code.¹⁵ The B3LYP hybrid functional^16–19 was used for the exchange–correlation functional with the 6-311++G(d,p) basis set.²⁰ The electronic convergence criterion was set to 1 × 10⁻⁵ eV in terms of the total energy of the molecule for all ORCA calculations. After the molecular structure was calculated, the hydrogen atoms were removed from the –SH groups, as is expected when the molecule is placed in a junction, and the PLY-2S molecule was then placed in between two Au (111) electrodes with lattice constants of 2.94 Å determined by VASP relaxations as shown in Fig. 1.

Fig. 1 A schematic showing the two-probe electrode-molecule-electrode system. Gold-colored atoms represent Au, yellow atoms represent S, brown atoms represent C, and white atoms represent H. Atoms shaded in red are held frozen and all other atoms are free to move. Inset: Chemical structure of the phenalenyl radical.

Structural relaxation and ab initio molecular dynamics (AIMD) calculations were performed using the Vienna ab initio Simulation Package (VASP)²¹ and employed a plane-wave basis set, with ion-electron interactions approximated by projector augmented wave potentials^21–25 and the exchange–correlation functional provided by the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation.^26,27 Van der Waals forces were considered using the DFT-D3 method with Becke–Johnson damping, which is an empirical correction to the interatomic forces.²⁸

Structural relaxation calculations were performed on the two-probe structure shown in Fig. 1 to identify the optimal inter-electrode separation and binding configuration. Structural relaxations employed a mixture of the blocked-Davidson and RMM-DIIS algorithms and utilized the Γ k-point^29,30 with a 400-eV plane-wave energy cutoff. The relaxed two-probe structure was then used as the initial structure for AIMD trajectories. Each AIMD trajectory employed the Γ k-point, a 1-fs timestep, and a 400-eV plane-wave energy cutoff.^29,30 The electronic convergence criterion was set to 1 × 10⁻⁶ eV in terms of the total energy of the two-probe system for both AIMD and structural relaxation calculations. AIMD trajectories were performed for 5000 fs at temperatures of 300 K, 500 K, and 700 K which were fixed via a Nosé–Hoover thermostat.³¹

2.2. Non-equilibrium green's function-density functional theory

Structural data generated by the AIMD trajectories were extracted and used to calculate electron transport throughout the trajectories, at a sampling rate of once per 100 fs due to the high computational cost of electron transport calculations. This approach has been utilized previously by our group to combine AIMD and NEGF-DFT calculations.³² Over the course of the AIMD trajectories, structures with greater energetic favorability will be visited more often. This sampling at a set time interval mimics the stochastic nature of break-junction experiments.^33,34 Sampling more data points from the AIMD trajectories would provide a clearer picture of the processes occurring in the system but is prohibitively computationally expensive.

As the outermost layers of the two-probe system in Fig. 1 were held frozen, the two-probe system in Fig. 1 could be joined to Au (111) electrodes, which repeat to infinity on the left and right sides of what is now called the central cell.³⁵ Electron transport properties were calculated via the non-equilibrium Green's function technique in conjunction with density functional theory, as implemented in the NanoDCAL code.^35,36 These calculations employed the Γ k-point for the central cell, the PBE exchange–correlation functional, Troullier–Martins norm-conserving pseudopotentials³⁷ and a 50-Hartree real space grid cutoff energy, consistent with previous studies.⁶ We note that, while the PBE-D3 approach was employed in the AIMD trajectories to obtain reliable geometries of the molecule in-junction, van der Waals corrections are not suitable for obtaining the electronic structure used in transport calculations. The electronic convergence criterion was set to 2.7 × 10⁻³ eV in terms of the Hamiltonian and density matrices of the two-probe system. A double zeta polarized (DZP) basis set was used for C, H, and S atoms, while a single zeta polarized (SZP) basis set was utilized for Au. This approach was used to significantly save on computational cost and is consistent with previous studies.^6,38,39 Additionally, these ballistic transport calculations were performed at a temperature of 0 K. This means that transport behavior is dictated solely by structural changes induced via thermal fluctuations, as opposed to, e.g., electron–phonon scattering. We performed these calculations at 0 K due to the short length of the molecule and the results from a previous study by Topolnicki et al.¹² The short length of the molecule (0.7 nm) implies a lack of electron–phonon coupling,⁴⁰ and the previous study by Topolnicki et al. showed that including the effects of temperature via Fermi level smearing only produced minute static shifts in the transmission spectra, far outweighed by the structural changes induced via MD methods.

Calculation of transport properties starts from the calculation of the Green's function of the central cell as a function of energy E,

G(E) = [(E + iη)S − H − Σ₁ − Σ₂]⁻¹. (1)

In eqn (1), H and S are the Hamiltonian and overlap matrices for the scattering region, respectively, η is a positive infinitesimal, and Σ_1,2 are the self-energies of the left and right electrodes.^35,41 The self-energies describe the effects of coupling the electrodes to the scattering region and contain real and imaginary components. The real component represents a shift of the energy levels in the scattering region due to charge transfer from the electrodes, and the imaginary component represents the broadening of these levels corresponding to the broadening matrix seen in eqn (2),

The self-energy is calculated by an iterative technique.⁴² The electronic density matrix is then obtained as,

where μ_1,2 are the electrochemical potentials of the left and right electrodes and f(E,μ) is the Fermi–Dirac function which describes the occupation of an energy level given the energy and electrochemical potential, μ. With the Green's function obtained from the self-consistent, iterative process calculated using eqn (1)–(3),^35,36,42 we can now obtain the transmission function T(E) viaeqn (4), which gives the probability that an electron with a given energy E will transmit from one electrode through the scattering region into the other electrode as,^{35,36,41,43,44}

T(E) = Tr(Γ₁GΓ₂G^†). (4)

The electric current can then be obtained by integrating the transmission function over all energy points from the electrochemical potential of the left electrode to the electrochemical potential of the right electrode, and is given by:⁴⁵

which can be approximated at T = 0 K as,

For a given bias voltage, spin filter efficiency (SFE) is calculated as,

where I_α and I_β represent current through spin up and spin down channels (assigned arbitrarily), and G_α and G_β are the conductances of the α and β channels which can be recovered by dividing the calculated current by the bias voltage. All transport calculations were performed at a user-defined bias voltage of V_b = 0 V. In the limit of zero bias, the conductance across the simulation cell can be approximated as the value of the transmission function at the Fermi energy.^41,43,44 However, taking the SFE and conductance from the Fermi energy alone can make these results overly sensitive to slight shifts in T(E). To avoid this oversensitivity, we integrated over the transmission spectra in the energy window of [−0.1,0.1] eV, where E_F = 0, and divided by 0.2 to recover conductance. This ensures that fluctuations in T(E) near the Fermi level do not have an outsized effect on conductance and SFE. Accordingly, SFE and conductance are calculated as,andrespectively, where T_α(E) and T_β(E) are the α and β components of the transmission spectrum.

3. Results and discussion

The electronic and physical structures of PLY-2SH are shown in Fig. 2 and were obtained via ORCA code calculations. Using B3LYP as the exchange–correlation functional, the highest-occupied molecular orbital (HOMO) and lowest-unoccupied molecular orbital (LUMO) exhibited minimal splitting between their α and β spin components, with the HOMO orbital splitting by 0.16 eV and the LUMO orbital splitting by 0.05 eV.

Fig. 2 (a) Energies of the highest occupied (HOMO), singularly occupied (SOMO), singularly unoccupied (SUMO), and lowest unoccupied (LUMO) molecular orbitals in the PLY-2SH radical, calculated with B3LYP (left) and PBE (right). (b) Eigenfunctions of the orbitals plotted in (a).

The singularly occupied and unoccupied molecular orbitals (SOMO and SUMO), grouped by their shared spatial function, are split by 1.87 eV.⁴⁶ When the PBE exchange–correlation functional was utilized the HOMO and LUMO orbitals were split by 0.11 eV and 0.07 eV, respectively, while the SOMO and SUMO orbitals were split by 0.61 eV. While the exact value of the split between the SOMO and SUMO orbitals differs significantly between the B3LYP and PBE functionals, the orbital splitting values for the HOMO and LUMO levels are similar between the two functionals. Additionally, utilizing the PBE functional gives a SOMO–SUMO split value remarkably similar to that obtained via transport calculations performed in NanoDCAL, with NanoDCAL calculating a SOMO–SUMO split value of 0.52 eV, vide infra. We have also performed calculations on the isolated molecule using a DZP basis set, and these data are included in Fig. S19 of the SI. This result is in agreement with others reported in literature.⁶ The spin-split SOMO and SUMO indicate that spin-resolved transport is possible in the junction.

To that end, transport calculations were performed on the optimized metal–molecule–metal junction, and the resulting α and β spin transmission spectra are shown in Fig. 3(a). Clearly, the spin-split characteristic from the isolated molecule maintains itself in the junction as evidenced by the spin-split spectra in Fig. 3(a), which is in agreement with previous studies performed on PLY-2SH.⁶ The HOMO−1, HOMO, SOMO, and SUMO peaks were identified via scattering state calculations, and more detail on this process can be found in Section S2 of the SI. The β and α peaks at −0.05 eV were assigned to HOMO−1 and HOMO orbitals, respectively, and are indicated by the red and orange arrows in Fig. 3(a). The T(E) peak at −0.02 eV was assigned to the SOMO orbital and is indicated by the green arrow in Fig. 3(a). Finally, the peak at 0.5 eV was assigned to the SUMO orbital and is indicated by the blue arrow in Fig. 3(a). After the static transmission spectra were calculated, the Au-PLY-Au structure shown in Fig. 1 was employed for 5000 fs AIMD trajectories as described in Section 2.

Fig. 3 Transmission spectra calculated for (a) the static, geometrically optimized metal–molecule–metal junction and structures extracted every 100 fs from the (b) 300 K, (c) 500 K, and (d) 700 K AIMD trajectories. The association of HOMO−1, HOMO, SOMO, and SUMO orbitals with specific transmission resonances are indicated on the static transmission spectra by red, orange, green, and blue arrows, respectively.

Remarkably, PLY maintained its robust bond to the Au electrodes in each trajectory, never detaching from either electrode even at 700 K. Snapshots are plotted every 1000 fs from each AIMD trajectory in Section S1 of the SI, demonstrating the attachment of the PLY molecule to the electrode throughout each trajectory. From these AIMD trajectories, structures were extracted every 100 fs and employed for use in NEGF-DFT calculations. The resulting transmission spectra calculated from the 300 K, 500 K, and 700 K AIMD trajectories are plotted in Fig. 3(b–d). Clearly the thermal motion affects the transmission spectra, although different peaks in the spectra are affected to varying extents.

The HOMO and HOMO−1 peaks near −0.5 eV do react to thermal fluctuations, with their energetic positions exhibiting shifts of 0.1–0.5 eV compared to the optimized structure. The spread of the HOMO and HOMO−1 peaks is proportional to temperature, with far more variation in HOMO/HOMO−1 positions occurring in the 700 K trajectory than in the 500 K trajectory, which in turn exhibits more variation than the 300 K trajectory. Although the spectra exhibit large changes in shape, the energetic positions of the transmission peaks near the Fermi energy are relatively fixed. The SOMO peak near −0.02 eV is remarkably stable against thermal fluctuations in the 300 K trajectory, never deviating more than 0.1 eV away from the optimal structure peak position. This peak changes its energetic position by larger amounts as the temperature is increased, although it is much less varied than the HOMO/HOMO−1 or SUMO peaks. We suspect this is due to the extension of its tail across E_F, which implies that the SOMO level is pinned to the Fermi level.⁴⁷ The SUMO peak, near 0.5 eV, is more varied, with fluctuations changing the peak position by 0.2 eV in some cases. Again, the variation in peak position appears to be proportional to temperature, as high temperatures have larger spreads in SUMO energetic positions. Remarkably, while thermal-induced fluctuations in the junction do shift the energetic position of orbital resonances in the transmission spectra, spin degeneracy is never restored throughout the trajectory as evidenced by the spin-split nature of the transmission spectra.

From the generated transmission spectra, conductance and SFE were calculated as described in Section 2. The conductance and SFE of the PLY junction are plotted over the 300 K, 500 K, and 700 K AIMD trajectories in Fig. 4.

Fig. 4 Conductance and spin filter efficiency calculated every 100 fs over (a) the 300 K, (b) the 500 K, and (c) the 700 K AIMD trajectories. Spin filter efficiency is shown in blue, while conductance is shown in orange. Shaded purple circles denote SFE dips caused by a reduction in the SOMO peak.

Remarkably, the PLY radical maintains a strong spin filtering behavior in all three trajectories. The SFE is calculated to remain above 80% in the majority of the structures generated via the 300 K, 500 K, and 700 K trajectories. Dips occur in the SFE, notably at 200 and 2300 fs in the 300 K trajectory, at 1900 and 4100 fs in the 500 K trajectory, and at 600, 1400, 2100, and 4100 fs in the 700 K trajectory. Transmission spectra from these SFE dips are included in Section S3 of the SI. Interestingly, the SFE dips at 200 and 2300 fs in the 300 K trajectory, as well as those at 1900 fs and 4200 fs in the 500 K trajectory and 1400 and 4100 fs in the 700 K trajectory, are caused by a dramatic reduction in the SOMO peak compared to the preceding and succeeding transmission spectra. The dips caused by this feature are highlighted by shaded purple circles in Fig. 4.

The reduction in SOMO peak height is not attributable to a single physical cause, as analyses of structural data around the SFE dips in the 300 K trajectory and charge data over the first 3300 fs of the 300 K AIMD trajectory revealed no correlation between structural factors, charge, and SFE behavior (see Section S4 in the SI and discussion below). In some cases, the β HOMO resonance of the lowest SFE structure is closer to the Fermi energy (700 K, 1400 fs), while in others the β HOMO resonance is positioned further from the Fermi energy in comparison to preceding and succeeding spectra (700 K, 4100 fs). Still in other cases, the α transmission spectrum at the dip is an order of magnitude higher than the surrounding α spectra, in addition to a reduction in the height of the SOMO peak at the dip (500 K, 1900 fs).

The conductance, like the SFE, is more varied with increasing temperature as shown in Fig. 4. Dips and peaks in the conductance are short-lived, with the system returning to intermediate conductance values (0.05–0.15 G₀) within 200 fs. Notably, the changes in conductance do not occur at the same time as those in the SFE, and crucially, there is no obvious correlation or anti-correlation between conductance and SFE, as seen in Fig. 4 and Fig. S16 of the SI. This complements a previous study performed on a series of PLY derivatives,⁶ in which the PLY derivatives could either perform as efficient spin filters or high conductors. That anti-correlation between SFE and conductance was identified across a series of PLY derivatives, whereas this finding demonstrates a lack of anti-correlation between conductance and SFE in the same molecule moving in a junction.

Additionally, the amount of charge present on the molecule does not change significantly during these SFE dips, indicated by a lack of shift in the SOMO position, and no singular structural factor (i.e., dihedral angle) was identified to dominantly determine the conductance and SFE behavior. Regression analysis was performed to demonstrate correlation between dihedral angle and SFE or conductance around the 200 fs and 2300 fs SFE dips in the 300 K trajectory and is presented in Fig. S17 of the SI. Specifically, dihedral angle values were quantified by the sum of dihedral sin² products, which gives the alignment between the π-system of the molecule and the Au–S bonds connecting the molecule to the electrodes.^48,49 While there is some correlation between dihedral angles and conductance, there is no correlation between dihedral angles and SFE. Regression analysis was also performed between the amount of excess charge on the molecule and the SFE and conductance over the first 3300 fs of the 300 K trajectory and is presented in Fig. S18 of the SI. An excess charge of ca. −0.45e sits on the molecule in-junction, with minor variations throughout the trajectory, yet no correlation exists between charge and either SFE or conductance. Accordingly, changes in conductance and SFE in the junction are attributed to a multitude of structural factors such as bond distances, bond angles, and overall tip-molecule geometry. Importantly, the SFE is never destroyed in any of the trajectories. Momentary dips in the SFE are followed by a return to high filter efficiency behavior within 200 fs in all cases.

In simulations with higher temperatures, the SFE and conductance values become more disparate, as shown in Fig. 5. In all three trajectories, the conductance remains largely around 0.1 G₀, although the 700 K trajectory shows some higher conducting outliers.

Fig. 5 Spin filter efficiency (left) and conductance (right) binned into histograms of width 10% and 0.05 G₀ respectively, for the 300 K, 500 K, and 700 K trajectories.

The SFE also remains close to 100% for all three data sets, although the 500 K and 700 K data especially show low-filtering outliers. As expected, the number of structures which have either a conductance around 0.10 G₀ or an SFE near 100% decreases as the temperature of the simulation is increased. Notably, while higher temperature AIMD trajectories cause higher conductance values to appear, it has precisely the opposite effect on spin filtering behavior, causing lower SFE values to appear. This is in part due to the nearly 100% filtering behavior present in the PLY molecule at low temperatures. When disparate values appear, they can only decrease. The increase in conductance at higher temperatures is caused by motion of the SUMO towards the Fermi energy in some higher temperature structures, which also explains the appearance of lower SFE values in higher temperature trajectories. These results are in line with the conclusions drawn in a previous molecular electronics study in which minimal energy structures did not necessarily correspond to those with the highest conductance.¹³ Thermal motion allows the molecule to visit less stable structures that may conduct better than minimal energy structures. Additionally, this behavior is in line with conclusions drawn from studies on molecular break-junction histograms as well as previous MD studies on single-molecule junctions.^12,33

The combination of AIMD and NEGF-DFT methods allows for the construction of histograms like those shown in Fig. 5, providing a more complete connection to the actual results of molecular electronic experiments.⁴³ AIMD, through inclusion of thermal motion, generates a series of structures which have a certain amount of randomness due to the statistical nature of AIMD calculations.³¹ This is in-line with break-junction experiments performed on single molecules, which are highly stochastic in nature due to wide variety of bonding sites, tip geometry differences, and thermal energy present in the environment.^33,34 However, break-junction experiments have proven extremely difficult to perform on radicals due to their highly reactive nature,¹⁰ and we are unaware of any such published work to which we can directly compare our generated histograms. Instead, we hope our work will inspire future experiments and that our findings will inform future computation and experimental efforts and guide the design of spin filter molecules.

4. Conclusions

In summary, we have employed a brute-force method to evaluate the robustness of spin filtering efficiency in the PLY radical against thermal-induced structural fluctuations. We calculated T(E) from three sets of 50 unique structures obtained over AIMD trajectories performed at 300 K, 500 K, and 700 K, and from T(E) we obtained the conductance and spin filtering efficiency across the junction. Notably, conductance and spin filtering efficiency values were increasingly disparate with temperature, although high spin filtering efficiency was maintained throughout the majority of all trajectories. Importantly, no obvious correlation or anti-correlation was identified between conductance and spin filtering efficiency. These results demonstrate that SFE can indeed be robust against thermal-induced fluctuations, with the large dips in SFE establishing the variability in electronic behavior that can occur in a molecular junction, and highlight the need to consider environmental factors when calculating the in-junction behavior of a single molecule.

Author contributions

K. B. - conceptualization, data curation, formal analysis, writing: original draft, writing: review and editing. M. S. - conceptualization, writing: review and editing, supervision, project administration, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Additional data is available upon reasonable request from the author. Supplementary information: structural snapshots taken every 1000 fs from AIMD trajectories, scattering state calculations to perform peak assignments, transmission spectra corresponding to SFE dips in AIMD trajectories, different data analyses performed on generated SFE and conductance data, and calculations comparing DZP and 6-311++G(d,p) basis sets for isolated PLY molecule calculations. See DOI: https://doi.org/10.1039/d5cp03312h.

Acknowledgements

K. B. and M. S. were supported by funds from Binghamton University and the National Science Foundation under Grant No. DMR-2225370. NEGF-DFT calculations were performed on the Spiedie cluster at Binghamton University. AIMD calculations were performed on the Expanse cluster, a part of the Extreme Science and Engineering Discovery Environment (XSEDE/ACCESS), which is supported by the National Science Foundation Grant No. 1928224, under allocation TG-DMR18009.

References

1. D. Waldron, P. Haney, B. Larade, A. MacDonald and H. Guo, Phys. Rev. Lett., 2006, 96, 166804.
2. D. Waldron, V. Timoshevskii, Y. Hu, K. Xia and H. Guo, Phys. Rev. Lett., 2006, 97, 226802.
3. A. Fert, Rev. Mod. Phys., 2008, 80, 1517–1530.
4. S. Sanvito, Chem. Soc. Rev., 2011, 40, 3336–3355.
5. Z. Yu, F. Ding and H. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 205446.
6. M. Smeu, O. Monti and D. McGrath, Phys. Chem. Chem. Phys., 2021, 23, 24106–24110.
7. S. Naghibi, S. Sangtarash, V. Kumar, J. Wu, M. Judd, X. Qiao, E. Gorenskaia, S. Higgins, N. Cox, R. Nichols, H. Sadeghi, P. Low and A. Vezzoli, Angew. Chem., Int. Ed., 2022, 61, e202116985.
8. C. Herrmann, G. Solomon and M. Ratner, J. Am. Chem. Soc., 2010, 132, 3682–3684.
9. M. Smeu and G. DiLabio, J. Phys. Chem. C, 2010, 115, 17874–17879.
10. L. Li, C. Prindle, W. Shi, C. Nuckolls and L. Venkataraman, J. Am. Chem. Soc., 2023, 145, 18182–18204.
11. E. Montes, W. Rojas and H. Vazquez, J. Phys. Chem. C, 2025, 129, 9947–9953.
12. R. Topolnicki, R. Kucharczyk and W. Kaminski, J. Phys. Chem. C, 2021, 125, 19961–19968.
13. K. Batzinger, Q. Zhou, X. Ye, E. Borguet, S. Xiao and M. Smeu, Chem. – Eur. J., 2024, 30, e202400422.
14. P. Yasini, S. Shepard, T. Albrecht, M. Smeu and E. Borguet, J. Phys. Chem. C, 2020, 124, 9460–9469.
15. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 71–78.
16. A. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
17. C. Lee, W. Yang and R. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785.
18. S. Vosko, L. Wilk and M. Nusair, Can. J. Phys., 1980, 58, 1200–1211.
19. P. Stephens, F. Devlin, C. Chabalowski and M. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
20. R. Krishnan, J. Binkley, R. Seeger and J. Pople, J. Chem. Phys., 1980, 72, 650–654.
21. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
22. G. Kresse and J. Hafner, J. Phys.: Condens. Matter, 1994, 6, 8245–8257.
23. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
24. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
25. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
26. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
27. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
28. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
29. H. Monkhorst and J. Pack, Phys. Rev. B, 1976, 13, 5188.
30. J. Pack and H. Monkhorst, Phys. Rev. B, 1977, 16, 1748.
31. S. Nosé, Mol. Phys., 1984, 52, 255–268.
32. K. Batzinger, D. Liepinya and M. Smeu, Phys. Chem. Chem. Phys., 2024, 26, 5218–5225.
33. L. Mejía, P. Cossio and I. Franco, Nat. Commun., 2023, 14, 7646.
34. D. Dyer and O. Monti, J. Phys. Chem. Lett., 2023, 14, 10935–10942.
35. J. Taylor, H. Guo and J. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 245407.
36. J. Taylor, PhD thesis, McGill University, Montréal, Québec, 2000.
37. N. Troullier and J. Martins, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 1993–2006.
38. J. Ivie, N. Bamberger, K. Parida, S. Shepard, D. Dyer, A. Saraiva-Souza, R. Himmelhuber, D. McGrath, M. Smeu and O. Monti, ACS Appl. Mater. Interfaces, 2021, 13, 4267–4277.
39. P. Yasini, S. Shepard, M. Smeu and E. Borguet, J. Phys. Chem. B, 2023, 127, 9771.
40. X. Zhu, Surf. Sci. Rep., 2004, 56, 1–83.
41. S. Datta, Quantum Transport: Atom to Transistor, Cambridge University Press, Cambridge, 1st edn, 2005.
42. M. Lopez Sancho, J. Lopez Sancho and J. Rubio, J. Phys. F, 1984, 14, 1205–1215.
43. J. C. Cuevas and E. Scheer, Molecular Electronics: An Introduction to Theory and Experiment, World Scientific Publishing Co. Pte. Ltd, 2017.
44. S. Datta, Electronic transport in mesoscopic systems, Cambridge University Press, Cambridge, 1995.
45. M. Büttiker, Y. Imry, R. Landauer and S. Pinhas, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 6207.
46. S. Kasemthaveechok, L. Abella, J. Crassous, J. Autschbach and L. Favereau, Chem. Sci., 2022, 13, 9833–9847.
47. S. Rodrigues-Gonzalez, Z. Xie, O. Galangau, P. Selvanathan, L. Norel, C. Van Dyck, K. Costuas, C. Frisbie, S. Rigaut and J. Cornil, J. Phys. Chem. Lett., 2018, 9, 2394–2403.
48. N. Bamberger, D. Dyer, K. Parida, T. El-Assaad, D. Pursell, D. McGrath, M. Smeu and O. Monti, J. Phys. Chem. C, 2022, 126, 6653–6661.
49. D. Vonlanthen, A. Mishchenko, M. Elbing, M. Neuburger, T. Wandlowski and M. Mayor, Angew. Chem., Int. Ed., 2009, 48, 8886–8890.

---