for voltage in numpy.arange(4.5, 5.01, 0.1):
two_probe_method = TwoProbeMethod(
electrode_parameters = (left_electrode_parameters,right_electrode_parameters),
exchange_correlation_type = exchange_correlation_type,
iteration_mixing_parameters = iteration_mixing_parameters,
electron_density_parameters = electron_density_parameters,
basis_set_parameters = basis_set_parameters,
iteration_control_parameters = iteration_control_parameters,
energy_contour_integral_parameters = energy_contour_integral_parameters,
two_center_integral_parameters = two_center_integral_parameters,
electrode_voltages = (voltage, 0)*Volt,
algorithm_parameters = two_probe_algorithm_parameters
)
if processIsMaster(): nlPrint(voltage)
import datetime
if processIsMaster(): print "Calculation started:", datetime.datetime.now().replace(microsecond=0)
runtime_parameters = runtimeParameters(
verbosity_level = 1,
checkpoint_filename = '/home/master/sweep/%.1fV.nc' % voltage
)
# Perform self-consistent field calculation
scf = executeSelfConsistentCalculation(
twoprobe_configuration,
two_probe_method,
runtime_parameters = runtime_parameters,
initial_calculation=scf
)
This works...
Cheers,
Serhan
This question highlights an important conceptual point about the way ATK calculated two-probe systems under bias, which I would like to take a second to comment on in more detail.
The simplest, and perhaps most common way to compute the I-V curve for a tunneling device is precisely like was attempted above: to first compute the zero-bias transmission, and then just integrate the transmission spectrum over an increasingly broad "bias window" to get the current. This works fine, in some sense, and is called the linear response current. However, it only gives accurate results under one crucial assumption: that the transmission spectrum is not influenced by the bias. In many cases this assumption is valid, in others it is just assumed, without proof (often because there is no simple way to check it).
There are however many situations where the transmission spectrum T(E) itself depends strongly on the applied bias. The uniqueness of ATK lies in its ability to compute T(E) self-consistently under an applied finite bias.
To do this you must, exactly like serhan shows, include the bias in the actual self-consistent calculation. In this way, the bias will be included in the self-consistent electrostatic potential, and consequently in the electron density and the density matrix, which in turn are used to compute the transmission.
Newer versions (later than 2008.02) of ATK actually also include a method to compute the linear response current, which works almost as you tried to do, artingchen. To use this function, you must however first import a module called ATK.Future, and instead of calling calculateCurrent() you use a different function. A rewritten version of your original attempt that utilizes this function would be as below.
I have made this is as a standalone analysis task, loading an already converged calculation from the checkpoint file.
from ATK.TwoProbe import *
import numpy
from ATK.Future import calculateLinearResponseCurrent
from ATK.MPI import processIsMaster
scf = restoreSelfConsistentCalculation ('my_calculation.nc')
my_voltages = numpy.linspace(-4.0, 4.0, 100)*Volt
lr_currents = calculateLinearResponseCurrent(
self_consistent_calculation = scf,
voltages = my_voltages,
brillouin_zone_integration_parameters = brillouinZoneIntegrationParameters((1, 1)),
green_function_infinitesimal = 1.0e-5*electronVolt,
number_of_points = 100
)
if processIsMaster():
file = VNLFile ('linear_response_current.vnl')
print 'Bias (V)\tCurrent (A)'
for i in range(len(my_voltages)):
print my_voltages[i].inUnitsOf(Volt), '\t', lr_currents[i].inUnitsOf(Ampere)
file.addToSample(lr_currents[i], 'twoprobe_configuration', 'Current at %g V bias' % my_voltages[i].inUnitsOf(Volt))
Note that all voltages are given at once, instead of looping over them. I have also made use of the very nice NumPy function linspace() instead of arange() which is a tricky function.
Oh, and I reduced the number of points between -4 and +4 V from 400 to 100. Note that if you want to (and you do want to!) go by serhan's approach, the best approach is always to start by 0 V bias,.
For more details, see the manual page for the function calculateLinearResponseCurrent() (http://quantumwise.com/documents/manuals/ATK-2008.10/ref.calculatelinearresponsecurrent.html).
Do note, however, that calculating the current in this way may not give the correct result, if, as mentioned above, the transmission does depend on the bias. And the only way to find out if it does it to compute it... Still, the linear response current can be useful as a first quick test to get an overview. It is however always crucial to also study the details of the transmission spectrum, also at finite bias.
You have only indented the definition of the two-probe method, not the actual execution of it.
You need
for voltage in numpy.arange(0.0, 5.01, 0.1):
two_probe_method = TwoProbeMethod(
electrode_parameters = (left_electrode_parameters,right_electrode_parameters),
exchange_correlation_type = exchange_correlation_type,
iteration_mixing_parameters = iteration_mixing_parameters,
electron_density_parameters = electron_density_parameters,
basis_set_parameters = basis_set_parameters,
iteration_control_parameters = iteration_control_parameters,
energy_contour_integral_parameters = energy_contour_integral_parameters,
two_center_integral_parameters = two_center_integral_parameters,
electrode_voltages = (voltage, 0)*Volt,
algorithm_parameters = two_probe_algorithm_parameters
)
if processIsMaster(): nlPrint(voltage)
runtime_parameters = runtimeParameters(
verbosity_level = 10,
checkpoint_filename = 'C:/Users/617/Desktop/lih2li1.nc'
)
# Using initial density from self consistent calculation
scf = executeSelfConsistentCalculation(
twoprobe_configuration,
two_probe_method,
initial_calculation = scf,
runtime_parameters = runtime_parameters
)
######################################################################
# Calculate physical properties
######################################################################
current = calculateCurrent(
self_consistent_calculation = scf,
brillouin_zone_integration_parameters = brillouinZoneIntegrationParameters((1, 1)),
green_function_infinitesimal = 1.0e-5*electronVolt,
number_of_points = 100
)
if processIsMaster(): nlPrint(current)
if processIsMaster(): file.addToSample(current, 'twoprobe_configuration', 'Current')