Volume 133, Issue 3, 21 July 2010

Developed for complex systems undergoing rare events involving many (meta)stable states, the multiple state transition path sampling aims to sample from an extended path ensemble including all possible trajectories between any pair of (meta)stable states. The key issue for an efficient sampling of the path space in this extended ensemble is sufficient switching between different types of trajectories. When some transitions are much more likely than others the collective sampling of the different path types can become difficult. Here we introduce a Wang–Landau based biasing approach to improve the sampling. We find that the biasing of the multiple state path ensemble does not influence the switching behavior, but does improve the sampling and thus the quality of the individual path ensembles.
 ARTICLES

 Theoretical Methods and Algorithms

On the efficiency of biased sampling of the multiple state path ensemble
View Description Hide DescriptionDeveloped for complex systems undergoing rare events involving many (meta)stable states, the multiple state transition path sampling aims to sample from an extended path ensemble including all possible trajectories between any pair of (meta)stable states. The key issue for an efficient sampling of the path space in this extended ensemble is sufficient switching between different types of trajectories. When some transitions are much more likely than others the collective sampling of the different path types can become difficult. Here we introduce a Wang–Landau based biasing approach to improve the sampling. We find that the biasing of the multiple state path ensemble does not influence the switching behavior, but does improve the sampling and thus the quality of the individual path ensembles.

Calculation of IR frequencies and intensities in electrical and mechanical anharmonicity approximations: Application to small water clusters
View Description Hide DescriptionWe present a method for automatic computation of infrared (IR) intensities using parallel variational multiple window configuration interactionwave functions ( algorithm). Inclusion of both mechanical and electrical anharmonic effects permits fundamental vibrational frequencies, including combinations and overtones, to be assigned. We use these developments to interpret the nearIR (NIR) and midIR (MIR) spectra of individual water clusters . Cyclic and linear systems are studied to provide equivalent reference theoretical data to investigate the structure of water as a function of density using NIR and MIR experimental spectra. Various density functional theory methods for generating the potential energy surface have been compared to reference results obtained at the CCSD(T) level [X. Huang et al., J. Chem. Phys.128, 034312 (2008)]. For cyclic clusters, the IR intensities and frequencies obtained using B1LYP/ccpVTZ are found to be in very good agreement with the available experimental values and of the same orders of magnitude as the reference theoretical values. These data are completed by the vibrational study of linear systems.

On the cluster composition of supercritical water combining molecular modeling and vibrational spectroscopic data
View Description Hide DescriptionThe present study is aimed at a detailed analysis of supercritical water structure based on the combination of experimental vibrational spectra as well as molecular modeling calculations of isolated water clusters. We propose an equilibrium cluster composition model where supercritical water is considered as an ideal mixture of small water clusters at the chemical equilibrium and the vibrational spectra are expected to result from the superposition of the spectra of the individual clusters, Thus, it was possible to extract from the decomposition of the midinfrared spectra the evolution of the partition of clusters in supercritical water as a function of density. The cluster composition predicted by this model was found to be quantitatively consistent with the near infrared and Raman spectra of supercritical wateranalyzed using the same procedure. We emphasize that such methodology could be applied to determine the portion of cluster in water in a wider thermodynamic range as well as in more complex aqueous supercritical solutions.

Monte Carlo simulations of single and multistep enzymecatalyzed reaction sequences: Effects of diffusion, cell size, enzyme fluctuations, colocalization, and segregation
View Description Hide DescriptionFollowing the discovery of slow fluctuations in the catalytic activity of an enzyme in singlemolecule experiments, it has been shown that the classical Michaelis–Menten (MM) equation relating the average enzymatic velocity and the substrate concentration may hold even for slowly fluctuating enzymes. In many cases, the average velocity is that given by the MM equation with timeaveraged values of the fluctuating rate constants and the effect of enzyme fluctuations is simply averaged out. The situation is quite different for a sequence of reactions. For colocalization of a pair of enzymes in a sequence to be effective in promoting reaction, the second must be active when the first is active or soon after. If the enzymes are slowly varying and only rarely active, the product of the first reaction may diffuse away before the second enzyme is active, and colocalization may have little value. Even for singlestep reactions the interplay of reaction and diffusion with enzyme fluctuations leads to added complexities, but for multistep reactions the interplay of reaction and diffusion, cell size, compartmentalization, enzyme fluctuations, colocalization, and segregation is far more complex than for singlestep reactions. In this paper, we report the use of stochastic simulations at the level of whole cells to explore, understand, and predict the behavior of single and multistep enzymecatalyzed reaction systems exhibiting some of these complexities. Results for singlestep reactions confirm several earlier observations by others. The MM relationship, with altered constants, is found to hold for singlestep reactions slowed by diffusion. For singlestep reactions, the distribution of enzymes in a regular grid is slightly more effective than a random distribution. Fluctuations of enzyme activity, with average activity fixed, have no observed effects for simple singlestep reactions slowed by diffusion. Twostep sequential reactions are seen to be slowed by segregation of the enzymes for each step, and results of the calculations suggest limits for cell size. Colocalization of enzymes for a twostep sequence is seen to promote reaction, and rates fall rapidly with increasing distance between enzymes. Low frequency fluctuations of the activities of colocalized enzymes, with average activities fixed, can greatly reduce reaction rates for sequential reactions.

First passage time distribution in stochastic processes with moving and static absorbing boundaries with application to biological rupture experiments
View Description Hide DescriptionWe develop and investigate an integral equation connecting the first passage time distribution of a stochastic process in the presence of an absorbing boundary condition and the corresponding Green’s function in the absence of the absorbing boundary. Analytical solutions to the integral equations are obtained for three diffusion processes in timeindependent potentials which have been previously investigated by other methods. The integral equation provides an alternative way to analytically solve the three diffusioncontrolled reactive processes. In order to help analyze biological rupture experiments, we further investigate the numerical solutions of the integral equation for a diffusion process in a timedependent potential. Our numerical procedure, based on the exact integral equation, avoids the adiabatic approximation used in previous analytical theories and is useful for fitting the rupture force distribution data from singlemolecule pulling experiments or molecular dynamics simulation data, especially at larger pulling speeds, larger cantilever spring constants, and smaller reaction rates. Stochastic simulation results confirm the validity of our numerical procedure. We suggest combining a previous analytical theory with our integral equation approach to analyze the kinetics of force induced rupture of biomacromolecules.

Logarithm secondorder manybody perturbation method for extended systems
View Description Hide DescriptionWe propose progressive downsampling of wave vectors in the Brillouin zone integrations occurring in the secondorder manybody or Møller–Plesset perturbation (MP2) method for extended systems of onedimensional periodicity. Higherlying unoccupied and lowerlying occupied Bloch orbitals are subject to downsampling by an exponentially increasing factor (with base ), making the total number of Bloch orbitals included in the MP2 lattice sums grow only logarithmically with respect to the number of basis functions per unit cell. Unlike the downsampling scheme proposed earlier, this scheme reduces the scaling of the computational cost and thus achieves a greater speedup as the unit cell size increases. Correct band indexing is essential for accuracy. Twoelectron integrals entering the MP2 energy and quasiparticle energy expressions must be multiplied by quadrature weights that are a function of the energy bands involved, and an algorithm to compute the weights is proposed. A combined use of the and schemes can speedup the calculation of polyacetylene typically by a factor of 20 with an error in the correlation energy within a few percent relative to the conventional calculation. Similar combinations can reproduce the MP2 quasiparticle energy bands accurately at a fraction of the usual computational cost.

A new basis set for molecular bending degrees of freedom
View Description Hide DescriptionWe present a new basis set as an alternative to Legendre polynomials for the variational treatment of bending vibrational degrees of freedom in order to highly reduce the number of basis functions. This basis set is inspired from the harmonic oscillator eigenfunctions but is defined for a bending angle in the range . The aim is to bring the basis functions closer to the final (ro)vibronic wave functions nature. Our methodology is extended to complicated potential energy surfaces, such as quasilinearity or multiequilibrium geometries, by using several free parameters in the basis functions. These parameters allow several density maxima, linear or not, around which the basis functions will be mainly located. Divergences at linearity in integral computations are resolved as generalized Legendre polynomials. All integral computations required for the evaluation of molecular Hamiltonian matrix elements are given for both discrete variable representation and finite basis representation. Convergence tests for the low energy vibronic states of , , and HCCS are presented.

Excited state geometry of photoactive yellow protein chromophore: A combined conductorlike polarizable continuum model and timedependent density functional study
View Description Hide DescriptionAnalytic gradient of the combined conductorlike polarizable continuum model (CPCM) and timedependent density functional theory method is derived and implemented. Due to the use of the fixed points with variable areas tessellation scheme, the excited state potential energy surfaces (PESs) are rigorously continuous and smooth. The CPCM/TDB3LYP method is used to study an analog of the photoactive yellow protein chromophore, anionic thiomethyl coumaric acid . Although CPCM/TDB3LYP method may not be accurate in predicting solvent effect on vertical excitation of , it may be used to predict redshiftings of emission maxima relative to absorption maxima with an accuracy of . We also found that the excited trans tends to form a single bond twisted structure in the gas phase but a double bond twisted structure in aqueous solution. The TDB3LYP minimum energyisomerization pathway shows a barrier of 3.6 kcal/mol in aqueous solution and 5.2 kcal/mol in the gas phase. The gas phase double bond twisted structure is trapped in a well of the excited state PES, with a depth of (0.88 eV), in good agreement with an experimental value of .

Firstprinciples theories for anharmonic lattice vibrations
View Description Hide DescriptionSizeextensive generalizations of the vibrational selfconsistent field (VSCF), vibrational Møller–Plesset perturbation (VMP), and vibrational coupledcluster (VCC) methods are made to anharmonic lattice vibrations of extended periodic systems on the basis of a quartic force field (QFF) in delocalized normal coordinates. Copious terms in the formalisms of VSCF that have nonphysical size dependence are identified algebraically and eliminated, leading to compact and strictly sizeextensive equations. This “quartic” VSCF method thus defined has no contributions from cubic force constants and alters only the transition energies of the underlying harmonicoscillator reference from a subset of quartic force constants. It also provides a way to evaluate an anharmonic correction to the lattice structure due to cubic force constants of a certain type. The secondorder VMP and VCC methods in the QFF based on the reference are shown to account for anharmonic effects due to all cubic and quartic force constants in a sizeextensive fashion. These methods can be readily extended to a higherorder truncated Taylor expansion of a potential energy surface in normal coordinates. An algebraic proof of the lack of sizeextensivity in the vibrational configurationinteraction method is also presented.

Firstprinciples calculations on anharmonic vibrational frequencies of polyethylene and polyacetylene in the approximation
View Description Hide DescriptionThe frequencies of the infrared and/or Ramanactive vibrations of polyethylene and polyacetylene are computed by taking account of the anharmonicity in the potential energy surfaces (PESs) and the resulting phononphonon couplings explicitly. The electronic part of the calculations is based on Gaussianbasisset crystalline orbital theory at the Hartree–Fock and secondorder Møller–Plesset (MP2) perturbation levels, providing one, two, and/or threedimensional slices of the PES (namely, using the socalled mode coupling approximation with ), which are in turn expanded in the fourthorder Taylor series with respect to the normal coordinates. The vibrational part uses the vibrational selfconsistent field, vibrational MP2, and vibrational truncated configurationinteraction (VCI) methods within the approximation, which amounts to including only phonons. It is shown that accounting for both electron correlation and anharmonicity is essential in achieving good agreement (the mean and maximum absolute deviations less than 50 and , respectively, for polyethylene and polyacetylene) between computed and observed frequencies. The corresponding values for the calculations including only one of such effects are in excess of 120 and , respectively. The VCI calculations also reproduce semiquantitatively the frequency separation and intensity ratio of the Fermi doublet involving the fundamental and first overtone in polyethylene.

Basis set effects on the hyperpolarizability of : Gaussiantype orbitals, numerical basis sets and realspace grids
View Description Hide DescriptionCalculations of the hyperpolarizability are typically much more difficult to converge with basis set size than the linear polarizability. In order to understand these convergence issues and hence obtain accurate ab initio values, we compare calculations of the static hyperpolarizability of the gasphase chloroform molecule using three different kinds of basis sets: Gaussiantype orbitals, numerical basis sets, and realspace grids. Although all of these methods can yield similar results, surprisingly large, diffuse basis sets are needed to achieve convergence to comparable values. These results are interpreted in terms of local polarizability and hyperpolarizability densities. We find that the hyperpolarizability is very sensitive to the molecular structure, and we also assess the significance of vibrational contributions and frequency dispersion.

Isomerization of nitrosomethane to formaldoxime: Energies, geometries, and frequencies from the parametric variational twoelectron reduceddensitymatrix method
View Description Hide DescriptionThe isomerization of nitrosomethane to transformaldoxime is treated with the parametric variational twoelectron reduceddensitymatrix (2RDM) method. In the parametric 2RDM method, the groundstateenergy is minimized with respect to a 2RDM that is parameterized to be both size extensive and nearly representable. The calculations were performed with an efficient version of the 2RDM method that we developed as an extension of the PSI3ab initio package. Details of the implementation, which scales like configuration interaction with single and double excitations, are provided as well as a comparison of two optimization algorithms for minimizing the energy functional. The conversion of nitrosomethane to transformaldoxime can occur by one of two pathways: (i) a 1,3sigmatropic hydrogen shift or (ii) two successive 1,2sigmatropic hydrogen shifts. The parametric 2RDM method predicts that the reaction channel involving two sequential 1,2shifts is about 10 kcal/mol more favorable than the channel with a single 1,3shift, which is consistent with calculations from other ab initio methods. We computed geometric parameters and harmonic frequencies for each stationary point on the reaction surfaces. Transitionstate energies, geometries, and frequencies from the 2RDM method are often more accurate than those from traditional wave function methods of a similar computational cost. Although electronicstructure methods generally agree that the 1,2shift is more efficient, the energy ordering of the reactant nitrosomethane and the 1,2shift intermediate formaldonitrone is unresolved in the literature. With an extrapolation to the completebasisset limit the parametric 2RDM method predicts formaldonitrone to be very slightly more stable than nitrosomethane.

Assigning quantum labels to variationally computed rotationalvibrational eigenstates of polyatomic molecules
View Description Hide DescriptionA procedure is investigated for assigning physically transparent, approximate vibrational and rotational quantum labels to variationally computed eigenstates. Pure vibrational wave functions are analyzed by means of normalmode decomposition (NMD) tables constructed from overlap integrals with respect to separable harmonic oscillator basis functions. Complementary rotational labels are determined from rigidrotor decomposition (RRD) tables formed by projecting rotationalvibrational wave functions onto products of symmetrized rigidrotor basis functions and previously computed vibrational eigenstates. Variational results for , HNCO, transHCOD, NCCO, and are presented to demonstrate the NMD and RRD schemes. The NMD analysis highlights several resonances at low energies that cause strong mixing and cloud the assignment of fundamental vibrations, even in such simple molecules. As the vibrational energy increases, the NMD scheme documents and quantifies the breakdown of the normalmode model. The RRD procedure proves effective in providing unambiguous rotational assignments for the chosen test molecules up to moderate values.

Distance and angular holonomic constraints in molecular simulations
View Description Hide DescriptionFinding the energy minima of systems with constraints is a challenging problem. We develop a minimization method based on the projection operator technique to enforce distance and angle constraints in minimization and reactionpath dynamics. The application of the projection operator alone does not maintain the constraints, i.e., they are slightly violated. Therefore, we use the SHAKEmethodology to enforce the constraints after each minimization step. We have extended SHAKE for bend angles and introduce SHAKE and SHAKE to constrain dihedral and outofplane angles, respectively. Two case studies are presented: (1) A mode analysis of unitedatom butane with various internal degrees of freedom kept frozen and (2) the minimization of chromene at a fixed approach toward the catalytic site of a (salen)Mn. The obtained information on energetics can be used to explain why specific enantioselectivity is observed. Previous minimization methods work for the free molecular case, but fail when molecules are tightly confined.

Quantum corrected Langevin dynamics for adsorbates on metal surfaces interacting with hot electrons
View Description Hide DescriptionWe investigate the importance of including quantized initial conditions in Langevin dynamics for adsorbates interacting with a thermal reservoir of electrons. For quadratic potentials the time evolution is exactly described by a classical Langevin equation and it is shown how to rigorously obtain quantum mechanical probabilities from the classical phase space distributions resulting from the dynamics. At short time scales, classical and quasiclassical initial conditions lead to wrong results and only correctly quantized initial conditions give a close agreement with an inherently quantum mechanical master equation approach. With CO on Cu(100) as an example, we demonstrate the effect for a system with ab initio frictional tensor and potential energy surfaces and show that quantizing the initial conditions can have a large impact on both the desorption probability and the distribution of molecular vibrational states.

QuasiNewton parallel geometry optimization methods
View Description Hide DescriptionAlgorithms for parallel unconstrained minimization of molecular systems are examined. The overall framework of minimization is the same except for the choice of directions for updating the quasiNewton Hessian. Ideally these directions are chosen so the updated Hessian gives steps that are same as using the Newton method. Three approaches to determine the directions for updating are presented: the straightforward approach of simply cycling through the Cartesian unit vectors (finite difference), a concurrent set of minimizations, and the Lanczos method. We show the importance of using preconditioning and a multiple secant update in these approaches. For the Lanczos algorithm, an initial set of directions is required to start the method, and a number of possibilities are explored. To test the methods we used the standard 50dimensional analytic Rosenbrock function. Results are also reported for the histidine dipeptide, the isoleucine tripeptide, and cyclic adenosine monophosphate. All of these systems show a significant speedup with the number of processors up to about eight processors.

Optimizing working parameters of the smooth particle mesh Ewald algorithm in terms of accuracy and efficiency
View Description Hide DescriptionWe construct an accurate error estimate for the root mean square force error of the smooth particle mesh Ewald (SPME) algorithm, which is often used for molecular dynamics simulations, where charge configurations under periodic boundary conditions require considerable amounts of CPU time. The error estimates are provided for the  as well as analytical force differentiation schemes, and their validity is tested for a random homogeneous sample system. Finally, we demonstrate that it is possible to straightforwardly and precisely determine the SPME parameters via the error estimates prior to the simulation for a predetermined accuracy. This can save precious computer and user time and allows an easy choice of a suitable parameter set for nearly optimal speed.

A harmonic transition state approximation for the duration of reactive events in complex molecular rearrangements
View Description Hide DescriptionMotivated by recent experimental efforts to measure the time a molecular system spends in transit between the reactants and the products of a chemical reaction, here we study the properties of the distribution of such transit times for the case of conservative dynamics on a multidimensional energy landscape. Unlike reaction rates, transit times are not invariant with respect to the order parameter (a.k.a. the experimental signal) used to monitor the progress of a chemical reaction. Nevertheless, such order parameter dependence turns out to be relatively weak. Moreover, for several model systems we find that the probability distribution of transit times can be estimated analytically, with reasonable accuracy, by assuming that the order parameter coincides with the direction of the unstable normal mode at the transition state. Although this approximation tends to overestimate the actual mean transit time measured using other order parameters, it yields asymptotically correct longtime behavior of the transit time distribution, which is independent of the order parameter.

Maximum caliber inference of nonequilibrium processes
View Description Hide DescriptionThirty years ago, Jaynes suggested a general theoretical approach to nonequilibrium statistical mechanics, called maximum caliber (MaxCal) [Annu. Rev. Phys. Chem.31, 579 (1980)]. MaxCal is a variational principle for dynamics in the same spirit that maximum entropy is a variational principle for equilibrium statistical mechanics. Motivated by the success of maximum entropyinference methods for equilibrium problems, in this work the MaxCal formulation is applied to the inference of nonequilibrium processes. That is, given some timedependent observables of a dynamical process, one constructs a model that reproduces these input data and moreover, predicts the underlying dynamics of the system. For example, the observables could be some timeresolved measurements of the folding of a protein, which are described by a fewstate model of the free energy landscape of the system. MaxCal then calculates the probabilities of an ensemble of trajectories such that on average the data are reproduced. From this probability distribution, any dynamical quantity of the system can be calculated, including population probabilities, fluxes, or waiting time distributions. After briefly reviewing the formalism, the practical numerical implementation of MaxCal in the case of an inference problem is discussed. Adopting various fewstate models of increasing complexity, it is demonstrated that the MaxCal principle indeed works as a practical method of inference: The scheme is fairly robust and yields correct results as long as the input data are sufficient. As the method is unbiased and general, it can deal with any kind of time dependency such as oscillatory transients and multitime decays.

Calculation of anharmonic OH phonon dispersion curves for the crystal
View Description Hide DescriptionAnharmonic OH phonon dispersion curves have been calculated for the crystal. A crystal Hamiltonian was set up for the vibrational problem, where the coordinates consists of the bond lengths of two hydroxide ions in the central unit cell. Its twodimensional potential energy surface was constructed from first principle calculations within the density functional theory approximation. Dispersion curves were calculated by diagonalizing the Hamiltonian in a basis of singly excited crystal functions. The single particle functions used to construct the crystal states were taken from a Morse oscillator basis set. These well chosen functions made it possible to restrict calculations to include only very few functions, which greatly contributed to a transparent presentation of the underlying theory. All calculations could be done analytically except for the calculation of a few integrals. We have compared our results with those of a series of harmonic lattice dynamics calculations and have found that the anharmonicity shifts the IR and Raman dispersion curves downward appreciably and slightly changes the energy differences between both curves. From an analysis of the harmonic results we conclude that incorporating the coupling between OH stretching motion and the motion of their centers of mass will appreciably change the overall features of the dispersion curves. Extension of the anharmonic model along these lines will cause no problem to the theoretical approach presented in this paper.