Volume 135, Issue 10, 14 September 2011
Index of content:

The structures and the electronic properties of two aluminumdoped boron clusters, AlB_{7} ^{−} and AlB_{8} ^{−}, were investigated using photoelectron spectroscopy and ab initio calculations. The photoelectron spectra of AlB_{7} ^{−} and AlB_{8} ^{−} are both broad, suggesting significant geometry changes between the ground states of the anions and the neutrals. Unbiased global minimum searches were carried out and the calculated vertical electron detachment energies were used to compare with the experimental data. We found that the Al atom does not simply replace a B atom in the parent B_{8} ^{−} and B_{9} ^{−} planar clusters in AlB_{7} ^{−} and AlB_{8} ^{−}. Instead, the global minima of the two dopedclusters are of umbrella shapes, featuring an Al atom interacting ionically with a hexagonal and heptagonal pyramidal B_{7} (C _{6v }) and B_{8} (C _{7v }) fragment, respectively. These unique umbrellatype structures are understood on the basis of the special stability of the quasiplanar B_{7} ^{3−} and planar B_{8} ^{2−} molecular wheels derived from double aromaticity.
 ARTICLES

 Theoretical Methods and Algorithms

NVU dynamics. I. Geodesic motion on the constantpotentialenergy hypersurface
View Description Hide DescriptionAn algorithm is derived for computer simulation of geodesics on the constantpotentialenergy hypersurface of a system of N classical particles. First, a basic timereversible geodesic algorithm is derived by discretizing the geodesic stationarity condition and implementing the constantpotentialenergy constraint via standard Lagrangian multipliers. The basic NVU algorithm is tested by singleprecision computer simulations of the LennardJones liquid. Excellent numerical stability is obtained if the force cutoff is smoothed and the two initial configurations have identical potential energy within machine precision. Nevertheless, just as for NVE algorithms, stabilizers are needed for very long runs in order to compensate for the accumulation of numerical errors that eventually lead to “entropic drift” of the potential energy towards higher values. A modification of the basic NVU algorithm is introduced that ensures potentialenergy and steplength conservation; centerofmass drift is also eliminated. Analytical arguments confirmed by simulations demonstrate that the modified NVU algorithm is absolutely stable. Finally, we present simulations showing that the NVU algorithm and the standard leapfrog NVE algorithm have identical radial distribution functions for the LennardJones liquid.

NVU dynamics. II. Comparing to four other dynamics
View Description Hide DescriptionIn the companion paper [T. S. Ingebrigtsen, S. Toxvaerd, O. J. Heilmann, T. B. Schrøder, and J. C. Dyre, “NVUdynamics. I. Geodesic motion on the constantpotentialenergy hypersurface,” J. Chem. Phys. (in press)] an algorithm was developed for tracing out a geodesic curve on the constantpotentialenergy hypersurface. Here, simulations of NVUdynamics are compared to results for four other dynamics, both deterministic and stochastic. First, NVUdynamics is compared to the standard energyconserving Newtonian NVEdynamics by simulations of the KobAndersen binary LennardJones liquid, its WCA version (i.e., with cutoff's at the pair potential minima), and the LennardJones Gaussian liquid. We find identical results for all quantities probed: radial distribution functions, incoherent intermediate scattering functions, and meansquare displacement as function of time. Arguments are presented for the equivalence of NVU and NVEdynamics in the thermodynamic limit; in particular, to leading order in 1/N these two dynamics give identical timeautocorrelation functions. In the final part of the paper, NVUdynamics is compared to Monte Carlodynamics, to a diffusive dynamics of smallstep random walks on the constantpotentialenergy hypersurface, and to NosHoover NVTdynamics. If time is scaled for the two stochastic dynamics to make singleparticle diffusion constants identical to that of NVEdynamics, the simulations show that all five dynamics are equivalent at low temperatures except at short times.

Quadratically convergent algorithm for orbital optimization in the orbitaloptimized coupledcluster doubles method and in orbitaloptimized secondorder MøllerPlesset perturbation theory
View Description Hide DescriptionUsing a Lagrangianbased approach, we present a more elegant derivation of the equations necessary for the variational optimization of the molecular orbitals (MOs) for the coupledcluster doubles (CCD) method and secondorder MøllerPlesset perturbation theory (MP2). These orbitaloptimized theories are referred to as OOCCD and OOMP2 (or simply “OD” and “OMP2” for short), respectively. We also present an improved algorithm for orbital optimization in these methods. Explicit equations for response density matrices, the MO gradient, and the MO Hessian are reported both in spinorbital and closedshell spinadapted forms. The NewtonRaphson algorithm is used for the optimization procedure using the MO gradient and Hessian. Further, orbital stability analyses are also carried out at correlated levels. The OD and OMP2 approaches are compared with the standard MP2, CCD, CCSD, and CCSD(T) methods. All these methods are applied to H_{2}O, three diatomics, and the molecule. Results demonstrate that the CCSD and OD methods give nearly identical results for H_{2}O and diatomics; however, in symmetrybreaking problems as exemplified by , the OD method provides better results for vibrational frequencies. The OD method has further advantages over CCSD: its analytic gradients are easier to compute since there is no need to solve the coupledperturbed equations for the orbital response, the computation of oneelectron properties are easier because there is no response contribution to the particle density matrices, the variational optimized orbitals can be readily extended to allow inactive orbitals, it avoids spurious secondorder poles in its response function, and its transition dipole moments are gauge invariant. The OMP2 has these same advantages over canonical MP2, making it promising for excited state properties via linear response theory. The quadratically convergent orbitaloptimization procedure converges quickly for OMP2, and provides molecular properties that are somewhat different than those of MP2 for most of the test cases considered (although they are similar for H_{2}O). Bond lengths are somewhat longer, and vibrational frequencies somewhat smaller, for OMP2 compared to MP2. In the difficult case of , results for several vibrational frequencies are significantly improved in going from MP2 to OMP2.

A molecular DebyeHückel theory and its applications to electrolyte solutions
View Description Hide DescriptionIn this report, a molecular DebyeHückel theory for ionic fluids is developed. Starting from the macroscopic Maxwell equations for bulk systems, the dispersion relation leads to a generalized DebyeHückel theory which is related to the dressed ion theory in the static case. Due to the multipole structure of dielectric function of ionic fluids, the electric potential around a single ion has a multiYukawa form. Given the dielectric function, the multiYukawa potential can be determined from our molecular DebyeHückel theory, hence, the electrostatic contributions to thermodynamic properties of ionic fluids can be obtained. Applications to binary as well as multicomponent primitive models of electrolyte solutions demonstrated the accuracy of our approach. More importantly, for electrolyte solution models with soft shortranged interactions, it is shown that the traditional perturbation theory can be extended to ionic fluids successfully just as the perturbation theory has been successfully used for shortranged systems.

Local pressure components and surface tension of spherical interfaces. Thermodynamic versus mechanical definitions. I. A mesoscale modeling of droplets
View Description Hide DescriptionWe report mesoscale simulations of spherical drops to investigate the surface tension and mechanical properties. The Monte Carlo simulations are performed with the multibody potential commonly used in the manybody dissipative particle dynamics simulations. We establish here the calculation of the local normal and transverse components of the pressure tensor via the perturbation volume within the thermodynamic route. The different profiles of these components are compared to those calculated using the mechanical approach. To complete the mesoscale modeling of drops, we investigate the curvature dependence of the surface tension in order to calculate the Tolman's length, which is found to be negative.

Enhanced sampling of particular degrees of freedom in molecular systems based on adiabatic decoupling and temperature or force scaling
View Description Hide DescriptionA method to enhance sampling of a small subset of N ^{ h } particular degrees of freedom of a system of N ^{ h } + N ^{ l } degrees of freedom is presented. It makes use of adiabatically decoupling these degrees of freedom by increasing their mass followed by either increasing their temperature or reducing their interaction or the force acting on them. The appropriate statisticalmechanical expressions for use of these methods in simulation studies are derived. As long as the subset of massincreased degrees of freedom is small compared to the total number of degrees of freedom of the system, sampling of this subset of degrees of freedom can be much enhanced at the cost of a slight perturbation of the configurational distribution. This is illustrated for a test system of 1000 SPC, simple point charge, water molecules at 300 K and a density of 997 kg m^{−3}. Various fractions N ^{ h }/(N ^{ h } + N ^{ l }) of water molecules were adiabatically decoupled to different degrees. The size of the diffusion coefficient of these decoupled water molecules was used as a measure for how much the sampling was enhanced and the average potential energy per water molecule was used as a measure of how much the configurational distribution of the system gets distorted. A variety of parameter values was investigated and it was found that for N ^{ h }/(N ^{ h } + N ^{ l }) ⩽ 0.1 the diffusion of the N ^{ h } molecules could be enhanced by factors up to 35 depending on the method, the ratio N ^{ h }/(N ^{ h } + N ^{ l }), the extent of adiabatic decoupling, and the temperature or force scaling factors, at the cost of a slight perturbation of the configurational distribution.

Ionorbital coupling in CarParrinello calculations of hydrogenbond vibrational dynamics: Case study with the NH_{3}–HCl dimer
View Description Hide DescriptionWe have performed CarParrinello molecular dynamics (CPMD) calculations of the hydrogenbonded NH_{3}–HCl dimer. Our main aim is to establish how ionicorbital coupling in CPMD affects the vibrational dynamics in hydrogenbonded systems by characterizing the dependence of the calculated vibrational frequencies upon the orbital mass in the adiabatic limit of CarParrinello calculations. We use the example of the NH_{3}–HCl dimer because of interest in its vibrational spectrum, in particular the magnitude of the frequency shift of the H–Cl stretch due to the anharmonic interactions when the hydrogen bond is formed. We find that an orbital mass of about 100 a.u. or smaller is required in order for the ionorbital coupling to be linear in orbital mass, and the results for which can be accurately extrapolated to the adiabatic limit of zero orbital mass. We argue that this is general for hydrogenbonded systems, suggesting that typical orbital mass values used in CPMD are too high to accurately describe vibrational dynamics in hydrogenbonded systems. Our results also show that the usual application of a scaling factor to the CPMD frequencies to correct for the effects of orbital mass is not valid. For the dynamics of the dimer, we find that the H–Cl stretch and the N–H–Cl bend are significantly coupled, suggesting that it is important to include the latter degree of freedom in quantum dynamical calculations. Results from our calculations with deuteriumsubstitution show that both these degrees of freedom have significant anharmonic interactions. Our calculated frequency for the H–Cl stretch using the Beckeexchange LeeYangParr correlation functional compares reasonably well with a previous secondorder MøllerPlesset calculation with anharmonic corrections, although it is low compared to the experimental value for the dimer trapped in a neonmatrix.

Relativistic JahnTeller effects in the photoelectron spectra of tetrahedral P_{4}, As_{4}, Sb_{4}, and Bi_{4}
View Description Hide DescriptionThe groupV tetrahedral cluster cations , , , and are known to exhibit exceptionally strong JahnTeller(JT)effects of electrostatic origin in their ^{2} Eground states and ^{2} T _{2}excited states. It has been predicted that there exist, in addition, JT couplings of relativistic origin (arising from the spinorbit (SO) operator) in ^{2} E and ^{2} T _{2} states of tetrahedral systems, which should become relevant for the heavier elements. In the present work, the JT and SO couplings in the groupV tetramer cations have been analyzed with ab initio relativistic electronic structure calculations. The vibronic line spectra and the band shapes of the photoelectron spectra were simulated with timedependent quantum wavepacket methods. The results provide insight into the interplay of electrostatic and relativistic JT couplings and SO splittings in the complex photoelectron spectra of these systems.

Single mode phonon energy transmission in functionalized carbon nanotubes
View Description Hide DescriptionAlthough the carbon nanotube(CNT) features superior thermal properties in its pristine form, the chemical functionalization often required for many applications of CNT inevitably degrades the structural integrity and affects the transport of energy carriers. In this article, the effect of the side wall functionalization on the phonon energy transmission along the symmetry axis of CNT is studied using the phononwave packet method. Three different functional groups are studied: methyl (–CH_{3}), vinyl (–C_{2}H_{3}), and carboxyl (–COOH). We find that, near Γ point of the Brillouin zone, acoustic phonons show ideal transmission, while the transmission of the optical phonons is strongly suppressed. A positive correlation between the energy transmission coefficient and the phonon group velocity is observed for both acoustic and optical phonon modes. On comparing the transmission due to functional groups with equivalent point mass defects on CNT, we find that the chemistry of the functional group, rather than its molecular mass, has a dominant role in determining phonon scattering, hence the transmission, at the defect sites.

Electron correlation via frozen Gaussian dynamics
View Description Hide DescriptionWe investigate the accuracy and efficiency of the semiclassical frozen Gaussian method in describing electron dynamics in real time. Modelsystems of two softCoulombinteracting electrons are used to study correlated dynamics under nonperturbative electric fields, as well as the excitation spectrum. The results show that a recently proposed method that combines exactexchange with semiclassical correlation to propagate the onebody densitymatrix holds promise for electron dynamics in many situations that either wavefunction or densityfunctional methods have difficulty describing. The results also however point out challenges in such a method that need to be addressed before it can become widely applicable.

A generalorder local coupledcluster method based on the clusterinmolecule approach
View Description Hide DescriptionA generalorder local coupledcluster (CC) method is presented which has the potential to provide accurate correlation energies for extended systems. Our method combines the clusterinmolecule approach of Li and coworkers [J. Chem. Phys.131, 114109 (2009)]10.1063/1.3218842 with the frozen natural orbital (NO) techniques widely used for the cost reduction of correlation methods. The occupied molecular orbitals (MOs) are localized, and for each occupied MO a local subspace of occupied and virtual orbitals is constructed using approximate Møller–Plesset NOs. The CC equations are solved and the correlation energies are calculated in the local subspace for each occupied MO, while the total correlation energy is evaluated as the sum of the individual contributions. The size of the local subspaces and the accuracy of the results can be controlled by varying only one parameter, the threshold for the occupation number of NOs which are included in the subspaces. Though our local CC method in its present form scales as the fifth power of the system size, our benchmark calculations show that it is still competitive for the CC singles and doubles (CCSD) and the CCSD with perturbative triples [CCSD(T)] approaches. For higher order CC methods, the reduction in computation time is more pronounced, and the new method enables calculations for considerably bigger molecules than before with a reasonable loss in accuracy. We also demonstrate that the independent calculation of the correlation contributions allows for a higher order description of the chemically more important segments of the molecule and a lower level treatment of the rest delivering further significant savings in computer time.

Amplitude equations for breathing spiral waves in a forced reactiondiffusion system
View Description Hide DescriptionBased on a multiple scale analysis of a forced reactiondiffusion system leading to amplitude equations, we explain the existence of spiral wave and its photoinduced spatiotemporal behavior in chlorine dioxideiodinemalonic acid system. When the photoillumination intensity is modulated, breathing of spiral is observed in which the period of breathing is identical to the period of forcing. We have also derived the condition for breakup and suppression of spiral wave by periodic illumination. The numerical simulations agree well with our analytical treatment.

Nonlocal continuum electrostatic theory predicts surprisingly small energetic penalties for charge burial in proteins
View Description Hide DescriptionWe study the energetics of burying charges, ion pairs, and ionizable groups in a simple proteinmodel using nonlocal continuum electrostatics. Our primary finding is that the nonlocal response leads to markedly reduced solvent screening, comparable to the use of applicationspecific proteindielectric constants. Employing the same parameters as used in other nonlocal studies, we find that for a sphere of radius 13.4 Å containing a single +1e charge, the nonlocal solvation free energy varies less than 18 kcal/mol as the charge moves from the surface to the center, whereas the difference in the local Poisson model is ∼35 kcal/mol. Because an ion pair (salt bridge) generates a comparatively more rapidly varying Coulomb potential, energetics for salt bridges are even more significantly reduced in the nonlocal model. By varying the central parameter in nonlocal theory, which is an effective length scale associated with correlations between solvent molecules, nonlocalmodel energetics can be varied from the standard local results to essentially zero; however, the existence of the reduction in chargeburial penalties is quite robust to variations in the proteindielectric constant and the correlation length. Finally, as a simple exploratory test of the implications of nonlocal response, we calculate glutamate pK _{a} shifts and find that using standard protein parameters (ε_{protein} = 2–4), nonlocal results match localmodel predictions with much higher dielectric constants. Nonlocality may, therefore, be one factor in resolving discrepancies between measured proteindielectric constants and the model parameters often used to match titration experiments. Nonlocal models may hold significant promise to deepen our understanding of macromolecular electrostatics without substantially increasing computational complexity.

Optimal diabatic bases via thermodynamic bounds
View Description Hide DescriptionDescribing kinetic processes within a perturbation theory approach such as Fermi's golden rule requires an understanding of the initial and final states of the system. A number of different methods have been proposed for obtaining these diabaticlike states, but a robust criterion for evaluating their accuracy has not been established. Here, we approach the problem of determining the most appropriate set of diabatic states for use in incoherent rate expressions. We develop a method that rotates an initial set of diabats into an optimized set beginning with a zerothorder diabatic Hamiltonian and choosing the rotation that minimizes the effect of nondiabatic terms on the thermodynamic free energy. The GibbsBogoliubov (GB) bound on the Helmholtz free energy is thus used as the diabatic criterion. We first derive the GB free energy for a two site system and then find an expression general for any electronic system Hamiltonian. Efficient numerical methods are used to perform the minimization subject to orthogonality constraints, and we examine the resulting diabats for system Hamiltonians in various parameter regimes. The transition from localized to delocalized states is clearly seen in these calculations, and some interesting features are discussed.

Mechanisms of kinetic trapping in selfassembly and phase transformation
View Description Hide DescriptionIn selfassembly processes, kinetic trapping effects often hinder the formation of thermodynamically stable ordered states. In a model of viral capsid assembly and in the phase transformation of a lattice gas, we show how simulations in a selfassembling steady state can be used to identify two distinct mechanisms of kinetic trapping. We argue that one of these mechanisms can be adequately captured by kinetic rate equations, while the other involves a breakdown of theories that rely on cluster size as a reaction coordinate. We discuss how these observations might be useful in designing and optimising selfassembly reactions.

Electric field effects on nuclear magnetic shielding of the 1:1 and 2:1 (homo and heterochiral) complexes of XOOX′ (X, X′ = H, CH_{3}) with lithium cation and their chiral discrimination
View Description Hide DescriptionThe set of 1:1 and 2:1 complexes of XOOX′ (X, X′ = H, CH_{3}) with lithium cation has been studied to determine if they are suitable candidates for chiral discrimination in an isotropic medium via nuclear magnetic resonance spectroscopy. Conventional nuclear magnetic resonance is unable to distinguish between enantiomers in the absence of a chiral solvent. The criterion for experimental detection is valuated by the isotropic part of nuclear shielding polarisabilitytensors, related to a pseudoscalar of opposite sign for two enantiomers. The study includes calculations at coupled HartreeFock and density functional theory schemes for ^{17}O nucleus in each compound. Additional calculations for ^{1}H are also included for some compounds. A huge static homogeneous electric field, perpendicular to the magnetic field of the spectromer, as big as ≈1.7 × 10^{8} V m^{−1} should be applied to observe a shift of ≈1 ppm for ^{17}O magnetic shielding in the proposed set of complexes.

Using swarm intelligence for finding transition states and reaction paths
View Description Hide DescriptionWe describe an algorithm that explores potential energy surfaces (PES) and finds approximate reaction paths and transition states. A few (≈6) evolving atomic configurations (“climbers”) start near a local minimum M1 of the PES. The climbers seek a shallow ascent, low energy, path toward a saddle point S12, cross over to another valley of the PES, and climb down to a new minimum M2 that was not known beforehand. Climbers use both energy and energy derivatives to make individual decisions, and they use relative fitness to make teambased decisions. In sufficiently long runs, they keep exploring and may go through a sequence M1–S12–M2–S23–M3 … of minima and saddle points without revisiting any of the critical points. We report results on eight small test systems that highlight advantages and disadvantages of the method. We also investigated the PES of Li_{8}, , Ag_{7}, and Ag_{2}NH_{3} to illustrate potential applications of this new method.
 Gas Phase Dynamics and Structure: Spectroscopy, Molecular Interactions, Scattering, and Photochemistry

Valence isoelectronic substitution in the B_{8} ^{−} and B_{9} ^{−} molecular wheels by an Al dopant atom: Umbrellalike structures of AlB_{7} ^{−} and AlB_{8} ^{−}
View Description Hide DescriptionThe structures and the electronic properties of two aluminumdoped boron clusters, AlB_{7} ^{−} and AlB_{8} ^{−}, were investigated using photoelectron spectroscopy and ab initio calculations. The photoelectron spectra of AlB_{7} ^{−} and AlB_{8} ^{−} are both broad, suggesting significant geometry changes between the ground states of the anions and the neutrals. Unbiased global minimum searches were carried out and the calculated vertical electron detachment energies were used to compare with the experimental data. We found that the Al atom does not simply replace a B atom in the parent B_{8} ^{−} and B_{9} ^{−} planar clusters in AlB_{7} ^{−} and AlB_{8} ^{−}. Instead, the global minima of the two dopedclusters are of umbrella shapes, featuring an Al atom interacting ionically with a hexagonal and heptagonal pyramidal B_{7} (C _{6v }) and B_{8} (C _{7v }) fragment, respectively. These unique umbrellatype structures are understood on the basis of the special stability of the quasiplanar B_{7} ^{3−} and planar B_{8} ^{2−} molecular wheels derived from double aromaticity.

Characterization of the potential minimum of the F ^{ ′}0_{ u } ^{ + }(^{1} D _{2}) ionpair state of Cl_{2} using (1 + 2^{′}) opticaloptical double resonance excitation and massresolved ion detection
View Description Hide DescriptionVibrational levels of the F ^{′}0_{ u } ^{ + }(^{1} D _{2}), (^{3} P _{0}), and D0_{ u } ^{ + }(^{3} P _{2}) ionpair states of ^{35}Cl_{2} and ^{35}Cl^{37}Cl in the range 62 500–67 600 cm^{−1} have been observed using (1 + 2^{′}) opticaloptical double resonance excitation with massresolved ion detection. The strong F ^{′}0_{ u } ^{ + }(^{1} D _{2})/(^{3} P _{0}) coupling has been modelled by a coupled twostate calculation. An optimized fit of the experimental data used an F ^{′}0_{ u } ^{ + }(^{1} D _{2}) state potential with a T _{ e } of 65 177 cm^{−1} and an R _{ e } of ≈2.636 Å with a coupling constant of ≈430 cm^{−1}. The calculation assigns the first observed members of the F ^{′}0_{ u } ^{ + }(^{1} D _{2}) state progression of ^{35}Cl_{2} and ^{35}Cl^{37}Cl at 64 998 and 65 094 cm^{−1}, respectively, as transitions to v = 0.

The visible spectrum of zirconium dioxide, ZrO_{2}
View Description Hide DescriptionThe electronic spectrum of a cold molecular beam of zirconium dioxide, ZrO_{2}, has been investigated using laser induced fluorescence(LIF) in the region from 17 000 cm^{−1} to 18 800 cm^{−1} and by massresolved resonance enhanced multiphoton ionization (REMPI) spectroscopy from 17 000 cm^{−1}–21 000 cm^{−1}. The LIF and REMPI spectra are assigned to progressions in the (ν_{1}, ν_{2}, ν_{3}) ← (0, 0, 0) transitions. Dispersed fluorescence from 13 bands was recorded and analyzed to produce harmonic vibrational parameters for the state of ω _{1} = 898(1) cm^{−1}, ω _{2} = 287(2) cm^{−1}, and ω _{3} = 808(3) cm^{−1}. The observed transition frequencies of 45 bands in the LIF and REMPI spectra produce origin and harmonic vibrational parameters for the state of T _{e} = 16 307(8) cm^{−1}, ω _{1} = 819(3) cm^{−1}, ω _{2} = 149(3) cm^{−1}, and ω _{3} = 518(4) cm^{−1}. The spectra were modeled using a normal coordinate analysis and FranckCondon factor predictions. The structures, harmonic vibrational frequencies, and the potential energies as a function of bending angle for the and states are predicted using timedependent density functional theory, complete active space selfconsistent field, and related firstprinciple calculations. A comparison with isovalent TiO_{2} is made.