Vibration analysis

You can calculate the vibrational modes of a ListOfAtoms in the harmonic approximation using the Using the

For example, assuming you have a ListOfAtoms called atoms with an attached calculator, you can use the module like this:

>>> from ASE.Utilities.VibrationalCalculation import *
>>> vib=VibrationalCalculation(filetoken='CO_vib',
...                            atoms=loa,
...                            freeatoms=[0,1],
...                            displacements=[0.05] * 2,
...                            method=0)
>>> vib.RunCalculations()
>>> vib.PrintFrequencies()
>>> print 'Zero-point energy = %1.2f eV' % vib.GetZeroPointEnergy()

filetoken is a string that is prefixed to the names of all the files created. atoms is a ListOfAtoms that is either at a fully relaxed ground state or at a saddle point. freeatoms is a list of the atoms which the vibrational modes will be calculated for, the rest of the atoms are considered frozen. displacements is a list of displacements, one for each free atom that are used in the finite difference method to calculate the Hessian matrix. method is -1 for backward differences, 0 for centered differences, and 1 for forward differences.


Using the dacapo caculator you must make sure that the symmetry program in dacapo finds the same number of symmetries for the displaced configurations in the vibrational modules as found in the ground state used as input. This is because the wavefunction is reused from one displacement to the next. One way to ensure this is to tell dacapo not to use symmetries.

This will show op as a python error 'Frames are not aligned'. This could be the case for other calculators as well.

You can get a NetCDF trajectory corresponding to a specific mode by using:

>>> mode=0
>>> vib.CreateModeTrajectory(mode=mode,scaling=5)

This will create a NetCDF trajectory file CO_vib_mode_0.traj, corresponding to the highest frequency mode. scaling is an option argument, that will give the amplitude of the mode, default is 10.

Maximally localized Wannier functions

Wannier functions are maximally localized orbitals constructed from a set of extended eigenstates. This page describes how to construct the Wannier orbitals using the class Wannier.

Band structure and orbital analysis, based on the set of localized Wannier orbitals are discussed in the chapter band structure analysis.

For the theory behind this method see the paper Partly Occupied Wannier Functions, Thygesen, Hansen and Jacobsen, Phys. Rev. Lett, Vol.94, 26405 (2005).

The class Wannier

A simple initialization of the class Wannier looks like this:

>>> from ASE.Utilities.Wannier import Wannier
>>> atoms = Calculator.ReadAtoms('')
>>> calc = atoms.GetCalculator()
>>> wannier = Wannier(numberofwannier=6,calculator=calc)

The first argument to the Wannier constructor must be the number of Wannier functions to be constructed, in this case 6. The next argument is a calculator containing all the information about the electronic structure calculation, i.e. number of bands, k-points, wavefunctions and so on.

For examples of how to use the Wannier class, see the Wannier tutorial.


For a calculation using k-points the relevant unit cell for eg. visualization of the Wannier orbitals is not the original unit cell, but rather a larger unit cell defined by repeating the original unit cell by the number of k-points in each direction. We will refer to this unit cell as the large unit cell. Note that for a Γ-point calculation the large unit cell coinsides with the original unit cell. The large unitcell defines also the periodicity of the Wannier orbitals.


For calculations using k-points, make sure that the Γ-point is included in the k-point grid. Moreover you must shift all k-points by a small amount (but not less than 2e-5 in) in e.g. the x direction, before performing the Dacapo calculation. If this is not done the symmetry program in Dacapo will use time-reversal symmetry to reduce the number of k-points by a factor 2. The shift can be performed like this:

kpoints = calc.GetBZKPoints()
kpoints[[,0]] += 2e-5

Below is the full list of keyword arguments:

numberofwannier: Number of Wannier orbitals.
The number of Wannier orbitals to be constructed must be <= numberofbands.
calculator: Calculator holding the eigenstates, the unit cell and
providing the method GetWannierLocalizationMatrix.
numberofbands: Total number of bands defining the external
localization space. This number is optional. It must be <= the total number of bands in the DFT calculation. Default is the total number of bands used in the DFT calculation.

spin: Select specific spin for a spin-polarized calculation (0 or 1).

occupationenergy: Eigenstates below this energy are included in
the internal localization space. In practice this means that all eigenstates below this energy can be exactly reproduced in terms of the resulting Wannier functions. The energy is relative to the Fermi-level, and 0 is default (all occupied states are included in the localization space).
numberoffixedstates: Fix the number of states to be included in
the internal localization space starting from the lowest eigenstate. This keyword provides another way of specifying how many states should be included, and overrides occupationenergy if this is also set. Default is None meaning that the number of fixed states is set by the occupationenergy keyword.
initialwannier: Setup an initial set of Wannier orbitals.

initialwannier can set up a starting guess for the Wannier functions. This is important to speed up convergence in particular for large systems For transition elements with d electrons you will always find 5 highly localized d-orbitals centered at the atom. Placing 5 d-like orbitals with a radius of 0.4 Angstroms and center at atom no. 7, and 3 p-like orbitals with a radius of 0.4 Angstroms and center at atom no. 27 looks like this:

initialwannier = [[[7],2,0.4],[[27],1,0.4]]

Placing only the l=2, m=-2 and m=-1 orbitals at atom no. 7 looks like this:

initialwannier = [[[7],2,-2,0.4],[[7],2,-1,0.4]]

I.e. if you do not specify the m quantum number all allowed values are used. Instead of placing an orbital at an atom, you can place it at a specified position. For example the following:

initialwannier = [[[0.5,0.5,0.5],0,0.5]]

places an s orbital with radius 0.5 Angstroms at the position (0.5,0.5,0.5) in scaled coordinates of the unit cell.


For the ethylene molecule (12 valence electrons) we assume that we have the Dacapo output file

>>> atoms = Calculator.ReadAtoms('')
>>> calc = atoms.GetCalculator()
>>> wannier = Wannier(numberofwannier=6,calculator=calc)

This will construct 6 Wannier orbitals using the 6 occupied valence states, corresponding to the 12 valence electrons.

Below is a list of the most important methods of the class Wannier:

Localize(step=0.5,tolerance=1.0e-08): Perform the localization of the
Wannier orbitals. This method will localize the Wannier orbitals, i.e will try to maximize the localization functional. If the functional value is not increasing decrease the step size from the default value 0.5.
GetWannierFunctionDOS(n,energies,width): Get projected density of states for WF.
Returns the projected density of states (PDOS) for Wannier function n. The calculation is performed over the energy grid specified in energies. The PDOS is produced as a sum of Gaussians centered at the points of the energy grid and with the specified width.
Same as GetWannierFunctionDOS, but writes the output to a NetCDF file.
GetElectronicState(wannierindex,repeat=None): Returns an ElectronicState instance
corresponding to the Wannier orbital with index wannierindex. The keyword repeat can be a list of 3 integers [n1,n2,n3], specifying how many times the unit cell is repeated along the unit cell basis vectors.
GetCentersAsAtoms: Returns a ListOfAtoms object with the Wannier centers.
The chemical element is set to 'X'.
TranslateAllWannierFunctionsToCell(cell): Move all Wannier orbitals to a specific unit cell.
There exists an arbitrariness in the positions of the Wannier orbitals relative to the unit cell. This method can move all orbitals to the unit cell specified by cell. For a Γ-point calculation, this has no effect. For a k-point calculation the periodicity of the orbitals are given by the large unit cell defined by repeating the original unitcell by the number of k-points in each direction. In this case it is usefull to move the orbitals away from the boundaries of the large cell before plotting them. For a bulk calculation with, say 10x10x10 k points, one could move the orbitals to the cell [2,2,2]. In this way the periodic boundary conditions will not be noticed.
WriteCube(wannierindex,filename,repeat=(7,7,7),real=False): Write a Cube formatted file.

A Cube formatted file is written for the given wannier index. repeat can be used to repeat the unitcell, this is only relevant for calculations using k-points. In this case repeat, will default be the number of k-points in each directions, i.e for a 11x11x11 k-point set, repeat will be (11x11x11). This cell size represents the periodicity of the Wannier orbitals.

Localized Wannier functions can often be chosen to be real. If the keyword real is set to True, the complex Wannier function will be transformed into a real one by multiplication be a suitable phase factor. In VMD you can use this to add two isosurfaces using +- isosurface value, to get an approximation for the sign of the Wannier function.

Save/ReadZIBlochMatrix(filename): Save and read ZI bloch matrix.
These methods save and restore the localization matrix generated from the initial set of bloch function. This can save time, since the ZI matrix must be provided each time a localization is performed. If a ZI matrix is not read from file, it will be calculated.
These methods can be used to save and restore the unitary matrices used to produce a set of Wannier functions. The method produces two files: filename_rot.pickle and filename_coeff.pickle.


You can save your localized Wannier orbital like this:

>>> wannier = Wannier(..)
>>> wannier.SaveZIBlochMatrix('fe_bloch.pickle')
>>> wannier.Localize(tolerance=0.000001)
>>> wannier.SaveRotation('fe')

and read them in again like this:

>>> wannier = Wannier(..)
>>> wannier.ReadZIBlochMatrix('fe_bloch.pickle')
>>> wannier.ReadRotation('fe')
>>> wannier.Localize(tolerance=0.000001)

Localize should now converge in one step.

Band structure and orbital analysis

The class HoppingParameters can generate a band structure using the set of Wannier orbitals.

An instance of HoppingParameters is initialized like this:

>>> from ASE.Utilities.Wannier import HoppingParameters
>>> hop = HoppingParameters(wannier,cutoff)

The cutoff distance truncates the Wannier orbitals at the specified distance. This distance should be smaller than half the length of large unitcell. The truncation is necessary because the Wannier functions will always be periodic (with a periodicity given by the large cell), and thus in order to describe completely localized orbitals the WFs must be truncated.

HoppingParameters have the following methods:

Returns the matrix element <n,0|H|m,R>, where (n,0) is Wannier function number n in unit cell 0=[0,0,0], and (m,R) and m is Wannier function number m in unit cell R=[n1,n2,n3] where n1,n2,n3 are integers.
WriteBandDiagramToNetCDFFile(filename,npoints,kpt1,kpt2): Write a band diagram to file.
A band structure plot is written to file filename. There will be npoints k-points distributed uniformly along the line connecting kpt1 and kpt2 in the BZ. Each coordinate of kpt1 and kpt2 should be between -0.5 and 0.5.
The Hamiltonian matrix in the basis of the Wannier orbitals are returned. We will refer to this Hamiltonian as H in that follows. The Hamiltonian refers to the large unit cell, and its dimension is therefore (N_w*N_k)x(N_w*N_k), where N_w is the number of Wannier functions in a unit cell and N_k is the number of k points. Periodic boundary conditions are imposed on the boundaries of the large cell.

The module HamiltonianTools have a number of useful methods for analysing problems in terms of the Wannier functions and the Hamiltonian matrix H. Definition and physical meaning of the term group-orbital (see below) can be found in the paper PRL 94,036807 (2005). The module is imported like this:

>>> from ASE.Utilities.Wannier import HamiltonianTools

The methods are described below:

H_rot,U,eigenvalues = HamiltonianTools.SubDiagonalize(h,listofindices):
This methods diagonalize the Hamiltonian h within the subspace spanned by the basis functions (Wannier functions) speficied in the list listofindices. This can be used to e.g. to obtain renormalized molecular orbitals for a molecule adsobed on a surface, by diagonalizing h within the subspace spanned by the Wannier functions centered at the molecule. H_rot will be the transformed Hamiltonian matrix, U is the unitary matrix that relates H_rot to the original h, and eigenvalues are the eigenvalues (energies) in the diagonalized subspace.
Returns the coupling matrix element between a basis function (Wannier function or renormalized orbital - see method above) and its so-called group orbital.
Returns the matrix h with all couplings involving the basis functions specified in the list indexlist set to zero.
Returns the projected density of states (PDOS) for the orbitals specified in listoforbitals. Each entity in listoforbitals can be an integer (the index of a basis function) or a normalized list of coordinates, depending on whether one wants the PDOS for a specific basis function or a linear- combination of such. hamiltonian is a Hamiltonian matrix, listofenergies is a Python array with an energy grid on which the PDOS is returned, and width sets the smearing scale of the PDOS.
The GetWannierLocalizationMatrix calculator method

The GetWannierLocalizationMatrix returns the matrix:


Here n and m runs over the number of bands and G_a a set of reciprocal lattice vectors of the unit cell, specified as one of:

GI = GI × R

G_I is one of:


and R is the reciprocal unit cell.

For the case of a cubic unitcell with length L, only 3 G_I values will be used:

G_x = 2pi*(1,0,0)/L
G_y = 2pi*(0,1,0)/L
G_z = 2pi*(0,0,1)/L
Γ-point calculations

The calculator should implement the method:


For the Γ-point calculations kpoint, nextkpoint , and dirG can be ignored.

Uniform k-point sampling

Under construction ...


The STMTool class generates Tersoff-Hamann STM topographs.

Here is how you could make a STM picture from the GPAW code, assuming you have allready generated the the Al100.gpw GPAW restart file, (see the STM tutorial for a full description):

>>> from ASE.Utilities.ElectronicStates import ElectronicStates
>>> from ASE.Utilities.STMTool import STMTool
>>> from ASE.Visualization.VTK import VTKPlotSTM
>>> loe = ElectronicStates(filename = 'Al100.gpw')
>>> stm = STMTool(loe,contourvalue=1.0e-4, channel="s",
...               normalaxis=2, smoothfactor=0.1)    # available channels: s,px,py,pz
>>> stmplot = VTKPlotSTM(stm)
>>> stmplot.SetColorRange((3.325,3.55))      # from surface span color range
>>> stmplot.Render()                         # updates picture

From the planewave code dacapo, the dacapo codes own version of the ElectronicStates class must be used:

>>> from Dacapo.ElectronicStates import ElectronicStates
>>> from ASE.Utilities.STMTool import STMTool
>>> from ASE.Visualization.VTK import VTKPlotSTM
>>> loe = ElectronicStates(filename='')

From this point the code is identical to the example above.

The list of keywords for the STMTool object is:

The ElectronicStates object holding information about each eigenstate of the slab.
The density of states at the Fermi level for which the topograph is generated. This argument is mandatory, as no universal (simple) default is natural. Often used values are [1.0e-6 - 1.0e-3], but is very system dependent.
Selects whether the topograph is generated by the value or gradient isosurface of the Fermi level DOS. This adjusts for the orbital character of the DOS of the actual STM tip in the probing energy window, i.e. the STM tip states that couples to the substrate. Available options are "s", "px", "py" and "pz". Default is "s".
The surface normal direction, (0,1,2) for (x,y,z) respectively. normalaxis = 2 is default.
The thermal convolution applied to the slab states. It is necessary to apply this due to finite k-point sampling.

STM linescans

The class STMLineScan can be used to simulate an STM linescan. It must be initialized with an instance of the STMTool class.

Note that the linescans are limited to be along one of the axis of the unit cell. The origin of the linescan can only assume discrete values corresponding to the FFT grid.

The class is used like:

linescan = STMLineScan(stm,fftnumber=12,axis=1)

See more details in the linescan section of the STM tutorial.

Bader Analysis

Henkelman et al. have implemented a fast and robust algorithm for calculating the electronic charges on individual atoms in molecules or crystals with the Bader scheme (G. Henkelman, A. Arnaldsson, and H. Jónsson, A fast and robust algorithm for Bader decomposition of charge density, Comput. Mater. Sci. (in press, 2005)). In that method the electron density is analyzed and a so-called zero-flux surfaces are used to divide a system into atoms. A more detailed description about the method can be found at Richard Bader's homepage at McMaster University and more details about the algorithm (how to use, about the output files, discussion forum and an article about the algorithm) can be found at Graeme Henkelman's homepage at the University of Texas at Austin. A presentation about the Bader method and this algorithm can be found here.

This algorithm suites very well for large solid state physical systems as well as large biomolecular systems. The computational time depends only on the size of the 3D grid used in the calculations and typically it takes less than a minute to do the analysis. To get more accurate results in the charge analysis it is recommended to use a finer grid than the default values, however, usually it is enough to calculate the electronic structure with a finer grid after the optimization has been done with the default values.

To be able to use this algorithm one has to create cube files from the NetCDFiles of e.g. a dacapo calculation. These cube files contain the size of the unitcell, coordinates of the atoms, and the electronic charge density on a grid. This can be done with a simple python script:

>>> from Dacapo import Dacapo
>>> from ASE.IO.Cube import *
>>> atoms = Dacapo.ReadAtoms('')
>>> calc = atoms.GetCalculator()
>>> dens = calc.GetDensityArray()
>>> density =  dens * (0.529177)**3     # changes from A^-3 units to bohr^-3 units
>>> WriteCube(atoms, density, 'filename.cube')

After the cube files have been created, one can run the bader code on the cube files with four possible options (1, 2, 3 or 4) where the choose of options is a choose of output files from the analysis as explained on the code's homepage. The bader code has been installed on all camp's desktop and will also be installed on niflheim in near future. Examples of using the program in the command line are:

bader filename.cube 4

where no Bader volumes are written at all. It is recommended to use this option first to see if the analysis is correct because writing the Bader volumes take a longer time. If one uses e.g.

bader filename.cube 3

then all the Bader volumes around each Bader maxima are written in the same cube file. If option 1 is used then all Bader volumes are written in separate files for each Bader region.

Few things should be mentioned for the people that wants to use this analysis tool;

  • The unitcell used in the calculations of the electronic structure has to be rectangular.
  • You can only create the cube files when your calculations are finished, because it is not until then that the charge density is written on the grid.
  • There can be a problem with O-H bonds (and possibly C-H bonds too) when the electron density is coming from pseudopotential code like dacapo.

For questions about the Bader algorithm in general should be directed to the code's discussion forum but questions of how to use the program with cube files written from NetCDFiles in e.g. Dacapo should be directed to