Volume 144, Issue 5, 07 February 2016

We present a new algorithm to solve the polarizable continuum modelequation in a framework compatible with the strategy previously developed by us for the conductorlike screening model based on Schwarz’s domain decomposition method (ddCOSMO). The new discretization is systematically improvable and is fully consistent with ddCOSMO so that it reproduces ddCOSMO results for large dielectric constants.
 COMMUNICATIONS


Communication: Librational dynamics in water, sI and sII clathrate hydrates, and ice Ih: Moleculardynamics insights
View Description Hide DescriptionEquilibrium moleculardynamics simulations have been performed for liquid water, and on metastable sI and sII polymorphs of empty hydrate lattices, in addition to ice Ih, in order to study the dynamical properties of librational motion (rotation oscillation) depicted by protons in water molecules. In particular, hydrate lattices were found to display prominent “bifurcated” features, or peaks, at circa 70 and 8095 meV (or ∼560 and 640760 cm^{−1}, respectively), also displayed by ice, in essentially quantitative agreement with experimental neutronscattering data. However, observed differences in dispersion between these librational modes between these two structures (both hydrate polymorphs visàvis ice), owing primarily to density effects, have been decomposed into contributions arising from angularvelocity dynamics about axes in the local molecular frame of water molecules, with inplane “wagging” and “twisting” rationalising one mode at ∼70 meV, and outofplane motion for the higherfrequency band. This was confirmed explicitly by a type of de facto normalmode analysis, in which only immediate layers of water molecules about the one under consideration were allowed to move. In contrast, liquid water displayed no marked preference for such local in or outofplane modes characterising librational motion, owing to the marked absence of rigid, pentamers or hexamers therein.

Communication: Consistent interpretation of molecular simulation kinetics using Markov state models biased with external information
View Description Hide DescriptionMolecular simulations can provide microscopic insight into the physical and chemical driving forces of complex molecular processes. Despite continued advancement of simulation methodology, model errors may lead to inconsistencies between simulated and reference (e.g., from experiments or higherlevel simulations) observables. To bound the microscopic information generated by computer simulations within reference measurements, we propose a method that reweights the microscopic transitions of the system to improve consistency with a set of coarse kinetic observables. The method employs the welldeveloped Markov state modeling framework to efficiently link microscopic dynamics with longtime scale constraints, thereby consistently addressing a wide range of time scales. To emphasize the robustness of the method, we consider two distinct coarsegrained models with significant kinetic inconsistencies. When applied to the simulated conformational dynamics of small peptides, the reweighting procedure systematically improves the time scale separation of the slowest processes. Additionally, constraining the forward and backward rates between metastable states leads to slight improvement of their relative stabilities and, thus, refined equilibrium properties of the resulting model. Finally, we find that difficulties in simultaneously describing both the simulated data and the provided constraints can help identify specific limitations of the underlying simulation approach.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

A new discretization for the polarizable continuum model within the domain decomposition paradigm
View Description Hide DescriptionWe present a new algorithm to solve the polarizable continuum modelequation in a framework compatible with the strategy previously developed by us for the conductorlike screening model based on Schwarz’s domain decomposition method (ddCOSMO). The new discretization is systematically improvable and is fully consistent with ddCOSMO so that it reproduces ddCOSMO results for large dielectric constants.

Efficient linearscaling secondorder MøllerPlesset perturbation theory: The divide–expand–consolidate RIMP2 model
View Description Hide DescriptionThe Resolution of the Identity secondorder MøllerPlesset perturbation theory (RIMP2) method is implemented within the linearscaling DivideExpandConsolidate (DEC) framework. In a DEC calculation, the full molecular correlated calculation is replaced by a set of independent fragment calculations each using a subset of the total orbital space. The number of independent fragment calculations scales linearly with the system size, rendering the method linearscaling and massively parallel. The DECRIMP2 method can be viewed as an approximation to the DECMP2 method where the RI approximation is utilized in each fragment calculation. The individual fragment calculations scale with the fifth power of the fragment size for both methods. However, the DECRIMP2 method has a reduced prefactor compared to DECMP2 and is wellsuited for implementation on massively parallel supercomputers, as demonstrated by test calculations on a set of mediumsized molecules. The DEC error control ensures that the standard RIMP2 energy can be obtained to the predefined precision. The errors associated with the RI and DEC approximations are compared, and it is shown that the DECRIMP2 method can be applied to systems far beyond the ones that can be treated with a conventional RIMP2 implementation.

Critical analysis of fragmentorbital DFT schemes for the calculation of electronic coupling values
View Description Hide DescriptionWe present a critical analysis of the popular fragmentorbital densityfunctional theory (FODFT) scheme for the calculation of electronic coupling values. We discuss the characteristics of different possible formulations or “flavors” of the scheme which differ by the number of electrons in the calculation of the fragments and the construction of the Hamiltonian. In addition to two previously described variants based on neutral fragments, we present a third version taking a different route to the approximate diabatic state by explicitly considering charged fragments. In applying these FODFT flavors to the two molecular test sets HAB7 (electron transfer) and HAB11 (hole transfer), we find that our new scheme gives improved electronic couplings for HAB7 (−6.2% decrease in mean relative signed error) and greatly improved electronic couplings for HAB11 (−15.3% decrease in mean relative signed error). A systematic investigation of the influence of exact exchange on the electronic coupling values shows that the use of hybrid functionals in FODFT calculations improves the electronic couplings, giving values close to or even better than more sophisticated constrained DFT calculations. Comparing the accuracy and computational cost of each variant, we devise simple rules to choose the best possible flavor depending on the task. For accuracy, our new scheme with chargedfragment calculations performs best, while numerically more efficient at reasonable accuracy is the variant with neutral fragments.

Timedependent wave packet averaged vibrational frequencies from femtosecond stimulated Raman spectra
View Description Hide DescriptionFemtosecond stimulated Raman spectroscopy (FSRS) on the Stokes side arises from a third order polarization,P^{(3)}(t), which is given by an overlap of a first order wave packet, , prepared by a narrow band (ps) Raman pump pulse, Epu(t), on the upper electronic e2potential energy surface (PES), with a second order wave packet, , that is prepared on the lower electronic e1 PES by a broadband (fs) probe pulse, Epr(t), acting on the firstorder wave packet. In offresonant FSRS, resembles the zeroth order wave packet on the lower PES spatially, but with a force on along the coordinates of the reporter modes due to displacements in the equilibrium position, so that will oscillate along those coordinates thus giving rise to similar oscillations in P^{(3)}(t) with the frequencies of the reporter modes. So, by recovering P^{(3)}(t) from the FSRS spectrum, we are able to deduce information on the timedependent quantummechanical wave packet averaged frequencies, , of the reporter modes j along the trajectory of . The observable FSRS Raman gain is related to the imaginary part of P^{(3)}(ω). The imaginary and real parts of P^{(3)}(ω) are related by the KramersKronig relation. Hence, from the FSRS Raman gain, we can obtain the complex P^{(3)}(ω), whose Fourier transform then gives us the complex P^{(3)}(t) to analyze for . We apply the theory, first, to a twodimensional model system with one conformational mode of low frequency and one reporter vibrational mode of higher frequency with good results, and then we apply it to the timeresolved FSRS spectra of the cistrans isomerization of retinal in rhodopsin [P. Kukura et al., Science 310, 1006 (2005)]. We obtain the vibrational frequency upshift time constants for the C12H wagging mode at 216 fs and for the C10H wagging mode at 161 fs which are larger than for the C11H wagging mode at 127 fs, i.e., the C11H wagging mode arrives at its final frequency while the C12H and C10H wagging modes are still upshifting to their final values, agreeing with the findings of Yan et al. [Biochemistry 43, 10867 (2004)].

Variational tensor approach for approximating the rareevent kinetics of macromolecular systems
View Description Hide DescriptionEssential information about the stationary and slow kinetic properties of macromolecules is contained in the eigenvalues and eigenfunctions of the dynamical operator of the molecular dynamics. A recent variational formulation allows to optimally approximate these eigenvalues and eigenfunctions when a basis set for the eigenfunctions is provided. In this study, we propose that a suitable choice of basis functions is given by products of onecoordinate basis functions, which describe changes along internal molecular coordinates, such as dihedral angles or distances. A sparse tensor product approach is employed in order to avoid a combinatorial explosion of products, i.e., of the basis set size. Our results suggest that the highdimensional eigenfunctions can be well approximated with relatively small basis set sizes.

Selfconsistent secondorder Green’s function perturbation theory for periodic systems
View Description Hide DescriptionDespite recent advances, systematic quantitative treatment of the electron correlation problem in extended systems remains a formidable task. Systematically improvable Green’s function methods capable of quantitatively describing weak and at least qualitatively strong correlations appear as promising candidates for computational treatment of periodic systems. We present a periodic implementation of temperaturedependent selfconsistent 2ndorder Green’s function (GF2) method, where the selfenergy is evaluated in the basis of atomic orbitals. Evaluating the realspace selfenergy in atomic orbitals and solving the Dyson equation in kspace are the key components of a computationally feasible algorithm. We apply this technique to the onedimensional hydrogen lattice — a prototypical crystalline system with a realistic Hamiltonian. By analyzing the behavior of the spectral functions, natural occupations, and selfenergies, we claim that GF2 is able to recover metallic, band insulating, and at least qualitatively Mott regimes. We observe that the iterative nature of GF2 is essential to the emergence of the metallic and Mott phases.

Phase space barriers and dividing surfaces in the absence of critical points of the potential energy: Application to roaming in ozone
View Description Hide DescriptionWe examine the phase space structures that govern reactiondynamics in the absence of critical points on the potential energy surface. We show that in the vicinity of hyperbolic invariant tori, it is possible to define phase space dividing surfaces that are analogous to the dividing surfaces governing transition from reactants to products near a critical point of the potential energy surface. We investigate the problem of capture of an atom by a diatomic molecule and show that a normally hyperbolic invariant manifold exists at large atomdiatom distances, away from any critical points on the potential. This normally hyperbolic invariant manifold is the anchor for the construction of a dividing surface in phase space, which defines the outer or loose transition state governing capture dynamics. We present an algorithm for sampling an approximate capture dividing surface, and apply our methods to the recombination of the ozone molecule. We treat both 2 and 3 degrees of freedom models with zero total angular momentum. We have located the normally hyperbolic invariant manifold from which the orbiting (outer) transition state is constructed. This forms the basis for our analysis of trajectories for ozone in general, but with particular emphasis on the roaming trajectories.

Natural occupation numbers in twoelectron quantum rings
View Description Hide DescriptionNatural orbitals (NOs) are central constituents for evaluating correlation energies through efficient approximations. Here, we report the closedform expression of the NOs of twoelectron quantum rings, which are prototypical finiteextension systems and new starting points for the development of exchangecorrelation functionals in density functional theory. We also show that the natural occupation numbers for these twoelectron paradigms are in general nonvanishing and follow the same power law decay as atomic and molecular twoelectron systems.

Energy landscapes and persistent minima
View Description Hide DescriptionWe consider a coarsegraining of highdimensional potential energy landscapes based upon persistences, which correspond to lowest barrier heights to lowerenergy minima. Persistences can be calculated efficiently for local minima in kinetic transition networks that are based on stationary points of the prevailing energy landscape. The networks studied here represent peptides, proteins, nucleic acids, an atomic cluster, and a glassy system. Minima with high persistence values are likely to represent some form of alternative structural morphology, which, if appreciably populated at the prevailing temperature, could compete with the global minimum (defined as infinitely persistent). Threshold values on persistences (and in some cases equilibrium occupation probabilities) have therefore been used in this work to select subsets of minima, which were then analysed to see how well they can represent features of the full network. Simplified disconnectivity graphs showing only the selected minima can convey the funnelling (including any multiplefunnel) characteristics of the corresponding full graphs. The effect of the choice of persistence threshold on the reduced disconnectivity graphs was considered for a system with a hierarchical, glassy landscape. Sets of persistent minima were also found to be useful in comparing networks for the same system sampled under different conditions, using minimum oriented spanning forests.

A simple quasidiabatization scheme suitable for spectroscopic problems based on oneelectron properties of interacting states
View Description Hide DescriptionWe present a simple quasidiabatization scheme applicable to spectroscopic studies that can be applied using any wavefunction for which oneelectron properties and transition properties can be calculated. The method is based on rotation of a pair (or set) of adiabatic states to minimize the difference between the given transition property at a reference geometry of high symmetry (where the quasidiabatic states and adiabatic states coincide) and points of lower symmetry where quasidiabatic quantities are desired. Compared to other quasidiabatization techniques, the method requires no special coding, facilitates direct comparison between quasidiabatic quantities calculated using different types of wavefunctions, and is free of any selection of configurations in the definition of the quasidiabatic states. On the other hand, the method appears to be sensitive to multistate issues, unlike recent methods we have developed that use a configurational definition of quasidiabatic states. Results are presented and compared with two other recently developed quasidiabatization techniques.

Accurate molecular dynamics and nuclear quantum effects at low cost by multiple steps in real and imaginary time: Using density functional theory to accelerate wavefunction methods
View Description Hide DescriptionThe development and implementation of increasingly accurate methods for electronic structure calculations mean that, for many atomistic simulation problems, treating light nuclei as classical particles is now one of the most serious approximations. Even though recent developments have significantly reduced the overhead for modeling the quantum nature of the nuclei, the cost is still prohibitive when combined with advanced electronic structure methods. Here we present how multiple time step integrators can be combined with ringpolymer contraction techniques (effectively, multiple time stepping in imaginary time) to reduce virtually to zero the overhead of modelling nuclear quantum effects, while describing interatomic forces at high levels of electronic structure theory. This is demonstrated for a combination of MP2 and semilocal DFT applied to the Zundel cation. The approach can be seamlessly combined with other methods to reduce the computational cost of path integral calculations, such as highorder factorizations of the Boltzmann operator or generalized Langevin equation thermostats.

Ab initio molecular dynamics with nuclear quantum effects at classical cost: Ring polymer contraction for density functional theory
View Description Hide DescriptionPath integral molecular dynamics simulations, combined with an ab initio evaluation of interactions using electronic structure theory, incorporate the quantum mechanical nature of both the electrons and nuclei, which are essential to accurately describe systems containing light nuclei. However, path integral simulations have traditionally required a computational cost around two orders of magnitude greater than treating the nuclei classically, making them prohibitively costly for most applications. Here we show that the cost of path integral simulations can be dramatically reduced by extending our ring polymer contraction approach to ab initiomolecular dynamics simulations. By using density functional tight binding as a reference system, we show that our ring polymer contraction scheme gives rapid and systematic convergence to the full path integral density functional theory result. We demonstrate the efficiency of this approach in ab initio simulations of liquid water and the reactive protonated and deprotonated water dimer systems. We find that the vast majority of the nuclear quantum effects are accurately captured using contraction to just the ring polymer centroid, which requires the same number of density functional theory calculations as a classical simulation. Combined with a multiple time step scheme using the same reference system, which allows the time step to be increased, this approach is as fast as a typical classical ab initiomolecular dynamics simulation and 35× faster than a full path integral calculation, while still exactly including the quantum sampling of nuclei. This development thus offers a route to routinely include nuclear quantum effects in ab initiomolecular dynamics simulations at negligible computational cost.

Characterizing metastable states beyond energies and lifetimes: Dyson orbitals and transition dipole moments
View Description Hide DescriptionThe theoretical description of electronic resonances is extended beyond calculations of energies and lifetimes. We present the formalism for calculating Dyson orbitals and transition dipole moments within the equationofmotion coupledcluster singles and doubles method for electronattached states augmented by a complex absorbing potential (CAPEOMEACCSD). The capabilities of the new methodology are illustrated by calculations of Dyson orbitals of various transient anions. We also present calculations of transition dipole moments between transient and stable anionic states as well as between different transient states. Dyson orbitals characterize the differences between the initial neutral and final electronattached states without invoking the meanfield approximation. By extending the molecularorbital description to correlated manyelectron wave functions, they deliver qualitative insights into the character of resonance states. Dyson orbitals and transition moments are also needed for calculating experimental observables such as spectra and cross sections. Physically meaningful results for those quantities are obtained only in the framework of nonHermitian quantum mechanics, e.g., in the presence of a complex absorbing potential (CAP), when studying resonances. We investigate the dependence of Dyson orbitals and transition moments on the CAP strength and illustrate how Dyson orbitals help understand the properties of metastable species and how they are affected by replacing the usual scalar product by the socalled cproduct.
 Atoms, Molecules, and Clusters

Lattice constants and expansivities of gas hydrates from 10 K up to the stability limit
View Description Hide DescriptionThe lattice constants of hydrogenated and deuterated CH4, CO2, Xe (clathratestructure type I) and N2hydrates (clathratestructure type II) from 10 K up to the stability limit were established in neutron and synchrotron diffraction experiments and were used to derive the related thermal expansivities. The following results emerge from this analysis: (1) The differences of expansivities of structure type I and II hydrates are fairly small. (2) Despite the larger guestsize of CO2 as compared to methane, CO2hydrate has the smaller lattice constants at low temperatures, which is ascribed to the larger attractive guesthost interaction of the CO2water system. (3) The expansivity of CO2hydrate is larger than for CH4hydrate which leads to larger lattice constants for the former at temperatures above ∼150 K; this is likely due to the higher motional degrees of freedom of the CO2 guest molecules. (4) The cage occupancies of Xe and CO2hydrates affect significantly the lattice constants. (5) Similar to ice Ih, the deuterated compounds have generally slightly larger lattice constants which can be ascribed to the somewhat weaker Hbonding. (6) Compared to ice Ih, the high temperature expansivities are about 50% larger; in contrast to ice Ih and the empty hydrate, there is no negative thermal expansion at low temperature. (7) A comparison of the experimental results with lattice dynamical work, with models based on an Einstein oscillator model, and results from inelastic neutron scattering suggest that the contribution of the guest atoms’ vibrational energy to thermal expansion is important, most prominently for CO2 and Xehydrates.

On the mean kinetic energy of the proton in strong hydrogen bonded systems
View Description Hide DescriptionThe mean atomic kinetic energies of the proton, Ke(H), and of the deuteron, Ke(D), were calculated in moderate and strongly hydrogen bonded (HB) systems, such as the ferroelectric crystals of the KDP type (XH2PO4, X = K, Cs, Rb, Tl), the DKDP (XD2PO4, X = K, Cs, Rb) type, and the X3H(SO4)2 superprotonic conductors (X = K, Rb). All calculations utilized the simulated partialphonon density of states, deduced from density functional theory based firstprinciple calculations and from empirical lattice dynamics simulations in which the Coulomb, short range, covalent, and van der Waals interactions were accounted for. The presently calculated Ke(H) values for the two systems were found to be in excellent agreement with published values obtained by deep inelastic neutron scattering measurements carried out using the VESUVIO instrument of the Rutherford Laboratory, UK. The Ke(H) values of the M3H(SO4)2 compounds, in which the hydrogen bonds are centrosymmetric, are much lower than those of the KDP type crystals, in direct consistency with the oxygenoxygen distance ROO, being a measure of the HB strength.

Could the description on polynuclear superhalogens by DFT be comparable with highlevel ab initio results? A comparison between DFT and CCSD(T)
View Description Hide DescriptionA systematic density functional theory study including 17 exchangecorrelation functionals was performed on different types of superhalogens with high level coupledcluster single double including perturbative triple excitations (CCSD(T)) results as the reference. The superhalogens selected here cover the ranges from mononuclear to polynuclear structures and from structures with halogenatom ligands to those with nonhalogen ligands, e.g., [MgX3]^{−}, [Mg2X5]^{−}, and [Mg3X7]^{−} (X = F, Cl, CN). It is clearly indicated that three doublehybrid functionals B2TPLYP, B2GPPLYP, B2KPLYP as well as the rangeseparated hybrid functional ωB97X are capable of providing results which approach the accuracy at the CCSD(T) level. The basis set effect is usually moderate and, in most cases, it is enough to utilize the basis set of tripleξ quality, e.g., Def2TZVP. In addition, the results of the HF and MP2 method are also acceptable here, especially for polynuclear superhalogens where CCSD(T) is probably unpractical.

Ab initio 3D potential energy and dipole moment surfaces for the CH4–Ar complex: Collisioninduced intensity and dimer content
View Description Hide DescriptionWe present new threedimensional potential energy surface (PES) and dipole momentsurfaces (DMSs) for the CH4–Ar van der Waals system. Ab initio calculations of the PES and DMS were carried out using the closedshell single and doubleexcitation coupled cluster approach with noniterative perturbative treatment of triple excitations. The augmented correlationconsistent augccpVXZ (X = D,T,Q) basis sets were employed, and the energies obtained were then extrapolated to the complete basis set limit. The dipole momentsurface was obtained using augccpVTZ basis set augmented with midbond functions for better description of exchange interactions. The second mixed virial coefficient was calculated and compared to available experimental data. The equilibrium constant for true dimer formation was calculated using classical partition function based on the knowledge of ab initio PES. Temperature variations of the zeroth spectral moment of the rototranslational collisioninduced band as well as its true dimer constituent were traced with the use of the Boltzmannweighted squared induced dipole properly integrated over respective phase space domains. Height profiles for N2–N2, N2–H2, CH4–N2, (CH4)2, and CH4–Ar true bound dimers in Titan’s atmosphere were calculated with the use of reliable ab initio PESs.

Timeresolved photoelectron spectroscopy of a dinuclear Pt(II) complex: Tunneling autodetachment from both singlet and triplet excited states of a molecular dianion
View Description Hide DescriptionTimeresolved pumpprobe photoelectron spectroscopy has been used to study the relaxation dynamics of gaseous [Pt2(μP2O5H2)4 + 2H]^{2−} after population of its first singlet excited state by 388 nm femtosecond laser irradiation. In contrast to the fluorescence and phosphorescence observed in condensed phase, a significant fraction of the photoexcited isolated dianions decays by electron loss to form the corresponding monoanions. Our transient photoelectron data reveal an ultrafast decay of the initially excited singlet ^{1}A2u state and concomitant rise in population of the triplet ^{3}A2u state, via subpicosecond intersystem crossing (ISC). We find that both of the electronically excited states are metastably bound behind a repulsive Coulomb barrier and can decay via delayed autodetachment to yield electrons with characteristic kinetic energies. While excited statetunneling detachment (ESETD) from the singlet ^{1}A2u state takes only a few picoseconds, ESETD from the triplet ^{3}A2u state is much slower and proceeds on a time scale of hundreds of nanoseconds. The ISC rate in the gas phase is significantly higher than in solution, which can be rationalized in terms of changes to the energy dissipation mechanism in the absence of solvent molecules. [Pt2(μP2O5H2)4 + 2H]^{2−} is the first example of a photoexcited multianion for which ESETD has been observed following ISC.