# Stochastic Simulation of smFRET Dynamic Shift from Langevin Dynamics

David Frost, Keisha Cook, Hugo Sanabria

Abstract

  * FRET

##  1 Introduction

An understanding of protein conformational dynamics has been a key feature in
biological research. In this search experimental tools from fluorescence
microscopy have been used. One of these tools is single molecule Förster (or
fluorescence) resonance energy transfer (smFRET). This tool helps quantify the
distance between two points on a given molecule. Moreover, it has been shown
that smFRET experiments can give insight into the dynamic behavior of the
underlying molecule. This is done by measuring the so-called dynamic shift in
FRET efficiency vs Lifetime plots. This dynamic shift is based on the
conditional joint distribution of the two random variables. A natural skew in
the lifetime distribution, which is conditional based on the distance,
presents a skew in the joint distribution. The theory of FRET lines details
how this skew can be related to dynamic switching in the underlying molecule.

The goal of this paper is to relate the dynamic shift to continuous state
dynamics that are based on fundamental physics. From this insight can be
gained in how the dynamic shift presents in FRET vs lifetime plots in various
physical scenarios and relate physical parameters to the dynamic shift. This
generalizes the approach mentioned as the dynamics have a continuous state
space rather than discrete. To this end we simulate a confocal smFRET
experiment while modeling the protein and dyes as evolving according to
Langevin equations. We do this in several orders of complexity while examining
the parameters used in each order. From this we see the influence of dye
linker size and dye radius as well as the influence of reaction rate dynamics
in protein conformation via the height of an energy barrier. Moreover, we also
see the influence of flexability, length, and dye position on the FRET vs
lifetime plots in the case of a randomly evolving string modelling a
disordered protein.

From this research we present a more comprehensive view of the dynamic shift
seen in smFRET experiments while also providing a wealth of observational
information on possible dynamic models for protein conformation and their
relation to actual experiments.

###  1.1 Continuous Time Markov Chains

Continuous time Markov Chains (CTMCs) are used for a wide variety of
applications to model random processes. In compartmental epidemiology models,
CTMCs model the spread of a disease or pathogen by separating populations into
compartments, studying the transitions times and rates associated with
individuals moving between the compartments of the system. The rates and
probabilities of transitioning between compartments can be used to predict
future behavior of the system. Applications related to operations research use
CTMCs to model waiting times associated with processes such as hospital
waiting rooms, production management, assembly lines, and any other service or
engineering applications. CTMCs are also used in cell biology to describe the
states involved in the process of the cell cycle, the transport of cells or
organelles inside the cell, interactions between cell populations, and many
other random processes that occur over time intervals. In essence, CTMCs can
be used to describe state transitions of random processes.

CTMCs are frequently used in the physical sciences to model dynamics of
populations, in chemistry they are invoked to find the concentrations of
different chemicals in a chemical reaction, using the so -called master
equation approach. In this work we develop a CTMC model for the FRET energy
transfer process. Using this derivations for lifetime and intensity based FRET
efficiencies can be obtained. We utilize this CTMC in simulations of a
confocal single molecule FRET experiment.

Roughly, a CTMC is a stochastic process in which transfers between states
occur at exponentially distributed times. The probability of transfer between
specific states is given by a transition matrix. The theory of CTMCs is widely
known and their properties are advantageous in computing the stationary
distribution of the process as well as hitting times and analysis of the
transition probabilities.

In this work we make use of CTMCs in two different ways. Firstly, we model the
energy transfer process involved in FRET as a CTMC, based on the fact that
fluorescent decays are exponentially distributed as well as the energy
transfer rate. Secondly, we make use of Markov chain based approximations to
simulate stochastic differential equations as is mentioned in more detail in
subsequent sections.

###  1.2 Förster (Fluorescence) Resonance Energy Transfer (FRET)

Förster (or Fluorescence) Resonance Energy Transfer (FRET) is a physical
process in which energy is transferred between a “donor” molecule and an
“acceptor” molecule nonradiatively due to electric dipole coupling. A detailed
account and derivation of FRET can be found in [Clegg, lackowickz].

Figure 1: Jablonski Diagram for FRET

We shall refer to the molecules as dyes throughout the paper. The rate of
energy transfer depends on the distance between the dipoles to the sixth
power. Due to this strong dependence on distance, FRET is often used as a
“molecular ruler”.

Consider the case in which the two dyes are completely static and the distance
between the two dyes is $r>0$. The energy transfer rate is then given by the
following

$k_{ET}=k_{D}\bigg{(}\frac{R_{0}}{r}\bigg{)}^{6}$

where $k_{D}$ is the radiative decay time for the donor dye, $R_{0}$ is a
physical parameter called the Förster radius where the probability of energy
transfer is 0.5. One can see that the energy transfer rate increases
dramatically as the distance decreases and vice versa as the distance
increases. From this sensitivity if one can estimate the energy transfer rate,
one can estimate the distance between the dyes.

The energy transfer rate can be difficult to directly estimate using lifetime
measurements as the equipement can be expensive. In replacement one can
calculate the FRET efficiency:

$\mathcal{E}=\frac{1}{\big{(}\frac{r}{R_{0}}\big{)}^{6}+1}$

this is the probability one sees a photon from the acceptor dye when the donor
is excited.

This relationship between distance and FRET has been used in biophysics and
structural biology to determine lengths between points of interest on proteins
in order to understand the conformation of the protein. Further, time resolved
FRET has been used to understand protein conformational dynamics.

###  1.3 Stochastic Differential Equations

Stochastic differential equations (SDE), also called Langevin equations, are a
type of differential equation that describes dynamics subject to noise. This
can be heuristically written as

$\frac{dx}{dt}=F(x)+\frac{d\psi(t)}{dt}$

where $\psi$ is Gaussian white noise with mean $0$ and variance $\sigma dt$.
These types of equations have long been used in physics to describe the motion
of particles subject to thermal fluctuations, in finance to model stock
prices, biology as population models, and epidemiology to study the dynamics
of a pandemic.

The area of mathematics that studies these types of equations is known as
stochastic analysis. The rigorous foundations for the expression above can be
met by an infinite amount of integrals depending on the evaluation point used
in the limit definition of the integral. The common integral used in physics
is the Stratonovich integral, rather than the Itô integral, due to its
adherence to the normal rules of calculus and the fact that it can be seen as
the limit of a Wong - Zakai limit of a particle subject to noise where the
correlation between the noise vanishes. However, it should be noted that since
we only consider the additive noise case the Itô and Stratonovich integral
will yield the same results. The expression above, interpreted as a
Stratonovich SDE, will be denoted

$dX_{t}=F(X_{t})dt+\sigma\circ dB_{t}$

as is common in stochastic analysis.

The case that is of mose importance to the current work are SDEs of the
following form:

$dX_{t}=-\nabla V(X_{t})dt+\sigma\circ dB_{t}$

where $V:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is some potential function.
These types of SDEs are called gradient systems and have advantageous
properties, they are reversible, and possess a simple stationary distribution,

$p(x)=\frac{1}{Z}e^{\frac{-V(x)}{\sigma}}$

where

$Z=\int_{\mathbb{R}^{n}}e^{\frac{-V(x)}{\sigma}}dx$

is the normalization constant. This stationary distribution is frequently
called the Gibbs distribution for the SDE.

Moreover these SDEs will exhibit metastable behavior. This is a phenomena in
which a SDE exhibits a switching behavior between multiple apparent stable
states. In the case of the gradient systems we work with the stable states
will be given by the local minima of the potential function. These local
minima are referred to as stable energy wells. An example sample path of a $3$
dimensional SDE exhbiting metastable transitions is given below.

Figure 2: Example of metastable transitions in $3$ dimensions between stable
vectors $(0,0,0)$ and $(1,0,0)$. Note the time spent in the transition state
between the stable wells.

Metastability is a key feature in protein conformational dynamcics and are an
interesting area of mathematics by itself. It is well known and shown
rigorously in the weak noise limit in[Freidlin&Wentzell] and was originally
estimated by Kramers[Kramers] that the transition time between two local
minima of the potential function will depend on the height of the barrier
between them, called the activation energy. More specifically, if
$x_{1},x_{2}$ are two stable states and $z$ is an index one saddle point
between the stable states then the approximate transition rate between the
stable point $x_{1}$ to $x_{2}$ is given by[weinan2021applied]

$\frac{|\lambda_{s}|}{2\pi}\sqrt{\frac{\det(H(x_{1}))}{|\det(z)|}}\exp\bigg{(}-\frac{\delta\Psi}{\sigma}\bigg{)}$

where $\det H(x)$ represents the determinant of the Hessian matrix of $\Psi$
evaluated at $x$, $\lambda_{s}$ is the unique negative eigenvalue of the
Hessian matrix at the saddle point, $\sigma$ is the noise present in the
gradient system, lastly $\delta\Psi=\Psi(z)-\Psi(x_{1})$ is the activation
energy.

These types of systems can be approximated by CTMCs located at the stable
points of the potential function with transition rates given by the Kramer’s
rate equation. However, in this work we consider the full continuous dynamics
as when the transition rate is large the system will spend more time in the
transition state between the wells. We believe that the continuous dynamics
used will more accurately capture these non-stable states only observed in
transitions.

###  1.4 Stochastic Dynamics

We consider the impact of several different types of dynamics on data
collected from smFRET experiments. These dynamics can be roughly ordered in
terms of complexity of the motion.

### Orders of Complexity

We consider three orders of complexity as follows

  * •

Static Case: This is the case wherein the dyes are the only parts of the
structure that exhibit dynamic behavior. In this case the dyes are model
according to the laws of an Ornstein - Uhlenbeck process. This implies that
the stationary distribution of each dye is given by a multidimensional
Gaussian. In this situation the underlying measure of interdye displacement
vector is again Gaussian distributed with known parameters. We then relate
this case to experimental situations which are applicable and we test five
variations of common dye configurations via the radius and tether length.

  * •

Metastable Case: In this case the protein exhibits metastable behavior and
transfers between a finite number of conformations randomly. This is modeled
as a nonlinear Langevin equation with equilibrium conformations given by the
minima of a potential function. The distance between the dyes is not bound by
a finite number of distances, as in the aforementioned case, but the is given
the freedom to sample the area around the equilibrium points and transfer due
to the thermal fluctuations. While this model is the closest related to the
finite state switching model wherein the dynamic shift is known, the dynamics
allow for a wider variety of behavior as well as direct physical
interpretation.

  * •

Disordered Case: The final case models a disordered protein. We do this by
simulating a discretized stochastic partial differential equation for the
evolution of a random string. This is done by simulating a large number of
connected springs subject to thermal noise. This polymer model is known as the
Rouse model. From this model the protein can sample a wide range of interdye
distances. We investigate cases in which the protein is comprised of a smaller
number of links and a larger number. Further, we consider variations in the
spring constant as well as dye placement along the string.

##  2 Methods

###  2.1 The FRET CTMC

Throughout we model the energy transfer process in FRET as a continuous time
Markov chain (CTMC) with the following rate matrix. Let $P$ be a CTMC on state
space $S=\\{D,A,F_{D},F_{A}\\}$ with rate matrix

$\begin{pmatrix}-(k_{D}+k_{ET})&k_{ET}&k_{D}&0\\\ 0&-k_{A}&0&k_{A}\\\
0&0&0&0\\\ 0&0&0&0\end{pmatrix}$

with graph representation given below.

$D$$A$$F_{D}$$F_{A}$$k_{ET}$$k_{D}$$k_{A}$

The state of the CTMC keeps track of where the excitation energy is:

  * •

State D represents an excited donor

  * •

State A represents an excited acceptor

  * •

State $F_{D}$ represents donor fluorescence.

  * •

State $F_{A}$ represents acceptor fluorescence.

The transition rates $k_{D},k_{A}$ represent the fluorescence decay times for
the donor and acceptor respectively. The transition rate $k_{ET}$ represents
the FRET energy transfer rate and is therefore given by

$k_{ET}=k_{D}\bigg{(}\frac{R_{0}}{r}\bigg{)}^{6}.$

Note that if $r=0$ then the CTMC is treated as a two state CTMC between states
$A$ and $F_{A}$ with rate $k_{A}$.

Using this CTMC model we can derive the FRET efficiency and two means of
approximating this efficiency. The FRET efficiency is the probability that the
acceptor fluoresces when the donor is excited. One can note that this can only
happen if energy transfer occurs, i.e. if in the CTMC we have a $D\rightarrow
A$ transition. This means that in a race of exponential random variables, the
$D\rightarrow A$ transition wins. So, denoting $\tau_{D}$ and $\tau_{ET}$ as
exponential random variables with rates $k_{D}$ and $K_{ET}$ respectively the
FRET efficiency is

$\displaystyle\mathcal{E}$
$\displaystyle=\mathbb{P}(\min(\tau_{D},\tau_{ET})=\tau_{ET})=\frac{k_{ET}}{k_{D}+k_{ET}}=\frac{k_{D}\big{(}\frac{R_{0}}{r}\big{)}^{6}}{k_{D}+k_{D}\big{(}\frac{R_{0}}{r}\big{)}^{6}}=\frac{1}{\big{(}\frac{r}{R_{0}}\big{)}^{6}+1}.$

The goal in most FRET experiments is to approximate this value. This can be
done in two different ways. The first is the experimental method of intensity
based FRET.[lakowicz] This method repeatedly excites the donor and measures
whether an acceptor of donor photon is observed. Effectively, this treats the
experiment as measuring the probability of success of a binomial random
variable with $p=\mathcal{E}$. It is well known that the estimator for this is
given by the number of successes observed divided by the total number of
trials. In this experiment this becomes

$\mathcal{E}_{I}=\frac{I_{A}}{I_{A}+I_{D}}$

where $I_{A}$, $I_{D}$ are the intensities (photon counts) for the acceptor
and donor respectively.

The second method is known as lifetime based FRET and comes from the
following.

$\displaystyle\mathbb{P}(\min(\tau_{D},\tau_{ET})=\tau_{ET})$
$\displaystyle+\mathbb{P}(\min(\tau_{D},\tau_{ET})=\tau_{D})=1$
$\displaystyle\rightarrow\mathcal{E}+\frac{k_{D}}{k_{D}+k_{ET}}=1$
$\displaystyle\rightarrow\mathcal{E}=1-\frac{k_{D}}{k_{D}+k_{ET}}.$

Noting that
$\mathbb{P}(\tau_{D}>t|\min(\tau_{D},\tau_{ET})=\tau_{D})\sim\text{exp}(k_{D}+k_{ET})$
[Ross] we write the following

$\mathcal{E}=1-\frac{\tau_{D}^{\prime}}{\tau_{D}}$

where $\tau_{D}^{\prime}=(k_{D}+k_{ET})^{-1}$ is the lifetime of the donor in
presence of the acceptor and $\tau_{D}=k_{D}^{-1}$ is the lifetime in absence
of acceptor. In this way one may measure the mean lifetime of the donor and
estimate the FRET efficiency, by approximating the mean of an exponential
random variable.

It is then apparent that the lifetime based approach can be easily extended to
the dynamic case, as has been done in PUT IN CITATIONS. This is done by noting
that $\tau_{ET}=\tau_{ET}(r)$ where $r$ is a random variable. Then one can
condition on $r=r_{0}$ and represent the FRET efficiency in terms of the
radial distribution. In this way one can develop the theory of FRET lines. It
is with this in mind that we conduct the following simulations. In this way
the FRET efficiency and lifetime distributions are compound distributions
conditioned on the underlying dynamics, meaning if $G_{r}$ is the measure of
the distance between the dyes as dictated by the associated Gibbs measure of
the dynamics on a state space $X$ then the FRET efficiency is simply the
expected FRET efficiency averaged over all possible states

$\mathcal{E}=\int_{X}\mathcal{E}(x)dG_{r}(x)$

further, the distribution for the lifetimes will be given by the compound
distribution of an exponential distribution with rate $k_{ET}(x)$ and the
Gibbs distribution,

$f_{\tau}(t)=\int_{X}k_{ET}(x)e^{-k_{ET}(x)t}dG_{r}(x).$

From this one can see that the difference in Gibbs distribution will influence
the FRET efficiency as well as the distribution of the lifetime. As noted, in
the case when the Gibbs distribution is a delta distribution the relationship
between FRET efficiency and the lifetime is a straightforward linear
relationship, but in the presence of a more complicated Gibbs distribution the
relationship becomes nonlinear. In the simplest case of two discrete states
with a Gibbs distribution given by

$G_{r}(x)=(\pi,1-\pi)$

this relationship has been calculated. For more other distributions such as a
Gaussian or Gaussian mixture models computational tools are generally used.

###  2.2 smFRET Experiments

As mentioned in the previous section there are two ways in which one can
estimate the FRET efficiency, intensity based and lifetime based. Both require
observing the absorbing state outcome of the energy transfer CTMC described.
In this section we will discuss the experimental set up for collection of the
data used to estimate the FRET efficiency and sources of noise and data loss
in the experiment. These considerations will be taken into account for our
simulations.

In practice one usually wishes to estimate the distance between two points on
a larger molecule, such as a protein. To achieve this one attaches fluorescent
dyes, such as Alexa 488 or Rhodamine B seen in the figure below to the
molecule of interest.

(a) Alexa 488

(b) Rhodamine B

Figure 3: Alexa 488 and Rhodamine B are two common fluorescent dyes. These
dyes are attached to a molecule via a sulfur point mutation and a methyl
chain.[Alexa][Rhodamine]

These are small molecules that are attached to the target molecule by a point
mutation and a methyl bridge consisting of several C-C links. Once the dyes
are in place one dilutes the dyed molecules in a solution called a buffer.
This is the reason for this part of the single molecule part of smFRET and why
this is done will be explained shortly. Once one has the diluted dyed molecule
one needs to repeatedly excite the donor molecule.

The excitation is done by pulsing a laser directed towards the diluted
solution of a wavelength such that the photons emitted by this laser will
excite the donor molecule.[george] These pulses happen approximately every
$50$ nanoseconds or $5\times 10^{-8}$ seconds and is called the pulse time.
The light is then passed through a linear polarizer and directed towards the
objective. The beam passes through a dichroic beamsplitterm or dichroic.

Once the light has passed through this dichroic it is focused using a convex
lens. This focusing creates an illumination profile in the sample. The
confocal volume is the region in which the illumination profile has, within a
cut off, the highest intensity and it represents the region that is “in
focus”. The confocal volume is well approximated by a $3$ dimensional Gaussian
distribution truncated at $e^{-2}$ given by[FCS]

$I=I_{0}\exp\bigg{(}\frac{-(x^{2}+y^{2}+z^{2}/k)}{w^{2}}\bigg{)}$

where $I_{0}$ is the maximum intensity, $z$ is the direction of propagation of
the light, $w$ is related to the beam radius in the $xy$ plane and $k$ is a
constant dependent on index of refraction of the solution and the numerical
aperture of the lens. This means that the confocal volume is symmetric along
the plane perpendicular to the axis of propagation of the light and asymmetric
along the axis of propagation.

The dyed molecule solution is diluted so that on average only a single
molecule of interest is within the confocal volume at a time. The molecules
then freely diffuse in and out of the confocal volume, making the time spent
in the illumination profile random. Let $V_{c}\subset\mathbb{R}^{3}$ be the
confocal volume. Each molecule in the solution undergoes a pure Brownian
motion, as it is subjected to the heat bath of the
solution.[lavenda2019nonequilibrium] So, the amount of time spent in the
confocal volume is dependent on the diffusion of the molecule. More
rigorously, consider a Brownian motion $B_{t}$ and define
$\tau_{1}=\inf\\{t>0:B_{t}\in V_{c}\\}$ and
$\tau_{2}=\inf\\{t>\tau_{1}:B_{t}\notin V_{c}\\}$ then the amount of time
spent in the confocal is given by $\tau_{2}-\tau_{1}$. The time spent in the
confocal volume is called the burst time.[george] This is due to the spike in
photons when the molecule enters the confocal volume when observing the photon
counts over time.

Provided a molecule is in the confocal volume it will be excited by an
incoming photon with a probability based on its orientation relative to the
polarization of the photon. Let $A$ be the event that the photon is absorbed
by the donor and let $\hat{\mu}_{D}$ be the electric dipole moment of the
donor and $\hat{I}$ the unit vector corresponding to the polarization of the
incident photon. Then

$\mathbb{P}(A)=\langle\hat{\mu}_{D},\hat{I}\rangle^{2}=\cos^{2}(\theta)$

where $\theta$ is the angle between the two
vectors.[kliger2012polarized][lakowicz] Notice that this distribution is
symmetric about the planes perpendicular to the polarization of the incident
photon. In this way we can see that the absorption probability increases the
more co-linear the dipole moment is with the polarization of the photon with a
maximum probability of absorption when the two are parallel.

Once excited the donor may transfer its energy via FRET, fluorescence, or
decay nonradiatively. This nonradiative decay is representative of other
energy pathways the absorbed energy may take, such as increasing vibration
between two bonds in the molecule, increasing rotation, or loss of energy due
to interactions with other molecules. The probability of emitting a photon
once excited is called the quantum yield of the molecule.[lakowicz] Denote
this by $Q_{D}$ and $Q_{A}$ for the donor and acceptor quantum yield
respectively.

If either dye fluoresces the photon will be emitted in a direction with a
distribution dependent on the electric dipole moment of the emitting dye. Let
$\hat{e}$ be the unit vector in the direction of propagation of the emitted
photon. Then the density of of the random vector $\hat{e}$ is given by

$p(e)=\frac{1}{\pi^{2}}(1-\langle\hat{\mu},\hat{e}\rangle^{2})$

where $\hat{\mu}$ is the dipole moment of the emitting dye, either the
acceptor or donor.[kliger2012polarized] Note that the distribution for the
direction of emission is symmetric about planes perpendicular to the electric
dipole moment. This shows that the emitted photon is more likely to travel in
a direction perpendicular to the dipole moment of the dye. In this way we can
see that the emitted photon will have a polarization representative of the
dipole moment of the fluorescing dye. This is the basis for fluorescence
anisotropy.[lakowicz][kliger2012polarized]

To see a photon it must travel in the direction of the detector. So, a great
loss of photons occurs due to the stochastic nature of the emitted direction.
Provided a photon travels in the appropriate direction, it is reflected by the
dichroic mirror the incident light passed through earlier. Once reflected
towards the detector it is split based on wavelength and polarization. Once
split the photon travels towards a photon detector which may or may not
register a photon based on the efficiency of the detector.[george] We refer to
the detector efficiency by the following probabilities of detection, $S_{A}$
and $S_{D}$ for the acceptor and donor detectors respectively. Once a photon
is registered is counted and time stamped with the time after the pulse and
the time in the over all experiment. This is then repeated to gather counts of
acceptor, and donor photons as well as their lifetimes.

A schematic of the confocal set up is shown below in figure 1.2.

Figure 4: An example of a confocal microscopy set up for smFRET experiments.
Taken from [george], credit to George Hamiltion.

There is a great loss of data in these experiments; approximately $99\%$ of
the photons emitted are lost due to the effects of quantum yield,
orientational factors, or detection efficiency. However, these experiments
create a large data pool with around $50,000$ photons counted during an
experiment. Once these photons are counted, one can estimate the FRET
efficiency in each burst and observe any dynamics in the underlying molecule
of interest, such as conformational changes in a protein.

###  2.3 Static Stochastic Dynamics

####  2.3.1 Spring Model

In the static case (reference section 1.4), we consider the case when the
underlying molecule is stationary in long-time scales. The dynamics in place
arise from the motion of the fluorescent dyes attached to the molecule. The
dyes are tethered to a point on the molecules via a strand of $CH_{2}$ bonds,
and their motion will influence FRET measurements.[dye1] The only forces
acting on the dye are the thermal noise and the force of the linker tethering
the dye to the protein. Therefore, we consider the dye as a mass attached to
an isotropic spring subject to Gaussian white noise. The equation of motion
for such a system in this environment may be formally written as the following
Langevin equation:

$m\ddot{X}(t)=-\kappa(X(t)-X_{eq})-\gamma\dot{X}(t)+\xi(t)$ (1)

where $\kappa,\gamma\in\mathbb{R}^{+}$ are constants representing the spring
constant of the dye linker, dependent on the number of $CH_{2}$ links used,
and friction on the dye, which is dependent on the viscosity of the
surrounding fluid, and the radius of the dye respectively. The constant
$X_{eq}\in\mathbb{R}^{3}$ is the equilibrium vector for the process and
represents the mean position of the dye. The thermal noise, $\xi(t)$, is taken
to be a white noise with variance $K_{B}T\gamma=\sigma^{2}$ where $K_{B}$ is
the Boltzmann constant and $T$ is the temperature.

Due to the low Reynolds number experienced in fluids at the molecular level,
inertial effects will be negligible. Consequently, the so-called overdamped
limit will be used. In this limit, changes in inertia, i.e., forces cancel
out, meaning
$F=m\ddot{X}=0$.[lavenda2019nonequilibrium][pottier2009nonequilibrium] In this
case, we get the following expression (after converting to an Itô SDE)

$\dot{X}(t)=-\frac{\kappa}{\gamma}(X_{eq}-X(t))+\frac{1}{\gamma}\xi(t).$

From here, we convert to a Stratonovich SDE and retrieve the equation

$dX_{t}=-\frac{\kappa}{\gamma}(X_{eq}-X_{t})dt+\frac{\sigma}{\gamma}\circ
dB_{t}.$

This type of process is well known and is called the Ornstein - Uhlenbeck
process. [Oskendal][Durrett]

Note that this is the general model for each dye, meaning the parameters in
this model are dependent on the dye used. Moreover, the equation above
considers the equilibrium position of the dye as a fixed position. This model
is useful when the dyes are fixed to a rigid molecule, such as DNA. We will
use this molecule in the first part of our simulations to determine the effect
of the physical characteristics of a dye on the dynamic shift in FRET vs.
lifetime plots.

We consider five different cases:

Table 1: Parameter table

Number of Links | Hydrodynamic Radius
---|---
3 | 0.15 nm
7 | 0.3 nm
11 | 0.3 nm
7 | 1.1 nm
11 | 1.1 nm

In this way, we can investigate the timescale separation of the dye motion in
three common scenarios, small dyes with rigid tethers and larger dyes with
more flexible tethers.

####  2.3.2 Elastic Pendulum Model

We consider a second model to contrast the more constrained spring model above
- an elastic pendulum model. Wherein the dyes are tethered at the anchor
point, and the dynamics are determined by the radial Hookean potential given
by

$V(r)=\frac{1}{2}K(r-r_{eq})^{2},$

where $K$ is the spring constant and $r_{eq}$ is the equilibrium radius. This
is a gradient system, and thus the Gibbs distribution is given by

$\frac{1}{Z}e^{-V(r)}$

a radially symmetric distribution, with Gaussian tails off the sphere of
radius $r_{eq}$. A visual example of the Gibbs distribution of this dye model
in $2$ dimensions is given in figure $5$ below.

Figure 5: $2$D stationary distribution for the radial potential above.

In this model, the dyes are constrained by the stiffness of the linker to
explore an area around a spherical region centered at the anchor point. The
angular freedom in the model introduces more dynamic behavior in the dyes.
This addition of dynamical freedom is proposed to be the source of the static-
dynamic shift. A sample trajectory of the dye motion is given in figure 6
below.

Figure 6: Sample path of the dye motion. Note how the dye mostly freely
diffuses along the sphere of radius equal to $r_{eq}$ but does not sample the
entire sphere during the burst time.

The equilibrium radius for the elastic pendulum model is determined by the
number of C-C bonds in the linker. Using $1.39\AA$ as the C-C bond length in
the methyl chain and $109.4^{\circ}$ as the bond angle, one can use the law of
cosines to deduce that a single link in the linker is $1.135\AA$. Therefore,
the equilibrium length will be $r_{eq}=1.135n\AA$, where $n$ is the number of
links in the linker.

###  2.4 Metastable Case

In the metastable case we introduce dynamics to the underlying protein. This
is done by modelling the motion of the equilibrium position in the static
model as a stochastic process. We simplify the complex protein dynamics by
considering two reaction coordinates, the anchor points for the dyes. In this
way, the dynamics of the protein are projected onto the dynamics of these two
locations. Consider when these reaction coordinates denoted $P^{d}$ and
$P^{A}$ evolve according to a nonlinear SDE with additive noise and drift
according to the gradient of a potential function. [Energy_landscape]

The reaction coordinate processes are given by

$\displaystyle dP^{D}_{t}=-\nabla\Phi(P^{D}_{t})dt+\sigma_{P^{D}}dB_{t}$
$\displaystyle dP^{A}_{t}=-\nabla\Psi(P^{A}_{t})dt+\sigma_{P^{A}}dB_{t}$

where $\Phi:\mathbb{R}^{3}\rightarrow\mathbb{R}$ and
$\Psi:\mathbb{R}^{3}\rightarrow\mathbb{R}$ are potentials representative of
the free energy landscape of the protein.

In this work we consider several cases of a double well potential wherein the
system exhibits transitions between two metastable states. In all cases we
consider a simplified system of reaction coordinates in which the
metastablilty occurs only along the x - axis in $\mathbb{R}^{3}$. This is
without loss of generality as we can always choose a coordinate system such
that the line connecting the two metastable points is defined as the x-axis.
Moreover, we only consider metastability occurring for the acceptor anchor.
The donor behaves as a noisy spring similar to the dye model. In this way
there are only two possible stable locations for the system corresponding to
the acceptor anchor transition along the x-axis.

The three cases of potential landscape geometry are shown below along with
their formulas.

  * •

The symmetric Double Well:

$\Psi(x,y,z)=E((x-W)^{2}-W)^{2}+y^{2}+z^{2}$

this gives a double well potential in which the height of each well is the
same, $E$, and the locations of the wells are given by $0$ and $W$. This
represents a case in which the motion in each well is the same, and the system
spends equal amounts of time in each well.

  * •

Asymmetric Height Double Well:

$\Psi(x,y,z)=E((x-W)^{2}-W)^{2}+Cx+y^{2}+z^{2}$

this gives a double well potential such that the energy barrier height differs
depending on the well the system is in. This represents a system where the
dynamics stay in one stable location longer than another, but still exhibits
similar behavior in each well.

  * •

Asymmetric Width Double Well: A cubic spline approximation to a asymmetric
width double well is used.

Figure 7: Example of the asymmetric width double well.

While both wells are the same height, one well is wider and hence will have a
smaller hessian value at the local minimum, and moreover the dynamics can
sample a wider area in one well as opposed to the other. Lastly, in this
situation oscillations occur at a fast rate in the well near one rather than
the well around 0.

This potential does give problems however, since we use splines to construct
the landscape, the Kramers rule is not readily available and the transition
state analysis is not as readily available as in the other cases. Due to this
the energy barrier heights differ from the previous two cases.

###  2.5 Disordered Case

In the disordered case the underlying protein is considered as a flexible
polymer. We model this polymer by considering a system of $C$ overdamped
springs with equal spring constant $k$. Each spring represents a part of the
protein, and the location of this part is given by a SDE governed by the
neighbors in the chain. We consider each spatial coordinate to evolve
independent of the others.

Figure 8: A simple diagram of the Rouse model with 13 links NEED TO MAKE A
DIFFERENT ONE! Each blue ball represents a $U\in\mathbb{R}^{3}$ where each
spring has a spring constant $k$.

Consequently the $i^{th}$ coordinate of each component of the chain is
governed by a linear system of $C$ SDEs

$dU^{i}_{t}=\mathbf{K}U^{i}_{t}dt+\boldsymbol{\sigma}dB_{t},~{}~{}i=1,2,3$

where $\mathbf{K}\in\mathbb{R}^{C\times C}$ is the spring constant matrix,
which is tridiagonal with $-k$ on the diagonal and $k$ on the off diagonal.
The noise component is driven by independent thermal fluctuations and
therefore $\boldsymbol{\sigma}=\sigma I$ where $\sigma$ is the standard
thermal noise coefficient and $I$ is the identity matrix.

Since this in each component this is a linear SDE in $\mathbb{R}^{C}$ driven
by Gaussian noise the solution is readily available in terms of a matrix
exponential, and the stationary distribution will be a multivariate Gaussian
with covariance given by the solution to the Lyaponov equation

$\mathbf{K}\boldsymbol{\omega}+\boldsymbol{\omega}\mathbf{K}^{T}=\boldsymbol{\sigma}\boldsymbol{\sigma}^{T}.$

Since $\boldsymbol{K}$ is symmetric and
$\boldsymbol{\sigma}\boldsymbol{\sigma}^{T}=\sigma^{2}I$ the solution to this
equation is given by $\boldsymbol{\omega}=-\sigma^{2}\boldsymbol{K}$.
Consequently, the stationary distribution for the system of equations is a
multivariate Gaussian with mean vector $0$ and covariance matrix
$-\sigma\boldsymbol{K}$. Note that this is in $C$ dimensional space, this does
not mean that the mean length of the chain is $0$.

###  2.6 Computational Methods

The goal of this project is to simulate smFRET experiments to relate first
principles physics to the dynamic shift in lifetimes. In this section, we will
explain the methods by which we conduct our simulations and generate synthetic
smFRET data.

To simulate the random nature of the number of bursts seen in a smFRET
experiment we define the total time in seconds the experiment is conducted and
generate a Poisson distributed random number for each of these seconds. We use
$\lambda=1$ to parameterize the arrival rates of bursts per second. The total
number of bursts, $N_{B}$, is defined as the sum of these outcomes.

For each of the $N_{B}$ bursts a sample path of the dynamics presented in the
previous section is simulated with an iteration step $\Delta$ which is
representative of the pulse time. Each individual sample path is simulated for
a total time $T$ where $T$ is a Gamma distributed random variable meant to
approximate the burst time and generated for each sample path. At each
iteration, corresponding to pulse time in the experiment, a stopped version of
the FRET CTMC is simulated to represent the absorption, energy transfer,
emission, and collection of a photon from a single pulse of the laser. The
result of this CTMC determines if the photon was collected, the color of the
photon, and the lifetime associated with the photon. The results from the CTMC
for each pulse are collected for the total burst. Using the results from each
burst intensity-based FRET and lifetime-based FRET estimations can be
conducted.

To simulate the dye and protein dynamics the Euler-Maruyama method is used. A
pseudocode for this algorthim is given below for simulating the general SDE

$dX_{t}=b(X_{t})dt+\sigma(X_{t})dW_{t}$

Algorithm 1 Euler - Maruyama Method for SDEs

$\Delta$ \- Time step, $T$ \- Total time

$X(0)\leftarrow x$ $\triangleright$ $x\sim F$ for some distribution $F$ is
also acceptable.

$\delta W(0)\leftarrow 0$ $\triangleright$ Set to $0$ since initial
distribution is in the previous step.

for $n$ from $1$ to $N$ do

$\delta W(n)\sim\mathcal{N}(0,\Delta)$

$X(n+1)\leftarrow X(n)+b(X(n))\Delta+\sigma(X(n))\delta W(n)$

end for

This method converges strongly, i.e. in mean square meaning,

$\max_{t\in[0,T]}\mathbb{E}[|X_{t}^{\Delta}-X_{t}|^{2}]\leq C\Delta$

where $X_{t}^{\Delta}$ is the approximating process.[weinan2021applied] This
means that the algorithm converges in mean square with rate
$\mathcal{O}(\Delta)$. Moreover the algorithm converges in expectation
meaning,

$\max_{t\in[0,T]}|\mathbb{E}[f(X_{t}^{\Delta})]-\mathbb{E}[f(X_{t})]|\leq
C_{f}\Delta$

for all $f\in C^{\infty}_{b}(\mathbb{R}^{d})$. So, in both senses, the
algorithm converges with a rate $\mathcal{O}(\Delta)$.[weinan2021applied]

In each order of complexity the initial positions of the dye-protein system
are chosen at random from the stationary distribution. In the static spring
setting a standard Gaussian is used. For the elastic spring model, and the
metastable models a multivariate accept-reject method is used to simulate data
from the Gibbs distribution. The dye processes in the metastable cases are
initialized following their Gaussian Gibbs distribtuion with mean at their
respective anchor process initial value. For the disordered case, the initial
configuration is simulated from a $3C$ dimensional Gaussian, with covariance
given by $\Sigma=\sigma^{2}\mathbf{K}^{-1}$. The components are then summed to
generate the intitial configuration of the polymer.

With each iteration, the processes dependent on these SDEs are also simulated.
Namely, at each iteration the vector displacement $L_{t}^{\Delta}$, the scalar
distance $r_{t}^{\Delta}$, the energy transfer rate $k_{ET}^{\Delta}(t)$ and
the FRET efficiency $\mathcal{E}_{t}^{\Delta}$ are calculated. Using these
processes a simulation of the energy transfer process and detection process is
conducted using the simulated $k_{ET}$. This process incorporates the loss of
photons via the quantum yield, the detector efficiency and loss by
orientational factors. The pseudocode for this is given below. This is a
modified version of the CTMC code in [dobrow2016introduction]

Algorithm 2 Energy Transfer Simulation

$\Delta$, $t$ $k_{ET}^{\Delta}(t)$ \- time step in algorithm and energy
transfer rate at a multiple of . $k_{D}$, $k_{A}$ \- fluorescence rates for
donor and acceptor resp. $Q_{D}$, $Q_{A}$ \- Quantum yield of donor and
acceptor resp. $S_{D}$, $S_{A}$ \- Detector efficiency for donor and acceptor
channels. $P_{ND}$ \- a probability based on the photon loss due to
orientation. state $\leftarrow D$

while $t<\Delta$ state $\neq ND,F_{D},F_{A}$ do111$ND$ represents non-
detection, $F_{D},F_{A}$ represents a donor or acceptor observation
respectively.

if state $=D$ then

Generate two Bernoulli random variables with probabilities $(1-Q_{D})$ and

$P_{ND}$ if either is $1$ state $\leftarrow ND$ and $\tau=0$.

Generate $\tau_{DF}\sim\exp(k_{D})$ and $\tau_{ET}\sim\exp(k_{ET}(t))$

end if

if $min(\tau_{D},\tau_{ET})=\tau_{D}$ then

$\tau=\tau_{D}$

state $\leftarrow F_{D}$, break

else

$\tau=\tau_{ET}$

state $\leftarrow A$

end if

if state $=F_{D}$ then

Generate Bernoulli random variable with probability of success $(1-S_{D})$

if $1$ then state $\leftarrow ND$.

end if

if state $=A$ then

Generate Bernoulli random variable with probability of success $(1-Q_{A})$

if $1$ then state $\leftarrow ND$.

Generate $\tau_{AF}\sim\exp(k_{A})$.

$\tau=\tau+\tau_{AF}$

State $\leftarrow F_{A}$, break

end if

if state $=F_{A}$ then

Generate Bernoulli random variable with probability of success $(1-S_{A})$

if $1$ then state $\leftarrow ND$.

end if

end while

Collect the outcome:

color[t] = state

lifetime[t] = $\tau$

For each time step $t$ the output is a lifetime and color for the associated
photon with lifetimes of $0$ corresponding to non-detection. From here,
conduct intensity-based and lifetime-based FRET analysis on the data
generated.

##  3 Results

In this section, the resulting FRET vs Lifetime plots generated from the
simulated data are discussed.

###  3.1 Static Case

The standard Gaussian spring model for the dye fluctuation yields surprising
results. From the simulations, as seen in figure 9, it is apparent that the
motion of the dyes is not variable enough to produce a noticeable dynamic
shift. The spring constant for the linker is too strong to generate enough
movement around the mean position. This can be seen as an “averaging out” of
the dynamics. Due to the highly oscillatory motion induced by the high spring
constants, the dye moves between extremes of its distribution within the same
burst. Due to this, the asymmetry to long lifetimes is lost. This prevents any
large dynamic shift.

(a) Dye with 11 Links and 1.1 nm in radius

(b) Dye with 7 links and 1.1 nm radius

(c) Dye with 7 links and 0.3 nm radius

(d) Dye with 11 links and 0.3 nm radius

(e) Dye with 3 links and 0.1 nm radius

Figure 9: FRET vs Lifetime plots for the five cases of dye configuration

Moreover, simulations were conducted using anisotropic springs to induce an
asymmetry in the motion. However, this produced similar results. Further, it
should be noted that the backbone noise has also been taken into account.
Using molecular dynamics simulations it has been shown that the average
oscillation of the backbone when considered as DNA is less than that of the
dye. When coupling these dynamics a dynamic shift is still not observed.

(a) Slow transition $\approx 10$ms

(b) Fast transition $\approx 0.1$ms

Figure 10: plots of….

This suggests something is wrong with the standard Gaussian model for dye
motion. When considering the system’s physical parameters, there is not enough
variability in the dye motion to produce a noticeable dynamic shift.

###  3.2 Elastic Pendulum

The elastic pendulum model is suggested to replace the aforementioned Gaussian
spring model. While the dyes are still constrained in their radial motion by
the rigidity of the linker, their angular motion is free. This allows the dye
to explore its accessible volume without the influence of intense
oscillations. Due to this, the dye motion will not average out as the dye will
explore only a region of the accessible volume during a burst time and not the
whole volume.

### Metastable Case

The metastable case mostly adds to the current understanding of FRET lines.
The overall shape appears the same as the discrete state approximations but
with some added structure due to the continuous dynamics.

####  3.2.1 Symmetric

(a) Slow transition $\approx 10$ms

(b) Moderate transition $\approx 1$ms

(c) Fast transition $\approx 0.1$ms

Figure 11: plots of….

#### Asymmetric Height

(a) Slow transition $\approx 10$ms

(b) Moderate transition $\approx 1$ms

(c) Fast transition $\approx 0.1$ms

Figure 12: plots of….

####  3.2.2 Asymmetric Width

(a) Slow transition $\approx 10$ms

(b) Moderate transition $\approx 1$ms

(c) Fast transition $\approx 0.1$ms

Figure 13: plots of….

###  3.3 Disordered Case

(a) 100 links, $k/\sigma\approx 1$

(b) 25 links, $k/\sigma\approx 0.1$

Next we investigate how different placements of the dyes along the polymer
impact the FRET vs lifetime distributions. Below we have two different
placement schemes, one where the donor is placed in the middle of the polymer
and the acceptor is placed at an end. For the other scheme both dyes are
placed at quarters of the polymer length. In this way we can see the impact of
having a dye attached at a free end and the other is more constrained, as well
as the impact of both dyes being constrained.

(a) 100 links $k/\sigma\approx 1$, donor - 50 and acceptor - 100

(b) 100 links $k/\sigma\approx 1$, Donor - 25 and acceptor - 75

## Conclusion


# Unraveling the Interplay between Quantum Transport and Geometrical
Conformations in Monocyclic Hydrocarbons Molecular Junctions

A. Martinez-Garcia  Departamento de Física Aplicada and Instituto
Universitario de Materiales de Alicante (IUMA), Universidad de Alicante,
Campus de San Vicente del Raspeig, E-03690 Alicante, Spain.  T. de Ara
Departamento de Física Aplicada and Instituto Universitario de Materiales de
Alicante (IUMA), Universidad de Alicante, Campus de San Vicente del Raspeig,
E-03690 Alicante, Spain.  L. Pastor-Amat  Departamento de Física Aplicada and
Instituto Universitario de Materiales de Alicante (IUMA), Universidad de
Alicante, Campus de San Vicente del Raspeig, E-03690 Alicante, Spain.  C.
Untiedt  Departamento de Física Aplicada and Instituto Universitario de
Materiales de Alicante (IUMA), Universidad de Alicante, Campus de San Vicente
del Raspeig, E-03690 Alicante, Spain.  E. B. Lombardi  Department of Physics,
Florida Science Campus, University of South Africa, Florida Park, Johannesburg
1710, South Africa  W. Dednam  Department of Physics, Florida Science Campus,
University of South Africa, Florida Park, Johannesburg 1710, South Africa  C.
Sabater  Electronic mail address<EMAIL_ADDRESS>Departamento de Física
Aplicada and Instituto Universitario de Materiales de Alicante (IUMA),
Universidad de Alicante, Campus de San Vicente del Raspeig, E-03690 Alicante,
Spain.

###### Abstract

In the field of molecular electronics, particularly in quantum transport
studies, the orientation of molecules plays a crucial role. This orientation,
with respect to the electrodes, can be defined through the cavity of ring-
shaped monocyclic hydrocarbon molecules. In this manuscript, we unveil the
geometrical conformation of these molecules when they are trapped between two
atomically sharp electrodes through a combination of dynamic simulations,
electronic transport calculations based on density functional theory, and
break junction experiments under room conditions. Moreover, we present a novel
criterion for determining the molecular orientation of benzene, toluene,
(aromatic) and cyclohexane (aliphatic) solvents. Our findings for the
identification of the molecular orientations on gold metal nanocontacts and
their associated transport properties, can improve the understanding of
molecular electronics using more complex cyclic hydrocarbons.

##  I Introduction

All the instruments and techniques used to measure the electronic transport of
atomic-sized contacts were developed over more than three decades agoAgraït
_et al._ (1993); Krans  _et al._ (1993). At the beginning of that period the
electronic transport of atomic contacts made of metalsKrans  _et al._ (1995),
semi-metals Krans and van Ruitenbeek (1994), and superconductorsScheer  _et
al._ (1998) were studied, leading to an immediate interest in molecular
junctions. The molecular electronics field initially considered very simple
molecules like diatomic hydrogen and deuterium Smit  _et al._ (2002), which
paved the way for obtaining knowledge from an experimental point of view. For
various reasons and based on attempts at trying to identify exotic properties
in the electronic transport field, innumerable types of different molecules
have been studied Venkataraman  _et al._ (2006); Li  _et al._ (2022); van der
Zant  _et al._ (2006); Reznikova  _et al._ (2021); Yelin  _et al._ (2013);
Taniguchi  _et al._ (2011); Yoshida  _et al._ (2015); Tewari  _et al._ (2019);
Bopp  _et al._ (2017); Cuevas and Scheer (2017); Agraït  _et al._ (2003);
Evers  _et al._ (2020); de Ara  _et al._ (2022). The creation of a molecular
junction basically depends on its physical state and the method of deposition
of the molecules, which can be grouped into three categories. The first one
corresponds to the gas phase molecules, which are associated with the unique
possibility of blowing the gas close to the electrodes. The second group
comprises molecules in the solid state that can be deposited by thermal
evaporation onto the junctions Yelin  _et al._ (2021). The final group is
composed of molecules that can be delivered in solution, and then, the usual
way to capture the target molecules is to wait for the solvent to evaporate.

It is still widely believed that the solvent disappears or has a negligible
effect on electronic transport after evaporation van Veen  _et al._ (2022).
However, recent works de Ara  _et al._ (2022) have demonstrated that the
solvent’s presence is indeed detectable via Scanning Tunnelling Microscopy
(STM) topography and STM-Break Junction (BJ) experiments. In our previous
published experimental results de Ara  _et al._ (2022), we unmasked the
molecular electronic signature of the aromatic (benzene, toluene) and
aliphatic (cyclohexane) moleculesQuek  _et al._ (2009); Venkataraman  _et al._
(2019). However, a systematic study in solvents exploring the relationship
between electronic transport and the geometric orientation of the molecules
has yet to be conducted.

Recent works based on atomically precise binding conformations Pan  _et al._
(2022) have inspired us to clearly identify the relationship between
electronic transport and molecular orientation by using Classical Molecular
Dynamics (CMD) simulations and Density Functional Theory (DFT) electronic
transport calculations de Ara  _et al._ (2022). Based on our simulations, we
have classified the typical final contact for benzene, cyclohexane, and
toluene. Furthermore, we have calculated the electronic transport using DFT
for several scenarios provided by CMD. Finally, we compare our simulations and
calculations with electronic transport experiments. The good agreement of this
comparison helps us to unmask the relationship between the conductance and the
orientation of molecules between the gold electrodes.

##  II Materials and methods

###  II.1 Identifying the orientation of the cyclic molecules

From a mathematical point of view, the normal vector is the simplest way to
identify a plane. Thanks to the fact that the aromatic solvents (such as
benzene and toluene) and aliphatic solvents (such as cyclohexane) analyzed in
our manuscript are monocyclic in nature, we can use their internal cavity as
an orientation plane. Its normal vector can help to identify if we are in
parallel or in perpendicular configuration with respect to the alignment
vector defined by the electrodes. In Fig. 1 blue arrows indicate the normal
vector over the plane defined by the cavity of the benzene molecule. The red
line indicates the direction of the vector defined by the alignment of the
electrodes, in this case, dominated by the center of the apex atoms of the
upper and bottom gold electrodes. As an adopted criterion we have defined
parallel configuration when the normal vector is in the same direction of the
alignment of the electrodes (see left illustration of Fig. 1). On the other
hand, when the angle formed by the normal vector with respect to the alignment
vector direction is between 16 and 90 degrees, we say that we are in the
perpendicular configuration (see right illustration of Fig. 1, where the
extreme case of 90 degrees between both lines is depicted).

Figure 1: Ilustration of“parallel”and “perpendicular” configurations of
benzene molecule between two gold electrodes

###  II.2 Calculations of conductance-based on Density Functional Theory and
Classical Molecular Dynamics simulations

To calculate the conductance, we employed DFT combined with the non-
equilibrium Green’s function (NEGF) approach to quantum scattering Louis  _et
al._ (2003). The electronic transport calculations were performed using the
well-established code ANT.GAUSSIAN Palacios  _et al._ (2001, 2002, ). This
code is built on top of and interfaces with GAUSSIAN09. To ensure the high
quality Zotti  _et al._ (2022) of our calculations, we used the LANL2DZ basis
setWadt and Hay (1985) for the atoms of the molecules and a few of the
adjacent metal layers on either side of the molecules, while the outer metal
layers were described by the smaller tight-binding-like CRENBS basis set Ross
_et al._ (1990).

In this work, we have calculated the electronic transport in two scenarios:
(i) For idealized initial test structures such as those shown in Fig. 2 and
(ii) for electrodes stretched by CMD simulations as described below. For each
of these scenarios, the DFT quantum transport methodology remained the same.

Figure 2: The upper left panel shows the type of gold point contact used. The
orange arrow indicates the direction of the pull-off of the upper electrode.
The blue rectangle next to it shows the initial structures of the benzene,
cyclohexane and toluene junctions for one of the two configurations studied in
this work. The starting distance between the two gold apex atoms of the tips
in the “parallel” configuration is depicted with dashed lines and is $2$ Å.
Three of the “perpendicular” configurations have a starting distance between
the gold tip apex atoms, indicated by dashed lines, of $6$ Å. In the fourth
“perpendicular” configuration, the toluene molecule with the methyl group
attached to the upper electrode has gold tip atoms that start apart $7$ Å. The
“perpendicular” configuration of the toluene is labelled as L.R. and U.D.
which indicates if the methyl group is located left-right or up-down.

A simulation based on classical molecular dynamics essentially solves Newton’s
second law in order to obtain the trajectories of all the atoms involved in
the phenomena it is meant to describe Allen and Tildesley (1989); Frenkel and
Smit (2002). We have traced the evolution of a collection of atoms at a very
fine time resolution ($\sim 1$ fs). To simulate the rupture-formation cycles
of the gold nanowire with the target molecules, we used the LAMMPS
codePlimpton (1995); Plimpton  _et al._ with a ReaxFF potential Järvi  _et
al._ (2011); Senftle  _et al._ (2016), which reproduces quite well the
mechanical and catalytic behaviour of metal-organic compounds Soria  _et al._
(2018); Gomzi  _et al._ (2021). In order to reproduce the experimental
conditions, the simulations were performed at a constant temperature of 300 K
with a Nosé-Hoover thermostat. The simulated nanowire is created with an area
of narrower cross-section ($\approx$ 17 x 17 x 25 Å3) grown along the (001)
crystallographic direction Sabater  _et al._ (2018), where the rupture will
take place. In this area of constriction, we place four molecules on each of
the four square crystallographic facets of the gold nanowire (see the ”Input”
panels in Fig. 4); hence, sixteen molecules are deposited in total for each
initial condition. This is repeated for each of Benzene, Cyclohexane, or
Toluene in turn.

###  II.3 STM-Break Junction experiments

The molecular electronics experiments were performed by STM-BJ approach under
ambient conditions. As electrodes, we used two gold wires (0.5 mm in diameter
and of 99.99% purity as supplied by Goodfellow Goodfellow (city)). The voltage
applied in all the experiments was 100 mV. The current through our molecular
junctions was converted into volts in three stages of the current-voltage
amplifier with respective gains of ${10^{6}}$, ${10^{8}}$, and ${10^{9}}$
$V/A$, in order to be recorded by our DAC system. This homemade instrument was
tested and used in a previous publication de Ara  _et al._ (2022). Knowing the
applied bias voltage and current, we can easily calculate the conductance,
which is expressed in quantum of conductance ($G_{0}=2e^{2}/h$), where factor
two comes from the spin degeneracy, $e$ is the charge of the electron and $h$
is the Planck’s constant. Historically, the representation of conductance vs
the relative displacement of the electrodes has been called a “traces of
conductance”. These traces can be labelled as “rupture” or “formation”,
depending on whether the electrodes are pulled apart or pushed together,
respectively. In this manuscript, we have only studied the rupture traces. As
internal protocol, we always use ultra-clean gold electrodes, and we verify
the extent of “cleanness” via a logarithmic histogram, which for clean
electrodes is characterised by the non-existence of any peak between the
atomic contact ($\approx$ 1 $G_{0}$) or any level of noise in the I-V
converter. Only after obtaining ultra-clean gold, the organic solvent is
deposited via drop casting over electrodes de Ara  _et al._ (2022).

##  III Results and discussion

To study the change in conductance with relative displacement and compare the
orientation of molecules between the electrode tips, we assumed that the
molecule and electrode atoms are static and that the only movement allowed is
the bulk displacement of the upper electrode (as shown in Fig. 2, where the
gold contact and its orange arrow indicating the direction of movement are
present). We have studied the evolution of the conductance vs the relative
displacement of the upper electrode for benzene, cyclohexane and toluene
molecules in the parallel and perpendicular configurations (see Fig. 3). In
all the cases we have moved the upper electrode in step intervals of $0.2$ Å
without relaxation. At every step, we have calculated the electronic transport
(conductance) by DFT+NEGF. In Fig. 3, the three panels show the conductance in
units of $G_{0}$ versus displacement up to $2.5-3.6$ Å (depending on the
molecule). This range has been selected to exclude the tunnelling regime
($2.5-3.6$ Å to $6-7$ Å) due to lack of information of interest (see for full
range and details sup. inf. Fig. S1.).

Fig. 3 upper, middle and bottom panels show the evolution of the conductance
for benzene, cyclohexane and toluene, respectively. In each panel, spheres-
lines represent parallel configurations and triangles-lines perpendicular
configurations. Moreover, the colour code used in Fig. 2 for every molecule
also carries over to these panels a, b, and c. In direct correspondence with
Fig. 2, the “parallel” and “perpendicular” structures are indicated by blue
and grey shaded areas, respectively. The height of this rectangle indicates
the maximum and minimum values attributable to each configuration. This colour
code of blue and grey facilitates understanding of the experimental data.

Figure 3: Calculated electronic transport vs relative displacement of the
initial structures of benzene, cyclohexane, and toluene, labelled as a,b, and
c, respectively. Coloured markers represent the calculated values of benzene
(black), cyclohexane (pink) and toluene (green).

From Fig. 3, one of the first observations we can make is that “parallel”
configurations of benzene and toluene exhibit very similar conductance traces
during stretching. Additionally, the calculated conductance traces for the
perpendicular cases of benzene and toluene “perpendicular-U.D.” show the same
behaviour. These results suggest that the electronic transport in these
aromatic molecules is relatively similar. By contrast, the aliphatic molecule
cyclohexane shows different electronic transport values of around $\sim
10^{-3}G_{0}$ for a displacement of $2.5$ Å. These “toy-model” conductance’s
calculations of idealized structures by DFT+NEGF are not representative of the
vast quantity of possible positions that can be adopted by the molecule during
the experiment. For this reason, we have decided to simulate via CMD
simulations the different types of contacts that can likely be produced.

Fig. 4 shows a summary of the snapshot of the cycles of rupture and formation
obtained by CMD simulations. Each row corresponds to a different molecule, as
indicated by the colour code (black, purple, and green for benzene,
cyclohexane, and toluene, respectively). At the top of each snapshot, we
indicate the simulation cycle and step in units of ksteps (where 1 step equals
1 femtosecond). The first column depicts the initial input structure used. For
each molecule, multiple junction rupture and formation cycles are repeated.
Representative snapshots are taken from different cycles at various steps in
the cycle, illustrating the different positions that the molecules can take.

Figure 4: Snapshots of the simulation during the process of rupture and
formation cycles. Each row corresponds to a different molecule: benzene (top),
cyclohexane (middle) and toluene (bottom).

We analyzed only the rupture cycles thanks to experimental data, that only is
based on rupture traces. In particular, we have focused on the last point
immediately before rupture and we have created the following classification
hierarchy.

Figure 5: Classification hierarchy of molecular contacts used in this work. To
the left and right, we show example structures that illustrate each type of
last contact.

At the top of this hierarchy, the last contact before rupture can be purely
metallic or purely molecular or a combination of both. The red blocks shown at
the top in Fig. 5 represent these three common types of final molecular
contacts. The structures on the left and right of the red blocks depict final
contacts composed of pure metal and a mixture of metal and molecules,
respectively. When the final contact is purely molecular, two possibilities
can occur: either a single molecule or multiple molecules are held by the
nanojunction. This secondary classification of the pure molecular last contact
is shown in the second row of Fig. 5 and is represented by yellow blocks in
the diagram. Moreover, the structures to the left and right of these yellow
blocks illustrate single or multiple molecular junctions respectively. As Fig.
2 illustrated, a single molecule can be captured either in a “parallel” or a
“perpendicular” configuration by the junction, as indicated by the purple
blocks in the final and lowest level of the classification hierarchy. The
structure on the left of this last row illustrates the parallel orientation,
while that on the right corresponds to the perpendicular one.

Thanks to our classification, we can determine how many times a given type of
last molecular contact occurs in our simulations (see table 1). The first
column indicates the molecule used in the simulation. The second column
corresponds to the total number of rupture cycles analyzed. The third column
represents the number of events in which we found a purely metallic contact at
the moment of last contact. The fourth column is the number of events that we
have identified as last contact involving a molecular contact (independently
of whether the contacts contain single or multiple molecules). Finally, the
fifth column shows the number of events in which we observed combined
molecular and metallic contacts.

Table 1: Classification of the last type of contact during the simulated process of the rupture. Molecule | Cycles | Metallic | Molecular | Met+Mol
---|---|---|---|---
Benzene | $15$ | $3$ | $4$ | $8$
Cyclohexane | $11$ | $0$ | $11$ | $0$
Toluene | $16$ | $1$ | $10$ | $5$
Figure 6: Percentage calculated from table 1 for the different types of last
rupture in the benzene, cyclohexane and toluene.

Fig. 6 represents in a bar graph the percentage of the events obtained in
table 1, the colour code used is the same as in previous plots. From this
figure, we observe clearly that the most probable last contact for the
cyclohexane and toluene is the single-molecule junction, and for the case of
benzene is the combination of the metallic-molecular junction. However, from
the results shown in table 1 and Fig. 6, it is difficult to deduce a refined
classification of the type of pure molecular junction that is possible. For
this reason, we present table 2 where the first column indicates the molecule,
the second and third columns correspond to the single parallel and single
perpendicular configurations, respectively, and the fourth column corresponds
to the multiple molecular junctions that can be found at the moment of last
contact.

Table 2: Refined classification of the pure molecular junctions of our simulation results. The number of events of Single Parallel (Sin. ara.), Single Perpendicular (Sin. Perp.), and Multiple Molecular (Mul. Mol.) junctions. Molecule | Sin. Para. | Sin. Perp. | Mul. Mol.
---|---|---|---
Benzene | $2$ | $1$ | $1$
Cyclohexane | $0$ | $7$ | $4$
Toluene | $3$ | $3$ | $4$

In order to ascertain the accuracy and reliability of our simulations and
conductance calculations, we compare the simulation results with the
experimental data. In each experiment, we first study bare gold for which we
usually collect over 2000 trace files as reference. Once we have finished the
analysis of pure gold, we drop cast the molecules over the bare gold and start
the process of breaking and reforming of the junction, whilst at the same time
collecting up to 10000 traces of conductance. Typically, molecular signatures
are observed under the 1 G0 region, so a logarithmic scale is usually most
appropriate for analyzing the data statistically. A representative collection
of experimental single rupture traces as a function of the stretching of the
junction is shown in Fig. 7 panels a, b, and c. Each trace exhibits a
different initial plateau-like behaviour followed by a plateau at lower
conductance values or in the tunnelling regime, related to the bridging of the
electrodes by the molecules. Moreover, here we also display within the plots
the blue and gray shading obtained from the DFT calculations shown in Fig. 3.
These shadings help us identify if the plateau aligns with parallel or
perpendicular configurations.

Figure 7: In panels a, b, and c, we show rupture traces for benzene,
cyclohexane and toluene, respectively. Grey and blue shaded areas represent
the range of values of the theoretical traces calculated via DFT based on the
different geometric configurations tested. Bottom panels d,e, and f show one-
dimensional experimental histograms of conductance for benzene, cyclohexane
and toluene.

When we plot conductance histograms in the linear scale, we also observe humps
or shoulders around the conductance of the metallic contact (1 $G_{0}$). Fig.
7 panels d, e, and f show the experimental histograms for each molecule in
dark grey and for bare gold in dark gold colour. The signals around 1 G0
conductance can be ascribed to two different scenarios. The first could be to
a parallel configuration of the molecules compressed between the electrodes,
so the coupling may be stronger and electrons hopping between the electrodes
are also allowed, and hence the measured conductances are higher. The second
scenario can be ascribed to the measurements of the gold atom with the
molecule alongside it, providing greater stability to the junction. The data
in table I and Fig. 6 show that the latter is the most likely scenario.

Figure 8: Experimental conductance-displacement density plots for benzene,
cyclohexane, and toluene are shown in panels a, b, and c, respectively. Inside
each density plot, a coloured rectangle demarcates the area where the traces
of the theoretical model overlap with the experimental data. A zoom-in of the
calculated traces of conductance and an illustration of the molecular contacts
are shown in insets to the right of each density plot. Here, a1 and a2
correspond to benzene, b1 and b2 to cyclohexane and finally c1 and c2, to
toluene.

In Fig. 3 idealized scenarios were shown, whose results can not reproduce
accurately a dynamic process as our experiments. Hence, to identify the
possible binding configurations, we have computed the conductance from some
scenarios obtained via CMD, and compared it with the experimental data using
2D histograms (see Fig. 8). In constructing the experimental 2D histograms we
have included only traces involving molecular junctions under 0.7 G0, aligning
them with the value of 0.5 G0 to enhance solely the molecular contribution.
Overlaid on the 2D histograms, are the computed conductance values obtained as
the systems evolve, i.e. as the distance between the electrodes increases. The
colour maps on the logarithmic scale show the most repeated occurrences of the
plateaus. Rectangular frames facilitate focusing on the range where
theoretical values are obtained. Two insets also accompany each panel to show
the evolution of the junction with the separation distance and the different
atomic arrangements of the molecules between the gold electrodes. The insets’
edges are coloured in concordance with the theoretical trace. The three panels
a,b, and c of Fig. 8, provide insight of the detailed molecular orientations
of the molecules in the junction during the junction break process, with a
high level of agreement between experimental data and the results obtained
from our theoretical model consisting of CMD and DFT+NEGF.

##  IV Conclusions

Through a combination of electronic transport experiments under ambient
conditions, molecular dynamics simulations and DFT transport calculations of
the conductance, we have untangled the relation between electrical conductance
and binding conformations of aromatic and aliphatic molecular junctions.

In the first set of DFT transport calculations we performed, our idealized
DFT-based models reproduced the experimentally observed range of conductance
versus displacement values for molecular junctions made of benzene,
cyclohexane and toluene. Subsequently, classical Molecular Dynamics
simulations aided in the identification and classification of the type of last
contact that is likely formed immediately before rupture in the experiments.
Moreover, CMD simulations generated different scenarios of the contacts that
have allowed us to compute the electronic transport and obtain statistics with
which to make a direct comparison with the experimental data. Consequently,
the results obtained by means of the combination of CMD with DFT are in
complete agreement with the experimental data, as Fig.8 shows. We further
believe that the classification that we have established for the different
types of molecular contacts can be applied in molecular electronics involving
all manner of ring-shaped or planar molecules.

The combination of our previous findingsde Ara  _et al._ (2022) and the
results presented in this manuscript reveal that the solvents benzene,
cyclohexane and toluene are never fully evaporated and remain adsorbed on the
electrodes. Our simulations, calculations and experiments thus allow us to
obtain the characteristic conductance values of these solvents, such that when
they are employed alongside other target molecules in future experiments, we
will be able to clearly distinguish the conductance signatures of the solvents
and their geometric relationships with the electrodes.

##  V Acknowledgment

This work was supported by the Generalitat Valenciana through CIDEXG/2022/45,
CDEIGENT/2018/028, and PROMETEO/2021/017 and the Spanish government through
PID2019-109539-GB-C41. This study forms part of the Advanced Materials program
and was supported by the Spanish MCIN with funding from European Union
NextGenerationEU and by Generalitat Valenciana through MFA/2022/045. The
theoretical modelling was performed on the high-performance computing
facilities of the University of South Africa and the University of Alicante.

## References

  * Agraït  _et al._ (1993) N. Agraït, J. G. Rodrigo, and S. Vieira, Phys. Rev. B 47, 12345 (1993).
  * Krans  _et al._ (1993) J. M. Krans, C. J. Muller, I. K. Yanson, T. C. M. Govaert, R. Hesper, and J. M. van Ruitenbeek, Phys. Rev. B 48, 14721 (1993).
  * Krans  _et al._ (1995) J. M. Krans, J. M. van Ruitenbeek, V. V. Fisun, I. K. Yanson, and L. J. de Jongh, Nature 375, 767 (1995).
  * Krans and van Ruitenbeek (1994) J. M. Krans and J. M. van Ruitenbeek, Phys. Rev. B 50, 17659 (1994).
  * Scheer  _et al._ (1998) E. Scheer, N. Agraït, J. C. Cuevas, A. L. Yeyati, B. Ludoph, A. Martín-Rodero, G. R. Bollinger, J. M. van Ruitenbeek, and C. Urbina, Nature 394, 154 (1998).
  * Smit  _et al._ (2002) R. H. M. Smit, Y. Noat, C. Untiedt, N. D. Lang, M. C. van Hemert, and J. M. van Ruitenbeek, Nature 419, 906 (2002).
  * Venkataraman  _et al._ (2006) L. Venkataraman, J. E. Klare, C. Nuckolls, M. S. Hybertsen, and M. L. Steigerwald, Nature 442, 904 (2006).
  * Li  _et al._ (2022) L. Li, J. Z. Low, J. Wilhelm, G. Liao, S. Gunasekaran, C. R. Prindle, R. L. Starr, D. Golze, C. Nuckolls, M. L. Steigerwald, F. Evers, L. M. Campos, X. Yin, and L. Venkataraman, Nature Chemistry 14, 1061 (2022).
  * van der Zant  _et al._ (2006) H. S. J. van der Zant, E. A. Osorio, M. Poot, and K. O’Neill, physica status solidi (b) 243, 3408 (2006).
  * Reznikova  _et al._ (2021) K. Reznikova, C. Hsu, W. M. Schosser, A. Gallego, K. Beltako, F. Pauly, H. S. J. van der Zant, and M. Mayor, Journal of the American Chemical Society 143, 13944 (2021).
  * Yelin  _et al._ (2013) T. Yelin, R. Vardimon, N. Kuritz, R. Korytár, A. Bagrets, F. Evers, L. Kronik, and O. Tal, Nano Letters 13, 1956 (2013).
  * Taniguchi  _et al._ (2011) M. Taniguchi, M. Tsutsui, R. Mogi, T. Sugawara, Y. Tsuji, K. Yoshizawa, and T. Kawai, Journal of the American Chemical Society 133, 11426 (2011).
  * Yoshida  _et al._ (2015) K. Yoshida, I. V. Pobelov, D. Z. Manrique, T. Pope, G. Mészáros, M. Gulcur, M. R. Bryce, C. J. Lambert, and T. Wandlowski, Scientific Reports 5, 9002 (2015).
  * Tewari  _et al._ (2019) S. Tewari, C. Sabater, and J. van Ruitenbeek, Nanoscale 11, 19462 (2019).
  * Bopp  _et al._ (2017) J. M. Bopp, S. Tewari, C. Sabater, and J. M. van Ruitenbeek, Low Temp. Phys. 43, 905 (2017), https://doi.org/10.1063/1.5001283 .
  * Cuevas and Scheer (2017) J. C. Cuevas and E. Scheer,  _Molecular Electronics_, 2nd ed. (WORLD SCIENTIFIC, 2017) https://www.worldscientific.com/doi/pdf/10.1142/10598 .
  * Agraït  _et al._ (2003) N. Agraït, A. L. Yeyati, and J. M. Van Ruitenbeek, Phys. Rep. 377, 81 (2003).
  * Evers  _et al._ (2020) F. Evers, R. Korytár, S. Tewari, and J. M. van Ruitenbeek, Rev. Mod. Phys. 92, 035001 (2020).
  * de Ara  _et al._ (2022) T. de Ara, C. Sabater, C. Borja-Espinosa, P. Ferrer-Alcaraz, B. C. Baciu, A. Guijarro, and C. Untiedt, Materials Chemistry and Physics 291, 126645 (2022).
  * Yelin  _et al._ (2021) T. Yelin, S. Chakrabarti, A. Vilan, and O. Tal, Nanoscale 13, 18434 (2021).
  * van Veen  _et al._ (2022) F. H. van Veen, L. Ornago, H. S. J. van der Zant, and M. El Abbassi, The Journal of Physical Chemistry C 126, 8801 (2022).
  * Quek  _et al._ (2009) S. Y. Quek, H. J. Choi, S. G. Louie, and J. B. Neaton, Nano Letters 9, 3949 (2009).
  * Venkataraman  _et al._ (2019) A. Venkataraman, P. Zhang, and C. Papadopoulos, AIP Advances 9, 035122 (2019), https://doi.org/10.1063/1.5087413 .
  * Pan  _et al._ (2022) X. Pan, C. Qian, A. Chow, L. Wang, and M. Kamenetska, The Journal of Chemical Physics 157, 234201 (2022), https://doi.org/10.1063/5.0103642 .
  * Louis  _et al._ (2003) E. Louis, J. A. Vergés, J. J. Palacios, A. J. Pérez-Jiménez, and E. SanFabián, Phys. Rev. B 67, 155321 (2003).
  * Palacios  _et al._ (2001) J. J. Palacios, A. J. Pérez-Jiménez, E. Louis, and J. A. Vergés, Phys. Rev. B 64, 115411 (2001).
  * Palacios  _et al._ (2002) J. J. Palacios, A. J. Pérez-Jiménez, E. Louis, E. SanFabián, and J. A. Vergés, Phys. Rev. B 66, 035322 (2002).
  * (28) J. J. Palacios _et al._ , “Computer code ANT.Gaussian, with SOC corrections,” Available from https://github.com/juanjosepalacios/ANT.Gaussian.
  * Zotti  _et al._ (2022) L. A. Zotti, W. Dednam, E. B. Lombardi, and J. J. Palacios, Nanomaterials 12 (2022), 10.3390/nano12071234.
  * Wadt and Hay (1985) W. R. Wadt and P. J. Hay, The Journal of Chemical Physics 82, 284 (1985), https://doi.org/10.1063/1.448800 .
  * Ross  _et al._ (1990) R. B. Ross, J. M. Powers, T. Atashroo, W. C. Ermler, L. A. LaJohn, and P. A. Christiansen, The Journal of Chemical Physics 93, 6654 (1990), https://doi.org/10.1063/1.458934 .
  * Allen and Tildesley (1989) M. P. Allen and D. J. Tildesley,  _Computer simulation of liquids_ (Oxford university press, 1989).
  * Frenkel and Smit (2002) D. Frenkel and B. Smit,  _Understanding Molecular Simulation: From Algorithms to Applications_ (Academic Press, 2002).
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * (35) S. Plimpton _et al._ , ”computer code LAMMPS, publicly available at http://lammps.sandia.gov”.
  * Järvi  _et al._ (2011) T. T. Järvi, A. C. T. van Duin, K. Nordlund, and W. A. Goddard, The Journal of Physical Chemistry A 115, 10315 (2011).
  * Senftle  _et al._ (2016) T. P. Senftle, S. Hong, M. M. Islam, S. B. Kylasa, Y. Zheng, Y. K. Shin, C. Junkermeier, R. Engel-Herbert, M. J. Janik, H. M. Aktulga, T. Verstraelen, A. Grama, and A. C. T. van Duin, npj Computational Materials 2, 15011 (2016).
  * Soria  _et al._ (2018) F. A. Soria, W. Zhang, P. A. Paredes-Olivera, A. C. T. van Duin, and E. M. Patrito, The Journal of Physical Chemistry C 122, 23515 (2018).
  * Gomzi  _et al._ (2021) V. Gomzi, I. M. Šapić, and A. Vidak, The Journal of Physical Chemistry A 125, 10649 (2021), pMID: 34883013, https://doi.org/10.1021/acs.jpca.1c06939 .
  * Sabater  _et al._ (2018) C. Sabater, W. Dednam, M. R. Calvo, M. A. Fernández, C. Untiedt, and M. J. Caturla, Phys. Rev. B 97, 075418 (2018).
  * Goodfellow (city) Goodfellow, (description and publicity), ”metal trading supplier, publicly available at https://www.goodfellow.com”.

##  VI Supplementary information

Figure S1: Electronic transport vs relative displacement calculated for the
initial structures in panel (a) benzene, (b) cyclohexane and (c) toluene.
Coloured markers represent the calculated values of gold (orange), benzene
(black), cyclohexane (pink) and toluene (green).

Fig.S1 shows the evolution of the calculated conductance vs the relative
displacement in the range $0$ Å to $6$ Å. This figure is the same as Fig. 3
but without the subtraction of the tunnelling regime. In all the cases we have
moved the upper electrode in step intervals of $0.2$ Å without relaxation and
we have calculated the conductance by DFT+NEGF. The three panels show the
conductance in units of $G_{0}$ for benzene, cyclohexane and toluene molecules
in parallel and perpendicular configurations. The red lines indicate the slope
of the tunnelling regime. Yellow areas highlight the first value that stat to
be separated from the red trend line. We used the value of the relative
displacement highlighted to create Fig. 3.