Adaptive-Partitioning Multilayer Dynamics Simulations: 2. Implementations of the Permuted and Interpolated Adaptive-Partitioning Gradients.

Recently, an adaptive-partitioning multilayer Q1/Q2/MM method was proposed, where Q1 and Q2 denote, respectively, two distinct quantum-mechanical levels of theory and MM, the molecular-mechanical force fields. Such a multilayer model resembles the ONIOM (our own N-layered integrated molecular orbital and molecular mechanics) model by Morokuma and co-workers, but it is distinguished by on-the-fly reclassifying atoms to be Q1, Q2, or MM in dynamics simulations. To smoothly blend the levels of descriptions of the atoms, buffer zones are introduced between adjacent layers, and the energy is smoothly interpolated. In particular, the Q1/Q2 interaction energy was expressed in two different formalisms: permuted and interpolated adaptive-partitioning (PAP and IAP), respectively. While the PAP energy is based on a weighted many-body expansion, the IAP energy is derived via alchemical quantum calculations with interpolated Fock and overlap matrices. In this article, we examine in-depth the irregularities in the IAP energy near the boundary between the buffer and Q2 zones, which were found prominent in some calculations. These irregularities are due to basis-set linear dependencies, which can be effectively suppressed using a cutoff for the weighted atomic orbital coefficients. Furthermore, we derived and implemented the gradients for both PAP and IAP. Test calculations on a series of water cluster models show perfectly smooth gradients in PAP, while a minor discontinuity occurs in IAP gradients at the buffer/Q2 boundary. The energy and gradient discontinuities in IAP become smaller when moving the buffer/Q2 boundary further away from the Q1 center and when increasing the size of the basis sets used. Overall, those discontinuities are controllable, and possible ways to further diminish them are discussed.

Quasi-atomic orbital analysis of halogen bonding interactions.

A quasi-atomic orbital analysis of the halogen bonded NH3⋯XF complexes (X = F, Cl, Br, and I) is performed to gain insight into the electronic properties associated with these σ-hole interactions. It is shown that significant sharing of electrons between the nitrogen lone pair of the ammonia molecule and the XF molecule occurs, resulting in a weakening of the X-F bond. In addition, the N-X bond shows increasing covalent character as the size of the halogen atom X increases. While the Mulliken outer complex NH3⋯XF appears to be overall the main species, the strength of the covalent interaction of the N-X bond becomes increasingly similar to that of the N-X bond in the [NH3X]+ cation as the size of X increases.

The General Atomic and Molecular Electronic Structure System (GAMESS): Novel Methods on Novel Architectures.

The primary focus of GAMESS over the last 5 years has been the development of new high-performance codes that are able to take effective and efficient advantage of the most advanced computer architectures, both CPU and accelerators. These efforts include employing density fitting and fragmentation methods to reduce the high scaling of well-correlated (e.g., coupled-cluster) methods as well as developing novel codes that can take optimal advantage of graphical processing units and other modern accelerators. Because accurate wave functions can be very complex, an important new functionality in GAMESS is the quasi-atomic orbital analysis, an unbiased approach to the understanding of covalent bonds embedded in the wave function. Best practices for the maintenance and distribution of GAMESS are also discussed.

Adaptive-Partitioning Multilayer Dynamics Simulations: 1. On-the-Fly Switch between Two Quantum Levels of Theory.

We propose to generalize the previously developed two-layer permuted adaptive-partitioning quantum-mechanics/molecular-mechanics (QM/MM), which reclassifies atoms as QM or MM on-the-fly in dynamics simulations, to multilayer adaptive-partitioning algorithms that enable multiple levels of theory. In this work, we formulate two new algorithms that smoothly interpolate the energy between two QM (Q1 and Q2) levels of theory. The first "permuted adaptive-partitioning" scheme is based on the weighted many-body expansion of the potential, as in the adaptive-partitioning QM/MM. Unconventional and potentially more efficient, the second "interpolated adaptive-partitioning" method employs alchemical QM calculations with Q1/Q2-mixed basis sets, Fock matrices, and overlap matrices. To our knowledge, this is the first time that such alchemical calculations are performed in QM, although they are routinely done in MM. Test calculations on water-cluster models show that both new algorithms indeed yield smooth energy curves when water molecules shift between Q1 and Q2.

Computation of host-guest binding free energies with a new quantum mechanics based mining minima algorithm.

A new method called QM-VM2 is presented that efficiently combines statistical mechanics with quantum mechanical (QM) energy potentials in order to calculate noncovalent binding free energies of host-guest systems. QM-VM2 efficiently couples the use of semi-empirical QM (SEQM) energies and geometry optimizations with an underlying molecular mechanics (MM) based conformational search, to find low SEQM energy minima, and allows for processing of these minima at higher levels of ab initio QM theory. A progressive geometry optimization scheme is introduced as a means to increase conformational sampling efficiency. The newly implemented QM-VM2 is used to compute the binding free energies of the host molecule cucurbit[7]uril and a set of 15 guest molecules. The results are presented along with comparisons to experimentally determined binding affinities. For the full set of 15 host-guest complexes, which have a range of formal charges from +1 to +3, SEQM-VM2 based binding free energies show poor correlation with experiment, whereas for the ten +1 complexes only, a significant correlation (R2 = 0.8) is achieved. SEQM-VM2 generation of conformers followed by single-point ab initio QM calculations at the dispersion corrected restricted Hartree-Fock-D3(BJ) and TPSS-D3(BJ) levels of theory, as post-processing corrections, yields a reasonable correlation with experiment for the full set of host-guest complexes (R2 = 0.6 and R2 = 0.7, respectively) and an excellent correlation for the +1 formal charge set (R2 = 1.0 and R2 = 0.9, respectively), as long as a sufficiently large basis set (triple-zeta quality) is employed. The importance of the inclusion of configurational entropy, even at the MM level, for the achievement of good correlation with experiment was demonstrated by comparing the calculated ΔE values with experiment and finding a considerably poorer correlation with experiment than for the calculated free energy ΔE - TΔS. For the complete set of host-guest systems with the range of formal charges, it was observed that the deviation of the predicted binding free energy from experiment correlates somewhat with the net charge of the systems. This observation leads to a simple empirical interpolation scheme to improve the linear regression of the full set.

Why is Si2H2 Not Linear? An Intrinsic Quasi-Atomic Bonding Analysis.

The molecular energy of Si2H2 geometric structures increases in the order dibridged < trans-bent < linear, in contrast to acetylene, C2H2, for which the linear structure is the global minimum. In this study, the intra-atomic (antibonding) and bonding contributions to the total molecular energy of these valence isoelectronic molecules are computed by expressing the density matrices of the full valence space multiconfiguration self-consistent field wave function in terms of quasi-atomic orbitals. The analysis shows that the intra-atomic contributions to the molecular energy become less favorable in the order dibridged → trans-bent → linear for both C2H2 and Si2H2. By contrast, the inter-atomic bonding contributions become energetically more favorable in that order for both C2H2 and Si2H2. The two systems differ as follows. For Si2H2, the antibonding intra-atomic energy changes that occur when the dibridged molecule reconstructs into the trans-bent and linear structures prevail over the interatomic interactions that induce bond formation. In contrast, for C2H2, the interatomic interactions that create bonds prevail over the intra-atomic energy changes that occur when the dibridged molecule reconstructs into the trans-bent and linear structures. The intra-atomic energy changes that occur in these systems are related to the hybridization of the heavy atoms in an analogous manner to the hybridization of C in CH4 from (2s)2(2p)2 to sp3 hybrid orbitals.

Recent developments in the general atomic and molecular electronic structure system.

A discussion of many of the recently implemented features of GAMESS (General Atomic and Molecular Electronic Structure System) and LibCChem (the C++ CPU/GPU library associated with GAMESS) is presented. These features include fragmentation methods such as the fragment molecular orbital, effective fragment potential and effective fragment molecular orbital methods, hybrid MPI/OpenMP approaches to Hartree-Fock, and resolution of the identity second order perturbation theory. Many new coupled cluster theory methods have been implemented in GAMESS, as have multiple levels of density functional/tight binding theory. The role of accelerators, especially graphical processing units, is discussed in the context of the new features of LibCChem, as it is the associated problem of power consumption as the power of computers increases dramatically. The process by which a complex program suite such as GAMESS is maintained and developed is considered. Future developments are briefly summarized.

Quantum Mechanical Modeling of the Interactions between Noble Metal (Ag and Au) Nanoclusters and Water with the Effective Fragment Potential Method.

Explicit solvent interactions can significantly alter the physical and chemical properties of noble metal (e.g., gold and silver) nanoclusters. In order to compute these solvent interactions at a reasonable computational cost, a quantum mechanical (QM)/molecular mechanics (MM) approach, where the metal nanocluster is treated with full QM and the water molecules are treated with a MM force field, can be used. However, classical MM force fields were typically parameterized using molecules containing main group elements as the reference. The accuracy of noble metal-solvent interactions obtained with these force fields therefore remains unpredictable. The effective fragment potential (EFP) force field, designed to model explicitly solvated systems, represents an attractive method to simulate solvated noble metal nanoclusters because it is derived from first principles and contains few or no fitted parameters, depending on implementation. At the density functional theory-optimized geometries, good correlation is obtained between the nanocluster-water interaction energies computed with EFP and those computed with the reference coupled cluster singles, doubles, and perturbative triples method. It is shown that the EFP method gives qualitatively accurate interaction energies at medium-large intermolecular distances for various molecular configurations. In order to achieve higher quantitative accuracy, the first solvation shell should be treated with full QM, if possible. EFP is therefore a promising method for the QM modeling of explicitly solvated silver and gold nanoclusters.

Accuracy of the PM6 and PM7 Methods on Bare and Thiolate-Protected Gold Nanoclusters.

Semiempirical quantum mechanical (SEQM) methods offer an attractive middle ground between fully ab initio quantum chemistry and force-field simulations, allowing for a quantum mechanical treatment of the system at a relatively low computational cost. However, SEQM methods have not been frequently utilized in the study of transition metal systems, mostly due to the difficulty in obtaining reliable parameters. This paper examines the accuracy of the PM6 and PM7 semiempirical methods to predict geometries, ionization potentials, and HOMO-LUMO energy gaps of several bare gold clusters (Aun) and thiolate-protected gold nanoclusters (AuSNCs). Contrary to PM6, the PM7 method can predict qualitatively correct geometries and ionization potentials when compared to DFT. PM6 fails to predict the characteristic gold core and gold-sulfur ligand shell (staple motifs) of the AuSNC structures. Both the PM6 and PM7 methods overestimate the HOMO-LUMO gaps. Overall, PM7 provides a more accurate description of bare gold and gold-thiolate nanoclusters than PM6. Nevertheless, refining the gold parameters could help achieve better quantitative accuracy.

Benchmarking of the R-7 Anisotropic Dispersion Energy Term on the S22 Dimer Test Set.

The effects of including the anisotropic E7 term to the dispersion energy in addition to the leading E6 term are examined by using the effective fragment potential (EFP) method on the S22 test set. In this study, the full anisotropic E7 term is computed whereas the isotropic and spherical approximations are used for the E6 term. It is found that the E7 term is positive for hydrogen-bonded complexes and has a magnitude that can be as large as 50% of E6, giving rise to larger intermolecular distances than those obtained with E6 alone. The large positive value of E7 is analyzed for the hydrogen-bonded uracil dimer; it is found to originate from the large magnitude of the dynamic polarizability tensors as well as the proximity of the LMOs involved in hydrogen bonding. Conversely, E7 tends to be negative for dispersion-dominant complexes, and it has a very small magnitude for such complexes. The optimized geometries for these systems are therefore not greatly affected by the presence of the E7 term. For the mixed systems in the S22 test set, an intermediate behavior is observed. Overall, the E7 term is most important for systems with hydrogen bonding interactions and mixed systems. A full anisotropic treatment of the E6 term and higher order terms may need to be included to obtain more accurate interaction energies and geometries.

Benchmarking the Effective Fragment Potential Dispersion Correction on the S22 Test Set.

The usual modeling of dispersion interactions in density functional theory (DFT) is often limited by the use of empirically fitted parameters. In this study, the accuracies of the popular empirical dispersion corrections and the first-principles derived effective fragment potential (EFP) dispersion correction are compared by computing the DFT-D and HF-D equilibria interaction energies and intermolecular distances of the S22 test set dimers. Functionals based on the local density approximation (LDA) and generalized gradient approximation (GGA), as well as hybrid functionals, are compared for the DFT-D calculations using coupled cluster CCSD(T) at the complete basis set (CBS) limit as the reference method. In general, the HF-D(EFP) method provides accurate equilibrium dimerization energies and intermolecular distances for hydrogen-bonded systems compared to the CCSD(T)/CBS reference data without using any empirical parameters. For dispersion-dominant and mixed systems, the structures and interaction energies obtained with the B3LYP-D(EFP) method are similar to or better than those obtained with the other DFT-D and HF-D methods. Overall, the first-principles derived -D(EFP) correction presents a robust alternative to the empirical -D corrections when used with the B3LYP functional for dispersion-dominant and mixed systems or with Hartree-Fock for hydrogen-bonded systems.

Dispersion Interactions in Water Clusters.

The importance of dispersion forces in water clusters is examined using the effective fragment potential (EFP) method. Since the original EFP1 water potential does not include dispersion, a dispersion correction to the EFP1 potential (EFP1-D) was derived and implemented. The addition of dispersion to the EFP1 potential yields improved geometries for water clusters that contain 2-6 molecules. The importance of the odd E7 contribution to the dispersion energy is investigated. The E7 dispersion term is repulsive for all of the water clusters studied here and can have a magnitude that is as large as half of the E6 value. The E7 term therefore contributes to larger intermolecular distances for the optimized geometries. Inclusion of many-body effects and/or higher order terms may be necessary to further improve dispersion energies and optimized geometries.

Derivation and Implementation of the Gradient of the R(-7) Dispersion Interaction in the Effective Fragment Potential Method.

The dispersion interaction energy may be expressed as a sum over R(-n) terms, with n ≥ 6. Most implementations of the dispersion interaction in model potentials are terminated at n = 6. Those implementations that do include higher order contributions commonly only include even power terms, despite the fact that odd power terms can be important. Because the effective fragment potential (EFP) method contains no empirically fitted parameters, the EFP method provides a useful vehicle for examining the importance of the leading R(-7) odd power term in the dispersion expansion. To fully evaluate the importance of the R(-7) contribution to the dispersion energy, it is important to have analytic energy first derivatives for all terms. In the present work, the gradients of the term E7 ∼ R(-7) are derived analytically, implemented in the GAMESS software package, and evaluated relative to other terms in the dispersion expansion and relative to the total EFP interaction energy. Periodic boundary conditions in the minimum image convention are also implemented. A more accurate dispersion energy contribution can now be obtained during molecular dynamics simulations.

Time-dependent density functional theory study of the luminescence properties of gold phosphine thiolate complexes.

The origin of the emission of the gold phosphine thiolate complex (TPA)AuSCH(CH3)2 (TPA = 1,3,5-triaza-7-phosphaadamantanetriylphosphine) is investigated using time-dependent density functional theory (TDDFT). This system absorbs light at 3.6 eV, which corresponds mostly to a ligand-to-metal transition with some interligand character. The P-Au-S angle decreases upon relaxation in the S1 and T1 states. Our calculations show that these two states are strongly spin-orbit coupled at the ground state geometry. Ligand effects on the optical properties of this complex are also discussed by looking at the simple AuP(CH3)3SCH3 complex. The excitation energies differ by several tenths of an electronvolt. Excited state optimizations show that the excited singlet and triplet of the (TPA)AuSCH(CH3)2 complex are bent. On the other hand, the Au-S bond breaks in the excited state for the simple complex, and TDDFT is no longer an adequate method. The excited state energy landscape of gold phosphine thiolate systems is very complex, with several state crossings. This study also shows that the formation of the [(TPA)AuSCH(CH3)2]2 dimer is favorable in the ground state. The inclusion of dispersion interactions in the calculations affects the optimized geometries of both ground and excited states. Upon excitation, the formation of a Au-Au bond occurs, which results in an increase in energy of the low energy excited states in comparison to the monomer. The experimentally observed emission of the (TPA)AuSCH(CH3)2 complex at 1.86 eV cannot be unambiguously assigned and may originate from several excited states.

Dispersion correction derived from first principles for density functional theory and Hartree-Fock theory.

The modeling of dispersion interactions in density functional theory (DFT) is commonly performed using an energy correction that involves empirically fitted parameters for all atom pairs of the system investigated. In this study, the first-principles-derived dispersion energy from the effective fragment potential (EFP) method is implemented for the density functional theory (DFT-D(EFP)) and Hartree-Fock (HF-D(EFP)) energies. Overall, DFT-D(EFP) performs similarly to the semiempirical DFT-D corrections for the test cases investigated in this work. HF-D(EFP) tends to underestimate binding energies and overestimate intermolecular equilibrium distances, relative to coupled cluster theory, most likely due to incomplete accounting for electron correlation. Overall, this first-principles dispersion correction yields results that are in good agreement with coupled-cluster calculations at a low computational cost.

Quantum mechanical origin of the plasmon: from molecular systems to nanoparticles.

The surface plasmon resonance (SPR) of noble metal nanoparticles is reviewed in terms of both classical and quantum mechanical approaches. The collective oscillation of the free electrons responsible for the plasmon is well described using classical electromagnetic theory for large systems (from about 10 to 100 nm). In cases where quantum effects are important, this theory fails and first principle approaches like time-dependent density functional theory (TDDFT) must be used. In this paper, we give an account of the current understanding of the quantum mechanical origin of plasmon resonances. We provide some insight into how the discrete absorption spectrum of small noble metal clusters evolves into a strong plasmon peak with increasing particle size. The collective character of the plasmon is described in terms of the constructive addition of single-particle excitations. As the system size increases, the number of single-particle excitations increases as well. A configuration interaction (CI) approach can be applied to describe the optical properties of particles of all shapes and sizes, providing a consistent definition of plasmon resonances. Finally, we expand our analysis to thiolate-protected nanoparticles and analyze the effects of ligands on the plasmon.

Quantum coherent plasmon in silver nanowires: a real-time TDDFT study.

A plasmon-like phenomenon, arising from coinciding resonant excitations of different electronic characteristics in 1D silver nanowires, has been proposed based on theoretical linear absorption spectra. Such a molecular plasmon holds the potential for anisotropic nanoplasmonic applications. However, its dynamical nature remains unexplored. In this work, quantum dynamics of longitudinal and transverse excitations in 1D silver nanowires are carried out within the real-time time-dependent density functional theory framework. The anisotropic electron dynamics confirm that the transverse transitions of different electronic characteristics are collective in nature and oscillate in-phase with respect to each other. Analysis of the time evolutions of participating one-electron wave functions suggests that the transverse transitions form a coherent wave packet that gives rise to a strong plasmon resonance at the molecular level.

Plasmon resonance analysis with configuration interaction.

Dipole plasmon resonances are described quantum mechanically using configuration interaction (CI). A fictitious system of three interacting configurations is considered, which yields three excited states. Excited states energies and oscillator strengths are derived from the eigenvalues and eigenvectors of the CI matrix, where the diagonal elements α(i) (i = 1, 2, 3) correspond to the interacting one-electron transition energies and the off-diagonal elements ß(ij) correspond to the coupling between these configurations. The plasmonic state is easily identified by its higher energy and much larger oscillator strength. This high oscillator strength is due to a constructive addition of the eigenvectors contributing to this state. The maximum oscillator strength enhancement of the plasmon peak is equal to the number of configurations in the CI matrix (three here), which occurs in the ideal case where all α elements are equal and all ß elements are equal. When the transitions involved in the CI have different energies (different α values), the oscillator strength of the plasmon peak becomes smaller and its energy is shifted in comparison to the ideal case. Increasing all the coupling values from 0 up to the point where the coupling is similar in magnitude to the difference in α values leads to a rapid rise of the plasmon peak oscillator strength whereas its energy slightly blue-shifts. A further increase of the coupling values does not affect the oscillator strength of the plasmon peak, which remains near its maximum enhancement value, but drastically affects its energy, which rapidly rises. The plasmonic behavior of noble metal model systems is successfully described using configuration interaction.

Theoretical analysis of the optical excitation spectra of silver and gold nanowires.

The excitation spectra of linear atomic chains of silver and gold with various sizes have been calculated using time-dependent density functional theory. Silver chains show longitudinal and transverse peaks as well as a low-intensity d-band. The longitudinal peak, corresponding to the HOMO-LUMO transition (along the main axis of the chain), shifts linearly to the red as the length of the system increases, consistent with the particle-in-a-box model. The transverse peak remains at approximately constant energy for all systems studied and corresponds to ∑(m)→Π(m) transitions in the xy plane perpendicular to the chain. As the chain grows, transitions arising from d orbitals contribute to the transverse peak, which affects its oscillator strength. Contrary to silver, gold chains display a strong d-band that converges to a distinct pattern at a chain length of about twelve atoms. The transitions involved in the d-band originate from localized d-orbitals with a d(z(2)) character since they have the right symmetry to give transitions into the LUMO, LUMO + 1, …, which have ∑ symmetry. Transitions arising from these localized d-orbitals also affect the position of the longitudinal peak and generate a wide transverse band. Although the majority of the transitions involved in the transverse band have a d∑â†’Π or dΠ→∑ character, they are hidden by much stronger excitations of dΠâ†’Π character in gold nanowires.

Development of a charge-perturbed particle-in-a-sphere model for nanoparticle electronic structure.

The complex surface structure of gold-thiolate nanoparticles is known to affect the calculated density functional theory (DFT) excitation spectra. However, as the nanoparticle size increases, it becomes impractical to calculate the excitation spectrum using DFT. In this study, a new method is developed to determine the energy levels of the thiolate-protected gold nanoparticles [Au(25)(SR)(18)](-), Au(102)(SR)(44) and Au(144)(SR)(60). A 3 nm thiolate-protected nanoparticle is also modeled. The particle-in-a-sphere model is used to represent the core while the ligands are treated as point charge perturbations. The electronic structures obtained with this model are qualitatively similar to DFT results. The symmetry of the arrangement of the perturbations around the core plays a major role in determining the splitting of the orbitals. The radius chosen to represent the core also affects the orbital splitting. Increasing the number of perturbations around the core shifts the orbitals to higher energies but does not significantly change the band gaps and orbital splitting as long as the symmetrical arrangement of the perturbations is conserved. This model can be applied to any gold nanoparticle with a spherical core, regardless of its size or the nature of the ligands, at very low computational cost.