THE JOURNAL OF CHEMICAL PHYSICS 141, 014109 (2014)

Spin dynamics simulation of electron spin relaxation in Ni2+ (aq ) Jyrki Rantaharju,a) Jiˇrí Mareš,b) and Juha Vaarac) NMR Research Group, Department of Physics, University of Oulu, P.O. Box 3000, Oulu, FIN-90014, Finland

(Received 28 March 2014; accepted 13 June 2014; published online 3 July 2014) The ability to quantitatively predict and analyze the rate of electron spin relaxation of open-shell systems is important for electron paramagnetic resonance and paramagnetic nuclear magnetic resonance spectroscopies. We present a combined molecular dynamics (MD), quantum chemistry (QC), and spin dynamics simulation method for calculating such spin relaxation rates. The method is based on the sampling of a MD trajectory by QC calculations, to produce instantaneous parameters of the spin Hamiltonian used, in turn, to numerically solve the Liouville-von Neumann equation for the time evolution of the spin density matrix. We demonstrate the approach by simulating the relaxation of electron spin in an aqueous solution of Ni2+ ion. The spin-lattice (T1 ) and spin-spin (T2 ) relaxation rates are extracted directly from the simulations of the time dependence of the longitudinal and transverse magnetization, respectively. Good agreement with the available, indirectly obtained experimental data is obtained by our method. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4885050] I. INTRODUCTION

The rate of electron spin relaxation is of fundamental interest to electron paramagnetic resonance (EPR)1, 2 and paramagnetic nuclear magnetic resonance (pNMR)3, 4 spectroscopies. Solving the stochastic Liouville equation5 (SLE) has constituted a state-of-the-art method for electron spin relaxation calculations, in which the dynamics of the spatial degrees of freedom of the system is described by a diffusion operator. The Bloch-Redfield-Wangsness (BRW) relaxation theory,4, 6–8 widely used in the conventional NMR of diamagnetic substances, is rarely applicable to paramagnetic systems. The BRW theory assumes that the contributing interactions are weak as compared to the strength of the average Zeeman interaction with the external magnetic field B, and that the correlation time of the interactions is short as compared to the relaxation time. There exist a number of papers reporting models for electron spin relaxation in the literature.9–14 In this paper, we introduce a direct spin dynamics (SD) simulation method for investigating spin relaxation rates, using the electron spin of an aqueous solution of the Ni2+ ion as the example system. We simulate individual spins and average over the trajectories of sufficiently many of them, to approximate an ensemble. Each individual spin experiences a different realization of the interactions contributing to the relaxation of the ensemble. The time-dependent interactions that drive the time evolution of the electron spin, are obtained from a molecular dynamics (MD) simulation trajectory of Ni2+ (aq), with snapshots sampled using quantumchemical (QC) electronic structure calculations. We consider a fluctuating electron spin Hamiltonian that includes the timedependent g and zero-field splitting (ZFS) tensors, which parameterize the Zeeman interaction of the two unpaired electrons of Ni2+ with B and the spin-orbit-induced splitting of a) [email protected] b) [email protected] c) [email protected]

0021-9606/2014/141(1)/014109/6/$30.00

the electronic ground state, respectively. Related simulations were published earlier for Ni2 + (aq) by Odelius, Ribbing, and Kowalewski.15 In contrast to that work, our MD snapshot sampling method naturally includes all the motional modes of the characteristic hexagonal Ni(H2 O)2+ 6 moiety prevalent in this system. Furthermore, the currently available QC methodology to calculate spin Hamiltonian parameters enables reaching a very good agreement with the available experimental spin relaxation rates, which are obtained indirectly in pNMR relaxation experiments.11, 16, 17 At present, there also exist efficient codes for SD simulations, such as SpinDynamica18 and Spinach.19 An example of the application of SD simulations for parameterized, non-first-principles models used for the spin Hamiltonian and molecular dynamics, is provided by Ref. 20.

II. SPIN DYNAMICS

The two electrons of Ni2+ constitute an effective spin-1 system. In the spin Hamiltonian,2 we include the g and ZFS tensors (g and D) ˆ Hˆ (t) = μB Sˆ · g(t) · B + Sˆ · D(t) · S,

(1)

where μB is the Bohr magneton and Sˆ the effective spin operˆ can be expressed as a sum of the time-independent ator. H(t) isotropic and time-dependent anisotropic (interaction) parts Hˆ 0 = μB gB Sˆz = ¯ω0 Sˆz , ˆ Hˆ I (t) = μB Sˆ · [g(t) − g1] · B + Sˆ · D(t) · S,

(2)

respectively, where g is the isotropic part of g, which we calculate as 1/3 of the time average of the trace of g. QC calculations of MD snapshots produce a piecewise constant, time-ordered series of spin Hamiltonians Hˆ 1 , Hˆ 2 , Hˆ 3 , . . . , Hˆ l−1 , Hˆ l , with the time step τ . We transform ˆˆ ˆˆ ˆˆ ˆˆ ˆˆ the {Hˆ i } to superpropagators eL1 τ , eL2 τ , eL3 τ , . . . , eLl−1 τ , eLl τ ,

141, 014109-1

© 2014 AIP Publishing LLC

014109-2

Rantaharju, Mareš, and Vaara

J. Chem. Phys. 141, 014109 (2014)

ˆˆ are Liouvillians defined as21 where the L j ˆˆ = i [·, Hˆ ]. L j j ¯ The density operator for pure state is defined as ρ(t) ˆ = |(t)(t)|,

(3)

(4)

where |(t) is the state of the spin system. ρˆ is the solution of the Liouville-von Neumann equation,21, 22 d ρ(t) ˆ ˆˆ ρ(t), = L(t) ˆ (5) dt which means that we can approximate ρ(t) ˆ at the time instant t ˆˆ τ L ˆˆ ˆˆ τ L L n n−1 τ 1 . . . e ρ(0). ˆ Expressed in the ro= nτ as ρ(nτ ˆ )=e e ˆˆ nτ ˆ ˆˆ is r −L 0 ˆ ρ(nτ ˆ ) ≡ L(n)ρ(0), ˆ where L tating frame, ρˆ (nτ ) = e 0 the Liouvillian form of Hˆ 0 . We consider an ensemble of spin systems, with each member experiencing interactions caused by a different realization of the series of Hamiltonians, with the assumption that all the members of the ensemble are initially in the same state. Then, the ensemble average of ρˆ r can ˆˆ be calculated as ρˆ r (nτ ) = L(n) ρ(0). ˆ Assuming ergodicity and that the series {Hˆ i } is long ˆˆ is obtained as enough, the ensemble average of L(n) m ˆˆ 1 ˆˆ  Lˆˆ i+n τ Lˆˆ i+n−1 τ ˆˆ L(n) = e−L0 nτ e e . . . e Li τ , m i=1

(6)

where m is the number of the sub-series Hˆ i , . . . , Hˆ i+n that can be extracted from {Hˆ i }. m depends on both n and the length l of the series {Hˆ i }, m(l, n) = l − n + 1. For the values of n at ˆˆ which we calculate L(n), m(l, n) should be large enough to m ˆ ˆ render (1/m) i=1 H i (nτ ) a good approximation of Hˆˆ 0 . ρˆ r lives in a nine-dimensional space and its components can be expressed in the following orthonormal basis consistˆ ing of the shift operators and the z component of S: {Bˆ i } =

Sˆ− Sˆ− Sˆ− Sˆz + Sˆz Sˆ− Sˆ− Sˆz 3Sˆz Sˆz − 21ˆ , , ,√ , , √ 2 2 2 6 2 1ˆ Sˆ+ Sˆz + Sˆz Sˆ+ Sˆ+ Sˆ+ Sˆ+ , , , √ , 2 2 2 3

(7)

where the shift operators are Sˆ+ = Sˆx + iSˆy and Sˆ− = Sˆx − iSˆy . This is the built-in operator basis of the SpinDynamica software, which we use presently. ρˆ can now be expressed ˆˆ as a ket vector and L(n) as a matrix operator, in the ninedimensional space |ρ(0)) ˆ =

9 

9  ˆˆ cj |Bˆ j ) ; L(n) = Lj k (n)|Bˆ j )(Bˆ k |.

j =1

Mzr (nτ ) = L44 (n); Mzr (0)

Mxr (nτ ) L33 (n) + L88 (n) = , Mxr (0) 2 (10)

respectively. The correlation functions of the normalized components of Sˆ are defined as (Sˆμ |Sˆν (nτ )) r (nτ ) =  . Cμν (Sˆμ |Sˆμ )(Sˆν |Sˆν )

(11)

We calculate the Liouvillian form of the Hamiltonians {Hˆ i } as well as the matrix representations of the Liouvillians with the help of the routines of SpinDynamica. Furthermore, the matrix representations of the exponential propagators are obtained with the MatrixExp routine of Mathematica.23 We implemented a Mathematica routine for computing the evoluˆˆ tion of L(n) from the list of matrix representations of the exponential propagators. In Eq. (6), the sub-series starting with adjacent Hamiltonians are statistically correlated to the extent that much fewer than the maximum number of sub-series could be used in generating the SD trajectories, without any significant loss in the accuracy of the average magnetization ˆˆ decay. However, the propagation of L(n) is very fast and memory efficient so that we presently do not need to save on the number of SD trajectories simulated. The Ni(H2 O)2+ 6 entity has an isotropic rotational distribution. Also the dynamical behaviour of the g and ZFS tensors is independent of the direction of B. Thus, we gain threefold improved statistics by rotating these molecular property tensors so that effectively the magnetic field is, in turn, along each of the x, y, and z axes of the MD simulation box. The SD trajectories of the longitudinal and transverse magnetization obtained with the different choices are averaged over.

(8)

j,k=1

III. MOLECULAR DYNAMICS SIMULATION

Here, the elements ˆˆ Lj k (n) = (Bˆ j |Bˆ k (nτ )) = (Bˆ j |L(n)| Bˆ k )

180◦ with a π -pulse and the relaxation rate of the z component, Mz , is measured, as it decays back to equilibrium. In a T2 experiment, M is rotated by 90◦ with a π /2-pulse to the x axis, after which the relaxation rate of Mx is measured. In the present simulation method, the feedback from the individual members of the spin ensemble to the surrounding lattice is neglected and, as a result, all components of M decay to zero instead of the nonvanishing equilibrium value of Mz in a real, interacting spin system. At high temperatures, this model gives accurate results. The expressions for the normalized Cartesian z component of M in a T1 experiment, as well as the x component in a T2 experiment (where M is initially aligned along the x axis) are, at time nτ ,

(9)

are correlation functions of the members Bˆ i of the shift and z operator basis. The inner product (Bˆ j |Bˆ k (nτ )) is defined as † Tr[Bˆ j Bˆ k (nτ )]. In thermal equilibrium, the magnetization M is aligned with the Cartesian z axis. In a T1 experiment, M is flipped by

In order to perform a MD simulation of the Ni2+ ion solvated in water, we used high-level ab initio coupledcluster singles, doubles, and perturbative triples [CCSD(T)] calculations24 to obtain an applicable parameterization for the AMOEBA polarizable forcefield within the Tinker 5.0 molecular modeling package.25, 26 A simulation box containing 464 water molecules and one Ni2+ ion (under periodic boundary conditions) was first equilibrated by running MD

014109-3

Rantaharju, Mareš, and Vaara

in the standard laboratory conditions of constant pressure (1 atm) and temperature (300 K) with the Berendsen scaling algorithm.27 Subsequently, the production trajectory was obtained with a 0.5 fs integration time step in a microcanonical (NV E) run in a cubic box of size (24.0444 Å)3 , which is the average dimension of the box in an equilibrated constantpressure and temperature simulation. The long-range electrostatic interactions were treated by the particle-mesh Ewald algorithm,28 where the value of 0.3413 was used for the Ewald coefficient, together with splines of order 10 and a 40 × 40 × 40 charge grid. The cut-off lengths of the non-bonded interactions such as van der Waals, short-range electrostatics, and the direct-space part of the Ewald summation, were set to 11 Å. The N V E ensemble was used in order to avoid any artificial fluctuations of time-dependent properties. Analysis of the MD simulation30 resulted in very good structural parameters of the Ni2+ (aq) complex and, as compared to an earlier first-principles MD simulation,29 much superior dynamical properties of water in the Ni2+ (aq) solution. The most important parameters in the present context are the Ni–O distance of 2.05 Å as measured from the radial distribution function, water tilt angle of 36.4◦ in the Ni2+ complex, second-order rotational correlation time 2.0 ps of the water molecules in the bulk, and the translational diffusion constant of 0.25 Å2 /ps. Since all these parameters are in close agreement with experimental data,30 the MD trajectory can be deemed suitable for subsequent calculations. IV. QUANTUM-CHEMICAL CALCULATIONS

To obtain the series of Hamiltonians {Hˆ i }, instantaneous g and ZFS tensors were calculated from the MD snapshots. For this purpose, clusters containing the Ni2+ ion together with six water molecules constituting the first solvation shell were extracted by selecting the water molecules in which any of the atoms is closer than 2.5 Å from the Ni2+ ion. The gtensor calculations were carried out in the ORCA31–34 program at the unrestricted Kohn-Sham density-functional theory level using the PBE functional,35 together with the resolutionof-identity approximation36 and the def2-TZVP basis set.37 The ZFS tensors were calculated at the n-electron valence-state perturbation theory (NEVPT2) level38–40 applied on top of complete active space self-consistent field (CASSCF) calculations as implemented in ORCA, using the def2-TZVP basis set. In particular, only the spin-orbit contribution is presently included in the ZFS tensors calculated using quasi-degenerate perturbation theory. The omitted spinspin contribution is negligible based on calculations at the density-functional theory level. The active orbital space used in the CASSCF calculations consisted of the eight d-electrons of the metal ion in the five ligand-field orbitals of the 3d Ni2 + ion at an octahedral site. All 10 triplet and 15 singlet states accessible to the spin-orbit interaction operator from the triplet ground state with this active orbital space, were included in state-average fashion in the CASSCF. Various methods for calculating the ZFS tensor for Ni2+ systems were recently investigated by Kubica et al.41 The length of the production part of the MD trajectory was 750 ps and it was sampled with the time step of 48 fs to

J. Chem. Phys. 141, 014109 (2014)

obtain a piecewise constant series of l = 15 625 spin Hamiltonians {Hˆ i }. Whereas the calculation of the MD trajectory of this length and system size could be completed on a single CPU core within a few months, the QC snapshot calculations consumed altogether around 8 years of CPU time. Such resources could only be extracted from a nation-wide grid computing facility. To gain understanding of the required timescale of the simulations, the autocorrelation functions of the relevant interaction tensors were calculated in the form of scalar quantities A(0) : A(τ ), where A is either the ZFS or g-tensor. Both correlation functions show a sharp drop to less than 20% of their respective zero-offset values within the Hamiltonian sampling interval (48 fs). Thereafter a slower, approximately monoexponential decay takes place with the characteristic time constant of 4.6 ps and 4.3 ps for the ZFS and g-tensors, respectively. V. RESULTS AND DISCUSSION

In principle, it would be possible to investigate the individual effects of the different motional modes of the hexagonal Ni(H2 O)2+ 6 moiety on the relaxation. This corresponds to using quantum-chemically preparametrized property surfaces for the g and D tensors, expressed in the basis of the motional modes. This would be interesting and valuable particularly for selecting the relevant motional modes to be included in semianalytical relaxation models. In our direct simulation approach, all the modes are simultaneously included, however, and it would require adopting an entirely different computational strategy than the snapshot-based method that we have used for this work. References 42 and 43 present possible methods for such investigation carried out in terms of individual modes. Our choice of method is also related to our long-term goals of investigating more complicated processes than the present electron spin relaxation problem. An example of such applications is provided by the coupled spin dynamics of the unpaired electron and nuclear spins, resulting in 1 H pNMR relaxation in the present Ni2+ (aq) system. QC preparametrization of the hyperfine coupling terms of, say, beyond 1st solvation-shell protons, would involve very many structural coordinates, rendering the individual mode-based analysis impractical for such problems. Figure 1 shows the simulated magnetization components [Eq. (10)] and the cross-correlation function Czx (t), at B = 4.5 T. The plot contains 2000 SD simulation steps and has, hence, the length of 96 ps. Both the normalized z and x components of M are seen to decay single-exponentially according to Mzr (t) = e−R1 t ; Mzr (0)

Mxr (t) = e−R2 t , Mxr (0)

(12)

where R1 = 1/T1 and R2 = 1/T2 are the corresponding spinlattice and spin-spin relaxation rates, respectively. We extract T1, 2 by fitting Eq. (12) to the simulated M, with the FindFit routine of Mathematica. A fluctuation of the cross-correlation function Czx (t) can be observed. Also M appears slightly noisy. These findings arise due to the limited simulation

014109-4

Rantaharju, Mareš, and Vaara

J. Chem. Phys. 141, 014109 (2014)

FIG. 1. Electron spin magnetisation decay Mzr (t)/Mzr (0) (blue), Mxr (t)/Mxr (0) (red) in T1 and T2 relaxation simulations, respectively, for Ni2+ (aq) at 300 K. Also shown are the respective single-exponential fits, r (t) plotted with dashed lines, as well as the cross-correlation function Czx (brown) in a 4.5 T magnetic field along the z axis.

statistics due to the finite sampling frequency and length of the series {Hˆ i }. The T1 and T2 relaxation times extracted from this particular trajectory are 9.8 ps and 3.8 ps, respectively. Examples of simulated M at different field strengths are provided in Figures S1– S7 in the supplementary material.44 With B below 3 T there is a clearly visible fluctuation of the diagonal and off-diagonal elements of the ensemble-averaged superpropagator. The amplitude of the fluctuations becomes larger as B diminishes. The fluctuation amplitude depends on the length of the series of Hamiltonians {Hˆ i } used to propagate the SD trajectories; the amplitude increases as the length of {Hˆ i } shortens and the statistical quality of the simulations is reduced. Mz and Mx decay single-exponentially with B above 3 T. However as B increases, the simulated M curves become increasingly noisy. At fields above 12 T, the statistics are clearly not good enough for accurate simulations. We tested the dependence of the T1 and T2 relaxation times on the length l of the series of Hamiltonians {Hˆ i } in the range from 1563 to 15 625, thereby affecting the sampling statistics. According to the findings (Figures S8–S12 in the supplementary material44 ), the relaxation times are relatively stable in that range of l. Figure 2 shows the dependence of the simulated relaxation rates on B in the range from 0.001 T to 100 T. The rates follow the equation R1,2 (B) = Y1,2 +

A1,2 , U1,2 + B 2

(13)

FIG. 2. Simulated R1 (blue open circles) and R2 (red open circles) relaxation rates of electron spin in Ni2+ (aq) at 300 K. Also shown are their fits to Eq. (13) as blue and red solid lines, respectively. The data are presented as functions of the strength of the magnetic field.

FIG. 3. Simulated T1 (blue) and T2 (red) electron spin relaxation times in Ni2+ (aq) at 300 K as a function of magnetic field. The open circles are the simulation results and the solid lines are fitted. The filled circles correspond to simulations done without the fluctuating part of the g tensor. The stars represent the available experimental results: T1 = T2 = 2.9 ps (blue) at B = 0 T and temperature of 298 K (Ref. 11), T1 = T2 = 5.7 ps (blue, assuming T1 = T2 ) at 1.39 T and 243 K (Ref. 16), T1 = 3.4 ps (blue, assuming T1  T2 ), and T1 = T2 = 1.9 ps (red, assuming T1 = T2 ) at 2.11 T and 243 K (Ref. 17).

where Y1, 2 , A1, 2 , and U1, 2 are fitted constants. The functional form is that of the spectral density arising from BRW theory for single-exponential decay of the appropriate correlation functions.6 Figure 3 shows the field dependence of T1 and T2 extracted from the simulations and fits and Table I lists the corresponding numerical data. Simulations were also performed at a few field strengths without considering the fluctuating g tensor (replaced by its static average). Figure 3 shows that the ZFS interaction dominates relaxation in our range of B values. Experimental values from the literature11, 16, 17 are also displayed. These are indirect data obtained from the analysis of the pNMR relaxation rates of the 1 H and 17 O nuclei of water molecules, at different experimental conditions. Our simulation results are in a good agreement with these data; the order of magnitude of experimental T1 and T2 is in all cases reproduced. Taken the facts that our non-empirical method combines three molecular modeling paradigms, as well as that the analysis of the pNMR relaxation data involves various approximations, the level of agreement is remarkable. Besides scarce numerical data, the experimental paper by Friedman, Holz, and Hertz11 includes a figure of the electron T1 and T2 relaxation times as functions of magnetic field in the range between 0 T and 14 T, extracted from analytical nuclear relaxation models. The paper uses a statistical model for the ZFS interaction and fits the theory to data obtained from pNMR experiments and other works cited in Ref. 11. The differences between our direct, non-empirical simulation and the empirical approach used in Ref. 11 imply differences between the obtained results. The T2 results from Ref. 11 and our simulated data are in a good agreement. Below the field of 10 T, the simulated magnitude of T1 also agrees with the data by Friedman et al.11 However, at higher fields, the experimental T1 is larger than our simulated T1 .

014109-5

Rantaharju, Mareš, and Vaara

J. Chem. Phys. 141, 014109 (2014)

TABLE I. Simulated T1 and T2 electron spin relaxation times (in ps) of Ni2 + (aq) at 300 K, as well as the results of fits to a single-exponential magnetization decay, as a function of the magnetic field strength. B (T)

0

0.11

0.61

1.0

2.1

3.5

4.5

5.8

7.4

9.5

12

16

20

26

33

43

55

70

90

Simulated T1 Fitted T1

2.30 2.30

2.31 2.30

2.52 2.54

2.94 2.94

4.99 4.87

8.16 8.16

9.78 10.5

12.6 13.1

17.7 15.6

15.8 17.8

18.3 19.5

20.6 20.9

22.3 21.7

21.7 22.3

24.9 22.6

27.1 22.8

18.5 23.0

21.4 23.1

32.4 23.1

Simulated T2 Fitted T2

2.30 2.30

2.30 2.30

2.35 2.35

2.44 2.43

2.82 2.82

3.42 3.45

3.83 3.86

4.35 4.30

4.72 4.69

5.09 5.02

5.16 5.26

5.51 5.46

5.50 5.57

5.66 5.65

5.77 5.69

5.80 5.72

5.64 5.74

5.74 5.75

5.70 5.76

We tested the dependence of the simulation results on the chosen time step used for sampling the Hamiltonians, which is also the time step of the SD simulation. The test was carried out with the spacing between the adjacent Hamiltonians chosen in the range from 48 to 480 fs. The results presented in Figures S13–S17 of the supplementary material44 indicate that especially at high magnetic fields, the dependence on the time step is relatively strong. The magnitude of T2 results is comparatively stable with respect to the sampling interval, whereas T1 increases towards more denser sampling. Particularly at low fields, the order of magnitude of the indirect experimental results is safely reproduced, however. The presently chosen sampling interval of 48 fs represents a compromise, taking into account the computational cost of the QC snapshot calculations. VI. CONCLUSIONS

Accurate modeling of spin dynamics and relaxation is fundamentally important for magnetic resonance and its applications. We have developed a novel, coupled MD and SD method to simulate spin relaxation processes in molecular systems. The method is applied for the external magneticfield dependence of the electronic spin-spin and spin-lattice relaxation times of Ni2+ (aq). The SD calculations are timeand memory-efficient and the results are in a good agreement with the available experimental data. The method and its implementation are also applicable to more complex problems involving a combination of spin and molecular dynamics. In ongoing work, we use the method in first-principles investigations of the paramagnetic relaxation enhancement in pNMR spectroscopy. SD-based simulations of electronic and nuclear relaxation and spin exchange processes can now be performed for, in principle, any system for which the well-established MD and QC methods are feasible. Particularly, the computational cost of the QC snapshot calculations sets a practical upper limit for the sampling frequency of the set of Hamiltonians used to drive the propagation of the spin density matrix. This compromises the accuracy of the present simulation results slightly. ACKNOWLEDGMENTS

We are grateful to Malcolm H. Levitt, Pär Håkansson (Southampton), Michael C. D. Tayler (Radboud), and Jozef Kowalewski (Stockholm) for useful discussions. Financial support has been obtained from the Magnus Ehrnrooth Foundation (J.R.), the Exactus doctoral program of the University of Oulu graduate school (J.R.), the Academy of Finland, the

European Union Seventh Framework Programme (FP7/20072013) under Grant Agreements Nos. 254552 (J.M.) and 317127 (J.R., J.M., and J.V.), University of Oulu (J.V.), and the Tauno Tönning Foundation (J.V.). Computational resources due to CSC (Espoo, Finland) and the Finnish Grid Initiative project were used. 1 A.

Schweiger and G. Jeschke, Principles of Pulse Electron Paramagnetic Resonance (Oxford University Press, Oxford, 2001). 2 J. E. Harriman, Theoretical Foundations of Electron Spin Resonance (Academic Press, New York, 1978). 3 I. Bertini, C. Luchinat, and G. Parigi, Solution NMR of Paramagnetic Molecules: Applications to Metallobiomolecules and Models (Elsevier Science B.V., Amsterdam, 2001). 4 J. Kowalewski, D. Kruk, and G. Parigi, Adv. Inorg. Chem. 57, 41 (2005). 5 G. Moro and J. H. Freed, J. Phys. Chem. 84, 2837 (1980). 6 J. Kowalewski and L. Mäler, Nuclear Spin Relaxation in Liquids: Theory, Experiments, and Applications (Taylor and Francis, London, 2006). 7 A. G. Redfield, IBM J. Res. Dev. 1, 19 (1957). 8 R. K. Wangsness and F. Bloch, Phys. Rev. 89, 728 (1953). 9 D. T. Pegg and D. M. Doddrell, Austral. J. Chem. 29, 1869 (1976). 10 N. Benetis, J. Kowalewski, L. Nordenskiöld, H. Wennerström, and P.-O. Westlund, Mol. Phys. 48, 329 (1983). 11 H. L. Friedman, M. Holz, and H. G. Hertz, J. Chem. Phys. 70, 3369 (1979). 12 S. Rast, P. H. Fries, E. Belorizky, A. Borel, L. Helm, and A. E. Merbach, J. Chem. Phys. 115, 7554 (2001). 13 K. Åman and P.-O. Westlund, Mol. Phys. 102, 1085 (2004). 14 P.-O. Westlund and H. Wennerström, Phys. Chem. Chem. Phys. 12, 201 (2010). 15 M. Odelius, C. Ribbing, and J. Kowalewski, J. Chem. Phys. 104, 3181 (1996). 16 A. M. Chmelnick and D. Fiat, J. Am. Chem. Soc. 93, 2875 (1971). 17 J. Granot, A. M. Achlama, and D. Fiat, J. Chem. Phys. 61, 3043 (1974). 18 M. H. Levitt, SpinDynamica, with contributions by J. Rantaharju, A. Brinkmann, and S. S. Roy (see www.SpinDynamica.soton.ac.uk). 19 H. J. Hogben, M. Krzystyniak, G. T. P. Charnock, P. J. Hore, and I. Kuprov, J. Magn. Reson. 208, 179 (2011). 20 N. Schaefle and R. Sharp, J. Chem. Phys. 121, 5387 (2004). 21 J. Jeener, Adv. Magn. Reson. 10, 1 (1982). 22 J. J. Sakurai, in Modern Quantum Mechanics, 2nd ed., edited by J. Napolitano (Addison-Wesley, 2010). 23 Mathematica, version 8.0.4 (Wolfram Research, Inc., Champaign, IL, 2011). 24 Forcefield parameters available on request from the authors. 25 P. Ren and J. W. Ponder, J. Comput. Chem. 23, 1497 (2002). 26 P. Ren and J. W. Ponder, J. Phys. Chem. B 107, 5933 (2003). 27 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984). 28 T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993). 29 J. Mareš, H. Liimatainen, T. O. Pennanen, and J. Vaara, J. Chem. Theory Comput. 7, 3248 (2011). 30 J. Mareš and J. Vaara, “Solvation Structure and Dynamics of Ni2+ (aq) from a polarizable forcefield,” Chem. Phys. (submitted). 31 F. Neese, J. Chem. Phys. 122, 034107 (2005). 32 F. Neese, J. Chem. Phys. 127, 164112 (2007). 33 F. Neese, ORCA, version 2.9, 2012. 34 F. Neese, J. Chem. Phys. 118, 3939 (2003). 35 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996); 78, 1396(E) (1997).

014109-6 36 F.

Rantaharju, Mareš, and Vaara

Neese, J. Comput. Chem. 24, 1740 (2003). Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005). 38 C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger, and J.-P. Malrieu, J. Chem. Phys. 114, 10252 (2001). 39 C. Angeli, R. Cimiraglia, and J.-P. Malrieu, Chem. Phys. Lett. 350, 297 (2001). 40 C. Angeli, R. Cimiraglia, and J.-P. Malrieu, J. Chem. Phys. 117, 9138 (2002). 41 A. Kubica, J. Kowalewski, D. Kruk, and M. Odelius, J. Chem. Phys. 138, 064304 (2013). 37 F.

J. Chem. Phys. 141, 014109 (2014) 42 M.

Odelius, C. Ribbing, and J. Kowalewski, J. Chem. Phys. 103, 1800 (1995). 43 D. Kruk, J. Kowalewski, and P.-O. Westlund, J. Chem. Phys. 121, 2215 (2004). 44 See supplementary material at http://dx.doi.org/10.1063/1.4885050 for illustrations of M trajectories simulated at different external field strengths, the dependence of the calculated T1 and T2 relaxation times on the length of the used series of Hamiltonians {Hi }, and the dependence of the calculated T1 and T2 relaxation times on the sampling interval of {Hi }.

The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp

Spin dynamics simulation of electron spin relaxation in Ni²⁺(aq).

The ability to quantitatively predict and analyze the rate of electron spin relaxation of open-shell systems is important for electron paramagnetic re...
369KB Sizes 0 Downloads 3 Views