Large Piezoelectric Response of van der Waals Layered Solids
Abstract
The bulk piezoelectric response, as measured by the piezoelectric modulus tensor (d), is determined by a combination of charge redistribution due to strain and the amount of strain produced by the application of stress (stiffness). Motivated by the notion that less stiff materials could exhibit large piezoelectric responses, herein we investigate the piezoelectric modulus of van der Waalsbonded quasi2D ionic compounds using firstprinciples calculations. From a pool of 869 known binary and ternary quasi2D materials, we have identified 135 noncentrosymmetric crystals of which 48 systems are found to have d components larger than the longitudinal piezoelectric modulus of AlN (a common piezoelectric for resonators), and three systems with the response greater than that of PbTiO, which is among the materials with largest known piezoelectric modulus. None of the identified materials have previously been considered for piezoelectric applications. Furthermore, we find that large d components always couple to the deformations (shearing or axial) of van der Waals “gaps” between the layers and are indeed enabled by the weak intralayer interactions.
I Introduction
Coupling between the mechanical degrees of freedom and the electric polarization in a solid, a hallmark of piezoelectric materials, has found use in applications that range from sensors, resonators, motors, and actuators, to highresolution ultrasound devices and miniature filters for cellular communications.jaffe2012a; uchino1996a; liu2018double Interestingly, only about 10 piezoelectric materials, including SiO (quartz), LiTaO, LiNbO, PZT (lead zirconate titanate)based, BaTiObased, (K,Na)NbObased, BiTiObased, AlN, and ZnO, are technologically relevant and cover virtually all of these applications.ieeeweigel; wu2016piezotronics; uchino1998piezoelectric Expanding the pool of known piezoelectric materials would help broaden the range of applications and allow earthabundant and nontoxicsaito2004lead replacements for materials that are presently in use. Furthermore, it would offer more cost and performanceeffective choices beyond a relatively limited set of materials that are at present used for their piezoelectric properties.
The need for new piezoelectrics was recognized before and was the main motivation behind recent computational efforts in highthroughput screening of inorganic materials for piezoelectric performance. These include the work on perovskite alloys,armiento_PRB:2014; armiento2011screening and creation of a database of piezoelectric properties of compounds.de2015database Materials with reduced dimensionality, such as 2D mono and multilayers have also been investigated recently.acsnano2017; blonsky2015ab; Li2015; yin2017giant; duerloo2012intrinsic; hinchet2018piezoelectric In virtually all of these works, the quantity of interest was the piezoelectric coefficient tensor , which relates the polarization to strain ,nye1985physical; eijdij (in Voigt notation,voigt2014lehrbuch ) and is typically obtained from firstprinciples calculations.
In this work, we also apply firstprinciples calculations to screen for candidate piezoelectrics, but instead of the piezoelectric coefficient tensor we consider the piezoelectric modulus tensor , which relates the induced polarization to the applied stress, .nye1985physical; eijdij The advantage of using is in that it represents an important figure of merit in a wide range of technological applicationsjaffe2012a; uchino1996a and is a more commonly measured piezoelectric property, especially for materials in the thin film form. Because the induced polarization depends on the applied stress through a combination of the charge redistribution due to strain and the amount of strain that is produced by the applied stress, the piezoelectric modulus tensor depends on both the piezoelectric coefficient tensor and the elastic tensor as:
(1) 
From equation (1), we can infer that large piezoelectric modulus can be expected in materials with large large , but also in systems with low stiffness . Softer (less stiff) materials can exhibit large piezoelectric response, as measured by , compared to stiffer materials with similar . This is illustrated in Figure 2, where we notice less stiff materials overall exhibit larger piezoelectric responses in the corresponding direction.
For this reason, we concentrate our investigation on layered (quasi2D) van der Waals bonded solids, which can be expected to have relatively low stiffness in the outofplane direction due to the presence of weak van der Waals (vdW) interactions between the layers. The questions we are addressing are: (a) whether materials that belong to this class and exhibit strong piezoelectric response () can be found, and (b) if this is true, what is the role of vdW interactions in enabling the strong response. To our knowledge, these questions have not been previously addressed in a systematic way.
Our results confirm the expectations. The search has revealed a number of quasi2D materials with relatively large components. Out of 869 considered binary and ternary layered VdW systems, we have identified 135 noncentrosymmetric crystals. Out of those we find more than one third (48 compounds) exhibit piezoelectric moduli greater than that of AlN (=5.5 pC/N), a commonly used piezoelectric material in resonator applications. In addition, we found three layered systems with piezoelectric moduli even larger than that of PbTiO (=119 pC/N ), another established piezoelectric material known for its very large response. It is important to note that none of these vdW systems have been considered previously for piezoelectric applications.
After performing a thorough analysis of the coupling between various stress and components, we find that in all of these systems, large piezoelectric response is always coupled to the stress components that imply deformations (axial or shear) of the van der Waals “gaps” between the layers. This is consistent with the fact that the softest elastic constants are related to these deformations. Ultimately, our results point to the layered vdW systems as a rich chemical space for finding new piezoelectric materials, and introduce elastic properties as additional design criteria for finding materials with large piezoelectric modulus.
Ii Computational Methodology
In this section, we discuss the details of the computational methodology used in our calculations, broadly divided into three main subsections. The first subsection describes the procedure for identifying quasi2D structures from the Inorganic Crystal Structure Database (ICSD).bergerhoff1983inorganic; belsky2002new Next, we describe the drawbacks of GGA exchange correlation functionals for predicting the properties of layered materials and our approach to overcome this issue. Then, we provide detailed descriptions of our workflow for evaluating the piezoelectric modulus tensors (d) of quasi2D solids.
ii.1 Automated Identification of Quasi2D Materials
An essential component of this work is the identification of layered (quasi2D) vdW materials from ICSD. To accomplish this, we extend the applications of our previouslydeveloped proceduregorai2016computational for automated identification of the binary quasi2D materials from the ICSD to include ternary chemistries. Similar algorithms have been developed and used for the purpose of broad identification of vdW bonded layered systems by others, including the work of Ashton et al.,ashton2017topology Mounet et al.,mounet2016novel and Cheon et al.cheon2017data
Our procedure relies on a slab cutting routine and bond counting. In the first step, we cut out stoichiometric slabs of a certain thickness for all symmetry inequivalent sets of Miller indices within a certain range (). Next, for each slab we find the terminations of its surfaces that minimize the number of broken bonds by translating surface atoms from one side of the slab to the other using appropriate lattice vectors. We then count the (minimal) number of broken bonds, i.e., the undercoordination, of the surface atoms. The condition of quasi lowdimensional crystals then implies the existence of directions for which the corresponding slabs do not have any undercoordinated atom relative to their bulk coordination in the first shell. If there is exactly one such , the material is a layered material with relatively large spatial gaps separating individual layers. If the number of directions is larger than one, then the corresponding systems are of lower dimensionality, quasi1D for two such directions and molecular crystals for larger than two. If there are no such directions, the structure is then a connected 3D structure without large spatial gaps.
In Ref. gorai2016computational we demonstrated the success of our algorithm in searching complex quasi2D materials, including those with layer stacking in oblique directions, materials with corrugated, accordionlike layers, and those with individual layers composed of multiple atomic layers. Using the described automated algorithm, we have in this work considered 3500 binary and 8000 additional ternary compounds from the ICSD and classified them into layered (quasi2D) and not layered materials. We restricted our search to stoichiometric and ordered systems that do not contain rare earth elements and have 50 of less atoms in the unit cell. From our calculations, we have identified 426 binary and 443 ternary quasi2D compounds. A full list of these materials can be found in supplementary information with their ICSD id.
ii.2 Calculating Piezoelectric Properties of Quasi2D Materials
In quasi2D structures considered in this work, the individual layers are held together by relatively weak vdW interactions. The standard exchangecorrelation functionals typically employed in density functional theory (DFT) calculations, including the calculations of elastic and piezoelectric properties, are known to fail to describe the vdW interactions. This is evident from relatively large errors in the outofplane lattice constants and the associated elastic properties.lebegue2013two To overcome this issue, we employed a vdWcorrected functional (optB86) as implemented in VASP (Vienna Abinitio Simulation Package) codekresse1993ab; kresse1996efficient to calculate the lattice parameters, elastic, and piezoelectric properties of quasi2D materials.PhysRevB.83.195131; klimevs2009chemical To evaluate the piezoelectric coefficient tensors, we utilize the VASP implementation of the density functional perturbation theory (DFPT)baroni2001phonons; baroni1987green; gonze1995adiabatic calculations. A relatively large plane wave cutoff energy of 540 eV is used for structural relaxation, calculation of elastic tensors, and piezoelectric coefficient tensors. A dense kpoint grid, defined by , where is number of atoms in the primitive cell and is the number of kpoints, is employed. In all our calculations, a very high tolerance of 10 eV for energy convergence is used, which is an important consideration for conducting DFPT calculations.de2015database For calculation of elastic tensors we use a finite difference method. Here, the full elastic tensor is calculated by conducting six finite distortions of the lattice and obtaining elastic constants () from the stressstrain relationship.le2002symmetry; PRB2010elastic
The importance of incorporating vdWcorrections is illustrated in Figure 3, where we notice significant improvement in predicting elastic constants, particularly and , with vdWcorrected functional (optB86) PhysRevB.83.195131; klimevs2009chemical compared to standard GGAPBE functional.perdew1996generalized The GGAPBE predicted elastic constants are sourced from Ref.jain2013commentary. A more detailed analysis of the data presented in Figure 3 reveals that the GGAPBE is still better in predicting inplane elastic coefficients and , but the error in reproducing and is a factor of 10 or larger. This is due the fact that these particular two elastic constants are directly related to the deformations of the spatial gaps between the layers and the failure of GGAPBE in reproducing relatively weak vdW interactions. The comparison of predicted properties with GGAPBE and vdWcorrected functional is limited only to the elastic constants (Figure 3); the lack of experimental data on piezoelectric properties of quasi2D materials prevented us from making similar comparisons for predicted piezoelectric properties. In the supplemental information Table S1, we have provided the experimental details (e.g. measurement techniques, measuring temperatures, etc.) for each compound shown in Figure 3. The comparison of calculated and experimental values for the piezoelectric modulus of a few commercially important piezoelectric materials (including AlN and PbTiO) are shown supplementary information. Predicted piezoelectric moduli are found to be in good agreement with experimental values.
ii.3 Workflow for identifying quasi2D piezoelectrics
A complete workflow we developed for identifying promising quasi2D piezoelectric materials is illustrated in Figure 4. The binary and ternary crystal structures from the ICSD database are first screened using the automated algorithm for identifying quasi2D structures. Then, we filter out all centrosymmetric structures based on the space group assigned in ICSD. Out of 11500 binary and ternary materials we find 869 layered systems, out of which 135 are identified as having noncentrosymmetric structures. Next, the noncentrosymmetric, structures are relaxed using the previously described firstprinciples calculations employing a vdWcorrected functional. As the piezoelectric materials need to have sizable band gaps for their properties to not be screened by the existence of free charge carriers we next employ a bandgap filter. As suggested in the previous works,gorai2016computational; marom2009describing for electronic structure calculation, we perform selfconsistent GGAPBE calculations on the vdWrelaxed structures using dense kpoint grids. Because of the known band gap error in DFT calculations we use a relatively generous band gap cutoff of 0.1 eV. Fifty additional materials with their band gap smaller than 0.1 eV are discarded as a result. Finite difference calculations are performed with vdWcorrected functional to obtain the elastic constants (C). Five materials with elastic tensors with negative eigenvalues are also discarded. According to Born stability criteria,born1940stability the elastically stable materials always have positive eigenvalues of stiffness matrix. This means that an elastically stable materials always have positive elastic energy for arbitrary homogeneous deformation by an infinitesimal strain.mouhat2014necessary With the remaining 80 candidates, we proceed our calculations by performing Density Functional Perturbation Theory (DFPT) baroni2001phonons; baroni1987green; gonze1995adiabatic calculations to assess their piezoelectric coefficient tensor, e. Out of the remaining 80 candidates, 38 materials contain transition metal elements. For these, we perform a limited search of the magnetic ground state by enumerating all magnetic configurations in the unit cell and calculating their total energies. The elastic and piezo calculations are then performed only on the lowest energy spin state. This step is necessary because of the strong dependence of the electronic structure and other properties on spin configuration discussed in more details in Ref. gorai2016thermoelectricity. The automated DFT calculations including initial file generation, calculating properties, data extraction and data handling are performed with the help of PyLada,pylada a Python framework for highthroughput firstprinciples calculations.
Iii Results and Discussions
iii.1 Promising Quasi2D Piezoelectric Materials
Based on the calculated piezoelectric modulus tensor, a number of candidate materials with relatively large components have emerged. They are shown in Figure 5 with the full list together with the corresponding , and values provided in the supplementary information. In addition, the top 20 most promising systems, based on the largest component, are listed in Table 1.
The piezoelectric modulus tensor (d) is a third rank tensor, and any isotropic averaging scheme will yield zero.li2000effective To rank these materials based on the merit of their piezoelectric response, we define as the largest element of the absolute matrix. Then the is plotted against the space group number in Figure 5. The two reference lines have been drawn for categorizing these candidates – one representing calculated of PbTiO and the other representing calculated of AlN (see supplementary information for the benchmark agains experiments).
Please note that PbTiO is also a ferroelectric and here we are only using the value for its piezoelectric response. More precisely, the PbTiO reference line corresponds to the bulk piezoelectric modulus corresponding to the singlecrystal singledomain samples.
[b]
Compound  SG  (eV)  (pC/N)  max  (C/m)  max  (C/m) 

InTe  9  0.7  351.7  2.6  2.6  
PbS^{†}  186  0.2  161.4  8.3  8.3  
GeTe  160  0.6  148.4  3.3  3.3  
CuVO  1  0.9  106.9  0.8  0.2  
SnO  31  1.6  67.1  1.1  1.1  
BiInO  33  2.8  56.1  4.7  4.7  
BiWO^{‡}  41  1.7  54.1  3.9  2.9  
NaIO  81  2.8  48.4  0.7  0.6  
NaN  12  1.4  40.7  0.3  0.1  
KSnF  143  3.0  40.5  0.3  0.3  
MoVO  35  0.8  40.4  2.9  2.9  
TlBrO  160  3.0  38.8  1.0  1.0  
NaSnN  186  1.1  36.6  0.6  0.6  
CsTe^{‡}  36  0.5  31.3  0.6  0.4  
BiMoO  61  1.7  28.7  1.6  1.3  
AgI  186  1.3  27.8  0.3  0.1  
SbFCl  79  1.5  24.4  0.2  0.1  
MgCl  115  4.9  23.5  0.2  0.2  
PbRbO  36  1.3  22.2  0.7  0.7  
BiGeO  9  2.3  21.1  3.0  3.0 

hypothetical structures

dynamically unstable
The materials shown in the left panel of Figure 5 can broadly be divided in three categories. The first category is comprised of quasi2D compounds with larger than the longitudinal piezoelectric modulus of PbTiO (=119 pC/N)jaffe2012a – the key end member of most commercial highstrain piezoelectrics. We found three materials (InTe, PbS, and GeTe) in this category. Among them PbS is not in its ground state rocksalt phase, but in the hypothetical distorted NiAs structure which has found its way into the ICSD.PhysRevB.84.045206 The other two compounds have previously been experimentally synthesized,chattopadhyay1987neutron; sutherland1976indium but their piezoelectric moduli have not been reported so far.
In the second category we group all compounds which have larger than the longitudinal piezoelectric modulus of AlN (=5.5 pC/N)dubois1999properties and lower than the longitudinal piezoelectric modulus of PbTiO. The majority (48 compounds) of the piezoelectric candidates from our study fall in this category revealing that overall the vdW bonded quasi2D systems indeed exhibit a propensity toward large piezoresponse. This group is composed from oxides and other chalcogenides such as CuVO, SnO, BiInO, MoVO, SnS, InSe, CsTe and other (total of 30); halides such as KSnF, AgI, MgCl, and PbI (total of 4); pnictides such as NaSnN, NaSnP, KSnAs, and NaN (total of 4). Also, a number of materials in this group (10) are the mixed anion systems, e.g., NaIO, SbFCl. Not surprisingly, the more ionic systems like oxides, halides and nitrides are more frequently found closer to the top of the range. Finally, two compounds CsTe and BiWO are found in calculations to be dynamically unstable, but both have been experimentally synthesized (likely hightemperature phases).chuntonov1982synthesis; wolfe1969crystal
The last group is composed of materials with the piezoelectric response lower than the longitudinal piezoelectric modulus of AlN. Though these candidates exhibit low piezoelectric response, they could still be useful as the calculated moduli are comparable with that of quartz (=2.27 pC/N).bottom1970measurement We found a total of 12 compounds which fall in this category. Examples include: WS, RhF, ZrCl, and GaInS. The moduli d, e and C with other informations such as band gap, space group of these compounds are provided in the supplementary information.
We also observe that the piezoelectric compounds in the quasi2D family of solids are clustered mainly in three specific space groups, i.e., space group no. 36 (Cmc2), 160 (R3m), and 186 (P6mc). This is mainly a reflection of the population bias, as these are the three most frequently occurring noncentrosymmetric space groups in the quasi2D family of crystals.
Table 1 shows that in 9 of the highresponse quasi2D piezoelectric compounds, the component appears as . Note that because of the freedom in choosing the inplane axes, and are virtually indistinguishable (see the discussion section). This component corresponds to the thickness shearing deformation where the material shears like a deck of cards in the inplane direction, with no change in the other dimension. Materials with large can be used in a variety of applications including: sensors, actuators, accelerometer, material testing structural health monitoring, nondestructive testing (NDT), and nondestructive evaluation (NDE).boivin2017torsional The components of appearing as usually coincide with the components of appearing as . The distribution of different components of (and ) appearing as (and ) are discussed in more detail in the next section.
Another quantity that could influence the piezoelectric response is the band gap of the material. In this work the band gaps are calculated at the DFT level and are also shown in Table 1 and supplementary information. Only about 1/3 of the studied materials are found to have DFT band gaps below 1 eV. Given the well known underestimation of band gaps in DFT based methods we do not think that materials with DFT band gaps larger than 1 eV would suffer from problems related to the existence of free charge carriers due to thermal fluctuations. However, for those with smaller gaps thermal fluctuation may cause sufficient number of free charge carriers, which may lower the polarization upon straining these materials despite having large piezoelectric moduli. Of course, provided that the real band gap is sufficiently close to the DFT one. For these materials, a more accurate assessment of the electronic structure might be needed before they are considered for applications.
In relation to the chemical composition and toxicity it is also important to not that currently, the most widely used piezoelectric material is lead zirconate titanate (PbZrTiO or PZT).shrout2007lead However, PZT causes significant environmental problems because of its high lead content.shrout2007lead Hence, significant efforts have been made to develop leadfree piezoelectric materials. shrout2007lead; takenaka2005current; baettig2005theoretical In our work, we have identified 44 candidates that do not contain any toxic elements including Pb. Out of 44 candidates, 33 of them have their piezoelectric modulus larger than AlN.
iii.2 Role of van der Waals Interactions
In order to understand the role of van der Waals interactions on the piezoelectric response of quasi2D materials, we analyze the relationships of and to the corresponding strain and stress components, respectively. A histogram showing the number of compounds with a given component appearing as is shown in Figure 6 (a). We observe that the most frequent is . This indicates that in the majority of quasi2D materials the largest piezoelectric response, as measured by the , is along the layer stacking direction. The reason for this behavior is that the relatively weak vdW interactions allow large charge redistribution in the layering direction upon straining the system.
On the other hand, the piezoelectric modulus tensor d relates relates stress to polarization and combines two types of effects: (1) amount of strain due to application of stress, and (2) amount of charge redistribution (polarization) due to resultant strain produced by the applied stress. A similar histogram of the components appearing as is shown in Figure 6(b). We divide the components into three groups depending on the deformation types (stress components) and the polarization direction. The schematics of the polarization directions and the associated stress components is shown in Figure 6(c) together with the components connecting the two. Every component of represents a separate piezoelectric operating mode. Schematics of all possible piezoelectric modes are provided in supplementary information.
Group I: Applied stress deforms van der Waals bonds and the measured polarization coincides with the direction of the deformations. The longitudinal mode () and shear modes ( and ) fall in this class. As evident from the histogram in Figure 6(b) these are the most frequently appearing components. In our considered materials, the modes and are indistinguishable because of the arbitrariness of the choice of axes ‘1’ and ‘2’, while the axis ‘3’ is fixed by layer stacking directions. In both of these modes, the same stress component () is responsible for the deformation, which implies shearing of the van der Waals gaps. On the other hand, in the , the applied stress axially deforms (stretches or compresses) van der Waals gaps. Hence, the large piezoelectric responses are achieved by deforming (axial or shear) the relatively soft van der Waals bonds. The barheights of and in the histogram of in Figure 6(b) are larger compared to . This is because the shearing resistance values () of quasi2D materials are lower compared to their axial resistance values () (refer to Figure 3(b)).
Group II: Applied stress deforms van der Waals bonds, but the measuring polarization directions are different from their deformation directions. The face shear modes ( and ) and the thicknessextension modes ( and ) fall in this class. The schematics of and are provided in the supplementary information. Usually, the direction of polarization is facilitated by the direction of deformation. In such modes, the deformation directions are different from their measured polarization directions. This is the reason behind their low occurrence in the histogram of though the van der Waals bond are deformed by the applied stress.
Group III: The van der Waals bond does not deform by the applied stress in these modes. This is why these modes do not appear frequently in the histogram of . This includes length or width extension modes ( or ) and shearing modes of type and .
The above discussion implies that the large piezoelectric response is always accompanied and caused by the stress that deforms the “soft” van der Waals gaps either through stretching, compression, or shearing. In addition, the analysis of the maximal piezoelectric moduli and the associated stress components provided guidance to experimentalists of how thin films should be grown to utilize large . It also describes what kind of mechanical actuation is necessary to achieve large piezoelectric response. This also helps in the design of new devices to take advantage of large . Finally, these findings can open up a wide variety of devices based on their operation modes or based on new materials for nonconventional modes. For example, materials with large faceshear () mode response will be attractive for for torsional applications, like novel gyroscopic sensors or highprecision torsional MEMS actuators.boivin2017torsional
iii.3 Axial Piezoelectric and Elastic Anisotropy
In addition to revealing new candidate piezoelectric materials and explaining the origin of strong piezoelectric response, we have also investigated the anisotropy in axial piezoelectric response in relation to the axial elastic anisotropy. We analyze how the response in the outofplane direction compares to the inplane responses. The inplane and outofplane directions in quasi2D materials are trivial to define and are illustrated in Figure 1. We define the axial anisotropy in both elastic and piezoelectric responses by the ratio of the outofplane component to the inplane component. Here, the inplane component is defined by the arithmetic mean of the ‘11’ and ‘22’components (invariant to the choice of the inplane axes), whereas the outofplane response is solely defined by the ‘33’component. Hence, the axial anisotropy of d, e, and C can be expressed as , , and respectively. If the axial anisotropy equals or nearly equals 1 then the material is considered to be isotropic in the corresponding quantity responses with respect to the inplane and outofplane directions. If the axial anisotropy is greater (or lower) than 1 the response of a material to axial deformation is dominated in outofplane (inplane) direction.
All possible correlations among the axial anisotropy of d, e, and C have been investigated. The correlation between axial anisotropy of d with respect to the axial anisotropy of e and C are shown in Figure 7 (a) and (b) respectively. From the comparative studies between Figure 7 (a) and (b), we see that the axial anisotropy of d is mainly dictated by the axial anisotropy of e not by the axial anisotropy of C. The plot between axial anisotropy of C and e is provided in the supplementary information.
From Figure 7 (a) and (b), we observe that most quasi2D materials have axial piezoelectric anisotropy parameter (both in d and e) 1 and axial elastic anisotropy parameter 1. This implies that the quasi2D piezoelectric materials are dominant in outofplane piezoelectric responses but elastically they are dominant inplane directions. This result corroborates the correlation between elastic softness and large piezoelectric response, i.e., the large piezoelectric responses are observed in elastically softer directions. We found quasi2D piezoelectric compounds such as CuVO and LiSbO are dominated by inplane piezoelectric response. Compounds such as BiWO, AlZnS, GeZnO, and GeTe are nearly isotropic in their axial elastic responses but highly anisotropic in axial piezoelectric responses.
Iv Conclusions
In conclusion, we performed a largescale computational (firstprinciples) assessment of the bulk piezoelectric properties of layered (quasi2D), vdW bonded materials. In our study we concentrate on the piezoelectric modulus as the measure of the piezoelectric response, which relates mechanical stress and electric polarization and depends on a combination of charge redistribution due to strain and the amount of strain produced by the stress. Overall, out of 135 noncentrosymmetric quasi2D binary and ternary structures from ICSD we have discovered 51 materials with piezoelectric response larger than that of AlN, a wellknown piezoelectric materials used in applications. Out of these 51 systems, we find three with the piezoelectric modulus even larger than that of PbTiO that has the piezoelectric modulus among the largest known. More importantly, 33 out of the 51 layered compounds do not contain any toxic elements including Pb. Our results also reveal that the large piezoelectric modulus in vdW systems is directly enabled by the vdW interactions between layers as in majority of compounds the large components of the piezoelectric modulus tensor couple to the stress components that imply deformations (both shear and axial) of the “soft” vdW bonds between layers. Our results suggest that quasi2D layered materials are a rich structural space for discovering new piezoelectric materials.
Acknowledgement The authors gratefully acknowledge the support of the National Science Foundation through Grant No. DMREF1534503. The calculations were performed using the highperformance computing facilities at National Renewable Energy Laboratory (NREL) and at Colorado School of Mines (the Golden Energy Computing Organization) .