Volume 140, Issue 17, 07 May 2014

The Frenkel exciton Hamiltonian is at the heart of many simulations of excitation energy transfer in molecular aggregates. It separates the aggregate into Coulombcoupled monomers. Here it is shown that the respective parameters, i.e., monomeric excitation energies and Coulomb couplings between transition densities can be efficiently calculated using timedependent tightbindingbased density functional theory (TDDFTB). Specifically, Coulomb couplings are expressed in terms of selfconsistently determined Mulliken transition charges. The approach is applied to two dimer systems. First, formaldehyde oxime for which a detailed comparison with standard DFT using the B3LYP and the PBE functionals as well as with SCSCC2 is provided. Second, the Coulomb coupling is explored in dependence on the intermolecular coordinates for a perylene bisimide dimer. This provides structural evidence for the previously observed biphasic aggregation behavior of this dye.
 COMMUNICATIONS


Communication: Quantum molecular dynamics simulation of liquid parahydrogen by nuclear and electron wave packet approach
View Description Hide DescriptionLiquid parahydrogen (pH2) is a typical quantum liquid which exhibits strong nuclear quantum effects (NQEs) and thus anomalous static and dynamic properties. We propose a realtime simulation method of wave packet (WP) molecular dynamics (MD) based on nonempirical intra and intermolecular interactions of nonspherical hydrogen molecules, and apply it to condensedphase pH2. The NQEs, such as WP delocalization and zeropoint energy, are taken into account without perturbative expansion of prepared model potential functions but with explicit interactions between nuclear and electron WPs. The developed MD simulation for 100 ps with 1200 hydrogen molecules is realized at feasible computational cost, by which basic experimental properties of pH2 liquid such as radial distribution functions, selfdiffusion coefficients, and shear viscosities are all well reproduced.

Communication: Helium nanodroplet isolation and rovibrational spectroscopy of hydroxymethylene
View Description Hide DescriptionHydroxymethylene (HCOH) and its d1isotopologue (HCOD) are isolated in low temperature helium nanodroplets following pyrolysis of glyoxylic acid. Transitions identified in the infrared spectrum are assigned exclusively to the transconformation based on previously reported anharmonic frequency computations [P. R. Schreiner, H. P. Reisenauer, F. C. Pickard, A. C. Simmonett, W. D. Allen, E. Mátyus, and A. G. Császár, Nature453, 906 (2008); L. Koziol, Y. M. Wang, B. J. Braams, J. M. Bowman, and A. I. Krylov, J. Chem. Phys.128, 204310 (2008)]. For the OH(D) and CH stretches, a and btype transitions are observed, and when taken in conjunction with CCSD(T)/ccpVTZ computations, lower limits to the vibrational band origins are determined. The relative intensities of the a and btype transitions provide the orientation of the transition dipole moment in the inertial frame. The He nanodroplet data are in excellent agreement with anharmonic frequency computations reported here and elsewhere, confirming an appreciable Armatrix shift of the OH and OD stretches and strong anharmonic resonance interactions in the highfrequency stretch regions of the midinfrared.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

A new efficient method for calculation of Frenkel exciton parameters in molecular aggregates
View Description Hide DescriptionThe Frenkel exciton Hamiltonian is at the heart of many simulations of excitation energy transfer in molecular aggregates. It separates the aggregate into Coulombcoupled monomers. Here it is shown that the respective parameters, i.e., monomeric excitation energies and Coulomb couplings between transition densities can be efficiently calculated using timedependent tightbindingbased density functional theory (TDDFTB). Specifically, Coulomb couplings are expressed in terms of selfconsistently determined Mulliken transition charges. The approach is applied to two dimer systems. First, formaldehyde oxime for which a detailed comparison with standard DFT using the B3LYP and the PBE functionals as well as with SCSCC2 is provided. Second, the Coulomb coupling is explored in dependence on the intermolecular coordinates for a perylene bisimide dimer. This provides structural evidence for the previously observed biphasic aggregation behavior of this dye.

Ionic size effects to molecular solvation energy and to ion current across a channel resulted from the nonuniform sizemodified PNP equations
View Description Hide DescriptionIonic finite size can impose considerable effects to both the equilibrium and nonequilibrium properties of a solvated molecular system, such as the solvation energy, ionic concentration, and transport in a channel. As discussed in our former work [B. Lu and Y. C. Zhou, Biophys. J.100, 2475 (2011)], a class of sizemodified PoissonBoltzmann (PB)/PoissonNernstPlanck (PNP) models can be uniformly studied through the general nonuniform sizemodified PNP (SMPNP) equations deduced from the extended free energy functional of Borukhov et al. [I. Borukhov, D. Andelman, and H. Orland, Phys. Rev. Lett.79, 435 (1997)] This work focuses on the nonuniform size effects to molecular solvation energy and to ion current across a channel for real biomolecular systems. The main contributions are: (1) we prove that for solvation energy calculation with nonuniform size effects (through equilibrium SMPNP simulation), there exists a simplified approximation formulation which is the same as the widely used one in PB community. This approximate form avoids integration over the whole domain and makes energy calculations convenient. (2) Numerical calculations show that ionic size effects tend to negate the solvation effects, which indicates that a higher molecular solvation energy (lower absolute value) is to be predicted when ionic size effects are considered. For both calculations on a protein and a DNA fragment systems in a 0.5M 1:1 ionic solution, a difference about 10 kcal/mol in solvation energies is found between the PB and the SMPNP predictions. Moreover, it is observed that the solvation energy decreases as ionic strength increases, which behavior is similar as those predicted by the traditional PB equation (without size effect) and by the uniform sizemodified PoissonBoltzmann equation. (3) Nonequilibrium SMPNP simulations of ion permeation through a gramicidin A channel show that the ionic size effects lead to reduced ion current inside the channel compared with the results without considering size effects. As a component of the current, the drift term is the main contribution to the total current. The ionic size effects to the total current almost come through the drift term, and have little influence on the diffusion terms in SMPNP.

Analytical gradients of complete active space selfconsistent field energies using Cholesky decomposition: Geometry optimization and spinstate energetics of a ruthenium nitrosyl complex
View Description Hide DescriptionWe present a formulation of analytical energy gradients at the complete active space selfconsistent field (CASSCF) level of theory employing density fitting (DF) techniques to enable efficient geometry optimizations of large systems. As an example, the ground and lowest triplet state geometries of a ruthenium nitrosyl complex are computed at the DFCASSCF level of theory and compared with structures obtained from density functional theory (DFT) using the B3LYP, BP86, and M06L functionals. The average deviation of all bond lengths compared to the crystal structure is 0.042 Å at the DFCASSCF level of theory, which is slightly larger but still comparable with the deviations obtained by the tested DFT functionals, e.g., 0.032 Å with M06L. Specifically, the rootmeansquare deviation between the DFCASSCF and best DFT coordinates, delivered by BP86, is only 0.08 Å for S0 and 0.11 Å for T1, indicating that the geometries are very similar. While keeping the mean energy gradient errors below 0.25%, the DF technique results in a 13fold speedup compared to the conventional CASSCF geometry optimization algorithm. Additionally, we assess the singlettriplet energy vertical and adiabatic differences with multiconfigurational secondorder perturbation theory (CASPT2) using the DFCASSCF and DFT optimized geometries. It is found that the vertical CASPT2 energies are relatively similar regardless of the geometry employed whereas the adiabatic singlettriplet gaps are more sensitive to the chosen triplet geometry.

Solidstate dimer method for calculating solidsolid phase transitions
View Description Hide DescriptionThe dimer method is a minimum mode following algorithm for finding saddle points on a potential energy surface of atomic systems. Here, the dimer method is extended to include the cell degrees of freedom for periodic solidstate systems. Using this method, reaction pathways of solidsolid phase transitions can be determined without having to specify the final state structure or reaction mechanism. Example calculations include concerted phase transitions between CdSe polymorphs and a nucleation and growth mechanism for the A15 to BCC transition in Mo.

Mixed quantum classical calculation of proton transfer reaction rates: From deep tunneling to over the barrier regimes
View Description Hide DescriptionWe present mixed quantum classical calculations of the proton transfer (PT) reaction rates represented by a double well system coupled to a dissipative bath. The rate constants are calculated within the so called nontraditional view of the PT reaction, where the proton motion is quantized and the solvent polarization is used as the reaction coordinate. Quantization of the proton degree of freedom results in a problem of nonadiabatic dynamics. By employing the reactive flux formulation of the rate constant, the initial sampling starts from the transition state defined using the collective reaction coordinate. Dynamics of the collective reaction coordinate is treated classically as over damped diffusive motion, for which the equation of motion can be derived using the path integral, or the mixed quantum classical Liouville equation methods. The calculated mixed quantum classical rate constants agree well with the results from the numerically exact hierarchical equation of motion approach for a broad range of model parameters. Moreover, we are able to obtain contributions from each vibrational state to the total reaction rate, which helps to understand the reaction mechanism from the deep tunneling to over the barrier regimes. The numerical results are also compared with those from existing approximate theories based on calculations of the nonadiabatic transmission coefficients. It is found that the twosurface LandauZener formula works well in calculating the transmission coefficients in the deep tunneling regime, where the crossing point between the two lowest vibrational states dominates the total reaction rate. When multiple vibrational levels are involved, including additional crossing points on the free energy surfaces is important to obtain the correct reaction rate using the LandauZener formula.

Predicting the influence of longrange molecular interactions on macroscopicscale diffusion by homogenization of the Smoluchowski equation
View Description Hide DescriptionThe macroscopic diffusion constant for a charged diffuser is in part dependent on (1) the volume excluded by solute “obstacles” and (2) longrange interactions between those obstacles and the diffuser. Increasing excluded volume reduces transport of the diffuser, while longrange interactions can either increase or decrease diffusivity, depending on the nature of the potential. We previously demonstrated [P. M. KekenesHuskey et al. , Biophys. J.105, 2130 (2013)] using homogenization theory that the configuration of molecularscale obstacles can both hinder diffusion and induce diffusional anisotropy for small ions. As the density of molecular obstacles increases, van der Waals (vdW) and electrostatic interactions between obstacle and a diffuser become significant and can strongly influence the latter's diffusivity, which was neglected in our original model. Here, we extend this methodology to include a fixed (timeindependent) potential of mean force, through homogenization of the Smoluchowski equation. We consider the diffusion of ions in crowded, hydrophilic environments at physiological ionic strengths and find that electrostatic and vdW interactions can enhance or depress effective diffusion rates for attractive or repulsive forces, respectively. Additionally, we show that the observed diffusion rate may be reduced independent of nonspecific electrostatic and vdW interactions by treating obstacles that exhibit specific binding interactions as “buffers” that absorb free diffusers. Finally, we demonstrate that effective diffusion rates are sensitive to distribution of surface charge on a globular protein, Troponin C, suggesting that the use of molecular structures with atomisticscale resolution can account for electrostatic influences on substrate transport. This approach offers new insight into the influence of molecularscale, longrange interactions on transport of charged species, particularly for diffusioninfluenced signaling events occurring in crowded cellular environments.

Model reduction for slow–fast stochastic systems with metastable behaviour
View Description Hide DescriptionThe quasisteadystate approximation (or stochastic averaging principle) is a useful tool in the study of multiscale stochastic systems, giving a practical method by which to reduce the number of degrees of freedom in a model. The method is extended here to slow–fast systems in which the fast variables exhibit metastable behaviour. The key parameter that determines the form of the reduced model is the ratio of the timescale for the switching of the fast variables between metastable states to the timescale for the evolution of the slow variables. The method is illustrated with two examples: one from biochemistry (a fastspeciesmediated chemical switch coupled to a slower varying species), and one from ecology (a predator–prey system). Numerical simulations of each model reduction are compared with those of the full system.

Compressible generalized hybrid Monte Carlo
View Description Hide DescriptionOne of the most demanding calculations is to generate random samples from a specified probability distribution (usually with an unknown normalizing prefactor) in a highdimensional configuration space. One often has to resort to using a Markov chain Monte Carlo method, which converges only in the limit to the prescribed distribution. Such methods typically inch through configuration space step by step, with acceptance of a step based on a Metropolis(Hastings) criterion. An acceptance rate of 100% is possible in principle by embedding configuration space in a higher dimensional phase space and using ordinary differential equations. In practice, numerical integrators must be used, lowering the acceptance rate. This is the essence of hybrid Monte Carlo methods. Presented is a general framework for constructing such methods under relaxed conditions: the only geometric property needed is (weakened) reversibility; volume preservation is not needed. The possibilities are illustrated by deriving a couple of explicit hybrid Monte Carlo methods, one based on barrierlowering variablemetric dynamics and another based on isokinetic dynamics.

Graphics processing units accelerated semiclassical initial value representation molecular dynamics
View Description Hide DescriptionThis paper presents a Graphics Processing Units (GPUs) implementation of the Semiclassical Initial Value Representation (SCIVR) propagator for vibrational molecular spectroscopy calculations. The timeaveraging formulation of the SCIVR for power spectrum calculations is employed. Details about the GPU implementation of the semiclassical code are provided. Four molecules with an increasing number of atoms are considered and the GPUcalculated vibrational frequencies perfectly match the benchmark values. The computational time scaling of two GPUs (NVIDIA Tesla C2075 and Kepler K20), respectively, versus two CPUs (Intel Core i5 and Intel Xeon E52687W) and the critical issues related to the GPU implementation are discussed. The resulting reduction in computational time and power consumption is significant and semiclassical GPU calculations are shown to be environment friendly.

A Gibbsensemble based technique for Monte Carlo simulation of electric double layer capacitors (EDLC) at constant voltage
View Description Hide DescriptionCurrent methods for molecular simulations of Electric Double Layer Capacitors (EDLC) have both the electrodes and the electrolyte region in a single simulation box. This necessitates simulation of the electrodeelectrolyte region interface. Typical capacitors have macroscopic dimensions where the fraction of the molecules at the electrodeelectrolyte region interface is very low. Hence, large systems sizes are needed to minimize the electrodeelectrolyte region interfacial effects. To overcome these problems, a new technique based on the Gibbs Ensemble is proposed for simulation of an EDLC. In the proposed technique, each electrode is simulated in a separate simulation box. Application of periodic boundary conditions eliminates the interfacial effects. This in addition to the use of constant voltage ensemble allows for a more convenient comparison of simulation results with experimental measurements on typical EDLCs.

Calculating vibrational spectra with sum of product basis functions without storing fulldimensional vectors or matrices
View Description Hide DescriptionWe propose an iterative method for computing vibrational spectra that significantly reduces the memory cost of calculations. It uses a direct product primitive basis, but does not require storing vectors with as many components as there are product basis functions. Wavefunctions are represented in a basis each of whose functions is a sum of products (SOP) and the factorizable structure of the Hamiltonian is exploited. If the factors of the SOP basis functions are properly chosen, wavefunctions are linear combinations of a small number of SOP basis functions. The SOP basis functions are generated using a shifted block power method. The factors are refined with a rank reduction algorithm to cap the number of terms in a SOP basis function. The ideas are tested on a 20D model Hamiltonian and a realistic CH3CN (12 dimensional) potential. For the 20D problem, to use a standard direct product iterative approach one would need to store vectors with about 10^{20} components and would hence require about 8 × 10^{11} GB. With the approach of this paper only 1 GB of memory is necessary. Results for CH3CN agree well with those of a previous calculation on the same potential.

An algorithm for nonrelativistic quantummechanical finitenuclearmass variational calculations of nitrogen atom in L = 0, M = 0 states using allelectrons explicitly correlated Gaussian basis functions
View Description Hide DescriptionAn algorithm for quantummechanical nonrelativistic variational calculations of L = 0 and M = 0 states of atoms with an arbitrary number of s electrons and with three p electrons have been implemented and tested in the calculations of the ground ^{4}S state of the nitrogen atom. The spatial part of the wave function is expanded in terms of allelectrons explicitly correlated Gaussian functions with the appropriate preexponential Cartesian angular factors for states with the L = 0 and M = 0 symmetry. The algorithm includes formulas for calculating the Hamiltonian and overlap matrix elements, as well as formulas for calculating the analytic energy gradient determined with respect to the Gaussian exponential parameters. The gradient is used in the variational optimization of these parameters. The Hamiltonian used in the approach is obtained by rigorously separating the centerofmass motion from the laboratoryframe allparticle Hamiltonian, and thus it explicitly depends on the finite mass of the nucleus. With that, the mass effect on the total groundstate energy is determined.

Strong field ionization rates simulated with timedependent configuration interaction and an absorbing potential
View Description Hide DescriptionIonization rates of molecules have been modeled with timedependent configuration interaction simulations using atom centered basis sets and a complex absorbing potential. The simulations agree with accurate gridbased calculations for the ionization of hydrogen atom as a function of field strength and for charge resonance enhanced ionization of as the bond is elongated. Unlike gridbased methods, the present approach can be applied to simulate electron dynamics and ionization in multielectron polyatomic molecules. Calculations on HCl^{+} and HCO^{+} demonstrate that these systems also show charge resonance enhanced ionization as the bonds are stretched.

Equationofmotion coupled cluster perturbation theory revisited
View Description Hide DescriptionThe equationofmotion coupled cluster (EOMCC) framework has been used for deriving a novel series of perturbative corrections to the coupled cluster singles and doubles energy that formally converges towards the full configuration interaction energy limit. The series is based on a MøllerPlesset partitioning of the Hamiltonian and thus size extensive at any order in the perturbation, thereby remedying the major deficiency inherent to previous perturbation series based on the EOMCC ansatz.

Quantum mechanical/molecular mechanical/continuum style solvation model: Second order MøllerPlesset perturbation theory
View Description Hide DescriptionA combined quantum mechanical/molecular mechanical/continuum (QM/MM/C) style second order MøllerPlesset perturbation theory (MP2) method that incorporates induced dipole polarizable force field and induced surface charge continuum solvation model is established. The Zvector method is modified to include induced dipoles and induced surface charges to determine the MP2 response density matrix, which can be used to evaluate MP2 properties. In particular, analytic nuclear gradient is derived and implemented for this method. Using the Assisted Model Building with Energy Refinement induced dipole polarizable protein force field, the QM/MM/C style MP2 method is used to study the hydrogen bonding distances and strengths of the photoactive yellow protein chromopore in the wild type and the Glu46Gln mutant.
 Advanced Experimental Techniques

Atomically precise (catalytic) particles synthesized by a novel cluster deposition instrument
View Description Hide DescriptionWe report a new high vacuum instrument which is dedicated to the preparation of welldefined clusters supported on model and technologically relevant supports for catalytic and materials investigations. The instrument is based on deposition of size selected metallic cluster ions that are produced by a high flux magnetron cluster source. The throughput of the apparatus is maximized by collecting and focusing ions utilizing a conical octupole ion guide and a linear ion guide. The size selection is achieved by a quadrupole mass filter. The new design of the sample holder provides for the preparation of multiple samples on supports of various sizes and shapes in one session. After cluster deposition onto the support of interest, samples will be taken out of the chamber for a variety of testing and characterization.
 Atoms, Molecules, and Clusters

Homogenous nucleation rates of npropanol measured in the Laminar Flow Diffusion Chamber at different total pressures
View Description Hide DescriptionNucleation rates of npropanol were investigated in the Laminar Flow Diffusion Chamber. Nucleation temperatures between 270 and 300 K and rates between 10^{0} and 10^{6} cm^{−3} s^{−1} were achieved. Since earlier measurements of nbutanol and n‑pentanol suggest a dependence of nucleation rates on carrier gas pressure, similar conditions were adjusted for these measurements. The obtained data fit well to results available from literature. A small positive pressure effect was found which strengthen the assumption that this effect is attributed to the carbon chain length of the nalcohol [D. Brus, A. P. Hyvärinen, J. Wedekind, Y. Viisanen, M. Kulmala, V. Ždímal, J. Smolík, and H. Lihavainen, J. Chem. Phys.128, 134312 (2008)] and might be less intensive for substances in the homologous series with higher equilibrium vapor pressure. A comparison with the theoretical approach by Wedekind et al. [Phys. Rev. Lett.101, 12 (2008)] shows that the effect goes in the same direction but that the intensity is much stronger in experiments than in theory.

Pairing and unpairing electron densities in organic systems: Twoelectron three center through space and through bonds interactions
View Description Hide DescriptionTwoelectron threecenter bonding interactions in organic ions like methonium ( ), ethonium ( ), and protonated alkanes isomers (butonium cations) are described and characterized within the theoretical framework of the topological analysis of the electron density decomposition into its effectively paired and unpaired contributions. These interactions manifest in some of this type of systems as a concentration of unpaired electron cloud around the bond paths, in contrast to the well known paradigmatic boron hydrids in which it is not only concentrated close to the atomic nucleus and the bond paths but out of them and over the region defined by the involved atoms as a whole. This result permits to propose an attempt of classification for these interactions based in such manifestations. In the first type, it is called as interactions through bonds and in the second type as interactions through space type.