job types
Here we will review different job types that we typically execute for our thermo-kinetic computations. Each job type could be run at a different level of theory.
Before we discuss the levels of theory and the computational chemistry methods themselves, let's first consider what types of jobs we need to properly characterize a species or a reaction transition state (TS).
Consider ethanol (CH3CH2OH), a relatively "simple" species. As users, we demand from ARC that our input will be convenient, for example as SMILES (see the entry on 2D molecular representation). The SMILES for ethanol are CCO. Since the electronic structure calculations are done on a collection of atoms in 3D space, the first task we need to perform (or ARC automatically does for us) is to understand how our species looks like in 3D. But wait! the problem is that ethanol has several 3D stable configurations in space (called "conformers"):
Which of the conformers above is the correct one? How do we define "correct"? The most apparent difference between these conformers is that the "OH" group is rotating - we call this internal rotation and the rotation angle is called a torsional angle. In addition to H number 9 moving around, another difference is in the CH3 (methyl) group - see how the indices sometimes change place. They are always in the same order, though. This is another internal rotor, but since its symmetric (a symmetry number of 3), we don't relate to it much when considering conformers.
We should consider that some of the above conformers will have a different potential energies since the forces between the atoms are different (but some will not - can you see which ones?) It's likely that most of the species population will be most of the time in the lowest potential energy state, and once in a while overcome the energetic barrier for internal rotation. So the conformer we care about the most is the lowest energy one. That's definitely not the entire story, since we need to somehow consider the conformers the species is in for less amount of the time. More on that later.
Typically, the first computational chemistry "jobs" we'll need to run are optimizing the 3D geometry of different conformers and computing their respective energies. Below we discuss these job types (refereed to as "opt" and "sp").
A single-point (sp) energy calculation is "simple" in essence: it computes the potential energy (repulsion and attraction energies between all particles in a molecule) for some arrangement of atoms in space. The underlying molecule could be at any state, i.e., an energy well, an energy peak, a saddle point, or any other point. Many of the following job types rely on performing many single-point energy calculations (opt, scan, IRC).
The energy computed using the sp job is called E0.
A geometry optimization job uses many sp jobs and an optimization algorithm to identify the localized minimal energy well on the potential energy surface (PES). Note that the optimization could also be constrained if desired: A user may request to perform a constrained optimization job in which some degrees of freedom are kept constant (constrained), and all other degrees of freedom are allowed to "relax" and be optimized. For chemical species, we will optimize by minimization of the energy to obtain the local minimum. For transition states (TSs) we optimize to a saddle point, i.e., the optimization algorithm will search for a minimum energy is all modes of motion but one (the reaction coordinate) which should be a maximum. See some more discussion about TS below under Frequencies. When optimizing a TS, we will have to tell the electronic structure software that this is a TS by passing it a certain keyword.
You can think of a frequencies job as identifying the frequencies of masses connected through springs as in classical mechanics, just to get an intuition. This just gives us some intuition, of course the computation we're actually doing utilizes quantum chemical principles and is not classic at all. By computing frequencies we can understand the molecular vibrations.
Each molecule has 3N degrees of freedom in 3D space (where N is the number of atoms). However, here we care about internal degrees of freedom, so we exclude translation in X, Y, Z of the molecule as a whole, as well as rigid rotation in 3 principal directions. So overall, a non-linear multi-atom molecule has 3N - 6 internal degrees of freedom, and will therefore have 3N - 6 frequencies, each corresponds to a different mode of motion. Note that linear molecules have 3N - 5 internal degrees of freedom, and there is one less external degree of freedom which is the rotation around the main axis. Trivially, molecules that consist of a single atom have no internal degrees of freedom and no frequencies. The number of frequencies a molecule has corresponds to the number of degrees of freedom it has.
Why do we care about computing frequencies for a molecule? Because even at 0K (absolute zero temperature) molecules are not still and they vibrate. We need to know the vibration frequency of each mode of motion to know how much energy the molecule has at 0K above the E0 value computed using the sp job type. In other words, if you molecule has more than one atom, computing sp only is not enough to know its energy. Another reason to compute frequencies is that we want to know the heat capacity of the molecule, and how much energy can be absorbed by each mode of motion.
The Morse potential for diatomic molecule is schematically drawn in the figure above for HCl. This motion mode is not harmonic, but we model it using an harmonic oscillator approximation, justified by the fact the the lower part of the potential does resemble a parabola. We call this model the rigid-rotor harmonic-oscillator (RRHO). As the system temperature increases, the molecule will have a higher kinetic energy which means it has a higher velocity and will rotate around its axis faster; moreover, at higher temperatures the molecule will have a higher vibrational quantum number v. Note that a molecule cannot be in a state where it has no vibrations at all - it will always vibrate in each of its motion modes, even if it's the smallest vibration allowed (the smallest vibrational quantum number, v=0). This means that computing the sp energy (E0) is not enough to determine the energy of a species, since each vibrational mode has some amount of additional energy. The sum of these energy contributions is called the zero-point energy correction, and the total energy of the molecule is then called the zero-point energy (ZPE). So we can write:
E0 + zero-point energy correction = ZPE
The electronic structure software (e.g., Gaussian, QChem, Psi4, etc.) will compute a Hessian matrix, which is the matrix of second order partial derivatives of the energy with respect to the position of the atoms. The Hessian matrix is sometimes called the "force matrix" and is often used in geometry optimization (energy minimization). The eigenvalues from the Hessian matrix are then converted to the vibrational frequencies. We know that all second derivatives must be positive for a minimum (well), corresponding to positive eigenvalues. At a saddle point (a TS), one eigenvalue is negative, giving an imaginary frequency (a complex number). Such saddle points with exactly one imaginary frequencies represent reaction transition states.
Each frequency has a "normal mode", that is an X-Y-Z direction in space for each atom in the molecule along which this particular vibration causes motions. We call this the "displacement normal mode", and this is how we can visuallize the vibrations.
As noted above under "conformers", we need to consider that some population of the species is at a different conformation than the lowest energy one. We will account for that by performing scans of anharmonic modes and considering these scans in our statistical mechanics processing (in the thermodynamic molecular partition function). What are anharmonic modes? Any internal molecular mode of motion that is not well described by an harmonic oscillator (does not behave like a parabola). These modes include internal rotations (torsion or dihedral angles), ring puckering, and inversions.
Mathematically, a dihedral angle (or a torsion) is defined as the angle between two planes. Three points in space which are not on a straight line form a plane. If we have four atoms, A-B-C-D, atoms A-B-C form plane 1, and atoms B-C-D form plane 2 (if they’re not on a straight line). Thus these four atoms define a dihedral angle. This pure mathematical point of view is not always intuitive, and it’s easier to think of a torsion in a molecule as two “handles” that one can grab in both hands and turn them in different directions. For example, consider atoms 9, 1, 2, 3 in the ethanol conformers depicted above. To attain conformers number 2 and 3 starting from conformer number 1 (counting from left), one needs to hold the bond between 9-1 and the bond between 2-3 to keep the rest of the molecule in place, and rotate them relative to each other. The actual motion is around bond 1-2 in the middle. This describes the “OH” internal rotor. The forth ethanol conformer above is attained by rotating the CH3 internal rotor (e.g., atoms 1-2-3-6). See the animation below of a rotating internal CH3 rotor (it is shown in a transition state, but that’s not important for now).
Below are examples for internal rotations in a species. This species is a transition state (can you spot the reactive region?), and we can see two different methyl (-CH3) internal rotation modes. Every single-bond could potentially be considered as an axis for a hindered rotation motion, we only present two here.
Next, we’d like to plot the potential energy of the species vs. the torsional dihedral angle. Below are several examples of such internal rotation scans. We call these 1D scans, since only one torsion was modifies. Here are examples of two separate 1D scans for ethylamine (SMILES CCN):
These graphs shoe the potential energy (in kcal/mol) vs. the dihedral angle in radians (0-360 degrees). The scan on the left shows the NH2 rotor scan, and the scan on the right shows the CH3 rotor scan (how can we tell that? look at the symmetry). The black symbols are the actual calculations. The red and blue curves are fits using cosines or a Fourier series, respectively. See how the cosine fit fails to capture the middle peak on the left (the NH2 rotor) - the cosine fit only works for symmetric rotors. In these figures we used a resolution of 10 degrees (36 points), which is widely practiced in the literature. However, it is better to use a higher resolution of 8 degrees (45 points) to assist the algorithms in automatically identifying the symmetries of the peaks (although it a more expensive in terms of the computations needed).
Sometimes (often?) two or more internal rotors are “conjugated”, i.e., there are significant internal molecular forces between parts of the molecule (usually attraction forces such as hydrogen bonds) that cause another internal rotor to rotate along with the one we scan. When rotors are conjugated, we need to scan them together, e.g., in 2D. ARC can do these scans for us, but they are currently not taken into account when computing thermodynamic properties (work in progress). Here's an example for a 2D scan of ethanol (CH3CH2OH):
In this 2D scan, the X axis represents the CH3 rotor, and the Y axis represents the OH rotor. See how the data on the X axis repeats itself three times: indeed the CH3 rotor has an internal symmetry number of 3.
We call these internal rotors "hindered rotors", meaning that there is some barrier for rotation, and the motion could be back and forth within a certain energy well, and once in a while the rotor could accumulate enough energy to pass the barrier. There's another kind of rotor called "free rotor" where the barrier for rotation is negligible.
Intrinsic reaction coordinate (IRC) is a calculation that slightly perturbes a TS in both directions of the displacement normal mode corresponding to the imaginary frequency and optimizes each of the two resulting geometries into wells. We then get the result of the lowest energy path that connects the TS to the reactants and products well. This method allows us to make sure that our TS geometry really corresponds to the reaction we care about. Here are two visualizations of TSs:
This is a check that verifies whether the electronic structure is stable. This calculation will try to find whether an electron can be populated in a lower molecular orbital. If it finds one, this means that our original electronic structure wasn't the ground state we usually try to capture, but rather an electronically excited state.
A composite method is a recipe that specifies which level of theory to use per job type. Only the Gaussian software supports composite methods. For example, a common and well-known composite method is CBS-QB3, you can read about it here and here.
Specifying CBS-QB3 means the following:
- B3LYP/6-311G(2d,d,p) geometry optimization
- MP22(FC)/6-31G geometry optimization
- UMP2/6-311+G(3d2f,2df,2p) energy and complete basis set (CBS) extrapolation
- MP4(SDQ)/6-31+G(d(f),p) sp energy
- QCISD(T)/6-31+G sp energy
We try to explain what these levels mean under methods and levels of theory.