Comparison of different methods for calculating the paramagnetic relaxation enhancement of nuclear spins as a function of the magnetic field

The enhancement of the spin-lattice relaxation rate for nuclear spins in a ligand bound to a paramagnetic metal ion [known as the paramagnetic relaxation enhancement (PRE)] arises primarily through the dipole-dipole (DD) interaction between the nuclear spins and the electron spins. In solution, the DD interaction is modulated mostly by reorientation of the nuclear spin-electron spin axis and by electron spin relaxation. Calculations of the PRE are in general complicated, mainly because the electron spin interacts so strongly with the other degrees of freedom that its relaxation cannot be described by second-order perturbation theory or the Redfield theory. Three approaches to resolve this problem exist in the literature: The so-called slow-motion theory, originating from Swedish groups [Benetis et al., Mol. Phys. 48, 329 (1983); Kowalewski et al., Adv. Inorg. Chem. 57, (2005); Larsson et al., J. Chem. Phys. 101, 1116 (1994); T. Nilsson et al., J. Magn. Reson. 154, 269 (2002)] and two different methods based on simulations of the dynamics of electron spin in time domain, developed in Grenoble [Fries and Belorizky, J. Chem. Phys. 126, 204503 (2007); Rast et al., ibid. 115, 7554 (2001)] and Ann Arbor [Abernathy and Sharp, J. Chem. Phys. 106, 9032 (1997); Schaefle and Sharp, ibid. 121, 5387 (2004); Schaefle and Sharp, J. Magn. Reson. 176, 160 (2005)], respectively. In this paper, we report a numerical comparison of the three methods for a large variety of parameter sets, meant to correspond to large and small complexes of gadolinium(III) and of nickel(II). It is found that the agreement between the Swedish and the Grenoble approaches is very good for practically all parameter sets, while the predictions of the Ann Arbor model are similar in a number of the calculations but deviate significantly in others, reflecting in part differences in the treatment of electron spin relaxation. The origins of the discrepancies are discussed briefly.     

Relative energy of the high-(5T2g) and low-(1A1g) spin states of the ferrous complexes [Fe(L)(NHS4)]: CASPT2 versus density functional theory

High-level ab initio calculations using multiconfigurational perturbation theory [complete active space with second-order perturbation theory (CASPT2)] were performed on the transition energy between the lowest high-spin (corresponding to (5T2g) in Oh) and low-spin (corresponding to 1A1g in Oh) states in the series of six-coordinated Fe(II) molecules [Fe(L)(NHS4)], where NHS4 is 2,2'-bis(2-mercaptophenylthio)diethylamine dianion and L=NH3, N2H4, PMe3, CO, and NO+. The results are compared to (previous and presently obtained) results from density functional theory (DFT) calculations with four functionals, which were already shown previously by Casida and co-workers [Fouqueau et al., J. Chem. Phys. 120, 9473 (2004); Ganzenmuller et al., ibid. 122, 234321 (2005); Fouqueau et al., ibid. 122, 044110 (2005); Lawson Daku et al., ChemPhysChem 6, 1393 (2005)] to perform well for the spin-pairing problem in these and other Fe(II) complexes, i.e., OLYP, PBE0, B3LYP, and B3LYP*. Very extended basis sets were used both for the DFT and CASPT2 calculations and were shown to be necessary to obtain quantitative results with both types of method. This work presents a sequel to a previous DFT/CASPT2 study of the same property in the complexes [Fe(H2O)6]2+, [Fe(NH3)6]2+, and [Fe(bpy)3]2+ [Pierloot et al., J. Chem. Phys. 125, 124303 (2006)]. The latter work was extended with new results obtained with larger basis sets and including the OLYP functional. For all considered complexes, the CASPT2 method predicts the correct ground state spin multiplicity. Since experimental data for the actual quintet-singlet (free) energy differences are not available, the performance of the different DFT functionals was judged based on the comparison between the DFT and CASPT2 results. From this, it was concluded that the generalized gradient OLYP functional performs remarkably well for the present series of ferrous compounds, whereas the success of the three hybrid functionals varies from case to case.     

Modeling of HeN+ clusters. II. Calculation of He3+ vibrational spectrum

We have computed the vibrational spectrum of the helium ionized trimer He(3)(+) using three different potential energy surfaces [D. T. Chang and G. L. Gellene, J. Chem. Phys. 119, 4694 (2003); E. Scifoni et al., ibid. 125, 164304 (2006); I. Paidarova et al., Chem. Phys. 342, 64 (2007)]. Differences in the details of these potential energy surfaces induce discrepancies between bound state energies of the order of 0.01 eV. The effects of the geometric phase induced by the conical intersection between the ground electronic potential energy surface and the first excited one are studied by computing vibrational spectra with and without this phase. The six lowest vibrational bound states are negligibly affected by the geometric phase. Indeed, they correspond to wavefunctions localized in the vicinity of the linear symmetric configurations and can be assigned well defined vibrational quantum numbers. On the other hand, higher excited states are delocalized, cannot be assigned definite vibrational quantum numbers, and the geometric phase shifts their energies by approximately 0.005 eV.     

Hybrid one-electron/many-electron methods for ionized states of molecular clusters

To describe singly-ionized states of molecular clusters we devised an effective Hamiltonian approach that combines (1) accurate monomer ionization potentials from many-electron wave functions with (2) polarization shifts and (3) effective monomer couplings obtained from a simple one-electron approach (the superposition-of-fragment-states (SFS) method [Valeev et al., J. Am. Chem. Soc., 2006, 128, 9882]). The accuracy of the intermolecular coupling parameters evaluated with SFS Hartree-Fock (HF) and Density-Functional-Theory (DFT) variants was evaluated for several weakly-bound dimers and compared against the state-of-the-art equation-of-motion ionization-potential coupled-cluster singles and doubles (EOM-IP-CCSD) data of Krylov and co-workers. The SFS-HF method produces coupling integrals accurate to a few percent, whereas SFS-DFT predictions are substantially worse. A hybrid approach combining SFS-HF couplings and shifts with EOM-IP-CCSD ionization potentials of monomers (denoted as SFS-EOM-IP-CCSD) was applied to ionized states of two conformers of a benzene dimer and ten representative DNA base pairs. The 16 considered SFS-EOM-IP-CCSD ionization potentials of the benzene dimer differed from the reference EOM-IP-CCSD IPs of Krylov and co-workers [Pieniazek et al., J. Chem. Phys. 2007, 127, 044317; Bravaya et al., Phys. Chem. Chem. Phys. 2010, 12, 2261] by less than 0.1 eV on average, and at most by 0.2 eV. For the DNA base pairs the mean absolute (median) deviation of the SFS-EOM-IP-CCSD IPs was 0.27 (0.23) eV; several deviations for non-Koopmans states were as large as 0.9 eV. The SFS-EOM-IP-CCSD method can be readily applied to large molecular clusters with computational effort scaling cubically with the size of the cluster.     

Molecular simulation of polymer assisted protein refolding

Protein refolding in vitro, the formation of the tertiary structure that enables the protein to display its biological function, can be significantly enhanced by adding a polymer of an appropriate hydrophobicity and concentration into the refolding buffer. A molecular simulation of the refolding of a two-dimensional simple lattice protein was presented. A protein folding map recording the occurrence frequency of specified conformations was derived, from which the refolding thermodynamics and kinetics were interpreted. It is shown that, in the absence of polymer, the protein falls into the "energy trapped" conformations characterized by a high intramolecular hydrophobic interaction, denoted as HH contact, and a high magnitude of the structure overlap function, chi. This makes it difficult for the protein to fold to the native state. The polymer with a suitable chain length, concentration, and hydrophobicity has formed complex with partially folded protein and created diversified intermediates with low chi. This gives more pathways for the protein to fold to the native state. At a given hydrophobicity, the short chain polymer has a broader concentration range where it assists protein folding than those of long chains. The above simulation agrees well with the experimental results reported elsewhere [Cleland et al., J. Biol. Chem. 267, 13327 (1992); ibid., Bio/Technology 10, 1013 (1992); Chen et al., Enzyme Microb. Technol. 32, 120 (2003); Lu et al., Biochem. Eng. J. 24, 55 (2005); ibid., J. Chem. Phys. 122, 134902 (2005); ibid., Biochem. Eng. J. (to be published)] and is of fundamental importance for the design and application of polymers for protein refolding.     

A global 12-dimensional ab initio potential energy surface and dynamical studies for the SiH4+H-->SiH3+H2 reaction

A global 12-dimensional ab initio interpolated potential energy surface (PES) for the SiH(4)+H-->SiH(3)+H(2) reaction is presented. The ab initio calculations are based on the unrestricted quadratic configuration interaction treatment with all single and double excitations together with the cc-pVTZ basis set, and the modified Shepard interpolation method of Collins and co-workers [K. C. Thompson et al., J. Chem. Phys. 108, 8302 (1998); M. A. Collins, Theor. Chem. Acc. 108, 313 (2002); R. P. A. Bettens and M. A. Collins, J. Chem. Phys. 111, 816 (1999)] is applied. Using this PES, classical trajectory and variational transition state theory calculations have been carried out, and the computed rate constants are in good agreement with the available experimental data.     

On the stability of glycine-water clusters with excess electron: implications for photoelectron spectroscopy

Calculations are presented for the glycine-(H(2)O)(n) (-) (n=0-2) anionic clusters with excess electron, with the glycine core in the canonical or zwitterion form. A variety of conformers are predicted, and their relative energy is examined to estimate thermodynamic stability. The dynamic (proton transfer) pathways between the anionic clusters with the canonical and the zwitterion glycine core are examined. Small barrier heights for isomerization from the zwitterion glycine-(H(2)O)(2) (-) anion to those with canonical glycine core suggest that the former conformers may be kinetically unstable and unfavorable for detection of neutral glycine zwitterion-(H(2)O)(n) (n=1,2) clusters by photodetachment, in accordance with the photoelectron spectroscopic experiments by Bowen and co-workers [Xu et al., J. Chem. Phys. 119, 10696 (2003)]. The calculated stability of the glycine-(H(2)O)(n) (-) anion clusters with canonical glycine core relative to those with zwitterion core indicates that the observation of the anionic conformers with the canonical glycine core would be much more feasible, as revealed by Johnson and co-workers [Diken et al. J. Chem. Phys. 120, 9902 (2004)].     

Accuracy of recent potential energy surfaces for the He-N2 interaction. II. Molecular beam scattering and bulk gas relaxation phenomena

A new semiempirical exchange-Coulomb model potential energy surface for the N(2)-He interaction was reported recently [A. K. Dham et al., J. Chem. Phys. 127, 054302 (2007)] and, using it, the temperature dependence of bulk gas properties of N(2)-He mixtures, such as the second virial coefficient and traditional transport phenomena, most of which depend primarily on the isotropic component of the interaction potential energy surface, was determined. Values of these properties, along with values calculated using two high-quality ab initio potential energy surfaces [C.-H. Hu and A. J. Thakkar, J. Chem. Phys. 104, 2541 (1996); K. Patel et al., ibid 119, 909 (2003)] were compared critically to available experimental data. The present paper reports on the ability of the same three potential energy surfaces to predict state-to-state and total differential cross sections, total integral cross sections, and the temperature dependence of bulk gas relaxation phenomena (including magnetic field effects on transport coefficients). While all three potential energy surfaces give total differential and higher speed integral scattering results that fall within the experimental uncertainties, integral scattering results and state-to-state differential cross section measurements consistently exceed the calculated values. All three surfaces give similar agreement with the relaxation properties of N(2)-He binary mixtures, with the semiempirical exchange-Coulomb model potential energy surface giving slightly better overall agreement with experiment than the two ab initio potential energy surfaces.     

Electronic structure and spectrum of UO2 2+ and UO2Cl4 2-

A theoretical study is presented of the electronic spectra of the UO(2) (2+) and UO(2)Cl(4) (2-) ions, based on multiconfigurational perturbation theory (CASSCF/CASPT2), combined with a recently developed method to treat spin-orbit coupling [P.-A. Malmqvist et al., Chem. Phys. Lett. 357, 230 (2002); B. O. Roos and P.-A. Malmqvist, Phys. Chem. Chem. Phys. 6, 2919 (2004)]. The results are compared to the experimental spectroscopic data obtained for uranyl ions in Cs(2)UO(2)Cl(4) crystals from Denning [Struct. Bonding (Berlin) 79, 215 (1992)] and to previous theoretical calculations performed using a combined configuration-interaction spin-orbit treatment [Z. Zhang and R. M. Pitzer, J. Phys. Chem. A 103, 6880 (1999); S. Matsika and R. M. Pitzer, J. Phys. Chem. A. 105, 637 (2001)]. As opposed to the latter results, the calculations performed in this work point to a significant effect of the weakly bound equatorial chlorine ligands on the excitation energies.     

Adding salt to an aqueous solution of t-butanol: is hydrophobic association enhanced or reduced?

Recent neutron scattering experiments on aqueous salt solutions of amphiphilic t-butanol by Bowron and Finney [Phys. Rev. Lett. 89, 215508 (2002); J. Chem. Phys. 118, 8357 (2003)] suggest the formation of t-butanol pairs, bridged by a chloride ion via O-H...Cl- hydrogen bonds, leading to a reduced number of intermolecular hydrophobic butanol-butanol contacts. Here we present a joint experimental/theoretical study on the same system, using a combination of molecular dynamics (MD) simulations and nuclear magnetic relaxation measurements. Both MD simulation and experiment clearly support the more classical scenario of an enhanced number of hydrophobic contacts in the presence of salt, as it would be expected for purely hydrophobic solutes. [T. Ghosh et al., J. Phys. Chem. B 107, 612 (2003)]. Although our conclusions arrive at a structurally completely distinct scenario, the molecular dynamics simulation results are within the experimental error bars of the Bowron and Finney data.     

The new insight into dynamic crossover in glass forming liquids from the apparent enthalpy analysis

One of the most intriguing phenomena in glass forming systems is the dynamic crossover (T(B)), occurring well above the glass temperature (T(g)). So far, it was estimated mainly from the linearized derivative analysis of the primary relaxation time τ(T) or viscosity η(T) experimental data, originally proposed by Stickel et al. [J. Chem. Phys. 104, 2043 (1996); J. Chem. Phys. 107, 1086 (1997)]. However, this formal procedure is based on the general validity of the Vogel-Fulcher-Tammann equation, which has been strongly questioned recently [T. Hecksher et al. Nature Phys. 4, 737 (2008); P. Lunkenheimer et al. Phys. Rev. E 81, 051504 (2010); J. C. Martinez-Garcia et al. J. Chem. Phys. 134, 024512 (2011)]. We present a qualitatively new way to identify the dynamic crossover based on the apparent enthalpy space (H(a) (')=dlnτ/d(1/T)) analysis via a new plot lnH(a) (') vs. 1∕T supported by the Savitzky-Golay filtering procedure for getting an insight into the noise-distorted high order derivatives. It is shown that depending on the ratio between the "virtual" fragility in the high temperature dynamic domain (m(high)) and the "real" fragility at T(g) (the low temperature dynamic domain, m = m(low)) glass formers can be splitted into two groups related to f < 1 and f > 1, (f = m(high)∕m(low)). The link of this phenomenon to the ratio between the apparent enthalpy and activation energy as well as the behavior of the configurational entropy is indicated.     

Pair-correlated speed distributions for the OH+CH4/CD4 reactions: Further remarks on their classical trajectory calculations in a quantum spirit

Ten years ago, Liu and co-workers measured pair-correlated speed distributions for OH+CH₄/CD₄ reactions by means of velocity map imaging (VMI) techniques at a collision energy of ∼10 kcal/mol [B. Zhang, W. Shiu, J. J. Lin and K. Liu, J. Chem. Phys 122, 131102 (2005); B. Zhang, W. Shiu and K. Liu, J. Phys. Chem. A 109, 8989 (2005)]. Recently, two of us could semi-quantitatively reproduce these measurements by performing full-dimensional quasi-classical trajectory calculations in a quantum spirit on an ab-initio potential energy surface of their own [J. Espinosa-Garcia and J. C. Corchado, Theor. Chem. Acc. (2015) 134: 6; J. Phys. Chem. B 120, 1446 (2016)]. The goal of the present work is to show that these results can be significantly improved by adding a few more constraints in order to better comply with the restrictions imposed by VMI. Overall, the level of agreement between theory and experiment is remarkable owing to the large dimensionality of the reactions under scrutiny. This is an encouraging result considering the computational challenges of quantum scattering calculations for such large processes.     

Active-space coupled-cluster methods through connected quadruple excitations

Coupled-cluster methods that include just a subset of all connected triple, quadruple, or both excitation amplitudes, according to the ansatz of and Adamowicz co-workers [Int. Rev. Phys. Chem. 12, 339 (1993); J. Chem. Phys. 99, 1875 (1993); 100, 5792 (1994)] and Piecuch et al. [J. Chem. Phys. 110, 6103 (1999)], have been implemented into parallel execution programs. They are applicable to closed- and open-shell species and they take advantage of real Abelian point-group symmetry. A symbol manipulation program has been invoked to automate the implementation. These methods have been applied to the singlet-triplet separations of five triatomic hydrides (CH2, NH2+, SiH2, PH2+, and AsH2+) with consideration of scalar relativistic effects. They have been shown to be remarkably effective with errors arising from the use of a very small subset of higher-order excitations being no more than a few tenths of 1 kcal/mol.     

Modeling collective behavior of molecules in nanoscale direct deposition processes

We present a theoretical model describing the collective behavior of molecules in nanoscale direct deposition processes such as dip-pen nanolithography. We show that strong intermolecular interactions combined with nonuniform substrate-molecule interactions can produce various shapes of molecular patterns including fractal-like structures. Computer simulations reveal circular and starlike patterns at low and intermediate densities of preferentially attractive surface sites, respectively. At large density of such surface sites, the molecules form a two-dimensional invasion percolation cluster. Previous experimental results showing anisotropic patterns of various chemical and biological molecules correspond to the starlike regime [P. Manandhar et al., Phys. Rev. Lett. 90, 115505 (2003); J.-H. Lim and C. A. Mirkin, Adv. Mater. (Weinheim, Ger.) 14, 1474 (2002); D. L. Wilson et al., Proc. Natl. Acad. Sci. U.S.A. 98, 13660 (2001); M. Su et al., Appl. Phys. Lett. 84, 4200 (2004); R. McKendry et al., Nano Lett. 2, 713 (2002); H. Zhou et al., Appl. Surf. Sci. 236, 18 (2004); G. Agarwal et al., J. Am. Chem. Soc. 125, 580 (2003)].     

New theoretical considerations in polymer rheology: elastic breakdown of chain entanglement network

Recent experimental evidence has motivated us to present a set of new theoretical considerations and to provide a rationale for interpreting the intriguing flow phenomena observed in entangled polymer solutions and melts [P. Tapadia and S. Q. Wang, Phys. Rev. Lett. 96, 016001 (2006); 96, 196001 (2006); S. Q. Wang et al., ibid. 97, 187801 (2006)]. Three forces have been recognized to play important roles in controlling the response of a strained entanglement network. During flow, an intermolecular locking force f(iml) arises and causes conformational deformation in each load-bearing strand between entanglements. The chain deformation builds up a retractive force f(retract) within each strand. Chain entanglement prevails in quiescence because a given chain prefers to stay interpenetrating into other chains within its pervaded volume so as to enjoy maximum conformational entropy. Since each strand of length l(ent) has entropy equal to k(B)T, the disentanglement criterion is given by f(retract)>f(ent) approximately k(B)Tl(ent) in the case of interrupted deformation. This condition identifies f(ent) as a cohesive force. Imbalance among these forces causes elastic breakdown of the entanglement network. For example, an entangled polymer yields during continuous deformation when the declining f(iml) cannot sustain the elevated f(retract). This opposite trend of the two forces is at the core of the physics governing a "cohesive" breakdown at the yield point (i.e., the stress overshoot) in startup flow. Identifying the yield point as the point of force imbalance, we can also rationalize the recently observed striking scaling behavior associated with the yield point in continuous deformation of both shear and extension.     

Vibrational calculation for the HOCO radical and the cis-HOCO anion

We present numerically exact vibrational transitions for trans-HOCO, cis-HOCO, and cis-HOCO(-) for the quartic force fields of Fortenberry et al. [J. Chem. Phys. 135, 134301 (2011); ibid. 135, 214303 (2011)], obtained by means of a computational strategy based on the discrete variable representation. Several adiabatic projection schemes have been employed to characterize the vibrational levels and to study the relevance of the intermode coupling (vibrational mixing). Our results help to clear up a large discrepancy between previously reported vibrational perturbation theory and vibrational configuration interaction predictions for the torsional frequency.     

Effect of the overall rotation on the cis-trans isomerization of HONO induced by an external field

Rovibrational eigenenergies of HONO are computed and compared to experimental energies available in the literature. For their computation, we use a previously developed potential energy surface (PES) and a newly derived exact kinetic energy operator (KEO) including the overall rotation for a tetra-atomic molecule in non-orthogonal coordinates. In addition, we use the Heidelberg Multi-Configuration Time-Dependent Hartree (MCTDH) package. We compare the experimental rovibrational eigenvalues of HONO available in the literature with those obtained with MCTDH and a previously developed potential energy surface (PES) [F. Richter et al., J. Chem. Phys., 2004, 120, 1306.] for the cis geometry. The effect of the overall rotation on the process studied in our previous work on HONO [F. Richter et al., J. Chem. Phys., 2007, 127, 164315.] leading to the cis→trans isomerization of HONO is investigated. This effect on this process is found to be weak.     

Speed dependent rotational angular momentum polarization of the O2 (a 1Deltag) fragment following ozone photolysis in the wavelength range 248-265 nm

The translational anisotropy and rotational angular momentum polarization of a selection of rotational states of the O2 (a 1Deltag; v=0) photofragment formed from ozone photolysis at 248, 260, and 265 nm have been determined using the technique of resonance enhanced multiphoton ionization in combination with time of flight mass spectrometry. At 248 nm, the dissociation is well described as impulsive in nature with all rotational states exhibiting similarly large, near-limiting values for the bipolar moments describing their angular momentum alignment and orientation. At 265 nm, however, the angular momentum polarization parameters determined for consecutive odd and even rotational states exhibit clear differences. Studies at the intermediate wavelength of 260 nm strongly suggest that such a difference in the angular momentum polarization is speed dependent and this proposal is consistent with the angular momentum polarization parameters extracted and reported previously for longer photolysis wavelengths [G. Hancock et al., Phys. Chem. Chem. Phys. 5, 5386 (2003); S. J. Horrocks et al., J. Chem. Phys. 126, 044308 (2007)]. The alternation of angular momentum polarization for successive odd and even J states may be a consequence of the different mechanisms leading to the formation of the two O2 (a 1Deltag) Lambda doublets. Specifically, the involvement of out of plane parent rotational motion is proposed as the origin for the observed depolarization for the Delta- relative to the Delta+ state.     

Local structure of Pt and Pd ions in Ce(1-x)TixO2: x-ray diffraction, x-ray photoelectron spectroscopy, and extended x-ray absorption fine structure

Ce(1-x-y)Ti(x)Pt(y)O(2-delta) (x=0.15; y=0.01) and Ce(1-x-y)Ti(x)Pd(y)O(2-delta) (x=0.25; y=0.02 and 0.05) are found to be good CO oxidation catalysts [T. Baidya et al., J. Phys. Chem. B 110, 5262 (2006); T. Baidya et al., J. Phys. Chem. C 111, 830 (2007)]. A detailed structural study of these compounds has been carried out by extended x-ray absorption fine structure along with x-ray diffraction and x-ray photoelectron spectroscopy. The gross cubic fluorite structure of CeO(2) is retained in the mixed oxides. Oxide ion sublattice around Ti as well as Pt and Pd ions is destabilized in the solid solution. Instead of ideal eight coordinations, Ti, Pd, and Pt ions have 4+3, 4+3, and 3+4 coordinations creating long and short bonds. The long Ti-O, Pd-O, and Pt-O bonds are approximately 2.47 A (2.63 A for Pt-O) which are much higher than average Ce-O bonds of 2.34 A.     

The role of the excited electronic states in the C+ + H2O reaction

The electronic excited states of the [COH2]+ system have been studied in order to establish their role in the dynamics of the C+ + H2O-->[COH]+ +H reaction, which is a prototypical ion-molecule reaction. The most relevant minima and saddle points of the lowest excited state have been determined and energy profiles for the lowest excited doublet and quartet electronic states have been computed along the fragmentation and isomerization coordinates. Also, nonadiabatic coupling strengths between the ground and the first excited state have been computed where they can be large. Our analysis suggests that the first excited state could play an important role in the generation of the formyl isomer, which has been detected in crossed beam experiments [D. M. Sonnenfroh et al., J. Chem. Phys. 83, 3985 (1985)], but could not be explained in quasiclassical trajectory computations [Y. Ishikawa et al., Chem. Phys. Lett. 370, 490 (2003); J. R. Flores, J. Chem. Phys. 125, 164309 (2006)].