Fast and accurate quantum Monte Carlo for molecular crystals

Beyond Single-Reference Fixed-Node Approximation in Ab Initio Diffusion Monte Carlo Using Antisymmetrized Geminal Power Applied to Systems with Hundreds of Electrons

Diffusion Monte Carlo (DMC) is an exact technique to project out the ground state (GS) of a Hamiltonian. Since the GS is always bosonic, in Fermionic systems, the projection needs to be carried out while imposing antisymmetric constraints, which is a nondeterministic polynomial hard problem. In practice, therefore, the application of DMC on electronic structure problems is made by employing the fixed-node (FN) approximation, consisting of performing DMC with the constraint of having a fixed, predefined nodal surface. How do we get the nodal surface? The typical approach, applied in systems having up to hundreds or even thousands of electrons, is to obtain the nodal surface from a preliminary mean-field approach (typically, a density functional theory calculation) used to obtain a single Slater determinant. This is known as single reference. In this paper, we propose a new approach, applicable to systems as large as the C60 fullerene, which improves the nodes by going beyond the single reference. In practice, we employ an implicitly multireference ansatz (antisymmetrized geminal power wave function constraint with molecular orbitals), initialized on the preliminary mean-field approach, which is relaxed by optimizing a few parameters of the wave function determining the nodal surface by minimizing the FN-DMC energy. We highlight the improvements of the proposed approach over the standard single-reference method on several examples and, where feasible, the computational gain over the standard multireference ansatz, which makes the methods applicable to large systems. We also show that physical properties relying on relative energies, such as binding energies, are affordable and reliable within the proposed scheme.

First-principles spectroscopy of aqueous interfaces using machine-learned electronic and quantum nuclear effects

Vibrational spectroscopy is a powerful approach to visualising interfacial phenomena. However, extracting structural and dynamical information from vibrational spectra is a challenge that requires first-principles simulations, including non-Condon and quantum nuclear effects. We address this challenge by developing a machine-learning enhanced first-principles framework to speed up predictive modelling of infrared, Raman, and sum-frequency generation spectra. Our approach uses machine learning potentials that encode quantum nuclear effects to generate quantum trajectories using simple molecular dynamics efficiently. In addition, we reformulate bulk and interfacial selection rules to express them unambiguously in terms of the derivatives of polarisation and polarisabilities of the whole system and predict these derivatives efficiently using fully-differentiable machine learning models of dielectric response tensors. We demonstrate our framework's performance by predicting the IR, Raman, and sum-frequency generation spectra of liquid water, ice and the water-air interface by achieving near quantitative agreement with experiments at nearly the same computational efficiency as pure classical methods. Finally, to aid the experimental discovery of new phases of nanoconfined water, we predict the temperature-dependent vibrational spectra of monolayer water across the solid-hexatic-liquid phases transition.

Affinity of Telluronium Chalcogen Bond Donors for Lewis Bases in Solution: A Critical Experimental-Theoretical Joint Study

Telluronium salts [Ar2 MeTe]X were synthesized, and their Lewis acidic properties towards a number of Lewis bases were addressed in solution by physical and theoretical means. Structural X-ray diffraction analysis of 21 different salts revealed the electrophilicity of the Te centers in their interactions with anions. Telluroniums' propensity to form Lewis pairs was investigated with OPPh3 . Diffusion-ordered NMR spectroscopy suggested that telluroniums can bind up to three OPPh3 molecules. Isotherm titration calorimetry showed that the related heats of association in 1,2-dichloroethane depend on the electronic properties of the substituents of the aryl moiety and on the nature of the counterion. The enthalpies of first association of OPPh3 span -0.5 to -5 kcal mol-1 . Study of the affinity of telluroniums for OPPh3 by state-of-the-art DFT and ab-initio methods revealed the dominant Coulombic and dispersion interactions as well as an entropic effect favoring association in solution. Intermolecular orbital interactions between [Ar2 MeTe]+ cations and OPPh3 are deemed insufficient on their own to ensure the cohesion of [Ar2 MeTe ⋅ Bn ]+ complexes in solution (B=Lewis base). Comparison of Grimme's and Tkatchenko's DFT-D4/MBD-vdW thermodynamics of formation of higher [Ar2 MeTe ⋅ Bn ]+ complexes revealed significant molecular size-dependent divergence of the two methodologies, with MBD yielding better agreement with experiment.

Toward accurate modeling of structure and energetics of bulk hexagonal boron nitride

Materials that exhibit both strong covalent and weak van der Waals interactions pose a considerable challenge to many computational methods, such as DFT. This makes assessing the accuracy of calculated properties, such as exfoliation energies in layered materials like hexagonal boron nitride (h-BN) problematic, when experimental data are not available. In this paper, we investigate the accuracy of equilibrium lattice constants and exfoliation energy calculation for various DFT-based computational approaches in bulk h-BN. We contrast these results with available experiments and reference fixed-node diffusion quantum Monte Carlo (QMC) results. From our reference QMC calculation, we obtained an exfoliation energy of - 33 ± 2 meV/atom (-0.38 ± 0.02 J/m2 ).

Many-Body Methods for Surface Chemistry Come of Age: Achieving Consensus with Experiments

The adsorption energy of a molecule onto the surface of a material underpins a wide array of applications, spanning heterogeneous catalysis, gas storage, and many more. It is the key quantity where experimental measurements and theoretical calculations meet, with agreement being necessary for reliable predictions of chemical reaction rates and mechanisms. The prototypical molecule-surface system is CO adsorbed on MgO, but despite intense scrutiny from theory and experiment, there is still no consensus on its adsorption energy. In particular, the large cost of accurate many-body methods makes reaching converged theoretical estimates difficult, generating a wide range of values. In this work, we address this challenge, leveraging the latest advances in diffusion Monte Carlo (DMC) and coupled cluster with single, double, and perturbative triple excitations [CCSD(T)] to obtain accurate predictions for CO on MgO. These reliable theoretical estimates allow us to evaluate the inconsistencies in published temperature-programed desorption experiments, revealing that they arise from variations in employed pre-exponential factors. Utilizing this insight, we derive new experimental estimates of the (electronic) adsorption energy with a (more) precise pre-exponential factor. As a culmination of all of this effort, we are able to reach a consensus between multiple theoretical calculations and multiple experiments for the first time. In addition, we show that our recently developed cluster-based CCSD(T) approach provides a low-cost route toward achieving accurate adsorption energies. This sets the stage for affordable and reliable theoretical predictions of chemical reactions on surfaces to guide the realization of new catalysts and gas storage materials.

H2O2(s) and H2O2·2H2O(s) crystals compared with ices: DFT functional assessment and D3 analysis

The H2O and H2O2 molecules resemble each other in a multitude of ways as has been noted in the literature. Here, we present density functional theory (DFT) calculations for the H2O2(s) and H2O2·2H2O(s) crystals and make selected comparisons with ice polymorphs. The performance of a number of dispersion-corrected density functionals-both self-consistent and a posteriori ones-are assessed, and we give special attention to the D3 correction and its effects. The D3 correction to the lattice energies is large: for H2O2(s) the D3 correction constitutes about 25% of the lattice energy using PBE, much more for RPBE, much less for SCAN, and it primarily arises from non-H-bonded interactions out to about 5 Å.The large D3 corrections to the lattice energies are likely a consequence of several effects: correction for missing dispersion interaction, the ability of D3 to capture and correct various other kinds of limitations built into the underlying DFT functionals, and finally some degree of cell-contraction-induced polarization enhancement. We find that the overall best-performing functionals of the twelve examined are optPBEvdW and RPBE-D3. Comparisons with DFT assessments for ices in the literature show that where the same methods have been used, the assessments largely agree.

Physical binding energies using the electron localization function in 4-hydroxyphenylboronic acid co-crystals with aza donors

Binding energies are traditionally simulated using cluster models by computation of each synthon for each individual co-crystal former. However, our investigation of the binding strengths using the electron localization function (ELF) reveals that these can be determined directly from the crystal supercell computations. We propose a new modeling protocol for the computation of physical binding energies directly from bulk simulations using ELF analysis. In this work, we establish a correlation between ELF values and binding energies calculated for co-crystals of 4-hydroxyphenylboronic acid (4HPBA) with four different aza donors using density functional theory with varying descriptions of dispersion. Boronic acids are gaining significant interest in the field of crystal engineering, but theoretical studies on their use in materials are still very limited. Here, we present a systematic investigation of the non-covalent interactions in experimentally realized co-crystals. Prior diffraction studies on these complexes have shown the competitive nature between the boronic acid functional group and the para-substituted phenolic group forming heteromeric interactions with aza donors. We determine the stability of the co-crystals by simulating their lattice energies, and the different dispersion descriptions show similar trends in lattice energies and lattice parameters. Our study bolsters the experimental observation of the boronic acid group as a competitive co-crystal former in addition to the well-studied phenolic group. Further research on correlating ELF values for physical binding could potentially transform this approach to a viable alternative for the computation of binding energies.

New Density Matrix Renormalization Group Approaches for Strongly Correlated Systems Coupled with Large Environments

Thanks to the high compression of the matrix product state (MPS) form of the wave function and the efficient site-by-site iterative sweeping optimization algorithm, the density matrix normalization group (DMRG) and its time-dependent variant (TD-DMRG) have been established as powerful computational tools in accurately simulating the electronic structure and quantum dynamics of strongly correlated molecules with a large number (10^1-2) of quantum degrees of freedom (active orbitals or vibrational modes). However, the quantitative characterization of the quantum many-body behaviors of realistic strongly correlated systems requires a further consideration of the interaction between the embedded active subsystem and the remaining correlated environment, e.g., a larger number (10^2-3) of external orbitals in electronic structure or infinite condensed-phase phononic modes in nucleus dynamics. To this end, we introduced three new post-DMRG and TD-DMRG approaches, namely (1) DMRG2sCI-MRCI and DMRG2sCI-ENPT by the reconstruction of selected configuration interaction (sCI) type of compact reference function from DMRG coefficients and the use of externally contracted MRCI (multireference configuration interaction) and Epstein-Nesbet perturbation theory (ENPT), without recourse to the expensive high order n-electron reduced density matrices (n-RDMs). (2) DMRG combined with RR-MRCI (renormalized residue-based MRCI), which improves the computational accuracy and efficiency of internally contracted (ic) MRCI by renormalizing the contracted bases with small-sized buffer environment(s) of a few external orbitals as probes based on quantum information theory. (3) HM (hierarchical mapping)-TD-DMRG in which a large environment is reduced to a small number of renormalized environmental modes (which accounts for the most vital system-environment interactions) through stepwise mapping transformation. These advances extend the efficacy of highly accurate DMRG/TD-DMRG computations to the quantitative characterization of the electronic structure and quantum dynamics in realistic strongly correlated systems coupled with large environments and are reviewed in this paper.

Quantifying the Impact of Halogenation on Intermolecular Interactions and Binding Modes of Aromatic Molecules

Halogenation of aromatic molecules is frequently used to modulate intermolecular interactions with ramifications for optoelectronic and mechanical properties. In this work, we accurately quantify and understand the nature of intermolecular interactions in perhalogenated benzene (PHB) clusters. Using benchmark binding energies from the fixed-node diffusion Monte Carlo (FN-DMC) method, we show that generalized Kohn-Sham semicanonical projected random phase approximation (GKS-spRPA) plus approximate exchange kernel (AKX) provides reliable interaction energies with mean absolute error (MAE) of 0.23 kcal/mol. Using the GKS-spRPA+AXK method, we quantify the interaction energies of several binding modes of PHB clusters ((C₆X₆)_n; X = F, Cl, Br, I; n = 2, 3). For a given binding mode, the interaction energies increase 3-4 times from X = F to X = I; the X-X binding modes have energies in the range of 2-4 kcal/mol, while the π-π binding mode has interaction energies in the range of 4-12 kcal/mol. SAPT-DFT-based energy decomposition analysis is then used to show that the equilibrium geometries are dictated primarily by the dispersion and exchange interactions. Finally, we test the accuracy of several dispersion-corrected density functional approximations and show that only the r2SCAN-D4 method has a low MAE and correct long-range behavior, which makes it suitable for large-scale simulations and for developing structure-function relationships of halogenated aromatic systems.

Discriminating and understanding molecular crystal polymorphism

Polymorph discrimination for a molecular crystal has long been a challenging task, which, nonetheless, is a major concern in the pharmaceutical industry. In this work, we have investigated polymorph discrimination on three different molecular crystals, tetrolic acid, oxalic acid, and oxalyl dihydrazide, covering both packing polymorphism and conformational polymorphism. To gain more understanding, we have performed energy decomposition analysis based on many-body expansion, and have compared the results from the XO-PBC method, that is, the eXtended ONIOM method (XO) with the periodic boundary condition (PBC), with those from some commonly used dispersion corrected density functional theory (DFT-D) methods. It is shown here that, with the XYG3 doubly hybrid functional chosen as the target high level to capture the intra- and short-range intermolecular interactions, and the periodic PBE as the basic low level to take long range interactions into account, the XO-PBC(XYG3:PBE) method not only obtains the correct experimental stability orderings, but also predicts reasonable polymorph energy ranges for all three cases. Our results have demonstrated the usefulness of the present theoretical methods, in particular XO-PBC, while highlighted the importance of a better treatment of different kinds of interactions to be beneficial to polymorph control.

The first-principles phase diagram of monolayer nanoconfined water

Water in nanoscale cavities is ubiquitous and of central importance to everyday phenomena in geology and biology. However, the properties of nanoscale water can be substantially different from those of bulk water, as shown, for example, by the anomalously low dielectric constant of water in nanochannels¹, near frictionless water flow² or the possible existence of a square ice phase³. Such properties suggest that nanoconfined water could be engineered for technological applications in nanofluidics⁴, electrolyte materials⁵ and water desalination⁶. Unfortunately, challenges in experimentally characterizing water at the nanoscale and the high cost of first-principles simulations have prevented the molecular-level understanding required to control the behaviour of water. Here we combine a range of computational approaches to enable a first-principles-level investigation of a single layer of water within a graphene-like channel. We find that monolayer water exhibits surprisingly rich and diverse phase behaviour that is highly sensitive to temperature and the van der Waals pressure acting within the nanochannel. In addition to multiple molecular phases with melting temperatures varying non-monotonically by more than 400 kelvins with pressure, we predict a hexatic phase, which is an intermediate between a solid and a liquid, and a superionic phase with a high electrical conductivity exceeding that of battery materials. Notably, this suggests that nanoconfinement could be a promising route towards superionic behaviour under easily accessible conditions.

BIGDML-Towards accurate quantum machine learning force fields for materials

Machine-learning force fields (MLFF) should be accurate, computationally and data efficient, and applicable to molecules, materials, and interfaces thereof. Currently, MLFFs often introduce tradeoffs that restrict their practical applicability to small subsets of chemical space or require exhaustive datasets for training. Here, we introduce the Bravais-Inspired Gradient-Domain Machine Learning (BIGDML) approach and demonstrate its ability to construct reliable force fields using a training set with just 10-200 geometries for materials including pristine and defect-containing 2D and 3D semiconductors and metals, as well as chemisorbed and physisorbed atomic and molecular adsorbates on surfaces. The BIGDML model employs the full relevant symmetry group for a given material, does not assume artificial atom types or localization of atomic interactions and exhibits high data efficiency and state-of-the-art energy accuracies (errors substantially below 1 meV per atom) for an extended set of materials. Extensive path-integral molecular dynamics carried out with BIGDML models demonstrate the counterintuitive localization of benzene-graphene dynamics induced by nuclear quantum effects and their strong contributions to the hydrogen diffusion coefficient in a Pd crystal for a wide range of temperatures.

A complete description of thermodynamic stabilities of molecular crystals

Predicting stable polymorphs of molecular crystals remains one of the grand challenges of computational science. Current methods invoke approximations to electronic structure and statistical mechanics and thus fail to consistently reproduce the delicate balance of physical effects determining thermodynamic stability. We compute the rigorous ab initio Gibbs free energies for competing polymorphs of paradigmatic compounds, using machine learning to mitigate costs. The accurate description of electronic structure and full treatment of quantum statistical mechanics allow us to predict the experimentally observed phase behavior. This constitutes a key step toward the first-principles design of functional materials for applications from photovoltaics to pharmaceuticals.

Predictions of relative stabilities of (competing) molecular crystals are of great technological relevance, most notably for the pharmaceutical industry. However, they present a long-standing challenge for modeling, as often minuscule free energy differences are sensitively affected by the description of electronic structure, the statistical mechanics of the nuclei and the cell, and thermal expansion. The importance of these effects has been individually established, but rigorous free energy calculations for general molecular compounds, which simultaneously account for all effects, have hitherto not been computationally viable. Here we present an efficient “end to end” framework that seamlessly combines state-of-the art electronic structure calculations, machine-learning potentials, and advanced free energy methods to calculate ab initio Gibbs free energies for general organic molecular materials. The facile generation of machine-learning potentials for a diverse set of polymorphic compounds—benzene, glycine, and succinic acid—and predictions of thermodynamic stabilities in qualitative and quantitative agreement with experiments highlight that predictive thermodynamic studies of industrially relevant molecular materials are no longer a daunting task.

A pathway toward high-throughput quantum Monte Carlo simulations for alloys: A case study of two-dimensional (2D) GaSxSe1-x

The study of alloys using computational methods has been a difficult task due to the usually unknown stoichiometry and local atomic ordering of the different structures experimentally. In order to combat this, first-principles methods have been coupled with statistical methods such as the cluster expansion formalism in order to construct the energy hull diagram, which helps to determine if an alloyed structure can exist in nature. Traditionally, density functional theory (DFT) has been used in such workflows. In this paper, we propose to use chemically accurate many-body variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC) methods to construct the energy hull diagram of an alloy system due to the fact that such methods have a weaker dependence on the starting wavefunction and density functional, scale similarly to DFT with the number of electrons, and have had demonstrated success for a variety of materials. To carry out these simulations in a high-throughput manner, we propose a method called Jastrow sharing, which involves recycling the optimized Jastrow parameters between alloys with different stoichiometries. We show that this eliminates the need for extra VMC Jastrow optimization calculations and results in significant computational cost savings (on average 1/4 savings of total computational time). Since it is a novel post-transition metal chalcogenide alloy series that has been synthesized in its few-layer form, we used monolayer GaS_xSe_1-x as a case study for our workflow. By extensively testing our Jastrow sharing procedure for monolayer GaS_xSe_1-x and quantifying the cost savings, we demonstrate how a pathway toward chemically accurate high-throughput simulations of alloys can be achieved using many-body VMC and DMC methods.

Monte Carlo MP2-F12 for Noncovalent Interactions: The C60 Dimer

A scalable stochastic algorithm is presented that can evaluate explicitly correlated (F12) second-order many-body perturbation (MP2) energies of weak, noncovalent, intermolecular interactions. It first transforms the formulas of the MP2 and F12 energy differences into a short sum of high-dimensional integrals of Green's functions in real space and imaginary time. These integrals are then evaluated by the Monte Carlo method augmented by parallel execution, redundant-walker convergence acceleration, direct-sampling autocorrelation elimination, and control-variate error reduction. By sharing electron-pair walkers across the supermolecule and its subsystems spanned by the joint basis set, the statistical uncertainty is reduced by one to 2 orders of magnitude in the MP2 binding energy corrected for the basis-set incompleteness and superposition errors. The method predicts the MP2-F12/aug-cc-pVDZ binding energy of 19.1 ± 4.0 kcal mol^-1 for the C₆₀ dimer at the center distance of 9.748 Å.

Importance of Nuclear Quantum Effects for NMR Crystallography

The resolving power of solid-state nuclear magnetic resonance (NMR) crystallography depends heavily on the accuracy of computational predictions of NMR chemical shieldings of candidate structures, which are usually taken to be local minima in the potential energy. To test the limits of this approximation, we systematically study the importance of finite-temperature and quantum nuclear fluctuations for ¹H, ¹³C, and ¹⁵N shieldings in polymorphs of three paradigmatic molecular crystals: benzene, glycine, and succinic acid. The effect of quantum fluctuations is comparable to the typical errors of shielding predictions for static nuclei with respect to experiments, and their inclusion improves the agreement with measurements, translating to more reliable assignment of the NMR spectra to the correct candidate structure. The use of integrated machine-learning models, trained on first-principles energies and shieldings, renders rigorous sampling of nuclear fluctuations affordable, setting a new standard for the calculations underlying NMR structure determinations.

Interactions between large molecules pose a puzzle for reference quantum mechanical methods

Quantum-mechanical methods are used for understanding molecular interactions throughout the natural sciences. Quantum diffusion Monte Carlo (DMC) and coupled cluster with single, double, and perturbative triple excitations [CCSD(T)] are state-of-the-art trusted wavefunction methods that have been shown to yield accurate interaction energies for small organic molecules. These methods provide valuable reference information for widely-used semi-empirical and machine learning potentials, especially where experimental information is scarce. However, agreement for systems beyond small molecules is a crucial remaining milestone for cementing the benchmark accuracy of these methods. We show that CCSD(T) and DMC interaction energies are not consistent for a set of polarizable supramolecules. Whilst there is agreement for some of the complexes, in a few key systems disagreements of up to 8 kcal mol^-1 remain. These findings thus indicate that more caution is required when aiming at reproducible non-covalent interactions between extended molecules.

Crystal Structure Prediction Methods for Organic Molecules: State of the Art

The prediction of the crystal structures that a given organic molecule is likely to form is an important theoretical problem of significant interest for the pharmaceutical and agrochemical industries, among others. As evidenced by a series of six blind tests organized over the past 2 decades, methodologies for crystal structure prediction (CSP) have witnessed substantial progress and have now reached a stage of development where they can begin to be applied to systems of practical significance. This article reviews the state of the art in general-purpose methodologies for CSP, placing them within a common framework that highlights both their similarities and their differences. The review discusses specific areas that constitute the main focus of current research efforts toward improving the reliability and widening applicability of these methodologies, and offers some perspectives for the evolution of this technology over the next decade.

r2SCAN-3c: A "Swiss army knife" composite electronic-structure method

The recently proposed r²SCAN meta-generalized-gradient approximation (mGGA) of Furness and co-workers is used to construct an efficient composite electronic-structure method termed r²SCAN-3c. To this end, the unaltered r²SCAN functional is combined with a tailor-made triple-ζ Gaussian atomic orbital basis set as well as with refitted D4 and geometrical counter-poise corrections for London-dispersion and basis set superposition error. The performance of the new method is evaluated for the GMTKN55 database covering large parts of chemical space with about 1500 data points, as well as additional benchmarks for non-covalent interactions, organometallic reactions, and lattice energies of organic molecules and ices, as well as for the adsorption on polar salt and non-polar coinage-metal surfaces. These comprehensive tests reveal a spectacular performance and robustness of r²SCAN-3c: It by far surpasses its predecessor B97-3c at only twice the cost and provides one of the best results of all semi-local density-functional theory (DFT)/QZ methods ever tested for the GMTKN55 database at one-tenth of the cost. Specifically, for reaction and conformational energies as well as non-covalent interactions, it outperforms prominent hybrid-DFT/QZ approaches at two to three orders of magnitude lower cost. Perhaps, the most relevant remaining issue of r²SCAN-3c is self-interaction error (SIE), owing to its mGGA nature. However, SIE is slightly reduced compared to other (m)GGAs, as is demonstrated in two examples. After all, this remarkably efficient and robust method is chosen as our new group default, replacing previous composite DFT and partially even expensive high-level methods in most standard applications for systems with up to several hundreds of atoms.

Data-efficient machine learning for molecular crystal structure prediction

The combination of modern machine learning (ML) approaches with high-quality data from quantum mechanical (QM) calculations can yield models with an unrivalled accuracy/cost ratio. However, such methods are ultimately limited by the computational effort required to produce the reference data. In particular, reference calculations for periodic systems with many atoms can become prohibitively expensive for higher levels of theory. This trade-off is critical in the context of organic crystal structure prediction (CSP). Here, a data-efficient ML approach would be highly desirable, since screening a huge space of possible polymorphs in a narrow energy range requires the assessment of a large number of trial structures with high accuracy. In this contribution, we present tailored Δ-ML models that allow screening a wide range of crystal candidates while adequately describing the subtle interplay between intermolecular interactions such as H-bonding and many-body dispersion effects. This is achieved by enhancing a physics-based description of long-range interactions at the density functional tight binding (DFTB) level-for which an efficient implementation is available-with a short-range ML model trained on high-quality first-principles reference data. The presented workflow is broadly applicable to different molecular materials, without the need for a single periodic calculation at the reference level of theory. We show that this even allows the use of wavefunction methods in CSP.

Biomolecular QM/MM Simulations: What Are Some of the "Burning Issues"?

QM/MM simulations have become an indispensable tool in many chemical and biochemical investigations. Considering the tremendous degree of success, including recognition by a 2013 Nobel Prize in Chemistry, are there still "burning challenges" in QM/MM methods, especially for biomolecular systems? In this short Perspective, we discuss several issues that we believe greatly impact the robustness and quantitative applicability of QM/MM simulations to many, if not all, biomolecules. We highlight these issues with observations and relevant advances from recent studies in our group and others in the field. Despite such limited scope, we hope the discussions are of general interest and will stimulate additional developments that help push the field forward in meaningful directions.

Identification of synthesisable crystalline phases of water – a prototype for the challenges of computational materials design

We discuss the identification of experimentally realisable crystalline phases of water to outline and contextualise some of the diverse building blocks of a computational materials design process.

Perspectives on van der Waals Density Functionals: The Case of TiS2

The van der Waals interaction is of foundational importance for a wide variety of physical systems. In particular, van der Waals forces lie at the heart of potential device technologies that may be realized from the functional organization of layered two-dimensional (2D) nanomaterials. For intermediate to large-scale applications modeling, van der Waals density functionals have become the de facto choice for first-principles calculations. In particular, the vdW-DF family of functionals have provided a systematic approach to this theoretically challenging problem. While much progress has been made, there remains room for improvement in the microscopic description of vdW forces from these density functionals. In this work, we compute benchmark results for the binding energy and the electronic density response to binding in TiS₂ via accurate diffusion quantum Monte Carlo calculations. We compare these benchmark data to results obtained from local, semilocal, and van der Waals functionals. In particular, we gauge the quality of the original vdW-DF/vdW-DF2 functionals, as well as updated variants such as vdW-DF-C09, vdW-DF-optB88, vdW-DF-optB86b, and vdW-DF2-B86R. We find a close relationship between the accuracy of predicted interlayer separation distances and binding energies for TiS₂, with the vdW-DF-optB88 functional performing very well in terms of both quantities. In general, the more recently developed functionals are systematic improvements over older ones. However, when considering the response of the electron density to binding, we find that local-density approximation (LDA) and PBEsol generally outperform the vdW-DF functionals in describing the interlayer charge accumulation with vdW-DF-C09 variants performing the best overall.

Hydration of NH_{4}^{+} in Water: Bifurcated Hydrogen Bonding Structures and Fast Rotational Dynamics

Understanding the hydration and diffusion of ions in water at the molecular level is a topic of widespread importance. The ammonium ion (NH_{4}^{+}) is an exemplar system that has received attention for decades because of its complex hydration structure and relevance in industry. Here we report a study of the hydration and the rotational diffusion of NH_{4}^{+} in water using ab initio molecular dynamics simulations and quantum Monte Carlo calculations. We find that the hydration structure of NH_{4}^{+} features bifurcated hydrogen bonds, which leads to a rotational mechanism involving the simultaneous switching of a pair of bifurcated hydrogen bonds. The proposed hydration structure and rotational mechanism are supported by existing experimental measurements, and they also help to rationalize the measured fast rotation of NH_{4}^{+} in water. This study highlights how subtle changes in the electronic structure of hydrogen bonds impacts the hydration structure, which consequently affects the dynamics of ions and molecules in hydrogen bonded systems.

KITE: high-performance accurate modelling of electronic structure and response functions of large molecules, disordered crystals and heterostructures

We present KITE, a general purpose open-source tight-binding software for accurate real-space simulations of electronic structure and quantum transport properties of large-scale molecular and condensed systems with tens of billions of atomic orbitals (N ∼ 1010). KITE's core is written in C++, with a versatile Python-based interface, and is fully optimized for shared memory multi-node CPU architectures, thus scalable, efficient and fast. At the core of KITE is a seamless spectral expansion of lattice Green's functions, which enables large-scale calculations of generic target functions with uniform convergence and fine control over energy resolution. Several functionalities are demonstrated, ranging from simulations of local density of states and photo-emission spectroscopy of disordered materials to large-scale computations of optical conductivity tensors and real-space wave-packet propagation in the presence of magneto-static fields and spin-orbit coupling. On-the-fly calculations of real-space Green's functions are carried out with an efficient domain decomposition technique, allowing KITE to achieve nearly ideal linear scaling in its multi-threading performance. Crystalline defects and disorder, including vacancies, adsorbates and charged impurity centres, can be easily set up with KITE's intuitive interface, paving the way to user-friendly large-scale quantum simulations of equilibrium and non-equilibrium properties of molecules, disordered crystals and heterostructures subject to a variety of perturbations and external conditions.

Extension and evaluation of the D4 London-dispersion model for periodic systems

We present an extension of the DFT-D4 model [J. Chem. Phys., 2019, 150, 154122] for periodic systems. The main new ingredients are additional reference polarizabilities for highly-coordinated group 1-5 elements derived from pseudo-periodic electrostatically-embedded cluster calculations. To illustrate the performance of the updated method, several test cases are considered, for which we compare D4 to its predecessor D3(BJ), as well as to a comprehensive set of other dispersion-corrected methods. The largest improvements are observed for solid-state polarizabilities of 16 inorganic salts, where the D4 model achieves an unprecedented accuracy, surpassing its predecessor as well as other, computationally much more demanding approaches. For cell volumes and lattice energies of two sets of chemically diverse molecular crystals, the accuracy gain is less pronounced compared to the already excellently performing D3(BJ) method. For the challenging adsorption energies of small organic molecules on metallic as well as on ionic surfaces, DFT-D4 provides values in good agreement with experimental and/or high-level references. These results suggest the application of the proposed D4 model as a physically improved yet computationally efficient dispersion correction for standard DFT calculations as well as low-cost approaches like semi-empirical or even force-field models.

Efficient and accurate description of adsorption in zeolites

Accurate theoretical methods are needed to correctly describe adsorption on solid surfaces or in porous materials. The random phase approximation (RPA) with singles corrections scheme and the second order Møller-Plesset perturbation theory (MP2) are two schemes, which offer high accuracy at affordable computational cost. However, there is little knowledge about their applicability and reliability for different adsorbates and surfaces. Here, we calculate adsorption energies of seven different molecules in zeolite chabazite to show that RPA with singles corrections is superior to MP2, not only in terms of accuracy but also in terms of computer time. Therefore, RPA with singles is a suitable scheme for obtaining highly accurate adsorption energies in porous materials and similar systems.

Excited States with Selected Configuration Interaction-Quantum Monte Carlo: Chemically Accurate Excitation Energies and Geometries

We employ quantum Monte Carlo to obtain chemically accurate vertical and adiabatic excitation energies, and equilibrium excited-state structures for the small, yet challenging, formaldehyde and thioformaldehyde molecules. A key ingredient is a robust protocol to obtain balanced ground- and excited-state Jastrow-Slater wave functions at a given geometry, and to maintain such a balanced description as we relax the structure in the excited state. We use determinantal components generated via a selected configuration interaction scheme which targets the same second-order perturbation energy correction for all states of interest at different geometries, and fully optimize all variational parameters in the resultant Jastrow-Slater wave functions. Importantly, the excitation energies as well as the structural parameters in the ground and excited states are converged with very compact wave functions comprising few thousand determinants in a minimally augmented double-ζ basis set. These results are obtained already at the variational Monte Carlo level, the more accurate diffusion Monte Carlo method yielding only a small improvement in the adiabatic excitation energies. We find that matching Jastrow-Slater wave functions with similar variances can yield excitation energies compatible with our best estimates; however, the variance-matching procedure requires somewhat larger determinantal expansions to achieve the same accuracy, and it is less straightforward to adapt during structural optimization in the excited state.

Ion-Hydroxyl Interactions: From High-Level Quantum Benchmarks to Transferable Polarizable Force Fields

Ion descriptors in molecular mechanics models are calibrated against reference data on ion-water interactions. It is then typically assumed that these descriptors will also satisfactorily describe interactions of ions with other functional groups, such as those present in biomolecules. However, several studies now demonstrate that this transferability assumption produces, in many different cases, large errors. Here we address this issue in a representative polarizable model and focus on transferability of cationic interactions from water to a series of alcohols. Both water and alcohols use hydroxyls for ion-coordination, and, therefore, this set of molecules constitutes the simplest possible case of transferability. We obtain gas phase reference data systematically from "gold-standard" quantum Monte Carlo and CCSD(T) methods, followed by benchmarked vdW-corrected DFT. We learn that the original polarizable model yields large gas phase water → alcohol transferability errors - the RMS and maximum errors are 2.3 and 5.1 kcal/mol, respectively. These errors are, nevertheless, systematic in that ion-alcohol interactions are overstabilized, and systematic errors typically imply that some essential physics is either missing or misrepresented. A comprehensive analysis shows that when both low- and high-field responses of ligand dipole polarization are described accurately, then transferability improves significantly - the RMS and maximum errors in the gas phase reduce, respectively, to 0.9 and 2.5 kcal/mol. Additionally, predictions of condensed phase transfer free energies also improve. Nevertheless, within the limits of the extrathermodynamic assumptions necessary to separate experimental estimates of salt dissolution into constituent cationic and anionic contributions, we note that the error in the condensed phase is systematic, which we attribute, at least, partially to the parametrization in long-range electrostatics. Overall, this work demonstrates a rational approach to boosting transferability of ionic interactions that will be applicable broadly to improving other polarizable and nonpolarizable models.

Physisorption of Water on Graphene: Subchemical Accuracy from Many-Body Electronic Structure Methods

Wet carbon interfaces are ubiquitous in the natural world and exhibit anomalous properties, which could be exploited by emerging technologies. However, progress is limited by lack of understanding at the molecular level. Remarkably, even for the most fundamental system (a single water molecule interacting with graphene), there is no consensus on the nature of the interaction. We tackle this by performing an extensive set of complementary state-of-the-art computer simulations on some of the world's largest supercomputers. From this effort a consensus on the water-graphene interaction strength has been obtained. Our results have significant impact for the physical understanding, as they indicate that the interaction is weaker than predicted previously. They also pave the way for more accurate and reliable studies of liquid water at carbon interfaces.

Cost-Effective Quantum Mechanical Approach for Predicting Thermodynamic and Mechanical Stability of Pure-Silica Zeolites

Several computational techniques for solid-state applications have recently been proposed to enlarge the scope of computer simulations of large molecular systems. In this contribution, we focused on two of these, namely, HF-3c and PBEh-3c. They were recently proposed by the Grimme's group, as "low-cost" ab initio-based techniques for electronic structure calculation of large systems and were proved to be effective essentially for organic molecules. HF-3c is based on a Hartree-Fock Hamiltonian with a minimal Gaussian quality basis set, whereas PBEh-3c is a density functional theory (DFT) based method with a hybrid functional and a medium-quality basis set. Both HF-3c and PBEh-3c account for dispersion (London) interactions and are free from the basis set superposition error due to limited basis set size, through several pairwise semiempirical corrections. To the best of our knowledge, despite the promising results on the cost-accuracy side of molecular simulations of organic molecules, these methods have been used only in few cases for solid-state applications. In this contribution, we studied the performance of HF-3c and PBEh-3c for predicting the properties of inorganic crystals to enlarge the applicability of these cheap and fast methodologies. As a testing ground, we have chosen a well-known class of material, e.g., microporous all-silica zeolites. We benchmarked geometries, formation energies, vibrational features, and mechanical properties by comparing the results with literature data from both experiment and computer simulation. For structures, HF-3c is extremely accurate in predicting the zeolites cell volume, albeit we do not include any vibrational contribution, neither zero point nor thermal, on the computed volumes, which may introduce small variations in the predicted values. For the energetic, the relative stability of the zeolites using the DFT//HF-3c approach allows predictions within the experimental error for most of the cases taken into consideration when the experimental enthalpies were corrected back to electronic energies by using the HF-3c thermodynamic contributions computed in the harmonic approximation. This strategy is particularly convenient, as the slow step (geometry optimization) is carried out with the cheapest HF-3c method, whereas the fast step (single point energy evaluation) is carried out with costly DFT methods. In this sense, the use of the DFT//HF-3c approach results to be a promising one to predict the stability and structure of microporous materials. Finally, the HF-3c method predicts the mechanical properties of the zeolite set in reasonable agreement with respect to those computed with the state-of-the-art DFT simulations, indicating the HF-3c method as a possible technique for the mechanical stability screenings of microporous materials.

Understanding non-covalent interactions in larger molecular complexes from first principles

Non-covalent interactions pervade all matter and play a fundamental role in layered materials, biological systems, and large molecular complexes. Despite this, our accumulated understanding of non-covalent interactions to date has been mainly developed in the tens-of-atoms molecular regime. This falls considerably short of the scales at which we would like to understand energy trends, structural properties, and temperature dependencies in materials where non-covalent interactions have an appreciable role. However, as more reference information is obtained beyond moderately sized molecular systems, our understanding is improving and we stand to gain pertinent insights by tackling more complex systems, such as supramolecular complexes, molecular crystals, and other soft materials. In addition, accurate reference information is needed to provide the drive for extending the predictive power of more efficient workhorse methods, such as density functional approximations that also approximate van der Waals dispersion interactions. In this perspective, we discuss the first-principles approaches that have been used to obtain reference interaction energies for beyond modestly sized molecular complexes. The methods include quantum Monte Carlo, symmetry-adapted perturbation theory, non-canonical coupled cluster theory, and approaches based on the random-phase approximation. By considering the approximations that underpin each method, the most accurate theoretical references for supramolecular complexes and molecular crystals to date are ascertained. With these, we also assess a handful of widely used exchange-correlation functionals in density functional theory. The discussion culminates in a framework for putting into perspective the accuracy of high-level wavefunction-based methods and identifying future challenges.

Theory and practice of modeling van der Waals interactions in electronic-structure calculations

The accurate description of long-range electron correlation, most prominently including van der Waals (vdW) dispersion interactions, represents a particularly challenging task in the modeling of molecules and materials. vdW forces arise from the interaction of quantum-mechanical fluctuations in the electronic charge density. Within (semi-)local density functional approximations or Hartree-Fock theory such interactions are neglected altogether. Non-covalent vdW interactions, however, are ubiquitous in nature and play a key role for the understanding and accurate description of the stability, dynamics, structure, and response properties in a plethora of systems. During the last decade, many promising methods have been developed for modeling vdW interactions in electronic-structure calculations. These methods include vdW-inclusive Density Functional Theory and correlated post-Hartree-Fock approaches. Here, we focus on the methods within the framework of Density Functional Theory, including non-local van der Waals density functionals, interatomic dispersion models within many-body and pairwise formulation, and random phase approximation-based approaches. This review aims to guide the reader through the theoretical foundations of these methods in a tutorial-style manner and, in particular, highlight practical aspects such as the applicability and the advantages and shortcomings of current vdW-inclusive approaches. In addition, we give an overview of complementary experimental approaches, and discuss tools for the qualitative understanding of non-covalent interactions as well as energy decomposition techniques. Besides representing a reference for the current state-of-the-art, this work is thus also designed as a concise and detailed introduction to vdW-inclusive electronic structure calculations for a general and broad audience.

Accurate Predictions of Electron Binding Energies of Dipole-Bound Anions via Quantum Monte Carlo Methods

Neutral molecules with sufficiently large dipole moments can bind electrons in diffuse nonvalence orbitals with most of their charge density far from the nuclei, forming so-called dipole-bound anions. Because long-range correlation effects play an important role in the binding of an excess electron and overall binding energies are often only on the order of 10s-100s of wave numbers, predictively modeling dipole-bound anions remains a challenge. Here, we demonstrate that quantum Monte Carlo methods can accurately characterize molecular dipole-bound anions with near-threshold dipole moments. We also show that correlated sampling Auxiliary Field Quantum Monte Carlo is particularly well-suited for resolving the fine energy differences between the neutral and anionic species. These results shed light on the fundamental limitations of quantum Monte Carlo methods and pave the way toward using them for the study of weakly bound species that are too large to model using traditional electron structure methods.

Simplified DFT methods for consistent structures and energies of large systems

Kohn-Sham density functional theory (DFT) is routinely used for the fast electronic structure computation of large systems and will most likely continue to be the method of choice for the generation of reliable geometries in the foreseeable future. Here, we present a hierarchy of simplified DFT methods designed for consistent structures and non-covalent interactions of large systems with particular focus on molecular crystals. The covered methods are a minimal basis set Hartree-Fock (HF-3c), a small basis set screened exchange hybrid functional (HSE-3c), and a generalized gradient approximated functional evaluated in a medium-sized basis set (B97-3c), all augmented with semi-classical correction potentials. We give an overview on the methods design, a comprehensive evaluation on established benchmark sets for geometries and lattice energies of molecular crystals, and highlight some realistic applications on large organic crystals with several hundreds of atoms in the primitive unit cell.

QMCPACK: an open source ab initio quantum Monte Carlo package for the electronic structure of atoms, molecules and solids

QMCPACK is an open source quantum Monte Carlo package for ab initio electronic structure calculations. It supports calculations of metallic and insulating solids, molecules, atoms, and some model Hamiltonians. Implemented real space quantum Monte Carlo algorithms include variational, diffusion, and reptation Monte Carlo. QMCPACK uses Slater-Jastrow type trial wavefunctions in conjunction with a sophisticated optimizer capable of optimizing tens of thousands of parameters. The orbital space auxiliary-field quantum Monte Carlo method is also implemented, enabling cross validation between different highly accurate methods. The code is specifically optimized for calculations with large numbers of electrons on the latest high performance computing architectures, including multicore central processing unit and graphical processing unit systems. We detail the program's capabilities, outline its structure, and give examples of its use in current research calculations. The package is available at http://qmcpack.org.