To calculate an I-V curve, a self-consistent calculation for each voltage is required. The following script perform such a calculations for voltages, 0.2, 0.4, ..., 3.0 V.

#read in the 0 V configuration
device_configuration = nlread("au_dtb_au.nc",DeviceConfiguration)[0]
calculator = device_configuration.calculator()

# Define bias voltages
voltage_list= 0.2 *numpy.arange(1,16)*Volt
#make loop
for voltage in voltage_list:
    # Set new calculator with modified electrode voltages on the configuration
    # use the self consistent state of  the old calculation as starting input.
    device_configuration.setCalculator(
          calculator(electrode_voltages=(-0.5*voltage, 0.5*voltage)),
          initial_state=device_configuration)

    # Calculate the transmission spectrum
    transmission_spectrum = TransmissionSpectrum(
        configuration=device_configuration,
        energies=numpy.linspace(-3,3,100)*eV,
        kpoints=MonkhorstPackGrid(3,3),
        )

    #save the results
    nlsave('au_dtb_au.nc', device_configuration)
    nlsave('au_dtb_au.nc', transmission_spectrum)

ivscan.py

The script loops over a list of bias voltages and performs a self-consistent calculation for each bias. The self-consistent state of the previous calculation is used as initial guess for the next higher bias. The full calculation of the I-V curve will take about 20 times longer than the zero-bias calculation. If available, it is recommended to use a parallel computer. If you do not have access to a parallel computer, you may start the script on your local machine using the Job Manager and run it over-night.

In this chapter you will calculate the I-V curve from the transmission spectrum at each bias.

Launch the Custom Analyzer tool from the toolbar in the main VNL window.

Start the built-in I-V curve analyzer by selecting Analyzers → I-V Curve from the Custom Analyzer top menu bar.

Next, drag and drop the file au_dtb_au.nc onto the NC file drop zone.

Tick off symmetrize I-V and you should see the following plot

The left-most plots show the I-V and dI/dV-V curves. The curves are obtained by importing each transmission spectrum in the nc file and integrating the transmission coefficient over the bias window. In the input field to the left it is possible to change the electron temperature in the electrodes, and thereby the Fermi distribution used, when performing the integration over the bias windows.

The right-most plot shows the transmission spectra for different bias voltages, the spectra have been displaced vertically for clarity. The wedge signifies the part of the transmission spectra which are inside the bias window.

In the following you will analyze the transport properties at 1.6 V in a little more detail. To visualize the transmission spectrum, locate the file au_dtb_au.nc in the result browser, and open the transmission spectrum with ID gID018.

Note that the zero-bias peak at -1.5 eV appears to have split into two peaks at -0.8 and -2 eV. To further investigate the electronic structure for an applied bias of 1.6 V you will now calculate the LDDOS.

Using the same approach as in the section called “ Energy dependent LDDOS”, calculate LDDOS(E,z) of the system at 1.6 V. Below is shown the script for calculating the LDDOS.

configuration = nlread('au_dtb_au.nc', DeviceConfiguration)[8]

# make list of energies to be used for ldos
energies = numpy.linspace(-3,3,31)

# calculate ldos for each energy in the list
for e in energies:
    local_device_density_of_states = LocalDeviceDensityOfStates(
        configuration=configuration,
        energy=e*eV,
        kpoints=MonkhorstPackGrid(3,3),
        contributions=All,
        energy_zero_parameter=AverageFermiLevel,
        infinitesimal=1e-06*eV,
        self_energy_calculator=KrylovSelfEnergy(),
        spin=Spin.Sum,
        )
    nlsave(r'lddos16.nc', local_device_density_of_states)

lddos16.py

Compared with Figure 7 the resonance at -1.5 eV has indeed split into two resonances, at -0.8 and -2 eV, as we saw already in the transmission spectrum above. From the LDDOS plot we see that the resonance at -0.8 is located towards the left electrode, while the resonance at -2 eV is located towards the right electrode. The splitting arises from the drop in the electrostatic potential through the molecule, as will be investigated in the next section.

You will now calculate the voltage drop and induced density for an applied bias of 1.6 V. Carry out the analysis by dropping the following script named voltage_drop.py on the Job Manager.

# Read the configurations
configuration_list = nlread('au_dtb_au.nc', DeviceConfiguration)

#calculate the electrostatic potentials
potential_0V = ElectrostaticDifferencePotential(configuration_list[0])
potential_16V =ElectrostaticDifferencePotential(configuration_list[8])

# Calculate the voltage drop.
voltage_drop = potential_16V-potential_0V

# Save the voltage drop to voltage_drop.nc
nlsave('voltage_drop.nc', voltage_drop, object_id='drop16V')

#calculate the difference densities
density_0V = ElectronDifferenceDensity(configuration_list[0])
density_16V =ElectronDifferenceDensity(configuration_list[8])

# Calculate the difference
density_diff = density_16V-density_0V

# Save the density difference to voltage_drop.nc
nlsave('voltage_drop.nc', density_diff, object_id='diff16V')

voltage_drop.py

The result of this script is to save an ElectrostaticDifferencePotential object and an ElectronDifferenceDensity object into voltage_drop.nc. This potential corresponds to the potential drop across the device at 1.6 V, and the density is the induced density due to the applied bias. Visualizing the voltage drop as a contour plot and the induced density as an isosurface, the following plot can be obtained

The plot was obtained by setting the following properties for the widgets