Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package.

This article summarizes technical advances contained in the fifth major release of the Q-Chem quantum chemistry program package, covering developments since 2015. A comprehensive library of exchange-correlation functionals, along with a suite of correlated many-body methods, continues to be a hallmark of the Q-Chem software. The many-body methods include novel variants of both coupled-cluster and configuration-interaction approaches along with methods based on the algebraic diagrammatic construction and variational reduced density-matrix methods. Methods highlighted in Q-Chem 5 include a suite of tools for modeling core-level spectroscopy, methods for describing metastable resonances, methods for computing vibronic spectra, the nuclear-electronic orbital method, and several different energy decomposition analysis techniques. High-performance capabilities including multithreaded parallelism and support for calculations on graphics processing units are described. Q-Chem boasts a community of well over 100 active academic developers, and the continuing evolution of the software is supported by an "open teamware" model and an increasingly modular design.

Multi-state Approach to Chemical Reactivity in Fragment Based Quantum Chemistry Calculations.

We introduce a multistate framework for Fragment Molecular Orbital (FMO) quantum mechanical calculations and implement it in the context of protonated water clusters. The purpose of the framework is to address issues of nonuniqueness and dynamic fragmentation in FMO as well as other related fragment methods. We demonstrate that our new approach, Fragment Molecular Orbital Multistate Reactive Molecular Dynamics (FMO-MS-RMD), can improve energetic accuracy and yield stable molecular dynamics for small protonated water clusters undergoing proton transfer reactions.

Improving Generalized Born Models by Exploiting Connections to Polarizable Continuum Models. I. An Improved Effective Coulomb Operator.

We investigate the generalized Born (GB) implicit solvation model in comparison with polarizable continuum models (PCMs). We show that the GB model is intimately connected to the conductor-like PCM (C-PCM), a method that is accurate for high-dielectric solvents but less so for weakly polar and nonpolar solvents. The formal connection between C-PCM and the GB model suggests that C-PCM calculations place a limit on the accuracy that one should expect from GB models but also demonstrates that comparison of GB and C-PCM calculations directly interrogates the accuracy of the effective Coulomb operator that is used in the pairwise GB energy expression. We introduce a simple alternative to the "canonical" pairwise interaction operator first proposed by Still et al. and show that this alternative reduces the cost of the pairwise GB energy summation by as much as a factor of 3. At the same time, the new operator reduces statistical errors in solvation energies (as compared to C-PCM benchmarks) by 0.3% with respect to the canonical operator that exhibits an error of roughly 1.0%.

Improving Generalized Born Models by Exploiting Connections to Polarizable Continuum Models. II. Corrections for Salt Effects.

A previous analytical investigation of the generalized Born (GB) implicit solvation model is extended to solvents of nonzero ionic strength. The GB model with salt effects (GB-SE) is shown to resemble the Debye-Hückel-like screening model (DESMO), a polarizable continuum model (PCM) that we have recently developed for salty solutions. DESMO may be regarded either as a generalization of the conductor-like PCM (C-PCM) that extends C-PCM to electrolyte solutions or alternatively as a generalization of Debye-Hückel theory to arbitrary cavity shapes. The connection between GB-SE and DESMO suggests how the former can be modified to account for the exclusion of mobile ions from the cavity interior, an effect that is typically absent in GB-SE models. We propose two simple GB-SE models that are exact for a point charge in a spherical cavity and that introduce the ability to account, albeit approximately, for the finite size of the mobile ions. The accuracy of these new models is demonstrated by applications to both model systems and real proteins. These tests also demonstrate the accuracy of the DESMO approach, as compared to more sophisticated PCMs developed for electrolyte solutions.

A simple polarizable continuum solvation model for electrolyte solutions.

We propose a Debye-Hückel-like screening model (DESMO) that generalizes the familiar conductor-like screening model (COSMO) to solvents with non-zero ionic strength and furthermore provides a numerical generalization of the Debye-Hückel model that is applicable to non-spherical solute cavities. The numerical implementation of DESMO is based upon the switching/Gaussian (SWIG) method for smooth cavity discretization, which we have recently introduced in the context of polarizable continuum models (PCMs). This approach guarantees that the potential energy is a smooth function of the solute geometry and analytic gradients for DESMO are reported here. The SWIG formalism also facilitates analytic implementation of two other PCMs that are based on a screened Coulomb potential: the "integral equation formalism" (IEF-PCM) and the "surface and simulation of volume polarization for electrostatics" [SS(V)PE] method. Fully analytic implementations of these screened PCMs are reported here for the first time. Numerical results, for model systems where an exact solution of the linearized Poisson-Boltzmann equation is available, demonstrate that these screened PCMs are highly accurate. In realistic test cases, they are as accurate as the best available three-dimensional finite-difference methods. In polar solvents, DESMO is nearly as accurate as more sophisticated screened PCMs, but is significantly simpler and more efficient.

A smooth, nonsingular, and faithful discretization scheme for polarizable continuum models: the switching/Gaussian approach.

Polarizable continuum models (PCMs) are a widely used family of implicit solvent models based on reaction-field theory and boundary-element discretization of the solute/continuum interface. An often overlooked aspect of these theories is that discretization of the interface typically does not afford a continuous potential energy surface for the solute. In addition, we show that discretization can lead to numerical singularities and violations of exact variational conditions. To fix these problems, we introduce the switching/Gaussian (SWIG) method, a discretization scheme that overcomes several longstanding problems with PCMs. Our approach generalizes a procedure introduced by York and Karplus [J. Phys. Chem. A 103, 11060 (1999)], extending it beyond the conductor-like screening model. Comparison to other purportedly smooth PCM implementations reveals certain artifacts in these alternative approaches, which are avoided using the SWIG methodology. The versatility of our approach is demonstrated via geometry optimizations, vibrational frequency calculations, and molecular dynamics simulations, for solutes described using quantum mechanics and molecular mechanics.

Both intra- and interstrand charge-transfer excited states in aqueous B-DNA are present at energies comparable to, or just above, the (1)pipi* excitonic bright states.

Vertical electronic excitations in model systems representing single- and double-stranded B-DNA are characterized using electronic structure theory, including both time-dependent density functional theory (TD-DFT) and correlated wave function techniques. Previous TD-DFT predictions of charge-transfer (CT) states well below the optically bright (1)pipi* states are shown to be artifacts of the improper long-range behavior of standard density-functional exchange approximations, which we rectify here using a long-range correction (LRC) procedure. For nucleobase dimers (hydrogen-bonded or pi-stacked), TD-LRC-DFT affords vertical excitation energies in reasonable agreement with the wave function methods, not only for the (1)npi* and (1)pipi* states but also for the CT states, and qualitatively reproduces well-known base-stacking effects on the absorption spectrum of DNA. The emergence of (1)pipi* Frenkel exciton states, localized on a single strand, is clearly evident, and these states (rather than low-energy CT states) are primarily responsible for the fact that DNA's absorption spectrum exhibits a red tail that is absent in monomer absorption spectra. For B-DNA in aqueous solution, the low-energy tail of the CT band (representing both intra- and interstrand CT states) appears at energies comparable to those of the optically bright (1)pipi* exciton states. In systems with more than one base pair, we also observe the emergence of delocalized, interstrand CT excitations, whose excitation energies may be significantly lower than the lowest CT excitation in a single base pair. Together, these observations suggest that a single Watson-Crick base pair is an inadequate model of the photophysics of B-DNA.

Charge-transfer excited states in a pi-stacked adenine dimer, as predicted using long-range-corrected time-dependent density functional theory.

The lowest few electronic excitations of a pi-stacked adenine dimer in its B-DNA geometry are investigated, in the gas phase and in a water cluster, using a long-range-corrected version of time-dependent density functional theory (TD-DFT) that asymptotically incorporates Hartree-Fock exchange. Long-range correction is shown to eliminate the catastrophic underestimation of charge-transfer (CT) excitation energies that plagues conventional TD-DFT, at the expense of introducing one adjustable parameter, mu, that determines the length scale on which Hartree-Fock exchange is turned on. This parameter allows us to interpolate smoothly between hybrid density functionals and time-dependent Hartree-Fock theory. Excitation energies for CT states (in which an electron is transferred from one adenine molecule to the other) are found to increase dramatically as a function of mu. Uncorrected hybrid functionals underestimate the CT excitation energies, placing them well below the valence excitations, while time-dependent Hartree-Fock calculations place these states well above the valence states. Values for mu determined from certain benchmark calculations place the CT states well above the valence pipi* and npi* states at the Franck-Condon point.