Future directions of chemical theory and computation

Density functional theory

Density functional theory (DFT) is a quantum mechanical method that studies the electronic structure of multi-electron systems. It has been widely used in chemistry, physics, and materials science [21], [22], [23]. In the framework of the Kohn-Sham formulation of DFT, the most difficult many-body problem due to electron interactions is reduced to a problem of independent electrons moving in an effective potential field [24]. This effective potential field includes the effects of the external potential and Coulombic interactions between electrons including exchange and correlation interactions.

Since 1970, DFT has been widely used in solid-state physics calculations. Density functional theory using local density approximation (LDA) with Slater exchange and VWN correlation from a uniform electron gas usually leads to very satisfactory results with very low computational cost especially when compared with other methods dealing with the many-body problems in quantum mechanics [25]. However, LDA offered less satisfactory results for a molecule than for solids due to the inhomogeneity in electron density. The Parr-Yang book is a milestone in DFT with faster development of better exchange correlation functions suitable for chemical problems [23]. The development of an exchange-correlation function began with LDA and then moved to the general gradient approximation (GGA) such as PBE and PW91. The field progressed with meta-GGA in consideration of orbital kinetic energy and exact exchange such as TPSS and M06-L. Hybrid GGA is now the most widely employed and is represented by B3LYP, hybrid meta-GGA such as TPSSh and M06-2x, and double hybrid-GGA containing unoccupied orbital contributions such as XYG3 and B2PLYP. Savin first proposed to separate long and short range Coulombic interactions [26] that were long ignored by the community but later found to be important in describing the charge transfer excitation and response properties [27]. Range separation for the functional with an optimal ω (the range separation parameter) is now seen in daily practice for modeling the molecular energy materials (often with charge transfer excitations).

One of the main challenges of DFT is to maintain its simplicity while improving its accuracy or expanding its functionals. This is also a problem faced by the entire field of computer science. It is important to keep the theory and calculation concise and to ensure the computational costs are within an acceptable range because this leads to proper density functional approximations. The great success and popularity of DFT lie in that some simple approximations work very well for the calculations of the structural and thermodynamic properties of molecules and solids [28], [29], [30]. However, there are still many qualitative failures of DFT which result from the inappropriate approximations of the exchange-correlation functionals. The core issue that how to construct the universally applicable functionals has been elusive, which is full of challenges.

To describe the chemical phenomena and processes more comprehensively, it is necessary to jump out of the equilibrium geometries of molecules and consider the weak interactions between molecules and the transition states in chemical reactions. It remains a great challenge to describe van der Waals forces, covalent bonds, and reaction barriers correctly and efficiently.

It is known that all semilocal density functionals fail to describe the long-range dispersion interactions [31]. Various approaches were proposed to include dispersion corrections into DFT, such as DFT-D-type methods [32], [33], [34], [35], [36], vdW-DFs [37], [38], 1ePOT [39], [40], [41], and so on. Among them, the DFT-D method proposed by Grimme (DFT-D1 [32], DFT-D2 [33], DFT-D3 [35], and DFT-D4 [36]) is the most widely-used one. The DFT-D method solves the dispersion problem of DFT in a rather general way and the correction can be coupled in a simple form to standard density functionals. The DFT-D3 has been refined with higher accuracy, applies to a broader range of applicability, and adopts less empiricism [35]. It is just now replacing DFT-D2 as the worldwide de facto standard in dispersion corrected DFT calculations. The DFT-D method opens completely new possibilities for the application of DFT in the areas of condensed matter, materials science, and biochemistry where dispersion effects are often of utmost importance [32], [33], [34], [35], [36].

The empirical elements are embedded in almost all of the current dispersion corrections in different ways. Hence, it is essential to perform a systemical benchmark on experimental or reliable theoretical data [42]. For the intramolecular and intramolecular interactions, the gold standard is the WF-based singles and doubles coupled-cluster method with perturbative triples (CCSD(T)). The method can provide accurate results, whose error is less than 1 kcal/mol for typical chemical reactions [43], [44]. The domain-based local pair natural orbital coupled-cluster method, abbreviated as DLPNO-CCSD(T), is an approximation to CCSD(T) and it has the following advantages: (1) Accurate. DLPNO-CCSD(T) recovers more than 99.9 % of the CCSD(T) correlation energy. Reaction energies are calculated with a mean deviation of 0.3 kcal/mol for 12 test reactions of medium-sized molecules [45]. (2) The computational cost of DLPNO-CCSD(T) is comparable to DFT but scales linearly with the system size [45], [46]. (3) DLPNO-CCSD(T) operates like a black box without adjusting any complicated parameters.

The majority of the calculation failures of DFT result from its underestimation of the reaction barrier, dissociation energy of ionic molecules, energy gaps, excited states, intermolecular interactions, and charge transfer excitation. Meanwhile, DFT overestimates the binding energies of charge transfer complexes and the response to an electric field in molecules and materials [47]. Both the underestimation and overestimation result from the same origin, i.e., the delocalization error of the approximated functionals. Due to the discrete nature of electrons, the exact energy of the atom as a function of the charge should be a straight-line interpolation between the integers [48]. Yang et al. found that there were delocalization errors when the fractional charge deviates from a linear relationship. They further found that the deviation of fractional spins is responsible for static correlation-related errors, e.g., a strong correlation effect in oxide and magnetic systems, degenerate states, and bond-breaking [49].

Langreth-Perdew’s adiabatic connection-fluctuation-dissipation (ACFD) theorem is derived from density fluctuations and has become an important foundation for systematic improvements in functionals from many-body theory, e.g., connecting with the random-phase approximation (RPA-DFT). Parallel to density fluctuation-based ACFD, Yang et al. proposed another adiabatic connection based on a pairing matrix fluctuation that led to the particle-particle RPA formulation of DFT and ppRPA-DFT. The ppRPA-DFT demonstrates the nearly linear fractional charge behavior and thus the delocalization errors are minimized: the ppRPA-DFT correctly describes the van der Waals interaction and has almost no static correlation error for single-bond systems [50]. Both the occupied and unoccupied Kohn-Sham orbitals should be considered. Thus, learning from wavefunction theory in dealing with the correlation problem is critical to developing radically different functionals to advance DFT.

Post-Hartree-Fock and beyond

The wavefunction based methods are increasingly popular because of their systematic accuracy improvements and the development of efficient low-scaling computational techniques. Wavefunctions also have value in fields requiring high accuracy such as chemical reaction dynamics. Wavefunction theory starts with the self-consistent field Hartree-Fock (HF) mean-field theory, which is the origin of the modern ab initio quantum chemistry. Later, many post-Hartree-Fock methods have been proposed to treat the correlation effects including configuration interaction (CI), perturbation (MPn), coupled-cluster (CC), multiconfiguration self-consistent field (MCSCF), complete active space self-consistent field (CASSCF), and CAS perturbation. Multireference CC remains very time-consuming.

The explicit correlation methods (R12 and F12) can incorporate the interelectron distance r₁₂ into the wavefunction ansatz. This idea was originally developed by Hylleraas already in 1929 for the helium atom. It can correctly describe the cusp behavior of the wavefunction [51]. R12 was first introduced in MP2 and then made popular in coupled cluster theory. Later, a Slater-type correlation F12 form was postulated, which could be fitted through a linear combination of a few Gaussian-type functions. It is still a challenge to formulate an analytical gradient with an explicit correlation form.

The introduction of quantum Monte Carlo (QMC) into traditional wavefunction theory by the Prof. Ali Alavi of Cambridge into perturbation theory was a major step forward. There was later full configuration interaction QMC (FCIQMC) [52]. The imaginary time Schroedinger equation becomes a diffusion equation in QMC. The required antisymmetric property is violated because the lowest energy solution is generally nodeless and symmetric. This has long been a significant problem for QMC. Alavi suggested performing a long-time integration in the space of a Slater determinant; there was a propagation step via population dynamics. The “walker” in the simulation represents an instantaneous wavefunction that carries a sign and a pair of walkers coinciding with the same determinant but with different signs; these are removed from the simulation. Such a new QMC algorithm in the Slater determinant space can efficiently converge to a full CI energy. This approach was quickly applied to coupled cluster even with periodic conditions for solids with great success [53]. The results offered unprecedented accuracy.

The density matrix renormalization group (DMRG) method was initially proposed in condensed matter physics for strongly correlated models [54] and has been quickly applied to quantum chemistry [55], [56]. DMRG can now handle a much larger active space than conventional CASSCF [57]. DMRG iteratively optimizes the eigenvectors of the reduced density matrix, which are used as approximate natural orbitals for constructing a full CI-like but much lower dimension Hilbert space. Recently, DMRG ansatz is understood as a matrix product state (MPS) formalism from applied mathematics and quantum information theory. The optimization procedures are actually consecutive singular value decompositions and maintain a few important singular values in each step. The International Academy of Quantum Molecular Science has twice awarded the annual (single) medal for this field. In 2010, Garnet KL Chan was awarded the IAQMS medal with citation “for his outstanding contribution to the density matrix renormalization group theory of molecular system”, and Takeshi Yanai in 2013 “for his development of novel approaches to incorporate dynamical correlation into DMRG using canonical transformation theory”.

The most recent remarkable advances are made for the time-dependent MPS formalism for finite-temperature problems, which is quickly emerging as a promising method to deal with quantum dynamics for complex systems [58]. In fact, the MPS ansatz along with the more recent tensor tree network state can grasp the essential quantum entanglements, thus offering proper spaces for time evolution. The applications are envisioned for non-adiabatic dynamics and electronic processes in complex systems such as carrier transport, various optical spectroscopies, singlet fission, and organic light-emitting diodes.

Molecular design

Molecular design is a hot research topic especially for materials design and drug design. The design and optimization of these two kinds of molecules can bring great benefits and have great potential value. Therefore, related studies have dramatically stimulated the application of deep learning in related fields. However, the process of molecular design is still very challenging because it is expensive and time-consuming considering [49], [59], [60].

Most of today’s technologies such as batteries, aerospace, and renewable energy strongly rely on the synthesis and application of advanced materials. Artificial intelligence has only recently begun to affect the field of materials design, and the potential impact of such data-driven materials science is tremendous [61], [62], [63]. Computer-aided materials design can significantly reduce the typical cycle for the development and commercialization of new materials [64]. Specifically, the integration of artificial intelligence algorithms could help address the inverse chemical paradigm, automatically discovering the property-biased molecular structures, and thus expediting the design of novel useful compounds [65]. However, the current material design faces the following challenges: (1) obscure terminology in material informatics makes it difficult for typical material scientists to see how data-driven methods are applied to their own work; (2) limited access to structured data and a lack of data standards; (3) absence of a mature data-sharing mechanism [66].

Synthesis planning

Synthesis planning can be divided into three parts: retrosynthesis, reaction prediction, and reaction optimization. In retrosynthesis, the product is known and the initial reagents from which it can be made are determined. In reaction prediction, the reagents are known, and the main products are predicted. In reaction optimization, both reagents and products are known, and the yield or efficiency of this chemical reaction can be maximized by changing the reaction conditions [67], [68].

The application of machine learning and artificial intelligence technologies in synthesis planning has unprecedented opportunities. Since the machine learning models could discover the hidden relationship between the dataset and the experimental results, efficiently and relatively low-priced, it would promote the development of retrosynthesis and replace the expensive quantum calculations. However, there are still some urgent problems to be solved: (1) Detailed reaction data are limited. (2) The literature prefers to report only successful reactions rather than failed ones, and access to negative reaction data is necessary for more efficient synthesis planning [67].

Challenges in theory and simulations of polymer science

The challenges faced by polymer theory and simulation include the following aspects: The need for design rules for a specific functional polymeric material including hierarchical composites and biomaterials; understanding the transport process of electrons, ions, photons, and energy in polymer systems (solar cells, living cells, batteries, sensors, and membranes); and exploring more efficient and stable processing strategy for polymer materials with complex structures [69], [70].

Polymer processing often involves a strong flow field, temperature gradients, and other external conditions that make polymer systems far away from the equilibrium state. It is a great challenge to predict the properties of polymer materials at the micro, meso, and macro time or length scales. It is impractical to develop a theory considering all variables, and thus it is very important to make proper approximations and reasonable coarse graining.

Electrically charged polymeric systems

Charged polymers involve a series of synthetic materials and consist of the basic unit of most functional biomaterials. The polyelectrolytes in polymer solutions are strongly correlated in terms of both topologies and electrostatic interactions. Moreover, there are lots of factors influencing the polyelectrolyte behaviors, such as the charge distribution on the polyelectrolytes, the electrostatic correlations, size effects, and the polarizability effects around charged polymers [71]. Despite the theoretical efforts devoted to the scaling theories, the analytical theories basing on the field theory, and the liquid-state theories, there are still two major challenges that need to be done. The first challenge associated with electrically charged polymeric systems is that the dielectric properties are spatially heterogeneous, which means that the mean-field theory fails to describe such a system. The second challenge is the theoretical description of molecular conformations or structure factors of charged polymers in saline/salt-free solutions: the correlation length of two charged units is large enough to even match the size of the molecular conformation, thus coupling with the molecular conformations on larger length scales [69], [71].

Crystallization and glass transition

Although substantial progress has been made in these classic problems in polymer physics, there are still significant challenges in understanding the following situations. (1) The interplay between conformational entropy (favoring disorder) and the short-range attraction between segments (favoring order) results in extremely complex behaviors of polymer crystallization and glass transition. Such processes are far away from the thermodynamic equilibrium and usually take place in supercooled metastable states. The current available theoretical models succeed in describing the polymer crystallization and glass transition processes, but fail to predict these behaviors. Therefore, a universal theoretical framework is urgently needed to describe the kinetics of crystallization and glass transition [72], [73]. (2) It is still far from understood about the origin of the temperature and pressure dependences of the equilibrium dynamics and the corresponding quantitative description of the behaviors remains to be explored. Meanwhile, the origin of the Kauzmann paradox and the relationship between the evolution of the glass structure and physical aging are still unanswered [74]. (3) How the chain connectivity and uncrossability affect the mechanical properties of the entangled polymer system? (4) How to theoretically explain the special nonequilibrium state generated by the thermal history also needs to be further understood, such as the brittle-ductile transition and glass transition [73], [74].

Nonlinear rheological behaviors of entangled polymer fluids

Nonlinear rheology is to study the mechanical behavior of materials under a large strain or rapid deformation conditions. The stress response of materials depends on the deformation amplitude, rate, and flow field type. Owing to the failure of the Boltzmann superposition principle, the nonlinear behavior cannot be predicted by linear viscoelastic behavior. The development of rheological measurement technologies and structural characterization methods has led to the development of nonlinear rheological experimental research. The field moved from focusing on the study of steady-state shear properties to focusing on transient/dynamic shear, transient tension, and even more complex flow field rheological response. The description of nonlinear rheological behavior has also developed from empirical continuum mechanics model to molecular model. However, there is still no universal constitutive model to describe all nonlinear behaviors. Therefore, a key goal of nonlinear rheology is to understand the structural evolution mechanism and dynamic behavior under the condition of complex flow, and to establish a relationship with rheological properties [75], [76], [77], [78], [79].

Accelerated computational models

Utilizing physics-based calculations, the computational models of chemistry are adopted to investigate various properties and behaviors of the modeling system [68]. There are two different ways to apply deep learning in these computational models. The first one is to integrate deep learning methods with physics-based approaches directly, which involves network training to predict the key property. The second method is to predict properties by establishing the relationship between the molecular structure and particular properties. The corresponding examples of the integrated physics-based approach include the prediction of force fields [84], potential energy surfaces [85], [86], [87], and corrections of ab initio calculations [88]. Such method is more flexible due to the physical basis, however, it has a slower speed when compared with the second method.

Quantitative structure property/activity relationships

In computational chemistry, there is an alternative approach to realize the deep learning by directly mapping a simple representation of the considered molecules to the target property [68]. Such an approach can be broadly divided into THE quantitative structure property relationship (QSPR) and quantitative structure activity relationship (QSAR). Generally speaking, QSPR predicts the properties of molecular systems. However, QSAR focuses on molecule activity. Both methods are employed to improve the prediction accuracy [89], [90]. The training data available directly determines the properties which can be predicted. Therefore, the choice of the available databases is essential, which are summarized in the excellent review of Butler. It should be noted that some properties can be readily computed with DFT, such as the dipole moments, ionization energies, ground state energies, and so on. The inclusion of computational data sets greatly maximizes the speed and enriches the available data as much as possible.

Challenges in large-scale molecular dynamics

Progress in theoretical and computational chemistry originates from technical research and developments in the field of materials science. Related research can directly solve practical problems leading to the so-called social demand-driven research. Owing to the understanding of fundamental problems, researchers inevitably have to solve these problems when they try to explain the mechanism of experimental phenomena and processes, i.e., problem-driven research. Although the specific problems and systems studied change with time, the current hot research topics will be gradually replaced with new problems. Researchers always have a strong interest in performing larger and larger simulations processing more complex systems, simulating longer dynamics, analyzing more microscopic details, and increasing the diversity of the studied systems [66], [91], [92]. However, large-scale molecular dynamics simulations require access to multicore clusters or supercomputers that are not always available to all researchers and more and more researchers explore the potential applications of GPUs consequently.

Applications in reaction prediction and optimization

The development of technologies for chemical synthesis traditionally required expert knowledge and special laboratory practice. Artificial intelligence might solve these problems to create new knowledge more effectively [15]. By putting a huge amount of experimental and related theoretical data into the computer program, deep learning and other artificial intelligence technologies would help experimental chemists find new important reaction routes [17]. Automation of retrosynthetic analysis has already been put into practice and different learning methods for the prediction of reaction products are also being developed. Using the same logic, one could also optimize reaction conditions and design natural product drugs. Further application of artificial intelligence technology in those fields must accelerate related research and industry [93]. Although reinforcement learning and machine learning has been used to optimize reaction conditions, access to high-quality, interpretable, and standardized data sets suitable for machine learning is a current bottleneck as the literature on chemical reactions is often unstructured, exists in multiple formats, sometimes behind paywalls, and was collected on different reaction setups [94], [95], [96], [97].