Name

TotalEnergy — Class for calculating the total energy

Synopsis

Namespace: NanoLanguage
TotalEnergy(configuration)

Description

Constructor for the TotalEnergy object.

TotalEnergy Arguments

configuration

Configuration with a calculator that supports total energy calculations

Type: A DeviceConfiguration, BulkConfiguration, or MoleculeConfiguration

Default: None

TotalEnergy Methods

A TotalEnergy object provides the following methods:

  • components(): Return dict with the energy components of the total energy.

  • evaluate(): Return the total energy as a physical quantity.

  • nlprint(stream): Print a string with the value of the total energy and its components.

    stream

    The stream the total energy should be written to.

    Type: A stream that supports strings being written to using 'write'.

    Default: sys.stdout

Usage Examples

Calculate the total energy of a hydrogen Molecule

#setup H-H geometry
elements = [ Hydrogen, Hydrogen]
cartesian_coordinates = [[ 0.0,0.0,0.0],
                         [ 0.0,0.0,0.75]]*Angstrom
molecule_configuration = MoleculeConfiguration(
    elements=elements,
    cartesian_coordinates=cartesian_coordinates
    )

#define calculator
calculator = LCAOCalculator()
molecule_configuration.setCalculator(calculator)

#calculate total energy
total_energy = TotalEnergy(molecule_configuration)
components = total_energy.components()
for term in components.keys():
    print '%-30s = %12.4f eV' % (term,components[term].inUnitsOf(eV))
print '----------------------------------------'
energy = total_energy.evaluate().inUnitsOf(eV)
print 'Total energy                   = %12.4f eV' % energy

h2_energy.py

Script for Calculating the equilibrium bond length and dissociation energy of a hydrogen molecule

#define calculator
calculator = LCAOCalculator()

#initialize arrays
energies = numpy.array([])
old_configuration = None

#setup distances where h2 energy is evaluated
distances = numpy.linspace(0.7,0.8,6)

#calculate total energy as function of distance
for d in distances:
    #setup H-H geometry
    elements = [ Hydrogen, Hydrogen]
    cartesian_coordinates = [[ 0.0,0.0,0.0],
                       [ 0.0,0.0,d]]*Angstrom
    molecule_configuration = MoleculeConfiguration(
        elements=elements,
        cartesian_coordinates=cartesian_coordinates
        )
    # set the calculator
    if old_configuration == None:
        molecule_configuration.setCalculator(calculator())
    else:
        molecule_configuration.setCalculator(calculator(),
                initial_state=old_configuration)
    #store the scf state
    old_configuration = molecule_configuration
    #calculate total energy
    total_energy = TotalEnergy(molecule_configuration)
    #append the total energy to energies
    energies = numpy.append(energies,total_energy.evaluate().inUnitsOf(eV))

#get the energy of atomic hydrogen
elements = [ Hydrogen]
cartesian_coordinates = [[ 0.0,0.0,0.0]]*Ang
molecule_configuration = MoleculeConfiguration(
        elements=elements,
        cartesian_coordinates=cartesian_coordinates
        )
molecule_configuration.setCalculator(calculator())
h_energy = TotalEnergy(molecule_configuration)

#calculate energies relative to atomic hydrogen
energies = energies - 2.*h_energy.evaluate().inUnitsOf(eV)

#fit polynomium to data
(a,b,c) = numpy.polyfit(distances,energies,2)

#print out the results
print 'distance(Angstrom) Energy (eV) polynomial_fit (eV)'
for i in range(len(distances)):
    d=distances[i]
    energy_fit = a*d*d+b*d+c
    print ' %8.2f         %8.4f         %8.4f '% (d,energies[i],energy_fit)

#calculate the bonding distance (experimental 0.742* Angstrom)
d_min = -b/(2.*a)
print 'bonding distance     = %8.3f Angstrom'% d_min

# calculate vibrational frequency (experimental 4395 cm-1)
k = 2.*a*eV/(Ang*Ang)
effective_mass = Hydrogen.atomicMass()/math.sqrt(2)
omega = (k/effective_mass)**(0.5)
wave_number = (omega/speed_of_light/2./math.pi).inUnitsOf(Meter**-1)/100.
print 'vibrational frequency = %8.f cm-1 '% (wave_number)

#calculate H-H bond dissociation energy (experimental 436 kJ/mol)
e_min = (a*d_min*d_min+b*d_min+c)*eV
#add zero point energy
e_min = e_min + 0.5 * hbar * omega
#convert to kj/mol
e_min_kjmol = (e_min*avogadro_number).inUnitsOf(Joule/Units.Mol)/1000.
print 'dissociation energy   = %8.f kJ/mol'% e_min_kjmol

h2.py

[Note] Note
To obtain a more accurate dissociation energy of the hydrogen molecule, we must perform a spin polarized calculation of the reference energy of the hydrogen atom

Notes LCAOCalculator

In ATK-DFT, the TotalEnergy object returns the free energy given by

\displaystyle

    E[n] = T[n] + E^{xc}[n]+E^{H}[n]+E^{ext}[n]- k_b T S

where n is the electron density

The terms in this equation are

  • T[n] is the kinetic energy of the Kohn-Sham orbitals and obtained with the method E.components()['Kinetic']

  • E^{xc}[n] is the exchange-correlation energy and obtained with the method E.components()['Exchange-Correlation']

  • E^{H}[n] is all the electro-static terms, i.e. the Hartree energy and the interaction energy with the pseudo potential ions. The sum of these terms is obtained with the method E.components()['Electrostatic']

  • E^{ext}[n] is the interaction energy with an external field. The value of this term is obtained with the method E.components()['External-Field']

  • - k_b T S is the Mermin free energy term, where k_b is the Boltzmann constant,  T is the electron temperature, and S = \sum_i [ f_i ln f_i + (1-f_i) ln(1-f_i) ] is the electronic entropy, given by the eigenstate occupation factors,f_i. The value of this term is obtained with the method E.components()['Entropy-Term']

To print the data of the total energy, use the method nlprint.

Notes DeviceLCAOCalculator

For a device configuration, the TotalEnergy object will in a DFT calculation contain the free energy functional given by

\displaystyle

    E[n] = T[n] + E^{xc}[n]+E^{H}[n]+E^{ext}[n]- e \, \Delta N_L \, \mu_L - e \, \Delta N_R \, \mu_R

where n is the electron density. Compared to the bulk system the electronic entropy term is missing,  - T S . The entropy term is not included since it is usually very insignificant, and it is very difficult and time consuming to calculate for a DeviceConfiguration.

The last two terms in the above equation arise from charge flow between the central region and the electrode reservoirs[22]. The symbol  \mu_L ( \mu_R ) is the electro-chemical potential of the left (right) electrode and  \Delta N is the excess number of electrons from the electrode. If the central region is charge neutral, we will have   \Delta N_L + \Delta N_R = 0 . As detailed below, for finite bias systems, these reservoir terms are calculated approximately, and the use of the total energy object is therefore approximate for DeviceConfigurations at finite bias, and the calculated force may in this case not be fully consistent with the derivative of the total energy.

The terms in the equation are

  • T[n] is the kinetic energy of the Kohn-Sham orbitals and obtained with the method E.components()['Kinetic']

  • E^{xc}[n] is the exchange-correlation energy and obtained with the method E.components()['Exchange-Correlation']

  • E^{H}[n] is all the electro-static terms, i.e. the Hartree energy and the interaction energy with the pseudo potential ions. The sum of these terms is obtained with the method E.components()['Electrostatic']

  • E^{ext}[n] is the interaction energy with an external field. This term is obtained with the method E.components()['External-Field']

  • -  e \, \Delta N_L \, \mu_L is the energy of the electrons which have entered the central region from the left reservoir. The number of electrons \, \Delta N_L is calculated approximately as the Mulliken charge of the atoms closest to the left electrode. This term is obtained with the method E.components()['Reservoir-Left']

  • -  e \, \Delta N_R \, \mu_R is the energy of the electrons which have entered the central region from the right reservoir. The number of electrons \, \Delta N_R is calculated approximately as the Mulliken charge of the atoms closest to the right electrode. This term is obtained with the method E.components()['Reservoir-Right']

To print the data of the total energy, use the method nlprint.

Notes HuckelCalculator

In ATK-SE, the TotalEnergy object returns a free energy containing only attractive total energy terms

\displaystyle

    E[n] =  E_{1el}^0 +E^{H}[n]+E^{ext}[n]- k_b T S

where n is the electron density. The total energy is missing a repulsive pair potential term. The energy can therefore only be used to compare the total energy of systems with the same distances between atoms. This, is for instance the case for the calculation of charging energies, i.e. the difference in the total energy as function of the charge on a molecule[24].

The terms in the equation are

  • E_{1el}^0 = E_{1el} - \int V^H({\bf r}) n({\bf r}) d{\bf r} is the one-electron energy, subtracted electro-static double counting terms. The term may also be written as  E_{1el}^0 = Tr[H^0 D] , where  H^0 is the non-self-consistent part of the Huckel Hamiltonian and D is the density matrix. The value of the term is obtained with the method E.components()['One-Electron']

  • E^{H}[n] is the Hartree energy and obtained with the method E.components()['Hartree']

  • E^{ext}[n] is the interaction energy with an external field. The value of this term is obtained with the method E.components()['External-Field']

  • - k_b T S is the Mermin free energy term[21], and the value of this term is obtained with the method E.components()['Entropy-Term']

To print the data of the total energy, use the method nlprint.

Notes DeviceHuckelCalculator

For a device configuration, the TotalEnergy object will in a Huckel calculation contain the terms

\displaystyle

       E[n] =  E_{1el}^0 +E^{H}[n]+E^{ext}[n]- e \, \Delta N_L \, \mu_L - e \, \Delta N_R \, \mu_R

where n is the electron density. Compared to the bulk system the DeviceConfiguration total energy object is missing the electronic free energy term,  - k_b T S . This term is omitted for convenience, since it is very difficult and time consuming to calculate for a DeviceConfiguration. For most systems this term is insignificant.

There are two new terms arising from charge flow between the central region and the electrode reservoirs. The reservoir terms are calculated approximately, and the use of the total energy object is therefore approximate for DeviceConfigurations.

The terms in the equation are

  • E_{1el}^0 = E_{1el} - \int V^H({\bf r}) n({\bf r}) d{\bf r} is the one-electron energy, subtracted electro-static double counting terms. The term may also be written as  E_{1el}^0 = Tr[H^0 D] , where  H^0 is the non-self-consistent part of the Huckel Hamiltonian and D is the density matrix. The value of the term is obtained with the method E.components()['One-Electron']

  • E^{H}[n] is the Hartree energy and obtained with the method E.components()['Hartree']

  • E^{ext}[n] is the interaction energy with an external field. The value of this term is obtained with the method E.components()['External-Field']

  • -  e \, \Delta N_L \, \mu_L is the energy of the electrons which have entered the central region from the left reservoir. The number of electrons \, \Delta N_L is calculated approximately as the Mulliken charge of the atoms closest to the left electrode. This term is obtained with the method E.components()['Reservoir-Left']

  • -  e \, \Delta N_R \, \mu_R is the energy of the electrons which have entered the central region from the right reservoir. The number of electrons \, \Delta N_R is calculated approximately as the Mulliken charge of the atoms closest to the right electrode. This term is obtained with the method E.components()['Reservoir-Right']

To print the data of the total energy, use the method nlprint.