Volume 134, Issue 24, 28 June 2011

Studies on confined water are important not only from the viewpoint of scientific interest but also for the development of new nanoscale devices. In this work, we aimed to clarify the properties of confined water in the cylindrical pores of singlewalled carbon nanotubes(SWCNTs) that had diameters in the range of 1.46 to 2.40 nm. A combination of xray diffraction(XRD),nuclear magnetic resonance, and electrical resistance measurements revealed that water inside SWCNTs with diameters between 1.68 and 2.40 nm undergoes a wetdry type transition with the lowering of temperature; below the transition temperature T _{wd}, water was ejected from the SWCNTs.T _{wd} increased with increasing SWCNT diameter D. For the SWCNTs with D = 1.68, 2.00, 2.18, and 2.40 nm, T _{wd} obtained by the XRD measurements were 218, 225, 236, and 237 K, respectively. We performed a systematic study on finite length SWCNT systems using classical molecular dynamics calculations to clarify the effect of open ends of the SWCNTs and water content on the waterstructure. It was found that icestructures that were formed at low temperatures were strongly affected by the bore diameter, a = D − σ_{OC}, where σ_{OC} is gap distance between the SWCNT and oxygen atom in water, and the number of water molecules in the system. In small pores (a < 1.02 nm), tubule ices or the socalled icenanotubes(ice NTs) were formed irrespective of the water content. On the other hand, in larger pores (a > 1.10 nm) with small water content, filled water clusters were formed leaving some empty space in the SWCNT pore, which grew to fill the pore with increasing water content. For pores with sizes in between these two regimes (1.02 < a < 1.10 nm), tubule ice also appeared with small water content and grew with increasing water content. However, once the tubule ice filled the entire SWCNT pore, further increase in the water content resulted in encapsulation of the additional water molecules inside the tubule ice. Corresponding XRD measurements on SWCNTs with a mean diameter of 1.46 nm strongly suggested the presence of such a filled structure.
 COMMUNICATIONS


Communication: Universality of the melting curves for a wide range of interaction potentials
View Description Hide DescriptionWe demonstrate that the melting curves of various model systems of interacting particles collapse to (or are located very close to) a universal master curve on a plane of appropriately chosen scaled variables. The physics behind this universality is discussed. An equation for the emerging “universal melting curve” is proposed. The obtained results can be used to approximately predict melting of various substances in a wide range of conditions.

Communication: A new spectroscopic window on hydroxyl radicals using UV + VUV resonant ionization
View Description Hide DescriptionA 1 + 1^{′}multiphoton ionization (MPI) detection scheme for OH radicals is presented. The spectroscopic approach combines initial excitation on the wellcharacterized A^{2}Σ^{+}–X^{2}Π band system with vacuum ultraviolet (VUV)ionization via autoionizing Rydberg states that converge on the OH^{+} A^{3}Π ion state. Jetcooled MPI spectra on the (1,0) and (2,0) bands show anomalous rotational line intensities, while initial excitation on the (0,0) band does not lead to detectable OH^{+} ions. The onset of ionization with the (1,0) band is attributed to an energetic threshold; the combined UV + VUVphoton energies are above the first member of the autoionizing (A^{3}Π)nd Rydberg series. Comparison of the OH 1 + 1^{′} MPI signal with that from single photonVUVionization of NO indicates that the cross section for photoionization from OH A^{2}Σ^{+}, v ^{′} = 1 is on the order of 10^{−17} cm^{2}.

Communication: Linearexpansion shooting techniques for accelerating selfconsistent field convergence
View Description Hide DescriptionBased on the corrected HohenbergKohnSham total energy density functional[Y. A. Zhang and Y. A. Wang, J. Chem. Phys.130, 144116 (2009)]10.1063/1.3104662, we have developed two linearexpansion shooting techniques (LIST)– direct LIST (LISTd) and indirect LIST (LISTi), to accelerate the convergence of selfconsistent field (SCF) calculations. Case studies show that overall LISTi is the most robust and efficient algorithm for accelerating SCF convergence, whereas LISTd is advantageous in the early stage of an SCF process. More importantly, LISTi outperforms Pulay's direct inversion in the iterative subspace (DIIS) [P. Pulay, J. Comput. Chem.3, 556 (1982)]10.1002/jcc.540030413 and its two recent improvements, energyDIIS [K. N. Kudin, G. E. Scuseria, and E. Cancès, J. Chem. Phys.116, 8255 (2002)]10.1063/1.1470195 and augmented RoothaanHall energyDIIS [X. Hu and W. Yang, J. Chem. Phys.132, 054109 (2010)]10.1063/1.3304922.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

A multiconfigurational timedependent HartreeFock method for excited electronic states. I. General formalism and application to openshell states
View Description Hide DescriptionThe solution of the timedependent Schrödinger equation for systems of interacting electrons is generally a prohibitive task, for which approximate methods are necessary. Popular approaches, such as the timedependent HartreeFock (TDHF) approximation and timedependent density functional theory (TDDFT), are essentially singleconfigurational schemes. TDHF is by construction incapable of fully accounting for the excited character of the electronic states involved in many physical processes of interest; TDDFT, although exact in principle, is limited by the currently available exchangecorrelation functionals. On the other hand, multiconfigurational methods, such as the multiconfigurational timedependent HartreeFock (MCTDHF) approach, provide an accurate description of the excited states and can be systematically improved. However, the computational cost becomes prohibitive as the number of degrees of freedom increases, and thus, at present, the MCTDHF method is only practical for fewelectron systems. In this work, we propose an alternative approach which effectively establishes a compromise between efficiency and accuracy, by retaining the smallest possible number of configurations that catches the essential features of the electronic wavefunction. Based on a timedependent variational principle, we derive the MCTDHF working equation for a multiconfigurational expansion with fixed coefficients and specialise to the case of general openshell states, which are relevant for many physical processes of interest.

A multiconfigurational timedependent HartreeFock method for excited electronic states. II. Coulomb interaction effects in single conjugated polymer chains
View Description Hide DescriptionConjugated polymers have attracted considerable attention in the last few decades due to their potential for optoelectronic applications. A key step that needs optimisation is charge carrier separation following photoexcitation. To understand better the dynamics of the exciton prior to charge separation, we have performed simulations of the formation and dynamics of localised excitations in single conjugated polymer strands. We use a nonadiabatic molecular dynamics method which allows for the coupled evolution of the nuclear degrees of freedom and of multiconfigurational electronic wavefunctions. We show the relaxation of electronhole pairs to form excitons and oppositely chargedpolaron pairs and discuss the modifications to the relaxation process predicted by the inclusion of the Coulomb interaction between the carriers. The issue of charge photogeneration in conjugated polymers in dilute solution is also addressed.

Quantum effects in energy and charge transfer in an artificial photosynthetic complex
View Description Hide DescriptionWe investigate the quantum dynamics of energy and charge transfer in a wheelshaped artificial photosynthetic antennareaction center complex. This complex consists of six lightharvesting chromophores and an electronacceptor fullerene. To describe quantum effects on a femtosecond time scale, we derive the set of exact nonMarkovian equations for the Heisenberg operators of this photosynthetic complex in contact with a Gaussian heat bath. With these equations we can analyze the regime of strong systembath interactions, where reorganization energies are of the order of the intersite exciton couplings. We show that the energy of the initially excited antenna chromophores is efficiently funneled to the porphyrinfullerene reaction center, where a chargeseparated state is set up in a few picoseconds, with a quantum yield of the order of 95%. In the singleexciton regime, with one antenna chromophore being initially excited, we observe quantum beatings of energy between two resonant antenna chromophores with a decoherence time of ∼100 fs. We also analyze the doubleexciton regime, when two porphyrin molecules involved in the reaction center are initially excited. In this regime we obtain pronounced quantum oscillations of the charge on the fullerene molecule with a decoherence time of about 20 fs (at liquid nitrogen temperatures). These results show a way to directly detect quantum effects in artificial photosynthetic systems.

A comparison of accelerators for direct energy minimization in electronic structure calculations
View Description Hide DescriptionWe compare three different methods for direct energy minimization in electronic structure calculations where the gradient of the energy functional with respect to the molecular orbitals is available. These methods make use of the preconditioned gradient to increase robustness. An orbital transformation is used to ensure that the orthogonality constraint on the orbitals remains satisfied when using standard minimization methods. In addition, we propose an adaptive scheme for estimating the curvature of the energy functional to increase the performance of a line search free quasiNewton method. We show that the performance of all methods is similar when robustness of the methods is ensured.

Index k saddles and dividing surfaces in phase space with applications to isomerization dynamics
View Description Hide DescriptionIn this paper, we continue our studies of the phase space geometry and dynamics associated with index k saddles (k > 1) of the potential energy surface. Using PoincaréBirkhoff normal form (NF) theory, we give an explicit formula for a “dividing surface” in phase space, i.e., a codimension one surface (within the energy shell) through which all trajectories that “cross” the region of the index k saddle must pass. With a generic nonresonance assumption, the normal form provides k (approximate) integrals that describe the saddle dynamics in a neighborhood of the index k saddle. These integrals provide a symbolic description of all trajectories that pass through a neighborhood of the saddle. We give a parametrization of the dividing surface which is used as the basis for a numerical method to sample the dividing surface. Our techniques are applied to isomerizationdynamics on a potential energy surface having four minima; two symmetry related pairs of minima are connected by low energy index 1 saddles, with the pairs themselves connected via higher energy index 1 saddles and an index 2 saddle at the origin. We compute and sample the dividing surface and show that our approach enables us to distinguish between concerted crossing (“hilltop crossing”) isomerizing trajectories and those trajectories that are not concerted crossing (potentially sequentially isomerizing trajectories). We then consider the effect of additional “bath modes” on the dynamics, by a study of a four degreeoffreedom system. For this system we show that the normal form and dividing surface can be realized and sampled and that, using the approximate integrals of motion and our symbolic description of trajectories, we are able to choose initial conditions corresponding to concerted crossing isomerizing trajectories and (potentially) sequentially isomerizing trajectories.

Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems
View Description Hide DescriptionPadé spectrum decomposition is an optimal sumoverpoles expansion scheme of Fermi function and Bose function [J. Hu, R. X. Xu, and Y. J. Yan, J. Chem. Phys.133, 101106 (2010)]10.1063/1.3484491. In this work, we report two additional members to this family, from which the best among all sumoverpoles methods could be chosen for different cases of application. Methods are developed for determining these three Padé spectrum decomposition expansions at machine precision via simple algorithms. We exemplify the applications of present development with optimal construction of hierarchical equationsofmotion formulations for nonperturbative quantum dissipation and quantum transport dynamics. Numerical demonstrations are given for two systems. One is the transient transport current to an interacting quantumdots system, together with the involved highorder cotunneling dynamics. Another is the nonMarkovian dynamics of a spinboson system.

Dynamical reweighting: Improved estimates of dynamical properties from simulations at multiple temperatures
View Description Hide DescriptionDynamical averages based on functionals of dynamical trajectories, such as timecorrelation functions, play an important role in determining kinetic or transport properties of matter. At temperatures of interest, the expectations of these quantities are often dominated by contributions from rare events, making the precise calculation of these quantities by molecular dynamics simulation difficult. Here, we present a reweighting method for combining simulations from multiple temperatures (or from simulated or parallel tempering simulations) to compute an optimal estimate of the dynamical properties at the temperature of interest without the need to invoke an approximate kinetic model (such as the Arrhenius law). Continuous and differentiable estimates of these expectations at any temperature in the sampled range can also be computed, along with an assessment of the associated statistical uncertainty. For rare events, aggregating data from multiple temperatures can produce an estimate with the desired precision at greatly reduced computational cost compared with simulations conducted at a single temperature. Here, we describe use of the method for the canonical (NVT) ensemble using four common models of dynamics (canonical distribution of Hamiltonian trajectories, Andersen thermostatting, Langevin, and overdamped Langevin or Brownian dynamics), but it can be applied to anythermodynamic ensemble provided the ratio of path probabilities at different temperatures can be computed. To illustrate the method, we compute a timecorrelation function for solvated terminallyblocked alanine peptide across a range of temperatures using trajectories harvested using a modified parallel tempering protocol.

Optimal use of data in parallel tempering simulations for the construction of discretestate Markov models of biomolecular dynamics
View Description Hide DescriptionParallel tempering (PT) molecular dynamics simulations have been extensively investigated as a means of efficient sampling of the configurations of biomolecular systems. Recent work has demonstrated how the short physical trajectories generated in PT simulations of biomolecules can be used to construct the Markovmodels describing biomolecular dynamics at each simulated temperature. While this approach describes the temperaturedependent kinetics, it does not make optimal use of all available PT data, instead estimating the rates at a given temperature using only data from that temperature. This can be problematic, as some relevant transitions or states may not be sufficiently sampled at the temperature of interest, but might be readily sampled at nearby temperatures. Further, the comparison of temperaturedependent properties can suffer from the false assumption that data collected from different temperatures are uncorrelated. We propose here a strategy in which, by a simple modification of the PT protocol, the harvested trajectories can be reweighted, permitting data from all temperatures to contribute to the estimated kinetic model. The method reduces the statistical uncertainty in the kinetic model relative to the single temperature approach and provides estimates of transition probabilities even for transitions not observed at the temperature of interest. Further, the method allows the kinetics to be estimated at temperatures other than those at which simulations were run. We illustrate this method by applying it to the generation of a Markovmodel of the conformational dynamics of the solvated terminally blocked alanine peptide.

Ab initio study of excited state electronic circular dichroism. Two prototype cases: Methyl oxirane and R(+)1,1^{′}bi(2naphthol)
View Description Hide DescriptionA computational approach to the calculation of excited state electronic circular dichroism (ESECD) spectra of chiral molecules is discussed. Frequency dependent quadratic response theory is employed to compute the rotatory strength for transitions between excited electronic states, by employing both a magnetic gauge dependent and a (velocitybased) magnetic gauge independent approach. Application is made to the lowest excited states of two prototypical chiral molecules, propylene oxide, also known as 1,2epoxypropane or methyl oxirane, and R(+)1,1^{′}bi(2naphthol), or BINOL. The dependence of the rotatory strength for transitions between the lowest three excited states of methyl oxirane upon the quality and extension of the basis set is analyzed, by employing a hierarchy of correlation consistent basis sets. Once established that basis sets of at least triple zeta quality, and at least doubly augmented, are sufficient to ensure sufficiently converged results, at least at the HartreeFock selfconsistent field (HFSCF) level, the rotatory strengths for all transitions between the lowest excited electronic states of methyl oxirane are computed and analyzed, employing HFSCF, and density functional theory(DFT) electronic structure models. For DFT, both the popular B3LYP and its recently highly successful CAMB3LYP extension are exploited. The strong dependence of the spectra upon electron correlation is highlighted. A HFSCF and DFT study is carried out also for BINOL, a system where excited states show the typical pairing structure arising from the interaction of the two monomeric moieties, and whose conformational changes following photoexcitation were studied recently with via timeresolved CD.

Benchmark of density functional theory methods on the prediction of bond energies and bond distances of noblegas containing molecules
View Description Hide DescriptionWe have tested three pure density functional theory(DFT) functionals, BLYP, MPWPW91, MPWB95, and ten hybrid DFT functionals, B3LYP, B3P86, B98, MPW1B95, MPW1PW91, BMK, M052X, M062X, B2GPPLYP, and DSDBLYP with a series of commonly used basis sets on the performance of predicting the bond energies and bond distances of 31 small neutral noblegas containing molecules. The reference structures were obtained using the CCSD(T)/augccpVTZ theory and the reference energies were based on the calculation at the CCSD(T)/CBS level. While in general the hybrid functionals performed significantly better than the pure functionals, our tests showed a range of performance by these hybrid functionals. For the bond energies, the MPW1B95/6311+G(2df,2pd), BMK/augccpVTZ, B2GPPLYP/augccpVTZ, and DSDBLYP/augccpVTZ methods stood out with mean unsigned errors of 2.0−2.3 kcal/mol per molecule. For the bond distances, the MPW1B95/6311+G(2df,2pd), MPW1PW91/6311+G(2df,2pd), and B3P86/6311+G(2df,2pd), DSDBLYP/6311+G(2df,2pd), and DSDBLYP/augccpVTZ methods stood out with mean unsigned errors of 0.008−0.013 Å per bond. The current study showed that a careful selection of DFT functionals is very important in the study of noblegas chemistry, and the most recommended methods are MPW1B95/6311+G(2df,2pd) and DSDBLYP/augccpVTZ.

Comparison of two adaptive temperaturebased replica exchange methods applied to a sharp phase transition of protein unfoldingfolding
View Description Hide DescriptionTemperaturebased replica exchange (TReX) enhances sampling of molecular dynamics simulations by autonomously heating and cooling simulation clients via a Metropolis exchange criterion. A pathological case for TReX can occur when a change in state (e.g., folding to unfolding of a protein) has a large energetic difference over a short temperature interval leading to insufficient exchanges amongst replica clients near the transition temperature. One solution is to allow the temperature set to dynamically adapt in the temperature space, thereby enriching the population of clients near the transition temperature. In this work, we evaluated two approaches for adapting the temperature set: a method that equalizes exchange rates over all neighbor temperature pairs and a method that attempts to induce clients to visit all temperatures (dubbed “current maximization”) by positioning many clients at or near the transition temperature. As a test case, we simulated the 57residue SH3 domain of alphaspectrin. Exchange rate equalization yielded the same unfoldingfolding transition temperature as fixedtemperature ReX with much smoother convergence of this value. Surprisingly, the current maximization method yielded a significantly lower transition temperature, in close agreement with experimental observation, likely due to more extensive sampling of the transition state.

Analytical evaluation of firstorder electrical properties based on the spinfree DiracCoulomb Hamiltonian
View Description Hide DescriptionWe report an analytical scheme for the calculation of firstorder electrical properties using the spinfree DiracCoulomb (SFDC) Hamiltonian, thereby exploiting the welldeveloped densitymatrix formulations in nonrelativistic coupledcluster (CC) derivative theory. Orbital relaxation effects are fully accounted for by including the relaxation of the correlated orbitals with respect to orbitals of all types, viz., frozencore, occupied, virtual, and negative energy state orbitals. To demonstrate the applicability of the presented scheme, we report benchmark calculations for firstorder electrical properties of the hydrogen halides, HX with X = F, Cl, Br, I, At, and a first application to the iodo(fluoro)methanes, CH_{ n }F_{3 − n }I, n = 0–3. The results obtained from the SFDC calculations are compared to those from nonrelativistic calculations, those obtained via leadingorder direct perturbation theory as well as those from full DiracCoulomb calculations. It is shown that the full inclusion of spinfree (SF) relativistic effects is necessary to obtain accurate firstorder electrical properties in the presence of fifthrow elements. The SFDC scheme is also recommended for applications to systems containing lighter elements because it introduces no extra cost in the ratedetermining steps of a CC calculation in comparison to the nonrelativistic case. On the other hand, spinorbit contributions are generally small for firstorder electrical properties of closedshell molecules and may be handled efficiently by means of perturbation theory.

Brueckner doubles coupled cluster method with the polarizable continuum model of solvation
View Description Hide DescriptionWe present the theory and implementation for computing the (free) energy and its analytical gradients with the Brueckner doubles (BD) coupled cluster method in solution, in combination with the polarizable continuum model of solvation (PCM). The complete model, called PTED, and an efficient approximation, called PTE, are introduced and tested with numerical examples. Implementation details are also discussed. A comparison with the coupledcluster singles and doubles CCSDPCMPTED and CCSDPCMPTE schemes, which use HartreeFock (HF) orbitals, is presented. The results show that the two PTED approaches are mostly equivalent, while BDPCMPTE is shown to be superior to the corresponding CCSD scheme when the HF reference wave function is unstable. The BDPCMPTE scheme, whose computational cost is equivalent to gas phase BD, is therefore a promising approach to study molecular systems with complicated electronic structure in solution.

Decoherence and surface hopping: When can averaging over initial conditions help capture the effects of wave packet separation?
View Description Hide DescriptionFewestswitches surfacehopping (FSSH) is a popular nonadiabaticdynamics method which treats nuclei with classical mechanics and electrons with quantum mechanics. In order to simulate the motion of a wave packet as accurately as possible, standard FSSH requires a stochastic sampling of the trajectories over a distribution of initial conditions corresponding, e.g., to the Wigner distribution of the initial quantum wave packet. Although it is wellknown that FSSH does not properly account for decoherence effects, there is some confusion in the literature about whether or not this averaging over a distribution of initial conditions can approximate some of the effects of decoherence. In this paper, we not only show that averaging over initial conditions does not generally account for decoherence, but also why it fails to do so. We also show how an apparent improvement in accuracy can be obtained for a fortuitous choice of model problems, even though this improvement is not possible, in general. For a basic set of onedimensional and twodimensional examples, we find significantly improved results using our recently introduced augmented FSSH algorithm.

A general formulation for the efficient evaluation of nelectron integrals over products of Gaussian charge distributions with Gaussian geminal functions
View Description Hide DescriptionIn this work, we present a general formulation for the evaluation of manyelectron integrals which arise when traditional one particle expansions are augmented with explicitly correlated Gaussian geminal functions. The integrand is expressed as a product of charge distributions, one for each electron, multiplied by one or more Gaussian geminal factors. Our formulation begins by focusing on the quadratic form that arises in the general nelectron integral. Using the Rys polynomial method for the evaluation of potential energy integrals, we derive a general formula for the evaluation of any nelectron integral. This general expression contains four parameters ω, θ, v, and h, which can be evaluated by an examination of the general quadratic form. Our analysis contains general expressions for any nelectron integral over stype functions as well as the recursion needed to build up arbitrary angular momentum. The general recursion relation requires at most n + 1 terms for any nelectron integral. To illustrate the general method, we develop explicit expressions for the evaluation of two, three, and four particle electron repulsion integrals as well as two and three particle overlap and nuclear attraction integrals. We conclude our exposition with a discussion of a preliminary computational implementation as well as general computational requirements. Implementation on parallel computers is briefly discussed.

Mixed quantumclassical simulations of charge transport in organic materials: Numerical benchmark of the SuSchriefferHeeger model
View Description Hide DescriptionThe electronphonon coupling is critical in determining the intrinsic charge carrier and exciton transport properties in organic materials. In this study, we consider a SuSchriefferHeeger (SSH) model for molecular crystals, and perform numerical benchmark studies for different strategies of simulating the mixed quantumclassical dynamics. These methods, which differ in the selection of initial conditions and the representation used to solve the time evolution of the quantum carriers, are shown to yield similar equilibrium diffusion properties. A hybrid approach combining molecular dynamics simulations of nuclear motion and quantumchemical calculations of the electronic Hamiltonian at each geometric configuration appears as an attractive strategy to modelchargedynamics in large size systems “on the fly,” yet it relies on the assumption that the quantum carriers do not impact the nuclear dynamics. We find that such an approximation systematically results in overestimated chargecarriermobilities, with the associated error being negligible when the roomtemperature mobility exceeds ∼4.8 cm^{2}/Vs (∼0.14 cm^{2}/Vs) in onedimensional (twodimensional) crystals.

Development and application of the analytical energy gradient for the normalized elimination of the small component method
View Description Hide DescriptionThe analytical energy gradient of the normalized elimination of the small component (NESC) method is derived for the first time and implemented for the routine calculation of NESC geometries and other first order molecular properties. Essential for the derivation is the correct calculation of the transformation matrix U relating the small component to the pseudolarge component of the wavefunction. The exact form of is derived and its contribution to the analytical energy gradient is investigated. The influence of a finite nucleus model and that of the picture change is determined. Different ways of speeding up the calculation of the NESC gradient are tested. It is shown that first order properties can routinely be calculated in combination with HartreeFock, density functional theory(DFT),coupled clustertheory, or any electron correlation corrected quantum chemical method, provided the NESC Hamiltonian is determined in an efficient, but nevertheless accurate way. The general applicability of the analytical NESC gradient is demonstrated by benchmark calculations for NESC/CCSD (coupled cluster with all single and double excitation) and NESC/DFT involving up to 800 basis functions.