The Computational 2D Materials Database: High-Throughput Modeling and Discovery of Atomically Thin Crystals

We introduce the Computational 2D Materials Database (C2DB), which organises a variety of structural, thermodynamic, elastic, electronic, magnetic, and optical properties of around 1500 two-dimensional materials distributed over more than 30 different crystal structures. Material properties are systematically calculated by state-of-the art density functional theory and many-body perturbation theory (G$_0\!$W$\!_0$ and the Bethe-Salpeter Equation for $\sim$200 materials) following a semi-automated workflow for maximal consistency and transparency. The C2DB is fully open and can be browsed online or downloaded in its entirety. In this paper, we describe the workflow behind the database, present an overview of the properties and materials currently available, and explore trends and correlations in the data. Moreover, we identify a large number of new potentially synthesisable 2D materials with interesting properties targeting applications within spintronics, (opto-)electronics, and plasmonics. The C2DB offers a comprehensive and easily accessible overview of the rapidly expanding family of 2D materials and forms an ideal platform for computational modeling and design of new 2D materials and van der Waals heterostructures.


Introduction
Over the past decade, atomically thin two-dimensional (2D) materials have made their way to the forefront of several research areas including batteries, (electro-) catalysis, electronics, and photonics [1,2]. This development was prompted by the intriguing and easily tunable properties of atomically thin crystals and has been fueled by the constant discovery of new 2D materials and the emergent concepts of lateral [3] and vertical [4] 2D heterostructures, which opens completely new possibilities for designing materials with tailored and superior properties.
In the search for new materials with tailored properties or novel functionalities, first-principles calcul ations are playing an increasingly important role. The continuous increase in computing power and significant advancements of theoretical methods and numerical algorithms have pushed the field to a point where first-principles calculations are comparable to experiments in terms of accuracy and greatly The Computational 2D Materials Database: high-throughput modeling and discovery of atomically thin crystals 2 S Haastrup et al surpass them in terms of speed and cost. For more than a century, experimental databases on e.g. structural, thermal, and electronic properties, have been a cornerstone of materials science, and in the past decade, the experimental data have been augmented by an explosion of computational data obtained from firstprinciples calculations. Strong efforts are currently being focused on storing and organising the computational data in open repositories [20,21]. Some of the larger repositories, together containing millions of material entries, are the Materials Project [22], the Automatic Flow for Materials Discovery (AFLOW-LIB) [23], the Open Quantum Materials Database (OQMD) [24,25], and the Novel Materials Discovery (NOMAD) Repository [26].
The advantages of computational materials databases are many. Most obviously, they facilitate open sharing and comparison of research data whilst reducing duplication of efforts. In addition, they underpin the development and benchmarking of new methods by providing easy access to common reference systems [27]. Finally, the databases enable the application of machine learning techniques to identify deep and complex correlations in the materials space and to use them for designing materials with tailored properties and for accelerating the discovery of new mat erials [28][29][30]. Among the challenges facing the computational databases is the quality of the stored data, which depends both on the numerical precision (e.g. the employed k-point grid and basis set size) and the accuracy of the employed physical models (e.g. the exchange-correlation functional). Most of the existing computational databases store results of standard density functional theory (DFT) calculations. While such methods, when properly conducted, are quite reliable for ground state properties such as structural and thermodynamic properties, they are generally not quantitatively accurate for excited state properties such as electronic band structures and optical absorption spectra.
Compared to databases of bulk materials, databases of 2D materials are still few and less developed. Early work used DFT to explore the stability and electronic structures of monolayers of group III-V honeycomb lattices [31,32] and the class of MX 2 trans ition metal dichalcogenides and oxides [33]. Later, by datafiltering the inorganic crystal structure database (ICSD), 92 experimentally known layered crystals were identified and their electronic band structures calculated at the DFT level [34]. Another DFT study, also focused on stability and band structures, explored around one hundred 2D materials selected from different structure classes [35]. To overcome the known limitations of DFT, a database with many-body G 0 W 0 band structures for 50 semiconducting TMDCs was established [36]. Very recently, data mining of the Materials Project and experimental crystal structure databases in the spirit of [34], led to the identification of close to one thousand experimentally known layered crystals from which single layers could potentially be exfoliated [16][17][18][19]. These works also computed basic energetic, structural and electronic properties of the monolayers (or at least selected subsets) at the DFT level.
In this paper, we introduce the open Computational 2D Materials Database (C2DB) which organises a variety of ab initio calculated properties for more than 1500 different 2D materials. The key characteristics of the C2DB are: • Materials: the database focuses entirely on 2D materials, i.e. isolated monolayers, obtained by combinatorial lattice decoration of known crystal structure prototypes. • Consistency: all properties of all materials are calculated using the same code and parameter settings following the same workflow for maximum transparency, reproducibility, and consistency of the data. • Properties: the database contains a large and diverse set of properties covering structural, thermodynamic, magnetic, elastic, electronic, dielectric and optical properties. • Accuracy: Hybrid functionals (HSE06) as well as beyond-DFT many-body perturbation theory (G 0 W 0 ) are employed to obtain quantitatively accurate band structures, and optical properties are obtained from the random phase approximation (RPA) and Bethe-Salpeter equation (BSE). • Openness: the database is freely accessible and can be directly downloaded and browsed online using simple and advanced queries.
The systematic combinatorial approach used to generate the structures in the database inevitably produces many materials that are unstable and thus unrealistic and impossible to synthesise in reality. Such 'hypothetical' structures may, however, still be useful in a number of contexts, e.g. for method development and benchmarking, testing and training of machine learning algorithms, identification of trends and structureproperty relationships, etc. For this reason we map out the properties of all but the most unstable (and thus chemically unreasonable) compounds. Nevertheless, the reliable assessment of stability and synthesisability of the candidate structures is an essential issue. Using the 55 materials in the C2DB, which have been experimentally synthesised in monolayer form, as a guideline, we set down the criteria that a hypothesised 2D material should fulfill in order for it to be 'likely synthesisable'. On the basis of these criteria, we introduce a simple stability scale to quantify a candidate material's dynamic and thermodynamic stability. Out of an initial set of around 1900 monolayers distributed over 32 different crystal structures, we find 350 in the most stable category. In addition to the 55 experimentally synthesised monolayers, this set also includes around 80 mono layers from experimentally known vdW layered bulk materials, and thus around 200 completely new and potentially synthesisable 2D materials.
In section 2, we describe the computational workflow behind the database. The structure and properties of the materials are calculated using well established state-of-the-art methodology. Technical descriptions of the different steps in the workflow are accompanied by illustrative examples and comparisons with literature data. Since documentation and validation is the main purpose of the section, we deliberately focus on well known 2D materials like the Xenes and transition metal dichalcogenides where plenty of both computational and experimental reference data is available. It should be clear that the novelty of the present work does not lie in the employed methodology nor in the type of materials properties that we calculate (we note, however, that to the best of our knowledge the present compilation of GW and BSE calculations represents the largest of its kind reported so far). The significance of our work is rather reflected by the fact that when large and consistently produced data sets are organised and made easily accessible, new scientific opportunities arise. As outlined below, this paper presents several examples of this effect.
In section 3 we give an overview of the materials and the data contained in the C2DB and provide some specific examples to illustrate its use. Using an extensive set of many-body G 0 W 0 calculations as a reference, we establish the performance of various DFT xcfunctionals for predicting band gaps, band edge positions, and band alignment at hetero-interfaces, and we propose an optimal strategy for obtaining accurate band energies at low computational cost. Similarly, the 250 BSE calculations allow us to explore trends in exciton binding energies and perform a statistically significant and unbiased assessment of the accuracy and limitations of the widely used Mott-Wannier model for 2D excitons. From the data on more than 600 semiconductor monolayers, we present strong empirical evidence against an often employed relation between effective masses and band gaps derived from k · p perturbation theory. Inspired by the potential of using 2D materials as building blocks for plasmonics and photonics, we propose a model to predict the plasmon dispersion relations in 2D metals from the (intraband) plasma frequency and the onset of interband trans itions and use it to identify 2D metals with plasmons in the optical frequency regime. We propose several new magnetic 2D materials (including both metals and semiconductors) with ferromagnetic or anti-ferromagnetic ordering and significant out-ofplane magnetic anisotropy. Finally, we point to new high-mobility 2D semiconductors including some with band gaps in the range of interest for (opto)electronic applications.
In section 4 we provide our conclusions together with an outlook discussing some opportunities and possible future directions for the C2DB.

Workflow
The workflow used to generate the data in the C2DB is illustrated in figure 1. It consists of two parts: In the first part (left panel) the unit cell and atom positions are optimised for different magnetic configurations: non-magnetic (NM), ferro-magnetic (FM) and antiferro-magnetic (AFM). Materials satisfying certain stability and geometry criteria (indicated by green boxes) are subject to the second part (right panel) where the different properties are computed using DFT and many-body methods. The G 0 W 0 band structure and BSE absorbance calculations have been performed only for semiconducting materials with up to four atoms in the unit cell. Per default, properties shown in the online database include spin-orbit coupling (SOC); however, to aid comparison with other calculations, most properties are also calculated and stored without SOC.
All DFT and many-body calculations are performed with the projector augmented wave code GPAW [37] using plane wave basis sets and PAW potentials version 0.9.2. The workflow is managed using the Python based atomic simulation environment (ASE) [38]. We have developed a library of robust and numerically accurate (convergence verified) ASE-GPAW scripts to perform the various tasks of the workflow, and to create the database afterwards. The library is freely available, under a GPL license.
Below we describe all steps of the workflow in detail. As the main purpose is to document the workflow, the focus is on technical aspects, including numerical convergence and benchmarking. An overview of the most important parameters used for the different calculations is provided in table 1.

Structure relaxation
The workflow is initiated with a crystal structure defined by its unit cell (Bravais lattice and atomic basis). The crystal lattice is typically that of an experimentally known prototype (the 'seed structure') decorated with atoms picked from a subset of the periodic table, see figure 2. We refer to materials by the chemical formula of their unit cell followed by the crystal structure. The latter is indicated by a representative material of that prototype, as described in section 3.1. For example, monolayer MoS 2 in the hexagonal H and T phases are denoted MoS 2 -MoS 2 and MoS 2 -CdI 2 , respectively. Now, MoS 2 is in fact not stable in the T phase, but undergoes a 2 × 1 distortion to the so-called T′ phase. Because the T′ phase is the thermodynamically stable phase of WTe 2 , we denote MoS 2 in the distorted T phase by Mo 2 S 4 -WTe 2 . In the following, we shall refer to the unit cell with which the workflow is initiated, i.e. the unit cell of the seed structure, as the primitive cell or the 1 × 1 cell, even if this cell is not dynamically stable for the considered material (see section 2.4).
The unit cell and internal coordinates of the atoms are relaxed in both a spin-paired (NM), ferromagnetic (FM), and anti-ferromagnetic (AFM) configuration. Calculations for the AFM configuration are performed only for unit cells containing at least two metal atoms.
The symmetries of the initial seed structure are kept during relaxation. All relevant computational details are provided in table 1.
After relaxation, we check that the structure has remained a covalently connected 2D material and not disintegrated into 1D or 0D clusters. This is done by defining clusters of atoms using the covalent radius [39] + 30% as a measure for covalent bonds between atoms. The dimensionality of a cluster is obtained from the scaling of the number of atoms in a cluster upon repetition of the unit cell following the method described by Ashton et al [16]. Only materials containing exactly one cluster of dimensionality 2 are given further consideration (an exception is made for the metal-organic perovskites (prototype PbA 2 I 4 ) for which the metal atom inside the octahedron represents a 0D cluster embedded in a 2D cluster). To illustrate the effect of the covalent radius + 30% threshold, figure 3 shows the distribution of the candidate structures in the database as a function of the covalent factor needed to fully connect the structure. Most materials have a critical covalent factor below 1.3 and fall in the green shaded region. There is, however, a tail of around 100 disconnected materials (red region); these materials are not included in the database (see first green box in figure 1).
We also check that the material is not already contained in the database (second green box in figure 1). This is done by measuring the root mean square distance (RMSD) [40] relative to all other materials in the C2DB with the same reduced chemical formula. A threshold of 0.01 Å is used for this test.
In case of multiple metastable magnetic configurations (in practice, if both a FM and AFM ground state are found), these are regarded as different phases of the same material and will be treated separately throughout the rest of the workflow. To indicate the magnetic phase we add the extensions 'FM' or 'AFM' to the material name. The total energy of the spinpaired ground state is always stored, even when it is not the lowest. If the energy of the non-magnetic state is higher than the most stable magnetic state by less than 10 meV/atom, the workflow is also performed for the non-magnetic state. This is done in recognition of the finite accuracy of DFT for predicting the correct energetic ordering of different magnetic states.
We have compared the lattice constants of 29 monolayers with those reported in [41], which were obtained with the VASP code using PBE and very similar numerical settings and find a mean absolute deviation of 0.024 Å corresponding to 0.4%. The small yet finite deviations are ascribed to differences in the employed PAW potentials.

Crystal structure classification 2.2.1. Symmetry
To classify the symmetries of the crystal structure the 3D space group is determined using the crystal symmetry library Spglib [42] on the 3D supercell with a tolerance of 10 −4 Å.

Prototypes
The materials are classified into crystal structure prototypes based on the symmetry of the crystals. For two materials to belong to the same prototype, we require that they have the same space group, the same stoichiometry, and comparable thicknesses. The last requirement is included to distinguish between materials with the same symmetry and stoichiometry but with different number of atomic layers, see for example monolayer BN and GaS in figure 4. Each prototype is labelled by a specific representative material. For prototypes which have been previously investigated, we comply with the established conventions. However, since the field of 2D materials is still young and because C2DB contains a large number of never-synthesised materials, some of the considered crystal structures fall outside the known prototypes. In these cases we have chosen the representative material to be the one with the lowest energy with respect to the convex hull. Some of the crystal structure prototypes presently contained in the C2DB are shown in figure 4.

Thermodynamic stability
The heat of formation, ∆H , is defined as the energy of the material with respect to the standard states of its constituent elements. For example, the heat of formation per atom of a binary compound, A x B y , is Table 1. Overview of the methods and parameters used for the different steps of the workflow. If a parameter is not specified at a given step, its value equals that of the last step where it was specified. ; PW cutoff = 50 eV; truncated Coulomb interaction; at least 4 occupied and 4 empty bands a For the cases with convergence issues, we set a k-point density of 9.0 and a smearing of 0.02 eV.
where E(A x B y ) is the total energy of the material A x B y , and E(A) and E(B) are the total energies of the elements A and B in their standard state. When assessing the stability of a material in the C2DB, it should be kept in mind that the accuracy of the PBE functional for the heat of formation is only around 0.2 eV/atom on average [43]. Other materials databases, e.g. OQMD, Materials Project, and AFLOW, employ fitted elementary reference energies (FERE) [44] and apply a Hubbard U term [45] for the rare earth and transition metal atoms (or a selected subset of them). While such correction schemes in general improve ∆H they also introduce some ambiguity, e.g. the dataset from which the FERE are determined or the exact form of the orbitals on which the U term is applied. Thus in order not to compromise the transparency and reproducibility of the data we use the pure PBE energies.
For a material to be thermodynamically stable it is necessary but not sufficient that ∆H < 0. Indeed, thermodynamic stability requires that ∆H be negative not only relative to its pure elemental phases but relative to all other competing phases, i.e. its energy must be below the convex hull [46]. We stress, however, that in general, but for 2D materials in particular, this definition cannot be directly applied as a criterion for stability and synthesisability. The most important reasons for this are (i) the intrinsic uncertainty on the DFT energies stemming from the approximate xc-functional (ii) substrate interactions or other external effects that can stabilise the monolayer (iii) kinetic barriers that separate the monolayer from other lower energy phases rendering the monolayer (meta)stable for all practical purposes.
We calculate the energy of the 2D material relative to the convex hull of competing bulk phases, ∆H hull . The convex hull is currently constructed from the 2836 most stable binary bulk compounds which were obtained from the OQMD [24]. The energies of the bulk phases were recalculated with GPAW using the PBE xc-functional and the same numerical settings as applied for the 2D materials (but the structure was not re-optimised). Because the bulk reference structures from OQMD were optimised with the VASP code and with Hubbard U corrections for materials containing 3d elements, and because the PBE misses attractive vdW interaction, the bulk energies could be slightly overestimated relative to the monolayers. As a consequence, monolayers that also exist in a layered bulk phase could have ∆H hull < 0, even if the layered bulk phase is part of the convex hull and thus should be energetically more stable than the monolayer. Comparing our ∆H hull values for 35 compounds with the exfoliation energies calculated in [18] employing vdW compliant xc-functionals for both bulk and monolayer, we estimate the errors in the convex hull energies to be below 0.1 eV/atom.
As an example, the convex hull for Fe x Se 1−x is shown in figure 5. The convex hull as defined by the bulk binaries is indicated by the blue lines. The labels for the 2D materials refer to the crystal prototype and magnetic order. Clearly, most 2D materials lie above the convex hull and are thus predicted to be thermodynamically unstable in freestanding form under standard conditions. However, as mentioned above, depending on the material, errors on the PBE formation energies can be sizable and thus the hull diagram should only be taken as guideline. Nevertheless, in the present example we find that FeSe (which is itself a prototype) with anti-ferromagnetic ordering lies slightly below the convex hull and is thus predicted to be thermodynamically stable. This prediction is consistent with the recent experimental observation that monolayer FeSe deposited on SrTiO 3 exhibits AFM order [47].

Phonons and dynamic stability
Due to the applied symmetry constraints and/or the limited size of the unit cell, there is a risk that the structure obtained after relaxation does not represent a local minimum of the potential energy surface, but only a saddle point. We test for dynamical stability by calculating the Γ-point phonons of a 2 × 2 repeated cell (without re-optimising the structure) as well as the elastic constants (see section 2.5). These quantities represent second-order derivatives of the total energy with respect to atom displacements and unit cell lengths, respectively, and negative values for either quantity indicate a structural instability.
The Γ-point phonons of the 2 × 2 supercell are obtained using the finite displacement method [48]. We displace each atom in the primitive cell by ±0.01 Å, and calculate the forces induced on all the atoms in the supercell. From the forces we construct the dynamical matrix, which is diagonalised to obtain the Γ-point phonons of the 2 × 2 cell (or equivalently the Γ-point and zone boundary phonons of the primitive cell). The eigenvalues of the dynamical matrix correspond to the square of the mass-renormalised phonon frequencies, ω. Negative eigenvalues are equivalent to imaginary frequencies and signal a saddle point. Our procedure explicitly tests for stability against local distortions of periodicities up to 2 × 2 and thus provides a necessary, but not sufficient condition for dynamic stability. We stress, however, that even in cases where a material would spontaneously relax into a structure with periodicity larger than 2 × 2, the Γ-point dynamical matrix of the 2 × 2 cell could exhibit negative eigenvalues. Our test is thus more stringent than it might seem at first glance. In principle, a rigorous test for dynamic stability would require the calculation of the full phonon band structure. Mathematically, the instabilities missed by our approach are those that result in imaginary phonons in the interior of the BZ but not at the zone boundary. Physically, such modes could be out of plane buckling or charge density wave-driven reconstructions with periodicities of several unit cells. In general, however, these types of instabilities are typically rather weak (as measured by the magnitude of the imaginary frequency) as compared to more local distortions such as the T to T′ distortion considered below. Moreover, they could well be a special property of the isolated monolayer and become stabilised by the ubiquitous interactions of the 2D material with its environment, e.g. substrates. This is in fact supported by the full phonon calculations by Mounet et al for ∼ 250 isolated monolayers predicted to be easily exfoliable from experimentally known layered bulk phases [18]. Indeed, most of the instabilities revealed by their calculations are of the type described above and would thus be missed by our test. However, these instabilities cannot be too critical as the monolayers are known to be stable in the vdW bonded layered bulk structure.
As an example, figure 6 compares the dynamical stability of a subset of transition metal dichalcogenides and -oxides in the T and T′ phases (CdI 2 and WTe 2 prototypes). The two upper panels show the smallest eigenvalue of the Γ-point dynamical matrix of the 2 × 2 cell. Only materials above the dashed line are considered dynamically stable (for this example we do not consider the sign of the elastic constants which could further reduce the set of dynamically stable materials). Since the unit cell of the T′ phase contains that of the T phase it is likely that a material initially set up in the T′ phase relaxes back to the T phase. To identify these cases, and thereby avoid the presence of duplicates in the database, the third panel shows the root mean square distance (RMSD) between the structures obtained after relaxations starting in the T-and T′ phase, respectively. Structures below the dashed line are considered identical. The color of each symbol refers to the four different potential energy surfaces illustrated at the bottom of the figure.

Stability criteria
To assess the stability of the materials in the C2DB, we turn to the set of experimentally synthesised/exfoliated monolayers. For these materials, the calculated energy above the convex hull and minimum eigenvalue of the dynamical matrix are shown on figure 7. It is clear that all but five known monolayers have a hull energy below 0.2 eV/atom, and three of these have only been synthesised on a metal substrate. Turning to the dynamical stability, all but one of the experimentally known monolayers have a minimum eigenvalue of the dynamical matrix above −2 eVÅ −2 , and 70% have a minimum eigenvalue above −1 × 10 −5 eVÅ −2 .
Guided by these considerations, we assign each material in the C2DB a stability level (low, medium or high) for both dynamical and thermodynamic stability, as illustrated in table 2. For ease of reference, we also define the overall stability level of a given mat erial as the lower of the dynamical and thermodynamic stability levels. If a material has 'low' overall stability (marked by bold in the table), we consider it unstable and do not carry out the rest of the workflow. Mat erials with 'high' overall stability are considered likely to be stable and thus potentially synthesisable. Mat erials in the 'medium' stability category, while unlikely to be stable as freestanding monolayers, cannot be discarded and might be metastable and possible to synthesise under the right conditions. For example, freestanding silicene has a heat of formation of 0.66 eV/ atom, but can be grown on a silver substrate. Likewise, the T′ phase of MoS 2 (WTe 2 prototype) has an energy of 0.27 eV/atom higher than the thermodynamically stable H phase, but can be stabilised by electron doping.

Elastic constants
The elastic constants of a material are defined by the generalised Hooke's law, where σ ij , C ijkl and ε kl are the stress, stiffness and strain tensors, respectively, and where we have used the Einstein summation convention. In two dimensions, the stress and strain tensors have three independent components, namely planar stress/ strain in the x and y directions, as well as shear stress/ strain. The stiffness tensor is a symmetric linear map between these two tensors, and therefore has up to six independent components. Disregarding shear deformations, the relationship between planar strain and stress is For all materials in the C2DB, we calculate the planar elastic stiffness coefficients C 11 , C 22 , and C 12 . These are calculated using a central difference approximation to the derivative of the stress tensor: the material is strained along one of the coordinate axes, x or y, and the stress tensor is calculated after the ions have relaxed. We use strains of ±1% which we have found to be sufficiently large to eliminate effects of numerical The convex hull as defined by the bulk phases is represented by the blue lines. Blue squares denote bulk binary reference phases while orange triangles represent 2D materials. The labels for the 2D materials refer to the crystal prototype and magnetic order.
noise and sufficiently small to stay within the linear response regime. Table 3 shows the calculated planar stiffness coefficients of a set of 2D materials. As can be seen the values from the C2DB are in very good agreement with previously published PBE results. For the isotropic mat erials MoS 2 , WSe 2 and WS 2 , C 11 and C 22 should be identical, and we see a variation of up to 0.6%. This provides a test of how well converged the values are with respect to numerical settings.

Magnetic anisotropy
The energy dependence on the direction of magnetisation, or magnetic anisotropy (MA), arises from spin-orbit coupling (SOC). According to the magnetic force theorem [96] this can be evaluated from the eigenvalue differences such that the correction to the energy becomes where εn kn and f (εn kn ) are the eigenenergies and occupation numbers, respectively, obtained by diagonalising the Kohn-Sham Hamiltonian including SOC in a basis of collinear spinors aligned along the direction n, while ε 0 kn and f (ε 0 kn ) are the bare Kohn-Sham eigenenergies and occupation numbers without SOC.
For all magnetic materials we have calculated the energy difference between out-of-plane and inplane magnetisation E MA (i) = ∆E(ẑ) − ∆E(i), (i =x,ŷ). Negative values of E MA (i) thus indicate that there is an out-of-plane easy axis of magnet isation.
Calculations for the ground state have been performed with plane-wave cutoff and energetic convergence threshold set to 800 eV and 0.5 meV/atom respectively. For all calculations we have used a Γ-centered Monkhorst-Pack k-point with a density of 20/Å −1 . The SOC contribution is introduced via a non-self-consistent diagonalisation of the Kohn-Sham Hamiltonian evaluated in the projectoraugmented wave formalism [97].

Projected density of states
The projected density of states (PDOS) is a useful tool for identifying which atomic orbitals comprise a band. It is defined as where ψ kn are the Kohn-Sham wave functions with eigenvalues ε kn and φ a l,m are the spin-paired Kohn-Sham orbitals of atomic species S with angular momentum l (s, p, d, f ). We sum over all atoms belonging to species S so every atomic species has one T T T T T T T T Figure 6. Dynamical stability of a set of transition metal dichalcogenides and -oxides in the T and Tʹ phases (CdI 2 and WTe 2 prototypes), respectively. The first and second panels show the minimum eigenvalue of the Γ-point dynamical matrix of the 2 × 2 unit cell (containing 12 and 24 atoms for the T and T′ phase, respectively. The lower panel shows the root mean square distance (RMSD) between the relaxed structures. The color indicates whether the material is dynamically stable in the T phase (black), the T′ phase (blue), both phases (orange) or neither of the phases (green). entry per angular momentum channel. In the PAW formalism this can be approximated as where ψ kn are the pseudo wave functions and p a l,m are the PAW projectors associated with the atomic orbitals φ a l,m . The PDOS is calculated from equation (6) using linear tetrahedron interpolation [98] (LTI) of energy eigenvalues obtained from a ground state calculation with a k-point sampling of 36 . In contrast to other techniques for calculating the PDOS using smearing, the PDOS yielded by the LTI method returns exactly zero at energies with no states. Examples of PDOS Figure 7. The calculated energy above the convex hull and minimum eigenvalue of the dynamical matrix (evaluated at the Γ-point for the 2 × 2 cell) for the 55 materials in the C2DB that have been synthesised or exfoliated in monolayer form, see [6,9,10,12,. The three materials highlighted in red have only been synthesised on metallic substrates. The black dashed lines indicate the thresholds used to categorise the thermodynamic and dynamic stability of materials in the C2DB. are shown in figure 9 (right) for respectively the ferromagnetic metal VO 2 and the semiconductor WS 2 in the H phase (MoS 2 prototype).

Band structures
Electronic band structures are calculated along the high symmetry paths shown in figure 8 for the five different types of 2D Bravais lattices. The band energies are computed within DFT using three different xc-functionals, namely PBE, HSE06, and GLLBSC. These single-particle approaches are complemented by many-body G 0 W 0 calculations for materials with a finite gap and up to four atoms in the unit cell (currently around 250 materials). For all methods, SOC is included by non-selfconsistent diagonalisation in the full basis of Kohn-Sham eigenstates. Band energies always refer to the vacuum level defined as the asymptotic limit of the Hartree potential, see figure 12. Below we outline the employed methodology while section 3.2.1 provides an overview and comparison of the band energies obtained with the different methods.

PBE band structure
The electron density is determined self-consistently on a uniform k-point grid of density 12.0/Å −1 . From this density, the PBE band structure is computed non-selfconsistently at 400 k-points distributed along the band path (see figure 8). Examples of PBE band structures are shown in figure 9 for the ferromagnetic metal VO 2 and the semiconductor WS 2 both in the MoS 2 prototype structure. The expectation value of the out-of-plane spin component, χ nkσ |Ŝ z |χ nkσ , is evaluated for each spinorial wave function, χ nkσ = (ψ nk↑ , ψ nk↓ ), and is indicated by the color of the band. For materials with inversion symmetry, the SOC cannot induce band splitting, meaning that χ nkσ |Ŝ z |χ nkσ is ill-defined and no color coding is used. The band structure without SOC is indicated by a dashed grey line. We have compared our PBE + SOC band gaps of 29 different monolayers with those obtained with the VASP code in [41] and find a mean absolute deviation of 0.041 eV.

HSE band structure
The band structure is calculated non-selfconsistently using the range-separated hybrid functional HSE06 [99] on top of a PBE calculation with k-point density 12.0/Å −1 and 800 eV plane wave cutoff. We have checked for selected systems that the HSE band structure is well converged with these settings. The energies along the band path are obtained by spline interpolation from the uniform k-point grid. As an example, the HSE band structure of WS 2 is shown in the left panel of figure 10 (black line) together with the PBE band structure (grey dashed). The PBE band gap increases from 1.52 eV to 2.05 eV with the HSE06 functional in good agreement with earlier work reporting band gaps of 1.50 eV (PBE) and 1.90 eV (HSE) [100] and 1.55 eV (PBE) and 1.98 eV (HSE) [101], respectively. A more systematic comparison of our results with the HSE + SOC band gaps obtained with the VASP code in [41] for 29 monolayers yield a mean absolute deviation of 0.14 eV. We suspect this small but non-zero deviation is due to differences in the employed PAW potentials and the non-selfconsistent treatment of the HSE in our calculations. Table 2. The materials in the C2DB distributed over the nine stability categories defined by the three levels (high, medium and low) of dynamical stability (columns) and thermodynamic stability (rows). The overall stability of the materials is defined as the lower of the two separate stability scales. Materials with low overall stability (bold) are considered unstable.
Thermodynamic stability (eV/atom) Dynamic stability (eVÅ −2 )  Table 3. Planar elastic stiffness coefficients (in N m −1 ) calculated at the PBE level. The results of this work are compared to previous calculations from the literature and the mean absolute deviation (MAD) is shown.

GLLBSC fundamental gap
For materials with a finite PBE band gap, the fundamental gap (i.e. the difference between the ionisation potential and electron affinity) also sometimes referred to as the quasiparticle gap, is calculated self-consistently using the GLLBSC [102] xc-functional with a Monkhorst-Pack k-point grid of density 12.0/Å −1 . The GLLBSC is an orbitaldependent exact exchange-based functional, which evaluates the fundamental gap as the sum of the Kohn-Sham gap and the xc-derivative discontinuity, E gap = ε KS gap + ∆ xc . The method has been shown to yield excellent quasiparticle band gaps at very low computational cost for both bulk [102,103] and 2D semiconductors [36].
In the exact Kohn-Sham theory, ε KS v should equal the exact ionisation potential and thus ∆ xc should be used to correct only the conduction band energies [104]. Unfortunately, we have found that in practice this procedure leads to up-shifted band energies (compared with the presumably more accurate G 0 W 0 results, see figure 20). Consequently, we store only the fundamental gap and ∆ xc in the database. However, as will be shown in section 3.2.1 the center of the gap is in fact reasonably well described by PBE suggesting that efficient and fairly accurate predictions of the absolute band edge energies can be obtained by a symmetric GLLBSC correction of the PBE band edges.

G 0 W 0 band structure
For materials with finite PBE band gap the quasiparticle (QP) band structure is calculated using the G 0 W 0 approximation on top of PBE following our earlier work [105,106]. Currently, this resource demanding step is performed only for materials with up to four atoms in the unit cell. The number of plane waves and the number of unoccupied bands included in the calculation of the non-interacting density response function and the GW self-energy are always set equal. The individual QP energies are extrapolated to the infinite basis set limit from calculations at plane wave cutoffs of 170, 185 and 200 eV, following the standard 1/N G dependence [107,108], see figure 11 (right). The screened Coulomb interaction is represented on a non-linear real frequency grid ranging from 0 eV to 230 eV and includes around 250 frequency points. The exchange contribution to the self-energy is calculated using a Wigner-Seitz truncation scheme [109] for a more efficient treatment of the long range part of the exchange potential. For the correlation part of the selfenergy, a 2D truncation of the Coulomb interaction is used [110,111]. We stress that the use of a truncated Coulomb interaction is essential to avoid unphysical screening from periodically repeated layers which otherwise leads to significant band gap reductions.
Importantly, the use of a truncated Coulomb interaction leads to much slower k-point convergence because of the strong q-dependence of the 2D di electric function around q = 0. We alleviate this problem by using an analytical expression for the screened interaction when performing the BZ int egral around q = 0 [106]. This allows us to obtain well converged results with a relatively low k-point density of 5.0/Å −1 (corre sponding to Figure 9. Band structure (left) and projected density of states (right) for VO 2 (top) and WS 2 (bottom) in the MoS 2 prototype, calculated with the PBE xc-functional. The z-component of the spin is indicated by the color code on the band structure.
12 × 12 k-points for MoS 2 ). For example, with this setting the G 0 W 0 band gap of MoS 2 is converged to within 0.05 eV, see figure 11 (left). In comparison, standard BZ sampling with no special treatment of the q = 0 limit, requires around 40 × 40 k-points to reach the same accuracy. Figure 10 (right) shows the PBE and G 0 W 0 band structures of WS 2 (including SOC). The G 0 W 0 selfenergy opens the PBE band gap by 1.00 eV and the HSE gap by 0.47 eV, in good agreement with previous stud-ies [112]. We note in passing that our previously published G 0 W 0 band gaps for 51 monolayer TMDCs [36] are in good agreement with the results obtained using the workflow described here. The mean absolute error between the two data sets is around 0.1 eV and can be ascribed to the use of PBE rather than LDA as starting point and the use of the analytical expression for W around q = 0.
A detailed comparison of our results with previously published G 0 W 0 data is not meaningful because   of the rather large differences in the employed implementations/parameter settings. In particular, most reported calculations do not employ a truncated Coulomb interaction and thus suffer from spurious screening effects, which are then corrected for in different ways. Moreover, they differ in the amount of vacuum included the supercell, the employed k-point grids and basis sets, the in-plane lattice constants, and the DFT starting points. For example, published values for the QP band gap of monolayer MoS 2 vary from from 2.40 to 2.90 eV [113][114][115][116][117][118][119][120] (see [119] for a detailed overview). The rather large variation in published GW results for 2D materials is a result of the significant numerical complexity of such calculations and underlines the importance of establishing large and consistently produced benchmark data sets like the present.
For bulk materials, self-consistency in the Green's function part of the self-energy, i.e. the GW 0 method, has been shown to increase the G 0 W 0 band gaps and improve the agreement with experiments [121]. The trend of band gap opening is also observed for 2D materials [106,120,120,122], however, no systematic improvement with respect to experiments has been established [122]. For both bulk and 2D materials, the fully self-consistent GW self-energy systematically overestimates the band gap [121,122] due to the neglect of vertex corrections [122,123]. In G 0 W 0 the neglect of vertex corrections is partially compensated by the smaller band gap of the non-interacting Kohn-Sham Green's function compared to the true interacting Green's function. In this case, the vertex corrections corrections will affect mainly the absolute position of the bands relative to vacuum while the effect on the band gap is relatively minor [122].
In table 4 we compare calculated band gaps from C2DB with experimental band gaps for three monolayer TMDCs and phosphorene. The exper imental data has been corrected for substrate interactions [122,124], but not for zero-point motion, which is expected to be small (<0.1 eV). The G 0 W 0 results are all within 0.2 eV of the experiments. A further (indirect) test of the G 0 W 0 band gaps against exper imental values is provided by the comparison of our BSE spectra against experimental photoluminescence data in table 7, where we have used a G 0 W 0 scissors operator. Finally, we stress that the employed PAW potentials are not norm-conserving, and this can lead to errors for bands with highly localised states (mainly 4f and 3d orbitals), as shown in [108]. Inclusion of vertex corrections and use of norm conserving potentials will be the focus of future work on the C2DB.

Band extrema
For materials with a finite band gap, the positions of the valence band maximum (VBM) and conduction band minimum (CBM) within the BZ are identified together with their energies relative to the vacuum level. The latter is defined as the asymptotic value of the electrostatic potential, see figure 12. The PBE electrostatic potential is used to define the vacuum level in the non-selfconsistent HSE and G 0 W 0 calculations. For materials with an out-of-plane dipole moment, a dipole correction is applied during the selfconsistent DFT calculation, and the vacuum level is defined as the average of the asymptotic electrostatic potentials on the two sides of the structure. The PBE vacuum level shift is also stored in the database.

Fermi surface
The Fermi surface is calculated using the PBE xcfunctional including SOC for all metallic compounds in the database. Based on a ground state calculation with a k-point density of at least 20/Å −1 , the eigenvalues are interpolated with quadratic splines and plotted within the first BZ. Figure 13 (left) shows an example of the Fermi surface for VO 2 -MoS 2 with color code indicating the out-of-plane spin projection S z . The band structure refers to the ferromagnetic ground state of VO 2 -MoS 2 , which has a magnetic moment of 0.70 µ B per unit cell, characterised by alternating spinpolarised lobes with S z = ±1.

Effective masses
For materials with a finite PBE gap, the effective electron and hole masses are calculated from the PBE eigenvalues; initially these are calculated on an ultrafine k-point mesh of density 45.0/Å −1 uniformly distributed inside a circle of radius 0.015 Å −1 centered at the VBM and CBM, respectively. The radius is chosen to be safely above the noise level of the calculated eigenvalues but still within the harmonic regime; it corresponds to a spread of eigenvalues of about 1 meV within the circle for an effective mass of 1 m 0 . For each band within an energy window of 100 meV above/below the CBM/VBM, the band Even though the masses represent the second derivative of the band energies, we have found that the inclusion of 3rd order terms stabilises the fitting procedure and yields masses that are less sensitive to the details of the employed k-point grids. For each band the mass tensor is diagonalised to yield heavy and light masses in case of anisotropic band curvatures. The masses (in two directions) and the energetic splitting of all bands within 100 meV of the band extremum are calculated both with and without SOC and stored in the database. Other approaches exist for calculating effective masses, such as k · p perturbation theory (see e.g. [128] and references therein); the present scheme was chosen for its simplicity and ease of application to a wide range of different materials.
In addition to the effective masses at the VBM and CBM, the exciton reduced mass is calculated by applying the above procedure to the direct valence-conduction band transition energies, For direct band gap materials the exciton reduced mass is related to the electron and hole masses by 1/µ ex = 1/m * e + 1/m * h , but in the more typical case of indirect band gaps, this relation does not hold.
As an example, figure 14 shows a zoom of the band structure of SnS-GeSe around the VBM and CBM (upper and lower panels). The second order fits to the band energies (extracted from the fitted 3rd order polynomial) are shown by red dashed lines. It can be seen that both the conduction and valence bands are anisotropic leading to a heavy and light mass direction (left and right panels, respectively). The valence band is split by the SOC resulting in two bands separated by ∼ 10 meV and with slightly different curvatures. The conduction band exhibits a non-trivial band splitting in one of the two directions. The peculiar band splitting is due to a Rashba effect arising from the combination of spin-orbit coupling and the finite perpend icular electric field created by the permanent dipole of the SnS structure where Sn and S atoms are displaced in the out of plane direction leading to a sizable vacuum level difference of 1.13 eV, see figure 12. Table 5 shows a comparison between selected effective masses from the C2DB and previously published data also obtained with the PBE xc-correlation functional and including SOC. Overall, the agreement is very satisfactory.

Work function
For metallic compounds, the work function is obtained as the difference between the Fermi energy and the asymptotic value of the electrostatic potential in the vacuum region, see figure 12. The work function is determined for both PBE and HSE band structures (both including SOC) on a uniform k-point grid of density 12.0/Å −1 . Since the SOC is evaluated non-selfconsistently, the Fermi energy is adjusted afterwards based on a charge neutrality condition.

Deformation potentials
For semiconductors, the deformation potentials quantify the shift in band edge energies (VBM or CBM) upon a linear deformation of the lattice. The uniaxial absolute deformation potential along axis i (i = x, y) is defined as [129,130] where ∆E α is the energy shift upon strain and ε ii are the strains in the i-directions. The deformation potentials are important physical quantities as they provide an estimate of the strength of the (acoustic) electron-phonon interaction, see section 3.2.2. Moreover, they are obviously of interest in the context of strain-engineering of band gaps and they can be used to can be used to infer an error bar on the band gap or band edge positions due to a known or estimated error bar on the lattice constant.
The calculation of D α ii is based on a central difference approximation to the derivative. A strain of ± 1% is applied separately in the x and y directions and the ions are allowed to relax while keeping the unit cell fixed. Calculations are performed with the PBE xcfunctional, a plane wave cutoff of 800 eV, and a k-point The change in band energy, ∆E α , is measured relative to the vacuum level. In cases with nearly degenerate bands, care must be taken to track the correct bands as different bands might cross under strain. In this case, we use the expectation value Ŝ z to follow the correct band under strain. Figure 15 shows how the band structure of MoS 2 changes as a function of strain. Both the VBM and the CBM shift down (relative to the vacuum level) when tensile strain is applied in the x direction, but the conduction band shows a much larger shift, leading to an effective band gap closing under tensile strain. Table 6 shows a comparison between the deformation potentials in the C2DB, and literature values obtained using similar methods. There is generally   Phosphorene (zig-zag) Γ 1.24 1.24 [95] 6.56 6.48 [95] Phosphorene (armchair) Γ 0.14 0.13 [95] 0.13 0.12 [95] MAD 0.02 -0.03 -good agreement, and part of the discrepancy can be ascribed to the fact that, in contrast to [131], our numbers include spin-orbit coupling.

Plasma frequencies
The dielectric response of a 2D material is described by its 2D polarisability, α 2D (see section 2.15 for a general introduction of this quantity). For metals, it can be separated into contributions from intraband and interband transitions, i.e. α 2D = α 2D,intra + α 2D,inter . We have found that local field effects (LFEs) are negligible for the intraband component, which consequently can be treated separately and evaluated as an integral over the Fermi surface. Specifically, this leads to the Drude expression for the polarisability in the long wave length limit α 2D,intra (ω) = −ω 2 P,2D /(2πω 2 ) where ω P,2D is the 2D plasma frequency, which in atomic units is given by where v snk = snk|p/m 0 |snk is a velocity matrix element (with m 0 the electron mass), q = q/q is the polarisation direction, s, n, k denote spin, band and momentum indices, and A is the supercell area. The 2D plasma frequency is related to the conventional 3D plasma frequency by ω 2 P,2D (ω) = ω 2 P,3D (ω)L/2 where L is the supercell height.
The plasma frequency defined above determines the intraband response of the 2D metal to external fields.
In particular, it determines the dispersion relation of plasmon excitations in the metal. The latter are defined by the condition ε 2D (ω P ) = 1 + 2πqα 2D (ω P ) = 0, where q is the plasmon wave vector. Neglecting interband transitions (the effect of which is considered in section 3.2.4), the 2D plasmon dispersion relation becomes The plasma frequencies, ω P,2D , for polarisation in the x and y directions, respectively, are calculated for all metals in the C2DB using the linear tetrahedron method [98] to interpolate matrix elements and eigenvalues based on a PBE band structure calculation with a k-point density of 20/Å −1 .

Electronic polarisability
The polarisability tensor α ij is defined by where P i is the i'th component of the induced polarisation averaged over a unit cell and E j is the j'th component of the macroscopic electric field. Using that where ij is the dielectric function. In contrast to the dielectric function, whose definition for a 2D material is not straightforward [119], the polarisability allows for a natural generalisation to 2D by considering the induced dipole moment per unit area, Since the P i is a full unit cell average and P 2D i is integrated in the direction orthogonal to the slab, we have P 2D i = LP i and α 2D ij = Lα ij , where L is the length of the unit cell in the direction orthogonal to the slab.
In the following, we focus on the longitudinal components of the polarisability and dielectric tensors, which are simply denoted by α and ε. These are related to the density-density response function, χ, via the relations where v c is the Coulomb interaction and where Cell is the supercell with volume V . The response function, χ, satisfies the Dyson equation [132] χ = χ irr + χ irr v c χ, where χ irr is the irreducible density-density response function. In the random phase approximation (RPA) χ irr is replaced by the non-interacting response function, χ 0 , whose plane wave representation is given by [133,134] where G, G are reciprocal lattice vectors and Ω is the crystal volume. For all materials in the database, we calculate the polarisability within the RPA for both in-plane and out-of-plane polarisation in the optical limit q → 0. For metals, the interband contribution to the polarisability is obtained from equation (15) while the intraband contribution is treated separately as described in section 2.14. The single-particle eigenvalues and eigenstates used in equation (15) are calculated with PBE, a k-point density of 20/Å −1 (corresponding to a k-point grid of 48 × 48 for MoS 2 and 60 × 60 for graphene), and 800 eV plane wave cutoff. The Dyson equation is solved using a truncated Coulomb potential [105,111] to avoid spurious interactions from neighboring images. We use the tetrahedron method to interpolate the eigenvalues and eigenstates and a peak broadening of η = 50 meV. Local field effects are accounted for by including G-vectors up to 50 eV. For the band summation we include 5 times as many unoccupied bands as occupied bands, which roughly corresponds to an  In figure 16 we show the real and imaginary part of α 2D for the semiconductor MoS 2 . The PBE band gap of this material is 1.6 eV and we see the onset of dissipation at that energy. We also see that the initial structure of Im α is a constant, which is exactly what would be expected from the density of states in a 2D mat erial with parabolic dispersion. Finally, we note that the static polarisability Re α| ω=0 ≈ 6Å, which can easily be read off the figure. The polarisability is also shown for the metallic 1T-NbS 2 where we display the real part with and without the intraband Drude contribution ω 2 P,2D /( ω + iη) 2 .

Optical absorbance
The power absorbed by a 2D material under illumination of a monochromatic light field with polarisation ê is quantified by the dimensionless absorbance: where c is the speed of light. In section 2.15 we gave a prescription for evaluating α 2D in the RPA. However, absorption spectra of 2D semiconductors often display pronounced excitonic effects, which are not captured by the RPA. The Bethe-Salpeter equation (BSE) is a well-known method capable of describing excitonic effects and has been shown to provide good agreement with experimental absorption spectra for a wide range of materials [135]. For materials with finite band gap and up to four atoms per unit cell, we have calculated the RPA and the BSE absorption spectra for electric fields polarised parallel and perpendicular to the layers. The calculations are performed on top of PBE eigenstates and eigenvalues with spin-orbit coupling included and all unoccupied band energies shifted by a constant in order to reproduce the G 0 W 0 quasiparticle gap (the scissors operator method). If the G 0 W 0 gap is not available we use the GLLBSC gap for non-magnetic materials and the HSE gap for magnetic materials (since GLLBSC is not implemented in GPAW for spin-polarised systems). The screened interaction entering the BSE Hamiltonian is calculated within the RPA using a non-interacting response function evaluated from equation (15) with local field effects (i.e. G-vectors) included up to 50 eV and 5 times as many unoccupied bands as occupied bands for the sum over states. We apply a peak broadening of η = 50 meV and use a truncated Coulomb interaction. The BSE Hamiltonian is constructed from the four highest occupied and four lowest unoccupied bands on a k-point grid of density of 20/Å −1 , and is diagonalised within the Tamm-Dancoff approximation. We note that the Tamm-Dancoff approximation has been found to be very accurate for bulk semiconductors [136]. For monolayer MoS 2 we have checked that it reproduces the full solution of the BSE, but its general validity for 2D materials, in particular those with small band gaps, should be more thoroughly tested.
In figure 17 we show the optical absorption spectrum of MoS 2 obtained with the electric field polarised parallel and perpendicular to the layer, respectively. Both RPA and BSE spectra are shown (the in-plane RPA absorbance equals the imaginary part of the RPA polarisability, see figure 16 (left), apart from the factor 4πω and the scissors operator shift). The low energy part of the in-plane BSE spectrum is dominated by a double exciton peak (the so-called A and B excitons) and is in excellent agreement with experiments [55].
In general, calculations of electronic excitations of 2D materials converge rather slowly with respect to k-points due to the non-analytic behavior of the dielectric function in the vicinity of q = 0 [119,137,138]. In figure 18 we show the k-point dependence of the binding energy of the A exciton in MoS 2 obtained as the difference between the direct band gap and the position of the first peak in figure 17. We observe a strong overestimation of the exciton binding energy at small k-point samplings, which converges slowly to a value of ∼0.5 eV at large k-point samplings. For 48 × 48 k-points, corresponding to the k-point sampling used for the BSE calculations in the database, the exciton binding energy is 0.53 eV, whereas a 1/N 2 k extrapolation to infinite k-point sampling gives 0.47 eV (see inset in figure 18). In general, the exciton binding energy decreases with increasing k-point sampling, and thus the exciton binding energies reported in the C2DB might be slightly overestimated. However, since G 0 W 0 band gaps also decrease when the k-point sampling is increased (see figure 11) the two errors tend to cancel such that the absolute position of the absorption peak from BSE-G 0 W 0 conv erges faster than the band gap or exciton binding energy alone. The BSE-G 0 W 0 method has previously been shown to provide good agreement with experimental photoluminescence and absorption measurements on 2D semiconductors. In table 7 we show that our calculated position of the first excitonic peak agree well with experimental observations for four different TMDCs and phosphorene. Experimentally, the monolayers are typically supported by a substrate, which may alter the screening of excitons. However the resulting decrease in exciton binding energies is largely cancelled by a reduced quasiparticle gap such that the positions of the excitons are only slightly red-shifted as compared with the case of pristine monolayers [139].

Database overview
Having described the computational workflow, we now turn to the content of the database itself. We first present a statistical overview of all the materials in the C2DB (i.e. without applying any stability filtering) by displaying their distribution over crystal structure prototypes and their basic properties. We also provide a short list with some of the most stable materials, which to our knowledge have not been previously studied. Next, the predicted stability of the total set of materials is discussed and visualised in terms of the descriptors for thermodynamic and dynamic stability introduced in section 2.4.1. In section 3.2 we analyse selected properties in greater detail focusing on band gaps and band alignment, effective masses and mobility, magnetic properties, plasmons, and excitons. Throughout the sections we explore general trends and correlations in the data and identify several promising materials with interesting physical properties.

Materials
In table 8 we list the major classes of materials currently included in the database. The materials are grouped according to their prototype, see section 2.2.2. For each prototype we list the corresponding space group, the total number of materials, and the number of materials satisfying a range of different conditions. The atomic structure of some of the different prototypes were shown in figure 4. The vast majority of the 2D materials that have been experimentally synthesised in monolayer form are contained in the C2DB (the 55 materials in figure 7 in addition to seven metal-organic perovskites). These materials are marked in the database and a literature reference is provided. Additionally, 80 of the monolayers in the C2DB could potentially be exfoliated from experimentally known layered bulk structures [16][17][18][19]. These materials are also marked and the ID of the bulk compound in the relevant experimental crystal structure database is provided.
To illustrate how all the materials are distributed in terms of stability, we show the energy above the convex hull plotted against ω 2 min in figure 19. It can be seen that the structures naturally sort themselves into two clusters according to the dynamic stability. The points have been colored according to the three levels for dynamic stability introduced in section 2.4. The lower panel shows the distribution of the materials in the grey region on a linear scale. While most of the experimentally known materials (red and black dots) have high dynamic stability, a significant part of them fall into the medium stability category. The marginal distributions on the plot show that the more dynamically stable materials are also more thermodynamically stable. The mean energy above the convex hull is 0.12 eV for the materials with high dynamical stability, while it is 0.25 eV for the others.
In table 9 we show the key properties of a selected set of stable materials, distributed across a variety of different crystal structure prototypes. To our knowledge, these materials are not experimentally known, and they are therefore promising candidates further study and experimental synthesis.

Properties: example applications
In the following sections we present a series of case studies focusing on different properties of 2D materials including band gaps and band alignment, effective masses and mobility, magnetic order, plasmons and excitons. The purpose is not to provide an in-depth nor material specific analysis, but rather to explore trends and correlations in the data and showcase some potential applications of the C2DB. Along the way, we report some of the novel candidate materials revealed by this analysis, which could be interesting to explore closer in the future.

Band gaps and band alignment
The band gaps and band edge positions of all semiconductors and insulators in the C2DB have been calculated with the PBE, HSE06, and GLLBSC xc-functionals while G 0 W 0 calculations have been performed for the ∼250 simplest materials. The relatively large size of these datasets and the high degree of consistency in the way they are generated (all calculations performed with the same code using same PAW potentials and basis set etc) provide a unique opportunity to benchmark the performance of the different xcfunctionals against the more accurate G 0 W 0 method. Figure 20 compares the size and center of the band gaps obtained with the density functionals to the G 0 W 0 results. Relative to G 0 W 0 the PBE functional underestimates the gaps by 45%, i.e. on average the PBE values must be scaled by 1.83 to reproduce the G 0 W 0 results. The HSE06 band gaps are closer to G 0 W 0 but are nevertheless systematically underesti-mated by more than 20%. On the other hand, GLLBSC shows very good performance with band gaps only 2% smaller than G 0 W 0 on average. Table 10 shows the mean absolute deviations of the DFT methods relative to G 0 W 0 . We note that although GLLBSC provides an excellent description of the G 0 W 0 band gaps on average the spread is sizable with a mean absolute deviation of 0.4 eV.
We note a handful of outliers in figure 20 with large HSE band gaps compared to PBE and G 0 W 0 . For one of these, namely the ferromagnetic CoBr 2 -CdI 2 , we obtain the band gaps: 0.  agreement with the C2DB results. It should be interesting to explore the reason for the anomalous behavior of the HSE band gap in these materials. Compared to the band gaps, the gap centers predicted by PBE and HSE06 are in overall better agreement with the G 0 W 0 results. This implies that, on average, the G 0 W 0 correction of the DFT band energies is symmetric on the valence and conduction band. In contrast, the GLLBSC predicts less accurate results for the gap center. This suggests that an accurate and efficient prediction of absolute band energies is obtained by combining the GLLBSC band gap with the PBE band gap center.
Next, we consider the band alignment at the interface between different 2D materials. Assuming that the bands line up with a common vacuum level and neglecting hybridisation/charge transfer at the interface, the band alignment is directly given by the VBM and CBM positions relative to vacuum.
We focus on pairs of 2D semiconductors for which the G 0 W 0 band alignment is either Type II (∆E > 0) or Type III (∆E < 0), see figure 21 (left). Out of approximately 10 000 bilayers predicted to have Type II band alignment by G 0 W 0 , the PBE and HSE06 functionals predict qualitatively wrong band alignment (i.e. Type III) in 44% and 21% cases, respectively (grey shaded areas). In particular, PBE shows a sizable and systematic underestimation of ∆E as a direct consequence of the underestimation of the band gaps in both monolayers.

Effective masses and mobilities
The carrier mobility relates the drift velocity of electrons or holes to the strength of an applied electric field and is among the most important parameters for a semiconductor material. In general, the mobility is a sample specific property which is highly dependent on the sample purity and geometry, and (for 2D materials) interactions with substrate or embedding layers. Here we consider the phonon-limited mobility, which can be considered as the intrinsic mobility of the material, i.e. the mobility that would be measured in the absence of any sample specific-or external scattering sources.
The effective masses of the charge carriers have been calculated both with and without SOC for ∼600 semiconductors. Figure 22 shows the electron mass plotted against the hole mass. The data points are scat- tered, with no clear correlation between the electron and hole masses. Overall, the electron masses are generally slightly smaller than the hole masses. The mean electron mass is 0.9 m 0 , while the mean hole mass is 1.1 m 0 , and 80% of the electron masses are below m 0 while the fraction is only 65% for the holes. This is not too surprising, since, on average, the energetically lower valence band states are expected to be more localised and thus less dispersive than the conduction band states.
The right panel of figure 22 shows the effective mass for electrons and holes plotted as a function of the inverse band gap. It can be seen that there is no clear correlation between the two quantities, which is confirmed by calculating the cross-correlation coefficient: for both electrons and holes it is less than 0.02. This provides empirical evidence against the linear relation between effective masses and inverse band gaps derived from k · p perturbation theory. The relation is based on the assumption that the perturbative expansion is dominated by the conduction and valence band and that the momentum matrix element between these states, u c p u v , does not vary too much as function of the considered parameter (here the type of material). These assumptions clearly do not hold across a large set of different semiconductors. If we focus on a specific class of materials, e.g. sulfides in the MoS 2 structure indicated by the highlighted symbols, we see a slightly improved trend but still with significant fluctuations. Table 9. Key properties of selected stable materials in the C2DB, which have not been previously synthesised. The calculated properties are the magnetic state, formation energy, energy above the convex hull, work function, PBE gap and and the nature of the gap (direct or indirect). If one assumes energetically isolated and parabolic bands, the intrinsic mobility limited only by scattering on acoustic phonons can be estimated from the Takagi formula [147], Here i refers to the transport direction, ρ is the mass density, v i is the speed of sound in the material, m * i is the carrier mass, m * d is the equivalent isotropic densityof-state mass defined as m * d = m * x m * y , and D i is the deformation potential. We stress that the simple Takagi formula is only valid for temperatures high enough that the acoustic phonon population can be approximated by the Rayleigh-Jeans law, n ≈ ω ac /k B T , but low enough that scattering on optical phonons can be neglected.
For the semiconductors in the C2DB we have found that the denominator of equation (17) varies more than the numerator. Consequently, a small product of deformation potential and effective mass is expected to correlate with high mobility. Figure 23 shows the deformation potential plotted against the carrier mass for the valence and conduction bands, respectively. The shaded area corresponds, somewhat arbitrarily, to the region for which m * i D i < m 0 (1 eV). The 2D semiconductors which have been synthesised in monolayer form are indicated with orange symbols while those which have been used in field effect transistors are labeled. Consistent with experimental findings, phosphorene (P) is predicted to be among the materials with the highest mobility for both electrons and holes.
Interestingly, a number of previously unknown 2D materials lie in this shaded region and could be candidates for high mobility 2D semiconductors. Table 11 lists a few selected materials with high intrinsic mobility according to equation (17), which all have 'high' overall stability (see section 2.4.1). In the future, it will be interesting to explore the transport properties of these candidate materials in greater detail.
To put the numbers in table 11 to scale, we consider the well studied example of MoS 2 . For this material we obtain an electron mobility of 240 cm 2 V −1 s −1 while a full ab initio calculation found a phonon-limited mobility of 400 cm 2 V −1 s −1 (in good agreement with experiments on hBN encapsulated MoS 2 [148]), with the acoustic phonon contribution corre sponding to a mobility of 1000 cm 2 V −1 s −1 . Similarly, for the series MX 2 (M = W, Mo, X = S, Se), we calculate roomtemper ature electron mobilities between 200 cm 2 V −1 s −1 and 400 cm 2 V −1 s −1 , which are all within 50% of the ab initio results [149]. Presumably, as in the case for MoS 2 , the good quantitative agreement is partly a result of error cancellation between an overestimated acoustic phonon scattering and the neglect of optical phonon scattering. Importantly, however, the relative ordering of the mobilities of the four MX 2 mono layers is correctly predicted by equation (17) for all but one pair (MoS 2 and WSe 2 ) out of the six pairs. These results illustrate that equation (17) should only be used for 'order of magnitude' estimates of the mobility but that relative comparisons of mobilities in different materials are probably reliable.

Magnetic properties
Recently, a single layer of CrI 3 was reported to exhibit ferromagnetic order with a Curie temperature of 45 K [12]. This comprises the first example of a pure 2D material exhibiting magnetic order and there is    currently an intense search for new 2D materials with magnetic order persisting above room temperature [150][151][152].
For 2D materials, magnetic order will only persist at finite temperatures in the presence of magnetic anisotropy (MA). Indeed, by virtue of the Mermin-Wagner theorem, magnetic order is impossible in 2D unless the rotational symmetry of the spins is broken [153]. A finite MA with an out of plane easy axis breaks the assumption behind the Mermin-Wagner theorem and makes magnetic order possible at finite temperature. The critical temperature for magnetic order in 2D materials will thus have a strong dependence on the anisotropy.
The MA originates from spin-orbit coupling and is here defined as the energy difference between inplane and out-of-plane orientation of the magnetic moments, see equation (4). With our definition, a negative MA corresponds to an out-of-plane easy axis. We note that most of the materials in the C2DB are nearly isotropic in-plane. Consequently, if the easy axis lies in the plane, the spins will exhibit an approximate in-plane rotational symmetry, which prohibits magnetic order at finite temperatures. Since spinorbit coupling becomes large for heavy elements, we generally expect to find larger MA for materials containing heavier elements. In general the magnitude of the MA is small. For example, for a monolayer of CrI 3 with a Curie temper ature of 45 K [12] we find a MA of -0.85 meV per Cr atom in agreement with previous calculations [154]. Although small, the MA is, however, crucial for magnetic order to emerge at finite temperature.
In figure 24 we show the magnitude of the magn etic anisotropy (red triangles) and the magnetic moment per metal atom (blue triangles) averaged over all materials with a given chemical composition. The plot is based on data for around 1200 materials in the medium to high stability categories (see table 2) out of which around 350 are magnetic. It is interesting to note that while the magn etic moment is mainly determined by the metal atom, the MA depends strongly on the non-metal atom. For example, the halides (Cl, Br, I) generally exhibit much larger MAs than the chalcogenides (S, Se, Te). Overall, iodine (I) stands out as the most significant element for a large MA while the 3d metals Cr, Mn, Fe, Co are the most important elements for a large magnetic moment. Since the MA is driven by spin-orbit coupling (SOC) and the spin is mainly located on the metal atom, one would expect a large MA to correlate with a heavy metal atom. However, it is clear from the figure that it is not essential that the spin-carrying metal atom should also host the large SOC. For example, we find large MA for several 3d metal-iodides despite of the relatively weak SOC on the 3d metals. This shows that the MA is governed by a rather complex interplay between the spins, orbital hybridisation and crystal field.
A selection of materials predicted to have high overall stability (see section 2.4.1) and high outof-plane magnetic anisotropy ( MA < −0.3 meV/ Table 11. Key transport properties of selected materials with high intrinsic room-temperature mobility according to equation (17). All the materials shown have 'high' overall stability as defined in section 2.4.1. µ high is the larger value of the mobilities in the x or y directions, m * is the corresponding effective mass, and µ high /µ low is the ratio of the mobilities in the two directions.  Figure 24. Absolute magnetic anisotropy (red) and magnetic moment (blue) averaged over all materials in the C2DB with a given composition. The two red boxes highlight the halides and 3d metals, respectively. magn etic atom) is listed in table 12. We find several semiconductors with anisotropies comparable to CrI3 and some metals with higher values. If we also look at materials with medium overall stability, we find semiconductors with anisotropies up to 2 meV/atom. It is likely that some of these materials will have Curie temper atures similar to, or even higher than, CrI 3 .
In addition to the MA, the critical temperature depends sensitively on the magnetic exchange couplings. We are presently developing a workflow for systematic calculation of exchange coupling constants, which will allow us to estimate the Curie temperature of all the magnetically ordered 2D materials. The database contains several 2D materials with antiferromagnetic order. As a note of caution, we mention that the magnetic interactions in AFM materials typically arise from the super-exchange mechanism, which is poorly described by PBE and requires a careful verification using a PBE + U scheme [155].

Plasmons
The unique optical properties of 2D materials make them highly interesting as building blocks for nanophotonic applications [156,157]. Many of these applications involve electron rich components which can capture, focus, and manipulate light via plasmons or plasmon-polaritons. Graphene sheets can host plasmons that are long lived, can be easily tuned via electrostatic or chemical doping, and can be used to confine light to extremely small volumes [158].  Figure 25. (left) Plasmon dispersion relations for the unscreened (i.e. intraband) and true plasmons, ω P and ω true P , respectively, for NbS 2 in the H phase (the MoS 2 crystal structure prototype). This is compared to the full first principles calculations of the plasmons in NbS 2 by Andersen et al (data points) [159]. (right) The in-plane averaged true plasmon frequency versus the unscreened plasmon frequency for all metals in C2DB at a plasmon wavelength of λ P = 50 nm corresponding to q 0 in the left panel. The data points are colored by the overall stability level as defined in section 2.4.1, and the straight line corresponds to ω P = ω true P .
However, due to the limited charge carrier density achievable in graphene, its plasmons are limited to the mid-infrared regime. Here we show that some metallic monolayers support plasmons with significantly higher energies than graphene and could potentially push 2D plasmonics into the optical regime. Figure 25 (left) shows the plasmon dispersion for monolayer NbS 2 in the MoS 2 crystal structure. The effect of interband transitions on the plasmon is significant as can be seen by comparison to the pure intraband plasmon ( ω P ). The true plasmon energies are obtained from the poles of the (long wave length limit) dielectric function including the interband transitions, = 1 + 2πq(α 2D,intra + α 2D,inter ), yielding ω true P = ω P,2D q 1/2 [1 + 2πqα 2D,inter (ω true P )] −1/2 . For simplicity we ignore the frequency dependence of the interband polarisability, i.e. we set α 2D,inter (ω true P ) ≈ α 2D,inter (ω = 0), which should be valid for small plasmon energies (far from the onset of interband transitions). The validity of this approximation is confirmed by comparing to the full ab initio calculations of Andersen et al (blue dots) which include the full qand ω-dependence [159]. The figure shows that interband screening generally reduces the plasmon energy and becomes increasingly important for larger q. Figure 25 (right) shows the in-plane averaged true plasmon energy of all metals in the C2DB plotted against the intraband plasmon energy at a fixed plasmon wavelength of λ P = 50 nm (corresponding to q 0 at the dashed vertical line in the left panel). For comparison, the plasmon energy of freestanding graphene at this wavelength and for the highest achievable doping level (E F = ±0.5 eV relative to the Dirac point) is around 0.4 eV. The data points are colored according to the overall stability level as defined in section 2.4.1. Table 13 shows a selection of the 2D metals with 'high' overall stability (see section 2.4.1) and large plasmon frequencies. We briefly note the interesting band structures of the metals in the FeOCl prototype (not shown) which exhibits band gaps above or below the partially occupied metallic band(s), which is likely to lead to reduced losses in plasmonic applications [160].
A detailed study of the plasmonic properties of the lead candidate materials will be published elsewhere. However, from figure 25 (right) it is already clear that several of the 2D metals have plasmon energies around 1 eV at λ P = 50 nm, which significantly exceeds the plasmon energies in highly doped graphene.

Excitons
Two-dimensional materials generally exhibit pronounced many-body effects due to weak screening and strong quantum confinement. In particular, exciton binding energies in monolayers are typically an order of magnitude larger than in the corresponding layered bulk phase and it is absolutely crucial to include excitonic effects in order to reproduce experimental absorption spectra.
The electronic screening is characterised by the in-plane 2D polarisability, see equation (12). For a strictly 2D insulator, the screened Coulomb potential can be written as W 2D (q) = v 2D c (q)/ 2D (q) with 2D (q) = 1 + 2πα 2D q and v 2D c (q) = 2π/q is the 2D Fourier transform of the Coulomb interaction. The q-dependence of 2D indicates that the screening is non-local, i.e. it cannot be represented by a qindependent dielectric constant, and that Coulomb interactions tend to be weakly screened at large distances (small q vectors) [119,161,162]. This is in sharp contrast to the case of 3D semiconductors/ insulators where screening is most effective at large distances where its effect can be accounted for by a qindependent dielectric constant. For a two-band model with isotropic parabolic bands, the excitons can be modeled by a 2D Hydrogen-like (Mott-Wannier) Hamiltonian where the Coulomb interaction is replaced by W = 1/ r and the electron mass is replaced by a reduced excitonic mass µ ex derived from the curvature of conduction and valence bands. It has previously been shown that the excitonic Rydberg series of a 2D semiconductor can be accurately reproduced by this model if the dielectric constant, ε, is obtained by averaging 2D (q) over the extent of the exciton in reciprocal space [163]. For the lowest exciton (n = 1), the binding energy can then be expressed as It has furthermore been demonstrated that this expression gives excellent agreement with a numerical solution of the Mott-Wannier model employing the full q-dependent dielectric function, 2D (q) = 1 + 2πα 2D q, for 51 transition metal dichalcogenides [163]. We note that the previous calculations were based on LDA and we generally find that the PBE values for α 2D obtained in the present work are 10-20% smaller compared with LDA.
In figure 26, we compare the binding energy of the lowest exciton obtained from BSE-PBE with G 0 W 0 scissors operator and the Mott-Wannier model equation (18), respectively. Results are shown for the 194 Table 13. Selection of 2D metals with high plasmon energies ω true P for a plasmon wavelength of λ P = 50 nm. The interband screening α 2D,inter at ω = 0 and the in-plane averaged 2D plasma frequency ω P,2D , which are required to reproduce ω true P , are also reported. non-magnetic semiconductors out of the total set of ∼ 250 materials for which BSE calculations have been performed. We focus on the optically active zeromomentum excitons and compute the exciton masses by evaluating the curvature of the band energies at the direct gap, see section 2.11. For anisotropic materials we average the heavy and light exciton masses as well as the x and y components of the polarisability, α 2D , to generate input parameters for the isotropic model equation (18). Although a clear correlation with the BSE results is observed, it is also evident that the Mott-Wannier model can produce significant errors. The mean absolute deviation between BSE and the model is 0.28 eV for all materials and 0.20 eV for the subset of transition metal dichalcogenides (TMDCs). Furthermore, the Mott-Wannier model seems to overestimate E B for more strongly bound excitons while the opposite trend is seen for weakly bound excitons. As explained below these trends are a consequence of systematic errors in the Mott-Wannier model which can be traced to two distinct sources.
(i) Weak screening: If α 2D is small (on the order of 1 Å), the exciton becomes strongly localised and the orbital character of the states comprising the exciton plays a significant role. In general, the Mott-Wannier model tends to overestimate the exciton binding energy in this case as can be seen from the relatively large deviation of points with model binding energies >2.0 eV in figure 26. The overestimated binding energy results from the homogeneous electron and hole distributions implicitly assumed in the Mott-Wannier model. In reality, the short range variation of the electron and hole distributions is determined by the shape of the conduction and valence band states. In general these will differ leading to a reduced spatial overlap of the electron and hole and thus a lower Coulomb interaction. Materials with small band gaps often exhibit hyperbolic rather than parabolic band structures in the vicinity of the band gap. This typically happens in materials with small band gaps such as BSb in the BN prototype. In figure 26 these materials can be identified as the cluster of points with model binding energies <0.25 eV and BSE binding energies >0.25 eV. A similar situation occurs if the conduction and valence bands flatten out away from the band gap region. In both of these cases, the excitons tend to be delocalised over a larger area in the Brillouin zone than predicted by the parabolic band approximation of the Mott-Wannier model. Typically, such delocalisation will result in larger binding energies than predicted by the model. For example, FeI 2 in the CdI 2 prototype exhibits shallow band minima in a ring around the Γ-point and has a BSE binding energy of 1.1 eV and a model binding energy of 0.5 eV because the model assumes that the exciton will be located in the vicinity of the shallow minimum (and thus more delocalised in real space). A detailed inspection reveals that such break down of the parabolic band approximation is responsible for most of the cases where the model underestimates the binding energy.
Other sources of errors come from contributions to the exciton from higher/lower lying bands, i.e. break down of the two-band approximation, and anisotropic exciton masses not explicitly accounted for by equation (18).
Based on this comprehensive and unbiased assessment of the Mott-Wannier model, we conclude that while the model can be useful for understand- ing trends and qualitative properties of excitons, its quanti tative accuracy is rather limited when applied to a broad set of materials without any parameter tuning.
For quantitative estimates α 2D should not be too small (certainly not less than 2 Å) and the the validity of the effective mass approximation should be carefully checked by inspection of the band structure. It has been argued that there should exist a robust and universal scaling between the exciton binding energy and the quasiparticle band gap of 2D semiconductors, namely E B ≈ E gap /4 [164]. This scaling relation was deduced empirically based on BSE-GW calculations for around 20 monolayers and explained from equation (18) and the relation α 2D ∝ 1/E gap from k · p perturbation theory. Another work observed a similar trend [165] but explained it from the 1/E gap dependence of the exciton effective mass expected from k · p perturbation theory. Based on our results we can completely refute the latter explanation (see figure 22 (right)). In figure 26 (right) we show the exciton binding energy plotted versus the direct PBE and G 0 W 0 band gaps, respectively. While there is a correlation, it is by no means as clear as found in [164].

Conclusions and outlook
The C2DB is an open database with calculated properties of two-dimensional materials. It currently contains more than 1500 materials distributed over 32 different crystal structures. A variety of structural, elastic, thermodynamic, electronic, magnetic and optical properties are computed following a high-throughput, semi-automated workflow employing state of the art DFT and many-body methods. The C2DB is growing continuously as new structures and properties are being added; thus the present paper provides a snapshot of the present state of the database. The C2DB can be browsed online using simple and advanced queries, and it can be downloaded freely at https://c2db.fysik.dtu.dk/ under a Creative Commons license.
The materials in the C2DB comprise both experimentally known and not previously synthesised structures. They have been generated in a systematic fashion by combinatorial decoration of different 2D crystal lattices. The full property workflow is performed only for structures that are found to be dynamically stable and have a negative heat of formation. We employ a liberal stability criterion in order not to exclude potentially interesting materials that could be stabilised by external means like substrate interactions or doping even if they are unstable in freestanding form. As an important and rather unique feature, the C2DB employs beyond-DFT methods, such as the many-body GW approximation, the random phase approx imation (RPA) and the Bethe-Salpeter equation (BSE). Such methods are essential for obtaining quantitatively accurate descriptions of key properties like band gaps and optical spectra. This is par ticularly important for 2D materials due the weak dielectric screening in reduced dimensions, which tends to enhance many-body effects. For maximal transparency and reproducibility of the data in the C2DB, all relevant parameters have been provided in this paper. Additionally, all scripts used to generate the data are freely available for download under a GPL license.
Beyond its obvious use as a look-up table, the C2DB offers access to numerous well documented, high-quality calculations, making it ideally suited for benchmarking and comparison of different codes and methodologies. The large set of different available properties makes the C2DB interesting as a playground for exploring structure-property relations and for applying and advancing machine learning approaches in materials science. Moreover, the C2DB should be useful as a stepping stone towards the development of theoretical models for more complex 2D structures such as van der Waals heterostructures (see below).
As reported in this work, based on the combinatorial screening approach, we have identified a number of new, potentially synthesisable 2D materials with interesting properties including ferromagnets with large magnetic anisotropy, semiconductors with high intrinsic carrier mobility, and metals with plasmons in the visible frequency range. These predictions are all based on the computed properties of the perfect crystalline materials. While the pristine crystal constitutes an important baseline reference it remains an idealised model of any real material. In the future, it would be interesting to extend the database to monolayers with adsorbed species and/or point defects. Not only would this allow for a more realistic assessment of the magnetic and (opto)electronic properties, it would also facilitate the design and discovery of 2D materials for e.g. battery electrodes and (electro)catalysis [166,167].
The C2DB should also be useful as a platform for establishing parametrisations of computationally less expensive methods such as tight-binding models [168] and k · p perturbation theory [128]. Such methods are required e.g. for device modeling, description of magn etic field effects, and van der Waals heterostructures. The database already provides band structures, spin orbit-induced band splittings, and effective masses, which can be directly used to determine model parameters. It would be straightforward to complement these with momentum matrix elements at band extrema for modeling of optical properties and construction of full k · p Hamiltonians. Similarly, the spread functional required as input for the construction of Wannier functions e.g. by the ASE [38] or the Wannier90 [169] packages, could be easily and systematically produced. This would enable direct construction of minimal basis set Hamiltonians and would allow for the calculation of Born charges and piezo electric coefficients as well as certain topological invariants [170]. A workflow to calculate exchange couplings of magnetic 2D materials is currently being developed with the aim of predicting magnetic phase transitions and critical temperatures.
Of specific interest is the modeling of the electronic and optical properties of vdW heterostructures. Due to lattice mismatch or rotational misalignment between stacked 2D layers, such structures are difficult or even impossible to treat by conventional ab initio techniques. Different simplified models have been proposed to describe the electronic bands, including tight-binding Hamiltonians derived from strained lattice configurations [171] and perturbative treatments of the interlayer coupling [172]. In both cases, the data in the C2DB represents a good starting point for constructing such models. The effect of dielectric screening in vdW heterostructures can be incorporated e.g. by the quantum electrostatic heterostructure (QEH) model [173] which computes the dielectric function of the vdW heterostructure from the polarisabilities of the isolated monolayers. The latter are directly available in the C2DB, at least in the long wavelength limit.
Finally, it would be relevant to supplement the current optical absorbance spectra by other types of spectra, such as Raman spectra, infrared absorption or XPS, in order to assist experimentalists in characterising their synthesised samples. The automatic firstprinciples calculation of such spectra is, however, not straightforward and will require significant computational investments.

Acknowledgments
The Center for Nanostructured Graphene is sponsored by the Danish National Research Foundation, Project DNRF103. The project also received funding from the European Unions Horizon 2020 research and innovation programme under Grant Agreement No. 676580 with The Novel Materials Discovery (NOMAD) Laboratory, a European Center of Excellence. This work was also supported by a research grant (9455) from VILLUM FONDEN. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 773122, LIMA).