We start from the around 80000 inorganic compounds in the 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:
- Standard conventional unit cell was used for all the compounds.
- 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.
- 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.
- Make a 3×3×3 supercell and classify all the atoms in the supercell into clusters again. If the number of clusters in the supercell was three times the ones in the unit cell, the structure was tagged as layered.
- 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.
- 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 steps. A 20 Å vacuum along the 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.
The bottom-up 2D materials are generated by chemical substitution of the unary and binary top-down compounds by the following process.
- All the elements of the periodic table (from H to Bi) are categorized according to their group number. Radioactive elements, lanthanide and actinide with f electrons are excluded except La which is included in group 3 elements.
- Systematically replace each element in a known 2D material by one other element in the same category. For instance, 24 new materials can be generated from BN by replacing B with [B, Al, Ga, In, Tl] and N with [N, P, As, Sb, Bi] systematically.
- Structure matcher tools from Pymatgen were used to remove duplicates.
The standard workflow developed by the 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.
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 1.0x10-4
. 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.
The exfoliation energy of 2D materials is calculated by Eexf
, where E2D
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 undergo 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.
The following three cases are considered to the find the layered bulk counterpart of the 2D materials. If a 2D material has only one corresponding layered bulk in Materials Project, this bulk material is tagged as its layered bulk counterpart. If there are multiple layered bulks in MP corresponding to one 2D material, only the one with the lowest energy is tagged as its layered bulk counterpart. For 2D materials obtained through the bottom-up approach, there is often no corresponding layered bulk material in existence. In such a case, we also construct a “layered bulk” counterpart by stacking the monolayers in the same way the initial structure is stacked in its bulk counterpart. It is noted that for a few special cases like Silicene, the “exfoliation energy” is provided only for reference value. The energetic stability of a 2D material is estimated by the decomposition energy as discussed below.
The decomposition energy is defined as the energy required to separate a compound into its components. Here, we applied a modified version of the energy_above_hull ((https://hackingmaterials.com/2014/05/08/pdcomic/
), and define the decomposition energy as the energy required (per atom) to decompose a given material into the set of most stable materials at this chemical composition. This definition of decomposition energy excludes the corresponding layered bulk of the given 2D material from its competing phases. The energy difference between a 2D material and its layered bulk is given by the exfoliation energy.
The competing entries are queried from Materials Project, which are calculated without the Van der Waals correction. In order to make direct comparison, we also calculate the total energies of 2D materials without the Van der Waals correction. This energy and the energies of all entries belonging to the same chemical system in Materials Project are used to generate a phase diagram. The decomposition energy is obtained from the difference between the energy of the 2D material and those of its competing phases.