Structural, electronic, vibration and elastic properties of the layered AgInP₂S₆ semiconducting crystal – DFT approach

Abstract

Detailed first principles calculations of the structural, electronic and vibrational properties of the AgInP₂S₆ crystal are reported. The energy band spectra of the mentioned material using DFT/GGA-D methodology with the PBE functional was calculated for the first time. Stability of the AgInP₂S₆ crystals in contrast to Cu-containing representatives of M1M2P₂X₆ materials family (M1, M2 – metal, X – chalcogen) has been explained in the framework of the second-order Jahn–Teller effect. The high covalence of the Ag–[P₂S₆] bonds and strong hybridization of the 4d- and 5s-orbitals of the Ag atoms are responsible for the stability of the considered crystal. The calculated vibrational properties were compared with the available experimental data derived from Raman scattering spectroscopy and their good agreement was demonstrated. The electronic and vibration properties within the framework of a group theory approach were studied. Also elastic properties of the AgInP₂S₆ crystal were modeled and analyzed for the first time.

---

1 Introduction

Recently, hexachalcogenohypodiphosphates have attracted great attention due to their unique physical and chemical properties. Layered crystals belonging to the M1M2P₂X₆ (M1, M2 – metal cations, X – chalcogen) family are promising materials for functional electronics, because many of them exhibit ferroelectric, ferro- or antiferromagnetic as well as piezoelectric properties.¹ They exhibit mixed electron–ionic conductivity^2–4 and promising optical^5,6 and thermoelectric properties.⁷ Therefore, these layered crystals are interesting from the fundamental point of view and also they are important materials for practical as well as for technological applications. In this context, it is important to explain the microscopic origin of the wide variety of above-mentioned properties and their dependence on morphological and componential peculiarities of particular representatives of the M1M2P₂X₆ crystal family.

Recently, an exploration of functional 2D materials, besides graphene, evolved rapidly thanks to the synthesis of the variety of new materials and their functionalities.^8–17 In the recent works^8–23 the main interest is primarily focused on the unique electronic and optical properties of the bulk metal phosphorus trichalcogenides. Especially it applies to the quaternary compounds belonging to the M_1/2^IM_1/2^IIIPX₃ type crystals (mainly M_1/2^I = Cu, Ag; M^III = Cr, V, Al, Ga, In, Bi, Sc, Er, and Tm; X = S and Se).²⁴ In addition, the nonlinear effects as second harmonic generation (SHG) and chiral electroluminescence appearing in thin atomic layers have expanded the functionalities and potential applications of the 2D materials. The mentioned effects disappear in the centrosymmetric bulk counterparts.^25–30 Also, the broken inversion symmetry in monolayer 2D materials has shown the viability of optical valley control.^31,32 Furthermore, the external modulation opens a path to explore the physical and technological applications of the functional 2D materials. One of the most important result that was obtained in recent years is the fact that the CuInP₂S₆ ferrielectric is a promising functional material with polarization switching not only in bulk samples but even in thin flakes of few tens nanometers thickness.^33,34

Consideration of materials with different chemical diversity and structural complexity gives the possibility to investigate their comprehensive properties. Recently, the creation of artificial materials (heterostructures) based on components with the same structural properties is very popular. The heterostructure constructed from ferrielectric CuInP₂S₆ and paraelectric In_4/3P₂S₆ crystals may be here mentioned.³⁵ Semiconducting AgInP₂S₆ crystal discussed in present work also have a similar structure. Thus, one can assume the possibility to create an artificial materials on its basis. However, application of the AgInP₂S₆ crystal in the future, first of all needs their physical properties investigations.

The AgInP₂S₆ crystal was chosen as material representing the M1M2P₂X₆ family. For the first time, the mentioned semiconductor was synthesized in 1987 and its crystalline structure was reported in the work of Ouili and co-workers.³⁶ Being isostructural with the paraelectric phase of the CuInP₂Se₆ crystal, the AgInP₂S₆ structure does not manifest the presence of the dipole ordering and the ferroelectric phase transition.⁶ Detailed microscopic investigations of an origin of the AgInP₂S₆ crystal stability can give also a valuable insight to explain an intriguing behavior of other representatives of the above-mentioned crystal family. Also, the stable atomic thin flacks offers a platform to manipulate the structure of an artificial magnetic or ferroelectric samples by their doping using relevant atoms. Obviously, a synthesis condition of such new functional materials could be elaborated only based on known and clearly described physical properties of the host crystal.

In this case, the present investigations are focused on physical properties modeling of the AgInP₂S₆ crystal using quantum chemical methods. To our knowledge, information about properties of the pure AgInP₂S₆ crystals are almost absent in the literature. All these facts motivated us to study the mentioned unexplored system. We present the quantum chemical calculation of the structural, electronic, vibrational and elastic properties of the AgInP₂S₆ crystals, their analysis, and comparison with experimental data.

The presented paper is organized as follows: in the next section, the structure of the AgInP₂S₆ crystal accompanied by symmetry considerations and method of physical properties calculations are described. After that, the paper continues with presenting the results of the structural, electronic, vibrational and elastic properties calculations accompanied by their analysis. The work is ended by conclusions with a short summary.

2 Crystal structure of the AgInP₂S₆ crystal and method of calculations

According to ref. 36, the AgInP₂S₆ crystallizes in rhombohedral structure with the P3̄1c space group (no. 163). The lattice parameters and atomic reduced coordinates are collected in the Table 1. The AgInP₂S₆ crystal has prominent layered structure with two nonequivalent packet layers in the unit cell. They may be superimposed by the 180 degrees rotation around the Oz crystallographic axis and are separated by the so-called van der Waals (vdW) gap. Each packet layer consists of the [P₂S₆] anionic complex and two metallic cations (Ag and In) located at the center of sulfur near-octahedral polyhedrons connected one with the other by edges (Fig. 1). Every Ag-centered polyhedron is connected with three In-centered polyhedrons (and vice versa) meeting the threefold rotational symmetry of the P3̄1c space group.

Table 1 Atomic positions in the hexagonal AgInP₂S₆ crystal (reduced coordinates)

Atoms	x	y	z	Wyckoff position	Site symm.
Ag	2/3	1/3	1/4		3.2
In	0.0	0.0	1/4		3.2
P	1/3	2/3	0.16376		3
S	0.32303	0.34647	0.11906	12i(x, y, z)	1
a [Å]	6.18200
b [Å]	6.18200
c [Å]	12.9570
V [Å^3]	428.837

---

Fig. 1 The unit cell of the AgInP₂S₆ crystal structure and topology of the Ag- and In-centered polyhedrons.

The Fig. 2 presents k-path involving a set of high-symmetry points in the hexagonal Brillouin zone (BZ) used to calculate the electronic and phonon properties of the AgInP₂S₆ crystal structure. Additionally, as an initial input for all performed calculations of the structural, electronic and vibrational properties of the mentioned crystal the structural data from ref. 36 were taken. The CASTEP, module³⁷ of the Materials Studio program package, was used to carry out planned calculations. They were performed applying the generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof functional³⁸ or with the CA-PZ local functional based on the Ceperley and Alder data³⁹ parameterized by Perdew and Zunger.⁴⁰ The AgInP₂S₆ crystal possesses layered structure characterized by the presence of so-called vdW gap. Therefore, the extended DFT methodology is needed to calculate electronic properties of the mentioned material. In this case, one may propose to use the DFT-D method taking into account dispersion interaction elaborated by Grimme.⁴¹ The ultra-soft pseudopotential⁴² was used to perform calculations for Ag – 4d¹⁰5s¹5p⁰, In – 4d¹⁰5s²5p¹, P – 3s²3p³, S – 3s²3p⁴ atomic configurations.

Fig. 2 The hexagonal Brillouin zone and k-path used for electron and phonon spectra calculations. Coordinates of the high-symmetry points are the following: Γ(0, 0, 0), , , , , and .

The plane-wave basis set cut-off was chosen to be equal to 600 eV. The Monkhorst–Pack k-points grid⁴³ sampling was set at 12 × 12 × 3 points for the Brillouin zone. The convergence tolerance parameters were as follows: energy 5 × 10⁻⁶ eV, force 0.01 eV Å⁻¹; stress 0.02 GPa; displacement 0.05 Å. The total energy convergence criterion was assumed to be fulfilled when the self-consistent field (SCF) tolerance reach the value 10⁻⁷ eV per atom.

3 Results of calculations

3.1. Structural properties and symmetry analysis

As it was mentioned above, symmetry of the AgInP₂S₆ crystal lattice is defined by the P3̄1c (D_3d²) space group. Symmetry operations (h_j – rotational part, α_j – nontrivial translation, ) suitable for mentioned structure are as follows:

Following the rules contained in the book of Kovalev⁴⁴ the possible non-trivial translations in relative coordinates and rotations denoted in the oblique coordinate system can be classified into six classes:

corresponding to the identity operation (E), two threefold rotations (C₃), three horizontal twofold screw rotations (C′₂), inversion (I), inversion rotations (S₃) and vertical glide planes (m_v), respectively. Table of characters of the irreducible representation (irrep) of little groups of the D_3d² group in high-symmetry points in the BZ and corresponding basis functions are presented in Table 2. The characters of the irreps presented in Table 2 can be used to establish the selection rules for the direct dipole optical transition for example in the center of BZ. It is well known that matrix element of the dipole transition is nonzero when the direct product of irreps Γ^*_α × Γ_β × Γ_γ contains the Γ₁ irrep. In the mentioned expressions the Γ_α corresponds to the irrep of the initial state, Γ_γ describes the symmetry of a final state and Γ_β is an irrep describing transformation rules of a radius-vector r⃑(x, y, z) which coincides with a dipole operator. One may assume that the initial state has Γ₆ symmetry, final state transforms according to Γ₁ and the polarization of incident light is in XOY plane (symmetry is Γ₆). Consequently, the direct product is as follows Γ^*₆ × Γ₆ × Γ₁ = Γ₁ + Γ₃ + Γ₅. As far as the product contains Γ₁ the above transition is allowed. Summary of all possible initial and final states and all possible polarizations is presented in Table 3. In particular, the direct Γ₆ → Γ₁ transition is allowed only for E ⊥ OZ polarization.

Table 2 Characters of irreps and the basis functions of little groups of the D_3d² space group in high-symmetry points of BZ. Here (x, y, z) and (J_x, J_y, J_z) are components of polar and axial vectors, respectively. Representations in curly brackets are combined due to the time inversion symmetry. Shaded cells represent the symmetry operations absent in corresponding little groups

---

Table 3 Summary of all possible initial and final states and all possible polarizations

Initial state	Allowed final states
	E ‖ OZ (Γ4)	E ⊥ OZ (Γ6)
Γ1	Γ4	Γ6
Γ2	Γ3	Γ5
Γ3	Γ2	Γ6
Γ4	Γ1	Γ5
Γ5	Γ6	Γ2, Γ4, Γ6
Γ6	Γ5	Γ1, Γ3, Γ5

---

As it was mentioned above, the AgInP₂S₆ crystal structure has two nonequivalent layers in the unit cell. In the quasiparticle spectra of such crystals the small splitting, known as Davydov's splitting, take place in dispersion curves.⁴⁵ Group theory allows to analyze these states basing on the knowledge of the single layer symmetry. In the considered structure, the symmetry group of the single layer is D₃¹ and it contains the following six operations:

Irreducible representations of the mentioned subgroup in the center of the BZ are collected in Table 4. Symmetry analysis allows to establish relations between irreducible representations of the D_3d² group and its subgroup D₃¹ giving the symmetry description of Davydov's splitting. Results of this relations are summarized in the following diagrams:

Γ^layer₁ → (Γ₁, Γ₂); Γ^layer₂ → (Γ₃, Γ₄); Γ^layer₃ → (Γ₅, Γ₆). (4)

Table 4 Characters of irreps of the symmetry group corresponding to a single layer of AgInP₂S₆

	E (g1)	C3 (g2, g3)	C′2 (g4, g5, g6)	f(r)
Γ^layer1 (A1)	1	1	1	1, x^2 + y^2, z^2
Γ^layer2 (A2)	1	1	−1	z, Jz
Γ^layer3 (E)	2	−1	0	(x, y), (xz, yz), (x^2 − y^2, xy), (Jx, Jy)

---

The above relations mean that all branches of the AgInP₂S₆ crystal band spectrum will be the symmetrically split doublets defining transformation rules by above-written couples of irreps.

Knowing the site symmetry of atomic positions (see Table 1) and thus the localization of the valence charge the energy band structure can be predicted by elementary energy bands concept.^46,47 Using the induction procedure^48,49 the representations of the crystal space group induced by irreps of the site symmetry groups for particular Wyckoff positions could be obtained. These induced representations are the elementary energy bands building the whole valence band. The structure of these elementary energy bands originating from atomic positions are present below:

position

A₁: Γ₁ + Γ₂ − M₁ + M₂ − L₁ − A₃ − K₁ + K₂ − H₁ + H₂;

A₂: Γ₃ + Γ₄ − M₃ + M₄ − L₁ − A₃ − K₁ + K₂ − H₁ + H₂;

E: Γ₅ + Γ₆ − M₁ + M₂ + M₃ + M₄ − 2L₁ − {A₁ ⊕ A₂} − 2K₃ − 2H₃;

position and position

A₁: Γ₁ + Γ₂ − M₁ + M₂ − L₁ − A₃ − K₃ − H₃;

A₂: Γ₃ + Γ₄ − M₃ + M₄ − L₁ − A₃ − K₃ − H₃;

E: Γ₅ + Γ₆ − M₁ + M₂ + M₃ + M₄ − 2L₁ − {A₁ ⊕ A₂} − K₁ + K₂ + K₃ − H₁ + H₂ + H₃;

position

A₁: Γ₁ + Γ₂ + Γ₃ + Γ₄ − M₁ + M₂ + M₃ + M₄ − 2L₁ − 2A₃ − 2K₃ − 2H₃;

E: Γ₅ + Γ₆ − M₁ + M₂ + M₃ + M₄ − 2L₁ − {A₁ ⊕ A₂} − K₁ + K₂ + K₃ − H₁ + H₂ + H₃.

For general position i(x, y, z) all corresponding irreps are present in each k-point:

A: Γ₁ + Γ₂ + Γ₃ + Γ₄ + 2Γ₅ + 2Γ₆ − 3(M₁ + M₂ + M₃ + M₄) − 6L₁ − 2({A₁ ⊕ A₂} + A₃) − 2(K₁ + K₂ + 2K₃) − 2(H₁ + H₂ + 2H₃).

The induced representations for position located in the middle of P–P bond are also listed above. Their structure completely coincides with those induced by irreps for the local position . These relations will be used to analyze the structure of the calculated energy band spectrum and to investigate the spatial valence charge density distribution depending on the symmetry of the structure.

The electronic properties calculations were performed for fully optimized AgInP₂S₆ crystal structure. It means that the atomic positions and lattice parameters were relaxed according to the total energy minimization. The optimization was performed using the GGA/PBE and GGA/PBE-D methodology. It turned out that the lattice parameters obtained by using the GGA/PBE methodology are somewhat greater than the experimental ones. The data obtained using the GGA/PBE-D method are better matched to the experimental results. These results are collected in Table 5. Good compatibility of experimental data with theoretical ones allows to conclude that the applied calculation methodology can be used to predict an electronic and vibrational properties of the AgInP₂S₆ crystal presented in the next sections.

Table 5 Experimental and calculated structural parameters obtained by GGA/PBE-D methodology for hexagonal AgInP₂S₆ crystal

Bond identification	Bond length
	Exp.^34 [Å]	GGA/PBE-D [Å]
Ag1–S1	2.7784(10)	2.74837
Ag1–In1	3.5692(70)	3.55313
Ag1–P1	3.7387(80)	3.72616
In1–P1	3.7387(80)	3.72616
P1–P1	2.2260(19)	2.24459
In1–In1	6.1820(11)	6.15420
In1–S1	2.6457(90)	2.67644
S1–S1	3.3715(15)	3.36013
Ag1–Ag1	7.3966(10)	7.41372
a = b	6.182(2)	6.154207
c	12.957(2)	13.013575
V, [Å^3]	428.8(4)	426.846

---

3.2. Electronic properties

The energy band structure and the partial density of states (pDOS) calculated by GGA/PBE-D methodology for the AgInP₂S₆ crystal with optimized geometry as mentioned above are presented in Fig. 3. Constructed diagrams allow to conclude that the investigated crystal is indirect band gap semiconductor with the gap edges localized at the H and Γ high symmetry points (see Fig. 3b). The calculated value of the indirect energy gap E^calc_gi ∼ 1.072 eV is significantly smaller than available experimental estimation of the optical gap which is equal to 2.43 eV.⁵⁰ In general, this is consistent with the DFT methodology which systematically underestimates the value of the energy band gap of semiconductors. However, such essential discrepancy needs to be explained. The experimental measurements⁵⁰ have been carried out for AgInP₂S₆ powder samples at room temperature by diffuse light reflection method taking into consideration only direct transitions contribution. According to presented theoretical results, the studied AgInP₂S₆ crystal is an indirect-gap material, moreover the lowest direct transition (Γ₅ → Γ₁) is forbidden by symmetry for each polarization of an incident light (see Table 3). Additionally, as it can be seen in Fig. 3b, the lowest energy level of the conduction band is formed by an energetically isolated couple of branches which give the fairly low contribution to reflection coefficient at room temperature.⁵⁰

Fig. 3 Energy band structure of the AgInP₂S₆ crystal and corresponding partial DOS (pDOS) calculated by GGA/PBE-D methodology: overall view (a) and the vicinity of the band gap (b).

The extra validation of reliability of our results can be found by comparison of calculated (on the basis of band structure Fig. 3) direct and indirect optical gap (see Fig. 4) and experimental results.⁵⁰ The Fig. 4 presents the calculated energy dependence of the light absorption coefficient α (see insertion), the (α × hν)² and (α × hν)^1/2. Based on the presented data and the simple parabolic band model the rough estimations for indirect and direct optical gaps are equal to E^ind_g ∼ 1.75 eV and E^dir_g = 1.91 eV, respectively. The indirect gap is defined by order of magnitude of the optical phonon energy in the H point (∼±0.01 eV). These values obviously are in qualitative and better quantitative agreement with results of.⁵⁰

Fig. 4 Calculated energy dependence of light absorption coefficient α (insertion) of the AgInP₂S₆ crystal and (α × hν)², and (α × hν)^1/2 curves in low energy region.

Additionally, in semiconducting alloys, the variations of components composition strongly affects the band gap. Ferreira et al. has shown that with temperature increasing of the PbSnTe the energy gap also increases.⁵¹ The bandgap variation with temperature measured for the CsSnI₃ crystals shows the same uncommon temperature dependence.⁵² In present work the energy gap was calculated at the 0 K temperature and it seems that the same uncommon temperature dependence is observed.

From the fundamental solid state point of view, the calculated dispersion of the band structure is principal to calculate the effective masses of charge carriers. In present work the effective mass of electrons (m^*_e) was evaluated from the curvature (energy derivatives) of the bottom of conduction band in Γ point of the BZ and the effective mass of holes (m^*_h) was calculated analysing the curvature of the top of valence band in H point of the BZ. The diagonal elements of the effective mass tensor for electrons and holes were calculated following the equation:

In consequence, the effective mass of electrons and holes was determined by fitting the conduction and valence band, respectively, to a parabolic function. The values of reduced effective masses of charge carriers in the AgInP₂S₆ crystal are as follows: ; ; ; . Using these data the ratio of m_‖/m_⊥ = 1.28 for electrons and ratio of m_‖/m_⊥ = 1.53 for holes were obtained. The relatively low values of reported electrons effective masses suggest that the mobility of charges in investigated AgInP₂S₆ single crystal should be relatively high. Dziaugys et al.⁴ measured that the conductivity of AgInP₂S₆ crystals is ionic (Ag⁺ ions) and is equal to 10⁻⁷ S m⁻¹ at room temperature. It is in agreement to the high value of the effective mass of holes because in ionic semiconductors the mobility of holes is very low.

In Fig. 5 the calculated projected density of states (prDOS) for [P₂S₆] anion complexes and cation atoms are presented. Analyzing Fig. 3 and 5 one may say that the valence band can be separated into nine subbands denoted by D1…D9 (see Fig. 5). The energy band structure can be analyzed in the framework of elementary energy bands concept⁴⁷ and, in consequence, the nature of the mentioned subbands may be explained.

Fig. 5 Projected densities of states (prDOS) for anion complexes [P₂S₆] (a) and cations (b) as well as crystal orbital Hamilton population (COHP)^53,54 of Ag–S (c) and In–S (d) orbitals overlapping in AgInP₂S₆ crystal calculated by GGA/PBE-D methodology.

Wave functions analysis of symmetry transformation properties in all high-symmetry points of the BZ defined for the AgInP₂S₆ crystal can designate all valence and set of the lowest conduction states. Corresponding sequences of irreps in high symmetry points Γ(0, 0, 0) and defined as the band structure extremes are as follows:

in Γ(0, 0, 0) point:

(Γ₁ Γ₂) (Γ₃ Γ₄ Γ₆ Γ₅ Γ₆ Γ₅ Γ₁ Γ₂) (Γ₅ Γ₆ Γ₅ Γ₆) (Γ₁ Γ₂) (Γ₃ Γ₄) (Γ₁ Γ₂) (Γ₆ Γ₅ Γ₆ Γ₅ Γ₂ Γ₁) (Γ₆ Γ₃ Γ₅ Γ₁ Γ₆ Γ₅ Γ₂ Γ₄ Γ₁ Γ₅ Γ₃ Γ₆ Γ₄ Γ₂ Γ₅ Γ₆ Γ₅ Γ₆) (Γ₆ Γ₅) ↓ (Γ₁ Γ₂) (Γ₃ Γ₄ Γ₆ Γ₄ Γ₁ Γ₃ Γ₅ Γ₂)…

in point:

(H₃) (H₃ H₃ H₃ H₁ H₂ H₃ H₃) (H₁ H₂ H₃ H₂ H₁ H₃) (H₃) (H₃) (H₁ H₂) (H₃ H₁ H₁ H₃ H₃ H₂ H₂) (H₃ H₁ H₂ H₃ H₃ H₁ H₃ H₁ H₃ H₂ H₁ H₂ H₃ H₃ H₂ H₃ H₃ H₃) (H₁ H₂ H₃) ↓ (H₁ H₂) (H₃ H₃ H₃ H₁ H₃ H₂)…

Also, the symmetry of the wave functions in all other high symmetry points of the Brillouin zone was determined (see Table 2). These band symmetry designators combined into sets correspond to subbands denoted by D1…D9 (see Fig. 5). As a result, the decomposition of the whole valence band into elementary energy bands was obtained:

D1: Γ₁ + Γ₂ − M₁ + M₂ − L₁ − A₃ − K₃ − H₃ → or ;

D2: Γ₃ + Γ₄ − M₃ + M₄ − L₁ − A₃ − K₃ − H₃ → or ;

Γ₁ + Γ₂ − M₁ + M₂ − L₁ − A₃ − K₃ − H₃ → or ;

Γ₅ + Γ₆ − M₁ + M₂ + M₃ + M₄ − 2L₁ − {A₁ ⊕ A₂} − 2K₃ − 2H₃ → ;

Γ₅ + Γ₆ − M₁ + M₂ + M₃ + M₄ − 2L₁ − {A₁ ⊕ A₂} − K₁ + K₂ + K₃ − H₁ + H₂ + H₃ → or or ;

D3: 2(Γ₅ + Γ₆) − 2(M₁ + M₂ + M₃ + M₄) − 2L₁ − 2{A₁ ⊕ A₂} − 2(K₁ + K₂ + K₃) − 2(H₁ + H₂ + H₃) → or or/and ;

D4: Γ₁ + Γ₂ − M₁ + M₂ − L₁ − A₃ − K₃ − H₃ → or ;

D5: Γ₃ + Γ₄ − M₃ + M₄ − L₁ − A₃ − K₃ − H₃ → or ;

D6: Γ₁ + Γ₂ − M₁ + M₂ − L₁ − A₃ − K₁ + K₂ − H₁ + H₂ → ;

D7: Γ₁ + Γ₂ − M₁ + M₂ − L₁ − A₃ − K₃ − H₃ → or ;

2(Γ₅ + Γ₆) − 2(M₁ + M₂ + M₃ + M₄) − 4L₁ − 2{A₁ ⊕ A₂} − 2(K₁ + K₂ + K₃) − 2(H₁ + H₂ + H₃) → or or/and ;

D8: 2(Γ₁ + Γ₂ + Γ₃ + Γ₄ + 2Γ₅ + 2Γ₆) − 6(M₁ + M₂ + M₃ + M₄) − 12L₁ − 4({A₁ ⊕ A₂} + A₃) − 4(K₁ + K₂ + 2K₃) − 4(H₁ + H₂ + 2H₃) → ;

Γ₅ + Γ₆ − M₁ + M₂ + M₃ + M₄ − 2L₁ − {A₁ ⊕ A₂} − K₁ + K₂ + K₃ − H₁ + H₂ + H₃ → or or ;

D9: Γ₅ + Γ₆ − M₁ + M₂ + M₃ + M₄ − 2L₁ − {A₁ ⊕ A₂} − K₁ + K₂ + K₃ − H₁ + H₂ + H₃ → or or .

Above presented underlined positions are those whose contribution is dominant, as it was shown in calculated prDOS (see Fig. 5).

According to electronic configuration of atoms, the total number of occupied energy states forming the valence band is equal to 70 in each k-point. Analyzing the distribution of energy bands and the densities of states (pDOS, prDOS) one may conclude that the width of the valence band is 16 eV and is formed by nine distinguished bundles of branches with different genesis (see Fig. 3 and 5).

The lowest energy states (from −16 eV to −15.5 eV), labeled as D1, are built up by 3s orbitals of P and S atoms (see Fig. 5). Next valence subband D2 located from −15.5 eV up to −14 eV has an essential contribution of localized 4d electrons of indium, which are rather chemically inert in this crystal. This subband has the contribution of 3s and 3p orbitals of P and S atoms as well. In the energy region from −13.5 eV to −12.5 eV, the subband D3 formed by overlapping of S 3s, P 3s, and S 3p states and admixture of In 4d orbitals is placed. Next three couples of energy bands combined into three subbands D4, D5 and D6 are located from ∼−11 to −6.5 eV in energy scale. These bands, except the D6, are formed by 3s and 3p orbitals of S and P atoms mixed in different proportions. The D6 subband also contains the strong contribution of In 5s states. These slightly split doublets are the prominent evidence of the layered structure of the investigated crystal and represent above-mentioned Davydov's splitting,⁴⁵ typical for crystals with translationally nonequivalent complexes in the unit cell.

The valence band branches, defined as D7 and located between −6.5 and −4 eV arises mainly from overlapping of the S 3p, P 3p, and S 3s orbitals with admixing of In 5p and 5s as well as Ag 4d and 5s states. In this energy range, the essential mixing of [P₂S₆] molecular orbitals with different cation orbitals (s and p for In, and s and d for Ag) take place. Both, the broadest D8 and rather narrow D9 subbands spread in energy from −4 to −1 eV and from −1 to 0 eV, respectively, are originated from 4d-states of Ag atoms and 3p-orbitals of sulfur with noticeable contribution of the P 3p and P 3s states (mainly in the D8 subband). Also, in the considered energy range (D8) one can find In 5p, Ag 5s, In 5s, and Ag 5p orbitals contribution (in order of magnitude decrease) which also indicate the complicated mixing of [P₂S₆]' molecular orbitals and electronic states of both metal cations. The topmost Ag 4d states are split into two separated ranges. The same situation is seen for the deep 4d states of indium. Such splitting appears due to the partial lifting of a degeneracy of atomic d-orbitals in the local crystal field. It is in agreement with the site symmetry of the In and Ag atoms located in the center of sulfur polyhedrons. It should be noted, that for ML_n-type molecules (M – transition metal, L – ligands) the d-orbitals split in accordance with the scheme (t_2g + e_g) for octahedral coordination or (e_g + t_2g) for tetragonal coordination.⁵⁵ However in considered structure with D_3d site symmetry the d-states combined into two higher isoenergetic pairs as (E′_g (d_xy, d_x²−y²) and E′′_g (d_xz, d_yz)) and one lower state with A_1g (d_z²) symmetry (see Fig. 6). The interaction of metallic cation states with molecular orbitals of the whole [P₂S₆] anion complexes leads to mixing of the s and d states of Ag, as well as the s and p states of In (whose d orbitals are lying much lower in energy). Mentioned interaction is weaker than in pure ML_n-molecules but involves almost all states and significantly raises the covalence of M–[P₂S₆] bonding.

Fig. 6 Calculated projected densities of states of d-orbitals of cations in the AgInP₂S₆ crystal.

The bottom of the conduction band is mainly formed by the In 5s and S 3p orbitals. As one can easily see from Fig. 5 and 6, the roles of d orbitals of two cations of Ag and In are completely distinct. It seems also important to mention that s and d orbitals of In atoms are significantly separated in energy (∼16 eV), while s and d orbitals of Ag atoms are separated by only ∼3 eV gap. Therefore, the d electrons of the In atoms form more chemically inert (no overlapping) and localized charge distribution, while self-overlapping Ag orbitals make them more labile and facilitate their delocalization as well as indirect exchange with neighboring molecular complexes. These features of cationic electrons redistribution in the AgInP₂S₆ crystal can be illustrated by maps of squared modulus of wave functions distribution in the plane (110).

It visualizes the charge distributions within atom layer packages (see Fig. 7). In addition, to confirm of the obtained results, the Mulliken charge population was calculated using GGA/PBE-D approach. Obtained results for s, p, d orbitals of In, Ag, P, S atoms are collected in Table 6. The presented results exhibit the cationic nature of the Ag and In atoms. It proved the localization of the In electrons and dispersed character of the Ag electrons.

Fig. 7 Maps of the valence energy states formed by D1…D9 subbands calculated by GGA/PBE-D methodology for AgInP₂S₆ crystal in Γ point. The visualization was done in the (110) plane.

Table 6 Mulliken charges in AgInP₂S₆ crystal calculated by the GGA/PBE-D method

	s	p	d	Total	Charge
Ag	0.50	0.30	9.98	10.78	0.22
In	0.94	1.27	10.0	12.21	0.79
P	1.49	2.96	—	4.45	0.55
S	1.89	4.47	—	6.35	−0.35

---

Discussing a stability of the AgInP₂S₆ seems right to compare the studied materials with a structure of other representatives of this materials family. In this case the CuInP₂S(Se)₆ crystal was chosen. It is well known that instability of Cu-contained crystals, successfully explained by the second-order Jahn–Teller effect (SOJT),⁵⁶ can destabilize the closed-shell ions such as Cu¹⁺, In³⁺ or Ag¹⁺.

Main factor, which affect the SOJT effect realization in these crystals is the covalency of bonds and possibility of hybridization of the s and d orbitals of cation.⁵⁷ Electronic configuration of the Cu atom, unlike the configuration of Ag or In atoms (discussed above) is characterised by the extremely close appearance of 3d- and 4s-orbitals, which are rather shallow on the energy scale. The d-orbitals of both considered atoms (Ag and Cu) can hybridize with molecular orbitals of [P₂S₆] anion complexes,⁵⁸ but while Ag 4d states spread in energy in rather wide range (∼−7 to 0 eV) and the 3d states of Cu are more localized in energy (∼−5 to 0 eV).

Thus, the s–d hybridization is more efficient in the AgInP₂S₆ crystal than in paraelectric phase of the Cu-containing isostructures. Situation changes at Cu displacement towards trigonal coordination on ferrielectric phase transition when symmetry lowering allows more effective mixing of the s and d states of Cu stabilizing the system and making distort crystalline structure energetically preferable (this results will be published thereafter). Also, we can conclude that distinguish in cations electronic structure lead to lowering the covalence in Cu–[P₂S(Se)₆] bonding in comparison with Ag–[P₂S₆] ones. It follows from the set of experimental facts. In particular, the difference between sum of ionic radii for Cu and S ions (r_ion(Cu) = 0.73 Å, r_ion(S) = 1.84 Å (ref. 59)) and calculated length of Cu–S bonds (∼2.54 Å) is significantly less than the one for Ag and S ions (r_ion(Ag) = 1.15 Å, r_ion(S) = 1.84 Å (ref. 59)), which is about of ∼0.2 Å, that also indicates the considerable covalent contribution in Ag–S bonding. Above estimations as well as experimental values of the energy gaps (E_g = 2.43 eV for AgInP₂S₆,⁵⁰ and E_g = 2.65 eV for CuInP₂S₆ (in paraelectric phase)⁶⁰) also indicates the lowering of covalence on Ag to Cu changing. This tendency also reflected in vibration spectrum of above-mentioned crystals. Namely, the comparison of the lowest optical vibration frequencies of CuInP₂S₆ and AgInP₂S₆ demonstrate decrease from 25 cm⁻¹ to 21 cm⁻¹ (ref. 61) and corresponding squared frequency lowers to about 20% while the atomic mass growths in about 40%. This also indicates the growth of Ag–S bonds rigidity in comparison to Cu–S ones. Thus, we can state that according to our numerical modeling the AgInP₂S₆ crystal should be stable in contrast to isostructural Cu-containing compounds according to known experimental facts. Concerning the role of the second cation (In³⁺) in considered structures it is worth to highlight that the main contribution in chemical bonding are made by its s- and p-orbitals, while relatively low lying d-states are almost inert. Being more delocalized the In 5s- and 5p-electrons strongly mix with molecular orbitals of the anion complexes and make their contribution in almost each subband of the energy band spectrum. The In³⁺ ion plays here the role of covalently (stable) bonded electron donors. It reflects in larger charge transfer from its s- and p-orbitals (5s^0.945p^1.27) in comparison to the Ag ion (5s^0.55p^0.3).

3.3. Vibrational and elastic properties

In the present work also the dynamics of the AgInP₂S₆ crystal lattice was investigated and reported. The phonon frequency of the AgInP₂S₆ crystal was measured by polarized Raman spectroscopy. The experiment was performed in Z(XX)Z geometry at room temperature with 514.5 nm Ar laser excitation lines. In the high-frequency region the spectra were satisfactorily fitted to Lorentzian profiles but in the low-frequency range, they were fitted by Voight profiles (Lorentzian convoluted with a Gaussian of appropriate width). Obtained results are depicted in the Fig. 8.

Fig. 8 Raman scattering spectrum of AgInP₂S₆ crystal measured experimentally in Z(XX)Z geometry at room temperature.

To determine the nature of the experimental Raman scattering spectra they were compared to the eigenfrequencies of phonon modes calculated by GGA/PBE-D methodology (see Table 7). The calculated frequencies are also depicted in Fig. 8 by red and pink arrows. The agreement between experimental and theoretically obtained data is very good. It proves that the GGA/PBE-D methodology optimize well the structure of the AgInP₂S₆ crystal.

Table 7 Vibration modes of AgInP₂S₆ lattice calculated by GGA/PBE-D methodology and their comparison to experimentally measured Raman frequencies

Mode N	Frequency (cm^−1)	Irrep.	Raman active	Raman Z(XX)Z [exp.]
1, 2	0.022501	Eu (x, y)	N
3	0.013105	A2u (z)	N
4	20.926421	A2g (Jz)	N
5	26.107787	A2u (z)	N
6, 7	35.979029	Eg (xz, yz)	Y	—
8	61.671083	A2g (Jz)	N
9, 10	75.514009	Eu (x, y)	N
11, 12	78.935074	Eg (xz, yz)	Y	65
13, 14	93.954836	Eu (x, y)	N
15, 16	100.369898	Eg (xz, yz)	Y	100
17	117.278994	A2u (z)	N
18	127.060009	A2g (Jz)	N
19, 20	135.102694	Eu (x, y)	N
21, 22	138.545071	Eg (xz, yz)	Y	133
23	179.280872	A2g (Jz)	N
24	179.576487	A2u (z)	N
25	180.786519	A1u (y(3x^2 − y^2))	N
26	182.083987	A1g (z^2)	Y	—
27, 28	184.492954	Eg (xz, yz)	Y	182
29, 30	185.042316	Eu (x, y)	N
31, 32	215.374008	Eu (x, y)	N
33, 34	215.601394	Eg (xz, yz)	Y	216
35	244.040140	A1u (y(3x^2 − y^2))	N
36	250.638375	A1g (z^2)	Y	—
37, 38	254.658638	Eu (x, y)	N
39	258.888355	Eg (xz, yz)	Y	252
41, 42	272.692033	Eu (x, y)	N
43, 44	275.028018	Eg (xz, yz)	Y	282
45	306.988799	A2u (z)	N
46	314.703934	A2g (Jz)	N
47	366.465808	A1g (z^2)	Y	374
48	367.398277	A1u (y(3x^2 − y^2))	N
49	446.468206	A2u (z)	N
50	448.103233	A2g (Jz)	N
51	538.402358	A1g (z^2)	Y	—
52	539.244632	A1u (y(3x^2 − y^2))	N
53, 54	547.301483	Eg (xz, yz)	Y	—
55, 56	548.380855	Eu (x, y)	N
57, 58	549.779358	Eg (xz, yz)	Y	560
59, 60	550.875484	Eu (x, y)	N

---

The bulk phonon dispersions along several symmetry lines together with the corresponding phonon density of states are presented in Fig. 9. Analyzing the corresponding atoms contribution into vibrational modes the frequencies may be divided into four regions: the lowest-energy region below 35 cm⁻¹, first middle region (35–175 cm⁻¹), second middle region (175–400 cm⁻¹) and the high-energy excitation region above 450 cm⁻¹. The low-frequency range mainly consists of Ag atom vibrations. The Raman modes in this region are largely affected by the vibration of metal atoms. No Raman peak could be resolved below 50 cm⁻¹ in AgInP₂S₆ while the peak around 25 cm⁻¹ is observed in the spectrum of the CuInP₂S₆.⁵⁹ The first middle region consists mostly of the In, S and P atom vibrations.

Fig. 9 Phonon dispersion curves and phonons partial DOS of calculated by GGA/PBE-D method for AgInP₂S₆ crystal along the symmetry directions of the BZ.

The second middle region includes only vibrations of the (P₂S₆)⁴⁻ anion complexes. It is characteristic for all hexachalcogenohypodiphosphates.⁶² In the high-frequency region, vibrations of phosphorus atoms forming P–S bonds are seen. The same results are observed for another chalcogenide crystals.^24,61

The theoretical results presented in Table 7 indicate that there is no peak allowed by symmetry located around 450 cm⁻¹ for the AgInP₂S₆ crystal but in our work (see Fig. 8) as well as in the work⁶ it has been detected. A presence of the mentioned peak can be an evidence of slight imperfection of investigated AgInP₂S₆ crystals. Wang and co-authors⁶ had attributed this peak to the stretch vibration of the P–P bond.

The elastic constants of AgInP₂S₆ crystal were also calculated using the GGA/PBE-D methodology. Complete set of these results is collected in Table 8. No reports on the elastic properties of the AgInP₂S₆ crystal have been found so far. Therefore, comparison of these results with experimentally or theoretical obtained data is not possible. The experimental data that can be used to validate the results presented in Table 8 are reported in the work of Samulionis and coauthors.⁶³ The ultrasonic sound velocity measured there for the CuInP₂S₆ crystal along the c-axis is equal to 4450 m s⁻¹ at 295 K. This velocity increases with decreasing temperature.

Table 8 The elastic modulus B and coefficients of elastic stiffness C_ij calculated by GGA/PBE-D method for AgInP₂S₆ crystal. All values are in GPa units except for the non-dimensional Poisson ratios ε_j (i, j = x, y, z)

Property	Value
C11	108.92750 ± 4.897
C33	39.95850 ± 9.673
C44	18.55615 ± 1.647
C12	35.81315 ± 2.981
C13	17.46785 ± 3.575
C15	5.15943 ± 0.354
Bulk modulus B	33.42266 ± 4.567
Young's modulus Ex, Ey	91.07583
Young's modulus Ez	35.74233
εxy = εyx	0.2965
εxz = εyz	0.3075
εzx = εzy	0.1207

---

In this case, the sound velocity for AgInP₂S₆ crystal was calculated using Christoffel equations⁶⁴ and obtained results were compared to data measured for the Cu(In,Cr)P₂S₆ crystal.^26,63

The calculated elastic stiffness constants C_iklj were used to construct the Christoffel matrix Γ_ij:

here n_i = cos(k⃑, a⃑_i) define directional cosine defining the propagation direction of the sound wave.

The system of linear eqn (6) can be solved when the corresponding secular equation meets dependency:

|Γ_ij − ρv²δ_ij| = 0 (7)

A solution of the secular eqn (7) allows to obtain three values of sound velocities corresponding to one longitudinal (LA) and two transverses (TA₁, TA₂) sound waves. The sound propagation velocity in (001) and (100) plane is presented in Fig. 10.

Fig. 10 The calculated sound velocity indicatrices in the plane of the atomic layer (001) (top panel) and in the perpendicular plane (100) (bottom panel).

The maximum and minimum values of sound velocity in the plane of the atomic layer (001) and in the perpendicular plane (100) are collected in Table 9. The obtained results are in good agreement with the value measured for the CuInP₂S₆ and CuCrP₂S₆ crystals.^63,65 It indirectly confirms that the elastic constants of the AgInP₂S₆ crystal were calculated properly using GGA/PBE-D methodology. The acoustic wave propagates in the (001) plane with the same velocity in all directions. It means that the AgInP₂S₆ crystal in (001) plane is isotropic. Perpendicularly to this plane the mentioned crystal is anisotropic. The ratio V_‖/V_⊥ = 1.65 can serve as a mechanical anisotropy magnitude indicator. Also the ratio of electron effective masses m_‖/m_⊥ = 1.28 confirm the anisotropy character of the AgInP₂S₆ crystal.

Table 9 Sound velocities in AgInP₂S₆ crystal calculated by GGA/PBE-D method using Christoffel' equations

	(100) plane		(001) plane
	Vmax, m s^−1	Vmin, m s^−1	Vmax, m s^−1	Vmin, m s^−1
LA	5430.2	3288.9	5437.5	5430.2
TA1	3204.4	2241.7	3204.4	3145.8
TA2	2429.7	2156.7	2223.4	2156.7

---

Despite absence in literature, the experimental results on sound velocity measurements for the AgInP₂S₆ crystal as mentioned above, can estimate these parameters by comparison to known values for the similar crystals as the CuInP₂S₆ and CuCrP₂S₆. According to ref. 63 sound velocity in [001] direction is equal to v₀ = 4460 ± 50 m s⁻¹ for CuInP₂S₆ and v₀ = 4160 ± 50 m s⁻¹ for CuCrP₂S₆, at low temperature (polar phase), and is equal to v₀ = 3875 ± 50 m s⁻¹ and v₀ = 3975 ± 50 m s⁻¹ at high temperature (centrosymmetric phase). Our results obtained for centrosymmetric AgInP₂S₆ crystal (v₀ = 3288 m s⁻¹) are close enough to above listed experimental values for related structures.

4 Conclusions

In the present work, the first principles calculations of physical properties of the chalcogenide AgInP₂S₆ crystal was performed. First of all the crystal structure was optimized applying the GGA/PBE-D methodology. Then the electronic, vibrational and elastic constants were also calculated using the mentioned DFT functional. Comparing the theoretically obtained structural data with the experimental results was shown that for layered chalcogenide crystals the DFT-D approach should be used. The detailed structural information defined by DFT-D approach give opportunity to explain chemical bonding character of the AgInP₂S₆ crystal.

The formation of subbands in the energy band spectrum of the layered crystal was analyzed. In this case the reason of the AgInP₂S₆ crystal stability in comparison with other (almost) isostructural representatives of the considered materials family (e.g. CuInP₂S₆ or CuInP₂Se₆ crystals) can be elucidated. Having the electronic structure calculated for the AgInP₂S₆ crystal one can assume that stability of this compound in comparison with Cu-containing structures, in particular, is based on a rather big covalence of the Ag–[P₂S₆] bonds and therefore more rigid Ag–S polyhedra (this is an important factor for realization or absence of the SOJT effect mechanism, which is often employed in analysis of phase transitions in this type of crystals).

For more detailed understanding physical properties of the AgInP₂S₆ crystal, the lattice dynamics was investigated within the GGA/PBE-D approach. The phonon dispersion curves and partial density of the AgInP₂S₆ crystal phonon states were calculated. Theoretically obtained results were compared with experimental data. Satisfactory agreement between calculated and experimental vibrational modes was achieved.

For the first time, the elastic modulus, coefficients of elastic stiffness of AgInP₂S₆ crystal and Poisson ratios were calculated in presented work. These data were used to calculate the sound velocity in the plane of the (001) atomic layer and in the perpendicular plane (100). Obtained results are in good agreement with the experimental results.

The presented work confirm that the electronic and structural properties of the layered chalcogenide structure can be characterized theoretically applying the DFT methodology at the GGA level of approximation adding the van der Waals correction. It seems to be reasonable for the layered structure characterized by van der Waals interlayer gap.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The calculations have been performed in Wroclaw Center for Networking and Supercomputing (Grant No. 171). The MATERIALS STUDIO package was used under POLAND COUNTRY-WIDE LICENSE. The financial support of National Science Centre, Poland, within the project No. 2016/21/N/ST3/00461 is also acknowledged.

Notes and references

1. A. Dziaugys, J. Banys, J. Macutkevic, R. Sobiestianskas and Y. Vysochanskii, Phys. Status Solidi A, 2010, 207, 1960.
2. V. Maisonneuve, V. B. Cajipe, A. Simon, R. VonderMuhl and J. Ravez, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, 10186.
3. Y. M. Vysochanskii, A. A. Molnar, M. I. Gurzan and V. B. Cajipe, Ferroelectrics, 2001, 257, 147.
4. A. Dziaugys, J. Banys, J. Macutkevic and Y. Vysochanskii, Integr. Ferroelectr., 2009, 103, 52.
5. I. P. Studenyak, O. A. Mykajlo, V. O. Stephanovich, M. I. Gurzan, Y. M. Vysochanskii and V. B. Cajipe, Phys. Status Solidi A, 2003, 198, 487.
6. X. Wang, K. Du, W. Liu, P. Hu, X. Lu, W. Xu, C. Kloc and Q. Xiong, Appl. Phys. Lett., 2016, 109, 123103.
7. Y. Vysochanskii, K. Glukhov, M. Maior, K. Fedyo, A. Kohutych, V. Betsa, I. Prits and M. Gurzan, Ferroelectrics, 2011, 418, 124.
8. A. K. Geim and I. V. Grigorieva, Nature, 2013, 499(7459), 419.
9. Z. Sun and H. Chang, ACS Nano, 2014, 8(5), 4133.
10. O. V. Yazyev and Y. P. Chen, Nat. Nanotechnol., 2014, 9(10), 755.
11. F. Xia, H. Wang, D. Xiao, M. Dubey and A. Ramasubramaniam, Nat. Photonics, 2014, 8(12), 899.
12. F. Bonaccorso, L. Colombo, G. Yu, M. Stoller, V. Tozzini, A. C. Ferrari, R. S. Ruoff and V. Pellegrini, Science, 2015, 347(6217), 1246501.
13. X. Luo, X. Lu, G. K. Koon, A. H. Castro Neto, B. Ozyilmaz, Q. Xiong and S. Y. Quek, Nano Lett., 2015, 15(6), 3931.
14. Y. Zhao, X. Luo, H. Li, J. Zhang, P. T. Araujo, C. K. Gan, J. Wu, H. Zhang, S. Y. Quek, M. S. Dresselhaus and Q. Xiong, Nano Lett., 2013, 13(3), 1007.
15. X. Lu, M. I. B. Utama, J. Lin, X. Luo, Y. Zhao, J. Zhang, S. T. Pantelides, W. Zhou, S. Y. Quek and Q. Xiong, Adv. Mater., 2015, 27(30), 4502.
16. X. Huang and H. Zhang, Sci. China Mater., 2015, 58(1), 5.
17. Z. Sun, D. Popa, T. Hasan, F. Torrisi, F. Wang, E. J. R. Kelleher, J. C. Travers, V. Nicolosi and A. C. Ferrari, Nano Res., 2010, 3(9), 653.
18. Y. Zhao, M. de la Mata, R. L. J. Qiu, J. Zhang, X. Wen, C. Magen, X. P. A. Gao, J. Arbiol and Q. Xiong, Nano Res., 2014, 7(9), 1243.
19. K. F. Mak, C. Lee, J. Hone, J. Shan and T. F. Heinz, Phys. Rev. Lett., 2010, 105(13), 136805.
20. X. Yin, Z. Ye, D. A. Chenet, Y. Ye, K. O'Brien, J. C. Hone and X. Zhang, Science, 2014, 344(6183), 488.
21. M. Gruning and C. Attaccalite, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89(8), 081102.
22. Y. Li, Y. Rao, K. F. Mak, Y. You, S. Wang, C. R. Dean and T. F. Heinz, Nano Lett., 2013, 13(7), 3329.
23. H. Yu, D. Talukdar, W. Xu, J. B. Khurgin and Q. Xiong, Nano Lett., 2015, 15(8), 5653.
24. K. Du, X. Wang, Y. Liu, P. Hu, M. I. B. Utama, C. K. Gan, Q. Xiong and C. Kloc, ACS Nano, 2016, 10(2), 1738.
25. X. Yin, Z. Ye, D. A. Chenet, Y. Ye, K. O'Brien, J. C. Hone and X. Zhang, Science, 2014, 344(6183), 488.
26. M. Gruning and C. Attaccalite, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89(8), 081102.
27. Y. Li, Y. Rao, K. F. Mak, Y. You, S. Wang, C. R. Dean and T. F. Heinz, Nano Lett., 2013, 13(7), 3329.
28. W. Kleemann, V. V. Shvartsman, P. Borisov, J. Banys and Y. M. Vysochanskii, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84(9), 094411.
29. Y. J. Zhang, T. Oka, R. Suzuki, J. T. Ye and Y. Iwasa, Science, 2014, 344(6185), 725.
30. E. Peeters, M. P. T. Christiaans, R. A. J. Janssen, H. F. M. Schoo, H. P. J. M. Dekkers and E. W. Meijer, J. Am. Chem. Soc., 1997, 119(41), 9909.
31. H. Zeng, J. Dai, W. Yao, D. Xiao and X. Cui, Nat. Nanotechnol., 2012, 7(8), 490.
32. K. F. Mak, K. He, J. Shan and T. F. Heinz, Nat. Nanotechnol., 2012, 7(8), 494.
33. M. Chyasnavichyus, M. A. Susner, A. V. Ievlev, E. A. Eliseev, S. V. Kalinin, N. Balke, A. N. Morozovska, M. A. McGuire and P. Maksymovych, Appl. Phys. Lett., 2016, 109, 171903.
34. A. Belianinov, Q. He, A. Dziaugys, P. Maksymovych, E. Eliseev, A. Borisevich, A. Morozovska, J. Banys, Y. Vysochanskii and S. V. Kalinin, Nano Lett., 2015, 15(6), 3808.
35. M. A. Susner, M. Chyasnavichyus, A. A. Puretzky, Q. He, B. S. Conner, Y. Ren, D. A. Cullen, P. Ganesh, D. Shin, H. Demir, J. W. McMurray, A. Y. Borisevich, P. Maksymovych and M. A. McGuire, ACS Nano, 2017, 11, 7060.
36. Z. Ouili, A. Leblanc and P. Colombet, J. Solid State Chem., 1987, 66(1), 86.
37. S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, K. Refson and M. C. Payne, Z. Kristallogr., 2005, 220, 567.
38. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
39. D. M. Ceperley and B. J. Alder, Phys. Rev. Lett., 1980, 45, 566.
40. J. P. Perdew and A. Zunger, Phys. Rev. B, 1981, 23, 5048.
41. S. J. Grimme, Comput. Chem., 2006, 27, 1787.
42. D. Vanderbilt, Phys. Rev. B: Condens. Matter Mater. Phys., 1990, 41, 8412.
43. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13(12), 5188–5192.
44. O. V. Kovalev, Representations of the Crystallographic Space Groups. Irreducible Representations, Induced Representations and Corepresentations, Brigham Young University, Provo, 1993.
45. A. Davydov, Theory of Absorption of Light by Molecular Crystals, Naukova Dumka, Kiev, 1951.
46. L. Michel and J. Zak, Europhys. Lett., 2000, 50(4), 519.
47. D. Bercha, K. Glukhov and M. Sznajder, Phys. Status Solidi A, 2007, 244(4), 1318.
48. J. Zak, Phys. Rev. Lett., 1980, 45, 1025.
49. J. Zak, Phys. Rev. B: Condens. Matter Mater. Phys., 1982, 25, 1344.
50. M. V. Potoriy, S. F. Motrya, I. P. Prits and P. M. Milyan, Nonorganic Chemistry, Sci. Herald of Volyn National Univ., 2009, vol. 24, p. 20.
51. S. Ferreira, E. Abramof, P. Motisuke, P. Rappl, H. Closs, A. Ueta, C. Boschetti and I. Bandeira, Braz. J. Phys., 1999, 29, 771.
52. C. Yu, Z. Chen, J. J. Wang, W. Pfenninger, N. Vockic, J. T. Kenney and K. Shum, J. Appl. Phys., 2011, 110, 063526.
53. R. Dronskowski and P. E. Blöchl, J. Phys. Chem., 1993, 97, 8617.
54. V. L. Deringer, A. L. Tchougreeff and R. Dronskowski, J. Phys. Chem. A, 2011, 115, 5461.
55. M. A. Halcrow, Chem. Soc. Rev., 2011, 40, 4119.
56. Y. Fagot-Revurat, X. Bourdon, F. Bertran, V. B. Cajipe and D. Malterre, J. Phys.: Condens. Matter, 2003, 15, 595.
57. J. K. Burdett and O. Eisentein, Inorg. Chem., 1992, 31, 1758.
58. K. Glukhov, K. Fedyo, J. Banys and Y. Vysochanskii, Int. J. Mol. Sci., 2012, 13, 14356.
59. R. D. Shannon, Acta Crystallogr., Sect. A: Found. Crystallogr., 1976, 32, 751.
60. I. P. Studenyak, V. V. Mitrovcij, G. S. Kovacs, M. I. Gurzan, O. A. Mykajlo, Y. M. Vysochanskii and V. B. Cajipe, Phys. Status Solidi B, 2003, 236, 678.
61. Y. M. Vysochanskii, V. A. Stephanovich, A. A. Molnar and V. B. Cajipe, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58(14), 9119.
62. Y. Mathey, R. Clement, J. P. Audiere, O. Poizat and C. Sourisseau, Solid State Ionics, 1983, 9(10), 459.
63. V. Samulionis, Ultragarsas, 2002, 4(45), 7.
64. L. Landau and E. M. Lifshitz, Theory of Elasticity, England: Butterworth Heinemann, Oxford, 3rd edn, 1986.
65. V. Samulionis, Y. Vysochanskii and V. Cajipe, Ferroelectrics, 2003, 295, 21.

---