Volume 137, Issue 13, 07 October 2012

We introduce a sequencedependent parametrization for a coarsegrained DNA model [T. E. Ouldridge, A. A. Louis, and J. P. K. Doye, J. Chem. Phys.134, 085101 (2011)] originally designed to reproduce the properties of DNA molecules with average sequences. The new parametrization introduces sequencedependent stacking and basepairing interaction strengths chosen to reproduce the melting temperatures of short duplexes. By developing a histogram reweighting technique, we are able to fit our parameters to the melting temperatures of thousands of sequences. To demonstrate the flexibility of the model, we study the effects of sequence on: (a) the heterogeneous stacking transition of single strands, (b) the tendency of a duplex to fray at its melting point, (c) the effects of stacking strength in the loop on the melting temperature of hairpins, (d) the forceextension properties of single strands, and (e) the structure of a kissingloop complex. Where possible, we compare our results with experimental data and find a good agreement. A simulation code called oxDNA, implementing our model, is available as a free software.
 COMMUNICATIONS


Communication: Explicitly correlated fourcomponent relativistic secondorder MøllerPlesset perturbation theory
View Description Hide DescriptionWe propose explicitly correlated Ansatz for fourcomponent relativistic methods within the framework of the nopair approximation. Kinetically balanced geminal basis is derived to satisfy the cusp conditions in the nonrelativistic limit based on the LévyLeblendlike equation. Relativistic variants of strongorthogonality projection operator (Ansätze 2α and 2β) suitable for practical calculations are introduced by exploiting the orthogonal complement of the largecomponent basis. A pilot implementation is performed for the second order MøllerPlesset perturbation theory.

Communication: On the isotope anomaly of nuclear quadrupole coupling in molecules
View Description Hide DescriptionThe dependence of the nuclear quadrupole coupling constants (NQCC) on the interaction between electrons and a nucleus of finite size is theoretically analyzed. A deviation of the ratio of the NQCCs obtained from two different isotopomers of a molecule from the ratio of the corresponding bare nuclear electric quadrupole moments, known as quadrupole anomaly, is interpreted in terms of the logarithmic derivatives of the electric field gradient at the nuclear site with respect to the nuclear charge radius. Quantum chemical calculations based on a Diracexact relativistic methodology suggest that the effect of the changing size of the Au nucleus in different isotopomers can be observed for Aucontaining molecules, for which the predicted quadrupole anomaly reaches values of the order of 0.1%. This is experimentally detectable and provides an insight into the charge distribution of nonspherical nuclei.

Communication: Restoring full size extensivity in internally contracted multireference coupled cluster theory
View Description Hide DescriptionThe reason for the lack of size extensivity in the valence space in current implementations of internally contracted multireference coupled clustertheories is the procedure used to eliminate redundant components from the cluster operator. We present a simple way to restore full size extensivity by performing this critical step in a basis of excitation operators that are normal ordered with respect to the multiconfigurational reference function.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

The Ehrenfest force field: Topology and consequences for the definition of an atom in a molecule
View Description Hide DescriptionThe Ehrenfest force is the force acting on the electrons in a molecule due to the presence of the other electrons and the nuclei. There is an associated force field in threedimensional space that is obtained by the integration of the corresponding Hermitian quantum force operator over the spin coordinates of all of the electrons and the space coordinates of all of the electrons but one. This paper analyzes the topology induced by this vector field and its consequences for the definition of molecular structure and of an atom in a molecule. Its phase portrait reveals: that the nuclei are attractors of the Ehrenfest force, the existence of separatrices yielding a dense partitioning of threedimensional space into disjoint regions, and field lines connecting the attractors through these separatrices. From the numerical point of view, when the Ehrenfest force field is obtained as minus the divergence of the kinetic stress tensor, the induced topology was found to be highly sensitive to choice of Gaussian basis sets at long range. Even the use of large split valence and highly uncontracted basis sets can yield spurious critical points that may alter the number of attraction basins. Nevertheless, at short distances from the nuclei, in general, the partitioning of threedimensional space with the Ehrenfest force field coincides with that induced by the gradient field of the electron density. However, exceptions are found in molecules where the electron density yields results in conflict with chemical intuition. In these cases, the molecular graphs of the Ehrenfest force field reveal the expected atomic connectivities. This discrepancy between the definition of an atom in a molecule between the two vector fields casts some doubts on the physical meaning of the integration of Ehrenfest forces over the basins of the electron density.

Increasing the applicability of density functional theory. III. Do consistent KohnSham density functional methods exist?
View Description Hide DescriptionThe concept of a “consistent,” KohnSham (KS) density functional theory(DFT) is discussed, where the functional is able to provide good total energies and its selfconsistent potential is such that the KS eigenvalues correspond to accurate approximations to the principal ionization potentials for the molecule. Today, none of the vast number of DFT approximations show this property. The one exception is the ab initiodft method built upon the optimized effective potential strategy for exchange and correlation. This qualifies as a DFT method because it represents the correlated density as a single determinant and by imposing that condition, generates local exchange and correlation operators which are used in selfconsistent solutions of the orbitals and eigenvalues. Such a “consistent” DFT shares many of the properties of the Dyson equation, but without its frequency dependence and associated complications. The relationship between ab initiodft based on MBPT2 functional and GW method is discussed. Ab initiodft provides a selfconsistent, frequency independent, effective independent particle alternative with a local correlation potential.

Density functional theory for molecular multiphoton ionization in the perturbative regime
View Description Hide DescriptionA general implementation of the lowest nonvanishing order perturbation theory for the calculation of molecular multiphoton ionization cross sections is proposed in the framework of density functional theory. Bound and scattering wave functions are expanded in a multicentric basis set and advantage is taken of the full molecular point group symmetry, thus enabling the application of the formalism to mediumsize molecules. Multiphoton ionization cross sections and angular asymmetry parameters have been calculated for the two and fourphoton ionization of the molecule, for linear and circular light polarizations. Both fixed and random orientations of the target molecule have been considered. To demonstrate the efficiency of the proposed methodology, the twophoton cross section and angular asymmetry parameters for the HOMO and HOMO1 orbital ionization of benzene are also presented.

A generalized Irving–Kirkwood formula for the calculation of stress in molecular dynamics models
View Description Hide DescriptionIn nonequilibrium molecular dynamics simulations,continuum mechanics quantities can be computed from the position and momentum of the particles based on the classical Irving–Kirkwood formalism. For practical purposes, the implementations of Irving–Kirkwood formulas often involve a spatial averaging using a smooth kernel function. The resulting formula for the stress has been known as Hardy stress. Usually results obtained this way still need to be further processed to reduce the fluctuation, e.g., by ensemble or time averaging. In this paper we extend Hardy's formulas by systematically incorporating both spatial and temporal averaging into the expression of continuum quantities. The derivation follows the Irving–Kirkwood formalism, and the average quantities still satisfy conservation laws in continuum mechanics. We will discuss the selection of kernel functions and present several numerical tests.

Energy conserving, linear scaling BornOppenheimer molecular dynamics
View Description Hide DescriptionBornOppenheimer molecular dynamics simulations with longterm conservation of the total energy and a computational cost that scales linearly with system size have been obtained simultaneously. Linear scaling with a low prefactor is achieved using density matrix purification with sparse matrix algebra and a numerical threshold on matrix elements. The extended Lagrangian BornOppenheimer molecular dynamics formalism [A. M. N. Niklasson, Phys. Rev. Lett.100, 123004 (2008)10.1103/PhysRevLett.100.123004] yields microcanonical trajectories with the approximate forces obtained from the linear scaling method that exhibit no systematic drift over hundreds of picoseconds and which are indistinguishable from trajectories computed using exact forces.

Structural optimization of molecular clusters with density functional theory combined with basin hopping
View Description Hide DescriptionIdentifying the energy minima of molecular clusters is a challenging problem. Traditionally, search algorithms such as simulated annealing, genetic algorithms, or basin hopping are usually used in conjunction with empirical force fields. We have implemented a basin hopping search algorithm combined with density functional theory to enable the optimization of molecular clusters without the need for empirical force fields. This approach can be applied to systems where empirical potentials are not available or may not be sufficiently accurate. We illustrate the effectiveness of the method with studies on water, methanol, and water + methanol clusters as well as protonated water and methanol clusters at the B3LYP+D/631+G* level of theory. A new lowest energy structure for H^{+}(H_{2}O)_{7} is predicted at the B3LYP+D/631+G* level. In all of the protonated mixed water and methanol clusters, we find that H^{+} prefers to combine with methanol rather than water in the lowestenergy structures.

Semiclassical evaluation of kinetic isotope effects in 13atomic system
View Description Hide DescriptionThe semiclassical instanton approach discussed by Kryvohuz [J. Chem. Phys.134, 114103 (2011)10.1063/1.3565425] is applied to calculate kinetic H/D isotope effect (KIE) of intramolecular hydrogen transfer in cis1,3pentadiene. All 33 vibrational degrees of freedom are treated quantum mechanically with semiclassical approximation. Nuclear quantum effects such as tunneling under the barrier and zeropoint energy are automatically incorporated in the theory, and are shown to be responsible for the observed appreciable kinetic isotope effect in cis1,3pentadiene. Over the barrier passage is also automatically included. Numerical calculations are performed on an empirical valence bond potential energy surface and compared with the previous experimental and theoretical studies. An estimation of heavyatom ^{12}C/^{13}C KIE in the same system is also provided and the factors contributing to it are discussed.

Efficient and accurate solver of the threedimensional screened and unscreened Poisson's equation with generic boundary conditions
View Description Hide DescriptionWe present an explicit solver of the threedimensional screened and unscreened Poisson's equation, which combines accuracy, computational efficiency, and versatility. The solver, based on a mixed planewave/interpolating scaling function representation, can deal with any kind of periodicity (along one, two, or three spatial axes) as well as with fully isolated boundary conditions. It can seamlessly accommodate a finite screening length, nonorthorhombic lattices, and charged systems. This approach is particularly advantageous because convergence is attained by simply refining the real space grid, namely without any adjustable parameter. At the same time, the numerical method features scaling of the computational cost (N being the number of grid points) very much like planewave methods. The methodology, validated on model systems, is tailored for leadingedge computer simulations of materials (including ab initioelectronic structure computations), but it might as well be beneficial for other research domains.

Comparison of some dispersioncorrected and traditional functionals with CCSD(T) and MP2 ab initio methods: Dispersion, induction, and basis set superposition error
View Description Hide DescriptionWe compare dispersion and induction interactions for noble gas dimers and for Ne, methane, and 2butyne with HF and LiF using a variety of functionals (including some specifically parameterized to evaluate dispersion interactions) with ab initio methods including CCSD(T) and MP2. We see that inductive interactions tend to enhance dispersion and may be accompanied by chargetransfer. We show that the functionals do not generally follow the expected trends in interaction energies, basis set superposition errors (BSSE), and interaction distances as a function of basis set size. The functionals parameterized to treat dispersion interactions often overestimate these interactions, sometimes by quite a lot, when compared to higher level calculations. Which functionals work best depends upon the examples chosen. The B3LYP and X3LYP functionals, which do not describe pure dispersion interactions, appear to describe dispersion mixed with induction about as accurately as those parametrized to treat dispersion. We observed significant differences in highlevel wavefunction calculations in a basis set larger than those used to generate the structures in many of the databases. We discuss the implications for highly parameterized functionals based on these databases, as well as the use of simple potential energy for fitting the parameters rather than experimentally determinable thermodynamic state functions that involve consideration of vibrational states.

Maintain rigid structures in Verlet based Cartesian molecular dynamics simulations
View Description Hide DescriptionAn algorithm is presented to maintain rigid structures in Verlet based Cartesian molecular dynamics (MD) simulations. After each unconstrained MD step, the coordinates of selected particles are corrected to maintain rigid structures through an iterative procedure of rotation matrix computation. This algorithm, named as SHAPE and implemented in CHARMM program suite, avoids the calculations of Lagrange multipliers, so that the complexity of computation does not increase with the number of particles in a rigid structure. The implementation of this algorithm does not require significant modification of propagation integrator, and can be plugged into any Cartesian based MD integration scheme. A unique feature of the SHAPE method is that it is interchangeable with SHAKE for any object that can be constrained as a rigid structure using multiple SHAKE constraints. Unlike SHAKE, the SHAPE method can be applied to large linear (with three or more centers) and planar (with four or more centers) rigid bodies. Numerical tests with four model systems including two proteins demonstrate that the accuracy and reliability of the SHAPE method are comparable to the SHAKE method, but with much more applicability and efficiency.

Improved coarsegraining of Markov state models via explicit consideration of statistical uncertainty
View Description Hide DescriptionMarkov state models (MSMs)–or discretetime master equation models–are a powerful way of modeling the structure and function of molecular systems like proteins. Unfortunately, MSMs with sufficiently many states to make a quantitative connection with experiments (often tens of thousands of states even for small systems) are generally too complicated to understand. Here, I present a Bayesian agglomerative clustering engine (BACE) for coarsegraining such Markovmodels, thereby reducing their complexity and making them more comprehensible. An important feature of this algorithm is its ability to explicitly account for statistical uncertainty in model parameters that arises from finite sampling. This advance builds on a number of recent works highlighting the importance of accounting for uncertainty in the analysis of MSMs and provides significant advantages over existing methods for coarsegraining Markov state models. The closedform expression I derive here for determining which states to merge is equivalent to the generalized JensenShannon divergence, an important measure from information theory that is related to the relative entropy. Therefore, the method has an appealing information theoretic interpretation in terms of minimizing information loss. The bottomup nature of the algorithm likely makes it particularly well suited for constructing mesoscale models. I also present an extremely efficient expression for Bayesian model comparison that can be used to identify the most meaningful levels of the hierarchy of models from BACE.

Finitetemperature electronic simulations without the BornOppenheimer constraint
View Description Hide DescriptionThe adiabatic approximation, typically assumed when performing standard BornOppenheimer (BO) molecular dynamics, can become unreliable at finite temperature, and specifically when the temperature is larger than the electronic energy gap between the ground state and the lowlying excited states. In this regime, relevant for many important chemical processes, the nonadiabatic couplings between the electronic energy states can produce finite temperature effects in several molecular properties, such as the geometry, the vibrational frequencies, the binding energy, and several chemical reactions. In this work, we introduce a novel finitetemperature nonadiabatic molecular dynamics based on a novel covariant formulation of the electronic partition function. In this framework, the nuclei are not constrained to move in a specific electronic potential energy surface. Then, by using a rigorous variational upper bound to the free energy, we are led to an approximate partition function that can be evaluated numerically. The method can be applied to any technique capable to provide an energy value over a given wave function ansatz depending on several variational parameters and atomic positions. In this work, we have applied the proposed method within a quantum Monte Carlo (QMC) scheme. In particular, we consider in this first application only classical ions, but we explicitly include an electronic correlation (Jastrow) term in the wave function, by extending in this way the standard variational QMC method, from ground state to finite temperature properties. We show that our approximation reduces correctly to the standard groundstate BornOppenheimer (gsBO) at zero temperature and to the correct high temperature limit. Moreover, at temperatures large enough, this method improves the upper bound of the free energy obtained with a single BO energy surface, since within our approach it is possible to estimate the electron entropy of a correlated ansatz in an efficient way. We test this new method on the simple hydrogen molecule, where at low temperature we recover the correct gsBO low temperature limit. Moreover, we show that the dissociation of the molecule is possible at a temperature much smaller than the one corresponding to the gsBO energy surface, in good agreement with experimental evidence. Several extensions of the proposed technique are also discussed, as for instance the inclusion of quantum effects for ions and the calculation of critical (magnetic, superconducting) temperatures.

Exploring constrained quantum control landscapes
View Description Hide DescriptionThe broad success of optimally controlling quantum systems with external fields has been attributed to the favorable topology of the underlying control landscape, where the landscape is the physical observable as a function of the controls. The control landscape can be shown to contain no suboptimal trapping extrema upon satisfaction of reasonable physical assumptions, but this topological analysis does not hold when significant constraints are placed on the control resources. This work employs simulations to explore the topology and features of the control landscape for purestate population transfer with a constrained class of control fields. The fields are parameterized in terms of a set of uniformly spaced spectral frequencies, with the associated phases acting as the controls. This restricted family of fields provides a simple illustration for assessing the impact of constraints upon seeking optimal control. Optimization results reveal that the minimum number of phase controls necessary to assure a high yield in the target state has a special dependence on the number of accessible energy levels in the quantum system, revealed from an analysis of the first and secondorder variation of the yield with respect to the controls. When an insufficient number of controls and/or a weak control fluence are employed, trapping extrema and saddle points are observed on the landscape. When the control resources are sufficiently flexible, solutions producing the globally maximal yield are found to form connected “level sets” of continuously variable control fields that preserve the yield. These optimal yield level sets are found to shrink to isolated points on the top of the landscape as the control field fluence is decreased, and further reduction of the fluence turns these points into suboptimal trapping extrema on the landscape. Although constrained control fields can come in many forms beyond the cases explored here, the behavior found in this paper is illustrative of the impacts that constraints can introduce.
 Advanced Experimental Techniques

Recoupling of chemical shift anisotropy by Rsymmetry sequences in magic angle spinning NMR spectroscopy
View Description Hide Description^{13}C and ^{15}N chemical shift (CS) interaction is a sensitive probe of structure and dynamics in a wide variety of biological and inorganic systems, and in the recent years several magic angle spinning NMR approaches have emerged for residuespecific measurements of chemical shift anisotropy (CSA) tensors in uniformly and sparsely enriched proteins. All of the currently existing methods are applicable to slow and moderate magic angle spinning (MAS) regime, i.e., MAS frequencies below 20 kHz. With the advent of fast and ultrafast MAS probes capable of spinning frequencies of 40–100 kHz, and with the superior resolution and sensitivity attained at such high frequencies, development of CSA recoupling techniques working under such conditions is necessary. In this work, we present a family of Rsymmetry based pulse sequences for recoupling of ^{13}C/^{15}N CSA interactions that work well in both natural abundance and isotopically enriched systems. We demonstrate that efficient recoupling of either firstrank (σ 1) or secondrank (σ 2) spatial components of CSA interaction is attained with appropriately chosen γencoded RN n ^{ v } symmetry sequences. The advantage of these γencoded RN n ^{ v }symmetry based CSA (RNCSA) recoupling schemes is that they are suitable for CSA recoupling under a wide range of MAS frequencies, including fast MAS regime. Comprehensive analysis of the recoupling properties of these RN n ^{ v } symmetry sequences reveals that the σ 1CSA recoupling symmetry sequences exhibit large scaling factors; however, the partial homonuclear dipolar Hamiltonian components are symmetry allowed, which makes this family of sequences suitable for CSA measurements in systems with weak homonuclear dipolar interactions. On the other hand, the γencoded symmetry sequences for σ 2CSA recoupling have smaller scaling factors but they efficiently suppress the homonuclear dipoledipole interactions. Therefore, the latter family of sequences is applicable for measurements of CSA parameters in systems with strong homonuclear dipolar couplings, such as uniformly^{13}C labeled biological solids. We demonstrate RNCSA NMR experiments and numerical simulations establishing the utility of this approach to the measurements of ^{13}C and ^{15}N CSA parameters in model compounds, [^{15}N]Nacetylvaline (NAV), [U^{13}C, ^{15}N]alanine, [U^{13}C,^{15}N]histidine, and present the application of this approach to [U^{13}C/^{15}N]Tyr labeled Cterminal domain of HIV1 CA protein.
 Atoms, Molecules, and Clusters

State and species selective energy flow in gas ensembles containing vibrationally excited O_{2}
View Description Hide DescriptionStatetostate, collisioninduced, energy transfer is followed to equilibrium through sequences of collision cycles in gas ensembles containing vibrationally excited oxygen molecules (v = 8 and 1) in several different atomic and molecular bath gases. Quantum state distributions for each of the constituent species are available at each stage of the ensemble's evolution and enable the dominant energy exchange mechanisms to be identified. Equilibration is generally a complex process that evolves through several phases of inter and intramolecular events, each with their characteristic response rate to collisions. The results suggest that single quantum state population loss rate constants, however precisely determined, may miss key features of the overall equilibration process.

Full dimensional quantummechanical simulations for the vibronic dynamics of difluorobenzene radical cation isomers using the multilayer multiconfiguration timedependent Hartree method
View Description Hide DescriptionFull dimensional multilayer multiconfiguration timedependent Hartree (MLMCTDH) calculations of the dynamics of the three difluorobenzene cationic isomers in five lowestlying doublet electronic states using the ab initio multistate multimode vibronic coupling Hamiltonian (MMVCH) model are carried out using the Heidelberg MCTDH package. The same dynamical problems, but treated with the MCTDH scheme and using a reduced dimensional ab initio MMVCH model, have been previously reported [S. Faraji, H.D. Meyer, and H. Köppel, “Multistate vibronic interactions in difluorobenzene radical cations. II Quantum dynamical simulations,” J. Chem. Phys.129, 074311 (2008)10.1063/1.2958918]. For easy comparison with the reduced dimensional results, 11D or 10D MLMCTDH calculations are also performed. Extensive MLMCTDH test calculations are performed to find appropriate MLMCTDH wavefunction structures (MLtrees), and the convergence of the MLMCTDH calculations are carefully checked to ensure accurate results. Based on the appropriate MLtrees, the photoelectron (PE) spectrum and the mass analyzed threshold ionization (MATI) spectrum are simulated, analyzed, and compared with corresponding experimental spectra. Because of its efficient simulation capability for large systems, MLMCTDH calculations save a considerable amount of central processing unit (CPU)time, even when a reduced dimensional MMVCH is used, i.e., the same reduced model as in the corresponding MCTDH calculations. Simulations of the experimental PE spectra by full dimensional MLMCTDH calculations reproduced main peaks, which originate from different electronic states. The agreement is improved as compared to the reduced dimensionality calculations. Unfortunately, the experimental PE spectra are not very well resolved. Therefore, we compare our calculations additionally with highly resolved MATI spectra, which, however, are only available for the state. Based on a series of MLMCTDH simulations with longer propagation time for , a number of vibrational modes, including fundamentals, their combinations, and overtones are simulated and assigned by comparing with the experimental assignments and the ab initio frequencies. Excellent correlation between the experimental and full dimensional MLMCTDH results show that MLMCTDH is accurate and very efficient and that the ab initio MMVCH model is very suitable for MLMCTDH calculations.

Al_{6}H_{18}: A baby crystal of γAlH_{3}
View Description Hide DescriptionUsing globalminima search methods based on the density functional theory calculations of (AlH_{3})_{ n } (n = 1–8) clusters, we show that the growth pattern of alanes for n ≥ 4 is dominated by structures containing hexacoordinated Al atoms. This is in contrast to the earlier studies where either linear or ring structures of AlH_{3} were predicted to be the preferred structures in which the Al atoms can have a maximum of fivefold coordination. Our calculations also reveal that the Al_{6}H_{18} cluster, with its hexacoordination of the Al atoms, resembles the unitcell of γ AlH_{3}, thus Al_{6}H_{18} is designated as the “baby crystal.” The fragmentation energies of the (AlH_{3})_{ n } (n = 2–8) along with the dimerization energies for even n clusters indicate an enhanced stability of the Al_{6}H_{18} cluster. Both covalent (hybridization) and ionic (charge) contribution to the bonding are the driving factors in stabilizing the isomers containing hexacoordinated Al atoms.