We have used different computational approaches to carry out the different calculations for the three topics developed in these studies. In the lines of research regarding the interaction of small molecules based on phen derivatives and Mo complexes, including phen with dDNA and GQ, we used several models along with different levels of calculation. That is, in the case of sandwich models [49], we used the M06-2X/6-31+G(d,p) level of calculation [50,51,52] with Gaussian09 in order to study the intercalation of phen derivatives with dDNA. Ref. [53] The M06-2X functional is recommended to study systems in which weak interactions such as dispersion are important to explain their behavior (as the intercalation of small molecules between bps of DNA). In the case of the ring models [49] including sugar, phosphates, and Na⁺ counterions, and for the study of the groove binding interactions with the d(GTCGAC)₂ hexamer, we used the semi-empirical Hamiltonian PM6-DH2 [54] including dispersion effects with the MOPAC software. Ref. [55] In the case of the ring model 2 ps semi-empirical MD, simulations with the PM6-DH2 Hamiltonian were also performed with a time step of 1fs. Explicit water molecules were considered by means of the TIP3PBOX solvent model for the simulations by using a rectangular box with edges no closer than 5 Å to any atom of the solute. Finally, in the case of the studies on the interaction of the Mo[(η³-C₃H₅)Br(CO)₂(phen)] complex with dDNA and GQ, we took into account the d(AGACGTCT)₂ octamer for the dDNA, coming from the 1n37 PDB structure, and the GQ coming from the 2jwq PDB structure to study the interaction of the Mo[(η³-C₃H₅)Br(CO)₂(phen)] with this non-canonical structure. These structures were studied at the LMKLL/DZDP level, [56,57] which includes van der Waals corrections. The core electrons were substituted by norm-conserving pseudopotentials. Refs. [58,59] Such LS-DFT computations with ~500 atoms for the dDNA and ~1000 atoms for the GQ were carried out with the Spanish Initiative for Electronic Simulations of Thousands of Atoms (SIESTA) method and associated software. Ref. [60] The cut-off radii for the atomic orbitals of each element were obtained for an energy shift [60] of 30 meV. The tolerances used for the optimizations were 10⁻⁵ eV for the energetics, whereas the tolerance for the forces was 0.02 eV/Å for the dDNA and 0.1 eV/Å for the GQ. In order to gain deeper insight into the interaction between the small molecules and DNA substrates, we performed the Energy Decomposition Analysis (EDA). Refs. [61,62] To carry out the EDA, we mainly used the B3LYP-D3/TZP level of calculation [63,64,65,66] since the B3LYP-D3 functional includes an explicit Grimme’s D3 correction for dispersion and, therefore, an additional ΔE_disp term appears in the EDA for this functional. For this reason, we thought that the discussion was simpler to visualize better the trends of the intercalated systems. It must be said that the M06-2X/TZP and M06-L/TZP levels were also checked for the EDA. In any case, the three functionals led to results comparable to the MP2/6-31G* (0.25) level of theory already used by Řeha et al. [67] which gave results comparable to the benchmark CCSD(T) data. Ref. [68] These EDA were carried out with the ADF software. Refs. [69,70,71] Another way to gain deeper insight into the interaction between the studied small molecules and DNA structures is the topological analysis of the structures. Two kinds of approaches were used, namely the classical QTAIM developed by Bader et al. [72] and the most recent approach developed by Johnson et al. [73] based on the NCI. The latter provides a rich, 3D representation for the non-covalent interactions with surfaces based on the peaks that appear in the reduced density gradient at low values of ρ. Such isosurfaces are mapped according to values of the sign of the second Hessian eigenvalue, (λ₂), and while negative values (i.e., stabilizing interactions) are depicted in blue and pale green, positive values (i.e., destabilizing interactions) are represented in yellow and red. QTAIM and NCI computations were performed with the AIM2000 [74] and AIMALL [75] software with the wave functions generated at M06-2X/6-31+G(d,p) level with Gaussian09 in all cases with the exception of the interactions of [Pt(en)(phen)]²⁺ derivatives via groove binding with dDNA in which the wave function was obtained at B3LYP/6-31G(d,p) level with Gaussian09. Finally, for some of these systems in which small molecules are interacting with DNA, we studied the polarization/charge transfer for the different modes of interaction with several charge schemes arising from different approaches (Mulliken, [76] APT, [77] Hirshfeld [78] and NPA/NBO [79,80]).

In the study of the activity of Mo-oxo species as promoters and catalysts for the hydrolysis of the phosphoester bond, we carried out the DFT calculations by using ADF [69,70,71] and the Gaussian 16 Revision A.03 software. Ref. [81] In the case of ADF, geometry optimizations were carried out with the BP86 functional, which uses the Vosko-Wilk-Nusair exchange-correlation potential [82] with the generalised gradient approximation exchange correction reported by the Becke (1998) exchange functional [83] and Perdew (1986) correlation correction [84], as well as Grimme dispersion corrections (BP86-D3). Ref. [66] Relativistic effects were treated with the zero order regular approximation (ZORA) Hamiltonian. Refs. [85,86] The frozen core approximation and triple-ζ Slater-type orbitals (STO) were used to describe the valence shells of C and N (2s and 2p). One polarization function was added to C, N, O, and Mo (single-ζ, 3d, 4f). Triple-ζ STOs were used to describe the valence shells of H (1s) augmented with one polarization function (single-ζ, 2s, 2p). Solvent effects were included with the COSMO [87] and standard parameters (Water, ε = 78.39). Analytical frequencies were calculated to characterise the obtained stationary points and calculate the Gibbs free energies (standard state T = 298.15 K, p = 1 atm). Transition states were followed after a fractional displacement of the imaginary vibrational mode to both the reactant(s) and product(s). On the other hand, Gaussian calculations were carried out with the B3LYP hybrid functional. Refs. [63,64,65] The LANL2DZ effective core potential with the associated double zeta basis set supplemented with f polarisation functions was used for Mo atoms [88,89,90,91,92] and the 6-31+G(d,p) basis set [52] for the rest of the atoms. Dispersion effects were included by using the third version of Grimme dispersion with the Becke–Johnson damping approach. Ref. [93] Stationary points were characterised by means of frequency calculations, and the intrinsic reaction coordinate was followed in order to obtain the geometries of the reactants and products followed by unconstrained optimisations. Refs. [94,95] For some mechanistically relevant stationary points the energies were refined with the 6-311++G(3df,2p) basis set. Ref. [96] For these calculations with Gaussian 16 Revision A.03 the Polarisable Continuum Model (PCM) was used to take into account the solvent in an implicit way [97,98].