Benjamin G. Levine

Conical intersections and double excitations in time-dependent density functional theory

There is a clear need for computationally inexpensive electronic structure theory methods which can model excited state potential energy surfaces. Time-dependent density functional theory (TDDFT) has emerged as one of the most promising contenders in this context. Many previous tests have concentrated on vertical excitation energies, which can be compared to experimental absorption maxima. Here, we focus attention on more global aspects of the resulting potential energy surfaces, especially conical intersections which play a key role in photochemical mechanisms. We introduce a new method for minimal energy conical intersection (MECI) searches which does not require knowledge of the non-adiabatic coupling vector. Using this new method, we compute MECI geometries with multi-state complete active space perturbation theory (MS-CASPT2) and TDDFT. We show that TDDFT in the linear response and adiabatic approximations can predict MECI geometries and energetics quite accurately, but that there are a number of qualitative deficiencies which need to be addressed before TDDFT can be used routinely in photochemical problems. †Dedicated to Professor M. A. Robb on the occasion of his 60th birthday.

Structure and mechanism of proton transport through the transmembrane tetrameric M2 protein bundle of the influenza A virus

The M2 proton channel from influenza A virus is an essential protein that mediates transport of protons across the viral envelope. This protein has a single transmembrane helix, which tetramerizes into the active channel. At the heart of the conduction mechanism is the exchange of protons between the His37 imidazole moieties of M2 and waters confined to the M2 bundle interior. Protons are conducted as the total charge of the four His37 side chains passes through 2+ and 3+ with a pKa near 6. A 1.65 Å resolution X-ray structure of the transmembrane protein (residues 25–46), crystallized at pH 6.5, reveals a pore that is lined by alternating layers of sidechains and well-ordered water clusters, which offer a pathway for proton conduction. The His37 residues form a box-like structure, bounded on either side by water clusters with well-ordered oxygen atoms at close distance. The conformation of the protein, which is intermediate between structures previously solved at higher and lower pH, suggests a mechanism by which conformational changes might facilitate asymmetric diffusion through the channel in the presence of a proton gradient. Moreover, protons diffusing through the channel need not be localized to a single His37 imidazole, but instead may be delocalized over the entire His-box and associated water clusters. Thus, the new crystal structure provides a possible unification of the discrete site versus continuum conduction models.

Ab initio multiple spawning dynamics using multi-state second-order perturbation theory

We have implemented multi-state second-order perturbation theory (MS-CASPT2) in the ab initio multiple spawning (AIMS) method for first-principles molecular dynamics including nonadiabatic effects. The nonadiabatic couplings between states are calculated numerically using an efficient method which requires only two extra energy calculations per time step. As a representative example, we carry out AIMS-MSPT2 calculations of the excited state dynamics of ethylene. Two distinct types of conical intersections, previously denoted as the twisted-pyramidalized and ethylidene intersections, are responsible for ultrafast population transfer from the excited state to the ground state. Although these two pathways have been observed in prior dynamics simulations, we show here that the branching ratio is affected by dynamic correlation with the twisted-pyramidalized intersection overweighting the ethylidene-like intersection during the decay process at the AIMS-MSPT2 level of description.

Functional studies and modeling of pore-lining residue mutants of the influenza A virus M2 ion channel

The A/M2 protein of influenza A virus forms a tetrameric proton-selective pH-gated ion channel. The H(37)xxxW(41) motif located in the channel pore is responsible for its gating and proton selectivity. Channel activation most likely involves protonation of the H37 residues, while the conductive state of the channel is characterized by two or three charged His residues in a tetrad. A/M2 channel activity is inhibited by the antiviral drug amantadine. Although a large number of functional amantadine-resistant mutants of A/M2 have been observed in vitro, only a few are observed in highly transmissible viruses in the presence or absence of amantadine. We therefore examined 49 point mutants of the pore-lining residues, representing both natural and nonnatural variants. Their ion selectivity, amantadine sensitivity, specific activity, and pH-dependent conductance were measured in Xenopus oocytes. These measurements showed how variations in the sequence lead to variations in the proton conduction. The results are consistent with a multistep mechanism that allows the protein to fine-tune its pH-rate profile over a wide range of proton concentrations, hypothesized to arise from different protonation states of the H37 tetrad. Mutations that give native-like conductance at low pH as well as minimal leakage current at pH 7.0 were surprisingly rare. Moreover, the results are consistent with a location of the amantadine-binding site inside the channel pore. These findings have helped to define the set of functionally fit mutants that should be targeted when considering the design of novel drugs that inhibit amantadine-resistant strains of influenza A virus.

Ab initio multiple spawning dynamics of excited butadiene: role of charge transfer.

Ab initio multiple spawning simulations of the photochemical reaction dynamics of s-trans-1,3-butadiene were performed. It is found that nonadiabatic events involving two low-lying excited states begin as early as 10 fs after excitation, resulting in the population being split between the bright 1(1)B(u) state and the dark 2(1)A(g) state. The molecule subsequently twists about a terminal carbon-carbon bond regardless of whether it is on the 1(1)B(u) or 2(1)A(g) electronic state. This twisting motion leads to conical intersections between S(1) and S(0). Several regions of the intersection seam involving states of differing character are accessed. The regions of the seam involving intersection between a state of charge-transfer character and a state of covalent character dominate the quenching dynamics, but intersections between two covalent states are also accessed a small percentage of the time. The existence and relative energies of these intersections are validated by optimization at the multistate complete active space second-order perturbation level of theory (MS-CASPT2). Our results point to a new mechanism for photoisomerization of butadiene that emphasizes the role of charge-transfer states.

Self-assembly of coarse-grained ionic surfactants accelerated by graphics processing units

Due to the relatively long time scales inherent to ionic surfactant self-assembly (>μs), an aggressive computational approach is needed to obtain converged data on micellar solutions. This work presents a study of micellization using a coarse-grained (CG) model of aqueous ionic surfactants in replicated molecular dynamics (MD) simulations run on graphics processing unit hardware. The performance of our implementation of the CG model with electrostatics into the HOOMD-Blue GPU-accelerated MD software package is comparable to that of a modest sized cluster running a highly optimized parallel CPU code. From 0.36 ms of cumulative trajectory data, we are able to predict equilibrium thermodynamic and morphological properties of ionic surfactant micellar solutions. Estimating the critical micelle concentrations (CMC) from the free monomer (ρ1) and premicellar concentrations obtained from simulations of sodium hexyl sulfate (S6S, CMC of 460 ± 6 mM) at high (1 M) concentration, a value in good agreement with experimental results is obtained; however, the same method applied to simulations of sodium nonyl sulfate (S9S, ρ1 of 2.4 ± 0.01 mM) and sodium dodecyl sulfate (SDS, ρ1 of 0.02 ± 0.01 mM) at the same total concentration systematically underestimates the CMCs. An alternative method for calculating the CMC is presented, where the free monomer concentration computed from high concentration CG-MD data is used as the input to a simple theoretical model which can be used to extrapolate to a more accurate prediction of the CMC. Better agreement between the empirical and predicted CMC is obtained from this theory for S9S (28.7 ± 0.3 mM) and SDS (3.32 ± 0.04 mM), though the CMC for S6S is slightly underestimated (304 ± 3 mM). We also present statistically converged morphological data, including aggregation number distributions and the principal components of the gyration tensor. This data suggest a transition from spherical micelles to rod-like at a specific aggregation number, which increases with increasing hydrocarbon length.

Micellization Studied by GPU-Accelerated Coarse-Grained Molecular Dynamics

The computational design of advanced materials based on surfactant self-assembly without ever stepping foot in the laboratory is an important goal, but there are significant barriers to this approach, because of the limited spatial and temporal scales accessible by computer simulations. In this paper, we report our work to bridge the gap between laboratory and computational time scales by implementing the coarse-grained (CG) force field previously reported by Shinoda et al. [Shinoda, W.; DeVane, R.; Klein, M. L. Mol. Simul. 2007, 33, 27-36] into the HOOMD-Blue graphical processing unit (GPU)-accelerated molecular dynamics (MD) software package previously reported by Anderson et al. [Anderson, J. A.; Lorenz, C. D.; Travesset, A. J. Comput. Phys. 2008, 227, 5342-5359]. For a system of 25 750 particles, this implementation provides performance on a single GPU, which is superior to that of a widely used parallel MD simulation code running on an optimally sized CPU-based cluster. Using our GPU setup, we have collected 0.6 ms of MD trajectory data for aqueous solutions of 7 different nonionic polyethylene glycol (PEG) surfactants, with most of the systems studied representing ∼1 000 000 atoms. From this data, we calculated various properties as a function of the length of the hydrophobic tails and PEG head groups. Specifically, we determined critical micelle concentrations (CMCs), which are in good agreement with experimental data, and characterized the size and shape of micelles. However, even with the microsecond trajectories employed in this study, we observed that the micelles composed of relatively hydrophobic surfactants are continuing to grow at the end of our simulations. This suggests that the final micelle size distributions of these systems are strongly dependent on initial conditions and that either longer simulations or advanced sampling techniques are needed to properly sample their equilibrium distributions. Nonetheless, the combination of coarse-grained modeling and GPU acceleration marks a significant step toward the computational prediction of the thermodynamic properties of slowly evolving surfactant systems.

Nanoscale Multireference Quantum Chemistry: Full Configuration Interaction on Graphical Processing Units

Methods based on a full configuration interaction (FCI) expansion in an active space of orbitals are widely used for modeling chemical phenomena such as bond breaking, multiply excited states, and conical intersections in small-to-medium-sized molecules, but these phenomena occur in systems of all sizes. To scale such calculations up to the nanoscale, we have developed an implementation of FCI in which electron repulsion integral transformation and several of the more expensive steps in σ vector formation are performed on graphical processing unit (GPU) hardware. When applied to a 1.7 × 1.4 × 1.4 nm silicon nanoparticle (Si72H64) described with the polarized, all-electron 6-31G** basis set, our implementation can solve for the ground state of the 16-active-electron/16-active-orbital CASCI Hamiltonian (more than 100,000,000 configurations) in 39 min on a single NVidia K40 GPU.

On the extent and connectivity of conical intersection seams and the effects of three-state intersections.

We discuss the connectivity of intersection spaces and the role of minimal energy points within these intersection spaces (minimal energy conical intersections or MECIs) in promoting nonadiabatic transitions. We focus on malonaldeyde as a specific example, where there is a low-lying three-state conical intersection. This three-state intersection is the global minimum on the bright excited electronic state, but it plays a limited role in population transfer in our ab initio multiple spawning (AIMS) simulations because the molecule must traverse a series of two-state conical intersections to reach the three-state intersection. Due to the differences in seam space dimensionality separating conventional (two-state) and three-state intersections, we suggest that dynamical effects arising directly from a three-state intersection may prove difficult to observe in general. We also use a newly developed method for intersection optimization with geometric constraints to demonstrate the connectivity of all the stationary points in the intersection spaces for malonaldehyde. This supports the conjecture that all intersection spaces are connected, and that three-state intersections play a key role in extending this connectivity to all pairs of states, e.g. the S1/S0 and S2/S1 intersection spaces.

Evaluation of the Time-Derivative Coupling for Accurate Electronic State Transition Probabilities from Numerical Simulations.

Spikes in the time-derivative coupling (TDC) near surface crossings make the accurate integration of the time-dependent Schrödinger equation in nonadiabatic molecular dynamics simulations a challenge. To address this issue, we present an approximation to the TDC based on a norm-preserving interpolation (NPI) of the adiabatic electronic wave functions within each time step. We apply NPI and two other schemes for computing the TDC in numerical simulations of the Landau-Zener model, comparing the simulated transfer probabilities to the exact solution. Though NPI does not require the analytical calculation of nonadiabatic coupling matrix elements, it consistently yields unsigned population transfer probability errors of ∼0.001, whereas analytical calculation of the TDC yields errors of 0.0-1.0 depending on the time step, the offset of the maximum in the TDC from the beginning of the time step, and the coupling strength. The approximation of Hammes-Schiffer and Tully yields errors intermediate between NPI and the analytical scheme.