Volume 142, Issue 3, 21 January 2015
Index of content:

To perform its barrier function, the lipid bilayer membrane requires a robust resistance against pore formation. Using a selfconsistent field (SCF) theory and a molecularly detailed model for membranes composed of charged or zwitterionic lipids, it is possible to predict structural, mechanical, and thermodynamical parameters for relevant lipid bilayer membranes. We argue that the edge energy in membranes is a function of the spontaneous lipid monolayer curvature, the mean bending modulus, and the membrane thickness. An analytical Helfrichlike model suggests that most bilayers should have a positive edge energy. This means that there is a natural resistance against pore formation. Edge energies evaluated explicitly in a twogradient SCF model are consistent with this. Remarkably, the edge energy can become negative for phosphatidylglycerol (e.g., dioleoylphosphoglycerol) bilayers at a sufficiently low ionic strength. Such bilayers become unstable against the formation of pores or the formation of lipid disks. In the weakly curved limit, we study the curvature dependence of the edge energy and evaluate the preferred edge curvature and the edge bending modulus. The latter is always positive, and the former increases with increasing ionic strength. These results point to a small window of ionic strengths for which stable pores can form as too low ionic strengths give rise to lipid disks. Higher order curvature terms are necessary to accurately predict relevant pore sizes in bilayers. The electric double layer overlap across a small pore widens the window of ionic strengths for which pores are stable.
 ARTICLES

 Theoretical Methods and Algorithms

On the edge energy of lipid membranes and the thermodynamic stability of pores
View Description Hide DescriptionTo perform its barrier function, the lipid bilayer membrane requires a robust resistance against pore formation. Using a selfconsistent field (SCF) theory and a molecularly detailed model for membranes composed of charged or zwitterionic lipids, it is possible to predict structural, mechanical, and thermodynamical parameters for relevant lipid bilayer membranes. We argue that the edge energy in membranes is a function of the spontaneous lipid monolayer curvature, the mean bending modulus, and the membrane thickness. An analytical Helfrichlike model suggests that most bilayers should have a positive edge energy. This means that there is a natural resistance against pore formation. Edge energies evaluated explicitly in a twogradient SCF model are consistent with this. Remarkably, the edge energy can become negative for phosphatidylglycerol (e.g., dioleoylphosphoglycerol) bilayers at a sufficiently low ionic strength. Such bilayers become unstable against the formation of pores or the formation of lipid disks. In the weakly curved limit, we study the curvature dependence of the edge energy and evaluate the preferred edge curvature and the edge bending modulus. The latter is always positive, and the former increases with increasing ionic strength. These results point to a small window of ionic strengths for which stable pores can form as too low ionic strengths give rise to lipid disks. Higher order curvature terms are necessary to accurately predict relevant pore sizes in bilayers. The electric double layer overlap across a small pore widens the window of ionic strengths for which pores are stable.

The abinitio density matrix renormalization group in practice
View Description Hide DescriptionThe abinitio density matrix renormalization group (DMRG) is a tool that can be applied to a wide variety of interesting problems in quantum chemistry. Here, we examine the density matrix renormalization group from the vantage point of the quantum chemistry user. What kinds of problems is the DMRG wellsuited to? What are the largest systems that can be treated at practical cost? What sort of accuracies can be obtained, and how do we reason about the computational difficulty in different molecules? By examining a diverse benchmark set of molecules: πelectron systems, benchmark maingroup and transition metal dimers, and the Mnoxosalen and Feporphine organometallic compounds, we provide some answers to these questions, and show how the density matrix renormalization group is used in practice.

Hybrid pathwise sensitivity methods for discrete stochastic models of chemical reaction systems
View Description Hide DescriptionStochastic models are often used to help understand the behavior of intracellular biochemical processes. The most common such models are continuous time Markov chains (CTMCs). Parametric sensitivities, which are derivatives of expectations of model output quantities with respect to model parameters, are useful in this setting for a variety of applications. In this paper, we introduce a class of hybrid pathwise differentiation methods for the numerical estimation of parametric sensitivities. The new hybrid methods combine elements from the three main classes of procedures for sensitivity estimation and have a number of desirable qualities. First, the new methods are unbiased for a broad class of problems. Second, the methods are applicable to nearly any physically relevant biochemical CTMC model. Third, and as we demonstrate on several numerical examples, the new methods are quite efficient, particularly if one wishes to estimate the full gradient of parametric sensitivities. The methods are rather intuitive and utilize the multilevel Monte Carlo philosophy of splitting an expectation into separate parts and handling each in an efficient manner.

Towards an exact theory of linear absorbance and circular dichroism of pigmentprotein complexes: Importance of nonsecular contributions
View Description Hide DescriptionA challenge for the theory of optical spectra of pigmentprotein complexes is the equal strength of the pigmentpigment and the pigmentprotein couplings. Treating both on an equal footing so far can only be managed by numerically costly approaches. Here, we exploit recent results on a normal mode analysis derived spectral density that revealed the dominance of the diagonal matrix elements of the excitonvibrational coupling in the exciton state representation. We use a cumulant expansion technique that treats the diagonal parts exactly, includes an infinite summation of the offdiagonal parts in secular and Markov approximations, and provides a systematic perturbative way to include nonsecular and nonMarkov corrections. The theory is applied to a model dimer and to chlorophyll (Chl) a and Chl b homodimers of the reconstituted watersoluble chlorophyllbinding protein (WSCP) from cauliflower. The model calculations reveal that the nonsecular/nonMarkov effects redistribute oscillator strength from the strong to the weak exciton transition in absorbance and they diminish the rotational strength of the exciton transitions in circular dichroism. The magnitude of these corrections is in a few percent range of the overall signal, providing a quantitative explanation of the success of timelocal convolutionless density matrix theory applied earlier. A close examination of the optical spectra of Chl a and Chl b homodimers in WSCP suggests that the opening angle between Qy transition dipole moments in Chl b homodimers is larger by about 9^{∘} than for Chl a homodimers for which a crystal structure of a related WSCP complex exists. It remains to be investigated whether this change is due to a different mutual geometry of the pigments or due to the different electronic structures of Chl a and Chl b.

Shortrange stabilizing potential for computing energies and lifetimes of temporary anions with extrapolation methods
View Description Hide DescriptionThe energy of a temporary anion can be computed by adding a stabilizing potential to the molecular Hamiltonian, increasing the stabilization until the temporary state is turned into a bound state, and then further increasing the stabilization until enough bound state energies have been collected so that these can be extrapolated back to vanishing stabilization. The lifetime can be obtained from the same data, but only if the extrapolation is done through analytic continuation of the momentum as a function of the square root of a shifted stabilizing parameter. This method is known as analytic continuation of the coupling constant, and it requires—at least in principle—that the boundstate input data are computed with a shortrange stabilizing potential. In the context of molecules and ab initio packages, longrange Coulomb stabilizing potentials are, however, far more convenient and have been used in the past with some success, although the error introduced by the longrang nature of the stabilizing potential remains unknown. Here, we introduce a softVoronoi box potential that can serve as a shortrange stabilizing potential. The difference between a Coulomb and the new stabilization is analyzed in detail for a onedimensional model system as well as for the ^{2}Π u resonance of CO , and in both cases, the extrapolation results are compared to independently computed resonance parameters, from complex scaling for the model, and from complex absorbing potential calculations for CO . It is important to emphasize that for both the model and for CO , all three sets of results have, respectively, been obtained with the same electronic structure method and basis set so that the theoretical description of the continuum can be directly compared. The new softVoronoiboxbased extrapolation is then used to study the influence of the size of diffuse and the valence basis sets on the computed resonance parameters.

Sublinear scaling for timedependent stochastic density functional theory
View Description Hide DescriptionA stochastic approach to timedependent density functional theory is developed for computing the absorption cross section and the random phase approximation (RPA) correlation energy. The core idea of the approach involves timepropagation of a small set of stochastic orbitals which are first projected on the occupied space and then propagated in time according to the timedependent KohnSham equations. The evolving electron density is exactly represented when the number of random orbitals is infinite, but even a small number (≈16) of such orbitals is enough to obtain meaningful results for absorption spectrum and the RPA correlation energy per electron. We implement the approach for silicon nanocrystals using realspace grids and find that the overall scaling of the algorithm is sublinear with computational time and memory.

Deviations from piecewise linearity in the solidstate limit with approximate density functionals
View Description Hide DescriptionIn exact density functional theory, the total groundstate energy is a series of linear segments between integer electron points, a condition known as “piecewise linearity.” Deviation from this condition is indicative of poor predictive capabilities for electronic structure, in particular of ionization energies, fundamental gaps, and charge transfer. In this article, we take a new look at the deviation from linearity (i.e., curvature) in the solidstate limit by considering two different ways of approaching it: a large finite system of increasing size and a crystal represented by an increasingly large reference cell with periodic boundary conditions. We show that the curvature approaches vanishing values in both limits, even for functionals which yield poor predictions of electronic structure, and therefore cannot be used as a diagnostic or constructive tool in solids. We find that the approach towards zero curvature is different in each of the two limits, owing to the presence of a compensating background charge in the periodic case. Based on these findings, we present a new criterion for functional construction and evaluation, derived from the sizedependence of the curvature, along with a practical method for evaluating this criterion. For large finite systems, we further show that the curvature is dominated by the selfinteraction of the highest occupied eigenstate. These findings are illustrated by computational studies of various solids, semiconductor nanocrystals, and long alkane chains.

A finiteelement visualization of quantum reactive scattering. II. Nonadiabaticity on coupled potential energy surfaces
View Description Hide DescriptionThis is the second in a series of papers detailing a MATLAB based implementation of the finite element method applied to collinear triatomic reactions. Here, we extend our previous work to reactions on coupled potential energy surfaces. The divergence of the probability current density field associated with the two electronically adiabatic states allows us to visualize in a novel way where and how nonadiabaticity occurs. A twodimensional investigation gives additional insight into nonadiabaticity beyond standard onedimensional models. We study the F(^{2}P) + HCl and F(^{2}P) + H2 reactions as model applications. Our publicly available code (http://www2.chem.umd.edu/groups/alexander/FEM) is general and easy to use.

On the accuracy of coherent modified Redfield theory in simulating excitation energy transfer dynamics
View Description Hide DescriptionIn this study, we investigate the accuracy of a recently developed coherent modified Redfield theory (CMRT) in simulating excitation energy transfer (EET) dynamics. The CMRT is a secular nonMarkovian quantum master equation that is derived by extending the modified Redfield theory to treat coherence dynamics in molecular excitonic systems. Herein, we systematically survey the applicability of the CMRT in a large EET parameter space through the comparisons of the CMRT EET dynamics in a dimer system with the numerically exact results. The results confirm that the CMRT exhibits a broad applicable range and allow us to locate the specific parameter regimes where CMRT fails to provide adequate results. Moreover, we propose an accuracy criterion based on the magnitude of secondorder perturbation to characterize the applicability of CMRT and show that the criterion summarizes all the benchmark results and the physics described by CMRT. Finally, we employ the accuracy criterion to quantitatively compare the performance of CMRT to that of a small polaron quantum master equation approach. The comparison demonstrates the complementary nature of these two methods, and as a result, the combination of the two methods provides accurate simulations of EET dynamics for the full parameter space investigated in this study. Our results not only delicately evaluate the applicability of the CMRT but also reveal new physical insights for factors controlling the dynamics of EET that should be useful for developing more accurate and efficient methods for simulations of EET dynamics in molecular aggregate systems.

Accuracy of maximum likelihood estimates of a twostate model in singlemolecule FRET
View Description Hide DescriptionPhoton sequences from singlemolecule Förster resonance energy transfer (FRET) experiments can be analyzed using a maximum likelihood method. Parameters of the underlying kinetic model (FRET efficiencies of the states and transition rates between conformational states) are obtained by maximizing the appropriate likelihood function. In addition, the errors (uncertainties) of the extracted parameters can be obtained from the curvature of the likelihood function at the maximum. We study the standard deviations of the parameters of a twostate model obtained from photon sequences with recorded colors and arrival times. The standard deviations can be obtained analytically in a special case when the FRET efficiencies of the states are 0 and 1 and in the limiting cases of fast and slow conformational dynamics. These results are compared with the results of numerical simulations. The accuracy and, therefore, the ability to predict model parameters depend on how fast the transition rates are compared to the photon count rate. In the limit of slow transitions, the key parameters that determine the accuracy are the number of transitions between the states and the number of independent photon sequences. In the fast transition limit, the accuracy is determined by the small fraction of photons that are correlated with their neighbors. The relative standard deviation of the relaxation rate has a “chevron” shape as a function of the transition rate in the loglog scale. The location of the minimum of this function dramatically depends on how well the FRET efficiencies of the states are separated.

Selfconsistent continuum solvation for optical absorption of complex molecular systems in solution
View Description Hide DescriptionWe introduce a new method to compute the optical absorption spectra of complex molecular systems in solution, based on the Liouville approach to timedependent densityfunctional perturbation theory and the revised selfconsistent continuum solvation model. The former allows one to obtain the absorption spectrum over a whole wide frequency range, using a recently proposed Lanczosbased technique, or selected excitation energies, using the Casida equation, without having to ever compute any unoccupied molecular orbitals. The latter is conceptually similar to the polarizable continuum model and offers the further advantages of allowing an easy computation of atomic forces via the HellmannFeynman theorem and a ready implementation in periodicboundary conditions. The new method has been implemented using pseudopotentials and planewave basis sets, benchmarked against polarizable continuum model calculations on 4aminophthalimide, alizarin, and cyanin and made available through the Quantum ESPRESSO distribution of opensource codes.

Stabilized quasiNewton optimization of noisy potential energy surfaces
View Description Hide DescriptionOptimizations of atomic positions belong to the most commonly performed tasks in electronic structure calculations. Many simulations like global minimum searches or characterizations of chemical reactions require performing hundreds or thousands of minimizations or saddle computations. To automatize these tasks, optimization algorithms must not only be efficient but also very reliable. Unfortunately, computational noise in forces and energies is inherent to electronic structure codes. This computational noise poses a severe problem to the stability of efficient optimization methods like the limitedmemory Broyden–Fletcher–Goldfarb–Shanno algorithm. We here present a technique that allows obtaining significant curvature information of noisy potential energy surfaces. We use this technique to construct both, a stabilized quasiNewton minimization method and a stabilized quasiNewton saddle finding approach. We demonstrate with the help of benchmarks that both the minimizer and the saddle finding approach are superior to comparable existing methods.

Fundamental aspects of recoupled pair bonds. I. Recoupled pair bonds in carbon and sulfur monofluoride
View Description Hide DescriptionThe number of singly occupied orbitals in the groundstate atomic configuration of an element defines its nominal valence. For carbon and sulfur, with two singly occupied orbitals in their ^{3}P ground states, the nominal valence is two. However, in both cases, it is possible to form more bonds than indicated by the nominal valence—up to four bonds for carbon and six bonds for sulfur. In carbon, the electrons in the 2s lone pair can participate in bonding, and in sulfur the electrons in both the 3p and 3s lone pairs can participate. Carbon 2s and sulfur 3p recoupled pair bonds are the basis for the tetravalence of carbon and sulfur, and 3s recoupled pair bonds enable sulfur to be hexavalent. In this paper, we report generalized valence bond as well as more accurate calculations on the a ^{4}Σ^{−} states of CF and SF, which are archetypal examples of molecules that possess recoupled pair bonds. These calculations provide insights into the fundamental nature of recoupled pair bonds and illustrate the key differences between recoupled pair bonds formed with the 2s lone pair of carbon, as a representative of the early pblock elements, and recoupled pair bonds formed with the 3p lone pair of sulfur, as a representative of the late pblock elements.

Fundamental aspects of recoupled pair bonds. II. Recoupled pair bond dyads in carbon and sulfur difluoride
View Description Hide DescriptionFormation of a bond between a second ligand and a molecule with a recoupled pair bond results in a recoupled pair bond dyad. We examine the recoupled pair bond dyads in the a ^{3}B1 states of CF2 and SF2, which are formed by the addition of a fluorine atom to the a ^{4}Σ^{−} states of CF and SF, both of which possess recoupled pair bonds. The two dyads are very different. In SF2, the second FS–F bond is very strong (De = 106.3 kcal/mol), the bond length is much shorter than that in the SF(a ^{4}Σ^{−}) state (1.666 Å versus 1.882 Å), and the three atoms are nearly collinear (θe = 162.7°) with only a small barrier to linearity (0.4 kcal/mol). In CF2, the second FC–F bond is also very strong (De = 149.5 kcal/mol), but the bond is only slightly shorter than that in the CF(a ^{4}Σ^{−}) state (1.314 Å versus 1.327 Å), and the molecule is strongly bent (θe = 119.0°) with an 80.5 kcal/mol barrier to linearity. The a ^{3}B1 states of CF2 and SF2 illustrate the fundamental differences between recoupled pair bond dyads formed from 2s and 3p lone pairs.

NonMarkovian Quantum State Diffusion for temperaturedependent linear spectra of light harvesting aggregates
View Description Hide DescriptionNonMarkovian Quantum State Diffusion (NMQSD) has turned out to be an efficient method to calculate excitonic properties of aggregates composed of organic chromophores, taking into account the coupling of electronic transitions to vibrational modes of the chromophores. NMQSD is an open quantum system approach that incorporates environmental degrees of freedom (the vibrations in our case) in a stochastic way. We show in this paper that for linear optical spectra (absorption, circular dichroism), no stochastics is needed, even for finite temperatures. Thus, the spectra can be obtained by propagating a single trajectory. To this end, we map a finite temperature environment to the zero temperature case using the socalled thermofield method. The resulting equations can then be solved efficiently by standard integrators.

Twocomponent hybrid timedependent density functional theory within the TammDancoff approximation
View Description Hide DescriptionWe report the implementation of a twocomponent variant of timedependent density functional theory (TDDFT) for hybrid functionals that accounts for spinorbit effects within the TammDancoff approximation (TDA) for closedshell systems. The influence of the admixture of HartreeFock exchange on excitation energies is investigated for several atoms and diatomic molecules by comparison to numbers for pure density functionals obtained previously [M. Kühn and F. Weigend, J. Chem. Theory Comput. 9, 5341 (2013)]. It is further related to changes upon switching to the local density approximation or using the full TDDFT formalism instead of TDA. Efficiency is demonstrated for a comparably large system, Ir(ppy)3 (61 atoms, 1501 basis functions, lowest 10 excited states), which is a prototype molecule for organic lightemitting diodes, due to its “spinforbidden” tripletsinglet transition.

Arbitrary order permanent Cartesian multipolar electrostatic interactions
View Description Hide DescriptionRecently, there has been a concerted effort to implement advanced classical potential energy surfaces by adding higher order multipoles to fixed point charge electrostatics in a bid to increase the accuracy of simulations of condensed phase systems. One major hurdle is the unwieldy nature of the expressions which in part has limited developers mostly to including only dipoles and quadrupoles. In this paper, we present a generalization of the Cartesian formulation of electrostatic multipolar interactions that enables the specification of an arbitrary order of multipoles. Specifically, we derive formulas for arbitrary order implementation of the particle mesh Ewald method and give a closed form formula for the stress tensor in the reciprocal space. In addition, we provide recurrence relations for common electrostatic potentials employed in molecular simulations, which allows for the generalization to arbitrary order and guarantees a computational cost that scales as O(p ^{3}) for Cartesian multipole interactions of order p.

Adaptive hybrid simulations for multiscale stochastic reaction networks
View Description Hide DescriptionThe probability distribution describing the state of a Stochastic Reaction Network (SRN) evolves according to the Chemical Master Equation (CME). It is common to estimate its solution using Monte Carlo methods such as the Stochastic Simulation Algorithm (SSA). In many cases, these simulations can take an impractical amount of computational time. Therefore, many methods have been developed that approximate sample paths of the underlying stochastic process and estimate the solution of the CME. A prominent class of these methods include hybrid methods that partition the set of species and the set of reactions into discrete and continuous subsets. Such a partition separates the dynamics into a discrete and a continuous part. Simulating such a stochastic process can be computationally much easier than simulating the exact discrete stochastic process with SSA. Moreover, the quasistationary assumption to approximate the dynamics of fast subnetworks can be applied for certain classes of networks. However, as the dynamics of a SRN evolves, these partitions may have to be adapted during the simulation. We develop a hybrid method that approximates the solution of a CME by automatically partitioning the reactions and species sets into discrete and continuous components and applying the quasistationary assumption on identifiable fast subnetworks. Our method does not require any user intervention and it adapts to exploit the changing timescale separation between reactions and/or changing magnitudes of copynumbers of constituent species. We demonstrate the efficiency of the proposed method by considering examples from systems biology and showing that very good approximations to the exact probability distributions can be achieved in significantly less computational time. This is especially the case for systems with oscillatory dynamics, where the system dynamics change considerably throughout the timeperiod of interest.

Molecular quantum mechanical gradients within the polarizable embedding approach—Application to the internal vibrational Stark shift of acetophenone
View Description Hide DescriptionWe present an implementation of analytical quantum mechanical molecular gradients within the polarizable embedding (PE) model to allow for efficient geometry optimizations and vibrational analysis of molecules embedded in large, geometrically frozen environments. We consider a variational ansatz for the quantum region, covering (multiconfigurational) selfconsistentfield and Kohn–Sham density functional theory. As the first application of the implementation, we consider the internal vibrational Stark effect of the C=O group of acetophenone in different solvents and derive its vibrational linear Stark tuning rate using harmonic frequencies calculated from analytical gradients and computed local electric fields. Comparisons to PE calculations employing an enlarged quantum region as well as to a nonpolarizable embedding scheme show that the inclusion of mutual polarization between acetophenone and water is essential in order to capture the structural modifications and the associated frequency shifts observed in water. For more apolar solvents, a proper description of dispersion and exchange–repulsion becomes increasingly important, and the quality of the optimized structures relies to a larger extent on the quality of the LennardJones parameters.

Timedependent nonequilibrium dielectric response in QM/continuum approaches
View Description Hide DescriptionThe Polarizable Continuum Models (PCMs) are some of the most inexpensive yet successful methods for including the effects of solvation in quantummechanical calculations of molecular systems. However, when applied to the electronic excitation process, these methods are restricted to dichotomously assuming either that the solvent has completely equilibrated with the excited solute charge density (infinitetime limit), or that it retains the configuration that was in equilibrium with the solute prior to excitation (zerotime limit). This renders the traditional PCMs inappropriate for resolving timedependent solvent effects on nonequilibrium solute electron dynamics like those implicated in the instants following photoexcitation of a solvated molecular species. To extend the existing methods to this nonequilibrium regime, we herein derive and apply a new formalism for a general timedependent continuum embedding method designed to be propagated alongside the solute’s electronic degrees of freedom in the time domain. Given the frequencydependent dielectric constant of the solvent, an equation of motion for the dielectric polarization is derived within the PCM framework and numerically integrated simultaneously with the timedependent Hartree fock/density functional theory equations. Results for small molecular systems show the anticipated dipole quenching and electronic state dephasing/relaxation resulting from outofphase charge fluctuations in the dielectric and embedded quantum system.