Effective electronic-only Kohn-Sham equations for the muonic molecules
1 Effective electronic-only Kohn-Sham equations for the muonic molecules
Milad Rayka , Mohammad Goli and Shant Shahbazian Department of Physics and Department of Physical and Computational Chemistry, Shahid Beheshti University, G. C., Evin, Tehran, Iran, 19839, P.O. Box 19395-4716. School of Nano Science, Institute for Research in Fundamental Sciences (IPM), Tehran 19395-5531, Iran
E-mails: Mohammad Goli : [email protected] Shant Shahbazian: [email protected] * Corresponding authors 2
Abstract
A set of effective electronic-only Kohn-Sham (EKS) equations are derived for the muonic molecules (containing a positively charged muon), which are completely equivalent to the coupled electronic-muonic Kohn-Sham equations derived previously within the framework of the Nuclear-Electronic Orbital density functional theory (NEO-DFT). The EKS equations contain effective non-coulombic external potentials depending on parameters describing muon’s vibration, which are optimized during the solution of the EKS equations making muon’s KS orbital reproducible. It is demonstrated that the EKS equations are derivable from a certain class of effective electronic Hamiltonians through applying the usual Hohenberg-Kohn theorems revealing a “duality” between the NEO-DFT and the effective electronic-only DFT methodologies. The EKS equations are computationally applied to a small set of muoniated organic radicals and it is demonstrated that a mean effective potential maybe derived for this class of muonic species while an electronic basis set is also designed for the muon. These computational ingredients are then applied to muoniated ferrocenyl radicals, which had been previously detected experimentally through adding muonium atom to ferrocene. In line with previous computational studies, from the six possible species the staggered conformer, where the muon is attached to the exo position of the cyclopentadienyl ring, is deduced to be the most stable ferrocenyl radical.
Keywords
Positively charged muon; muonic molecules; Nuclear-electronic orbital methodology; non-Coulombic density functional theory; effective theory 3
I. Introduction
The bound states of the muonic molecules, derived from the attachment of the positively charged muon, , to usual molecules, have been in focus in recent years through the spin rotational/relaxation/resonance (μSR) spectroscopy with vast applications from condensed matter physics to chemistry and molecular biology [1-8]. The basic measured quantities in the μSR spectroscopy are the hyperfine coupling constants which are used to locate the molecular site where the is attached [1,3]. However, the assignment procedure of the hyperfine coupling constants is not always an easy task since in complex molecules there are several sites potentially capable of trapping . Thus, it is desirable to derive theoretically both preferred sites for addition and their corresponding hyperfine coupling constants from quantum mechanical calculations [9-28]. The commonly-used molecular quantum mechanical methods are based on the clamped nucleus paradigm conceiving electrons as quantum particles and the nuclei as point charges [29,30]. Nevertheless, it is not generally evident that to what extent a light-mass , which is virtually one-ninth of proton’s mass, could be properly treated as a point charge [22-28]. One possible response to this concern is treating as a quantum particle like electron and trying to incorporate simultaneously the kinetic energy operators of electrons and into the Hamiltonian of the muonic molecules for quantum mechanical calculations. Recently, we have employed this strategy to consider several muonic molecules within the context of the Nuclear-Electronic Orbital (NEO) ab initio methodology [31-37]. As expected, one may extend such studies to more complex muonic molecules and try to develop more accurate ab initio computational procedures to achieve 4 experimental accuracy. Nevertheless, the recent “effective” reformulation of the NEO-Hartree-Fock (NEO-HF) equations [38,39], based on a simplified wavefunction proposed by Auer and Hammes-Schiffer [40], opens a new door to design “muon-specific” ab initio procedures. The present study is also a continuation of the same path trying to incorporate electron-electron (ee) correlation, which is absent in the effective NEO-HF (EHF) method, into the effective NEO theory. For this purpose, one may pursue two independent paths; trying to incorporate the ee correlation using various wavefunction-based post-NEO-HF procedures or using the NEO-density functional theory (NEO-DFT). In this report, the latter possibility is considered through introducing effective “electronic-only” Kohn-Sham (KS) equations for the muonic molecules while the former option will be considered in a subsequent report. The idea of extending DFT to systems containing multiple quantum particles, i.e., electrons plus at least one other type of quantum particles, is not new and have been tried for the electron-hole and the electron-positron systems many decades ago [41-43]. However, the seminal paper by Parr and coworkers is usually perceived as the first rigorous formulation of the multi-component DFT [44], though since then various extensions have been proposed [45-49]. Particularly, the recent interest in developing orbital-based ab initio non-Born-Oppenheimer procedures treating both electrons and nuclei as quantum particles from outset [31-33,50-53], triggered a renewed interest in practical reformulation and computational implementation of the multi-component DFT. After the pioneering studies [48,54], a number of papers have appeared dealing with various aspects of computational implementation and the functional design of the multi-component DFT in non-Born-Oppenheimer realm [55-69]. The NEO-DFT as proposed and extended by 5 Hammes-Schiffer and coworkers in recent years belongs to this category [56,62-68]. In the NEO-DFT, in contrast to the usual electronic DFT [70-73], one is faced with the coupled KS equations for each type of quantum particles, and therefore, each equation contains its own effective KS potential energy with some unique terms. These terms appear from multitude of the exchange-correlation functionals and while the previously designed electronic exchange-correlation functionals are used to describe the ee correlation, new functionals for other types of quantum particles and their interactions must be designed. In the case of the muonic molecules the basic types of correlations are the ee and eµ correlations and based on these correlations the muonic NEO-DFT may be divided into NEO-DFT(ee) and NEO-DFT(ee+eµ) categories. The latter theory is, in principle, exact though needs the proper introduction of the eµ correlation functional while the former theory basically neglects the eµ correlation. For systems containing quantum protons there have been several attempts to introduce electron-proton correlation functionals [48,49,59-67]. However, no systematic comparative study has been conducted yet to compare the relative merits of the proposed functionals though many (but not all) of them seem to be inspired in some way from the original Colle-Salvetti formula for the ee correlation [74]. In this study, we will neglect the eµ correlation and an effective formulation of the NEO-DFT(ee) is proposed, similar to that proposed for the NEO-HF theory [38,39]. The paper is organized as follows. In section II the effective NEO-DFT(ee) and corresponding KS equations are discussed while our computational details are provided in Section III. In section IV, the optimized electronic basis sets used for are introduced through considering a small but representative set of organic molecules. In addition, the computational implementation of the effective equations for large molecules with an 6 example is studied to demonstrate the efficiency of the proposed theory. Finally, our conclusions are offered in Section V. II. Theory
The Hamiltonian for a system containing N electrons, a single massive quantum positively charged quantum particle (PCP) with mass m (assuming e m m ) and q clamped nuclei, is as the following in the atomic units ( e m ): ˆ ˆ q qtotal NEO Z ZH H R R
N N N NNEO i exti i i j ii i j
H m Vr r r r q qNext i i
Z ZV R r R r (1) This is a two-component NEO Hamiltonian that includes the basic physics of the muonic molecules if the mass of the PCP to be fixed at the mass of ( e m m ). Since the formalism of the NEO-DFT has been extensively discussed for the general multi-component systems, herein, only the main points are discussed briefly and the interested reader may consult the original literature for details [56,62-67]. In principle, not only but also all nuclei may be treated as quantum particles and a DFT formalism for resulting “self-bound” muonic system can be devised [75-78], however, this generalized formalism is beyond the scope of this paper. It is possible to demonstrate that the external potential, ext V , is uniquely determined by the following one-particle densities: ... e N spins r N dr dr dr , ... N spins r dr dr (2) 7 In these equations and its complex conjugate stand for the muonic ground state wavefunction while the resulting universal functional in the Levy’s constrained search procedure depends on both densities and the mass of the PCP (for a very brief but lucid discussion see section 9.6 in [72]). Assuming the KS reference system to be a non-interacting system of electrons and the PCP the KS wavefunction is a product of a determinant composed of the KS spin-orbitals for electrons and a single spatial orbital for the PCP. The total energy of the system for a typical closed-shell electronic system, neglecting eµ correlation, is the following: ˆ ˆ2 NNEO i i i i i ii
E dr r h r dr r h r - eee e exc e r rJ E dr dr r r ˆ 1 2 qi i i Zh R r , ˆ 1 2 q Zh m R r e eee e r rJ dr dr r r (3) In the preceding equations exc E , the exchange-correlation energy functional, stands for all so-called “non-classical” effects resulting from the electron exchange, ee correlation and residual electronic kinetic energy beyond the KS reference system [71-73]. If the energy is varied with respect to the electronic and muonic spatial orbitals the following coupled KS equations are derived:
21 1 1 1 eKS j j j v r r r j N KS m v r r r q eeKS exc rZ rv r dr dr v rr r r rR r q eKS Z rv r dr r rR r exc eexc e Ev (4) These differential equations are usually transformed into algebraic equations by expanding all orbitals through known basis sets and then solving algebraic equations employing the self-consistent field procedure (extension to open-shell electronic systems is also straightforward and is not considered herein). As discussed in the previous communications [38,39], since the PCP is well-localized due to its large mass relative to that of electron, the orbital of the PCP may be approximated with a wavefunction describing a 3D quantum oscillator. The simplest example is an isotropic harmonic oscillator with the following ground state wavefunction: c r Exp r R , where is the width and c R is the center of the Gaussian-type function (GTF), which are the standard parameters of a quantum oscillator [38]. It is possible to use more complicated anharmonic and anisotropic oscillator models with a large set of parameters instead, as discussed in detail elsewhere [39], however, the isotropic harmonic oscillator may be employed as an illustrative example of the formulation of the effective NEO-DFT(ee) (for a similar idea see [79]). Incorporating this s-type GTF in equation (3) yields the following expression for the energy: ˆ2 2 N eNEO i i i i i i ci c rE dr r h r dr erf r Rr R ee e exc e J E U
3, , 22 qc cc
ZU R R erf R Rm R R (5) In this equation, the PCP disappears as a quantum particle and instead novel non-coulombic potential energy terms appear containing the parameters of the original GTF. Particularly, U can be conceived as a classical potential energy term similar to the nuclear repulsion term in equation (1) and does not need to be considered in subsequent functional variation. Thus, the total and effective KS energies are now: total eff KS classical E E E ˆ2 2 N eeff KS i i i i i i ci c rE dr r h r dr erf r Rr R ee e exc e J E , , q qclassical c
Z ZE R R U R R (6) Upon variation of eff KS E with respect to the electronic orbitals the following effective electronic-only KS equations, hereafter briefly called the EKS equations, arise:
21 1 1 1 eeff KS j j j v r r r j N q ee eeff KS eff exc Z rv r v r dr v rr rR r eeff cc v r erf R rR r (7) Formally, solving the EKS equations is equivalent to the solution of the coupled equations offered in equation (4) with a single GTF as a basis set. In practice the price that has been 10 payed is adding two new parameters of the GTF, i.e. and c R , in the optimization procedure of classical E along the usual geometry optimization of the clamped nuclei. Like the case of equations (4) these differential equations may be transformed into algebraic equations by expanding the electronic KS orbitals in GTFs. The one-electron integrals resulting from the electron-PCP interaction potential energy term, eeff v r , are available analytically [80]. Still, even in the absence of analytical formulas, the same numerical integration procedure employed to derive the integrals associated to exc v r may also be used to evaluate this type of integrals [81-84]. At this stage of development, the algebraic EKS equations can be solved by using any known electronic exchange-correlation functionals and basis sets for the study of muonic system. However, before considering computational implementation, let us try to grasp certain ramifications of the effective formulation. Although equations (7) were derived assuming the two-component Hamiltonian and associated NEO-DFT, it is possible to reverse the procedure and try to construct an effective electronic Hamiltonian that yields equations (7) directly through its own DFT (for a comprehensive discussion on “building” DFTs for “model” Hamiltonians see [85]). The following Hamiltonian is a proper candidate:
3ˆ ˆ 22 q q qeff elec eff cc
Z Z ZH H erf R RmR R R R
1ˆ 1 2
N N Nelec eff i exti i j i i j
H Vr r qN Next c ii ii c i
ZV erf R rR r R r (8) 11 It is straightforward to demonstrate that the whole machinery of the electronic DFT is equally applicable, and equations (7) arise assuming a non-interacting electronic KS reference system for the effective electronic Hamiltonian. It is interesting to try to generalize this result and introduce a picture not tied to specific vibrational models used for the PCP (for a discussion on complicated vibrational models beyond harmonic oscillator see [39]). Accordingly, there is a correspondence between the effective potential and the Gaussian or Slater type basis sets (or any other well-designed mathematical function) used to expand the orbital of the PCP. This correspondence originates from integrating the kinetic energy integral of the PCP and the electron-PCP interaction integrals which leads to the remaining of only basis function parameters, denoted as k c . Thus, the resulting correspondence of the PCP orbital and the effective interaction potential is: ; ; , ek eff k k r c v r c U c [39]. Hence, the most general form of the EKS equations, reiterating equations (5) to (7), is:
21 1 1 1 eeff KS j j j v r r r j N ; q ee eeff KS eff k exc Z rv r dr v r c v rr rR r , q qclassical k k Z ZE c R U c R R (9) In these equations k c must be optimized like the nuclear geometry and using the optimized parameters, one may reconstruct the orbital of the PCP and corresponding one-particle density, which describes the vibrational motion of the PCP. The corresponding effective Hamiltonian may be written as the following: 12 ˆ ˆ q qeff elec eff k Z ZH H U cR R
1ˆ 1 2
N N Nelec eff i exti i j i i j
H Vr r ; qN N eext eff i ki ii ZV v r cR r (10) Equations (9) and (10) are the heart of the effective NEO-DFT for the muonic systems though it can be used also for systems containing a single quantum proton or any other heavier particle and may easily be extended to the multi-component cases with more than two types of quantum particles as well.
III. Computational details
Based on the proposed effective DFT in order to start solving the EKS equations, at first step, the effective potentials must be introduced. In this study the used potentials are constructed from a fully-optimized single s-type GTF, [1s], and from a scaled [5s5p] Gaussian basis set; the former has been explicitly given in equation (7). It must be emphasized that instead of deriving the potential separately for each basis set, an automated algorithm may be constructed to produce the potential after determining the type and number of GTFs used to expand the muonic orbital (Goli and Shahbazian, under preparation). For describing electronic distribution, 6-311++g(d,p) basis set was placed on the clamped nuclei [86-88], while for the muon a [4s1p] electronic basis set was placed at a banquet atom; in the [1s] associated effective potential, this center is denoted by c R in equations (7). In all ab initio calculations, the B3LYP exchange-correlation hybrid functional was employed for electrons without further modifications or re-optimization of 13 its parameters [89-91] (we leave the possibility of reparametrizing this functional for the muonic systems to a future study). The whole computational level is termed EKS-B3LYP/[6-311++g(d,p)/4s1p]. The optimized parameter of the effective potential, , in equations (7) as well as the energy-optimized exponents of [4s1p] electronic basis set of a representative set of organic molecules ( vide infra ) were determined through a full optimization of the EKS equations using a non-linear numerical optimization procedure [39]. In the case of [5s5p] associated effective potential partial optimization was done and only half of parameters, the linear coefficients ( vide infra ), were determined through direct optimization as discussed in the next section. Besides, the energy-optimized exponents of [4s1p] electronic basis derived from the [1s] associated effective potential were used without further optimization for the EKS calculations with [5s5p] associated effective potential. The geometry of the clamped and banquet nuclei was optimized using the analytical gradients during the geometry optimization procedure. Since all considered muonic species are odd-electron systems, the unrestricted (U) as well as the restricted-open (RO) versions of the algebraic EKS equations were utilized for ab initio calculations. A modified version of the GAMESS package was used for all ab initio calculations and the original implemented numerical integration algorithm for the exchange-correlation integrals was employed without any modifications [92,93]. Throughout the calculations, the used masses for , proton (H), deuterium (D) and tritium (T) are 206.76828, 1836.15267, 3670.48296 and 5496.92153, respectively, in atomic units. IV. Results and discussion A.
Designing the effective potential and the electronic basis set
14 In the previous EHF study on a large set of species it was demonstrated that the parameters of the effective potential as well as the exponents of the electronic basis set are mainly determined by the mass of the PCP and relatively insensitive to the chemical environment [39]. To have a clear picture of the mass-dependence of the effective potential in the EKS equations a comparative study was done on
XCN ( , , , X H D T ) species, where the carbon and nitrogen nuclei were considered as clamped point charges while X was treated as a quantum particle. Figure 1 offers the final optimized results of the EKS calculations revealing that the effective potential associated to the heavier particle has less deviations from the point charge coulombic potential far from the banquet center while its one-particle density is more localized. These observations are in line with the expectation that a massive quantum particle behaves more like a clamped point charge than a lighter one. To have a quantified picture of these variations the effective potentials in equations (5) and (7) are rewritten using the known mass-dependence of from the isotropic harmonic oscillator model [94]: eeff i i ci c v r erf km r Rr R
3, , , 4 qc cc
ZkU k m R R erf km R Rm R R (11) In these potentials k stands for the force constant which is mainly characteristic of the environment and theoretically, independent from the mass. Thus, the model is to be taken seriously only when slight variation of this effective force constant is seen upon isotopic substitution in the EKS calculations. To test the reliability of the model, Figure 1 compares the mass-dependence of the optimized and k values demonstrating that the latter is 15 indeed much less sensitive to the mass variations, and may even be treated as a constant in the case of the hydrogen isotopes. This observation reveals that equations (11) are reliable for modeling the intrinsic mass-dependence of the effective potential while the smaller variations of the potentials induced by the variations of the environment are all absorbed in the effective force constant. If ab initio EKS calculations demonstrate that the variations of the effective force constant in a certain set of muonic species are also small, then it is legitimate to introduce a “mean” and corresponding mean effective potential for the corresponding set. The introduction of the mean effective potentials particularly bypasses the costly non-linear optimization of for each species separately. Seven small organic molecules: diazene, acetylene, methenamine, hydrogen cyanide, formamide, fomaldehyde, and ethylene, where their atom types and bonds are typical to many organic molecules used in the experimental μSR studies [1,3], were selected as the representative set. Since in the corresponding experiments muonium atom ( plus an electron) is attached to organic molecules, we have also attached muonium atom to these seven molecules and the resulting set of open-shell muonic species was employed in order to evaluate the numerical values of . Practically, to have an initial geometry, a hydrogen atom with a clamped proton was first attached to the molecules of the representative set and from various possible conformers the lower energy minima were extracted after the geometry optimization at B3LYP/6-311++g(d,p) level. In next step, a muonium atom replaced the hydrogen atom, i.e. a banquet atom and corresponding electronic [4s1p] basis functions were added instead of hydrogen atom. Subsequently, the position of the banquet atoms, and the exponents of [4s1p] electronic basis set were simultionsly optimized while the geometry of the remaining nuclei was held fixed. The 16 resulting eleven low-energy muoniated structures are depicted in Figure 2 while the corresponding final optimized geometries have been given in the supporting information. Table 1 offers the optimized values at both EKS-UB3LYP/[6-311++g(d,p)/4s1p] and EKS-ROB3LYP/[6-311++g(d,p)/4s1p] levels of calculations; the resulting optimized values are distributed narrowly, EKS , and insensitive to the U or the RO versions of the EKS calculations. This is not far from the mean derived from the EHF calculations in the previous study on a completely different set of closed-shell muonic species [39], EHF , and confirms that the idea of the mean effective potential is reliable enough to be used in the EKS calculations. On the other hand, in the previous EHF study the mean exponents of the basis functions in the extended muonic [2s2p2d] basis set were narrowly distributed around the mean:
EHF EHF [39]. In present study the exponents of [5s5p] muonic basis set were also scaled around the mean:
EKS EKS EKS EKS EKS . Then, they are used in the construction of corresponding extended effective potential, which in contrast to the simple effective potentials in equations (5) and (7), includes the anharmonicity and the anisotropy of ’s vibrations (for a thorough discussion see [39]). Hence, only the linear coefficients in the extended effective potential were optimized during the EKS calculations (see the supporting information of [39] for details of the procedure). The exponents of the electronic [4s1p] basis set were also optimized both at EKS-UB3LYP/[6-311++g(d,p)/4s1p] and EKS-ROB3LYP/[6-311++g(d,p)/4s1p] levels and the final results have been gathered in Table 2. Once again, it is evident from this table that the optimized exponents are relatively insensitive to chemical environment and the mean values are proper representatives for the optimized exponents. Therefore, the mean values derived at 17 the EKS-UB3LYP/[6-311++g(d,p)/4s1p] level will be used in the rest of this paper. The computed mean exponents are comparable to the mean exponents derived in the previous study using the EHF equations: s p [39], and are insensitive to the U and the RO versions of the EKS calculations as well. B. Muoniated ferrocenyl radicals
The idea of adding muonium atoms to organic molecules is not new [1], by the way, more recently muonium atoms have been attached also to carbenes, their organosilicon analogs and organometallic molecules [95,96]. One of the interesting organometallic targets considered both experimentally and computationally is the iconic ferrocene molecule though the corresponding experimental μSR spectrum of its muoniated radical is not yet conclusively assigned [97-99]. Taking the size of ferrocene molecule, all previous computational studies employed DFT to model this system where instead of the muonium atom, a hydrogen atom with a clamped nucleus has been added to ferrocene [97,99]. In this section the EKS-UB3LYP/[6-311++g(d,p)/4s1p] method is used in conjunction with the extended effective potential developed in the previous section in order to study the muoniated ferrocenyl radicals. In order to start the calculations, we reoptimized the ferrocenyl radical structures reported by McKenzie at UB3LYP/6-311++g(d,p) level [99]. The four optimized structures include ferrocenyl radicals after hydrogen addition to cyclopentadienyl ring or the iron atom while considering the relative configuration of the cyclopentadienyl (Cp) rings that could be staggered or eclipsed [99]. In the next step, the hydrogen atom was eliminated from the structures and a banquet atom with a [4s1p] basis set was added instead. In the case of radicals originating from adding hydrogen atom to the Cp ring, 18 banquet atom could be placed at exo, i.e. with the farthest distance to the iron center, or endo, i.e. with the closest distance to the iron center, positions. Therefore, six distinct muoniated structures were prepared and the EKS-UB3LYP/[6-311++g(d,p)/4s1p] method along with the geometry optimization were applied. Figure 3 depicts the final optimized structures and the used nomenclature while supporting information contains the optimized coordinates. In the case of exo-Cp-eclipsed, endo-Cp-eclipsed and Fe-staggered structures, the full geomtry optimization did not yield stable structures thus they were derived by imposing a plane of symmetry as a contriant. Indeed, the hydorgenic analogs of these structures are saddle points on the corresponding energy hypersurface as reproted by McKenzie and independently confirmed in the present study [99]. Figure 3 also contains the “mean “distance between and neighboring nucleus; it is important to stress that this mean distance is distinct from the distance between banquet center and neighboring nucleus and is the expectation value of ’s position operator where the neighboring nucleus is at the center of the coordinate system. The mean C- distances, 1.16 - 1.18 Å, if compared to usual C-H distances, 1.08 - 1.10 Å, are quite longer revealing the nonnegligible elongations of the bond lengths upon isotopic substitution. This trend is completely absent in the conventional ab initio calculations based on using a clamped hydrogen atom to model muonium addition to molecular systems. This observation is in line with one of our previous studies where substitution of with one of the hydrogen atoms of malonaldehyde varied the conformation drastically [37]. All observed traits point to the fact that even at the structural level the EKS results contain novel information that are absent if the clamped nucleus model be used instead. Table 3 includes the total and relative energies for the considered species revealing that exo-Cp-staggered radical is the 19 most stable species in line with McKenzie’s results [99]. As stressed before, (exo/endo)-Cp-ecllipesed radicals are not stable structures, and threfore, the next stable structure is the exo-Cp-staggered radical, which is only ~3 kJ/mol less stable than the correspodning exo radical whereas the Fe-ecplised radical is ~53 kJ/mol less table than the exo. Thus, it seems reasonable to assume that only the Cp-staggered conformer has a major contribution in the muoniation in gas phase though as disucssed in detail by McKenzie [99], in the solid-state uncertainties emerge due to the crystal packing effects and possible structural deformations of ferreocne molecule itself. V. Conclusion
The present study demonstrates that the basic idea of the effective theory, proposed recently [38,39], is extendable both theoretically and computationally to the DFT of muonic systems. Nevertheless, this is just a primary result and it is desirable to include ee correlation at the wavefunction-based levels of the effective theory as well as devising the effective version of the NEO-DFT(ee+eµ), both will be discussed in subsequent reports. In the meantime, the EKS equations provide a framework to start studying the structural and energetic aspects of large muonic systems though without including eµ correlation the quantitative prediction of the μSR spectrum remains yet elusive. The EKS equations may also yield the required information, e.g. one-particle densities, for “atoms in molecules” study of muonic species as considered in some previous reports [34-37]. This path is also now under scrutiny in our laboratory and the results will be discussed in future communications.
Conflicts of interest
There are no conflicts of interest to declare. 20
Acknowledgments
The authors are grateful to Cina Foroutan-Nejad for detailed reading of this paper.
References [1] E. Roduner,
The Positive Muon as a Probe in Free Radical Chemistry (Lecture Notes in Chemistry, Vol. 49, Springer, Berlin, 1988). [2] S. J. Blundell, Contemp. Phys. , 175 (1999). [3] K. Nagamine, Introductory Muon Science (Cambridge University Press, Cambridge, 2003). [4] S. J. Blundell, Chem. Rev. , 5717 (2004). [5] J. H. Brewer, Phys. Proc. , 2 (2012). [6] L. Nuccio, L. Schulz, and A. J. Drew, J. Phys. D: Appl. Phys. , 473001 (2014). [7] J. Phys. Soc. Jap. No. 9 (2016), Special Topics: Recent Progress of Muon Science in Physics, 091001-091016. [8] K. Wang, P. Murahari, K. Yokoyama, J. S. Lord, F. L. Pratt, J. He, L. Schulz, M. Willis, J. E. Anthony, N. A. Morley, L. Nuccio, A. Misquitta, D. J. Dunstan, K. Shimomura, I. Watanabe, S. Zhang, P. Heathcote and A. J. Drew, Nature Mat. , 467 (2017). [9] C. G. Van de Walle, Phys. Rev. Lett. , 669 (1990). [10] A. R. Porter, M. D. Towler, and R. J. Needs, Phys. Rev. B , 13534 (1999). [11] D. Cammarere, R. H. Scheicher, N. Sahoo, T. P. Das, and K. Nagamine, Physica B , 636 (2000). [12] C. P. Herrero and R. Ramírez, Phys. Rev. Lett. , 205504 (2007). [13] H. Maeter, H. Luetkens, Yu. G. Pashkevich, A. Kwadrin, R. Khasanov, A. Amato, A. A. Gusev, K. V. Lamonova, D. A. Chervinskii, R. Klingeler, C. Hess, G. Behr, B. Büchner, and H.-H. Klauss, Phys. Rev. B , 094524 (2009). [14] E. V. Sheely, L. W. Burggraf, P. E. Adamson, X. F. Duan, and M. W. Schmidt, J. Phys.: Conf. Ser. , 012049 (2010). [15] E. L. Silva, A. G. Marinopoulos, R. C. Vilão, R. B. L. Vieira, H. V. Alberto, J. Piroto Duarte, and J. M. Gil, Phys. Rev. B , 165211 (2012). 21 [16] F. Bernardini, P. Bonfâ, S. Massidda, and R. De Renzi, Phys. Rev. B , 115148 (2013). [17] P. Bonfà, F. Sartori and R. De Renzi, J. Phys. Chem. C , 4278 (2015). [18] R. B. L. Vieira, R. C. Vilão, A. G. Marinopoulos, P. M. Gordo, J. A. Paixão, H. V. Alberto, J. M. Gil, A. Weidinger, R. L. Lichti, B. Baker, P. W. Mengyan, and J. S. Lord, Phys. Rev. B , 115207 (2016). [19] S. J. Blundell, J. S. Möller, T. Lancaster, P. J. Baker, F. L. Pratt, G. Seber, and P. M. Lahti, Phys. Rev. B , 064423 (2013). [20] F. R. Foronda, F. Lang, J. S. Möller, T. Lancaster, A. T. Boothroyd, F. L. Pratt, S. R. Giblin, D. Prabhakaran, and S. J. Blundell, Phys. Rev. Lett. , 017602 (2015). [21] P. Bonfà and R. De Renzi, J. Phys. Soc. Jap. , 091014 (2016). [22] R. M. Valladares, A. J. Fisher, and W. Hayes, Chem. Phys. Lett. , 1 (1995). [23] M. I. J. Probert, and A. J. Fisher, Chem. Phys. Lett. , 271 (1996). [24] M. I. J. Probert, and A. J. Fisher, J. Phys.: Condens. Matter , 3241 (1997). [25] Y. K. Chen, D. G. Fleming, and Y. A. Wang, J. Phys. Chem. A , 2765 (2011). [26] D. G. Fleming, D. J. Arseneau, M. D. Bridges, Y. K. Chen, and Y. A. Wang, J. Phys. Chem. C , 16523 (2013). [27] K. Yamada, Y. Kawashima, and M. Tachikawa, J. Chem. Theory Comput. , 2005 (2014). [28] Y. Oba, T. Kawatsu, and M. Tachikawa, J. Chem. Phys. , 064301 (2016). [29] A. Szabo and N. S. Ostlund,
Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Publications Inc., New York, 1996). [30] T. Helgaker, P. Jørgenson, and J. Olsen,
Molecular Electronic-Structure Theory (John Wiley & Sons, New York, 2000). [31] S. P. Webb, T. Iordanov, and S. Hammes-Schiffer, J. Chem. Phys. , 4106 (2002). [32] T. Iordanov and S. Hammes-Schiffer, J. Chem. Phys. , 9489 (2003). [33] C. Swalina, M. V. Pak, and S. Hammes-Schiffer, Chem. Phys. Lett. 404, (2005). [34] M. Goli and Sh. Shahbazian, Phys. Chem. Chem. Phys . , 6602 (2014). [35] M. Goli and Sh. Shahbazian, Phys. Chem. Chem. Phys . , 245 (2015). 22 [36] M. Goli and Sh. Shahbazian, Phys. Chem. Chem. Phys . , 7023 (2015). [37] M. Goli and Sh. Shahbazian, Chem. Eur. J . , 2525 (2016). [38] M. Gharabaghi and Sh. Shahbazian, Phys. Lett. A , 3983 (2016). [39] M. Ryaka, M. Goli, and Sh. Shahbazian, Phys. Chem. Chem. Phys .
20, 4466 (2018). [40] B. Auer and S. Hammes-Schiffer, J. Chem. Phys. , 084110 (2010). [41] L. M. Sander, H. B. Shore and L. J. Sham, Phys. Rev. Lett. , 533 (1973). [42] R. K. Kalia and P. Vashishta, Phys. Rev. B , 2655 (1978). [43] R. M. Nieminen, E. Boroński and L. J. Lantto, Phys. Rev. B , 1377 (1985). [44] J. F. Capitani, R. F. Nalewajski and R. G. Parr, J. Chem. Phys. , 568 (1983). [45] H. Englisch and R. Englisch, Phys. Stat. Sol. , 711 (1984). [46] E. S. Kryachko, E. V. Ludeńa and V. Mujica, Int. J. Quantum Chem. , 589 (1991). [47] N. Gidopoulos, Phys. Rev. B , 2146 (1998). [48] T. Kreibich and E. K. U. Gross, Phys. Rev. Lett. , 2984 (2001). [49] T. Kreibich, R. van Leeuwen and E. K. U. Gross, Phys. Rev. A , 022501 (2008). [50] M. Tachikawa, K. Mori, H. Nakai and K. Iguchi, Chem. Phys. Lett. , 437 (1998). [51] H. Nakai, Int. J. Quantum Chem. , 2849 (2007). [52] T. Ishimoto, M. Tachikawa and U. Nagashima, Int. J. Quantum Chem. , 2677 (2009). [53] R. Flores-Moreno, E. Posada, F. Moncada, J. Romero, J. Charry, M. Díaz-Tinoco, S. A. González, N. F. Aguirre and A. Reyes, Int. J. Quantum Chem. , 50 (2014). [54] Y. Shigeta, H. Takahashi, S. Yamanaka, M. Mitani, H. Nagao and K. Yamaguchi, Int. J. Quantum Chem. , 659 (1998). [55] T. Udagawa and M. Tachikawa, J. Chem. Phys. , 244105 (2006). [56] M. V. Pak, A. Chakraborty, and S. Hammes-Schiffer, J. Phys. Chem. A , 4522 (2007). [57] Y. Kita and M. Tachikawa, Comput. Theor. Chem. , 9 (2011). [58] Y. Kita, H. Kamikubo, M. Kataoka and M. Tachikawa, Chem. Phys. , 50 (2013). 23 [59] T. Udagawa, T. Tsuneda and M. Tachikawa, Phys. Rev. A , 052519 (2014). [60] Y. Imamura, H. Kiryu and H. Nakai, J. Comput. Chem. , 735 (2008). [61] Y. Imamura, Y. Tsukamoto, H. Kiryu and H. Nakai, Bull. Chem. Soc. Jpn. , 1133 (2009). [62] A. Chakraborty, M. V. Pak, and S. Hammes-Schiffer, Phys. Rev. Lett. , 153001 (2008). [63] A. Chakraborty, M. V. Pak, and S. Hammes-Schiffer, J. Chem. Phys. , 124115 (2009). [64] A. Sirjoosingh, M. V. Pak, and S. Hammes-Schiffer, J. Chem. Theory Comput. , 2689 (2011). [65] A. Sirjoosingh, M. V. Pak, and S. Hammes-Schiffer, J. Chem. Phys. , 174114 (2012). [66] T. Culpitt, K. R. Brorsen, M. V. Pak, and S. Hammes-Schiffer, J. Chem. Phys. , 044106 (2016). [67] Y. Yang, K. R. Brorsen, T. Culpitt, M. V. Pak and S. Hammes-Schiffer, J. Chem. Phys. , 114113 (2017). [68] K. R. Brorsen, Y. Yang, and S. Hammes-Schiffer, J. Phys. Chem. Lett. , 3488 (2017). [69] R. Requist and E. K. U. Gross, Phys. Rev. Lett. , 193001 (2016). [70] P. Hohenberg and W. Kohn, Phys. Rev. B , 864 (1964). [71] W. Kohn and L. J. Sham, Phys. Rev. A , 1133 (1965). [72] R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules (Oxford University Press, Oxford, 1989). [73] E. Engel and R. M. Dreizler,
Density Functional Theory: An advanced Course (Springer, Heidelberg, 2011). [74] R. Colle and O. Salvetti, Theoret. Chim. Acta (Berl.) , 329 (1975). [75] J. Engel, Phys. Rev. C , 014306 (2007). [76] N. Barnea, Phys. Rev. C , 067302 (2007). [77] J. Messud, M. Bender and E. Suraud, Phys. Rev. C , 054314 (2009). [78] J. Messud, Phys. Rev. C , 052113 (2011). [79] T. Shimazaki and M. Kubo, Chem. Phys. Lett. , 134 (2012). 24 [80] P. M. W. Gill and R. D. Adamson, Chem. Phys. Lett. , 105 (1996). [81] A. D. Becke, J. Chem. Phys. , 2547 (1988). [82] A. D. Becke, J. Chem. Phys. , 2993 (1988). [83] C. W. Murray, N. C. Handy, and G. J. Laming, Mol. Phys. , 997 (1993). [84] O. Treutler and R. Ahlrichs, J. Chem. Phys. , 346 (1995). [85] K. Capelle and V. L. Campo, jr., Phys. Rep. , 91 (2013). [86] W. J. Hehre, R. Ditchfield, and J. A. Pople, J. Chem. Phys. , 2257 (1972). [87] P. C. Hariharan and J. A. Pople, Theoret. Chim. Acta , 213 (1973). [88] T. Clark, J. Chandrasekhar, G.W. Spitznagel, and P. v. R. Schleyer. J. Comp. Chem. , 294 (1983). [89] A. D. Becke, Phys. Rev. A , 3098 (1988). [90] C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B , 785 (1988). [91] A. D. Becke, J. Chem. Phys. , 5648 (1993). [92] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis and J. A. Montgomery, J. Comput. Chem. , 1347 (1993). [93] M. S. Gordon and M. W. Schmidt, Theory and Applications of Computational Chemistry: the first forty years , pp. 1167-1189, Eds. C. E. Dykstra, G. Frenking, K. S. Kim, G. E. Scuseria (Elsevier, Amsterdam, 2005). [94] M. Goli and Sh. Shahbazian, Theor. Chem. Acc. , 1362 (2013). [95] R. West and P. W. Percival, Dalton Trans. , 9209 (2010). [96] R. West, K. Samedov and P. W. Percival, Chem. Eur. J. , 9184-9190 (2014). [97] R. M. Macrae, Physica B , 307 (2006). [98] U. A. Jayasooriya, R. Grinter, P. L. Hubbard, G. M. Aston, J. A. Stride, G. A. Hopkins, L. Camus, I. D. Reid, S. P. Cottrell and S. F. J. Cox, Chem. Eur. J. , 2266 (2007). [99] I. McKenzie, Phys. Chem. Chem. Phys. , 10600 (2014). 25 Tables:
Table 1- The optimized values (in atomic units) derived at EKS-UB3LYP/[6-311++g(d,p)/4s1p] and EKS-ROB3LYP/[6-311++g(d,p)/4s1p] levels for the eleven representative set of the muoniated species depicted in Figure 2. backbone U RO Acetylene Diazene
Ethylene
C-Formaldehyde
O-Formaldehyde
C-Formamide
O-Formamide
C-Hydrogen cyanide
N-Hydrogen cyanide
C-Methenamine
N-Methenamine mean
U S S S S P Acetylene
Diazene
Ethylene
C-Formaldehyde
O-Formaldehyde
C-Formamide
O-Formamide
C-Hydrogen cyanide
N-Hydrogen cyanide
C-Methenamine
N-Methenamine mean
RO S S S S P Acetylene
Diazene
Ethylene
C-Formaldehyde
O-Formaldehyde
C-Formamide
O-Formamide
C-Hydrogen cyanide
N-Hydrogen cyanide
C-Methenamine
N-Methenamine mean
Radicals E ∆E exo-Cp-eclipsed -1651.35270 2.4 endo-Cp-eclipsed -1651.35153 5.4
Fe-eclipsed -1651.33346 52.9 exo-Cp-staggered -1651.35360 0.0 endo-Cp-staggered -1651.35246 3.0
Fe-staggered -1651.33227 56.0 28 Figure 1- The one-particle density and the optimized effective potential (see equation (7) in the text) of
XCN ( , , , X H D T ) series of species obtained from the EKS-B3LYP/[6-311++g(d,p)/4s1p] calculations. The black line in the effective potential panel is the point charge coulombic potential (proportional to the inverse of distance). The optimized parameter of the effective potential, i.e. the exponent of the s-type GTF, and the corresponding effective force constants are: , H , D , T , and k , H k , D k , T k , respectively. All quantities are given in atomic units. 29 Figure 2- The structures of muoniated (a) diazene, (b) acetylene, (c) N-methenamine, (d) C-methenamine, (e) N-hydrogen cyanide, (f) C-hydrogen cyanide, (g) O-formamide, (h) C-formamide, (i) O-fomaldehyde, (j) C-fomaldehyde, and (k) ethylene adducts (N, C and O symbols are used to descriminate muon’s attachment site). Blue, grey, red, white and green spheres indicate the locations of nitrogen, carbon, oxygen, hydrogen and muonium banquet atoms, respectively. 30 Figure 3- The equilibriumm structures of ferrocene muoniated adducts (a) exo-Cp-eclipsed, (b) endo-Cp-eclipsed, (c) Fe-eclipsed, (d) exo-Cp-staggered, (e) endo-Cp-staggered and (f) Fe-staggered. The mean inter-nuclear distance between muonium (banquet) atom and its binding site (neigboring nucleus) for each adduct has been given over the corresponding bond (in angstroms). Orange, grey, white and green spheres indicate the locations of iron, carbon, hydrogen and muonium (banquet) atoms, respectively. 1 Supporting Information
Effective electronic-only Kohn-Sham equations for the muonic molecules
Milad Rayka , Mohammad Goli and Shant Shahbazian Department of Physics and Department of Physical and Computational Chemistry, Shahid Beheshti University, G. C., Evin, Tehran, Iran, 19839, P.O. Box 19395-4716. School of Nano Science, Institute for Research in Fundamental Sciences (IPM), Tehran 19395-5531, Iran
E-mails: Mohammad Goli : [email protected] Shant Shahbazian: [email protected] * Corresponding authors 2
Table of contents
Pages 3 - 8 : The Cartesian coordinates of the optimized muoniated structures depicted in Figure 2 in the main text.
Pages 9 - 14 : The Cartesian coordinates of the optimized muoniated ferrocenyl radicals depicted in Figure 3 in the main text.: The Cartesian coordinates of the optimized muoniated ferrocenyl radicals depicted in Figure 3 in the main text.