Frequencies¶

Analytical Frequencies¶

The frequencies are calculated analytically by specifying the AnalyticalFreq block keyword (see below for more details). The analytical frequencies are as accurate as the numerical frequencies for the same integration accuracy, but can be up to 3 to 5 times quicker to compute, depending on the molecule, integration grid parameters, and choice of basis set. The analytical frequencies are fully parallelized and linearly scaled.

Calculating the analytical frequencies requires the solution of the Coupled Perturbed Kohn-Sham (CPKS) equations, which is an iterative process. This part of the process is of order 3 x number of atoms, and is generally the main bottle neck in calculating the frequencies. (The immediate result of the solution of the CPKS equations is the U1 matrix, the components of which are closely related to the derivatives of the MO coefficients. One of the adjustable parameters in the input of an analytical frequencies calculation can be used to control the accuracy of the U1 matrix components.)

One disadvantage in calculating analytical frequencies is that the range of exchange-correlation functionals is limited. (This is because derivative formulas have to be derived for each exchange-correlation functional in ADF, which is not an straight forward task). Here are the currently available functionals:

LDA: XONLY, VWN, STOLL, PW92

Exchange GGA: Becke88, OPTx, PBEx, rPBEx, revPBEx

Correlation GGA: LYP, Perdew86, PBEc

XC GGA shortcuts: BP86, PBE, RPBE, revPBE, BLYP, OLYP, OPBE

Any functional not mentioned above is not implemented, including PW91 and Hartree-Fock.

To calculate the frequencies analytically, include the block key word ANALYTICALFREQ. Subkeys are available, but in general, to calculate the frequencies, no subkeys are required, and including the following in your run file is sufficient:

AnalyticalFreq
End

Unlike the numerical frequencies, the analytical frequencies can be computed immediately after a geometry optimization by including both block keywords in the same input file:

Geometry
  ... Geometry optimization options here ...
End

AnalyticalFreq
End

A note of caution: For accurate frequencies it is especially important to also have an accurately optimized geometry. During a geometry optimization the integration accuracy is set by default to “Normal”, and so the resulting frequencies will also have this level of integration accuracy while it may be desirable to have frequencies computed with a higher accuracy. One might consider using Good NumericalQuality (or BeckeGrid quality) and set the convergence criteria for the geometry optimization tighter.

The format for the AnalyticalFreq block key is:

AnalyticalFreq
    PRINT {eigs} {u1} {parts} {raw_freq}
    DEBUG {fit} {hessian} {b1} {densities} {numbers} {symmetry} {all}
    MAX_CPKS_ITERATIONS Niter
    CHECK_CPKS_FROM_ITERATION  N
    U1_ACCURACY  x
    NUC N1 N2 ... Nk
End

An explanation of the subkeys follow.

PRINT

This is primarily for debugging purposes. Choosing EIGS results in the print out of the MO eigenvectors, while U1 results in the print out of the U1 matrices. Except for small molecules this will result in a lot of data being output, and so they are not recommended. Choosing PARTS results in the print out of various sub-hessians that add up to give the final analytical hessian. RAW_FREQ gives the eigenvalues of the initial force matrix, which are essentially the frequencies before rotational and translational degrees of freedom have been removed from the force matrix.

DEBUG

This is for debugging purposes. The choice FIT results in the print out of information related to the calculation of the fit coefficients and their derivatives. HESSIAN results in the printing out of many of the sub-Hessians that add ups to give the final Hessian. Many more Hessians are printed out with this option that with the print parts subkey option (mentioned above). Choosing B1 gives data related to the frozen core orthogonalization coefficients and their derivatives. DENSITIES gives the integrals and moments of various densities computed during the calculation of the frequencies. Including NUMBERS results in the print out of numbers of basis functions, fit functions etc, as well as various integer arrays that are crucial to the calculation of the analytical second derivatives. SYMMETRY results in symmetry information and symmetry matrices being printed out. ALL can be used to print out all debug information.

MAX_CPKS_ITERATIONS  Niter

Calculating the analytical frequencies requires the solution of the Coupled Perturbed Kohn-Sham (CPKS) equations, which is an iterative process. For most systems tested so far, convergence to the required accuracy in the U1 matrix is achieved within Niter=20 iterations, which is the default. If convergence is not achieved (a warning will be printed in the output if this is the case) then this subkey can be used to increase the number of iterations, although convergence is not guaranteed. The user required accuracy of the U1 matrix, as well as the ADF integration accuracy, can effect the rates of convergence.

CHECK_CPKS_FROM_ITERATION N

Solution of the CPKS equations is an iterative process, and convergence is achieved if the difference between U1 matrix of successive iterations falls below a certain threshold. This key can be used to determine at which iteration the checking should start taking place. The default is 1.

U1_ACCURACY x

Solution of the CPKS equations is an iterative process, and convergence is achieved if the difference between U1 matrix of successive iterations falls below a certain threshold. This subkey can be used to set the threshold. The accuracy of the U1 will be 10**(-x). So, the higher the number the more accurate the U1 will be. The default is 4. While this parameter effects the accuracy of the frequencies, other factors also effect the accuracy of the frequencies, especially the ADF integration accuracy.

NUC N1 N2 ... Nk

By default, when calculating the frequencies analytically, the derivatives of the energy with respect to all nuclei are calculated. This gives a complete Hessian (second derivative) matrix, from which the vibrational frequencies of the molecule can be calculated. However, there may be certain cases where only derivatives with respect to a subset of all the nuclei are required. In this case it is a considerable saving in time if only a partial Hessian is calculated. With this subkey, a list of the nuclei for which the derivatives are required can be specified. However, the frequencies in this case are not the vibrational frequencies of the molecule, but may be used in guiding certain transition state searches.

Restarting Analytical Frequency jobs

Analytical frequency jobs can be restarted if a previous job did not finish. A restart can be done if the file TAPE21 has been saved without any corruption to the file. Upon doing a restart for analytical frequencies, the ADF program will check for the presence of an incomplete Hessian and/or for the presence of an incomplete U1 matrix, then attempt to figure out what more needs to be done.

The restart option may also be useful in combination with the atom selection option (i.e. by specifying a list of atoms following the keyword nuc in the analyticalfreq block key, see previous notes in this section). So for instance, you could calculate the partial Hessian for a subset of atoms of a molecule, and at a later time add another subset of atoms, or the rest of the atoms of the molecule to complete the Hessian.

Simplified interface for isotopic shift calculations in ADF

It is now possible to calculate an IR isotopic shift without modifying the atom type. The alternative mass is specified using the IsotopicShift keyword as follows:

IsotopicShift n1=m1 n2=m2 ...

n1=m1, n2=m2

Here, n1, n2, etc. are sequential numbers of atoms in the input geometry and m1, m2, etc. are their alternative masses.

All frequencies and intensities are calculated as usual with the original atomic masses taken from the corresponding fragment files and the isotopic shifts are calculated and printed afterwards. The Hessian is not broken down into irreducible representations and the rigid motions are not projected out prior to the normal mode analysis. Normal modes (eigenvectors) corresponding to the original and the shifted frequencies are matched against each other based on their overlap, which allows to calculate an isotopic shift for each frequency. This analysis may, however, fail for highly symmetric molecules with degenerate modes.

This feature currently works in conjunction with analytical frequencies only.

Numerical Frequencies¶

Input options

Calculation of the numerical frequencies is specified using the FREQUENCIES keyword in the GEOMETRY block. Most of the subkeys in the geometry block are meaningless for the calculation of frequencies. Indeed, a Frequencies calculation is not a variation on optimization, but rather a sequence of Single Point runs for the equilibrium geometry and a series of slightly different geometries. By comparison of the computed gradients the force constants and hence the frequencies are computed (in the harmonic approximation of the energy surface).

GEOMETRY
 Frequencies {Symm} {Allowed} {Numdif=Numdif} {Disrad=drad}
             {Disang=dang} {SCANALL} {NOSCAN}
 iterations Niter
end

Symm

This switch requests that frequencies are calculated in symmetric displacements. During such a calculation first symmetric atomic displacements are constructed. The number of such displacements in each irreducible representation corresponds to the number of frequencies with the corresponding symmetry. All displaced geometries within one representation have the same symmetry, which enables us to use it to speed up the computation significantly. This is a new option and for now it only works with geometries specified in Cartesian coordinates. This option does not work correctly with a restart file. This option does not work correctly when symmetry is explicitly specified in the input file.

Allowed

Another advantage of the symmetric displacements is that only a subset of frequencies can be calculated. The ALLOWED option requests computation of only IR-visible frequencies. This option is only useful for symmetric molecules where it can be a big time-saver.

Numdif

Must have the value between 1 and 4 and specifies the type of numerical differentiation that is applied to compute the force constants from gradients in slightly displaced geometries: 1-, 2-, 3-, or 4-point numerical differentiation. In the case of 1-point differentiation the gradients of the displaced geometry are compared with the gradients at the input (equilibrium) geometry. In 2-point case both a negative and a positive displacement are applied, yielding much more accurate results but at the expense of more computations. This option is the default. In certain cases the 3-point differentiation method gives better result than 2-point because it also takes gradients in the middle point into account. This is the case when geometry has not completely converged and the residual gradients are not quite close to zero. In this method, a formula is used that interpolates the second derivative (i.e. force constant) at the zero-force point. This way, the error due to a small deviation from the minimum geometry is decreased. The requirement is that the residual forces are small enough, more precisely, less that forces at displaced geometries (that is, using numdif=3 for arbitrary geometries is a bad idea). When Numdif=4 is specified, force constants matrix will be computed by making two displacements in each direction, the standard (see drad, dang below) and twice as short. The force constant is then computed using the Romberg formula that reduces the higher-order and noise components: H(tot)=(4*H(dx/2)-H(dx))/3. Although this method requires twice as many single-point evaluations, one can probably get reliable results using lower integration accuracy, which might be faster than the default.

dang and drad

The displacements of the coordinates that will be varied. Dang applies to angles (bond and dihedral) in degrees and drad applies to Cartesian (x, y, z) coordinates and to bond lengths, in angstrom. Defaults: 1 degree and 0.01 angstrom.

Niter

In a calculation of frequencies it is the total number of (displaced) geometries for which gradients are computed. By default this is internally determined such that the calculation of frequencies can be completed. If you reduce it, the run will only partially build the matrix of force constants and a restart is required to complete the computation.

WARNING: you cannot combine a Frequencies calculation with the QM/MM feature.

SCANALL and NOSCAN

ADF can scan some or all normal modes after a frequency calculation to verify the corresponding frequencies. By default, the normal modes corresponding all found imaginary frequencies are scanned. This can be switched off by specifying NOSCAN. Specifying SCANALL will tell ADF to scan along all normal modes. These options apply only to calculations in atomic displacements, that is when SYMM is not specified. These two options are mutually exclusive, the one specified last taking precedence.

Cartesian versus Z-matrix displacements

Cartesian displacements yield usually a higher accuracy than Z-matrix displacements because in the former case cancellation of numerical integration errors between the different geometries is (almost always) larger.

If Z-matrix coordinates are used as the displacement variables, then make sure that no bond angles of 180 (or zero) degrees are among them. They will very probably be treated incorrectly. If your molecule has such bond angles, use dummies to redefine the coordinates or use Cartesian displacements.

Frequencies and GEOVAR keyword

The use of the GEOVAR keyword in combination with a Frequencies run implies that constraints may be applied to the displacements, even if no coordinates are explicitly frozen. If different coordinates are connected to the same variable in the GEOVAR block, only combined displacements of the atoms will be allowed that correspond to a small change in the GEOVAR variable. For this reason, the combination of the GEOVAR and Frequencies keywords is to be handled with extreme caution. If no constraints are intended, it is recommended not to use the GEOVAR keyword, but to use DEFINE instead, or to specify the coordinates explicitly

Mobile Block Hessian (MBH)¶

The mobile block Hessian (MBH) method [282, 283] is useful when calculating vibrational frequencies of a small part of a very large system (molecule or cluster). Calculation of the full spectrum of such a system may be inefficient and is unnecessary if one is interested in one particular part. Besides, it may be difficult to extract normal modes related to the interesting sub-system out of the whole spectrum. Using MBH it is possible to treat parts of the system as rigid blocks. Each block will usually have only six frequencies related to its rigid motions compared to 3*N for when each atom of the block is treated separately.

The calculation of frequencies using mobile blocks is invoked by specifying FREQUENCIES and MBH keywords at the same time in the GEOMETRY input block:

GEOMETRY
  FREQUENCIES
  MBH blockname1 blockname2 ...
End

The names of blocks must correspond to the ones specified in the b= parameter in the ATOMS input block.

The second derivatives with respect to block motions are calculated by numerical differentiation. Since the number of degrees of freedom is reduced, the number of second derivatives is reduced as well. Therefore the MBH can realize a speed-up in the calculation of the Hessian compared to a full numerical frequency calculation.

MBH for partially optimized structures

MBH is suitable to calculate frequencies in partially optimized structures. Assume a geometry optimization is performed with the BLOCK subkey in the CONSTRAINTS section [see constrained geometry optimizations]. During the geometry optimization, the shape of the block is not changed. The internal geometry of the block is kept fixed, but the block as a whole can still translate or rotate.

At the end of such a partial geometry optimization, the position and orientation of the block are optimized. The total force on the block is zero. However, there might be still some residual forces within a block, since those degrees of freedom were not optimized.

A traditional frequency calculation performed on this partially optimized structure might result in non-physical imaginary frequencies without a clear interpretation. Therefore one should use an adapted formulation of normal mode analysis: the Mobile Block Hessian method. The MBH does not consider the internal degrees of freedom of the block (on which residual forces) apply, but instead uses the position/orientation of the block as coordinates. In the resulting normal mode eigenvectors, all atoms within the same block move collectively.

Of course, MBH can also be applied on a fully optimized structure.

Accuracy

The second derivatives with respect to Cartesian displacements (3 translations) of the free atoms (atoms not belonging to any block) and those with respect to block motions (3 block translation/3 block rotations) are calculated by numerical differentiation of the gradient. The accuracy of the second derivatives is mainly influenced by the accuracy of the gradient evaluation (e.g. accuracy numerical integration) and the step size in the numerical differentiation. The parameters DISRAD and DISANG can be specified to set the step size for Cartesian displacements (translations) and block rotations respectively. The step size for angles is automatically scaled with the block size.

FREQUENCIES {DISRAD=drad} {DISANG=dang}

The default values for drad and dang are the same as in the case of a standard numerical frequency run.

MBH Notes

Blocks consisting of 1 or 2 atoms are not supported and are ignored. This means that each atom of such a block is treated separately.

At this moment, it is not possible to calculate IR intensities with the MBH method.

The printed output of the MBH normal mode analysis lists all frequencies, including the ones corresponding to rotations and translation of the molecule as the whole. Note that while frequencies corresponding to translations should always be close to zero, the magnitude of the ones corresponding to rotations depends on how well the geometry is optimized.

Debug information can be obtained if one includes DEBUG MBHNormalModes.

Thermodynamics¶

At the end of a completed Frequencies calculation, a survey is given of thermodynamic properties: Heat Capacity, Internal Energy, Entropy. The computed results assume an ideal gas, and electronic contributions are ignored. The latter is a serious omission if the electronic configuration is (almost) degenerate, but the effect is small whenever the energy difference with the next state is large compared to the vibrational frequencies.

THERMO {P=pressure} {T=temp1 {temp2}} {nT=nT}

pressure

The Pressure in atmospheres. Default value: 1.0. A zero or negative pressure is adjusted by the program to a (very) small number 1e-4

temp1,temp2

The endpoints of the Temperature range (in K), for which the results are computed. By default only room temperature is considered (298.15 K). If the option T= is used and only one value supplied (temp1), then temp2 is assumed to be equal to temp1. A zero or negative temperature is adjusted by the program to a (very) small number 1e-4

nT

The number of steps by which the temperature interval is scanned. By default it is computed by the program from the temperature range (temp1, temp2), such that the step size is as close as possible to 10 K. Note that the number of temperatures for which the calculations are done is one more than the number of temperature steps.

The thermal analysis is based on the temperature dependent partition function. The energy of a (non-linear) molecule is (if the energy is measured from the zero-point energy)

The summation is over all harmonic \(\nu_j\), \(h\) is Planck’s constant and \(D\) is the dissociation energy

Contributions from low (less than 20 1/cm) frequencies to entropy, heat capacity and internal energy have always been excluded from the total values for thermodynamical quantities listed before. In ADF2013 they are also listed separately so the user can add them if they wish.

                                              hartree              eV         kcal/mol           kJ/mol
                                 --------------------     -----------       ----------      -----------

 Total Bonding Energy:             -0.744356253597793        -20.2550          -467.091         -1954.31


Zero-Point Energy :      0.033172 a.u.
===================      0.902646 eV

     Temp                                                       Transl     Rotat    Vibrat     Total
     ----                                                       ------     -----    ------     -----
     298.15   Entropy (cal/mole-K):                             34.441    11.494     0.126    46.061
              Internal Energy (Kcal/mole):                       0.889     0.889    20.847    22.624
              Constant Volume Heat Capacity (cal/mole-K):        2.981     2.981     0.533     6.495

Accuracy¶

Accuracy is a crucial aspect in the computation of frequencies, in particular for modes with low frequencies: the gradients at the geometries displaced along that mode will hardly change - analytically - from their equilibrium values, so numerical integration noise may easily affect the reliability of the computed differences in gradients. It is worthwhile to consider carefully the size of the displacements. At one hand they should be small in order to suppress the effect of higher order (anharmonic) terms in the energy surface around the minimum, at the other hand they should be large enough to get significant differences in gradients so that these are computed reliably.

High precision calculations where low frequency modes are involved may require high integration settings [14]. It is sometimes advisable to increase the numerical accuracy of the calculation:

NumericalQuality {Good|VeryGood}

Using 2-point differentiation rather than 1-point differentiation implies two-sided displacements of the atoms. This doubles the computational effort but in the so-computed force constants all anharmonic terms of odd order are eliminated. Since in general the lowest anharmonicity is third order this eliminates the first anharmonicity. Again, this is a feature directed primarily at obtaining highly accurate and reliable results.

The 3- and 4-point methods are intended to assist in special cases and as an extra check when the results obtained with the 2-point formula are not satisfactory. The 3-point formula should be used when residual forces after geometry optimization are between 0.01 and 0.0001 a.u./angstrom. In this case frequencies obtained with the 3-point formula are much closer to those that would be computed at the exact optimum geometry.

If a Frequencies calculation is carried out only to construct a good start-up Hessian for a TS search (see the restart key), accurate results are not crucial. The most important thing in such a situation is to get a fair guess for the negative eigenvalue and its associated mode, and to avoid spurious additional negative eigenvalues. We recommend to avoid the rather time-consuming standard Hessian-computing preparation run for a TS search and to lower the precision of the Frequencies run. A reasonable value should be 4.0.

Isotope Shifts of Vibrational Frequencies¶

To calculate isotopic shifts using ADF do the following:

• Calculate frequencies and save TAPE21 with a different name, say result.t21
• Modify the input file as follows:
  • add “RESTART result.t21” anywhere in the input file
  • create new fragment file with different mass
  • specify the fragment file in FRAGMENTS section
• Run ADF with the new input

Alternatively one can use the the key ATOMPROPS to change the masses of the atoms.

Please note that if you change the fragment file for an atom that has symmetry-equivalent ones then the new fragment file will be applied to all of the atoms. Example: first calculate the NH3 frequencies in the C(3v) symmetry and then change H to D. This will mean that one calculates the frequencies of a ND3 molecule and not of NH2D as one might want to do. If one wants to calculate the frequencies of NH2D one first has to do a calculation with lower symmetry, say C(s), to be able to change isotope of only one of the hydrogens.

Scanning a Range of Frequencies¶

In ADF2006 it was already possible to request a full scan of all frequencies obtained by finite differences. This was originally done to help identify spurious imaginary frequencies that sometimes appear where one would expect a very low (nearly zero) frequency. Most frequently this happens when there is a barrier-free rotation of, for example, methyl groups.

Starting from ADF2007.01, it is possible to scan any range of frequencies calculated in the same run or found in a restart file. Note that one should not use the key SCANFREQ in combination with the key SYMMETRY. The input keyword used to request the scan is as follows:

SCANFREQ low high {NUM=num DISRAD=disrad}

low, high

Two values defining an interval of frequencies to scan. Frequencies that fall within the interval will be recalculated by numerical differentiation of the gradient along the respective frequency’s normal mode. This means that 2**N* single-point calculations with gradients will be performed, where N is the number of frequencies within the range. Imaginary frequencies are specified using negative values, which is consistent with the notation adopted within ADF.

num

Num is an integer number specifying how many points are to be used for numerical differentiation: 2, 4, or 6. The default value is 2.

disrad

Disrad specifies the step size (in Angstrom) to be made in each direction for 2-point differentiation. For the 4-point differentiation the maximum deviation from the equilibrium geometry will be twice as large as for the 2-point one and so on. The default value for disrad is the same as used for numerical frequencies.

The main advantage of this method is that single-point calculations used to obtain a force constant are performed within the same symmetry and, usually, with the same numerical integration grid, which significantly reduces the level of numerical noise and thus increases accuracy of the calculated frequency.

Moments of inertia¶

In case frequencies are calculated in ADF, ADF also reports the moments of inertia of the molecule in units of amu bohr² (amu = atomic mass unit).