### Bayesian Error Estimate (BEE)

For a given observable calculated within DFT, say an atomization energy of a molecule, the
variation of the calculated value of this observable within an
ensemble of density functionals provides a useful estimate of how
large the error of the value is compared to experiment. The
`ASE.Utilities.BEE` module contains a finely tuned ensemble of such
functionals. Suppose you want to estimate the errorbar on the
dissociation energy of molecule AB, and you have done calculations for
systems A, B and AB, then the following script will do the job:

from Numeric import sum from Dacapo import Calculator from ASE.Utilities.BEE import GetEnsembleEnergies # Read in calculations for A, B and AB: a = Calculator.ReadAtoms('A.nc') b = Calculator.ReadAtoms('B.nc') ab = Calculator.ReadAtoms('AB.nc') # Get ensembles of energies: ea = GetEnsembleEnergies(a) eb = GetEnsembleEnergies(b) eab = GetEnsembleEnergies(ab) # ea, eb and eab are now arrays of 1000 energies # Calculate ensemble of binding energies: e = ea + eb - eab n = len(e) e0 = sum(e) / n sigma = (sum((e - e0)**2) / n)**0.5 print 'Best fit:', e0, '+-', sigma, 'eV'

The `GetEnsembleEnergies` function takes a `ListOfAtoms`-argument
and returns a Numeric array of length 1000 containing energies from
the ensemble. The technique is described in more detail here.

The observable that you want to calculate the BEE for must be a
function of a number of total energies for a set of different systems,
like in the example above (`ea + eb - eab`).
Adsorption energies, cohesive energies and all other energy
differences are good examples, but also bond lengths and lattice
constanst are possible observables.

Note

The calculation of a BEE is based on non-selfconstistent PBE calculations When the energy is evaluated for all points in the ensemble of xc-functionals the electron density is fixed, and so are the eigenvalues and the effective potential. It is therefore not possible to calculate the BEE for quantities like band gaps and work functions.

#### BEE for hydrogen bond length

In this example we calculate the bond length of a hydrogen molecule by doing three calculations at three different bond lengths:

import Numeric as num from gpaw import Calculator from ASE import ListOfAtoms, Atom from ASE.Utilities.BEE import GetEnsembleEnergies a = 5.0 # Size of unit cell (Angstrom) h2 = ListOfAtoms([Atom('H'), Atom('H')], cell=(a, a, a)) h2.SetCalculator(Calculator(h=0.2, xc='PBE', out=None)) b0 = 0.74 # Experimental bond length db = 0.05 e_pbe = [] e_bee = [] for b in [b0 - db, b0, b0 + db]: h2[0].SetCartesianPosition([a/2, a/2, a/2 - b/2]) h2[1].SetCartesianPosition([a/2, a/2, a/2 + b/2]) e_pbe.append(h2.GetPotentialEnergy()) e_bee.append(GetEnsembleEnergies(h2)) def f(e): """Fit to a parabola.""" # e(b) = a0 + a1*(b - b0) + 0.5*a2*(b - b0)^2 a1 = (e[2] - e[0]) / (2 * db) a2 = (e[0] - 2 * e[1] + e[2]) / db**2 return b0 - a1 / a2 b_pbe = f(e_pbe) b_bee = f(e_bee) n = len(b_bee) b_ave = num.sum(b_bee) / n sigma = (num.sum((b_bee - b_ave)**2) / n)**0.5 print 'b =', b_pbe, '+-', sigma, 'Angstrom'

### The GeometricTransforms module

The geometric transforms module contains functions that allow to translate and rotate a ListOfAtoms object. In general, the functions do not return anything and change the coordinates of the ListOfAtoms that is given to them. For example:

>>> from ASE import Atom, ListOfAtoms >>> atoms=ListOfAtoms([Atom('Mo', (1.0, 1.0, 1.0)), ... Atom('Mo', (-2.0, -2.0, -2.0))], ... cell=(10,10,10)) >>> from ASE.Utilities.GeometricTransforms import * >>> Translate(atoms, (1.0, 1.0, 1.0), 'cartesian')

This translates `atoms` about the vector `(1.0,1.0,1.0)` in cartesian coordinates,
use the argument `scaled` to do this in scaled coordinates. The function changes the
coordinates in `atoms`, if you want to keep an original unchanged ListOfAtoms, copy
it with

>>> newatoms=atoms.Copy()

before applying the transformations. All functions take a ListOfAtoms (specified by `atoms`)
as the first argument.

`CenterOfMass(atoms,mode='scaled')`:- Returns the center of mass of
`atoms`. Default are scaled coordinates, but the argument`cartesian`specifies that the center of mass is returned in cartesian coordinates. `Translate(atoms, trans_vector=(0.0,0.0,0.0), mode='scaled')`:- Translates
`atoms`about the vector`trans_vector`. The argument`scaled`specifies that`trans_vector`is in scaled coordinates, with`cartesian`it is assumed that the vector is in cartesian coordinates. `TranslateCOM(ListOfAtoms, target_vector=(0.0,0.0,0.0), mode='scaled'):`:- Translates the center of mass of
`atoms`to a point specified by`target_vector`. The argument`scaled`indicates that`target_vector`is assumed to be in scaled coordinates, with`cartesian`, it is assumed to be in cartesian coordinates. `TranslateAtom(ListOfAtoms,atom,target_vector=(0.0,0.0,0.0),mode='scaled')`:- Translates
`atoms`such that the atom specified by`atom`is located at`target_vector`. The argument`scaled`or`cartesian`specifies whether`target_vector`is given in scaled or cartesian coordinates. The atom`atom`is specified by its number in the ListOfAtoms, the counting starts at 0. Example:

>>> from ASE import Atom, ListOfAtoms >>> atoms=ListOfAtoms([Atom('Mo', (1.0, 1.0, 1.0)), ... Atom('Mo', (-2.0, -2.0, -2.0))], ... cell=(10,10,10)) >>> from ASE.Utilities.GeometricTransforms import * >>> TranslateAtom(atoms,1, (1.0, 1.0, 1.0), 'cartesian')

This example translates `atoms` such that atom number 1 (the second atom in the list) is
located at `(1.0,1.0,1.0)` in cartesian coordinates.

`RotateAboutAxis(atoms,axis=(1,0,0),theta=0):`:- Rotates the list of atoms about an axis specified by
`axis`and the origin by an angle specified by theta (in degrees). As the axis goes through the origin and the point`axis`you have to take care of that the atoms are centered around an appropriate origin, e.g. by translating them with the functions above. `RotateCOMAboutX(atoms,theta):`:- Rotates the list of atoms around the x axis through its center of mass by an angle
specified by
`theta`(in degrees). The function works such that first the center of mass is translated to the origin, then the rotation around the x axis is performed and then the center of mass is translated back. `RotateCOMAboutY(atoms,theta):`:- Rotates the list of atoms around the y axis through its center of mass, see function above.
`RotateCOMAboutZ(atoms,theta):`:- Rotates the list of atoms around the z axis through its center of mass, see function above.
`RotateUnitCellAboutX(atoms,theta):`:- Rotates the unit cell about the x axis by an angle specified by
`theta`in degrees. The positions of the atoms are changed with the unit cell so that their relative positions are the same as before. `RotateUnitCellAboutY(atoms,theta):`:- Rotates the unit cell about the y axis, see function above.
`RotateUnitCellAboutZ(atoms,theta):`:- Rotates the unit cell about the z axis, see function above.
`MoveAllAtomsIntoUnitCell(atoms):`:- Translates all atoms by unit cell vectors until they are all in the unit cell.
`RattleLOA(atoms,amount=0.01):`:- Changes the atomic positions by a small random amount within the interval specified by
`amount`. This function is designed to give some random shifts in atomic positions to break any symmetry present. `GetStrainTensor(strainvec):`:- Given a strain vector by
`strainvec`, this function returns the corresponding strain tensor. `ApplyStrain(atoms,epsilon):`:- Strains the unit cell of atoms according to the strain vector
`epsilon`. The strain vector`epsilon`is specified by`[exx, eyy, ezz, gamma_yz, gamma_xz, gamma_xy]`where gamma is the engineering strain. i.e.`e_xy = 0.5*gamma_xy`. The atoms are moved so that their relative positions are the same as before. `SetUnitCellVolume(ListOfAtoms,volume):`:- Scales the unit cell volume to the number given by
`volume`.

Apart from the functions the module contains some self-tests. These are probably useful when adding new functions.

### The EquationOfState module

This module provides tools to fit volumetric equations of state to calculate volumes and energies. Several equations of state are supported, including Murnaghan, Birch, BirchMurnaghan, Vinet, PoirerTarantola, AntonSchmidt and Taylor series expansion.

>>> from ASE.Utilities.EquationOfState import EquationOfState >>> volumes = [vol1, vol2, vol3, vol4, vol5] >>> energies = [E1, E2, E3, E4, E5] >>> eos = EquationOfState('Murnaghan',volumes,energies) >>> print eos >>> g = eos.GetPlot() >>> eos.SavePlot('murn.png') #save the figure as a png >>> print eos.GetPressure(somevolume) >>> print eos.GetBulkModulus(somevolume) >>> print eos.GetEnergy(somevolume) >>> print eos.GetReference()

This fits the Murnaghan equation of state to the volumes and energies provided. It then prints the volume, bulk modulus (in GPa), energy and pressure at the minumum found. A gnuplot window will pop up with the fit and will be saved as a png format. You can also find the pressure, bulk modulus and energy at other volumes if desired. Finally, you can print the reference where the equation came from. Typical output looks like:

Murnaghan: PRB 28, 5480 (1983) V0(A^3) B(Gpa) E0 (eV) pressure(GPa) 16.60 78.27 -56.4688 -0.00 chisq = 0.0000

A typical plot might look like:

`eos = EquationOfState(eos,volumes,energies):`:- creates an
EquationOfState class,
`volumes`and`energies`should each be a list of numbers.`eos`must be one of these strings "Murnaghan", "Birch", "BirchMurnaghan", "Vinet", "PoirerTarantola", "AntonSchmidt", "Taylor" `print eos:`:- prints a string containing all the fit information
`eos.GetReference():`:- returns literature reference for equation used.
`eos.GetV0():`:- returns volume of minimum energy
`eos.GetEnergy(volume):`:- returns the fitted energy at
`volume` `eos.GetPressure(volume):`:- returns the pressure at
`volume` `eos.GetBulkmodulus(volume):`:- returns the bulk modulus at
`volume` `eos.GetPlot():`:- pops up a gnuplot of raw data, fit and results. Returns the gnuplot object.
`eos.SavePlot(filename):`:- saves the gnuplot in
`filename`. Format is determined by the filename extension, png, eps, tex, and pdf(requires new gnuplot and libpdf) are supported.

### Crystal structures

Modules for creating crystal structures are found in
`ASE.Utilities.Lattice`. Most Bravais lattices are implemented, as
are a few important lattice with a basis. The modules can create
lattices with any orientation (see below).

#### Available crystal lattices

The following modules are currently available (the * mark lattices with a basis):

`ASE.Utilities.Lattice.Cubic``SimpleCubic``FaceCenteredCubic``BodyCenteredCubic``Diamond`(*)

`ASE.Utilities.Lattice.Tetragonal``SimpleTetragonal``CenteredTetragonal`

`ASE.Utilities.Lattice.Orthorhombic``SimpleOrthorhombic``BaseCenteredOrthorhombic``FaceCenteredOrthorhombic``BodyCenteredOrthorhombic`

`ASE.Utilities.Lattice.Monoclinic``SimpleMonoclinic``BaseCenteredMonoclinic`

`ASE.Utilities.Lattice.Triclinic``Triclinic`

`ASE.Utilities.Lattice.Hexagonal``Hexagonal``HexagonalClosedPacked`(*)`Graphite`(*)

The rhombohedral (or trigonal) lattices are not implemented. They will be implemented when the need arises (and if somebody can tell me the precise definition of the 4-number Miller indices - I only know that they are "almost the same as in hexagonal lattices").

`ASE.Utilities.Lattice.Compounds`Lattices with more than one element. These are mainly intended as examples allowing you to define new such lattices. Currenly, the following are defined

`B1`=`NaC``l = ``Rocksalt``B2`=`CsCl``B3`=`ZnS`=`Zincblende``L1_2`=`AuCu3``L1_0`=`AuCu`

Warning

Asap users should import from `Asap.Setup.Lattice.XXX` instead of
from `ASE.Utilities.Lattice.XXX` in order to produce Asap list of
atoms.

#### Usage

The lattice objects are called with a number of aguments specifying e.g. the size and orientation of the lattice. All arguments should be given as named arguments.

`directions`and/or`miller`:Specifies the orientation of the lattice as the Miller indices of the three basis vectors of the supercell (

`directions=...`) and/or as the Miller indices of the three surfaces (`miller=...`). Normally, one will specify either three directions or three surfaces, but any combination that is both complete and consistent is allowed, e.g. two directions and two surface miller indices (this example is slightly redundant, and consistency will be checked). If only some directions/miller indices are specified, the remaining should be given as`None`. If you intend to generate a specific surface, and prefer to specify the miller indices of the unit cell basis (`directions=...`), it is a good idea to give the desired Miller index of the surface as well to allow the module to test for consistency. Example:>>> atoms = BodyCenteredCubic(directions=[[1,-1,0],[1,1,-1],[0,0,1]], ... miller=[None, None, [1,1,2]], ...)

`size`:- A tuple of three numbers, defining how many times the fundamental repeat unit is repeated. Default: (1,1,1). Be aware that if high-index directions are specified, the fundamental repeat unit may be large.
`element`:- The element, specified by the atomic number (an integer), by the atomic symbol (i.e. "Au") or by an object returned by ASE.ChemicalElements.Element(). For compounds, a tuple or list of elements should be given.
`latticeconstant`:The lattice constant. If no lattice constant is specified, one is extracted from ASE.ChemicalElements provided that the element actually has the crystal structure you are creating. Depending on the crystal structure, there will be more than one lattice constant, and they are specified by giving a dictionary or a tuple (a scalar for cubic lattices). Distances are given in Angstrom, angles in degrees.

Structure Lattice constants Cubic 'a'

a

Tetragonal 'a', 'c' or 'c/a'

(a, c)

Orthorhombic 'a', 'b' or 'b/a', 'c' or 'c/a'

(a, b, c)

Triclinic 'a', 'b' or 'b/a', 'c' or 'c/a', 'alpha', 'beta', 'gamma'

(a, b, c, alpha, beta, gamma)

Monoclinic 'a', 'b' or 'b/a', 'c' or 'c/a', 'alpha'

(a, b, c, alpha)

Hexagonal 'a', 'c' or 'c/a'

(a, c)

Example:

>>> atoms = Monoclinic( ... , latticeconstant={'a':3.06, ... 'b/a': 0.95, 'c/a': 1.07, 'alpha'=74})

`debug`:- Controls the amount of information printed. 0: no info is printed. 1 (the default): The indices of surfaces and unit cell vectors are printed. 2: Debugging info is printed.

#### Defining new lattices

Often, there is a need for new lattices - either because an element crystallizes in a lattice that is not a simple Bravais lattice, or because you need to work with a compound or an ordered alloy.

All the lattice generating objects are instances of a class, you generate new lattices by deriving a new class and instantiating it. This is best explained by an example. The diamond lattice is two interlacing FCC lattices, so it can be seen as a face-centered cubic lattice with a two-atom basis. The Diamond object could be defined like this:

from ASE.Utilities.Lattice.Cubic import FaceCenteredCubicFactory, \ MakeListOfAtoms class DiamondFactory(FaceCenteredCubicFactory): "A factory for creating diamond lattices." xtal_name = "diamond" bravais_basis = [[0,0,0], [0.25, 0.25, 0.25]] Diamond = DiamondFactory(MakeListOfAtoms)

The `MakeListOfAtoms` is available from all the crystal structure
module. It is needed to support ASE and Asap list of atoms in the
same module. Asap users should import from `Asap.Setup.Lattice.XXX`
and should replace `MakeListOfAtoms` with `MakeAsapListOfAtoms` in
order to produce the right kind of list of atoms.

##### Lattices with more than one element

Lattices with more than one element is made in the same way. A new
attribute, `element_basis`, is added, giving which atoms in the
basis are which element. If there are four atoms in the basis, and
element_basis is (0,0,1,0), then the first, second and fourth atoms
are one element, and the third is the other element. As an example,
the AuCu3 structure (also known as L1_2) is defined as:

# The L1_2 structure is "based on FCC", but is really simple cubic # with a basis. class AuCu3Factory(SimpleCubicFactory): "A factory for creating AuCu3 (L1_2) lattices." bravais_basis = [[0, 0, 0], [0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0]] element_basis = (0, 1, 1, 1) AuCu3 = L1_2 = AuCu3Factory(MakeListOfAtoms)

Sometimes, more than one crystal structure can be used to define the crystal structure, for example the Rocksalt structure is two interpenetrating FCC lattices, one with one kind of atoms and one with another. It would be tempting to define it as

class NaClFactory(FaceCenteredCubicFactory): "A factory for creating NaCl (B1, Rocksalt) lattices." bravais_basis = [[0, 0, 0], [0.5, 0.5, 0.5]] element_basis = (0, 1) B1 = NaCl = Rocksalt = NaClFactory(MakeListOfAtoms)

but if this is used to define a finite system, one surface would be covered with one type of atoms, and the opposite surface with the other. To maintain the stochiometry of the surfaces, it is better to use the simple cubic lattice with a larger basis:

# To prevent a layer of element one on one side, and a layer of # element two on the other side, NaCl is based on SimpleCubic instead # of on FaceCenteredCubic class NaClFactory(SimpleCubicFactory): "A factory for creating NaCl (B1, Rocksalt) lattices." bravais_basis = [[0, 0, 0], [0, 0, 0.5], [0, 0.5, 0], [0, 0.5, 0.5], [0.5, 0, 0], [0.5, 0, 0.5], [0.5, 0.5, 0], [0.5, 0.5, 0.5]] element_basis = (0, 1, 1, 0, 1, 0, 0, 1) B1 = NaCl = Rocksalt = NaClFactory(MakeListOfAtoms)

More examples can be found in the file
`ASE/Utilities/Lattice/Compounds.py`.