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:

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.

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.

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.

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.

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 will get the conventional surface unit cell. For the primitive surface unit cell, select → To Primitive Cell. To get the conventional surface unit cell, select → To Conventional 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.

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).

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.

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.

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.

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

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.

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.

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

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 VASP: TiO₂ 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 Solid-liquid interface tutorial.

To build solid/liquid interfaces with Python scripting, see the PLAMS Packmol example.