Start up VNL and launch the Builder tool by clicking the icon on the Toolbar.

Click Add → From Database..., and when the database window opens up choose Databases → Molecules from the menu. In the search field, type benzene. Then click the plus icon in the lower right-hand corner to add the molecule to the Builder stash.

Click the Send To icon in the lower right-hand corner of the Builder and select Script Generator to transfer the benzene geometry to the Script Generator tool.

In the Script Generator:

The Script Generator window should now look like this.

Next double-click the MolecularEnergySpectrum icon that has appeared in the right panel in order to change the default parameters. Set the energy zero to Absolute Energy such that energies are reported relative to the Vacuum level.

Select File → Save and save the simulation script as benzene0.py.

Now drop the file benzene0.py on the Job Manager and start processing the queue.

After a few minutes the calculation has converged, and you can inspect the log file to determine the HOMO and the LUMO levels of the molecule. In the print-out of the molecular energy levels you should find

   14  -5.907277e+00   2.000000e+00
   15  -7.292105e-01   6.413524e-44

The occupations show that level 14 is the Highest Occupied Molecular Orbital (HOMO), while level 15 is the Lowest Unoccupied Molecular Orbital (LUMO). Thus,

It is instructive to compare these values with the ionization and affinity energy of benzene. The ionization energy, , is the cost of removing a single electron from the molecule and should be equal to . The affinity energy, is the energy gain of adding an additional electron, and should be equal to .

The experimental value for the ionization energy is = 9.25 eV [2], and the value is thus 3.35 eV off − a rather poor approximation.

The experimental affinity energy of benzene is negative, i.e. benzene cannot bind an additional electron and there is therefore no experimental data for the affinity energy. The experimental optical excitation energy of benzene is 10.54 eV [3] which gives an upper bound for the affinity energy of . The value is thus at least 2 eV off, again a very poor result.

All this is however to be expected: the HOMO and LUMO energies give a poor approximation for the ionization and affinity energies of benzene, since when the molecule becomes charged, there will be an additional charging energy which cannot be described by a calculation of the neutral molecule. In order to properly describe the charging energy, you need to perform self-consistent calculations for a charged molecule, and this is the topic of the next section.

Re-open the script generator tool, and change the default output file to benzene1.nc.

Next Double-left-click the New Calculator block and set the charge to 1

Now execute the script using the Job manager.

In the VNL main window, select benzene1.nc and left-click the “Total energy” and you should see the following

Similarly, inspect benzene0.nc and find the total energies

.

From these values, we can calculate the ionization energy

.

which is in excellent agreement with the experimental value.

Based on the calculated addition energies of benzene in the gas phase, you will now model how the molecule will function in a SET setup, where it is connected with gold electrodes. To use the formalism in the section called “The energy balance in the weak coupling regime”, you will also need the work function of gold [4],

.

Furthermore, you will need the dependence of the addition energies on the gate voltage, . To first order there is a linear dependence and the linear coefficient is called the gate coupling constant .

At this point the value of is unknown and it is in the following analysis set to 1 for simplicity. In the next chapter you will determine the gate coupling constant from a more detailed analysis of the fully self-consistent charge stability diagram. As it turns out, for this system is indeed close to 1, and the linear approximation is quite accurate.

The condition for transmission through a given molecular charge states is then given by

For a given value of and , you can now calculate the number of charge states in the bias window corresponding to the charge stability diagram. Below is a Python script, which uses the addition energies of the isolated benzene molecule (stored in the variable called delta_e, defined on line 7) to calculate the charge stability diagram.

import numpy
import pylab

#------------ Parameter definitons ---------------#

# Addition energies, already computed
delta_e = numpy.array([-15.73,-9.15,2.34,8.39])

# Workfunction for the assumed electrodes
w = 5.28

# Gate coupling constant
gate_coupling = 1.0

# Gate bias interval
v_g_interval = numpy.linspace(-15,15,151)

# Source-drain bias interval
v_sd_interval = numpy.linspace(-30,30,301)

#---------- End of parameter definitons ----------#
# test                                           test                         test                test               test                 test
# Calculate the number of charge states in the bias window
def conductionChannels(v_g,v_sd):
    return numpy.sum( abs( delta_e+w+gate_coupling*v_g) <= abs(v_sd/2) )

# Generate the mesh points of the contour plot
X, Y = numpy.meshgrid(v_g_interval,v_sd_interval)

# Evaluate the number of charge states for each mesh point
Z = [ conductionChannels(X[i,j],Y[i,j])
     for i in range(numpy.shape(X)[0])
          for j in range(numpy.shape(X)[1])]

Z = numpy.array(Z).reshape(numpy.shape(X))

# Make the plot
pylab.contour(X,Y,Z)
pylab.contourf(X,Y,Z)
pylab.xlabel("Gate voltage (Volt)")
pylab.ylabel("Source-Drain bias (Volt)")
pylab.show()

charge_stability_diagram.py

The so-obtained diagram is shown in Figure 5 (note the plot has been zoomed). As will be shown in the next chapter, the properties of the benzene molecule change, however, when it is placed in the SET environment, and the charge stability will be modified accordingly.