New Insights into the Catalytic Activity of Cobalt Orthophosphate Co3(PO4)2 from Charge Density Analysis

Abstract An extensive characterization of Co3(PO4)2 was performed by topological analysis according to Bader‘s Quantum Theory of Atoms in Molecules from the experimentally and theoretically determined electron density. This study sheds light on the reactivity of cobalt orthophosphate as a solid‐state heterogeneous oxidative‐dehydration and ‐dehydrogenation catalyst. Various faces of the bulk catalyst were identified as possible reactive sites given their topological properties. The charge accumulations and depletions around the two independent five‐ and sixfold‐coordinated cobalt atoms, found in the topological analysis, are correlated to the orientation and population of the d‐orbitals. It is shown that the (011) face has the best structural features for catalysis. Fivefold‐coordinated ions in close proximity to advantageously oriented vacant coordination sites and electron depletions suit the oxygen lone pairs of the reactant, mainly for chemisorption. This is confirmed both from the multipole refinement as well as from density functional theory calculations. Nearby basic phosphate ions are readily available for C−H activation.


Introduction
The synthesis of today's bulk chemicals heavily relies on catalysis, reducing energy consumption and preventing waste in atom-economically tailored reaction sequences. [1,2] The design of appropriate reagents precisely adapted to the specific demandi so ne of the fundamentalo bjectives in chemistry.A lthought he exact synthetic targets might differ,t he aim is always the same:s elective replacement of individual atoms or residues within am aterialt of inetune its properties and by these means its reactivity or functionality.A samajor drawback, it is still not alwaysp redictable which alterations cause the desired effects. Hence the benchmarking of working modelsa longe stablished concepts is required. The properu nderstanding of the structure-property relationship of ac hemi-cal system is vital for the improvement of the materials development process. Electron-density (ED) investigationsp rovide an ideal tool for understanding these interactions. [3] Alreadyi n 1932,P auling [4] proposed an interrelation between the structure of ac ompound and its properties: "The properties of ac ompound depend on two main factors, the nature of the bonds between the atoms, and the natureo ft he atomic arrangement.
[…] The satisfactory descriptiono ft he atomic arrangementi nacrystal or molecule necessitates the complete determination of the positiono ft he atomsr elative to one another." From the knowledge of the distances at the atomic level and the arrangement in the solid phase, many properties both at the molecular and macroscopic scale can be deduced. However,t he most basic concepts, such as the chemical bond and reactivity,a re still strongly discussable.
Single-crystals tructural analyses based on the independentatom model only provide the positions of the centroids of the atomsa nd the distances between the atoms. In contrast, the multipole model from high-resolution datap rovides physically meaningful properties far beyonds imple geometrical considerations. In this paperw ei nvestigate the ED in catalytically active cobalto rthophosphate Co 3 (PO 4 ) 2 (Scheme 1) derived from both high-resolution diffractiond ata and density functionalt heory (DFT) calculations, to get some insight in the various properties contributing to that activity.N ot only the bonding is studied but also the ED distribution at the two different metals ites is determined andl inked to various substrates.

Results and Discussion
Cobalt orthophosphate crystallizes in the monoclinic space group P2 1 /n,w hereas in earlierp ublications [6,7] sometimes the other setting P2 1 /c was chosen.F urthermore, am etastable phase [8] of cobalt orthophosphate is known in the nonstandard setting P2 1 /b. Theu nit cell contains two formula units of Co 3 (PO 4 ) 2 .O ne of the two independent Co II atoms is sixfoldc oordinated in an approximately ideal octahedral fashion (Figure 1a). The second is fivefold coordinated,g iving as everely distortedt rigonal bipyramidal coordination polyhedron (Figure 1b). The octahedrally coordinated Co ion is O-corner-sharing connected to six different phosphate ions;t he Co inside the distorted bipyramidal polyhedron is only connected to four different anions ( Figure 2). Three of these are corner-sharing and one is edge-sharing connected to the metal. The asymmetricu nit contains only one half of the formula unit, hence the octahedrally coordinated Co atom is located on a symmetry centera nd the fivefold coordinated Co atom is on a general position. To judge whether the latter polyhedron is closer to as quare-pyramidal or to at rigonal-bipyramidals tructure, the t parameter proposed by Addison [9] was determined.
It is defined as t = (bÀa) = 608 where b is the largest and a the second largestl igand-metal-ligand angle. For an ideal trigonalbipyramidal geometry, t is unity and for an ideal square pyramidal geometry zero. Here, t adopts av alue of 0.67, which suggestsapredominatelyt rigonal-bipyramidal environment.
Detailed bond lengths and angles are provided in the Supporting Information. For the sake of simplicity,t he different Co atoms will from now-on be termed Co [6o] and Co [5by] ,r espectively,a ccording to Machatschki's terminology. [10] Figure 2a shows the single-crystal structure composed from alternating layers of phosphate tetrahedra colored in magenta and Co polyhedra in bronze. The latter are connected through corners and edges to their adjacent Co polyhedra (Figure 2b).

Topological analysis
To model the ED more appropriately than feasible from the independent-atom model (IAM), the multipole model (MM), based on the Hansen and Coppensf ormalism, [11] was employed. It allows for the description of aspherical ED by spherical harmonics, in which the total density is partitioned in a spherical core density,aspherical valence, and an aspherical deformation valencedensity.T odescribe the aspherical density, local coordinate systems are necessary.W eu sed DFT calculations, together with manual judgment, to determinet he pertinent local coordinate systems, because the actual coordination polyhedra around the Co ions in Co 3 (PO 4 ) 2 are distortedf rom the ideal polyhedra (for details see the SupportingI nformation).T he definition of the local coordinate systemsi sf ound in Figure3aa nd Figure3b. The thus derived ED can be investigated by ac omplete topological analysisa ccording to QTAIM, [12] based on the experimental diffraction data (MM-Expt, MM = multipole model)a sw ell as the structure factor derived from DFTcalculations (MM-DFT). In the Bader methodology,t he presence of bond-critical points (BCP) in the molecular graph is ap rerequisite for an interaction between two atoms. The features at the BCP are also an indication for the bonding type and the strength of ab ond, which might correlate with the distance. [12] Ac onsistent set of critical points was found for cobalt phosphate. They obey Morse's relation [13] nÀb+ +rÀc = 0, where n, b, r,a nd ca re the numberso fn uclei, bond, ring, and cage-critical points, respectively. Scheme1.Catalytic reaction of butan-2-ol with oxygen to the products butene and butenone. [5] Figure 1. Differente nvironments of a) Co [6o] and b) Co [5by] .The coordination polyhedron of phosphorus is colored in magenta andthat of Co in bronze.

Phosphate
Regardingt he nature of the PÀOb onds in phosphate, there is still al ong-standing discussion about the presence or absence of double-bond character,h ence also about the presence or absence of valence expansion at the phosphorus atom. Different theoretical calculations led to variousb inding models. Molina, Sundberg et al. [14] investigated the bondings ituation of severalt heoretically calculated and geometry optimized struc-tures by meanso fQ TAIM. They showed that the PÀOb ond in phosphane oxide is best described with ah ighly polar s-bond with ac onsiderable electrostatic contribution. The same is valid for iminophosphoranes. [15] Durrant [16] calculated the number of valence electrons at the phosphorus atom in PS 4 3À to be 7.94 anda tt he sulfur atom in sulfate to be 4.34. He, amongm any others, regarded the PÀOb ond in phosphate as very polar.F or potassium sulfate, the high polarity of the SÀO bond wasa lready established in an experimental charge-density investigation. [17] In this study,w ea lso show for phosphate that the bond between phosphorus and oxygen is ah ighly polar covalentbond with no double-bond character. [18] The evidence for this stems from the high ED and apositive Laplacian, the trace of the Hessian matrix of the second derivative of the ED function, at the BCP.A se xpected, for all four PÀOb onds a high ED was found at the BCP,w hich suggests ac ovalent bond ( Table 1). The values vary from 1.54 to 1.61 e À3 and are in ag ood agreement with the theoretical data. For the Laplacian at the BCP, positive values were found in the range of 2-5e À5 ,w hich are slightly highert han those from theory (0-3e À5 ). Given that the BCP is located in ar ampant edge of the Laplace function, even as mall variationi nt he position of the BCP will cause at remendous change in the numerical value and sometimes even ac hange of sign (Table 1, FigureS4-2, Supporting Information). The ellipticity describes the deviation of the charged ensity around ab ond path from an ideal cylindrical symmetry: e = l 1 = l 2 À1, where l 1 and l 2 correspond to the eigenvectorso ft he Hessian matrix.N ot surprisingly,w e found e values close to zero for all PÀOb onds, corresponding to s single bonds. FigureS5-1 (Supporting Information) depicts the distances of the BCP from the oxygen-atom positionf rom theory (MM-DFT) plotted above those from the experiment (MM-Expt).
In the experiment, the BCP is slightly more shiftedt owards the electropositive phosphorusa tom, indicating am arginally higher ionic character in the experimental than in DFT data. Cremera nd Kraka [19] proposed to use the total energy density H to classify bonds. For covalentb onds, the potential energy dominates, resulting in an egative total energy compared with closed shell interactions. In addition, as introduced by Espinosa et al., [20] the ratio of potential tok inetic-energy density j V j = G can be used to distinguish between shared, intermediate and closed-shell interactions. For ac losed-shell interaction, the ratio j V j = Ga tt he bond-critical point is less than unity,f or a shared interaction greater than two,a nd for an intermediate interaction between one and two. The ratios j V j =Gfor the PÀ Ob onds ranging from 1.82 to 1.93 [1.90-2.01] are much more on the side of the shared interaction. If the ratio of the kineticenergy density per electron density at the BCP is greater than unity,this is indicative for polarity in ab ond.
For the PÀOb onds in phosphate, we found a G(r bcp ) = 1(r bcp ) ratio greatert han unity (1.07 to 1.25). Maxima of the negative Laplacian represent Valence-Shell Charge Contributions (VSCC). VSCCs in the nonbonding regions are indicative for lone pairs, frequently providing the link to the classical Lewis model and the VSEPR theory. [21] For the phosphate oxygen atoms, four tetrahedrally oriented VSCCs would be expectedd ue to the sp 3 Figure 3. a) and b) local coordinate systems for Co [6o] and Co [5by] ;c)and d) showt he deformation density with the iso-surface valuesofAE 0.65e À3 for negative charge accumulation in light blue and positive charge in orange; e) and f) showt he coordination sites;g )a nd h) depict the Laplacian with iso-surfaces of 240 e À5 (red)a nd À1250e À5 (dark blue) for the experimental, i) and k) theoreticaldata with iso-surfaces of 260 e À5 (red) and À1260 e À5 (dark blue). hybridization. [14] However,t he number of VSCCs could neither be resolved from MM-Expt nor from MM-DFT.Ase xpected, one bondingV SCC is located near the PÀObond axis. The VSCCs in the nonbonding region pointing away from that bond are fused to at oroidal orientation. Even with full refinement of all feasible multipole parameters from theoretical data, it was impossible to resolve them. We refrained from that endeavor because it clearly indicated overfitting. [22] Local mirror symmetry had to be taken into account. Additionally,t he multipole population at all oxygen atoms had to be set to be the same.

CoÀOPolyhedra
Compared with the covalent PÀOb onds, the CoÀOb onds show as ignificantly lower ED in the range of approximately 0.32 to 0.52 e À3 at the BCP.T he Laplacian at the BCPi ss ignificantly more positive and ranges from 4t o1 0e À5 .T he ED increase correlates with the reduction in distance. In the theoretical data set, 1(r BCP )i ss lightly largerf or all CoÀOb onds, showing al arger diversity at the Co [6o] ÀOb onds than the Co [5by] ÀO bonds. The found trend also reflects how the shared character increases with the decrease of the coordination number. Except for one, all Co [5by] ÀOb onds are shortera nd have a greater value for 1(r BCP )t han the Co [6o] ÀOb onds ( Figure S5-2, Supporting Information). This is also valid for many other MÀO bonds. The value for r 2 1(r BCP )i ncreases with the decreasing bond distance, which has also been observed for other transition-metal-oxygen interactions and nonmetal-oxygen interactions. [23] Comparing the differenceb etween the experimental and the theoretical data set, the values for r 2 1(r BCP )i nt he CoÀ Ob ond are almost identical. Ther atios j V j = Gf or the CoÀO bondso f0 .97-1.05 [1.02-1.08 from theory]i ndicate rather a closed-shell interaction. Ah int towards slightly sharedi nteraction in the CoÀOb ond is the negative total energy density, H (r BCP ).

Bader charges
The Bader charges (Table S4- [5by] ). As expected, the positive chargef rom both approaches is slightly higherf or the six-than for the fivefold coordinated metal atom.

Laplacian
The The vertices of the charge-depletion cube do not fully coincide with the ion-ligand bonds. In contrast, from the experimental multipole refinement (Figure 3g), there are only six charge depletions, that form an octahedron with vertices approximately along the six ion-ligandbonds.

Density of states and d-orbitalp opulation
DFT calculations found the lowest energy when all Co ions were in high-spin configurations. Figure 4s hows the calculated electronic density of states,D OS, for Co 3 (PO 4 ) 2 .T he high-spin complexes cause the spin-up (majority spin) and spin-down (minority spin) components of the DOS to be quite different from each other.T he spin-down DOS features well-localized Co d-states near the valence-bandm aximum (as seen from the partial DOS projected onto the individual Co ions, solid orange line and dashed brown line in Figure 4), with ao ne-electron dd-gap calculated to be around2 .5 eV.T he d-d transition energy is largely consistentw ith the experimentally measured absorption maximum for Co 3 (PO 4 ) 2 at 590 nm (2.10 eV); [24] the HOMO-LUMO gap is here assumedt ob eo nly af irst-order approximation to the excitation energy.I nc ontrast, the spin-up valence band has significant contributions from both of the two inequivalent Co ions in the range from À6t o0eV (where the zero level is set to the Fermi energy at 0K,E F ).
The d-orbital populations were estimated using the local coordinate systems around Co [5by] and Co [6o] from the multipole model.T he DOS was projected onto the individuald -orbitals and the populations were found by integrating the DOS up to the Fermi level. The resulting d-orbital populations are given in Ta ble 3a sp ercentage of the total d-populations on the respective atoms. The table also gives the estimated d-orbital populations from multipole refinement on the DFT-calculated structure factors( MM-DFT), as well as from the experimentally obtained structure factors (MM-Expt). Table 3s hows overall satisfactory agreement between both populations, confirming that the high-spin state for both Co ions found in the DFT calculations also corresponds to the experiment.
The relative populations behave accordingly to simple crystal-field theory:f or the distorted trigonal-bipyramidal complex   ,w e note that the populations obtained from the two DFT-based methods have an oticeably higher population on d xy (25 %) as compared to d xz (20 %), whereas this is not the case from the experimental data (in which the orbitals have roughlyt he same fractional population of 22 %). We explaint his discrepancy as follows:g iven the 0, d yz and d xy orbitals at Co [6o] are occupied by only five electrons in a high-spin state, there may be severaln early degenerate solutions to the Schrçdingere quation and the presentD FT calculations have converged to one of them.T oo btain at ruly accurate electronic description of the Co [6o] ion, it might therefore be necessary to make use of more advanced computational methods, such as multi-referencem ethods. Still, even with this limitation for Co [6o] ,t he agreement between the DFT calculations and the experiment is very good for Co [5by] .

Maximally localized Wannier functions (MLWFs)
To visualize real-space analogs of the electronic bands from the periodicD FT calculations,w ec alculated maximally localized Wannier functions centered on the Co atoms. The resulting MLWFs had clear connections to the d-orbital populations, as wella st he locationso ft he negative charge-depletion regions and VSCCsa round the Co ions. Given that all five spinup d-orbitalsa re occupied, they form ar oughly spherical ED. The two spin-down orbitals add to the ED in particulard irections, as can be seen from ac omparison between the spindown MLWFs and the VSCCs. Figure 5a depictst he negative charge-depletion points (red cube) at Co [5by] .F igure 5b-c show the spin-down MLWFs, pointingt owards the directions of the VSCCs.F igure 5d-h show the fives pin-up MLWFs centered on Co [5by] .I nterestingly, there is ac lear connection betweent he charge-depletion points and the orientation of the spin-up MLWFs. In particular, the two MLWFs in Figure5d-g, which have shapes that are reminiscent of atomicd z 2 orbitals, have one of their lobes pointingt owards the charge-depletion point along one of the equatorial CoÀOb onds, and the other lobe pointing towards one of the charge-depletion points that does not lie along any CoÀOb ond. In fact, this latter typeo fc harge depletion points towardsasmallv oid in the crystal. The MLWFs for Co [6o] are provided in the Supporting Information.
We also consider the charge-depletionr egions and MLWFs aroundt wo edge-sharing Co [5by] ions, and correlate these regions to the different orientation at variousC o 3 (PO 4 ) 2 faces. The two (011) and( 110) faces of cobalt orthophosphate expose nearby Co [5by] ions. Figure 6a shows as ide view of the Co 3 (PO 4 Figure 5d have lobes pointingout of the face, towards negative charge-depletionregions.C onsidering cooperativity between two nearby Co ions as for example, stabilizing low oxidation states [25] or molecular  conductors of electricity, [26] this pictures would suggest the crystal plane (011) as clearly catalytically more active than any other.

Consideringc atalyticabilities from ED
Phosphates,w hich contain cobalt and other transition metals, are knownt oc atalyzeo rganic reactions and are therefore in the spotlight of current research. The oxidation of styrene, for example, was realized using an ickel-cobalt catalyst. [27] The reductiono fp-nitrophenol to p-aminophenolc an be effectively catalyzed by nanosized solid solutions of copper/cobalt molybdate and chromium phosphate. [28] For ad etailed study of the cobalt-substituted calcium phosphate catalyst,w hich catalyzes oxidative dehydration and dehydrogenation,L egrouri et al. [29] used propane-2-ol as ab enchmark substrate for characterizing acid and basic properties of solid-state catalysts (Scheme 1, see above). The maximum catalytic effect was found at ac obalt content of 0.32, which also represents the maximum solubility of Co in calcium phosphate. The reactivity fore thane, propane, butan-2-ol, and propanol was studied by Aaddane et al. [5] and plausibility for the reactivity was provided.
Considering the reactiono fb utan-2-ol with oxygen to butene and butenone, at least one of the two reactants must be chemisorbed at the catalyst's surface and thus be activated. If the butanol is adsorbed over severala toms on the surface due to its size, the negative oxygen atom most likely coordinates through the two oxygen lone pairs to an electronically depleteds ite of cobalt atoms and the positivelyp olarized methyl hydrogen atoms to an earby basic phosphate site. If we assume differences in the heterogeneousc atalytic ability of the solid-statef aces in cobaltp hosphate we need to establish criteria to identify effective catalytic sites. [30] Given that the crystal structure of cobalt phosphate shows two differently coordinated metal sites, the first question to be asked is whether both Co atomsa re equally active, and if not, what causes the differences. If the plain coordination number of the cobalt atoms is considered, the saturated Co [6o] most likely is not tremendously active as it lacks av acant site for the substrate to interactw ith. Figures 3g and3hc learly show that there are no free sites at Co [6o] ,w hereas Co [5by] holds three of them. However,i th as to be kept in mind that Co-atoms exposed at the surfacei nevitably render free coordination sites. Depending on how exposed they are on the surface, the number of free coordinations ites also varies, which affects the probability of chemisorption. However,s tatistically speaking, Co [5by] would still have to have more free coordination sites than Co [6o] at the surface, because they already have three more vacant sites in the bulk. In addition, the present distortion of the geometry will promote catalysis. [31] When for example, CÀCb onds are activated in as ubstrate by rhodium and nickelc atalysts, it is proventhat this is astructurally sensitiver eaction that takes place on coordinatively unsaturated atom pairs of the d-metals on corners or edges of minute crystallites. [2] Thinking in terms of concepts of homometallicc ooperativity, [32] short distances between the metal centers are needed. In the crystal structure of cobalt phos-phate, there are such short distances between two fivefold coordinated edge-sharingC o-ions( 3.032 ) ( Figure 2b)a nd between af ivefold and as ixfold Co-ion (3.145 ). According to the Bravais,F riedel,D onnay, and Harker theory, [33] Figure 7 shows the four most likely low-index crystal faces with Co···Co cooperative sites. The layer is chosen such that no covalent bonds are broken in the phosphatea nd cobalt atoms appear on the surface. The positions of electronically depleted areas found from the multipoler efinedL aplacians are superimposed onto the surface of the cobalt atoms as black areas. At the (011) face only Co [5by] and Oa toms are visible,w here the Co atomso ccur isolated or in pairs with as hort interatomic distance of 3.032 .T he Co [5by] atoms face each other in pairs, each exposing three and four vacant coordinations ites. The pairs at the edges even have up to six vacant coordination sites per Co atom. They are perfectlya rranged to serve as adsorptions ites for the two oxygen lone pairs of the butanol molecules, which then s-bridges two Co atoms.
Looking at the other lattice planes from as imilar point of view,p lane (110) shows similar properties to (011), but all Coatoms are slightly deeper embedded in the surfacet hus they are much more difficult to be reached by the butanol oxygen atom. At the (10-1) face exclusively five-and sixfold coordinated cobalt atoms are exposed. Due to the solemnly electronically depleted sites at this surface, the chemisorption of am olecule with both, electron-rich and electron-poor regionsl ike butan-2-ol is much more difficult (Figure 7).The slightly positive s-hydrogen atom H 3 (Scheme 1) would probably not find a more basic nearbyp hosphate site to be activated. The (101) plane also displays Co-atom pairs with short metal-metald istancesb ut they occur between five-and sixfold coordinated Co atom. Furthermore, they are not as well arranged as at the (011) plane.Hence, we conclude that (011) plane is the catalytically most active. Finally,w en ote that the above arguments were made based on the ED for the bulk crystal. To verify that the predicted adsorptiono ff or example butan-2-ol onto the Co [5by] at the (011) surfacei sf avorable, we used DFT to calculate the adsorption energy of butan-2-ol ontos everal different adsorption sites for two different surfacet erminations of the (011) surface and found the predicteda dsorption site to be indeed the most favorable (see the Supporting Information) as expected, confirming that the bulk ED can be used as at ool to qualitativelyp redict adsorptionb ehavior at the surface. Nevertheless, the surface structure of ar eal catalystm ay deviate from bulk-truncated cuts, and so based on this analysis we cannotf ully exclude the presenceo fa dditional relevant active sites and surface facets. However,i nsofara st he catalyst surface structure resembles any of the low-index surfaces, the (011) surface hast he right structuralfeatures andE Df or catalytic conversion.

Conclusions
An extensive characterization of cobalt phosphate using topological analysis and theoretical methods fore xperimental and theoreticale lectron-density distribution was performed. The PÀOb onds were identified as polar single bonds in which the BCP is considerably shiftedt ot he phosphorus atom. The CoÀO bonds can be described as closed-shell interaction. The calculated MLWFs have ac lear connection to the d-orbital occupation, as wella st he position of the charge-depletion centers. In the solid state, the two independentc obalt atomse xhibit different numbers of charge depletion sites in the valence shell. Co [5by] has, in addition to the charge depletions directed towards the negativelyc harged phosphate oxygen atoms,t hree additional vacant coordinations ites. This suggests that the nucleophilic attack on Co [5by] is much more likely.T hese depletions, congruent with the unoccupied or partially occupied dorbitals, were analyzed on the face of differents olid-state surfaces. It hasb een shown that the (011) face is the most suitable for catalytic oxidative dehydrogenation or dehydrogenation, showingC o [5by] ···Co [5by] cooperativity.T herea re sufficient acidic and basic sites for the chemisorption of butan-2-ol available. On this face Co [5by] ions and phosphate ions are the most accessible;t herefore, this surface is best equipped to activate both, the CÀHa nd OÀHb onds. Furthermore, the lobes of the MLWFs of both cooperative Co [5by] atoms face each other to address both oxygen lone pairs of the substrate. The investigations show that the concerted but independente xperimental and theoretical determination of charge-density distribution is av aluable tool to rationalize catalytic active materials.

Experimental Section
High-resolution structure determination:C obalt orthophosphate was purchased as powder from Alfa Aesar.G iven that it is not soluble in any conventional organic solvents, it was crystallized from the melt. The solid was heated to 1180 8Cu nder air atmosphere in ac hamber oven, then cooled to 500 8C, tempered for one hour and finally cooled down to room temperature. The crystallization approach known from the Ref. [7] in slightly modified form was also tested, but did not produce better single crystals. Av iolet single crystal of 45 90 100 mmi ns ize was mounted at room temperature on the goniometer using inert perfluorinated polyether oil and cooled to 100 K.
High-resolution X-ray diffraction data were collected on aB ruker SMART APEX II diffractometer based on D8 three-circle goniometer system using an Incoatec microfocus AgKa source (ImS) and Incoatec QUAZAR mirror optics. [34] AB ruker APEX II CCD detector was used to record the diffracted intensities at l = 0.56086 .D ata reduction was carried out with SAINT [35] (version 8.30C) from the APEX2 [36] program package in which the integration box sizes were refined for every run using as tandard procedure. The data were scaled and equivalent reflections were merged with SADABS [37] (version 2016/2). An empirical absorption correction for as trong absorber was also performed with SADABS. Afterwards the structure was solved with SHELXT [38] using the direct methods and refined by full-matrix least square against F 2 using SHELXL [39] (version 2018/1) by means of the graphical user interface SHELXle. [40] To pological analysis according to the Quantum Theory of Atoms in Molecules [12] (QTAIM) was performed using the XDPROP and TOPXD programs included in the XD package. [41] Crystallographic details are provided in Ta ble S1-1, Supporting Information.
DFT calculations:P eriodic DFT calculations were performed on the experimentally obtained Co 3 (PO 4 ) 2 crystal structure using the WIEN2k program package version 17.1. [42] The composition of the conventional unit cell is Co 6 (PO 4 ) 4 .S pin-polarized calculations were performed using the PBE + Uf unctional (U eff = UÀJ = 3eV). The value of U eff was chosen on the basis of previous studies for cobalt oxides and hydroxides. [43] Although the value of U eff affects the calculated (one-electron) band gap and dd transition level, we have verified that other choices (U eff = 1eV, U eff = 4eV) do not affect the qualitative results of this study,i np articular as pertains to the shape and orientation of the maximally localized Wannier functions [44] (MLWFs). The calculations were performed using the linearized augmented plane wave (LAPW) method with RK max = 8.5. The Brillouin zone was sampled with a1 9 11 13 k-point grid, corresponding to 744 irreducible k-points;t he total energy was then converged within 3meV as compared to a9 5 6k-point grid. The lowest-energy structure was found when all Co ions were in high-spin configurations. We do not consider any magnetic ordering or relative alignment of spins, because the experiments in the current study were conducted at T = 100 Ka tw hich temperature the compound is paramagnetic. MLWFs were calculated using the WIEN2k Wannier interface to the wannier90 program [45] (i)byp rojecting the 12 highest occupied spin-down bands (two occupied bands per Co ion) onto the Co d-orbitals, and (ii)byp rojecting the highest 62 occupied spin-up bands onto the 30 Co d-orbitals in the conventional unit cell. We attempted several different rotations for the initial projections, but found that all of them converged to the same final result. The MLWFs were plotted using VMD. [46] The X-ray structure factors used as input for the multipole model were calculated for reflections up to sin q = l < 1.61 À1 using the WIEN2k lapw3 program. We also performed periodic surface DFT calculations using Quantum ESPRESSO, [47] and molecular DFT calculations using ADF. [48] The details of those calculations are given in the Supporting Information.