Unrestricted Floating Orbitals for the Investigation of Open Shell Systems

The floating orbital molecular dynamics approach treats the basis functions' centers in ab initio molecular dynamics simulations variationally optimized in space rather than keeping them strictly fixed on nuclear positions. An implementation of the restricted theory for closed shell systems is already available (Perlt et al., Phys. Chem. Chem. Phys., 2014, 16, 6997–7005). In this article, the extension of the methodology to the unrestricted theory in order to cover open shell systems is introduced. The methyl radical serves as a test system to prove the correctness of the implementation and to demonstrate the scope of this method. The available spin density plots and vibrational spectra are compared to those obtained from atom-centered bases. Finally, more complex systems as well as further properties to be studied in future investigations by floating orbitals are suggested.


Introduction
In recent years, the investigation of more and more complex systems by multi-scalar methods became feasible (Masson, Laino, Röthlisberger, & Hutter, 2009) (Ihrig, Schiffmann, & Sebastiani, 2011) (Kurzbach, Sharma, Sebastiani, Klinkhammer, & Hinderberger, 2011) (Golze, Iannuzzi, Nguyen, Passerone, & Hutter, 2013). This progress was only possible by introducing alternatives to standard methods. An example for such an alternative is the usage of floating orbitals which means that the centers of the basis functions for the construction of the wave function are not necessarily located on the nuclear positions but are optimized in space in order to minimize the total energy. A closed shell implementation of the floating orbital molecular dynamics (FOMD) approach has been presented in detail elsewhere (Perlt, Brüssel, & Kirchner, 2014). The general idea is to distinguish between centers of basis functions and nuclear coordinates when formulating the molecular integrals, which arise when applying the Hamiltonian to basis functions for the construction of the self-consistent field (SCF) equations. This is in contrast to conventional methods, where basis functions are centered on their respective nuclei. Having introduced another set of position variables, another gradient vector containing the derivative of the total energy with respect to the centers of the basis functionsthe electronic gradient -is available. This one is now used to minimize the total energy with respect to the basis functions' centers which is termed floating. Afterwards the molecular gradientthe derivative of the total energy with respect to nuclear coordinatesis evaluated as usual in order to determine forces which can be applied in the frame of an integration algorithm to perform molecular dynamics simulations. FOMD is a promising tool for several reasons. At first floating functions are capable of replacing additional local functions like bond functions (Neisius & Verhaegen, Bond functions for ab initio calculations on polyatomic molecules. Hydrocarbons, 1979), (Neisius & Verhaegen, Bond functions for ab initio calculations on polyatomic molecules. Molecules containing C, N, O and H., 1981) or augmented basis sets in the case of chemically complex systems (Pitarch-Ruiz, Evangelisti, & Maynau, 2005). Furthermore, floating orbitals are potentially helpful to describe non-nuclear attractors. These are maxima in electron density which cannot be assigned to a nucleus and have been observed for solvated electron clusters from theory (Timerghazin & Peslherbe, 2007) as well as experimentally by x-ray spectroscopy (Platts, Overgaard, Jones, Iversen, & Stasch, 2011). Finally, in molecular dynamics (MD) simulations the unwanted Pulay forces arise due to incomplete local basis sets and can be eliminated by the usage of floating orbitals (Marx & Hutter, 2009).

Method
At this point, a short note on restricted floating orbitals is given. As already stated, the centers of the basis functions are no longer fixed on nuclear positions but optimized in space. In ab initio MD simulations, the electronic structure is determined at each MD step, where for FOMD simulations the optimization of the basis functions' centers is additionally carried out. In order to optimize the basis positions, a gradient is calculated as the first derivative of the total energy with respect to the basis functions' centers. After successful optimization, the force on the nuclei is obtained from the first derivative of the total energy with respect to their respective position coordinates. The latter one is then applied to displace the atoms according to the integration algorithm. A simplified scheme of this approach is given in Fig. 1. A detailed derivation and further references can be found in (Perlt, Brüssel, & Kirchner, 2014). The extension of floating orbital molecular dynamics (FOMD) to open shell systems is now straightforward. Instead of considering one spatial orbital to comprise two electrons with α and β spin, respectively, spin orbitals are defined in the unrestricted Hartree-Fock (UHF) theory: where the spin variable ω has been introduced and separated from the spatial variable r. The variable x contains both spatial and spin coordinates, similarly denotes a spin orbital which is factorized to a spatial orbital and a spin value α or β. For either value of the spin variable, the one-electron eigenvalue equation for the one-electron Fock operator ̂ can be formulated: with denoting the energy eigenvalues. An analogous formulation for β spin is conceived, easily. Due to normalization of spin orbitals, the spin contribution in Eq. 2 can be eliminated by multiplication with the complex conjugate and integration. Thus for both spin coordinates, a set of eigenvalue equations determining the spatial orbitals are obtained, of which one is exemplarily given for an electron with index 1 in the orbital with index j possessing α spin: The one-electron Fock operator ̂ is composed of one-electron contributions ĥ , the exchange interaction with electrons of the same spin K and two-electron Coulomb repulsive interactions J with all remaining electrons: where and denote the number of α and β electrons, respectively. The formulation of the analogous operator for β spin is again straightforward. The spatial orbitals may now be expanded in a set of spin independent basis functions , whereas the expansion coefficients and are spin dependent: Finally, the well-known matrix SCF equations in their unrestricted formulation read where the matrix of overlap integrals S is spin independent. The Fock matrices and contain the energy integrals as obtained by the Fock operator as usual for electrons with alpha and beta spin, respectively. The vectors and contain the orbital energy eigenvalues. Eqs. 7 and 8 are mutually dependent on each other because of the dependence of both Fock operators on both spin variables. Consequently, they are solved in an iterative manner in accordance with the usual unrestricted Hartree-Fock approach. In analogy to the coefficient matrices, their product matricesthe density matrices and are spin dependent, as well. Their sum denotes the total density matrix tot , their difference, however, the spin density matrix . Finally, the total electronic energy 0 is obtained from To perform FOMD simulations, two distinct gradients are required: the molecular gradient for the force evaluation which is required for the integration algorithm in order to move the atoms in the molecular dynamics framework, and the electronic gradient. The latter one contains the first derivatives of the energy with respect to each basis function's center coordinates and is essential for the optimization within the floating orbital procedure. The components of both gradients can be formulated and rearranged to be composed of derivatives of molecular integrals. Those integrals comprise the components of the energy operator as well as the Gaussian type basis functions as defined by the basis set. Consequently, their derivatives are evaluated analogously to the closed shell case and afterwards combined with the spin dependent coefficient matrix elements. The formulation of the molecular gradient is textbook knowledge (Szabo & Ostlund, 1996) and the derivation of the electronic gradient is readily done according to the closed shell case. The optimization scheme used to minimize the energy with respect to orbital positions is adopted from the restricted case. That means, that as usual, the same basis set, i.e. identical functions at the same positions, is used for both α and β spin orbitals. Nevertheless, there is the possibility to separate the spin centers by occupying orbitals in different locations. This effect is obviously more pronounced for floated orbitals than for atom centered ones. In the following, the correctness of the implementation is demonstrated and certain properties are evaluated for the methyl radical.

Single Point Calculations
In order to demonstrate the correctness of the implementation, different calculations have been performed on the methyl radical which possesses nine electrons and therefore is necessarily an open shell system. Despite its small size, a number of theoretical studies investigating the structure and properties (Pacansky, Restricted and unrestricted Hartree-Fock calculations of the methyl radical, 1982), (Chipman, 1983) as well as spectra (Pacansky, Koch, & Miller, Analysis of the structures, infrared spectra, and Raman spectra for methyl, ethyl, isopropyl and tert-butyl radicals, 1991) of the methyl radical are available using Hartree-Fock methods. Another investigation considers the methyl radical to assess the accuracy of the spin-projected unrestricted Hartree-Fock method by investigating electron and spin densities (Glaser & Choy, 1993). Furthermore, reactions with other small molecules like carbon monoxide (Das & Lee, 2014), the recombination reaction with OH radicals (Oliveira & Bauerfeldt, 2012), and even the role in methanol synthesis (Zakharov, Ijagbuji, Tselishtev, Loriya, & Fedotov, 2015) have been studied recently. Due to its high reactivity, the methyl radical is involved in reactions which comprise the abstraction of a hydrogen atom so that many recent studies deal with the kinetic properties of these reactions (Mendes, Zhou, & Curran, 2014), (Saheb, 2015), (Tan, Yang, Krauter, Ju, & Carter, 2015). Firstly, our own FORTRAN 2003 code FORPLEASURE (Brüssel, Perlt, & Kirchner, 2012-2016 has been used to perform single point calculations with atom-centered basis functions, denoted with FP HF , as well as floating orbital single point energy calculations, denoted with FP FO . Both values are compared to the reference value as obtained with the program ORCA V. 3.0 (Neese, 2012) denoted as Orca HF . The molecular geometry has been the same for all calculations with all − bond lengths amounting to 108.9pm and all ∢( ) angles being 120°. The energies as obtained by UHF calculations are given in Table 1, where also the applied basis sets are given. The electron density plot has been generated with vmd (Humphrey, Dalke, & Schulten, 1996)  It is apparent from Table 1 that the present implementation in FORPLEASURE is capable of reproducing the results of the established and widely used program ORCA. Within the chosen SCF convergence criterion of 10 −7 ℎ , identical values for the total SCF energy are obtained. The application of a floating orbital optimization to the setup is able to reduce the SCF energy by 0.01 ℎ which may be considered as a significant stabilization. As the spin density is expected to change upon orbital displacement, spin density plots for the methyl radical as obtained with both basis sets, each with both atom-centered and floated orbitals, are provided in Figure 2.  Figure 2 shows that the spin density as obtained by floating orbitals is shifted to larger distances from the carbon nucleus for both basis sets. Especially for systems possessing lower symmetry and more polar bonds, the differences observed for the spin density are expected to be even larger. Finally, molecular dynamics simulations have been carried out to study the impact of the displaced basis functions on the molecular gradient and the resulting trajectory. The computational details of the simulations are given in Table 2. All simulations are carried out using the program FORPLEASURE. For the HFMD simulations, the 6-311G** basis set has been augmented by polarization functions so that the effect of floating orbitals can be compared to the effect of those additional functions.

UFO Molecular Dynamics Simulations
Simulations are carried out in the microcanonical ensemble at constant particle number N, volume V and total energy E. Velocities are initialized according to a Boltzmann distribution at the temperature given in Table 2. In order to study the influence of floating orbitals on dynamic properties, power and infrared (IR) spectra are studied in the following. Power spectra are obtained by TRAVIS (Brehm & Kirchner, 2011) as the Fourier transform of the velocity autocorrelation function. Therefore, the obtained data reveal the indirect effect of the floating orbitals on the trajectory of the system. Fig. 3 shows the signals at lower and higher wave numbers as well as a group of small signals www.ccsenet.org/ijc International Journal of Chemistry Vol. 8, No. 1;2016 around 400cm −1 in the inset. The positions of the peak maxima are given in Table 3, where reference values as obtained by static HF calculations with a 6-311G** basis set (Pacansky, Koch, & Miller, Analysis of the structures, infrared spectra, and Raman spectra for methyl, ethyl, isopropyl and tert-butyl radicals, 1991) as well as experimental data (Snelson, 1970) are given. The observed peak intensities in the power and IR spectra depend on various parameters such as temperature, simulation time and the amplitude of the oscillation, i.e. the changes in polarizability and dipole moment, respectively. Since the simulation times are not identical for all simulations, the resulting peak intensities will not be discussed in the following. The wavenumbers, at which the vibrational signals occur, denote the frequency of the oscillation. Considering the classical analogs, this is dependent on the reduced mass and the force constant of the oscillation. As the masses are not affected by the basis set, the observed shifts are due to different forces. Those arise from different basis setsby using larger basis sets, including polarization functions or applying floating basis functions. Consequently, different forces result in different trajectories showing different vibrational frequencies, which will be discussed in the following. Figure 3. Power spectra as obtained by FOMD and HFMD simulations with different basis sets. The inset shows the wavenumber range from 300cm −1 to 600cm −1 enlarged. Table 3. Positions of the peak maxima in the power spectrum of the methyl radical as observed by the different MD methods applied in this study as well as SCF (Pacansky, Koch, & Miller, Analysis of the structures, infrared spectra, and Raman spectra for methyl, ethyl, isopropyl and tert-butyl radicals, 1991) and experimental (Snelson, 1970)  Considering the two signals above 3000cm −1 , a slight shift towards smaller wavenumbers is observed by increasing the basis set for HFMD simulations from 3-21G to 6-311G. Both the inclusion of polarization functions in the 6-311G** basis set as compared to the 6-311G basis and floating of the orbitals do not change the positions of the signals significantly. The observed differences are not significant within the spectral resolution which is given by TRAVIS. However, this is not astonishing since the two signals can be assigned to C-H stretching vibrations. For the description of this nonpolar bond neither polarization functions nor floating orbitals are essential so that no significant effect is observed here. This is also the case for the in-plane bending vibration observed at ≈ 1500cm −1 . Although it is observed that improving the basis set by more basis functions, polarization, or floating, yields a shift to smaller wave numbers, the differences are hardly significant. Considering the signals below 1000cm −1 however, large differences are obtained for the peak positions. The underlying molecular vibration is the out-of-plane bending vibration and due to the presence of electron density above and below the molecular plane, the different basis sets yield different results. In that case, the effect observed by including polarization functions in the 6-311G basis set is similar to the effect obtained by floating the orbitals of the smaller basis: The peaks are shifted to smaller wavenumbers (Δ ≈ 80 − 100cm −1 ). This is in agreement with the observations from previous studies, where floating orbitals have also been found to cover the effect of polarization functions (Perlt, Brüssel, & Kirchner, 2014).
Furthermore, IR spectra are obtained from the autocorrelation function of the molecular dipole moment (Thomas, Brehm, Fligg, Vöhringer, & Kirchner, 2013). The dipole moment as a property directly related to the separation of charges has been found to be improved by the application of floating orbitals (Huber, 1981). Accordingly, the IR spectra are affected not only by the dynamics and their behavior with floating orbitals, but also by the dipole moment and its dependence on the basis functions' centers. The resulting IR spectra are given in Fig. 4, where again different intensities result from different simulation times and therefore will not be discussed.  For the IR spectra, results are similar to those obtained by the power spectra, see Table 4. All IR active modes observed in Refs. (Pacansky, Koch, & Miller, Analysis of the structures, infrared spectra, and Raman spectra for methyl, ethyl, isopropyl and tert-butyl radicals, 1991) and (Snelson, 1970) can be found by simulations and the corresponding peak positions are equal to those observed in the power spectra in Table 3. If systems with a permanent dipole moment are investigated, as planned in future studies, the dipole moment and consequently the IR spectra are expected to be improved by using floating orbitals. It is summarized that in agreement with the observations for the closed shell FOMD formulation (Perlt, Brüssel, & Kirchner, 2014), for all observed vibrational frequencies, the shifts as obtained by floating orbitals are similar to those obtained by additional polarization functions.

Summary and Outlook
The extension of floating orbitals to the unrestricted theory has been presented and the application to the methyl radical has been demonstrated. This rather small system has been chosen to rapidly evaluate the correctness of the implementation. However, in future studies, more complex systems will be examined. The great advantage of unrestricted floating orbital molecular dynamics is the possibility to investigate centers of charge, which do not necessarily need to be attributed to a nuclear center. Well-known problems of this type are solvated electrons. The floating orbital approach allows to assign a separate basis function to a solvated electron without placing it on a nucleus and optimize it. Thereby, solvated electrons may be treated by molecular dynamics and results can be compared to other theoretical studies (Uhlig & Jungwirth, Embedded Cluster Models for Reactivity of the Hydrated Electron, 2013), (Uhlig, Herbert, Coons, & Jungwirth, 2014) and experiments (Siefermann & Abel, 2011), (Abel, Buck, Sobolewski, & Domcke, 2012) which will be subject to a future study.