Real-time Quantum Chemistry
Significant progress in the development of efficient and fast algorithms for quantum chemical calculations has been made in the past two decades. The main focus has always been the desire to be able to treat ever larger molecules or molecular assemblies—especially linear and sub-linear scaling techniques are devoted to the accomplishment of this goal. However, as many chemical reactions are rather local, they usually involve only a limited number of atoms so that models of about two hundred (or even less) atoms embedded in a suitable environment are sufficient to study their mechanisms. Thus, the system size does not need to be enlarged, but remains constant for reactions of this type that can be described by less than two hundred atoms. The question then arises how fast one can obtain the quantum chemical results. This question is not directly answered by linear-scaling techniques. In fact, ideas such as haptic quantum chemistry or interactive quantum chemistry require an immediate provision of quantum chemical information which demands the calculation of data in “real time”. In this perspective, we aim at a definition of real-time quantum chemistry, explore its realm and eventually discuss applications in the field of haptic quantum chemistry. For the latter we elaborate whether a direct approach is possible by virtue of real-time quantum chemistry.
Submitted for publication as a perspective article in Int. J. Quantum Chem.
The large number of possible nuclear configurations puts the curse of dimensionality on studies of the chemical reactivity of large molecular systems like metalloenzymes or transition metal complexes in homogeneous catalysis. To find exactly those configurations which correspond to the reaction path searched for is very difficult. Although there exist several methods to sample the configuration space of the nuclear positions efficiently and in an unbiased way (see, e.g., Refs. 2, 3, 4, 5), a full ab initio treatment of the molecular system renders also these approaches very time demanding or even unfeasible. An ab initio treatment is, however, mandatory for the study of chemical reactivity involving the forming and breaking of chemical bonds.
The overwhelming amount of data generated by computational methods calls for new approaches to access it. In conventional methods, a change of nuclear coordinates of a reactive or functional molecular system changes the potential energy
i.e., the electronic energy defined in the Born–Oppenheimer (BO) approximation. Once this energy has been calculated by a computer program, the results are collected and can be evaluated. However, if the results were accessible immediately after applying the configurational change, the perception of the results of the calculations would be much more efficient since a real interaction between the scientist and the system under study could be established.
In this perspective article we discuss a point of view on such a challenge that we may summarize under the term Real-time Quantum Chemistry. “Real-time” in this context means that the results of quantum chemical calculations are basically instantaneously available, i.e., on a time scale that a specific human sense would interpret as instantaneous. Quantum chemical information available in real time then offers a very efficient way to interactively study the reactivity of molecular systems.
The structure of this work is as follows. We start with a definition of Real-time Quantum Chemistry and discuss its principles. Then, we review and evaluate existing techniques to accelerate quantum chemical calculations with respect to their potential usage in Real-time Quantum Chemistry. Finally, we discuss the feasibility of a direct approach to Haptic Quantum Chemistry[6, 7] which implements the ideas of Real-time Quantum Chemistry. In order to illustrate what can already be done with currently available quantum chemistry software we present some exploratory calculations.
2 Principles of Real-time Quantum Chemistry
A theoretical study of chemical reactivity requires the exploration of the potential energy surface. Hence, first the wave function of the molecular system has to be calculated, from which then the energy and the gradients can be obtained. Other properties like molecular orbital representations, the electron density, polarizabilities or partial charges could also be considered, but are secondary targets compared to energies and gradients. We thus define
Definition Real-time Quantum Chemistry shall denote the very fast calculation of the quantum mechanical response of a reactive molecular system in terms of the wave function and the corresponding energy and gradient due to a user-driven manipulation of the system’s molecular structure such that the operator experiences the information flow without any time delay.
Since finding the optimized wave function and calculating the energy and the gradient is essential to all reactivity studies, we define their instantaneous calculation as the core objective of Real-time Quantum Chemistry. Obtaining additional properties in real time could then be considered an extension of (core) Real-time Quantum Chemistry.
It is clear that the fast calculation of the core quantities is central also to ab initio molecular dynamics (AIMD). Accordingly, the challenges which one faces in the context of Real-time Quantum Chemistry can also be found in the area of AIMD, where the fast calculation of the nuclear gradients (called ionic forces) is of paramount importance. Therefore, many techniques developed in AIMD to speed up the calculation of gradients are also important in a Real-time Quantum Chemistry framework. Although the forces are required as quickly as possible, AIMD simulation results are interpreted after the simulation has been carried out. In Real-time Quantum Chemistry, however, we want to seamlessly merge the calculation and the perception of the calculated information.
The seamless integration of complex quantum chemical information allows the interactive exploration of chemical reactivity. An approach which utilizes a force-feedback device as an input and an output device to transmit the quantum mechanical information to the operator is Haptic Quantum Chemistry (HQC) introduced by us in 2009[6, 7]. A force-feedback device as in HQC is also utilized in the interactive molecular dynamics framework, which implements a classical treatment of the forces in molecular systems. Another example is the recently presented semi-empirical interactive implementation for optimizing molecular structures of hydrocarbons during editing the structures[10, 11].
As it has been noted in Ref. 11 a real-time visual experience requires at least ten frames per second (i.e. to update the system’s state), if a frame represents a step in the shift of nuclear positions during structure optimization or reaction. In Haptic Quantum Chemistry, however, a smooth tactile experience requires an update rate of 1000 frames per second. Accordingly, the central question for real-time quantum chemical approaches is the following: Is it—in principle—possible to perform sufficiently accurate ab initio quantum chemical calculations in such a short time, i.e., in the millisecond range for a reasonably sized molecular model of a reactive system?
The issue of accelerating quantum chemical calculations has been treated thoroughly before—mostly in the area of linear and sub-linear scaling techniques[12, 13]. However, it is important to note that there is a fundamental difference between the problem formulated in the question above and the problem of achieving linear scaling. In Real-time Quantum Chemistry, we do not ask for methods which allow us to treat larger and larger systems but rather how fast we can perform a calculation for a particular system of constant size and target accuracy. This implies that one also has to focus on the prefactors and onsets of quantum chemical methods and not only on their overall scaling behavior. In other words, the actual execution time for a given system is the prime target.
The reason why we focus on reactive systems with a (large but) constant number of atoms is that we concentrate on rather local chemical events, which are ubiquitous in chemistry. For example, they can be mediated by transition-metal centers as in enzymatic reactions or in homogeneous catalysis. Reactions at transition metals are usually restricted to a limited spatial region containing the metal center(s), its (their) ligand environment, the incoming reactant, co-factors, other reactants that may lead to side reactions and some proper model of the close environment in such a way that models on the order of atoms are sufficient for a meaningful description[14, 15].
3 Methods to Accelerate Electronic Structure Calculations
Almost all acceleration techniques for electronic structure calculations developed so far aim at a linear scaling behavior especially for very large molecules ( atoms). A multitude of methods has been developed as discussed in many excellent reviews[12, 13, 16]. Here, we need to evaluate the existing methods from a different point of view. The overall goal is to decide, whether there is a certain lower limit of computation time for mid-sized molecular systems ( atoms), and to identify approaches that are likely to be important for a real-time calculation of gradients and energies.
In view of the fact that single-determinant models like Hartree–Fock theory and, most importantly, density functional theory (DFT) are likely to be the best candidates for quantitative real-time reactivity exploration, we first need to consider the solution of the Roothaan–Hall equations. The two most important steps in a self-consistent solution of the Roothaan–Hall equations are the construction of the Fock matrix and the subsequent calculation of the density matrix. Of course the size of the basis set chosen is also very important, since it determines the size of the system.
We skip the derivation of the one-particle mean-field Hartree–Fock theory and of Kohn–Sham (KS) density functional theory[18, 19, 20] which lead both to an effective one-particle operator equation. For the sake of simplicity, the following equations are given for spin-restricted calculations. The generalization to the spin-unrestricted case is straightforward. The Fock operator for Hartree–Fock and Kohn–Sham models can be written as[21, 13]
where and are two parameters that define the electronic structure model. For a Hartree–Fock calculation, one chooses and , whereas for a ’pure’ Kohn–Sham calculation and with a non-vanishing hold. With and hybrid approaches are described.
The introduction of a finite basis set ,
is the most convenient way to solve the spatial integro-differential self-consistent-field (SCF) equations that result when is operating on an orbital . This yields the well-known Roothaan–Hall equations,
where is the matrix of the expansion coefficients defined in Eq. (3), is the overlap matrix and is the matrix of the orbital energies of the orbitals . The elements of the Fock matrix in the chosen basis are given by
Here, is the matrix representation of the one-electron operator, is the two-electron Coulomb matrix operator, is the two-electron exchange matrix operator and is the two-electron exchange–correlation matrix operator. For the calculation of the matrices and the two-electron repulsion integrals are contracted with the elements of the density matrix which is for real orbitals in closed-shell systems defined as
and thus the reason for the self-consistent iterative solution of the Roothaan–Hall equation.
The elements of the matrix are derived from the exchange–correlation density functional which can be a functional of the electron density in the local density approximation, of the density and the gradient of the density in the generalized gradient approximation (GGA), or of the density, the gradient of the density, and the kinetic energy density in meta-GGA functionals. Since the matrix elements cannot be evaluated analytically the integration is usually performed on a mesh of grid points constructed from merged atomic grids
where is the number of points in the grid of atom and are the corresponding weights. is the function that splits the molecular grid into atomic sub-grids. The electron density at a grid-point is given by
Eq. (4) has to be solved iteratively, since depends on the elements of . After reaching self consistency the total energy of the system is calculated from the converged density matrix elements and the Fock matrix elements,
where represents the Coulombic pair interaction of all atomic nuclei. From this equation also the gradient with respect to the nuclear coordinates can be derived. Within a finite basis of, e.g., Cartesian Gaussian functions one obtains after a few rearrangements the expression[22, 17]
It includes the Pulay forces that are due to the origin dependence of atom-centered basis functions. In case of plane-wave basis functions these forces vanish and the expression can be simplified. For a more compact notation an energy weighted density matrix ,
has been introduced. The derivatives of the exchange–correlation matrix elements depend on the specific nature of the exchange–correlation potential . General formulations of the energy gradient can be found in Refs. 23, 24, 25, 26 for the local (spin) density approximation and in Ref. 27 for gradient-corrected functionals.
3.1 Basis Sets
Very important for any acceleration technique is the choice of the basis set in which the Kohn–Sham or Hartree–Fock orbitals are expanded. The larger the basis set the more accurate the calculations are, but also the more time they consume. Reducing their size is, therefore, a very seductive means to reduce computation time. In addition certain basis sets speed-up certain parts of the Fock matrix calculation but may have disadvantages in other parts.
The two most widely employed explicit basis sets are linear combinations of Gauss-type atomic orbitals and plane waves. Plane waves are especially suited for the calculation of periodic and homogeneous systems whereas the Gauss-type orbitals are usually employed in molecular systems. Density matrices in Gaussian basis sets are therefore mostly sparse, i.e. band diagonal and thus promote linear-scaling techniques.
Plane-wave basis sets are completely delocalized in the direct space but are very localized in reciprocal space, which is why they are usually applied in solid state physics. In molecular systems, however, a huge number of plane waves would be needed to obtain the same accuracy as with localized Gauss-type orbitals. On the other hand, some of the integrals can be calculated very fast by Fast Fourier Transforms (FFT) within a plane-wave basis set. Also the calculation of the forces on the nuclei is computationally less expensive since the basis functions are not position dependent. To increase the accuracy and limit the number of basis functions in molecular systems, hybrid codes combine Gaussians and plane waves [28, 29, 30]. These developments allow for the calculation of molecules with up to a million atoms.
Pseudopotentials[32, 33] can be employed for both types of basis sets in order to reduce the number basis functions. How many electrons are treated implicitly by the pseudopotential can, of course, be varied and determines the accuracy achieved. In plane-wave calculations pseudopotentials are crucial to avoid the description of the nodal structure of orbitals so that the number of plane-waves can be limited to a reasonably small value. Pseudopotentials are also employed for heavy elements to account for relativistic effects.
A less common alternative to the classes of basis sets described above are the so-called wavelets[35, 36, 37], which aim at combining the best of both worlds. They are localized in both direct and reciprocal space and the integrals can be calculated with very fast methods similar to Fast Fourier Transforms.
It is obvious, that for Real-time Quantum Chemistry the number of basis functions needs to be as small as possible. Since high accuracy does not need to be the main goal, comparatively small basis sets can be employed. For the same reason pseudopotentials are beneficial. Since plane waves and basis sets of plane waves mixed with Gaussians have successfully been applied in AIMD calculations, they are also useful for Real-time Quantum Chemistry.
3.2 Fock Matrix Calculation
The calculation of the Fock matrix can be divided into two different parts. The calculation of the elements of the one-electron and two-electron matrix operators. The calculation of the exchange–correlation matrix elements occurring in KS-DFT calculations will be discussed separately. In almost all ab initio electronic structure calculations the Fock matrix construction is the most time consuming part, because of the evaluation of all integrals. Therefore, many sophisticated algorithms have been devised to speed up their evaluation.
The calculation of the one-electron Hamiltonian matrix requires the evaluation of matrix elements, where denotes the number of basis functions. The matrix elements consist of two terms, the kinetic energy term and the electron–nuclei interaction part. If localized basis functions are employed, the number of integrals for the kinetic term will increase only linearly with system size. For the Coulomb interaction between the electrons and the nuclei, multipole expansions can be applied to reduce the number of integrals and achieve linear scaling. Although the evaluation of these matrix elements has a quadratic scaling, their contribution to the overall execution time is for mid-sized molecular systems very small and hence no problem for Real-time Quantum Chemistry. However, to even speed up this part, program parallelization can be employed efficiently as the integrals can be evaluated independently from one another.
The two-electron repulsion integrals (ERI) in the calculation of the Coulomb- and the exchange matrix elements are in principle four-index quantities and thus their evaluation formally requires operations. The first and most effective way to reduce their computation time is to discard all elements whose contribution is below a certain threshold. The most widely employed technique for such an integral screening is based on Schwarz-inequality integral estimates.[38, 39, 40] Another approach to screen for negligible matrix elements are the multipole-based integral estimates. They consider in addition the decay behavior between two charge distributions. In non-direct SCF calculations the integrals are pre-calculated and integral screening can only be done at the level of the integrals. In direct SCF methods also the density matrix elements with which the integrals are contracted can be taken into account for screening. Therefore, also large integrals can be neglected, if they have only a small weight assigned by the density matrix elements.
The calculation of the surviving integrals can then be accelerated. The first general approach is to fit the densities occurring in the Coulomb and exchange integrals with auxiliary basis sets, thus reducing the four index ERIs to two index quantities. This approach to fit the densities accelerates the evaluation of the Coulomb matrix elements by an order of magnitude. To obtain the auxiliary basis sets two different methods are employed. One is to determine them by a fitting procedure[42, 43, 44, 45], which is performed before the actual electronic structure calculation for each atom type and basis set and results in additional libraries of auxiliary basis sets. This density-fitting approach is also known as the resolution of the identity (RI)[46, 47, 48, 49]. Another method is to employ a Cholesky decomposition (CD) algorithm to determine the auxiliary basis set for each calculation separately. [50, 51, 52, 53, 54] CD methods are computationally more demanding, but generate an auxiliary basis set, which is “the best” for a given basis and they do not depend on pre-fitted parameters. Therefore, the error introduced by the technique is controllable and no extra auxiliary basis-set libraries are needed.
The calculation of Coulomb matrix elements can also be made more efficient by employing hierarchical multipole methods like the fast multipole method[55, 56, 57] or the quantum chemical tree code. Here, the problem of calculating the Coulomb interaction between many electrons is approximated by a truncated multipole expansion. A combination of a multipole expansion method and auxiliary basis sets can be applied to reduce the computation time even further.
For the exchange matrix elements special methods exist to obtain a linear-scaling behavior. The RI technique, which is very efficient for the Coulomb matrix elements, can also be applied for the exchange matrix elements not with the same efficiency though. In addition several methods ave been developed specifically for the exchange matrix elements. Examples are the -Exchange method (where denotes the number of basis functions and is identical to our ), the local K algorithm, the LinK method or the auxiliary density matrix methods for Hartree–Fock-type exchange.
In KS-DFT calculations the Fock matrix contains additional matrix elements from the exchange–correlation functional. The matrix elements are evaluated by numerical integration of the functional derivative on a grid. The computational cost depends, of course, on the size of the grid. The overall molecular grid is constructed from atomic grids, which are merged to obtain the molecular grid. A very common choice is the Becke atomic partitioning scheme. Since the calculation for each atomic grid can be done independently the overall scaling behavior can be made linear.[65, 66, 67, 68] The exchange–correlation matrix elements are usually considered as requiring a negligible amount of computation time compared with the Coulomb and exchange matrix elements. However, after accelerating the two latter by, e.g., RI methods their evaluation becomes the most time consuming part in the Fock matrix construction.
Considering the potential of the above-mentioned techniques to accelerate the Fock matrix construction, certainly the fast calculation of ERIs is of prime importance. Based on the current state of the art a rigorous screening of the integrals and density-fitting techniques are the methods of choice. Since the integrals can be evaluated independently they are prone for a parallelized evaluation, which will be discussed later. The exchange–correlation matrix elements require efficient numerical techniques for an accelerated calculation. Straightforwardly, the coarseness of the grids employed can be increased to decrease the computational effort within pre-defined accuracy bounds.
3.3 Density-matrix Construction
In conventional SCF-type electronic structure calculations the density matrix is calculated by a full diagonalization of the Fock matrix. This approach has the advantage that the complete eigenspace can be obtained, i.e., including all virtual orbitals. A disadvantage is, however, that the scaling with respect to system size is cubic (and can be easily prohibitive in plane-wave calculations). Due to the small pre-factor of the diagonalization procedure for Gauss-type orbitals it mainly poses a problem for systems with several thousands of atoms.
Almost all methods discussed here have been developed for problems where the construction of the Fock matrix is very fast and therefore the construction of the density matrix becomes the bottleneck of the calculations. An excellent review over methods, which avoid the diagonalization in semi-empirical calculations can be found in Ref. 69. An assessment of density-matrix methods for self-consistent-field calculations by comparing purification and minimization methods has recently been published by Rudberg et al.
The diagonalization of the Fock matrix does not directly yield the density matrix but the coefficient matrix which is then contracted to obtain the density matrix . Although the density matrix is a local quantity due to the nearsightedness of the electrons, the coefficient matrix is not. Therefore, a great reduction of computation time would be achieved calculating the density matrix directly from of the Fock matrix.
The so-called energy minimization techniques exploit the fact that the correct density matrix minimizes the expression . This minimization, however, has to be done under the constraints that the density matrix fulfills the idempotency condition and the trace condition, which requires that the trace of the density matrix yields the number of electrons. By contrast, the diagonalization procedures explicitly fulfill the idempotency condition. Li, Nunes, and Vanderbilt proposed the functional
as a method to compute the density matrix by implicitly fulfilling the idempotency condition. Here, denotes the chemical potential and the identity matrix. Although this functional requires the density matrix to be set up in an orthogonal basis, there is no major problem to obtain a formulation for non-orthogonal basis sets, for which the energy functional is modified to include also the overlap matrix. Independently, Daw proposed a similar approach in 1993 that applied additionally steepest-decent iterations to minimize the functional. Alternatively, a conjugate-gradient approach can be chosen, which is then called the conjugate-gradient density-matrix search method[74, 75]. Then, the Fock and the density matrix are transformed into an orthonormal basis. Another advantage of this approach is that the chemical potential does not need to be known in advance. Other techniques employ for example curvy steps or include second derivatives to obtain faster convergence. All these methods employ a purification of the matrix proposed by McWeeny in order to fulfill the idempotency condition of the density matrix. The so-called sign matrix methods[79, 31] follow a different way to achieve this by expressing the density matrix in terms of the sign matrix function, which can be computed by iterative schemes. Also possible is the introduction of a penalty functional for violating the idempotency condition. However, the penalty functional is difficult to employ since it cannot be minimized by standard methods. Linear scaling of these methods is obtained by weakening the idempotency condition through restricting the minimization to localized density matrices with density matrix elements corresponding to a distance larger than a certain threshold forced to be zero. Other methods, which are not that common in ab initio electronic structure calculations, shall just be mentioned here for the sake of completeness. Examples are the Fermi Operator Expansion[80, 81] and the Fermi Operator Projection method.
Another possibility to minimize the energy without explicitly imposing the orthonormality condition, as it is done by full diagonalization, is to minimize with respect to orbitals only with an implicit orthonormalization constraint. The Orbital Minimization[83, 84, 85, 86] or the Optimal Basis Density Matrix Minimization[87, 88] methods are typical examples. However, they are mainly applied in large tight-binding or semi-empirical calculations. A linear-scaling behavior can be achieved in these methods by carrying out the minimization with respect to localized orbitals (for example, by searching only over functions which non-zero values inside a specified region). This approach would not introduce any approximation if the localized orbitals could be obtained by a unitary transformation of the occupied eigenstates of the Roothaan–Hall equation.
3.4 Acceleration of SCF Convergence
Assuming that the build-up of the Fock matrix and the subsequent calculation of the density matrix are fast enough for Real-time Quantum Chemistry, the SCF procedure poses severe obstacles. The sheer number of SCF iterations strongly depends on the atomic configuration and is hardly predictable. Therefore, it cannot be guaranteed that the wave function and the gradients are available in time. From this point of view, for Real-time Quantum Chemistry it would be desirable to have a method which completely avoids any iterative methods. However, in all true ab initio electronic structure calculation an iterative procedure is unavoidable, since the Fock matrix depends on the density matrix elements. Consequently, to be able to obtain a real-time calculation, one has to focus on reducing the number of iterations to a minimum. Performing the structural manipulations in small steps is therefore the key for a working Real-time Quantum Chemistry implementation.
The steps discussed in the previous sections (Fock-matrix assembly and density-matrix calculation) are both part of a single step in the self-consistent-field procedure. The convergence of SCF iterations strongly depends on the first guess for the molecular orbitals and on the nuclear configuration. In a reactivity study, however, it may happen that configurations of the atoms occur, for which the SCF procedure converges only slowly or not at all. It is thus of utmost importance to have methods at hand, which yield stable and fast converging SCF iterations.
Direct inversion in the iterative subspace (DIIS)[90, 91] is a widely employed technique to accelerate and stabilize SCF iterations. Error vectors from previous iterations are calculated and minimized in a least-squares sense. Accordingly, previous iterations are utilized to extrapolate the Fock matrix in the next iteration. Quite closely related to the DIIS method are techniques called Fock matrix dynamics or electron density extrapolation  which are common in the field of BO molecular dynamics. Instead of accelerating the SCF algorithm itself the whole single-point SCF calculation is accelerated by extrapolating the information from previous time steps of the simulation. It is assumed that in between two steps the nuclear coordinates change only little so that the molecular orbitals, Fock matrix elements, or the density do not differe much from one step to the next and hence are a good starting point for the electronic structure optimization for the new nuclear configuration. Although SCF acceleration schemes have a long history, significant improvements can still be made as highlighted by the augmented Roothaan–Hall method.
The pseudo-diagonalization technique is based on the observation that only for the first and the last SCF iteration a full diagonalization is necessary. In between it is sufficient to eliminate all Fock matrix elements connecting the occupied and virtual molecular orbitals by unitary transformations. As a consequence, the diagonalization of the Fock matrix blocks corresponding to the virtual–virtual and the occupied–occupied molecular orbitals can be avoided, which greatly reduces the time of each SCF iteration.
In cases of a too narrow gap between the highest occupied and the lowest unoccupied molecular orbital, slow convergence or even divergent SCF iterations can occur. In such cases level shifting can be applied to enlarge the gap and therefore avoid a mixing. With this procedure the problem of divergence, slow convergence, or oscillating behavior can be cured in many cases.
In all techniques discussed above the construction of the density matrix and the optimization of the molecular orbitals in the SCF iterations were independent processes. However, for methods based on minimization of an energy functional for constructing the density matrix, methods have been developed that combine the density matrix optimization and the self-consistent-field iterations in one single optimization loop.[97, 98, 99, 100] The idea was developed for the coupled electron-nuclei problem (Car–Parrinello molecular dynamics) and is therefore often called the “molecular dynamics” method, but it can also be applied in situations were the nuclei are kept fixed. The difference to conventional matrix diagonalization procedures to solve for the eigenstates is that the variational principle is applied in a dynamic fashion and all eigenstates are determined simultaneously. The dynamic variables are here the coefficients of the basis functions with a fictitious electon mass. To exploit this methodology for a Real-time Quantum Chemistry implementation one could require that the manipulations are done in a continuous fashion, i.e., along a ’trajectory’ with sufficiently small configuration-change steps.
3.5 Gradient Calculation
As outlined in the beginning of this section the calculation of the energy gradient with respect to the position of the nuclear coordinates involves contributions from each term in the electronic energy. An efficient force evaluation for large molecular systems in the framework of pure DFT has recently been proposed by Reine et al. by combining screening with a fast multipole method. In addition the calculation of the Coulomb contributions were accelerated by employing a density-fitting scheme[24, 25] with auxiliary basis sets. There are also efficient gradient implementations available that do not not employ density fitting[102, 103].
3.6 Subsystem Approaches
To divide a molecular system under study into smaller subsystems is a key to the molecular-model approach discussed in the first two sections of this work and offers the possibility to reduce the computational effort further. As through a magnifying glass, the reactive part/region of a molecular assembly can be embedded in a spectator background and this magnifying lens can even be moved around in the whole system. So-called combined quantum mechanics/molecular mechanics (QM/MM)[105, 106] approaches are widely employed to enable calculations of large enzymatic systems. In QM/MM the reactive part is treated quantum mechanically and the surrounding environment is modeled by classical force fields. By contrast, QM/QM methods apply the laws of quantum mechanics to all subsystem but may treat them with different methodologies.
Depending on how the subsystems are embedded into each other a different level of acceleration can be achieved. Completely independent subsystem methods—also called Divide-and-Conquer (D&C) approaches—facilitate a massively parallelized calculation. In the original formulation of the D&C approach the independently calculated density matrices of the subsystems were merged to yield the full density matrix. Although the computation of the density matrices is done independently, the surrounding of a fragment is accounted for by buffer regions around the fragment. This approach does not only accelerate the calculation of the density matrix but also the calculation of the Fock matrices in the SCF calculations, which is the reason why this methodology is treated in this separate section and not together with other density matrix construction schemes. One difficulty for Real-time Quantum Chemistry is, however, that the calculation of the D&C force for the full system is not well defined. But there are cures available. The partitioning can also be carried out for other quantities than the density-matrix. For instance, in the fragment molecular orbital theory the fragmentation is done at the level of the molecular orbitals. For a recent and comprehensive review on fragmentation methods we refer to Ref. 111
Subsystem techniques allow for an in principle exact embedding of the subsystems thus retaining only the approximations introduced by choosing different electronic-structure methods for the subsystems. In the framework of density functional theory such methods have been proposed and are widely employed, for instance, to account for solvent effects[112, 113, 114, 115, 116, 117, 118, 119].
To assess the potential acceleration gained by any subsystem approach, one has to consider two aspects: (i) how small can the subsystems be chosen and (ii) how are they connected to each other, if at all. This affects not only the computation but also the accuracy. But since we are interested in local phenomena, the fragmentation is often already implied by the structure of the molecular system itself. Clearly, the smaller the subsystems can be chosen and the less one has to account for embedding effects, the better for Real-time Quantum Chemistry.
From the Real-time Quantum Chemistry point of view, Divide-and-Conquer and density embedding approaches are appealing, since they allow the greatest reduction in computation time if the whole system needs to be treated quantum mechanically. If, however, large parts of the molecular system are not directly involved in a reaction, but rather serve as a dielectric environment, then QM/MM methods are most suitable.
3.7 Technical Aspects: Parallelization and Special Hardware
Essentially all of the most popular quantum chemistry codes can be run in parallel. However, these software approaches usually have a rather high overhead making them not the best candidates for Real-time Quantum Chemistry applications. But electronic structure calculations have not only been accelerated by developing increasingly efficient algorithms but also by employing latest and even specialized hardware.
Although graphical processing units (GPUs) were originally designed for the fast rendering of three dimensional graphics, they can be employed for the acceleration of quantum chemical calculations. In the past years the development and application of special algorithms for quantum chemical calculations that exploit the computing power of GPUs has gained considerable attention.[120, 121, 122, 123, 124, 125] This trend is due to a significant effort undertaken to program quantum chemical algorithms specifically designed for GPUs but also due to new graphic cards produced with a focus on scientific calculations. For example, the widely used GAMESS US package or the Terachem program are able to efficiently exploit the advantages of GPUs for electronic structure calculations. Compared to calculations performed only on the central processing unit (CPU) of a computer they are able to achieve considerable speed-ups in the order of one magnitude. Also in the field of semi-empirical calculations GPUs have attracted some attention.
Massive parallelization is also possible on processor architectures other than GPUs. For Kohn–Sham DFT calculation, for instance, the application of processors from ClearSpeed Technology Ltd. has been reported to accelerate electronic structure and gradient calculations. [129, 130] Besides these developments for specific processor types, there are also some more general considerations about how to exploit the advantages of emerging new processor types available in the literature.[131, 132] A comprehensive overview of special processors and their potential for electronic structure calculation can be found in Ref. 133. In the field of classical molecular dynamics simulations the so-called ANTON processor developed by Shaw et al. is a successful attempt to build such a special purpose processor. Also for electronic structure calculations special processors have been designed; namely ERIC, the ERI Calculation specific chip-multiprocessor or the molecular orbital calculation specific embedded high performance computing (EHPC) system.
Almost all special processor types mentioned here accelerate the electronic structure calculations by parallelization of the ERI calculation. This is the reason for the success of GPUs but also of high performance clusters like, for instance, the IBM Blue Gene series[137, 138]. The calculations of ERIs is in almost all quantum chemical calculation the bottleneck because of their sheer number. Specialized hardware which allows for a massive parallelization can directly accelerate the calculations and not only improves the scaling behavior. Thus, exploiting these techniques is certainly imperative for Real-time Quantum Chemistry.
4 Direct Haptic Quantum Chemistry
After having elaborated on the available and future quantum chemical methods for Real-time Quantum Chemistry implementations, we shall now discuss their benefits for Haptic Quantum Chemistry[6, 7] and subsequently discuss their capabilities in an out-of-the-box application presented in the next section.
The concept and implementation of Haptic Quantum Chemistry as presented in Refs. 6, 7 is to employ a force-feedback device as depicted in Fig. 1 as an input and an output tool allowing for an intuitive manipulation of molecular structures while feeling the gradients on the manipulated atoms as forces. Such approaches are referred to as haptic enabled interactive molecular visualizations systems in the literature. There are a few such methods already available, but they only employ classical force fields to calculate the forces rendered by the device.[140, 141, 142, 143] Thus, they prohibit the study of chemical reactions, which would require the ability tp form and break chemical bonds.
In our current set-up, the haptic device is a pen-like pointer which allows the user to manipulate objects in a virtual reality framework and, at the same time, to physically experience a force feedback (cf. Fig. 1). Hence, one is able to feel the curvature of the potential energy surface of the manipulated nuclear coordinates in the whole system. The visual presentation of the structure, the gradients on the atoms not manipulated with the device, the orbitals, the electron density and other properties allow one to perceive complex information in a very intuitive way compared with the sole visual presentation. Hence, the haptic quantum chemical approach immerses the scientist more into the scientific problem. Even deeper immersion can of course be achieved by employing 3D displays or 3D glasses and other techniques from the field of virtual reality.
Since the human haptic sense is much more sensitive than the visual sense, haptic devices usually have an update rate of about . By contrast, to create the illusion of a smooth movement for the visual sense only are necessary. To circumvent this problem, in Haptic Quantum Chemistry[6, 7] so far the forces are not calculated directly but are obtained from interpolating single-point gradients . The force acting on an atom is calculated from the interpolated gradient by
where the gradient of the energy is given by
Here, denotes the three Cartesian derivatives with respect the nuclear coordinate of atom . The computationally inexpensive interpolation facilitates the very fast calculation of the forces necessary for the high update rate of the haptic device. The main drawback is, however, that a coarse-grained gradient field of the molecular system needs to be calculated in advance (though it can be refined during haptic exploration as more data points can be calculated in the background).
The pre-calculation of single-point gradients becomes, however, more and more demanding if parts of the molecular structure are allowed to relax. The molecular structure at one specific position in the configuration space of the mobile part is no longer unique and depends on the trajectory leading to it. As a consequence, one has to design algorithms that keep track of the history of the actual haptic exploration run. However, employing a Real-time Quantum Chemistry framework would allow us to circumvent this problem as the quantum chemical data is always immediately available so that no history (trajectory) needs to be stored. We may call this approach Direct Haptic Quantum Chemistry (D-HQC).
At first glance, the high update rate (more than 1 kHz) of the haptic device requires an update rate of a millisecond if the gradients from the electronic structure calculation were directly rendered by the device. But as it has been shown in the field of interactive fluid dynamics simulations, the servo loop of the haptic device and the simulation loop can be separated and can operate with different update rates. Both loops are connected by a shared memory from where the servo loop constantly reads the current force written by the simulation loop as soon as the new force is available. To reduce the occurring artifacts the stepwise force update is smoothed by a force filter. With this technique the update rate of the molecular system can be lowered below . Accordingly, D-HQC would require that within a few hundred millisecond a new gradient on the manipulated atoms has to be available and the other atom positions have to be relaxed (Fig. 2).
This almost instantaneous relaxation of the whole structure is, however, not always wanted. By contrast, it might be even desired to be able to alter the molecular structure faster than the system relaxes in order to simulate non-adiabatic changes. Therefore, the relaxation does not need to be instantaneous as long as the force change due to the relaxation can be rendered smoothly. For almost adiabatic changes the structure alterations have to be slowed down or done in small steps so that the structure of the whole system can relax. Altering the structure only in very small increments also speeds up the electronic structure optimizations, since the molecular orbitals change only very little. The program could force the user to perform only small changes (slow movements) by applying re-stalling forces on the haptic device.
D-HQC can be described as probing the potential energy surface in the configuration subspace spanned by the manipulated atoms of the molecular system. For a more formal description of the force calculation in D-HQC the set of nuclei is partitioned into nuclei controlled by the haptic device and the remaining nuclei . The force on a nucleus is then calculated as the negative spatial derivative of the total energy, which is minimized with respect to the basis set and the nuclear coordinates of the remaining nuclei .
with the total energy functional is given by
Many of the above-mentioned developments of linear and sub-linear scaling methods are available in standard program packages, though the ultimate package for real time quantum chemistry has not been developed yet—mostly for the reason that each package has been designed to serve a certain purpose. For example, huge molecular models have already been studied with the CP2K [28, 29, 30] program that employs mixed Gaussian and plane wave basis sets in AIMD simulations. Here, we choose the very efficient DFT modules of the Turbomole package[145, 146] combined with our D-HQC setup to demonstrate that such studies are in reach for molecular models of relevant size. The calculations exploit density fitting, effective core potentials and small basis sets. It is clear that the resulting accuracy then does not necessarily live up to the current standard. However, this is also not decisive as a reaction pathway recorded during a D-HQC exploration can always be relaxed on a more accurate potential energy surface.
5.1 D-HQC for a bromine molecule approaching an ethene
As depicted in Fig. 3 in this exemplary study a bromine molecule was pushed onto an ethene molecule to probe its reactivity. In this setup the position of the bromine atom closest to the ethene was dragged towards the ethene molecule so that the relaxation at each step had the bromine to ethene distance as a constraint. In this way a trajectory in the subspace spanned by the position of the bromine atom was generated. The distance was incrementally decreased by 0.5 Bohr. The resulting distances in Ångstroms are given in Fig. 3, which shows three exemplary points from the trajectory. In addition, the average time per structure relaxation step, i.e., the time for updating the system’s structure, as well as the overall time needed to converge the structure (second line) are depicted.
Following the results of the discussion about basis sets we chose for the bromine atoms the Stuttgart ECP-28-MWP pseudo potential and the def2-SV(P) basis set, for the carbon atom the def2-SV(P) basis set and for the hydrogen atoms the STO-3G HONDO basis set[149, 150] in order to obtain a very fast calculation of single-point energies and gradients. The calculations were performed employing the BP86 exchange–correlation density functional[151, 152, 153, 154, 155] on a coarse numerical grid. In addition, also the resolution-of-identity technique (RI) was applied to accelerate the calculations.
The execution times in Fig. 3 show that the time needed to update the system is almost constant during the trajectory, but the structure optimization time increases when the reactants get closer, which indicates that it needs more steps to converge. As it was outlined before, the update rate is the important quantity for a real-time experience. The execution times per update step can be shortened by sampling the trajectory in smaller steps, which means that the SCF procedure can converge faster since the wave function does not change too much from step to step.
The computations here were performed by running the individual Turbomole modules sequentially. Note that there is still room for improving the efficiency in terms of passing the information from one calculation to the next and avoid read/write accesses to the hard disk as these processes have not yet been optimized for a D-HQC implementation in standard programs like Turbomole.
5.2 D-HQC for an S2 reaction of fluoride with chloromethane
Another example with a more pronounced effect of structural relaxation of the remaining atom positions is expected to demonstrate how this influences the systems update rate: a S2 type reaction. In this example a fluoride ion approaches a chloromethane molecule and replaces the chlorine atom. The trajectory was recorded by moving the fluoride anion in steps of 0.1 Bohr towards the C atom. At each step the structure of the remaining atoms was optimized. The intermediate electronic structure optimizations were carried out again with small basis sets (def2-SVP  for the C, F, and Cl atoms and STO-3G HONDO [149, 150] for the H atoms) in combination with effective core potentials (ECP-10-MWB for Cl and ECP-2-SDF for F) and BP86/RI on a coarse numerical grid. To reduce the number of SCF cycles in each geometry optimization cycle the molecular orbitals from the preceding point of the trajectory were taken.
In Figure 4 three intermediate points of the trajectory are shown. The average time per electronic structure calculation in the structure optimization is almost constant and on the order of several hundred milliseconds. The overall time for the structure optimizer to reach convergence is, however, not constant and increases significantly when the attacking fluoride approaches its position in the transition state. The detailed timings show again that not the electronic structure optimization but the increased number of cycles in the structure optimization give rise to the increased overall time. As already discussed in the previous section this is not a severe issue as only the system updates have to be fast, which is this case in this example. The user would have to slow down the movement of the fluoride when approaching the C atom in order to obtain a reasonable minimum energy path. For a non-adiabatic simulation the movement can be faster although the remaining atoms of the system are not able to relax in time.
The possibility and necessity of obtaining the result of quantum chemical calculations in real time is beneficial in many respects for studying the reactivity of chemical systems and may change the way how quantum chemistry is done in the future. Not only the ever increasing amount of information provided by calculations but also the inherent complexity of chemical problems calls for new approaches like (Direct) Haptic Quantum Chemistry that require a Real-time Quantum Chemistry framework. The overview of currently available techniques to accelerate calculations provided here clearly showed that Real-time Quantum Chemistry is in reach and will be possible for relevant system sizes in the near future.
The evaluation of existing algorithms and technology for Real-time Quantum Chemistry also demonstrated, however, that a paradigm change is needed. Almost all techniques presented here were not specifically designed to allow quantum chemical calculations of energies and gradients in real time. The aim of their development was the overall scaling behavior to allow the treatment of ever larger molecules or molecular systems. For Real-time Quantum Chemistry the focus needs to be on reducing the actual execution time for a fixed system size to around . Although most of the currently available program packages in quantum chemistry have not been developed to allow the ultra-fast calculation of molecular systems consisting of atoms, the greatest potential for achieving considerable speed-up towards real time lies most probably in the activation of specialized hardware. Already in reach are calculations on GPUs which show a promising potential, but also completely new specialized hardware is desirable for Real-time Quantum Chemistry.
The overwhelming amount and the complexity of the data generated by current quantum chemical calculations already limits their fast and intuitive evaluation. Haptic Quantum Chemistry offers a new approach to tackle this problem. The instantaneous availability of the wave functions and the gradients offered by Real-time Quantum Chemistry allows an even more convenient way of studying chemical reactivity, as we have discussed for the Direct Haptic Quantum Chemistry variant. The exploitation of the human haptic sense to present scientific data more intuitively is only a first step. A deeper immersion by employing techniques already developed in the field of virtual reality would be the ultimate goal of any development in this direction.
We gratefully acknowledge financial support through TH grant ETH-08 11-2.
-  Bellman, R. Dynamic Programming; Princeton University Press, 1957.
-  Huber, T.; Torda, A. E.; Gunsteren, W. F. J. Comput. Aid. Mol. Des. 1994, 8, 695–708.
-  Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D. J. Chem. Phys. 1998, 108, 1964–1977.
-  Iannuzzi, M.; Laio, A.; Parrinello, M. Phys. Rev. Lett. 2003, 90, 238302.
-  Wales, D. J. Int. Rev. Phys. Chem. 2006, 25, 237–282.
-  Marti, K. H.; Reiher, M. J. Comput. Chem. 2009, 30, 2010–2020.
-  Haag, M. P.; Marti, K. H.; Reiher, M. ChemPhysChem 2011, 12, 3204–3213.
-  Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods; Cambridge University Press, 2009.
-  Grayson, P.; Tajkhorshid, E.; Schulten, K. Biophys. J. 2003, 85, 36 – 48.
-  Bosson, M.; Grudinin, S.; Bouju, X.; Redon, S. J. Comput. Phys. 2011, 231, 2581–?2598.
-  Bosson, M.; Richard, C.; Plet, A.; Grudinin, S.; Redon, S. J. Comput. Chem. 2012, 33, 779–790.
-  Goedecker, S. Rev. Mod. Phys. 1999, 71, 1085–1123.
-  Rubensson, E. H.; Rudberg, E.; Salek, P. In Linear-Scaling Techniques in Computational Chemistry and Physics; Zalesny, R.; Papadopoulos, M. G.; Mezey, P. G.; Leszczynski, J., Eds.; Springer Netherlands, 2011, Vol. 13, pp 263–300.
-  Siegbahn, P. E.; Himo, F. WIREs Comput. Mol. Sci. 2011, 1, 323–336.
-  Podewitz, M.; Stiebritz, M. T.; Reiher, M. Faraday Discuss. 2011, 148, 119–135.
-  Ochsenfeld, C.; Kussmann, J.; Lambrecht, D. S. In Linear-Scaling Methods in Quantum Chemistry; John Wiley & Sons, Inc., 2007; pp 1–82.
-  Szabó, A.; Ostlund, N. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; Dover Books on Chemistry Series; Dover Publications, 1996.
-  Hohenberg, P.; Kohn, W. Phys. Rev. 1964, 136, B864–B871.
-  Kohn, W.; Sham, L. J. Phys. Rev. 1965, 140, A1133–A1138.
-  Parr, R.; Yang, W. Density-Functional Theory of Atoms and Molecules; International Series of Monographs on Chemistry; Oxford University Press, 1994.
-  Rudberg, E.; Rubensson, E. H.; Sa?ek, P. J. Chem. Theory Comput. 2011, 7, 340–350.
-  Pulay, P. In Analytical Derivative Methods in Quantum Chemistry; John Wiley & Sons, Inc., 1987; pp 241–286.
-  Satoko, C. Phys. Rev. B 1984, 30, 1754–1764.
-  Versluis, L.; Ziegler, T. J. Chem. Phys. 1988, 88, 322–328.
-  Fournier, R.; Andzelm, J.; Salahub, D. J. Chem. Phys. 1989, 90, 6371–6377.
-  Delley, B. J. Chem. Phys. 1991, 94, 7245–7250.
-  Pople, J. A.; Gill, P. M.; Johnson, B. G. Chem. Phys. Lett. 1992, 199, 557 – 560.
-  Lippert, G.; Parrinello, M.; Hutter, J. Mol. Phys. 1997, 92, 477–488.
-  Lippert, G.; Hutter, J.; Parrinello, M. Theor. Chem. Acc. 1999, 103, 124–140.
-  VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Comput. Phys. Commun. 2005, 167, 103 – 128.
-  VandeVondele, J.; Borštnik, U.; Hutter, J. J. Chem. Theory Comput. 2011, in press, DOI: 10.1021/ct200897x.
-  Schwerdtfeger, P. ChemPhysChem 2011, 12, 3143–3155.
-  Dolg, M.; Cao, X. Chem. Rev. 2012, 112, 403–480.
-  Reiher, M.; Wolf, A. Relativistic Quantum Chemistry; Wiley-VCH: Weinheim, 2009.
-  Cho, K.; Arias, T. A.; Joannopoulos, J. D.; Lam, P. K. Phys. Rev. Lett. 1993, 71, 1808–1811.
-  Han, S.; Cho, K.; Ihm, J. Phys. Rev. B 1999, 60, 1437–1440.
-  Genovese, L.; Neelov, A.; Goedecker, S.; Deutsch, T.; Ghasemi, S. A.; Willand, A.; Caliste, D.; Zilberberg, O.; Rayson, M.; Bergman, A.; Schneider, R. J. Chem. Phys. 2008, 129, 014109.
-  Almlöf, J.; Faegri, K.; Korsell, K. J. Comput. Chem. 1982, 3, 385–399.
-  Häser, M.; Ahlrichs, R. J. Comput. Chem. 1989, 10, 104–111.
-  Cremer, D.; Gauss, J. J. Comput. Chem. 1986, 7, 274–282.
-  Lambrecht, D. S.; Ochsenfeld, C. J. Chem. Phys. 2005, 123, 184101.
-  Whitten, J. L. J. Chem. Phys. 1973, 58, 4496–4501.
-  Baerends, E.; Ellis, D.; Ros, P. Chem. Phys. 1973, 2, 41 – 51.
-  Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. International Journal of Quantum Chemistry 1977, 12, 81–87.
-  Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. J. Chem. Phys. 1979, 71, 3396–3402.
-  Vahtras, O.; Almlöf, J.; Feyereisen, M. Chem. Phys. Lett. 1993, 213, 514 – 518.
-  Eichkorn, K.; Treutler, O.; Öhm, H.; Häser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 240, 283–290.
-  Skylaris, C.-K.; Gagliardi, L.; Handy, N.; Ioannou, A.; Spencer, S.; Willetts, A. Comp. Theor. Chem. 2000, 501?502, 229 – 239.
-  Weigend, F. Phys. Chem. Chem. Phys. 2002, 4, 4285–4291.
-  Beebe, N. H. F.; Linderberg, J. Int. J. Quantum Chem. 1977, 12, 683–705.
-  Koch, H.; de Meras, A. S.; Pedersen, T. B. J. Chem. Phys. 2003, 118, 9481–9484.
-  Boström, J.; Aquilante, F.; Pedersen, T. B.; Lindh, R. J. Chem. Theory Comput. 2009, 5, 1545–1553.
-  Aquilante, F.; Lindh, R.; Pedersen, T. B. J. Chem. Phys. 2007, 127, 114107.
-  Aquilante, F.; Boman, L.; Boström, J.; Koch, H.; Lindh, R.; Merás, A. S.; Pedersen, T. B. In Linear-Scaling Techniques in Computational Chemistry and Physics; Zalesny, R.; Papadopoulos, M. G.; Mezey, P. G.; Leszczynski, J., Eds.; Springer Netherlands, 2011, Vol. 13, pp 301–343.
-  Greengard, L.; Rokhlin, V. J. Comput. Phys. 1987, 73, 325–348.
-  White, C. A.; Head-Gordon, M. J. Chem. Phys. 1994, 101, 6593–6605.
-  Strain, M. C.; Scuseria, G. E.; Frisch, M. J. Science 1996, 271, 51–53.
-  Challacombe, M.; Schwegler, E. J. Chem. Phys. 1997, 106, 5526–5536.
-  Sierka, M.; Hogekamp, A.; Ahlrichs, R. J. Chem. Phys. 2003, 118, 9136–9148.
-  Schwegler, E.; Challacombe, M. J. Chem. Phys. 1996, 105, 2726–2734.
-  Aquilante, F.; Pedersen, T. B.; Lindh, R. J. Chem. Phys. 2007, 126, 194106.
-  Ochsenfeld, C.; White, C. A.; Head-Gordon, M. J. Chem. Phys. 1998, 109, 1663–1669.
-  Guidon, M.; Hutter, J.; VandeVondele, J. J. Chem. Theory Comput. 2010, 6, 2348–2364.
-  Becke, A. D. J. Chem. Phys. 1988, 88, 2547–2553.
-  van Wüllen, C. Chem. Phys. Lett. 1994, 219, 8 – 14.
-  PÃ©rez-Jorda, J.; Yang, W. Chem. Phys. Lett. 1995, 241, 469 – 476.
-  Stratmann, R. E.; Scuseria, G. E.; Frisch, M. J. Chem. Phys. Lett. 1996, 257, 213–223.
-  Burow, A. M.; Sierka, M. J. Chem. Theory Comput. 2011, 7, 3097–3104.
-  Daniels, A. D.; Scuseria, G. E. J. Chem. Phys. 1999, 110, 1321–1328.
-  Kohn, W. Phys. Rev. Lett. 1996, 76, 3168–3171.
-  Li, X.-P.; Nunes, R. W.; Vanderbilt, D. Phys. Rev. B 1993, 47, 10891–10894.
-  Nunes, R. W.; Vanderbilt, D. Phys. Rev. B 1994, 50, 17611–17614.
-  Daw, M. S. Phys. Rev. B 1993, 47, 10895–10898.
-  Millam, J. M.; Scuseria, G. E. J. Chem. Phys. 1997, 106, 5569–5577.
-  Challacombe, M. J. Chem. Phys. 1999, 110, 2332–2342.
-  Helgaker, T.; Larsen, H.; Olsen, J.; Jørgensen, P. Chem. Phys. Lett. 2000, 327, 397–403.
-  Ochsenfeld, C.; Head-Gordon, M. Chem. Phys. Lett. 1997, 270, 399 – 405.
-  McWeeny, R. Rev. Mod. Phys. 1960, 32, 335–369.
-  Beylkin, G.; Coult, N.; Mohlenkamp, M. J. J. Comput. Phys. 1999, 152, 32 – 54.
-  Goedecker, S.; Teter, M. Phys. Rev. B 1995, 51, 9455–9464.
-  Goedecker, S.; Colombo, L. Phys. Rev. Lett. 1994, 73, 122–125.
-  Goedecker, S. J. Comput. Phys. 1995, 118, 261 – 268.
-  Kim, J.; Mauri, F.; Galli, G. Phys. Rev. B 1995, 52, 1640–1648.
-  Mauri, F.; Galli, G. Phys. Rev. B 1994, 50, 4316–4326.
-  Mauri, F.; Galli, G.; Car, R. Phys. Rev. B 1993, 47, 9973–9976.
-  Ordejón, P. Comp. Mater. Sci. 1998, 12, 157 – 191.
-  Hernández, E.; Gillan, M. J. Phys. Rev. B 1995, 51, 10157–10160.
-  Hierse, W.; Stechel, E. B. Phys. Rev. B 1994, 50, 17811–17819.
-  Galli, G.; Parrinello, M. Phys. Rev. Lett. 1992, 69, 3547–3550.
-  Pulay, P. Chem. Phys. Lett. 1980, 73, 393 – 398.
-  Pulay, P. J. Comput. Chem. 1982, 3, 556–560.
-  Pulay, P.; Fogarasi, G. Chem. Phys. Lett. 2004, 386, 272 – 278.
-  Niklasson, A. M. N.; Tymczak, C. J.; Challacombe, M. Phys. Rev. Lett. 2006, 97, 123001.
-  Høst, S.; Olsen, J.; Jansik, B.; Thøgersen, L.; JÃ¸rgensen, P.; Helgaker, T. J. Chem. Phys. 2009, 129, 124106.
-  Stewart, J. J. P.; Császár, P.; Pulay, P. J. Comput. Chem. 1982, 3, 227–228.
-  Saunders, V. R.; Hillier, I. H. Int. J. Quantum Chem. 1973, 7, 699–705.
-  Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471–2474.
-  Štich, I.; Car, R.; Parrinello, M.; Baroni, S. Phys. Rev. B 1989, 39, 4997–5004.
-  Teter, M. P.; Payne, M. C.; Allan, D. C. Phys. Rev. B 1989, 40, 12255–12263.
-  Payne, M. C.; Teter, M. P.; Allan, D. C.; Arias, T. A.; Joannopoulos, J. D. Rev. Mod. Phys. 1992, 64, 1045–1097.
-  Reine, S.; Krapp, A.; Iozzi, M. F.; Bakken, V.; Helgaker, T.; Pawlowski, F.; Salek, P. J. Chem. Phys. 2010, 133, 044102.
-  Burant, J. C.; Strain, M. C.; Scuseria, G. E.; Frisch, M. J. Chem. Phys. Lett. 1996, 248, 43 – 49.
-  Shao, Y.; White, C. A.; Head-Gordon, M. J. Chem. Phys. 2001, 114, 6572–6577.
-  Gascon, J. A.; Leung, S. S. F.; Batista, E. R.; Batista, V. S. J. Chem. Theory Comput. 2006, 2, 175–186.
-  Senn, H.; Thiel, W. Top. Curr. Chem. 2007, 268, 173–290.
-  Lin, H.; Truhlar, D. Theor. Chem. Acc. 2007, 117, 185–199.
-  Humbel, S.; Sieber, S.; Morokuma, K. J. Chem. Phys. 1996, 105, 1959–1967.
-  Yang, W. Phys. Rev. Lett. 1991, 66, 1438–1441.
-  Zhao, Q.; Yang, W. J. Chem. Phys. 1995, 102, 9598–9603.
-  Fedorov, D. G.; Kitaura, K. J. Phys. Chem. A 2007, 111, 6904–6914.
-  Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Chem. Rev. 2012, 112, 632–672.
-  Wesolowski, T. A.; Warshel, A. J. Phys. Chem. 1993, 97, 8050–8053.
-  Wesolowski, T. A. In Computational Chemistry: Reviews of Current Trends; Leszczynski, J., Ed.; World Scientific: Singapore, 2006, Vol. 10, p 82.
-  Iannuzzi, M.; Kirchner, B.; Hutter, J. Chem. Phys. Lett. 2006, Chem. Phys. Lett., 16–20.
-  Fux, S.; Kiewisch, K.; Jacob, C. R.; Neugebauer, J.; Reiher, M. Chem. Phys. Lett. 2008, 461, 353–359.
-  Kiewisch, K.; Eickerling, G.; Reiher, M.; Neugebauer, J. J. Chem. Phys. 2008, 128, 044114.
-  Jacob, C. R.; Neugebauer, J.; Visscher, L. J. Comput. Chem. 2008, 29, 1011–1018.
-  Neugebauer, J. Phys. Rep. 2010, 489, 1–87.
-  Severo Pereira Gomes, A.; Jacob, C. R. Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 2012, 108, 222–277.
-  Yasuda, K. J. Chem. Theory Comput. 2008, 4, 1230–1236.
-  Yasuda, K. J. Comput. Chem. 2008, 29, 334–342.
-  Ufimtsev, I. S.; Martinez, T. J. J. Chem. Theory Comput. 2008, 4, 222–231.
-  Ufimtsev, I. S.; Martinez, T. J. J. Chem. Theory Comput. 2009, 5, 2619–2628.
-  Ufimtsev, I. S.; Martinez, T. J. J. Chem. Theory Comput. 2009, 5, 1004–1015.
-  Luehr, N.; Ufimtsev, I. S.; Marti?nez, T. J. J. Chem. Theory Comput. 2011, 7, 949–954.
-  Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347–1363.
-  Stone, J. E.; Hardy, D. J.; Ufimtsev, I. S.; Schulten, K. J. Mol. Graph. Model. 2010, 29, 116 – 125.
-  Wu, X.; Koslowski, A.; Thiel, W. J. Chem. Theory Comput. 2012, in press.
-  Brown, P.; Woods, C.; McIntosh-Smith, S.; Manby, F. R. J. Chem. Theory Comput. 2008, 4, 1620–1626.
-  Brown, P.; Woods, C. J.; McIntosh-Smith, S.; Manby, F. R. J. Comput. Chem. 2010, 31, 2008–2013.
-  Ramdas, T.; Egan, G. K.; Abramson, D.; Baldridge, K. K. Comput. Phys. Commun. 2008, 178, 817 – 834.
-  Ramdas, T.; Egan, G. K.; Abramson, D.; Baldridge, K. K. Comput. Phys. Commun. 2009, 180, 1221 – 1229.
-  Ramdas, T.; Egan, G.; Abramson, D.; Baldridge, K. Theor. Chem. Acc. 2008, 120, 133–153.
-  Shaw, D. E. et al. SIGARCH Comput. Archit. News 2007, 35, 1–12.
-  Nakamura, K.; Honda, H.; Inoue, K.; Sato, H.; Uehara, M.; Komatsu, H.; Umeda, H.; Inadomi, Y.; Araki, K.; Sasaki, T. et al. Proceedings of the workshop on unique chips and systems (UCAS), 2005; pp 87–94.
-  Umeda, H.; Inadomi, Y.; Honda, H.; Nagashima, U. J. Comput. Chem. 2009, 30, 826–831.
-  IBM Blue Gene team, IBM J. Res. Dev. 2008, 52, 199 –220.
-  Fletcher, G. D.; Fedorov, D. G.; Pruitt, S. R.; Windus, T. L.; Gordon, M. S. J. Chem. Theory Comput. 2012, 8, 75–79.
-  Marchese, F. T. In Trends in Interactive Visualization; Liere, R.; Adriaansen, T.; Zudilova-Seinstra, E., Eds.; Advanced Information and Knowledge Processing; Springer London, 2009, pp 251–267.
-  Bayazit, O.; Song, G.; Amato, N. Robotics and Automation, 2001. Proceedings 2001 ICRA. IEEE International Conference on, 2001; pp 954 – 959.
-  Lai-Yuen, S.; Lee, Y.-S. Computer Aided Design and Computer Graphics, 2005. Ninth International Conference on, 2005; p 6 pp.
-  Bidmon, K.; Reina, G.; Bos, F.; Pleiss, J.; Ertl, T. EuroHaptics Conference, 2007 and Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems. World Haptics 2007. Second Joint, 2007; pp 537 –542.
-  Wollacott, A. M.; Jr., K. M. M. J. Mol. Graph. Model. 2007, 25, 801 – 805.
-  Baxter, W.; Lin, M. C. Proceedings of Graphics Interface 2004, School of Computer Science, University of Waterloo, Waterloo, Ontario, Canada, 2004; pp 81–88.
-  Treutler, O.; Ahlrichs, R. J. Chem. Phys. 1995, 102, 346–354.
-  von Arnim, M.; Ahlrichs, R. J. Comput. Chem. 1998, 19, 1746–1757.
-  Schwerdtfeger, P.; Dolg, M.; Schwarz, W. H. E.; Bowmaker, G. A.; Boyd, P. D. W. J. Chem. Phys. 1989, 91, 1762–1774.
-  Schäfer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571–2577.
-  Hehre, W. J.; Ditchfield, R.; Stewart, R. F.; Pople, J. A. J. Chem. Phys. 1970, 52, 2769–2773.
-  Hehre, W. J.; Stewart, R. F.; Pople, J. A. J. Chem. Phys. 1969, 51, 2657–2664.
-  Dirac, P. A. M. Proc. R. Soc. Lond. A 1929, 123, 714–733.
-  Slater, J. C. Phys. Rev. 1951, 81, 385–390.
-  Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200–1211.
-  Becke, A. D. Phys. Rev. A 1988, 38, 3098–3100.
-  Perdew, J. P. Phys. Rev. B 1986, 33, 8822–8824.