Skip to content

Step 1 — Building and Optimizing the Reactants

Basic tutorial

This step assumes familiarity with the Geometry Optimization + Frequency Calculation tutorial. Only the parts specific to this reaction are covered here.

Building the initial structures

1,3-butadiene and ethene were built in their s-cis / planar minima and exported as Cartesian coordinates. Since the two molecules do not react until they approach each other, they were placed in the same input file, far apart and without a connectivity bond between them:

%nprocshared=1
%mem=2GB
%chk=diels-alder-opt-freq.chk
# opt freq b3lyp/6-31g geom=connectivity

Title Card Required

0 1
 C   -0.05107701   2.62262333  -2.41495425
 C    0.34811294   3.26894028  -3.53722374
 ...(1,3-butadiene, 10 atoms)...
 C    1.63529934   0.78409688   1.86498813
 C    2.98818596   0.77382567   1.94346957
 ...(ethene, 6 atoms)...

 1 2 2.0 9 1.0 10 1.0
 2 3 1.0 8 1.0
 3 4 2.0 7 1.0
 ...
 11 12 2.0 13 1.0 14 1.0
 12 15 1.0 16 1.0

The idea was to obtain both optimized reactant fragments — and their combined SCF energy — from a single job, so the resulting energy could be compared directly against the transition state and product without bookkeeping two separate output files.

What actually happened

The job did not converge. The Gaussian log shows the Berny optimizer oscillating for the full 86-step limit before stopping:

 Step number  85 out of a maximum of   86
 Maximum Force            0.008567     0.000450     NO

 Step number  86 out of a maximum of   86
 Maximum Force            0.022472     0.000450     NO

 Optimization stopped.
    -- Number of steps exceeded,  NStep=  86
 Error termination request processed by link 9999.

The Maximum Force column alternates between roughly 0.0086 and 0.0225 Hartree/Bohr for the last several steps instead of decreasing — a classic sign that the optimizer is stuck oscillating between two nearby geometries rather than converging.

Common pitfall — optimizing non-bonded fragments together

When two molecules are placed in the same input file with no bond between them, the redundant internal coordinate scheme that the Berny optimizer builds has to describe the relative translation and rotation of one fragment with respect to the other using Cartesian-like coordinates. Because the fragments are far enough apart that they barely interact electronically, the gradient along these "inter-fragment" coordinates is extremely small and numerically noisy. The optimizer has almost nothing to converge against on that mode, so it can oscillate indefinitely instead of settling down — exactly what happened here.

The fix: optimize each reactant in its own, independent job (its own .com file, its own opt freq calculation). This is also what the DFT workflow means by "Optimize Reactant" — one species per job. Combine the two converged electronic energies (E(diene) + E(dienophile)) afterward for the energetics table in Step 4.

Using the result anyway (with a caveat)

Because the optimization never reached a stationary point, the frequency job that was meant to follow it never ran either — there is no zero-point energy or thermal correction available for this attempt. The last few SCF energies from the oscillating tail of the job are:

Step E(RB3LYP) / Hartree
84 −234.525788496
85 −234.526136060
86 −234.525788718

These values are close to each other (within 0.0004 Hartree, or about 0.2 kcal/mol) even though the geometry is still oscillating, which suggests the two fragments are already near their individual equilibrium geometries — the residual instability is coming from the poorly defined inter-fragment coordinate, not from real chemistry. For the purpose of this worked example, −234.5261 Hartree is carried forward in Step 4 as an illustrative, non-rigorous reactant energy, clearly flagged as such.

What to do for a real project

Rerun butadiene and ethene as two separate opt freq jobs at the same level of theory, confirm each has zero imaginary frequencies, and use the sum of their converged electronic (and, if needed, thermally corrected) energies as the reactant reference. This is the number that should ultimately be used in place of the illustrative value above.


← Back to Overview · Continue to Step 2 — Transition State Search →