Calculations Guide

Layered-material screening

We are starting from the around 60000 inorganic compounds in Materials Project database. A geometry-based algorithm [Phys. Rev. Lett. 118, 106101 (2017)] was used to identify layered structures among these compounds by the following steps:
  1. Standard conventional unit cell was uses for all the compounds.
  2. Check whether two atoms in the unit cell are bonded or not. Use a threshold as the sum of the covalent radii [Dalton Trans., 2008, 2832-2838 DOI: 10.1039/B801115J] of two elements with a small tolerance. If the distance of two atoms smaller than the threshold: boned, else: not bonded.
  3. Classify all the atoms in the conventional unit cell into different clusters by checking if they are bonded. Count the number of clusters in the unit cell.
  4. Make a 2×2×2 supercell and classify all the atoms in the supercell into clusters again. If the number of clusters in the supercell was twice of the ones in the unit cell, the structure was tagged as layered.
  5. A set of tolerances from 0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, to 0.4 was used and only the ones identified as layered by at least two tolerances are left.
  6. We exclude the compounds with more than 4 elements or 40 atoms in the unit cell.

Theoretical exfoliation for 2D top-down materials

Two-dimensional (2D) materials were theoretically exfoliated by extracting one cluster in the standard conventional unit cell of the layered structures screened in the above step. A 20 Å vacuum along c axis was imposed to minimize the interactions of image slabs by periodic condition. Structure matcher tools from Pymatgen were used to find duplicates of the exfoliated 2D materials. Among the duplicates, only the ones that their layered bulk counterparts have the lowest energy_above_hull are used. There are totally 2274 2D materials generated by this process.

High-throughput calculations

The standard workflow developed by Materials Project was used to perform high-throughput calculations for all the layered bulk and 2D materials screened in this project. The calculations were performed by density functional theory as implemented in the Vienna Ab Initio Simulation Package (VASP) software with Perdew-Burke-Ernzerhof (PBE) approximation for the exchange-correlation functional and the frozen-core all-electron projector-augmented wave (PAW) method for the electron-ion interaction. The cutoff energy for the plane wave expansion was set to 520 eV.

Structure optimization

For structure optimization of 2D materials, both cell shape and volume was allowed to change (ISIF = 3) while the c axis was kept fixed. The energy difference for ionic convergence is set to 5.0e-05 per atom in the cell. The dispersion-corrected vdW-optB88 exchange-correlation functional was applied for structure relaxation.

DOS and band calculations

A static run with uniform (gamma centered) k-point grid for the relaxed structure is performed first to generate the charge density. Non-self-consistence runs based on the charge density with an energy range from EF (Fermi level)-10 eV to EF+10 eV with 2000 intervals were used for DOS simulation. The band structure computation uses line-mode k-point grid along the high symmetry points [Setyawan, W.; Curtarolo, S. Computational Materials Science 2010, 49, 299-312.] of the 2D Brillouin zone. The band structure is displayed with full lines for up spin and dashed lines for down spin. A colour map was used to indicate the contribution of each element. The total and projected DOS are shown with line in different colours. Please note that the DOS data and band structure can be slightly different as the kpoint grid is not the same. It is noted that Van der Waals corrections are not applied for these steps.

Exfoliation Energy

The exfoliation energy of 2D materials is calculated by Eexf = E2D-Ebulk, where E2D and Ebulk are the total energy per atom of the 2D and its corresponding bulk layered materials, respectively. To get the total energy, both the 2D and bulk materials undergoes a structure optimization and static calculation. The total energy of the static run is used to calculate the exfoliation energy. The dispersion-corrected vdW-optB88 exchange-correlation functional was applied to all the runs for exfoliation energy.


Energy_above_hull is defined by the energy of decomposition of this material into the set of most stable materials at this chemical composition. [Chem. Mater. 20, 1798 (2008)]. The following steps were applied to compute the energy_above_hull (in meV/atom) of a 2D material: The energy_above_hull is the energy difference between the 2D material and the convex hull formed by all stable entries at each composition. (

Decomposition stability

In order to estimate the decomposition stability of a 2D material, we also applied a modified version of energy_above_hull, which excludes the corresponding layered bulk form of this 2D material from its competing phases. The energy difference between a 2D material and its layered bulk can be found in the exfoliation energy.

Bottom-up 2D materials

Among the 2274 2D materials from top-down approach, there are 560 binary compounds. The bottom-up 2D materials are generated by chemical substitution of these binary compounds by the following process.
  1. All the elements in same column of the periodic table (from H to Bi) are categorized in one group. Then, radioactive elements, lanthanide and actinide are removed except that La is added to group 3. H and O are also excluded due to their significant difference with the other elements in the group.
  2. Replace the elements in the top-down 2D materials with all the other elements in the same group. For instance, replace B with [B, Al, Ga, In, Tl] and N with [N, P, As, Sb, Bi] of BN and in principle 24 new materials can be generated. By such a process, 5407 new binary compounds are obtained.
  3. Remove duplicates of these 5407 within bottom-up compounds and with top-down 2D materials. 3196 unique new 2D binary compounds are generated.
  4. Standard high-throughput calculations with structure optimization->static->DOS, band structure runs are performed the 3196 compounds. Vdw-optB88 is applied for the structure optimization and static steps but not for DOS and band structure runs.