Magnetic Energy Landscape of Dimolybdenum Tetraacetate on a Bulk Insulator Surface

The magnetic states and the magnetic anisotropy barrier of a transition metal molecular complex, dimolybdenum tetraacetate, are investigated via density functional theory (DFT). Calculations are performed in the gas phase and on a calcite (10.4) bulk insulating surface, using the Generalized-Gradient Approximation (GGA)-PBE and the Hubbard-corrected DFT + U and DFT + U + V functionals. The molecular complex (denoted MoMo) contains two central metallic molybdenum atoms, embedded in a square cage of acetate groups. Recently, MoMo was observed to form locally regular networks of immobile molecules on calcite (10.4), at room conditions. As this is the first example of a metal-coordinated molecule strongly anchored to an insulator surface at room temperature, we explore here its magnetic properties with the aim to understand whether the system could be assigned features of a single molecule magnet (SMM) and could represent the basis to realize stable magnetic networks on insulators. After an introductory review on SMMs, we show that, while the uncorrected GGA-PBE functional stabilizes MoMo in a nonmagnetic state, the DFT + U and DFT + U + V approaches stabilize an antiferromagnetic ground state and several meta-stable ferromagnetic and ferrimagnetic states. Importantly, the energy landscape of magnetic states remains almost unaltered on the insulating surface. Finally, via a noncollinear magnetic formalism and a newly introduced algorithm, we calculate the magnetic anisotropy barrier, whose value indicates the stability of the molecule’s magnetic moment.


Introduction
Single-molecule magnets (SMMs) have been the objects of intense investigation for almost three decades. SMMs carry a magnetic moment (M) at the molecular level and are characterized by a large magnetic relaxation time at low temperatures and by the presence of a magnetic hysteresis loop. They differ from magnetic nanoparticles where M results from a 3D arrangement of ≈1000-100,000 atoms (superparamagnets), as SMMs are comparatively much smaller. Their behaviour also differs from the one of paramagnets, where the magnetic moments are carried by single atoms. Thus, as they are small and not 3D entities, but still composed of many atoms, SMMs are said to have a superparamagneticlike behaviour [1].
The first molecule identified as a single magnet was a dodecanuclear Mn complex exhibiting a magnetic anisotropy energy [1] (MAE), i.e., an energy barrier separating two anti-parallel magnetic moments along an "easy" magnetization axis [1,3,12] (EMA).
Thermal fluctuations allow these systems to overcome the magnetic energy barrier. However, thermal-assisted resonant tunnelling also allows to switch between magnetic ground states, as observed by the presence of steps in SMMs hysteresis loops [25]. As it is often difficult to observe these steps, a more practical signature of SMM features is given by the frequency dependence of the magnetic susceptibility imaginary part when employing external ac magnetic fields [20].
Important challenges for SMMs applications have been the possibility to operate at blocking temperatures (T B ) above liquid-helium and, currently, above the boiling point of nitrogen (77 K). Recently, a magnetic bistability up to T B = 60 K [26] and even to T B = 80 K [27] was found in different SMMs containing dysprosium, which represent dramatic improvements over the T B ≈ 3 K of the archetypal [Mn 12 ] compound. In order to increase coherence times, encapsulation techniques with endohedral fullerenes have been employed, with the fullerene cage used to shield external sources of decoherence affecting the magnetic ions spins [6].
Since the early 2000s, SMMs synthesized in the gas-phase and as bulk materials have been deposited on substrates [9,11,28], such as Au(111), where SMMs self-assemble, forming molecular chains on terraces [29]. Examples include on-surface lanthanide and transition-metal phthalocyanines and porphyrins, which showed a high magnetic tunability by changing ligands, metal-centers, end-groups and underlying substrates [11,28]. These complexes exhibit a high thermal stability and versatility [11], and besides being synthesised in gas-phase as catalysts [30], they were sandwiched on thin films [31], deposited on metallic substrates [32] or used as thin films on silicon surfaces [33,34]. Metallophthanocyanines magnetic properties on substrates exhibit both a ferromagnetic-like, collective behaviour (due to exchange interactions) and an individual, SMM behavior with no collective order [28,35].
A relevant issue is that SMMs, when deposited on conductive nonmagnetic substrates, might loose their magnetic properties, making their magnetic remanence vanishingly small and exhibiting narrow hysteresis loops [43]. In general, the magnetization and the hysteresis of a molecule depend on the chemical environment surrounding the magnetic ion(s). This includes the crystal field of the molecular complex but also the presence of the substrate, via electronic screening, molecule-surface electronic hybridization and charge transfer [44,45]. Substrates-molecules combinations show cases where the molecular magnetic moment is unaffected or partially/completely quenched. Interestingly, sometimes this is not only related to the ion magnetic moment, but to an overall change in magnetization, where also other atoms (e.g., N or C) contribute to the total moment [46].
Strategies to preserve the SMMs properties include the decoupling between the SMM and a metallic surface, through the insertion of nonmagnetic insulator films [43]. Clearly then, the successful attempts in keeping the properties of SMMs when adsorbed on insulating films such as MgO [43,47], pose the question as whether the adsorption on bulk insulating surfaces represents a valid alternative to preserve the molecular magnetic properties of SMMs.
Recently, a coordination complex, dimolybdenum tetraacetate Mo 2 (O 2 CCH 3 ) 4 , denoted in the following as "MoMo"), was deposited on the bulk insulating calcite (10.4) surface (CaCO 3 ) [48] and other metallic and insulating surfaces [49]. Despite the wellknown challenges in having stable molecular self-assemblies on insulators, MoMo was observed to form stable and locally regular arrays at high coverage [48], adopting an adsorption geometry on calcite that allows the molecule to be immobile at room temperature, as also confirmed by ab initio DFT calculations [48].
Realizing stable and ordered molecular structures that spontaneously assemble on bulk insulators is an important step forward in the field of molecular electronics [50][51][52][53]. The ordered MoMo arrays open up this possibility, with the extra feature of containing transition metals, regularly spaced by the acetate groups belonging to the molecule [48]. Due to the TM presence, a relevant question is whether these complexes can carry magnetic moments, resulting from stable or meta-stable electronic configurations. A related issue is whether MoMo can be assigned SMM features. Di-molybdenum complexes are known to feature triple or quadruple bonds between the Mo centers (depending on the strength of the π-donation from their ligands) [54,55] which imply the pairing of the involved electrons and inhibit the formation of magnetic moments. To the best of our knowledge, however, the electronic structure of the tetraacetate complex considered in this work was never directly addressed; it received, instead, a precise structural characterization as reported in Ref. [56].
Magnetic networks on bulk insulating surfaces are expected to have advantages over metallic substrates: for example, magnetic moments and hysteresis cycles might not be affected or quenched by the substrate-molecule (S-M) electronic hybridization [43,47,50,51]. A second advantage is that due the weaker reactivity of insulators as compared to metals [50][51][52][53], the S-M interaction prevents drastic changes in the structural, and thus electronic and magnetic properties of the molecule [50,51]. Thus, provided that an S-M anchoring mechanism is established (as certainly is in MoMo [48]), promising SMM features of a molecule in the gas-phase can be, in principle, easily transferred on insulators.
In this work we perform a computational study of the electronic structure of MoMo with a particular focus on its magnetic states and properties. Results obtained for the molecule in gas phase are contrasted with those for the system adsorbed on the calcite (10.4) insulating surface. In spite of the stability of the expected multiple Mo-Mo bonds, our goal here is to capture the energy landscape made of (meta-)stable magnetic configurations, and to evaluate the energy cost to deviate the magnetization of the system in particularly interesting FM states from its ground-state orientation. This is done here for the isolated molecule and, for a selection of interesting magnetic states, also for the system adsorbed on calcite (10.4), using the adsorption geometry determined in Ref. [48]. After modeling the system with the GGA-PBE functional, in order to improve the description of localized Mo d electrons, we resort to Hubbard corrected functionals, able to stabilize several meta-stable magnetic states. The gas phase results are then contrasted with the ones on the calcite surface. Finally, we evaluate the MAE which gives an indication about the possible SMM properties of this system. The long-term objectives are to understand whether a magnetic state can be stabilized in gas phase and on the insulating surface and to assess, through an evaluation of its MAE, whether it could be considered a viable prospect SMM.

Results
The electronic configuration of the Mo atom is [Mo] = [Kr]4d 5 5s 1 . However, based on our calculations, the Mo atoms in this molecule retain about 4.5-4.6 electrons on its d states, four of which are almost fully occupied (according to Refs [57,58], having four almost occupied d states corresponds to a 2+ oxidation state for the Mo ion). With a less-than-half-filled d shell, the maximum magnetization achievable is thus 4.5 µ B per Mo atom; however, this limit is never reached and calculations predict a maximum atomic magnetization of about 3.4 µ B per Mo atom. In molecular systems with multiple TM centers, the magnetic state also depends on the relative alignment of atomic moments, and the possible interactions between them and/or with their chemical environment. The MoMo structure (Figure 1a), precisely characterized in Ref. [56], shows that each Mo is coordinated with four oxygens, forming an almost planar square centered on the TM, which sits slightly off plane. The two MoO 4 plaquettes are overlayed on top of each other and linked by the external C-CH 3 groups. Importantly, there is no oxygen directly bridging the two Mo atoms nor charge transfer emerges between them, which prevents to use the super-exchange theory to draw direct conclusions on the coupling between their magnetic moments. Similar compounds featuring Mo-Mo central units are reported to form quadruple bonds between the Mo centers [55,56]. In these systems, however, the planar plaquettes formed by each Mo center and its four ligands do not feature external bridges and are only linked together by the Mo-Mo bond. At variance with these systems, in our system there are four C-CH 3 units that externally bridge the two MoO 4 units. The structural stability of these external links allows to promote the Mo d electrons from bonding to antibonding states, without compromising the structural integrity of the molecule. This is, in fact, an important feature that allows to avoid the spin pairing of the Mo d electrons and to consider magnetic configurations. In order to study the MoMo magnetic configurations and to investigate whether (meta)-stable magnetic ground states (GS) can be achieved when the molecule is adsorbed on the calcite surface, we will first start to characterize the isolated molecule. Structural optimizations in the gas-phase will be performed for different magnetic configurations to determine their relative stability. Next, MoMo will be deposited on calcite (10.4) and re-relaxed for some of the lowest-energy magnetic configurations. Comparing the relative stability of various MoMo magnetic configurations in gas-phase and on-surface will indicate the role of the surface in stabilizing these states. In this Section, we present DFT calculations performed with spin-resolved generalizedgradient approximation (σ-GGA) for the exchange-correlation (xc) functional, with the PBE prescription [59]. A selection of input and output files are available as Supplementary Materials.
Using this functional, independent of the starting value of the Mo's spin polarization, we always obtain a nonmagnetic (NM) ground state, with vanishing magnetic moments on each Mo center (i.e., an equal distribution of d electrons in the two spin channels). This result is consistent with the formation of a quadruple bond between Mo centers; in fact, a Löwdin population analysis of their d states reveals that four of them are approximately half-occupied for each spin in each Mo, which reflects the pairing of the electrons spin on the bonding states being formed. Consistent with the available literature [55] the d states involved are: the z 2 (σ bond), the xz and yz (π bonds) and the xy (δ bond), where the z axis is taken along the Mo-Mo axis. However, the well known tendency of approximate xc functionals to over-delocalize d electrons, might impair the stabilization of magnetic ground states that could result from the prevalence of the exchange interaction on atomic orbitals (first Hund's rule). To further probe the system's behavior, as an attempt to force a magnetic configuration, we imposed an imbalance between the two spin populations, by fixing the occupations of the Kohn-Sham (KS) spin states during the electronic selfconsistent solution. In this way, several magnetic states could be achieved, with different values of the total magnetization M. During the procedure, we relaxed the atomic positions for each value of M. Table 1 reports the total energy E (relative to the nonmagnetic GS) for all the achieved states, along with the total magnetization M, the distance between the Mo centers and their atomic magnetization (while M is the unbalance between the number of electrons in the two spin channels and has an integer value in a system with gap, n ↑ -n ↓ are obtained from projecting KS states on atomic d orbitals, thus yielding fractional numbers. Note, however, that the sum of the atomic magnetization approaches the value of M). Table 1. GGA-PBE results for the MoMo molecule in the gas-phase. For each state, we report the total energy E (relative to the nonmagnetic ground state energy), the total magnetization M, the Mo-Mo distance d Mo−Mo and the magnetization on the d states of each Mo, n ↑ -n ↓ . Double entries in n ↑ -n ↓ refer to the Mo atoms, when they present different magnetizations. In the first row, the numbering at the end of the label refers to the corresponding value of M. As evident from the table, the Mo-Mo distance in the NM GS is in perfect agreement with the experimental value of 2.09 Å [56]. The table also illustrates that our approach allowed us to achieve some M > 0 ferromagnetic (FM) states.

NM
As revealed by a Löwdin analysis of d states populations, these FM states, of increasing energy, are obtained by promoting electrons from the Mo-Mo bonding to the corresponding anti-bonding states of progressively higher energy. This gradual disruption of the Mo-Mo bonds implies the corresponding increase of the Mo-Mo distance without, however, compromising the integrity of the molecule, which is still held together by the external C-CH 3 groups. From a technical standpoint, the FM states were obtained through constraining the KS occupations and by re-optimizing the molecular structure each time. However, with the exception of the M = 6 µ B state, they are all excited states; in fact they did not result in stable configurations and a structural relaxation with no constraint on the occupations evolves the system back to its non-magnetic GS. This suggests that the self-consistent solution of constrained KS equations did not succeed in creating a self-consistent potential able to stabilize the corresponding M state. As for the M = 6 µ B state, it represents a meta-stable configuration where the system remains trapped even after constraints on KS occupations are lifted. The energy of this state, however, is significantly (4.08 eV) higher than the NM ground state and this thus represents the only stable state for MoMo when modeled with GGA-PBE functionals.

Results from Hubbard-Corrected Functionals
Standard xc functionals, as the GGA-PBE, are affected by electronic self-interaction and miss most of the effects of strong correlations. Localized d or f electrons are mostly affected by these flaws as they are typically over-delocalized, especially in presence of degeneracy. Among the consequences, one of the most typical is the suppression of magnetism. In order to ascertain whether some results in the previous Section are in fact artifacts of the GGA-PBE functional (in particular, the stable Mo-Mo multiple bond) and whether a magnetic GS can actually be stabilized by alleviating the self-interaction of Mo d states, we performed calculations with the Hubbard corrected DFT + U method [60][61][62][63][64]. A selection of input and output files are available as Supplementary Materials. Specifically we use two different formulations: the traditional DFT + U with only on-site interactions (U), and a generalized corrective functional based on the extended Hubbard model containing both on-site (U) and inter-site (V) interactions, known as DFT + U + V [64,65]. While mostly employed for crystals, these corrective functionals have also been used to study molecular systems containing TM magnetic centers, for example to assess the equilibrium geometry and dissociation energies [66][67][68]; the relative stability of various magnetic configurations and their implications on chemistry [69]; the potential energy surface along relevant reaction pathways [70,71]; the energetics of low-lying excited states using a Slater-1/2 approach [72]. In many cases the results compare favorably with those of more sophisticated quantum chemical approaches. Here the Hubbard U is applied to the Mo d states, while the inter-site V operates between the Mo d and O p states and between Mo d states. In general, the extended DFT + U + V scheme improves the description of electronic localization, especially in presence of significant hybridization between neighbor sites. Therefore, this method should be quite suitable for the system considered here. Quite surprisingly, we found that the two Hubbard corrections give almost equivalent results and will be collectively indicated as "DFT+Hubbard".
DFT+Hubbard calculations were conducted following a similar scheme as for the GGA-PBE functional. Table 2 reports the Hubbard parameters obtained for the considered states of the system; these quantities were calculated using a self-consistent procedure which involves the re-optimization of both the electronic structure and molecular geometry (see the Materials and Methods section). The Table only shows results for the molecule in gas phase. For the DFT+U + V calculation on the calcite surface we used the same value obtained in the gas phase. The Hubbard parameters reflect the substantial structural equivalence of the two Mo atoms, except for those magnetic configurations (AFiM2 and FiM4) where they assume quite different magnetic moments. In any case, the dependence of the Hubbard couplings on the magnetic configuration turns out to be relatively small (a negative value is obtained for the FM6 state) with a more pronounced decrease only for the FM8 state. Note also that the V between the two Mo centers is always much smaller than the one between each Mo and the O atoms in its ligands shell. This result is what makes it easy for d electrons on each Mo to align their spin and give rise to spin-polarized ground states (see the discussion below Table 3).
The main results of Hubbard corrected calculations are summarized in Table 3, which illustrates the relative stability of all the meta-stable magnetic states obtained with DFT + U and DFT + U + V for MoMo in the gas phase and with the molecule adsorbed on calcite (denoted DFT + U + V ads ). The scenario that emerges from DFT+Hubbard calculations is radically different from that of GGA-PBE. First of all, the nonmagnetic solution is not the GS state of the system (in fact, it is one of the highest in energy). While still characterized by a zero total magnetization, the GS shows instead an AFM order, with the magnetic moments of the two Mo atoms pointing in opposite directions. This is a first indication that, due to a more pronounced localization of their d states, the Mo's retain much of their atomic identity. It is important to note that the AFM GS features two almost completely spin-polarized Mo atoms, essentially as if all eight electrons forming the quadruple bond in the NM state were now localized on either Mo atom, depending on their spin. This outcome might seem quite surprising in light of the fact that other AFM could in principle be obtained by breaking one of the four bonds at a time, and localizing the two electrons of opposite spin on either Mo atom. Indeed solutions of this type (intermediate between AFM and NM) could be stabilized using values of V Mo−Mo in the 1.5-3.5 eV range, much higher than the self-consistent values reported in Table 2. Note that the FM2 state is in fact of this kind, with three remaining Mo-Mo bonds and two electrons (forming the δ bond in the NM state) segregated on either Mo with parallel spin. This state features a Mo-Mo distance which is very similar to the experimental value [56]. The stabilization of a NM and of "intermediate" AFM states is prevented by the low value of the self-consistent V Mo−Mo which is not effective in establishing a strong hybridization between Mo centers. It is notable anyway that the molecule does not break apart thanks to the external C-CH 3 linker groups. Whether this picture is an artifact of our linear response-based self-consistent evaluation of the Hubbard parameters, remains to be established.
Compared to the GGA-PBE results shown in Table 1, the Mo-Mo distance is almost exactly the same for states of corresponding M, indicating that the geometry of the molecule mostly depends on the number of broken Mo-Mo bonds rather than on the degree of d states localization.
The Hubbard correction has instead a role in destabilizing odd-valued magnetization (OVM) states. In fact, no OVM can be found in Table 3. These solutions could not be stabilized with DFT+Hubbard by constraining M, not even at the lowest odd-value (M = 1 µ B ). A OVM state corresponds to fractional transfers of electrons between the two spin channels, i.e., to fractional occupations of some Mo d states. Fractional occupations are made energetically unfavorable by the Hubbard correction, which operates on these states and favours integer occupations, compatible with the even-valued magnetization (EVM) states shown in Table 3.
The energy comparison between the NM and AFM solutions indicates that, at variance with GGA-PBE, the total energy depends not only on the total magnetization (equal to zero in both AFM and NM) but also, and more significantly, on the relative spin orientation of the two Mo atoms. By increasing the total magnetization to M = 2 µ B and M = 4 µ B , the total energy increases. As evident from the atomic magnetization of each Mo (n ↑ -n ↓ ), for both M = 2 µ B and M = 4 µ B , we were able to achieve two kinds of states: a ferromagnetic one (FM2 and FM4), with both atomic moments essentially equal; and a ferrimagnetic state (AFiM2 and FiM4) with very different moments. In fact, while for both M values the system was initialized with the two spin in opposite directions, only for M = 2 µ B an "antiferrimagnetic" (AFiM2) ground state was stabilized, while for M = 4 µ B the system preferred a ferrimagnetic (FiM4) configuration. The AFiM2 and FiM4 states are only possible through a reduction of molecular symmetry, that makes the two Mo atoms loose their equivalence. Note that, independent of the value of M, the AFiM2 and FiM4 solutions have energies lower than the FM ones (FM2 and FM4) with the same M. This stems from the fact that when the equivalence of the Mo's is lost, one of them can reach a saturation magnetization (3.36 µ B ), approaching the Hund's rule limit (highest S z ).
AFiM2 is one of the states whose energy is closest to the GS energy (E ≈ 1 eV). Importantly, a comparable energy can be obtained even for FM configurations of higher magnetizations, such as M = 6 and 8 µ B (FM6 and FM8). The total energy of FM6 (E ≈ 1.3 eV) is lower than FiM4, while FM8 presents an energy even slightly lower than AFiM2. Note that for M > 4 µ B only FM configurations can be stabilized, a fact related to the atomic nature of the magnetization, unlike in GGA-PBE. In fact, with ≈4 electrons in the d shell of Mo 2+ ions, an AFiM configuration with M > 4µ B would require maximum magnetization on one Mo and vanishing magnetization on the other or even more exotic configurations entailing, for example, the spin polarization of O 2p states. The electronic structure of the single Mo atoms is also probably responsible for the energy decrease of FM6 and FM8 configurations as compared to FM4. In fact, these magnetizations require both (equivalent) Mo atoms to assume atomic magnetizations of ≈3 and 4 µ B , respectively, which are within reach and allow to satisfy Hund's first rule.
Importantly, unlike those achieved with the GGA-PBE functional, all the DFT+Hubbard magnetic configurations are (meta-)stable, i.e., energy minima with vanishing atomic forces. This is related to the Mo atomic character and is a consequence of the finite gap appearing in the molecular spectrum, giving the system enhanced stability. In other words, the finite energy cost of transferring electrons between states of opposite spin prevents the system to change M (and the Mo atoms to change their moments), avoiding the evolution towards the AFM ground state (E = 0 eV).
The DFT + U and DFT + U + V results were presented together, as the two calculational schemes provide qualitatively similar results, including atomic magnetizations and Mo-Mo distances. However, there are some differences that are worth pointing out. First, one can observe a slightly different energy ordering. While AFM and FM8 are in both schemes the GS and the first excited state (maintaining, remarkably, the exact same energy difference), the next higher state for DFT + U is the ferromagnetic FM6 configuration, while it is the AFiM2 for DFT + U + V. A second difference concerns the FM4 state. In spite of its initialization with equal Mo atomic moments, within DFT + U the FM4 is unstable towards the saturation of one of the Mo magnetization and an antiferrimagnetic solution practically identical to FiM4 (except for the inversion of the two Mo atoms), while in DFT + U + V the two moments are very similar, with an energy significantly higher than FiM4. The imbalance of about two electrons in the population of the two spin for the FM4 state with DFT + U + V suggests that there are two residual bonds between Mo, as also confirmed by the shorter Mo-Mo distance ( 2.4 Å) compared to the one obtained with DFT + U (2.7 Å). This is analogous to what has been already discussed for the FM2 state, featuring three residual Mo-Mo bonds. At the same time FM6 corresponds to a configuration with only one bond is surviving.

Magnetic States of MoMo Adsorbed on the Calcite (10.4) Surface
We present now calculations for the MoMo molecule adsorbed on the calcite (10.4) surface (see Figure 1b). The results, shown in Table 3 (DFT + U + V ads ), can be directly compared with the ones for the isolated molecule. Given the significant increase in the computational cost, on-surface calculations were performed exclusively with DFT + U + V and only for a subset of configurations, namely the three lowest-energy configurations, AFM, AFiM2 and FM8, with the addition of the NM one to ascertain that the MoMo magnatization is not eroded in presence of the surface. For each magnetic configuration, a structural relaxation of MoMo on calcite was performed starting from the equilibrium configuration obtained in the gas-phase. As evident from Table 3, besides the NM configuration being pushed to higher energy than in the gas-phase, the main effect of the adsorption is the inversion of the energy order of the two lowest excited states, AFiM2 and FM8. While in the gas-phase these are almost degenerate, for the adsorbed molecule AFiM2 results 0.35 eV lower than FM8. At the same time, for all the considered states, the Mo-Mo distance slightly decreases with respect to the gas phase, with a concomitant reduction of the Mo atomic magnetic moments for both the AFM and the AFiM2 states.
In summary, we conclude that the adsorption on the calcite surface favors, for the MoMo molecule, states with lower magnetization. We point out, however, that the adsorption on a nonmagnetic insulating surface like calcite (10.4) does not alter significantly the gas-phase magnetic landscape of MoMo. In particular, as the AFM is still the GS of the system, the presence of the surface does not suppress magnetism, as often observed with metallic surfaces [11,43]. While a detailed chemical analysis of these results is left for future investigations, this represents an indication of the importance to use bulk insulating substrates to preserve the magnetic properties of potential SMMs.

Magnetic Anisotropy Energy
Considering one of the magnetic metastable state (we chose the lowest energy FM8) we now evaluate the magnetic anisotropy energy (MAE), i.e., the energy needed to change the orientation of the total magnetization with respect to the structure of the molecule, from its ground state position (defining the easy magnetization axis, EMA) up to an orthogonal direction (see Figure 2). The energy cost is associated with a finite spin-orbit coupling (SOC), which binds the magnetic moment to the geometrical structure of the molecule. The SOC was included in the calculations by updating the Mo pseudopotential to a fullyrelativistic version. We define the orientation of M via the polar angle θ, with θ = 0 in correspondence of the Mo-Mo direction, Figure 2a. Given the substantial σ h symmetry of MoMo (only broken by the extremal H atoms) a symmetric behavior for angles θ > 90 • is assumed (magnetic bistability).
Our SOC calculations used a non-collinear magnetic formalism and have been performed with DFT + U for the isolated module. We assumed the lowest-energy FM configuration, FM8, as the reference state. The MAE is then obtained from the dependence of the total energy on the polar angle θ between the total magnetization and the EMA, which, we anticipate, our calculations demonstrated to coincide with the Mo-Mo axis. Examples of input and output files of MAE calculations are available as Supplementary Materials. The results are shown in Figure 3, while the technical details of these calculations are discussed in the Material and Methods Section. The plot compares two sets of results, that were obtained for two different values of the azimuthal angle φ of the magnetization around the EAM. In fact, given the lack of cylindrical symmetry of the molecule around the Mo-Mo axis, with four acetate arms protruding along perpendicular directions in the x-y plane (with the EAM taken as z), we assumed a (light) dependence of the energy on this variable as well. The results shown in the figure, obtained for two values of φ corresponding to one of the acetate arms (black diamonds) and with a bisecting direction (solid blue line) show a very minor effect of φ on the total energy. The energy cost associated to the deviation of the magnetization vector from its minimum energy direction is thus almost entirely due to the polar angle θ it forms with the EMA. The figure shows that the energy barrier is about 9.3 meV. This value is compatible with the one found in the [Mn 12 ] compound (5.2 meV) [1] but lower than the MAE observed in single lanthanide ion SSMs [7]. However, to the best of our nowledge, no experimental data exists for the considered system.
It is important to remark that the MAE of the system is entirely due to (and so, dependent on) the chemical environment of the Mo atoms in the molecule. In fact repeating the same spin-orbit calculations on the bare Mo dimer resulted in a flat energy profile.

Materials and Methods
DFT calculations were performed using the publicly available, plane-wavepseudopotential QUANTUM ESPRESSO package [73,74]. The xc functional was constructed using a generalized gradient approximation (GGA) with the PBE parametrization [59]. In molecule on-surface calculations, van der Waals interactions were included through the Grimme D2 scheme [75]. Ultrasoft pseudopotentials, taken from the PSlibrary [76] were employed for all species. The electronic wavefunctions and charge density were expanded up to a kinetic energy cut off of 80 and 640 Ry, respectively.
To make the energy comparisons consistent, all calculations were performed using the simulation cell resulting from the structural optimization of the calcite (10.4) surface. The cell has an almost perfect orthorhombic symmetry with in-plane lattice parameters of 15.13 and 16.21 Å (the angle between them being 89.996 • ) and a perpendicular axis of 26.0 Å. The Brillouin zone was sampled with the Γ point.
The calcite (10.4) substrate was modeled with a periodically repeated slab of three layers, allowing a vacuum gap between the adsorbed molecule and the slab replica of 18 Å. The force threshold for atomic relaxations was set to 10 −3 eV/Å. Only the molecule and the slab top-layer were allowed to relax.
To capture the Mo d electron localization, both DFT + U [60][61][62][63][64] and the extended DFT + U + V [64,65] Hubbard corrective functionals were adopted. These schemes have been successfully used on a broad variety of cases and proved very effective in improving the structural [77], vibrational [64,[78][79][80][81][82][83][84], and electrochemical [85] properties as well as the phase-stability and transitions of d and f electron systems, notably involving magnetism [86][87][88][89]. Like in many of the works cited above the effective Hubbard interactions (i.e., the on-site U for Mo d states and inter-site V between Mo d and nearest-neighbor O p states) were determined from first-principles linear-response theory [77] using a density functional perturbation theory (DFPT) implementation [90,91]. The calculation of U and V was done for the molecule in gas phase in a self-consistent way that ensures the consistency of their values with both the electronic structure and the geometry of the system [64,85,89,91]. The obtained U and V values were also adopted for the molecule adsorbed on calcite, assuming negligible changes. Hubbard corrections are based on orthogonalized atomic orbitals and the same basis is adopted for the forces calculation (thanks to a recent extension [92]), thus improving the consistency between the molecular geometry and the electronic ground state.
Most calculations were performed with fixed KS occupations, possibly with an imbalance between the two spin populations to achieve the desired total magnetization.
As mentioned in Section 2.3 the calculation with the SOC required updating the pseudopotential of Mo to a fully relativistic version. In order to ascertain that this pseudopotential provides an equivalent description as the non relativistic one a preliminary set of calculations with the SOC was performed with DFT + U for the collinear (magnetic moments all along the Mo-Mo direction) AFM and FM8 configurations. These calculations recovered the same results (in particular, the same energy differences) as those shown in Table 3 for the collinear case.
The calculation of the MAE requires changing the polar angle θ between the total magnetization M and the easy axis of magnetization (EMA). Unfortunately, from a technical standpoint, acting directly on the magnetization vector (e.g., through the coupling with an auxiliary magnetic field pointing to the desired direction) is not straightforward and the small entity of the SOC with respect to other energy terms makes these calculations even more difficult to perform and delicate to converge.
In order to circumvent these problems we designed an ad hoc algorithm that uses an auxiliary system and the Hubbard correction as an effective constraint functional. This novel procedure is illustrated by the following steps:

1.
For each angle θ, perform a non-collinear calculation on an isolated Mo-Mo dimer (the auxiliary system) without SOC, imposing (from input) the desired direction to M. Use the same Hubbard U and Mo-Mo distance as in the molecule.

2.
Save the d states occupations of the two Mo atoms.

3.
Perform a preliminary DFT + U + SOC calculation on the molecule, starting from the saved d states occupations, which are kept fixed for ≈10 electronic iterations.

4.
From the potential achieved at the previous step, proceed with the DFT + U + SOC calculation, this time letting the occupations free to evolve and the molecule reach its self-consistent GS, at the given θ.
The calculation on the isolated dimer, point 1, is an expedient to achieve, at any angle θ, the corresponding atomic occupations. In fact, in absence of SOC, the energy of the system does not depend on θ and its value can be tuned at will very effectively. In addition, calculations on the bare dimer are much less intensive than in the molecule and this first step can be completed very efficiently. At this point the atomic occupations of Hubbard atoms (the two Mo here responsible for the magnetization) are compatible with a magnetic moment pointing in the desired direction. The calculation in point 3 is necessary to train the self-consistent potential until it becomes consistent with the desired deviation θ of the magnetization from the original direction. This "training" is enforced by the Hubbard potential (a finite U is of course needed) constructed with the atomic occupations obtained in point 1. Once the potential is "trained", the system is set free to converge to the self-consistent ground state. The obtained magnetization is typically within 1 degree from the desired direction. Of course the above procedure has to be repeated for all the directions of magnetization one desires to explore.

Discussion and Conclusions
We perform a detailed computational study of the magnetic states of the dimolybdenum tetraacetate molecule (MoMo), both in gas-phase and adsorbed on the calcite (10.4) insulating surface. The long-term goal of our investigation is to assess the possibility to use MoMo as a single molecule magnet (SMM), which would open up the possibility to have magnetic molecular networks on a bulk insulating substrate. The work is based on the recent observations that regular MoMo assemblies hold firmly on calcite (10.4) at ambient conditions [48]. MoMo seems particularly interesting for the vertical alignment of the Mo-Mo axis, that favors the symmetry breaking between the two Mo on the surface (a favorable condition for the onset of magnetic states [93]) and the external C-CH 3 linkers that prevent the molecule dissociation, even when all four Mo-Mo bonds are broken.While standard GGA-PBE functionals predict a nonmagnetic (NM) state with quadruple Mo-Mo bond and excited magnetic states at significantly higher energies (growing with the number of broken bonds and ultimately with the magnetization M), Hubbard corrections stabilize magnetic configurations (FM or AFM) with Mo atoms near their high-spin configurations and atomic magnetizations near their Hund's saturation limit of four. More importantly, the corrections render many FM states metastable. Although the obtained AFM ground state does not seem consistent with available literature on similar systems (but notably missing external linkers between the Mo centers), our calculations show that MoMo might have magnetic states provided that the system is excited to metastable configurations with Mo atoms pulled apart from each other.
In the gas-phase, the FM state at the lowest energy (denoted FM8) is about 1 eV above the ground state. It features parallel spins on the Mo atoms (close to maximum atomic magnetization), with a total magnetization of M = 8 µ B . On the surface, the MoMo second lowest-energy state features an antiferrimagnetic (AFiM) order with the atomic moments of the two Mo atoms having opposite direction and different length, still resulting in M = 2 µ B . Importantly, we find that the presence of the insulating surface does not quench the magnetic ordering found in the gas-phase. This is a necessary condition for having stable and regular SMMs arrays on an insulating substrate at room conditions. From this point of view, insulators seem an advantage over metallic surfaces, where the substrate-molecule electronic hybridization could alter or completely quench the magnetic properties of the isolated molecule.
For the lowest energy ferromagnetic state (FM8) we calculated the magnetic anisotropy energy (MAE) in the gas-phase, using a noncollinear formalism. We found that the molecule possesses an easy magnetization axis (EMA) that coincides with the Mo-Mo direction and we estimated a MAE barrier around 9.3 meV. This value is compatible with the one found in the [Mn 12 ] compound (5.2 meV) [1] although, to the best of our knowledge, no experimental magnetic measurements of MAE exist for our system.
The system could also represent a prospect SMM worth of further investigations, if a state with a non-vanishing total M could be experimentally stabilized on the calcite surface. In this perspective, a first direction of investigation concerns the calculation of the MAE for the molecule adsorbed on the surface, to ascertain that a finite magnetic moment is not too volatile. Other aspects of interest include the assessment of other insulating surfaces (e.g., fluorite) and the characterization of the behavior of molecular systems on metallic surfaces, aimed at verifying the survival of finite magnetic moments [49]. It is also worth extending the present investigation to other TM-based molecules that can be adsorbed firmly on insulating surfaces. These include systems with non-planar TM substructures (like MoMo, with the two Mo atoms aligned along a vertical direction) in which the symmetry breaking upon adsorption may favor the stabilization of magnetic ground states and other interesting properties.
From a methodological point of view, a relevant novelty introduced in this work is the use of the Hubbard corrective functional as a constraint on the magnetization M, to rotate its direction from the system EMA. The algorithm developed here for this purpose has proved both reliable and effective and will represent a valuable tool in future calculations of this type. Acknowledgments: M.C. is grateful to A. Lascialfari for useful discussions. We would like to thank A. Kühnle for discussions and for carefully reading the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: