Can DFT simulate chemical reactions?

Question

I was wondering, can you simulate the chemical reaction between two or more molecules using DFT? If not, what numerical method is usually used for this?

Answer 1

DFT can be used for many things. One of them is to calculate an approximation of the ground state energy of a molecule. By calculating the ground state energies of all reactants, and all products, you can gain valuable information about a chemical reaction between them, such as an estimate of how much energy would be released if the chemical reaction were to occur, or how much extra energy would need to be given to the system in order for the reaction to become likely. DFT can also be used to map out the potential energy surface for a chemical reaction.

However, other aspects of a reaction such as, how must faster it will be at some temperature/pressure compared to another temperature/pressure, may be better modeled using molecular-dynamics. I have given some explanations about how molecular dynamics works in my answers to:

• How to calculate the volumetric energy density of a molecule?,
• How do I simulate the interaction between two atoms?, and again here,
• How would one find a material's equilibrium structure at any specific temperature?.

My answer to the last of those questions also discusses how to do a more accurate chemical dynamics simulation by solving the time-dependent Schroedinger equation (which is in practice quite computationally demanding and often seen as overkill for most applications today).

Answer 2

In this answer I would like to supplement the existing, concise answers with some detailed contexts.

I assume that you are interested in either the thermodynamics or kinetics, or both, of a molecular reaction, and want to derive this information from a computational study. There are a few ways to do this:

1. The most straightforward idea is to put the reactant molecules together and simulate the dynamics of the system frame by frame (typically you would need a frame for every 0.5~2 femtoseconds), until the molecules react; then, from the average time needed for a reaction to occur, one can deduce the rate constant, and if you allow the reaction to go back and forth multiple times, and count how many frames the system is in the product state, you can read out the equilibrium constant. You can also get plenty information about how the reaction occurs, for example which sequence of formation, distortion and breakage of bonds are involved in the reaction. This is called a molecular dynamics (MD) simulation.
2. When the reaction takes longer than maybe a microsecond, a direct MD simulation is not affordable, as the total number of frames needed to cover the desired time scale is too large. In this case you may speed up the simulation by encouraging the molecules to participate in a reaction, say by applying forces to them to bring them together; this will of course lead to a drastic overestimation of the rate constant, but this overestimation can be mathematically corrected after the simulation. Thus you can, for example, simulate a reaction that takes milliseconds by a simulation of merely a nanosecond. Depending on the exact technique used, this may be called a metadynamics simulation, or it may have some other fancy names.
3. Sometimes, even metadynamics is not affordable, either because the required simulation time is still too long, or because the reaction cannot be easily promoted by an artificial force. In this case one can optimize the transition state of the reaction, and calculate the rate constant from the activation Gibbs free energy using the Eyring equation, which is less accurate (maybe introducing an error of about an order of magnitude to the rate constant) but is considerably cheaper, with a computational cost usually comparable to no more than a picosecond of MD simulation. This is called a transition state theory approach, and is probably much more popular than MD in DFT-based studies, but some may refuse to call it a simulation since it does not involve a frame-by-frame simulation of a real process.
4. Of course, when only the equilibrium constant is desired, one does not even need to find the transition state; simply calculating the Gibbs free energies of the reactant(s) and the product(s) is enough. This can be done in two ways, either by a conformational sampling (e.g. using MD or Monte Carlo), or by Taylor expanding the potential energy surface (PES) at the equilibrium geometries (which encompasses the extremely commonplace approach of estimating the Gibbs free energy from a harmonic approximation, as well as approaches that incorporate anharmonic corrections, e.g. VPT2). Not everyone would call these simulations, although if a conformation sampling is involved and it is done by MD, it can usually be called a simulation.

As you can see, the question may be further clarified by specifying whether the word "simulate" encompasses all computational chemistry studies (including transition state theory calculations, where essentially only three "checkpoints" of the reaction are calculated), or only those computational studies where the reaction is simulated frame by frame. Fortunately, the answer to the question does not depend much on this.

If one works under the Born-Oppenheimer approximation and is studying a ground state reaction, then the electrons in the molecules always minimize their energy instantly as soon as the atomic nuclei move to a new set of positions. This means that, in a MD simulation one only needs to simulate the dynamics of the nuclei, which (apart from the dependency on the initial condition and the environment) is completely determined by the adiabatic potential energy surface (PES) of the system. In a transition state theory approach, the activation Gibbs free energy is used instead, and this is also completely determined by the adiabatic PES. As DFT is an exact theory for the ground state electronic energy of a system under any given nuclear configuration, it follows that we can do an exact Born-Oppenheimer MD (BOMD) simulation of an arbitrary reaction, or an exact Gibbs free energy calculation of any species, if we have the exact functional at hand. In this sense the answer to your question is affirmative. Of course we only have access to approximate functionals, which make the simulation approximate as well, sometimes even qualitatively wrong, but you can blame the approximate functionals rather than the DFT theory itself for that.

When the reaction is an excited state reaction, DFT cannot be directly used because it only gives ground state energies. However you can calculate excitation energies using TDDFT (through the response functions of the ground state), which means that you can perform excited state MD simulations and Gibbs free energy calculations using TDDFT. Thus, if you consider a TDDFT simulation as a type of DFT simulation (as some people do), the answer in this case is also affirmative. Again, the simulation is only exact when you have access to the exact functional (and the exact frequency-dependent XC kernel), which you don't have; and if you use approximate functionals and XC kernels, you may sometimes experience qualitative failures, for example the electronic energy may become a complex number at some time during a MD simulation or a geometry optimization, effectively forcing the simulation to abort. Again, you may blame the approximate functionals and XC kernels instead of blaming the DFT theory itself.

When the reaction involves at least two electronic states (e.g. when it involves a spin crossover, or an internal conversion from an excited state, or involves a spontaneous excitation to an excited state), the Born-Oppenheimer approximation fails, and one should perform a non-adiabatic MD (NAMD) simulation instead. A NAMD simulation requires not only the electronic energies of the electronic states, but also the non-adiabatic coupling matrix elements (NACMEs) between the electronic states. It was not until the last decade that it was shown that NACMEs can be calculated exactly by TDDFT (you may have a look at a review on this topic by me and coworkers). Therefore, one can exactly simulate a reaction (at least under the assumption that the nuclei are non-relativistic particles) using TDDFT even if the Born-Oppenheimer approximation does not hold, and this is a non-trivial statement! Again, with approximate functionals you only get approximate results, and this time there is an additional kind of possible qualitative failure, namely the approximate NACMEs may become infinite even when the exact NACMEs are finite. In the review linked above, I mentioned a way of getting rid of this divergence at little loss of accuracy, but the resulting equation is probably no longer exact when one uses the exact functional and XC kernel.

In summary, the conclusion is: (1) DFT can be used to simulate adiabatic, ground-state reactions; (2) DFT can only be used to simulate non-adiabatic and excited-state reactions if you think a TDDFT simulation can be called a DFT simulation; (3) if you ever find the opposite during an actual simulation, this is due to the additional approximations involved, and is not a problem of (TD)DFT itself - that being said, if the problem is universal among all approximate functionals available to you, this may still force you to abandon DFT and use, e.g. wavefunction methods instead.

Answer 3

Can you simulate the chemical reaction between two or more molecules using DFT?

The short answer is no. In essence, DFT is solving the Kohn-Sham equation, which is a time-independent Schrodinger equation and hence can't capture the dynamical processes in chemical reactions.

If not, what numerical method is usually used for this?

The time-dependent density functional theory (TDDFT) maybe can be used for that purpose.