Specifying a temperature in a DRC

The behavior of a system at a specified temperature can be modeled by use of KINETIC.  This keyword allows extra kinetic energy to be added to the system.  In order to determine how much extra energy to add, an understanding of the issues involved is essential.

Energy issues involved in a dynamic system

The total energy of a system can be expressed as three sums:

(A) The heat of formation of the system, ΔH_f0.  This is an irreducible minimum, and represents the energy of the system at equilibrium.

(B) The potential energy of the system.  This is the heat of formation of the system with its current geometry, ΔH_f, minus ΔH_f0. It represents the energy of distortion from the equilibrium geometry.

(C) The kinetic energy of the system.  This is the sum of the vibrational energies of motion of all the atoms in the system.

Energy term (A) is a constant, regardless of temperature. Energy terms (B) plus (C) represent the internal energy (enthalpy), U, of the system.  At zero Kelvin, U is zero. At any other temperature, T, the enthalpy represents the integral of the heat capacity from zero to T. In a DRC calculation, any desired temperature can be specified by defining the associated internal energy.  The internal energy can be calculated using THERMO, and specifying the temperature to be used.  In the output of a THERMO calculation, the enthalpy needed is given at the intersection of ENTHALPY and VIB.  In the following example, this would be 578.5 cal/mol, or 0.578 kcal/mol.

    TEMP. (K)   PARTITION FUNCTION   H.O.F.    ENTHALPY   HEAT CAPACITY  ENTROPY
                                    KCAL/MOL   CAL/MOLE    CAL/K/MOL   CAL/K/MOL


 298.00  VIB.     0.2342D+01                   578.4823     4.2357     3.6320
         ROT.     0.3897D+04                   888.2813     2.9808    19.4109
         INT.     0.9125D+04                  1466.7636     7.2165    23.0429
         TRA.     0.1515D+27                  1480.4688     4.9680    36.0322
         TOT.                        10.536   2947.2324    12.1845    59.0751

To set up a run at a given temperature, the  ΔH_f of the optimized system is needed.  Calculate this first, then do a THERMO calculation to get the enthalpy at the desired temperature.  The next step might be unexpected.  Distort the geometry of the system slightly, and re-calculate the heat of formation at the distorted geometry.  Make sure that it has increased by at least 0.2 kcal/mol, and preferably by a large fraction of the enthalpy. A non-equilibrium starting geometry is needed in a DRC calculation because otherwise the atoms would not be moving, and adding in excess kinetic energy would not be possible (You can't make the atoms move faster if they were not moving originally).

Work out how much extra kinetic energy would need to be added to equal the vibrational enthalpy at the desired temperature.

Let the calculated equilibrium heat of formation be -100.000 kcal/mol.
Let the desired internal vibrational energy be 0.578 kcal/mol.
Let the heat of formation of the distorted geometry be -99.500 kcal/mol.

Then the extra kinetic energy needed would be 0.078 kcal/mol.

Set up the DRC calculation using the distorted geometry and KINETIC=0.078.  Run this system.  In the output, the starting heat of formation will be that of the distorted system, and the kinetic energy will be zero or almost zero.  Over the next few femtoseconds, the heat of formation will become more negative, and the kinetic energy will rise.  When it exceeds 0.2 kcal/mol, the velocity vector for the system will be well defined, and the extra kinetic energy will then be added.  This will temporarily confuse the predictor-corrector error reduction function, so errors in the calculation will rise for a while, but once they become small the total energy of the system (heat of formation of the current geometry plus kinetic energy) will be equal to ΔH_f0  plus the desired enthalpy.

A legitimate question is, why do it this way, why not add the velocity vector "by hand"?  The problem is, the kinetic energy must not include rotational or translational terms. Translational contributions to the enthalpy are irrelevant for a single molecule, and the rotational contributions are conserved, so they, too, can be ignored (or can they?), all that matters is the dynamics of the non-translating, non-rotating system.  Defining a velocity vector that achieves this is definitely non-trivial.