Article pubs.acs.org/JPCB
Atomic Spectral Methods for Ab Initio Molecular Electronic Energy Surfaces: Transitioning From Small-Molecule to BiomolecularSuitable Approaches Jeffrey D. Mills,† Michal Ben-Nun,‡ Kyle Rollin,§ Michael W. J. Bromley,∥ Jiabo Li,⊥ Robert J. Hinde,# Carl L. Winstead,∇ Jeffrey A. Sheehy,○ Jerry A. Boatz,† and Peter W. Langhoff*,◆ †
Air Force Research Laboratory, 10 East Saturn Boulevard, Edwards AFB, California 93524-7680, United States Predictive Science, Inc., 9990 Mesa Rim Road #170, San Diego, California 92121, United States § Northrup Grumman Corporation, 1 Rancho Carmel Drive, San Diego, California 92128, United States ∥ School of Mathematics and Physics, The University of Queensland, Brisbane, Queensland 4072, Australia ⊥ Accelrys Inc., 10188 Telesis Court #100, San Diego, California 92121-4779, United States # Department of Chemistry, University of Tennessee, Knoxville, Tennessee 37996-1600, United States ∇ A.A. Noyes Laboratory of Chemical Physics, California Institute of Technology, Pasadena, California 91125, United States ○ NASA Headquarters, 300 E Street SW, Suite 5R30, Washington, DC 201546, United States ◆ Department of Chemistry and Biochemistry, University of California, San Diego, 9500 Gilman Drive, MS 0365, La Jolla, California 92093-0365, United States ‡
ABSTRACT: Continuing attention has addressed incorportation of the electronically dynamical attributes of biomolecules in the largely static first-generation molecular-mechanical force fields commonly employed in molecular-dynamics simulations. We describe here a universal quantummechanical approach to calculations of the electronic energy surfaces of both small molecules and large aggregates on a common basis which can include such electronic attributes, and which also seems well-suited to adaptation in ab initio molecular-dynamics applications. In contrast to the more familiar orbital-product-based methodologies employed in traditional small-molecule computational quantum chemistry, the present approach is based on an “ex-post-facto” method in which Hamiltonian matrices are evaluated prior to wave function antisymmetrization, implemented here in the support of a Hilbert space of orthonormal products of many-electron atomic spectral eigenstates familiar from the van der Waals theory of long-range interactions. The general theory in its various forms incorporates the early semiempirical atoms- and diatomics-in-molecules approaches of Moffitt, Ellison, Tully, Kuntz, and others in a comprehensive mathematical setting, and generalizes the developments of Eisenschitz, London, Claverie, and others addressing electron permutation symmetry adaptation issues, completing these early attempts to treat van der Waals and chemical forces on a common basis. Exact expressions are obtained for molecular Hamiltonian matrices and for associated energy eigenvalues as sums of separate atomic and interaction-energy terms, similar in this respect to the forms of classical force fields. The latter representation is seen to also provide a long-missing general definition of the energies of individual atoms and of their interactions within molecules and matter free from subjective additional constraints. A computer code suite is described for calculations of the many-electron atomic eigenspectra and the pairwise-atomic Hamiltonian matrices required for practical applications. These matrices can be retained as functions of scalar atomic-pair separations and employed in assembling aggregate Hamiltonian matrices, with Wigner rotation matrices providing analytical representations of their angular degrees of freedom. In this way, ab initio potential energy surfaces are obtained in the complete absence of repeated evaluations and transformations of the one- and two-electron integrals at different molecular geometries required in most ab inito molecular calculations, with large Hamiltonian matrix assembly simplified and explicit diagonalizations avoided employing partitioning and Brillouin−Wigner or Rayleigh−Schrödinger perturbation theory. Illustrative applications of the important components of the formalism, selected aspects of the scaling of the approach, and aspects of “onthe-fly” interfaces with Monte Carlo and molecular-dynamics methods are described in anticipation of subsequent applications to biomolecules and other large aggregates.
1. INTRODUCTION
Special Issue: J. Andrew McCammon Festschrift
Molecular-dynamics simulations of the structures and dynamical properties of biomolecules continue to provide one of the most useful and widely employed methodologies generally © XXXX American Chemical Society
Received: February 26, 2016 Revised: May 4, 2016
A
DOI: 10.1021/acs.jpcb.6b02021 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
providing limited quantitative information relevant to other systems of interest. The dependence upon orbital representations of wave-function-based approaches, and the resulting consequences, arguably present intrinsic computational limitations to the aforementioned methods that can not easily be overcome, even by the ever-increasing capabilities of computer resources available for their application. An attractive alternative approach, potentially capable of overcoming the aforementioned limitations, is based on enforcing the antisymmetry requirements of Pauli’s principle “ex-post-facto” by unitary transformation of a Hamiltonian matrix which is considerably simpler to evaluate than those of conventional developments, also offering additional flexibility in choice of forms of product representations employed for this purpose.36−48 In the case of simple atomic- or molecular-orbital representations in Hartree-product forms, for example, it has been shown that results identical with conventional use of Slater determinants are obtained, employing either finite or arbitrarily large orbital basis sets,44 and adding confidence to the approach, but not providing any particular computational advantage in orbital-product cases. By contrast, use of Hartreelike products of universal many-electron atomic eigenstates as an orthonormal representation for purposes of analysis and for support suitable for molecular computations, as described in the original formulation of the method,36,37 has led to a continuing series of investigations in clarification of aspects of the theory,39,40,42−44,46 and for purposes of implementation.38,41,45,47,48 These developments also incorporate in a unified ab initio formalism earlier attempts to employ atomicbased representations in combined descriptions of van der Waals and chemical forces,49 semiempirical “atoms-in-molecules” calculations of chemical bond strengths,50 adoption of atomic-pair energies in semiempirical “diatomics-in-molecules” approximations to potential energy surfaces,51 and clarification of the role of electron permutation symmetry in long-range atomic interactions.52,53 The present contribution to this J. Andrew McCammon Festschrift describes the ex-post-facto methodology in the “atomic spectral-product” form, and a computer code suite for implementation of electronic structure calculations potentially suitable for applications to biomolecules and other large aggregates. A brief recapitulation of the essential features of the underlying general theory, detailed description of the methodology and code suite, selected calculations performed to illustrate the components of the code, and issues related to combined electronic energy surface calculations and Monte Carlo and molecular-dynamics applications are described. The theoretical formalism is reported in section 2 employing partitionings of the molecular Hamiltonian operator in forms suited to atoms-in-molecules and diatomics-in-molecules versions of the development; conditions are indicated under which the aggregate Hamiltonian matrix factors in a manner that requires only enforcement of atomic pairwise antisymmetry for calculations of potential energy surfaces in orthonormal atomic spectral-product representations, and the specific forms of the molecular Hamiltonian matrices obtained are reported. Section 3 describes the computational methodology devised and the code suite assembled for performing atomic pairwise spectral-product calculations. Symmetric-group methods are reported which have been employed to provide so-called phaseconsistent standard atomic tableau functions and molecular representations in subduced atomic-product forms,52,53 which are well-suited to calculations of the atomic and atomic-pair
available to the computational biochemistry community. Although the earliest of such studies employed largely “static” molecular-mechanics potential energy surfaces,1−6 continuing improvements and inovations have led to incorporation of additional selected physical and chemical attributes of biomolecules,7,8 largely through adoption of combinations of ad hoc procedures and ab initio quantum-mechanical methods.9−15 It spite of these promising advances, a generally accepted and widely available “ab inito molecular dynamics” methodology for treating even the lowest-lying vibronic states of an entire protein or other macromolecule, including the effects of electronic charge polarization, fluctuation, and transfer, dynamical reprotonation of titratable residues, selected nonadiabatic behaviors, and other modifications of firstgeneration force fields, does not yet appear to be in evidence.16 Although continuing advances in computational hardware and software for this purpose will play an important role in further improvements,17 new possibly ameliorating approaches are always welcome. Small-molecule electronic structure approaches employing largely conventional methods based on mature computational algorithms would seem to provide a significant arsenal of candidates to choose from in addressing their transitions to biomolecular applications. Choices of methodologies,18−35 aspects of which have already been adopted to some degree for biomolecular applications, include highly efficient Hartree− Fock and density-functional methods, which continue to provide first approximations to chemical structures,18,19 as well as a range of ab initio wave function treatments of electron correlation in ever-larger molecules.20−23 Representative examples of these are so-called divide-and-conquer,24,25 Nrepresentable reduced density matrix,26,27 density matrix renormalization group,28 QM/MM,29 fragment-based,30,31 and embedding32,33 approaches, as well as promising classical methods that can incorporate quantum corrections,34 and many other methods that aspired to both accuracy and linearscaling O(N) performance.35 It is commonly asserted that adoption of well-developed quantum-mechanical small-molecule methods for biomolecular calculations is computationally prohibitive in view of the large number of combined electronic and atomic degrees of freedom presented by even small proteins. Since most atoms in proteins contain less than 10 electrons, it would seem, however, that an increase in the number of degrees of freedom by less than an order of magnitude in combined ab initio electronic and molecular-dynamics calculations is most likely not, by itself, a prohibitive factor. Of course, combinations of some form of accurate ab initio methods with molecular-dynamics or Monte Carlo simulations applicable to large molecules would be highly desirable. Conventional wave-function-based quantum treatments of the potential energy surfaces of small molecules have in common a reliance on atomic and/or molecular orbitals in forming totally antisymmetric many-electron configurational state functions to provide support appropriate for the desired Schrödinger eigenstates of an individual large molecule or cluster. Such approaches entail repeated evaluations, transformations, and storage of large numbers of orbital-based oneand two-electron integrals at different geometries, and related Hamiltonian matrix assembly and diagonalization, requiring significant computational expenditure. Moreover, such calculations must generally be started all over again when a new molecule is addressed, with previous calculations generally B
DOI: 10.1021/acs.jpcb.6b02021 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
corresponding energy eigenvalues collected in the diagonal matrix E(R). Sets of electrons r ≡ (1, 2, ..., n) have been arbitrarily assigned to particular atoms (α = 1,2, ..., N), with R = (R1, R2, ···RN) referring collectively to all N atomic-position coordinate vectors (Rα). The row vector Φ(α)(i) includes the discrete and continuum orthonormal antisymmetric eigenstates of constituent atom α,57 specified by the quantum numbers (E, L, ML, S, MS, P)α,58 with i referring to the vector coordinates of the nα electrons measured relative to a coordinate system located at the atomic center Rα having space-fixed laboratoryframe orientation. The universal eigenstates associated with both the discrete and essential portions of the atomic eigenspectra can be determined variationally in appropriate L2 basis-set representations,59 and need not be point-wise exact eigensolutions of the corresponding nonrelativistic atomic Hamiltonians to ensure the validity of the ensuing theoretical development. The subscript O indicates adoption of an odometer ordering convention for the sequence of N-term atomic spectral products in the row vector Φ(r: R), in which the enumeration of the atomic eigenstates appearing later in eq 2 is advanced prior to those appearing earlier. This convention has a number of consequences, including particularly the specific forms of matrix representatives of Hamiltonian operators in the spectralproduct basis, with additional aspects of notational conventions described in the sequel and reported in great detail elsewhere.36−48 The basis of eq 2 is known to span only once the totally antisymmetric irreducible representation, or Pauli irrep, of the aggregate electron symmetric group (Snt),42 in spite of transforming under the tensor-product subgroup Sn1 ⊗ Sn2 ⊗···SnN of Snt,39 and to also span selected nontotally antisymmetric or non-Pauli irreps of Snt.60−64 As a consequence, the molecular Schrödinger eigenstates Ψ(r: R) constructed in the representation of eq 2 can include both unphysical nonPauli and physically significant Pauli solutions.42,65−67 When the Hilbert space of eq 2 provides a suitable closure limit56,68
information required in a universally applicable implementation. Additionally, aggregate energy determinations employing previously reported matrix partitioning and perturbation methods that avoid explicit diagonalizations of large Hamiltonian matrices are described. In section 4, calculations are reported in illustration of the components devised and assembled to implement the methodology. Convergence of the atomic-product representation in the absence of explicit electronic charge-transfer, effective compression of large-scale configuration-interaction calculations in the atomic-based representation, and equivalence of results obtained with standard quantum-mechanical calculations are illustrated, and aspects of scaling issues described. Additionally, calculations of the promotion energies of individual atoms and their pairwise interaction energies within a molecule are reported, apparently for the first time on the basis of predictions made directly from the molecular Schrödinger equations in the absence of additional subjective constraints.54 Calculations of large numbers of electronic states of polyatomic molecules are compared with conventional valence-bond results to demonstrate the accuracy and convergence of the method, and aspects of on-the-fly construction of potential energy surfaces in combination with Monte Carlo and moleculardynamics calculations are described in anticipation of applications to larger systems to be reported subsequently elsewhere. Discussion and concluding remarks are provided in section 5, and detailed connections are made with the aforementioned early combined studies of van der Waals and chemical interactions, and with semiempirical atomic- and diatomicrepresentation-based approximations to molecular potential energy surfaces.49−52,55
2. ATOMIC SPECTRAL-PRODUCT FORMULATION The spectral-product methodology is described in section 2A in a form particularly appropriate to the atomic pairwise development considered here, and notational conventions are established and the pairwise aggregate Hamiltonian matrix in the spectral-product basis reported. Unitary transformation procedures for constructing the correctly antisymmetric exactpair form of the aggregate Hamiltonian matrix are presented in section 2B, and computational aspects of the approach in finite subspace representations are described in section 2C; conditions are described in section 2D under which such representations can be employed in forms which do not require overall aggregate antisymmetry for molecular energy calculations. 2A. Spectral-Product Representation. Solutions of the Schrödinger equation56 Ĥ (r: R)Ψ(r: R) = Ψ(r: R) ·E(R)
Ĥ (r: R)Φ(r: R) → Φ(r: R) ·H(R)
a faithful Hermitian representative of the aggregate Hamiltonian matrix H(R) is obtained, with the basis supporting both Pauli and non-Pauli Schrödinger eigenstates and energies in this limit. Attempts to partition molecular Hamiltonian operators into meaningful “fragment” component operators encounter issues related to electron indistinguishability.69 Specifically, assignments of electrons to particular nuclei in a molecule are thought to generally have no validity in molecular quantum theory,70,71 since such assignments are not invariant under arbitrary electron permutations. Coulombic interactions between electrons and nuclei, for example, cannot be assigned as either intra- or interatomic terms under arbitrary electron permutations, seemingly preventing a meaningful partitioning of molecular Hamiltonian operators into fragment components. Moreover, assignments of indistinguishable electrons to particular nuclei or fragment operators are seemingly made meaningless upon integrations over the totally symmetric probability distributions obtained from products of totally antisymmetric basis states, providing total energy expressions in the familiar form of integrals over a single one- or two-electron operator and a corresponding one- or two-electron reduced
(1)
for a molecule or other atomic aggregate can be constructed variationally employing a Hilbert space of outer or tensor products (⊗) of many-electron atomic eigenstates, providing an orthonormal spectral-product (or “van der Waals”) representation of the universal form37 Φ(r: R) = {Φ(1)(1) ⊗ Φ(2)(2) ⊗ . . . Φ(N )(n)}O
(3)
(2)
Here, Ĥ (r: R) is the nonrelativistic many-electron (nt) Hamiltonian operator for a molecule having N atoms,56 and Ψ(r: R) is a row vector of molecular Schrödinger eigenstates to be represented in the row-vector basis Φ(r: R), with C
DOI: 10.1021/acs.jpcb.6b02021 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B density matrix constructed from the electronic distribution,54 casting doubt on the early atomic- and pair-based approaches employing explicitly antisymmetric wave functions.50,51 When made in conjunction with the spectral-product representation of eq 2, however, in which particular sets of electrons are “permanently” assigned to specific nuclei, corresponding assignments of electrons to individual atomic operators, and to the interaction potentials between them, can be made, and Hermitian matrix representatives of such operators evaluated in the spectral-product representation. In this way, an atomic-pairwise form of the total molecular Hamiltonian operator can be written as a sum over essentially self-adjoint atomic-pair operators, as well as over pairs of atomic operators which correct for their overcounting in the sums over pairs, as N
Ĥ (r: R) ≡ N−1
=
∑ α=1
−
N−1
N
∑ Ĥ (α)(i) + ∑ ∑ α=1 N
V̂
constructing a partitioned matrix representative of the total Hamiltonian operator of eq 4. Employing the formally complete orthonormal spectralproduct basis and the Coulombic Hamiltonian operator of eq 4 in a linear variational solution of eq 1 gives the familiar matrix Schrödinger equation56 H(R) ·UH(R) = UH(R) ·E(R)
where the columns of the unitary matrix UH(R) provide the desired Schrödinger eigenfunctions in the form Ψ(r: R) = Φ(r: R) ·UH(R)
E(R) = UH(R)† ·H(R) ·UH(R)
(α , β)
(i , j: R αβ)
α=1 β=α+1
Φ(α , β)(i , j) ≡ {Φ(α)(i) ⊗ Φ(β)(j)}
⎫ (N − 2) ̂ (α) (β) {H (i) + Ĥ (j)}⎬ (N − 1) ⎭
(4)
(5)
N
=
where the atomic operator Ĥ (i)= ⎧
Ze ℏ 2 ∇i − α ⎬ + riα ⎭ ⎩ 2m
∑ ⎨− i
ZαZβe 2 R αβ
nα
nα
2
∑ ∑ i
i ′= i + 1
e rii ′
−
∑ i
(6)
nβ
−
∑ j
Zαe 2 + r jα
nα
nβ
∑∑ i
j
e2 rij
∑ ∑
V (α , β)(R αβ)
⎧ (α , β) ⎨H (R αβ) β=α+1 ⎩ N
∑ ∑ −
⎫ (N − 2) (α) {H + H(β)}⎬ (N − 1) ⎭
(12)
Here, V(α,β) (Rαβ) = H(α,β) (Rαβ) − H(α) − H(β) is the pair interaction-energy Hamiltonian matrix, whereas the individual atomic-pair energy matrices in the spectral-product basis are
(i; j: R αβ)=
riβ
+
N
α=1 β=α+1
α=1
(α , β)
Zβe 2
∑H N−1
=
is a sum over appropriate terms for the electrons (i) assigned to the atom α, with a similar expression employed for the atom β with i → j, and the pairwise interaction term V̂
N−1 (α)
α=1
(α)
nα − 1
(11)
⟨Φ(r: R)|Ĥ (r: R)|Φ(r: R)⟩
H(R) ≡
(α , β) (α) (β) Ĥ (i) + Ĥ (j) + V̂ (i ; j: R αβ)
α
