{'fraction_non_alphanumeric': 0.026530405062434436, 'fraction_numerical': 0.011742411424317282, 'mean_word_length': 5.085667215815486, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
First-principles interpretation of electron transport through single-molecule junctions using molecular dynamics of electron attached states

The electron transport through the single-molecule junction of benzene-1,4-diamine (BDA) is modelled using ab initio quantum-classical molecular dynamics of electron attached states. Observations on the nature of the process are made by time-resolved analysis of energy differences, non-adiabatic transition probabilities and the spatial distribution of the excess electron. The role of molecular vibrations that facilitate the transport by being responsible for the periodic behaviour of these quantities is shown using normal mode analysis. The results support a mechanism involving the electron's direct hopping between the electrodes, without its presence on the molecule, with the prime importance of the bending vibrations that periodically alter the molecule–electrode interactions. No relevant differences are found between results provided by the ADC(2) and SOS-ADC(2) excited state models. Our approach provides an alternative insight into the role of nuclear motions in the electron transport process, one which is more expressive from the chemical perspective. GRAPHICAL ABSTRACT


Introduction
Molecular electron transport is the focus of scientific interest in different fields from material science to biochemistry. Research in this area was spectacularly boosted in the past two decades [1][2][3][4], in The theoretical modelling of these transport processes was also developed at the same time and now has an extensive literature [3,4,[10][11][12][13][14][15][16][17][18][19][20], allowing a reasonable modelling of such systems. The majority of the investigations on SMJs are based on the Landauer approach. The key quantity of this model is the transmission function, usually calculated from scattering theory using Green's function techniques. The corresponding model Hamiltonian is usually a simple one-electron representation of the molecular states, with parameters from density functional theory (DFT) calculations. [11,15,16,18,[21][22][23].
Extension of this theory by considering nuclear motion effects is also possible [14][15][16][17] but rarely done. One way to do this is via so called electron-phonon couplings, [11,14,16,19,24] usually described using a Hamiltonian of the typê whereĤ 0 is the (mean-field) electronic Hamiltonian of the molecule and the terminal region of the metallic contacts (these combined forming the so-called 'extended system') andV M−out is the interaction of this extended system with the bulk electrodes. The linear coupling of the electronic and nuclear degrees of freedom is described by the last term of Equation (1), withd i being electronic creation-annihilation operators andQ being the molecule's vibrational coordinate operators -usually also expressed in a second quantised formalism. The coupling elements with C being the transformation matrix between cartesian and normal mode coordinates, are very much reminiscent to the respective terms of the linear vibronic coupling (LVC) model used widely in chemistry [25][26][27][28][29] and can be regarded as a perturbation on the system's electronic structure by the molecule's vibrations. The above listed methods are relatively simple to apply, work acceptably in practice and their different forms are thus widely adopted by the physics community. However, from the quantum chemist's perspective this methodology is not quite satisfactory from the following points of view: (a) they do not provide a picturesque explanation of the mechanism of the electron transport (not distinguishing e.g. so-called inelastic tunnelling or hopping), and the role of the molecule's electronic structure is hidden as it plays the role of a 'scattering centre'; (b) for complex systems where the number of nuclear degrees of freedom included in Equation (1) is limited by computational complexity, the selection of considered normal modes is quite difficult without a more in-depth knowledge of their nature and impact; c) the model assumes a meanfield Hamiltonian of the system not allowing to check the importance of electron correlation or, in general, the dependence on the quantum chemical description of the electronic structure.
The aim of this proof-of-concept paper is to present an alternative description of the electron transport process, one which is more familiar from the chemical perspective. We employ a toolset proven to be effective for studying and interpreting time-resolved electronic-nuclear interactions: molecular dynamics. Our approach should shed light on the role of molecular motions and unveil the nature of structural distortions that facilitate the translocation of an attached electron, thereby providing a chemically more meaningful picture of electronic conduction. This model easily incorporates all normal mode vibrations of the investigated molecule and makes it possible to identify the most relevant ones. This methodology readily allows the use of high level electron correlation methods, i.e. it does not rely on the mean-field descriptions of the electronic structure that are quite widespread in the field despite their often controversial accuracy.

Model system
The physical picture behind this study is that an electron attached to the natural model system -imagined as having come from the bulk conductors -could move across the junction via finite-probability transitions between electron attached states localised on different sites of the system.
We chose to investigate the single-molecule junction of benzene-1,4-diamine (BDA), simulating the circumstances of the respective experiments. The conductance properties of this molecule have been extensively studied many times both experimentally [30][31][32] and theoretically [33][34][35]. The geometry of the BDA junction is rather well-defined due to the selectivity of the gold-amine bonding, resulting in the narrow spread of the measured conductance values. Another advantage of this system is that the energy of the highest occupied molecular orbital (HOMO) of BDA lies relatively close to the Fermi level of the electrodes, allowing the electron transport to be observed at a low bias.
Following the concept given in [34], our model system consists of a BDA molecule placed between two gold electrodes, modelled as 6-atomic clusters in a distorted octahedral configuration, as shown on Figure 1. The size of the chosen gold clusters is an important trade-off between physical accuracy and computational complexity: the number of tractable explicit gold atoms is limited by the scaling of the electronic structure method, and the number and density of electronic states -approaching a band structure in the asymptotic limit -also increases with the cluster size [36]. According to the ionisation properties of gold clusters of various sizes [36][37][38][39] the chosen six gold atoms provide a reasonable model of the terminal region of gold electrodes, approaching the IP and EA values of larger clusters to an extent which is satisfactory for the aim of this study.
The equilibrium structure of the BDA molecule was obtained at the MP2 level using the def2-SV(P) basis [40]. In the gold cluster the gold-gold bond lengths were set to the bulk value of 2.89 Å and the bond angles to 60 • , 90 • and 120 • , based on that given in [36]. The BDA was placed between the gold clusters so that the distance between the terminal gold and amino-nitrogen atoms is 2.38 Å, while the Au-N-C and Au-N-C-C angles are 112 • and 87.5 • , respectively. The structure of the model system is available in the Supplemental Material.
An advantage of the chosen model system is its relatively modest size, allowing us to employ more sophisticated electron correlation methods, surpassing the rather compromised reliability of various mean-field approaches found in the existing literature.

Ab initio methods
The electronic structure of the system was calculated using the Algebraic Diagrammatic Construction Second Order (ADC(2)) [41][42][43] method, with the core electrons excluded from the correlation treatment. The TURBOMOLE [44,45] program system was used for the calculations with the resolution-of-identity approximation, offering an attractive computational performance suited for the molecular dynamics calculations explained below. To also address the concerns raised by the rather controversial performance of ADC(2) in some cases [46], results with the spin-component-scaled variant SOS-ADC(2) were also obtained as comparison.
Recent studies by our group [47][48][49] revealed that the unsatisfactory performance of ADC(2) is remedied to a considerable extent by spin-component-scaling, while an agreement with the parent variant is often an indication of a good accuracy of both techniques. The double-ζ quality def2-SV(P) basis of Schäfer and coworkers [40] was used in the calculations with a def-ECP effective core potential (ECP) present on the gold atoms.

Continuum orbital strategy
The electron attached states were calculated using the continuum orbital strategy, an elegant idea first proposed by Stanton [50] which allows the use of excitation energy methodology to describe such states. This technique was successfully used later by us [49,51,52] and others [53,54]. In this approach the system is augmented by a continuum orbital occupied by two excess electrons and the electron attached states of the neutral system are obtained as excited states of this doubly charged anion, dominated by an excitation from the continuum orbital to one of the virtual MOs. Technically, the continuum orbital is included by adding a dummy atomic centre to the calculation, equipped with a single, gaussian type basis function of infinite size (in practice, an exponent value of 10 −8 was used) which, having no overlap with any other basis function, results in a standalone molecular orbital with an orbital energy of precisely zero. A great advantage of this idea is that any program code able to calculate excited electronic states can be made capable of providing ionised or electron attached states in a closed-shell formalism [51]. This not only has benefits on the performance side but also eliminates complications associated with open-shell wave functions. Another important aid of the continuum orbital technique is that all electron attached states are treated on an equal footing even with single-reference methods, due to the fact that they are all described with respect to a common reference -in contrast to, e.g. as excited states of a singly charged anionic ground state.

Analysis of the spatial extent of the wave function
The electron attached states can be conveniently analysed using the TheoDORE toolset of Plasser and co-workers [55][56][57]. This approach is based on a population analysis of the transition density integrated over the fragments that build up the system. In our case these are trivially chosen as the molecule (Fragment No. 1) as well as the left (Fragment No. 2) and right (Fragment No. 3) electrode, while the elements of the transition density essentially describe the weight of each virtual orbital in the electron attachment. The Natural Transition Orbitals (NTOs, see [55,56]) of the electron attached states and the respective diagonal elements of the socalled -matrix [56][57][58] thus give an expressive picture of the excess electron's spatial distribution. A quantity very useful for investigating electron attached states is the POS f descriptor [57] calculated from the -matrix, which tells on which pre-defined fragment the electron is attached. For a delocalised EA state the POS f takes a value between the involved fragments' numbers.  Table 1. ADC(2) isosurfaces with the contour value of 0.02 are shown.

Molecular dynamics
The impact of the molecule's nuclear motions was explored within a nonadiabatic molecular dynamics framework using the Newton-X software package of Barbatti et al. [59]. Relying on the concept of trajectory surface hopping (TSH), the excited state process is generally modelled by a large set of Newtonian nuclear trajectories, allowed to jump from one electronic state into another in case the associated 'hopping probability' is sufficiently large. For a detailed introduction to the concept, limitations and applications of the TSH approach, the reader is advised to check out the excellent reviews of Barbatti and co-workers [59][60][61].
To obtain a picture of the conduction process, the molecular dynamics is started in an electron attached state localised on the electrodes (called electrode states hereafter, see below), and the potential hopping events towards other states are investigated from the electron transport perspective. In order to mimic the solid-state nature of the electrodes, the atoms of the gold clusters were not allowed to move, i.e. the trajectory was exclusively driven by the BDA molecule's internal motions. A total of 70 trajectories were propagated for a duration of 100 fs with a nuclear integration step of 0.5 fs. The initial conditions were generated for the isolated molecule, following the uncorrelated harmonic oscillator Wigner distribution model described in [62]. The ground state vibrational frequencies and normal modes were evaluated at the DFT/def2-SV(P) level, using the B3-LYP functional. The initial structures were then placed between the gold clusters. Due to the metallic structure of the gold clusters, the number and density of electrode states is relatively large, which, if trajectories were allowed to actually perform surface hopping, would result in a high frequency of such events, requiring a very large number of trajectories to be analysed. To avoid the complexity arising from this phenomenon, the trajectories were kept in a single, initially chosen electron attached state while the hopping probabilities among all states were calculated in each dynamics step using the wave function overlap described in [63,64]. The calculations were performed on the ATLASZ computer cluster of Eötvös Loránd University.

Electron attached states
The first eight electron attached states of the model system are summarised in Table 1 and depicted on Figure 2 using their NTOs. One can see that the lowest six states are all localised on the gold clusters (electrode states (El x )), associated with a negative electron affinity. In the ADC(2) model the first two states show a very low value of less than −2.0 eV, followed by a rather dense set of another four states between −0.70 and −0.62 eV. The symmetric structure of the system has a consequence that in the absence of a non-symmetric perturbation (e.g. electric field or a distortion of the molecular structure) the electrode states show up as combinations of those localised on the left and right electrodes. The first electron attached state that is localised on the BDA molecule (Mol x ) is the 7th in the row with a considerably higher electron affinity (0.11 eV), dominantly attached to the molecule's lowest lying π * orbital (LUMO). Considerably higher in energy is the second molecule state found, also of the π * type with a minor component on the amino group nitrogen atoms' lone pair.
Compared to the ADC(2) data, the SOS-ADC(2) results show very similar features with all predicted electron affinities higher by 0.35 (± 0.11) eV, with no relevant differences in the states' character. This is in agreement with our previous studies [46,49] where the parent ADC(2) and CC2 methods were found to underestimate electron affinities, whereas their spin-component scaled variants predicted very accurate results as compared to CCSD and higher level references. Since the present study focuses on the energy differences and interactions of EA states, this consistent up-shift of the energies does not affect the results considerably. Therefore, in the following we present the ADC(2) results in this paper, while some SOS-ADC(2) data are provided in the Supplemental Material.

Nuclear trajectories
The sixth electrode state El 6 was chosen as the state the dynamics was running in. This is a somewhat arbitrary choice, however, as discussed in Section 2.3, in practice the Newtonian nuclear trajectories are driven by the forces on the BDA molecule's atoms and this, with the system remaining in one of the electrode states throughout the entire trajectory, means that the changes observed during the dynamics simulation are to a large extent associated with the neutral BDA molecule's vibrations and their capability to perturb the electronic structure of the molecule and the electrodes. The electrode states El 3 -El 6 being very close in energy, these latter effects can be assumed to be similar in any of them, while the highest lying El 6 state should be the most interesting from the viewpoint of its interaction with the Mol 1 state and thus, unveil a potential explicit role of the latter in the electron transport process. Nevertheless, transition probabilities were monitored between all potentially relevant pairs of states via the respective surface hopping probability. The nominal values of these dimensionless quantities are of limited physical relevance per se. Their relative change, in particular their growth by orders of magnitude at certain regions, however indicate a strong mixing and a likely nonadiabatic transition between the two states in question.
A comprehensive statistics on the trajectories are shown in Table 2 with data presented for the starting point (t = 0 fs) and two selected (50 and 100 fs) timesteps, as well as for the entire trajectory -the latter showing values averaged over all trajectories. One can observe that the state energy differences undergo a considerable variation with the molecule's structure, dropping to a fraction of the equilibrium value in certain points. Generally, the approaching state energies can be expected to increase the corresponding transition probability, and the opposite happens if the energy gap increases. These changes are reflected by the large alteration of the p values, from a practical disappearance of the transition probability up to two orders of magnitude larger than the mean starting value.
The general agreement of the mean E and p data with the starting point and the selected timesteps, as well as the mean entire trajectory is an indication of the molecule's practically unperturbed periodic vibration around the reference point, the phases of which are randomly sampled by the set of initial conditions.
No pronounced directions of structural deformations could be visually identified within the ADC(2) set of trajectories, meaning that the presence of the excess electron with no external voltage bias does not cause the system to undergo any particular structural changesinstead, an MD representation of a quasi-stationary system is observed. In such system a natural behaviour of the state energies along the trajectories is a periodic change with a frequency corresponding to the dominant normal modes -i.e. the internal motions that affect the energies to the largest extent. However, depending on the number and relative starting phase of the dominant vibrations -chosen randomly and independently for each trajectory -the change of the energies during the trajectory can either indeed show a characteristic periodicity or a much less systematic behaviour, eventually even nearly no variation. In the following, for the simplicity of the interpretation, we discuss selected trajectories that show the above periodical changes, keeping in mind that other trajectories do not possess any conceptually different behaviour either. (The presented data are available for all trajectories in the Supplemental Material.) By evaluating the displacements along the molecule's normal modes it is possible to analyse the importance of individual motions in the change of transition probabilities and state energies. This is done by transforming the dynamics snapshots' displacement vectors into the reference structure's normal coordinate representation and correlating the displacements with the investigated quantities.

Transitions to the molecule state Mol 1
A sample trajectory is depicted on Figure 3, with the energy difference of the Mol 1 and El 6 state shown, together with the associated hopping probabilities towards the Mol 1 state, as well as the POS f values of the states. One can observe that the energy gap between the states varies considerably with the molecule's vibrations, going as low as 20 percent of the equilibrium value up to 155 percent of it. This behaviour is typical for all trajectories, the average minimum and maximum values being 0.29 and 1.04 eV, respectively.
By taking a look at the correlation between the increase of the normal coordinate values and the energy gap, depicted on the top panel of Figure 4, one can see that the coordinates most responsible for altering the Mol 1 -El 6 energy gap are the 31-34 normal modes which are the joint bending motions of the BDA's ring hydrogens and the amino groups' scissoring (see the normal modes illustrated in the Supplemental Material), the latter being more significant. The respective plots of the correlation with the individual states' electron affinities on the middle and bottom panels indicate that this periodic variation of the energy gap is practically attributed to the Mol 1 state energy alone -the El 6 state is not articulately affected by these motions. This means that although the distortion of the aromatic ring and the amino groups' displacements towards the gold electrodes are inseparable in these relevant vibrations, the HCC bendings are likely to be the more important due to their destabilisation effect on the aromatic ring and thus, the reduction of the Mol 1 state energy.
The reduced energy gap is normally expected to be accompanied by an increased hopping probability and some very pronounced peaks in the probability are indeed found (see the bottom panel of Figure 3). However, the precise position of these peaks does not coincide with the energy gap minima. This is probably due to the fact that the molecule and electrode states have very different wave functions with a low overlap, which can thus remain insignificant even with the energies falling close to each other. This is true for other electrode states as well, as none of the El 3 -El 6 states was found to have a considerably higher overlap with Mol 1 . Nevertheless, the repetitive occurrence of points with an extraordinarily high transition probability (up to 86 times the mean value) is an indication that the molecule's vibrations indeed create probability for the transition to take place, even though the absolute values remain minuscule throughout the entire trajectory. This is partially explained by the absence of delocalised states, as indicated by the excess electron's positions (middle panel of Figure 3): in the electrode state the excess electron is essentially hopping back and forth between the left and right electrodes (fragments 2 and 3), remaining on one of them for longer periods of about 10 fs, while the Mol 1 state is constantly localised on the molecule (fragment 1).
The preserved quasi-orthogonality of the electrode and molecule state wave functions can be regarded as an indication that the explicit role of the El 6 →Mol 1 transition in the electron transport process is, despite how low the energy separation can become, rather unlikely.
It has to be noted that at the SOS-ADC(2) level, a slow rotation of the BDA molecule between the electrodes can also be observed during the dynamics, presumably caused by the large Coulomb repulsion of the excess electron and the molecule's electron density. Simultaneously, the energy gap of the Mol1 and El6 states slowly increases while the electrode states approach each other further with the dynamics timesteps. This effect, a possible artefact caused by our choice of modelling the external electron in the system, is found negligible at the ADC(2) level.

Transitions to other electrode states
The six states localised on the electrodes are, as discussed previously, close in energy, so the transition probabilities between them should be expected to be significant. The states are, due to the symmetry of the system, delocalised in the reference configuration, but become localised on one of the electrodes if the molecule has a non-symmetric structure during the dynamics, i.e. most of the time. Figures 5-7 show the same trajectory as Figure 3 for the El 6 / El 5 , El 6 / El 4 and El 6 / El 3 state pairs, respectively. For the former it is clearly seen that the 0.05 eV energy gaps at the reference structure reduce considerably at many steps during the trajectory, down to a quasidegeneracy of less than 0.01 eV. For the other two pairs the gaps barely drop below the equilibrium value, but show well defined minima with a characteristic periodicity. It is common to all states that they are localised on one of the electrodes most of the time, jumping to the other  with ∼ 10 fs periods, spending very short times in the delocalised range. While El 6 and El 5 jump almost simultaneously between the left and right electrodes, the El 6 / El 4 and El 6 / El 3 pairs show a pattern of periodic swapping. In the latter cases, the locations of the energy gap minima precisely coincide with the delocalised points. The peaks of the hopping probability are also predominantly located at these regions, and orders of magnitude larger than towards the Mol 1 state, indicating that with both states delocalised on both electrodes the nonadiabatic transition is facilitated. In the adiabatic picture, the back and forth oscillation of the excess electron's locality corresponds to the electron transport between the electrodes. Naturally, in the absence of an external voltage bias these swapping events show a more or less symmetric pattern whose periods correspond to those of the internal motions that drive the process. The oscillations' period can be conveniently determined by fitting a periodic sine function to the POS f curves of the states. The results, shown in Table 3, are similar for all electrode states under investigation: considering all trajectories we obtain mean periods of 24-31 fs. This corresponds to that of the 21-28 normal modes, which are also combinations of the in-plane ring hydrogen bendings and the twisting/rocking and wagging of the amino groups.
Due to the fact that all electrode states appear to be almost equally affected, the swapping effect is likely to be attributed to the periodic oscillation of the molecule-electrode interaction. This is reasonably driven by the regular reduction and elongation of the amino groups' distance from the electrodes via its wagging and rocking motions.

Conclusions
The main objective of the current study is to test a new methodology providing an alternative for the description of electron transfer in Single-Molecule Junctions, based on a picture that is more adapted to the chemical perspective. The Newtonian molecular dynamics of electron attached states simulate the fate of an excess electron in the system, allowing to see many details of the transport process, including the nature of internal motions that influence the conduction.
Our results do not support a transport mechanism involving the excess electrons localisation on the molecule as the transition probabilities are seen to be very low to and from the electrodes. The translocation from one electrode to another is more likely to happen in the electrode states directly: the periodic swapping of the excess electron's position in the adiabatic electrode states is a characteristic behaviour driven by the molecule's vibrations and their effect on the electrodes. This picture is in line with the inelastic tunnelling formulation of the electron transport in the physicist community.
We found that in this respect, the terminal groups' bending motions that alter the molecule-electrode interactions in a periodic manner are of great importance. This can potentially be exploited in the application of models where the number of considered vibrations is limited by computational tractability and the selection of relevant motions is critical.
Our approach also presents an alternative to the widespread single-orbital mean field descriptions based on density functional theory as it allows electron correlation to be incorporated at a higher level via many-body wave functions. The agreement of the ADC(2) and SOS-ADC(2) models employed in this work is, based on earlier experience, reassuring that both provide a reasonable description of electron attached states in this model and can be recommended choices for studying chemically similar single-molecule junctions.
It should be emphasised that the MD approach used in this paper is not intended to provide a quantitative model for electrical conductance. It is clear that for real electron transport events to be observed, the application of an external voltage bias would be indispensable. A molecular dynamics study of this kind could eventually be even more expressive as it would reflect an asymmetry in the observed oscillations. Therefore such calculations will be performed in our group in the near future.
Nevertheless, the present results still provide a chemically meaningful interpretation of vibronic effects in such processes, thereby facilitating the harmonisation of different viewpoints and formulations in the scientific community.