the NLTypeError and how to find the bandgap??

[rosen]

I am just start to learn to use ATK, when I try to do the mulliken_population and the following Transmission calculation in a two-probe method, sometimes it will fail and give the following feedback:

Traceback (most recent call last):
  File "c:/docume~1/admini~1/locals~1/temp/tmpikz0uw.nl", line 271, in ?
    mulliken_population = calculateMullikenPopulation(self_consistent_calculation = scf)
NLTypeError: The population must contain non-negative numbers only.
Terminated Abnormally


Besides that, if I want to see the relationship of the insulator band gap  and the voltage in the two probe system, how to find the band gap (the exact LUMO and HOMO) value? It seems that the DOS is fluctuate dramatically and is hard to derive that band gap value.

[Nordland]

The reason for the error is that the system has converged to some unphysical state.
This is either a geometry issue or insufficient parameters. Perhaps if you posted your
script, I could help in trying to fix this issue.

Regarding the band gap - is it a molecule or is it bulk you are trying to determine the band gap in?

[rosen]

the whole script  :

from ATK.TwoProbe import *
from ATK.MPI import processIsMaster

# Generate time stamp
if processIsMaster():
    import platform, time
    print '#',time.ctime()
    print '#',platform.node(),platform.platform()+'\n'

# Opening vnlfile
if processIsMaster(): file = VNLFile('Oxy1v.vnl')

# Scattering elements and coordinates
scattering_elements = [Platinum, Platinum, Platinum, Platinum,
                       Platinum, Platinum, Platinum, Platinum,
                       Oxygen,   Oxygen,   Silicon,  Oxygen,   
                       Oxygen,   Silicon,  Oxygen,   Oxygen,   
                       Oxygen,   Silicon,  Oxygen,   Oxygen,   
                       Platinum, Platinum, Platinum, Platinum,
                       Platinum, Platinum, Platinum, Platinum]
scattering_coordinates = [[  0.00000000e+00,   0.00000000e+00,   6.79656744e+00],
                          [ -1.38734353e+00,   2.40294933e+00,   6.79656744e+00],
                          [  2.77468705e+00,   0.00000000e+00,   6.79656744e+00],
                          [  1.38734353e+00,   2.40294933e+00,   6.79656744e+00],
                          [ -1.38734353e+00,   4.00491571e+00,   9.06208992e+00],
                          [  0.00000000e+00,   1.60196626e+00,   9.06208992e+00],
                          [  1.38734353e+00,   4.00491571e+00,   9.06208992e+00],
                          [  2.77468705e+00,   1.60196626e+00,   9.06208992e+00],
                          [  1.20103330e+00,   2.32279668e+00,   1.81736437e+01],
                          [  1.58815643e+00,   4.62207968e+00,   1.64930061e+01],
                          [  1.80105513e+00,   2.87943831e+00,   1.66156304e+01],
                          [  3.51000495e+00,   2.48519392e+00,   1.64697903e+01],
                          [ -1.02802235e+00,   2.28316448e+00,   1.51466260e+01],
                          [  6.81155340e-01,   1.88944240e+00,   1.50025398e+01],
                          [  8.94502430e-01,   1.46872210e-01,   1.51234308e+01],
                          [  1.28141397e+00,   2.44658245e+00,   1.34461835e+01],
                          [ -3.96967870e-01,   2.88425220e+00,   1.11567331e+01],
                          [  1.20751363e+00,   2.39587521e+00,   1.16897585e+01],
                          [  1.53630998e+00,   7.57744760e-01,   1.11369934e+01],
                          [  2.40856772e+00,   3.49402221e+00,   1.10200002e+01],
                          [ -1.38734353e+00,   4.00491571e+00,   2.01299996e+01],
                          [  0.00000000e+00,   1.60196626e+00,   2.01299996e+01],
                          [  1.38734353e+00,   4.00491571e+00,   2.01299996e+01],
                          [  2.77468705e+00,   1.60196626e+00,   2.01299996e+01],
                          [  8.88178420e-16,   3.20393252e+00,   2.23955221e+01],
                          [  1.38734353e+00,   8.00983131e-01,   2.23955221e+01],
                          [  2.77468705e+00,   3.20393252e+00,   2.23955221e+01],
                          [  4.16203070e+00,   8.00983131e-01,   2.23955221e+01]]*Angstrom
       

electrode_elements = [Platinum, Platinum, Platinum]
electrode_coordinates = [[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
                         [  4.44089210e-16,   1.60196629e+00,   2.26552246e+00],
                         [  1.38734350e+00,   8.00983146e-01,   4.53104491e+00]]*Angstrom

electrode_cell = [[ 2.77468701,  0.        ,  0.        ],
                  [-1.3873435 ,  2.40294944,  0.        ],
                  [ 0.        ,  0.        ,  6.79656737]]*Angstrom

# Set up electrodes
electrode_configuration = PeriodicAtomConfiguration(
    electrode_cell,
    electrode_elements,
    electrode_coordinates
    )

# Set up two-probe configuration
twoprobe_configuration = TwoProbeConfiguration(
    (electrode_configuration,electrode_configuration),
    scattering_elements,
    scattering_coordinates,
    electrode_repetitions=[[2,2],[2,2]],
    equivalent_atoms=([0,0],[2,27])
    )
if processIsMaster(): nlPrint(twoprobe_configuration)
if processIsMaster(): file.addToSample(twoprobe_configuration, 'Oxy1v')

######################################################################
# Central region parameters
######################################################################
exchange_correlation_type = GGA.PBE

iteration_mixing_parameters = iterationMixingParameters(
    algorithm = IterationMixing.Pulay,
    diagonal_mixing_parameter = 0.1,
    quantity = IterationMixing.Hamiltonian,
    history_steps = 6
)

electron_density_parameters = electronDensityParameters(
    mesh_cutoff = 100.0*Rydberg
)

basis_set_parameters = basisSetParameters(
    type = DoubleZetaPolarized,
    radial_sampling_dr = 0.001*Bohr,
    energy_shift = 0.01*Rydberg,
    delta_rinn = 0.8,
    v0 = 40.0*Rydberg,
    charge = 0.0,
    split_norm = 0.15
)

iteration_control_parameters = iterationControlParameters(
    tolerance = 1e-005,
    criterion = IterationControl.TotalEnergy,
    max_steps = 100
)

electrode_voltages = (0.0,1.0)*Volt

two_probe_algorithm_parameters = twoProbeAlgorithmParameters(
    electrode_constraint = ElectrodeConstraints.RealSpaceDensity,
    initial_density_type = InitialDensityType.EquivalentBulk
)

energy_contour_integral_parameters = energyContourIntegralParameters(
    circle_points = 30,
    integral_lower_bound = 3*Rydberg,
    fermi_line_points = 10,
    fermi_function_poles = 4,
    real_axis_infinitesimal = 0.01*electronVolt,
    real_axis_point_density = 0.02*electronVolt
)

two_center_integral_parameters = twoCenterIntegralParameters(
    cutoff = 2500.0*Rydberg,
    points = 1024
)

######################################################################
# Left electrode parameters
######################################################################
left_electrode_electron_density_parameters = electronDensityParameters(
    mesh_cutoff = 100.0*Rydberg
)

left_electrode_iteration_control_parameters = iterationControlParameters(
    tolerance = 1e-005,
    criterion = IterationControl.TotalEnergy,
    max_steps = 100
)

left_electrode_brillouin_zone_integration_parameters = brillouinZoneIntegrationParameters(
    monkhorst_pack_parameters = (2, 2, 100)
)

left_electrode_iteration_mixing_parameters = iterationMixingParameters(
    algorithm = IterationMixing.Pulay,
    diagonal_mixing_parameter = 0.1,
    quantity = IterationMixing.Hamiltonian,
    history_steps = 6
)

left_electrode_eigenstate_occupation_parameters = eigenstateOccupationParameters(
    temperature = 300.0*Kelvin
)

######################################################################
# Collect left electrode parameters
######################################################################
left_electrode_parameters = ElectrodeParameters(
    brillouin_zone_integration_parameters = left_electrode_brillouin_zone_integration_parameters,
    electron_density_parameters = left_electrode_electron_density_parameters,
    eigenstate_occupation_parameters = left_electrode_eigenstate_occupation_parameters,
    iteration_mixing_parameters = left_electrode_iteration_mixing_parameters,
    iteration_control_parameters = left_electrode_iteration_control_parameters
)

######################################################################
# Right electrode parameters
######################################################################
right_electrode_electron_density_parameters = electronDensityParameters(
    mesh_cutoff = 100.0*Rydberg
)

right_electrode_iteration_control_parameters = iterationControlParameters(
    tolerance = 1e-005,
    criterion = IterationControl.TotalEnergy,
    max_steps = 100
)

right_electrode_brillouin_zone_integration_parameters = brillouinZoneIntegrationParameters(
    monkhorst_pack_parameters = (2, 2, 100)
)

right_electrode_iteration_mixing_parameters = iterationMixingParameters(
    algorithm = IterationMixing.Pulay,
    diagonal_mixing_parameter = 0.1,
    quantity = IterationMixing.Hamiltonian,
    history_steps = 6
)

right_electrode_eigenstate_occupation_parameters = eigenstateOccupationParameters(
    temperature = 300.0*Kelvin
)

######################################################################
# Collect right electrode parameters
######################################################################
right_electrode_parameters = ElectrodeParameters(
    brillouin_zone_integration_parameters = right_electrode_brillouin_zone_integration_parameters,
    electron_density_parameters = right_electrode_electron_density_parameters,
    eigenstate_occupation_parameters = right_electrode_eigenstate_occupation_parameters,
    iteration_mixing_parameters = right_electrode_iteration_mixing_parameters,
    iteration_control_parameters = right_electrode_iteration_control_parameters
)

######################################################################
# Initialize self-consistent field calculation
######################################################################
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 = electrode_voltages,
    algorithm_parameters = two_probe_algorithm_parameters
)
if processIsMaster(): nlPrint(two_probe_method)


runtime_parameters = runtimeParameters(
    verbosity_level = 10,
    checkpoint_filename = 'Oxy1v.nc'
)

# Using initial density from self consistent calculation
scf = executeSelfConsistentCalculation(
    twoprobe_configuration,
    two_probe_method,
    runtime_parameters = runtime_parameters
)
######################################################################
# Calculate physical properties
######################################################################
atomic_forces = calculateAtomicForces(self_consistent_calculation = scf)
if processIsMaster(): nlPrint(atomic_forces)
if processIsMaster(): file.addToSample(atomic_forces, 'Oxy1v', 'Atomic Forces')

current = calculateCurrent(
    self_consistent_calculation = scf,
    brillouin_zone_integration_parameters = brillouinZoneIntegrationParameters((2, 2)),
    green_function_infinitesimal = 1.0e-5*electronVolt,
    number_of_points = 100
)
if processIsMaster(): nlPrint(current)
if processIsMaster(): file.addToSample(current, 'Oxy1v', 'Current')

import numpy
density_of_states = calculateDensityOfStates(
    self_consistent_calculation = scf,
    energies = numpy.arange(-10.0, 10.0+0.004, 0.04)*electronVolt,
    brillouin_zone_integration_parameters = brillouinZoneIntegrationParameters((2, 2)),
    green_function_infinitesimal = 1.0e-5*electronVolt
)
if processIsMaster(): nlPrint(density_of_states)
if processIsMaster(): file.addToSample(density_of_states, 'Oxy1v', 'Density of states')

mulliken_population = calculateMullikenPopulation(self_consistent_calculation = scf)
if processIsMaster(): nlPrint(mulliken_population)
if processIsMaster(): file.addToSample(mulliken_population, 'Oxy1v', 'Mulliken Population')

# Set verbosity level so that all energy components are printed
import ATK
verbosity_level=ATK.verbosityLevel()
ATK.setVerbosityLevel(10)
total_energy = calculateTotalEnergy(self_consistent_calculation = scf)
ATK.setVerbosityLevel(verbosity_level)

if processIsMaster(): nlPrint(total_energy,'Total energy')
if processIsMaster(): file.addToSample(total_energy, 'Oxy1v', 'Total energy')

import numpy
transmission_spectrum = calculateTransmissionSpectrum(
    self_consistent_calculation = scf,
    energies = numpy.arange(-10.0, 10.0+0.004, 0.04)*electronVolt,
    brillouin_zone_integration_parameters = brillouinZoneIntegrationParameters((2, 2)),
    green_function_infinitesimal = 1.0e-5*electronVolt
)
if processIsMaster(): nlPrint(transmission_spectrum)
if processIsMaster(): file.addToSample(transmission_spectrum, 'Oxy1v', 'Transmission Spectrum')




In the output files, after the Density of States (some of its value is negative), there is the error message:

Traceback (most recent call last):
  File "Oxy1v.py", line 269, in ?
    mulliken_population = calculateMullikenPopulation(self_consistent_calculation = scf)
NLTypeError: The population must contain non-negative numbers only.


Before the calculation, the structure didn't do the structure optimization, I think that is the mainly reason of these error. But with the bias of 0.1V, the calculation (also without optimize) can finish without any error report.


the band gap I want to calculate is about the molecule between the two proble system. For example, in the above script it is three SiO2 chain, one may expect that with the increasing voltage, the band gap of this insulator chain embeded in the electrode will decrease (from 5eV to 1eV or even zero), I just want to verify it and don't know how to get its bandgap value? because the DOS value fluctuate without a smooth gap.

[Nordland]

Hey Rosen.

I have posted a script for calculating the band gap in bulk systems.
You can find it here - ( Band gap for bulk systems )

For molecule it is straight forward, you simply perform a calculation on the desired molecule, and
calculate the molecular spectrum. The substract the first positive energy from the last negative and you get the homo-lumo:

molecular_energy_spectrum = calculateMolecularEnergySpectrum(self_consistent_calculation = scf)
nlPrint(molecular_energy_spectrum)

For the H2 molecule this gives the following output:

# -----------------------------------------------------------------------------
# Energy Spectrum
# -----------------------------------------------------------------------------
# Energy (eV)
       -10.17
         2.84
         8.41
        12.92
        12.92
        16.42
        31.75
        32.56
        32.56
        70.26

So in this case the band gap is 2.84 - (-10.17) = 13.01 eV.

[rosen]

Dear Nordland
Thank you for your script, is it possible to get the bandgap of the insulator embeded in a two-probe system?

[Nordland]

Hej Rosen.

You can either perform a bulk calculation of the twoprobe system as described here (Equivalent Bulk Calculation of TwoProbe) and the use my script for finding the bandgap.

However if it is a molecular junction you are modelling, then I suggest making a seperate calculation of the molecule alone, and use the molecular spectrum to determine the band gap.

[rosen]

OK, thank you very much, it will be a better way to determine its bandgap as a separate molecule, and it won't have the influence of the electrode.

[Anders Blom]

A tip: What's really interesting is to compare the calculation made on the isolated molecule (the "intrinsic" homo-lumo gap) with that of the molecule embedded in the two-probe system. The molecular energy levels are typically shifted (and broadened) by the coupling to the electrodes, and often this can have a profound influence on the transmission properties. The hall-mark example of this is the famous Li-H2-Li system.

To evaluate the "HOMO" and "LUMO" levels (or what corresponds to it) of the molecule embedded in the electrode, use the function calculateProjectedHamiltonianEnergySpectrum() (see the online manual) and make sure to designate all central region atoms except those in the molecule as "surface atoms".

The broadening of the levels shows up if you calculate the (surface) density of states for the two-probe system, provided that the levels are coupled to the electrodes and not localized.

[Anders Blom]

Second, regarding the problem of converging the system to a physical state: I suggest using at least 6 layers of Pt atoms in the electrode instead of 3. Cf. the discussion here: http://quantumwise.com/forum/index.php?topic=7.msg22#msg22

Also, the system is quite small in the XY plane, and hence a larger k-point sampling than (2x2) is advisable.

Finally, unless you've already done so, of course, always calculate the zero-bias case first, before attempting the finite bias. Then, use the zero bias state to initialize the finite bias calculation, otherwise it's often very hard to converge it, especially at such a relatively high bias at 1 V.