302
Acta Chim. Slov. 2014, 61, 302-307
Scientific paper
Rigid- and Polarizable-ion Potentials for Modeling Ru-polyoxometalate Catalysts for Water Oxidation
Changru Ma,1 Simone Piccinin1 and Stefano Fabris1
1 CNR-IOM Istituto Officina dei Materiali, Centro DEMOCRITOS and SISSA, Via Bonomea 265, Trieste, Italy * Corresponding author: E-mail: fabris@democritos.it
Received: 23-12-2013
Paper based on a presentation at the 4th RSE-SEE 2013 Symposium on Electrochemistry in Ljubljana, Slovenia
Abstract
This work assesses the predictive power and capabilities of classical interatomic potentials for describing the atomistic structure of a fully inorganic water-oxidation catalyst in the gas phase and in solution. We address a Ru-polyoxometala-te molecule (Ru-POM) that is presently one of the most promising catalysts for water oxidation due to its efficiency and stability under reaction conditions. The Ru-POM molecule is modeled with two interatomic potentials, the rigid ion model and the shell model potentials, which are used to perform molecular dynamics simulations. The predictions of these two approaches are discussed and compared to the available ab-initio data. These results allow us to establish the suitable level of theory to model complex heterogeneous interfaces between the Ru-POM and electrodes in solution.
Keywords: Artificial photosynthesis, homogeneous catalysis, polarizable interatomic potentials, molecular dynamics simulations
1. Introduction
The sunlight-driven splitting of water into molecular H2 and O2 allows for storing solar energy into chemical fuels.1-4 The water oxidation half-reaction (or oxygen evolution reaction) is by far the most challenging step in the electrochemical water splitting and represents a bottleneck for the development of efficient artificial photosynthesis devices for the production of solar fuels.5-7
The most efficient oxygen-evolving (OE) catalysts are metallorganic molecules or metal-oxide surfaces based on Ru, Ir, Fe, Co or Mn. These two classes of homogeneous (metallorganic molecules) and heterogeneous (metal oxides surfaces) catalysts presents different advantages and disadvantages in terms of the chemical reactivity, catalytic efficiency, and structural stability during opera-tion.4 The development of alternative catalysts that could combine the best properties of these catalysts classes has drawn much attention during the past years.
A fully inorganic Ru-based molecular complex has recently emerged as a very promising candidate for the anodic reaction in artificial photosynthetic devices.8'9 The Ru-polyoxometalate complex (Ru-POM) is displayed in Fig. 1. It is a 10- ion and consists of a tetraruthenium-oxo core [Ru4O4(OH)2(H2O)4]6+ sandwiched between two [SiW10O36]8- POM units. Ru-POM has been reported to
promote water oxidation with low over-potential (0.35 V), high turn over frequency (> 450 cycles per hour) and no deactivation. It is considered to be one of the best OE catalysts for water oxidation reported to date.
Figure 1. Molecular structure for the Ru-POM complex. Green, red, white, yellow and pink spheres represent Ru, O, H, Si and W atoms, respectively.
X-ray diffraction measurements show that the active core of the Ru-POM catalyst is formed by four Ru atoms at the vertexes of a tetrahedron and that they are linked by two ^-hydroxo and four ^-oxo bridges.8 9 The Ru-POM molecule has D2d symmetry. One water ligand coordinates each Ru center, whose oxidation state at zero applied po-
Ma et al.: Rigid- and Polarizable-ion Potentials for Modeling
302	Acta Chim. Slov. 2014, 61, 303-307

tential is proposed to be Ru(IV). Together with the two POM caps, the Ru-POM anion has a 10- charge.
In a previous work we have investigated the structural, electronic and thermodynamic properties of this molecule in the gas phase and in solution by means of Density Functional Theory (DFT).10 These calculations suggested that the frontier orbitals of the catalyst are localized on the tetraruthenium-oxo core. Assuming four proton-coupled electron-transfer (PCET) oxidation steps, this study demonstrated that the free energy of the Ru(IV)-H2O / Ru(V)-OH couple (3.38 eV with B3LYP and 4.00 eV with HSE06) is significantly lower than the thermodynamical limit for water oxidation (4.92 eV in experiment and 4.56 eV in B3LYP calculations). This demonstrates that higher oxidation states of the Ru-POM would be involved to promote the oxidation of water. More recent calculations suggest that the catalyst is activated by the formation of a Ru-oxo moiety.11 Ab initio metadynamics simulations predict that the formation of O-O bond is triggered via the nuc-leophilic attack of a solvent water molecule.
To anchor the catalyst to a conductive substrate and hence to fabricate an anodic material based on this molecular catalyst, Toma et al.12 as well as Guo coworkers13 designed a specific support that binds Ru-POM and that preserves the catalytic efficiency of this molecular catalyst. The support consists of multiwalled carbon nanotu-be functionalized with polyamidoamine ammonium den-drimers (MWCNT-dend). A schematic model of this complex catalyst is displayed in Fig. 2. Resonant Raman spec-troscopy indicates that the structure of Ru-POM is preserved during the assembly process. STEM images clearly demonstrate that the catalyst binds to the nanotube surface. Small angle X-ray scattering diffraction spectra show that the inorganic catalyst approaches single-molecule homogeneous behavior. On the basis of cyclic voltammetry it is concluded that the Ru-POM/MWCNT-dend assembly preserves the electro-catalytic water oxidation of the unsupported Ru-POM molecule and that this nano structured electrode is as stable as the homogeneous complex, suggesting that the catalytic mechanism of water oxidation is preserved even in the presence of MWCNT-dend / Ru-POM interaction. In conclusion, the conductive MWCNT-dend provides heterogeneous support for the Ru-POM cluster, allows controlling the material morphology, and facilitates the sequential electron transfer to the electrode. The Ru-POM/MWCNT-dend system provides a promising electrode for the water oxidation in artificial photosynthesis device. A detailed understanding of the structural, electronic, and chemical properties of the Ru-POM/MWCNT-dend interface is still missing. To this end, numerical materials modeling can provide useful insight, however the size and complexity of the system (Ru-POM, MWCNT, and polyamidoamine ammonium dendri-mers, see Fig. 2) together with the presence of the solvent water molecules call for combining different levels of theory and computational techniques. This requires develo-
Figure 2. Molecular structure for the full Ru4-POM complex @ functionalized graphene from the top view a) and the side view b). Red, pink, white, yellow and green spheres represent O, W, H, Si and Ru atoms, respectively. Red, white, blue, cyan sticks represent O, H, N, C atoms, respectively. Part of graphene layer and all water molecules are not displayed for simplicity.
ping empirical interatomic potentials for the Ru-POM/MWCNT-dend interface that would allow performing molecular dynamics on the complete electrode, also including the presence of the solvent. Such simulations have the potential to provide key information on how the Ru-POM catalyst binds to the functionalized MWCNT electrode through the organic dendrimers. As a first step towards this goal, in this work we explore the predictive power of classical interatomic potentials for describing the structure of the Ru-POM/MWCNT-dend interface.
2. Methods
In this work we consider two levels of approximation for describing the Ru-POM system in solution, namely the rigid ion model (RIM) and the shell model (SM) potentials. Both approaches are empirical and rely on a set of parameters, which we partially fit to our ab initio calcu-lations10 and partially adapt from previous works. The parameter sets will be described in the following section, here we introduce the functional form of these interatomic potentials and define the parameters involved.
The RIM was used to describe the interactions between the solvent water molecules and the Ru-POM molecule, which was treated as a rigid body. The functional form of the RIM consists of a Lennard-Jones (LJ) potential and of a Coulomb charge-charge electrostatic term:
\ t "ij \ <iW (1)
As a refinement to the rigid ion model, we also consider the SM potential that includes higher multipoles in the electrostatic expansion up to dipole-dipole terms. In the SM,14 the partial polarization of the ionic bonding is included by describing the atoms with a charged core connected by a harmonic spring to a massless charged shell.
Ma et al.: Rigid- and Polarizable-ion Potentials for Modeling
302
Acta Chim. Slov. 2014, 61, 304-307
The core and the shell interact through a harmonic potential, whose strength is one of the parameters of the model:
where ki is the spring constant and 8ri is the distance between the shell center and the core of atom i.
The equilibrium position of the core and of the shell are obtained with a self-consistent procedure, which leads to an effective electronic polarization of the atom. In addition, the SM consists of the Coulomb interactions and of the short-range Buckingham two-body potentials.
where Atj, ptj and Cij are parameterized constants specific to each pair of atomic species i and j (i ^ j) at a distance r.
Water solution is modeled as the widely used TIP3P model (transferable intermolecular potential 3P).15 The TIP3P model specifies a 3-site water molecule with charges and Lennard-Jones parameters assigned to each of the 3 atoms. O-H bonds (r) and H-O-H angles (6) are controlled via harmonic potentials:
Since the outer surface of the POM caps comprises of oxygen atoms, we follow the approach used in Ref. 10 and limit the interaction of the POM with the solvent to these exposed O atoms only. In particular, their interaction with the TIP3P water molecules and Na+ counter ions was modeled with the same LJ parameters of the oxygens in the TIP3P model.
The calculations are performed with the General Utility Lattice Program (GULP)1617 and the DL_POLY Molecular Simulation Package.18
Table I. Shell model potential parameters
	qcre	qshell	k(eVÄ~2)
OPOM	0.848	-2.848	74.9204
OruO	0.513	-2.513	20.53
Si	4.0		
W	6.0		
H	1.0		
			
	AjfeV)		C..(eVÄ"6)
Si-Oshdl	1283.907	0.32052	10.66158
W_osheU	1305.22	0.375	0.0
H-Oshdl	227.046	0.27620	0.0
Ru-OsheU	13733.40	0.2259	0.49
°RuO"°RuO	25.41	0.6937	32.32
°POM-°POM	22764.00	0.149	27.879
3. Results and Discussion
3. 1. Shell Model Parameterization and Simulations
There are no SM parameters available in the literature for the specific case of the Ru- POM molecule. All the O atoms of the Ru-POM molecule are modeled as a core of charge qcore and a shell of charge qsheU, while all the other atoms (W, Si, Ru, H) are represented only by point charges. For the electrostatic and core-shell interactions, the O atoms in the POM caps and in the tetraruthenium-oxo core are described as different species, OPOM and ORu-O, respectively.
Our starting point for the parameterization of the SM for the Ru-POM system is the studies of Ru oxide,19 ZrW2O820 and SiO2.21 In particular, we take the parameters for the Ru-OsheU interaction from the work of Battle et al., The W-Oshell and OPOM parameters from Pryde et al.,20 while those for Si-OsheU from Sanders et al.21 The parameter set employed for the SM simulations of Ru-POM molecule is reported in Table I.
These parameters were used to perform a molecular dynamics simulation with the SM potential for the Ru-POM molecule in the gas phase at room temperature. The MD simulation was performed in the NVT ensemble with a time step of 1 fs. The equations of motion were integrated with the Verlet algorithm for 10 ns. We report in Fig. 3a) the radial distribution functions for the W-O, Si-O, Ru-O and O-O pairs (red lines). The SM results are benchmarked against ab initio density functional theory (DFT) MD simulations of the Ru-POM in vacuum performed as reported in Ref. 10. The radial distribution functions obtained from these ab-initio MD are also displayed in Fig.3 (green lines). For all the atom pairs, both the heights and positions of the first and second peaks in SM MD simulations are in good agreement with the DFT results. We therefore conclude that our parameterization of the SM potential provides a reliable description of the stability and structure of the Ru-POM molecule in the gas phase.
These sets of SM parameters were then used to study the solvation of the Ru-POM molecule in water. Although the present SM describes well the structural properties of the catalyst in the gas phase, the interaction with the solvent described explicitly with the TIP3P model led to structural instabilities. In particular, the SM description of Ru-POM solvation in water dismantled the octahedral coordination around the Ru centers in the active core of the catalyst and opened up the POM caps (see Fig. 3b). Even by reparameterizing the potential, we could not identify a set of SM parameters capable to predict a stable and physically meaningful catalyst/solvent interface.
Given the high stiffness of the Ru-POM frame, in the following we explore a different approach for modeling Ru-POM molecules in water solution.
Ma et al.: Rigid- and Polarizable-ion Potentials for Modeling
302	Acta Chim. Slov. 2014, 61, 305-307

a)
A	DFT
fi __	
A _	
	
Table II. Parameters for the RIM potential of the Ru-POM molecule
Distance [Angstrom]
b)
Figure 3. a) Comparison of the SM and DFT radial distribution functions (RDF) of the W-O (top), Ru-O (middle), Si-O (middle) and O-O (bottom) pairs for the Ru-POM molecule in the gas phase. b) Equilibrated structure of the solvated Ru-POM resulting from MD simulations with the SM potential.
3. 2. Rigid Ion Model Simulations
In this approach, the structure of the Ru-POM molecule is kept fixed and we fit the LJ parameters of the wa-ter/Ru-POM interaction so as to obtain a structure of the solvation shell that is compatible with the one resulting from QM/MM calculations.
The RIM parameters obtained for the Ru-POM molecule and for its interaction with the solvent water molecules are reported in Table II. The charges on each atom, qp are derived by fitting the total DFT charge density according to densities and derived atomic point charges (DDAP) method of Bloechl22 and then averaged over all the atoms with the same type. We note that fitting the value of these charges on the electrostatic potential is not possible in this case because the definition RESP (Restrained ElectroStatic Potential) charges are not well defined for transition metals (Ru and W in our case).23 In principle one should expect two different types of O atoms, one for the RuO core (ORuO) and one in the POM caps (OPOM). Although we initially assigned two different types of oxygen, we found that the atomic charges for these two types are quite similar: -0.914 for OPOM and -0.902 for ORuO. Thus we simplified the model by working with only one type of O atoms for what concerns electrostatics. It turns out that these values are in line with the charge of the O atom in the TIP3P potential of water, -0.834. In summary,
Atoms	E (eV)	O (A)	Atomic Charge
OPOM	0.006596	3.8000	-0.912
ORu	0.006596	3.1507	-0.912
H	0.0000	0.0000	+0.589
W	0.009588	2.3400	+2.437
Si	0.02094	3.5800	+1.821
Ru	0.001561	1.2040	+1.883
Runter	0.1691	1.9477	
the fitting procedure allows for identifying a set of point charges that describe the electrostatic potential generated by Ru-POM.
The Lennard-Jones parameters for the interaction of the atomic species on the POM surface and water were selected as follow. The parameters for the W-Owater interaction were selected by following the work of Lopez et al.,24 who addressed the case of POM molecules in solution. The parameters for the Si atoms were taken from the work of Tang et al.,25 who simulated silicon monocrystals in water. The LJ parameters for the interaction between the OPOM and the Owater were set equal to those of the Owater in the TIP3P water. The LJ parameters for ORu and Ru were
a)
b)
Figure 4. a) Snapshot from the MD simulation of the Ru-POM molecule solvated in water described with the RIM potential. The inset shows the binding of a water molecule at the Ru sites of the tetrarut-henium-oxo core. b) Comparison between the classical and QM/MM
radial distribution functions of the Ru-O and Opom-H pairs
Ma et al.: Rigid- and Polarizable-ion Potentials for Modeling
302
Acta Chim. Slov. 2014, 61, 302-307
obtained by fitting the geometrical properties in QM/MM simulations of Ru-POM in solution10 calculated by CP2K, where the catalyst was treated quantum mechanically and the solution was described by TIP3P model as reference.
We show in Fig. 4a) a snapshot of the equilibrated structure of the Ru-POM molecule in solution obtained from the MD simulation with the RIM model. This level of theory provides a reasonable description of the Ru-POM/solvent interface, as shown by the good agreement with QM/MM results. In particular, the inset of Fig. 4a) show the Ru-water moiety. During the first MD steps of the Ru-POM in solution, four solvent water molecules bind to the four Ru sites of the active tetraruthenium-oxo core. The resulting Ru-Owater bonds at the Ru sites are approximately equal, with an average distance of «2.1 A predicted by the classical potential, and of «2.2 A by the QM/MM simulations. Also the other Ru-Owater and OPOM-Hwater peaks in the RDF (Fig. 4b) display a good agreement between the classical force field and DFT simulations.
4. Conclusions
In conclusion, we have explored the predictive power of different classical interatomic potentials for describing the structure of the Ru-POM molecule in solution. The SM potential has been successfully used to model metal-oxide nanoparticles, also based on Ru, W and Si. Building from this literature, we have parameterized a SM potential for the Ru-POM molecule and have tested the resulting equilibrium structure in the gas phase. The results show good agreement with the available DFT calculations, thus suggesting that the SM potential allows for capturing the stability and the main features of the intramolecular interactions. We could however not identify a set of parameters capable to correctly describe the Ru-POM/water interaction. Ru-POM solvation at the SM level resulted in structural instabilities, which are not observed in our QM/MM simulations.
The failure of the SM potentials for the Ru-POM/solvent interaction calls for a simplified description that nevertheless could provide an acceptable prediction of this solid/liquid interface. To this end, we constrained the geometry of the Ru-POM molecule, described its potential with a RIM, and captured its short-range interaction with the TIP3P solvent via a LJ potential. This approach describes the catalyst/solvent structure in reasonable agreement with the available QM/MM simulations, and constitutes the starting point for studying complex heterogeneous catalysts in solution such as the one displayed in Fig.2. In a forthcoming work,26 we employ this classical interatomic potential to study the possible binding modes of the solvated Ru-POM catalyst to the MW-CNT-dend support, thus addressing the complex structure of this functional interface for artificial photosynthesis.
5. Acknowledgements
We thank A. Sartorel, M. Bonchio, M. Prato, and G. Scoles for useful discussions. This work was partially funded by the EU through the FP7 Marie Curie IRG Program (Project H2OSPLIT, Grant PIRG04-GA-2008-239199). Computational resources were made available through PRACE Project ENCORASOL (2012061002) and Artifi-cialLeaf (2012061080). SF acknowledges the financial support from the EU-FP7 COST action CM1104.
6. References
1.	H. B. Gray, Nat. Chem. 2009, 1, 7-7.
2.	N. S. Lewis, D. G. Nocera, Proc. Natl. Acad. Sci. U.S.A. 2005, 103, 15729-15735.
3.	V. Balzani, A. Credi, M. Venturi, ChemSusChem. 2008, 1, 26-58.
4.	M. G. Walter, E. L. Warren, J. R. McKone, S. W. Boettcher, Q. Mi, E. A. Santori, N. S. Lewis, Chem. Rev. 2010, 110, 6446-6473.
5.	H. Yamazaki, A. Shouji, M. Kajita, M. Yagi, Coord. Chem. Rev. 2010, 254, 2483-2491.
6.	X. Sala, I. Romero, M. Rodriguez, L. Escriche, A. Llobet, Angew. Chem. Int. Ed. 2009, 48, 2842-2852.
7.	N. Armaroli, V. Balzani, Angew. Chem. Int. Ed. 2006, 46, 52-66.
8.	A. Sartorel, M. Carraro, G. Scorrano, R. De Zorzi, S. Gere-mia, N. D. McDaniel, S. Bernhard, M. Bonchio, J. Am. Chem. Soc. 2008, 130, 5006-5007.
9.	Y. V. Geletii, B. Botar, P. Kogerler, D. A. Hillesheim, D. G. Musaev, C. L. Hill, Angew. Chem. Int. Ed. 2008, 47, 38963899.
10.	S. Piccinin, S. Fabris, Phys. Chem. Chem. Phys. 2011, 13, 7666-7674.
11.	S. Piccinin, A. Sartorel, G. Aquilanti, A. Goldoni, M. Bonc-hio, S. Fabris, Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 4917-4922.
12.	F. Toma, A. Sartorel, M. Iurlo, M. Carraro, P. Parisse, Nat. Chem. 2010, 2, 826-831.
13.	S. X. Guo, Y. Liu, C. Y. Lee, A. M. Bond, J. Zhang, Y. V. Geletii, C. L. Hill, Energy Environ. Sci. 2013, 6, 2654-2663.
14.	B. Dick, A. Overhauser, Phys. Rev. 1958, 112, 90-103.
15.	W. Jorgensen, J. Chandrasekhar, J. Madura, R. Impey, M. Klein, J. Chem. Phys. 1983, 79, 926.
16.	J. Gale, J. Chem. Soc., Faraday Trans. 1997, 93, 629-637.
17.	J. Gale, A. Rohl, Mol. Simulat. 2003, 29, 291-341.
18.	I. T. Todorov, W. Smith, K. Trachenko, M. T. Dove, J. Mater. Chem. 2006, 16, 1911-1918.
19.	P. Battle, T. Bush, C. Catlow, J. Am. Chem. Soc. 1995, 117, 6292-6296.
20.	A. Pryde, K. Hammonds, M. Dove, V. Heine, J. Gale, M. Warren, J. Phys.: Condens. Matter. 1996, 8, 10973-10982.
21.	M. Sanders, M. Leslie, C. Catlow, J. Chem. Soc.: Chem. Commun. 1984, 19, 1271-1273.
Ma et al.: Rigid- and Polarizable-ion Potentials for Modeling
302	Acta Chim. Slov. 2014, 61, 95-307

22.	P. E. Blochl, J. Chem. Phys. 1995, 103, 7422.
23.	J. Wang, W. Wang, P. A. Kollman, D. A. Case, J. Mol. Graphics Modell. 2006, 25, 247-260.
24.	X. Lopez, C. Nieto-Draghi, C. Bo, J. B. Avalos, J. M. Pöblet, J. Phys. Chem. A 2005, 109, 1216-1222.
25.	C. Y. Tang, L. C. Zhang, Nanotechnology 2004, 16, 15-20.
26.	C. Ma, S. Piccinin, S. Fabris, Phys. Chem. Chem. Phys. 2014, 16, 5333-5341.
Povzetek
V tem prispevku smo ocenili napovedno moč in zmožnosti klasičnih medatomskih potencialov za opis atomistične strukture anorganskega katalizatorja oksidacije vode v plinski fazi in v raztopini. Obravnavali smo molekulo rutenijeve-ga polioksimetalata (Ru-POM), ki je trenutno zaradi svoje učinkovitosti in stabilnosti pri reakcijskih pogojih eden izmed najbolj obetavnih katalizatorjev za oksidacijo vode. Molekulo Ru-POM smo modelirali z dvema medatomskima potencialoma, rigidnim ionskim modelom in modelom lupine. Ta dva potenciala smouporablili pri simulacijah molekulske dinamike. Obravnavali smo napovedno moč teh dveh pristopov in ju primerjali z razpoložljivimi ab-initio podatki. Na podlagi dobljenih rezultatov smo postavili ustrezno teorijo za modeliranje kompleksnih heterogenih faznih mej med Ru-POM in elektrodami v raztopini.
Ma et al.: Rigid- and Polarizable-ion Potentials for Modeling