Ab initio property characterisation of thousands of previously unexplored 2D materials

We perform extensive density functional theory calculations to determine the stability and elementary properties of 4249 previously unexplored monolayer crystals. The monolayers comprise the most stable subset (energy within 0.1 eV/atom of the convex hull) of a larger portfolio of two-dimensional (2D) materials recently discovered using a deep generative model and systematic lattice decoration schemes. The relaxed 2D structures are run through the basic property workflow of the Computational 2D Materials Database (C2DB) to evaluate the dynamical stability and obtain the stiffness tensor, piezoelectric tensor, deformation potentials, Born and Bader charges, electronic band structure, effective masses, plasma frequency, Fermi surface, projected density of states, magnetic moments, magnetic exchange couplings, magnetic anisotropy, topological indices, optical- and infrared polarisability. We provide statistical overviews of the property data and highlight a few specific examples of interesting materials. Our work exposes previously unknown parts of the 2D chemical space and provides a basis for the discovery of 2D materials with specific properties. All data is available in the C2DB.


I. INTRODUCTION
The combination of powerful ab initio modeling codes and artificial intelligence (AI) models is opening new opportunities in materials science.Today, it is possible to determine the atomic and electronic structure of even fairly complex materials within a few hours on a single compute node implying that thousands of materials can be scrutinised in a time span of weeks on standard highperformance computing clusters.Such data can then be curated and organised in databases [1][2][3][4] on which AI models can be trained to establish structure-property relations [5][6][7] .Over the past few years numerous examples of such applications have been demonstrated, and AI has already become an integral part of the computational materials scientist's toolbox when it comes to modeling of structures and predicting their properties [8][9][10][11][12] .Compared to such supervised learning tasks, a far more intricate and challenging application of AI is the unsupervised generation of new materials (composition and structure) with prescribed properties -the most basic such property being the thermodynamic stability of the material.
Within the past couple of years, computational materials generation projects (some of them partly AI-driven) have led to the discovery of hundreds of thousands of previously unknown inorganic bulk crystals with high thermodynamic stability, i.e. at or very close to the socalled convex hull [13][14][15] .It must be expected that many of these compounds can be synthesised, and thus they represent an enormous reservoir of candidate materials some of which could be used to improve the performance of existing technologies or even enable new ones.The in-silico expansion of the set of known inorganic crystals by almost an order of magnitude within a couple of years is a tremendous intellectual achievement.However, for practical purposes knowledge of the structure and composition of stable materials is not very useful in itself because the decision to synthesise and deploy a given material usually requires some presumptions about the material's properties.Therefore, to make computational materials discovery relevant for experiments and applications, the determination of stable crystal structures must be complemented by a characterisation of the most basic materials properties.Two-dimensional (2D) materials represent an emerging class of materials whose unique and unconventional properties make them interesting for both fundamental science and technological applications.There exist around 800 layered bulk materials in experimental crystal structure databases whose layer-layer bonds are predicted by ab initio calculations to be sufficiently weak that single layers can be exfoliated 16,17 .The vast majority of these materials, including the around 100 that have already been produced in mono-or few-layer form 18 , are contained in the Computational 2D Materials Database (C2DB) since 2021 19 .
Recently, we trained a crystal diffusion variational autoencoder (CDVAE) 20 on the most stable crystal structures in the C2DB and used it to generate a large set of new, thermodynamically stable 2D materials 21 .The set of CDVAE-generated crystals was complemented by 2D crystals generated by a more traditional lattice decoration approach where the atoms in the known structures were substituted by chemically similar ones.After removing duplicates and non-2D crystal structures, this resulted in a portfolio of 11.630 previously unknown 2D crystals, which were subsequently relaxed using density functional theory (DFT) calculations.Concurrently with our work, Wang et al 22 discovered around 6500 lowenergy 2D materials by systematically occupying all the arXiv:2402.02783v2[cond-mat.mtrl-sci]17 Jun 2024 Wyckoff positions of selected layered space groups by all possible atoms.
In the present work, we calculate the elementary physical properties of the most stable subset of monolayers resulting from our own structure generation project 21 as well those found by Wang et al 22 .Specifically, we consider all monolayers with an energy above the convex hull of less than 100 meV/atom.This amounts to a total of 4249 materials of which 629 originate from the CDVAE, 2702 were produced by lattice decoration while the remaining 918 come from Wang et al.As a testimony to the good thermodynamic stability and chemical validity of the structures, we find that the majority of them (about 70%) are dynamically stable, i.e. stable against small perturbations of the atom positions and unit cell shape.For this subset of 2759 materials we employ the computational workflow behind the C2DB to compute a wide variety of properties, leaving the more computationally demanding steps of the workflow (e.g.GW quasiparticle band structures, absorbance spectra including excitonic effects, and nonlinear optical properties) for future work.All the structures and properties are available on on the C2DB website where they can be browsed or downloaded.
Prior to the present work, the C2DB (which is currently the largest 2D materials database) contained 1345 monolayers with convex hull energy below 0.1 eV/atom.Thus the new set of monolayers characterised in this work, triples the number of (theoretically) known stable 2D materials.

II. RESULTS
In this section we first introduce the 2D materials investigated in this work and describe the different types of properties calculated.A few selected materials with specific or particularly interesting properties are highlighted along the way for illustrative purposes.Out of the 4249 materials with formation energies within 0.1 eV/atom of the convex hull, 2759 are found to be dynamically stable.Further property calculations are performed for this subset only.A complete list of the properties explored and the number of materials for which each property has been determined, is provided in Table I.The table does not comprise all properties of the C2DB workflow.The remaining properties will be computed at a later stage.
For consistency, all calculations presented in this work employ the xc-functionals PBE and HSE06 (only for band structures); see Method section for more details.It is well known that materials with highly localised, partially filled states (e.g open d-shells) can exhibit strong correlation effects, which are poorly described by PBE (and HSE06).For this reason, we have also performed systematic PBE+U calculations for all the materials containing one or more of the 3d transition metals V, Cr, Mn, Fe, Co, Ni.An analysis of these results will be presented elsewhere and all the results will be available in C2DB.

A. Generation of materials
Recently, we used a crystal diffusion variational autoencoder (CDVAE) to generate 5003 monolayer structures 21 .The CDVAE model was trained on a set of 2615 monolayers from the Computational 2D Materials Database (C2DB) 18,24 with energy above the convex hull, ∆H hull , below 0.3 eV/atom.We shall refer to this set of materials as the "seed structures".The CDVAEgenerated materials were complemented by 14192 monolayers obtained by replacing atoms in the seed structures by chemically similar ones.We shall refer to this procedure as the lattice decoration protocol (LDP).The generated structures were subsequently relaxed using DFT and their heat of formation and energy above the convex hull calculated.In the end 3073 (8579) unique 2D structures were obtained using the CDVAE (LDP) generation models.
Importantly, we found that both the CDVAE and LDP generated crystals inherited the good stability properties of the seed structures.As a side remark we mention that the CDVAE-generated structures showed high degree of structural and chemical diversity extending beyond that of the seed structures and the lattice decorated structures.For more details on the CDVAE and LDP generation processes, along with a comprehensive comparison we refer to Ref. 21 .
For the present work we have selected the 3331 most stable subset of the 2D materials generated by the CDVAE and LDP, namely those with ∆H hull < 0.1 eV/atom.
The materials from Wang et al 22 are all binary and ternary compounds and are primarily generated by a symmetry-based approach where the different Wyckoff positions of a given layer group are occupied by up to three different types of elements.This approach is free of any structural or chemical bias originating from a set of seed structures (as is the case for the LDP and to some extent the CDVAE).However, this symmetry-based generation of materials does also not contain any bias towards producing stable materials (as is the case for both LDP and CDVAE).Consequently, Wang et al. invoke several screening criteria to remove unstable materials before the final DFT relaxation.These include conditions on charge neutrality and electronegativity as well as a pre-relaxation of the structure using a machine learning universal interatomic potential.Additionally, Wang et al. performed atom substitution with all elements from the periodic table for a subset of the generated materials and used a crystal-graph neural network to screen the resulting structures for stability.All the structures generated by Wang et al. with ∆H hull < 0.1 eV/atom have been re-relaxed using the same computational framework as used for the CDVAE and LDP structures (see Methods).In addition, we have removed any structures categorized as non-2D by our dimensionality analysis 24 and 81 structures which were already present in the C2DB or the set of CDVAE/LDP-generated materials.Table I: Properties calculated by the C2DB workflow.The computational method and the criteria used to decide whether the property should be evaluated for a given material is also shown.The symbol '*' indicates that spin-orbit coupling (SOC) is included.All calculations are performed with the GPAW electronic structure code using a plane wave basis set.The variations in count numbers among properties sharing identical criteria stem from convergence problems in the DFT calculations.DS: Dynamical stable.E g,half : Half-metal gap.
The combined set of new materials (presenting no overlap with the original structures in C2DB) comprises 4249 unique crystals with ∆H hull < 0.1 eV/atom.This set of crystals contains 3610 unique (reduced) chemical formulas of which only 356 are present in the original C2DB.In the following we explore the properties of the 4249 new materials and compare them to the subset of structures in the C2DB with ∆H hull < 0.1 eV/atom, referred to as the set of stable original materials.

B. Crystal symmetry
The layer group (the analogue of the space group for a crystal with a 2D lattice) has been determined for all the materials using the algorithm in Ref. 23 .In Fig. 1 we show the histograms of layer group numbers for each of the four sets of materials, namely the original stable materials of the C2DB (original), the materials generated by the deep generative crystal diffusion model (CDVAE), the materials generated by the lattice decoration protocol (LDP), and the materials generated by the symmetrybased approach of Wang et al. (Wang et al.).Note the logarithmic scale.Only materials with formation energies within 0.1 eV/atom of the convex hull are included in the plot.
Not surprisingly, the layer group distribution of the LDP materials is very similar to that of the original materials from which they were generated.The CDVAE structures contain relatively few layer groups while the

C. Thermodynamic stability
All the monolayers considered have formation energies within 0.1 eV/atom of the convex hull.In this work we use a convex hull defined by a reference database consisting of the most stable 1,2, and 3-element solid phases from the Open Quantum Materials Database (OQMD) 25 amended by the monolayers themselves (this ensures that the energy above the hull is always non-negative).The threshold value of 0.1 eV/atom has been chosen to account for the finite accuracy of the PBE xc-functional, and to allow for inclusion of meta-stable crystal structures.There are dozens of examples of meta-stable 2D crystal structures that have been experimentally characterised, e.g. the T'-phases of Mo-based dichalcogenides with ∆H hull ranging from 0.18 to 0.02 eV/atom.
The distribution of heat of formation, ∆H, and energy above the convex hull, ∆H hull , for the new materials and stable original materials is shown in Fig. 2.An example of an LDP-generated monolayer on the convex hull (∆H hull = 0 eV/atom) is MoF 3 as seen on the right panel.This monolayer phase has a formation energy of 0.04 eV/atom below the most stable bulk phase of the same composition.We note that MoF 3 (like all the monolayers studied here) does not have a known layered bulk counterpart and thus could not have been found by analysing bulk crystal structure databases.
In total 735 monolayers are predicted to lie exactly on the convex hull.We stress that in reality a monolayer would never lie exactly on the convex hull because its layered bulk form would always have a lower energy due to the attractive interlayer forces.Such layered bulk crystals are not included in the reference database of bulk crystals used to calculate the convex hull, and therefore the monolayers can end up be exactly on-hull.For monolayers with high thermodynamic stability, like the ones studied here, interlayer interactions will be weak and mostly of the van der Waals type (see next section).Therefore, we can safely assume that the monolayers will remain close to the convex hull even if the layered bulk form was included in the convex hull.On a more technical note, a correct description of the interlayer interactions in vdW bulk crystals require the use of an xcfunctional that can account for dispersive forces.In our work we use the PBE functional that does not capture such forces.Thus, even if the layered bulk form of a given monolayer would be included in the convex hull, it would only be predicted to have a lower energy than the monolayer if the layers form stronger chemical (covalent/ionic) bonds.

D. Exfoliability
In general, there are two main routes to the production of monolayer crystals, namely chemical synthesis (e.g. via atomic layer deposition, chemical vapor deposition, or pulsed laser deposition) or exfoliation from a bulk crystal (either mechanical or via selective chemical etching).The two approaches may be viewed as bottom-up or topdown, respectively.
Common to all the monolayers studied in this work is that they do not have a known (naturally occurring or previously synthesised) layered bulk analogue.Consequently, they would have to be synthesised bottom-up.Here it should be noted that chemical synthesis might result in a film containing several monolayers.In such cases, an exfoliation step is required to isolate a single monolayer.
An important descriptor for evaluating the synthesisability of a monolayer crystal is its exfoliation energy, i.e. the difference in total energy between the isolated monolayer and the layered bulk (per layer) normalised by the surface area.When the parent bulk structure is known, the calculation of the exfoliation energy is in principle straightforward 17 .When this is not the case, the calculation is more involved as various stacking configurations must be explored in order to find the most stable layered structure 26 .Based on an analysis of the exfoliation energies calculated using DFT-vdW for a set of known materials, Mounet et al. 17 concluded that an exfoliation energy of 30 meV/ Å2 is a reasonable upper limit for "easily exfoliable" materials, and that materials with exfoliation energies up to 150 meV/ Å2 are "potentially exfoliable".
The size of the exfoliation energy is largely determined by the chemical inertness of the monolayer (e.g.absence of dangling bonds).Thus a high thermodynamic stability of the monolayer should be indicative of low exfoliation energy.This expectation is supported by our recent work 26 where we perform DFT-vdW calculations of the exfoliation energies of 1052 monolayers in the original C2DB dataset with ∆H hull = 0.2 eV/atom and found that the vast majority of the monolayers had exfoliation energies below 50 meV/ Å2 .Based on these findings it is very likely that the monolayers of the present study will have comparably low exfoliation energies.

E. Dynamical stability
The dynamical stability is assessed by calculating the phonons at the high-symmetry points of the Brillouin zone corresponding to the q-points (0,0), (0,0.5),(0.5,0), and (0.5,0.5) in fractional coordinates.An imaginary frequency signals a phonon instability.Moreover, the stiffness tensor is calculated and diagonalized.A negative eigenvalue signals an instability of the shape of the unit cell.A material is termed dynamically stable if all phonon frequencies and stiffness tensor eigenvalues are real and positive.More details including justification of the scheme can be found in Ref. 7 .
Figure 3 shows the distribution of materials according to the dynamical stability (phonons and stiffness, respectively).The materials have been subdivided according to their origin, i.e.LDP, CDVAE, and Wang et al..For all groups of materials, phonon instabilities occur more frequently than stiffness instabilities.In particular, very few materials that are stable with respect to phonons show stiffness instability (0.2%, 0.8%, and 1.1% for the three groups of materials).A slightly higher percentage (76%) of the LDP materials are found to be dynamically stable with respect to phonons as compared to the CDVAE materials (70%), which in turn have a higher phonon sta-  The percentage of dynamically stable structures produced by the LDP and CDVAE schemes (70-76%) is similar to that of the subset of seed structures with ∆H hull < 0.1 eV/atom (76%) suggesting that the LDP/CDVAE models learn to generate dynamically stable crystals from the seed/training structures.
In total 2759 of the new materials are found to be dynamically stable.The remaining part of the property workflow is limited to this set of materials.

F. Electronic band structure
The electronic band structure has been calculated for all the dynamically stable materials using the PBE xcfunctional.For materials with a finite band gap, the band structure has also been obtained with the HSE06 hybrid functional.In both cases spin-orbit coupling (SOC) is included.
The distribution of the PBE and HSE06 band gaps is shown in Fig. 4 (metals not included).About 30% the new materials generated by LDP and CDVAE are metallic, which is very similar to that of the seed structures (31 %).The fraction of metallic compounds is significantly lower in the set from Wang  ing on the 1982 non-metallic compounds, we find that 26% (CDVAE), 42% (LDP), and 47% (Wang) of the new structures have a direct band gap, see pie charts on Fig. 4.Here we classify band gaps as direct if the direct gap is below or within 5 meV of the indirect gap.Ultra wide band gap materials are important for applications such as insulating dielectrics and ultraviolet photonics.We find 18 new 2D materials with very large band gaps, here defined as the HSE06 band gap exceeding 8 eV.In comparison, the largest band gap in the original set of stable materials in C2DB is 7.49 eV (MgB 2 H 8 ).17 of the new large band gap materials are fluorides and most of them have chemical formula XYF 6 or YF 3 , where X is an alkali metal and Y is a transition metal.Only one of the new large band gap materials B 4 O 6 does not contain fluoride.
In Fig. 5 we show the PBE and HSE06 band structures for four non-metallic compounds selected from the set of new monolayers: AgHgISe, MoF 3 , AlF 3 , and CrFO 2 .All four materials have a direct band gap.
With a band gap of 10.8 eV (HSE06), AlF 3 has the largest band gap among all the new materials.Such large band gap 2D insulators are relevant for several reasons, including as gate dielectrics in field effect transistors employing 2D semiconductors as the channel material 27,28 .As an example of a simple, previously unknown monolayer we show the band structure of MoF 3 .This material lies on the convex hull (∆H hull = 0) and is predicted to have a direct band gap of 1.0 eV (HSE06), which is close to a widely used telecommunication wavelength band (Oband) 29 , making it interesting for photonics and optoelectronic applications in the visible frequency range.As an example of a more complex material we show the band structure of AgHgISe.This material has an energy slightly above the convex hull of (∆H hull = 0.02 eV/atom) and a direct band gap of 1.88 eV (HSE06).Finally, as an example of a new magnetic material, we show the band structure of CrFO 2 .This material also lies on the convex hull and has a direct band gap of 3.65 eV (HSE06).

G. Magnetic materials
All the DFT structure relaxations are performed with spin polarisation starting from an initial ferromagnetic spin configuration.If the final structure has absolute magnetic moments on all atoms below 0.1µ B , the materials is considered non-magnetic and all subsequent property evaluations are performed on-top of a spin-paired ground state calculation.For magnetic materials, a nearest neighbor exchange coupling, J, is derived from the total energy of one specific anti-ferromagnetic configuration.If J > 0, the material is classified as 'ferromagnetic'.If J < 0, the material is classified as 'magnetic'.Note that in the case of J < 0, we do not classify the ma-  terial as 'anti-ferromagnetic', because (1) we have only considered one of many possible anti-ferromagnetic configurations and (2) the true magnetic ground state could be more complex, e.g. a spin spiral.Fig. 6 shows three pie charts depicting the fraction of new materials that are magnetic/non-magnetic and metallic/insulating.Magnetic materials with semiconducting properties are of particular interest for spintronics applications, e.g.spin transistors 30 .Interestingly, the CDVAE-generated structures contain a significantly larger fraction of magnetic, non-metallic materials (16%) than both the LDP (7%) and Wang et al. (6%) structures.
Another interesting class of materials within spintronics are half-metals, which are metallic in one of the spin channels while the other has a band gap 31 .Calculating the half-metal gap was not previously part of the C2DB workflow but we have recently added it and calculated the half-metal gap for all monolayers in C2DB, which are magnetic, metallic, and dynamically stable.In total 254 monolayers satisfy these conditions and of these 74 have a non-zero gap in one of the spin-channels when using the PBE xc-functional.Among these 74 monolayers demonstrating a PBE-based half-metallic gap, further analysis was conducted employing the more accurate HSE06 xc-functional.In total 51 of these monolayers remains half-metallic at the HSE06 level.We note that spin orbit coupling (SOC) may mix the two spin channels making the very concept of a half-metal somewhat blurred.For this reason the calculation of the half-metal gap does not include SOC.The distribution of the halfmetal gaps for both PBE and HSE06 is show in Fig. 7.The half-metallic band gap of three of the monolayers, namely FeCl 2 , FeBr 2 , and FeI 2 , was previously investigated at the HSE06 level in Ref. 31 .Our findings align well with theirs; they reported the (HSE06) half-metallic band gap to be 6.4 eV, 5.5 eV, and 4.0 eV, respectively, while our calculations yield values of 6.1 eV, 5.2 eV, and 4.0 eV for the same materials.In addition to the total magnetic moment and the nearest neighbor exchange coupling, the magnetic anisotropy is also calculated by the workflow for materials with one or two magnetic atoms in the unit cell.It represents the total energy difference between spins aligned in the in-plane (x and y) and out-of-plane (z) directions, respectively.Thus the sign of the magnetic anisotropy determines whether the magnetic material has an easy axis or easy plane.
The magnetic properties of a material can be described using a Heisenberg model (assuming a single magnetic site per unit cell and nearest neighbor interactions) Here J is the nearest neighbor exchange coupling, A  is the single-ion anisotropy (out-of-plane) and B is the anisotropic exchange (out-of-plane).The sums over i and j is restricted to nearest neighbors.The parameters of the Heisenberg model can be obtained from DFT total energies as described in Ref. 32 .Using spin wave theory, the spin wave gap, ∆, i.e. the smallest energy required for a magnetic excitation, is given by 32 where N nn is the number of nearest neighbors in the given crystal and S is the spin.∆ > 0 is required for for out-of-plane magnetism in 2D materials.By combining the exchange coupling and spin wave gap with simple structural parameters, e.g.number of nearest neighbors and lattice type, it is possible to estimate the critical temperature as described in Ref. 32 .In general for a high critical temperature, both the exchange coupling and the spin wave gap should be as large as possible.
In Fig. 8 the exchange coupling, J, is plotted against the spin wave gap, ∆, for the magnetic and non-metallic materials.The two CDVAE generated materials V 2 I 2 SSe and V 2 BrIS 2 exhibit notably high values of J and ∆.These two share the same crystal structure and belong to layer group 11.

H. Dielectric screening
The frequency-dependent 2D optical polarisability, α opt 2D (ω), in the long wave lenght limit (q → 0), is calculated within the random phase approximation (RPA).For metals, a Drude term, is added to the interband polarisability to account for intraband transitions.In this expression the plasma frequency is obtained as an integral of ⟨nk|q•∇|nk⟩ over the Fermi surface, where q is a unit vector in the direction of the electric field 18 In materials with a finite band gap, the atomic lattice can also contribute to the polarisability at frequencies below or comparable to the maximum phonon frequency.To include the contribution to the polarisability from infrared (IR) active phonons, we first determine the atomic Born charges describing the change in the macroscopic polarisation, P, due to displacement of the atom, The IR polarisability can then be determined by combining the Born charges with the eigenmodes and frequencies of the optical phonons at the Γ-point.The total polarisability is then obtain as the sum All quantities in the above equation are tensors.It has been proposed that the dielectric constant of a layered van der Waals (vdW) bulk crystal can be obtained from the polarisability of its constituent monolayers via the relations 33 where d is taken as an effective thickness of the 2D material.We estimate d as the distance between the two outermost atoms of the monolayer plus the vdW radii 34 of the outermost atoms.We have checked that this approximation gives a good qualitative agreement with the DFT calculated distance between monolayers in vdW bilayers, albeit slightly underestimating the thickness 26 .As previously mentioned, 2D vdW materials with good dielectric properties are being actively sought for due to the potential application in 2D electronics, e.g. as gate insulators 27,28,35 .Good field effect transistor gate dielectrics are characterised by a large electronic band gap (to limit leakage currents) and a large out-of-plane static dielectric constant (to minimise gate thickness and threshold voltage).To screen the new 2D materials for gate dielectric candidates we plot the fraction α ⊥ 2D /d against the HSE06 band gap, see Fig. 9.While our results show that Eq. ( 6) is generally quite accurate, we have found that Eq. ( 7) can lead to a diverging or even negative out-of-plane dielectric constant -even when the layer thickness (d) is derived from more accurate DFT calculations or experimental interlayer distances.For this reason we show α ⊥ 2D /d, which directly expresses the ability of the individual 2D layer to screen an electric field.We note in passing that the electronic contribution to the total static polarisability is expected to scale as 1/E gap .).These observations agree with previous findings 33 .We ascribe this different behavior to the larger influence of local field effects on the out-of-plane polarisability.
A few of the new materials with particularly large values of the key quantity E gap α ⊥ 2D /d are highlighted and their atomic structure shown.In particular, CaNaI

I. Piezoelectric tensor
The piezoelectric tensor has been calculated for 858 of the new materials that are dynamically stable, noncentrosymmetric, and have a finite band gap.
The piezoelectric tensor, c, of an insulating crystal is a rank-3 tensor relating the macroscopic polarization, P , to an applied strain.It is non-zero only for crystals lack-  ing an inversion center.In Voigt notation, c is expressed as a 3 × N matrix relating the (x, y, z) components of the macroscopic polarizability to the N independent components of the strain tensor.The piezoelectric tensor is evaluated as a finite difference of the polarization under three independent strains of the unit cell with the atom positions fully relaxed.The polarisation in the periodic directions is calculated as an integral over Berry phases.
The polarization in the non-periodic direction is obtained by direct evaluation of the first moment of the electron density.For 2D materials, the strain tensor has three independent components comprising two linear (xx and yy) and one shear (xy) component.Thus strain can be represented as a 3-vector, and the Piezoelectric tensor as a 3 × 3 matrix, where i = x, y, z and j = xx, yy, xy.As a validation of the computational methodology we mention that the calculated piezoelectric coupling of freestanding MoS 2 is 0.35 nC/m in good agreement with the experimental value of 0.3 nC/m.See Ref. 24 for further details on the computational method.Piezoelectric 2D materials could find applications in nanoscale electro-mechanical devices, mechanicalelectrical energy conversion, and sensing 36,37 .For a high mechanical-electrical energy conversion, the straininduced polarisation should be as large as possible, and the elastic energy as small as possible.The strain direction of maximum polarisation is the eigenvector (e max ) corresponding to the largest eigenvalue of the matrix c † c.The stiffness of the material in this direction is then e † max Ce max .In Fig. 10 we plot the maximum polarisation against the inverse stiffness along the direction of the maximum polarisation for the different materials groups.Some of the materials that look particularly promising mechanical-electrical energy conversion, are highlighted.

J. Topological invariants
The Berry phase spectrum is calculated for 340 of the new materials with a band gap in the range 0 < E PBE gap < 0.7 eV and with less than 13 atoms in the unit cell.The Berry phase spectrum gives the Berry phase, γ n ( k1 ), of an occupied band, n, along the path k2 = 0 → k2 = 2π, where a specific k point in reciprocal space is given by ( k1 /2π)b 1 + ( k2 /2π)b 2 .The topological indices, in particular the Chern number (C), the mirror Chern number (C M ), and the Z 2 invariant (ν), can be determined by inspection of the Berry phase spectrum.We refer to Refs. 19nd 38 for more details on the methodology.
An example of a quantum spin Hall insulator is the CDVAE-generated Ta 2 Te 2 S. The Berry phase spectrum of this materials is shown in Fig. 11 (right) while the electronic band structures calculated with PBE and HSE06 are shown in the left and middle panels, respectively.From the band structure we note that Ta 2 Te 2 S, like graphene 39 and silicene 40 , hosts a Dirac cone at the K point.However, the band gap of Ta 2 Te 2 S (150 meV) is much larger than that of pristine graphene (no gap) and silicene (1.55 meV 41 ), making it more suitable for applications.
The only example among the new materials of an anomalous Hall insulator is the CDVAE-generated magnetic monolayer Mn 2 Br 2 O 3 .Its Berry phase spectrum and PBE electronic band structure is shown in Fig. 12.

III. SUMMARY
We have performed an extensive computational characterisation of 2759 two-dimensional (2D) crystals all of which are predicted to be dynamically and thermodynamically stable, but never have been explored before.Using state of the art ab initio calculations we have determined a variety of basic properties of this previously unknown set materials and made them available in the Computational 2D Materials Database (C2DB).As a result, this work increases the number of stable materials contained in the C2DB by almost a factor of three.
We identify 633 monolayers with a band gap in the semiconductor energy range 0.5-2 eV, of which 156 are direct gaps (results obtained with the HSE06 xc-functional including spin-orbit coupling).We find 406 materials with a magnetic ground state of which 216 have a finite band gap.Among these several exhibit a large magnetic anisotropy, which is an essential requirement for magnetic order at finite temperatures in the 2D limit.In particular, we find a number of semiconducting materials (e.g. the chalcohalides V 2 BrIS 2 and V 2 I 2 SSe) with simultaneously large exchange coupling and spin wave gaps making them candidates for 2D magnets with high Curie temperature.
Our calculations provide a complete characterisation of the linear dielectric properties of the non-metallic monolayers.Specifically, both the electronic and phononic contributions to the polarisability are calculated as function of frequency for in-plane and out-of-plane polarisation directions.Focusing on the static out-of-plane polarisability, which is of greatest relevance for dielectric applications in 2D electronics, we identify a number of promising materials combining a high polarisability with a large band gap, e.g.CaNaI 3 and SrNaI 3 .Not surprisingly, these materials are characterised by a relatively large contribution to the polarisability from out-of-plane optical phonons.
The piezoelectric tensor is calculated for more than 800 non-centrosymmetric and insulating monolayers.Based on a combined analysis of the piezoelectric tensor and the stiffness tensor, we identified a number of promising materials for mechanical-electrical energy conversion.
Finally, we have calculated the Berry phase spectrum of 340 materials with small band gaps and used to identify materials with topologically non-trivial band structures.Here we highlighted Ta 2 Te 2 S, which we predict to be a quantum spin Hall insulator with a gapped Dirac cone, and Mn 2 Br 2 O 3 , which is classified as a ferromagnetic anomalous Hall insulator.

IV. METHOD
The employed C2DB computational workflow was constructed within the Atomic Simulation Recipes (ASR) Python framework 19 and executed using the MyQueue 42 task scheduler.The workflow performs calculations using the GPAW 43 electronic structure code and the Atomic Simulation Environment (ASE) Python library 44 .
All GPAW DFT calculations were performed using an 800 eV plane wave cut-off and a k-point grids with density between 4 and 12 Ådepending on the property calculated.The PBE exchange-correlation functional 45 was used in all calculations, except for the HSE06 46 band structure calculations.Spin-orbit coupling was included in calculations of single-particle band energies and the magnetic anisotropy.An out-of-plane vacuum region of at least 12 Å is included in the calculations to eliminate any spurious interactions between the periodic images of the 2D layers.All the property calculations has previously been shown to be well converged for the given computational setting used 18,24 .The precise computational settings used for each type of property calculation can be found in Ref. 24 .

V. DATA AVAILABILITY
The data that support the findings of this study are openly available.All the crystal structures and their properties are available as a part of C2DB (https: //cmr.fysik.dtu.dk/c2db/c2db.html)

Figure 1 :Figure 2 :
Figure 1: Histograms showing the distribution of the 2D crystals investigated in this work according to their layer group number.Note the logarithmic scale.Histograms are shown separately for the four distinct groups of materials: The original stable materials of the C2DB (original), the materials generated by the deep generative crystal diffusion model (CDVAE), the materials generated by lattice decoration of the original C2DB materials (LDP), and the materials generated by the symmetry-based approach of Wang et al. (Wang et al.).The layer group has been determined with the algorithm in Ref. 23 using a symmetry precision threshold of 0.1 Å.

Figure 4 :
Figure 4: Left: Stacked histograms showing the distribution of band gaps (metals not included) for the newly generated materials and the stable original materials in the C2DB.Band gaps obtained using the PBE and HSE06 xc-functionals are shown in the left and right histograms, respectively.Right: The pie charts show the fraction of the newly generated materials with a direct and indirect PBE band gap, respectively.All band energies include spin-orbit coupling.

Figure 5 :
Figure 5: Examples of band structures.The band structure calculated using the PBE xc-functional with and without spin orbit coupling as well as HSE06 hybrid functional is shown for the four materials AgHgISe, MoF 3 , AlF 3 , and CrFO.The four selected materials have direct band gaps of varying sizes, are dynamically stable, and have energy above the convex hull of at most 0.03 eV/atom (MoF 3 and CrFO lie on the convex hull).CrFO has a ferromagnetic ground state.

Figure 6 :
Figure6: Distribution of the new materials according to their metallic/non-metallic and magnetic/non-magnetic nature.

Figure 7 :
Figure 7: Histogram showing the distribution of half-metal gaps for materials which are magnetic, metallic, and have a gap in one spin channel.The half-metal gaps are calculated using PBE and HSE06, both without spin orbit coupling.

Figure 8 :
Figure 8: Nearest neighbor exchange coupling, J, versus the spin wave gap, ∆, for the magnetic and non-metallic materials.Only materials with ∆ > −1 meV are shown, as a positive spin wave gap is a requirement for magnetic order in 2D.

Figure 9 :
Figure 9: Static out-of-plane polarisability normalized by the estimated van der Waals thickness of the monolayer, α ⊥ 2D /d, versus the HSE06 band gap (metals not included).Results are shown for the new materials from the different datasets and the original set of stable materials in the C2DB.The absence of any correlation between α 2D and the size of the band gap is ascribed to two effects, namely, the contribution from infrared active phonons and the influence of local field effects, which is particularly strong for out-of-plane polarisation.
3 and SrNaI 3 , which shares structure prototype, have exceptionally high α ⊥ 2D /d and a reasonably large band gap of 4.4 eV.This is due to a very large phonon contribution to the polarisability, e.g. in the case of CaNaI 3 , α ⊥,IR 2D = 1.20 Åwhile α ⊥,opt 2D = 0.39 Å. Hexagonal BN, which is commonly used in experimental studies, is also highlighted as is AlF 3 whose band structure is shown in Fig. 5.In general, the materials from Wang et al. contain several candidates with large phonon contributions to α ⊥ 2D .

Figure 10 :
Figure 10: Maximum polarisation, P max , versus the inverse stiffness, 1/C.The stiffness is calculated in the direction of the maximum polarisation.

Figure 11 :Figure 12 :
Figure 11: Electronic band structure calculated with the PBE (left) and HSE06 (middle) xc-functionals, and the Berry phase spectrum (right) of the CDVAE-generated monolayer Ta 2 Te 2 S. The material is a quantum spin Hall insulator (Z 2 index ν = 1) with a HSE06 band gap of 0.15 eV at the K-point.
Overview of the dynamical stability of the new materials for the three material groups.Green (red) colors denote stable (unstable) materials.The test for dynamical stability involves a separate test of the phonon frequencies and the eigenvalues of the stiffness tensor.The diagrams show the result of both tests (inner and outer rims).It can be noted that materials are much more likely to be dynamically unstable due to phonons (distortions of the atomic structure) than the stiffness tensor (instability of the unit cell shape).