Quantum chemical methods can be used to calculate heats of reactions. The Heat of formation (HOF) is a special type of reaction.
For example, to calculate the bond separation energy of
Experience has shown that the quality of the calculated heats of reaction depends on both the method used to calculate the energy (E[..]) and the type of reaction.
Our guide book explains that the "best" type of reaction is an isodesmic reaction. (An isodesmic reaction is one in which the number and type of bonds remain the same during the reaction. For an example see the discussion on basicities.) The book also lists other types of reactions and gives quantitative examples of the accuracy one can expect for different theories and basis sets.
In principle the heat of formation (HOF) is calculated like any other reaction, see the discussion above. Sadly, the experimental HOF consists of some of the most difficult reactions to calculate: breaking all bonds and creating atomically pure compounds such as H2 for hydrogen and graphite for carbon.
Thus, good HOF calculations usually require methods better than QCISD with very large basis sets. These can be extremely time-consuming, even for molecules of only a few heavy atoms! The G3 method is a good place to start if you are interested in accurate Heats of Formation.
Chances are that if you think you want the HOF formation, you may only need to compare the relative differences in energy of certain reactions, and can use a less time-consuming approach than the methods required for Heat of formation calculations. The guide book includes many examples of HF/6-31G* without vibrational corrections performing well (against experiment). The discussion on finding the *best* energy also covers this and other common theory levels.
The G3 method (the most used of the Gx series) is a time consuming recipe for calculating accurate "Heats of Formation" (HOF) from first principles. G3 consists of 10 calculations:
All these energies are combined, along with similar calculations done for each atomic species. The result is then subtracted by experimental energies of the "standard state atomically pure" systems to arrive at a good approximation for the actual HOF. A quick overview of how these steps are used can be found in our discussion of calculating the exact energy.
Note that each of these steps is relatively time consuming and practical only for very small systems. The G3(MP2) method is often more than twice as fast, but still prohibitively time consuming for all but small molecules.
The reference to the G3 papers can be found at the end of this FAQ.
The T1 method is a simpler recipe to calculate the absolute heat of formation. It is similar in concept to the G3 method but replaces the post MP2 calculations with empirical relations based on Mulliken bond order and implements a faster MP2 approximation. Thus it is much faster than G3 or G3(MP2) and can be applied to much larger molecules. The steps of the calculation are:
The quantum mechanical calculations in Spartan assume that molecules are isolated, with a temperature of zero Kelvin, and with stationary nuclei. Real experiments are carried out with vibrating molecules at finite temperature (often 298.15 K). In order to correct for these differences Spartan can use calculated frequency data to determine a set of normal-mode vibrational frequencies (v_{i}). (These are the same data used to generate IR spectra.)
Using this approximation, energy corrections are fairly straightforward. Splitting the Enthalpy correction into four (4) parts to include the 'Zero point energy' (Zp), the temperature correction (Hv), enthalpy due to translation (Ht) , and enthalpy due to rotation (Hr).
Note that for small molecules this is a very good approximation. However as molecules get larger the 'normal-mode' approximation becomes less valid. For example:
These problems are magnified if one is examining the entropy or free energies of larger molecules. For these reasons low frequency vibrations are often systematically ignored.
The entropy in Spartan is calculated from the same data as vibrational enthalpy (Hv). The equation for entropy (S) can be written as follows:
The constants u_{i} are proportional to the
vibrational frequencies and v_{A} are inversely
proportional to the moment of inertia of the rigid body. It is
important to note that these equations ignore any entropy due to
conformational flexibility, which may dominate for
larger molecules. Also note that low frequency vibrations will
dominate the calculation of entropy. It is these
modes which break the normal mode assumption, thus making
entropy calculations of larger molecules suspect.
Return to Top
You can use the isodesmic proton transfer reaction relative to a standard to calculate the relative basicities. For example if we choose the ammonia molecule as a standard we would write:
where B is the compound of interest. Then the relative proton affinity is:
Relative proton affinity is a measure of the basicity which can easily be converted to pH. In our methods book (table 6-17) we show the accuracy of this calculation for a number of methods and basis sets for series of Nitrogen Bases. A regression analysis shows that even simple HF/6-31G* methods provide results with an accuracy of less than 4 kcal/mol on a system covering 125 kcal/mol range.
There is also an interesting example of predicting pKa's from the graphical map of the electrostatic potential. This can be found on page 478 of A Guide to Molecular Mechanics and Quantum Chemical Calculations.
The pH of a compound should be related the energy (enthalpy) of the reaction described above
If you know the pH of a few standards you can use these in a regression fit (using regression tools in Spartan's spreadsheet). This will result in a simple linear relation correlating 'energy' with 'pH'.
While pK_{a}'s can be calculated using methods described above for pH and basicity it is often possible for to use empirical relationships. For example the Groups of Organic Molecules tutorial (menu > Tutorials) has an example of determining the pK_{a} of carboxylic acids from the electrostatic potential of the acidic hydrogen. With a simple linear regression one can get rather good prediction for a rather limited case of molecules. As is shown in the first plot below which is the result of going through the above mentioned tutorial.
For many cases empirical rules such as this work "good enough" and are so quick to calculate that this is the recommended approach. In principle we calculate the pK_{a} from first principles using the deceptively simple equation
However, using the fact that the change in energy should be proportional to the energy of de-protonation one can do some quick calculations to approximate pK_{a}'s. The second plot in figure above is a linear fit of the same carboxylic acids but calculate the energy difference between the neutral acid, and the de-protonated acid. In this example the geometry used was the same as the neutral case, except the acid hydrogen was removed. In some cases this kind of fit might be "good enough". The third fit is a similar calculation but done in solvent.
Science is a series of approximations. The trick is to use the largest simplification one can and still get reasonably 'good' predictions. This said, we always feel more comfortable with a theory if a recipe exists to get the *exact* answer. In QM calculations this is so prohibitively time consuming as to be impossible, but it can be described. We describe this 'exact' calculation in the next section.
One definition of the *best* energy is the energy which answers one's questions with enough accuracy, yet consumes the minimum amount of time and resources. As mentioned in the discussion on calculating heats of reactions, this is often achieved by writing a near-isodesmic reaction and calculating an energy difference. (The key here is the 'energy difference'. It allows cancellation of systematic errors.)
In order to find what is *best* for your question, we suggest starting with a question/reaction in which you know the answer and that is similar to your reaction of interest. Try this reaction with HF/6-31G* or HF/3-21G*. If that isn't accurate enough then add the zero-point energy (ZPE) corrections. (scaled by .90 for HF or by .99 for DFT) If that fails, try MP2 or DFT energies with larger basis sets. If that fails attempt single point CCSD with ZPE corrections. Finally try the T1, G3(MP2), or G3 methods.
Conversely, if HF/6-31G* is good enough, you may want to try faster methods to see if they are also acceptable. HF/3-21G* is a popular basis set and much faster than 6-31G*. Semi-Empirical is the simplest method using quantum mechanics approaches while requiring an order of magnitude fewer resources and time. Molecular mechanics methods, while limited to ground states of common chemistry groups, is extremely quick and in the areas in which it is parameterized can be fairly accurate.
If you are new to molecular modeling,
Wavefunction's guide book,
discusses many types of reactions and molecular properties,
quantitatively comparing many different theory levels.
Return to Top
The complete Enthalpy (H) and the Gibb's Free Energy (G) of a given conformer of a molecule in the dilute gas phase can be written
Calculating all of these terms is impossible for the foreseeable future for any real/non-trivial molecule. However a number of methods and recipes have been proposed that attempt to approach this limit. Below we describe such an algorithm using the G3 method. The G3 method is a specific thermochemical recipe available in Spartan that attempts to produce accurate results in a relatively short amount of time. By studying it we can generalize it to imagine a method which should give the *exact* results even if such a recipe is impractical.
Below is the edited result of a G3 run on a water molecule, used as a reference:
Energy HF/6-31G(d) -76.009809 1 Energy MP2/6-31G(d) -76.196848 2 Energy MP4/6-31G(d) -76.207326 3 perturbation correction -0.000565 4 polarization correction -0.074488 5 diffuse correction -0.012979 6 basis set correction -0.081653 7 spin-orbit correction 0.000000 8 electron pair HLC -0.025544 9 ---------------------- ------------------ G3 energy (Ee) -76.402556 HF/6-31G(d) zero point 0.020515 11 G3 energy (Eo) -76.382040 HF/6-31G(d) Vibrational 0.020518 12 HF/6-31G(d) Translation 0.001416 13 HF/6-31G(d) Rotational 0.001416 RT constant 0.000944 G3 kT energy (H298) -76.378260 Atomic Energy -76.032990 14 dH(298) -0.091381 dH(298) -57.34 kCal/mol
A description of successive improvements to the energy is provided. The details of each step are not discussed here, as the main point is to provide an overview of what is required for an *exact* answer.
HF methods ignore electron correlation energy. The goal of DFT and other post HF methods is to account for this electron correlation energy.
One weakness in the approach described so far is the basis set. There are a number of possible improvements that can be made to the basis set:
If we carry all the above corrections correctly to their infinite limit we should have the *exact* electronic energy. We label this energy Ee. But this would still not be good enough. We've made the assumption that the nuclei do not move. Quantum mechanics tells us that everything is moving, even at absolute zero; the so called "zero point" energy (Zo). By doing a normal mode analysis we can calculate this zero point energy.
The above energy would be correct at zero Kelvin, but chances are that the system we are studying exists at a non-zero temperature. We can make normal mode corrections to the energy, which we now call enthalpy and label it H298. (298.15 K being the most common temperature of interest.)
At this point we notice that the energy is in some "strange" units called Hartrees. In order to compare with other results it is useful to do a conversion.
This final conversion step is not always necessary. What is necessary in QM calculations is to ensure that when comparing energies, the energies are not only in the same units, but define the same zero. (Typically that means the same basis set and theory level, and molecules with the exact same constituents and number of electrons.)
Often, when examining reactions we require even more; the Gibbs free energy (G) which is calculated by subtracting entropy: G = H - TS + PV.
The above steps should cover all the physics of a completely isolated molecule. (i.e. In a vacuum or a dilute gas.) Adding a real environment (such as a solvent or a crystal lattice) is vastly more difficult.
Combining all of these corrections terms "correctly" (for which there is currently no 'good' systematic procedure) one should arrive at the "exact" answer.
Thankfully, a careful examination of the types of chemistry questions one typically asks shows that one does not need to do all this work. (See our guide book for a thorough discussion of this.) Surprisingly, in many cases, HF/6-31G(d) is a good approximation to reality. This happy accident makes quantum chemistry possible on today's computers.
Spartan gives users the ability to choose from a number of advanced correlated methods, and implements the G3 recipe, which is a standard approach for improving the accuracy of calculated energies.
Spartan also provides the flexibility to work with a wide variety of theoretical models:.
All three approaches have a different definition of zero, and use different units (by convention).
Given this definition one cannot compare molecular mechanics energies of different molecules. This restriction holds true even for different isomers. One can compare energies of different conformations. In fact this is one of the most popular and appropriate uses of molecular mechanics.
Semi-Empirical methods use the same definition as ab initio, but a constant is subtracted for each atom in the molecule. The goal of this 'shift' is to approximate a heat of formation (HOF). In reality, this approximation is not very accurate and should only be used as a first guess at the HOF. Semi-Empirical methods are most useful for predicting rough geometries and transition states quickly.
Ab initio methods (including HF and DFT) compare the energy of the combined molecule with the energy required to remove every electron, and separate every nuclei. While this reaction is not a very realistic experiment, it is theoretically simple and can be used to calculate energy differences of conformers, isomers, and entirely different molecules. See the discussion on calculating the heat of formation for a more detailed description of using this energy.
By convention ab initio results are returned in units of Hartrees. See the units section to convert from Hartrees to more common units
Each method is different, but an overly short answer would be:
For an exact number one will need to examine the primary literature on these methods. However, slightly more detailed discussion of these methods presented below. More details can be found in Spartan's documentation.
Spartan can access a large number of DFT functionals. Only the most popular and well tested are available from the "Density Functional" pull-down menu of the Calculations dialogue. For example Truhlar group's M06-L function while not available by default from the menu can be accessed by clicking on "More..." at the bottom of the menu list. This leads to a hierarchical dialogue offering a wide range of functionals organized by classes. Additionally, specific functionals may be typed in the "Options" line of the Calculations dialogue. When the Enter/Return key is struck the specified functional will disappear from the Options line and will appear as the selected functional in the "Density Functional" pull-down menu.
It is possible to define your own functionals with the EXCHANGE= and CORRELATION= keywords. Each functional must be followed a colon (':') and then its weighting. Different functionals are then separated by commas. For example a functional consisting of 9 parts "Slater" exchange and 1 part "Becke" exchange would look like:
If any correlation functionals are required one would use the CORRELATION keyword. For example:
The values allowed exchange functionals are
For those who are familiar with the Q-Chem Exchange/Correlation notation; using EXCHANGE or CORRELATION without any colons sets the appropriate "REM" values. With colons included, an "XC_FUNCTIONAL" section is created.
We implement the main equation of the Vosko, Wilk, and Nusair paper. Sometimes, people will use equation #5 of this paper, which, in Spartan is called B3LYP5.
In some references the coefficients are written differently. For example we write the B3LYP functional as
Another subtle difference between different B3LYP calculations can be the selection of the numerical grid used to calculate the density. See "What size grid does Spartan use for DFT calculations?" for a discussion of this.
Recently there has been some large scale coverage of potential problems with grid sizes in DFT calculations. For example a paper by Bootsma and Wheeler and mentioned in a C&EN article.
We think our defaults are good for most applications, but it is "easy" to check by re-running with the HUGEGRID keyword on sample "energy differences" (such as the two end-points in a reaction) to confirm.
We determine default grid sizes based on what give consistent gradients, as opposed to "best energy", but regardless the grid error is nearly always less than the "incomplete basis set" error, both of which are presumably random for reaction energies (and other "relative" values, ie. not atomization energies)
Meta functionals like M06* and wB97XM-V do require bigger grids sets than wB97X-D and other "2nd generation" sets. (even though we use (70,302) similar to wB97X-D's default to follow Truhlar's paper)
wB97X-D which uses (75,302) by default requires more than B3LYP for which we use the SG-1/SG-0 grids (which themselves were optimized for the Pople basis sets 6-31G and 6-311G).
We contend that "grid errors" should be given as "error per atom" dEg/a (read as delta Energy-grid per atom). Assuming that this error it is not systematic, the "total" error/noise would then be [dEg/a] x [N^{1/2}]
The curious might want to examine the discussion of the grids we are most familiar with (SG1 and SG0) which we use for B3LYP and EDF2. Most modern papers covering "new" functionals will spend time discussing the preferred grid.
Of course in Spartan you can easily change the default gradi as discussed below.
For the "simpler" functionals, such as B3LYP and EDF2, Spartan employs a hybrid SG-0 and SG-1 grid as follows: SG-0 for H,N,C, and O, SG-1 for all other atoms. For the more "modern" functionals we use a larger grid (because they are more sensitive to small changes in grid). Specifically, we use the "BIGGRID" (70,302) for the Truhlar functionals and Head-Gordon functionals such as:
Following are some keywords which can be typed into the Options line of Spartan's Calculations dialogue to change the grid:
In the above notation the first number is the number of shells in the radial direction, the second number (i.e. 194, 434 etc.) is the number of Lebedev radial points.
One can also use the BIGGIRD keyword with an equal sign, to enter a Q-Chem like grid notation. Specifically a 12 digit number with the first six (counting leading/implied zeros) defining the number of shells in the radial direction and the next 6 defining the number of Lebedev radial points. i.e., BIGGRID and BIGGRID=70000302 both refer to (70,302)
Valid values for Lebedev grids are:
6, 18, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302,
350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702,
3074, 3470, 3890, 4334, 4802, 5294.
If you want to use Gauss-Legendre angular points (a=2N^2) instead of Lebedev numbers use BIGGRID=-rrraaaaaa (i.e. use a negative number).