Crystals and Surfaces¶

Files¶

AMSinput can import and export coordinates in many different file formats, including .cif and .pdb files. Select File → Import coordinates or File → Export coordinates.

When exporting to .xyz format, AMS also writes the lattice vectors into the file.

Unit cells¶

In AMS, a 3D-periodic crystal can be represented as a periodic repetition of some unit cell. This is known as periodic boundary conditions (PBC). The unit cell is defined by three lattice vectors.

There are infinitely many choices of unit cells. The two most important ones are:

• The primitive cell, which is the smallest possible unit cell (containing only one lattice point), and
• The conventional cell, which follows guidelines outlined by the International Union for Crystallography.

For example, copper is a cubic closed-packed (ccp) metal with a face-centered cubic (fcc) lattice. The primitive and conventional cells are:

Fig. 7 Primitive unit cell of Cu.

Fig. 8 Conventional unit cell of Cu.

Download and import Cu_ccp_primitive.xyz into AMSinput

In AMSinput you can switch between primitive and conventional cells. Select → To Conventional Cell, or → To Primitive Cell. View the lattice vectors and lattice parameters under Model → Lattice.

All other possible unit cells can be generated as supercells from the primitive cell. Select → Generate Super Cell.

You can also view some of the periodic images of the system, without generating a supercell. Click the periodic view tool.

For many crystals, the primitive and conventional unit cells are identical.

When to use primitive cell, conventional cell, or supercell¶

• Always use the primitive cell for electronic or phonon band structure calculations. (For the calculation of phonon band structures, you should use the primitive cell but also specify the size of the supercell used internally by AMS).
• Always use the conventional cell when creating a surface (slab) from Miller indices.
• Use some suitable supercell when you do not have a perfect crystal, for example if you want to introduce a defect like a vacancy. The size of the supercell then determines the defect concentration. Supercells can also be used for approximately modeling amorphous solids, if atoms are displaced from their crystal positions.

Lattice vectors and lattice parameters¶

A unit cell in AMS is represented by three lattice vectors \(\vec{a}\), \(\vec{b}\), and \(\vec{c}\).

In crystallography it is more common to describe crystal structures with six lattice parameters: \(a, b, c, \alpha, \beta, \gamma\), where

• \(a, b, c\) are the lengths of the three lattice vectors
• \(\alpha\) is the angle between \(\vec{b}\) and \(\vec{c}\)
• \(\beta\) is the angle between \(\vec{a}\) and \(\vec{c}\)
• \(\gamma\) is the angle between \(\vec{a}\) and \(\vec{b}\)

In AMSinput, go to Model → Lattice to see the lattice vectors and lattice parameters (see Fig. 7). You can also edit them, but make sure to tick the Adjust atoms when changing lattice vectors box if you want the fractional coordinates to stay the same (if you want the atoms to move when you change the lattice).

In AMSinput, you can also select View → Properties to see the density.

Fractional coordinates¶

The positions of the atoms can either be expressed in

• Cartesian coordinates in angstrom (default in AMS), or
• Fractional coordinates (in the basis of the lattice vectors), which are common in crystallography

Toggle between Cartesian and fractional coordinates in AMSinput, by going to the Model → Coordinates panel, and tick the Use fractional coordinates checkbox.

Fig. 9 Fractional coordinates of Cu atoms in the conventional Cu unit cell.

You may also specify fractional coordinates in the text input file.

Rotation of lattice vectors relative to cartesian axes¶

The components of the lattice vectors are cartesian (xyz) coordinates in angstroms.

The lattice vectors can be rotated in an infinite number of ways relative to the cartesian axes. In AMS it does in general not matter how the lattice vectors are rotated. However, some other programs or file formats, which use the lattice parameters \(a, b, c, \alpha, \beta, \gamma\) may require that some convention is followed. Two common conventions are

• \(\vec{c}\) is parallel to z, and \(\vec{b}\) lies in the yz plane.
• \(\vec{a}\) is parallel to x, and \(\vec{b}\) lies in the xy plane.

In AMSinput, you can rotate the lattice vectors to follow the above conventions under Model → Lattice. Tick the Adjust atoms when changing lattice vectors, and press the Align c to z-axis or Align a to x-axis buttons.

Fig. 10 The lattice of ZnO_wurtzite.xyz

Fig. 11 The lattice after aligning \(\vec{a}\) to x. The fractional coordinates were kept fixed by ticking the “Adjust atoms” checkbox.

Run a lattice optimization¶

It is common in computational chemistry to optimize the lattice vectors of a crystal before constructing a supercell or surface. It is necessary to first optimize the lattice before calculating phonons or the elastic tensor.

In AMSinput, select Task → Geometry optimization on the main panel, and then on the Details → Geometry optimization panel tick the Optimize lattice checkbox.

You can also set the maximum number of iterations.

Convergence criteria for lattice optimizations¶

Four convergence criteria are considered:

• The maximum gradient (atomic force) component
• The root-mean-square of the gradient (atomic force) components
• The maximum displacement from one step to the next
• The maximum stress energy per atom

The first three are defined as for normal geometry optimizations.

The maximum stress energy per atom is used as a convergence criterion for the components of the stress tensor. It is the maximum value of the matrix stress_tensor * cell_volume / number_of_atoms.

In AMSinput, change the convergence criteria on the Details → Geometry Optimization panel.

Tips for lattice optimizations¶

• Have dense (good) k-point grids, if you use BAND, DFTB, Quantum ESPRESSO, or VASP. The quality of a k-point grid depends on the lattice, and during a lattice optimization the lattice changes. Therefore make sure that the quality is better than normal.
• BAND can be quite slow for lattice optimizations. Consider optimizing the lattice with Quantum ESPRESSO, and then switch back to BAND.
• To only allow isotropic scaling of the unit cell, go to Model → Geometry Constraints and PES Scan, and tick the off-diagonal (xy, xz, and yz) checkboxes under Freeze Strain and the diagonal (xx, yy, and zz) checkboxes under Equal Strain.
• To optimize the lattice under external pressure, go to Details → Geometry Optimization and set the Pressure.

Example: lattice optimization of wurtzite ZnO with DFTB¶

This example uses the znorg-0-1 DFTB parameter set. You can also run the example with ForceField (UFF), BAND, or Quantum ESPRESSO.

ZnO crystallizes in the hexagonal wurtzite structure, for which a = b, α = β = 90°, and γ = 120°.

1. Switch to the panel.

2. Import ZnO_wurtzite.xyz into AMSinput: Select File → Import coordinates.

3. Set Model to SCC-DFTB.

4. Set Parameter directory to DFTB.org/znorg-0-1

5. Set K-space to Good.

6. Set Task to Geometry Optimization, and click the arrow next to it.

7. Select Optimize lattice: Yes.

8. Run the job with File → Run, and save with the name zno_lattopt.ams.

9. When the job has finished, AMSinput will ask whether to update the coordinates. Select Yes, New job.

10. View the new lattice vectors and lattice parameters under Model → Lattice.

11. Export your structure to .xyz format: File → Export coordinates → .xyz. Save it with the name ZnO_optimized_znorg.xyz.

12. Optional: switch to AMSjobs (SCM → Jobs), select your job, and build a spreadsheet summary with Tools → Build spreadsheet.

Table 3 gives the a and c lattice parameters for different k-space qualities. This is an example of a k-point convergence test.

Surface Miller indices¶

Surfaces are often defined based on Miller indices. A surface plane is denoted (hkl).

The Miller indices give the inverse fractional intercepts along the lattice vectors in the conventional unit cell. If an index is 0, the surface lateral plane is parallel to that direction. A negative index is denoted by a bar.

Example 1: The Cu(100) surface is parallel to the second and third lattice vectors of Cu in the conventional unit cell.

Example 2: The Cu(111) surface is formed by the plane slicing all lattice vectors of the Cu conventional unit cell at the maximum coordinates for the unit cell.

Example 3: The Cu(511) surface is formed by the plane slicing the second and third lattice vectors of the Cu conventional unit cell at the maximum coordinates for the unit cell, but slicing the first vector at one fifth. This forms a stepped surface.

Fig. 12 Side views of slabs for different Cu surfaces.

For hexagonal crystals (e.g. Zn or ZnO), four Miller indices (hkil) are often used. The third index (i) is then redundant: h + k + i = 0.

For hexagonal crystals with a rhombohedral lattice (e.g. α-Al₂O₃ or calcite CaCO₃), the Miller indices can refer to either the rhombohedral primitive cell or the hexagonal conventional cell. You need to know which convention is used; when publishing or communicating with others, you should explicitly state which convention you use.

Crystal directions¶

A direction in a crystal is written as [hkl], and denotes the linear combination \(h\vec{a} + k\vec{b} + l\vec{c}\).

This direction is not necessarily perpendicular to the plane (hkl). Example: in a monoclinic cell with α = γ = 90° and β > 90°, the [001] direction is parallel to \(\vec{c}\) which is not perpendicular to the (001) plane; the (001) plane is parallel to both \(\vec{a}\) and \(\vec{b}\).

For hexagonal lattices, crystal directions can either be given using three indices as above, or a special four-index [hkil] notation. How to convert between the two types of notation is described in materials science or crystallography text books.

AMSinput and the other GUI modules do not show crystal directions.

Polar surface stabilization mechanisms¶

A polar surface can be stabilized by

• A change in stoichiometry, e.g., vacancies or other defects
• Adsorption
• charge transfer from one side of the slab to the other ↗, potentially leading to surface metallization. This is common in electronic structure calculations but uncommon in experiment.
• for thin slabs, the entire crystal structure might change ↗.

Example: ZnO(0001) and ZnO(\(000\bar{1}\))

ZnO(0001) can be stabilized by, for example

• 1/4 surface Zn vacancies, or
• 1/2 ML OH adsorption

ZnO(\(000\bar{1}\)) can be stabilized by, for example,

• 1/4 surface O vacancies, or
• 1/2 ML H adsorption

1. Download and extract polar_stabilization.zip

2. Import clean_good.xyz in AMSinput

3. On the panel, set Task → Single Point and Model → SCC-DFTB

4. Set KSpace to Good (surface metallization requires dense k-space sampling)

5. Set Parameter directory to DFTB.org/znorg-0-1.

6. Run the job (it may require a few hundred SCC cycles)

7. View results in AMSdos: SCM → Dos

8. To view the partial DOS from only a few atom, select the atoms in the AMSdos window

9. View the DOS for the top layer, a middle layer, and the bottom layer

10. Repeat steps 2-9 for vacancies_good.xyz and adsorbates_good.xyz

Fig. 18 a shows the DFTB-calculated DOS for a clean ZnO(0001) slab (terminated as in Fig. 13 d). The positions of the bands shift when going from the top to the bottom layer. The conduction band at the Zn-terminated side overlaps with the valence band on the O-terminated side, giving rise to charge transfer and surface metallization. Because of this, the SCC (or SCF for DFT calculations) can be difficult to converge.

Fig. 18 b and c show the DFTB-calculated DOS after vacancy or adsorption stabilization. In these cases, there is no surface metallization and the ZnO slab remains a semiconductor.

Fig. 18 DOS for ZnO(0001)/ZnO(000-1) calculated with SCC-DFTB znorg-0-1. a) clean bulk-terminated surfaces, b) surfaces stabilized by surface vacancies, c) surfaces stabilized by OH and H adsorption.

If you use VASP via AMS, you can use pseudo-atoms with a fractional number of electrons by selecting the corresponding POTCAR files. Such pseudo-atoms are sometimes used to passivate polar surfaces, but it is also possible to use real elements, like H, with some suitable stoichiometry.

Create a conventional surface unit cell¶

When you generate a slab in AMSinput, you may not always get the conventional surface unit cell. To get the conventional surface unit cell, select → To Conventional Cell. For the primitive surface unit cell, select → To Primitive Cell.

Example: The Cu(111) surface

1. Import Cu_ccp_primitive.yz

2. → To Conventional Cell. Note: this generates the conventional cell of the bulk Cu crystal, which is needed for the slab generation.

3. → Generate Slab

4. Type Miller indices 1, 1, 1

5. Number of layers: 5

6. Generate slab. The generated slab has 4 atoms in the top layer. (This is not the conventional surface cell)

7. → To Conventional Cell. The slab now has only 1 atom in the top layer. This is the conventional Cu(111) surface unit cell for a 5-layers-thick slab. See also Fig. 19

Create a supercell¶

From the conventional surface cell, it is common to create supercells, either to accommodate surface reconstructions or overlayers (adlayers, adsorbed layers) of some particular structure.

Supercells are often described using either Wood’s notation (not supported by AMSinput) or matrix notation (supported by AMSinput).

Fig. 19 Top views of the a) Cu(100) and b) Cu(111) surfaces, with only the top layer atoms shown as balls. Surface supercells of the primitive surface unit cells are given in matrix notation. The primitive (same as conventional) supercells are shown in blue.

Wood’s notation gives supercell’s lattice vector lengths as multiples of the original lattice lengths, together with any applied rotation. This notation can only be used if the angle γ between the lattice vector does not change when creating the supercell.

Examples of Wood’s notation:

• Cu(100) (2×1)
• Cu(100) (√2×√2)-R45°
• Cu(111) (√3×√3)-R30°

In AMSinput, you cannot specify supercells using Wood’s notation but must instead use matrix notation, which works for any supercell (potentially allowing the angle γ to change). See Fig. 19. The above three examples would translate to

• Cu(100) \(\begin{pmatrix}2 & 0 \\ 0 & 1\end{pmatrix}\)
• Cu(100) \(\begin{pmatrix}1 & 1 \\ 1 & -1\end{pmatrix}\)
• Cu(111) \(\begin{pmatrix}2 & 1 \\ -1 & 1\end{pmatrix}\)

In AMSinput, select → Generate supercell and enter the supercell matrix.

For hexagonal surfaces, like Cu(111) or ZnO(0001), γ = 120°. It is often more convenient to work with a supercell in which γ = 90°, especially if the purpose is to construct a large slab for a molecular dynamics simulation.

Convert hexagonal (γ = 120°) surface unit cell to rectangular (γ = 90°) supercell:

• Cu(111) \(\begin{pmatrix}1 & 0 \\ 1 & 2\end{pmatrix}\)

Example: Surface energy of ZnO(\(10\bar{1}0\))¶

This ZnO(\(10\bar{1}0\)) surface energy example uses DFTB and the ZnO_optimized_znorg.xyz file from the lattice optimization example. You can also run these calculations with BAND, ForceField (UFF), ReaxFF, Quantum ESPRESSO, or VASP via AMS. Use the same method that you used when generating ZnO_optimized_znorg.xyz.

1. Import ZnO_optimized_znorg.xyz into AMSinput

2. On the panel, set Model → SCC-DFTB

3. Set Parameter directory to DFTB.org/znorg-0-1

4. Set Kspace to Normal

5. Set Task to Geometry Optimization. On the Details → Geometry Optimization panel, make sure that Optimize lattice is not activated.

6. → Generate Slab

7. Miller indices 1, 0, 0 (for hexagonal crystals, ignore the third Miller index \(\bar{1}\)).

8. Number of layers: 4 (in this case the number of “double”-layers)

9. Select the Zn(2) atom, to get the correct surface termination (Fig. 13 b).

10. Press the Generate Slab button.

11. On the Model → Lattice panel, note the surface area \(|\vec{a} \times \vec{b}|\) = 17.68 Ų.

11. Save and run the job.

12. Repeat the above steps for a few other thicknesses (number of layers).

Now, tabulate the slab energies \(E^{\text{slab}}\) for each thickness. To get the energy, either

• SCM → Jobs, select the job, Job → Edit Comments, in the Results dropdown menu choose Energy, click Set as Default, click Save, select View → Comments.
• SCM → Output, use the bottom search box to search for calculation results. The energy is printed directly below.
• SCM → KFbrowser, File → Related files → dftb.rkf, AMS results → Energy.
• SCM → Jobs, select all jobs for the different thicknesses, Tools → Build spreadsheet, click the Do It button.

Using method 1, fit a straight line through the datapoints, as in Fig. 20.

Fig. 20 Energies of fully relaxed slabs of ZnO(\(10\bar{1}0\)) (red dots) calculated with SCC-DFTB znorg-0-1, k-space quality Normal. The number of (double-)layers is equal to 1/4 of the number of atoms. Only the three thickest slabs were used for the linear fit (blue line).

The intercept = \(2AE^\text{surf}\), where A = 17.68 Ų for all slabs. The surface energy becomes

If method 2 is used for the thickest slab with Normal k-space sampling, \(E^{\text{surf}}\) becomes 1.04 J/m². When reoptimizing the thickest slab with Good k-space sampling, \(E^{\text{surf}}\) instead becomes 1.31 J/m², which is the same number we found above using method 1.

Surface energy vs. cleavage energy¶

The cleavage energy is the energy required to create a slab from the bulk. If the slab exposes two identical surfaces, the cleavage energy is double the surface energy.

For example, you can only calculate the cleavage energy of the combined Zn-terminated ZnO(0001) and O-terminated ZnO(\(000\bar{1}\)) surfaces, and not the individual surface energies.

Adsorption on surfaces in AMSinput¶

You can place an adsorbate on a surface (substrate) in arbitrary conformations by rotating and translating the molecule.

Some surfaces have high-symmetry adsorption sites, like ontop, bridge, or hollow sites.

Fig. 21 Top view of water molecules adsorbed on different adsorption sites at Cu(111). The second Cu layer is colored green, the third is colored yellow.

Example: ontop, bridge, fcc hollow and hcp hollow sites on Cu(111)¶

Download Cu_111_4l_H2O.xyz, or create a 4-layer Cu(111) (5x5) surface supercell from the conventional Cu unit cell and add a water molecule somewhere on top of it.

For easier visualization, we will color the 2nd Cu layer green, and the 3rd Cu layer yellow.

1. In AMSinput, View → View Direction → Along x-axis

2. Atoms → Details (Color, Radius, Mass, …)

3. Select an atom in the 2nd Cu layer. On the right hand side, click the brown “Color” rectangle for this atom.

4. Change the color to green

5. Select all atoms in the 2nd Cu layer

6. Right-click on the green rectangle from before. This colors all atoms in the 2nd layer green.

7. Repeat steps 3-6 for the 3rd Cu layer, and color it yellow.

8. View → View Direction → Along z-axis

This is what it looks like from the side and from the top:

1. Edit → Tune Geometry

2. Select the water molecule and click the ⬭ label to the right of Atoms to move.

3. Select three atoms in the top Cu layer and click the ⬭ label to the right of Plane.

4. Set Distance to plane to 2.5 angstrom, and choose to closest atom.

5.

• Ontop: Select one Cu atom in the top layer

• Bridge: Select two neighboring Cu atoms in the top layer

• fcc hollow: Select three neighboring Cu atoms in the top layer that surround a yellow (third-layer) atom when viewed from above.

• hcp hollow: Select three neighboring Cu atoms in the top layer that surround a green (second-layer) atom when viewed from above.

6. Click the ⬭ label to the right of Above atoms.

7. Click the Move button. This moves the adsorbate to the chosen adsorption site.

Fig. 22 Example setup for moving an H₂O adsorbate to the bridge position on Cu(111).

Adsorption energy¶

The adsorption energy \(\Delta E^\text{ads}\) is the binding energy of a molecule (adsorbate) on a surface:

All energies on the right hand side are calculated for geometry-optimized systems. With the above equation, a positive adsorption energy implies exothermic adsorption. A different sign convention can also be used:

With this equation, a negative adsorption energy implies exothermic adsorption.

When reporting adsorption energies, always specify the equation that you used to calculate them.

Fig. 23 a-c) The slab, mol, and slab+mol systems needed to calculate the adsorption energy. d) The slab and molecule widely separated, used for the in-cell approach

Basis set superposition error with BAND¶

In BAND, if your basis set is not large enough, the adsorption may become artificially stable from the basis set superposition error (BSSE). The best way to address BSSE is to increase the size of the basis set. Alternatively, you can also apply the counterpoise correction.

DFTB, Quantum ESPRESSO, VASP via AMS, and the other engines do not show any BSSE.

Example: water molecule adsorption on ZnO(\(10\bar{1}0\)).

The figure below shows calculated adsorption energies with and without counterpoise correction. Each geometry was optimized with the corresponding basis set, and the counterpoise correction was applied to final slab+mol system.

Fig. 24 Adsorption energies calculated with and without counterpoise correction using BAND with the PBE functional. A more negative value implies more exothermic adsorption. Download PLAMS scripts that set up and postprocess jobs calculating the counterpoise correction and which generate this figure.

In-cell approach for adsorption energies with ReaxFF¶

In ReaxFF, the electronegativity equalization method (EEM) is used to determine the charges on the atoms, which in turn determine the Coulomb energy contribution to the total energy.

In an isolated molecule, charge can only be redistributed among the atoms in the molecule. Similarly, for an isolated slab, charge can only be redistributed among the atoms in the slab.

In the combined slab+molecule system, there is more freedom to redistribute the charge, potentially lowering the energy and giving more stable adsorption. This issue is similar to the basis set superposition error.

For ReaxFF calculations, it is therefore an option to use the in-cell approach:

where the first term is calculated for a system where the slab and molecule are placed far apart in the same cell (“in-cell”), see Fig. 23 d. In this way, for both systems on the right hand side, the charges have the same freedom to redistribute.

A different sign convention can also be used:

Example: H₂O adsorption on ZnO(\(10\bar{1}0\)). (3×2) supercell, 5 double-layers thick, constructed from bulk lattice parameters a = 3.28 Å, c = 5.28 Å. Download a PLAMS python script calculating these numbers.

The difference in adsorption energy between the two approaches is 10 kcal/mol, which stems from the fact that

This energy difference gets larger for larger supercells and thicker slabs, see Fig. 25

Fig. 25 Energy difference between the “standard” and “in-cell” approaches for calculating adsorption energies. The data was calculated from single-point calculations for a 5 double-layer thick slab of ZnO(\(10\bar{1}0\)) using the force field from Raymand et al. (ZnOH.ff) Download a PLAMS python script and xyz file to generate this figure.

One-sided vs. two-sided adsorption¶

A slab exposes two surfaces, and adsorbates can either be added to both surfaces or only one.

For Quantum ESPRESSO and VASP, double-sided adsorption can prevent the formation of a dipole moment across the vacuum gap. The other engines support 2D periodicity for which there is no vacuum gap.

Fixing the central and/or bottom layers of the slab¶

It is common to fix the positions of either the central or the bottom layers of the slab.

When relaxing a clean (without adsorbates) slab, it is a good idea to relax either all atoms, or to keep the central layers fixed. In this way, both the top and bottom sides of the slab relax. This relaxed slab can then be used for

• one-sided adsorption (relaxing everything or fixing the central and bottom layers),
• double-sided adsorption (relaxing everything or fixing the central layers).

Even in the case of one-sided adsorption (on the top surface) it is good to have a relaxed bottom surface when using DFTB, BAND, Quantum ESPRESSO, or VASP. This will help SCF convergence and give a “cleaner” DOS (compare Fig. 16 and Fig. 17), especially for covalent or ionic materials. For metals the effect is smaller.

In AMSinput, select the atoms that you would like to keep fixed. On the Model → Geometry Constrained and PES Scan panel, click the plus button next to selected atoms (fix positions).

Example: See the TiO2 surface relaxation tutorial.

Solid-solid interfaces¶

A solid/solid interface (heterojunction) under periodic boundary conditions corresponds to a slab of material A adsorbed onto a slab of material B.

A crystalline interface can be written as A(hkl) || B(h’k’l’) with A[uvw] || B[u’v’w’], where the first part specifies which crystal planes that are parallel, and the second part specifies a crystal direction in A in the (hkl) plane that is parallel to a crystal direction in B in the (h’k’l’) plane.

The lattice vectors of the slab models of the two materials must match exactly (coincidence site lattice). Two materials usually have very different crystal structures, so to create such a model, supercells that are possibly strained or rotated need to be created, before adsorbing one slab on top of the other.

Solid-liquid interfaces¶

Solid/liquid interfaces can easily be created in AMSinput. Create a slab, and select Edit → Builder.

For more details, see the Water on an aluminum surface tutorial.