Comparison of Methods for Active Orbital Selection in Multiconfigurational Calculations.

Several methods of constructing the active orbital space for multiconfigurational wave functions are compared on typical moderately strongly or strongly correlated ground-state molecules. The relative merits of these methods and problems inherent in multiconfigurational calculations are discussed. Strong correlation in the ground electronic state is found typically in larger conjugated and in antiaromatic systems, transition states which involve bond breaking or formation, and transition metal complexes. Our examples include polyenes, polyacenes, the reactant, product and transition state of the Bergman cyclization, and two transition metal complexes: Hieber's anion [(CO)3FeNO]- and ferrocene. For the systems investigated, the simplest and oldest selection method, based on the fractional occupancy of unrestricted Hartree-Fock natural orbitals (the UNO criterion), yields the same active space as much more expensive approximate full CI methods. A disadvantage of this method used to be the difficulty of finding broken spin symmetry UHF solutions. However, our analytical method, accurate to fourth order in the orbital rotation angles (Tóth and Pulay J. Chem. Phys. 2016, 145, 164102.), has solved this problem. Two further advantages of the UNO criterion are that, unlike most other methods, it measures not only the energetic proximity to the Fermi level but also the magnitude of the exchange interaction with strongly occupied orbitals and therefore allows the estimation of the correlation strength for orbital selection in Restricted Active Space methods.

Analytical Energy Gradients for the Cluster-in-Molecule MP2 Method and Its Application to Geometry Optimizations of Large Systems.

An efficient analytical energy gradient algorithm for the cluster-in-molecule (CIM) second order Møller-Plesset perturbation theory (MP2) method is presented. In our algorithm, the gradient contributions from the nonseparable term of the two-body density matrix on a given atom is extracted from calculations on a cluster constructed for this atom. The other terms in the CIM-MP2 energy gradient expression are evaluated by constructing the density matrices of the whole system with the contributions from all clusters constructed. For basis sets with diffuse functions, tight CIM parameters are necessary to obtain accurate gradients. Benchmark calculations show that the CIM-MP2 method can accurately reproduce the conventional MP2 gradients and geometries for larger systems. The optimized structure of a 174-atom oligopeptide using the CIM-MP2 method with the cc-pVDZ basis set is in good agreement with the corresponding crystal structure. The present CIM-MP2 gradient program can be used for optimizing the geometries of large systems with hundreds of atoms on ordinary workstations.

Automatic Construction of the Initial Orbitals for Efficient Generalized Valence Bond Calculations of Large Systems.

We propose an efficient general strategy for generating initial orbitals for generalized valence bond (GVB) calculations which makes routine black-box GVB calculations on large systems feasible. Two schemes are proposed, depending on whether the restricted Hartree-Fock (RHF) wave function is stable (scheme I) or not (scheme II). In both schemes, the first step is the construction of active occupied orbitals and active virtual orbitals. In scheme I, active occupied orbitals are composed of the valence orbitals (the inner core orbitals are excluded), and the active virtual orbitals are obtained from the original virtual space by requiring its maximum overlap with the virtual orbital space of the same system at a minimal basis set. In scheme II, active occupied orbitals and active virtual orbitals are obtained from the set of unrestricted natural orbitals (UNOs), which are transformed from two sets of unrestricted HF spatial orbitals. In the next step, the active occupied orbitals and active virtual ones are separately transformed to localized orbitals. Localized occupied and virtual orbital pairs are formed using the Kuhn-Munkres (KM) algorithm and are used as the initial guess for the GVB orbitals. The optimized GVB wave function is obtained using the second-order self-consistent-field algorithm in the GAMESS program. With this procedure, GVB energies have been obtained for the lowest singlet and triplet states of polyacenes (up to decacene with 96 pairs) and the singlet ground state of two di-copper-oxygen-ammonia complexes. We have also calculated the singlet-triplet gaps for some polyacenes and the relative energy between two di-copper-oxygen-ammonia complexes with the block-correlated second-order perturbation theory based on the GVB reference.

Benchmark Relative Energies for Large Water Clusters with the Generalized Energy-Based Fragmentation Method.

The generalized energy-based fragmentation (GEBF) method has been applied to investigate relative energies of large water clusters (H2O)n (n = 32, 64) with the coupled-cluster singles and doubles with noniterative triple excitations (CCSD(T)) and second-order Møller-Plesset perturbation theory (MP2) at the complete basis set (CBS) limit. Here large water clusters are chosen to be representative structures sampled from molecular dynamics (MD) simulations of liquid water. Our calculations show that the GEBF method is capable of providing highly accurate relative energies for these water clusters in a cost-effective way. We demonstrate that the relative energies from GEBF-MP2/CBS are in excellent agreement with those from GEBF-CCSD(T)/CBS for these water clusters. With the GEBF-CCSD(T)/CBS relative energies as the benchmark results, we have assessed the performance of several theoretical methods widely used for ab initio MD simulations of liquids and aqueous solutions. These methods include density functional theory (DFT) with a number of different functionals, MP2, and density functional tight-binding (the third generation, DFTB3 in short). We find that MP2/aug-cc-pVDZ and several DFT methods (such as LC-ωPBE-D3 and ωB97XD) with the aug-cc-pVTZ basis set can provide satisfactory descriptions for these water clusters. Some widely used functionals (such as B3LYP, PBE0) and DFTB3 are not accurate enough for describing the relative energies of large water clusters. Although the basis set dependence of DFT is less than that of ab initio electron correlation methods, we recommend the combination of a few best functionals and large basis sets (at least aug-cc-pVTZ) in theoretical studies on water clusters or aqueous solutions.

Approximate Force Constants from Uncoupled Self-Consistent Field Perturbation Theory Using Nonhybrid Density Functional Theory.

Second derivatives of the molecular energy with respect to the nuclear coordinates (the nuclear Hessian or force constant matrix) are important for predicting infrared and Raman spectra, for calculating thermodynamic properties, for characterizing stationary states, and for guiding geometry optimization. However, their calculation for larger systems scales with molecular size one power higher than the calculation of the energy and the forces. The step responsible for the steep scaling of the nuclear Hessian is the coupled-perturbed self-consistent field (CP-SCF) iteration. This is omitted in the uncoupled SCF (UC-SCF) approximation. We have found that, though UC-SCF performs rather poorly at the Hartree-Fock and hybrid DFT levels, its performance for "pure" (non-hybrid) DFT is remarkably good. This is valid also for imaginary frequencies that characterize transition states. UC-SCF vibrational frequencies and normal modes are compared with coupled calculations for various exchange-correlation functionals including Hartree-Fock, and with basis sets ranging from simple to large for a variety of organic and some organometallic molecules. Their unexpectedly good performance makes them good candidates for calculating thermodynamic properties and for guiding difficult geometry optimizations, including the determination of transition states.

Finding symmetry breaking Hartree-Fock solutions: The case of triplet instability.

Determining the lowest unrestricted Hartree-Fock (UHF) solution is often difficult in even-electron systems. We have developed a deterministic method for locating approximately the UHF minimum using the restricted Hartree-Fock triplet instability matrix. The current method is truncated to fourth order. The minimum energy solution for this model can be determined by solving a small linear system of equations. This solution gives a suitable starting point to determine the exact UHF solution. This should be useful for the black-box determination of active spaces spanned by the fractionally occupied charge natural orbitals of the ground-state UHF wavefunction. The results can be generalized to higher (6th and 8th) degree expansions (odd expansion orders vanish by symmetry), and to other types of instability, including complex instability. The results are illustrated by calculations on ozone, benzene, nitrobenzene, butadiene, hexatriene, octatetraene, dichromium, and nickel porphine. Further examples are given in the supplementary material.

Selection of active spaces for multiconfigurational wavefunctions.

The efficient and accurate description of the electronic structure of strongly correlated systems is still a largely unsolved problem. The usual procedures start with a multiconfigurational (usually a Complete Active Space, CAS) wavefunction which accounts for static correlation and add dynamical correlation by perturbation theory, configuration interaction, or coupled cluster expansion. This procedure requires the correct selection of the active space. Intuitive methods are unreliable for complex systems. The inexpensive black-box unrestricted natural orbital (UNO) criterion postulates that the Unrestricted Hartree-Fock (UHF) charge natural orbitals with fractional occupancy (e.g., between 0.02 and 1.98) constitute the active space. UNOs generally approximate the CAS orbitals so well that the orbital optimization in CAS Self-Consistent Field (CASSCF) may be omitted, resulting in the inexpensive UNO-CAS method. A rigorous testing of the UNO criterion requires comparison with approximate full configuration interaction wavefunctions. This became feasible with the advent of Density Matrix Renormalization Group (DMRG) methods which can approximate highly correlated wavefunctions at affordable cost. We have compared active orbital occupancies in UNO-CAS and CASSCF calculations with DMRG in a number of strongly correlated molecules: compounds of electronegative atoms (F2, ozone, and NO2), polyenes, aromatic molecules (naphthalene, azulene, anthracene, and nitrobenzene), radicals (phenoxy and benzyl), diradicals (o-, m-, and p-benzyne), and transition metal compounds (nickel-acetylene and Cr2). The UNO criterion works well in these cases. Other symmetry breaking solutions, with the possible exception of spatial symmetry, do not appear to be essential to generate the correct active space. In the case of multiple UHF solutions, the natural orbitals of the average UHF density should be used. The problems of the UNO criterion and their potential solutions are discussed: finding the UHF solutions, discontinuities on potential energy surfaces, and inclusion of dynamical electron correlation and generalization to excited states.

Zundel-type H-bonding in biomolecular ions.

Using quantum chemical calculations and infrared multiphoton dissociation (IRMPD) spectroscopy in the fingerprint and X-H stretching regions, we demonstrate here that the all-Ala (b6) fragment ion features a macrocyclic structure with C(2) symmetry. For this structure, the ionizing proton is equally shared by the Ala(1) and Ala(4) amide oxygens in a Zundel-type symmetric (X…H(+)…X) H-bond.

The accuracy of quantum chemical methods for large noncovalent complexes.

We evaluate the performance of the most widely used wavefunction, density functional theory, and semiempirical methods for the description of noncovalent interactions in a set of larger, mostly dispersion-stabilized noncovalent complexes (the L7 data set). The methods tested include MP2, MP3, SCS-MP2, SCS(MI)-MP2, MP2.5, MP2.X, MP2C, DFT-D, DFT-D3 (B3-LYP-D3, B-LYP-D3, TPSS-D3, PW6B95-D3, M06-2X-D3) and M06-2X, and semiempirical methods augmented with dispersion and hydrogen bonding corrections: SCC-DFTB-D, PM6-D, PM6-DH2 and PM6-D3H4. The test complexes are the octadecane dimer, the guanine trimer, the circumcoronene…adenine dimer, the coronene dimer, the guanine-cytosine dimer, the circumcoronene…guanine-cytosine dimer, and an amyloid fragment trimer containing phenylalanine residues. The best performing method is MP2.5 with relative root mean square deviation (rRMSD) of 4 %. It can thus be recommended as an alternative to the CCSD(T)/CBS (alternatively QCISD(T)/CBS) benchmark for molecular systems which exceed current computational capacity. The second best non-DFT method is MP2C with rRMSD of 8 %. A method with the most favorable "accuracy/cost" ratio belongs to the DFT family: BLYP-D3, with an rRMSD of 8 %. Semiempirical methods deliver less accurate results (the rRMSD exceeds 25 %). Nevertheless, their absolute errors are close to some much more expensive methods such as M06-2X, MP2 or SCS(MI)-MP2, and thus their price/performance ratio is excellent.

A benchmark comparison of σ/σ and π/π dispersion: the dimers of naphthalene and decalin, and coronene and perhydrocoronene.

The stacking interaction between π systems is a well-recognized structural motif, but stacking between σ systems was long considered of secondary importance. A recent paper points out that σ stacking can reach the energy of chemical bonds and concludes that "σ/σ and π/π interactions are equally important" (Fokin, A. F.; Gerbig, D.; Schreiner, P. R. J. Am. Chem. Soc. 2011, 133, 20036). Our analysis shows that strong dispersion interaction requires rigid subsystems and good fits of their repulsive potential walls, conditions which are satisfied for both graphenes and larger graphanes (perhydrographenes). Comparison of the dimerization energies of decalin and perhydrocoronene with those of the naphthalene and coronene dimers at the coupled cluster (CC) CCSD(T) level confirms the substantial σ-stacking energies in graphanes but shows lower binding energies than do the B97D calculations of Fokin et al. Graphane dimerization energies are substantially lower at the CC level than the corresponding π-stacking energies: the value for perhydrocoronene is only 67% of the value for coronene, and the difference increases with system size. Our best estimate for the dimerization energy of naphthalene is 6.1 kcal/mol. Spin-component scaled MP2 is unbalanced: it gives only 70% of the CCSD(T) binding energy in σ dimers. The B3LYP-D3 method compares very well with CC for both σ and π dimers at the van der Waals minimum but underestimates the binding at larger distances. We used the largest possible atomic basis for these systems with 64-bit arithmetic, half-augmented-pVDZ, and the results were corrected for basis set incompleteness at the MP2 level.

The Ethidium-UA/AU Intercalation Site: Effect of Model Fragmentation and Backbone Charge State.

We report a systematic analysis of the intermolecular interactions of cationic ethidium intercalated into a UA/AU step of RNA for a single conformation based on crystallographic coordinates. Interaction energies at the MP2/6-31G** level were partitioned into electrostatic, exchange, delocalization, and correlation components. Various pairwise interaction models built from chemically intuitive fragments reproduce within a few percent values obtained when treating the intercalation site as a whole. Gas phase results are very sensitive to the charge state of the two phosphate groups, with the electrostatic term nearly tripling when the counterions are removed. But this is largely compensated by solvation, an effect represented here within the polarizable continuum model. In a few cases, more diffuse and larger basis sets as well as QCISD(T) corrections were applied in an effort to estimate plausible ethidium-nucleobase electron correlation effects.

A reliable and efficient first principles-based method for predicting pK(a) values. 2. Organic acids.

In part 1 of this series, we developed a protocol for the large-scale calculation of pK(a) values in aqueous solutions from first principles calculations, with the goal of striking a compromise between accuracy and computational efficiency. Following previous workers in the field, pK(a) values are calculated from a linear regression fit to deprotonation energies: pK(a)(f) = alpha(f)(E(A)(-) - E(HA)) + beta(f), where f denotes a family of functional groups. In this paper, we derive (alpha(f), beta(f)) values for the acidic functional groups -COOH, -POOH, alcoholic and phenolic -OH, -SH, -NHOH/ horizontal lineNOH, and -NROH, using a data set of 449 experimental pK(a) values. Several groupings of these functional groups were explored; our final recommended method uses five families (10 empirical parameters). Mean absolute deviations between our fits and experiment are 0.4 pK(a) units or less for each with a maximum error range of +/-1.5 pK(a) units. In certain subgroups, such as monocarboxylic acids, considerably better fits (mean absolute deviation approximately 0.20 pK(s) units) were obtained at the cost of more empirical parameters. Almost 70% of pK(a)'s calculated by our protocol lie within +/-0.4 pK(a) units and over 90% within +/-0.8 pK(a) units of the experimental reference value. Our results compare favorably with previous similar models which have greater computational cost.

A reliable and efficient first principles-based method for predicting pK(a) values. 1. Methodology.

We have developed an efficient and reliable protocol for the calculation of pK(a) values in aqueous solution from density functional calculations. We establish a standard linear regression fit using only calculated energies of deprotonation and experimental pK(a) values; all other factors, including most entropic effects, are absorbed into the fitting constants. In this article we fit a small training set of 34 experimentally well-characterized molecules to determine the best level of theory among those tested (i.e., the optimum compromise between efficiency and accuracy for the basis set, the exchange-correlation functional, the (continuum) solvation model and the level of geometry optimization). Our main findings are that a relatively modest basis set (6-311+G**) suffices for the calculation of the energy differences, with an even small basis set (3-21G*) sufficient for the preceding geometry optimization. Using a solvation model (COSMO in our case) throughout is essential to achieve reliable results. The exchange-correlation functional plays only a modest role; in particular, pure DFT functionals that allow the efficient calculation of the Coulomb term are perfectly adequate. The final protocol will be applied subsequently to data sets much larger than commonly used in such studies.

Comments on the molecular geometry of ferrocene: The dangers of using quantum chemistry programs as black boxes.

The dangers of using standard quantum chemistry programs as black boxes is illustrated by analyzing some results in a recent paper published in this journal (Zhang et al., J Comput Chem 2007, 28, 2260). The main danger is that nonlinear optimizations of both the wavefunction and the molecular geometry may converge to higher local minima or to saddle points, producing misleading results. For instance, some of the calculated molecular geometries of ferrocene in the aforementioned paper correspond to an SCF solution that converged to an excited state. This is the cause of the apparent large variation in the calculated iron-ring distance with the basis set. Another problem we noticed is that the source of the diffuse functions used in the earlier work in connection with the 6-31G and 6-311G basis sets for transition metals is not specified in the literature or the program manual. They are also a poor match for the 6-31G basis set. We re-emphasize that the 6-31G basis set used in this paper lacks the necessary diffuse d-type functions for the late first-row transition metals, and ought to be replaced by the m6-31G basis that offers a more balanced description of the atomic valence states.

Quantum chemistry in parallel with PQS.

This article describes the capabilities and performance of the latest release (version 4.0) of the Parallel Quantum Solutions (PQS) ab initio program package. The program was first released in 1998 and evolved from the TEXAS program package developed by Pulay and coworkers in the late 1970s. PQS was designed from the start to run on Linux-based clusters (which at the time were just becoming popular) with all major functionality being (a) fully parallel; and (b) capable of carrying out calculations on large-by ab initio standards-molecules, our initial aim being at least 100 atoms and 1000 basis functions with only modest memory requirements. With modern hardware and recent algorithmic developments, full accuracy, high-level calculations (DFT, MP2, CI, and Coupled-Cluster) can be performed on systems with up to several thousand basis functions on small (4-32 node) Linux clusters. We have also developed a graphical user interface with a model builder, job input preparation, parallel job submission, and post-job visualization and display.

Efficient Parallel Implementation of the CCSD External Exchange Operator and the Perturbative Triples (T) Energy Calculation.

A new, efficient parallel algorithm is presented for the most expensive step in coupled cluster singles and doubles (CCSD) energy calculations, the external exchange operator (EEO). The new implementation requires much less input/output than our previous algorithm and takes better advantage of integral screening. It is formulated as a series of matrix multiplications. Both the atomic orbital integrals and the corresponding CC coefficients are broken up into smaller blocks to diminish the memory requirement. Integrals are presorted to make their sparsity pattern more regular. This allows the simultaneous use of two normally conflicting techniques for speeding up the CCSD procedure: the use of highly efficient dense matrix multiplication routines and the efficient utilization of sparsity. We also describe an efficient parallel implementation of the perturbative triples correction to CCSD and related methods. Using the Array Files tool for distributed filesystems, parallelization is straightforward and does not compromise efficiency. Representative timings are shown for calculations with 282-1528 atomic orbitals, 68-228 correlated electrons, and various symmetries, C1 to C2h.

Parallel DFT gradients using the Fourier Transform Coulomb method.

The recently described Fourier Transform Coulomb (FTC) algorithm for fast and accurate calculation of Density Functional Theory (DFT) gradients (Füsti-Molnar, J Chem Phys 2003, 119, 11080) has been parallelized. We present several calculations showing the speed and accuracy of our new parallel FTC gradient code, comparing its performance with our standard DFT code. For that part of the total derivative Coulomb potential that can be evaluated in plane wave space, the current parallel FTC gradient algorithm is up to 200 times faster in total than our classical all-integral algorithm, depending on the system size and basis set, with essentially no loss in accuracy. Proposed modifications should further improve the overall performance relative to the classical algorithm.

New parallel algorithm for MP2 energy gradient calculations.

A new parallel algorithm has been developed for calculating the analytic energy derivatives of full accuracy second order Møller-Plesset perturbation theory (MP2). Its main projected application is the optimization of geometries of large molecules, in which noncovalent interactions play a significant role. The algorithm is based on the two-step MP2 energy calculation algorithm developed recently and implemented into the quantum chemistry program, GAMESS. Timings are presented for test calculations on taxol (C47H51NO14) with the 6-31G and 6-31G(d) basis sets (660 and 1032 basis functions, 328 correlated electrons) and luciferin (C11H8N2O3S2) with aug-cc-pVDZ and aug-cc-pVTZ (530 and 1198 basis functions, 92 correlated electrons). The taxol 6-31G(d) calculations are also performed with up to 80 CPU cores. The results demonstrate the high parallel efficiency of the program.

Accuracy of the three-body fragment molecular orbital method applied to Møller-Plesset perturbation theory.

The three-body energy expansion in the fragment molecular orbital method (FMO) was applied to the 2nd order Møller-Plesset theory (MP2). The accuracy of both the two and three-body expansions was determined for water clusters, alanine n-mers (alpha-helices and beta-strands) and one synthetic protein, using the 6-31G* and 6-311G* basis sets. At the best level of theory (three-body, two molecules/residues per fragment), the absolute errors in energy relative to ab initio MP2 were at most 1.2 and 5.0 mhartree, for the 6-31G* and 6-311G* basis sets, respectively. The relative accuracy was at worst 99.996% and 99.96%, for 6-31G* and 6-311G*, respectively. A three-body approximation was introduced and the optimum threshold value was determined. The protein calculation (6-31G*) at the production level (FMO2/2) took 3 h on 36 3.2-GHz Pentium 4 nodes and had the absolute error in the MP2 correlation energy of only 2 kcal/mol.