Modeling of Solvent Effects in the Electrical ... - ACS Publications

Article pubs.acs.org/JPCC

Modeling of Solvent Effects in the Electrical Response of π‑Stacked Molecular Junctions Tahereh Ghane,† Andrii Kleshchonok,† Rafael Gutierrez,*,† and Gianaurelio Cuniberti†,‡,§ Institute for Materials Science, ‡Dresden Center for Computational Materials Science, and §Center for Advancing Electronics Dresden, Dresden University of Technology, 01062 Dresden, Germany

Downloaded by UNIV OF MANITOBA on August 30, 2015 | http://pubs.acs.org Publication Date (Web): August 19, 2015 | doi: 10.1021/acs.jpcc.5b06867

†

ABSTRACT: We present theoretical modeling of the influence of THF solvents on the mechanical stability and the electrical response of two different π-stacked molecular junctions based on cysteamine conjugates of naphthalic anhydride and of pyrene. Combining molecular dynamics simulations and quantum transport calculations, we show that for junctions with a weaker π−π stackingas measured by the stacking energydynamical breaking of the stacking induced by the solvent can take place. However, contrary to what may be expected, the conductance of the system is not suppressed due to the emergence of an additional transport channel which bypasses the broken π overlap of the perylene cores. However, an additional gating-like effect in such a situation does reduce the low bias current when comparing with situations, where π-stacking is preserved.

■

INTRODUCTION Over the past two decades, electronic-transport properties of single molecular junctions have extensively been investigated both experimentally and theoretically in view of potential applications in nano- and subnanoelectronic devices.1−7 In many instances, measurements of the electrical response of molecular junctions are carried out at room temperature in nonconductive solvents,8−12 which can potentially influence the junction formation as well as its conductance by influencing the mechanically stable conformations of the molecules.13 To mention only one example, using a scanning tunneling microscope based point-contact technique, electrical transport measurements of 1,4-benzenediamine molecular junctions have been carried out in 13 different ambient-solvated environments, showing that the conductance of the junction could vary by more than 50% depending on the solvent.13 Clarifying the impact of solvents on the molecular junction conductance is thus important because it exemplifies the potential of the environment in controlling the electronic response of nanoelectronic devices. Theoretical works that have been carried out on the impact of the solvent on the conductance of molecular junctions usually address the influence of the solvent on the static14 and on the dynamical15−18 molecular junction geometries, the major focus lying in this latter case on aqueous environments. For example, ref 15 addressed the effects of water solvent on the transport properties of molecular nanojunctions, which resulted in a transmission spectrum shifted by about 0.6 eV with respect to that of the dry junction. This confirms the fact that using polar solvents like H2O can © 2015 American Chemical Society

significantly affect the transport properties because of the generation of local dipole fields. Recently, Kotiuga and co-workers19 developed a quantitative understanding and a general model of the effects of the solvent on the conductance of singlemolecule junctions. They demonstrated that an electrostatic model approximating the junction and its surroundings by an array of point-polarizable dipoles quantitatively captures these effects and can be extended to incorporate the influence of nearby solvent molecules on the local potential and the conductance.19 The use of tetrahydrofuran (THF) as a solvent has been recently extended into the area of nanoelectronics. In nonpolar solvents such as liquid Xe,20 liquid methane and low density liquid ethane21 excess electrons do not absorb in the visible or near-IR, which is consistent with the simulation picture of a highly delocalized ground-state wave function and no bound electronic excited states.22 However, in weakly polar fluids such as tetrahydrofuran (THF) the solvated electron has an absorption spectrum, but it is quite different from that in polar fluids. Results presented in ref 23 show that the peak of the solvated electron’s spectrum in THF occurs at much lower energies and is substantially broader than the solvated electron’s spectrum in water. For this reason, THF is considered as a proper solvent in electronic transports measurements. Received: July 16, 2015 Revised: August 12, 2015 Published: August 12, 2015 20201

DOI: 10.1021/acs.jpcc.5b06867 J. Phys. Chem. C 2015, 119, 20201−20209

Article

The Journal of Physical Chemistry C

Concerning the THF solvent, we have used force field parameters previusly reported in the literature.28 We chose the united atom model by Helfrich and Hentschke29 because this model performs best for simulating the THF molecules on Au(111) surfaces.28 In this model, Lennard-Jones and electrostatic interactions in the THF molecule are carried by the oxygen atom and the four carbon atoms, while the hydrogen atoms are only implicitly included. THF−THF and THF−molecule nonbonding Lennard-Jones interactions are described as


In this paper, we report calculations of the charge transport properties of π-stacked molecular systems sandwiched between Au(111) electrodes and embedded in a THF solvent. In particular, we focus on the effect of the THF solvent on the transmission function of the molecular junctions. We aimed at addressing the question of how the presence of the solvent influences the electrical transport signatures of the junction by dynamically modifying the accessible conformational states of the molecule. We deal with this problem by combining classical molecular dynamics simulations, electronic structure calculations, and quantum transport calculations. In particular, we show that solvent-induced breaking of the π-stacking does not necessarily suppress transport through the junctions due to the emergence of alternative charge transport pathways. This study partially relies on a previous joint experimental−theoretical work,24 where π-stacked molecular junctions were studied within a multiscale approach in vacuum to address their relevance for tuning charge transport in nanoparticle networks.

UijLJ

⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σij σij = 4ϵij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠

(1)

where i and j located on different molecules. The cross interaction potential parameters σij and ϵij are calculated from the Lorentz−Berthelot combining rules, if interacting sites i and j are on different species. The configurational information and the potential parameters of the rigid THF model are summarized in Table 1.

■

COMPUTATIONAL METHODOLOGIES A. Geometry Optimization: First-Principles Calculations. We have addressed two different molecules building π-stacked dimers: a cysteamine conjugate of naphthalic anhydride (CYS-NA) and a cysteamine conjugate of pyrene (CYS-PA).24 The corresponding molecular junctions consist of pairs of molecules of the same type with overlapping perylene cores (π-stacking) and attached on each side to metallic electrodes via thiol linkers. The side chains of the CYS-NA and CYS-PA molecules have different lengths: in the CYS-NA case the side chain contains two methylene groups, while in the CYS-PA case there are three methylene groups interrupted by an amino group. At the end of the side chains there is a thiol linker, which binds to the gold electrodes. Charge transport signatures of theses systems in vacuum were studied in a previous publication,24 but a more realistic setup mimicking mechanically controllable break junction experiments in solution needs to include solvent effects at least in an approximate way. Density-functional theory (DFT) calculations were first used for the geometry optimization and electronic structure calculations of isolated CYS-PA and CYS-NA molecules and π-stacked dimers chemisorbed on Au(111) surfaces via thiol linkers. Hereby, we have used combined plane-wave and atomic-orbital-based approaches as implemented in the cp2k code.25 Periodic boundary conditions were employed along all three unit cell directions and gold slabs three layers thick were used to model the Au(111) surfaces. DFT calculations used the generalized gradient approximation (GGA) and the PBE exchange-correlation energy functionals.26 The van der Waals (vdW) interactions were further taken into account through the standard Grimme27 parametrization. We carried out full geometry optimizations, in which all atoms of the molecules and the first layer of the Au(111) surface were allowed to relax. The relaxed geometries served as starting conformations for the classical molecular dynamics simulations described in the next section. B. Classical Molecular Dynamics. To describe the structural fluctuations of the π-stacked molecules through force field (FF) based classical molecular dynamics (MD) simulations, we have validated and improved in a previous study24 the OPLS/AA force field parameters for CYS-PA and CYS-NA molecular dimers in a Au−molecule−Au junction to properly describe the π-stacking interactions. We thus refer the reader to this study for details of the force field parametrization.

Table 1. Molecular Configuration and Potential Parameters for the THF Structure28 vdW carbon group O bond stretching

ϵ(kcal/mol)

r0 (Å)

0.090 0.200

4.300 3.200 r0 (Å)

C−O C−C angles

1.432 1.523 θ0 (deg)

C−O−C O−C−C C−C−C

109.6 108.0 101.4

The MD simulations of π-stacked single molecular junctions in a vacuum and in the presence of THF solvent were based on the Gromacs package.30 The simulations were carried out in the NVT ensemble, using the Nosé−Hoover thermostat to maintain the temperature at 300 K. Three-dimensional periodic boundary conditions were applied and Newton’s equations of motion were solved using the Leapfrog algorithm,31 with a time step of 1 fs. Particle mesh Ewald (PME) electrostatic summation was used with a real space cutoff at 11 Å, whereas a force-switched cutoff starting at 9 Å and ending at 10 Å was used for Lennard-Jones nonbonded interactions. The simulations in a vacuum have been reported in ref 24, and here we present the results of calculations in the THF solvent. For these simulations, a gold slab with a surface area of 48 × 55 Å2 and five atomic layers thick was used to simulate the metallic substrate. The size of the gold surface gives a gap of approximately 55 Å between the adsorbed molecules and the nearest periodic image of the slab to avoid the interaction of molecules with their images. The space between the two gold surfaces in the junction was filled with THF solvent molecules. To address the effect of the solvent concentration, we prepared different boxes with different concentrations. The solvent boxes were minimized using steepest descent until an energy gradient of 10 kJ mol−1 nm−1 was reached, and then the boxes were added to the junction. The solvated junction was minimized again by steepest descent algorithm with solute atoms harmonically constrained, followed by the minimization of the whole system using the same algorithm. To let the molecules 20202


Article



From the single-point calculations with the DFTB code, we obtain the Hamiltonian H and the overlap S matrix elements needed to build the necessary quantities in the quantum transport formalism. Since we assume coherent charge transport through the system (i.e., no dynamical coupling to vibrational degrees of freedom), we use the Landauer scheme to compute the charge transport characteristics.35 Here, the key quantity is the energy-dependent quantum mechanical transmission function T(E), which quantifies the degree of ”transparency” of the system to a charge propagating from one electrode to the other. The transmission can be computed using Green’s functions and is given by the following expression:

relaxing, short annealing MD simulations of up to 40 ps were first carried out. The temperature of the gold substrate was maintained at 300 K throughout, while that of the molecules and solvent was incrementally decreased by 1 K. The starting temperature (200 K) decreased to 100 and 1 K in 10 and 30 ps, respectively. After that, the system was allowed to equilibrate for 10 ps more. The MD simulations were then performed over 3 ns with snapshots saved every 0.5 ps. C. Transport Calculations. In general terms, when a charge is injected into a molecule, it can, during its propagation, modify the local electronic structure as well as the local atomic arrangement, leading in some cases to effects like polaron formation.32 Whether such effects are important or not depends on different parameters such as the strength of the electron− vibrational coupling, the length of the molecular region, and the ratio between the typical leakage time into the electrodes (roughly given by the inverse of the dissipation rate γ, see below) and typical time scales of molecular vibrations τvibr. In our study, we will assume an adiabatic-like approximation, where the basic assumption is that ℏγ−1 ≪ τvibr; i.e., the electronic coupling to the electrodes is large enough, so that a charge will only probe static molecular conformations during the tunneling process. In our case, this assumption is approximately valid, since the coupling between the thiol linkers and the gold atom to which it is covalently attached is strong (short leakage times). We might expect deviations for specific conformations where the π-stacking is fully broken, since charge localization effects could become important. However, as we will show later on, additional transport channels emerge in such cases reestablishing a rather large transmission. We therefore compute the transmission function along the MD trajectory without taking into account the potential back-reaction of the structural fluctuations on the charge. As a result, a time series of transmission functions is generated at each energy, encoding the global influence of the structural dynamics as well as the solvent effects We also notice that the influence of the solvent molecules manifests mostly in steric effects on the π-stacked molecules, indirectly modifying the transport signatures by changing the conformational space explored by the molecules during the MD simulation. This of course influences the molecular electronic structure and thus the transmission spectrum. In fact, at least in one of the studied junctions (the CYS-PA one), a gating-like effect of the solvent is found, leading to a rigid shift of the time-averaged transmission function. We also remark that charge transport is considered to be taking place along the molecular dimer and not through pathways leading through the solvent. Although the static junction geometries were structurally relaxed with a full first-principles approach, the quantum transport calculations for each snapshot extracted from the classical MD simulation were carried out with a density-functional parametrized tight-binding (DFTB) approach,33,34 since computing the transmission function only involved a single-point calculation of the electronic structure with no further relaxation. The DFTB code uses a minimal valence basis set to build a generalized (i.e., beyond nearest neighbors) nonorthogonal tight-binding Hamiltonian, with electronic matrix elements being precomputed according to a meanwhile well-established procedure.33,34 This approach is computationally very efficient for calculating the electronic structure of complex systems and allows for computations with up to 2000 atoms with a reasonable CPU time. This is especially relevant for charge transport calculations in molecular junctions, where not only the molecular system needs to be explicitly considered but also the mesoscopic electrodes.

T (E) = Tr(ΓL(E)Gr (E)ΓR(E)Ga(E))

(2)

Here, Gr(a) is the retarded (advanced) Green’s function of the extended molecular system,36,37 which is built in our case by including not only the molecular π-stacked dimer as such but also the first two layers of the gold electrodes. This ensures that all elastic scattering processes taking place in the system are appropriately accounted for. The interaction with the left and right leads is taken into account by means of the so-called self-energy functions ΣrL(R)(E), which are related to the spectral densities ΓL(R)(E) = i(ΣrL(R)(E) − ΣaL(R)(E)). The self-energies encode both the surface electronic structure of the electrodes and the electronic coupling between atoms on the electrode surface and atoms in the molecular system, in this particular case the Au−S interaction. The self-energies account for the fact that we are dealing with an open quantum system; they can be numerically computed in a standard way by using the iterative Lopez−Sancho procedure.38 The retarded (and also the advanced) Green’s function Gr satisfies a Dyson equation in the energy domain, given by [Gr(E)]−1 = (E + iη)S − H − ΣrL(E) − ΣrR(E), where clearly all quantities are matrices with the dimension of the corresponding electronic space of the extended molecule. Whenever required, the corresponding electrical current can be computed as I(V) = (2e/h)∫ dE T(E)( f L − f R), where f L and f R are the Fermi−Dirac distribution functions of the left and right electrodes with symmetrically applied bias ±eV/2. As mentioned above, both

Figure 1. Two studied π-stacked molecular junctions: (a) CYS-PA and (b) CYS-NA π-stacked dimers sandwiched between two gold electrodes. 20203


Article



Figure 2. Snapshots from a molecular dynamics run for the CYS-PA (top panel, a, b) and CYS-NA (low panel, c, d) molecular junctions for low concentration (a, c) and high concentration (b, d) THF boxes. In both cases, panels display molecular conformations taken at 3 ns. The CYS-PA molecule is presented in blue, and the CYS-NA molecule is presented in green.

Figure 3. Short-range vdW interaction for CYS-NA and CYS-PA molecular junctions in the low concentration box (upper panel) and high concentration box (lower panel). This analysis illustrates the strong THF−THF interaction when compared to THF−Au.

the transmission function and the current get an additional time labeling when computed along the MD trajectory, so that also time-average quantities will be computed later on. Notice that the spectral functions quantify the dissipation rate into the electrodes, and their matrix elements are related to the time scale γ mentioned at the beginning of this section.

■

RESULTS AND DISCUSSION Classical Molecular Dynamics Simulations. THF solvents have been used in experimental studies of alkanethiols on a Au surface.39,40 In this study we focus on the effect of the THF solvent on the dynamics of the π-stacked molecular junctions shown in Figure 1a,b. We have prepared three solvent boxes with different concentrations of THF molecules: 2 molecules/nm2 (model I), 3 molecules/nm2 (model II), and 4 molecules/nm2 (model III). After performing first trial MD

Figure 4. Upper panel: radius of gyration Rg calculated around the center of mass of the (a) CYS-PA and (b) CYS-NA molecules during 3 ns in two different boxes. The curves are smoothed over the range of 100 ps. Lower panel: statistics of the angle between the normal vectors of the two aromatic rings in (c) CYS-PA and (d) CYS-NA π-stacked molecules for the same simulation scenarios as in the upper panel. 20204


Article


between the two Au electrodes. Cluster formation is stronger in the case of the low concentration box (model I). Our simulations with different concentrations showed that because of the restricted space between the Au electrodes the concentration of THF solvent has an important effect on the formation of clusters. As Figure 2 shows, there are no distinct boundaries between the THF layers in both cases. Indeed, this finding has been confirmed by Monte Carlo simulations.28 Our energy analysis shows that in both cases THF−THF interactions are stronger than THF−Au interactions, which clarifies why the solvent molecules are clustering. Figure 3 shows the short-range van der Waals interactions for both systems in the low concentration and high concentration models. In all cases, the interaction between THF molecules is stronger than the interactions of the THF with the surface gold atoms. Molecular Dimers. Our analysis of the molecular linkers shows that the impact of the THF solvent on the CYS-PA molecule varies by changing the solvent concentration while solvent effects are considerably weaker in the CYS-NA junction. This is summarized in Figure 4. Figure 4a,b shows the radius of gyration Rg of the linkers around their center of mass during the simulation time. In the CYS-PA dimer, Figure 4a, the conformational changes of the molecule during the 3 ns dynamics are highly affected by the number of THF molecules. On the contrary, in the CYS-NA junction, Figure 4b, there is no noticeable difference in the dynamics of the molecule in the two different solvent models. Further analysis, presented in Figure 4c,d,


simulations and analyzing the corresponding MD trajectories, we found that the impact of the solvent on the conformational dynamics of the π-stacked molecules in cases I and II turned out to be rather similar. For this reason, we limit the following discussion to cases I and III. We divide the trajectory analysis into two parts: analysis of the solvent and analysis of the π-stacked molecules. THF Solvent. Figure 2a−d presents snapshots of the THF solvent for two different concentrations in the two studied molecular junctions. After a 3 ns MD run, we found that the THF molecules stick together rather than spreading uniformly

Figure 5. Radial distribution function (RDF) of THF solvent around the center of mass of (a) CYS-PA and (b) CYS-NA π-stacked molecules during a 3 ns MD run in two different solvent boxes.

Figure 6. Upper panel: transmission function through the CYS-NA molecular junction with (a) low (2 molecules/nm2) and (b) high (4 molecules/nm2) THF solvent concentration. Lower panel: similar transmission function plot for the CYS-PA junction with (c) low (2 molecules/nm2) and (d) high (4 molecules/nm2) THF solvent concentration. In all cases, the red curve shows the time-averaged transmission, while the bunches of blue curves are snapshots of the transmission function along the MD trajectory (shown are 81 snapshots in the time window between 1 and 2 ns), taken from the MD simulations. The insets show in all case typical snapshots within the chosen time window. The π-stacking is preserved for low solvent concentrations (a, c) as well as for the smaller molecule CYS-NA in higher solvent concentration (b). However, for the CYS-PA case, an increase of the THF solvent concentration leads to a breaking of the π-stacking (d). Nevertheless, the interaction between carbon and nitrogen atoms provides an alternative charge transport pathway. 20205


Article



larger, a fact related to a solvent-induced suppression of too strong conformational fluctuations in the current study. In this sense, the THF solvent acts as a stabilizing factor for the CYS-NA molecular dimer. The situation changes drastically when considering now the larger CYS-PA molecular dimer (see panels c and d in Figure 6). In the case of a low concentration solvent (Figure 6c), the π-stacking is preserved, and the transmission curves are qualitatively similar to the previous case of the CYS-NA junction, although obviously different spectral features arise. In particular, the (nonresonant) transmission within the HOMO−LUMO gap is effectively increased (Figure 6c) due to the presence of the broad spectral feature within the energy window ranging from −1.8 to −0.8 eV below the Fermi level. The tails of transmission resonances in this energy window enter the molecular HOMO−LUMO gap and increase the transmssion around the Fermi energy. However, if we increase the THF solvent concentration up to 4 molecules/nm2, the π-stacking is dynamically broken; i.e., there exist a large set of statistically relevant conformations where the π−π overlap of the two CYS-PA molecules in the dimer is strongly perturbed (see the inset of Figure 6d), which shows a typical conformation belonging to this class. In this case, one may naively expect that the conductance of the system should be suppressed over a broad energy range. Indeed, we find specific conformations where charge transport is fully suppressed due to the decoupling of the two CYS-PA molecules in the dimer. There is however a large number of trajectories, where despite an apparent breaking of the π-stacking, charge transport still takes place and a sizable conductance is found (see the transmission curves in Figure 6d). To understand this

reveals that in the CYS-PA junction enhancing the solvent concentration has a major effect on the π−π stacking. Here, Figure 4c,d displays the probability of finding the normal vector of one plane of carbon rings at a specific angle with respect to the normal vector of the second plane in the π-stacked structure in such a way that 0 and (±)90 correspond to perfect stacking and breaking of the stack, respectively. The results clearly illustrate the imapct of the solvent on the molecular conformational dynamics. As we can see, the original conformation (π-stacking) of the CYS-PA dimer is still preserved at low concentrations (model I), while the range of angle fluctuations of the aromatic rings in the high concentration box (model III) is considerably wider. Thus, the dynamics of the CYS-PA junction clearly displays solvent-dependent features, while the structural dynamics of the CYS-NA junction is considerably less affected by the presence of the solvent. For reference, the calculated stacking energy using the previously discussed classical force fields are −188.3 and −80.75 kJ/mol for CYS-NA and CYS-PA, respectively. This indicates that the π−π stacking in the latter dimer is considerably weaker and can thus be stronger perturbed in a higher concentration solvent. To further study the behavior of the solvent around the molecule, we additionally computed the radial distribution function of THF molecules around the center of mass of the linker (see Figure 5a,b). In the case of the CYS-PA dimer, shown in Figure 5a, there are some peaks very close to the center of mass which are due to the breaking of the molecular stack in the high concentration model. Electronic Structure and Transport Properties. To analyze the electronic transport properties of the system in the presence of the solvent, we computed the quantum mechanical transmission function along the lines described in section C for the CYS-PA and CYS-NA molecular dimers and for the cases of low (2 molecules/nm2, model I) and high (4 molecules/nm2, model III) solvent concentrations. To account for the influence of structural fluctuations, the transmission was calculated at 80 snapshots within a simulation time window between 1 and 2 ns, where we found the largest changes in the structural properties of the junctions. The main results are summarized in Figure 6a−d, where we plot the quantum mechanical transmission functions at different snapshots as well as their time average, as a function of the electron’s injection energy. Although the energy window relevant for linear transport lies around the Fermi energy of the system, we plot a broader energy range to better illustrate the impact of conformational fluctuations at different energies (both at resonance and within transmission gaps). As previously mentioned, the π−π stacking in the CYS-NA molecule is preserved during the MD simulation for both concentrations of the THF solvent. This behavior is clearly reflected in the transmission curves for this molecule, as shown in Figure 6a,b. For both solvent concentrations, models I and III, the π-stacking is well preserved, so that the time average transmission remains nearly the same. Notice also that the Fermi level is very close to a broad set of states corresponding to the unoccupied (LUMO) states of the CYS-NA dimer. Although the transmission at the Fermi energy is rather low, small variations in the position of the Fermi level could lead to a strong increase of the conductance by moving the LUMO states closer to EF. In general, the fluctuations around the time average can be rather large, covering e.g. 3 orders of magnitude within the HOMO−LUMO gap. A similar behavior has been previously obtained for simulations of the same molecule in a vacuum;24 the range of the fluctuations in the transmission was however

Figure 7. Top panel: a typical snapshot of the CYS-PA junction showing the broken π-stacking. The region encircled includes the nitrogen atom on one molecule and the closest carbon of the other one, at a distance of about 3.4 Å. Through this link, a transport pathways opens and rather large values of the transmission can be obtained. Lower panel: total and projected (on the C−N region) density of states (in arbitrary units) of the junction as well as the transmission function for the specific conformation shown in the upper panel. Notice that the C−N region appreciably contributes to a broad spectral range, thus explaining the high transmission despite the broken π-stacking. 20206


Article



Figure 8. Time-averaged transmission functions for the (a) CYS-PA and (b) CYS-NA molecular junctions with low solvent concentration (2 molecules/nm2, green curve) and high solvent concentration (4 molecules/nm2, magenta curve). In case (a) there is a gating-like effect leading ot a horizontal shift of the transmission, while in case (b) the two curves are almost identical.

consequent redistribution of charges within the molecules; as a result, the transmission function around the Fermi energy (±1 eV) is drastically changed when comparing to the low concentration case. In the CYS-NA dimer, on the contrary, the π-stacking is preserved, and the transmission curves for the two solvent concentrations are basically identical. We also plot in the insets of Figure 8a,b the calculated time-averaged I−V characteristics, which are experimentally measurable quantities. Here, the gating-like effect seen in the average transmission of the CYS-PA dimer manifests as a suppression of the current up to 2.5 V, after which conducting states enter the integration window. It is interesting to note that the solvent effect is clearly revealed in the I−V response, which encompasses both the gating effect and the specific energy-dependent features of the transmission. We finally remark that these results for the I−V characteristics are supposed to only display qualitative trends; a more complete calculation would require a full nonequilibrium treatment, including the bias dependence of the transmission as well as, possibly, the influence of the electric field in the junction on the solvent.

behavior, we have analyzed the projected electronic density of states, which is shown in Figure 7. First of all, we remark that there are configurations of the dimer where the distance between a C atom of one molecule and the N atom of the other one is ∼3.4 Å, which is comparable with the obtained typical distance of 3.8 Å, between the aromatic rings in the π-stacking. The relevant region is highlighted in the upper panel of Figure 7. The lower panel of Figure 7 shows the electronic density of states (DOS) for the solvent model III, projected on different atomic regions of the system: total DOS (red) and projected DOS on the molecular region (green) and on the C−N region (brown). For reference, the quantum mechanical transmission for the single structure (single snapshot) shown in the upper panel of Figure 7 is also included. First of all, a rather good correlation between the energetic position of the transmission peaks and the molecular states in the DOS can be clearly seen. More important is, however, the fact that the C−N pair contributes to the total DOS over a broad spectral range. This result explains the high transmission found in this junction despite the break of the π-stacking, since now an alternative, equally efficient electron transport pathway opens. In Figure 8a,b we compare the time-averaged transmission functions for the two molecular junctions and the two solvent models I and III. The average transmission in the CYS-PA junction (panel a) in the high concentration case is horizontally shifted (gated), a result related to the broken π-stacking and the

■

CONCLUSIONS In summary, we have investigated, by a combination of electronic structure calculations, classical molecular dynamics, and quantum transport methods, the influence of THF solvent on the electrical transport signatures of two π-stacked molecular 20207


Article


(11) Mishchenko, A.; Zotti, L. A.; Vonlanthen, D.; Burkle, M.; Pauly, F.; Cuevas, J. C.; Mayor, M.; Wandlowski, T. Single-Molecule Junctions Based on Nitrile-Terminated Biphenyls: A Promising New Anchoring Group. J. Am. Chem. Soc. 2011, 133, 184−187. (12) Bernard, L.; Kamdzhilov, Y.; Calame, M.; van der Molen, S. J.; Liao, J.; Schönenberger, C. Spectroscopy of Molecular Junction Networks Obtained by Place Exchange in 2D Nanoparticle Arrays. J. Phys. Chem. C 2007, 111, 18445−18450. (13) Fatemi, V.; Kamenetska, M.; Neaton, J. B.; Venkataraman, L. Environmental Control of Single-Molecule Junction Transport. Nano Lett. 2011, 11, 1988−1992. (14) Baldea, I.; Koppel, H.; Wenzel, W. Applying the Extended Molecule Approach to Correlated Electron Transport: Important Insight from Model Calculations. Phys. Chem. Chem. Phys. 2013, 15, 1918−1928. (15) Rungger, I.; Chen, X.; Schwingenschlogl, U.; Sanvito, S. FiniteBias Electronic Transport of Molecules in a Water Solution. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 235407−235415. (16) Tawara, A.; Tada, T.; Watanabe, S. Electrostatic and Dynamical Effects of an Aqueous Solution on the Zero-Bias Conductance of a Single Molecule: A First-principles Study. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 073409−073412. (17) Leary, E.; Hobenreich, H.; Higgins, S. J.; van Zalinge, H.; Haiss, W.; Nichols, R. J.; Finch, C. M.; Grace, I.; Lambert, C. J.; McGrath, R.; et al. Single-Molecule Solvation-Shell Sensing. Phys. Rev. Lett. 2009, 102, 086801−086804. (18) Cao, H.; Jiang, J.; Ma, J.; Luo, Y. Temperature-Dependent Statistical Behavior of Single Molecular Conductance in Aqueous Solution. J. Am. Chem. Soc. 2008, 130, 6674−6675. (19) Kotiuga, M.; Darancet, P.; Arroyo, C. R.; Venkataraman, L.; Neaton, J. B. Adsorption-Induced Solvent-Based Electrostatic Gating of Charge Transport through Molecular Junctions. Nano Lett. 2015, 15, 4498−4503. (20) Gallicchio, E.; Berne, B. J. On the Calculation of Dynamical Properties of Solvated Electrons by Maximum Entropy Analytic Continuation of Path Integral Monte Carlo Data. J. Chem. Phys. 1996, 105, 7064−7078. (21) Liu, Z.; Berne, B. J. Electron Solvation in Methane and Ethane. J. Chem. Phys. 1993, 99, 9054−9069. (22) Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. The Role of Solvent Structure in the Absorption Spectrum of Solvated Electrons: Mixed Quantum/Classical Simulations in Tetrahydrofuran. J. Chem. Phys. 2005, 122, 134506−1345016. (23) Gorti, S.; Plank, L.; Ware, B. R. Determination of Electrolyte Friction from Measurements of the Tracer Diffusion Coefficients, Mutual Diffusion Coefficients, and Electrophoretic Mobilities of Charged Spheres. J. Chem. Phys. 1984, 81, 909−914. (24) Ghane, T.; Nozaki, D.; Dianat, A.; Vladyka, A.; Gutierrez, R.; Chinta, J. P.; Yitzchaik, S.; Calame, M.; Cuniberti, G. Interplay between Mechanical and Electronic Degrees of Freedom in π-Stacked Molecular Junctions: From Single Molecules to Mesoscopic Nanoparticle Networks. J. Phys. Chem. C 2015, 119, 6344−6355. (25) CP2K: Open Source Molecular Dynamics (http://www.cp2k. org/), accessed August 04, 2015. (26) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (27) Grimme, S. Semiempirical GGA-type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (28) Zhao, X.; Leng, Y.; Cummings, P. T. Self-Assembly of 1,4Benzenedithiolate/Tetrahydrofuran on a Gold Surface: A Monte Carlo Simulation Study. Langmuir 2006, 22, 4116−4124. (29) Helfrich, J.; Hentschke, R. Molecular Dynamics Simulation of Macromolecular Interactions in Solution: Poly(.Gamma.-Benzyl Glutamate) in Dimethylformamide and Tetrahydrofuran. Macromolecules 1995, 28, 3831−3841. (30) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718.

junctions. We have found that the role played by the solvent can be completely different in dependence of its concentration but also on the molecular structure. Thus, while for low solvent concentrations, the π-stacking is as a rule, well preserved, increasing the solvent concentration can induce in the junction with a lower stacking energy a dynamical break of the π-stacking and a gating-like effect manifesting in a shift of the conductance along the energy axis. This, in turns, leads to clear differences in the finite bias electrical response of the junction. However, contrary to what it might have been expected, breaking of the π-stacking does not automatically imply full suppression of the charge transport through the system. Indeed, we show that additional transport pathways can emerge, which give a rather high transmission.

■

AUTHOR INFORMATION


Corresponding Author

*Tel +49 (0)351 46331419; e-mail rafael.gutierrez@ nano.tu-dresden.de (R.G.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank A. Dianat, D. Nozaki, M. Calame, and S. Yitzchaik for very fruitful discussions. This work was partly funded by the EU within the projects Synaptic Molecular Networks for Bioinspired Information Processing (SYMONE, project no. 318597). This work has also been partly supported by the German Research Foundation (DFG) within the Cluster of Excellence “Center for Advancing Electronics Dresden”. We acknowledge the Center for Information Services and High Performance Computing (ZIH) at TU Dresde for computational resources.

■

REFERENCES

(1) Reed, M. A.; Zhou, C.; Muller, C. J.; Burgin, T. P.; Tour, J. M. Conductance of a Molecular Junction. Science 1997, 278, 252−254. (2) Xiao, X.; Xu, B.; Tao, N. J. Measurement of Single Molecule Conductance: Benzenedithiol and Benzenedimethanethiol. Nano Lett. 2004, 4, 267−271. (3) Aradhya, S. V.; Venkataraman, L. Single-Molecule Junctions beyond Electronic Transport. Nat. Nanotechnol. 2013, 8, 399−410. (4) Natelson, D. Mechanical Break Junctions: Enormous Information in a Nanoscale Package. ACS Nano 2012, 6, 2871−2876. (5) Nitzan, A.; Ratner, M. A. Electron Transport in Molecular Wire Junctions. Science 2003, 300, 1384−1389. (6) Gschneidtner, T. A.; Diaz Fernandez, Y. A.; Moth-Poulsen, K. Progress in Self-Assembled Single-Molecule Electronic Devices. J. Mater. Chem. C 2013, 1, 7127−7133. (7) Kondo, H.; Kino, H.; Nara, J.; Ozaki, T.; Ohno, T. ContactStructure Dependence of Transport Properties of a Single Organic Molecule between Au Electrodes. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 235323−235332. (8) Xu, B.; Tao, N. J. Measurement of Single-Molecule Resistance by Repeated Formation of Molecular Junctions. Science 2003, 301, 1221− 1223. (9) Venkataraman, L.; Klare, J. E.; Nuckolls, C.; Hybertsen, M. S.; Steigerwald, M. L. Dependence of Single-Molecule Junction Conductance on Molecular Conformation. Nature 2006, 442, 904− 907. (10) Hybertsen, M. S.; Venkataraman, L.; Klare, J. E.; Whalley, A. C.; Steigerwald, M. L.; Nuckolls, C. Amine-Linked Single-Molecule Circuits: Systematic Trends Across Molecular Families. J. Phys.: Condens. Matter 2008, 20, 374115−374128. 20208


Article


The Journal of Physical Chemistry C (31) Hockney, R. W.; Goel, S. P.; Eastwood, J. Quiet HighResolution Computer Models of a Plasma. J. Comput. Phys. 1974, 14, 148−158. (32) Gollub, C.; Avdoshenko, S.; Gutierrez, R.; Berlin, Y.; Cuniberti, G. Can Propagating Charges Affect the Key Physical Quantities Controlling Their Motion? Isr. J. Chem. 2012, 52, 452−460. (33) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 7260−7268. (34) Pecchia, A.; Penazzi, G.; Salvucci, L.; Di Carlo, A. NonEquilibrium Green’s Functions in Density Functional Tight Binding: Method and Applications. New J. Phys. 2008, 10, 065022−065038. (35) Datta, S. Electronic Transport in Mesoscopic Systems; Cambridge University Press: New York, 1997. (36) Cuniberti, G., Fagas, G., Richter, K. Introducing Molecular Electronics; Springer: Berlin, 2005. (37) Cuevas, J. C.; Scheer, E. Molecular Electronics: An Introduction to Theory and Experiment; World Scientific: Singapore, 2010. (38) Sancho, M. L.; Sancho, J. L.; Sancho, J. L.; Rubio, J. Highly Convergent Schemes for the Calculation of Bulk and Surface Green Functions. J. Phys. F: Met. Phys. 1985, 15, 851−858. (39) Silin, V. I.; Wieder, H.; Woodward, J. T.; Valincius, G.; Offenhausser, A.; Plant, A. L. The Role of Surface Free Energy on the Formation of Hybrid Bilayer Membranes. J. Am. Chem. Soc. 2002, 124, 14676−14683. (40) Choo, H.; Cutler, E.; Shon, Y. S. Synthesis of Mixed MonolayerProtected Gold Clusters from Thiol Mixtures: Variation in the Tail Group, Chain Length, and Solvent. Langmuir 2003, 19, 8555−8559.

20209


Modeling of Solvent Effects in the Electrical ... - ACS Publications

Recommend Documents