Controlling piezoresistance in single molecules through the isomerisation of bullvalenes

Nanoscale electro-mechanical systems (NEMS) displaying piezoresistance offer unique measurement opportunities at the sub-cellular level, in detectors and sensors, and in emerging generations of integrated electronic devices. Here, we show a single-molecule NEMS piezoresistor that operates utilising constitutional and conformational isomerisation of individual diaryl-bullvalene molecules and can be switched at 850 Hz. Observations are made using scanning tunnelling microscopy break junction (STMBJ) techniques to characterise piezoresistance, combined with blinking (current-time) experiments that follow single-molecule reactions in real time. A kinetic Monte Carlo methodology (KMC) is developed to simulate isomerisation on the experimental timescale, parameterised using density-functional theory (DFT) combined with non-equilibrium Green’s function (NEGF) calculations. Results indicate that piezoresistance is controlled by both constitutional and conformational isomerisation, occurring at rates that are either fast (equilibrium) or slow (non-equilibrium) compared to the experimental timescale. Two different types of STMBJ traces are observed, one typical of traditional experiments that are interpreted in terms of intramolecular isomerisation occurring on stable tipped-shaped metal-contact junctions, and another attributed to arise from junction‒interface restructuring induced by bullvalene isomerisation.


Supplementary Note 1. STMBJ blinking approach experiments and correlation maps
In the blinking approach, the STM tip is fixed at a specific distance from the surface in the presence of a dilute concentration (5 μM) of the bullvalene molecules in 1,2,4−trichlorobenzene as the solvent.To determine the distance between the two Au electrodes, a specific set-point tunnelling current is first chosen through the STM feedback loop, and the tunnelling current decay β is determined by measuring the current decay for a defined distance.In the absence of molecules, the current decay values obtained were ≈7 nm −1 , and the electrode-electrode gap distance can be evaluated using the equation (1)  =  0  − where  is the distance separating the two electrodes and (2)  0 = 2 2 /ℎ is the quantum of conductance.Given this, an estimate can be made for the junction separation from a measured value of the set-point current at a specific voltage bias.
Errors in the determination of electrode-electrode distance could arise if a high concentration of the target molecule is used, because the decay rate of the current with increasing distance can be reduced if a self-assembled monolayer forms inside the junction.To prevent such a monolayer forming, ensuring that at most one isolated molecule is in the junction region at a time, a low concentration of 5 µM of bullvalene molecules was used in the blinking experiments.
The blinking method differs significantly from the STMBJ experiment in that no physical contact occurs between the STM tip and the Au(111) surface, a fact that minimises possible tip-induced artefacts.In these experiments, the electrode-electrode gap  is set using the STM control electronics to bring the STM tip close enough to the surface to achieve the desired tunnelling current.Once this distance is established, the feedback system is disabled.Occasionally, a bullvalene molecule bridges between the tip and the surface, and this occurrence is accompanied by a sudden jump in current referred to as 'the blink".
We built the correlation (conductance) map following previous procedures. 1An automated algorithm fragments the captured conductance range in different regions (y-axis bins).The algorithm analyses the captures inside the specific bin only and, therefore, only the counts of the current decays containing plateaus inside the region are accumulated and plotted in the xaxis.The key point of this methodology is that when a plateau is detected and the complete current trace contains more plateaus at conductance different to the one studied, they are accumulated in the 2D map, thus highlighting the interdependence between the occurrences of the different current plateaus.

Supplementary Note 2. Conductance trace simulations using the 4-atom tip model
The rate of isomerisaton from isomer X to Y is given by the Arrhenius equation: (3) where Δ  ‡ is the PBE-D3BJ calculated activation energy (difference in energies between those in Figure 5c and those in Figure 5b),  = 300 K is the temperature, and  is determined such that the analogously calculated rate for isomerisation of bullvalene in solution agrees with the observed value.Solution of the coupled equations for all 20 constitutional and conformational isomers considered yields the mole fractions shown in Figure 5f,g.These kinetics simulations were performed using a time stop of 11 ns, with faster reactions that this being taken to be in equilibrium.Also, the kinetics equations were solved applying a tip retraction rate of 0.5 Å ms -1 that corresponds to the value used in the experiments.
To simulate individual conductance traces, the mole fraction at any instant in time must be constrained to depict a single isomer, i.e., instantaneous mole fractions can be only 0 or 1.To achieve this, the probabilities of reaction at each instant in time are obtained from equations of the above form, and then a random number is chosen and used to determine the composition at the next time instant.In these simulations, each "instant in time" is taken to last for a timestep of 11 ns, and 3000 such time steps are then sampled and the current averaged.This average is therefore performed over the sampling time of 0.033 ms used in the experiments to make each individual current measurement.In this process, all isomerisation reactions that occur on a timescale faster than 0.033 ms are averaged over.In general, the observed conductance in each individual measurement is found to arise from averages over multiple isomers.Initially, the Amm and its symmetrically equivalent App conformer are predicted to dominate the conductance, but subsequently the small amount of Bpm found on average to be present in each 0.033 ms time interval is predicted to dominate (Figure 5f,h).
A large number of simulated conductance traces can then be averaged, resulting in average isomer mole fractions that converge towards the directly simulated results shown in Figure 5f,g.This is demonstrated in Supplementary Figure 6 in which the results from averaging over 3000 traces (Supplementary Figure 6b) are shown and compared to those from Figure 5 that are reproduced in Supplementary Figure 6a.In contrast, averaging the conductance using the easily determined average mole fractions from Figure 5f and 5g, shown in Supplementary Figure 7a, does not lead to the same conductance histogram as obtained from averaging the traces shown in Figure 5h (and reproduced in Supplementary Figure 7b).This is because utilisation of the average mole fraction in this way does not take into account the fact that the current is averaged only over 0.033 ms intervals during each measurement.Using the average mole fractions assumes that each individual trace represents complete chemical equilibrium, whereas many reactions are too slow for this to occur (see Figure 5d).
Results from conductance-trace simulations are also shown in Supplementary Figures 6-8 obtained by varying the calculated isomer energy differences and transition-state energies in simple ways.One set of results involves adding 0.1 eV to all transition-state energies.The second involves stabilizing all A conformers by 0.1 eV, with all related transition states stabilised by half this amount, as well as analogously destabilizing all C isomers and transition states by the same amounts.The third variation involves destabilizing all C and D conformers by 0.1 eV, with again the associated transition states destabilised by 0.05 eV.Supplementary Figure 8 shows that all three of these variations significantly inhibits reverse reactions, predicting conductance traces in much better agreement with the observed traces.Also, both the original simulated histogram and those arising from the three variations display two sharp peaks in the low-conductance region (Supplementary Figure 7), in good agreement with experiment (Figure 3c).In each case, the two peaks are attributed to variations in the small amount (Supplementary Figure 6) of the highly conductive (Figure 5e) isomer Bpm found in the dynamic equilibria established on the 0.033 ms experimental timescale.The simulations performed by increasing all transition-state energies by 0.1 eV differ from the others in that 68% of the traces show A conformers dissociating from the gold contacts before isomerisation to B occurs.This is in qualitative agreement with experiment, in which over 90% of the traces did not depict isomerisation to low-conductance forms; quantitative agreement with this experimental feature can be obtained by enhancing the transition-state energy increase to 0.13 eV.In these simulations, the A to D transition-state energy is increased such that isomerisation becomes slow on the ≈3 ms time scale need to go from separations at which the reaction can occur to separations at which dissociation can occur.

Figure 7 |
Conductance histograms evaluated by different methods.a. obtained by direct solution of the kinetics equations, hence making the assumption of extensive time averaging during the sampling time.b. obtained by averaging 3000 simulated conductance traces (see Supplementary Figure 8) with 0.033 ms sampling time, reproducing results from Figure 5h.c.From 3000 traces, increasing all activation energies by 0.1 eV.d.From 3000 traces, stabilising all A conformers by 0.1 eV and destabilizing all C conformers by 0.1 eV, and associated activation energies by half of these changes.e.From 3000 traces, destabilising all C and D conformers by 0.10 eV, and associated activation energies by 0.05 eV.Source data are provided as a Source Data file.

Figure 14 |
NEGF calculated transmissions through Au-diarylbullvalene-Au linkages.The plot shows transmissions as functions of energy away from the Fermi energy computed for different conformers and stretching distances d: a) Amm, b) Apm, c) Bmm, d) Bmp, e) Bpm, f) Cmm.

Supplementary Table 1 | Calculated properties of diaryl bullvalene constitutional and conformational isomers in the gas phase.
Source data are provided as a Source Data file.

Table 2 |
Calculated relative energies, in eV, for diaryl bullvalene isomers bound between gold contacts with 4-atom tips, at cluster-model apex separation .

Table 3 |
Calculated transition-state energies, in eV, relative to the isomer energies listed in Supplementary Table 2, for diaryl bullvalene isomers bound between gold contacts with 4-atom tips, at cluster-model apex separation .