Computational Investigation of Copper Phosphides as Conversion Anodes for Lithium-Ion Batteries

Using first-principles structure searching with density-functional theory (DFT), we identify a novel Fm3̅m phase of Cu2P and two low-lying metastable structures, an I4̅3d–Cu3P phase and a Cm–Cu3P11 phase. The computed pair distribution function of the novel Cm–Cu3P11 phase shows its structural similarity to the experimentally identified Cm–Cu2P7 phase. The relative stability of all Cu–P phases at finite temperatures is determined by calculating the Gibbs free energy using vibrational effects from phonon modes at 0 K. From this, a finite-temperature convex hull is created, on which Fm3̅m–Cu2P is dynamically stable and the Cu3–xP (x < 1) defect phase Cmc21–Cu8P3 remains metastable (within 20 meV/atom of the convex hull) across a temperature range from 0 to 600 K. Both CuP2 and Cu3P exhibit theoretical gravimetric capacities higher than contemporary graphite anodes for Li-ion batteries; the predicted Cu2P phase has a theoretical gravimetric capacity of 508 mAh/g as a Li-ion battery electrode, greater than both Cu3P (363 mAh/g) and graphite (372 mAh/g). Cu2P is also predicted to be both nonmagnetic and metallic, which should promote efficient electron transfer in the anode. Cu2P’s favorable properties as a metallic, high-capacity material suggest its use as a future conversion anode for Li-ion batteries; with a volume expansion of 99% during complete cycling, Cu2P anodes could be more durable than other conversion anodes in the Cu–P system, with volume expansions greater than 150%. The structures and figures presented in this paper, and the code used to generate them, can be interactively explored online using Binder.


■ INTRODUCTION
Graphite is the most commonly employed lithium-ion battery (LIB) anode but is inherently limited by a maximum theoretical capacity of 372 mAh/g upon the formation of LiC 6 . Phosphorus (black or red) has a significantly higher theoretical capacity of 2596 mAh/g due to the formation of Li 3 P; however, it suffers from capacity deterioration, primarily caused by deleterious volume expansion that occurs upon charging, which constrains the capacity to 350−500 mAh/g in a limited voltage window. 1 In addition, P and its lithiated phases have limited electrical conductivity, requiring dopants and additives to improve performance. By adding transition metals to P, through nanostructuring or synthesis, both electrical conductivity and stability during cycling can be enhanced. 2 Transition-metal phosphides (TMPs) provide a large design space in which to engineer such high-capacity, conversion anodes for LIBs. 3 High-throughput computational screening has previously identified TMPs with high capacities for LIB electrodes, including TiP, Co 2 P, Mn 2 P, and others. 4 As conversion anodes for LIBs, TMPs offer both the added gravimetric capacity (ranging from 500 to 2000 mAh/g 4 ) and stability against volume expansion over several battery cycles. 5 In addition to bulk or powdered TMPs being used as LIB conversion anodes, 6 nanostructured TMPs can often display improved electrochemical cycling performance. 7 Despite these efforts, TMPs have yet to be widely adopted as conversion anodes, given the large volume expansion (between 150 and 300% 4 ) exhibited by anodes with high P content, which limits their cyclability. Despite this drawback in volume expansion, TMPs show higher average voltages than graphite, which has an average voltage of 0.1 V. For example, CoP has an average voltage of 0.67 V, the ternary metal phosphide LiFeP has an average voltage of 0.4 V, and MnP has an average voltage of 0.62 V. 4 Higher average voltages give the metal phosphides improved safety while sacrificing energy density, making them an ideal choice for large-scale and long-term energy storage.
Several previously studied TMP anodes include FeP x=1,2,4 , 8 Fe 2 P nanoparticles, 9 Ni 2 P, 10 CuP 2 , 11,12 and Cu 3 P, 13,14 among others. Of the TMPs tested as conversion anodes, the copper phosphides (specifically CuP 2 and Cu 3 P) have shown promise for their cyclability and capacity. The copper phosphides offer additional benefits to the other TMPs, as Cu is already used as a common current collector, providing further cycling capability and resistance to degradation. 15 Cu 3 P prepared by hightemperature synthesis had a first-cycle capacity of 527 mAh/g, 13 and a porous Cu 3 P anode synthesized by facile chemical corrosion exhibited a capacity between 360 and 380 mAh/g over 70 cycles. 16 The capacity of high-temperaturesynthesized Cu 3 P exceeds that of graphite, and the cyclability of porous Cu 3 P is improved relative to other Cu 3 P anodes. 13,17 CuP 2 , on the other hand, delivers a higher initial capacity of 815 mAh/g, but can only be cycled stably 10 times before the capacity fades to 360 mAh/g. 12 The main factor in this degradation is the high concentration of P in the CuP 2 , which, while enabling high capacity, also contributes to the structural instability of CuP 2 during cycling as the lithium-rich Li 3 P phase forms. To optimize the trade-off between stability and capacity, it would be beneficial to discover a compound with higher P content than Cu 3 P to offer higher capacity, and with a Cu content higher than CuP 2 to aid in cyclability.
By performing crystal structure prediction, combining both ab initio random structure searching (AIRSS) and a genetic algorithm (GA), in addition to structural prototyping with known crystal structures of related chemistries, 18−20 we produced the compositional phase diagram of the copper phosphide system. We describe this approach to structure prediction and the application of open-source Python packages matador (v0.9), 21 for high-throughput first-principles calculations, and ilustrado (v0.3), 22 for computational structure prediction with GAs. Crystal structure prediction for battery anodes is a well-tested method, 23 used for identifying both novel anode materials 4 and unknown phases, which form during battery cycling. 24,25 AIRSS has been used previously to search for additional phases of Li−P and Na−P, which form during battery cycling. 26 The GA was also employed to search for new phases of Na−P, which were confirmed experimentally by solidstate NMR spectroscopy. 27 As applied here to Cu−P, these methods predict a novel metallic Fm3̅ m−Cu 2 P phase at 0 K, within the target composition range of Cu 1<x<3 P, for a highcapacity, low-volume expansion conversion anode; we compare its electronic structure to other TMPs to show a similarity to Fm3̅ m−Rh 2 P and Fm3̅ m−Ir 2 P. Two other phases, Cm−Cu 3 P 11 and I4̅ 3d−Cu 3 P, are identified as metastable, both bearing structural similarity to known copper phosphides. We calculate the convex hull of Cu−P at temperatures up to 600 K, confirming the dynamic and chemical stabilities of Cu 2 P across this temperature range. A ground-state voltage profile from density-functional theory (DFT) shows that Fm3̅ m−Cu 2 P undergoes the same lithiation process as P6 3 cm−Cu 3 P; however, Fm3̅ m−Cu 2 P has a higher capacity of 508 mAh/g, with an average voltage of 0.86 V versus Li/Li + (compared to 0.91 V for P6 3 cm−Cu 3 P).

■ METHODS
To search for novel copper phosphides, we first performed structural relaxations of the 13 structures from the Inorganic Crystal Structure Database (ICSD) 28 of Cu x P (0 < x < 1). The Python package matador 21 was used to query 1053 prototype binary structures from the Open Quantum Materials Database (OQMD), 19 with chemical compositions containing a pnictogen and a transition metal from the first two rows, namely, {Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Zr, Nb, Mo, Tc, Ru, Rh, Pd, Ag, Cd}−{P, As, Sb}; each composition was then transmuted to the corresponding stoichiometry of Cu−P, yielding 909 unique structures after geometry optimization. To extend this search beyond existing prototypes, two additional structure prediction steps were performed, namely, AIRSS 29 and an evolutionary search with the GA implemented in the ilustrado 22 package.
When performing AIRSS, one proceeds by generating random "sensible" (symmetry, density, and atomic separation constrained) trial cells and then geometry-optimizing them to their corresponding local minima. All relaxations can be performed concurrently, with no interdependence between calculations. New trial structures are generated until the ground state of each stoichiometry (within the constraints of the search) has been found multiple times.
We initially performed an exploratory AIRSS search consisting of around 5000 trial structures, with constraints on cell size, stoichiometry, and number of atoms in the cell. In this initial search, the total number of atoms in the cell was constrained to be ≤40, and the total number of formula units was randomized between 1 and 4, while still keeping the total number of atoms below 40. The number of atoms of Cu and P were randomized between 1 and 9 in each cell, and the cell volume (V) was constrained based on the total number of atoms in the cell (N) to be 8N Å 3 ≤ V ≤ 20N Å 3 , based on the average densities of Cu−P phases within the ICSD.
Structures from the searching and enumeration procedures were then used, with fitness weighted according to their distance from the convex hull, as the initial configurations for a GA implemented in the Python package ilustrado. 22 The ilustrado package uses a simple cut-and-splice crossover operation, supplemented by mutation operators (random noise, atomic permutations, vacancies, and adatoms). 30 To avoid stagnation, each trial structure was filtered for similarity (via pair distribution function, PDF, overlap) against existing structures in the population. Three independent GA runs were performed with 10 generations each, yielding a further 1049 relaxed structures. Finally, a directed AIRSS search of Cu x P y , where x + y < 8, was performed to create a final set of ∼20 000 structures within the Cu−P chemical space. In all cases, to constrain the search to physically reasonable structures, a minimum atomic separation of 1.5 Å was enforced and the maximum number of atoms in the cell was constrained to 10 for the initial ∼10 000 AIRSS searches and 40 atoms per cell for the final ∼3000 trials.
All calculations were performed using CASTEP (v18.1 and v19.1), the plane wave pseudopotential DFT package. 31 To maximize computational efficiency, the initial calculations were performed with loose convergence criteria that ensured formation energies converged to 10 meV/atom. The Perdew−Burke−Ernzerhof (PBE) exchange− correlation functional was used 32 with Vanderbilt ultrasoft pseudopotentials 33 that required a plane wave kinetic energy cutoff of 300 eV to converge energies to within 10 meV/atom. The Brillouin zone (BZ) was sampled with a Monkhorst−Pack grid k-point spacing finer than 2π × 0.05 Å −1 ; the grid was frequently recomputed to accommodate any changes in the cell shape and size during relaxation. Each structure was geometry-optimized at this accuracy to a force tolerance of 0.05 eV/Å. The structures with a formation energy within 50 meV of the convex hull were then further optimized once more using CASTEP's on-the-fly (OTF) "C18" library of ultrasoft pseudopotentials, a with a finer k-point sampling of 2π × 0.03 Å −1 and plane wave kinetic energy cutoff of 500 eV, which yielded formation energies converged to within 2.5 meV/atom. To predict the voltage profiles with the same convergence criteria (formation energies within 2.5 meV/atom), the relaxation of known Li−P structures required a higher plane wave cutoff of 700 eV. Therefore, to compare ternary phases of Cu−Li−P in the voltage profile, all Cu−Li−P phases were reoptimized at a plane wave kinetic energy cutoff of 700 eV.
To identify stable structures from this search, a convex hull of the copper phosphides was created. The formation energy E f of each structure Cu x P y was calculated using where E(Cu) is the DFT total energy of the Fm3̅ m−Cu structure, and E(P) is the energy of Cmca-P (black phosphorus). Black phosphorus was used as the P chemical potential instead of the lower-energy polymorph red phosphorus; as has been previously discussed in Mayo et al., 26 black phosphorus is commonly used when making electrochemical cells. 34 Electrochemical voltage profiles for Li insertion into the stable Cu−P phases were calculated from the computed formation energies from the ternary convex hull of Cu−Li−P. To calculate the voltage profiles shown in the Cu 2 P as a Li-Ion Battery Conversion Anode section, the voltage, V̅ , between two tie-lines in the ternary convex hull with compositions Li x 1 Cu n P and Li x 2 Cu n P was calculated using as stated by Urban et al. 35 In eq 2, E(Li x 1 Cu n P) and E(Li x 2 Cu n P) are the ground-state energies of two phases along the reaction pathway of the ternary convex hull, in which x 1 and x 2 are the relative amounts of Li in the starting and ending products at each point in the pathway. All phonon calculations were performed under the harmonic approximation with the PBE xc-functional in a 2 × 2 × 2 supercell (corresponding to a phonon q-point spacing of 2π × 0.046 Å −1 for Fm3̅ m−Cu 2 P) using the finite displacement method implemented in the CASTEP code. The dynamical matrix was then Fourierinterpolated onto the BZ path provided by the SeeK-path Python package 36,37 to compute the phonon dispersion and onto a fine Monkhorst−Pack grid to compute the phonon density of states.
The band structure for Cu 2 P was calculated using the higher accuracy parameters and pseudopotentials mentioned previously, and the electronic density of states was integrated and projected onto atomic orbitals using the OptaDOS code. 38,39 Vibrational properties of all stable phases were computed using the finite displacement method, with an added many-body dispersion (MBD denoted MBD* in CASTEP v19.0) correction 40

to account for interlayer interactions in black phosphorus.
The open-source Python package matador (v0.9) 21 was used to run the CASTEP calculations, perform the analysis, and create the plots found in this article. All of this analysis, as well as the underlying source code and data, can be explored interactively using Binder and is found on GitHub at harpaf13/data.copper-phosphides. The input and output files associated with our calculations have been deposited into the Cambridge Repository at https://doi.org/10.17863/CAM.54795.

■ RESULTS
From the search of ∼20 000 trial structures, there are 42 unique phases within 50 meV/atom of the convex hull. Previous computational structure searches have used a distance above the hull of 25 meV/atom, 41 and given the accuracy of PBE, 42 we chose to increase this cutoff to 50 meV/atom. Furthermore, the experimentally verified P6 3 cm−Cu 3 P structure 43 is 37 meV/atom above the convex hull tie-line, further justifying this cutoff. The uniqueness was determined by computing pairwise overlap integrals of the pair distribution functions of phases at each stoichiometry using matador. The set of 42 unique phases contains four experimentally reported copper phosphides from the ICSD: P1̅ −CuP 10 synthesized by a mineralization reaction 44 and C2/m−Cu 2 P 7 , 45 P2 1 /c− CuP 2 4545 and P6 3 cm−Cu 3 P 43,46 from high-temperature sintering.
Oloffson's experiments on single-crystal P6 3 cm−Cu 3 P synthesized at high temperature, and subsequent work by De Trizio et al., 47 show that Cu 3 P has several defects 43 with a range of stoichiometries between Cu 2.6 P and Cu 2.8 P. DFT studies of the Cu vacancies indicate that Cu 3 P is substoichiometric, 47 and to search this substoichiometric space, unit cells of P6 3 cm− Cu 18 P 6 were enumerated with 1, 2, and 3 Cu vacancies, resulting in 76 Cu 3−x P structures. The lowest-energy defect was a Cmc2 1 − Cu 8 P 3 (Cu 2.67 P) phase 26 meV/atom above the convex hull tieline, denoted as vacancy enumeration in Figure 1.
The convex hull of Cu−P, with points colored by the provenance of each structure, is presented in Figure 1; the experimentally identified phases, and a new Fm3̅ m−Cu 2 P phase, all lie on the convex hull tie-line and are each labeled with an arrow.
Details of the 24 structures, which are both negative in the formation energy relative to Cu and P and are within 50 meV/atom of the convex hull, are given in Table 1. Phases on the convex hull tie-line in Figure 1 are indicated with an asterisk in Table 1, and phases that are experimentally confirmed are denoted by their ICSD Collection Code as "ICSD #". Phases not reported previously, within 20 meV/atom of the convex hull tie-line, are bolded in Table 1. The provenance of each phase is given in the last column of Table 1. Phases that were found by swapping the elements of a prototype ICSD structure are denoted by the ICSD Collection Code of the prototype used as "prototype #".
Of the 24 binary structures in Table 1, 9 were discovered by AIRSS, 6 by the GA, 4 from structural prototyping, and 4 were previously known Cu−P structures from the ICSD. Of particular interest are three new phases, bolded in Table 1, Fm3̅ m−Cu 2 P, I4̅ 3d−Cu 3 P, and Cm−Cu 3 P 11 , which are all within 20 meV/atom of the convex hull and will be discussed further in the following sections.
Phosphorus-Rich Phase Cm−Cu 3 P 11 . Cm−Cu 3 P 11 is a new structure that was found by relaxing the prototype Ag 3 P 11 (ICSD 26563); it is 17 meV/atom from the hull tie-line and has structural similarity to the ICSD structure C2/m−Cu 2 P 7 (ICSD 35281), 45 as shown in Figure 2. Both of these structures have repeating chains of P atoms, as seen in the supercells in Figure 2, in which alternating patterns of Cu or Cu−P are connected to a zig-zag chain of P atoms. All known phases in the P-rich (Cu x P, where x > 1) region of the convex hull, namely, C2/m−Cu 2 P 7 (ICSD 35281), 45 P2 1 /c−CuP 2 (ICSD 35281), 45 and P1̅ −CuP 10 (ICSD 418805), 44 have long chains of P atoms, similar to the layered P2/c-P (ICSD 29273, 54 red P).
In the P-rich region, five new phases were identified within 50 meV/atom of the convex hull: Pmmn−CuP 3 , Cm−Cu 8 P 27 , Cm−Cu 3 P 11 , Cm−Cu 7 P 27 , and Cm−CuP 4 . Using the GA, it was possible to include structures with stoichiometries of P up to Figure 1. Convex hull of Cu−P phases from structure searching. Four structures lie on the convex hull: CuP 10 , Cu 2 P 7 , CuP 2 , and Cu 2 P. Structures are colored according to their provenance: either from a searching method (AIRSS, GA, prototyping, vacancy enumeration) or from an existing database (ICSD). Prototyping refers to using a prototype structure from the ICSD and replacing the atoms with Cu or P, as described in the Methods section. The vacancy enumeration phases are the phases optimized after adding Cu vacancies to P6 3 cm− Cu 3 P. Phases within 20 meV/atom of the convex hull lie within the shaded gray region, and Cu 3 P is labeled for reference.
Chemistry of Materials pubs.acs.org/cm Article 27 atoms in the unit cell and thus found structural variations on Cu 2 P 7 such as Cu 3 P 11 . To compare the new metastable Cm− Cu 3 P 11 structure with other P-rich structures, the pair distribution functions (PDFs) and calculated powder X-ray diffraction (PXRD) peaks of CuP 2 , Cu 2 P 7 , and Cu 3 P 11 were calculated and are compared in Figure 3. In all three cases, the initial sharp peak in the PDF between 2.20 and 2.24 Å shows, unsurprisingly, the same Cu−P and P−P distances shared by all three structures. The peaks at radii above 3 Å show the longerrange similarity between Cu 3 P 11 and Cu 2 P 7 , which is not shared by CuP 2 . Comparison of the PXRD patterns of C2/m−Cu 2 P 7 and Cm−Cu 3 P 11 shows that Cm−Cu 3 P 11 is distinguished by a peak at a 2θ value of 16°, where C2/m−Cu 2 P 7 has an indistinguishable peak at this point. Given the shared symmetry operations between Cm and C2/m, we expect to see peaks at the same 2θ values, but the intensities will vary between the structures. We deduce that these three phases could be verified using experimental PXRD, by using the peaks at 2θ < 30°to distinguish between the phases. Cu 3−x P Phases (x ≤ 1). Within the stoichiometry range Cu 3−x P (x ≤ 1), four unique Cu 3 P phases, Cu 17 P 6 , Cu 8 P 3 , Cu 7 P 3 , and Cu 2 P, were found. Of these, P6 3 cm−Cu 3 P was the only phase previously experimentally determined and had a formation energy of 37 meV/atom above the convex hull tieline. Olofsson identified the stoichiometry of P6 3 cm−Cu 3 P at 975 K to be between Cu 2.867 P and Cu 2.755 P due to Cu vacancies within the unit cell of P6 3 cm−Cu 18 P 6 (shown in Figure S1). 43 A study on low-temperature phases of Cu 3−x P proposes phases from Cu 2.3 P to Cu 2.9 P. 55 The lowest-energy Cu 3−x P (x ≤ 1) phases identified in Table 1, P1−Cu 17 P 6 (Cu 2.83 P), Cmc2 1 − Cu 8 P 3 (Cu 2.66 P), and P1−Cu 7 P 3 (Cu 2.33 P), are all defect structures of P6 3 cm−Cu 3 P with 1, 2, and 4 Cu vacancies respectively, from the P6 3 cm−Cu 18 P 6 unit cell of Cu 3 P. Of these three P6 3 cm−Cu 3 P defect structures, Cmc2 1 −Cu 8 P 3 (Cu 2.66 P) has the smallest distance from the hull (ΔE = 26 meV/atom). This corroborates the previous DFT calculations suggesting Cu 3 P has two Cu vacancies. 47 In addition to the ICSD phase of P6 3 cm−Cu 3 P (ΔE = 37 meV/atom), two other Cu 3 P phases were found, which are closer to the convex hull tie-line than P6 3 cm−Cu 3 P; these are the P2 1 /m−Cu 3 P (ΔE = 30 meV/atom) and I4̅ 3d−Cu 3 P phases (ΔE = 11 meV/atom). The P2 1 /m−Cu 3 P phase is structurally related to the Fm3̅ m−Cu 2 P (ΔE = 0 meV/atom) phase (discussed in the following section). These two phases are shown in Figure 4, in which the P2 1 /m−Cu 3 P can be described as a stacking of the Fm3̅ m−Cu 2 P phase. While the Fm3̅ m−Cu 2 P phase has not been observed experimentally, it is likely that the two phases could be distinguished, given their distinct PDF and PXRD patterns shown in Figure S2. The PXRD pattern for P2 1 /m−Cu 3 P has additional low-intensity peaks to the right of the 46°peaks and is distinct from the other low-energy phases of Cu 3 P as shown in Figure S2, which would further distinguish this phase in experiment.
The lowest-energy Cu 3 P phase is an I4̅ 3d phase 11 meV/atom above the tie-line, which was identified by relaxing the prototype I4̅ 3d−Cu 3 As structure (ICSD 64715 49,57 ). The I4̅ 3d−Cu 3 P structure is the highest-symmetry Cu 3 P phase and is the only cubic phase in the set of low-energy Cu 3 P structures. I4̅ 3d− Cu 3 P contains 8 formula units in the primitive unit cell and has 9-fold coordinated P atoms, whereas P6 3 cm−Cu 3 P has 8-fold coordinated P atoms. The resulting crystal structures, shown in Figure S1, show two different long-range orderings of the Cu subnetwork. P6 3 cm−Cu 3 P has only one, 8-fold coordinated, P site, which results in continuous zig-zag chains of Cu atoms surrounding the P, which are at the peaks of the buckles in the zig-zag. In I4̅ −Cu 3 P, there are two 9-fold coordinated sites: one site at the center of the surrounding Cu (seen in Figure S1) and one at the edges, which together form a hexagonal Cu cage surrounding the P atom in the center. While both phases have high-coordinated P atoms, the I4̅ 3d−Cu 3 P shows a network of Cu atoms surrounding a central P atom, where P6 3 cm−Cu 3 P contains infinite Cu chains in the c direction.
Another trigonal phase, P3̅ c1−Cu 3 P (ICSD 16841, 49 ΔE > 50 meV/atom), has the same structure as P3̅ c1−Cu 3 As (ICSD 16840 49 ); however, it is 82 meV/atom above the convex hull tieline. To the best of our knowledge, there are no reports of an I4̅ 3d−Cu 3 P phase, either experimentally or in a computational database. The PDF and PXRD patterns of I4̅ 3d−Cu 3 P given in Figure S2 show no relation to any other Cu 3 P phase, or the Fm3̅ m−Cu 2 P phase; thus, if energetically stable, it could be distinguished using PXRD in the experiment.
Fm3̅ m−Cu 2 P. The Fm3̅ m−Cu 2 P phase was found from the prototype Fm3̅ m−Ir 2 P (ICSD 640898). 58 Comparing the Cu 2 P phase to both Ir 2 P and Rh 2 P using PDFs in Figure S3 shows that the PDFs are identical between all three structures, and the  49 c Structure from single-crystal diffractometry. 43 d Prototype structure I4̅ −Cr 3 P by single-crystal X-ray diffraction. 50 e Prototype structure Fm3̅ m−Rh 2 P by X-ray diffraction. 51 f Structure from X-ray diffraction. 45 g Structure from X-ray diffraction. 45 h Prototype structure Cm−Ag 3 P 11 by single-crystal X-ray diffraction. 52 i Structure from single-crystal X-ray diffraction. 44 j Black phosphorus structure from powder X-ray diffraction. 53 Chemistry of Materials pubs.acs.org/cm Article PXRD plot of Cu 2 P has the same peaks, all shifted to slightly higher values of 2θ due to structural relaxations in the geometry optimization of Cu 2 P. Previously, a 2D structure of Cu 2 P was predicted theoretically as a buckled nonmagnetic material, 56 in which the magnetism expected was inhibited by the buckled layers. The buckled layers from the 2D phase are also present in the bulk Fm3̅ m−Cu 2 P, and the nonmagnetic nature was confirmed in the bulk phase by the lack of spin-polarization in the density of states shown in Figure  S4. The bulk Fm3̅ m−Cu 2 P structure described above has the same structural motifs as the 2D hexagonal phase predicted by Yang et al. 56 and has the same electronic properties.
Fm3̅ m−Cu 2 P lies on the convex hull tie-line and is energetically more stable than both the experimentally confirmed phase of P6 3 cm−Cu 3 P and its defect structure Cmc2 1 −Cu 8 P 3 . Figure 5 shows the phonon dispersion for the Fm3̅ m−Cu 2 P computed as mentioned in the Methods section. No imaginary phonon frequencies were present in the dynamical matrix (interpolated or otherwise), indicating that Fm3̅ m−Cu 2 P is dynamically stable.
The electronic structure of Fm3̅ m−Cu 2 P is related to the electronic structure of other Fm3̅ m TMPs, suggesting that it belongs to the same class of materials as Fm3̅ m−Ir 2 P and Rh 2 P. Of the TMPs in the Materials Project database, 59 21 are insulating and 68 are metallic with a high density of transitionmetal d-bands below the Fermi level. Figure 6 shows the electronic band structure and density of states of Fm3̅ m−Cu 2 P projected by species along the high-symmetry path from SeeKpath used previously and the density of states projected by the angular momentum channel on a fine Monkhorst−Pack grid. The band structure shows that Cu 2 P is a metal with P and Cu bands touching at the Γ point of ∼2.0 eV above the Fermi level. In addition, there is a characteristic high density of flat bands localized on the Cu ions that exhibit a d-character of around 2.5 eV below the Fermi level. Calculating this band structure using the HSE06 functional (shown in Figure S5), a hybrid functional designed to correct for band-gap underestimation, the gap at Γ between the Cu and P bands is closed. Figure 2. (a) Cu 2 P 7 unit cell, (b) Cu 3 P 11 unit cell, and (c) Cu 2 P 7 2 × 2 × 2 supercell in which the P−P connectivity is shown to highlight the P chains in the supercell structure. (d) Cu 3 P 11 2 × 2 × 2 supercell with the P−P connectivity to show P chains as in the Cu 2 P 7 supercell.  (ICSD 35282), and C2/m−Cu 2 P 7 (ICSD 35281) shows all three have a first peak between 2.20 and 2.24 Å, while CuP 2 has two peaks around 3.6 Å, where both Cu 2 P 7 and Cu 3 P 11 have one. All PDFs are artificially broadened with Gaussians of width 0.1 Å, and PXRDs are calculated using a Cu Kα source. (b) Simulated PXRD patterns of both Cu 2 P 7 and Cu 3 P 11 share peak positions, as is expected from their shared symmetries. The Cu 3 P 11 phase could be identified experimentally by the higher-intensity peaks at 2θ < 30°, including a distinct peak at 16°, not present in Cu 2 P 7 . Details of PXRD calculations can be found in the Supporting Information.
Chemistry of Materials pubs.acs.org/cm Article Many M 2 P phases (where M is a transition metal) have a structure similar to P6̅ 2m--Ni 2 P 60 and Fe 2 P, 61,62 in which the metal atoms sit in a cage of 3-fold coordinated P and 4-fold coordinated metal atoms. Fm3̅ m−Cu 2 P is most similar to the other Fm3̅ m TMPs, as it was derived from the prototype structure Fm3̅ m−Ir 2 P, and has a 4-fold coordinated Cu with an 8-fold coordinated P. The Cu 2 P band structure in Figure 6 is also similar to those of Ir 2 P and Rh 2 P. In Rh 2 P, there is a directionally opened gap of 1 eV above the Fermi level at the Γ point, not present in Cu 2 P or Ir 2 P (Cu 2 P, Ir 2 P, and Rh 2 P band structures calculated with spin−orbit coupling are given in Figure S6). . P2 1 /m−Cu 3 P found by the GA and Fm3̅ m−Cu 2 P from a swap with Fm3̅ m−Ir 2 P. Here Cu atoms are colored blue and P atoms are colored pink. P2 1 /m−Cu 3 P structure can be described as a stacking of Cu 2 P layers separated by Cu atoms. The stacking pattern, which is present in both structures, is indicated by the black dashed line surrounding the atoms in each structure. This is meant to guide the eye and show the similarity in the two structures. The unit cells of each structure are outlined with a thin gray box. Cu 2 P is predicted to be a stable two-dimensional (2D) phase, 56 which could be layered to produce the Cu 3 P phase shown in panel (a) here.  . Electronic band structure of Cu 2 P projected onto the Cu and P states and density of states projected onto the Cu s, p, d and P s and p states, for Cu 2 P using an energy cutoff of 500 eV with a 2π × 0.03 Å −1 kpoint spacing and C18 on-the-fly pseudopotentials. The projected band structure is produced by OptaDOS, 38 and band energies are calculated by CASTEP. 31 Chemistry of Materials pubs.acs.org/cm Article Both of these structures exhibit spin−orbit coupling due to their heavy metal ions, while Cu has negligible spin−orbit coupling effects. The Rh 2 P and Ir 2 P band structures are calculated, including spin−orbit coupling. Finite-Temperature Phase Stability. The temperaturedependent convex hull was constructed by calculating the finitetemperature Gibbs free energies by including vibrational effects at the harmonic level 63 of several related structures on or near the convex hull from Figure 1. All structures within 20 meV/atom of the hull at 0 K were included in the finitetemperature hull; these were Fm3̅ m−Cu 2 P, I4̅ 3d−Cu 3 P, Cmc2 1 −Cu 8 P 3 (the structure with 2 Cu vacancies from P6 3 cm−Cu 3 P discussed previously), CuP 2 , CuP 10 , Cu 2 P 7 , and Cu 3 P 11 .
The chemical potentials for this binary convex hull were Cmca-P (black phosphorus) and Fm3̅ m−Cu. Previously, Mayo et al. noted that the inclusion of semiempirical dispersion corrections for black phosphorus changed the energetics of the convex hull, 26 and therefore, it is not possible to combine optimized structures with and without dispersion corrections on the same convex hull. However, to obtain nonimaginary phonon frequencies of Cmca-P, it is necessary to account for dispersion. To account for this, the many-body dispersion correction (MBD) was applied during the geometry optimization and phonon calculation. 40 Using PBE, the distance between P chains in black phosphorus is 3.95 Å. By applying this correction, the P−P chain distance was reduced to 3.58 Å. To include the MBD black phosphorus on the convex hull in Figure 7a, we calculated the free energy of black phosphorus F(T) at a given temperature T as where H is the enthalpy without the dispersion correction, and ΔF MBD* (T) is the free-energy contribution at temperature T with the MBD dispersion correction, which includes the zeropoint energy. In this way, the energies of black phosphorus were referenced to the ground-state energy without dispersion. The SCAN functional accurately describes the phonon modes of black phosphorus without any added dispersion corrections (i.e., no imaginary modes are observed) and therefore is used as a comparison to the MBD-corrected PBE functional in Figure S7. Figure S7 shows that for any temperature T, both ΔF MBD (T) and ΔF SCAN (T) are on the same scale; only the zero-point energy is shifted (by 2.6 meV/atom) for the PBE + MBD calculation. Therefore, we expect the results of the PBE + MBD free energies to be comparable with nondispersion-corrected PBE free energies.
Using the MBD correction on black phosphorus in addition to the phonon modes of the previously mentioned phases of Cu−P, the hull in Figure 7a was constructed up to 600 K, above which no changes to stabilities are observed. A maximum value of 600 K was chosen so as not to approach the melting point of any phases, as the known phases of Cu−P typically have a melt between 800 and 1200 K. Furthermore, the harmonic approximation is a limited approach, and at higher temperatures, anharmonicity should be accounted for. Fm3̅ m−Cu 2 P remains on the hull at 600 K, suggesting that it could be synthesized at high temperature. The convex hull is a confirmation that the Cmc2 1 −Cu 8 P 3 phase formed from two Cu vacancies in P6 3 cm− Cu 3 P is the more stable phase at room temperature, as at 300 K, Cmc2 1 −Cu 8 P 3 is within 10 meV/atom of the convex hull, as shown in Figure 7b. In addition, the destabilization of I4̅ 3d− Cu 3 P at high temperatures, shown in Figure 7b, suggests that this phase is not experimentally realizable and provides an explanation as to why it has not yet been experimentally synthesized. We can clearly see that P6 3 cm−Cu 3 P is stabilized at higher temperatures, as shown in Figure 7b, in which it is within 10 meV/atom of the convex hull at temperatures higher than 450 K.
Previous work on the Cu 3−x P phases of P6 3 cm−Cu 3 P 47 confirms that the formation of two vacancies in Cu 3 P is energetically stabilizing. By enumerating all of the possible structures with two Cu vacancies using the vacancy enumeration procedure described in the Methods section, we have determined that the Cmc2 1 −Cu 8 P 3 phase with two Cu vacancies in the 6c Wyckoff positions is the lowest-energy vacancy phase. Given the large number of ways to introduce these vacancies into the structure, configurational entropy will further stabilize this phase at high temperatures. To fully understand the nature of vacancy formation in P6 3 cm−Cu 3 P, a full cluster expansion could be performed, which is beyond the scope of this paper.
Cu 2 P as a Li-Ion Battery Conversion Anode. Fm3̅ m− Cu 2 P was computationally predicted to be energetically stable both as a 2D material 56 and now, in this article, as a bulk phase. The previous sections predict the stability of Fm3̅ m−Cu 2 P at temperatures up to 600 K and characterize it as a metal with Figure 7. (a) Temperature-dependent convex hull with PBE xcfunctional and MBD correction on the black phosphorus chemical potential. (b) Distance from the hull for the four structures that were above the convex hull tie-line at 0 K. All structures except I4̅ 3d−Cu 3 P get closer to the tie-line as the temperature increases, suggesting that they are stabilized by temperature. This agrees with the experimental evidence for Cu 8 P 3 and P6 3 cm−Cu 3 P and suggests that Cu 3 P 11 may form experimentally.
Chemistry of Materials pubs.acs.org/cm Article dispersive bands and delocalized conduction states at the Fermi level. An intuitive choice of application for Cu 2 P lies in conversion anodes for Li-ion batteries, where previously both CuP 2 and Cu 3 P were used as anodes with gravimetric capacities between 300 and 800 mAh/g. 11−14 The crystal structure of P6 3 cm−Cu 3 P has a theoretical capacity of 363 mAh/g and experimentally has exhibited a range of capacities based upon the preparation method used. 13 The powdered Cu 3 P anodes prepared by Bichat et al. 13 ranged in the initial capacity from 272 mAh/g using high-temperature synthesis in a silica tube to 527 mAh/g using low-temperature solvothermal synthesis. In the solvothermal route, the Cu 3 P powders were prepared with copper chloride, water, and NH 4 OH with white phosphorus, which could have resulted in copper oxide impurities leading to the initial capacity that is above the theoretical capacity of crystalline Cu 3 P. Pfeiffer et al. synthesized Cu 3 P powder by a solid-state reaction with red P in an ethanol suspension and Cu foil, and the initial capacity of Cu 3 P was 415 mAh/g. 14 Energy-dispersive X-ray analysis showed that the stoichiometry was close to Cu 3 P (though not exact), suggesting that this initial structure could have been within the stoichiometric range of Cu 3−x P to achieve that initial capacity, in addition to the added capacity from likely oxide impurities.
In contrast, Cu 2 P has a theoretical capacity of 509 mAh/g, which is above that of graphite at 372 mAh/g. The metallic nature of Cu 2 P further enhances its use as a Li-ion battery anode, enabling fast electronic transfer through the electrode of the battery. In fact, Cu is already widely used as a current collector in contemporary Li-ion batteries, and previous studies on Cu 3 P nanorods suggest that Cu−P anodes create a synergistic chemical interface with the Cu-current collector, which promotes cyclability. 15 Furthermore, because of its comparatively lower P content, volume changes during cycling are reduced, and therefore, the degradation is likely to be less severe. 12 The volume expansion for a conversion anode with an overall conversion reaction where V(Cu a P b ) is the volume per formula unit of each phase in the conversion reaction. Using this equation, the volume expansion of Fm3̅ m−Cu 2 P is 99%. This is comparable to the calculated volume expansion of P6 3 cm−Cu 3 P, which is 86%, and far superior to the volume expansion of CuP 2 , which is 165%. The volume expansion for each binary Cu−P phase is shown in Figure S8, confirming that Cu 2 P has the lowest volume expansion of the four stable phases on the convex hull. Experimental reports on the cycling of ball-milled CuP 2 12 suggest that volume expansion occurs, as after cycling for 10 cycles, the capacity is reduced by 50%, although they give no estimate of the level of volume expansion in the cell. The expansion is partially mitigated through the use of nanostructuring, 11 which allows cycling for 200 cycles. However, there is still capacity fading in this case, which reiterates the need for a high-capacity conversion anode with low-volume expansion, so as to reduce the need for nanostructuring or other postprocessing techniques to mitigate volume expansion. As both Cu 3 P and Cu 2 P have lower predicted volume expansions, and the synthesized Cu 3 P shows no evidence of deleterious volume expansion, 13 it is likely that Cu 2 P would also have minimal volume expansion in the experiment.
Using the convex hull constructed in Figure 1 and the structures on the ternary hull of Cu−Li−P, a voltage profile was constructed from the DFT ground-state energies for both Fm3̅ m−Cu 2 P and P6 3 cm−Cu 3 P. All of the known ternary structures were included in this hull: P3̅ m1−Cu 2 LiP, I4/mmm− Cu 2 LiP 2 , Immm−Cu 4 Li 5 P 6 , and Cmcm−CuLi 2 P, as well as the binary Li−P structures Cmcm-Li 3 P, P2 1 /c-LiP, P2 1 2 1 2 1 -Li 3 P 7 , and I4 1 /acd-LiP 7 . A plane wave kinetic energy cutoff of 700 eV was used, and all structures on the hull were re-relaxed at this higher cutoff. The ternary hull is shown in Figure S9, in which the pathways from Cu−P to Li are also shown to depict how the voltage profiles for these Cu−P phases were calculated. The hull is shaded with a colormap to show the relative formation energy of phases on the hull, indicating that the Li−P phases have more negative formation energies (and thus create a deeper convex hull) than the Cu−P phases.
Although Cu−Li phases are predicted to be stable under the approximation of PBE, the formation energy of the predicted Cu 3 Li phase is only 26 meV/atom in the OQMD database 19,64 and no phases of Cu−Li are predicted at finite temperature in the experiment. 65 Furthermore, Cu is used as a current collector in Li-ion batteries, specifically for its properties in resisting Li intercalation, and dead Li is found during cycling rather than Cu−Li phases. 66 Therefore, no Cu−Li compounds were included in the convex hull.
There are three ternary compounds on the Cu−Li−P hull in Figure S9; these are I4/mmm−Cu 2 LiP 2 , Immm−Cu 4 Li 5 P 6 , and Cmcm−CuLi 2 P. Experiments suggest that a hexagonal LiCu 2 P phase forms 67 during cycling; however, the P3̅ m1−Cu 2 LiP (ICSD 659706) 46 phase of this structure is 39 meV/atom above the hull at a plane wave cutoff of 700 eV.
From this ternary hull, shown in Figure S9, the voltage profile shown in Figure 8 was constructed. This hull is calculated as usual, without incorporating vibrational effects at 0 K. As the Figure 8. Ground-state voltage profile for Fm3̅ m−Cu 2 P and P6 3 cm− Cu 3 P generated from the DFT ground-state structures of Cu−Li−P. The experimental profile was adapted from ref 67. The experimental onset voltage matches closely with Cu 3 P and shows a similar capacity to Cu 2 P, though this added capacity is likely due to oxide impurities in the experiment.
Chemistry of Materials pubs.acs.org/cm Article P3̅ m1−Cu 2 LiP phase suggested in the experiment 46,67 is 39 meV/atom above the convex hull, it cannot be in the voltage profile calculated in Figure 8. The 0 K voltage profile includes the I4/mmm−Cu 2 LiP 2 phase, which has been previously synthesized through a solid-state reaction, 68 and is a high T c pnictide superconductor. 69 The I4/mmm−Cu 2 LiP 2 phase has not, to our knowledge, been identified during cycling in Li-ion batteries previously. Both Cu 2 P and Cu 3 P have the same overall reaction mechanism given by eq 4, and the stable phases during the reaction from charging Cu 3 P are given in Table 2. These reactions show that Cu 2 P operates in a narrower voltage window than Cu 3 P and has a higher predicted gravimetric capacity. The predicted structures forming at each capacity and voltage are given in the 4th column of Table 2. Though Fm3̅ m−Cu 2 P undergoes the same lithiation process as P6 3 cm−Cu 3 P, Fm3̅ m− Cu 2 P has a higher capacity of 508 mAh/g and a higher average voltage of 0.86 V versus Li/Li + , while P6 3 cm−Cu 3 P has a capacity of 363 mAh/g and an average voltage of 0.91 V.
Each plateau in the ground-state voltage profile in Figure 8 represents a three-phase region of the ternary hull in which phases of Cu−Li−P are stable. Here, Cu 3 P was "stabilized" on the ternary hull by artificially excluding the CuP 2 −Cu 2 P 7 and Cu 2 P phases. This is an approximation of a convex hull in which Cu 3 P is on the tie-line, which does not affect the formation energy (and thus predicted voltages) of the other phases. The experimental voltage curve shown in Figure 8 from ref 67 exhibits a similar trend in phase transitions along the cycle as the theoretical curve for Cu 3 P.

■ CONCLUSIONS
Using four different computational crystal structure-searching techniques on the copper phosphides, several structures that lie close to the convex hull (within 20 meV/atom) were predicted, including Cm−Cu 3 P 11 , I4̅ 3d−Cu 3 P, and Fm3̅ m−Cu 2 P; the experimentally characterized P1̅ −CuP 10 , C2/m−Cu 2 P 7 , and P2 1 /c−CuP 2 were all on the convex hull tie-line. By calculating the phonon dispersion curves of all structures within 20 meV/atom of the Cu−P convex hull, we constructed a temperature-dependent convex hull that predicted Fm3̅ m− Cu 2 P to be stable up to 600 K, while I4̅ 3d−Cu 3 P was destabilized with increasing temperature. We have also shown that the Cmc2 1 −Cu 8 P 3 phase formed from two Cu vacancies at the 6c Wyckoff positions of P6 3 cm−Cu 3 P is stabilized with increasing temperature and is within 10 meV/atom of the convex hull above 300 K. Experimental diffractometry on single crystals of Cu 3 P suggests that the phase has a range of stoichiometries between Cu 2.6 P and Cu 2.8 P, 43,47 and Cu 8 P 3 (or Cu 2.67 P equivalently) is within these bounds.
In addition to confirming the stability of Cmc2 1 −Cu 8 P 3 , we also confirmed that Cm−Cu 3 P 11 remains metastable up to high temperatures, as shown in the temperature-dependent hull in Figure 7. While Cu 3 P 11 is unlikely to be used as a Li-ion battery anode, given its high P content, and therefore susceptibility to volume expansion, it could be a novel phase to consider within the Cu−P phase diagram. CuP 10 was identified experimentally by preparing Cu 2 P 7 in excess P; 70 given the structural similarity between Cm−Cu 3 P 11 and C2/m−Cu 2 P 7 shown in Figure 2, it is possible that Cm−Cu 3 P 11 is also formed in excess P. Using the PXRD patterns presented in Figure 3, it may be possible to distinguish the Cm−Cu 3 P 11 phase from C2/m−Cu 2 P 7 experimentally, from the change in the peak intensity at 16°and peak differences at 2θ values < 20°; however, further experimental analysis is likely required given the low intensity of this peak.
Finally, Fm3̅ m−Cu 2 P is the only phase identified through crystal structure prediction, which was found on the hull at 0 K and which remained on the convex hull at finite temperature, strongly suggesting that it is possible to synthesize Cu 2 P experimentally. Furthermore, its synthesis could provide a novel conversion anode, with favorable properties for Li-ion batteries. Hybrid functional calculations of the electronic properties of Cu 2 P predict it to be isostructural and qualitatively similar electronically to both Rh 2 P and Ir 2 P, which are also Fm3̅ m metals with dispersive bands at the Fermi level. This was confirmed using spin-polarized calculations, with both vector and scalar spin treatments, hybrid functional calculations using the HSE06 functional, and finally a projected band structure and density of states using PBE. This confirmation of the metallic nature of Cu 2 P using a wide range of functionals and spin treatments suggests that this could be a better choice for anode than Li−P, which are insulators with wide band gaps. Furthermore, the presence of such dispersive bands suggests high electron mobility within the anode, which would mitigate fast charge transfer between the Cu current collector and Liions. Finally, given its higher capacity (509 mAh/g) compared to Cu 3 P, Cu 2 P has potential as an experimentally realizable conversion anode, which has a capacity that is competitive with graphite, conductive to promote electronic transfer within the anode, and less vulnerable to degradation compared to high P content conversion anodes due to reduced levels of cyclic volume changes.

* sı Supporting Information
The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.chemmater.0c02054. The input and output files associated with our calculations have been deposited into the Cambridge Repository at https://doi. org/10.17863/CAM.54795 in the Cambridge Repository; and all analysis is on GitHub at harpaf13/data.copper-phosphides and Binder.
Details of the simulation of powder X-ray diffraction patterns and pair distribution functions, as well as the electronic band structure calculations using spin−orbit coupling (PDF)