PluginPlay: Enabling exascale scientific software one module at a time.

For many computational chemistry packages, being able to efficiently and effectively scale across an exascale cluster is a heroic feat. Collective experience from the Department of Energy's Exascale Computing Project suggests that achieving exascale performance requires far more planning, design, and optimization than scaling to petascale. In many cases, entire rewrites of software are necessary to address fundamental algorithmic bottlenecks. This in turn requires a tremendous amount of resources and development time, resources that cannot reasonably be afforded by every computational science project. It thus becomes imperative that computational science transition to a more sustainable paradigm. Key to such a paradigm is modular software. While the importance of modular software is widely recognized, what is perhaps not so widely appreciated is the effort still required to leverage modular software in a sustainable manner. The present manuscript introduces PluginPlay, https://github.com/NWChemEx-Project/PluginPlay, an inversion-of-control framework designed to facilitate developing, maintaining, and sustaining modular scientific software packages. This manuscript focuses on the design aspects of PluginPlay and how they specifically influence the performance of the resulting package. Although, PluginPlay serves as the framework for the NWChemEx package, PluginPlay is not tied to NWChemEx or even computational chemistry. We thus anticipate PluginPlay to prove to be a generally useful tool for a number of computational science packages looking to transition to the exascale.

Gaussian approximation potentials for accurate thermal properties of two-dimensional materials.

Two-dimensional materials (2DMs) continue to attract a lot of attention, particularly for their extreme flexibility and superior thermal properties. Molecular dynamics simulations are among the most powerful methods for computing these properties, but their reliability depends on the accuracy of interatomic interactions. While first principles approaches provide the most accurate description of interatomic forces, they are computationally expensive. In contrast, classical force fields are computationally efficient, but have limited accuracy in interatomic force description. Machine learning interatomic potentials, such as Gaussian Approximation Potentials, trained on density functional theory (DFT) calculations offer a compromise by providing both accurate estimation and computational efficiency. In this work, we present a systematic procedure to develop Gaussian approximation potentials for selected 2DMs, graphene, buckled silicene, and h-XN (X = B, Al, and Ga, as binary compounds) structures. We validate our approach through calculations that require various levels of accuracy in interatomic interactions. The calculated phonon dispersion curves and lattice thermal conductivity, obtained through harmonic and anharmonic force constants (including fourth order) are in excellent agreement with DFT results. HIPHIVE calculations, in which the generated GAP potentials were used to compute higher-order force constants instead of DFT, demonstrated the first-principles level accuracy of the potentials for interatomic force description. Molecular dynamics simulations based on phonon density of states calculations, which agree closely with DFT-based calculations, also show the success of the generated potentials in high-temperature simulations.

High-Accuracy Heats of Formation for Alkane Oxidation: From Small to Large via the Automated CBH-ANL Method.

It is generally challenging to obtain high-accuracy predictions for the heat of formation for species with more than a handful of heavy atoms, such as those of importance in standard combustion mechanisms. To this end, we construct the CBH-ANL approach and illustrate that, for a set of 194 alkane oxidation species, it can be used to produce ΔHf(0 K) values with 2σ uncertainties of 0.2-0.5 kcal mol-1. This set includes the alkanes, hydroperoxides, and alkyl, peroxy, and hydroperoxyalkyl radicals for 17 representative hydrocarbon fuels containing up to 10 heavy atoms with various degrees of branching in the alkane backbone. The CBH-ANL approach, automated in the QTC and AutoMech software suites, builds balanced chemical equations for the calculation of ΔHf(0 K), in which the reference species may be up to five heavy atoms. The high-level ANL0 and ANL1 reference ΔHf(0 K) values are further refined for even the largest of these reference species with a novel laddering approach. We perform a comprehensive quantification of the uncertainties for both the individual reference species (the largest of which is 0.15 kcal mol-1) and the propagation of those uncertainties when used in the calculation of ΔHf(0 K) for the 194 target species. We examine the sensitivity of the predicted ΔHf(0 K) values to (i) electronic energies from various methods, including ωB97X-D/cc-pVTZ, B2PLYP-D3/cc-pVTZ, CCSD(T)-F12b/cc-pVDZ-F12//B2PLYP-D3/cc-pVTZ, and CCSD(T)-F12b/cc-pVTZ-F12//B2PLYP-D3/cc-pVTZ; (ii) the zero-point vibrational energies (ZPVEs), where we consider harmonic ZPVEs as well as two scaling-based estimates of the anharmonic ZPVEs, all implemented for both ωB97X-D/cc-pVTZ and B2PLYP-D3/cc-pVTZ calculations; (iii) the particular CBH-ANL scheme employed; and (iv) the procedure for choosing the reference conformer for the analyses. The discussion concludes with a summary of the estimated overall uncertainty in the predictions and a validation of the predictions for the alkane subset.

Prediction of correlation energies using variational subspace valence bond.

In the variational subspace valence bond (VSVB) [G. D. Fletcher, J. Chem. Phys. 142, 134112 (2015)] method, the electronic orbitals comprising the wave function correspond to chemically meaningful objects, such as bonds, lone pairs, atomic cores, and so on. Selected regions of a molecule (for example, a single chemical bond, as opposed to all the bonds) can be modeled with different levels of basis set and possible methods for modeling correlation from the other regions. The interactions between the components of a molecule (say, a bond and a neighboring orbital) can then be studied in detail for their impact on a chemical phenomenon while avoiding the expense of necessarily applying the higher levels and methods to the entire molecule. This work presents the theoretical basis for modeling correlation effects between specific electron pairs by incorporating terms in the inter-electronic coordinates ("r12") into VSVB. The approach is validated with calculations on small systems using single-reference wave functions.

Exploring the anharmonic vibrational structure of carbon dioxide trimers.

Our previously developed mbCO2 potential [O. Sode and J. N. Cherry, J. Comput. Chem. 38, 2763 (2017)] is used to describe the vibrational structure of the intermolecular motions of the CO2 trimers: barrel-shaped and cyclic trimers. Anharmonic corrections are accounted for using the vibrational self-consistent field theory, vibrational second-order Møller-Plesset perturbation (VMP2) theory, and vibrational configuration interaction (VCI) methods and compared with experimental observations. For the cyclic structure, we revise the assignments of two previously observed experimental peaks based on our VCI and VMP2 results. We note that the experimental band observed near 13 cm-1 is the out-of-phase out-of-plane degenerate motion with E″ symmetry, while the peak observed at 18 cm-1 likely corresponds to the symmetric out-of-plane torsion A″ vibration. Since the VCI treatment of the vibrational motions accounts for vibrational mixing and delocalization, overtones and combination bands were also observed and quantified in the intermolecular regions of the two trimer isomers.

From NWChem to NWChemEx: Evolving with the Computational Chemistry Landscape.

Since the advent of the first computers, chemists have been at the forefront of using computers to understand and solve complex chemical problems. As the hardware and software have evolved, so have the theoretical and computational chemistry methods and algorithms. Parallel computers clearly changed the common computing paradigm in the late 1970s and 80s, and the field has again seen a paradigm shift with the advent of graphical processing units. This review explores the challenges and some of the solutions in transforming software from the terascale to the petascale and now to the upcoming exascale computers. While discussing the field in general, NWChem and its redesign, NWChemEx, will be highlighted as one of the early codesign projects to take advantage of massively parallel computers and emerging software standards to enable large scientific challenges to be tackled.

Understanding the anharmonic vibrational structure of the carbon dioxide dimer.

Understanding the vibrational structure of the CO2 system is important to confirm the potential energy surface and interactions in such van der Waals complexes. In this work, we use our previously developed mbCO2 potential function to explore the vibrational structure of the CO2 monomer and dimer. The potential function has been trained to reproduce the potential energies at the CCSD(T)-F12b/aug-cc-pVTZ level of electronic structure theory. The harmonic approximation, as well as anharmonic corrections using vibrational structure theories such as vibrational self-consistent field, vibrational second-order Møller-Plesset perturbation, and vibrational configuration interaction (VCI), is applied to address the vibrational motions. We compare the vibrational results using the mbCO2 potential function with traditional electronic structure theory results and to experimental frequencies. The anharmonic results for the monomer most closely match the experimental data to within 3 cm-1, including the Fermi dyad frequencies. The intermolecular and intramolecular dimer frequencies were treated separately and show good agreement with the most recent theoretical and experimental results from the literature. The VCI treatment of the dimer vibrational motions accounts for vibrational mixing and delocalization, such that we observe the dimer Fermi resonance phenomena, both in the intramolecular and intermolecular regions.

Valence : A Massively Parallel Implementation of the Variational Subspace Valence Bond Method.

This work describes the software package, Valence, for the calculation of molecular energies using the variational subspace valence bond (VSVB) method. VSVB is an ab initio electronic structure method based on nonorthogonal orbitals. Important features of practical value include high parallel scalability, wave functions that can be constructed automatically by combining orbitals from previous calculations, and ground and excited states that can be modeled with a single configuration or determinant. The open-source software package includes tools to generate wave functions, a database of generic orbitals, example input files, and a library build intended for integration with other packages. We also describe the interface to an external software package, enabling the computation of optimized molecular geometries and vibrational frequencies. © 2019 Wiley Periodicals, Inc.

SIESTA-SIPs: Massively parallel spectrum-slicing eigensolver for an ab initio molecular dynamics package.

Integration of Shift-and-Invert Parallel Spectral Transformation (SIPs) eigensolver (as implemented in the SLEPc library) into an ab initio molecular dynamics package, SIESTA, is described. The effectiveness of the code is demonstrated on applications to polyethylene chains, boron nitride sheets, and bulk water clusters. For problems with the same number of orbitals, the performance of the SLEPc eigensolver depends on the sparsity of the matrices involved, favoring reduced dimensional systems such as polyethylene or boron nitride sheets in comparison to bulk systems like water clusters. For all problems investigated, performance of SIESTA-SIPs exceeds the performance of SIESTA with default solver (ScaLAPACK) at the larger number of cores and the larger number of orbitals. A method that improves the load-balance with each iteration in the self-consistency cycle by exploiting the emerging knowledge of the eigenvalue spectrum is demonstrated. © 2018 Wiley Periodicals, Inc.

The influence of surface functionalization on thermal transport and thermoelectric properties of MXene monolayers.

The newest members of a two-dimensional material family, involving transition metal carbides and nitrides (called MXenes), have garnered increasing attention due to their tunable electronic and thermal properties depending on the chemical composition and functionalization. This flexibility can be exploited to fabricate efficient electrochemical energy storage (batteries) and energy conversion (thermoelectric) devices. In this study, we calculated the Seebeck coefficients and lattice thermal conductivity values of oxygen terminated M2CO2 (where M = Ti, Zr, Hf, Sc) monolayer MXene crystals in two different functionalization configurations (model-II (MD-II) and model-III (MD-III)), using density functional theory and Boltzmann transport theory. We estimated the thermoelectric figure-of-merit, zT, of these materials by two different approaches, as well. First of all, we found that the structural model (i.e. adsorption site of oxygen atom on the surface of MXene) has a paramount impact on the electronic and thermoelectric properties of MXene crystals, which can be exploited to engineer the thermoelectric properties of these materials. The lattice thermal conductivity κl, Seebeck coefficient and zT values may vary by 40% depending on the structural model. The MD-III configuration always has the larger band gap, Seebeck coefficient and zT, and smaller κl as compared to the MD-II structure due to a larger band gap, highly flat valence band and reduced crystal symmetry in the former. The MD-III configuration of Ti2CO2 and Zr2CO2 has the lowest κl as compared to the same configuration of Hf2CO2 and Sc2CO2. Among all the considered structures, the MD-II configuration of Hf2CO2 has the highest κl, and Ti2CO2 and Zr2CO2 in the MD-III configuration have the lowest κl. For instance, while the band gap of the MD-II configuration of Ti2CO2 is 0.26 eV, it becomes 0.69 eV in MD-III. The zTmax value may reach up to 1.1 depending on the structural model of MXene.

Benchmark Calculations of Energetic Properties of Groups 4 and 6 Transition Metal Oxide Nanoclusters Including Comparison to Density Functional Theory.

The heats of formation and the normalized clustering energies (NCEs) for the group 4 and group 6 transition metal oxide (TMO) trimers and tetramers have been calculated by the Feller-Peterson-Dixon (FPD) method. The heats of formation predicted by the FPD method do not differ much from those previously derived from the NCEs at the CCSD(T)/aT level except for the CrO3 nanoclusters. New and improved heats of formation for Cr3O9 and Cr4O12 were obtained using PW91 orbitals instead of Hartree-Fock (HF) orbitals. Diffuse functions are necessary to predict accurate heats of formation. The fluoride affinities (FAs) are calculated with the CCSD(T) method. The relative energies (REs) of different isomers, NCEs, electron affinities (EAs), and FAs of (MO2)n (M = Ti, Zr, Hf, n = 1-4) and (MO3)n (M = Cr, Mo, W, n = 1-3) clusters have been benchmarked with 55 exchange-correlation density functional theory (DFT) functionals including both pure and hybrid types. The absolute errors of the DFT results are mostly less than ±10 kcal/mol for the NCEs and the EAs and less than ±15 kcal/mol for the FAs. Hybrid functionals usually perform better than the pure functionals for the REs and NCEs. The performance of the two types of functionals in predicting EAs and FAs is comparable. The B1B95 and PBE1PBE functionals provide reliable energetic properties for most isomers. Long range corrected pure functionals usually give poor FAs. The standard deviation of the absolute error is always close to the mean errors, and the probability distributions of the DFT errors are often not Gaussian (normal). The breadth of the distribution of errors and the maximum probability are dependent on the energy property and the isomer.

Shift-and-invert parallel spectral transformation eigensolver: Massively parallel performance for density-functional based tight-binding.

The Shift-and-invert parallel spectral transformations (SIPs), a computational approach to solve sparse eigenvalue problems, is developed for massively parallel architectures with exceptional parallel scalability and robustness. The capabilities of SIPs are demonstrated by diagonalization of density-functional based tight-binding (DFTB) Hamiltonian and overlap matrices for single-wall metallic carbon nanotubes, diamond nanowires, and bulk diamond crystals. The largest (smallest) example studied is a 128,000 (2000) atom nanotube for which ∼330,000 (∼5600) eigenvalues and eigenfunctions are obtained in ∼190 (∼5) seconds when parallelized over 266,144 (16,384) Blue Gene/Q cores. Weak scaling and strong scaling of SIPs are analyzed and the performance of SIPs is compared with other novel methods. Different matrix ordering methods are investigated to reduce the cost of the factorization step, which dominates the time-to-solution at the strong scaling limit. A parallel implementation of assembling the density matrix from the distributed eigenvectors is demonstrated.

Fermi resonance in solid CO2 under pressure.

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.

Optimized coordinates for anharmonic vibrational structure theories.

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.

Size-extensive vibrational self-consistent field methods with anharmonic geometry corrections.

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.

Extensivity of energy and electronic and vibrational structure methods for crystals.

A pedagogical proof is presented for the extensivity of energies of metallic and nonmetallic crystals that proceeds by elucidating the asymptotic distance dependence of the effective chemical interactions: kinetic, Coulomb, exchange, and correlation. On this basis, a guideline for the size-consistent design of electronic and vibrational methods is proposed. This guideline underscores the significance of the distinct use of the intermediate and standard normalization of wave functions for extensive and intensive quantities, includes the extensive and intensive diagram theorems as the unambiguous criteria for determining size consistency of a method for extensive and intensive quantities, and introduces the extensive-intensive consistency theorem, which stipulates the precise balance between the determinant spaces reached by extensive and intensive operators. Electronic and vibrational methods for crystals are reviewed that are inspired by these formal analyses or developed in accordance with the guideline.

Size-extensive vibrational self-consistent field method.

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.

First-principles theories for anharmonic lattice vibrations.

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.

First-principles calculations on anharmonic vibrational frequencies of polyethylene and polyacetylene in the Gamma approximation.

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.