Search

Specialized computational chemistry packages have permanently reshaped the landscape of chemical and materials science by providing tools to support and guide experimental efforts and for the prediction of atomistic and electronic properties. In this regard, electronic structure packages have played a special role by using first-principle-driven methodologies to model complex chemical and materials processes. Over the past few decades, the rapid development of computing technologies and the tremendous increase in computational power have offered a unique chance to study complex transformations using sophisticated and predictive many-body techniques that describe correlated behavior of electrons in molecular and condensed phase systems at different levels of theory. In enabling these simulations, novel parallel algorithms have been able to take advantage of computational resources to address the polynomial scaling of electronic structure methods. In this paper, we briefly review the NWChem computational chemistry suite, including its history, design principles, parallel tools, current capabilities, outreach, and outlook.

Several closely related ab initio thermal mean-field theories for fermions, both well-established and new ones, are compared with one another at the formalism level and numerically. The theories considered are Fermi-Dirac theory; thermal Hartree-Fock (HF) theory; two modifications of the thermal single-determinant approximation of Kaplan and Argyres, Ann. Phys. 92, 1-24 (1975); and the first-order finite-temperature many-body perturbation theory based on a zero-temperature or thermal HF reference. Thermal full-configuration-interaction theory is used as the benchmark., (© 2024 Author(s). Published under an exclusive license by AIP Publishing.)

The widely used thermal Hartree-Fock (HF) theory is generalized to include the effect of electron correlation while maintaining its quasi-independent-particle framework. An electron-correlated internal energy (or grand potential) is postulated in consultation with the second-order finite-temperature many-body perturbation theory (MBPT), which then dictates the corresponding thermal orbital (quasiparticle) energies in such a way that all fundamental thermodynamic relations are obeyed. The associated density matrix is of a one-electron type, whose diagonal elements take the form of the Fermi-Dirac distribution functions, when the grand potential is minimized. The formulas for the entropy and chemical potential are unchanged from those of Fermi-Dirac or thermal HF theory. The theory thus stipulates a finite-temperature extension of the second-order Dyson self-energy of one-particle many-body Green's function theory and can be viewed as a second-order, diagonal, frequency-independent, thermal inverse Dyson equation. At low temperatures, the theory approaches finite-temperature MBPT of the same order, but it may outperform the latter at intermediate temperatures by including additional electron-correlation effects through orbital energies. A physical meaning of these thermal orbital energies is proposed (encompassing that of thermal HF orbital energies, which has been elusive) as a finite-temperature version of Janak's theorem., (© 2024 Author(s). Published under an exclusive license by AIP Publishing.)

A unified theory is presented for finite-temperature many-body perturbation expansions of the anharmonic vibrational contributions to thermodynamic functions, i.e., the free energy, internal energy, and entropy. The theory is diagrammatically size-consistent at any order, as ensured by the linked-diagram theorem proved in this study, and, thus, applicable to molecular gases and solids on an equal footing. It is also a basis-set-free formalism, just like its underlying Bose-Einstein theory, capable of summing anharmonic effects over an infinite number of states analytically. It is formulated by the Rayleigh-Schrödinger-style recursions, generating sum-over-states formulas for the perturbation series, which unambiguously converges at the finite-temperature vibrational full-configuration-interaction limits. Two strategies are introduced to reduce these sum-over-states formulas into compact sum-over-modes analytical formulas. One is a purely algebraic method that factorizes each many-mode thermal average into a product of one-mode thermal averages, which are then evaluated by the thermal Born-Huang rules. Canonical forms of these rules are proposed, dramatically expediting the reduction process. The other is finite-temperature normal-ordered second quantization, which is fully developed in this study, including a proof of thermal Wick's theorem and the derivation of a normal-ordered vibrational Hamiltonian at finite temperature. The latter naturally defines a finite-temperature extension of size-extensive vibrational self-consistent field theory. These reduced formulas can be represented graphically as Feynman diagrams with resolvent lines, which include anomalous and renormalization diagrams. Two order-by-order and one general-order algorithms of computing these perturbation corrections are implemented and applied up to the eighth order. The results show no signs of Kohn-Luttinger-type nonconvergence., (© 2023 Author(s). Published under an exclusive license by AIP Publishing.)

A second-order many-body perturbation correction to the relativistic Dirac-Hartree-Fock energy is evaluated stochastically by integrating 13-dimensional products of four-component spinors and Coulomb potentials. The integration in the real space of electron coordinates is carried out by the Monte Carlo (MC) method with the Metropolis sampling, whereas the MC integration in the imaginary-time domain is performed by the inverse-cumulative distribution function method. The computational cost to reach a given relative statistical error for spatially compact but heavy molecules is observed to be no worse than cubic and possibly quadratic with the number of electrons or basis functions. This is a vast improvement over the quintic scaling of the conventional, deterministic second-order many-body perturbation method. The algorithm is also easily and efficiently parallelized with 92% strong scalability going from 64 to 4096 processors.

A comprehensive and detailed account is presented for the finite-temperature many-body perturbation theory for electrons that expands in power series all thermodynamic functions on an equal footing. Algebraic recursions in the style of the Rayleigh-Schrödinger perturbation theory are derived for the grand potential, chemical potential, internal energy, and entropy in the grand canonical ensemble and for the Helmholtz energy, internal energy, and entropy in the canonical ensemble, leading to their sum-over-states analytical formulas at any arbitrary order. For the grand canonical ensemble, these sum-over-states formulas are systematically transformed to sum-over-orbitals reduced analytical formulas by the quantum-field-theoretical techniques of normal-ordered second quantization and Feynman diagrams extended to finite temperature. It is found that the perturbation corrections to energies entering the recursions have to be treated as a nondiagonal matrix, whose off-diagonal elements are generally nonzero within a subspace spanned by degenerate Slater determinants. They give rise to a unique set of linked diagrams-renormalization diagrams-whose resolvent lines are displaced upward, which are distinct from the well-known anomalous diagrams of which one or more resolvent lines are erased. A linked-diagram theorem is introduced that proves the size-consistency of the finite-temperature many-body perturbation theory at any order. General-order algorithms implementing the recursions establish the convergence of the perturbation series toward the finite-temperature full-configuration-interaction limit unless the series diverges. The normal-ordered Hamiltonian at finite temperature sheds light on the relationship between the finite-temperature Hartree-Fock and first-order many-body perturbation theories.

A scalable, stochastic algorithm evaluating the fourth-order many-body perturbation (MP4) correction to energy is proposed. Three hundred Goldstone diagrams representing the MP4 correction are computer generated and then converted into algebraic formulas expressed in terms of Green's functions in real space and imaginary time. They are evaluated by the direct (i.e., non-Markov, non-Metropolis) Monte Carlo (MC) integration accelerated by the redundant-walker and control-variate algorithms. The resulting MC-MP4 method is efficiently parallelized and is shown to display O(n 5.3 ) size-dependence of cost, which is nearly two ranks lower than the O(n 7 ) dependence of the deterministic MP4 algorithm. It evaluates the MP4/aug-cc-pVDZ energy for benzene, naphthalene, phenanthrene, and corannulene with the statistical uncertainty of 10 mE h (1.1% of the total basis-set correlation energy), 38 mE h (2.6%), 110 mE h (5.5%), and 280 mE h (9.0%), respectively, after about 10 9 MC steps.

In the Monte Carlo many-body perturbation (MC-MP) method, the conventional correlation-correction formula, which is a long sum of products of low-dimensional integrals, is first recast into a short sum of high-dimensional integrals over electron-pair and imaginary-time coordinates. These high-dimensional integrals are then evaluated by the Monte Carlo method with random coordinates generated by the Metropolis-Hasting algorithm according to a suitable distribution. The latter algorithm, while advantageous in its ability to sample nearly any distribution, introduces autocorrelation in sampled coordinates, which, in turn, increases the statistical uncertainty of the integrals and thus the computational cost. It also involves wasteful rejected moves and an initial "burn-in" step as well as displays hysteresis. Here, an algorithm is proposed that directly produces a random sequence of electron-pair coordinates for the same distribution used in the MC-MP method, which is free from autocorrelation, rejected moves, a burn-in step, or hysteresis. This direct-sampling algorithm is shown to accelerate second- and third-order Monte Carlo many-body perturbation calculations by up to 222% and 38%, respectively.

The use of many control variates is proposed as a method to accelerate the second- and third-order Monte Carlo (MC) many-body perturbation (MC-MP2 and MC-MP3) calculations. A control variate is an exactly integrable function that is strongly correlated or anti-correlated with the target function to be integrated by the MC method. Evaluating both integrals and their covariances in the same MC run, one can effect a mutual cancellation of the statistical uncertainties and biases in the MC integrations, thereby accelerating its convergence considerably. Six and thirty-six control variates, whose integrals are known a priori, are generated for MC-MP2 and MC-MP3, respectively, by systematically replacing one or more two-electron-integral vertices of certain configurations by zero-valued overlap-integral vertices in their Goldstone diagrams. The variances and covariances of these control variates are computed at a marginal cost, enhancing the overall efficiency of the MC-MP2 and MC-MP3 calculations by a factor of up to 14 and 20, respectively.

A finite-temperature many-body perturbation theory is presented, which expands in power series the electronic grand potential, chemical potential, internal energy, and entropy on an equal footing. Sum-over-states and sum-over-orbitals analytical formulas for the second-order perturbation corrections to these thermodynamic properties are obtained in a time-independent, nondiagrammatic, algebraic derivation, relying on the sum rules of the Hirschfelder-Certain degenerate perturbation energies in a degenerate subspace as well as nine algebraic identities for the zeroth-order thermal averages of one- through four-indexed quantities and products thereof. They reproduce numerically exactly the benchmark data obtained as the numerical derivatives of the thermal-full-configuration-interaction results for a wide range of temperatures.

A highly scalable stochastic algorithm is proposed and implemented for computing the basis-set-incompleteness correction to the diagonal, frequency-independent self-energy of the second-order many-body Green's function (GF2) theory within the explicitly correlated (F12) formalism. The 6-, 9-, 12-, and 15-dimensional integrals comprising the F12 correction are directly evaluated by the Monte Carlo method using appropriate weight functions for importance sampling. The method is naturally and easily parallelized, involves minimal memory space and no disk I/O, and can use virtually any mathematical form of a correlation factor. Its computational cost to correct all ionization energies (IEs) is observed to increase as the fourth power of system size, as opposed to the fifth power in the case of the deterministic counterparts. The GF2 calculations and their F12 corrections for the first IEs of C 60 and C 70 were executed on 128 graphical processing units (GF2) and 896 central processing units (F12), respectively, to reach the results with statistical errors of 0.04 eV or less. They showed that the basis-set-incompleteness (from aug-cc-pVDZ) accounts for only 50%-60% of the deviations from experiments, suggesting the significance of higher-order perturbation corrections.

A similarity-transformed equation-of-motion vibrational coupled-cluster (STEOM-XVCC) method is introduced as a one-mode theory with an effective vibrational Hamiltonian, which is similarity transformed twice so that its lower-order operators are dressed with higher-order anharmonic effects. The first transformation uses an exponential excitation operator, defining the equation-of-motion vibrational coupled-cluster (EOM-XVCC) method, and the second uses an exponential excitation-deexcitation operator. From diagonalization of this doubly similarity-transformed Hamiltonian in the small one-mode excitation space, the method simultaneously computes accurate anharmonic vibrational frequencies of all fundamentals, which have unique significance in vibrational analyses. We establish a diagrammatic method of deriving the working equations of STEOM-XVCC and prove their connectedness and thus size-consistency as well as the exact equality of its frequencies with the corresponding roots of EOM-XVCC. We furthermore elucidate the similarities and differences between electronic and vibrational STEOM methods and between STEOM-XVCC and vibrational many-body Green's function theory based on the Dyson equation, which is also an anharmonic one-mode theory. The latter comparison inspires three approximate STEOM-XVCC methods utilizing the common approximations made in the Dyson equation: the diagonal approximation, a perturbative expansion of the Dyson self-energy, and the frequency-independent approximation. The STEOM-XVCC method including up to the simultaneous four-mode excitation operator in a quartic force field and its three approximate variants are formulated and implemented in computer codes with the aid of computer algebra, and they are applied to small test cases with varied degrees of anharmonicity.

A thorough analytical and numerical characterization of the whole perturbation series of one-particle many-body Green's function (MBGF) theory is presented in a pedagogical manner. Three distinct but equivalent algebraic (first-quantized) recursive definitions of the perturbation series of the Green's function are derived, which can be combined with the well-known recursion for the self-energy. Six general-order algorithms of MBGF are developed, each implementing one of the three recursions, the ΔMPn method (where n is the perturbation order) [S. Hirata et al., J. Chem. Theory Comput. 11, 1595 (2015)], the automatic generation and interpretation of diagrams, or the numerical differentiation of the exact Green's function with a perturbation-scaled Hamiltonian. They all display the identical, nondivergent perturbation series except ΔMPn, which agrees with MBGF in the diagonal and frequency-independent approximations at 1≤n≤3 but converges at the full-configuration-interaction (FCI) limit at n=∞ (unless it diverges). Numerical data of the perturbation series are presented for Koopmans and non-Koopmans states to quantify the rate of convergence towards the FCI limit and the impact of the diagonal, frequency-independent, or ΔMPn approximation. The diagrammatic linkedness and thus size-consistency of the one-particle Green's function and self-energy are demonstrated at any perturbation order on the basis of the algebraic recursions in an entirely time-independent (frequency-domain) framework. The trimming of external lines in a one-particle Green's function to expose a self-energy diagram and the removal of reducible diagrams are also justified mathematically using the factorization theorem of Frantz and Mills. Equivalence of ΔMPn and MBGF in the diagonal and frequency-independent approximations at 1≤n≤3 is algebraically proven, also ascribing the differences at n = 4 to the so-called semi-reducible and linked-disconnected diagrams.

The explicitly correlated approach is one of the most important breakthroughs in ab initio electronic structure theory, providing arguably the most compact, accurate, and efficient ansatz for describing the correlated motion of electrons. Since Hylleraas first used an explicitly correlated wave function for the He atom in 1929, numerous attempts have been made to tackle the significant challenges involved in constructing practical explicitly correlated methods that are applicable to larger systems. These include identifying suitable mathematical forms of a correlated wave function and an efficient evaluation of many-electron integrals. R12 theory, which employs the resolution of the identity approximation, emerged in 1985, followed by the introduction of novel correlation factors and wave function ansätze, leading to the establishment of F12 theory in the 2000s. Rapid progress in recent years has significantly extended the application range of explicitly correlated theory, offering the potential of an accurate wave-function treatment of complex systems such as photosystems and semiconductors. This perspective surveys explicitly correlated electronic structure theory, with an emphasis on recent stochastic and deterministic approaches that hold significant promise for applications to large and complex systems including solids.

A stochastic algorithm is proposed and implemented that computes a basis-set-incompleteness (F12) correction to an ab initio second-order many-body perturbation energy as a short sum of 6- to 15-dimensional integrals of Gaussian-type orbitals, an explicit function of the electron-electron distance (geminal), and its associated excitation amplitudes held fixed at the values suggested by Ten-no. The integrals are directly evaluated (without a resolution-of-the-identity approximation or an auxiliary basis set) by the Metropolis Monte Carlo method. Applications of this method to 17 molecular correlation energies and 12 gas-phase reaction energies reveal that both the nonvariational and variational formulas for the correction give reliable correlation energies (98% or higher) and reaction energies (within 2 kJ mol -1 with a smaller statistical uncertainty) near the complete-basis-set limits by using just the aug-cc-pVDZ basis set. The nonvariational formula is found to be 2-10 times less expensive to evaluate than the variational one, though the latter yields energies that are bounded from below and is, therefore, slightly but systematically more accurate for energy differences. Being capable of using virtually any geminal form, the method confirms the best overall performance of the Slater-type geminal among 6 forms satisfying the same cusp conditions. Not having to precompute lower-dimensional integrals analytically, to store them on disk, or to transform them in a nonscalable dense-matrix-multiplication algorithm, the method scales favorably with both system size and computer size; the cost increases only as O(n 4 ) with the number of orbitals (n), and its parallel efficiency reaches 99.9% of the ideal case on going from 16 to 4096 computer processors.

Ice Ih displays several anomalous thermodynamic properties such as thermal contraction at low temperatures, an anomalous volume isotope effect (VIE) rendering the volume of D2O ice greater than that of H2O ice, and a pressure-induced transition to the high-density amorphous (HDA) phase. Furthermore, the anomalous VIE increases with temperature, despite its quantum-mechanical origin. Here, embedded-fragment ab initio second-order many-body perturbation (MP2) theory in the quasiharmonic approximation (QHA) is applied to the Gibbs energy of an infinite, proton-disordered crystal of ice Ih at wide ranges of temperatures and pressures. The quantum effect of nuclei moving in anharmonic potentials is taken into account from first principles without any empirical or nonsystematic approximation to either the electronic or vibrational Hamiltonian. MP2 predicts quantitatively correctly the thermal contraction at low temperatures, which is confirmed to originate from the volume-contracting hydrogen-bond bending modes (acoustic phonons). It qualitatively reproduces (but underestimates) the thermal expansion at higher temperatures, caused by the volume-expanding hydrogen-bond stretching (and to a lesser extent librational) modes. The anomalous VIE is found to be the result of subtle cancellations among closely competing isotope effects on volume from all modes. Consequently, even ab initio MP2 with the aug-cc-pVDZ and aug-cc-pVTZ basis sets has difficulty reproducing this anomaly, yielding qualitatively varied predictions of the sign of the VIE depending on such computational details as the choice of the embedding field. However, the temperature growth of the anomalous VIE is reproduced robustly and is ascribed to the librational modes. These solid-state MP2 calculations, as well as MP2 Born-Oppenheimer molecular dynamics, find a volume collapse and a loss of symmetry and long-range order in ice Ih upon pressure loading of 2.35 GPa or higher. Concomitantly, rapid softening of acoustic phonons is observed starting around 2 GPa. They constitute a computational detection of a mechanical instability in ice Ih and the resulting pressure-induced amorphization to HDA.

Diagrammatically size-consistent and basis-set-free vibrational coupled-cluster (XVCC) theory for both zero-point energies and transition frequencies of a molecule, the latter through the equation-of-motion (EOM) formalism, is defined for an nth-order Taylor-series potential energy surface (PES). Quantum-field-theoretical tools (the rules of normal-ordered second quantization and Feynman-Goldstone diagrams) for deriving their working equations are established. The equations of XVCC and EOM-XVCC including up to the mth-order excitation operators are derived and implemented with the aid of computer algebra in the range of 1 ≤ m ≤ 8. Algorithm optimizations known as strength reduction, intermediate reuse, and factorization are carried out before code generation, reducing the cost scaling of the mth-order XVCC and EOM-XVCC in an nth-order Taylor-series PES (m ≥ n) to the optimal value of O(N(m+⌊n/2⌋)), where N is the number of modes. The calculated zero-point energies and frequencies of fundamentals, overtones, and combinations as well as Fermi-resonant modes display rapid and nearly monotonic convergence with m towards the exact values for the PES. The theory with the same excitation rank as the truncation order of the Taylor-series PES (m = n) seems to strike the best cost-accuracy balance, achieving the accuracy of a few tenths of cm(-1) for transitions involving (m - 3) modes and of a few cm(-1) for those involving (m - 2) modes. The relationships between XVCC and the vibrational coupled-cluster theories of Prasad and coworkers and of Christiansen and coworkers as well as the size-extensive vibrational self-consistent-field and many-body perturbation theories are also elucidated.

A useful definition of orbital degeneracy—form-degeneracy—is introduced, which is distinct from the usual energy-degeneracy: Two canonical spatial orbitals are form-degenerate when the energy expectation value in the restricted Hartree-Fock (RHF) wave function is unaltered upon a two-electron excitation from one of these orbitals to the other. Form-degenerate orbitals tend to have isomorphic electron densities and occur in the highest-occupied and lowest-unoccupied molecular orbitals (HOMOs and LUMOs) of strongly correlated systems. Here, we present a mathematical proof of the existence of a triplet instability in a real or complex RHF wave function of a finite system in the space of real or complex unrestricted Hartree-Fock wave functions when HOMO and LUMO are energy- or form-degenerate. We also show that a singlet instability always exists in a real RHF wave function of a finite system in the space of complex RHF wave functions, when HOMO and LUMO are form-degenerate, but have nonidentical electron densities, or are energy-degenerate. These theorems provide Hartree-Fock-theory-based explanations of Hund's rule, a singlet instability in Jahn-Teller systems, biradicaloid electronic structures, and a triplet instability during some covalent bond breaking. They also suggest (but not guarantee) the spontaneous formation of a spin density wave (SDW) in a metallic solid. The stability theory underlying these theorems extended to a continuous orbital-energy spectrum proves the existence of an oscillating (nonspiral) SDW instability in one- and three-dimensional homogeneous electron gases, but only at low densities or for strong interactions.

One-dimensional (1D) solids exhibit a number of striking electronic structures including charge-density wave (CDW) and spin-density wave (SDW). Also, the Peierls theorem states that at zero temperature, a 1D system predicted by simple band theory to be a metal will spontaneously dimerize and open a finite fundamental bandgap, while at higher temperatures, it will assume the equidistant geometry with zero bandgap (a Peierls transition). We computationally study these unique electronic structures and transition in polyyne and all-trans polyacetylene using finite-temperature generalizations of ab initio spin-unrestricted Hartree-Fock (UHF) and spin-restricted coupled-cluster doubles (CCD) theories, extending upon previous work [He et al., J. Chem. Phys. 140, 024702 (2014)] that is based on spin-restricted Hartree-Fock (RHF) and second-order many-body perturbation (MP2) theories. Unlike RHF, UHF can predict SDW as well as CDW and metallic states, and unlike MP2, CCD does not diverge even if the underlying RHF reference wave function is metallic. UHF predicts a gapped SDW state with no dimerization at low temperatures, which gradually becomes metallic as the temperature is raised. CCD, meanwhile, confirms that electron correlation lowers the Peierls transition temperature. Furthermore, we show that the results from all theories for both polymers are subject to a unified interpretation in terms of the UHF solutions to the Hubbard-Peierls model using different values of the electron-electron interaction strength, U/t, in its Hamiltonian. The CCD wave function is shown to encompass the form of the exact solution of the Tomonaga-Luttinger model and is thus expected to describe accurately the electronic structure of Luttinger liquids.

A stochastic algorithm based on Metropolis Monte Carlo (MC) is presented for the size-extensive vibrational self-consistent field methods (XVSCF(n) and XVSCF[n]) for anharmonic molecular vibrations. The new MC-XVSCF methods substitute stochastic evaluations of a small number of high-dimensional integrals of functions of the potential energy surface (PES), which is sampled on demand, for diagrammatic equations involving high-order anharmonic force constants. This algorithm obviates the need to evaluate and store any high-dimensional partial derivatives of the potential and can be applied to the fully anharmonic PES without any Taylor-series approximation in an intrinsically parallelizable algorithm. The MC-XVSCF methods reproduce deterministic XVSCF calculations on the same Taylor-series PES in all energies, frequencies, and geometries. Calculations using the fully anharmonic PES evaluated on the fly with electronic structure methods report anharmonic effects on frequencies and geometries of much greater magnitude than deterministic XVSCF calculations, reflecting an underestimation of anharmonic effects in a Taylor-series approximation to the PES.

A normal-ordered second-quantized form of the Hamiltonian is derived for quantum dynamics in a bound potential energy surface expressed as a Taylor series in an arbitrary set of orthogonal, delocalized coordinates centered at an arbitrary geometry. The constant, first-, and second-order excitation amplitudes of this Hamiltonian are identified as the ground-state energy, gradients, and frequencies, respectively, of the size-extensive vibrational self-consistent field (XVSCF) method or the self-consistent phonon method. They display the well-defined size dependence of V(1-n/2), where V is the volume and n is the number of coordinates associated with the amplitudes. It is used to rapidly derive the equations of XVSCF and vibrational many-body perturbation methods with the Møller-Plesset partitioning of the Hamiltonian.

A new quantum Monte Carlo (QMC) method for anharmonic vibrational zero-point energies and transition frequencies is developed, which combines the diagrammatic vibrational many-body perturbation theory based on the Dyson equation with Monte Carlo integration. The infinite sums of the diagrammatic and thus size-consistent first- and second-order anharmonic corrections to the energy and self-energy are expressed as sums of a few m- or 2m-dimensional integrals of wave functions and a potential energy surface (PES) (m is the vibrational degrees of freedom). Each of these integrals is computed as the integrand (including the value of the PES) divided by the value of a judiciously chosen weight function evaluated on demand at geometries distributed randomly but according to the weight function via the Metropolis algorithm. In this way, the method completely avoids cumbersome evaluation and storage of high-order force constants necessary in the original formulation of the vibrational perturbation theory; it furthermore allows even higher-order force constants essentially up to an infinite order to be taken into account in a scalable, memory-efficient algorithm. The diagrammatic contributions to the frequency-dependent self-energies that are stochastically evaluated at discrete frequencies can be reliably interpolated, allowing the self-consistent solutions to the Dyson equation to be obtained. This method, therefore, can compute directly and stochastically the transition frequencies of fundamentals and overtones as well as their relative intensities as pole strengths, without fixed-node errors that plague some QMC. It is shown that, for an identical PES, the new method reproduces the correct deterministic values of the energies and frequencies within a few cm(-1) and pole strengths within a few thousandths. With the values of a PES evaluated on the fly at random geometries, the new method captures a noticeably greater proportion of anharmonic effects.

The structure, equation of state, IR, Raman, and inelastic neutron scattering (INS) spectra of high-pressure, proton-ordered phase VIII of ice are calculated by the second-order many-body perturbation and coupled-cluster singles and doubles methods. Nearly all the observed features of the pressure-dependence of the structures and spectra are reproduced computationally up to 60 GPa insofar as the anharmonic effects can be neglected. The calculations display no sign of the hypothetical isostructural transition in 2-3 GPa to phase VIII('), the existence of which has been a matter of controversy for over a decade, while they do not contradict the interpretation of the spectral anomaly at 10-14 GPa as a precursor of the VIII-X phase transition. The calculated INS spectra correct a systematic error in the peak positions of the observed spectra.

A stochastic algorithm is proposed that can compute the basis-set-incompleteness correction to the second-order many-body perturbation (MP2) energy of a polyatomic molecule. It evaluates the sum of two-, three-, and four-electron integrals over an explicit function of electron-electron distances by a Monte Carlo (MC) integration at an operation cost per MC step increasing only quadratically with size. The method can reproduce the corrections to the MP2/cc-pVTZ energies of H2O, CH4, and C6H6 within a few mEh after several million MC steps. It circumvents the resolution-of-the-identity approximation to the nonfactorable three-electron integrals usually necessary in the conventional explicitly correlated (R12 or F12) methods.

A new, alternative set of interpretation rules of Feynman-Goldstone diagrams for many-body perturbation theory is proposed, which translates diagrams into algebraic expressions suitable for direct Monte Carlo integrations. A vertex of a diagram is associated with a Coulomb interaction (rather than a two-electron integral) and an edge with the trace of a Green's function in real space and imaginary time. With these, 12 diagrams of third-order many-body perturbation (MP3) theory are converted into 20-dimensional integrals, which are then evaluated by a Monte Carlo method. It uses redundant walkers for convergence acceleration and a weight function for importance sampling in conjunction with the Metropolis algorithm. The resulting Monte Carlo MP3 method has low-rank polynomial size dependence of the operation cost, a negligible memory cost, and a naturally parallel computational kernel, while reproducing the correct correlation energies of small molecules within a few mEh after 10(6) Monte Carlo steps.

Finite-temperature extensions of ab initio Gaussian-basis-set spin-restricted Hartree-Fock (HF) and second-order many-body perturbation (MP2) theories are implemented for infinitely extended, periodic, one-dimensional solids and applied to the Peierls and charge-density-wave (CDW) transitions in polyyne and all-trans polyacetylene. The HF theory predicts insulating CDW ground states for both systems in their equidistant structures at low temperatures. In the same structures, they turn metallic at high temperatures. Starting from the "dimerized" low-temperature equilibrium structures, the systems need even higher temperatures to undergo a Peierls transition, which is accompanied by geometric as well as electronic distortions from dimerized to non-dimerized forms. The conventional finite-temperature MP2 theory shows a sign of divergence in any phase at any nonzero temperature and is useless. The renormalized finite-temperature MP2 (MP2R) theory is divergent only near metallic electronic structures, but is well behaved elsewhere. MP2R also predicts CDW and Peierls transitions occurring at two different temperatures. The effect of electron correlation is primarily to lower the Peierls transition temperature.

A simple, easily implemented, accurate, and efficient approximation of long-range electron-electron-repulsion and electron-nucleus-attraction integrals is proposed. It replaces each product of two atomic-orbital (AO) basis functions of an electron by a point charge centered at the midpoint of the two AO's. The magnitude of the point charge is equal to the overlap integral of the two AO's. Each integral is then rapidly evaluated in the direct algorithm as a Coulomb interaction between two point charges. This scheme is implemented in ab initio Hartree-Fock crystalline orbital theory and tested for one-, two-, and three-dimensional solids of metallic, semimetallic, and nonmetallic electronic structures, in which the lattice sums of the direct Coulomb and/or exchange interactions are expected to be slowly convergent. It is shown that this approximation reduces operation and/or memory costs by up to an order of magnitude to achieve converged lattice sums, although the scaling (size dependence) of operation cost is unchanged. An improved criterion for truncating the exchange lattice sum is also proposed.

Second-order many-body perturbation theories for anharmonic vibrational frequencies and zero-point energies of molecules are formulated, implemented, and tested. They solve the vibrational Dyson equation self-consistently by taking into account the frequency dependence of the Dyson self-energy in the diagonal approximation, which is expanded in a diagrammatic perturbation series up to second order. Three reference wave functions, all of which are diagrammatically size consistent, are considered: the harmonic approximation and diagrammatic vibrational self-consistent field (XVSCF) methods with and without the first-order Dyson geometry correction, i.e., XVSCF[n] and XVSCF(n), where n refers to the truncation rank of the Taylor-series potential energy surface. The corresponding second-order perturbation theories, XVH2(n), XVMP2[n], and XVMP2(n), are shown to be rigorously diagrammatically size consistent for both total energies and transition frequencies, yield accurate results (typically within a few cm(-1) at n = 4 for water and formaldehyde) for both quantities even in the presence of Fermi resonance, and have access to fundamentals, overtones, and combinations as well as their relative intensities as residues of the vibrational Green's functions. They are implemented into simple algorithms that require only force constants and frequencies of the reference methods (with no basis sets, quadrature, or matrix diagonalization at any stage of the calculation). The rules for enumerating and algebraically interpreting energy and self-energy diagrams are elucidated in detail.

Kohn and Luttinger [Phys. Rev. 118, 41 (1960)] showed that the conventional finite-temperature extension of the second-order many-body perturbation theory had the incorrect zero-temperature limit in metals and, on this basis, argued that the theory was incorrect. We show that this inconsistency arises from the noninclusion of the temperature effect in the energies of the zeroth-order eigenstates of the perturbation theory, which causes not only the Kohn-Luttinger conundrum but also another inconsistency with the zero-temperature many-body perturbation theory, namely, the different rates of divergence of the correlation energy in a homogeneous electron gas (HEG). We propose a renormalized many-body perturbation theory derivable from the finite-temperature extension of the normal-ordered second quantization applied to the denominators of the energy expression, which involves the energies of the zeroth-order states, as well as to the numerators. The renormalized theory is shown to have the correct zero-temperature limit and the same rate of divergence in a HEG as the zero-temperature counterpart, and is, therefore, the correct finite-temperature many-body perturbation theory.

A stochastic method is proposed that evaluates the second-order perturbation corrections to the Dyson self-energies of a molecule (i.e., quasiparticle energies or correlated ionization potentials and electron affinities) directly and not as small differences between two large, noisy quantities. With the aid of a Laplace transform, the usual sum-of-integral expressions of the second-order self-energy in many-body Green's function theory are rewritten into a sum of just four 13-dimensional integrals, 12-dimensional parts of which are evaluated by Monte Carlo integration. Efficient importance sampling is achieved with the Metropolis algorithm and a 12-dimensional weight function that is analytically integrable, is positive everywhere, and cancels all the singularities in the integrands exactly and analytically. The quasiparticle energies of small molecules have been reproduced within a few mEh of the correct values with 10(8) Monte Carlo steps. Linear-to-quadratic scaling of the size dependence of computational cost is demonstrated even for these small molecules.

The symmetric-stretching fundamental (ν1) and the bending first overtone (2ν2) of CO2, which are accidentally degenerate with the same symmetry, undergo a Fermi resonance and give rise to two Raman bands with a frequency difference of 107 cm(-1) and an intensity ratio of 2.1. Both the frequency difference and intensity ratio can be varied by pressure applied to CO2 in condensed phases, which has been utilized as a spectroscopic geobarometer for minerals with CO2 inclusion. This study calculates the pressure dependence of the Fermi dyad frequency difference and intensity ratio by combining the embedded-fragment second-order Mo̸ller-Plesset perturbation calculations of harmonic frequencies of solid CO2 under pressure and the coupled-cluster singles and doubles with noniterative triples and vibrational configuration-interaction calculations of anharmonic frequencies of molecular CO2. It reproduces frequency difference quantitatively and intensity ratio qualitatively up to 10 GPa. The analysis of the results is shown to render strong support for one particular order of unperturbed frequencies, ν1 > 2ν2, in both the gas and solid phases, which has been a matter of controversy for decades.

With the aid of the Laplace transform, the canonical expression of the second-order many-body perturbation correction to an electronic energy is converted into the sum of two 13-dimensional integrals, the 12-dimensional parts of which are evaluated by Monte Carlo integration. Weight functions are identified that are analytically normalizable, are finite and non-negative everywhere, and share the same singularities as the integrands. They thus generate appropriate distributions of four-electron walkers via the Metropolis algorithm, yielding correlation energies of small molecules within a few mE(h) of the correct values after 10(8) Monte Carlo steps. This algorithm does away with the integral transformation as the hotspot of the usual algorithms, has a far superior size dependence of cost, does not suffer from the sign problem of some quantum Monte Carlo methods, and potentially easily parallelizable and extensible to other more complex electron-correlation theories.

Ice Ih is arguably the most important molecular crystal in nature, yet our understanding of its structural and dynamical properties is still far from complete. We present embedded-fragment calculations of the structures and vibrational spectra of the three-dimensional, proton-disordered phase of ice Ih performed at the level of second-order many-body perturbation theory with a basis-set superposition error correction. Our calculations address previous controversies such as the one related to the O-H bond length as well as the existence of two types of hydrogen bonds with strengths differing by a factor of two. For the latter, our calculations suggest that the observed spectral features arise from the directionality or the anisotropy of collective hydrogen-bond stretching vibrations rather than the previously suggested vastly different force constants. We also report a capability to efficiently compute infrared and Raman intensities of a periodic solid. Our approach reproduces the infrared and Raman spectra, the variation of inelastic neutron scattering spectra with deuterium concentration, and the anomaly of heat capacities at low temperatures for ice Ih.

A procedure to determine optimal vibrational coordinates is developed on the basis of an earlier idea of Thompson and Truhlar [J. Chem. Phys. 77, 3031 (1982)]. For a given molecule, these coordinates are defined as the unitary transform of the normal coordinates that minimizes the energy of the vibrational self-consistent-field (VSCF) method for the ground state. They are justified by the fact that VSCF in these coordinates becomes exact in two limiting cases: harmonic oscillators, where the optimized coordinates are normal, and noninteracting anharmonic oscillators, in which the optimized coordinates are localized on individual oscillators. A robust and general optimization algorithm is developed, which decomposes the transformation matrix into a product of Jacobi matrices, determines the rotation angle of each Jacobi matrix that minimizes the energy, and iterates the process until a minimum in the whole high dimension is reached. It is shown that the optimized coordinates are neither entirely localized nor entirely delocalized (or normal) in any of the molecules (the water, water dimer, and ethylene molecules) examined (apart from the aforementioned limiting cases). Rather, high-frequency stretching modes tend to be localized, whereas low-frequency skeletal vibrations remain normal. On the basis of these coordinates, we introduce two new vibrational structure methods: optimized-coordinate VSCF (oc-VSCF) and optimized-coordinate vibrational configuration interaction (oc-VCI). For the modes that become localized, oc-VSCF is found to outperform VSCF, whereas, for both classes of modes, oc-VCI exhibits much more rapid convergence than VCI with respect to the rank of excitations. We propose a rational configuration selection for oc-VCI when the optimized coordinates are localized. The use of the optimized coordinates in VCI with this configuration selection scheme reduces the mean absolute errors in the frequencies of the fundamentals and the first overtones/combination tones from 104.7 (VCI) to 10.7 (oc-VCI) and from 132.4 (VCI) to 8.2 (oc-VCI) cm(-1) for the water molecule and the water dimer, respectively. It is also shown that the degree of coupling in the potential for ethylene is reduced effectively from four modes to three modes by the transformation from the normal to optimized coordinates, which enhances the accuracy of oc-VCI with low-rank excitations.

Can the zero-point vibrational energies (ZPVE) of molecular clusters and crystals be evaluated as sums of ZPVE of constituent molecular fragments embedded in the cluster or crystal electrostatic environment? What is the appropriate unit of fragmentation: monomers or overlapping dimers? Can the contributions of acoustic phonons, which are fundamentally delocalized, be recuperated at satisfactory accuracy? These questions are answered by this study applying embedded monomer- and dimer-fragmentation methods to the harmonic ZPVE of hydrogen fluoride clusters, hydrogen fluoride crystal, and water clusters. Our findings are as follows: (1) ZPVE are reproduced accurately by both fragmentation schemes within a few percents of exact values or a few tenths of 1 kcal mol(-1) per molecule even for crystalline hydrogen fluoride, which has acoustic phonons. (2) Both the monomer- and dimer-based fragmentation are nearly equally accurate and useful for the absolute values of ZPVE, but the latter is more reliable than the former in reproducing the relative ZPVE of cluster isomers of the same size. (3) The embedding field is essential as it renders nonzero frequencies to the translational and rotational motions of monomers and dimers, accounting for the pseudo-translational and librational motions of the entire clusters or crystals. (4) Some of these low-frequency modes of fragments are calculated to have imaginary frequencies because the fragments are not at their equilibrium geometries, causing ZPVE to be complex. The imaginary part of ZPVE, which is nonphysical and is guaranteed to vanish in the exact limit of the many-body expansion, is nonetheless a useful estimate of errors in the real part.

In the size-extensive vibrational self-consistent field (XVSCF) method introduced earlier [M. Keçeli and S. Hirata, J. Chem. Phys. 135, 134108 (2011)], only a small subset of even-order force constants that can form connected diagrams were used to compute extensive total energies and intensive transition frequencies. The mean-field potentials of XVSCF formed with these force constants have been shown to be effectively harmonic, making basis functions, quadrature, or matrix diagonalization in the conventional VSCF method unnecessary. We introduce two size-consistent VSCF methods, XVSCF(n) and XVSCF[n], for vibrationally averaged geometries in addition to energies and frequencies including anharmonic effects caused by up to the nth-order force constants. The methods are based on our observations that a small number of odd-order force constants of certain types can form open, connected diagrams isomorphic to the diagram of the mean-field potential gradients and that these nonzero gradients shift the potential minima by intensive amounts, which are interpreted as anharmonic geometry corrections. XVSCF(n) evaluates these mean-field gradients and force constants at the equilibrium geometry and estimates this shift accurately, but approximately, neglecting the coupling between these two quantities. XVSCF[n] solves the coupled equations for geometry corrections and frequencies with an iterative algorithm, giving results that should be identical to those of VSCF when applied to an infinite system. We present the diagrammatic and algebraic definitions, algorithms, and initial implementations as well as numerical results of these two methods. The results show that XVSCF(n) and XVSCF[n] reproduce the vibrationally averaged geometries of VSCF for naphthalene and anthracene in their ground and excited vibrational states accurately at fractions of the computational cost.

The vibrational self-consistent field (VSCF) method is a mean-field approach to solve the vibrational Schrödinger equation and serves as a basis of vibrational perturbation and coupled-cluster methods. Together they account for anharmonic effects on vibrational transition frequencies and vibrationally averaged properties. This article reports the definition, programmable equations, and corresponding initial implementation of a diagrammatically size-extensive modification of VSCF, from which numerous terms with nonphysical size dependence in the original VSCF equations have been eliminated. When combined with a quartic force field (QFF), this compact and strictly size-extensive VSCF (XVSCF) method requires only quartic force constants of the ∂(4)V/∂Q(i)(2)∂Q(j)(2) type, where V is the electronic energy and Q(i) is the ith normal coordinate. Consequently, the cost of a XVSCF calculation with a QFF increases only quadratically with the number of modes, while that of a VSCF calculation grows quartically. The effective (mean-field) potential of XVSCF felt by each mode is shown to be harmonic, making the XVSCF equations subject to a self-consistent analytical solution without matrix diagonalization or a basis-set expansion, which are necessary in VSCF. Even when the same set of force constants is used, XVSCF is nearly three orders of magnitude faster than VSCF implemented similarly. Yet, the results of XVSCF and VSCF are shown to approach each other as the molecular size is increased, implicating the inclusion of unnecessary, nonphysical terms in VSCF. The diagrams of the XVSCF energy expression and their evaluation rules are also proposed, underscoring their connected structures., (© 2011 American Institute of Physics)

A hybrid of the coupled-cluster singles and doubles (CCSD) and second-order Møller-Plesset perturbation (MP2) methods [M. Nooijen, J. Chem. Phys. 111, 10815 (1999); A. D. Bochevarov and C. D. Sherrill, ibid. 122, 234110 (2005); A. D. Bochevarov et al., ibid. 125, 054109 (2006)] is formulated and implemented for one-dimensional periodic extended systems, in which the excitation (T) amplitudes of active bands are determined iteratively by CCSD, while the T amplitudes of mixed active/inactive bands are held fixed at the first-order Møller-Plesset perturbation values. The occupied and virtual bands near the Fermi level, which can cause instability in MP2 when they are (quasi-)degenerate, are selected as active bands to be treated by CCSD, which can, in principle, resist such instability. Two contraction schemes of the T amplitudes (Contractions A and B) are considered. Contraction A is the one proposed for molecules and used also for extended systems because it is efficient for CCSD, but not necessarily so for the hybrid CCSD/MP2. Contraction B is introduced to be more optimally efficient for the hybrid CCSD/MP2 by maximizing the number of intermediate quantities made of the inactive T amplitudes and molecular integrals, which do not vary during CCSD iterations and are computed only once, stored, and reused. In an application to trans-polyacetylene, a smooth transition of the results of the hybrid CCSD/MP2 is observed toward those of CCSD and MP2 by increasing and decreasing, respectively, the number of active bands. With the smallest active space, the hybrid CCSD/MP2 with Contractions A and B achieves a speedup by a factor of 360 and 520, respectively, relative to CCSD. When all of the occupied bands and about half of the virtual bands are active, the hybrid CCSD/MP2 can recover 98% of the CCSD correlation energy or half of the difference between CCSD and MP2 at less than a tenth of the usual CCSD cost., (© 2011 American Institute of Physics)

Size-extensive generalizations of the vibrational self-consistent field (VSCF), vibrational Moller-Plesset perturbation (VMP), and vibrational coupled-cluster (VCC) methods are made to anharmonic lattice vibrations of extended periodic systems on the basis of a quartic force field (QFF) in delocalized normal coordinates. Copious terms in the formalisms of VSCF that have nonphysical size dependence are identified algebraically and eliminated, leading to compact and strictly size-extensive equations. This "quartic" VSCF method (qVSCF) thus defined has no contributions from cubic force constants and alters only the transition energies of the underlying harmonic-oscillator reference from a subset of quartic force constants. It also provides a way to evaluate an anharmonic correction to the lattice structure due to cubic force constants of a certain type. The second-order VMP and VCC methods in the QFF based on the qVSCF reference are shown to account for anharmonic effects due to all cubic and quartic force constants in a size-extensive fashion. These methods can be readily extended to a higher-order truncated Taylor expansion of a potential energy surface in normal coordinates. An algebraic proof of the lack of size-extensivity in the vibrational configuration-interaction method is also presented.

We propose progressive downsampling of wave vectors in the Brillouin zone integrations occurring in the second-order many-body or Møller-Plesset perturbation (MP2) method for extended systems of one-dimensional periodicity. Higher-lying unoccupied and lower-lying occupied Bloch orbitals are subject to downsampling by an exponentially increasing factor (with base n), making the total number of Bloch orbitals included in the MP2 lattice sums grow only logarithmically with respect to the number of basis functions per unit cell. Unlike the mod n downsampling scheme proposed earlier, this log n scheme reduces the scaling of the computational cost and thus achieves a greater speedup as the unit cell size increases. Correct band indexing is essential for accuracy. Two-electron integrals entering the MP2 energy and quasiparticle energy expressions must be multiplied by quadrature weights that are a function of the energy bands involved, and an algorithm to compute the weights is proposed. A combined use of the log n and mod n schemes can speedup the MP2/6-31G** calculation of polyacetylene typically by a factor of 20 with an error in the correlation energy within a few percent relative to the conventional calculation. Similar combinations can reproduce the MP2 quasiparticle energy bands accurately at a fraction of the usual computational cost.

The frequencies of the infrared- and/or Raman-active (k=0) vibrations of polyethylene and polyacetylene are computed by taking account of the anharmonicity in the potential energy surfaces (PESs) and the resulting phonon-phonon couplings explicitly. The electronic part of the calculations is based on Gaussian-basis-set crystalline orbital theory at the Hartree-Fock and second-order Møller-Plesset (MP2) perturbation levels, providing one-, two-, and/or three-dimensional slices of the PES (namely, using the so-called n-mode coupling approximation with n=3), which are in turn expanded in the fourth-order Taylor series with respect to the normal coordinates. The vibrational part uses the vibrational self-consistent field, vibrational MP2, and vibrational truncated configuration-interaction (VCI) methods within the Gamma approximation, which amounts to including only k=0 phonons. It is shown that accounting for both electron correlation and anharmonicity is essential in achieving good agreement (the mean and maximum absolute deviations less than 50 and 90 cm(-1), respectively, for polyethylene and polyacetylene) between computed and observed frequencies. The corresponding values for the calculations including only one of such effects are in excess of 120 and 300 cm(-1), respectively. The VCI calculations also reproduce semiquantitatively the frequency separation and intensity ratio of the Fermi doublet involving the nu(2)(0) fundamental and nu(8)(pi) first overtone in polyethylene.

A formalism for the second-order Møller-Plesset perturbation method employing basis functions that depend explicitly on electron-electron distances (the MP2-R12 or F12 method) is derived and implemented into computer codes for extended systems periodic in one dimension. The excitation amplitudes on these functions are held fixed at values that satisfy the first-order cusp condition. Necessary many-electron integrals over Gaussian-type functions involving Slater-type geminals are evaluated by means of the resolution-of-the-identity approximation with a complementary auxiliary basis set. These integrals and thus the final correlation energy are shown to have the correct size dependence. The valence MP2 correlation energy of polyethylene near the complete basis-set limit is obtained and shown to be considerably greater in magnitude than the value obtained without the R12 treatment.

Coupled-cluster singles and doubles (CCSD) or coupled-cluster singles, doubles, and triples (CCSDT) with noniterative, perturbation corrections for higher-order excitations have been extended to include the basis functions that explicitly depend on interelectronic distances (r(12)) in the wave function expansions with the aim of dramatically accelerating the basis-set convergence of correlation energies. The extension has been based on the so-called R12 (or F12) scheme and applied to a second-order triples correction to CCSD [CCSD(2)(T)-R12], a second-order triples and quadruples correction to CCSD [CCSD(2)(TQ)-R12], a third-order triples correction to CCSD [CCSD(3)(T)-R12], and a second-order quadruples correction to CCSDT [CCSDT(2)(Q)-R12]. A simplified R12 treatment suggested by Fliegl et al. [J. Chem. Phys. 122, 084107 (2005)] has been combined with some of these methods, introducing CCSD(2)(T)(R12) and CCSD(2)(TQ)(R12). The CCSD(T)-R12 method has also been developed as an approximation to CCSD(2)(T)-R12. These methods have been applied to dissociation of hydrogen fluoride and double dissociation of water. For the molecules at their equilibrium geometries, molecular properties predicted by these methods converge extremely rapidly toward the complete-correlation, complete-basis-set limits with respect to the cluster excitation rank, perturbation order, and basis-set size. Although the R12 scheme employed in this work does not improve the basis-set convergence of connected triples or quadruples corrections, the basis-set truncation errors in these contributions have roughly the same magnitude as small residual basis-set truncation errors in the connected singles and doubles contributions even in the dissociation of hydrogen fluoride. In the double dissociation of water, the basis-set truncation errors in the connected triples contribution can be a few times as great as those in the connected singles and doubles contributions.

Efficient computer codes for the explicitly correlated coupled-cluster (CC-R12 or F12) methods with up to triple (CCSDT-R12) and quadruple excitations (CCSDTQ-R12), which take account of the spin, Abelian point-group, and index-permutation symmetries and are based on complete diagrammatic equations, have been implemented with the aid of the computerized symbolic algebra SMITH. Together with the explicitly correlated coupled-cluster singles and doubles (CCSD-R12) method reported earlier [T. Shiozaki et al., J. Chem. Phys. 129, 071101 (2008)], they form a hierarchy of systematic approximations (CCSD-R12

A method for the routine first-principles determination of energies, structures, and phonons of molecular crystals by high-accuracy electron-correlation theories has been proposed. It approximates the energy per unit cell of a crystal by a sum of monomer and dimer energies in an embedding field of self-consistent (and, therefore, polarizable) atomic charges and dipole moments. First and second energy derivatives with respect to atom positions and lattice constants (useful for characterizing structures and phonons) have also been computed efficiently with a long-range electrostatic correction. The method has been applied to solid formic acid modeled as infinite one-dimensional hydrogen-bonded chains. Accurate energies (with corrections for basis-set superposition errors), structural parameters, and frequencies have been obtained for three polymorphic structures (beta(1), beta(2), and alpha) with second-order perturbation theory or higher. On this basis, reliable assignments of their infrared, Raman, and inelastic neutron scattering spectral bands have been proposed. The diffraction and spectroscopic data are shown to be consistent with the pristine beta(1) form and the hitherto-inexplicable infrared band splitting can be assigned to the in-phase and out-of-phase vibrations of adjacent hydrogen-bonded molecules rather than speculated polymorphism. Spectral features expected from the beta(2) and alpha forms have also been predicted and are found to be incompatible with the observed Raman and inelastic neutron scattering spectra in the low-frequency region.

The explicitly correlated coupled-cluster singles and doubles (CCSD-R12) and related methods-its linearized approximation CCSD(R12) and explicitly correlated second-order Moller-Plesset perturbation method-have been implemented into efficient computer codes that take into account point-group symmetry. The implementation has been largely automated by the computerized symbolic algebra SMITH that can handle complex index permutation symmetry of intermediate tensors that occur in the explicitly correlated methods. Unlike prior implementations that invoke the standard approximation or the generalized or extended Brillouin condition, our CCSD-R12 implementation is based on the nontruncated formalisms [T. Shiozaki et al., Phys. Chem. Chem. Phys. 10, 3358 (2008)] in which every diagrammatic term that arises from the modified Ansatz 2 is evaluated either analytically or by the resolution-of-the-identity insertion with the complementary auxiliary basis set. The CCSD-R12 correlation energies presented here for selected systems using the Slater-type correlation function can, therefore, serve as benchmarks for rigorous assessment of other approximate CC-R12 methods. Two recently introduced methods, CCSD(R12) and CCSD(2)(R12), are shown to be remarkably accurate approximations to CCSD-R12.

The anharmonic vibrational frequencies of FHF(-) were computed by the vibrational self-consistent-field, configuration-interaction, and second-order perturbation methods with a multiresolution composite potential energy surface generated by the electronic coupled-cluster method with various basis sets. Anharmonic vibrational averaging was performed for the bond length and nuclear magnetic resonance indirect spin-spin coupling constants, where the latter computed by the equation-of-motion coupled-cluster method. The calculations placed the vibrational frequencies at 580 (nu(1)), 1292 (nu(2)), 1313 (nu(3)), 1837 (nu(1) + nu(3)), and 1864 cm(-1) (nu(1) + nu(2)), the zero-point H-F bond length (r(0)) at 1.1539 A, the zero-point one-bond spin-spin coupling constant [(1)J(0)(HF)] at 124 Hz, and the bond dissociation energy (D(0)) at 43.3 kcal/mol. They agreed excellently with the corresponding experimental values: nu(1) = 583 cm(-1), nu(2) = 1286 cm(-1), nu(3) = 1331 cm(-1), nu(1) + nu(3) = 1849 cm(-1), nu(1) + nu(2) = 1858 cm(-1), r(0) = 1.1522 A, (1)J(0)(HF) = 124+/-3 Hz, and D(0) = 44.4+/-1.6 kcal/mol. The vibrationally averaged bond lengths matched closely the experimental values of five excited vibrational states, furnishing a highly dependable basis for correct band assignments. An adiabatic separation of high- (nu(3)) and low-frequency (nu(1)) stretching modes was examined and found to explain semiquantitatively the appearance of a nu(1) progression on nu(3). Our calculations predicted a value of 186 Hz for experimentally inaccessible (2)J(0)(FF).