Real and virtual polymorphism of titanium selenide with robust interatomic potentials †

The ﬁ rst successful pairwise potential for a layered material, TiSe 2 , has been parameterised to ﬁ t the experimental data, using a genetic algorithm as the optimisation tool for the parameters of the interatomic potential. This potential has been tested on a wide range of hypothetical isomorphous AX 2 metastable phases using ab initio derived data. From the initial survey, the ground state 1T – TiSe 2 structure remains the lowest enthalpy phase in a wide range of pressures (0 to 25 GPa), which leaves open questions about the nature of a reported unknown high-pressure phase.


Introduction
Titanium diselenide belongs to a class of highly topical transition metal dichalcogenide materials with exciting physical properties and a wide range of applications in advanced electronics and catalysis.Even though research focusing on TiSe 2 goes back many decades, [1][2][3][4][5][6][7] there has been a recent growing scientic interest in this material fuelled by an observed charge density wave (CDW) transition 8 and superconducting properties achieved upon copper intercalation, 9 when applying pressure 10 or under electrical gating. 11Under ambient conditions, the thermodynamically most stable TiSe 2 structure exhibits Se-Ti-Se layers repeated along the c-axis (see Fig. 1), which are bonded together with weak van der Waals (vdW) forces between Se ions of adjacent layers, similar to the CdI 2 structure.This layered compound undergoes a CDW transition at a critical temperature of T c z 200 K, below which the structure adopts a 2 Â 2 Â 2 periodicity with aw ea kla tt i cedi sto r ti on . 12,13he novel physical properties shown by TiSe 2 and other layered quasi-two-dimensional (quasi-2D) materials, e.g.graphene, MoS 2 , and hexagonal boron nitride (h-BN), have stimulated considerable interest in modelling and improving their properties.Modelling at the ab initio level provides a very accurate description of the atomic and electronic structure, but is unfeasible when dealing with large systems and long time-scale simulations.This limitation implies serious compromises when comparing with experiment.However, accurate force elds can provide a very good description of the structure and crystal properties of a material, including, e.g., thermal transport, mechanical deformation, defect formation energies, among others.Recently, for example, long time-scale molecular dynamics (MD) simulations using hundreds of thousands of atoms have been used to study a wide range of phenomena/materials, such as nucleation and microstructure evolution, 14 selfassembly in complex uids, 15 radiation damage, 16 and amorphisation of oxides. 17,18here is, therefore, a strong need for the development of reliable interatomic potentials for layered quasi-2D materials, which would allow us to explore scenarios that are computationally too expensive at the ab initio level, including exploration of the materials' potential energy landscapes, prediction of new structures, structure evolution on time scales of nanoseconds, global optimisation calculations, and atomic structures on the nano/micrometer scale.So far, force elds reported for quasi-2D materials have been based on many-body potentials for covalent materials, see e.g., MoS 2 , [19][20][21] black phosphorus 20 and others. 22The potential models for covalent materials include valence-force eld models and potentials of the Stillinger-Weber, 23 Tersoff 24,25 and Brenner 26 type.The computational cost of these models might, however, be much higher than the one of twobody interatomic potentials, while not providing a more accurate physical description of the interactions in the system.
In this paper, we use TiSe 2 as an example to demonstrate that interatomic interactions in layered materials can, indeed, be described with a polarisable shell model (PSM) two-body interatomic potential (PSM IP).This PSM IP can describe the, at ambient conditions, thermodynamically most stable TiSe 2 crystal structure and its elastic and phonon properties with good accuracy and, to our knowledge, it is the rst successful pairwise potential for a layered material.Moreover, the PSM IP proved to be effective in our exploration of the potential energy landscape of TiSe 2 ,where we nd a good agreement with ab initio results for the energies and possible hierarchy of local minima of TiSe 2 bulk polymorphs and clusters.The quality of the TiSe 2 PSM IP presented here opens up the possibility of studying this material's properties in situations where ab initio calculations are not feasible.
We describe the physical properties and atomic structure of TiSe 2 by a polarisable shell model using pairwise parameter-dependent analytical functions in the form of combined Buckingham, Morse and Lennard-Jones termssee Table 1.Unusually, the model includes short-range interactions between shells and cores of different anions, which allows us to account for strongly non-linear polarisation effects in electron-rich Se species (a local dipole dependence of the short-range interaction).Moreover, the soness of Se results in strong overlap/penetration effects between nearest-neighbour Ti and Se ions, which could not be described by a single simple potential form.A combination of the Morse and Buckingham forms allows for a more thorough search through the functional space when tting parameters of the PSM IP, but also it helps to separate explicitly the short-and long-range attraction and repulsion effects.The rapid increase in the Morse potential for shorter distances effectively alleviates the Buckingham catastrophe for the sensitive range of bonding distances in different coordination environments.We note that our PSM IP has not been tted to any explicit data about the vdW interactions in this system; we have relied on the tting to experimentally observed structure and physical properties.The dielectric properties and phonon dispersion have been key to xt h e spring constants in our tted set.We have also included repulsive Lennard-Jones terms to avoid spurious short interatomic distances caused by the divergence of the Cr À6 Table 1 Potential parameters for the TiSe 2 system.The forms of the interatomic potentials are: E Buckingham ¼ Ae Àr/r À C/r 6 , E Morse ¼ D e ((1 À e (Àa(rÀr0)) ) 2 À 1), E Lennard-Jones ¼ A/r 12 À B/r 6 , where r represents the distance between the ions in question.All other values are parameters of the potentials

Morse potentials
Range ( Å) D e (eV) a ( ÅÀ1 ) r 0 ( Å) Buckingham term.‡ The tting procedure includes the TiSe 2 crystal structure, its elastic constants and phonon frequencies (shown in Table 2) as given in ref.
12 and 13.In total, our PSM IP uses 14 parameters, which were tted to a set of experimental data as summarised in Table 1.
2 Computational methods

Fitting procedure
The tting of the TiSe 2 force eld was performed in three steps using a novel procedure.First, we used a genetic algorithm (GA) as implemented in the GULP code as a global optimisation tool. 27,28500 parallel GA runs were performed in this way, with 1000 generations and 10 000 congurations in the parametric space for each generation.The owchart for each GA run is shown in Fig. 2. The tournament and crossover probabilities were optimised in trial runs and xed at 0.8 and 0.4, respectively.Each GA run had a different random number seed to make sure that the starting point for each run was different.Next, the parameters obtained from the GA were optimised to their closest local minimum in the parametric space using the simplex algorithm 29 and least squares techniques as implemented in GULP.The resulting best set of parameters is presented in Table 1, whereas Table 2 shows the comparison between the observed crystal structure properties and those calculated with the optimised PSM IP.The total energy is calculated as the sum of Coulomb interactions and every force-eld term.We note that the TiSe 2 material has a relatively low polarity; therefore, we have tested models with different degrees of ionicity and found the best t to result from the Ti and Se charges of 1/4 of the formal ionic charge: +1.0e and À0.5e, respectively.We also note that, within the force constant formalism, Takaoka and Motizuki 12 have calculated the TiSe 2 phonon dispersion curves using a PSM with a similar Ti charge of 0.75e.

Monte Carlo search
Monte Carlo routines, as implemented in our in-house knowledge led master code (KLMC), [30][31][32] have been used to identify the Fig. 2 Schematic representation of the GA strategy used for fitting of the PSM IP. ‡ In the case of, e.g.MD and Monte Carlo (MC) simulations, we also anticipate that two ions of the same species can be close to each other, and in that case it is indispensable to counteract the Coulomb attraction of oppositely charged species, i.e. usually the core of one ion and the shell of another.As the short-range potentials are automatically calculated for any pair, we should also exclude the double counting of interactions between the core and the shell of the same ion, which is normally described by a spring potential with a default cutoff of 0.8 Å in GULP.In our PSM IP, this is done by having Lennard-Jones potentials between cores and shells of the same species.For forceelds where both the shell and the core are similarly charged, there is no need to add this potential -such ions will be repelled from each other by the Coulomb interaction.
lowest energy structures at ambient and high pressure.In these simulations, we used an orthorhombic extension of the hexagonal TiSe 2 primitive cell, which is commensurate with other hypothetical TiSe 2 phases (including the C2/m phase that has been proposed based on high pressure Raman data, see below).
Inside of this unit cell, a grid of lattice positions was evenly distributed with a separation of 1.3 Å between points, thus, creating a total of 891 lattice positions of which 18 were randomly occupied.For each run, 10 000 structures were created by KLMC and locally optimised with GULP using our new PSM IP.A good convergence of the energy distribution is already achieved when 2000 structures were used (see Fig. 3).We note that due to the restriction imposed by our settings (1.3 Å between points) and to help optimisation, structures were rst optimised using a rigid ion model and then fully optimised with the shells included.4][35] We report the results of KLMC MC runs at 0 and 25 GPa, with two independent runs for each pressure.Additionally, the effect of the initial volume was evaluated with a larger unit cell, totalling 60 000 KLMCgenerated and GULP-optimised structures.

Ab initio calculations
The reoptimisation of the PSM IP AX 2 structures was carried out by the periodic plane-wave DFT code VASP.Exchange and correlation were included using the generalized gradient approximation (GGA) PBEsol 36 functional with interactions between cores and valence electrons being described by the projector-augmented wave (PAW) 37,38 method.A good total energy convergence to 1 meV was achieved at an energy cutoff of 500 eV and a k-point spacing of 0.3 ÅÀ1 .The PBEsol functional was used for the hybrid calculations.

Results and discussion
In general, there is an excellent agreement for the calculated TiSe 2 lattice parameters when compared with experiment: the lattice parameters a and c are reproduced within ca.0.1% and 0.6%, respectively.There is also a good overall agreement with experiment for the calculated elastic constants and vibrational spectra (see Table 2).We note that the TiSe 2 PSM IP that we have developed are robust and can be used for MD calculations and exploration of the energy landscape, for example, using MC methods.This is achieved by including repulsive Lennard-Jones terms that heavily penalise spurious short distances between species that may result from high velocities or large displacements in various problems e.g.radiation damage or structure prediction, and ionic terms needed to establish the long-range ordering of the ions and their local coordination spheres, yielding a globally robust description of the energy landscape as we show below.
To check the robustness and transferability of our PSM IP we have investigated the TiSe 2 global energy landscape by modelling the geometry and energetics of the thermodynamically most stable TiSe 2 structure and other hypothetical sixteen AX 2 polymorphs, which were taken from well-known crystal structures with an AX 2 stoichiometry.Every structure was then fully reoptimised (including both lattice parameters and internal atomic coordinates) with the periodic plane-wave DFT code VASP (Vienna Ab initio Simulation Package) [39][40][41][42] keeping  symmetry constraints.The energy ranking (relative energies per TiSe 2 formula unit) and comparison between the PSM IP and DFT potential energy landscapes are shown in Fig. 4.
As we have noted above, the PSM IP calculations and experiment show good agreement for the ground state structure (Table 2).We also observe a very similar ranking of the local minima for PSM IP and DFT energies (Fig. 4).Moreover, during global searches in the low-energy region of the potential energy landscape with the threshold algorithm 43,44 as implemented in the G42+ code package 45 no structure with lower energy was observed, indicating that the PSM IP yields an accurate representation of the potential energy landscape; an exhaustive exploration of the TiSe 2 potential energy landscape is currently in progress.Fig. 4 shows that there are many similarities between the PSM IP and DFT energy landscapes around known thermodynamically metastable phases of isomorphous AX 2 crystalline compounds as described above.We note that the PSM IP energy differences are, however, lower than those at the DFT level.Similar high energy differences from a DFT approach have been observed in other systems, e.g., Cu/ZnO, 33 Ti clusters, 34 or Cu 4 Ag 4 clusters. 46Structurally, the optimised PSM IP structures are close to the corresponding DFT minima.The reliability of our new set of PSM IP is such that any observed structural differences between PSM IP and DFT calculations disappear on optimisation aer just a few small atomic displacements, pointing to the structural (geometric) similarity between the two landscapes.We observe in the absolute majority of cases that the optimised IP structures are very close to their DFT local minima.We have calculated the energy difference per TiSe 2 unit between the initial (DFT SP ) and nal structures obtained from the DFT optimisation process: for more than 70% of the structures this difference is less than 0.05 eV and in some cases (more than 17%) is less than 0.01 eV.This structural similarity is also reected in the small number of ionic steps needed to achieve convergence, which is fewer than 15 in about half of the cases.We note that, conversely, when the DFT AX 2 structures are re-optimised with the PSM IP, these structures (except ZnCl 2 ) revert to the ones shown in the PSM IP column in Fig. 4, indicating a close similarity of the shape of the PSM IP and the DFT landscape in the neighbourhood of the local minima.Only ZnCl 2 optimises to a higher energy conguration instead.On a closer inspection of the four thermodynamically accessible, lowest energy structures, we note that the mismatch between PSM IP and DFT for the lattice parameters is below 5%, with an average of 2.5%.As shown by the data markers (Fig. 4), the lowest energy structures have an octahedral environment.
Remarkably, the lowest energy conguration and the CdCl 2like structures are very close structurally and energetically on the PSM IP landscape.However, DFT enhances this small structural difference with an energetic difference per TiSe 2 unit of 0.014 eV between the two structures compared to the 2 Â 10 À4 eV calculated with the PSM IP.The three lowest energy structures are layered, whereas the fourth one is a distortion of the HgCl 2 -like phase, which retains its structure on DFT optimisation.Higher energy structures are expected to be difficult to access with conventional synthetic methods.Moreover, the new PSM IP predict that the thermodynamically accessible polymorphs are stable for a wide range of pressures with no phase transition (see Fig. 5 where the evolution of lattice enthalpy with pressure is shown for all lowest energy compounds investigated in the range of 0 to 25 GPa).§ Therefore, we conclude that there is a clear agreement between the new PSM IP and DFT, for different AX 2 crystal structures, at least as far as the stability and the local neighbourhood near the various minima, and the overall robustness far away from the minima, are concerned.The robustness of this PSM IP is in strong contrast with essentially all other (simple) interatomic potentials, which have been employed for the description of layered quasi-2D compounds in the past and have encountered serious problems when applied to regions of the energy landscape away from the atomic congurations to which they had been tted (e.g.Ref. 21).
The success of this approach in developing and parameterising interatomic potentials rests, we suggest, on the physically sound choice of a robust and transferable model.In particular, a simple chemical analysis of the bonding in TiSe 2 tells us that atoms are linked by simple heteropolar sigma-type bonds.In such systems, one can expect the dominant role of pairwise interactions and the need for an adequate description of Se polarisability.Both of these features are oen absent or underestimated in commonly applied many-body potential models.
As a rst application, we have also investigated the structure and properties of very small TiSe 2 clusters.This exercise is, however, purely illustrative and a brief discussion is given in the ESI, † whereas a more detailed study on TiSe 2 clusters will be published soon.In general, there is a good agreement between our PSM IP model and DFT, demonstrating, again, the predictive power of the new PSM IP.
We note a report of a possible pressure driven structural phase transition in the range of z15 À 19 GPa (ref.47) as shown by the evolution of the characteristic Raman signatures of fundamental optical phonon modes in TiSe 2 .It is proposed that as with 1T-TiTe 2 and IrTe 2 , 48,49 the TiSe 2 undergoes a transition to a non-layered monoclinic C2/m phase structure.We have checked this hypothesis using both our newly developed PSM IP and ab initio calculations in the range of relevant pressures.Both methods showed an excellent agreement with each other, but curiously no sign of a phase transition was found.The energy offset (in eV per TiSe 2 formula unit) between the two structures has been calculated to be 0.15 (0.17), 0.25 (0.39), 0.32 (0.48) for PSM IP, standard DFT and hybrid DFT, respectively at 0 GPa (25 GPa).Furthermore, we have performed an unbiased MC search for lowest enthalpy structures -at 0 and 25 GPausing orthorhombic supercells of the TiSe 2 structure (which is hexagonal in the primitive setting), commensurate with C2/m phase.Results of our search, which sampled over 60, 000 local minima, are illustrated in Fig. 6, where we observed that the layered TiSe 2 structure remains the ground state across the full range of pressures investigated; the main data mined AX 2 structures are rediscovered by the unbiased search; and a new type of metastable phase is predicted, which can be considered either as a defective layered structure, in which the layering pattern is regularly interrupted with Ti polyhedra displaced into the interlayer space connecting two adjacent layers,{ or alternatively as "generalization" of the C2/m structure with the polyhedral units extended in one or two dimensions.The hypothetical pressure induced phase transformation may be, therefore, a result of stoichiometry changes and will be affected by impurities present in the experimental material.Our PSM IP model, however, would need to be extended to describe such effects.In this context, we also note the multitude of possible polytypes, which can be generated by different close packings of the Se-atom layers and different distributions of the Ti atoms over the octahedral holes.It is planned to employ the so-called PCAE (Primitive Cell approach for Atom Exchange) method 51 to generate such structure candidates and explore their structural and electronic properties.A systematic investigation of the energy landscape of TiSe 2 will be reported in the near future, where we hope to be able to answer this experimental challenge.
To further check the PSM IP reliability and performance outside of the tting region, we have taken the lowest 1000 TiSe 2 structures from the distribution shown in Fig. 6.This set includes structures with energies as far as 2.5 kT from the global minimum.The energy for the structures was calculated using DFT SP (single point DFT), and the correlation between IP and DFT SP energies is evident from Fig. 7.In general, we see an excellent agreement for structures up to 1.75 kT higher in energy than the global minimum with R 2 values close to 0.9.We note that we compare the energies of alike structures instead of the ones fully optimised with each approach.Comparing the latter ones throws up more issues/aspects to consider, which are discussed in detail in, e.g., ref. 33.The greater energy dispersion at higher energies is due to the fact that the TiSe 2 structures  { Such structure elements have also been found in simulations of MgF 2 lm growth. 50ere randomly generated and have, most likely, no symmetry.They land, however, on parts of the energy landscape far away from the main basin and would be very difficult to access with conventional synthetic methods.Furthermore, the 1000 TiSe 2 structures from the previous test have been fully DFT reoptimised and their energies calculated with our PSM IP (IP SP , optimisation of shells; ions and lattice constants are kept xed).Thus, similar to above, we compare the energies (see Fig. 8) of alike structures, which in this case are the DFT optimised structures.In Fig. 8, we identied the structures with Se-Se distances (d Se-Se ) shorter than 3.0 Å.This bonding cutoff criterion (which is also widely used in molecular mechanics modelling) was chosen as 1.2 times the Se covalent radii (1.22 Å), such as d Se-Se ¼ 1.20 Â 1.22 Å + 1.20 Â 1.22 Å z 3.0 Å.We have excluded those structures because there is a strong indication of Se-Se bonds and, therefore, a change of oxidation states.We emphasise that the current PSM IP model does not perform well for oxidation state changes.Without the "outliers", there is an excellent agreement for structures up to 2 kT higher in energy than the global minimum with R 2 values above 0.9.

Conclusions
In summary, the new optimised TiSe 2 PSM IP, presented in Table 1, performs well in reproducing the physical properties and crystal structure of the thermodynamically most stable polymorph.This new PSM IP is also applicable for other hypothetical polymorphs of TiSe 2 , showing that it should be capable of identifying low energy structure candidates via a global search or during long MD simulations, for a subsequent renement on an ab initio level.Thus, this PSM IP constitutes both an accurate and globally applicable representation of the energy landscape of TiSe 2 , and could serve as a starting point for similar potentials to other layered quasi-2D compounds (e.g.Ti, Zr and Hf chalcogenides).

Fig. 1
Fig. 1 (a) Side and (b) top views of the thermodynamically most stable TiSe 2 crystal structure under ambient conditions.The unit cell is represented by the solid black lines.Gray and green colours are reserved for Ti and Se ions, respectively.

Fig. 3
Fig. 3 MC convergence.PSM IP total energy distribution as a function of the number of MC-generated structures.

Fig. 4
Fig. 4 Energy map ranking of AX 2 structures across the PSM IP and DFT energy landscapes.To facilitate the comparison, calculated DFT single point energies (DFT SP ) for the IP structures are included.DFT SP and DFT (optimised) energies are referenced to the ground state.The DFT energy as calculated for every fully reoptimised structure is linked by lines to the IP energy.Data markers are used to distinguish between structures with 6-coordinated Ti ions (circles) and those with ad i fferent coordination (squares).

Fig. 6
Fig. 6 Density distribution (with a Gaussian smearing of s ¼ 0.01 eV) of the PSM IP energies from the KLMC-generated TiSe 2 structures.In this plot, purple and red are reserved for 0 and 25 GPa results.The calculated energies from the data mining procedure are shown as comparison below the curves.

Fig. 7
Fig. 7 PSM IP and DFT SP energy comparison for the lowest 1000 PSM IP energy structures of our KLMC-MC approach.

Table 2
12,13rison of the TiSe 2 bulk properties -lattice parameters, elastic constants and phonon frequencies-as obtained by the new potentials and compared with experiment12,13