Volume 130, Issue 19, 21 May 2009

A new milestoning procedure using Voronoi tessellations is proposed. In the new procedure, the edges of Voronoi cells are used as milestones, and the necessary kinetic information about the transitions between the milestones is calculated by running molecular dynamics (MD) simulations restricted to these cells. Like the traditional milestoning technique, the new procedure offers a reduced description of the original dynamics and permits to efficiently compute the various quantities necessary in this description. However, unlike traditional milestoning, the new procedure does not require to reinitialize trajectories from the milestones, and thereby it avoids the approximation made in traditional milestoning that the distribution for reinitialization is the equilibrium one. In this paper we concentrate on Markovian milestoning, which we show to be valid under suitable assumptions, and we explain how to estimate the rate matrix of transitions between the milestones from data collected from the MD trajectories in the Voronoi cells. The rate matrix can then be used to compute mean first passage times between milestones and reaction rates. The procedure is first illustrated on testcase examples in two dimensions and then applied to study the kinetics of protein insertion into a lipid bilayer by means of a coarsegrained model.
 ARTICLES

 Theoretical Methods and Algorithms

Markovian milestoning with Voronoi tessellations
View Description Hide DescriptionA new milestoning procedure using Voronoi tessellations is proposed. In the new procedure, the edges of Voronoi cells are used as milestones, and the necessary kinetic information about the transitions between the milestones is calculated by running molecular dynamics (MD) simulations restricted to these cells. Like the traditional milestoning technique, the new procedure offers a reduced description of the original dynamics and permits to efficiently compute the various quantities necessary in this description. However, unlike traditional milestoning, the new procedure does not require to reinitialize trajectories from the milestones, and thereby it avoids the approximation made in traditional milestoning that the distribution for reinitialization is the equilibrium one. In this paper we concentrate on Markovian milestoning, which we show to be valid under suitable assumptions, and we explain how to estimate the rate matrix of transitions between the milestones from data collected from the MD trajectories in the Voronoi cells. The rate matrix can then be used to compute mean first passage times between milestones and reaction rates. The procedure is first illustrated on testcase examples in two dimensions and then applied to study the kinetics of protein insertion into a lipid bilayer by means of a coarsegrained model.

Twocomponent relativistic density functional method for computing nonsingular complex linear response of molecules based on the zeroth order regular approximation
View Description Hide DescriptionWe report the implementation of a frequencydependent twocomponent relativistic density functional theory method based on the zeroth order regular approximation (ZORA) for computations of complex linear response of molecules including spinorbit coupling. The implementation is based on Slatertype atomic orbital basis functions and makes extensive use of density fitting techniques. The complex response is obtained by applying damping in the computations. The method is validated by computations of the real and imaginary part of the static and dynamic polarizability of group 12 atoms, of a number of heavyatom diatomic molecules, of a range of two and threedimensional gold clusters, and of group 8 oxides and metallocenes. Simulated spectra—a plot of extinction coefficient as a function of frequency—obtained from the isotropic imaginary polarizability are compared to broadened spectra obtained from twocomponent ZORA excitation energies and oscillator strengths.

Revisiting the finite temperature string method for the calculation of reaction tubes and free energies
View Description Hide DescriptionAn improved and simplified version of the finite temperature string (FTS) method [W. E, W. Ren, and E. VandenEijnden, J. Phys. Chem. B109, 6688 (2005)] is proposed. Like the original approach, the new method is a scheme to calculate the principal curves associated with the Boltzmann–Gibbs probability distribution of the system, i.e., the curves which are such that their intersection with the hyperplanes perpendicular to themselves coincides with the expected position of the system in these planes (where perpendicular is understood with respect to the appropriate metric). Unlike more standard paths such as the minimum energy path or the minimum free energy path, the location of the principal curve depends on global features of the energy or the free energy landscapes and thereby may remain appropriate in situations where the landscape is rough on the thermal energy scale and/or entropic effects related to the width of the reaction channels matter. Instead of using constrained sampling in hyperplanes as in the original FTS, the new method calculates the principal curve via sampling in the Voronoi tessellation whose generating points are the discretization points along this curve. As shown here, this modification results in greater algorithmic simplicity. As a byproduct, it also gives the free energy associated with the Voronoi tessellation. The new method can be applied both in the original Cartesian space of the system or in a set of collective variables. We illustrate FTS on testcase examples and apply it to the study of conformational transitions of the nitrogen regulatory protein C receiver domain using an elastic network model and to the isomerization of solvated alanine dipeptide.

Exact ionization potentials from wavefunction asymptotics: The extended Koopmans’ theorem, revisited
View Description Hide DescriptionA simple explanation is given for the exactness of the extended Koopmans’ theorem, (EKT) for computing the removal energy of any manyelectron system to the lowestenergy ground state ion of a given symmetry. In particular, by removing the electron from a “removal orbital” of appropriate symmetry that is concentrated in the asymptotic region, one obtains the exact ionization potential and the exact Dyson orbital for the corresponding state of the ion. It is argued that the EKT is not restricted to manyelectron systems but holds for any finite manybody system, provided that the interaction vanishes for increasing interparticle distance. A necessary and sufficient condition for the validity of the EKT for any state (not just the lowestenergy states of a given symmetry) in terms of the thirdorder reduced density matrix is stated and derived.

Preparation of manybody states for quantum simulation
View Description Hide DescriptionWhile quantum computers are capable of simulating many quantum systems efficiently, the simulation algorithms must begin with the preparation of an appropriate initial state. We present a method for generating physically relevant quantum states on a lattice in real space. In particular, the present algorithm is able to prepare general pure and mixed manyparticle states of any number of particles. It relies on a procedure for converting from a secondquantized state to its firstquantized counterpart. The algorithm is efficient in that it operates in time that is polynomial in all the essential descriptors of the system, the number of particles, the resolution of the lattice, and the inverse of the maximum final error. This scaling holds under the assumption that the wave function to be prepared is bounded or its indefinite integral is known and that the Fock operator of the system is efficiently simulatable.

Band structures built by the elongation method
View Description Hide DescriptionA recently proposed approach for extracting band structures from finitecluster calculations is improved so that (avoided) band crossings can be handled and the problems related to socalled doublings and holes are reduced. In particular, we demonstrate how the method can be combined with the elongation method for the finitesystem calculations and apply it to extracting band structures for polymers from oligomer calculations. As illustrations of the approach we discuss a chain of water molecules, polyacetylene, polyethylene, and a BNnanotube without and with an impurity.

A study of cumulant approximations to electron valence multireference perturbation theory
View Description Hide DescriptionWe investigate the possibility of reducing the complexity of multireference perturbation theory through cumulant based approximations to the highorder density matrices that appear in such theories. Our test cases show that while the cumulant approximated forms are degraded in accuracy relative to the parent theory and exhibit intruder state problems that must be carefully handled, they may provide a route to a simple estimation of dynamic correlation when the parent perturbation theory is infeasible. Nonetheless, further work is clearly needed on better approximations to the denominators in the perturbation theory.

Parallel computation of coupledcluster hyperpolarizabilities
View Description Hide DescriptionStatic hyperpolarizabilities of molecules (water, acetonitrile, chloroform, and paranitroaniline) are calculated with large basis sets using coupledcluster response theory and compared to four common density functional theory methods. These results reveal which methods and basis sets are appropriate for nonlinear optical studies for different types of molecules and provide a means for estimating errors from the quantum chemical approximation when including vibrational contributions or solvent effects at the QM/MM level. The largest calculation reported, which was for 72 electrons in 812 functions at symmetry, took only a few hours on 256 nodes demonstrating that even larger calculations are quite feasible using modern supercomputers.

Magnetic anisotropy from density functional calculations. Comparison of different approaches: acetate as a test case
View Description Hide DescriptionMagnetic anisotropy is the capability of a system in a triplet or higher spin state to store magnetic information. Although the source of the magnetic anisotropy is the zerofield splitting of the ground state of the system, there is a difference between these two quantities that has to be fully rationalized before one makes comparisons. This is especially important for small spins such as triplets, where the magnetic anisotropy energy is only half of the zerofield splitting. Density functional calculations of magnetic anisotropy energies correspond to a highfield limit where the spins are aligned by the external magnetic field. Data are presented for the wellstudied molecular magnet acetate. Both perturbative and selfconsistent treatments, different quasirelativistic Hamiltonians (zeroth order regular approximation, Douglas–Kroll, effective core potentials) and exchangecorrelation functionals are compared. It is shown that some effects usually considered minor, such as the inclusion of the exchangecorrelation potential in the effective oneparticle spinorbit operator, lead to sizable differences when computing magnetic anisotropy energies. Higherorder contributions, that is, the difference between selfconsistent and perturbative results, increase the magnetic anisotropy energy somewhat but do not introduce sizeable quartic terms or an inplane anisotropy. In numerical experiments, on can switch off and on spinorbit coupling at individual atomic sites. This procedure yields singlesite contributions to the overall magnetic anisotropy energy that could be used as parameters in phenomenological spin Hamiltonians. If ferrimagnetic systems are treated with broken symmetry density functional methods where the Kohn–Sham reference function is not a spin eigenfunction, corrections are needed which depend on the size of the exchange couplings in the system and must therefore be evaluated case by case.

Nested variant of the method of moments of coupled cluster equations for vertical excitation energies and excitedstate potential energy surfaces
View Description Hide DescriptionIn this article we discuss the problem of proper balancing of the noniterative corrections to the ground and excitedstateenergies obtained with approximate coupled cluster (CC) and equationofmotion CC (EOMCC) approaches. It is demonstrated that for a class of excited states dominated by single excitations and for states with medium doubly excited component, the newly introduced nested variant of the method of moments of CC equations provides mathematically rigorous way of balancing the ground and excitedstate correlation effects. The resulting noniterative methodology accounting for the effect of triples is tested using its parallel implementation on the systems, for which iterative CC/EOMCC calculations with full inclusion of triply excited configurations or their most important subset are numerically feasible.

Approximate normal mode analysis based on vibrational subsystem analysis with high accuracy and efficiency
View Description Hide DescriptionNormal modeanalysis (NMA) has been proven valuable in modeling slow conformational dynamics of biomolecular structures beyond the reach of direct molecular simulations. However, it remains computationally expensive to directly solve normal modes for large biomolecular systems. In this study, we have evaluated the accuracy and efficiency of two approximate NMA protocols—one based on our recently proposed vibrational subsystem analysis (VSA), the other based on the rotation translation block (RTB), in comparison with standard NMA that directly solves a full Hessian matrix. By properly accounting for flexibility within blocks of residues or atoms based on a subsystemenvironment partition, VSAbased NMA has attained a much higher accuracy than RTB and much lower computing cost than standard NMA. Therefore, VSA enables accurate and efficient calculations of normal modes from allatom or coarsegrained potential functions, which promise to improve conformational sampling driven by lowfrequency normal modes.

Enhanced sampling in generalized ensemble with large gap of sampling parameter: Case study in temperature space random walk
View Description Hide DescriptionWe present an efficient sampling method for computing a partition function and accelerating configuration sampling. The method performs a random walk in the space, with being any thermodynamic variable that characterizes a canonical ensemble such as the reciprocal temperature or any variable that the Hamiltonian depends on. The partition function is determined by minimizing the difference of the thermal conjugates of (the energy in the case of ), defined as the difference between the value from the dynamically updated derivatives of the partition function and the value directly measured from simulation. Higherorder derivatives of the partition function are included to enhance the Brownian motion in the space. The method is much less sensitive to the system size, and to the size of window than other methods. On the two dimensional Ising model, it is shown that the method asymptotically converges the partition function, and the error of the logarithm of the partition function is much smaller than the algorithm using the Wang–Landau recursive scheme. The method is also applied to offlattice modelproteins, the AB models, in which cases many low energy states are found in different models.

Statistical mechanical theory for nonequilibrium systems. IX. Stochastic molecular dynamics
View Description Hide DescriptionThe general form for the probability density and for the transition probability of a nonequilibrium system is given. Maximization of the latter gives a generalized fluctuationdissipation theorem by providing a molecular basis for Langevin’s friction force that avoids continuum hydrodynamics. The result shows that the friction coefficient must be proportional to the variance of the stochastic equations of motion. Setting the variance to zero but keeping the friction coefficient nonzero reduces the theory to a Hoover thermostat without explicit constraint, although such a limit violates the physical requirement of proportionality between the dissipation and the fluctuation. A stochastic molecular dynamics algorithm is developed for both equilibrium and nonequilibrium systems, which is tested for steady heat flow and for a timevarying, driven Brownian particle.

On the accurate calculation of polarizabilities and second hyperpolarizabilities of polyacetylene oligomer chains using the CAMB3LYP density functional
View Description Hide DescriptionThe polarizability and second hyperpolarizability of polyacetylene oligomer chains of increasing size up to were investigated by means of the Coulombattenuating method (CAMB3LYP) using response theory. It was found that this longrange corrected density functional removes to large parts the overestimation observed for standard methods and in many cases provides results close to those of coupled cluster calculations. A direct comparison to experimentally observed dynamic hyperpolarizabilities is made to estimate the accuracy of the method. A basis set study revealed a noticeable contribution of diffuse orbitals to the hyperpolarizability also for larger oligomers. Furthermore, CAMB3LYP is also confirmed to provide molecular geometries close to experimentally observed structures, especially for longer chain lengths.

Selfassembly of nanocomponents into composite structures: Derivation and simulation of Langevin equations
View Description Hide DescriptionThe kinetics of the selfassembly of nanocomponents into a virus, nanocapsule, or other compositestructure is analyzed via a multiscale approach. The objective is to achieve predictability and to preserve key atomicscale features that underlie the formation and stability of the compositestructures. We start with an allatom description, the Liouville equation, and the order parameters characterizing nanoscale features of the system. An equation of Smoluchowski type for the stochastic dynamics of the order parameters is derived from the Liouville equation via a multiscale perturbation technique. The selfassembly of compositestructures from nanocomponents with internal atomic structure is analyzed and growth rates are derived. Applications include the assembly of a viral capsid from capsomers, a ribosome from its major subunits, and composite materials from fibers and nanoparticles. Our approach overcomes errors in other coarsegraining methods, which neglect the influence of the nanoscale configuration on the atomistic fluctuations. We account for the effect of order parameters on the statistics of the atomistic fluctuations, which contribute to the entropic and average forces driving order parameter evolution. This approach enables an efficient algorithm for computer simulation of selfassembly, whereas other methods severely limit the timestep due to the separation of diffusional and complexing characteristic times. Given that our approach does not require recalibration with each new application, it provides a way to estimate assembly rates and thereby facilitate the discovery of selfassembly pathways and kinetic deadend structures.
 Gas Phase Dynamics and Structure: Spectroscopy, Molecular Interactions, Scattering, and Photochemistry

Photodissociation dynamics of and excited by 527 nm photons
View Description Hide DescriptionThe photofragmentation dynamics of and clusters has been investigated at a 527 nm wavelength (2.35 eV) using a setup that allows simultaneous detection of the ionic and neutral fragments in a coincidence experiment. Measurement of positions and times of flight enables in principle a complete description of the fragmentation dynamics. The photofragmentation dynamics of clusters is similar to that of with, in addition, the ejection of a third fragment that can be neutral or ionized via a resonant electron capture. This is attributed to the triangular geometry of the ion.

VacuumUV negative photoion spectroscopy of , , and
View Description Hide DescriptionUsing synchrotron radiation, negative ions are detected by mass spectrometry following vacuumUV photoexcitation of trifluorochloromethane , trifluorobromomethane , and trifluoroiodomethane . The anions , , , , , , and are observed from all three molecules, where , Br, or I, and their ion yields recorded in the range of 8–35 eV. With the exception of and , the anions observed show a linear dependence of signal with pressure, showing that they arise from unimolecular ionpair dissociation.Dissociative electron attachment, following photoionization of and as the source of lowenergy electrons, is shown to dominate the observed and signals, respectively. Cross sections for ionpair formation are put onto an absolute scale by calibrating the signal strengths with those of from both and . These anion cross sections are normalized to vacuumUV absorption cross sections, where available, and the resulting quantum yields are reported. Anion appearance energies are used to calculate upper limits to 298 K bond dissociation energies for , which are consistent with literature values. We report new data for and . No ionpair formation is observed below the ionization energy of the parent molecule for and , and only weak signals (in both and ) are detected for . These observations suggest that neutral photodissociation is the dominant exit channel to Rydberg statephotoexcitation at these lower energies.

Timedependent wave packet and quasiclassical trajectory study of the reaction at the statetostate level
View Description Hide DescriptionThe first calculations of statetostate reaction probabilities and product stateresolved integral cross sections at selected collision energies (0.05, 0.1, 0.5, and 1.0 eV) for the title reaction on the ab initio potential energy surface of [Zanchet et al.J. Phys. Chem. A110, 12017 (2006)] with the OH reagent in selected rovibrational states have been carried out by means of the real wave packet (RWP) and quasiclassical trajectory(QCT) methods. Stateselected total reaction probabilities have been calculated for total angular momentum in a broad range of collision energies. Integral cross sections and statespecific rate coefficients have been obtained from the corresponding RWP reaction probabilities for initially selected rovibrational states by means of a capture model. The calculated RWP and QCT stateselected rate coefficients are practically temperature independent. Both RWP and QCTreaction probabilities, integral cross sections, and rate coefficients are almost independent of the initial rotational excitation. The RWP results are found to be in an overall good agreement with the corresponding QCT results. The present results have been compared with earlier wave packet calculations carried out on the same potential energy surface.

Highresolution spectroscopy of weak and shortlived bands of the transition of naphthalene
View Description Hide DescriptionRotationally resolved highresolution spectra and fluorescence decay curves have been observed for weak and shortlived vibronic bands of the transition of naphthalene. Fluorescence lifetime of the vibronic band with an excess energy of ( band) is remarkably shorter than that of other bands. Zeeman splitting of rotational lines is very small, so that the main radiationless process is not intersystem crossing to the triplet state but internal conversion to the ground state. The lifetime is thought to be governed by the strength of vibronic coupling between vibrational levels of the and states. As for the band, energy shifts were found in only a few rotational levels although the excess energy was higher than the threshold of intramolecular vibrational redistribution. We conclude that all of the rotational levels are mixed with other vibrational levels. The band spectrum is fairly complicated with numerous rotational lines, which is attributed to strong vibronic coupling with the state.

Theoretical study of complexes and transport of through RG
View Description Hide DescriptionWe present highlevel ab initio potential energy curves for barium cations and dications interacting with RG atoms . These potentials are employed to derive spectroscopic parameters for the and complexes, and also to derive the transport coefficients for and moving through a bath of the rare gas. The results are compared to the limited experimental data, which generally show reasonable agreement. We identify a large change in binding energy going from and to , which is not present in , and show that this is due to significant dispersion interactions in .