Volume 133, Issue 16, 28 October 2010

We present a new method for introducing stable nonequilibrium velocity and temperature gradients in molecular dynamics simulations of heterogeneous systems. This method extends earlier reverse nonequilibrium molecular dynamics (RNEMD) methods which use momentum exchange swapping moves. The standard swapping moves can create nonthermal velocity distributions and are difficult to use for interfacial calculations. By using nonisotropic velocity scaling (NIVS) on the molecules in specific regions of a system, it is possible to impose momentum or thermal flux between regions of a simulation while conserving the linear momentum and total energy of the system. To test the method, we have computed the thermal conductivity of model liquid and solid systems as well as the interfacial thermal conductivity of a metalwater interface. We find that the NIVSRNEMD improves the problematic velocity distributions that develop in other RNEMD methods.
 ARTICLES

 Theoretical Methods and Algorithms

A gentler approach to RNEMD: Nonisotropic velocity scaling for computing thermal conductivity and shear viscosity
View Description Hide DescriptionWe present a new method for introducing stable nonequilibrium velocity and temperature gradients in molecular dynamics simulations of heterogeneous systems. This method extends earlier reverse nonequilibrium molecular dynamics (RNEMD) methods which use momentum exchange swapping moves. The standard swapping moves can create nonthermal velocity distributions and are difficult to use for interfacial calculations. By using nonisotropic velocity scaling (NIVS) on the molecules in specific regions of a system, it is possible to impose momentum or thermal flux between regions of a simulation while conserving the linear momentum and total energy of the system. To test the method, we have computed the thermal conductivity of model liquid and solid systems as well as the interfacial thermal conductivity of a metalwater interface. We find that the NIVSRNEMD improves the problematic velocity distributions that develop in other RNEMD methods.

Identification of metastable states in peptide’s dynamics
View Description Hide DescriptionA recently developed spectral method for identifying metastable states in Markov chains is used to analyze the conformational dynamics of a fourresidue peptide valineprolinealanineleucine. We compare our results to empirically defined conformational states and show that the found metastable states correctly reproduce the conformational dynamics of the system.

Iterative Monte Carlo formulation of realtime correlation functions
View Description Hide DescriptionWe present an iterative Monte Carlo path integral methodology for evaluating thermally averaged realtime correlation functions. Standard path integral Monte Carlo methods are used to sample paths along the imaginary time contour. Propagation of the density matrix is performed iteratively on a grid composed of the end points of the sampled paths. Minimally oscillatory propagators are constructed using energy filtering techniques. A single propagation yields the values of the correlation function at all intermediate time points. Model calculations suggest that the method yields accurate results over several oscillation periods and the statistical error grows slowly with increasing propagation time.

Path integral based calculations of symmetrized time correlation functions. I.
View Description Hide DescriptionIn this paper, we examine how and when quantum evolution can be approximated in terms of (generalized) classical dynamics in calculations of correlation functions, with a focus on the symmetrized time correlation function introduced by Schofield. To that end, this function is expressed as a path integral in complex time and written in terms of sum and difference path variables. Taylor series expansion of the path integral’s exponent to first and second order in the difference variables leads to two original developments. The first order expansion is used to obtain a simple, path integral based, derivation of the socalled Schofield’s quantum correction factor. The second order result is employed to show how quantum mechanical delocalization manifests itself in the approximation of the correlation function and hinders, even in the semiclassical limit, the interpretation of the propagators in terms of sets of guiding classical trajectories dressed with appropriate weights.

Path integral based calculations of symmetrized time correlation functions. II
View Description Hide DescriptionSchofield’s form of quantum time correlation functions is used as the starting point to derive a computable expression for these quantities. The time composition property of the propagators in complex time is exploited to approximate Schofield’s function in terms of a sequence of short time classical propagations interspersed with path integrals that, combined, represent the thermal density of the system. The approximation amounts to linearization of the real time propagators and it becomes exact with increasing number of propagation legs. Within this scheme, the correlation function is interpreted as an expectation value over a probability density defined on the thermal and real path space and calculated by a Monte Carlo algorithm. The performance of the algorithm is tested on a set of benchmark problems. Although the numerical effort required is considerable, we show that the algorithm converges systematically to the exact answer with increasing number of iterations and that it is stable for times longer than those accessible via a brute force, path integral based, calculation of the correlation function. Scaling of the algorithm with dimensionality is also examined and, when the method is combined with commonly used filtering schemes, found to be comparable to that of alternative semiclassical methods.

Molecular dynamics at low time resolution
View Description Hide DescriptionThe internal dynamics of macromolecular systems is characterized by widely separated time scales, ranging from fraction of picoseconds to nanoseconds. In ordinary molecular dynamics simulations, the elementary time step used to integrate the equation of motion needs to be chosen much smaller of the shortest time scale in order not to cutoff physical effects. We show that in systems obeying the overdamped Langevin equation, it is possible to systematically correct for such discretization errors. This is done by analytically averaging out the fast molecular dynamics which occurs at time scales smaller than , using a renormalization group based technique. Such a procedure gives raise to a timedependent calculable correction to the diffusion coefficient. The resulting effective Langevin equation describes by construction the same longtime dynamics, but has a lower time resolution power, hence it can be integrated using larger time steps . We illustrate and validate this method by studying the diffusion of a pointparticle in a onedimensional toy model and the denaturation of a protein.

Spinstate splittings, highestoccupiedmolecularorbital and lowestunoccupiedmolecularorbital energies, and chemical hardness
View Description Hide DescriptionIt is known that the exact density functional must give groundstateenergies that are piecewise linear as a function of electron number. In this work we prove that this is also true for the lowestenergy excited states of different spin or spatial symmetry. This has three important consequences for chemical applications: the ground state of a molecule must correspond to the state with the maximum highestoccupiedmolecularorbital energy, minimum lowestunoccupiedmolecularorbital energy, and maximum chemical hardness. The beryllium, carbon, and vanadium atoms, as well as the and molecules are considered as illustrative examples. Our result also directly and rigorously connects the ionization potential and electron affinity to the stability of spin states.

The method of Gaussian weighted trajectories. V. On the 1GB procedure for polyatomic processes
View Description Hide DescriptionIn recent years, many chemical reactions have been studied by means of the quasiclassical trajectory(QCT) method within the Gaussian binning (GB) procedure. The latter consists of “quantizing” the final vibrational actions in Bohr spirit by putting strong emphasis on the trajectories reaching the products with vibrational actions close to integer values. A major drawback of this procedure is that if is the number of product vibrational modes, the amount of trajectories necessary to converge the calculations is larger than with the standard QCT method. Applying it to polyatomic processes is thus problematic. In a recent paper, however, Czakó and Bowman propose to quantize the total vibrational energy instead of the vibrational actions [G. Czakó and J. M. Bowman, J. Chem. Phys.131, 244302 (2009)], a procedure called 1GB here. The calculations are then only times more time consuming than with the standard QCT method, allowing thereby for considerable numerical saving. In this paper, we propose some theoretical arguments supporting the 1GB procedure and check its validity on model test cases as well as the prototype fouratom reaction.

Ab initio calculations of optical absorption spectra: Solution of the Bethe–Salpeter equation within density matrix perturbation theory
View Description Hide DescriptionWe describe an ab initio approach to compute the optical absorption spectra of molecules and solids, which is suitable for the study of large systems and gives access to spectra within a wide energy range. In this approach, the quantum Liouville equation is solved iteratively within first order perturbation theory, with a Hamiltonian containing a static selfenergy operator. This procedure is equivalent to solving the statically screened Bethe–Salpeter equation. Explicit calculations of single particle excited states and inversion of dielectric matrices are avoided using techniques based on density functional perturbation theory. In this way, full absorption spectra may be obtained with a computational workload comparable to ground state Hartree–Fock calculations. We present results for small molecules, for the spectra of a 1 nm Si cluster in a wide energy range (20 eV), and for a dipeptide exhibiting charge transfer excitations.

Strong electron correlation in the decomposition reaction of dioxetanone with implications for firefly bioluminescence
View Description Hide DescriptionDioxetanone, a key component of the bioluminescence of firefly luciferin, is itself a chemiluminescent molecule due to two conical intersections on its decomposition reactionsurface. While recent calculations of firefly luciferin have employed four electrons in four active orbitals [(4,4)] for the dioxetanone moiety, a study of dioxetanone [F. Liu et al., J. Am. Chem. Soc.131, 6181 (2009)] indicates that a much larger active space is required. Using a variational calculation of the twoelectron reduceddensitymatrix (2RDM) [D. A. Mazziotti, Acc. Chem. Res.39, 207 (2006)], we present the groundstatepotential energy surface as a function of active spaces from (4,4) to (20,17) to determine the number of molecular orbitals required for a correct treatment of the strong electron correlation near the conical intersections. Because the 2RDM method replaces exponentially scaling diagonalizations with polynomially scaling semidefinite optimizations, we readily computed large (18,15) and (20,17) active spaces that are inaccessible to traditional wave function methods. Convergence of the electron correlation with activespace size was measured with complementary RDMbased metrics, the von Neumann entropy of the oneelectron RDM as well as the Frobenius and infinity norms of the cumulant 2RDM. Results show that the electron correlation is not correctly described until the (14,12) active space with small variations present through the (20,17) space. Specifically, for active spaces smaller than (14,12), we demonstrate that at the first conical intersection, the electron in the orbital of the oxygenoxygen bond is substantially undercorrelated with the electron of the orbital and overcorrelated with the electron of the carbonyl oxygen's orbital. Based on these results, we estimate that in contrast to previous treatments, an accurate calculation of the strong electron correlation in firefly luciferin requires an active space of 28 electrons in 25 orbitals, beyond the capacity of traditional multireference wave function methods.

Frozen density embedding with hybrid functionals
View Description Hide DescriptionThe Kohn–Sham equations with constrained electron density are extended to hybrid exchangecorrelation (XC) functionals. We derive the frozen density embedding generalized Kohn–Sham (FDEGKS) scheme which allows to treat the nonlocal exactexchange in the subsystems. For practical calculations we propose an approximated version of the FDEGKS in which the nonadditive exchange potential is computed at a semilocal level. The proposed method is applied to compute the groundstate electronic properties of small test systems and selected DNA base pairs. The results of calculations employing the hierarchy of XC functionals BLYP/B3LYP/BHLYP and PBE/PBE0 are presented, in order to analyze the effect of nonlocal exchange contributions, and compared with reference coupledcluster singles and doubles results. We find that the use of hybrid functionals leads to a significant improvement in the description of groundstate electronic properties of the investigated systems. The semilocal version of the FDEGKS correctly reproduces the dipole and the electron density distribution of the exact GKS supramolecular system, with errors smaller than the ones obtained using conventional semilocal XC functionals.

Rangedependent adiabatic connections
View Description Hide DescriptionRecently, we have implemented a scheme for the calculation of the adiabatic connection linking the Kohn–Sham system to the physical, interacting system. This scheme uses a generalized Lieb functional, in which the electronic interaction strength is varied in a simple linear fashion, keeping the potential or the density fixed in the process. In the present work, we generalize this scheme further to accommodate arbitrary twoelectron operators, allowing the calculation of adiabatic connections following alternative paths as outlined by Yang [J. Chem. Phys.109, 10107 (1998)]. Specifically, we examine the errorfunction and Gaussianattenuated errorfunction adiabatic connections. It is shown that while the errorfunction connection displays some promising features, making it amenable to the possible development of new exchangecorrelation functionals by modeling the adiabatic connection integrand, the Gaussianattenuated errorfunction connection is less promising. We explore the highdensity and strong static correlation regimes for twoelectron systems. Implications of this work for the utility of rangeseparated schemes are discussed.

Unified interpretation of Hund’s first and second rules for and atoms
View Description Hide DescriptionA unified interpretation of Hund’s first and second rules for (C, N, O) and (Si, P, S) atoms is given by Hartree–Fock (HF) and multiconfiguration Hartree–Fock (MCHF) methods. Both methods exactly satisfy the virial theorem, in principle, which enables one to analyze individual components of the total energy , where , , and are the kinetic, the electronnucleus attraction, and the electronelectron repulsion energies, respectively. The correct interpretation for each of the two rules can only be achieved under the condition of the virial theorem by investigating how and interplay to attain the lower total potential energy . The stabilization of the more stable states for all the and atoms is ascribed to a greater that is caused by contraction of the valence orbitals accompanied with slight expansion of the core orbitals. The contraction of the valence orbitals for the two rules is a consequence of reducing the Hartree screening of the nucleus at short interelectronic distances. The reduced screening in the first rule is due to a greater amount of Fermi hole contributions in the state with the highest total spinangular momentum . The reduced screening in the second rule is due to the fact that two valence electrons are more likely to be on opposite sides of the nucleus in the state with the highest total orbitalangular momentum . For each of the two rules, the inclusion of correlation does not qualitatively change the HF interpretation, but HF overestimates the energy difference between two levels being compared. The magnitude of the correlation energy is significantly larger for the lower states than for the higher states since two valence electrons in the lower states are less likely to be on opposite sides of the nucleus. The MCHF evaluation of is in excellent agreement with experiment. The present HF and MCHF calculations demonstrate the above statements that were originally given by Katriel [Theor. Chem. Acta23, 309 (1972); 26, 163 (1972)]. We have, for the first time, analyzed the correlationinduced changes in the radial density distribution for the excited terms of the and atoms as well as for the ground term.
 Gas Phase Dynamics and Structure: Spectroscopy, Molecular Interactions, Scattering, and Photochemistry

N photoelectron angular distributions from fixedinspace molecules: Stereodynamics and symmetry considerations
View Description Hide DescriptionAngular distributions of N photoelectrons from fixedinspace molecules have been measured over the energy region of shape resonance and above. A multiplecoincidence velocitymap imaging technique for observation of molecular frame photoelectron angular distributions (MFPADs) has been extended to nonlinear molecular targets. Density functional theory calculations have also been conducted to elucidate the photoionizationdynamics and shape resonance in the N photoionization of . Results show that the N MFPADs exhibit strong shape variation as a function of both photoelectron kinetic energy and symmetries of final states, whereas asymmetry parameters of laboratory frame PADs show a local minimum around the shape resonance region and increase monotonically as the photon energy increases. Over the shape resonance, the spatial shape of the photoelectron wave function with symmetry closely resembles that of unoccupied molecular orbital of , although the MFPAD pattern for symmetry does not correspond directly to the orbital shape. At higher kinetic energy of 90 eV, MFPADs become less structured, but still show a significant dependence on the symmetry of final states.

Theoretical study of
View Description Hide DescriptionWe present the results of CCSD(T) calculations on the full set of complexes . Potential energy curves are calculated pointwise, employing the full counterpoise correction and basis sets of quadruple and quintuple quality, and then extrapolated to the complete basis set limit. Each curve has been employed to calculate rovibrational energy levels, from which spectroscopic parameters have been derived. These are compared to the available experimental data, and it is seen that there is excellent agreement with the values obtained from both Rydberg state extrapolations and highresolution laserinduced fluorescence studies. Finally, we have also used our potentials to calculate transport coefficients for moving through a bath of RG.

Complete experimental rovibrational eigenenergies of HNC up to above the ground state
View Description Hide DescriptionThe [H,C,N] system is one of the ideal candidate molecules to test new models aimed to calculate the manifold of the rotational, vibrational, and electronic states of a triatomic molecule. The isomerization reaction is one of the most important model systems for the study of unimolecular reactions. This paper reports on the experimental characterization of all 1191 eigenenergies up to relative to the ground state in the HNC part of the potential surface using high temperature hot gas emission spectroscopy. The spectroscopic constants for the first 27 vibrational states including highly excited bending vibrations up to are reported. The first 14 rotational perturbations have been identified and the perturbed eigenenergies were determined. The 3200 eigenenergies up to for the first 47 vibrational substates are included as supplement to this paper.

Ab initio adiabatic and quasidiabatic potential energy surfaces of lowest four electronic states of the system
View Description Hide DescriptionAb initio global adiabatic and quasidiabatic potential energy surfaces of lowest four electronic states of the system have been computed in the Jacobi coordinates using Dunning’s basis set at the internally contracted multireference (single and double) configuration interaction level of accuracy, which are relevant to the dynamics studies of inelastic vibrational and charge transfer processes observed in the scattering experiments. The computed equilibrium geometry parameters of the bound ion in the ground electronic state and other parameters for the transition state for the isomerization process, are in good quantitative agreement with those available from the high level ab initio calculations, thus lending credence to the accuracy of the potential energy surfaces. The nonadiabatic couplings between the electronic states have been analyzed in both the adiabatic and quasidiabatic frameworks by computing the nonadiabatic coupling matrix elements and the coupling potentials, respectively. It is inferred that the dynamics of energy transfer processes in the scattering experiments carried out in the range of 9.5–23 eV would involve all the four electronic states.

Ab initio characterization of the Mg–HF van der Waals complex
View Description Hide DescriptionThe equilibrium structure and the threedimensional potential energy surface of the Mg–HF van der Waals complex in its ground electronic state have been determined from accurate ab initio calculations using the coupledcluster method, CCSD(T), in conjunction with the basis sets of triple through quintuplezeta quality. The coreelectron correlation, highorder valenceelectron correlation, and scalar relativistic effects were investigated. The Mg–HF complex was confirmed to be linear at equilibrium, with a vibrationless dissociation energy (into Mg and HF) of . The vibrationrotation energy levels of two isotopologues, and , were predicted using the variational method. The predicted spectroscopic constants can be useful in a further analysis of highresolution vibrationrotation spectra of the Mg–HF complex.

Quasiclassical trajectory study of the postquenching dynamics of by on a global potential energy surface
View Description Hide DescriptionWe report fulldimensional, electronically adiabatic potential energy surfaces (PESs) for the ground state and excited state of . The PESs are permutationally invariant fits to roughly 23 000 electronic energies. Classical trajectory calculations of the postquenching dynamics of OH are carried out on the PES for and , at previously identified conical intersections (CoIs) [B. C. Hoffman and D. R. Yarkony, J. Chem. Phys.113, 10091 (2000)]. The initial momenta are sampled fully and partially microcanonically, corresponding to “adiabatic” and “diabatic” models of the dynamics, respectively. Branching ratios of reactive to nonreactive channels from separate , , and symmetries of CoIs are calculated, as are final rovibrational state distributions of OH and products. The rovibrational distributions of the OH and products, the D/Hatom translational energy distribution are calculated and compared to experimental ones. Agreement for these observable quantities is good. The branching between reactive and nonreactive quenching is sensitive to the momenta sampling; very good agreement with experiment is obtained using the diabatic sampling but not with the adiabatic sampling. The vibrational state distributions of and HOD (although not measured by experiment) are also presented.

Collisional quenching of by : Experimental and theoretical studies of the stateresolved product distribution and branching fraction
View Description Hide DescriptionWe report joint experimental and theoretical studies of outcomes resulting from the nonreactive quenching of electronically excited by . The experiments utilize a pumpprobe technique to detect the product state distribution under single collision conditions. The products are observed primarily in their lowest vibrational state with substantially less population in . The products are generated with a high degree of rotational excitation, peaking at with an average rotational energy of , and a strong propensity for populating the doublet component indicative of alignment of the halffilled orbital in the plane of OD rotation. Branching fraction measurements show that the nonreactive channel accounts for less than 20% of quenching outcomes. Complementary classical trajectory calculations of the postquenching dynamics are initiated from representative points along seams of conical intersections between the ground and excitedstate potentials of . Diabatic modeling of the initial momenta in the dynamical calculations captures the key experimental trends: products released primarily in their ground vibrational state with extensive rotational excitation and a branching ratio that strongly favors reactive quenching. The results are also compared with previous studies on the quenching of ; the two experimental studies show remarkably similar rotational energy distributions for the OH and radical products.