Volume 143, Issue 8, 28 August 2015
Index of content:

Molecular dynamics simulations and integral equation calculations of a simple equimolar mixture of diatomic molecules and monomers interacting via attractive and repulsive shortrange potentials show the existence of pattern formation (microheterogeneity), mostly due to depletion forces away from the demixing region. Effective sitesite potentials extracted from the pair correlation functions using an inverse Monte Carlo approach and an integral equation inversion procedure exhibit the features characteristic of a shortrange attractive and a longrange repulsive potential. When charges are incorporated into the model, this becomes a coarse grained representation of a room temperature ionic liquid, and as expected, intermediate range order becomes more pronounced and stable.
 COMMUNICATIONS


Communication: Conductivity enhancement in plasticcrystalline solidstate electrolytes
View Description Hide DescriptionFinding new ionic conductors that enable significant advancements in the development of energystorage devices is a challenging goal of current material science. Aside of material classes as ionic liquids or amorphous ion conductors, the socalled plastic crystals (PCs) have been shown to be good candidates combining high conductivity and favorable mechanical properties. PCs are formed by molecules whose orientational degrees of freedom still fluctuate despite the material exhibits a welldefined crystalline lattice. In the present work, we show that the conductivity of Li^{+} ions in succinonitrile, the most prominent molecular PC electrolyte, can be enhanced by several decades when replacing part of the molecules in the crystalline lattice by larger ones. Dielectric spectroscopy reveals that this is accompanied by a stronger coupling of ionic and reorientational motions. These findings, which can be understood in terms of an optimized “revolving door” mechanism, open a new path towards the development of better solidstate electrolytes.

Communication: Test of quantum chemistry in vibrationally hot hydrogen molecules
View Description Hide DescriptionPrecision measurements are performed on highly excited vibrational quantum states of molecular hydrogen. The v = 12, J = 0 − 3 rovibrational levels of H2 ( ), lying only 2000 cm^{−1} below the first dissociation limit, were populated by photodissociation of H2S and their level energies were accurately determined by twophoton Dopplerfree spectroscopy. A comparison between the experimental results on v = 12 level energies with the best ab initio calculations shows a good agreement, where the present experimental accuracy of 3.5 × 10^{−3} cm^{−1} is more precise than theory, hence providing a gateway to further test theoretical advances in this benchmark quantum system.

Communication: Activation energy of tensioninduced pore formation in lipid membranes
View Description Hide DescriptionTension plays a vital role in pore formation in biomembranes, but the mechanism of pore formation remains unclear. We investigated the temperature dependence of the rate constant of constant tension (σ)–induced pore formation in giant unilamellar vesicles of lipid membranes using an experimental method we developed. By analyzing this result, we determined the activation energy (U a) of tensioninduced pore formation as a function of tension. A constant (U 0) that does not depend on tension was found to contribute significantly to U a. Analysis of the activation energy clearly indicated that the dependence of U a on σ in the classical theory is correct, but that the classical theory of pore formation is not entirely correct due to the presence of U 0. We can reasonably consider that U 0 is a nucleation free energy to form a hydrophilic prepore from a hydrophobic prepore or a region with lower lateral lipid density. After obtaining U 0, the evolution of a prepore follows a classical theory. Our data provide valuable information that help explain the mechanism of tensioninduced pore formation in biomembranes and lipid membranes.

 ARTICLES

 Theoretical Methods and Algorithms

Thermal weights for semiclassical vibrational response functions
View Description Hide DescriptionSemiclassical approximations to response functions can allow the calculation of linear and nonlinear spectroscopic observables from classical dynamics. Evaluating a canonical response function requires the related tasks of determining thermal weights for initial states and computing the dynamics of these states. A class of approximations for vibrational response functions employs classical trajectories at quantized values of action variables and represents the effects of the radiationmatter interaction by discontinuous transitions. Here, we evaluate choices for a thermal weight function which are consistent with this dynamical approximation. Weight functions associated with different semiclassical approximations are compared, and two forms are constructed which yield the correct linear response function for a harmonic potential at any temperature and are also correct for anharmonic potentials in the classical mechanical limit of high temperature. Approximations to the vibrational linear response function with quantized classical trajectories and proposed thermal weight functions are assessed for ensembles of onedimensional anharmonic oscillators. This approach is shown to perform well for an anharmonic potential that is not locally harmonic over a temperature range encompassing the quantum limit of a twolevel system and the limit of classical dynamics.

Intrachain exciton dynamics in conjugated polymer chains in solution
View Description Hide DescriptionWe investigate exciton dynamics on a polymer chain in solution induced by the Brownian rotational motion of the monomers. Poly(paraphenylene) is chosen as the model system and excitons are modeled via the Frenkel exciton Hamiltonian. The Brownian fluctuations of the torsional modes were modeled via the Langevin equation. The rotation of monomers in polymer chains in solution has a number of important consequences for the excited state properties. First, the dihedral angles assume a thermal equilibrium which causes offdiagonal disorder in the Frenkel Hamiltonian. This disorder Anderson localizes the Frenkel exciton centerofmass wavefunctions into superlocalized local exciton ground states (LEGSs) and higherenergy more delocalized quasiextended exciton states (QEESs). LEGSs correspond to chromophores on polymer chains. The second consequence of rotations—that are lowfrequency—is that their coupling to the exciton wavefunction causes local planarization and the formation of an excitonpolaron. This torsional relaxation causes additional selflocalization. Finally, and crucially, the torsional dynamics cause the Frenkel Hamiltonian to be timedependent, leading to exciton dynamics. We identify two distinct types of dynamics. At low temperatures, the torsional fluctuations act as a perturbation on the polaronic nature of the exciton state. Thus, the exciton dynamics at low temperatures is a smalldisplacement diffusive adiabatic motion of the excitonpolaron as a whole. The temperature dependence of the diffusion constant has a linear dependence, indicating an activationless process. As the temperature increases, however, the diffusion constant increases at a faster than linear rate, indicating a second nonadiabatic dynamics mechanism begins to dominate. Excitons are thermally activated into higher energy more delocalized exciton states (i.e., LEGSs and QEESs). These states are not selflocalized by local torsional planarization. During the exciton’s temporary occupation of a LEGS—and particularly a quasiband QEES—its motion is semiballistic with a large group velocity. After a short period of rapid transport, the exciton wavefunction collapses again into an excitonpolaron state. We present a simple model for the activated dynamics which is in agreement with the data.

A new approach to the problem of bulkmediated surface diffusion
View Description Hide DescriptionThis paper is devoted to bulkmediated surface diffusion of a particle which can diffuse both on a flat surface and in the bulk layer above the surface. It is assumed that the particle is on the surface initially (at t = 0) and at time t, while in between it may escape from the surface and come back any number of times. We propose a new approach to the problem, which reduces its solution to that of a twostate problem of the particle transitions between the surface and the bulk layer, focusing on the cumulative residence times spent by the particle in the two states. These times are random variables, the sum of which is equal to the total observation time t. The advantage of the proposed approach is that it allows for a simple exact analytical solution for the double Laplace transform of the conditional probability density of the cumulative residence time spent on the surface by the particle observed for time t. This solution is used to find the Laplace transform of the particle mean square displacement and to analyze the peculiarities of its time behavior over the entire range of time. We also establish a relation between the double Laplace transform of the conditional probability density and the FourierLaplace transform of the particle propagator over the surface. The proposed approach treats the cases of both finite and infinite bulk layer thicknesses (where bulkmediated surface diffusion is normal and anomalous at asymptotically long times, respectively) on equal footing.

Construction of exchange repulsion in terms of the wave functions at QM/MM boundary region
View Description Hide DescriptionWe developed a simple method to calculate exchange repulsion between a quantum mechanical (QM) solute and a molecular mechanical (MM) molecule in the QM/MM approach. In our method, the size parameter in the Buckingham type potential for the QM solute is directly determined in terms of the oneelectron wave functions of the solute. The point of the method lies in the introduction of the exchange core function (ECF) defined as a Slater function which mimics the behavior of the exterior electron density at the QM/MM boundary region. In the present paper, the ECF was constructed in terms of the BeckeRoussel (BR) exchange hole function. It was demonstrated that the ECF yielded by the BR procedure can faithfully reproduce the radial behavior of the electron density of a QM solute. The size parameter of the solute as well as the exchange repulsion are, then, obtained using the overlap model without any fitting procedure. To examine the efficiency of the method, it was applied to calculation of the exchange repulsions for minimal QM/MM systems, hydrogenbonded water dimer, and H3O^{+}–H2O. We found that our approach is able to reproduce the potential energy curves for these systems showing reasonable agreements with those given by accurate full quantum chemical calculations.

The geometry of generalized force matching and related information metrics in coarsegraining of molecular systems
View Description Hide DescriptionUsing the probabilistic language of conditional expectations, we reformulate the force matching method for coarsegraining of molecular systems as a projection onto spaces of coarse observables. A practical outcome of this probabilistic description is the link of the force matching method with thermodynamic integration. This connection provides a way to systematically construct a local mean force and to optimally approximate the potential of mean force through force matching. We introduce a generalized force matching condition for the local mean force in the sense that allows the approximation of the potential of mean force under both linear and nonlinear coarse graining mappings (e.g., reaction coordinates, endtoend length of chains). Furthermore, we study the equivalence of force matching with relative entropy minimization which we derive for general nonlinear coarse graining maps. We present in detail the generalized force matching condition through applications to specific examples in molecular systems.

A general ansatz for constructing quasidiabatic states in electronically excited aggregated systems
View Description Hide DescriptionWe present a general method for analyzing the character of singly excited states in terms of charge transfer (CT) and locally excited (LE) configurations. The analysis is formulated for configuration interaction singles (CIS) singly excited wave functions of aggregate systems. It also approximately works for the secondorder approximate coupled cluster singles and doubles and the secondorder algebraicdiagrammatic construction methods [CC2 and ADC(2)]. The analysis method not only generates a weight of each character for an excited state, but also allows to define the related quasidiabatic states and corresponding coupling matrix elements. In the character analysis approach, we divide the target system into domains and use a modified PipekMezey algorithm to localize the canonical MOs on each domain, respectively. The CIS wavefunction is then transformed into the localized basis, which allows us to partition the wavefunction into LE configurations within domains and CT configuration between pairs of different domains. Quasidiabatic states are then obtained by mixing excited states subject to the condition of maximizing the weight of one single LE or CT configuration (localization in configuration space). Different aims of such a procedure are discussed, either the construction of pure LE and CT states for analysis purposes (by including a large number of excited states) or the construction of effective models for dynamics calculations (by including a restricted number of excited states). Applications are given to LE/CT mixing in πstacked systems, chargerecombination matrix elements in a heterodimer, and excitonic couplings in multichromophoric systems.

The raspberry model for hydrodynamic interactions revisited. I. Periodic arrays of spheres and dumbbells
View Description Hide DescriptionThe socalled “raspberry” model refers to the hybrid latticeBoltzmann and Langevin molecular dynamics scheme for simulating the dynamics of suspensions of colloidal particles, originally developed by Lobaskin and Dünweg [New J. Phys. 6, 54 (2004)], wherein discrete surface points are used to achieve fluidparticle coupling. This technique has been used in many simulation studies on the behavior of colloids. However, there are fundamental questions with regards to the use of this model. In this paper, we examine the accuracy with which the raspberry method is able to reproduce Stokeslevel hydrodynamic interactions when compared to analytic expressions for solid spheres in simplecubic crystals. To this end, we consider the quality of numerical experiments that are traditionally used to establish these properties and we discuss their shortcomings. We show that there is a discrepancy between the translational and rotational mobility reproduced by the simple raspberry model and present a way to numerically remedy this problem by adding internal coupling points. Finally, we examine a nonconvex shape, namely, a colloidal dumbbell, and show that the filled raspberry model replicates the desired hydrodynamic behavior in bulk for this more complicated shape. Our investigation is continued in de Graaf et al. [J. Chem. Phys. 143, 084108 (2015)], wherein we consider the raspberry model in the confining geometry of two parallel plates.

The Raspberry model for hydrodynamic interactions revisited. II. The effect of confinement
View Description Hide DescriptionThe socalled “raspberry” model refers to the hybrid latticeBoltzmann (LB) and Langevin molecular dynamics schemes for simulating the dynamics of suspensions of colloidal particles, originally developed by Lobaskin and Dünweg [New J. Phys. 6, 54 (2004)], wherein discrete surface points are used to achieve fluidparticle coupling. In this paper, we present a follow up to our study of the effectiveness of the raspberry model in reproducing hydrodynamic interactions in the Stokes regime for spheres arranged in a simplecubic crystal [Fischer et al., J. Chem. Phys. 143, 084107 (2015)]. Here, we consider the accuracy with which the raspberry model is able to reproduce such interactions for particles confined between two parallel plates. To this end, we compare our LB simulation results to established theoretical expressions and finiteelement calculations. We show that there is a discrepancy between the translational and rotational mobilities when only surface coupling points are used, as also found in Part I of our joint publication. We demonstrate that adding internal coupling points to the raspberry can be used to correct said discrepancy in confining geometries as well. Finally, we show that the raspberry model accurately reproduces hydrodynamic interactions between a spherical colloid and planar walls up to roughly one LB lattice spacing.

Largescale epitaxial growth kinetics of graphene: A kinetic Monte Carlo study
View Description Hide DescriptionEpitaxial growth via chemical vapor deposition is considered to be the most promising way towards synthesizing large area graphene with high quality. However, it remains a big theoretical challenge to reveal growth kinetics with atomically energetic and largescale spatial information included. Here, we propose a minimal kinetic Monte Carlo model to address such an issue on an active catalyst surface with graphene/substrate lattice mismatch, which facilitates us to perform large scale simulations of the growth kinetics over two dimensional surface with growth fronts of complex shapes. A geometrydetermined largescale growth mechanism is revealed, where the ratedominating event is found to be C 1attachment for concave growthfront segments and C 5attachment for others. This growth mechanism leads to an interesting timeresolved growth behavior which is well consistent with that observed in a recent scanning tunneling microscopy experiment.

Theoretical description of spinselective reactions of radical pairs diffusing in spherical 2D and 3D microreactors
View Description Hide DescriptionIn this work, we treat spinselective recombination of a geminate radical pair (RP) in a spherical “microreactor,” i.e., of a RP confined in a micelle, vesicle, or liposome. We consider the microreactor model proposed earlier, in which one of the radicals is located at the center of the micelle and the other one undergoes threedimensional diffusion inside the micelle. In addition, we suggest a twodimensional model, in which one of the radicals is located at the “pole” of the sphere, while the other one diffuses on the spherical surface. For this model, we have obtained a general analytical expression for the RP recombination yield in terms of the free Green function of twodimensional diffusion motion. In turn, this Green function is expressed via the Legendre functions and thus takes account of diffusion over a restricted spherical surface and its curvature. The obtained expression allows one to calculate the RP recombination efficiency at an arbitrary magnetic field strength. We performed a comparison of the two models taking the same geometric parameters (i.e., the microreactor radius and the closest approach distance of the radicals), chemical reactivity, magnetic interactions in the RP and diffusion coefficient. Significant difference between the predictions of the two models is found, which is thus originating solely from the dimensionality effect: for different dimensionality of space, the statistics of diffusional contacts of radicals becomes different altering the reaction yield. We have calculated the magnetic field dependence of the RP reaction yield and chemically induced dynamic nuclear polarization of the reaction products at different sizes of the microreactor, exchange interaction, and spin relaxation rates. Interestingly, due to the intricate interplay of diffusional contacts of reactants and spin dynamics, the dependence of the reaction yield on the microreactor radius is nonmonotonous. Our results are of importance for (i) interpreting experimental data for magnetic field effects on RP recombination in confined space and (ii) for describing kinetics of chemical reactions, which occur predominantly on the surfaces of biomembranes, i.e., lipid peroxidation reactions.

Electronic spectra from TDDFT and machine learning in chemical space
View Description Hide DescriptionDue to its favorable computational efficiency, timedependent (TD) density functional theory (DFT) enables the prediction of electronic spectra in a highthroughput manner across chemical space. Its predictions, however, can be quite inaccurate. We resolve this issue with machine learning models trained on deviations of reference secondorder approximate coupledcluster (CC2) singles and doubles spectra from TDDFT counterparts, or even from DFT gap. We applied this approach to lowlying singletsinglet vertical electronic spectra of over 20 000 synthetically feasible small organic molecules with up to eight CONF atoms. The prediction errors decay monotonously as a function of training set size. For a training set of 10 000 molecules, CC2 excitation energies can be reproduced to within ±0.1 eV for the remaining molecules. Analysis of our spectral database via chromophore counting suggests that even higher accuracies can be achieved. Based on the evidence collected, we discuss open challenges associated with datadriven modeling of highlying spectra and transition intensities.

Energy error bars in direct configuration interaction iteration sequence
View Description Hide DescriptionA computational scheme for approximate lower bound to eigenvalues of linear operators is elaborated, based on Löwdin’s bracketing function. Implementation in direct full configuration interaction algorithm is presented, generating essentially just input–output increase. While strict lower bound property is lost due to approximations, test calculations result lower bounds of the same order of magnitude, as the usual upper bound, provided by the expectation value. Difference of upper and lower bounds gives an error bar, characterizing the wavefunction at the given iteration step.

Auxiliary matrix formalism for interaction representation transformations, optimal control, and spin relaxation theories
View Description Hide DescriptionAuxiliary matrix exponential method is used to derive simple and numerically efficient general expressions for the following, historically rather cumbersome, and hard to compute, theoretical methods: (1) average Hamiltonian theory following interaction representation transformations; (2) BlochRedfieldWangsness theory of nuclear and electron relaxation; (3) gradient ascent pulse engineering version of quantum optimal control theory. In the context of spin dynamics, the auxiliary matrix exponential method is more efficient than methods based on matrix factorizations and also exhibits more favourable complexity scaling with the dimension of the Hamiltonian matrix.

Uniform electron gases. III. Lowdensity gases on threedimensional spheres
View Description Hide DescriptionBy combining variational Monte Carlo (VMC) and completebasisset limit HartreeFock (HF) calculations, we have obtained nearexact correlation energies for lowdensity samespin electrons on a threedimensional sphere (3sphere), i.e., the surface of a fourdimensional ball. In the VMC calculations, we compare the efficacies of two types of oneelectron basis functions for these strongly correlated systems and analyze the energy convergence with respect to the quality of the Jastrow factor. The HF calculations employ spherical Gaussian functions (SGFs) which are the curvedspace analogs of Cartesian Gaussian functions. At low densities, the electrons become relatively localized into Wigner crystals, and the natural SGF centers are found by solving the Thomson problem (i.e., the minimumenergy arrangement of n point charges) on the 3sphere for various values of n. We have found 11 special values of n whose Thomson sites are equivalent. Three of these are the vertices of fourdimensional Platonic solids — the hypertetrahedron (n = 5), the hyperoctahedron (n = 8), and the 24cell (n = 24) — and a fourth is a highly symmetric structure (n = 13) which has not previously been reported. By calculating the harmonic frequencies of the electrons around their equilibrium positions, we also find the firstorder vibrational corrections to the Thomson energy.

Efficient algorithms for HirshfeldI charges
View Description Hide DescriptionA new viewpoint on iterative Hirshfeld charges is presented, whereby the atomic populations obtained from such a scheme are interpreted as such populations which reproduce themselves. This viewpoint yields a selfconsistent requirement for the HirshfeldI populations rather than being understood as the result of an iterative procedure. Based on this selfconsistent requirement, much faster algorithms for HirshfeldI charges have been developed. In addition, new atomic reference densities for the HirshfeldI procedure are presented. The proposed reference densities are Nrepresentable, display proper atomic shell structure and can be computed for any charged species.

Quantum Monte Carlo calculation of the binding energy of the beryllium dimer
View Description Hide DescriptionThe accurate calculation of the binding energy of the beryllium dimer is a challenging theoretical problem. In this study, the binding energy of Be2 is calculated using the diffusion Monte Carlo (DMC) method, using single Slater determinant and multiconfigurational trial functions. DMC calculations using singledeterminant trial wave functions of orbitals obtained from density functional theory calculations overestimate the binding energy, while DMC calculations using HartreeFock or CAS(4,8), complete active space trial functions significantly underestimate the binding energy. In order to obtain an accurate value of the binding energy of Be2 from DMC calculations, it is necessary to employ trial functions that include excitations outside the valence space. Our best estimate DMC result for the binding energy of Be2, obtained by using configuration interaction trial functions and extrapolating in the threshold for the configurations retained in the trial function, is 908 cm^{−1}, only slightly below the 935 cm^{−1} value derived from experiment.

Theory of bimolecular association dynamics in 2D for accurate model and experimental parameterization of binding rates
View Description Hide DescriptionThe dynamics of association between diffusing and reacting molecular species are routinely quantified using simple rateequation kinetics that assume both wellmixed concentrations of species and a single rate constant for parameterizing the binding rate. In twodimensions (2D), however, even when systems are wellmixed, the assumption of a single characteristic rate constant for describing association is not generally accurate, due to the properties of diffusional searching in dimensions d ≤ 2. Establishing rigorous bounds for discriminating between 2D reactive systems that will be accurately described by rate equations with a single rate constant, and those that will not, is critical for both modeling and experimentally parameterizing binding reactions restricted to surfaces such as cellular membranes. We show here that in regimes of intrinsic reaction rate (ka) and diffusion (D) parameters ka/D > 0.05, a single rate constant cannot be fit to the dynamics of concentrations of associating species independently of the initial conditions. Instead, a more sophisticated multiparametric description than rateequations is necessary to robustly characterize bimolecular reactions from experiment. Our quantitative bounds derive from our new analysis of 2D ratebehavior predicted from Smoluchowski theory. Using a recently developed single particle reactiondiffusion algorithm we extend here to 2D, we are able to test and validate the predictions of Smoluchowski theory and several other theories of reversible reaction dynamics in 2D for the first time. Finally, our results also mean that simulations of reactive systems in 2D using rate equations must be undertaken with caution when reactions have ka/D > 0.05, regardless of the simulation volume. We introduce here a simple formula for an adaptive concentration dependent rate constant for these chemical kinetics simulations which improves on existing formulas to better capture nonequilibrium reaction dynamics from dilute to dense systems.