Dimensional classification and shape reliability in In_xSe_y clusters and 2D materials: a quasi-molecule perspective from subminimal-basis DFT

Abstract

The reliable prediction of shapes of molecules and molecular materials from compact basis sets remains a key issue in electronic-structure theory. We examine here the dimensional classification – linear, planar or otherwise – of 2D In_xSe_y clusters and materials using subminimal-based density-functional-theory (DFT) within the quasi-molecule framework. By focusing on shape classes rather than numerical bond lengths and angles, the extremely compact STO-3G basis set is shown to preserve the correct dimensional character of chemical-bound systems. Optimization of small clusters with split-valence and polarized basis sets confirms that the overall structural topology so obtained is remarkably robust to basis-set truncation. The results, apparently reported here for the first time for 2D In₂Se₅, demonstrate that the essential features of bonding and geometry are already encoded in the valence-orbital symmetry of the quasi-molecular system. The analysis clarifies the conceptual limits of geometry as an emergent property of the Born–Oppenheimer surface and supports the rational use of subminimal bases for cluster explorations. Attention is paid to neutral molecules that can provide insights into crystalline lattices even when conventionally viewed as independent or unrelated components within the crystal structure. Possible implications of the popular periodic boundary conditions model are tentatively suggested. While of no apparent impact on the results reported here, the role of including relativistic effects has also been tested for some of the title systems.

---

1 Introduction

The isolation of a single layer of graphene has triggered worldwide interest in 2D materials.^1,2 However, the lack of an intrinsic band gap in graphene was a major obstacle in developing devices of interest in electronics. Attention then turned to 2D materials that overcame such a limitation. Among those exfoliatable into individual layers with interesting structural and electronic properties are 2D transition-metal dichalcogenides.³ With a small band gap, they form a new family of layered materials with novel technological properties in the ultrathin limit.^4–8

Due to the progress in preparing high-quality 2D indium selenide materials (InSe and In₂Se₃) using techniques like chemical vapor deposition⁹ and physical vapor transport,¹⁰ the gas-phase indium–selenium systems became of compelling interest.¹¹ These methodologies address key challenges in producing large-area, phase-pure materials, which are difficult with exfoliation techniques. A prime example is the electrically induced solid-state amorphization in In₂Se₃, which would not have been discovered without the ability to synthesize the pristine, phase-pure, materials. This calls for studying the fundamental properties of gas-phase In–Se clusters, in particular the intermediate species and kinetic pathways that govern the gas-to-solid transition. Such studies are also key in quasi-molecule theory, where the goal is to understand the structure of large In–Se clusters and associated 2D lattices from the tiles and generalized tiles embedded on them. This is at the focal point in the current work.

The most applied model in studying the title 2D materials is the periodic boundary conditions (PBC) model, which approximates the bulk by replicating a unit cell. However, it can lead to known artifacts that affect the calculated properties. Suffice it to mention a few: artificial periodicity, finite-size effects, lack of long-range correlations, and non-periodic shapes. By forcing periodicity, the PBC model imposes (rather than predicts) a repeating structure on materials that in reality may not be perfectly periodic. It may then lead to skewed results when simulating non-periodic phenomena like defects, interfaces, or cracks. Although the PBC model may perfectly tile a three-dimensional space by repeating the unit cell, it naturally becomes unsuitable when studying systems with non-tiling geometries. Moreover, in a PBC simulation, the particles can only exhibit motions with wavelengths that fit within the simulation box. Regarding the unit cell size, increasing the simulation box reduces artifacts but requires additional significantly larger computational resources. Approaches that take into account the fairness or difficulties in the PBC model are then welcome, a major aim of this work.

2 Theory: scope and strategy

It is well established that atoms are not point-like, although such an assumption is key to explaining linearity or planarity of molecules. Accordingly, the possibility that a small molecule can be treated as a ‘pseudo-atom’ has been advanced,^12,13 and can be particularly valid if its internal structure and motions are (at least partly) ignorable. In fact, this assumption has been found to be fair in scenarios where such motions are assumed to be of little relevance for the envisaged problem, which to some extent is the case here also, as later discussed. Indeed, although the internal degrees of freedom are assumed not to play a key role, their geometries may get distorted (through deformations of their structures or twist motions around their axes), in a way judged similar to what happens with atoms when their assumed (or eventually true) sphericity is lost upon bond formation. Naturally, in isolation, the atomic electron clouds are spherically symmetric (especially the s-orbitals) in the absence of external (or most internal nuclear) fields distorting them. Conversely, rather than round, atoms in molecules take on distorted shapes (more like lobes or regions pointing toward the bonded atoms). Nevertheless, all degrees of freedom of such pseudo-atoms are later considered when calculating the forces in the optimization of the parent molecule. To our knowledge, such an approach has only recently been explicitly adopted.

Consider then: can the geometry of a molecule or molecular material be rationalized as pseudo-planar or otherwise from the geometries of the quasi-pseudo-tetra-atomic tiles (pseudo-tiles) it embeds? The answer is positive as implied by the

Lemma: if all tetrads of points (atoms or pseudo-atoms) in a set are on a plane, they are all on the same plane.

which is valid for atoms (not true points) and, in the present work, under the assumption that In₂Se₃ or In₂Se₅ can be treated as pseudo-atoms. This Lemma (called elsewhere^12,14–16 Lemma 2; Lemma 1 is for linearity,¹⁴ the tiles being triatomics¹⁷ or pseudo-triatomics¹²) has been discussed before,¹⁴ where the reader is addressed for details.

Of course, the concept of representing an atom as a point is key for classifying linearity and planarity of molecules, hence for applying the Lemma. The question is then, where is the point that represents the pseudo-atom located? Two possible candidates are the pseudo-atom center-of-mass and its center-of-charge. In the absence of an answer, rather than setting a choice from the beginning, we leave it open for discussion when analyzing the results in Section 5. Any of these choices is acceptable for atoms, since both the center-of-mass and center-of-charge coincide in this case with the geometric center (centroid).

Rather than rationalizing the results for the title systems from specific optimizations, the aim is to provide a background for predicting molecular geometries (linearity, planarity or otherwise) without resorting to direct calculations or experimental data. Of course, in the absence of the latter, direct ab initio or DFT calculations must be done to probe the predictions made from using quasi-molecule theory for the parent molecules.

Although the structure of a stable molecule is a minimum, predicting a saddle point may suffice in quasi-molecule theory. Uncommon to some extent, recall that one cannot distinguish these two stationary points in quasi-molecule theory (the Lemma involves no forces) before a push along a suitable off-planar coordinate. Of course, a single point calculation along such a coordinate suffices to disclose its nature. Recall that a saddle point may still behave as a minimum if vibrational stabilization is brought into play,^18,19 a topic of no concern in the present work.

A further remark to note is that one can disclose the nature of stationary points in quasi-molecule theory by employing a bisection method,^14–16 thus avoiding to perform ab initio or DFT calculations for all possible tiles (involving, in the case of planarity, four atoms or pseudo-atoms) embedded on the parent molecule. This greatly facilitates things, since the combinatorial law could otherwise make the task tedious or even unaffordable.

As noted above, for materials, quasi-molecule theory avoids the computational assumption of periodic boundary conditions, while scaling computationally efficiently with system size. Moreover, it offers chemically intuitive insights based on local bonding motifs. There are a few caveats though: (a) in cases where pseudo-atoms have significant internal asymmetry (e.g., In₂Se₃ and In₂Se₅), some deviations from planarity are to be expected; (b) unlike traditional molecular modeling, quasi-molecule theory does not compute forces or energy gradients for the global structure, with planarity being inferred from the geometric arrangement of the tiles, hence not enforced via optimization constraints nor resorting to empiricism [even if already stated, recall that stationary points can only be minima or saddle points of a given index; the former can be distinguished from the latter through a single-point calculation along an appropriate coordinate for off-planar motion]; and (c) the nature of the predicted stationary points of the potential energy surface^20,21 (PES) is commonly not (but could be) explored.

3 Quasi-molecules and quasi-pseudo-molecules

A conceptual subtlety is the definition of the electronic state of a quasi-molecule or quasi-pseudo-molecule, therein also often simply denoted as quasi-molecule. Of course, each molecule has a well-defined spin and spatial symmetry. In contrast, a quasi-molecule is a virtual construct: its existence precedes dissociation of the parent and hence it can formally occupy any state allowed by the spin–spatial coupling rules.^22,23 This idea draws inspiration from the quasi-atom concept of Ruedenberg and co-workers,^24–27 which highlights the utility of pre-fragmentation states for interpreting bonding and structure. In fact, unlike chemical reactions where reactants and products are infinitely separated, quasi-molecular tiles are defined within the parent structure, prior to dissociation,¹⁴ and hence cannot a priori be assigned an electronic state. The orbitals of the quasi-molecules may then be viewed as distorted molecular functions embedded in the parent wave function; artificial in the quasi-molecule, they become the molecular orbitals after fragmentation. Moreover, the quasi-molecules need not correspond to minima on the PES. In fact, saddle points may actually suffice for identifying meaningful geometric trends, provided vibrational stabilization or other entropic effects support their persistence^18,19 (even if repeatedly stated, recall that the Lemma invokes stationary points). Note further that the pseudo-atoms in a pseudo-molecule may upon optimization show an energy not rigorously relatable to the sum of the isolated pseudo-atom energies. For this reason, the popular account for the basis set superposition error²⁸ via counterpoise correction^29–31 can hardly be done, with extrapolation^32–37 of the uncorrected energies to the complete basis set (CBS) limit providing a more robust approach. Suffice it to recall now the multitude of stationary points (eventually minima) that the system may have, with the pseudo-atoms themselves evincing slightly distinct geometries.

Applying quasi-molecule theory to 2D In_xSe_y materials requires identifying all unique tetrads within a given In_xSe_y cluster. Although it is sufficient to warrant that they correspond to stationary structures, the calculation of forces to identify their topographical nature is done to probe the predictions reported here. Of course, experimental data can be used when available, thus making such calculations unnecessary. While the use of tiles would become combinatorially intensive and eventually unaffordable for large systems, such a task can be considerably simplified for large molecules by employing generalized tiles (planar leftovers) in conjunction with the bisection scheme,^14–16 which partitions efficiently the parent molecule while allowing evaluation of relevant dihedrals and planarity metrics.

4 Basis set selection: STO-3G, how?

The use of subminimal^32,33 (commonly known as minimal) basis sets, such as STO-3G, follows the molecular orbital construction introduced by Pople and co-workers.^38,39 Despite not having polarization functions, they reproduce the correct dimensional class of molecular geometries (linear, planar and non-planar), provided that the essential valence symmetries are included.^40–42 This qualitative robustness arises because the topology of the bonding framework is primarily determined by the symmetry of the occupied valence orbitals rather than by quantitative angular flexibility.^32,43,44 As noted by Woolley and Sutcliffe, molecular geometry is not an observable but an emergent concept derived from the Born–Oppenheimer (BO) PES.^45,46 When geometry is interpreted in this sense, the subminimal basis remains sufficient for reliable shape classification.^40,47,48 The heavier elements may occasionally require polarization functions to describe inner-valence angular distortions, but the shape topology of most other-element compounds remains unchanged.^6–8,41

As it is well established, subminimal basis sets include just enough functions to represent each occupied atomic orbital. They are about the smallest basis sets one can use to perform ab initio calculations, thus playing a key role when attempting to overcome the curses of size (physical- and basis-set-like) in quantum chemistry. Although polarization and diffuse functions are critical for an accurate energetics, one wonders whether the small size of STO-3G can be a serious problem, since bonding and molecular geometry are sensitive to how orbitals polarize and delocalize.

Subminimal basis sets are traditionally used when qualitative trends or very large systems are at stake, not to mention when using semiempirical or tight-binding methods.⁶⁰ They are also employed for scans before a higher-level optimization is done. So, they can give an idea of bonding but cannot be quantitatively accurate for determining a molecular structure. One must also recall the interplay between basis set limitations and method deficiencies: even with an exact method, a subminimal basis is likely unable to get the right orbital polarization. This is the basis set incompleteness error,^23,33,35 not a methodological one. There is then a hierarchy of limiting factors: the method controls energy accuracy and correlation, while the basis controls wavefunction flexibility and hence the molecular shape.

A few misleading cases are often invoked involving hypervalent species. An example is SF₆, where d-polarization is expected to produce the correct angular distribution. A subminimal basis may then be judged to be qualitatively problematic. A quick look is given to SF₆ by optimizing its STO-3G geometry. The result is a minimum with a qualitatively correct shape: an octahedral geometry of bond length 1.736 Å, which compares with 1.616 Å when using Dunning's^47,50 correlation consistent VDZ basis set, or 1.618 Å and 1.590 Å when employing its augmented congeners AVDZ and AVTZ, respectively. The difficulty may then arise with transition metal atoms. Because some are considered here, this will allow judgement on their study with the STO-3G basis.

It should be emphasized that no accurate determination of a molecular geometry is planned here, but instead our goal is to be able to predict whether a molecule is linear, planar or otherwise, the only classes into which shape can be clear-cut catalogued. We are then sticking to shape reliability, meaning bent vs. linear, planar vs. pyramidal, and square vs. tetrahedral, all this independently of whether the bond angles or distances are numerically highly accurate. It is then claimed that it is uncommon that subminimal bases predict a wrong structure even for small molecules. It may though be conceived to happen when the molecular structure is governed by subtle energy differences, since subminimal bases lack the flexibility required to mimic such differences. The STO-3G basis is then judged as giving the right shape for most molecules, except (but rarely) when the geometry depends critically on polarization, delocalization or subtle electronic correlation.

In support of the above, recall that shape comes from the qualitative features of the PES:^20,22,23 if this has a deep, well-defined minimum at a given geometry, even a subminimal basis will land there. If a nearly flat PES is at stake, any approximation (basis, method, integration grid) may shift the molecule's shape. So, for chemically bound systems, the shape is robust under a subminimal basis, while for soft or competing systems it is ambiguous until the basis and correlation (at least) are correctly treated. The STO-3G basis should then give the chemically correct shape for most systems because the topology of bonding is right. In fact, the molecular shape is physically ill-defined beyond that, because it is not an observable but an emergent concept.^43,44 To emphasize: if shape means linear, planar or otherwise, a subminimal basis should not fail for covalently bound molecules since it captures the correct dimensional class which arises from the symmetry and connectivity of the valence orbitals, not their fine polarization.

How can one be sure of the exact structural form in subtle cases? Possibly only experimentally, although experiments are themselves not error free. Or by solving Schrodinger's equation, but can it be done exactly? Because the wavefunction Ψ(r,R) depends on both the nuclear and electronic coordinates, the molecular geometry is just the position of the minimum in the BO PES. The best one can get are expectation values (bond lengths, angles) averaged over nuclear motion. So, even conceptually, there is no exact geometry, only R_e (equilibrium geometry of PES minimum) and 〈R〉 (average geometry) from the vibrational wavefunction. Experimentally, one infers geometries from microwave or rotational spectroscopy (effective bond lengths), electron or X-ray diffraction (average nuclear positions), and vibrational spectra (force constants). All refer to vibrationally averaged structures; even at 0 K, the molecule still vibrates (zero-point motion, ZPE), and hence the observed geometry is never the classical PES minimum. Moreover, for molecules with flat or multiple-well PESs (e.g. fluxional species,^51,52 Jahn–Teller distortions⁵³), the vibrational ground state delocalizes across both minima, which makes the notion of planar vs. pyramidal even more meaningless.

In summary, neither experiment nor theory gives a single, immutable structure. When classifying the latter as linear, planar or otherwise, one is describing the topology of the equilibrium PES. Having the above in mind, the suggested dimensional class is likely the most defensible level of structural inference. Indeed, the STO-3G basis captures all the qualitative features that are chemically meaningful. Possible exceptions (hypervalent, multicenter, or transition-metal systems) cannot be diagnosable with a low-level model as long as their bonding pattern is not defined by valence topology or the PES involves competing electronic configurations. Knowing which shape the molecule assumes would imply to know the exact PES, which is unaffordable for most many-electron species. Broadly stated, for ordinary molecules (99% of the chemistry one deals with), STO-3G plus a competent method can safely be used to infer the molecule's shape class. For exotic bonding motifs, even large basis plus correlation is not guaranteed to give an unambiguous geometry. To sum up, STO-3G offers an essentially correct reasoning that is conceptually defensible within the limits of what molecular shape can mean. Such a view is to some extent supported by the extensive use⁵⁴ of subminimal basis sets like STO-3G⁵⁵ or variants more recently developed^56–59 (and references therein).

5 Results and discussion

Since we are not aware of any theoretical studies for some of the In–Se clusters reported here, we survey the calculations here done from diatomics up to In₁₂Se₁₈ and In₁₀Se₂₅, the latter already involving 1200 and 1340 electrons, respectively. The present study employs computational tools to determine the lowest-energy geometric structures of all studied In_xSe_y clusters, revealing how atomic arrangement and composition influence stability and electronic properties. As noted elsewhere,^12,13 two points must be emphasized. First, for practical reasons, the title systems can only be studied using DFT. Second, to avoid numerical inaccuracies, the same atomic basis set should be employed when performing a comparative study of distinct members of the studied family. Although an empirical approach such as DFTxTB could offer a fast route,⁶⁰ the choice was to use a priori unrestricted DFT. As for the atomic basis sets, three have been mostly employed: two are subminimal (STO-3G, the simplest and computationally most advantageous, and 3-21G), basically the kind of basis sets affordable for the largest systems studied here, and the other is an extended triple-zeta basis set (ATZP) only employed for systems up to In₂Se₃. For testing purposes, the corresponding minimal basis sets (DZP and ADZP) are also occasionally employed.

5.1 Gas phase In–Se clusters up to tetratomics

To help decide on possible basis set limitations, further calculations of small clusters of In and Se atoms up to tetratomics have been carried out. Table 1 gathers the results so obtained. In their ground states, Se has the electronic configuration [Ar]3d¹⁰4s²4p⁴ (³P₂), while for In it is [Kr](4d¹⁰5s²5p¹) (²P^o_1/2). Although some values of 〈S²〉 deviate slightly more than 10% from their pure values for the singlet, doublet, triplet, and quartet states (0, 0.75, 2.0, and 3.75 in the same order), all spin contaminations were considered acceptable, except where indicated, having in mind that all calculations were performed at the unrestricted DFT level.

Table 1 Energies (from the code output, in E_h) for the In_xSe_y clusters studied in the present work up to tetratomics, dipole moment (in Debye), nature of the calculated stationary point^a, and spin state 〈S²〉

Cluster	State	STO-3G				ATZP
		Energy	μ	〈S^2〉	ΔE	Energy	μ	〈S^2〉	ΔE
a Unless specified otherwise, the reported structures are minima.b The linear doublet structure, shown in the first entry and predicted with STO-3G, corresponds to a linear saddle point with energy −17 053.93692647Eh and imaginary frequency −48.9 cm^−1, and is 37.3 kcal mol^−1 above the minimum.c Calculations including relativistic effects; both electronic states are linear minima. The relativistic ATZP-DKH basis set is unavailable; a similar result is obtained with TZP-DKH.d The reported reference doublet structure at the B3LYP/ATZP level is a saddle point with imaginary frequency −54.4 cm^−1.e The planar doublet structure (shown in the Graphical abstract) is a saddle point at the B3LYP/ATZP level, with energy −19 428.79120315Eh and imaginary frequency −357.3 cm^−1, and is 28.2 kcal mol^−1 above the minimum. Both doublets show spin contamination.
Se2	^1A	−4749.41331858	0.000	0.000	24.689	−4802.87116678	0.000	2.007	19.924
	^3A	−4749.45266373	0.000	1.999	0.000	−4802.90291774	0.000	2.007	0.000
InSe	^2A	−8059.38799360	3.162	0.756	0.000	−8139.02133293	2.795	0.757	0.000
	^4A	−8059.25277771	2.795	3.752	84.849	−8138.91886684	0.044	3.758	64.298
In2	^1A	−11 369.21831708	0.000	0.000	18.891	−11 475.09856430	0.000	0.000	0.000
	^3A	−11 369.24842161	0.000	2.002	0.000	−11 475.04291704	0.001	2.002	34.919
InSe2	^2A	−10 434.19947495	0.000	0.924	0.000	−10 540.51858189	0.000	0.776	0.000
	^4A	−10 434.10564193	1.949	3.767	58.881	−10 540.45823099	1.333	3.778	37.871
In2Se	^1A	−13 744.09995327	0.757	0.000	0.000	−13 876.66625649	1.516	0.000	0.000
	^3A	−13 744.06999771	1.697	2.005	18.797	−13 876.61320911	1.505	2.006	33.288
In2Se2	^1A	−16 118.95286916	1.236	0.000	0.000	−16 278.16157881	0.001	0.003	0.000
	^3A	−16 118.86241739	0.190	2.042	56.759	−16 278.09262164	3.103	2.015	43.271
Se3	^1A	−7124.19905644	0.000	0.000	0.000	−7204.35696461	0.000	0.000	0.000
	^3A	−7124.16574221	0.819	2.006	20.905	−7204.34186660	0.975	2.009	9.474
In3^b	^2A	−17 053.93692647	0.000	0.754	15.445	−17 212.70680490	0.000	0.758	11.964
	^4A	−17 053.96153933	0.458	3.755	0.000	−17 212.72587143	0.003	3.754	0.000
In3^c	^2A	−17 338.87506097	0.000	0.763	10.116	−17 545.84233303	0.000	0.751	6.536
	^4A	−17 338.89118155	0.000	3.755	0.000	−17 545.85274806	0.000	3.754	0.000
InSe3^d	^2A	−12 808.92887528	2.954	0.760	0.000	−12 941.98834772	0.649	0.757	0.000
	^4A	−12 808.83111482	1.614	3.764	61.346	−12 941.94218196	0.794	3.763	28.969
In3Se^e	^2A	−19 428.83607947	2.488	1.600	0.000	−19 614.26843707	2.330	1.280	0.000
	^4A	−19 428.83228305	1.693	3.762	2.382	−19 614.22153579	0.797	3.760	29.431

---

As shown in Table 1 and Fig. 1 and 2, the results calculated for both Se₂ and In₂ are mostly in agreement with their congeners in the corresponding upper rows of the periodic table, namely O₂ and Al₂, as well as the results reported in the NIST Chemistry WebBook. Regarding these homonuclear diatomics, the results in Table 1 are in agreement for Se₂, irrespective of the employed atomic basis set (STO-3G or ATZP), with the triplet being their ground state; for convenience, the electronic states are mostly classified in the C₁ symmetry point group. It turns out that the indium dimer shows a low vibrational frequency (11.0 cm⁻¹) even with the largest employed basis set. In this case, the triplet state is the lowest with the STO-3G basis, while the singlet is the ground state at the ATZP level. Due to shortage of available data, calculations have also been done at the high cost CCSD(T)/ATZP level of theory, but no attempt made to account for relativistic effects; however, see later. The calculated results are given in Table 2, and illustrated in Fig. 2. Because, to our knowledge, the cost-expensive harmonic vibrational frequencies for such clusters have not been reported, they are gathered in Table 3. As seen, the B3LYP/STO-3G results support the B3LYP/ATZP approach, although apparently suggesting some lack of accuracy in the bond distances and energetics.

Fig. 1 Optimized structures up to tetratomic In–Se ground-state (lowest energy) clusters using the STO-3G (left-hand-side) and ATZP (right) basis sets: (a) and (a′) Se₃; (b) and (b′) In₃ [the structure in panel (b) shows D_∞h symmetry when including relativistic effects; see the text]; (c) and (c′) InSe₂; (d) and (d′) In₂Se; (e) and (e′) In₂Se₂; (f) and (f′) InSe₃; and (g) and (g′) In₃Se. Structure (g) for In₃Se at the B3LYP/STO-3G level of theory is the only non-planar one, but 〈S²〉 = 1.600; a planar sp structure with the proper spin (〈S²〉 = 0.756) is also predicted, but it is 28.160 kcal mol⁻¹ above the tabulated minimum. Similarly, at the B3LYP/ATZP level, (g′) shows 〈S²〉 = 1.28. See Tables 1 and 2.

Fig. 2 Excitation energies for up to tetratomic In_xSe_y clusters when DFT/B3LYP optimized with the STO-3G and ATZP basis sets, and CCSD(T)/ATZP. When negative, the higher-spin state is the ground state. See also the text.

Table 2 Energies (from the code output, in E_h) for the In_xSe_y clusters studied in the present work up to tetratomics at the CCSD(T)/ATZP level of theory. Indicated in the seventh column is the excitation energy for the first electronic state

Cluster	State	Energy^a	Shape^b	Stat. pt.	〈S^2 − Sz^2 − Sz〉^c	ΔE
a Even if not all significant, decimal figures are kept as in code's output; “sp” stands for saddle point, “min” for minimum.b Approximate symmetry.c Spin contamination.d Weakly bound, with a bond length of 3.705 Å and a vibrational frequency of 62.9 cm^−1.e Bound, with three distinct sides.f Bound, L-shaped type.g Bound, open structure.h Highly cost-expensive (>35 d of CPU time for optimization with ATZP basis, ∼2 m with calculation of forces).i Assuming the doublet structure in panel (c) of Fig. 3.
Se2	^1A	−4799.93873623	D∞h	min	0.000	15.076
	^3A	−4799.96276095	D∞h	min	0.000	0.000
InSe	^2A	−8135.24007306	C∞v	min	0.000	0.000
	^4A	−8135.14363944	C∞v^d	min	0.000	60.513
In2	^1A	−11 470.48643867	D∞h	min	0.000	0.000
	^3A	−11 470.43460172	D∞h	min	0.000	32.528
InSe2	^2A	−10 535.26259194	D∞h	min	0.000	0.000
	^4A	−10 535.18314533	D∞h	sp	0.001	49.854
In2Se	^1A	−13 870.57919879	Cs^e	min	0.000	0.000
	^3A	−13 870.52477172	C2v	min	0.000	34.154
In2Se2	^1A	−16 270.60359598	D∞h	min	0.000	0.000
	^3A	−16 270.53837189	Cs^f	min	0.000	40.929
Se3	^1A	−7199.95556285	D3h	min	0.000	0.000
	^3A	−7199.92980891	C2v^g	min	0.000	16.161
In3	^2A	−17 205.76974412	D∞h	min	0.000	17.413
	^4A	−17 205.79749293	D∞h	min	0.000	0.000
InSe3	^2A	−12 935.24657531	C2v	sp^h	0.002	0.000
	^4A	−12 935.21744929	C2v	min	0.000	18.277^i
In3Se	^2A	−19 605.87090026	C2v	min	0.000	0.000
	^4A	−19 605.83317363	C2v	min	0.000	23.674

---

Table 3 Harmonic vibrational frequencies (ω_i, in cm⁻¹) of stationary points for up to tetratomic In_xSe_y clusters studied in the present work at the CCSD(T)/ATZP level of theory^a

System	State	ω1	ω2	ω3	ω4	ω5	ω6
a Empty spaces are not relevant: there are (3n − 6) vibrational frequencies, where n is the number of atoms.b Imaginary (negative) frequency: in this case, 77.0i cm^−1.
Se2	^1A	355.5
	^3A	380.0
InSe	^2A	222.2
	^4A	62.9
In2	^1A	128.6
	^3A	12.8
InSe2	^2A	86.6	90.2	219.0	354.2
	^4A	−77.0^b	−67.3	139.5	220.1
In2Se	^1A	28.1	204.0	220.4
	^3A	176.1	357.4	411.3
In2Se2	^1A	48.1	53.1	78.1	144.3	329.1	374.1
	^3A	27.0	55.4	75.4	148.3	232.0	337.1
Se3	^1A	247.7	247.7	331.4
	^3A	101.7	283.3	331.4
In3	^2A	61.6	116.4	237.9
	^4A	52.0	129.9	348.6
InSe3	^2A	−76.5	96.2	163.4	177.0	265.1	292.8
	^4A	60.9	70.6	98.0	182.9	280.4	306.0
In3Se	^2A	39.2	60.6	73.1	138.1	159.1	229.1
	^4A	61.5	105.6	116.4	158.7	186.3	233.3

---

Regarding Se₃, it is predicted with a cyclic (D_3h) geometry, but the literature reports it as structurally ambivalent: both open (C_2v) and cyclic (near D_3h symmetry) conformers lie within a few meV of each other, making it sensitive to computational methods or experimental conditions. Both DFT and CCSD(T) results here presented predict its optimum ground-state structure as D_3h. The predictions may then vary with theoretical methods, which underscores the importance of method selection for accurate results. This may be crucial, since dissociation of Se₃ makes its photoelectron spectra difficult to interpret, a well-documented issue.⁶¹

As for In₃, the commonly cited work is by Feng and Balasubramanian,⁶² who report its ground state to be . Modern theoretical or experimental confirmation appears to be absent or not easily discoverable in publicly indexed literature. To resolve the ambiguity, one would need to seek experimental work from main-group cluster spectroscopy or more recent high-level computational studies, e.g., using relativistic quantum chemistry. Unfortunately, both seem to be nonexistent. It turns out that its minimum is here predicted to be triangular at the B3LYP/STO-3G level, but linear with the larger ATZP basis set at both B3LYP and CCSD(T) levels of theory; see also later.

In the gas phase, indium selenides exist as small, discrete molecular clusters with electronic states determined by the specific In–Se bond and cluster geometry. Unlike bulk materials, which have continuous electronic bands, gas-phase clusters exhibit quantized energy levels. The exact electronic structure depends on factors like the number and arrangement of In and Se atoms in the cluster, leading to complex electronic transitions and optical properties. Feng and Balasubramanian⁶² investigated 7 low-lying electronic states of In₃ and reported the state to be the lowest, with lying somewhat higher in energy. They concluded that In₃ prefers a high-spin ground state, due to relatively weak p–p overlap and exchange stabilization. Later DFT/CCSD(T) studies (2000s onwards) predicted the doublet as the ground state, especially when including relativistic effects and better correlation. Somewhat surprising to some extent, using the same DFT approach (B3LYP), our calculations with the subminimal STO-3G basis set predict the doublet to be the lowest than the quartet by 0.985 eV. The reverse (linearity) is predicted with the extended ATZP basis set, which places the quartet −0.519 eV lower than the doublet. A similar result is observed at the CCSD(T)/ATZP level in Table 2, where the quartet is the lowest by 0.755 eV. Moreover, the energy gap between the quartet and doublet states is reported in the literature as relatively small (often <0.1 eV), and the possibility of switching the ordering is often advocated to depend on the level of theory. To our knowledge, direct experimental confirmation of the spin multiplicity of free gas-phase In₃ is lacking. In fact, most experimental work on In clusters focused on larger aggregates with photoelectron spectra and ionization energies, but not detailed spin assignments for In₃. To sum up, the current consensus from theory is early ab initio favouring a quartet ground state , while higher-level DFT/CC often favors the doublet ground state, sometimes the quartet. This is to some extent here replicated with B3LYP where linearity is predicted with the largest basis for both optimized structures [a result common to CCSD(T)] but a nearly equilateral triangular shape with the smallest basis set. Experimentally, there is no unambiguous direct confirmation yet. Hence, the ground state structure of In₃ does not appear to be settled: it is either a state or a nearly degenerate doublet , depending on the method and basis set (see also later).

Given the possibility of obtaining distinct results for the homonuclear di- and triatomic clusters depending on the methodology, suffice it to say from the results reported in Table 1 for the mixed clusters that the relative B3LYP predictions of the various electronic states do not seem to drastically change with the basis set. As summarized in Fig. 1, there is general agreement on the shapes of the predicted stationary points. Of course, differences arose in comparing the B3LYP results in Table 1 for In₃ [which is a quartet at B3LYP/ATZP and CCSD(T)/ATZP levels of theory but a doublet with B3LYP/STO-3G], besides difficulties due to spin contamination for the doublet In₃Se; see the following subsection for InSe₃ and Section 5.2 for InSe.

The incorporation of relativistic effects in addition to electron correlation is expected to be key for species involving heavy atoms. In fact, it is recognized to be mild for atomic numbers 40 ≤ Z ≤ 60, but crucial above that range and negligible below. Reviews on this topic include the early paper by Pyykko⁶³ and the more recent one by Nakajima and Hirao.⁶⁴ A very useful overview has also been reported by Truhlar and Li,⁶⁵ while the papers by Feng and Balasubramanian⁶² and Dolg et al.⁶⁶ should be cited in relation to indium chalcogenides, just to mention a few. Specifically, the paper by Feng and Balasubramanian is one of the earliest systematic theoretical studies of In–Se molecular species that points to equilibrium In–Se bond lengths in the range 2.6–2.7 Å, with relativistic electron core pseudo-potentials (ECPs) noted to be required for reliable bond distances. In turn, Dolg et al.⁶⁶ employed relativistic ECPs widely in the calculations of indium compounds known to cause the typical contraction of bond lengths. The question then arises whether the discrepancy reported for the B3LYP results of In₃(²A) in Table 3 [where the B3LYP method predicts a linear saddle point with the STO-3G basis (and a bent minimum) but only a linear minimum when ATZP is employed] is attributable to using the subminimal basis set or the lack of including relativistic effects (most likely, both). Although the prediction of a saddle or a minimum has no implications in the present approach (which focuses on stationary points), the In₃ results may already call for a clarification.

Molpro⁶⁷ supports relativistic effects for accurate calculations on heavy elements, including Douglas–Kroll–Hess (DKH) up to the 99th order, eXact-2-Component (X2C) Hamiltonians, and relativistic ECPs. Test calculations are therefore here also reported for In₃ employing the DKH method (activated with DKHO = 8) and the recommended exponential parametrization (DKHP = 2); see Section 5.1.1 for InSe₃ and Section 5.2 for the doublet of InSe. In fact, the inclusion of relativistic effects in all-electron DFT calculations via scalar relativistic corrections to the kinetic energy and the nuclear attraction energy has become rather common. The above referred bond-length contraction is observed (2.185 Å vs. 2.136 Å for the D_∞h saddle point of index 1 and a minimum obtained with the STO-3G basis set at the non-relativistic level, respectively). These results corroborate the non-relativistic ones in predicting stationary points for such structures except for the saddle at the B3LYP/STO-3G level that is now a minimum. Rather than using basis sets optimized for relativistic atoms,⁶⁵ the above traditional ones were used. The actual use of the former is tested for InSe and InSe₃ in Section 5.1.1. So, except when explicitly stated, the majority of the calculations reported here are non-relativistic.

5.1.1 InSe₃: a case study. Let us start by asking the question, is STO-3G preferable to the subminimal 3-21G split-valence (double-ζ) basis which is just slightly more cost-expensive? Recall that STO-3G has exactly one contracted Gaussian function per AO (each contracted function made of 3 primitive Gaussians chosen to best mimic a Slater-type orbital), and hence limited flexibility (the “shape” of an orbital cannot change much, since there is only one radial function per AO). Instead, the core orbitals in the 3-21G subminimal^32,33 basis set are represented by one contracted function built from 3 primitives, while the valence orbitals are split into two radial functions (one contracted from 2 primitives, the other from 1 primitive). So, the valence density can “breathe” in this case and adapt to bonding, thus being somewhat more flexible in the valence region. One then hopes 3-21G to yield better geometries, energies, and bonding descriptions. However, a further question is, since both are very far from complete, is the danger in making predictions with 3-21G increased when compared to STO-3G, which has the proper formal valence bonding? One is then led to suspect that the worst basis set errors may not be in subminimal sets like STO-3G, but in other small, variationally more flexible basis sets that cannot describe real correlation properly. In fact, 3-21G can often misbehave more than STO-3G simply because STO-3G's rigidity protects it from artifacts, thus generally predicting the formal valence structures, no artificial hyperconjugation, and less tendency to overpolarize/overdelocalize the electron density. Indeed, the larger flexibility can only warrant a more negative energy in the variational sense, not necessarily improved relative energies, etc. As a result, although STO-3G is far from a converged basis (hence open to quantitative errors), its subminimal nature will likely prevent the above misbehavior.

Since, for the many reasons pointed out in previous sections, one must pick a cost-effective basis set, STO-3G is expected to give stable (possibly realistic) structures. Additionally, because a cluster of one In and three Se atoms may be formed in many of the larger ones studied here, it is plausible to test the above basis sets for InSe₃ by considering both the DFT and CCSD methods, jointly with the flexible DZP, ADZP, and ATZP basis sets. Due to the high atomic number of In (ref. 49), both non-relativistic and scalar relativistic calculations are performed.

Fig. 3 evinces the results so obtained when non-relativistic bases are employed. Although other stationary points exist (no attempt is made to be exhaustive nor is it required for our main purposes; however, see later), the following gives per se an idea of how difficult it can be to extract a definite conclusion for the molecules (and materials) here at stake. To be as impartial as possible, all reported optimizations started from essentially the same geometry. The salient feature is that B3LYP/STO-3G predicts a minimum at the converged geometry, while all other basis sets yield a saddle point structure (one imaginary frequency for the twist motion of InSe and Se₂). Also visible is the slight tendency at the B3LYP level to diminish the imaginary frequency with basis set enhancement. This may suggest that a larger basis could even yield a smaller imaginary frequency, eventually a minimum. Because the STO-3G basis set performs well bond-length-wise when compared to the expensive ADZP and ATZP basis sets, this conveys a good reason for using it for larger clusters. Despite quite expensive, forces were also estimated for the above optimized structures at the CCSD(T)/ATZP level of theory; see Table 4.

Fig. 3 Optimized structures of InSe₃ and the lowest harmonic vibrational frequency (in red). All are planar, with energies in E_h and bond lengths in Å; for the entire set of frequencies at the CCSD(T)/ATZP level (hereinafter given in cm⁻¹), see the SI. Shown in italics in panels (c) and (e) are the bond distances, energies, and imaginary harmonic vibrational frequencies (only the lowest, when all are real) obtained when including relativistic effects; see the text.

Table 4 Center of mass planes

Pseudo-tile	Atoms #s	Centroid^a	Normal^b	(a,b,c,d) in^c ax + by + cz + d = 0
a x, y, z coordinates of the centroid.b x, y, z coordinates of the normal.c Parameters in equation ax + by + cz + d = 0 of fitted plane; see also the text.
In8Se12	In (5–8)	(−0.35209, −0.41417, +0.96003)	(0.41730, 0.84967, 0.32237)	(0.41730, 0.84967, 0.32237, 0.18935)
	In (13–16)	(0.66679, 0.51584, −1.08185)	(0.04762, 0.66695, −0.74358)	(0.04762, 0.66695, −0.74358, −1.18023)
	Se (1–4)	(−0.32354, −0.37970, 1.74791)	(0.18493, 0.80914, 0.55828)	(0.18493, 0.80914, 0.55828, −0.06797)
	Se (9–12)	(−0.22179, 0.02115, −0.75369)	(0.10265, 0.80123, 0.58952)	(0.10265, 0.80123, 0.58952, 0.10889)
	Se (17–20)	(−0.52012, 0.53883, −1.55037)	(0.02717, 0.85734, 0.51402)	(0.02717, 0.85734, 0.51402, −0.10522)
In8Se20	In (9–12)	(−0.68065, −0.79163, 1.75329)	(0.01929, −0.06327, 0.99781)	(0.01929, −0.06327, 0.99781, −1.78641)
	In (17–20)	(0.67466, 0.26254, −2.01685)	(−0.10347, 0.19146, 0.97603)	(−0.10347, 0.19146, 0.97603, 1.98805)
	Se (1–8)	(−1.26468, −0.58984, 3.07086)	(0.04773, −0.26698, 0.96252)	(0.04773, −0.26698, 0.96252, −3.05288)
	Se (13–16)	(0.58611, −0.46798, 0.05801)	(0.04715, 0.03123, 0.99840)	(0.04715, 0.03123, 0.99840, −0.07094)
	Se (21–28)	(0.97462, 1.08838, −2.96809)	(−0.01127, −0.03216, 0.99942)	(−0.01127, −0.03216, 0.99942, 3.01235)

---

Note that when starting the optimization of InSe₃ at a structure close to the optimized ones in Fig. 4, distinct results are predicted as shown in there. In fact, similar structures are also predicted with B3LYP with both STO-3G and the minimal DZP and ADZP basis sets. In turn, when starting at the CCSD(T)/ATZP geometry in panel (b) of Fig. 3, the optimized CCSD(T)/ADZP one is predicted to be a saddle point close to the ATZP one but with an imaginary harmonic vibrational frequency of −26.6 cm⁻¹. So, except for the nature of the stationary point (saddle point vs. minimum), the agreement is generally good with other structural parameters at the B3LYP/STO-3G level. With the latter, only a slight contraction is observed in the In–Se bonds while the Se–Se ones are slightly elongated. Conversely, the InSe bonds are slightly elongated and the Se–Se bonds shortened when using ADZP. Such variations are likely responsible for those in the harmonic vibrational frequencies at this level of theory. Upon inclusion of relativistic effects, the optimized stationary point becomes a third-order saddle point (index 3) at the ADZP level, with the dominant imaginary mode corresponding to out-of-plane folding and the remaining modes to in-plane distortion. Suffice it to add that similar calculations were also done with the M06-2X functional,⁶⁸ with the results generally similar to B3LYP.

Fig. 4 Optimized structures of InSe₃ and harmonic vibrational frequencies (in red) at a Y-shaped stationary point of InSe₃. All are planar, with energies in E_h and bond lengths in Å.

Moreover, relativistic B3LYP and CCSD(T) calculations have been performed at the stationary points of InSe₃ in Fig. 4 by employing both conventional and relativistic DZP-DKH basis sets (to our knowledge, the ADZP-DKH basis sets are unavailable, while the standard STO-3G basis set is inherently non-relativistic). The results agree in showing them all to be planar, but with differences concerning the nature of the stationary point (minimum vs. saddle point), hence distinct harmonic vibrational frequencies. These are shown to depend on both the basis set and method, but have no implications (all are planar stationary points) as far as the topic of the present work is concerned. With STO-3G, a slight contraction of In–Se bond distances is observed, which may partly explain the somewhat larger harmonic vibrational frequency at this level of theory. Since the unrestricted variant of DFT has been adopted, care was taken to check whether 〈S²〉 deviated less than up to very few percent from the desired value. Interestingly, the calculated B3LYP/STO-3G properties are found in good agreement with the rather cost-expensive CCSD(T) ones using both DZP and DZP-DKH basis sets.

5.2 2D In₂Se₃

Fig. 5 illustrates the optimized structures of the pseudo-atom In₂Se₃. It is a singlet-state structure, while the two molecular fragments into which it can be split are radicals, namely InSe and InSe₂, hence reactive species. Although a bulk InSe exists, it does not warrant a stable free diatomic molecule in the gas phase under ambient conditions. Known as an inorganic compound, it is here predicted at the B3LYP/STO-3G level with an equilibrium bond distance of 2.117 Å and harmonic vibrational frequency of 413.8 cm⁻¹ (the corresponding B3LYP/ATZP values are 2.330 Å and 270.1 cm⁻¹, while at the CCSD(T)/ATZP level they are 2.673 Å and 222.2 cm⁻¹); no published spectroscopic constants were found for the gas-phase InSe diatomic in Huber & Herzberg's constants of diatomic molecules. Its absence across this primary source strongly suggests that a free InSe diatomic has not been spectroscopically measured in the gas phase: if it forms at all (during high-T vaporization of InSe solids, laser ablation or plasmas), it is likely transient/low-abundance, and it has not received spectroscopic characterization. The present results predict it to be stable (i.e., a minimum). This is also the case when including relativistic effects and the ATZP basis set (ATZP-DKH is unavailable): E = −8274.70560991E_h with equilibrium geometry 2.330 Å, which is ΔE = 63.512 kcal mol⁻¹ lower than the quartet at its equilibrium geometry and same level of theory. A similar situation holds with DTZ-DKH basis sets, with the above values now reading E = −8305.65695989E_h at a slightly larger equilibrium distance 2.342 Å, and ΔE = 70.722 kcal mol⁻¹, in the same order. Regarding In₂, suffice it to add that B3LYP/STO-3G relativistic calculations have also been done for its singlet and triplet states, with the triplet still found to be lower in energy by Δ = 18.996 kcal mol⁻¹. The ATZP basis set is then the main reason for the singlet to be the lowest in energy; see Table 1.

Fig. 5 Optimized DFT structures (all are planar) for pseudo-atom In₂Se₃ using the B3LYP [(a) STO-3G; (b) 3-21G; (c) ADZP; (d) ATZP] and M06-2X [(e), ATZP] functionals. The bluish P implies planarity. Also indicated are the spin values 〈S²〉, dipole moment (in D), and energy (in E_h).

Notable from Fig. 5 is the fact that the predicted structures for the pseudo-atom agree well with each other when studied with both subminimal and extended basis sets. To test the functionals’ role, the geometry of In₂Se₃ has also been optimized with M06-2X and the ATZP basis set. The result in panel (e) agrees nicely with that from B3LYP in (d).

Interestingly, but perhaps not surprisingly, when seen top-down from Fig. 5, one observes Se–In–Se–In–Se, which one expects to be replicated in the sublayers of the 2D In₂Se₃ material. If In₂Se₃ is assumed to be a pseudo-atom, its center of mass (CoM) is somewhere inside the triangle formed by the Se atom number 2 and the In atoms 1 and 4. Given the near C_2v optimized symmetry of In₂Se₃, and the disparity of nuclear masses of the Se (78.97 a.m.u.) and indium (114.818 a.m.u.) atoms, we may for simplicity assume the CoM of the pseudo-atom to be close to the middle of the two In atoms. If so and the Lemma applies, the CoM of the four pseudo-atoms when forming the (In₂Se₃)₄ pseudo-tile should approximately sit in a plane. In fact, a similar simplified approach has been adopted¹³ in the study of the trilayer 2D MoS₂ material with good success.

Recalling next that a plane can be described by a point in the plane and a vector (normal) perpendicular to it, one can report the plane by its centroid and unit normal. This is particularly useful when the plane is an approximation of a set of data points, since the centroid represents the data's central location. Using this, the top In atoms (#s 5–8) in the pseudo-tile Se₈In₁₂ [shown in panel (c) of Fig. 6] have the centroid, normal, and best plane fit indicated in the top entry of Table 4. Suffice it to say that the fit has rmsd = 0.022 Å. Correspondingly, one gets for the bottom In atoms (#s 13–16) the results given in the second entry of Table 4, with rmsd = 0.132 Å. The angle between the above two planes is 69.7 deg, which shows that they are not parallel. As for the top file of Se atoms (#s 1–4) in the third row, rmsd = 0.330 Å, while for Se atoms 9–12, rmsd = 0.197 Å. Similarly, for Se atoms 17–20, the fit shows rmsd = 0.322 Å. The corresponding angles between the Se sublayer planes are 3.3 deg (top–middle), 5.4 deg (top–bottom), and 2.8 deg (middle–bottom). Hence, the three Se planes are mutually nearly parallel (within ∼5–6 deg). Relative to the In planes, for the top In vs. top/middle/bottom Se, they are 22.3 deg/19.9 deg/24.4 deg. In turn, for the bottom In vs. top/middle/bottom Se planes, the angles are 78.9 deg/79.6 deg/74.6 deg. Hence, the top In is moderately tilted relative to Se planes (∼20–24 deg), while the bottom In is ∼75–80 deg off (highly tilted).

Fig. 6 Two convenient views of the optimized (In₂Se₃)_x clusters up to the tetra-pseudo-atom tile: (a) x = 2; (b) x = 3; and (c) x = 4. Shown by red points on the right-hand-side of panel (c) are the CoM of the various pseudo-tiles in (In₂Se₃)₄ [involving In atom-pairs (5,13), (6,14), (7,15), and (8,16)], while the plane centroid of the four CoM and normal tip with 1.5 Å × n̂ are shown in white.

Regarding the center-of-charge (CoC) planes (based on the oxidation-state definition), since the net charge in Se₈In₁₂ is zero [8 × (+3) + 12 × (−2) = 0], a single global CoC is undefined. However, the separate CoC and principal planes for the cation sublattice (In³⁺) and anion sublattice (Se²⁻) can be computed. For the In³⁺ (cation) charge plane, the CoC, normal, and plane equation are shown in the first row of Table 5, with the fit showing a rmsd = 0.900 Å. In turn, the angles to the In planes are 16.1 deg (top In atoms) and 65.3 deg (bottom In). Similarly, the Se²⁻ (anion) charge plane has center of charge, normal n̂ and plane as defined in the second row of Table 5, with rmsd = 1.620 Å. Similarly, the angles to the Se planes are 27.3 deg (top Se atoms), 15.0 deg (middle Se), and 17.0 deg (bottom Se). Moreover, the angle between the cation and anion charge planes is 23.8 deg. So, the Se sublayers are nearly parallel to one another (≲6 deg spread). Conversely, the In sublayers are strongly non-parallel (≃70 deg apart); the upper In plane aligns moderately with the Se planes, while the lower In plane is almost orthogonal. In turn, the oxidation-state charge planes (In³⁺ vs. Se²⁻) are themselves tilted by ∼24 deg, and the In³⁺ charge plane tracks closely the top In plane (≃16 deg), whereas the Se²⁻ charge plane sits ∼15–27 deg from the Se planes.

Table 5 Center of charge (CoC) planes

Pseudo-tile	Atoms #s	Centroid^a	Normal^b	(a,b,c,d) in^c ax + by + cz + d = 0
a x, y, z coordinates of the centroid.b x, y, z coordinates of the normal.c Parameters in equation ax + by + cz + d = 0 of fitted plane; see also the text.
In8Se12	In^3+	(0.15735, 0.05083, −0.06091)	(0.15467, 0.94325, 0.29387)	(0.41730, 0.84967, 0.32237, 0.1893)
	Se^2−	(−0.02991, 0.09328, −0.18580)	(0.02222, 0.76608, 0.64237)	(0.02222, 0.76608, 0.64237, 0.04855)
In8Se20	In^3+	(−0.00300, −0.26455, −0.13178)	(−0.12538, 0.83930, 0.52902)	(−0.12538, 0.83930, 0.52902, 0.29137)
	Se^2−	(0.00120, 0.10582, 0.05271)	(0.70535, −0.65084, 0.28088)	(0.70535, −0.65084, 0.28088, 0.05322)

---

5.3 2D In₂Se₅

As in the previous case, Fig. 7 illustrates the optimized structures of pseudo-atom In₂Se₅. Note that it is a singlet-state, while the two smaller pseudo-atom units that make it are radicals, namely InSe₂ and InSe₃, hence reactive species. Once more, by top-down observation of pseudo-atom In₂Se₅ one wonders whether all layers have an equal number of Se and In atoms. Interestingly, one observes that the final optimized structures with the subminimal STO-3G and extended ATZP basis sets resemble much each other, offering support for employing hereinafter the affordable STO-3G basis set in all calculations.

Fig. 7 Optimized DFT structures for the pseudo-atom In₂Se₅ in the corresponding 2D material using the B3LYP [(a) STO-3G; (b) 3-21G; (c) 3-21G; (d) ATZP]) and M06-2X [(e) ATZP] functionals. The bluish N and P have the meaning assigned in Fig. 5, while the asterisk in N* implies mildly nonplanar. Also indicated at the bottom are the expectation values of 〈S²〉, dipole moment (in D), and energy (in E_h).

The corresponding CoM and CoC fitted planes to In and Se are at the bottom of Tables 4 and 5 and discussed next. Regarding the rmsd of the fits to In atoms #s 9–12 and 17–20 of In₈Se₂₀ [panel (c) of Fig. 8, a saddle point at this level of theory], one has 0.418 Å and 0.133 Å, respectively. Additionally, for the top (#s 1–8), middle (13–16), and bottom (21–28) Se planes, the rmsd are 0.611 Å, 0.007 Å, and 0.590 Å, respectively. Similarly, for the In³⁺ and Se²⁻ planes, the rmsd are 1.464 Å and 2.120 Å, respectively. In turn, the angles between the In vs. Se planes are for top In vs. top/middle/bottom Se 11.979 deg/5.647 deg/2.500 deg, while for the bottom In vs. top/middle/bottom Se they are 27.946 deg/12.691 deg/13.958 deg. Moreover, the angles of the In³⁺ plane to the In planes are top In 61.814 deg and bottom In 46.369 deg. Similarly, the angles of the Se²⁻ plane to the Se planes are top 61.460 deg, middle 72.941 deg, and bottom 72.921 deg.

Fig. 8 Two convenient views of the optimized clusters (In₂Se₅)_x up to the tetra-pseudo-atom tile: (a) x = 2; (b) x = 3; (c) x = 4 [involving In atom-pairs (9,17), (10,18), (11,19), (12,20)]. Shown by red points on the right-hand-side of panel (c) are the CoM of the various pseudo-tiles in (In₂Se₅)₄.

In summary, the two In planes are moderately tilted relative to each other (∼16 deg), unlike the ∼70 deg tilt observed in In₈Se₁₂. Moreover, the top In layer is nearly aligned with the bottom Se plane (only 2.5 deg off) and close to the middle Se plane (5.6 deg), thus suggesting a more symmetric stacking near the top side. Additionally, the oxidation-based charge planes are quite oblique to their respective atomic planes (≳60 deg), capturing the overall skew of the cation vs. anion frameworks across the 5 sublayers.

As noted above, with the exception of non-round cases, the centers of mass of the pseudo-atoms have been largely neglected. This may explain that one In plane usually shows a small rmsd in the fit, while the other is considerably larger. A better approach would be assuming it to be located at the midcenter of both In atoms of the pseudo-atom. The various mid-points so obtained are summarized in Table 6 for both In₈Se₁₂ and In₈Se₂₀. As before, the least-squares plane in In₈Se₁₂ through the above four average points can be defined by its centroid (0.15735, 0.05083, −0.0609) Å and normal n̂ = (0.07705, 0.08610, 0.99330), with the equation 0.07705x + 0.08610y + 0.99330z + 0.04400 = 0, and rmsd = 0.0121 Å. Clearly, there is excellent coplanarity. In turn, for In₈Se₂₀, the least squares fit plane has centroid (−0.00300, −0.26455, −0.13178) Å and normal n̂ = (−0.04978, 0.07881, 0.99565), having as equation −0.04978x + 0.07881y + 0.99565z + 0.15191 = 0, and rmsd = 0.1490 Å. Note that this “midpoint plane” lies close to the In³⁺ CoC plane, even sharing the same centroid to within the employed precision, and its normal is nearly along the +z axis (tilt only a few degrees), hence consistent with a mid-Se-like reference surface. Not surprisingly, there is good coplanarity but not as good as in the previous case.

Table 6 Center of charge (CoC) planes

Pseudo-tile	Atom-pair^a	Center (atom-pair)^b
a Averaged atom pair.b x, y, z centroid of the averaged pair.
In8Se12	(9,17)	(−4.60635, 1.25501, 0.20833)
	(10,18)	(−1.43983, −0.05806, 0.05782)
	(11,19)	(1.85114, −0.52701, −0.12452)
	(12,20)	(4.82443, −0.46660, −0.38527)
In8Se20	(9,17)	(−3.16090, −1.45857, −0.05025)
	(10,18)	(1.32916, −2.33308, −0.02736)
	(11,19)	(−1.05990, 1.51936, −0.49733)
	(12,20)	(2.87966, 1.21410, 0.04782)

---

One may possibly get a further improvement by employing the actual CoM of the various pseudo-atoms. This has also been done, with rmsd ≃ 0.067 Å which shows the four CoM of the In₂Se₅ units to be highly coplanar, markedly better than using just the In–In midpoints (thus ignoring the Se mass distribution). Indeed, the plane's normal is ∼12 deg away from the +z direction (since n̂·ẑ ∼ 0.981): the CoM plane is nearly horizontal relative to the employed coordinate frame. This matches the physical expectation: including the five nearby Se atoms for each unit balances the cation–anion mass about a local median layer, yielding a more coherent structural plane for the In₂Se₅ motifs.

The above corroborates that the four pseudo-atoms are coplanar when represented by their CoM. Although it also provides evidence supporting the PBC model, it cannot warrant that all In and Se planes lie parallel to each other as commonly assumed. Hence, one wonders whether keeping the CoM of the pseudo-tiles (or generalized pseudo-tiles) in a plane may be more realistic than the full planarity assumed in the PBC model. Of course, a proof of this would require considering more pseudo-tiles, a task beyond our goal in the present work. Nevertheless, it suggests that further flexibility may be deemed for the PBC model by allowing the orientation of the various fundamental building blocks (here assumed as tiles) to be optimized within the employed unit cell while keeping in a fixed plane their CoM. Testing whether it is so is out of the scope of the present work.

5.4 Evidence of quasi-planarity beyond tiles

The structural analysis was extended to reasonably large clusters by considering In₁₂S₁₈ and In₁₀Se₂₅ (see Fig. 8), which have been optimized for probing reasons at significant computational cost (weeks of uninterrupted computational work), when including the calculation of harmonic vibrational frequencies. These species exhibit slight out-of-plane distortions at their edges, but the interior In atoms remain generally well-aligned.

The quasi-optimized structure of In₁₂Se₁₈ is shown in panel (a) of Fig. 8, jointly with a visual plot of the six local centroids (one for each Se–In–Se–In–Se monomer stack). In turn, the best fit planes to the five sublayers are in the top part of Table 7. From the top Se plane to In, Se, In, and the bottom Se one, the rmsd were 1.165, 0.979, 0.746, 0.105, and 0.655 Å, respectively.

Table 7 Centers of mass planes which integrate the cluster

Cluster	Centroid
In12Se18	(−1.20139, −0.39519, 2.99952)
	(−1.24652, −0.63122, 1.08432)
	(0.24068, −0.77444, 0.86537)
	(1.03230, 0.33229, −1.03381)
	(0.92204, 0.13880, −2.31114)
In10Se25	(−3.15319, −1.94958, 2.08261)
	(−0.0900, −1.15781, −0.40223)
	(3.23413, −2.12588, −1.48170)
	(−2.6925, 1.08826, 0.26372)
	(1.70522, 1.58715, −0.27822)

---

In turn, a side-by-side detailed analysis of the In₂Se₅ CoM plane for the In₁₀Se₂₅ structure is now given. The CoM signed distances of the five In₂Se₅ pseudo-atoms are represented in Fig. 9 (marked for simplicity with “oxygen-atom” red dots, which should not be confused with true oxygen atoms). The best-fit plane to these CoM has centroid (−0.19927, −0.51157, 0.03683) Å, unit normal n̂ = (0.39692, 0.09349, 0.91308), and equation 0.396919x + 0.093490y + 0.913080z + 0.093287 = 0, with rmsd = 0.5219 Å. The signed distances of each CoM to the plane reveal a gentle “saddle” across the five units: one sits highest relative to the plane, another lowest, with the remaining three lying close to the plane. Regarding the reference planes, the CoM plane is closest to the Se middle/bottom planes (∼7 deg), indicating that the mass-weighted motif centers follow the chalcogen framework more closely than the individual In sublayers. The top and bottom In planes are nearly parallel (∼1.5 deg), as expected for a layered motif with mild overall tilt. The moderate CoM-plane rmsd (0.52 Å) over 5 units suggests mild long-wavelength puckering, not disorder.

Fig. 9 Optimized clusters up to penta- and hexa-pseudo-atomic tiles: (a) (In₂Se₃)₆; (b) (In₂Se₅)₅. In both cases, 〈s²〉 = 0.000 with the energy being excessively large (not quoted by default in MOLPRO).

A further remark to note is that all inter-sublayer contacts cluster in the 2.33–2.92 Å window, with means 2.54–2.66 Å consistent with typical In–Se bond lengths in layered indium selenides. The middle → 4th and top → 2nd sets show the widest spread (σ ≃ 0.18 deg), hinting at mild puckering/tilt across those interfaces.

Finally, it should be recalled that the stationary points reported here were obtained with GTO basis sets, while the popular Vienna ab initio simulation package^69,70 (VASP) or Quantum Expresso⁷¹ programs commonly employed in studying materials are plane-wave codes (an exception is TURBOMOLE^72,73 which employs periodic DFT with Gaussian basis functions). Hence, only the results of the former for a molecule placed in a large periodic cell with sufficient vacuum space to prevent interactions with its periodic images are comparable to those obtained with GTO-basis quantum chemistry software, particularly when using small basis sets.⁷⁴ So, the fair agreement for planarity of the title 2D materials provides some support to the assumed planarity in PBC given the disparity of employed methodologies.

5.5 Bond lengths in 2D In₂Se₃ and In₂Se₅ materials

Regarding typical In–Se bond lengths in layered 2D In₂Se₃, existing information points to 2.55–2.90 Å, depending on phase (α, β, γ, δ) and strain conditions.⁷⁵ A spread of In–Se bond distances in the range of 2.60–2.68 Å has also been reported,⁷⁶ while another suggests⁷⁶ 2.51–2.96 Å. The above three spreads are displayed for illustration in Fig. 10.

Fig. 10 In–Se bond distances in (In₂Se₃)_n, n = 1–6, compared with existing data from other sources. The various symbols represent the following: solid dots, Se (top sublayer)–In (2nd); open circles, In (2nd)–Se (middle); solid up-triangles, Se (middle)–In (4th); and solid down-triangles, In (4th)–Se (bottom). The dashed straight line in dark-green illustrates the growth pattern of the average In–Se bond distances with cluster size. See the text.

The author is not aware of any literature source for the inter-sublayer In–Se distances in 2D In₂Se₅ (penta-sublayer) structure which one could compare with the values here estimated. For example, by assuming the layering to match the z-ordering in our geometry, for In₈Se₂₀, the top Se atoms are #s 1–8, the In layer #s 9–12, the middle Se layer #s 13–16, the fourth In layer #s 17–20, and the bottom Se layer #s 21–28. With a bond cutoff of 3.0 Å, which is a conservative In–Se bonding distance (typical ∼ 2.6–2.8 Å), one finds the results evinced in Fig. 11. Given the currently optimized geometry for the title In₈Se₂₀ and the distances reported in Fig. 11, many of the present cross-layer In–Se distances are ∼2.6–2.7 Å. Hence, the tile is already “quasi-2D”, and one may expect that in the ideal infinite 2D limit those distances will be similar (perhaps slightly shorter if relaxation is unconstrained).

Fig. 11 In–Se bond distances in (In₂Se₅)_n, n = 1–5, compared with data from other sources. The symbols and line as in Fig. 10. See the text.

If the focus is on the cross-sublayer Se–In distances inside a quintuple-like stack, the best grounded range (based on InSe and α-In₂Se₃ monolayers) is Se (top) → In (2nd): ∼2.6–2.8 Å; In (2nd) → Se (middle): ∼2.6–2.9 Å (often the longest in α-In₂Se₃); Se (middle) → In (4th): ∼2.6–2.9 Å; and In (4th) → Se (bottom): ∼2.6–2.9 Å. Accordingly, ∼2.5–2.9 Å is a safe envelope for realistic, relaxed 2D indium–selenide quintuple-layer environments, with most bonds clustering near ∼2.65–2.80 Å (method/strain dependent). The exact split between the two inequivalent In–Se bonds (short vs. long) is polymorph- and relaxation-dependent, as the α-In₂Se₃ references show.⁷⁵

Rather than histogramming the data as in previous work,¹³ the geometric data have been analyzed by fitting the average of the In–Se distances for each cluster size to the sigmoidal function

d_MoS = a + b tanh(γn), (1)

where n is the cluster size. While the data for In₂Se₃ show a near-linear trend (a = 2.206, b = 0.495, γ = 0.147), for In₂Se₅ they show a rapidly varying pattern (a = −223 228, b = 223 231, γ = 7.17) that early attains the asymptote; in both cases, d outcomes are in Å. Clearly, Fig. 10 and 11 show that the calculated data match well the experimentally observed In–Se bond lengths in the corresponding 2D materials. Hence, as already noted, this suggests that the local bonding environments in the studied finite clusters evolve smoothly and fast toward the bulk-like behavior of the extended penta-sublayer sheets.

6 Concluding remarks

This work presents a quasi-molecular interpretation of the pseudo-planarity observed in 2D indium selenide compounds by modeling those materials as tilings of structurally meaningful molecular fragments made of pseudo-atomic-like In₂Se₃ and In₂Se₅ tiles, respectively In₈Se₁₂ and In₈Se₂₀, and even larger fragments. Through systematic DFT calculations, it was demonstrated that small clusters, although finite and only approximately planar, preserve the essential features of the extended 2D lattice. So, even if in striking contrast with the common view of 2D indium selenides as extended 2D crystals (or lattices) of In and Se atoms, the proposed perspective can be useful.

Notably, the novel results show that the structural and electronic identities of In₂Se₃ and In₂Se₅ are not entirely lost upon transitioning from the (gas phase) molecules to the monolayers of the corresponding 2D materials, themselves made up of five sub-layers: Se–In–Se–In–Se. The extended crystals are viewed as coherent networks of quasi-planar pseudo-tiles, whose collective arrangement gives rise to quasi-planarity, commonly assumed in PBC simulations. The present quasi-molecular perspective then bridges the gap between molecular and solid-state approaches, offering a complementary scheme to understand bonding, defect chemistry, and edge reactivity in 2D materials.

The present observations further validate to a large extent the underlying premise of the quasi-molecular tiling approach: local pseudo-planar geometries can collectively yield an extended, nearly planar structure. But, do such small deviations of planarity as suggested by the tilts of the In and Se planes fade away in the 2D In₂Se₃ and In₂Se₅ materials? If not, it would be valuable to test the implications (differences) relative to the popular PBC approach, a topic beyond the goal of the present study.

In conclusion, this study builds on our earlier proposal for 2D MoS₂ and demonstrates that the underlying conceptual framework is both insightful and broadly applicable to other 2D materials. Besides, by reporting novel work on the challenging 2D In₂Se₅ material, it also opens pathways for future work, including the development of scalable algorithms to systematically generate and assess tiling units, investigate how defects and doping influence local planarity, and include relativistic and dispersion effects in large-scale cluster calculations. Finally, recalling that atoms are not necessarily spherical, and even less pointwise, the Lemmas may be applicable to pseudo-atom-like^12,13 ones. One then wonders whether the approach could be used to explain the novel class of cyclic conjugated molecules composed of annelated n-membered rings^77,78 or even the aimed⁷⁹ free-standing planar 2D crystal of Sr₃Cr₂O₇, since this is known for its inherent 2D nature within the bulk structure. Even though significantly more cost-expensive, and apparently of no impact on the results reported in the present work, a comprehensive relativistic analysis of the title systems would also be valuable.

Conflicts of interest

The author declares no conflicts of interest.

Data availability

Data supporting this study are included within the article and/or supplementary information (SI). See DOI: https://doi.org/10.1039/d6cp00070c.

Acknowledgements

The author gratefully acknowledges support from the Double Hundred Experts Project of Shandong Province, China (Grant No. WSP2023008), the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Brazil – Finance Code 001, and the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Brazil. Support from the Chemistry Centre (CQC-IMS) at the University of Coimbra, Portugal, is also acknowledged, under the projects UIDB/00313/2025, UIDP/00313/2025, and LA/P/0056/2020, funded by the Fundação para a Ciência e a Tecnologia (FCT), Portugal.

Notes and references

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Electric Field Effect in Atomically Thin Carbon Films, Science, 2004, 306, 666–669.
2. J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth and S. Roth, The structure of suspended graphene sheets, Nature, 2007, 446, 60–63.
3. K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V. Khotkevich, S. V. Morozov and A. K. Geim, Two-dimensional atomic crystals, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 10451–10453.
4. K. F. Mak, C. Lee, J. Hone, J. Shan and T. F. Heinz, Atomically Thin MoS₂: A New Direct-Gap Semiconductor, Phys. Rev. Lett., 2010, 105, 136805.
5. L. Debbichi, O. Eriksson and S. Lebègue, Two-Dimensional Indium Selenides Compounds: An Ab Initio Study, J. Phys. Chem. Lett., 2015, 6, 3098–3103.
6. W. Ding, J. Zhu, Z. Wang, Y. Gao, D. Xiao, Y. Gu, Z. Zhang and W. Zhu, Prediction of intrinsic two-dimensional ferroelectrics in In2Se3 and other III2-VI3 van der Waals materials, Nat. Commun., 2017, 8, 14956.
7. M. Soleimani and M. Pourfath, Ferroelectricity and phase transitions in In2Se3 van der Waals material, Nanoscale, 2020, 12, 22688–22697.
8. X. Zheng, W. Han, K. Yang, L. W. Wong, C. S. Tsang, K. H. Lai, F. Zheng, T. Yang, S. P. Lau, T. H. Ly, M. Yang and J. Zhao, Phase and polarization modulation in two-dimensional In₂Se₃ via in situ transmission electron microscopy, Sci. Adv., 2022, 8, eabo0773.
9. Chemical Vapor Deposition: Principles and Applications, ed. M. L. Hitchman and K. F. Jensen, Academic Press, London, San Diego, 1993.
10. C. Paorici and G. Attolini, Vapour growth of bulk crystals by PVT and CVT, Prog. Cryst. Growth Charact. Mater., 2004, 48–49, 2–41.
11. D. W. Boukhvalov, B. Gürbulak, S. Duman, L. Wang, A. Politano, L. S. Caputi, G. Chiarello, A. Cupolillo, D. W. Boukhvalov, B. Gürbulak, S. Duman, L. Wang, A. Politano, L. S. Caputi, G. Chiarello and A. Cupolillo, The Advent of Indium Selenide: Synthesis, Electronic Properties, Ambient Stability and Applications, Nanomaterials, 2017, 7, 372.
12. A. J. C. Varandas, Calcium Carbides: Can Hexa-carbides Grow Unlimitedly? Theoretical Perspective and Issues That Oppose a Definite Answer, Phys. Chem. Chem. Phys., 2026, 28, 174–187.
13. A. J. C. Varandas, Quasi-Molecular Perspective on the Planarity of 2D Materials: The case of 2D MoS₂, Phys. Chem. Chem. Phys., 2026, 28, 5992–6005.
14. A. J. C. Varandas, From six to eight Pi-electron bare rings of group- XIV elements and beyond: can planarity be deciphered from the “quasi-molecules” they embed?, Phys. Chem. Chem. Phys., 2022, 24(2022), 8488–8507.
15. A. J. C. Varandas, Can the quasi-molecule concept help in deciphering planarity? The case of polycyclic aromatic hydrocarbons, Int. J. Quantum Chem., 2023, 123, e27036.
16. A. J. C. Varandas, Carbon-[n]Triangulenes and Sila-[n]Triangulenes: Which Are Planar?, J. Phys. Chem. A, 2023, 127, 5048–5064.
17. A. J. C. Varandas, Cn, CnH and their anions: Quest for linearity with n≤ 8 even versus odd, and beyond, Int. J. Quantum Chem., 2024, 124, e27287.
18. J. Manz, R. Meyer, E. Pollak and J. Römelt, A new possibility of chemical bonding: vibrational stabilization of IHI, Chem. Phys. Lett., 1982, 93, 184–187.
19. J. Manz, R. Meyer and J. Römelt, On vibrational bonding of IHI, Chem. Phys. Lett., 1983, 96, 607–612.
20. J. N. Murrell, S. Carter, S. C. Farantos, P. Huxley and A. J. C. Varandas, Molecular Potential Energy Functions, Wiley, Chichester, 1984.
21. A. J. C. Varandas, Intermolecular and Intramolecular Potentials: Topographical Aspects, Calculation, and Functional Representation via A Double Many-Body Expansion Method, in Advances in Chemical Physics, ed. I. Prigogine and S. A. Rice, Wiley, 1st edn, 1988, vol. 74, pp. 255–338.
22. G. Herzberg, Molecular Spectra and Molecular Structure. III Electronic Spectra and Electronic Structure of Polyatomic Molecules, Van Nostrand, New York, 1966.
23. T. Helgaker, P. Jørgensen and J. Olsen, Molecular Electronic-Structure Theory, Wiley, Chichester, 2000.
24. A. C. West, M. W. Schmidt, M. S. Gordon and K. Ruedenberg, A comprehensive analysis of molecule-intrinsic quasi-atomic, bonding, and correlating orbitals. I. Hartree-Fock wave functions, J. Chem. Phys., 2013, 139, 234107.
25. A. C. West, M. W. Schmidt, M. S. Gordon and K. Ruedenberg, A Comprehensive Analysis in Terms of Molecule-Intrinsic, Quasi-Atomic Orbitals. II. Strongly Correlated MCSCF Wave Functions, J. Phys. Chem. A, 2015, 119, 10360–10367.
26. A. C. West, M. W. Schmidt, M. S. Gordon and K. Ruedenberg, Intrinsic Resolution of Molecular Electronic Wave Functions and Energies in Terms of Quasi-atoms and Their Interactions, J. Phys. Chem. A, 2017, 121, 1086–1105.
27. E. B. Guidez, M. S. Gordon and K. Ruedenberg, Why is Si₂H₂ Not Linear? An Intrinsic Quasi-Atomic Bonding Analysis, J. Am. Chem. Soc., 2020, 142, 13729–13742.
28. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553.
29. A. J. C. Varandas, Extrapolation to the complete basis set limit without counterpoise. the pair potential of helium revisited, J. Phys. Chem. A, 2010, 114(2010), 8505–8516.
30. Ł. M. Mentel and E. J. Baerends, Can the Counterpoise Correction for Basis Set Superposition Effect Be Justified?, J. Chem. Theory Comput., 2014, 10, 252–267.
31. L. A. Burns, M. S. Marshall and C. D. Sherrill, Comparing Counterpoise-Corrected, Uncorrected, and Averaged Binding Energies for Benchmarking Noncovalent Interactions, J. Chem. Theory Comput., 2014, 10, 49–57.
32. A. J. C. Varandas, CBS extrapolation in electronic structure pushed to the end: A revival of minimal and subminimal basis sets, Phys. Chem. Chem. Phys., 2018, 20, 22084–22098.
33. A. J. C. Varandas, Straightening the Hierarchical Staircase for Basis Set Extrapolations: A Low-Cost Approach to High-Accuracy Computational Chemistry, Annu. Rev. Phys. Chem., 2018, 69, 177–203.
34. A. J. C. Varandas, CBS extrapolation of Hartree-Fock energy: Pople and Dunning basis sets hand-to-hand on the endeavour, Phys. Chem. Chem. Phys., 2019, 21, 8022–8034.
35. A. J. C. Varandas, Extrapolation in quantum chemistry: Insights on energetics and reaction dynamics, J. Theor. Comput. Chem., 2020, 19, 2030001.
36. A. J. C. Varandas, Canonical and explicitly-correlated coupled cluster correlation energies of sub-kJ mol-1 accuracy via cost-effective hybrid-post-CBS extrapolation, Phys. Chem. Chem. Phys., 2021, 23, 9571–9584.
37. A. J. C. Varandas, Scale-free-modeling (harmonic) vibrational frequencies: Assessing accuracy and cost-effectiveness by CBS extrapolation, J. Chem. Phys., 2022, 157, 174110.
38. W. J. Hehre, R. Ditchfield and J. A. Pople, Self-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules, J. Chem. Phys., 1972, 56, 2257–2261.
39. J. S. Binkley and J. A. Pople, Møller–Plesset theory for atomic ground state energies, Int. J. Quantum Chem., 1975, 9, 229–236.
40. E. R. Davidson and D. Feller, Basis set selection for molecular calculations, Chem. Rev., 1986, 86, 681–696.
41. P. Pulay, Ab initio calculation of force constants and equilibrium geometries in polyatomic molecules: I. Theory, Mol. Phys., 1969, 17, 197–204.
42. S. F. Boys, Construction of Some Molecular Orbitals to Be Approximately Invariant for Changes from One Molecule to Another, Rev. Mod. Phys., 1960, 32, 296–299.
43. R. G. Woolley and B. T. Sutcliffe, Molecular structure and the Born–Oppenheimer approximation, Chem. Phys. Lett., 1977, 45, 393–398.
44. B. T. Sutcliffe and R. G. Woolley, On the quantum theory of molecules, J. Chem. Phys., 2012, 137(22), 22A544.
45. M. Born and R. Oppenheimer, Zur Quantentheorie der Molekeln, Ann. Phys., 1927, 389, 457–484.
46. A. J. C. Varandas, J. da Providência and J. P. da Providência, Binding of muonated hydrogen molecules and Born-Oppenheimer approximation revisited, Can. J. Phys., 2020, 98, 379–384.
47. T. H. Dunning, Jr., Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen, J. Chem. Phys., 1989, 90, 1007–1023.
48. T. Li, Y. Wang, W. Li, D. Mao, C. J. Benmore, I. Evangelista, H. Xing, Q. Li, F. Wang, G. Sivaraman, A. Janotti, S. Law and T. Gu, Structural Phase Transitions between Layered Indium Selenide for Integrated Photonic Memory, Adv. Mater., 2022, 34, 2108261.
49. C. Bannwarth, S. Ehlert and S. Grimme, GFN2-xTB–An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions, J. Chem. Theory Comput., 2019, 15, 1652–1671.
50. D. E. Woon and T. H. Dunning, Jr., Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon, J. Chem. Phys., 1993, 98, 1358–1371.
51. A. J. C. Varandas, Tipping point to explain melamine's non-planarity or otherwise. Can quasi-molecules be “seen” in the parent molecule?, Chem. Phys. Lett., 2024, 852, 141493.
52. A. J. C. Varandas, A re-examination of the highly floppy melamine molecule and predictability of quasi-molecule theory, Chem. Phys. Lett., 2025, 873, 142137.
53. T. Mondal, A. Guerra-Barroso, J. Fang, J. Li and A. J. C. Varandas, What have the XH4+ (X = C, Si, Ge) ions in common? An updated summary, J. Chem. Phys., 2025, 162, 114309.
54. H. J. Kulik, N. Luehr, I. S. Ufimtsev and T. J. Martinez, Ab Initio Quantum Chemistry for Protein Structures, J. Phys. Chem. B, 2012, 116, 12501–12509.
55. W. J. Hehre, R. F. Stewart and J. A. Pople, Self-Consistent Molecular-Orbital Methods. I. Use of Gaussian Expansions of Slater-Type Atomic Orbitals, J. Chem. Phys., 1969, 51, 2657–2664.
56. K. Kapusta, E. Voronkov, S. Okovytyy, V. Korobov and J. Leszczynski, Reconstruction of STO-3G Family Basis Set for the Accurate Calculation of Magnetic Properties, Russ. J. Phys. Chem., 2018, 92, 2827–2834.
57. E. Keller, J. Morgenstein, K. Reuter and J. T. Margraf, Small basis set density functional theory method for cost-efficient, large-scale condensed matter simulations, J. Chem. Phys., 2024, 161, 074104.
58. R. Sarkar, Š. Budzák and D. Jacquemin, On the Search of a Modified STO-3G Basis Set Optimized for Molecules, Int. J. Quantum Chem., 2025, 125, e70125.
59. R. Pasquier, M. Graml and J. Wilhelm, Gaussian Basis Sets for All-Electron Excited-State Calculations of Large Molecules, J. Chem. Theory Comput., 2026, 22, 540–557.
60. C. Bannwarth, S. Ehlert and S. Grimme, GFN2-xTB–An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions, J. Chem. Theory Comput., 2019, 15, 1652–1671.
61. K. Kooser, D. T. Ha, E. Itälä, J. Laksman, S. Urpelainen and E. Kukk, Size selective spectroscopy of Se microclusters, J. Chem. Phys., 2012, 137, 044304.
62. P. Y. Feng and K. Balasubramanian, The low-lying electronic states of indium trimer, Chem. Phys., 1989, 138, 89–98.
63. P. Pyykko, Relativistic effects in structural chemistry, Chem. Rev., 1988, 88, 563–594.
64. T. Nakajima and K. Hirao, The Douglas-Kroll-Hess Approach, Chem. Rev., 2012, 112, 385–402.
65. D. G. Truhlar and X. Li, Introduction to Relativistic Electronic Structure Calculations, J. Phys. Chem. A, 2025, 129, 4301–4312.
66. M. Dolg, H. Stoll, A. Savin and H. Preuss, Energy-adjusted pseudopotentials for the rare earth elements, Theor. Chim. Acta, 1989, 75, 173–194.
67. M. W. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, W. Györffy, D. Kats, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, S. J. Bennie, A. Bernhardsson, A. Berning, D. L. Cooper and M. J. O. Deegan, MOLPRO, version 2010.1, a package of ab initio programs, 2019.
68. Y. Zhao and D. G. Truhlar, The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: Two new functionals and systematic testing of four M06-class functionals and 12 other functionals, Theor. Chem. Acc., 2008, 120, 215–241.
69. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
70. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
71. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter, 2017, 29, 465901.
72. Y. J. Franzke, C. Holzer, J. H. Andersen, T. Begušić, F. Bruder, S. Coriani, F. Della Sala, E. Fabiano, D. A. Fedotov, S. Fürst, S. Gillhuber, R. Grotjahn, M. Kaupp, M. Kehry, M. Krstić, F. Mack, S. Majumdar, B. D. Nguyen, S. M. Parker, F. Pauly, A. Pausch, E. Perlt, G. S. Phun, A. Rajabi, D. Rappoport, B. Samal, T. Schrader, M. Sharma, E. Tapavicza, R. S. Treß, V. Voora, A. Wodyński, J. M. Yu, B. Zerulla, F. Furche, C. Hättig, M. Sierka, D. P. Tew and F. Weigend, TURBOMOLE: Today and Tomorrow, J. Chem. Theory Comput., 2023, 19, 6859–6890.
73. M. Sharma, Y. J. Franzke, C. Holzer, F. Pauly and M. Sierka, Density Functional Theory for Molecular and Periodic Systems in TURBOMOLE: Theory, Implementation, and Applications, J. Phys. Chem. A, 2025, 129(39), 9062–9083.
74. J. Lee, X. Feng, L. A. Cunha, J. F. Gonthier, E. Epifanovsky and M. Head-Gordon, Approaching the basis set limit in Gaussian-orbital-based periodic calculations with transferability: Performance of pure density functionals for simple semiconductors, J. Chem. Phys., 2021, 155, 164102.
75. T. Nian, Z. Wang and B. Dong, Thermoelectric properties of α-In2Se3 monolayer, Appl. Phys. Lett., 2021, 118, 033103.
76. M. K. Horton, P. Huck, R. X. Yang, J. M. Munro, S. Dwaraknath, A. M. Ganose, R. S. Kingsbury, M. Wen, J. X. Shen, T. S. Mathis, A. D. Kaplan, K. Berket, J. Riebesell, J. George, A. S. Rosen, E. W. C. Spotte-Smith, M. J. McDermott, O. A. Cohen, A. Dunn, M. C. Kuner, G.-M. Rignanese, G. Petretto, D. Waroquiers, S. M. Griffin, J. B. Neaton, D. C. Chrzan, M. Asta, G. Hautier, S. Cholia, G. Ceder, S. P. Ong, A. Jain and K. A. Persson, Accelerated data-driven materials science with the Materials Project, Nat. Mater., 2025, 24, 1522–1532.
77. B. Esser, Theoretical analysis of [5.5.6]cyclacenes: Electronic properties, strain energies and substituent effects, Phys. Chem. Chem. Phys., 2015, 17, 7366–7372.
78. A. Somani, D. Gupta and H. F. Bettinger, Computational Studies of Dimerization of [n]-Cyclacenes, J. Phys. Chem. A, 2024, 128, 6847–6852.
79. T. Campos, P. Dally, S. Gbegnon, A. Blaizot, G. Trippé-Allard, M. Provost, M. Bouttemy, A. Duchatelet, D. Garrot, J. Rousset and E. Deleporte, Unraveling the Formation Mechanism of the 2D/3D Perovskite Heterostructure for Perovskite Solar Cells Using Multi-Method Characterization, J. Phys. Chem. C, 2022, 126, 13527–13538.

---