van der Waals forces control ferroelectric–antiferroelectric ordering in CuInP₂S₆ and CuBiP₂Se₆ laminar materials

Abstract

We show how van der Waals (vdW) forces outcompete covalent and ionic forces to control ferroelectric ordering in CuInP₂S₆ nanoflakes as well as in CuInP₂S₆ and CuBiP₂Se₆ crystals. While the self-assembly of these 2D layered materials is clearly controlled by vdW effects, this result indicates that the internal layer structure is also similarly controlled. Using up to 14 first-principles computational methods, we predict that the bilayers of both materials should be antiferroelectric. However, antiferroelectric nanoflakes and bulk materials are shown to embody two fundamentally different types of inter-layer interactions, with vdW forces strongly favouring one and strongly disfavouring the other compared to ferroelectric ordering. Strong specific vdW interactions involving the Cu atoms control this effect. Thickness-dependent significant cancellation of these two large opposing vdW contributions results in a small net effect that interacts with weak ionic contributions to control ferroelectric ordering.

---

Introduction

Of recent recognition is the general principle that van der Waals forces can compete against covalent and ionic forces to control chemical bonding.¹ A critical application of this is the stabilization of gold surfaces and nanoparticles by Au(0)-thiyl bonding rather than the Au(I)-thiolate bonding believed for decades to dominate in these systems.² This effect is not restricted to nanotechnology but also leads to Fe(II)-thiyl versus Fe(II)-thiolate valence tautomerisation in cytochrome-P450 model compounds.³ A classic catalytic system is the very strong binding of benzene to copper surfaces,⁴ a system traditionally believed to involve Cu–C covalent bonding that instead now forms a paradigm for the understanding of strong van der Waals forces.⁵ The van der Waals force is well known as being critical to the control of self assembly,^6,7 hence strongly affecting many conceived devices made from 2D materials. Here, we show that strong specific van der Waals interactions associated with soft^8,9 copper ions go beyond this to control internal ferroelectric ordering in polar laminar materials containing sheets of soft sulfur or selenium atoms. These forces prove to be sufficiently strong to, in certain cases, overturn the natural electrostatic forces driving ferroelectricity to make antiferroelectric structures – a somewhat rare effect that is only recently becoming understood.¹⁰

Bulk materials of the form ABP₂X₆ often take on a laminar shape and are of considerable current interest as possible ultrathin materials for use as ferroelectrics in memory storage and other devices. Bulk systems investigated in this context include: CuInP₂S₆,^11–17 CuBiP₂Se₆,¹⁸ AgBiP₂Se₆,¹⁸ AgBiP₂S₆,¹⁸ CuCrP₂S₆,¹⁹ and CuVP₂S₆,²⁰ with very closely related ferromagnetic systems including Cr₂Ge₂Te₆.^21–23 Specifically, materials of interest of the form ABP₂X₆ usually have X as a chalcogenide S, Se, Te, etc., and A/B as either a monovalent/trivalent or divalent/divalent combination. Their laminar shape comes about as the materials contain [(PX₃)–(PX₃)]⁴⁻ anions that assemble with their P–P bonds aligned like a bed of nails that is then covered by two parallel “sheets” of X atoms. Between the sheets, the A and B metal atoms embed into nominally octahedral holes formed between the X atoms of different anions. The property of interest is that the metal atoms may sit at the centre of these octahedral holes, making paraelectric structures with hexagonal symmetry that have no dipole moment, or else they can move to a side of the octahedral hole to become close to just one of the two X planes, creating local polarization orthogonal to the layer. Throughout a layer, these dipole moments can align ferroelectrically, antiferroelectrically, or randomly, giving rise to observed intra-layer phase transitions.^18,19

If a single layer is internally ferroelectric, then adjacent layers can align either ferroelectrically or antiferroelectrically with respect to it, and it is this process with which we are primarily concerned herein. It is relevant to both bulk materials as well as to ultrathin materials, substances often called “nanoflakes”. Nanoflakes can be made by exfoliating bulk materials¹⁷ and possibly also by self-assembly of individually prepared layers. In both cases, layer stacking will be influenced by both inter-layer electrostatic forces as well as by van der Waals forces, making nanoflakes susceptible to external manipulation. What are the dominant features that control this stacking? If a layer can be thought of in the same way as a laboratory bar magnet, then the electrostatic force will favour ferroelectric alignment. Does the van der Waals force also show preference?

For ABP₂X₆, recent experiments reveal results with quite different characteristics: CuInP₂S₆ bulk^11–13 and nanoflakes¹⁷ are ferroelectric whereas the closely related material CuBiP₂Se₆ forms an antiferroelectric bulk polymorph.¹⁸ In this work, we use first-principles calculations utilizing density-functional theory (DFT) to show how these features arise.

Methods

All calculations are performed using VASP 5.4.1,²⁴ where the valence electrons are separated from the core by use of projector-augmented wave pseudopotentials (PAW);²⁵ the used “POTCAR” files are dated: S, P, Se-06Sep2000, Cu-22Jun2005, and In, Bi-08Apr2002; the number of electrons included per atom are S, Se-6, P-5, Cu-11, In-3, and Bi-5. Test calculations performed including semicore electrons to the calculations resulted in changes to critical energy differences of at most 1 meV. The energy cut-off for the plane-wave basis functions was set at 500 eV for bulk solids but otherwise kept at the default values obtained using “PREC = HIGH”. The energy tolerance for the electronic structure determinations was set at 10⁻⁷ eV to ensure accuracy. We use a k-space grids of size: 4 × 4 × 1 for bilayers, 4 × 4 × 2 for two-layer solids, and 4 × 4 × 1 for six-layer solids. Test calculations involving expansion of these grids to double those sizes were shown to affect calculated energy differences by less than 1 meV. VASP performs periodic imaging in all three dimensions and so to model bilayers we introduced a vacuum region of ca. 15 Å between the periodic images in the direction normal to the bilayer. To minimize interactions between these periodic images we applied dipolar correction as implemented in VASP (by setting the “IDIPOL” to “TRUE”) in early calculations. However, the effect was found to be negligible (<1 meV) and so this sometimes troublesome procedure was no longer applied. As a second check on the calculation of ferroelectric interaction energies, calculations were performed on systems expanded three times in the z direction, with average energies also varying by less than 1 meV. In the related system CuInP₂Se₆, small changes associated with these calculation effects have been shown to have influence on ferroelectric structure.²⁶

Geometry optimizations were performed for all structures, terminating when the forces on all atoms fell below 0.01 eV Å⁻¹. For nanoflakes, in-plane hexagonal unit-cells of dimension 6.5532 Å for CuBiP₂Se₆ (ref. 18) and 6.088 Å for CuInP₂S₆,¹² are used, whereas all lattice vectors are fully optimized for bulk solids. The starting structures for geometry optimizations were based on the properties of adjacent layers observed in these materials, with test calculations performed for alternate structures using the PBE-D3 method not realizing lower-energy alternatives. Hence, whenever possible, hexagonal symmetry and/or inversion symmetry was enforced in the calculations, with molecular alignments chosen so that the symmetry of the structure paralleled symmetry elements in the plane-wave basis set. To further enhance numerical stability, the “ISYM = 2” flag was used to enforce density-matrix symmetry.

Many computational methods considered involve the raw PBE density functional plus this combined with various vdW corrections: Grimme's D2 empirical correction²⁷ (PBE-D2), Grimme's D3 empirical correction²⁸ in its original form without Becke–Johnson damping with PBE (PBE-D3) and revPBE²⁹ (revPBE-D3), the exchange-hole based correction of Steinmann and Corminboeuf³⁰ (PBE-dDsC), the Tkatchenko–Scheffler method³¹ (PBE-TS), that with self-consistent screening (SCS)³² (PBE-SCSTS), and that extended to make Tkatchenko's many-body dispersion method^33,34 (MBD@rsSCS) (VASP flag “IVDW = 202”). The other computational methods are all based around the vdW density-functional approach of Dion et al.³⁵ This vdW correction is applied with the revPBE density functional³⁶ (revPBEvdW), the optPBE density functional³⁶ (optPBEvdW) and the optB88 density functional³⁶ (optB88-vdW). Also, in a form modified by Lee et al., it is combined³⁷ with the BP86 density functional³⁸ (vdW-DF2). Finally, we also consider the raw HSE06 hybrid density functional³⁹ as well as that combined with D3.⁴⁰

All optimized coordinates are reported in ESI.† All monolayer and crystal energies are reported per cell per layer (i.e., per 10 atoms), whilst all bilayer energies are reported per cell (i.e., 20 atoms).

Results and discussion

A. Ferroelectric order in CuInP₂S₆ crystals and nanoflakes

Crystal structures of the bulk phase of CuInP₂S₆ depict internally ferroelectrically ordered layers that ferroelectrically self-assemble with each other, producing ferroelectric crystals that can be exfoliated to produce ferroelectric nanoflakes.^11–13,17 Its observed crystal structure is depicted in Fig. 1. Nominally hexagonal layers stack above each other in such a way that overall hexagonal symmetry is lost, with six symmetrically equivalent layers stacking vertically to form a unit that is periodic in the direction orthogonal to the layer plane. These six layers are named “a” to “f” in Fig. 1. While the In atoms remain close to the centres of their octahedral holes, the Cu atoms move to take on nearly trigonal-planar arrangements within a single S sheet.

Fig. 1 PBE-D3 optimized structures and relative energies for CuInP₂S₆ based on the observed ferroelectric structure of the bulk material: (b) this bulk material with layers named [abcdef] and an antiferroelectric structure made from it by pushing the Cu and In atoms to the opposite side of their octahedral holes in alternate layers to make structure [a′bc′de′f], ΔE per cell per layer; and (a) bilayers extracted from these solids, ΔE per cell. All six layers in the bulk ferroelectric structure are symmetrically equivalent to each other; similarly, the bilayer structures made by selecting any two adjacent layers in the bulk structure are all symmetrically equivalent, as are any (e.g., [b′c′]) made from them by moving both Cu and In atoms to the opposite sides. Relative energies per bilayer are indicated as calculated by PBE-D3, with, in parenthesis, analogous energies calculated by PBE. Top and sectioned side views are shown for each structure: P-orange, Se-yellow, Cu-green, Bi-purple.

No antiferroelectric arrangement of CuInP₂S₆ has been observed, so we consider a possible form made simply by moving the Cu atoms from one S sheet in a layer to the other, allowing the In and other atoms to relax accordingly. In Fig. 1, we apply this transformation to alternate layers, converting the original ferroelectric crystal with a structure depicted as [abcdef] to the antiferroelectric structure [a′bc′de′f]. The same procedure could also have been applied to make an alternate antiferroelectric structure [ab′cd′ef′], but this is symmetrically equivalent and hence not shown. All optimized antiferroelectric structures embody inversion symmetry linking adjacent layers, so these structures have no net dipolar polarization. By sliding adjacent layers with respect to each other, many possible variants of the shown antiferroelectric structure are conceivable. The energy differences between such structures are subtle and will be considered elsewhere.⁴¹ Here we focus on a key feature of all antiferroelectric structures: as Fig. 1 shows, two very different types of inter-layer interactions are apparent depending on whether or not the Cu atoms face each other across the interface.

To highlight this effect, we consider possible bilayers derived from the (observed) ferroelectric and (envisaged) antiferroelectric bulk structures of CuInP₂S₆. The ferroelectric bilayer is named “F” and comprised of the layers [bc] extracted from the crystal. Extracting any other pair of layers, or interchanging the layer order to make [cb], produces symmetrically equivalent bilayers, so only one unique bilayer can be produced from the crystal in this manner. However, two unique antiferroelectric structures can be conceived, with one, named “A_ii”, made by extracting layers [a′c], and the other, named “A_oo”, made by extracting layers [bc′] instead. The A_ii bilayers are so named as they have both Cu atoms adjacent to each other on the “inside” of the bilayer, whereas the A_oo bilayers have them on the “outside”, as far away from each other as is possible, see Fig. 1.

Fig. 1 indicates the energies ΔE of the different antiferroelectric bulk structures per cell per layer and bilayer structures per cell relative to that of the ferroelectric structures, obtained using two representative computational methods, (i) the PBE generalized-gradient approximation density functional,⁴² and (ii) this combined with the D3 dispersion correction in its original form (PBE-D3).²⁸ The critical properties of these computational methods of relevance is that they similarly include covalent and ionic effects but only PBE-D3 embodies a realistic description of the van der Waals attraction term.

Both the PBE and PBE-D3 computational approaches predict that the observed ferroelectric bulk polymorph is more stable than its envisaged antiferroelectric variant by ca. 10 meV per cell per layer, a small but significant quantity. That the two approaches predict very similar outcomes is naively suggestive that electrostatic forces rather than van der Waals forces are most important in determining the electronic structure of laminal materials.

Alternatively, results for the bilayers presents a completely different picture of the nature of the intrinsic processes involved in self-assembly. PBE predicts that the antiferroelectric A_ii and A_oo bilayers are each about 6 meV less stable per cell per layer than the ferroelectric bilayer F, with all structures barely stable to spontaneous layer separation. As it predicts the bulk ferroelectric phase to be more stable than alternating A_oo and A_ii bilayer-type arrangements by 9 meV per cell per layer, a long-range favourable ferroelectric electrostatic interaction of 3 meV can be deduced. However, PBE-D3 predicts that A_ii is more stable than F by 92 meV per cell whereas A_oo is less stable by 88 meV. Hence the van der Waals force is perceived as being very sensitive to the details of inter-layer arrangements, with the overall bulk structure emerging in part as a result of the way these contributions act to cancel each other out. Improved treatments of the electrostatic interactions²⁶ will not affect this analysis.

Next we consider how robust this bilayer result is likely to be considering possible variations and improvements in the computational methods used. Table 1 shows results for the A_ii bilayer energy with respect to the F bilayer energy obtained using the above two computational methods plus 12 others. In one case, the PBE generalize-gradient functional is replaced with the HSE06 (ref. 39) hybrid functional, providing a much enhanced description of electron exchange. The small PBE preference of 5 meV per cell for the ferroelectric bilayer compared to A_ii is reversed by a 15 meV preference for A_ii, but adding D3 provides an extra 44 meV preference. This change is only half that (97 meV) predicted when D3 is added to PBE, but the general scenario remains unchanged. The 10 other methods utilized in Table 1 all involve variations of the PBE density functional combined with different dispersion approaches. In terms of magnitude of the effect, these approaches show quite a range with the total preference for A_ii ranging from 27 to 105 meV per cell – a significant range but again the qualitative effect remains preserved. Understanding these method differences is a topic for ongoing research, but a noteworthy results is the net stabilization of 81 meV predicted by revPBE-D3,²⁹ a method recently highlighted as one of the most robust methods of its type across all types of chemical scenarios.⁴³ Another is the value of 43 meV predicted by the many-body dispersion^33,34 (MBD) method, an advanced treatment of dispersion that, to treat conductors better, does not assume pairwise additivity;⁴⁴ it is also known to be widely robust to chemical scenario.⁵

Table 1 Calculated bilayer structural energies differences ΔE_Aii = E_Aii − E_F per cell between the lowest-energy antiferroelectric arrangement A_oo and the ferroelectric structure F, in meV

Method	CuInP2S6	CuBiP2Se6
PBE	5	−1
PBE-D2	−42	−42
PBE-D3	−92	−86
revPBE-D3	−81	−63
PBE-DFTdDsC	−33	−36
PBE-MBDrsSCS	−43	−45
PBE-SCSTS	−100	−76
PBE-TS	−83	−87
revPBEvdW	−22	−26
optB88vdW	−40	−40
optPBEvdW	−30	−31
vdWDF2	−26	−23
HSE06	−15	−5
HSE06-D3	−59	−64

---

Lastly, we consider the robustness of the main result considering subtle effects not normally considered in energetics analyses of materials. Calculations of the zero-point energy difference between the two structures at the PBE-D3 level indicated a change of less than 1 meV. Also, calculations of the room-temperature Helmholtz free energy change⁴⁵ at the PBE-D3 level indicate corrections of just −3 meV per cell. Therefore the best-possible description from first principles calculations available at the moment is that bilayers of CuInP₂S₆ should be antiferroelectric.

To demonstrate that this prediction is consistent with the observation of ferroelectric nanoflakes of CuInP₂S₆, we consider what happens when a third layer is added to an existing A_ii bilayer, say by spontaneous self-assembly. As shown in Fig. 2, a third layer can be added in one of two ways, making structure A_iiF, in which the new layer adds ferroelectrically to its neighbour, as well as structure A_iiA_oo, in which it adds antiferroelectrically. PBE-D3 calculations indicate that A_iiF is more stable than A_iiA_oo by 88 meV, paralleling the previous result found for bilayer formation. Making a tetralayer system from this trilayer then proceeds analogously, with PBE-D3 predicting that A_iiFF is 87 meV per cell more stable than A_iiFA_oo. Taken together, these results suggest that the energies of antiferroelectric multi-layer systems can be usefully described simply in terms of sums of bilayer interactions. The key prediction here is that nanoflakes grown under such conditions favouring sequential layer deposition increasingly develop more and more ferroelectric character as the number of layers increases.

Fig. 2 Multilayer PBE-D3 optimized structures and relative energies ΔE per cell, made based on the observed crystal structure of ferroelectric CuInP₂S₆, for (a) trilayers and (b) tetralayers, labelled and cross-sectioned as in Fig. 1. P-brown, S-yellow, Cu-green, In-purple.

However, if instead two A_ii bilayers self-assemble to make a tetralayer, then the PBE-D3 calculations predict the resultant structure, A_iiA_ooA_ii, to be 8 meV per cell more stable than A_iiFF, making the antiferroelectric structure the thermodynamically most stable one. This result is also easily understood in terms of pairwise inter-layer interactions. In Fig. 1 and 2, all results for bilayers, trilayers, and tetralayers were obtained constraining the in-plane lattice vectors to reference values. This minimizes the variables changing between calculations, allowing focus to be placed on key features such as the change in the van der Waals interactions between structures. Alternatively, the shown bulk structures are at fully relaxed lattices obtained individually for each computational method used, allowing focus to be placed on the small energy differences that in the end determine which structure is observed and which is not. Calculations of the bulk materials at the reference lattice vectors indicate energy changes of only 5 meV, however, so lattice relaxation does not have a profound effect. Also, the energies of the ferroelectric crystals at these reference geometries are 26 meV more stable compared to the antiferroelectric structures than one would expect simply by adding bilayer energies, an effect arising from long-range attractive dipole–dipole interactions present in ferroelectric structures. Proper treatment of these effects,²⁶ plus use of the most accurate computational methods possible, will be critical to making robust predictions for properties of specific nanoflakes. For example, the calculated 8 meV energy difference between A_iiFF and A_iiA_ooA_ii tetralayers is too small and too susceptible to method variations for this to be regarded as a robust prediction of tetralayer properties.

Nevertheless, robust predictions from the calculations are that CuInP₂S₆ should be antiferroelectric, that under kinetic control subsequent layers will add ferroelectrically to an initial antiferroelectric bilayer, and that ferroelectric CuInP₂S₆ structures are intrinsically more stable than antiferroelectric ones for large nanoflakes. Given this, then at some point during sequential layer growth, the point will come at which kinetically produced structures transform into their thermodynamically more stable counterparts. In general, the issue of kinetic verses thermodynamics control in molecular self-assembled monolayer production is of considerable currently of interest, with new experimental and computational methods recently developed⁴⁶ and reviewed.¹ In the experimentally observed nanoflakes,¹⁷ layers were exfoliated from ferroelectric bulk material. Naively, one would expect such a high-energy exfoliation process to provide ample opportunity for the most thermodynamically stable products to be formed.

Through consideration of the details of the structures of the observed¹⁷ nanoflakes, the calculated results provide a basis for understanding. These nanoflakes were ca. 1 μm in size with ca. half their area being 3 layers high, a quarter 4 layers high, and the remaining quarter 2 layers high. Trilayer structures, and indeed all structures containing an odd number of internally ferroelectrically ordered layers, must have net dipole polarization, the issue being whether the net polarization corresponds to that of just a single layer or else multiples of this up to the full ferroelectric order of the nanoflake. When only a small amount of a fourth layer is added to a trilayer system, the energetics indicate that it can only add ferroelectrically. The fact that antiferroelectric structures come in two fundamental types and that these must be combined in multi-layer systems explains the experimental observations despite predictions that bilayers should be antiferroelectric.

B. Antiferroelectric order in CuBiP₂Se₆ crystal

This material is observed at low temperature to comprise an antiferroelectric combination of ferroelectric monolayers, see Fig. 3.¹⁸ It differs from CuInP₂S₆ primarily through the horizontal alignment of the layers, with the alignment in CuBiP₂Se₆ crystals such that the hexagonal symmetry of each layer is maintained throughout the crystal. Again, six layers are present in the smallest unit cell that shows periodicity in the direction normal to the layer plane, but for CuBiP₂Se₆ these consists of layer pairs that are replicated each three times. The layers in the observed crystal are labelled [ab′cd′ef′] using the same notation as that already applied in Fig. 1 and 2. By simply transferring Cu atoms within octahedral holes, a ferroelectric alternative is envisaged as [abcdef] shown in Fig. 3; this is symmetrically equivalent to the other structure makeable in this fashion, [a′b′c′d′e′f]. From the observed antiferroelectric structure, two fundamentally different bilayer types can be extracted as A_ii (e.g., [b′c] shown in Fig. 3) and as A_oo (e.g., [ab′] shown in Fig. 3). In contrast to bilayers built based on CuInP₂S₆, [bc′] is inequivalent to [ab′] while [a′b] is inequivalent to [b′c]; hence these other possibilities are also shown in Fig. 3 and labelled and , respectively. Also, two different types of ferroelectric bilayers can be simply constructed from the observed antiferroelectric crystal structure, with [ab] and [bc] shown as examples in Fig. 3 and named F and F′, respectively. While PBE-D3 calculations indicate that these ferroelectric variants have very similar energies, large differences are again predicted between antiferroelectric forms, and we focus on this effect.

Fig. 3 PBE-D3 optimized structures and relative energies for CuBiP₂Se₆ based on the observed antiferroelectric structure of the bulk material: (c) this bulk material with layers named [ab′cd′ef′] and a ferroelectric structure made from it by pushing the Cu and Bi atoms to the opposite side of their octahedral holes in alternate layers to make structure [abcdef], ΔE per cell per layer; (a) bilayers extracted from the top two layers of these solids, and (b) bilayers extracted from their second and third layers, with ΔE per cell. The shown ferroelectric bilayer structures [ab] and [bc] are symmetrically equivalent to the unshown structures [a′b′] and [b′c′], respectively, while for the bulk solid the unshown structure [a′b′c′d′e′f′] is equivalent to [abcdef]; also, in the bulk solids, the later layer pairs are all symmetrically equivalent to the first pair. Relative energies per bilayer are indicated as calculated by PBE-D3, with, in parenthesis, analogous energies calculated by PBE. Top and sectioned side views are shown for each structure: P-orange, Se-yellow, Cu-green, Bi-purple.

The results shown in Fig. 3 for PBE and PBE-D3, supported by those obtained using the other methods shown in Table 1, indicate that variations in covalent and ionic bonding poorly discriminate between bilayer structures, whereas dispersion again has a profound effect. As for CuInP₂S₆, the attractive arrangement of dipoles in a ferroelectric bulk material stabilizes this form by ca. 3 meV per cell per layer, a small effect compared to the magnitude of the differential van der Waals forces. However, the significant difference to CuInP₂S₆ is that for it the stabilization and destabilization of adjacent bilayer interfaces present in the bulk material cancel, the (PBE-D3) stabilization of A_ii by 86 meV per cell overpowers the destabilization of A_oo by 61 meV per cell to provide the driving force for the formation of the observed antiferroelectric crystal of 13 meV per cell per layer.

C. The van der Waals forces that control ferroelectric ordering

Table 2 explores the origins of the differential van der Waals bonding between the F, A_ii, and A_oo structures, listing the van der Waals contributions coming from the D3 correction to CuBiP₂Se₆ and CuInP₂S₆ bilayer energies. The D3 correction takes on a “pairwise additive” form involving summation of London dispersion energies between every pair of atoms in the system.^5,28,44 Hence, in principle it is possible to decompose this contribution into terms arising from constituent atoms, and here we do this considering contributions from Cu to Cu interactions, Cu to non-Cu interactions, and non-Cu to non-Cu interactions using:

Table 2 Analysis of the total D3-calculated dispersion energy E^D3, in eV per cell, for F, A_ii, and A_oo ABP₂X₆ bilayers, partitioned into contributions involving Cu atoms (Cu) and those not involving Cu atoms ; dispersion energy differences ΔE^D3, in meV, from the F structures are also listed, as well as the percentage of the total dispersion energy arising from contributions involving Cu

	E ^D3F
CuBiP 2 Se 6
Total	−4.04	−4.13	−3.96	−83	84
	−2.24	−2.21	−2.25	30	−16
	−1.78	−1.87	−1.69	−86	95
Cu–Cu	−0.03	−0.05	−0.02	−27	4
% Cu	45	47	43	136	119
CuInP 2 S 6
Total	−3.66	−3.75	−3.56	−95	99
	−2.57	−2.54	−2.59	30	−16
	−1.04	−1.12	−0.93	−72	110
Cu–Cu	−0.04	−0.09	−0.04	−53	5
% Cu	30	32	27	131	116

---

To do this, we evaluate the D3 energy for subsystems containing either only Cu atoms or else no Cu atoms and compare the results to the full D3 energy of the system. This procedure is somewhat approximate as D3 is strictly not purely pairwise additive, with the van der Waals properties of each atom determined in the context of its environment, but the results will be qualitatively indicative of the role of copper in controlling inter-layer van der Waals interactions.

Table 2 shows that the total D3 contributions to bilayer energies range from near −3.6 eV per cell for CuBiP₂Se₆ to near −4.0 eV per cell for CuBiP₂Se₆, the difference arising owing to the higher polarizabilities of Se and Bi compared to S and In. For CuBiP₂Se₆, ca. 45% of the dispersion energy arises through interactions with Cu atoms, reducing to ca. 30% for CuBiP₂Se₆. However, in terms of the D3 energy differences between the A_ii or A_oo structures and F, ΔE^D3, in both cases the Cu-involving contribution is ca. 125% of the whole, with the contribution opposing it but having only one fifth of the magnitude. The direct Cu–Cu van der Waals interaction contributes only one sixth to one third of this total, with the remainder coming from interactions. However, both types of Cu contributions are much more attractive across A_ii interfaces than they are across A_oo.

Why the copper van der Waals forces are strongest over A_ii junctions and weakest over A_oo ones is sketched in Fig. 4. The most important van der Waals interactions involve the Cu of one layer interacting with the Cu and the electron-rich part of its neighbour. The magnitude of such interactions clearly scale in the order A_ii > F > A_oo. Also in Fig. 4 the basic ionic inter-layer interactions are also sketched, clearly favouring F over A_ii and A_oo. It is the balance achieved between the weak ionic forces and the strong but opposing van der Waals forces that controls the ferroelectric ordering in ABP₂X₆ laminar materials.

Fig. 4 The primary ionic and van der Waals (vdW) interactions between layers of ABP₂X₆ laminar materials. The ionic forces are weak but provide net stabilization of ferroelectric polymorphs compared to antiferroelectric ones, whereas the differences in inter-layer van der Waals interactions focus on the interaction of the Cu atoms with the electron-rich parts of the neighbouring bilayer, strongly stabilizing A_ii bilayers while destabilizing A_oo bilayers compared to ferroelectric ones.

Conclusions

What is revealed for both ABP₂X₆ laminar materials considered is that antiferroelectric structures intrinsically involve two different types of inter-layer interactions with significantly different van der Waals forces, one type being much more attractive than the van der Waals force acting between ferroelectric structures while the other is much less attractive. The difference between these two van der Waals bonding effects compared to the weak dipole–dipole stabilization energy operative in ferroelectric structures then controls the properties of crystals and nanoflakes.

Realizing the nature of the balance of forces controlling the structure of ABP₂X₆ laminar materials will facilitate new synthetic strategies for designing desirable materials, as well as new measurements to exploit bilayer properties. How this works out for the bulk properties of CuInP₂S₆ and CuBiP₂Se has been developed elsewhere,⁴¹ with a focus added to consider also how these factors contribute to ferroelectric, paraelectric, and antiferroeelectric ordering within a single layer (as has been observed in these and related materials^18,19). This occurs as external electric fields can manipulate the dipole–dipole interaction, with controlling field strengths being determined by the differential van der Waals interaction.

In general, the role of van der Waals forces and the manifold ways they can compete with traditional ionic and covalent bonding scenarios is a topic of great current interest,¹ with the results found here for ABP₂X₆ laminar materials providing another example. van der Waals forces not only effect self-assembly but also can control localized chemical structure. Which computational methods provide the best description of the van der Waals forces in these systems is also the subject of future work,⁴¹ as is understanding how different methods perceive the competition with chemical forces and, as a consequence, how better methods can be developed in the future.⁴⁷

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by resources provided by the National Computational Infrastructure (NCI), and Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia, as well as Chinese NSF Grant #1167040630. This work is also supported by the Australian Research Council grant DP16010130.

Notes and references

1. J. R. Reimers, M. J. Ford, S. M. Marcuccio, J. Ulstrup and N. S. Hush, Nat. Rev. Chem., 2017, 1, 0017.
2. J. R. Reimers, M. J. Ford, A. Halder, J. Ulstrup and N. S. Hush, Proc. Natl. Acad. Sci. U.S.A., 2016, 113, E1424.
3. P. K. Das, S. Samanta, A. B. McQuarters, N. Lehnert and A. Dey, Proc. Natl. Acad. Sci. U.S.A., 2016, 113, 6611.
4. A. Bilic, J. R. Reimers, N. S. Hush, R. C. Hoft and M. J. Ford, J. Chem. Theor. Comput., 2006, 2, 1093.
5. J. R. Reimers, M. Li, D. Wan, T. Gould and M. J. Ford, in Noncovalent interactions in quantum chemistry and physics: theory and applications, ed. A. Otero de la Roza and G. DiLabio, Elsevier, Amsterdam, 2017, p. 13.
6. N. Nerngchamnong, L. Yuan, D.-C. Qi, J. Li, D. Thompson and C. A. Nijhuis, Nat. Nanotechnol., 2013, 8, 113.
7. A. K. Geim and I. V. Grigorieva, Nature, 2013, 499, 419.
8. R. G. Pearson, J. Am. Chem. Soc., 1963, 85, 3533.
9. R. G. Pearson, J. Chem. Educ., 1968, 45, 581.
10. A. K. Tagantsev, K. Vaideeswaran, S. B. Vakhrushev, A. V. Filimonov, R. G. Burkovsky, A. Shaganov, D. Andronikova, A. I. Rudskoy, A. Q. R. Baron, H. Uchiyama, D. Chernyshov, A. Bosak, Z. Ujma, K. Roleder, A. Majchrowski, J. H. Ko and N. Setter, Nat. Commun., 2013, 4, 2229.
11. A. Simon, J. Ravez, V. Maisonneuve, C. Payen and V. B. Cajipe, Chem. Mater., 1994, 6, 1575.
12. V. Maisonneuve, M. Evain, C. Payen, V. B. Cajipe and P. Molinié, J. Alloys Compd., 1995, 218, 157.
13. V. Maisonneuve, V. B. Cajipe, A. Simon, R. VonDerMuhll and J. Ravez, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, 10860.
14. D. M. Bercha, S. A. Bercha, K. E. Glukhov and M. Sznajder, Acta Phys. Pol., A, 2014, 126, 1143.
15. 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, 3808.
16. M. A. Susner, A. Belianinov, A. Borisevich, Q. He, M. Chyasnavichyus, H. Demir, D. S. Sholl, P. Ganesh, D. L. Abernathy, M. A. McGuire and P. Maksymovych, ACS Nano, 2015, 9, 12365.
17. F. Liu, L. You, K. L. Seyler, X. Li, P. Yu, J. Lin, X. Wang, J. Zhou, H. Wang, H. He, S. T. Pantelides, W. Zhou, P. Sharma, X. Xu, P. M. Ajayan, J. Wang and Z. Liu, Nat. Commun., 2016, 7, 12357.
18. M. A. Gave, D. Bilc, S. D. Mahanti, J. D. Breshears and M. G. Kanatzidis, Inorg. Chem., 2005, 44, 5293.
19. V. Maisonneuve, V. B. Cajipe and C. Payen, Chem. Mater., 1993, 5, 758.
20. E. Durand, G. Ouvrard, M. Evain and R. Brec, Inorg. Chem., 1990, 29, 4916.
21. V. Carteaux, D. Brunet, G. Ouvrard and G. Andre, J. Phys.: Condens. Matter, 1995, 7, 69.
22. H. Ji, R. A. Stokes, L. D. Alegria, E. C. Blomberg, M. A. Tanatar, A. Reijnders, L. M. Schoop, T. Liang, R. Prozorov, K. S. Burch, N. P. Ong, J. R. Petta and R. J. Cava, J. Appl. Phys., 2013, 114, 114907.
23. C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G. Louie, J. Xia and X. Zhang, Nature, 2017, 546, 265.
24. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
25. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
26. W. Song, R. Fei and L. Yang, Phys. Rev. B, 2017, 96, 235420.
27. S. Grimme, J. Comput. Chem., 2006, 27, 1787.
28. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
29. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413.
30. S. N. Steinmann and C. Corminboeuf, J. Chem. Phys., 2011, 134, 044117.
31. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett., 2009, 102, 073005.
32. A. Tkatchenko, R. A. DiStasio, R. Car and M. Scheffler, Phys. Rev. Lett., 2012, 108, 236402.
33. A. Ambrosetti, A. M. Reilly, R. A. Distasio Jr and A. Tkatchenko, J. Chem. Phys., 2014, 140, 18A508.
34. B. Tomáš, L. Sébastien, G. Tim and G. Á. János, J. Phys.: Condens. Matter, 2016, 28, 045201.
35. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
36. J. Klimeš, D. R. Bowler and A. Michaelides, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 195131.
37. K. Lee, A. K. Kelkkanen, K. Berland, S. Andersson, D. C. Langreth, E. Schröder, B. I. Lundqvist and P. Hyldgaard, J. Phys.: Condens. Matter, 2011, 84, 193408.
38. J. P. Perdew and Y. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8800.
39. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 224106.
40. J. Moellmann and S. Grimme, J. Phys. Chem. C, 2014, 118, 7615.
41. S. A. Tawfik, J. R. Reimers, C. Stampfl and M. J. Ford, J. Phys. Chem. C, 2018 DOI:10.1021/acs.jpcc.8b05349.
42. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
43. L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi and S. Grimme, Phys. Chem. Chem. Phys., 2017, 19, 32184.
44. J. F. Dobson, Int. J. Quantum Chem., 2014, 114, 1157.
45. J. W. Ochterski, Thermochemistry in GAUSSIAN, https://www.gaussian.com/g_whitepap/thermo/thermo.pdf.
46. J. R. Reimers, D. Panduwinata, J. Visser, Y. Chin, C. Tang, L. Goerigk, M. J. Ford, M. Sintic, T.-J. Sum, M. J. J. Coenen, B. L. M. Hendriksen, J. A. A. W. Elemans, N. S. Hush and M. J. Crossley, Proc. Natl. Acad. Sci. U.S.A., 2015, 112, E6101.
47. M. Li, J. R. Reimers, J. F. Dobson and T. Gould, 2018, arXiv:1806.07532.

---

Footnotes

† Electronic supplementary information (ESI) available: All optimized atomic coordinates. See DOI: 10.1039/c8sc01274a

‡ Current address: School of Science, RMIT University, Melbourne, 3001, Australia.

---