Crystals and Surfaces
*********************

BAND Quantum ESPRESSO DFTB Materials Science Crystal Surface Supercell
AMSinput

This tutorial covers, for **crystals**:

* Periodic crystals in AMSinput

* k-space sampling and convergence

* Lattice optimizations, with an example

* Overview of crystal properties in AMS

And for **surfaces**:

* Slabs vs. QM/MM in AMS

* Miller indices

* Create a surface slab in AMSinput

* Surface passivation, polar and nonpolar surfaces

* Surface relaxation, surface reconstruction

* Surface energy,

* Surface supercells, including matrix notation

* Adsorption sites and Adsorption energy calculations

Note:

  If you're not familiar with AMSinput and AMSjobs, first go through
  the Getting Started with AMS and Building Crystals and Slabs
  tutorials, as well as some tutorial for your compute engine, like

  * DFTB: Lattice Optimization and Phonons

  * BAND: DFT + Hubbard U, PDOS

  * Quantum ESPRESSO: Quantum ESPRESSO: Geometry and Lattice
    Optimization

  * VASP: VASP: TiO₂ surface relaxation

  * ReaxFF: Solid-liquid interface

  To run the PLAMS python scripts, first go through the Run Jobs with
  Python Scripts tutorial.In this tutorial, Zn atoms are colored brown
  for improved contrast (the default color is white). Learn how to
  change the default color of atoms.


3D-periodic crystals
====================

This tutorial covers only **3D-periodic crystals** and their surfaces.
Quasicrystals (aperiodic crystals) are not covered.


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:

+---------------------------+---------------------------+---------------------------+---------------------------+
| Cu unit cell              | # lattice points          | # atoms                   | angles between lattice    |
|                           |                           |                           | vectors                   |
|===========================|===========================|===========================|===========================|
| primitive                 | 1                         | 1                         | non-orthogonal            |
+---------------------------+---------------------------+---------------------------+---------------------------+
| conventional              | 4                         | 4                         | orthogonal (90°)          |
+---------------------------+---------------------------+---------------------------+---------------------------+

   [image]Primitive unit cell of Cu.

   [image]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 [image: CrystalTool][image] **→ To Conventional Cell**,
or [image: CrystalTool][image] **→ 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 [image: CrystalTool][image] **→ Generate Super
Cell**.

You can also view some of the periodic images of the system, without
generating a supercell. Click the [image: PeriodicViewTool][image]
periodic view tool.

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

Note:

  The conventional Cu unit cell contains 4 lattice points in a
  *crystallographic* sense. AMS does not recognize that all 4 Cu atoms
  at the lattice points are identical, but would treat them
  independently (at an increased computational cost compared to the
  primitive cell).

Tip:

  In AMSinput, select **View → Periodic → Show Lattice Vectors** and
  **Show Lattice Vector Labels** to see the lattice vectors. Select
  **View → Periodic → Show Unit Cell** to highlight the unit cell.


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. %s). 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 **Properties → Properties (including
estimated)** to see the **density**.

Tip:

  Whenever you change the lattice vectors or lattice parameters, use
  the [image: PeriodicViewTool][image] periodic view tool to double-
  check that the periodic boundary conditions are what you expect.


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.

   [image]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.

The below figures show the hexagonal unit cell of "wurtzite ZnO"

   [image]The lattice of "ZnO_wurtzite.xyz"

   [image]The lattice after aligning \vec{a} to *x*. The fractional
   coordinates were kept fixed by ticking the "Adjust atoms" checkbox.


k-space sampling
================

With DFTB, BAND, Quantum ESPRESSO, and VASP via AMS, you must specify
some k-point sampling for periodic calculations. ReaxFF, ForceField,
and MLPotential do not use k-points.

In **DFTB and BAND**, you can specify the k-space quality using the
values "GammaOnly", "Basic", "Normal", "Good", "VeryGood", and
"Excellent". You can also manually specify the number of k-points
along each reciprocal lattice vector, or for highly symmetric systems
use a Symmetric k-space grid.

In **Quantum ESPRESSO** and **VASP via AMS**, you must manually
specify the dimensions of the Monkhorst-Pack grid.

The required amount of k-points depends on the material and size of
the unit cell:

* **Metals** require many k-points

* **Insulators, semiconductors**, and **molecular crystals** typically
  require fewer k-points

* **Small unit cells** require many k-points

* **Large unit cells** and **supercells** require fewer k-points

To determine suitable k-point settings for your method and system, you
should do a **k-point convergence study** in which you systematically
increase the quality of the k-space sampling until the quantity that
you are interested in (e.g., energy or lattice parameter) is
converged. See Example: lattice optimization of wurtzite ZnO with DFTB
for an example.

Note:

  The "NumericalQuality" keyword for BAND affects the k-space
  sampling.

Tip:

  **If possible, avoid comparing energies calculated with nonidentical
  k-space sampling**. For example, if you create a supercell from the
  primitive cell and introduce an interstitial defect, calculate the
  insertion energy with respect to the energy of the nondefective
  *supercell*, not with respect to a multiple of the energy of the
  primitive cell. See also: Surface energy calculations.

See also: Recommendations for k-space


Lattice optimization
====================


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 [image: DFTBPanel][image] 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 [image: MoreBtn][image] 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, click the [image: UpdateMolecule][image] button in AMSinput.
   **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 %s gives the *a* and *c* lattice parameters for different
k-space qualities. This is an example of a **k-point convergence
test**.


Optimized lattice parameters for wurtzite ZnO with SCC-DFTB znorg-0-1.
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

+-----------------------------------+-----------------------------------+-----------------------------------+
| k-space                           | *a* (Å)                           | *c* (Å)                           |
|===================================|===================================|===================================|
| Basic                             | 3.281                             | 5.425                             |
+-----------------------------------+-----------------------------------+-----------------------------------+
| Normal                            | 3.295                             | 5.347                             |
+-----------------------------------+-----------------------------------+-----------------------------------+
| Good                              | 3.290                             | 5.375                             |
+-----------------------------------+-----------------------------------+-----------------------------------+
| Very Good                         | 3.289                             | 5.375                             |
+-----------------------------------+-----------------------------------+-----------------------------------+

Tip:

  k-point convergence tests can be scripted with the PLAMS Python
  module. Download a "PLAMS python script" and "xyz file" to generate
  the above table.


Solid-state properties
======================

* **Band gap** with BAND/DFTB/QE: The one-electron gap is printed
  towards the end of the logfile for DFTB and BAND calculations. The
  one-electron gap can also be extracted from band structure
  calculations, which include all high-symmetry k-points.

* **Electronic density of states (DOS)** with BAND/DFTB/QE: On the
  Main panel, tick the *Calculate DOS* checkbox. This also calculates
  the partial DOS, PDOS.

* **Electronic band structure** with BAND/DFTB/QE: On the Main panel,
  tick the *Calculate band structure* checkbox. Band structures are
  usually calculated for the primitive cell. To create a custom path
  through the Brillouin zone, go to *Properties → Band Structure*.

* **Fermi surface** with BAND: In AMSInput, go to *Properties → Band
  Structure* and tick *Calculate Fermi surface*. The calculated Fermi
  surface can be visualized with the Band Structure GUI module (SCM →
  Band Structure, then menu *View → Brillouin Zone*). See also the
  BAND documentation on Fermi surface.

* **Phonon dispersion curves and phonon DOS**: *Properties → Phonons
  and Elastic tensor → Phonons*. First optimize the lattice! See the
  Silicon lattice optimization and phonons tutorial.

* **Bulk modulus, Shear modulus, Young's modulus, Elastic tensor**:
  *Properties → Phonons and Elastic tensor → Elastic tensor*. First
  optimize the lattice with tight convergence criteria.

* **Stress-strain curve, yield point**: See the Snapping Polyacetylene
  Chain tutorial.

* **Thermal expansion coefficient**: Run well-equilibrated NPT
  molecular dynamics simulations at several different temperatures:
  See the Thermal expansion coefficients of thermoset polymers
  tutorial.

* **Atomization energy** and **lattice energy**: See Bond energy
  calculations.


Surface models in computational chemistry
=========================================

There are two main ways of representing surfaces in atomistic
modeling:

* A **slab model**, which is periodic in the lateral directions. It
  can also be periodic in the surface normal direction, in which case
  a **vacuum gap** separates the slabs (a repeated slab model). The
  AMS engines **BAND, DFTB,** and **ReaxFF** can be modeled with
  2-dimensional periodicity and do not require a vacuum gap. **Quantum
  ESPRESSO**  and **VASP via AMS** require a vacuum gap (3-dimensional
  periodicity).

* An **embedded cluster (QM/MM)**, in which a piece of the surface is
  modeled at a high level of theory (e.g. DFT), and the rest at a
  lower level of theory (force field or point charges). **ADF, BAND,
  and DFTB** can be used for this type of surface model. Combining
  different engines in QM/MM is done via the **Hybrid** engine.

The rest of this tutorial only covers **slab models**, which are the
most common type of surface model. For more details about QM/MM, see
the Inorganic linker in organic framework tutorial.

Tip:

  AMS includes a kinetics module for modeling reactions on surfaces.

+---------------------------+---------------------------+---------------------------+---------------------------+
| Engine                    | 2D slab                   | 3D slab (vacuum gap)      | QM/MM                     |
|===========================|===========================|===========================|===========================|
| ADF                       | ❌                        | ❌                        | ✔ (QM)                    |
+---------------------------+---------------------------+---------------------------+---------------------------+
| BAND                      | ✔                         | ✔                         | ✔ (QM)                    |
+---------------------------+---------------------------+---------------------------+---------------------------+
| DFTB                      | ✔                         | ✔                         | ✔ (QM, MM)                |
+---------------------------+---------------------------+---------------------------+---------------------------+
| ForceField                | ✔                         | ✔                         | ✔ (MM)                    |
+---------------------------+---------------------------+---------------------------+---------------------------+
| ReaxFF                    | ✔                         | ✔                         | ✔ (MM)                    |
+---------------------------+---------------------------+---------------------------+---------------------------+
| ML Potential              | ❌                        | ✔                         | ❌                        |
+---------------------------+---------------------------+---------------------------+---------------------------+
| Quantum ESPRESSO          | ❌                        | ✔                         | ❌                        |
+---------------------------+---------------------------+---------------------------+---------------------------+
| VASP via AMS              | ❌                        | ✔                         | ❌                        |
+---------------------------+---------------------------+---------------------------+---------------------------+


Miller indices
==============


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

   [image]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_2O_3
or calcite CaCO_3), 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.

Note:

  Crystal directions are almost always given with respect to the
  **conventional cell**.


Generate a slab with AMSinput
=============================

   **1.** Download "Cu_ccp_primitive.xyz", **File → Import coordinates**
   **2.** [image: CrystalTool][image] **→ Convert To Conventional Cell**
   **3.** [image: CrystalTool][image] **→ Generate Slab**
   **4.** Miller indices "(5, 1, 1)"
   **5.** Number of layers: 10
   **6.** **Generate slab**

This will generate a Cu(511) slab as in Fig. %s. See also the Building
Crystals and Slabs tutorial.

The slab is rotated such that the **surface normal direction** is
parallel to *z*.

If you have the [image: BANDPanel][image], [image: DFTBPanel][image],
[image: ReaxFFPanel][image], or [image: ForceFieldPanel][image] panel
active, the **periodicity** will be automatically changed to **Slab**
(2D). For the other engines, 3D periodicity (**Bulk**) is retained,
and a vacuum gap is introduced. View the lattice parameters under
**Model → Lattice**.

Note:

  When the periodicity is set to **Slab**, the \vec{a} and \vec{b}
  surface lattice vectors must lie in the *xy* plane. View the lattice
  vectors under **Model → Lattice**.

Tip:

  If the surface unit cell is hexagonal (examples: Cu(111) and
  ZnO(0001)), the angle \gamma between the \vec{a} and \vec{b} lattice
  vectors may become 60°. It is more common to use an angle of 120°.
  You can change the angle under **Model → Lattice**, or use the
  [image: CrystalTool][image] **→ To Conventional Cell** tool.


Surface terminations
====================

The above definition of (hkl) Miller indices was based on intercepts
in the conventional unit cell, but **any** plane that is parallel to
that plane is *also* a (hkl) surface, but with a potentially different
surface termination.

In general, you should carefully look at the generated surface
structure to see if it is reasonable.

**Example**. Import "ZnO_wurtzite.xyz" and generate 8 layers of a "(1,
0, -1, 0)" slab (you can also insert "(1, 0, 0)"; for hexagonal
systems, the third miller index is ignored when inputting four
indices).

Should the termination be like in (a) or (b) below? In (a), the
surface atoms become 2-coordinated, which is unfavorable. In (b), the
surface atoms become 3-coordinated, which is more stable.

**Starting with AMS2026 the slab builder automatically decides that
(b) is correct**.

Similarly, you can generate "good" and "bad" terminations of the
ZnO(0001) surface.

   [image]Side views of ZnO slabs. a) Bad ZnO(10\bar{1}0) surface
   termination (two-coordinated surface atoms). b) Good
   ZnO(10\bar{1}0) termination (three-coordinated surface atoms). c)
   Bad ZnO(0001) and ZnO(000\bar{1}) terminations (one-coordinated
   surface atoms). d) Good ZnO(0001) and ZnO(000\bar{1}) terminations
   (three-coordinated surface atoms). Note that the "good"
   terminations in (d) do not address surface polarity.

Tip:

  **Create a thicker slab than you need**, and manually remove atoms
  until the surface terminations are correct. **Double-check your
  stoichiometry** under *Properties → Properties*, or by **Select →
  Select All**.


What makes a good slab model
============================

In general when creating a slab, you would like to

* **Retain the stoichiometry**. Example: a ZnO slab should ideally
  contain the same number of Zn as O atoms (just like in the bulk).

* **Maximize the coordination numbers** of the surface atoms (see
  Surface terminations)

* **Have identical "top" and "bottom" sides**: This may not be
  possible for all crystals or surfaces (example: a slab exposing Zn-
  terminated ZnO(0001) at the "top" side will necessarily expose the
  O-terminated ZnO(000\bar{1}) surface at the bottom side, see Fig. %s
  d)

Sometimes it is very difficult to create ideal slabs, especially if
the crystal contains **composite ions like carbonate, nitrate,
sulfate, or phosphate**.


Surface passivation
===================

**Surface passivation** is often accomplished by adding **H atoms**
that are coordinated by the surface atoms. Such H atoms can passivate
**dangling bonds**.

**Example**: The diamond(100) surface is terminated by two-coordinated
sp^3 C atoms. This is very unstable. By passivating the surface with
H, the dangling bonds at the surface become passivated.

   **1.** **Builders → Crystal Structure → Cubic → Diamond**, click OK.
   **2.** [image: CrystalTool][image] **→ Generate Slab**.
   **3.** Click the **Convert To Conventional Cell** button.
   **4.** Miller indices "(1, 0, 0)"
   **5.** Click the **Generate Slab** button.
   **6.** **Bonds → Guess bonds**
   **7.** Either manually add 2 H atoms to each surface C atom, or use the **Atoms → Add Hydrogens** option.

Tip:

  Between steps 5 and 6 you can get the primitive surface unit cell
  with [image: CrystalTool][image] **→Convert To Primitive Cell**.

   [image]DOS for **unpassivated** diamond(100) calculated with a
   single-point GFN1-xTB calculation. The DOS is high at the Fermi
   level (-0.43 Ha), indicating a metallic solution.

   [image]DOS for **passivated** diamond(100) calculated with a
   single-point GFN1-xTB calculation. The Fermi level lies in the band
   gap (diamond is an insulator).

Note:

  The **Add Hydrogens** function does not add H atoms to the most
  energetically favored positions. Randomly move the H atoms a little,
  and perform a geometry optimization, to find more suitable positions
  for the H atoms.

Note:

  The **Add Hydrogens** function does not work for all chemical
  elements. For some elements, you need to add the H atoms manually.


Surface relaxation
==================

A geometry optimization of a surface will typically give the result
that the surface atoms relax away from the positions of a bulk-
truncated crystal. For example, the surface atoms may move closer to
the subsurface atoms.

Surface relaxation can also happen for atoms in the subsurface layers.

In AMS2024, the *Show Info: Axes and Planes* feature was introduced.

   SCM → New Input
   File → Import coordinates from file →  "ZnO_1010_optimized_znorg.xyz"
   Select the two top-most atoms Zn(23) and O(24)
   Properties → Show Info: Axes and Planes

This brings up a new window:

[image]

Here you can see information about the vector formed by the two atoms
and its projection onto various other axes and planes.

Example: the projection onto the z-axis is 0.40 Å. This means that the
Zn atom relaxed 0.40 Å further into the surface than the O atom.

Example: the angle between the vector and the ab-plane (the plane
spanned by the lattice vectors) is 12.35°

Tip:

  * You can change the selected atom while the *Axes and Planes*
    window is still open

  * The *Show Info: Axes and Planes* feature can also be used from
    AMSmovie to see how quantities change during optimization or
    dynamics

You can also do any custom analysis you like directly on the atomic
coordinates:

* in **AMSinput**, go to **Model → Coordinates** to see how the
  coordinates

* in **AMSmovie**, go to **View → Coordinates** to see how the
  coordinates of the atoms change during a geometry optimization,

* **SCM → KFbrowser**, select **File → Expert mode**, view the input
  coordinates (in bohr) under **InputMolecule → Coords**, and the
  final coordinates (in bohr) under **Molecule → Coords**.

* view the input and final coordinates in a spreadsheet summary
  (**Tools → Build spreadsheet** in **AMSjobs**).


Surface reconstruction
======================

A surface reconstruction is a significant structural change at the
surface compared to the bulk.

**Example 1**: The Si(100) surface forms a (2×1) reconstruction in
which the surface atoms form rows of dimers (not found in the bulk).
(2×2) and *c*(4×2) reconstructions are also possible. Download a
"PLAMS python script" to run the jobs shown below, and visualize the
DOS in AMSdos (**SCM → Dos**).

   [image]DOS for **non-reconstructed** Si(100) calculated with HSE06
   and BAND. The DOS at the Fermi level is high, signifying a metallic
   solution. This calculation required special SCF convergence tricks
   (an electronic temperature).

   [image]DOS for a **(2×1)-reconstructed** Si(100) surface calculated
   with HSE06 and BAND at the PBE-optimized geometry. The Fermi level
   is at the top of the valence band, indicating a semiconducting
   solution. This calculation did not require any special SCF
   convergence tricks.

**Example 2**: The polar Zn-terminated ZnO(0001) surface can form
quite large triangular pits and islands ↗, which help to stabilize the
surface because the stoichiometry of the pits and islands is different
from the bulk stoichiometry. See also: Polar and non-polar surfaces.


Polar and non-polar surfaces
============================

Nonpolar surfaces are the most common and are typically more stable
than polar surfaces. For ionic compounds, some surfaces are polar.

A surface is polar if the formal dipole moment perpendicular to the
slab grows with the slab thickness. Here, the *formal dipole moment*
is calculated based on the formal charges of the ions. The growing
dipole moment gives a diverging energy as the slab thickness
increases.

Nevertheless, polar surfaces can form in real life, as they are
stabilized in some manner.

In calculations, a polar surface slab exposes two different polar
surfaces, one on the "top" side and one on the "bottom" side.  For
example, the Zn-terminated ZnO(0001) at the top and the O-terminated
ZnO(000\bar{1}) at the bottom (Fig. %s d).


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 [image: DFTBPanel][image] panel, set **Task → Single Point** and **Model → SCC-DFTB**
   **4.** Set **KSpace** to **Good** (surface metallization requires dense k-space sampling)
   **5.** Set **Occupation** to **Fermi** and **Fermi temperature** to **10000** K
   **6.** Set **Parameter directory** to **DFTB.org/znorg-0-1**.
   **7.** Run the job (it may require a few hundred SCC cycles)
   **8.** View results in AMSdos: **SCM → Dos**
   **9.** To view the partial DOS from only a few atom, select the atoms in the AMSdos window
   **10.** View the DOS for the top layer, a middle layer, and the bottom layer
   **11.** Repeat steps for "vacancies_good.xyz" and "adsorbates_good.xyz", except you can use the default occupation method (**Auto**)

Fig. %s a shows the DFTB-calculated DOS for a **clean ZnO(0001) slab**
(terminated as in Fig. %s 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**, and it may help to use **Fermi occupation
smearing**.

Fig. %s 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**.

   [image]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.

Note:

  Ideal (not stabilized) polar surfaces can be difficult to converge
  in DFT and DFTB calculations. Always try to stabilize the surface
  via change in stoichiometry or 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.


Surface unit cells and supercells
=================================

Just like a crystal, a surface also has a primitive unit cell
(containing only one lattice point). In surface-science-style
calculations, it is common to use a surface supercell. By creating a
supercell, the coverage of adsorbates can be fine-tuned.


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 [image: CrystalTool][image] **→ To Primitive Cell**. To get the
conventional surface unit cell, select [image: CrystalTool][image] **→
To Conventional Cell**.

Note:

  For surfaces, the primitive and conventional cells only differ for
  centered rectangular lattices.

**Example: The Cu(111) surface**

   **1.** Import "Cu_ccp_primitive.yz"
   **2.** [image: CrystalTool][image] **→ To Conventional Cell**. Note: this generates the conventional cell of the bulk Cu crystal, which is needed for the slab generation.
   **3.** [image: CrystalTool][image] **→ Generate Slab**
   **4.** Type Miller indices "(1, 1, 1)"
   **5.** Number of layers: 5
   **6.** **Generate slab**.

Tip:

  On the **Model → Lattice** page, tick the **Adjust atoms when
  changing lattice** checkbox and rotate the surface lattice vectors
  so that \vec{a} is parallel to *x*.


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

   [image]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. %s. 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 [image: CrystalTool][image] **→ 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}

Tip:

  Hexagonal *bulk* crystals (e.g., wurtzite ZnO) can similarly be
  converted to an orthorhombic supercell. The matrix is then
  \begin{pmatrix}1 & 0 & 0 \\ 1 & 2 & 0 \\ 0 & 0 & 1\end{pmatrix}


Surface energy calculations
===========================

The **surface energy** E^{\text{surf}} gives the energy required to
form a surface, per unit surface area.

In a slab calculation with two identical surfaces, the slab energy
E^\text{slab} can be divided into a "bulk" contribution and a
"surface" contribution:

   E^{\text{slab}} = \frac{n^\text{slab}}{n^\text{bulk}}E^\text{bulk}
   + 2AE^\text{surf}

where n^\text{slab} is the number of formula units in the slab unit
cell, n^\text{bulk} is the number of formula units in the bulk unit
cell, *A* is the surface area, and E^{\text{surf}} is the surface
energy.  The factor 2 appears because the slab exposes two surfaces.

**Surface energy method 1 (recommended)**

* Plot E^\text{slab} vs. n^\text{slab} for progressively thicker
  slabs,

* Fit a straight line through the datapoints,

* Extract E^{\text{surf}} from the intercept with the y-axis (at 0
  "thickness", n^\text{slab} = 0).

In this way, E^{\text{bulk}} is not needed explicitly (but corresponds
to the slope of the plot), and the k-point sampling is identical for
each slab thickness. This means that the k-points in the surface
calculation do not need to be extremely accurately converged.

**Surface energy method 2**

The above equation can be rearranged to give

   E^{\text{surf}} = \frac{E^\text{slab} -
   \frac{n^\text{slab}}{n^\text{bulk}}E^\text{bulk}}{2A}

E^{\text{surf}} can be calculated using this equation for a thick
slab. This requires that E^{\text{bulk}} be calculated explicitly.
Both E^{\text{slab}} and E^{\text{bulk}} must be very well converged
with respect to k-point sampling, and E^{\text{bulk}} must be
calculated for a very carefully converged (optimized) lattice.

Note:

  It is easiest to calculate the surface energy for slabs constructed
  from the **optimized** bulk lattice parameters, so that there is no
  strain the surface-lateral direction.If your intention is to model a
  surface of a bulk crystal, then **do not optimize the lattice of the
  slab**. Optimize the lattice of the **bulk**, and use that to
  construct the slab.If your intention is to model a **thin film**,
  you can consider optimizing the lattice of the slab.


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 [image: DFTBPanel][image] 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.** [image: CrystalTool][image] **→ Generate Slab**
   **7.** **Miller indices** "(1, 0, -1, 0)" (in previous AMS version, ignore the third Miller index for hexagonal crystals and enter "1 0 0")
   **8.** **Number of layers**: 4 (in this case the number of "double"-layers)
   **9.** Press the **Generate Slab** button.
   **10.** **Bonds → Guess bonds** may help with visualization.
   **11.** Close the slab builder. 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.

[image]

Using **method 1**, fit a straight line through the datapoints, as in
Fig. %s.

   [image]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**

   E^\text{surf} = \frac{\text{intercept}}{2A} = \frac{0.106415\
   \text{Ha}}{2 \times 17.68\ \text{Å}^2} = 0.00301\ \text{Ha/Å}^2 =
   1.31\ \text{J/m}^2

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

Tip:

  See the Python scripting example for surface energies.

Tip:

  Calculate the surface energy of ZnO(11\bar{2}0) using method 1 and
  method 2 with SCC-DFTB znorg-0-1. Answer: E^{\text{surf}} = 1.39
  J/m^2.


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
==========


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.

   [image]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:

[image]

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

   [image]Example setup for moving an H_2O adsorbate to the bridge
   position on Cu(111).

Tip:

  In the  **Tune Geometry** window you can also

  * **Translate atoms, direction**: translate the *Atoms to move* a
    certain **distance** along a certain **direction**. Use the  **⬭**
    label to define the direction from selected atoms. Click the ◀ ▶
    buttons to move the atom.

  * **Translate atoms, target**: translate the *Atoms to move* towards
    a certain **target position**. Use the  **⬭** label to define the
    target position based on the currently selected atom. Click the ◀
    ▶ buttons to move the atom.

  * **Rotate atoms**: rotate the *Atoms to move* around a certain
    **axis** (vector) that goes through a **center** (position). For
    example, set *Angle* to 20°, *axis* to "0", "0", "1", and *center*
    to "0", "0", "0". to rotate the *Atoms to move* around the z axis
    in steps of 20°.

  If you place the mouse cursor over an input field, you will get more
  information in a help balloon.


Adsorption energy
-----------------

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

   ΔE^\text{ads} = E^\text{slab} + E^\text{mol} - E^\text{slab+mol}

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**:

   ΔE^\text{ads} = E^\text{slab+mol} - (E^\text{slab} + E^\text{mol})

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

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

   [image]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

See also:

  * Basis set superposition error with BAND

  * In-cell approach for adsorption energies with ReaxFF


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

[image]

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.

   [image]Adsorption energies calculated with and without counterpoise
   correction using BAND with the PBE functional. A more negative
   value implies more exothermic adsorption. See the Python scripting
   example for this workflow.

+----------------------------------------------------+----------------------------------------------------+
| Basis set                                          | Counterpoise correction for H2O on ZnO(10-10)      |
|                                                    | (kcal/mol)                                         |
|====================================================|====================================================|
| DZ                                                 | 3.69                                               |
+----------------------------------------------------+----------------------------------------------------+
| DZP                                                | 3.00                                               |
+----------------------------------------------------+----------------------------------------------------+
| TZP                                                | 1.31                                               |
+----------------------------------------------------+----------------------------------------------------+
| TZ2P                                               | 1.15                                               |
+----------------------------------------------------+----------------------------------------------------+
| QZ4P                                               | 0.14                                               |
+----------------------------------------------------+----------------------------------------------------+

See also:

  * Basis set superposition error tutorial Method 1,

  * BAND BSSE example.


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**:

   E^\text{ads} = E^\text{slab+mol widely separated} -
   E^\text{slab+mol}

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.
%s 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:

   E^\text{ads} = E^\text{slab+mol} - E^\text{slab+mol widely
   separated}

**Example:** H_2O 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 Å.

+----------------------------------------------------+----------------------------------------------------+
|                                                    | kcal/mol                                           |
|====================================================|====================================================|
| E^slab                                             | -10331.25                                          |
+----------------------------------------------------+----------------------------------------------------+
| E^mol                                              | -248.85                                            |
+----------------------------------------------------+----------------------------------------------------+
| E^slab + E^mol                                     | -10580.11                                          |
+----------------------------------------------------+----------------------------------------------------+
| E^slab+mol widely separated                        | -10590.11                                          |
+----------------------------------------------------+----------------------------------------------------+
| E^slab+mol                                         | -10606.43                                          |
+----------------------------------------------------+----------------------------------------------------+
| E^ads (standard)                                   | -26.3                                              |
+----------------------------------------------------+----------------------------------------------------+
| E^ads (in-cell)                                    | -16.3                                              |
+----------------------------------------------------+----------------------------------------------------+

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

   E^\text{slab} + E^\text{mol} \ne E^\text{slab+mol widely separated}

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

   [image]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) See the Python scripting example for this workflow.

Tip:

  Keep the in-cell approach in mind when creating training sets for
  ReaxFF parametrization.

See also:

  Python scripting example for comparing in-cell and "standard"
  adsorption energies in ReaxFF.


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. %s and Fig. %s), 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
[image: AddButton][image] plus button next to *selected atoms (fix
positions)*.

**Example**: See the VASP: TiO₂ surface relaxation tutorial.


Solid-solid and solid-liquid interfaces
=======================================


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

Note:

  The builder will create a solid/liquid interface under 3D periodic
  boundary conditions.

For more details, see the Solid-liquid interface tutorial.

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