TotalEnergy

class TotalEnergy(configuration)

Class for calculating the total energy

Parameters:configuration (MoleculeConfiguration | BulkConfiguration | DeviceConfiguration | SurfaceConfiguration) – Configuration with a calculator that supports total energy calculations.
alternativeEnergies()
Returns:A dictionary containing the values of other types of energies if available. If none are available it returns an empty dictionary.
Return type:dict
components()
Returns:A dictionary containing the energy components of the total free energy.
Return type:dict
evaluate()
Returns:The total free energy of the system.
Return type:PhysicalQuantity of type energy
metatext()
Returns:The metatext of the object or None if no metatext is present.
Return type:str | unicode | None
nlprint(stream=None)

Print a string containing an ASCII table useful for plotting the AnalysisSpin object.

Parameters:stream (python stream) – The stream the table should be written to.
Default: NLPrintLogger()
setMetatext(metatext)

Set a given metatext string on the object.

Parameters:metatext (str | unicode | None) – The metatext string that should be set. A value of “None” can be given to remove the current metatext.

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 and individual components
total_energy = TotalEnergy(molecule_configuration)
components = total_energy.components()

#print total energy and individual components
print() 
print('------------------------------------------------')
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

Calculate 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()
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

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

\[F[n] = T[n] + E^\mathrm{xc}[n]+E^{H}[n]+E^\mathrm{ext}[n]- \sigma 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^\mathrm{xc}[n]\) is the exchange-correlation energy and obtained with the method E.components()['Exchange-Correlation'].
  • \(E^{H}[n]\) includes all the electrostatic 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^\mathrm{ext}[n]\) is the interaction energy with an external field. The value of this term is obtained with the method E.components()['External-Field'].
  • \(- \sigma S\) is the entropy contribution due to smearing of the occupation function. Here \(\sigma\) is the broadening of the occupation function and \(S\) is the generalized entropy, see Occupation Methods. The value of this term is obtained with the method E.components()['Entropy-Term'].

It is possible to approximately correct for the effects introduced by using a smeared occupation function with a finite broadening, \(\sigma\), by the formula

\[E_{\sigma \to 0}(\sigma) = F_{\sigma \to 0}(\sigma) = E(\sigma) + \Delta E_{\sigma \to 0},\]

where the correction term, \(\Delta E_{\sigma \to 0}\) is determined by the specific occupation method used. This extrapolated zero-broadening energy is also written to output, but can also be obtained with the expression E.alternativeEnergies()['Zero-Broadening-Energy'].

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

Notes DeviceLCAOCalculator

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

\[F[n] = T[n] + E^\mathrm{xc}[n]+E^{H}[n]+E^\mathrm{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, \(- \sigma S\). The entropy term is not included since it usually is very insignificant and 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 [aST04]. 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 a DeviceConfiguration at finite bias. Therefore, 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^\mathrm{xc}[n]\) is the exchange-correlation energy and obtained with the method E.components()['Exchange-Correlation'].
  • \(E^{H}[n]\) includes all the electrostatic 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^\mathrm{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 TotalEnergy, use the method nlprint.

Notes HuckelCalculator

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

\[F[n] = E_\mathrm{1el}^0 +E^{H}[n]+E^\mathrm{ext}[n]- \sigma 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 [aKF08].

The terms in the equation are:

  • \(E_\mathrm{1el}^0 = E_\mathrm{1el} - \int V^H(\mathbf{r}) n(\mathbf{r}) d\mathbf{r}\) is the one-electron energy, subtracted electrostatic double counting terms. The term may also be written as \(E_\mathrm{1el}^0 = \mathrm{Tr}[H^0 D]\), where \(H^0\) is the non-self-consistent part of the Hückel 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^\mathrm{ext}[n]\) is the interaction energy with an external field. The value of this term is obtained with the method E.components()['External-Field'].
  • \(- \sigma S\) is the entropy contribution due to smearing of the occupation function. Here \(\sigma\) is the broadening of the occupation function and \(S\) is the generalized entropy, see Occupation Methods. The value of this term is obtained with the method E.components()['Entropy-Term'].

It is possible to approximately correct for the effects introduced by using a smeared occupation function with a finite broadening, \(\sigma\), by the formula

\[E_{\sigma \to 0}(\sigma) = F_{\sigma \to 0}(\sigma) = E(\sigma) + \Delta E_{\sigma \to 0},\]

where the correction term, \(\Delta E_{\sigma \to 0}\) is determined by the specific occupation method used. This extrapolated zero-broadening energy is also written to output, but can also be obtained with the expression E.alternativeEnergies()['Zero-Broadening-Energy'].

To print the data of the TotalEnergy, use the method nlprint.

Notes DeviceHuckelCalculator

For a DeviceConfiguration, the TotalEnergy object in a Hückel calculation will contain the terms

\[E[n] = E_\mathrm{1el}^0 +E^{H}[n]+E^\mathrm{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, \(- \sigma 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 a DeviceConfiguration.

The terms in the equation are:

  • \(E_\mathrm{1el}^0 = E_\mathrm{1el} - \int V^H(\mathbf{r}) n(\mathbf{r}) d\mathbf{r}\) is the one-electron energy, subtracted electrostatic double counting terms. The term may also be written as \(E_\mathrm{1el}^0 = \mathrm{Tr}[H^0 D]\), where \(H^0\) is the non-self-consistent part of the Hückel 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^\mathrm{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.

[aKF08]K. Kaasbjerg and K. Flensberg. Strong Polarization-Induced Reduction of Addition Energies in Single-Molecule Nanojunctions. Nano Letters, 8(11):3809–3814, 2008. doi:10.1021/nl8021708.
[aST04]A. P. Sutton and T. N. Todorov. A Maxwell relation for current-induced forces. Mol. Phys., 102(9-10):919–925, 2004. doi:10.1080/00268970410001703354.